Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /* 
0002  *
0003  *  \class EcalTestPulseAnalyzer
0004  *
0005  *  primary author: Julie Malcles - CEA/Saclay
0006  *  author: Gautier Hamel De Monchenault - CEA/Saclay
0007  */
0008 #include "TFile.h"
0009 #include "TTree.h"
0010 
0011 #include "EcalTestPulseAnalyzer.h"
0012 
0013 #include <sstream>
0014 #include <iomanip>
0015 
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 
0023 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0024 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0025 
0026 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h"
0027 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TSFit.h"
0028 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h"
0029 
0030 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h"
0031 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h"
0032 
0033 using namespace std;
0034 
0035 //========================================================================
0036 EcalTestPulseAnalyzer::EcalTestPulseAnalyzer(const edm::ParameterSet& iConfig)
0037     //========================================================================
0038     : iEvent(0),
0039       eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0040       eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0041       digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0042       digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0043       digiPNCollection_(iConfig.getParameter<std::string>("digiPNCollection")),
0044       rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0045       pnDiodeDigiToken_(consumes<EcalPnDiodeDigiCollection>(edm::InputTag(digiProducer_, digiPNCollection_))),
0046       mappingToken_(esConsumes()),
0047       // framework parameters with default values
0048       _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
0049       _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 3)),
0050       _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
0051       _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
0052       _samplemin(iConfig.getUntrackedParameter<unsigned int>("sampleMin", 3)),
0053       _samplemax(iConfig.getUntrackedParameter<unsigned int>("sampleMax", 9)),
0054       _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
0055       _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
0056       _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
0057       _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
0058       _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0059       _chi2max(iConfig.getUntrackedParameter<double>("chi2Max", 10.0)),
0060       _timeofmax(iConfig.getUntrackedParameter<double>("timeOfMax", 4.5)),
0061       _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0062       _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
0063       resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0064       nCrys(NCRYSEB),
0065       nMod(NMODEB),
0066       nGainPN(NGAINPN),
0067       nGainAPD(NGAINAPD),
0068       towerID(-1),
0069       channelID(-1),
0070       runType(-1),
0071       runNum(0),
0072       fedID(-1),
0073       dccID(-1),
0074       side(-1),
0075       iZ(1),
0076       phi(-1),
0077       eta(-1),
0078       event(0),
0079       apdAmpl(0),
0080       apdTime(0),
0081       pnAmpl(0),
0082       pnID(-1),
0083       moduleID(-1),
0084       channelIteratorEE(0)
0085 
0086 //========================================================================
0087 
0088 {
0089   if (_ecalPart == "EB") {
0090     ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0091   } else if (_ecalPart == "EE") {
0092     eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0093   }
0094 
0095   // Define geometrical constants
0096 
0097   if (_ecalPart == "EB") {
0098     nCrys = NCRYSEB;
0099   } else {
0100     nCrys = NCRYSEE;
0101   }
0102 
0103   iZ = 1;
0104   if (_fedid <= 609)
0105     iZ = -1;
0106 
0107   dccMEM = ME::memFromDcc(_fedid);
0108   modules = ME::lmmodFromDcc(_fedid);
0109   nMod = modules.size();
0110 }
0111 
0112 //========================================================================
0113 EcalTestPulseAnalyzer::~EcalTestPulseAnalyzer() {
0114   //========================================================================
0115 
0116   // do anything here that needs to be done at desctruction time
0117   // (e.g. close files, deallocate resources etc.)
0118 }
0119 
0120 //========================================================================
0121 void EcalTestPulseAnalyzer::beginJob() {
0122   //========================================================================
0123 
0124   // Define temporary file
0125 
0126   rootfile = resdir_;
0127   rootfile += "/TmpTreeTestPulseAnalyzer.root";
0128 
0129   outFile = new TFile(rootfile.c_str(), "RECREATE");
0130 
0131   for (unsigned int i = 0; i < nCrys; i++) {
0132     std::stringstream name;
0133     name << "Tree" << i;
0134 
0135     trees[i] = new TTree(name.str().c_str(), name.str().c_str());
0136 
0137     //List of branches
0138 
0139     trees[i]->Branch("iphi", &phi, "phi/I");
0140     trees[i]->Branch("ieta", &eta, "eta/I");
0141     trees[i]->Branch("side", &side, "side/I");
0142     trees[i]->Branch("dccID", &dccID, "dccID/I");
0143     trees[i]->Branch("towerID", &towerID, "towerID/I");
0144     trees[i]->Branch("channelID", &channelID, "channelID/I");
0145     trees[i]->Branch("event", &event, "event/I");
0146     trees[i]->Branch("apdGain", &apdGain, "apdGain/I");
0147     trees[i]->Branch("pnGain", &pnGain, "pnGain/I");
0148     trees[i]->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
0149     trees[i]->Branch("pnAmpl0", &pnAmpl0, "pnAmpl0/D");
0150     trees[i]->Branch("pnAmpl1", &pnAmpl1, "pnAmpl1/D");
0151 
0152     trees[i]->SetBranchAddress("ieta", &eta);
0153     trees[i]->SetBranchAddress("iphi", &phi);
0154     trees[i]->SetBranchAddress("side", &side);
0155     trees[i]->SetBranchAddress("dccID", &dccID);
0156     trees[i]->SetBranchAddress("towerID", &towerID);
0157     trees[i]->SetBranchAddress("channelID", &channelID);
0158     trees[i]->SetBranchAddress("event", &event);
0159     trees[i]->SetBranchAddress("apdGain", &apdGain);
0160     trees[i]->SetBranchAddress("pnGain", &pnGain);
0161     trees[i]->SetBranchAddress("apdAmpl", &apdAmpl);
0162     trees[i]->SetBranchAddress("pnAmpl0", &pnAmpl0);
0163     trees[i]->SetBranchAddress("pnAmpl1", &pnAmpl1);
0164   }
0165 
0166   // Initializations
0167 
0168   for (unsigned int j = 0; j < nCrys; j++) {
0169     iEta[j] = -1;
0170     iPhi[j] = -1;
0171     iModule[j] = 10;
0172     iTowerID[j] = -1;
0173     iChannelID[j] = -1;
0174     idccID[j] = -1;
0175     iside[j] = -1;
0176   }
0177 
0178   for (unsigned int j = 0; j < nMod; j++) {
0179     firstChanMod[j] = 0;
0180     isFirstChanModFilled[j] = 0;
0181   }
0182 
0183   // Define output results file name
0184 
0185   std::stringstream namefile;
0186   namefile << resdir_ << "/APDPN_TESTPULSE.root";
0187   resfile = namefile.str();
0188 
0189   // TP events counter
0190   TPEvents = 0;
0191 }
0192 
0193 //========================================================================
0194 void EcalTestPulseAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0195   //========================================================================
0196 
0197   ++iEvent;
0198 
0199   // Retrieve DCC header
0200   edm::Handle<EcalRawDataCollection> pDCCHeader;
0201   const EcalRawDataCollection* DCCHeader = nullptr;
0202   e.getByToken(rawDataToken_, pDCCHeader);
0203   if (!pDCCHeader.isValid()) {
0204     edm::LogError("nodata") << "Error! can't get the product  retrieving DCC header" << eventHeaderCollection_.c_str();
0205   } else {
0206     DCCHeader = pDCCHeader.product();
0207   }
0208 
0209   // retrieving crystal data from Event
0210   edm::Handle<EBDigiCollection> pEBDigi;
0211   const EBDigiCollection* EBDigi = nullptr;
0212   edm::Handle<EEDigiCollection> pEEDigi;
0213   const EEDigiCollection* EEDigi = nullptr;
0214   if (_ecalPart == "EB") {
0215     e.getByToken(ebDigiToken_, pEBDigi);
0216     if (!pEBDigi.isValid()) {
0217       edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0218     } else {
0219       EBDigi = pEBDigi.product();
0220     }
0221   } else {
0222     e.getByToken(eeDigiToken_, pEEDigi);
0223     if (!pEEDigi.isValid()) {
0224       edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0225     } else {
0226       EEDigi = pEEDigi.product();
0227     }
0228   }
0229 
0230   // Retrieve crystal PN diodes from Event
0231   edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
0232   const EcalPnDiodeDigiCollection* PNDigi = nullptr;
0233   e.getByToken(pnDiodeDigiToken_, pPNDigi);
0234   if (!pPNDigi.isValid()) {
0235     edm::LogError("nodata") << "Error! can't get the product " << digiCollection_.c_str();
0236   } else {
0237     PNDigi = pPNDigi.product();
0238   }
0239 
0240   // retrieving electronics mapping
0241   const auto& TheMapping = c.getData(mappingToken_);
0242   ;
0243 
0244   // ====================================
0245   // Decode Basic DCCHeader Information
0246   // ====================================
0247 
0248   for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0249        ++headerItr) {
0250     int fed = headerItr->fedId();
0251 
0252     if (fed != _fedid && _fedid != -999)
0253       continue;
0254 
0255     runType = headerItr->getRunType();
0256     runNum = headerItr->getRunNumber();
0257     event = headerItr->getLV1();
0258 
0259     dccID = headerItr->getDccInTCCCommand();
0260     fedID = headerItr->fedId();
0261 
0262     if (600 + dccID != fedID)
0263       continue;
0264 
0265     // Cut on runType
0266 
0267     if (runType != EcalDCCHeaderBlock::TESTPULSE_MGPA && runType != EcalDCCHeaderBlock::TESTPULSE_GAP &&
0268         runType != EcalDCCHeaderBlock::TESTPULSE_SCAN_MEM)
0269       return;
0270   }
0271 
0272   // Cut on fedID
0273 
0274   if (fedID != _fedid && _fedid != -999)
0275     return;
0276 
0277   // Count TP events
0278   TPEvents++;
0279 
0280   // ======================
0281   // Decode PN Information
0282   // ======================
0283 
0284   TPNFit* pnfit = new TPNFit();
0285   pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
0286 
0287   double chi2pn = 0;
0288   double ypnrange[50];
0289   double dsum = 0.;
0290   double dsum1 = 0.;
0291   double bl = 0.;
0292   double val_max = 0.;
0293   int samplemax = 0;
0294   unsigned int k;
0295   int pnG[50];
0296   int pngain = 0;
0297 
0298   std::map<int, std::vector<double> > allPNAmpl;
0299   std::map<int, std::vector<int> > allPNGain;
0300 
0301   for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end();
0302        ++pnItr) {  // Loop on PNs
0303 
0304     EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
0305 
0306     bool isMemRelevant = false;
0307     for (unsigned int imem = 0; imem < dccMEM.size(); imem++) {
0308       if (pnDetId.iDCCId() == dccMEM[imem]) {
0309         isMemRelevant = true;
0310       }
0311     }
0312 
0313     // skip mem dcc without relevant data
0314     if (!isMemRelevant)
0315       continue;
0316 
0317     for (int samId = 0; samId < (*pnItr).size(); samId++) {  // Loop on PN samples
0318       pn[samId] = (*pnItr).sample(samId).adc();
0319       pnG[samId] = (*pnItr).sample(samId).gainId();
0320 
0321       if (pnG[samId] != 1)
0322         edm::LogVerbatim("EcalTestPulseAnalyzer") << "PN gain different from 1 for sample " << samId;
0323       if (samId == 0)
0324         pngain = pnG[samId];
0325       if (samId > 0)
0326         pngain = TMath::Max(pnG[samId], pngain);
0327     }
0328 
0329     for (dsum = 0., k = 0; k < _presamplePN; k++) {
0330       dsum += pn[k];
0331     }
0332     bl = dsum / ((double)_presamplePN);
0333 
0334     for (val_max = 0., k = 0; k < _nsamplesPN; k++) {
0335       ypnrange[k] = pn[k] - bl;
0336 
0337       if (ypnrange[k] > val_max) {
0338         val_max = ypnrange[k];
0339         samplemax = k;
0340       }
0341     }
0342 
0343     chi2pn = pnfit->doFit(samplemax, &ypnrange[0]);
0344 
0345     if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
0346       pnAmpl = 0.;
0347     else
0348       pnAmpl = pnfit->getAmpl();
0349 
0350     allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
0351     allPNGain[pnDetId.iDCCId()].push_back(pngain);
0352   }
0353 
0354   // ===========================
0355   // Decode EBDigis Information
0356   // ===========================
0357 
0358   TSFit* pstpfit = new TSFit(_nsamples, 650);
0359   pstpfit->set_params(
0360       _nsamples, _niter, _presample, _samplemin, _samplemax, _timeofmax, _chi2max, _firstsample, _lastsample);
0361   pstpfit->init_errmat(10.);
0362 
0363   double chi2 = 0;
0364   double yrange[10];
0365   int adcgain = 0;
0366   int adcG[10];
0367 
0368   if (EBDigi) {
0369     for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
0370          ++digiItr) {  // Loop on EB crystals
0371 
0372       EBDetId id_crystal(digiItr->id());
0373       EBDataFrame df(*digiItr);
0374 
0375       int etaG = id_crystal.ieta();  // global
0376       int phiG = id_crystal.iphi();  // global
0377 
0378       int etaL;  // local
0379       int phiL;  // local
0380       std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0381 
0382       etaL = LocalCoord.first;
0383       phiL = LocalCoord.second;
0384 
0385       eta = etaG;
0386       phi = phiG;
0387 
0388       side = MEEBGeom::side(etaG, phiG);
0389 
0390       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0391 
0392       towerID = elecid_crystal.towerId();
0393       int strip = elecid_crystal.stripId();
0394       int xtal = elecid_crystal.xtalId();
0395       channelID = 5 * (strip - 1) + xtal - 1;  // FIXME
0396 
0397       int module = MEEBGeom::lmmod(etaG, phiG);
0398       int iMod = module - 1;
0399 
0400       assert(module >= *min_element(modules.begin(), modules.end()) &&
0401              module <= *max_element(modules.begin(), modules.end()));
0402 
0403       std::pair<int, int> pnpair = MEEBGeom::pn(module);
0404       unsigned int MyPn0 = pnpair.first;
0405       unsigned int MyPn1 = pnpair.second;
0406 
0407       unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0408 
0409       if (isFirstChanModFilled[iMod] == 0) {
0410         firstChanMod[iMod] = channel;
0411         isFirstChanModFilled[iMod] = 1;
0412       }
0413 
0414       iEta[channel] = eta;
0415       iPhi[channel] = phi;
0416       iModule[channel] = module;
0417       iTowerID[channel] = towerID;
0418       iChannelID[channel] = channelID;
0419       idccID[channel] = dccID;
0420       iside[channel] = side;
0421 
0422       // get adc samples
0423       //====================
0424 
0425       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0426         EcalMGPASample samp_crystal(df.sample(i));
0427         adc[i] = samp_crystal.adc();
0428         adcG[i] = samp_crystal.gainId();
0429 
0430         if (i == 0)
0431           adcgain = adcG[i];
0432         if (i > 0)
0433           adcgain = TMath::Max(adcG[i], adcgain);
0434       }
0435       // Remove pedestal
0436       //====================
0437       for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0438         dsum += adc[k];
0439         if (k < _presample - 1)
0440           dsum1 += adc[k];
0441       }
0442 
0443       bl = dsum / ((double)_presample);
0444 
0445       for (val_max = 0., k = 0; k < _nsamples; k++) {
0446         yrange[k] = adc[k] - bl;
0447         if (yrange[k] > val_max) {
0448           val_max = yrange[k];
0449         }
0450       }
0451 
0452       apdGain = adcgain;
0453 
0454       if (allPNAmpl[dccMEM[0]].size() > MyPn0)
0455         pnAmpl0 = allPNAmpl[dccMEM[0]][MyPn0];
0456       else
0457         pnAmpl0 = 0;
0458       if (allPNAmpl[dccMEM[0]].size() > MyPn1)
0459         pnAmpl1 = allPNAmpl[dccMEM[0]][MyPn1];
0460       else
0461         pnAmpl1 = 0;
0462 
0463       if (allPNGain[dccMEM[0]].size() > MyPn0)
0464         pnGain = allPNGain[dccMEM[0]][MyPn0];
0465       else
0466         pnGain = 0;
0467 
0468       // Perform the fit on apd samples
0469       //================================
0470 
0471       chi2 = pstpfit->fit_third_degree_polynomial(&yrange[0], ret_data);
0472 
0473       //Retrieve APD amplitude from fit
0474       //================================
0475 
0476       if (val_max > 100000. || chi2 < 0. || chi2 == 102) {
0477         apdAmpl = 0;
0478         apdTime = 0;
0479 
0480       } else {
0481         apdAmpl = ret_data[0];
0482         apdTime = ret_data[1];
0483       }
0484 
0485       trees[channel]->Fill();
0486     }
0487 
0488   } else {
0489     for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end();
0490          ++digiItr) {  // Loop on EE crystals
0491 
0492       EEDetId id_crystal(digiItr->id());
0493       EEDataFrame df(*digiItr);
0494 
0495       phi = id_crystal.ix();
0496       eta = id_crystal.iy();
0497 
0498       int iX = (phi - 1) / 5 + 1;
0499       int iY = (eta - 1) / 5 + 1;
0500 
0501       side = MEEEGeom::side(iX, iY, iZ);
0502 
0503       // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
0504 
0505       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0506 
0507       towerID = elecid_crystal.towerId();
0508       channelID = elecid_crystal.channelId() - 1;
0509 
0510       int module = MEEEGeom::lmmod(iX, iY);
0511       if (module >= 18 && side == 1)
0512         module += 2;  // Trick to fix endcap specificity
0513       int iMod = module - 1;
0514 
0515       assert(module >= *min_element(modules.begin(), modules.end()) &&
0516              module <= *max_element(modules.begin(), modules.end()));
0517 
0518       std::pair<int, int> pnpair = MEEEGeom::pn(module, _fedid);
0519 
0520       unsigned int MyPn0 = pnpair.first;
0521       unsigned int MyPn1 = pnpair.second;
0522 
0523       int hashedIndex = 100000 * eta + phi;
0524 
0525       if (channelMapEE.count(hashedIndex) == 0) {
0526         channelMapEE[hashedIndex] = channelIteratorEE;
0527         channelIteratorEE++;
0528       }
0529 
0530       unsigned int channel = channelMapEE[hashedIndex];
0531 
0532       if (isFirstChanModFilled[iMod] == 0) {
0533         firstChanMod[iMod] = channel;
0534         isFirstChanModFilled[iMod] = 1;
0535       }
0536 
0537       iEta[channel] = eta;
0538       iPhi[channel] = phi;
0539       iModule[channel] = module;
0540       iTowerID[channel] = towerID;
0541       iChannelID[channel] = channelID;
0542       idccID[channel] = dccID;
0543       iside[channel] = side;
0544 
0545       assert(channel < nCrys);
0546 
0547       // Get adc samples
0548       //====================
0549 
0550       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0551         EcalMGPASample samp_crystal(df.sample(i));
0552         adc[i] = samp_crystal.adc();
0553         adcG[i] = samp_crystal.gainId();
0554 
0555         if (i == 0)
0556           adcgain = adcG[i];
0557         if (i > 0)
0558           adcgain = TMath::Max(adcG[i], adcgain);
0559       }
0560 
0561       // Remove pedestal
0562       //====================
0563       for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0564         dsum += adc[k];
0565         if (k < _presample - 1)
0566           dsum1 += adc[k];
0567       }
0568 
0569       bl = dsum / ((double)_presample);
0570 
0571       for (val_max = 0., k = 0; k < _nsamples; k++) {
0572         yrange[k] = adc[k] - bl;
0573         if (yrange[k] > val_max) {
0574           val_max = yrange[k];
0575         }
0576       }
0577       apdGain = adcgain;
0578 
0579       int dccMEMIndex = 0;
0580       if (side == 1)
0581         dccMEMIndex += 2;  // Trick to fix endcap specificity
0582 
0583       if (allPNAmpl[dccMEM[dccMEMIndex]].size() > MyPn0)
0584         pnAmpl0 = allPNAmpl[dccMEM[dccMEMIndex]][MyPn0];
0585       else
0586         pnAmpl0 = 0;
0587       if (allPNAmpl[dccMEM[dccMEMIndex + 1]].size() > MyPn1)
0588         pnAmpl1 = allPNAmpl[dccMEM[dccMEMIndex + 1]][MyPn1];
0589       else
0590         pnAmpl1 = 0;
0591 
0592       if (allPNGain[dccMEM[dccMEMIndex]].size() > MyPn0)
0593         pnGain = allPNGain[dccMEM[dccMEMIndex]][MyPn0];
0594       else
0595         pnGain = 0;
0596 
0597       // Perform the fit on apd samples
0598       //=================================
0599 
0600       chi2 = pstpfit->fit_third_degree_polynomial(&yrange[0], ret_data);
0601 
0602       //Retrieve APD amplitude from fit
0603       //=================================
0604 
0605       if (val_max > 100000. || chi2 < 0. || chi2 == 102) {
0606         apdAmpl = 0;
0607         apdTime = 0;
0608 
0609       } else {
0610         apdAmpl = ret_data[0];
0611         apdTime = ret_data[1];
0612       }
0613 
0614       trees[channel]->Fill();
0615     }
0616   }
0617 
0618 }  // end of analyze
0619 
0620 //========================================================================
0621 void EcalTestPulseAnalyzer::endJob() {
0622   //========================================================================
0623 
0624   // Don't do anything if there is no events
0625   if (TPEvents == 0) {
0626     outFile->Close();
0627 
0628     // Remove temporary file
0629 
0630     std::stringstream del;
0631     del << "rm " << rootfile;
0632     system(del.str().c_str());
0633 
0634     edm::LogVerbatim("EcalTestPulseAnalyzer") << " No TP Events ";
0635     return;
0636   }
0637 
0638   edm::LogVerbatim("EcalTestPulseAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0639   edm::LogVerbatim("EcalTestPulseAnalyzer") << "\t+=+     Analyzing test pulse data: getting APD, PN  +=+";
0640 
0641   // Create output ntuples:
0642 
0643   //edm::LogVerbatim("EcalTestPulseAnalyzer")<< "TP Test Name File "<< resfile.c_str();
0644 
0645   resFile = new TFile(resfile.c_str(), "RECREATE");
0646 
0647   restrees = new TTree("TPAPD", "TPAPD");
0648   respntrees = new TTree("TPPN", "TPPN");
0649 
0650   restrees->Branch("iphi", &iphi, "iphi/I");
0651   restrees->Branch("ieta", &ieta, "ieta/I");
0652   restrees->Branch("dccID", &dccID, "dccID/I");
0653   restrees->Branch("side", &side, "side/I");
0654   restrees->Branch("towerID", &towerID, "towerID/I");
0655   restrees->Branch("channelID", &channelID, "channelID/I");
0656   restrees->Branch("moduleID", &moduleID, "moduleID/I");
0657   restrees->Branch("flag", &flag, "flag/I");
0658   restrees->Branch("gain", &gain, "gain/I");
0659   restrees->Branch("APD", &APD, "APD[6]/D");
0660 
0661   respntrees->Branch("pnID", &pnID, "pnID/I");
0662   respntrees->Branch("moduleID", &moduleID, "moduleID/I");
0663   respntrees->Branch("gain", &gain, "gain/I");
0664   respntrees->Branch("PN", &PN, "PN[6]/D");
0665 
0666   restrees->SetBranchAddress("iphi", &iphi);
0667   restrees->SetBranchAddress("ieta", &ieta);
0668   restrees->SetBranchAddress("dccID", &dccID);
0669   restrees->SetBranchAddress("side", &side);
0670   restrees->SetBranchAddress("towerID", &towerID);
0671   restrees->SetBranchAddress("channelID", &channelID);
0672   restrees->SetBranchAddress("moduleID", &moduleID);
0673   restrees->SetBranchAddress("flag", &flag);
0674   restrees->SetBranchAddress("gain", &gain);
0675   restrees->SetBranchAddress("APD", APD);
0676 
0677   respntrees->SetBranchAddress("pnID", &pnID);
0678   respntrees->SetBranchAddress("moduleID", &moduleID);
0679   respntrees->SetBranchAddress("gain", &gain);
0680   respntrees->SetBranchAddress("PN", PN);
0681 
0682   TMom* APDAnal[1700][10];
0683   TMom* PNAnal[9][2][10];
0684 
0685   for (unsigned int iMod = 0; iMod < nMod; iMod++) {
0686     for (unsigned int ich = 0; ich < 2; ich++) {
0687       for (unsigned int ig = 0; ig < nGainPN; ig++) {
0688         PNAnal[iMod][ich][ig] = new TMom();
0689       }
0690     }
0691   }
0692 
0693   for (unsigned int iCry = 0; iCry < nCrys; iCry++) {  // Loop on data trees (ie on cristals)
0694 
0695     for (unsigned int iG = 0; iG < nGainAPD; iG++) {
0696       APDAnal[iCry][iG] = new TMom();
0697     }
0698 
0699     // Define submodule and channel number inside the submodule (as Patrice)
0700 
0701     unsigned int iMod = iModule[iCry] - 1;
0702 
0703     moduleID = iMod + 1;
0704     if (moduleID >= 20)
0705       moduleID -= 2;  // Trick to fix endcap specificity
0706 
0707     for (Long64_t jentry = 0; jentry < trees[iCry]->GetEntriesFast(); jentry++) {
0708       // PN Means and RMS
0709 
0710       if (firstChanMod[iMod] == iCry) {
0711         PNAnal[iMod][0][pnGain]->addEntry(pnAmpl0);
0712         PNAnal[iMod][1][pnGain]->addEntry(pnAmpl1);
0713       }
0714 
0715       // APD means and RMS
0716 
0717       APDAnal[iCry][apdGain]->addEntry(apdAmpl);
0718     }
0719 
0720     if (trees[iCry]->GetEntries() < 10) {
0721       flag = -1;
0722       for (int j = 0; j < 6; j++) {
0723         APD[j] = 0.0;
0724       }
0725     } else
0726       flag = 1;
0727 
0728     iphi = iPhi[iCry];
0729     ieta = iEta[iCry];
0730     dccID = idccID[iCry];
0731     side = iside[iCry];
0732     towerID = iTowerID[iCry];
0733     channelID = iChannelID[iCry];
0734 
0735     for (unsigned int ig = 0; ig < nGainAPD; ig++) {
0736       APD[0] = APDAnal[iCry][ig]->getMean();
0737       APD[1] = APDAnal[iCry][ig]->getRMS();
0738       APD[2] = APDAnal[iCry][ig]->getM3();
0739       APD[3] = APDAnal[iCry][ig]->getNevt();
0740       APD[4] = APDAnal[iCry][ig]->getMin();
0741       APD[5] = APDAnal[iCry][ig]->getMax();
0742       gain = ig;
0743 
0744       // Fill APD tree
0745 
0746       restrees->Fill();
0747     }
0748   }
0749 
0750   // Get final results for PN and PN/PN
0751 
0752   for (unsigned int ig = 0; ig < nGainPN; ig++) {
0753     for (unsigned int iMod = 0; iMod < nMod; iMod++) {
0754       for (int ch = 0; ch < 2; ch++) {
0755         pnID = ch;
0756         moduleID = iMod;
0757         if (moduleID >= 20)
0758           moduleID -= 2;  // Trick to fix endcap specificity
0759 
0760         PN[0] = PNAnal[iMod][ch][ig]->getMean();
0761         PN[1] = PNAnal[iMod][ch][ig]->getRMS();
0762         PN[2] = PNAnal[iMod][ch][ig]->getM3();
0763         PN[3] = PNAnal[iMod][ch][ig]->getNevt();
0764         PN[4] = PNAnal[iMod][ch][ig]->getMin();
0765         PN[5] = PNAnal[iMod][ch][ig]->getMax();
0766         gain = ig;
0767 
0768         // Fill PN tree
0769         respntrees->Fill();
0770       }
0771     }
0772   }
0773 
0774   outFile->Close();
0775 
0776   // Remove temporary file
0777 
0778   std::stringstream del;
0779   del << "rm " << rootfile;
0780   system(del.str().c_str());
0781 
0782   // Save final results
0783 
0784   restrees->Write();
0785   respntrees->Write();
0786   resFile->Close();
0787 
0788   edm::LogVerbatim("EcalTestPulseAnalyzer") << "\t+=+    ...................................... done  +=+";
0789   edm::LogVerbatim("EcalTestPulseAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0790 }
0791 
0792 DEFINE_FWK_MODULE(EcalTestPulseAnalyzer);
[ 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!