Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-08 08:15:42 

0001 /* 
0002  *  \class EcalLaserAnalyzer
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/EcalLaserAnalyzer.h"
0017 
0018 #include <sstream>
0019 #include <fstream>
0020 #include <iomanip>
0021 
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025 
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0030 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0031 
0032 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h"
0033 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h"
0034 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h"
0035 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h"
0036 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TAPDPulse.h"
0037 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPN.h"
0038 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNPulse.h"
0039 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNCor.h"
0040 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMem.h"
0041 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0042 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h"
0043 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h"
0044 
0045 //========================================================================
0046 EcalLaserAnalyzer::EcalLaserAnalyzer(const edm::ParameterSet& iConfig)
0047     //========================================================================
0048     : iEvent(0),
0049       eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0050       eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0051       digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0052       digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0053       digiPNCollection_(iConfig.getParameter<std::string>("digiPNCollection")),
0054       rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0055       pnDiodeDigiToken_(consumes<EcalPnDiodeDigiCollection>(edm::InputTag(digiProducer_))),
0056       mappingToken_(esConsumes()),
0057       // Framework parameters with default values
0058       _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
0059       _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 2)),
0060       _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
0061       _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
0062       _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
0063       _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
0064       _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
0065       _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
0066       _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 2)),
0067       _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 9)),
0068       _timingquallow(iConfig.getUntrackedParameter<unsigned int>("timingQualLow", 3)),
0069       _timingqualhigh(iConfig.getUntrackedParameter<unsigned int>("timingQualHigh", 8)),
0070       _ratiomincutlow(iConfig.getUntrackedParameter<double>("ratioMinCutLow", 0.4)),
0071       _ratiomincuthigh(iConfig.getUntrackedParameter<double>("ratioMinCutHigh", 0.95)),
0072       _ratiomaxcutlow(iConfig.getUntrackedParameter<double>("ratioMaxCutLow", 0.8)),
0073       _presamplecut(iConfig.getUntrackedParameter<double>("presampleCut", 5.0)),
0074       _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0075       _fitab(iConfig.getUntrackedParameter<bool>("fitAB", false)),
0076       _alpha(iConfig.getUntrackedParameter<double>("alpha", 1.5076494)),
0077       _beta(iConfig.getUntrackedParameter<double>("beta", 1.5136036)),
0078       _nevtmax(iConfig.getUntrackedParameter<unsigned int>("nEvtMax", 200)),
0079       _noise(iConfig.getUntrackedParameter<double>("noise", 2.0)),
0080       _chi2cut(iConfig.getUntrackedParameter<double>("chi2cut", 10.0)),
0081       _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0082       _docorpn(iConfig.getUntrackedParameter<bool>("doCorPN", false)),
0083       _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
0084       _saveallevents(iConfig.getUntrackedParameter<bool>("saveAllEvents", false)),
0085       _qualpercent(iConfig.getUntrackedParameter<double>("qualPercent", 0.2)),
0086       _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
0087       resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0088       pncorfile_(iConfig.getUntrackedParameter<std::string>("pnCorFile")),
0089       nCrys(NCRYSEB),
0090       nPNPerMod(NPNPERMOD),
0091       nMod(NMODEE),
0092       nSides(NSIDES),
0093       runType(-1),
0094       runNum(0),
0095       fedID(-1),
0096       dccID(-1),
0097       side(2),
0098       lightside(2),
0099       iZ(1),
0100       phi(-1),
0101       eta(-1),
0102       event(0),
0103       color(-1),
0104       pn0(0),
0105       pn1(0),
0106       apdAmpl(0),
0107       apdAmplA(0),
0108       apdAmplB(0),
0109       apdTime(0),
0110       pnAmpl(0),
0111       pnID(-1),
0112       moduleID(-1),
0113       channelIteratorEE(0)
0114 
0115 //========================================================================
0116 
0117 {
0118   if (_ecalPart == "EB") {
0119     ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0120   } else if (_ecalPart == "EE") {
0121     eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0122   }
0123 
0124   // Geometrical constants initialization
0125 
0126   if (_ecalPart == "EB") {
0127     nCrys = NCRYSEB;
0128   } else {
0129     nCrys = NCRYSEE;
0130   }
0131   iZ = 1;
0132   if (_fedid <= 609)
0133     iZ = -1;
0134   modules = ME::lmmodFromDcc(_fedid);
0135   nMod = modules.size();
0136   nRefChan = NREFCHAN;
0137 
0138   for (unsigned int j = 0; j < nCrys; j++) {
0139     iEta[j] = -1;
0140     iPhi[j] = -1;
0141     iModule[j] = 10;
0142     iTowerID[j] = -1;
0143     iChannelID[j] = -1;
0144     idccID[j] = -1;
0145     iside[j] = -1;
0146     wasTimingOK[j] = true;
0147     wasGainOK[j] = true;
0148   }
0149 
0150   for (unsigned int j = 0; j < nMod; j++) {
0151     int ii = modules[j];
0152     firstChanMod[ii - 1] = 0;
0153     isFirstChanModFilled[ii - 1] = 0;
0154   }
0155 
0156   // Quality check flags
0157 
0158   isGainOK = true;
0159   isTimingOK = true;
0160 
0161   // PN linearity corrector
0162 
0163   pnCorrector = new TPNCor(pncorfile_);
0164 
0165   // Objects dealing with pulses
0166 
0167   APDPulse = new TAPDPulse(_nsamples,
0168                            _presample,
0169                            _firstsample,
0170                            _lastsample,
0171                            _timingcutlow,
0172                            _timingcuthigh,
0173                            _timingquallow,
0174                            _timingqualhigh,
0175                            _ratiomincutlow,
0176                            _ratiomincuthigh,
0177                            _ratiomaxcutlow);
0178   PNPulse = new TPNPulse(_nsamplesPN, _presamplePN);
0179 
0180   // Object dealing with MEM numbering
0181 
0182   Mem = new TMem(_fedid);
0183 
0184   // Objects needed for npresample calculation
0185 
0186   Delta01 = new TMom();
0187   Delta12 = new TMom();
0188 }
0189 
0190 //========================================================================
0191 EcalLaserAnalyzer::~EcalLaserAnalyzer() {
0192   //========================================================================
0193 
0194   // do anything here that needs to be done at destruction time
0195   // (e.g. close files, deallocate resources etc.)
0196 }
0197 
0198 //========================================================================
0199 void EcalLaserAnalyzer::beginJob() {
0200   //========================================================================
0201 
0202   // Create temporary files and trees to save adc samples
0203   //======================================================
0204 
0205   ADCfile = resdir_;
0206   ADCfile += "/APDSamplesLaser.root";
0207 
0208   APDfile = resdir_;
0209   APDfile += "/APDPNLaserAllEvents.root";
0210 
0211   ADCFile = new TFile(ADCfile.c_str(), "RECREATE");
0212 
0213   for (unsigned int i = 0; i < nCrys; i++) {
0214     std::stringstream name;
0215     name << "ADCTree" << i + 1;
0216     ADCtrees[i] = new TTree(name.str().c_str(), name.str().c_str());
0217 
0218     ADCtrees[i]->Branch("ieta", &eta, "eta/I");
0219     ADCtrees[i]->Branch("iphi", &phi, "phi/I");
0220     ADCtrees[i]->Branch("side", &side, "side/I");
0221     ADCtrees[i]->Branch("dccID", &dccID, "dccID/I");
0222     ADCtrees[i]->Branch("towerID", &towerID, "towerID/I");
0223     ADCtrees[i]->Branch("channelID", &channelID, "channelID/I");
0224     ADCtrees[i]->Branch("event", &event, "event/I");
0225     ADCtrees[i]->Branch("color", &color, "color/I");
0226     ADCtrees[i]->Branch("adc", &adc, "adc[10]/D");
0227     ADCtrees[i]->Branch("pn0", &pn0, "pn0/D");
0228     ADCtrees[i]->Branch("pn1", &pn1, "pn1/D");
0229 
0230     ADCtrees[i]->SetBranchAddress("ieta", &eta);
0231     ADCtrees[i]->SetBranchAddress("iphi", &phi);
0232     ADCtrees[i]->SetBranchAddress("side", &side);
0233     ADCtrees[i]->SetBranchAddress("dccID", &dccID);
0234     ADCtrees[i]->SetBranchAddress("towerID", &towerID);
0235     ADCtrees[i]->SetBranchAddress("channelID", &channelID);
0236     ADCtrees[i]->SetBranchAddress("event", &event);
0237     ADCtrees[i]->SetBranchAddress("color", &color);
0238     ADCtrees[i]->SetBranchAddress("adc", adc);
0239     ADCtrees[i]->SetBranchAddress("pn0", &pn0);
0240     ADCtrees[i]->SetBranchAddress("pn1", &pn1);
0241 
0242     nevtAB[i] = 0;
0243   }
0244 
0245   // Define output results filenames and shape analyzer object (alpha,beta)
0246   //=====================================================================
0247 
0248   // 1) AlphaBeta files
0249 
0250   doesABTreeExist = true;
0251 
0252   std::stringstream nameabinitfile;
0253   nameabinitfile << resdir_ << "/ABInit.root";
0254   alphainitfile = nameabinitfile.str();
0255 
0256   std::stringstream nameabfile;
0257   nameabfile << resdir_ << "/AB.root";
0258   alphafile = nameabfile.str();
0259 
0260   FILE* test;
0261   if (_fitab)
0262     test = fopen(alphainitfile.c_str(), "r");
0263   else
0264     test = fopen(alphafile.c_str(), "r");
0265   if (test == nullptr) {
0266     doesABTreeExist = false;
0267     _fitab = true;
0268   };
0269   delete test;
0270 
0271   TFile* fAB = nullptr;
0272   TTree* ABInit = nullptr;
0273   if (doesABTreeExist) {
0274     fAB = new TFile(nameabinitfile.str().c_str());
0275   }
0276   if (doesABTreeExist && fAB) {
0277     ABInit = (TTree*)fAB->Get("ABCol0");
0278   }
0279 
0280   // 2) Shape analyzer
0281 
0282   if (doesABTreeExist && fAB && ABInit && ABInit->GetEntries() != 0) {
0283     shapana = new TShapeAnalysis(ABInit, _alpha, _beta, 5.5, 1.0);
0284     doesABTreeExist = true;
0285   } else {
0286     shapana = new TShapeAnalysis(_alpha, _beta, 5.5, 1.0);
0287     doesABTreeExist = false;
0288     _fitab = true;
0289   }
0290   shapana->set_const(_nsamples, _firstsample, _lastsample, _presample, _nevtmax, _noise, _chi2cut);
0291 
0292   if (doesABTreeExist && fAB)
0293     fAB->Close();
0294 
0295   //  2) APD file
0296 
0297   std::stringstream nameapdfile;
0298   nameapdfile << resdir_ << "/APDPN_LASER.root";
0299   resfile = nameapdfile.str();
0300 
0301   // Laser events counter
0302 
0303   laserEvents = 0;
0304 }
0305 
0306 //========================================================================
0307 void EcalLaserAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0308   //========================================================================
0309 
0310   ++iEvent;
0311 
0312   // retrieving DCC header
0313   edm::Handle<EcalRawDataCollection> pDCCHeader;
0314   const EcalRawDataCollection* DCCHeader = nullptr;
0315   e.getByToken(rawDataToken_, pDCCHeader);
0316   if (!pDCCHeader.isValid()) {
0317     edm::LogError("nodata") << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str()
0318                             << " " << eventHeaderProducer_.c_str();
0319   } else {
0320     DCCHeader = pDCCHeader.product();
0321   }
0322 
0323   //retrieving crystal data from Event
0324   edm::Handle<EBDigiCollection> pEBDigi;
0325   const EBDigiCollection* EBDigi = nullptr;
0326   edm::Handle<EEDigiCollection> pEEDigi;
0327   const EEDigiCollection* EEDigi = nullptr;
0328   if (_ecalPart == "EB") {
0329     e.getByToken(ebDigiToken_, pEBDigi);
0330     if (!pEBDigi.isValid()) {
0331       edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0332     } else {
0333       EBDigi = pEBDigi.product();
0334     }
0335   } else if (_ecalPart == "EE") {
0336     e.getByToken(eeDigiToken_, pEEDigi);
0337     if (!pEEDigi.isValid()) {
0338       edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0339     } else {
0340       EEDigi = pEEDigi.product();
0341     }
0342   } else {
0343     edm::LogError("cfg_error") << " Wrong ecalPart in cfg file ";
0344     return;
0345   }
0346 
0347   // retrieving crystal PN diodes from Event
0348   edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
0349   const EcalPnDiodeDigiCollection* PNDigi = nullptr;
0350   e.getByToken(pnDiodeDigiToken_, pPNDigi);
0351   if (!pPNDigi.isValid()) {
0352     edm::LogError("nodata") << "Error! can't get the product " << digiPNCollection_.c_str();
0353   } else {
0354     PNDigi = pPNDigi.product();
0355   }
0356 
0357   // retrieving electronics mapping
0358   const auto& TheMapping = c.getData(mappingToken_);
0359 
0360   // ============================
0361   // Decode DCCHeader Information
0362   // ============================
0363 
0364   for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0365        ++headerItr) {
0366     // Get run type and run number
0367 
0368     int fed = headerItr->fedId();
0369     if (fed != _fedid && _fedid != -999)
0370       continue;
0371 
0372     runType = headerItr->getRunType();
0373     runNum = headerItr->getRunNumber();
0374     event = headerItr->getLV1();
0375 
0376     dccID = headerItr->getDccInTCCCommand();
0377     fedID = headerItr->fedId();
0378     lightside = headerItr->getRtHalf();
0379 
0380     // Check fed corresponds to the DCC in TCC
0381 
0382     if (600 + dccID != fedID)
0383       continue;
0384 
0385     // Cut on runType
0386 
0387     if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
0388         runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
0389       return;
0390 
0391     // Retrieve laser color and event number
0392 
0393     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0394     color = settings.wavelength;
0395     if (color < 0)
0396       return;
0397 
0398     std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0399     if (iter == colors.end()) {
0400       colors.push_back(color);
0401       edm::LogVerbatim("EcalLaserAnalyzer") << " new color found " << color << " " << colors.size();
0402     }
0403   }
0404 
0405   // Cut on fedID
0406 
0407   if (fedID != _fedid && _fedid != -999)
0408     return;
0409 
0410   // Count laser events
0411 
0412   laserEvents++;
0413 
0414   // ======================
0415   // Decode PN Information
0416   // ======================
0417 
0418   TPNFit* pnfit = new TPNFit();
0419   pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
0420 
0421   double chi2pn = 0;
0422   unsigned int samplemax = 0;
0423   int pnGain = 0;
0424 
0425   std::map<int, std::vector<double> > allPNAmpl;
0426   std::map<int, std::vector<double> > allPNGain;
0427 
0428   // Loop on PNs digis
0429 
0430   for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr) {
0431     EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
0432 
0433     if (_debug == 1)
0434       edm::LogVerbatim("EcalLaserAnalyzer")
0435           << "-- debug -- Inside PNDigi - pnID=" << pnDetId.iPnId() << ", dccID=" << pnDetId.iDCCId();
0436 
0437     // Skip MEM DCC without relevant data
0438 
0439     bool isMemRelevant = Mem->isMemRelevant(pnDetId.iDCCId());
0440     if (!isMemRelevant)
0441       continue;
0442 
0443     // Loop on PN samples
0444 
0445     for (int samId = 0; samId < (*pnItr).size(); samId++) {
0446       pn[samId] = (*pnItr).sample(samId).adc();
0447       pnG[samId] = (*pnItr).sample(samId).gainId();
0448       if (samId == 0)
0449         pnGain = pnG[samId];
0450       if (samId > 0)
0451         pnGain = int(TMath::Max(pnG[samId], pnGain));
0452     }
0453 
0454     if (pnGain != 1)
0455       edm::LogVerbatim("EcalLaserAnalyzer") << "PN gain different from 1";
0456 
0457     // Calculate amplitude from pulse
0458 
0459     PNPulse->setPulse(pn);
0460     pnNoPed = PNPulse->getAdcWithoutPedestal();
0461     samplemax = PNPulse->getMaxSample();
0462     chi2pn = pnfit->doFit(samplemax, &pnNoPed[0]);
0463     if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
0464       pnAmpl = 0.;
0465     else
0466       pnAmpl = pnfit->getAmpl();
0467 
0468     // Apply linearity correction
0469 
0470     double corr = 1.0;
0471     if (_docorpn)
0472       corr = pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
0473     pnAmpl *= corr;
0474 
0475     // Fill PN ampl vector
0476 
0477     allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
0478 
0479     if (_debug == 1)
0480       edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug -- Inside PNDigi - PNampl=" << pnAmpl << ", PNgain=" << pnGain;
0481   }
0482 
0483   // ===========================
0484   // Decode EBDigis Information
0485   // ===========================
0486 
0487   int adcGain = 0;
0488 
0489   if (EBDigi) {
0490     // Loop on crystals
0491     //===================
0492 
0493     for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr) {
0494       // Retrieve geometry
0495       //===================
0496 
0497       EBDetId id_crystal(digiItr->id());
0498       EBDataFrame df(*digiItr);
0499       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0500 
0501       int etaG = id_crystal.ieta();  // global
0502       int phiG = id_crystal.iphi();  // global
0503 
0504       std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0505 
0506       int etaL = LocalCoord.first;   // local
0507       int phiL = LocalCoord.second;  // local
0508 
0509       int strip = elecid_crystal.stripId();
0510       int xtal = elecid_crystal.xtalId();
0511 
0512       int module = MEEBGeom::lmmod(etaG, phiG);
0513       int tower = elecid_crystal.towerId();
0514 
0515       int apdRefTT = MEEBGeom::apdRefTower(module);
0516 
0517       std::pair<int, int> pnpair = MEEBGeom::pn(module);
0518       unsigned int MyPn0 = pnpair.first;
0519       unsigned int MyPn1 = pnpair.second;
0520 
0521       int lmr = MEEBGeom::lmr(etaG, phiG);
0522       unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0523       assert(channel < nCrys);
0524 
0525       setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);
0526 
0527       if (_debug == 1)
0528         edm::LogVerbatim("EcalLaserAnalyzer")
0529             << "-- debug -- Inside EBDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
0530             << " modules:" << modules.size();
0531 
0532       // APD Pulse
0533       //===========
0534 
0535       // Loop on adc samples
0536 
0537       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0538         EcalMGPASample samp_crystal(df.sample(i));
0539         adc[i] = samp_crystal.adc();
0540         adcG[i] = samp_crystal.gainId();
0541         adc[i] *= adcG[i];
0542         if (i == 0)
0543           adcGain = adcG[i];
0544         if (i > 0)
0545           adcGain = TMath::Max(adcG[i], adcGain);
0546       }
0547 
0548       APDPulse->setPulse(adc);
0549 
0550       // Quality checks
0551       //================
0552 
0553       if (adcGain != 1)
0554         nEvtBadGain[channel]++;
0555       if (!APDPulse->isTimingQualOK())
0556         nEvtBadTiming[channel]++;
0557       nEvtTot[channel]++;
0558 
0559       // Associate PN ampl
0560       //===================
0561 
0562       int mem0 = Mem->Mem(lmr, 0);
0563       int mem1 = Mem->Mem(lmr, 1);
0564 
0565       if (allPNAmpl[mem0].size() > MyPn0)
0566         pn0 = allPNAmpl[mem0][MyPn0];
0567       else
0568         pn0 = 0;
0569       if (allPNAmpl[mem1].size() > MyPn1)
0570         pn1 = allPNAmpl[mem1][MyPn1];
0571       else
0572         pn1 = 0;
0573 
0574       // Fill if Pulse is fine
0575       //=======================
0576 
0577       if (APDPulse->isPulseOK() && lightside == side) {
0578         ADCtrees[channel]->Fill();
0579 
0580         Delta01->addEntry(APDPulse->getDelta(0, 1));
0581         Delta12->addEntry(APDPulse->getDelta(1, 2));
0582 
0583         if (nevtAB[channel] < _nevtmax && _fitab) {
0584           if (doesABTreeExist)
0585             shapana->putAllVals(channel, adc, eta, phi);
0586           else
0587             shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0588           nevtAB[channel]++;
0589         }
0590       }
0591     }
0592 
0593   } else if (EEDigi) {
0594     // Loop on crystals
0595     //===================
0596 
0597     for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr) {
0598       // Retrieve geometry
0599       //===================
0600 
0601       EEDetId id_crystal(digiItr->id());
0602       EEDataFrame df(*digiItr);
0603       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0604 
0605       int etaG = id_crystal.iy();
0606       int phiG = id_crystal.ix();
0607 
0608       int iX = (phiG - 1) / 5 + 1;
0609       int iY = (etaG - 1) / 5 + 1;
0610 
0611       int tower = elecid_crystal.towerId();
0612       int ch = elecid_crystal.channelId() - 1;
0613 
0614       int module = MEEEGeom::lmmod(iX, iY);
0615       if (module >= 18 && side == 1)
0616         module += 2;
0617       int lmr = MEEEGeom::lmr(iX, iY, iZ);
0618       int dee = MEEEGeom::dee(lmr);
0619       int apdRefTT = MEEEGeom::apdRefTower(lmr, module);
0620 
0621       std::pair<int, int> pnpair = MEEEGeom::pn(dee, module);
0622       unsigned int MyPn0 = pnpair.first;
0623       unsigned int MyPn1 = pnpair.second;
0624 
0625       int hashedIndex = 100000 * eta + phi;
0626       if (channelMapEE.count(hashedIndex) == 0) {
0627         channelMapEE[hashedIndex] = channelIteratorEE;
0628         channelIteratorEE++;
0629       }
0630       unsigned int channel = channelMapEE[hashedIndex];
0631       assert(channel < nCrys);
0632 
0633       setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);
0634 
0635       if (_debug == 1)
0636         edm::LogVerbatim("EcalLaserAnalyzer")
0637             << "-- debug -- Inside EEDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
0638             << " modules:" << modules.size();
0639 
0640       // APD Pulse
0641       //===========
0642 
0643       if ((*digiItr).size() > 10)
0644         edm::LogVerbatim("EcalLaserAnalyzer") << "SAMPLES SIZE > 10!" << (*digiItr).size();
0645 
0646       // Loop on adc samples
0647 
0648       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0649         EcalMGPASample samp_crystal(df.sample(i));
0650         adc[i] = samp_crystal.adc();
0651         adcG[i] = samp_crystal.gainId();
0652         adc[i] *= adcG[i];
0653 
0654         if (i == 0)
0655           adcGain = adcG[i];
0656         if (i > 0)
0657           adcGain = TMath::Max(adcG[i], adcGain);
0658       }
0659 
0660       APDPulse->setPulse(adc);
0661 
0662       // Quality checks
0663       //================
0664 
0665       if (adcGain != 1)
0666         nEvtBadGain[channel]++;
0667       if (!APDPulse->isTimingQualOK())
0668         nEvtBadTiming[channel]++;
0669       nEvtTot[channel]++;
0670 
0671       // Associate PN ampl
0672       //===================
0673 
0674       int mem0 = Mem->Mem(lmr, 0);
0675       int mem1 = Mem->Mem(lmr, 1);
0676 
0677       if (allPNAmpl[mem0].size() > MyPn0)
0678         pn0 = allPNAmpl[mem0][MyPn0];
0679       else
0680         pn0 = 0;
0681       if (allPNAmpl[mem1].size() > MyPn1)
0682         pn1 = allPNAmpl[mem1][MyPn1];
0683       else
0684         pn1 = 0;
0685 
0686       // Fill if Pulse is fine
0687       //=======================
0688 
0689       if (APDPulse->isPulseOK() && lightside == side) {
0690         ADCtrees[channel]->Fill();
0691 
0692         Delta01->addEntry(APDPulse->getDelta(0, 1));
0693         Delta12->addEntry(APDPulse->getDelta(1, 2));
0694 
0695         if (nevtAB[channel] < _nevtmax && _fitab) {
0696           if (doesABTreeExist)
0697             shapana->putAllVals(channel, adc, eta, phi);
0698           else
0699             shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0700           nevtAB[channel]++;
0701         }
0702       }
0703     }
0704   }
0705 }
0706 // analyze
0707 
0708 //========================================================================
0709 void EcalLaserAnalyzer::endJob() {
0710   //========================================================================
0711 
0712   // Adjust channel numbers for EE
0713   //===============================
0714 
0715   if (_ecalPart == "EE") {
0716     nCrys = channelMapEE.size();
0717     shapana->set_nch(nCrys);
0718   }
0719 
0720   // Set presamples number
0721   //======================
0722 
0723   double delta01 = Delta01->getMean();
0724   double delta12 = Delta12->getMean();
0725   if (delta12 > _presamplecut) {
0726     _presample = 2;
0727     if (delta01 > _presamplecut)
0728       _presample = 1;
0729   }
0730 
0731   APDPulse->setPresamples(_presample);
0732   shapana->set_presample(_presample);
0733 
0734   //  Get alpha and beta
0735   //======================
0736 
0737   if (_fitab) {
0738     edm::LogVerbatim("EcalLaserAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0739     edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+     Analyzing data: getting (alpha, beta)     +=+";
0740     TFile* fAB = nullptr;
0741     TTree* ABInit = nullptr;
0742     if (doesABTreeExist) {
0743       fAB = new TFile(alphainitfile.c_str());
0744     }
0745     if (doesABTreeExist && fAB) {
0746       ABInit = (TTree*)fAB->Get("ABCol0");
0747     }
0748     shapana->computeShape(alphafile, ABInit);
0749     edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+    .................................... done  +=+";
0750     edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0751   }
0752 
0753   // Don't do anything if there is no events
0754   //=========================================
0755 
0756   if (laserEvents == 0) {
0757     ADCFile->Close();
0758     std::stringstream del;
0759     del << "rm " << ADCfile;
0760     system(del.str().c_str());
0761     edm::LogVerbatim("EcalLaserAnalyzer") << " No Laser Events ";
0762     return;
0763   }
0764 
0765   // Set quality flags for gains and timing
0766   //=========================================
0767 
0768   double BadGainEvtPercentage = 0.0;
0769   double BadTimingEvtPercentage = 0.0;
0770 
0771   int nChanBadGain = 0;
0772   int nChanBadTiming = 0;
0773 
0774   for (unsigned int i = 0; i < nCrys; i++) {
0775     if (nEvtTot[i] != 0) {
0776       BadGainEvtPercentage = double(nEvtBadGain[i]) / double(nEvtTot[i]);
0777       BadTimingEvtPercentage = double(nEvtBadTiming[i]) / double(nEvtTot[i]);
0778     }
0779     if (BadGainEvtPercentage > _qualpercent) {
0780       wasGainOK[i] = false;
0781       nChanBadGain++;
0782     }
0783     if (BadTimingEvtPercentage > _qualpercent) {
0784       wasTimingOK[i] = false;
0785       nChanBadTiming++;
0786     }
0787   }
0788 
0789   double BadGainChanPercentage = double(nChanBadGain) / double(nCrys);
0790   double BadTimingChanPercentage = double(nChanBadTiming) / double(nCrys);
0791 
0792   if (BadGainChanPercentage > _qualpercent)
0793     isGainOK = false;
0794   if (BadTimingChanPercentage > _qualpercent)
0795     isTimingOK = false;
0796 
0797   // Analyze adc samples to get amplitudes
0798   //=======================================
0799 
0800   edm::LogVerbatim("EcalLaserAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0801   edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+     Analyzing laser data: getting APD, PN, APD/PN, PN/PN    +=+";
0802 
0803   if (!isGainOK)
0804     edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1    +=+";
0805   if (!isTimingOK)
0806     edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ ............................ WARNING! TIMING WAS BAD        +=+";
0807 
0808   APDFile = new TFile(APDfile.c_str(), "RECREATE");
0809 
0810   int ieta, iphi, channelID, towerID, flag;
0811 
0812   for (unsigned int i = 0; i < nCrys; i++) {
0813     std::stringstream name;
0814     name << "APDTree" << i + 1;
0815 
0816     APDtrees[i] = new TTree(name.str().c_str(), name.str().c_str());
0817 
0818     //List of branches
0819 
0820     APDtrees[i]->Branch("event", &event, "event/I");
0821     APDtrees[i]->Branch("color", &color, "color/I");
0822     APDtrees[i]->Branch("iphi", &iphi, "iphi/I");
0823     APDtrees[i]->Branch("ieta", &ieta, "ieta/I");
0824     APDtrees[i]->Branch("side", &side, "side/I");
0825     APDtrees[i]->Branch("dccID", &dccID, "dccID/I");
0826     APDtrees[i]->Branch("towerID", &towerID, "towerID/I");
0827     APDtrees[i]->Branch("channelID", &channelID, "channelID/I");
0828     APDtrees[i]->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
0829     APDtrees[i]->Branch("apdTime", &apdTime, "apdTime/D");
0830     if (_saveallevents)
0831       APDtrees[i]->Branch("adc", &adc, "adc[10]/D");
0832     APDtrees[i]->Branch("flag", &flag, "flag/I");
0833     APDtrees[i]->Branch("flagAB", &flagAB, "flagAB/I");
0834     APDtrees[i]->Branch("pn0", &pn0, "pn0/D");
0835     APDtrees[i]->Branch("pn1", &pn1, "pn1/D");
0836 
0837     APDtrees[i]->SetBranchAddress("event", &event);
0838     APDtrees[i]->SetBranchAddress("color", &color);
0839     APDtrees[i]->SetBranchAddress("iphi", &iphi);
0840     APDtrees[i]->SetBranchAddress("ieta", &ieta);
0841     APDtrees[i]->SetBranchAddress("side", &side);
0842     APDtrees[i]->SetBranchAddress("dccID", &dccID);
0843     APDtrees[i]->SetBranchAddress("towerID", &towerID);
0844     APDtrees[i]->SetBranchAddress("channelID", &channelID);
0845     APDtrees[i]->SetBranchAddress("apdAmpl", &apdAmpl);
0846     APDtrees[i]->SetBranchAddress("apdTime", &apdTime);
0847     if (_saveallevents)
0848       APDtrees[i]->SetBranchAddress("adc", adc);
0849     APDtrees[i]->SetBranchAddress("flag", &flag);
0850     APDtrees[i]->SetBranchAddress("flagAB", &flagAB);
0851     APDtrees[i]->SetBranchAddress("pn0", &pn0);
0852     APDtrees[i]->SetBranchAddress("pn1", &pn1);
0853   }
0854 
0855   for (unsigned int iref = 0; iref < nRefChan; iref++) {
0856     for (unsigned int imod = 0; imod < nMod; imod++) {
0857       int jmod = modules[imod];
0858 
0859       std::stringstream nameref;
0860       nameref << "refAPDTree" << imod << "_" << iref;
0861 
0862       RefAPDtrees[iref][jmod] = new TTree(nameref.str().c_str(), nameref.str().c_str());
0863 
0864       RefAPDtrees[iref][jmod]->Branch("eventref", &eventref, "eventref/I");
0865       RefAPDtrees[iref][jmod]->Branch("colorref", &colorref, "colorref/I");
0866       if (iref == 0)
0867         RefAPDtrees[iref][jmod]->Branch("apdAmplA", &apdAmplA, "apdAmplA/D");
0868       if (iref == 1)
0869         RefAPDtrees[iref][jmod]->Branch("apdAmplB", &apdAmplB, "apdAmplB/D");
0870 
0871       RefAPDtrees[iref][jmod]->SetBranchAddress("eventref", &eventref);
0872       RefAPDtrees[iref][jmod]->SetBranchAddress("colorref", &colorref);
0873       if (iref == 0)
0874         RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplA", &apdAmplA);
0875       if (iref == 1)
0876         RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplB", &apdAmplB);
0877     }
0878   }
0879 
0880   assert(colors.size() <= nColor);
0881   unsigned int nCol = colors.size();
0882 
0883   // Declare PN stuff
0884   //===================
0885 
0886   for (unsigned int iM = 0; iM < nMod; iM++) {
0887     unsigned int iMod = modules[iM] - 1;
0888 
0889     for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
0890       for (unsigned int icol = 0; icol < nCol; icol++) {
0891         PNFirstAnal[iMod][ich][icol] = new TPN(ich);
0892         PNAnal[iMod][ich][icol] = new TPN(ich);
0893       }
0894     }
0895   }
0896 
0897   // Declare function for APD ampl fit
0898   //===================================
0899 
0900   PulseFitWithFunction* pslsfit = new PulseFitWithFunction();
0901   double chi2;
0902 
0903   for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
0904     for (unsigned int icol = 0; icol < nCol; icol++) {
0905       // Declare APD stuff
0906       //===================
0907 
0908       APDFirstAnal[iCry][icol] = new TAPD();
0909       IsThereDataADC[iCry][icol] = 1;
0910       std::stringstream cut;
0911       cut << "color==" << colors[icol];
0912       if (ADCtrees[iCry]->GetEntries(cut.str().c_str()) < 10)
0913         IsThereDataADC[iCry][icol] = 0;
0914     }
0915 
0916     unsigned int iMod = iModule[iCry] - 1;
0917     double alpha, beta;
0918 
0919     // Loop on events
0920     //================
0921 
0922     Long64_t nbytes = 0, nb = 0;
0923     for (Long64_t jentry = 0; jentry < ADCtrees[iCry]->GetEntriesFast(); jentry++) {
0924       nb = ADCtrees[iCry]->GetEntry(jentry);
0925       nbytes += nb;
0926 
0927       // Get back color
0928 
0929       unsigned int iCol = 0;
0930       for (unsigned int i = 0; i < nCol; i++) {
0931         if (color == colors[i]) {
0932           iCol = i;
0933           i = colors.size();
0934         }
0935       }
0936 
0937       // Retreive alpha and beta
0938 
0939       std::vector<double> abvals = shapana->getVals(iCry);
0940       alpha = abvals[0];
0941       beta = abvals[1];
0942       flagAB = int(abvals[4]);
0943       iphi = iPhi[iCry];
0944       ieta = iEta[iCry];
0945       towerID = iTowerID[iCry];
0946       channelID = iChannelID[iCry];
0947 
0948       // Amplitude calculation
0949 
0950       APDPulse->setPulse(adc);
0951       adcNoPed = APDPulse->getAdcWithoutPedestal();
0952 
0953       apdAmpl = 0;
0954       apdAmplA = 0;
0955       apdAmplB = 0;
0956       apdTime = 0;
0957 
0958       if (APDPulse->isPulseOK()) {
0959         pslsfit->init(_nsamples, _firstsample, _lastsample, _niter, alpha, beta);
0960         chi2 = pslsfit->doFit(&adcNoPed[0]);
0961 
0962         if (chi2 < 0. || chi2 == 102 || chi2 == 101) {
0963           apdAmpl = 0;
0964           apdTime = 0;
0965           flag = 0;
0966         } else {
0967           apdAmpl = pslsfit->getAmpl();
0968           apdTime = pslsfit->getTime();
0969           flag = 1;
0970         }
0971       } else {
0972         apdAmpl = 0;
0973         apdTime = 0;
0974         flag = 0;
0975       }
0976 
0977       if (_debug == 1)
0978         edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- apdAmpl=" << apdAmpl << ", apdTime=" << apdTime;
0979       double pnmean;
0980       if (pn0 < 10 && pn1 > 10) {
0981         pnmean = pn1;
0982       } else if (pn1 < 10 && pn0 > 10) {
0983         pnmean = pn0;
0984       } else
0985         pnmean = 0.5 * (pn0 + pn1);
0986 
0987       if (_debug == 1)
0988         edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- pn0=" << pn0 << ", pn1=" << pn1;
0989 
0990       // Fill PN stuff
0991       //===============
0992 
0993       if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
0994         for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
0995           PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
0996         }
0997       }
0998 
0999       // Fill APD stuff
1000       //================
1001 
1002       if (APDPulse->isPulseOK()) {
1003         APDFirstAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime);
1004         APDtrees[iCry]->Fill();
1005 
1006         // Fill reference trees
1007         //=====================
1008 
1009         if (apdRefMap[0][iMod + 1] == iCry || apdRefMap[1][iMod + 1] == iCry) {
1010           apdAmplA = 0.0;
1011           apdAmplB = 0.0;
1012           eventref = event;
1013           colorref = color;
1014 
1015           for (unsigned int ir = 0; ir < nRefChan; ir++) {
1016             if (apdRefMap[ir][iMod + 1] == iCry) {
1017               if (ir == 0)
1018                 apdAmplA = apdAmpl;
1019               else if (ir == 1)
1020                 apdAmplB = apdAmpl;
1021               RefAPDtrees[ir][iMod + 1]->Fill();
1022             }
1023           }
1024         }
1025       }
1026     }
1027   }
1028 
1029   delete pslsfit;
1030 
1031   ADCFile->Close();
1032 
1033   // Remove temporary file
1034   //=======================
1035   std::stringstream del;
1036   del << "rm " << ADCfile;
1037   system(del.str().c_str());
1038 
1039   // Create output trees
1040   //=====================
1041 
1042   resFile = new TFile(resfile.c_str(), "RECREATE");
1043 
1044   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1045     std::stringstream nametree;
1046     nametree << "APDCol" << colors[iColor];
1047     std::stringstream nametree2;
1048     nametree2 << "PNCol" << colors[iColor];
1049 
1050     restrees[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
1051     respntrees[iColor] = new TTree(nametree2.str().c_str(), nametree2.str().c_str());
1052 
1053     restrees[iColor]->Branch("iphi", &iphi, "iphi/I");
1054     restrees[iColor]->Branch("ieta", &ieta, "ieta/I");
1055     restrees[iColor]->Branch("side", &side, "side/I");
1056     restrees[iColor]->Branch("dccID", &dccID, "dccID/I");
1057     restrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1058     restrees[iColor]->Branch("towerID", &towerID, "towerID/I");
1059     restrees[iColor]->Branch("channelID", &channelID, "channelID/I");
1060     restrees[iColor]->Branch("APD", &APD, "APD[6]/D");
1061     restrees[iColor]->Branch("Time", &Time, "Time[6]/D");
1062     restrees[iColor]->Branch("APDoPN", &APDoPN, "APDoPN[6]/D");
1063     restrees[iColor]->Branch("APDoPNA", &APDoPNA, "APDoPNA[6]/D");
1064     restrees[iColor]->Branch("APDoPNB", &APDoPNB, "APDoPNB[6]/D");
1065     restrees[iColor]->Branch("APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D");
1066     restrees[iColor]->Branch("APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D");
1067     restrees[iColor]->Branch("flag", &flag, "flag/I");
1068 
1069     respntrees[iColor]->Branch("side", &side, "side/I");
1070     respntrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1071     respntrees[iColor]->Branch("pnID", &pnID, "pnID/I");
1072     respntrees[iColor]->Branch("PN", &PN, "PN[6]/D");
1073     respntrees[iColor]->Branch("PNoPN", &PNoPN, "PNoPN[6]/D");
1074     respntrees[iColor]->Branch("PNoPNA", &PNoPNA, "PNoPNA[6]/D");
1075     respntrees[iColor]->Branch("PNoPNB", &PNoPNB, "PNoPNB[6]/D");
1076 
1077     restrees[iColor]->SetBranchAddress("iphi", &iphi);
1078     restrees[iColor]->SetBranchAddress("ieta", &ieta);
1079     restrees[iColor]->SetBranchAddress("side", &side);
1080     restrees[iColor]->SetBranchAddress("dccID", &dccID);
1081     restrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1082     restrees[iColor]->SetBranchAddress("towerID", &towerID);
1083     restrees[iColor]->SetBranchAddress("channelID", &channelID);
1084     restrees[iColor]->SetBranchAddress("APD", APD);
1085     restrees[iColor]->SetBranchAddress("Time", Time);
1086     restrees[iColor]->SetBranchAddress("APDoPN", APDoPN);
1087     restrees[iColor]->SetBranchAddress("APDoPNA", APDoPNA);
1088     restrees[iColor]->SetBranchAddress("APDoPNB", APDoPNB);
1089     restrees[iColor]->SetBranchAddress("APDoAPDA", APDoAPDA);
1090     restrees[iColor]->SetBranchAddress("APDoAPDB", APDoAPDB);
1091     restrees[iColor]->SetBranchAddress("flag", &flag);
1092 
1093     respntrees[iColor]->SetBranchAddress("side", &side);
1094     respntrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1095     respntrees[iColor]->SetBranchAddress("pnID", &pnID);
1096     respntrees[iColor]->SetBranchAddress("PN", PN);
1097     respntrees[iColor]->SetBranchAddress("PNoPN", PNoPN);
1098     respntrees[iColor]->SetBranchAddress("PNoPNA", PNoPNA);
1099     respntrees[iColor]->SetBranchAddress("PNoPNB", PNoPNB);
1100   }
1101 
1102   // Set Cuts for PN stuff
1103   //=======================
1104 
1105   for (unsigned int iM = 0; iM < nMod; iM++) {
1106     unsigned int iMod = modules[iM] - 1;
1107 
1108     for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
1109       for (unsigned int icol = 0; icol < nCol; icol++) {
1110         PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
1111                                           PNFirstAnal[iMod][ich][icol]->getPN().at(1));
1112       }
1113     }
1114   }
1115 
1116   // Build ref trees indexes
1117   //========================
1118   for (unsigned int imod = 0; imod < nMod; imod++) {
1119     int jmod = modules[imod];
1120     if (RefAPDtrees[0][jmod]->GetEntries() != 0 && RefAPDtrees[1][jmod]->GetEntries() != 0) {
1121       RefAPDtrees[0][jmod]->BuildIndex("eventref");
1122       RefAPDtrees[1][jmod]->BuildIndex("eventref");
1123     }
1124   }
1125 
1126   // Final loop on crystals
1127   //=======================
1128 
1129   for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
1130     unsigned int iMod = iModule[iCry] - 1;
1131 
1132     // Set cuts on APD stuff
1133     //=======================
1134 
1135     for (unsigned int iCol = 0; iCol < nCol; iCol++) {
1136       std::vector<double> lowcut;
1137       std::vector<double> highcut;
1138       double cutMin;
1139       double cutMax;
1140 
1141       cutMin = APDFirstAnal[iCry][iCol]->getAPD().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1142       if (cutMin < 0)
1143         cutMin = 0;
1144       cutMax = APDFirstAnal[iCry][iCol]->getAPD().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1145 
1146       lowcut.push_back(cutMin);
1147       highcut.push_back(cutMax);
1148 
1149       cutMin = APDFirstAnal[iCry][iCol]->getTime().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1150       cutMax = APDFirstAnal[iCry][iCol]->getTime().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1151       lowcut.push_back(cutMin);
1152       highcut.push_back(cutMax);
1153 
1154       APDAnal[iCry][iCol] = new TAPD();
1155       APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),
1156                                      APDFirstAnal[iCry][iCol]->getAPD().at(1));
1157       APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),
1158                                         APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
1159       APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),
1160                                          APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
1161       APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),
1162                                          APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
1163       APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),
1164                                       APDFirstAnal[iCry][iCol]->getTime().at(1));
1165       APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut, highcut);
1166       APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut, highcut);
1167     }
1168 
1169     // Final loop on events
1170     //=======================
1171 
1172     Long64_t nbytes = 0, nb = 0;
1173     for (Long64_t jentry = 0; jentry < APDtrees[iCry]->GetEntriesFast(); jentry++) {
1174       nb = APDtrees[iCry]->GetEntry(jentry);
1175       nbytes += nb;
1176 
1177       double pnmean;
1178       if (pn0 < 10 && pn1 > 10) {
1179         pnmean = pn1;
1180       } else if (pn1 < 10 && pn0 > 10) {
1181         pnmean = pn0;
1182       } else
1183         pnmean = 0.5 * (pn0 + pn1);
1184 
1185       // Get back color
1186       //================
1187 
1188       unsigned int iCol = 0;
1189       for (unsigned int i = 0; i < nCol; i++) {
1190         if (color == colors[i]) {
1191           iCol = i;
1192           i = colors.size();
1193         }
1194       }
1195 
1196       // Fill PN stuff
1197       //===============
1198 
1199       if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
1200         for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
1201           PNAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
1202         }
1203       }
1204 
1205       // Get ref amplitudes
1206       //===================
1207 
1208       if (_debug == 1)
1209         edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- Last Loop event:" << event << " apdAmpl:" << apdAmpl;
1210       apdAmplA = 0.0;
1211       apdAmplB = 0.0;
1212 
1213       for (unsigned int iRef = 0; iRef < nRefChan; iRef++) {
1214         RefAPDtrees[iRef][iMod + 1]->GetEntryWithIndex(event);
1215       }
1216 
1217       if (_debug == 1)
1218         edm::LogVerbatim("EcalLaserAnalyzer")
1219             << "-- debug test -- Last Loop apdAmplA:" << apdAmplA << " apdAmplB:" << apdAmplB << ", event:" << event
1220             << ", eventref:" << eventref;
1221 
1222       // Fill APD stuff
1223       //===============
1224 
1225       APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
1226     }
1227 
1228     moduleID = iMod + 1;
1229 
1230     if (moduleID >= 20)
1231       moduleID -= 2;  // Trick to fix endcap specificity
1232 
1233     // Get final results for APD
1234     //===========================
1235 
1236     for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1237       std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
1238       std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
1239       std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
1240       std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
1241       std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
1242       std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
1243       std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
1244 
1245       for (unsigned int i = 0; i < apdvec.size(); i++) {
1246         APD[i] = apdvec.at(i);
1247         APDoPN[i] = apdpnvec.at(i);
1248         APDoPNA[i] = apdpn0vec.at(i);
1249         APDoPNB[i] = apdpn1vec.at(i);
1250         APDoAPDA[i] = apdapd0vec.at(i);
1251         APDoAPDB[i] = apdapd1vec.at(i);
1252         Time[i] = timevec.at(i);
1253       }
1254 
1255       // Fill APD results trees
1256       //========================
1257 
1258       iphi = iPhi[iCry];
1259       ieta = iEta[iCry];
1260       dccID = idccID[iCry];
1261       side = iside[iCry];
1262       towerID = iTowerID[iCry];
1263       channelID = iChannelID[iCry];
1264 
1265       if (!wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor] == 0) {
1266         flag = 0;
1267       } else
1268         flag = 1;
1269 
1270       restrees[iColor]->Fill();
1271     }
1272   }
1273 
1274   // Get final results for PN
1275   //==========================
1276 
1277   for (unsigned int iM = 0; iM < nMod; iM++) {
1278     unsigned int iMod = modules[iM] - 1;
1279 
1280     side = iside[firstChanMod[iMod]];
1281 
1282     for (unsigned int ch = 0; ch < nPNPerMod; ch++) {
1283       pnID = ch;
1284       moduleID = iMod + 1;
1285 
1286       if (moduleID >= 20)
1287         moduleID -= 2;  // Trick to fix endcap specificity
1288 
1289       for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1290         std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
1291         std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
1292         std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
1293         std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
1294 
1295         for (unsigned int i = 0; i < pnvec.size(); i++) {
1296           PN[i] = pnvec.at(i);
1297           PNoPN[i] = pnopnvec.at(i);
1298           PNoPNA[i] = pnopn0vec.at(i);
1299           PNoPNB[i] = pnopn1vec.at(i);
1300         }
1301 
1302         // Fill PN results trees
1303         //========================
1304 
1305         respntrees[iColor]->Fill();
1306       }
1307     }
1308   }
1309 
1310   // Remove temporary files
1311   //========================
1312   if (!_saveallevents) {
1313     APDFile->Close();
1314     std::stringstream del2;
1315     del2 << "rm " << APDfile;
1316     system(del2.str().c_str());
1317 
1318   } else {
1319     APDFile->cd();
1320     APDtrees[0]->Write();
1321 
1322     APDFile->Close();
1323     resFile->cd();
1324   }
1325 
1326   // Save results
1327   //===============
1328 
1329   for (unsigned int i = 0; i < nCol; i++) {
1330     restrees[i]->Write();
1331     respntrees[i]->Write();
1332   }
1333 
1334   resFile->Close();
1335 
1336   edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+    .................................................. done  +=+";
1337   edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
1338 }
1339 
1340 void EcalLaserAnalyzer::setGeomEB(
1341     int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr) {
1342   side = MEEBGeom::side(etaG, phiG);
1343 
1344   assert(module >= *min_element(modules.begin(), modules.end()) &&
1345          module <= *max_element(modules.begin(), modules.end()));
1346 
1347   eta = etaG;
1348   phi = phiG;
1349   channelID = 5 * (strip - 1) + xtal - 1;
1350   towerID = tower;
1351 
1352   std::vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1353   for (unsigned int iref = 0; iref < nRefChan; iref++) {
1354     if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1355       apdRefMap[iref][module] = channel;
1356     }
1357   }
1358 
1359   if (isFirstChanModFilled[module - 1] == 0) {
1360     firstChanMod[module - 1] = channel;
1361     isFirstChanModFilled[module - 1] = 1;
1362   }
1363 
1364   iEta[channel] = eta;
1365   iPhi[channel] = phi;
1366   iModule[channel] = module;
1367   iTowerID[channel] = towerID;
1368   iChannelID[channel] = channelID;
1369   idccID[channel] = dccID;
1370   iside[channel] = side;
1371 }
1372 
1373 void EcalLaserAnalyzer::setGeomEE(
1374     int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr) {
1375   side = MEEEGeom::side(iX, iY, iZ);
1376 
1377   assert(module >= *min_element(modules.begin(), modules.end()) &&
1378          module <= *max_element(modules.begin(), modules.end()));
1379 
1380   eta = etaG;
1381   phi = phiG;
1382   channelID = ch;
1383   towerID = tower;
1384 
1385   std::vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1386   for (unsigned int iref = 0; iref < nRefChan; iref++) {
1387     if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1388       apdRefMap[iref][module] = channel;
1389     }
1390   }
1391 
1392   if (isFirstChanModFilled[module - 1] == 0) {
1393     firstChanMod[module - 1] = channel;
1394     isFirstChanModFilled[module - 1] = 1;
1395   }
1396 
1397   iEta[channel] = eta;
1398   iPhi[channel] = phi;
1399   iModule[channel] = module;
1400   iTowerID[channel] = towerID;
1401   iChannelID[channel] = channelID;
1402   idccID[channel] = dccID;
1403   iside[channel] = side;
1404 }
1405 
1406 DEFINE_FWK_MODULE(EcalLaserAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!