Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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