CSCCondDB/src/CSCAFEBThrAnalysis.cc

0001 #include "OnlineDB/CSCCondDB/interface/CSCAFEBThrAnalysis.h"
0002 #include "OnlineDB/CSCCondDB/interface/CSCToAFEB.h"
0003 #include <OnlineDB/CSCCondDB/interface/CSCFitAFEBThr.h>
0004 #include "OnlineDB/CSCCondDB/interface/CSCOnlineDB.h"
0005 #include "CondFormats/CSCObjects/interface/CSCobject.h"
0006 #include "TMath.h"
0007
0008 class CSCFitAFEBThr;
0009
0010 CSCAFEBThrAnalysis::CSCAFEBThrAnalysis() {
0011   hist_file = nullptr;  // set to null
0012
0013   nmbev = 0;
0014   nmbev_no_wire = 0;
0015
0016   npulses = 100;
0017   vecDac.clear();
0018   BegDac = 1;
0019   EndDac = 29;
0020   EvDac = 1;
0021   StepDac = 1;
0022   indDac = 0;
0023
0024   vecDacOccup.assign(150, 0);
0025   m_wire_dac.clear();
0026   m_res_for_db.clear();
0027   mh_ChanEff.clear();
0028   mh_FirstTime.clear();
0029   mh_AfebDac.clear();
0030 }
0031
0032 void CSCAFEBThrAnalysis::setup(const std::string& histoFileName) {
0033   /// open the histogram file
0034   hist_file = new TFile(histoFileName.c_str(), "RECREATE");
0035   hist_file->cd();
0036 }
0037
0038 void CSCAFEBThrAnalysis::bookForId(int flag, const int& idint, const std::string& idstring) {
0039   hist_file->cd();
0040
0041   std::ostringstream ss;
0042
0043   if (flag == 100) {
0044     ss << idint << "_Anode_First_Time";
0045     mh_FirstTime[idint] = new TH2F(ss.str().c_str(), "", 675, 0.0, 675.0, 50, 0.0, 10.0);
0046     mh_FirstTime[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0047     mh_FirstTime[idint]->GetYaxis()->SetTitle("Anode First Time Bin");
0048     mh_FirstTime[idint]->SetOption("BOX");
0049     ss.str("");  // clear
0050   }
0051
0052   if (flag == 101) {
0053     ss << idint << "_Anode_Chan_Eff";
0054     mh_ChanEff[idint] = new TH1F(ss.str().c_str(), "", 675, 0.0, 675.0);
0055     mh_ChanEff[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0056     mh_ChanEff[idint]->GetYaxis()->SetTitle("Entries");
0057     ss.str("");  // clear
0058   }
0059
0060   if (flag == 200) {
0061     ss << idint << "_Anode_AfebDac";
0062     mh_AfebDac[idint] = new TH2F(ss.str().c_str(), "", 75, 0.0, 75.0, 50, 0.0, 50.0);
0063     mh_AfebDac[idint]->GetXaxis()->SetTitle("Threshold DAC");
0064     mh_AfebDac[idint]->GetYaxis()->SetTitle("AFEB Channel Occupancy");
0065     mh_AfebDac[idint]->SetOption("COL");
0066     ss.str("");  // clear
0067   }
0068
0069   if (flag == 300) {
0070     ss << idint << "_Anode_AfebThrPar";
0071     mh_AfebThrPar[idint] = new TH2F(ss.str().c_str(), "", 700, 0.0, 700.0, 50, 0.0, 50.0);
0072     mh_AfebThrPar[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0073     mh_AfebThrPar[idint]->GetYaxis()->SetTitle("AFEB Channel Threshold (DAC)");
0074     mh_AfebThrPar[idint]->SetOption("BOX");
0075     ss.str("");  // clear
0076   }
0077
0078   if (flag == 400) {
0079     ss << idint << "_Anode_AfebNoisePar";
0080     mh_AfebNoisePar[idint] = new TH2F(ss.str().c_str(), "", 700, 0.0, 700.0, 50, 0.0, 5.0);
0081     mh_AfebNoisePar[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0082     mh_AfebNoisePar[idint]->GetYaxis()->SetTitle("AFEB Channel Noise (DAC)");
0083     mh_AfebNoisePar[idint]->SetOption("BOX");
0084     ss.str("");  // clear
0085   }
0086
0087   if (flag == 500) {
0088     ss << idint << "_Anode_AfebNDF";
0089     mh_AfebNDF[idint] = new TH2F(ss.str().c_str(), "", 700, 0.0, 700.0, 25, -5.0, 20.0);
0090     mh_AfebNDF[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0091     mh_AfebNDF[idint]->GetYaxis()->SetTitle("AFEB Channel Fit NDF");
0092     mh_AfebNDF[idint]->SetOption("BOX");
0093     ss.str("");  // clear
0094   }
0095
0096   if (flag == 600) {
0097     ss << idint << "_Anode_AfebChi2perNDF";
0098     mh_AfebChi2perNDF[idint] = new TH2F(ss.str().c_str(), "", 700, 0.0, 700.0, 50, 0.0, 10.0);
0099     mh_AfebChi2perNDF[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
0100     mh_AfebChi2perNDF[idint]->GetYaxis()->SetTitle("AFEB Channel Fit Chi2/NDF");
0101     mh_AfebChi2perNDF[idint]->SetOption("BOX");
0102     ss.str("");  // clear
0103   }
0104 }
0105
0106 void CSCAFEBThrAnalysis::hf1ForId(std::map<int, TH1*>& mp, int flag, const int& id, float& x, float w) {
0107   std::map<int, TH1*>::iterator h;
0108   h = mp.find(id);
0109   if (h == mp.end()) {
0110     bookForId(flag, id, "");
0111     h = mp.find(id);
0112   }
0113   h->second->Fill(x, w);
0114 }
0115
0116 void CSCAFEBThrAnalysis::hf2ForId(std::map<int, TH2*>& mp, int flag, const int& id, float& x, float& y, float w) {
0117   std::map<int, TH2*>::iterator h;
0118   h = mp.find(id);
0119   if (h == mp.end()) {
0120     bookForId(flag, id, "");
0121     h = mp.find(id);
0122   }
0123   h->second->Fill(x, y, w);
0124 }
0125
0126 /* Analyze the hits */
0127 void CSCAFEBThrAnalysis::analyze(const CSCWireDigiCollection& wirecltn) {
0128   std::ostringstream ss;
0129   std::map<int, std::vector<int> >::iterator intIt;
0130   std::map<int, std::vector<std::vector<int> > >::iterator wiredacIt;
0131
0132   std::vector<int> vec;
0133   int afeb, ch;
0134   float x, y;
0135
0136   m_wire_ev.clear();
0137
0138   /// Store DAC
0139
0140   nmbev++;
0141   //indDac=(nmbev-1)/EvDac;
0142   float dac = BegDac + StepDac * indDac;
0143   if (vecDac.size() <= indDac)
0144     vecDac.push_back(dac);
0145
0146   //std::cout<<"  Event "<<nmbev;
0147   //std::cout<<"  "<<indDac<<" "<<vecDac[indDac]<<std::endl;
0148
0149   //Anode wires
0150
0151   CSCWireDigiCollection::DigiRangeIterator wiredetUnitIt;
0152   if (wirecltn.begin() == wirecltn.end())
0153     nmbev_no_wire++;
0154
0155   if (wirecltn.begin() != wirecltn.end()) {
0156     vecDacOccup[indDac] = vecDacOccup[indDac] + 1;
0157
0158     for (wiredetUnitIt = wirecltn.begin(); wiredetUnitIt != wirecltn.end(); ++wiredetUnitIt) {
0159       const CSCDetId& id = (*wiredetUnitIt).first;
0160
0161       const int idchamber = id.endcap() * 10000 + id.station() * 1000 + id.ring() * 100 + id.chamber();
0162       const int idlayer =
0163           id.endcap() * 100000 + id.station() * 10000 + id.ring() * 1000 + id.chamber() * 10 + id.layer();
0164
0165       //    std::cout<<idchamber<<" "<<idlayer<<std::endl;
0166
0167       const int maxwire = csctoafeb.getMaxWire(id.station(), id.ring());
0168       std::vector<int> wireplane(maxwire, 0);
0169
0170       const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
0171       for (CSCWireDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0172         int iwire = (*digiIt).getWireGroup();
0173         if (iwire <= maxwire) {
0174           if (wireplane[iwire - 1] == 0) {
0175             wireplane[iwire - 1] = (*digiIt).getBeamCrossingTag() + 1;
0176             ch = csctoafeb.getAfebCh(id.layer(), (*digiIt).getWireGroup());
0177             afeb = csctoafeb.getAfebPos(id.layer(), (*digiIt).getWireGroup());
0178
0179             /// Plot time bin of the first hit vs AFEB channels
0180
0181             x = (afeb - 1) * 16 + ch;
0182             y = wireplane[iwire - 1];
0183             hf2ForId(mh_FirstTime, 100, idchamber, x, y, 1.0);
0184
0185           }  // end if wireplane[iwire-1]==0
0186         }  // end if iwire<=csctoafeb.getMaxWire(id.station(),id.ring()
0187       }  // end of for digis in layer
0188
0189       /// Store time bin of the first hit into map
0190
0191       if (m_wire_ev.find(idlayer) == m_wire_ev.end())
0192         m_wire_ev[idlayer] = wireplane;
0193
0194     }  // end of cycle on detUnit
0195
0196     /// Accumulate hits into map of wires vs DAC
0197
0198     for (intIt = m_wire_ev.begin(); intIt != m_wire_ev.end(); ++intIt) {
0199       const int idwirev = (*intIt).first;
0200       const std::vector<int> wiretemp = (*intIt).second;  // What for?
0201       int nsize = EndDac - BegDac + 1;
0202       std::vector<int> zer(nsize, 0);
0203
0204       wiredacIt = m_wire_dac.find(idwirev);
0205       if (wiredacIt == m_wire_dac.end()) {
0206         for (unsigned int j = 0; j < wiretemp.size(); j++)
0207           m_wire_dac[idwirev].push_back(zer);
0208         wiredacIt = m_wire_dac.find(idwirev);
0209         //      std::cout<<idwirev<<"  "<<wiredacIt->second.size()<<"  "<<
0210         //                 wiredacIt->second[0].size()<<std::endl;
0211       }
0212       for (unsigned int i = 0; i < (*intIt).second.size(); i++)
0213         if ((*intIt).second[i] > 0)
0214           wiredacIt->second[i][indDac] = wiredacIt->second[i][indDac] + 1;
0215     }  // end of adding hits to the map of wire/DAC
0216   }  // end of if wire collection not empty
0217   indDac++;
0218   if (dac == (float)EndDac)
0219     indDac = 0;
0220 }  // end of void CSCAFEBThrAnalysis
0221
0222 void CSCAFEBThrAnalysis::done() {
0223   float x, y;
0224
0225   //This is for DB transfer
0226   CSCobject cn{};
0227
0228   std::map<int, std::vector<std::vector<int> > >::iterator mwiredacIt;
0229   std::map<int, std::vector<std::vector<float> > >::iterator mresfordbIt;
0230   std::vector<int>::iterator vecIt;
0231
0232   std::cout << "Events analyzed  " << nmbev << std::endl;
0233   std::cout << "Events no anodes " << nmbev_no_wire << std::endl << std::endl;
0234
0235   std::cout << "DAC occupancy" << std::endl;
0236   size_t ndacsize = EndDac - BegDac + 1;
0237   for (size_t i = 0; i < ndacsize; i++)
0238     std::cout << " " << vecDacOccup[i];
0239   std::cout << "\n\n" << std::endl;
0240
0241   std::vector<float> inputx;
0242   std::vector<float> inputy;
0243
0244   std::vector<float> mypar(2, 0.0);
0245   std::vector<float> ermypar(2, 0.0);
0246   float ercorr, chisq, edm;
0247   int ndf, niter;
0248
0249   int ch, afeb, idchmb;
0250
0251   CSCFitAFEBThr* fitAnodeThr;
0252
0253   std::vector<float> fitres(4, 0);
0254
0255   //  std::cout<<"Threshold curves:\n"<<std::endl;
0256
0257   for (mwiredacIt = m_wire_dac.begin(); mwiredacIt != m_wire_dac.end(); ++mwiredacIt) {
0258     int idwiredac = (*mwiredacIt).first;
0259
0260     int layer = idwiredac - (idwiredac / 10) * 10;
0261     idchmb = idwiredac / 10;
0262
0263     for (int unsigned iwire = 0; iwire < mwiredacIt->second.size(); iwire++) {
0264       afeb = csctoafeb.getAfebPos(layer, iwire + 1);
0265       ch = csctoafeb.getAfebCh(layer, iwire + 1);
0266       int afebid = (idwiredac / 10) * 100 + csctoafeb.getAfebPos(layer, iwire + 1);
0267
0268       indDac = 0;
0269       for (vecIt = mwiredacIt->second[iwire].begin(); vecIt != mwiredacIt->second[iwire].end(); ++vecIt) {
0270         x = vecDac[indDac];
0271         y = *vecIt;
0272         hf2ForId(mh_AfebDac, 200, afebid, x, y, 1.0);
0273
0274         ///  Plot "efficiency"  vs  AFEB channels
0275         if (indDac == 0) {
0276           x = (afeb - 1) * 16 + ch;
0277           y = *vecIt;
0278           hf1ForId(mh_ChanEff, 101, idchmb, x, y);
0279         }
0280
0281         inputx.push_back(x);
0282         inputy.push_back(y);
0283
0284         indDac++;
0285       }
0286       // end of DAC cycle to form vectors of input data (inputx,inputy)for fit
0287
0288       //           std::cout<<afebid<<" "<<ch<<std::endl;
0289       //           for(unsigned int i=0;i<inputx.size();i++)
0290       //              std::cout<<" "<<inputy[i];
0291       //           std::cout<<std::endl;
0292
0293       for (unsigned int i = 0; i < 2; i++) {
0294         mypar[i] = 0.0;
0295         ermypar[i] = 0.0;
0296       }
0297       ercorr = 0.0;
0298       chisq = 0.0;
0299       ndf = 0;
0300       niter = 0;
0301       edm = 0.0;
0302
0303       /// Fitting threshold curve
0304       fitAnodeThr = new CSCFitAFEBThr();
0305
0306       fitAnodeThr->ThresholdNoise(inputx, inputy, npulses, vecDacOccup, mypar, ermypar, ercorr, chisq, ndf, niter, edm);
0307       delete fitAnodeThr;
0308
0309       //  std::cout<<"Fit "<<mypar[0]<<" "<<mypar[1]<<" "<<ndf<<" "<<chisq
0310       //           <<std::endl<<std::endl;
0311
0312       /// Histogram fit results for given CSC vs wire defined as x=(afeb-1)*16+ch
0313
0314       x = (afeb - 1) * 16 + ch;
0315
0316       /// Threshold
0317       y = mypar[0];
0318       hf2ForId(mh_AfebThrPar, 300, idchmb, x, y, 1.0);
0319
0320       /// Noise
0321       y = mypar[1];
0322       hf2ForId(mh_AfebNoisePar, 400, idchmb, x, y, 1.0);
0323
0324       /// NDF
0325       y = ndf;
0326       hf2ForId(mh_AfebNDF, 500, idchmb, x, y, 1.0);
0327
0328       /// Chi2/NDF
0329       y = 0.0;
0330       if (ndf > 0)
0331         y = chisq / (float)ndf;
0332       hf2ForId(mh_AfebChi2perNDF, 600, idchmb, x, y, 1.0);
0333
0334       /// Prepare vector of fit results
0335       fitres[0] = mypar[0];
0336       fitres[1] = mypar[1];
0337       fitres[2] = ndf;
0338       fitres[3] = 0.0;
0339       if (ndf > 0)
0340         fitres[3] = chisq / (float)ndf;
0341
0342       /// Store fit results to map of wire vectors of vectors of results
0343
0344       mresfordbIt = m_res_for_db.find(idwiredac);
0345       if (mresfordbIt == m_res_for_db.end())
0346         m_res_for_db[idwiredac].push_back(fitres);
0347       else
0348         m_res_for_db[idwiredac].push_back(fitres);
0349
0350       inputx.clear();
0351       inputy.clear();
0352
0353     }  // end for(int iwire=0)
0354   }  // end of iteration thru m_wire_dac map
0355
0356   std::cout << "Size of map for DB " << m_res_for_db.size() << std::endl;
0357
0358   std::cout << "The following CSCs will go to DB" << std::endl << std::endl;
0359
0360   for (mresfordbIt = m_res_for_db.begin(); mresfordbIt != m_res_for_db.end(); ++mresfordbIt) {
0361     int idlayer = (*mresfordbIt).first;
0362     idchmb = idlayer / 10;
0363     int layer = idlayer - idchmb * 10;
0364     std::cout << "CSC " << idchmb << "  Layer " << layer << "  " << (*mresfordbIt).second.size() << std::endl;
0365   }
0366
0367   for (mresfordbIt = m_res_for_db.begin(); mresfordbIt != m_res_for_db.end(); ++mresfordbIt) {
0368     int idlayer = (*mresfordbIt).first;
0369
0370     //This is for DB transfer
0371     int size = (*mresfordbIt).second.size();
0372     cn.obj[idlayer].resize(size);
0373
0374     for (unsigned int i = 0; i < (*mresfordbIt).second.size(); i++) {
0375       std::cout << idlayer << " " << i + 1 << "    ";
0376
0377       for (int j = 0; j < 4; j++)
0378         std::cout << (*mresfordbIt).second[i][j] << " ";
0379       std::cout << std::endl;
0380
0381       //This is for DB transfer
0382       cn.obj[idlayer][i].resize(4);
0383       cn.obj[idlayer][i][0] = (*mresfordbIt).second[i][0];
0384       cn.obj[idlayer][i][1] = (*mresfordbIt).second[i][1];
0385       cn.obj[idlayer][i][2] = (*mresfordbIt).second[i][2];
0386       cn.obj[idlayer][i][3] = (*mresfordbIt).second[i][3];
0387     }
0388   }
0389
0390   //send data to DB
0391   //dbon->cdbon_last_run("afeb_thresholds",&run);
0392   //std::cout<<"Last AFEB thresholds run "<<run<<" for run file "<<myname<<" saved "<<myTime<<std::endl;
0393   //if(debug) dbon->cdbon_write(cn,"afeb_thresholds",run+1,myTime);
0394
0395   if (hist_file != nullptr) {  // if there was a histogram file...
0396     hist_file->Write();        // write out the histrograms
0397     delete hist_file;          // close and delete the file
0398     hist_file = nullptr;       // set to zero to clean up
0399     std::cout << "Hist. file was closed\n";
0400   }
0401
0402   std::cout << " End of CSCAFEBThrAnalysis" << std::endl;
0403 }