Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57:49 

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     for (Long64_t jentry = 0; jentry < ADCtrees[iCry]->GetEntriesFast(); jentry++) {  // Loop on events
0869       ADCtrees[iCry]->GetEntry(jentry);
0870 
0871       flagfit = 1;
0872       apdAmpl = 0.0;
0873       apdTime = 0.0;
0874       ieta = eta;
0875       iphi = phi;
0876 
0877       // Get back color
0878 
0879       unsigned int iCol = 0;
0880       for (unsigned int i = 0; i < nCol; i++) {
0881         if (color == colors[i]) {
0882           iCol = i;
0883           i = colors.size();
0884         }
0885       }
0886 
0887       // Amplitude calculation
0888 
0889       APDPulse->setPulse(adc);
0890       adcNoPed = APDPulse->getAdcWithoutPedestal();
0891 
0892       if (isSPRFine && APDPulse->isPulseOK()) {
0893         psfit->doFit(&adcNoPed[0]);
0894         apdAmpl = psfit->getAmpl();
0895         apdTime = psfit->getTime();
0896 
0897       } else {
0898         apdAmpl = 0;
0899         apdTime = 0;
0900         flagfit = 0;
0901       }
0902 
0903       if (_debug >= 1)
0904         edm::LogVerbatim("EcalLaserAnalyzer2")
0905             << "-- debug test -- endJob -- apdAmpl:" << apdAmpl << " apdTime:" << apdTime;
0906       double pnmean;
0907       if (pn0 < 10 && pn1 > 10) {
0908         pnmean = pn1;
0909       } else if (pn1 < 10 && pn0 > 10) {
0910         pnmean = pn0;
0911       } else
0912         pnmean = 0.5 * (pn0 + pn1);
0913 
0914       if (_debug >= 1)
0915         edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- endJob -- pnMean:" << pnmean;
0916 
0917       // Fill PN stuff
0918       //===============
0919 
0920       if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
0921         for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
0922           PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
0923         }
0924       }
0925 
0926       // Fill APD stuff
0927       //================
0928 
0929       if (apdAmpl != 0.0)
0930         APDFirstAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime);
0931       if (_debug >= 1)
0932         edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- endJob -- filling APDTree";
0933 
0934       APDtrees[iCry]->Fill();
0935 
0936       // Fill reference trees
0937       //=====================
0938 
0939       if (apdRefMap[0][iMod + 1] == iCry || apdRefMap[1][iMod + 1] == iCry) {
0940         apdAmplA = 0.0;
0941         apdAmplB = 0.0;
0942         eventref = event;
0943         colorref = color;
0944 
0945         for (unsigned int ir = 0; ir < nRefChan; ir++) {
0946           if (_debug >= 1)
0947             edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- ir:" << ir << " tt:" << towerID
0948                                                    << " refmap:" << apdRefMap[ir][iMod + 1] << " iCry:" << iCry;
0949 
0950           if (apdRefMap[ir][iMod + 1] == iCry) {
0951             if (_debug >= 1)
0952               edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- cut passed ";
0953             if (ir == 0)
0954               apdAmplA = apdAmpl;
0955             else if (ir == 1)
0956               apdAmplB = apdAmpl;
0957             if (_debug >= 1)
0958               edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- apdAmplA=" << apdAmplA;
0959             if (_debug >= 1)
0960               edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- apdAmplB=" << apdAmplB;
0961             if (_debug >= 1)
0962               edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- color=" << color << ", event:" << event
0963                                                      << ", ir:" << ir << " tt-1:" << towerID - 1;
0964 
0965             RefAPDtrees[ir][iMod + 1]->Fill();
0966 
0967             if (_debug >= 1)
0968               edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- tree filled" << event;
0969           }
0970         }
0971       }
0972     }
0973   }
0974 
0975   delete psfit;
0976 
0977   ADCFile->Close();
0978 
0979   if (_debug == 1)
0980     edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- endJob -- after apdAmpl Loop";
0981 
0982   // Remove temporary file
0983   //=======================
0984   stringstream del;
0985   del << "rm " << ADCfile;
0986   system(del.str().c_str());
0987 
0988   // Create output trees
0989   //=====================
0990 
0991   resFile = new TFile(resfile.c_str(), "RECREATE");
0992 
0993   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
0994     stringstream nametree;
0995     nametree << "APDCol" << colors.at(iColor);
0996     stringstream nametree2;
0997     nametree2 << "PNCol" << colors.at(iColor);
0998 
0999     restrees[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
1000     respntrees[iColor] = new TTree(nametree2.str().c_str(), nametree2.str().c_str());
1001 
1002     restrees[iColor]->Branch("iphi", &iphi, "iphi/I");
1003     restrees[iColor]->Branch("ieta", &ieta, "ieta/I");
1004     restrees[iColor]->Branch("side", &side, "side/I");
1005     restrees[iColor]->Branch("dccID", &dccID, "dccID/I");
1006     restrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1007     restrees[iColor]->Branch("towerID", &towerID, "towerID/I");
1008     restrees[iColor]->Branch("channelID", &channelID, "channelID/I");
1009     restrees[iColor]->Branch("APD", &APD, "APD[6]/D");
1010     restrees[iColor]->Branch("Time", &Time, "Time[6]/D");
1011     restrees[iColor]->Branch("APDoPN", &APDoPN, "APDoPN[6]/D");
1012     restrees[iColor]->Branch("APDoPNA", &APDoPNA, "APDoPNA[6]/D");
1013     restrees[iColor]->Branch("APDoPNB", &APDoPNB, "APDoPNB[6]/D");
1014     restrees[iColor]->Branch("APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D");
1015     restrees[iColor]->Branch("APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D");
1016     restrees[iColor]->Branch("ShapeCor", &ShapeCor, "ShapeCor/D");
1017     restrees[iColor]->Branch("flag", &flag, "flag/I");
1018 
1019     respntrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1020     respntrees[iColor]->Branch("pnID", &pnID, "pnID/I");
1021     respntrees[iColor]->Branch("PN", &PN, "PN[6]/D");
1022     respntrees[iColor]->Branch("PNoPN", &PNoPN, "PNoPN[6]/D");
1023     respntrees[iColor]->Branch("PNoPNA", &PNoPNA, "PNoPNA[6]/D");
1024     respntrees[iColor]->Branch("PNoPNB", &PNoPNB, "PNoPNB[6]/D");
1025 
1026     restrees[iColor]->SetBranchAddress("iphi", &iphi);
1027     restrees[iColor]->SetBranchAddress("ieta", &ieta);
1028     restrees[iColor]->SetBranchAddress("dccID", &dccID);
1029     restrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1030     restrees[iColor]->SetBranchAddress("towerID", &towerID);
1031     restrees[iColor]->SetBranchAddress("channelID", &channelID);
1032     restrees[iColor]->SetBranchAddress("APD", APD);
1033     restrees[iColor]->SetBranchAddress("Time", Time);
1034     restrees[iColor]->SetBranchAddress("APDoPN", APDoPN);
1035     restrees[iColor]->SetBranchAddress("APDoPNA", APDoPNA);
1036     restrees[iColor]->SetBranchAddress("APDoPNB", APDoPNB);
1037     restrees[iColor]->SetBranchAddress("APDoAPDA", APDoAPDA);
1038     restrees[iColor]->SetBranchAddress("APDoAPDB", APDoAPDB);
1039     restrees[iColor]->SetBranchAddress("ShapeCor", &ShapeCor);
1040     restrees[iColor]->SetBranchAddress("flag", &flag);
1041 
1042     respntrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1043     respntrees[iColor]->SetBranchAddress("pnID", &pnID);
1044     respntrees[iColor]->SetBranchAddress("PN", PN);
1045     respntrees[iColor]->SetBranchAddress("PNoPN", PNoPN);
1046     respntrees[iColor]->SetBranchAddress("PNoPNA", PNoPNA);
1047     respntrees[iColor]->SetBranchAddress("PNoPNB", PNoPNB);
1048   }
1049 
1050   // Set Cuts for PN stuff
1051   //=======================
1052 
1053   for (unsigned int iM = 0; iM < nMod; iM++) {
1054     unsigned int iMod = modules[iM] - 1;
1055 
1056     for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
1057       for (unsigned int icol = 0; icol < nCol; icol++) {
1058         PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
1059                                           PNFirstAnal[iMod][ich][icol]->getPN().at(1));
1060       }
1061     }
1062   }
1063 
1064   // Build ref trees indexes
1065   //========================
1066   for (unsigned int imod = 0; imod < nMod; imod++) {
1067     int jmod = modules[imod];
1068     if (RefAPDtrees[0][jmod]->GetEntries() != 0 && RefAPDtrees[1][jmod]->GetEntries() != 0) {
1069       RefAPDtrees[0][jmod]->BuildIndex("eventref");
1070       RefAPDtrees[1][jmod]->BuildIndex("eventref");
1071     }
1072   }
1073 
1074   // Final loop on crystals
1075   //=======================
1076 
1077   for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
1078     unsigned int iMod = iModule[iCry] - 1;
1079 
1080     // Set cuts on APD stuff
1081     //=======================
1082 
1083     for (unsigned int iCol = 0; iCol < nCol; iCol++) {
1084       std::vector<double> lowcut;
1085       std::vector<double> highcut;
1086       double cutMin;
1087       double cutMax;
1088 
1089       cutMin = APDFirstAnal[iCry][iCol]->getAPD().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1090       if (cutMin < 0)
1091         cutMin = 0;
1092       cutMax = APDFirstAnal[iCry][iCol]->getAPD().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1093 
1094       lowcut.push_back(cutMin);
1095       highcut.push_back(cutMax);
1096 
1097       cutMin = APDFirstAnal[iCry][iCol]->getTime().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1098       cutMax = APDFirstAnal[iCry][iCol]->getTime().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1099       lowcut.push_back(cutMin);
1100       highcut.push_back(cutMax);
1101 
1102       APDAnal[iCry][iCol] = new TAPD();
1103       APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),
1104                                      APDFirstAnal[iCry][iCol]->getAPD().at(1));
1105       APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),
1106                                         APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
1107       APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),
1108                                          APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
1109       APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),
1110                                          APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
1111       APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),
1112                                       APDFirstAnal[iCry][iCol]->getTime().at(1));
1113       APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut, highcut);
1114       APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut, highcut);
1115     }
1116 
1117     // Final loop on events
1118     //=======================
1119 
1120     for (Long64_t jentry = 0; jentry < APDtrees[iCry]->GetEntriesFast(); jentry++) {
1121       APDtrees[iCry]->GetEntry(jentry);
1122 
1123       double pnmean;
1124       if (pn0 < 10 && pn1 > 10) {
1125         pnmean = pn1;
1126       } else if (pn1 < 10 && pn0 > 10) {
1127         pnmean = pn0;
1128       } else
1129         pnmean = 0.5 * (pn0 + pn1);
1130 
1131       // Get back color
1132       //================
1133 
1134       unsigned int iCol = 0;
1135       for (unsigned int i = 0; i < nCol; i++) {
1136         if (color == colors[i]) {
1137           iCol = i;
1138           i = colors.size();
1139         }
1140       }
1141 
1142       // Fill PN stuff
1143       //===============
1144 
1145       if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
1146         for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
1147           PNAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
1148         }
1149       }
1150 
1151       // Get ref amplitudes
1152       //===================
1153 
1154       if (_debug >= 1)
1155         edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- LastLoop event:" << event << " apdAmpl:" << apdAmpl;
1156       apdAmplA = 0.0;
1157       apdAmplB = 0.0;
1158 
1159       for (unsigned int iRef = 0; iRef < nRefChan; iRef++) {
1160         RefAPDtrees[iRef][iMod + 1]->GetEntryWithIndex(event);
1161       }
1162 
1163       if (_debug == 1)
1164         edm::LogVerbatim("EcalLaserAnalyzer2")
1165             << "-- debug test -- LastLoop apdAmplA:" << apdAmplA << " apdAmplB:" << apdAmplB << ", event:" << event
1166             << ", eventref:" << eventref;
1167 
1168       // Fill APD stuff
1169       //===============
1170 
1171       APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
1172     }
1173 
1174     moduleID = iMod + 1;
1175 
1176     if (moduleID >= 20)
1177       moduleID -= 2;  // Trick to fix endcap specificity
1178 
1179     // Get final results for APD
1180     //===========================
1181 
1182     for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1183       std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
1184       std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
1185       std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
1186       std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
1187       std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
1188       std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
1189       std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
1190 
1191       for (unsigned int i = 0; i < apdvec.size(); i++) {
1192         APD[i] = apdvec.at(i);
1193         APDoPN[i] = apdpnvec.at(i);
1194         APDoPNA[i] = apdpn0vec.at(i);
1195         APDoPNB[i] = apdpn1vec.at(i);
1196         APDoAPDA[i] = apdapd0vec.at(i);
1197         APDoAPDB[i] = apdapd1vec.at(i);
1198         Time[i] = timevec.at(i);
1199         ShapeCor = shapeCorrection;
1200       }
1201 
1202       // Fill APD results trees
1203       //========================
1204 
1205       iphi = iPhi[iCry];
1206       ieta = iEta[iCry];
1207       dccID = idccID[iCry];
1208       towerID = iTowerID[iCry];
1209       side = iside[iCry];
1210       channelID = iChannelID[iCry];
1211 
1212       if (!wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor] == 0)
1213         flag = 1;
1214       else
1215         flag = 0;
1216 
1217       if (_debug >= 1)
1218         edm::LogVerbatim("EcalLaserAnalyzer2") << "-- debug test -- endJob -- APD[0]" << APD[0] << " APDoPN[0] "
1219                                                << APDoPN[0] << " APDoAPDA[0] " << APDoAPDA[0] << " flag " << flag;
1220       restrees[iColor]->Fill();
1221     }
1222   }
1223 
1224   // Get final results for PN
1225   //==========================
1226 
1227   for (unsigned int iM = 0; iM < nMod; iM++) {
1228     unsigned int iMod = modules[iM] - 1;
1229 
1230     side = iside[firstChanMod[iMod]];
1231 
1232     for (unsigned int ch = 0; ch < nPNPerMod; ch++) {
1233       pnID = ch;
1234       moduleID = iMod + 1;
1235 
1236       if (moduleID >= 20)
1237         moduleID -= 2;  // Trick to fix endcap specificity
1238 
1239       for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1240         std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
1241         std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
1242         std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
1243         std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
1244 
1245         for (unsigned int i = 0; i < pnvec.size(); i++) {
1246           PN[i] = pnvec.at(i);
1247           PNoPN[i] = pnopnvec.at(i);
1248           PNoPNA[i] = pnopn0vec.at(i);
1249           PNoPNB[i] = pnopn1vec.at(i);
1250         }
1251 
1252         if (_debug >= 1)
1253           edm::LogVerbatim("EcalLaserAnalyzer2")
1254               << "-- debug test -- endJob -- filling pn results'tree: PN[0]:" << PN[0] << " iModule:" << iMod
1255               << " iColor:" << iColor << " ch:" << ch;
1256 
1257         // Fill PN results trees
1258         //========================
1259 
1260         respntrees[iColor]->Fill();
1261       }
1262     }
1263   }
1264 
1265   // Remove temporary files
1266   //========================
1267   if (!_saveallevents) {
1268     APDFile->Close();
1269     stringstream del2;
1270     del2 << "rm " << APDfile;
1271     system(del2.str().c_str());
1272 
1273   } else {
1274     APDFile->cd();
1275     APDtrees[0]->Write();
1276 
1277     APDFile->Close();
1278     resFile->cd();
1279   }
1280 
1281   // Save results
1282   //===============
1283 
1284   for (unsigned int i = 0; i < nCol; i++) {
1285     restrees[i]->Write();
1286     respntrees[i]->Write();
1287   }
1288 
1289   resFile->Close();
1290 
1291   edm::LogVerbatim("EcalLaserAnalyzer2") << "\t+=+    .................................................. done  +=+";
1292   edm::LogVerbatim("EcalLaserAnalyzer2") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
1293 }
1294 
1295 //========================================================================
1296 bool EcalLaserAnalyzer2::getShapes() {
1297   //========================================================================
1298 
1299   // Get Pulse From Matacq Analysis:
1300   //================================
1301 
1302   bool IsMatacqOK = false;
1303 
1304   int doesMatFileExist = 0;
1305   int doesMatShapeExist = 0;
1306   FILE* test2;
1307   TProfile* laserShape = nullptr;
1308   test2 = fopen(matfile.c_str(), "r");
1309   if (test2)
1310     doesMatFileExist = 1;
1311 
1312   TFile* MatShapeFile;
1313   if (doesMatFileExist == 1) {
1314     MatShapeFile = new TFile(matfile.c_str());
1315     laserShape = (TProfile*)MatShapeFile->Get("shapeLaser");
1316     if (laserShape) {
1317       doesMatShapeExist = 1;
1318       double y = laserShape->Integral("w");
1319       if (y != 0)
1320         laserShape->Scale(1.0 / y);
1321     }
1322   } else {
1323     edm::LogError("file_not_found") << " ERROR! Matacq shape file not found !";
1324   }
1325   if (doesMatShapeExist)
1326     IsMatacqOK = true;
1327 
1328   // Get SPR from the average elec shape in SM6:
1329   //============================================
1330 
1331   int doesElecFileExist = 0;
1332   FILE* test;
1333   test = fopen(elecfile_.c_str(), "r");
1334   if (test)
1335     doesElecFileExist = 1;
1336 
1337   TFile* ElecShapesFile;
1338   TH1D* elecShape = nullptr;
1339 
1340   if (doesElecFileExist == 1) {
1341     ElecShapesFile = new TFile(elecfile_.c_str());
1342     stringstream name;
1343     name << "MeanElecShape";
1344     elecShape = (TH1D*)ElecShapesFile->Get(name.str().c_str());
1345     if (elecShape && doesMatShapeExist == 1) {
1346       double x = elecShape->GetMaximum();
1347       if (x != 0)
1348         elecShape->Scale(1.0 / x);
1349       isSPRFine = true;
1350     } else {
1351       isSPRFine = false;
1352     }
1353 
1354   } else {
1355     throw cms::Exception("file_not_found") << " ERROR! Elec shape file not found !";
1356   }
1357 
1358   if (IsMatacqOK) {
1359     ShapeFile = new TFile(shapefile.c_str(), "RECREATE");
1360 
1361     unsigned int nBins = int(laserShape->GetEntries());
1362     assert(nSamplesShapes == nBins);
1363     double elec_jj, laser_iiMinusjj;
1364     double sum_jj;
1365 
1366     if (isSPRFine == true) {
1367       unsigned int nBins2 = int(elecShape->GetNbinsX());
1368 
1369       if (nBins2 < nBins) {
1370         edm::LogError("cfg_error")
1371             << "EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins";
1372         isSPRFine = false;
1373       }
1374       assert(nSamplesShapes == nBins2);
1375 
1376       stringstream name;
1377       name << "PulseShape";
1378 
1379       PulseShape = new TProfile(name.str().c_str(), name.str().c_str(), nBins, -0.5, double(nBins) - 0.5);
1380 
1381       // shift shapes to have max close to the real APD max
1382 
1383       for (int ii = 0; ii < 50; ii++) {
1384         shapes[ii] = 0.0;
1385         PulseShape->Fill(ii, 0.0);
1386       }
1387 
1388       for (unsigned int ii = 0; ii < nBins - 50; ii++) {
1389         sum_jj = 0.0;
1390         for (unsigned int jj = 0; jj < ii; jj++) {
1391           elec_jj = elecShape->GetBinContent(jj + 1);
1392           laser_iiMinusjj = laserShape->GetBinContent(ii - jj + 1);
1393           sum_jj += elec_jj * laser_iiMinusjj;
1394         }
1395         PulseShape->Fill(ii + 50, sum_jj);
1396         shapes[ii + 50] = sum_jj;
1397       }
1398 
1399       double scale = PulseShape->GetMaximum();
1400       shapeCorrection = scale;
1401 
1402       if (scale != 0) {
1403         PulseShape->Scale(1.0 / scale);
1404         for (unsigned int ii = 0; ii < nBins; ii++) {
1405           shapesVec.push_back(shapes[ii] / scale);
1406         }
1407       }
1408 
1409       if (_saveshapes)
1410         PulseShape->Write();
1411     }
1412   }
1413   ShapeFile->Close();
1414 
1415   if (!_saveshapes) {
1416     stringstream del;
1417     del << "rm " << shapefile;
1418     system(del.str().c_str());
1419   }
1420 
1421   return IsMatacqOK;
1422 }
1423 
1424 void EcalLaserAnalyzer2::setGeomEB(
1425     int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr) {
1426   side = MEEBGeom::side(etaG, phiG);
1427 
1428   assert(module >= *min_element(modules.begin(), modules.end()) &&
1429          module <= *max_element(modules.begin(), modules.end()));
1430 
1431   eta = etaG;
1432   phi = phiG;
1433   channelID = 5 * (strip - 1) + xtal - 1;
1434   towerID = tower;
1435 
1436   vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1437   for (unsigned int iref = 0; iref < nRefChan; iref++) {
1438     if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1439       apdRefMap[iref][module] = channel;
1440     }
1441   }
1442 
1443   if (isFirstChanModFilled[module - 1] == 0) {
1444     firstChanMod[module - 1] = channel;
1445     isFirstChanModFilled[module - 1] = 1;
1446   }
1447 
1448   iEta[channel] = eta;
1449   iPhi[channel] = phi;
1450   iModule[channel] = module;
1451   iTowerID[channel] = towerID;
1452   iChannelID[channel] = channelID;
1453   idccID[channel] = dccID;
1454   iside[channel] = side;
1455 }
1456 
1457 void EcalLaserAnalyzer2::setGeomEE(
1458     int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr) {
1459   side = MEEEGeom::side(iX, iY, iZ);
1460 
1461   assert(module >= *min_element(modules.begin(), modules.end()) &&
1462          module <= *max_element(modules.begin(), modules.end()));
1463 
1464   eta = etaG;
1465   phi = phiG;
1466   channelID = ch;
1467   towerID = tower;
1468 
1469   vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1470   for (unsigned int iref = 0; iref < nRefChan; iref++) {
1471     if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1472       apdRefMap[iref][module] = channel;
1473     }
1474   }
1475 
1476   if (isFirstChanModFilled[module - 1] == 0) {
1477     firstChanMod[module - 1] = channel;
1478     isFirstChanModFilled[module - 1] = 1;
1479   }
1480 
1481   iEta[channel] = eta;
1482   iPhi[channel] = phi;
1483   iModule[channel] = module;
1484   iTowerID[channel] = towerID;
1485   iChannelID[channel] = channelID;
1486   idccID[channel] = dccID;
1487   iside[channel] = side;
1488 }
1489 
1490 DEFINE_FWK_MODULE(EcalLaserAnalyzer2);
[ 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!