EcalLaserAnalyzer/plugins/EcalABAnalyzer.cc

0001 /*
0002  *  \class EcalABAnalyzer
0003  *
0004  *  primary author: Julie Malcles - CEA/Saclay
0005  *  author: Gautier Hamel De Monchenault - CEA/Saclay
0006  */
0007
0008 #include <TAxis.h>
0009 #include <TH1.h>
0010 #include <TProfile.h>
0011 #include <TTree.h>
0012 #include <TChain.h>
0013 #include <TFile.h>
0014 #include <TMath.h>
0015
0016 #include "CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalABAnalyzer.h"
0017
0018 #include <sstream>
0019 #include <fstream>
0020 #include <iomanip>
0021
0022 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0023 #include <FWCore/Framework/interface/EventSetup.h>
0024 #include <FWCore/Framework/interface/Event.h>
0025 #include <FWCore/Framework/interface/MakerMacros.h>
0026 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0027 #include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
0028 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
0029
0030 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h>
0031 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
0032 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h>
0033 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPDPulse.h>
0034 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h>
0035 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h>
0036 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h>
0037
0038 using namespace std;
0039
0040 //========================================================================
0041 EcalABAnalyzer::EcalABAnalyzer(const edm::ParameterSet& iConfig)
0042     //========================================================================
0043     : iEvent(0),
0044       eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0045       eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0046       digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0047       digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0048       rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0049       mappingToken_(esConsumes()),
0050       // framework parameters with default values
0051       _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
0052       _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 2)),
0053       _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
0054       _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
0055       _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 2)),
0056       _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 9)),
0057       _timingquallow(iConfig.getUntrackedParameter<unsigned int>("timingQualLow", 3)),
0058       _timingqualhigh(iConfig.getUntrackedParameter<unsigned int>("timingQualHigh", 8)),
0059       _ratiomincutlow(iConfig.getUntrackedParameter<double>("ratioMinCutLow", 0.4)),
0060       _ratiomincuthigh(iConfig.getUntrackedParameter<double>("ratioMinCutHigh", 0.95)),
0061       _ratiomaxcutlow(iConfig.getUntrackedParameter<double>("ratioMaxCutLow", 0.8)),
0062       _presamplecut(iConfig.getUntrackedParameter<double>("presampleCut", 5.0)),
0063       _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0064       _alpha(iConfig.getUntrackedParameter<double>("alpha", 1.5076494)),
0065       _beta(iConfig.getUntrackedParameter<double>("beta", 1.5136036)),
0066       _nevtmax(iConfig.getUntrackedParameter<unsigned int>("nEvtMax", 200)),
0067       _noise(iConfig.getUntrackedParameter<double>("noise", 2.0)),
0068       _chi2cut(iConfig.getUntrackedParameter<double>("chi2cut", 100.0)),
0069       _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0070       _fedid(iConfig.getUntrackedParameter<int>("fedId", -999)),
0071       _qualpercent(iConfig.getUntrackedParameter<double>("qualPercent", 0.2)),
0072       _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
0073       resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0074       nCrys(NCRYSEB),
0075       runType(-1),
0076       runNum(0),
0077       fedID(-1),
0078       dccID(-1),
0079       side(2),
0080       lightside(2),
0081       iZ(1),
0082       phi(-1),
0083       eta(-1),
0084       event(0),
0085       color(-1),
0086       channelIteratorEE(0)
0087
0088 //========================================================================
0089
0090 {
0091   if (_ecalPart == "EB") {
0092     ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0093   } else if (_ecalPart == "EE") {
0094     eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0095   }
0096
0097   // Geometrical constants initialization
0098
0099   if (_ecalPart == "EB") {
0100     nCrys = NCRYSEB;
0101   } else {
0102     nCrys = NCRYSEE;
0103   }
0104   iZ = 1;
0105   if (_fedid <= 609)
0106     iZ = -1;
0107
0108   for (unsigned int j = 0; j < nCrys; j++) {
0109     iEta[j] = -1;
0110     iPhi[j] = -1;
0111     iTowerID[j] = -1;
0112     iChannelID[j] = -1;
0113     idccID[j] = -1;
0114     iside[j] = -1;
0115     wasTimingOK[j] = true;
0116     wasGainOK[j] = true;
0117     nevtAB[j] = 0;
0118   }
0119
0120   // Quality check flags
0121
0122   isGainOK = true;
0123   isTimingOK = true;
0124
0125   // Objects dealing with pulses
0126
0127   APDPulse = new TAPDPulse(_nsamples,
0128                            _presample,
0129                            _firstsample,
0130                            _lastsample,
0131                            _timingcutlow,
0132                            _timingcuthigh,
0133                            _timingquallow,
0134                            _timingqualhigh,
0135                            _ratiomincutlow,
0136                            _ratiomincuthigh,
0137                            _ratiomaxcutlow);
0138
0139   // Objects needed for npresample calculation
0140
0141   Delta01 = new TMom();
0142   Delta12 = new TMom();
0143   _fitab = true;
0144 }
0145
0146 //========================================================================
0147 EcalABAnalyzer::~EcalABAnalyzer() {
0148   //========================================================================
0149
0150   // do anything here that needs to be done at desctruction time
0151   // (e.g. close files, deallocate resources etc.)
0152 }
0153
0154 //========================================================================
0155 void EcalABAnalyzer::beginJob() {
0156   //========================================================================
0157
0158   //Calculate alpha and beta
0159
0160   // Define output results filenames and shape analyzer object (alpha,beta)
0161   //=====================================================================
0162
0163   // 1) AlphaBeta files
0164
0165   doesABTreeExist = true;
0166
0167   std::stringstream nameabinitfile;
0168   nameabinitfile << resdir_ << "/ABInit.root";
0169   alphainitfile = nameabinitfile.str();
0170
0171   std::stringstream nameabfile;
0172   std::stringstream link;
0173   nameabfile << resdir_ << "/AB.root";
0174
0175   FILE* test;
0176   test = fopen(nameabinitfile.str().c_str(), "r");
0177   if (test == nullptr) {
0178     doesABTreeExist = false;
0179     _fitab = true;
0180   } else {
0181     fclose(test);
0182   }
0183
0184   TFile* fAB = nullptr;
0185   TTree* ABInit = nullptr;
0186   if (doesABTreeExist) {
0187     fAB = new TFile(nameabinitfile.str().c_str());
0188   }
0189   if (doesABTreeExist && fAB) {
0190     ABInit = (TTree*)fAB->Get("ABCol0");
0191   }
0192
0193   // 2) Shape analyzer
0194
0195   if (doesABTreeExist && fAB && ABInit && ABInit->GetEntries() != 0) {
0196     shapana = new TShapeAnalysis(ABInit, _alpha, _beta, 5.5, 1.0);
0197     doesABTreeExist = true;
0198   } else {
0199     shapana = new TShapeAnalysis(_alpha, _beta, 5.5, 1.0);
0200     doesABTreeExist = false;
0201     _fitab = true;
0202   }
0203   shapana->set_const(_nsamples, _firstsample, _lastsample, _presample, _nevtmax, _noise, _chi2cut);
0204
0205   if (doesABTreeExist && fAB)
0206     fAB->Close();
0207
0208   if (_fitab) {
0209     alphafile = nameabfile.str();
0210   } else {
0211     alphafile = alphainitfile;
0212     link << "ln -s " << resdir_ << "/ABInit.root " << resdir_ << "/AB.root";
0213     system(link.str().c_str());
0214   }
0215
0216   // Define output results files' names
0217
0218   std::stringstream namefile;
0219   namefile << resdir_ << "/AB.root";
0220   alphafile = namefile.str();
0221 }
0222
0223 //========================================================================
0224 void EcalABAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0225   //========================================================================
0226
0227   ++iEvent;
0228
0229   // retrieving DCC header
0230   edm::Handle<EcalRawDataCollection> pDCCHeader;
0231   const EcalRawDataCollection* DCCHeader = nullptr;
0232   e.getByToken(rawDataToken_, pDCCHeader);
0233   if (!pDCCHeader.isValid()) {
0234     edm::LogError("nodata") << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str()
0235                             << " " << eventHeaderProducer_.c_str();
0236   } else {
0237     DCCHeader = pDCCHeader.product();
0238   }
0239
0240   //retrieving crystal data from Event
0241   edm::Handle<EBDigiCollection> pEBDigi;
0242   const EBDigiCollection* EBDigi = nullptr;
0243   edm::Handle<EEDigiCollection> pEEDigi;
0244   const EEDigiCollection* EEDigi = nullptr;
0245   if (_ecalPart == "EB") {
0246     e.getByToken(ebDigiToken_, pEBDigi);
0247     if (!pEBDigi.isValid()) {
0248       edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0249     } else {
0250       EBDigi = pEBDigi.product();
0251     }
0252   } else if (_ecalPart == "EE") {
0253     e.getByToken(eeDigiToken_, pEEDigi);
0254     if (!pEEDigi.isValid()) {
0255       edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0256     } else {
0257       EEDigi = pEEDigi.product();
0258     }
0259   } else {
0260     edm::LogError("cfg_error") << " Wrong ecalPart in cfg file";
0261     return;
0262   }
0263
0264   // retrieving electronics mapping
0265   const auto& TheMapping = c.getData(mappingToken_);
0266
0267   // =============================
0268   // Decode DCCHeader Information
0269   // =============================
0270
0271   for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0272        ++headerItr) {
0273     // Get run type and run number
0274
0275     int fed = headerItr->fedId();
0276     if (fed != _fedid && _fedid != -999)
0277       continue;
0278
0279     runType = headerItr->getRunType();
0280     runNum = headerItr->getRunNumber();
0281     event = headerItr->getLV1();
0282
0283     dccID = headerItr->getDccInTCCCommand();
0284     fedID = headerItr->fedId();
0285     lightside = headerItr->getRtHalf();
0286
0287     // Check fed corresponds to the DCC in TCC
0288
0289     if (600 + dccID != fedID)
0290       continue;
0291
0292     // Cut on runType
0293
0294     if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
0295         runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
0296       return;
0297
0298     // Retrieve laser color and event number
0299
0300     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0301     color = settings.wavelength;
0302     if (color < 0)
0303       return;
0304
0305     std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0306     if (iter == colors.end()) {
0307       colors.push_back(color);
0308     }
0309   }
0310
0311   // Cut on fedID
0312
0313   if (fedID != _fedid && _fedid != -999)
0314     return;
0315
0316   // ===========================
0317   // Decode EBDigis Information
0318   // ===========================
0319
0320   int adcGain = 0;
0321
0322   if (EBDigi) {
0323     // Loop on crystals
0324     //===================
0325
0326     for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
0327          ++digiItr) {  // Loop on crystals
0328
0329       // Retrieve geometry
0330       //===================
0331
0332       EBDetId id_crystal(digiItr->id());
0333       EBDataFrame df(*digiItr);
0334
0335       int etaG = id_crystal.ieta();  // global
0336       int phiG = id_crystal.iphi();  // global
0337
0338       int etaL;  // local
0339       int phiL;  // local
0340       std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0341
0342       etaL = LocalCoord.first;
0343       phiL = LocalCoord.second;
0344
0345       eta = etaG;
0346       phi = phiG;
0347
0348       side = MEEBGeom::side(etaG, phiG);
0349
0350       // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
0351
0352       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0353
0354       int towerID = elecid_crystal.towerId();
0355       int strip = elecid_crystal.stripId();
0356       int xtal = elecid_crystal.xtalId();
0357       int channelID = 5 * (strip - 1) + xtal - 1;
0358
0359       unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0360
0361       assert(channel < nCrys);
0362
0363       iEta[channel] = eta;
0364       iPhi[channel] = phi;
0365       iTowerID[channel] = towerID;
0366       iChannelID[channel] = channelID;
0367       idccID[channel] = dccID;
0368       iside[channel] = side;
0369
0370       // APD Pulse
0371       //===========
0372
0373       // Loop on adc samples
0374
0375       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0376         EcalMGPASample samp_crystal(df.sample(i));
0377         adc[i] = samp_crystal.adc();
0378         adcG[i] = samp_crystal.gainId();
0379         adc[i] *= adcG[i];
0380         if (i == 0)
0381           adcGain = adcG[i];
0382         if (i > 0)
0383           adcGain = TMath::Max(adcG[i], adcGain);
0384       }
0385
0386       APDPulse->setPulse(adc);
0387
0388       // Quality checks
0389       //================
0390
0391       if (adcGain != 1)
0392         nEvtBadGain[channel]++;
0393       if (!APDPulse->isTimingQualOK())
0394         nEvtBadTiming[channel]++;
0395       nEvtTot[channel]++;
0396
0397       // Fill if Pulse is fine
0398       //=======================
0399
0400       if (APDPulse->isPulseOK() && lightside == side) {
0401         Delta01->addEntry(APDPulse->getDelta(0, 1));
0402         Delta12->addEntry(APDPulse->getDelta(1, 2));
0403
0404         if (nevtAB[channel] < _nevtmax && _fitab) {
0405           if (doesABTreeExist)
0406             shapana->putAllVals(channel, adc, eta, phi);
0407           else
0408             shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0409           nevtAB[channel]++;
0410         }
0411       }
0412     }
0413
0414   } else if (EEDigi) {
0415     // Loop on crystals
0416     //===================
0417
0418     for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end();
0419          ++digiItr) {  // Loop on crystals
0420
0421       // Retrieve geometry
0422       //===================
0423
0424       EEDetId id_crystal(digiItr->id());
0425       EEDataFrame df(*digiItr);
0426
0427       phi = id_crystal.ix();
0428       eta = id_crystal.iy();
0429
0430       int iX = (phi - 1) / 5 + 1;
0431       int iY = (eta - 1) / 5 + 1;
0432
0433       side = MEEEGeom::side(iX, iY, iZ);
0434       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0435
0436       int towerID = elecid_crystal.towerId();
0437       int channelID = elecid_crystal.channelId() - 1;
0438
0439       int hashedIndex = 100000 * eta + phi;
0440
0441       if (channelMapEE.count(hashedIndex) == 0) {
0442         channelMapEE[hashedIndex] = channelIteratorEE;
0443         channelIteratorEE++;
0444       }
0445
0446       unsigned int channel = channelMapEE[hashedIndex];
0447
0448       assert(channel < nCrys);
0449
0450       iEta[channel] = eta;
0451       iPhi[channel] = phi;
0452       iTowerID[channel] = towerID;
0453       iChannelID[channel] = channelID;
0454       idccID[channel] = dccID;
0455       iside[channel] = side;
0456
0457       // APD Pulse
0458       //===========
0459
0460       if ((*digiItr).size() > 10)
0461         edm::LogVerbatim("EcalABAnalyzer") << "SAMPLES SIZE > 10!" << (*digiItr).size();
0462
0463       // Loop on adc samples
0464
0465       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0466         EcalMGPASample samp_crystal(df.sample(i));
0467         adc[i] = samp_crystal.adc();
0468         adcG[i] = samp_crystal.gainId();
0469         adc[i] *= adcG[i];
0470
0471         if (i == 0)
0472           adcGain = adcG[i];
0473         if (i > 0)
0474           adcGain = TMath::Max(adcG[i], adcGain);
0475       }
0476
0477       APDPulse->setPulse(adc);
0478
0479       // Quality checks
0480       //================
0481
0482       if (adcGain != 1)
0483         nEvtBadGain[channel]++;
0484       if (!APDPulse->isTimingQualOK())
0485         nEvtBadTiming[channel]++;
0486       nEvtTot[channel]++;
0487
0488       // Fill if Pulse is fine
0489       //=======================
0490
0491       if (APDPulse->isPulseOK() && lightside == side) {
0492         Delta01->addEntry(APDPulse->getDelta(0, 1));
0493         Delta12->addEntry(APDPulse->getDelta(1, 2));
0494
0495         if (nevtAB[channel] < _nevtmax && _fitab) {
0496           if (doesABTreeExist)
0497             shapana->putAllVals(channel, adc, eta, phi);
0498           else
0499             shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0500           nevtAB[channel]++;
0501         }
0502       }
0503     }
0504   }
0505 }  // analyze
0506
0507 //========================================================================
0508 void EcalABAnalyzer::endJob() {
0509   //========================================================================
0510
0511   edm::LogVerbatim("EcalABAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0512   edm::LogVerbatim("EcalABAnalyzer") << "\t+=+     Analyzing data: getting (alpha, beta)     +=+";
0513
0514   // Adjust channel numbers for EE
0515   //===============================
0516
0517   if (_ecalPart == "EE") {
0518     nCrys = channelMapEE.size();
0519     shapana->set_nch(nCrys);
0520   }
0521
0522   // Set presamples number
0523   //======================
0524
0525   double delta01 = Delta01->getMean();
0526   double delta12 = Delta12->getMean();
0527   if (delta12 > _presamplecut) {
0528     _presample = 2;
0529     if (delta01 > _presamplecut)
0530       _presample = 1;
0531   }
0532
0533   APDPulse->setPresamples(_presample);
0534   shapana->set_presample(_presample);
0535
0536   //  Get alpha and beta
0537   //======================
0538
0539   if (_fitab) {
0540     edm::LogVerbatim("EcalABAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0541     edm::LogVerbatim("EcalABAnalyzer") << "\t+=+     Analyzing data: getting (alpha, beta)     +=+";
0542     TFile* fAB = nullptr;
0543     TTree* ABInit = nullptr;
0544     if (doesABTreeExist) {
0545       fAB = new TFile(alphainitfile.c_str());
0546     }
0547     if (doesABTreeExist && fAB) {
0548       ABInit = (TTree*)fAB->Get("ABCol0");
0549     }
0550     shapana->computeShape(alphafile, ABInit);
0551
0552     // Set quality flags for gains and timing
0553
0554     double BadGainEvtPercentage = 0.0;
0555     double BadTimingEvtPercentage = 0.0;
0556
0557     int nChanBadGain = 0;
0558     int nChanBadTiming = 0;
0559
0560     for (unsigned int i = 0; i < nCrys; i++) {
0561       if (nEvtTot[i] != 0) {
0562         BadGainEvtPercentage = double(nEvtBadGain[i]) / double(nEvtTot[i]);
0563         BadTimingEvtPercentage = double(nEvtBadTiming[i]) / double(nEvtTot[i]);
0564       }
0565       if (BadGainEvtPercentage > _qualpercent) {
0566         wasGainOK[i] = false;
0567         nChanBadGain++;
0568       }
0569       if (BadTimingEvtPercentage > _qualpercent) {
0570         wasTimingOK[i] = false;
0571         nChanBadTiming++;
0572       }
0573     }
0574
0575     double BadGainChanPercentage = double(nChanBadGain) / double(nCrys);
0576     double BadTimingChanPercentage = double(nChanBadTiming) / double(nCrys);
0577
0578     if (BadGainChanPercentage > _qualpercent)
0579       isGainOK = false;
0580     if (BadTimingChanPercentage > _qualpercent)
0581       isTimingOK = false;
0582
0583     if (!isGainOK)
0584       edm::LogVerbatim("EcalABAnalyzer") << "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1    +=+";
0585     if (!isTimingOK)
0586       edm::LogVerbatim("EcalABAnalyzer") << "\t+=+ ............................ WARNING! TIMING WAS BAD        +=+";
0587
0588     edm::LogVerbatim("EcalABAnalyzer") << "\t+=+    .................................... done  +=+";
0589     edm::LogVerbatim("EcalABAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0590   }
0591 }
0592
0593 DEFINE_FWK_MODULE(EcalABAnalyzer);