EcalLaserAnalyzer/plugins/EcalPerEvtLaserAnalyzer.cc

0001 /*
0002  *  \class EcalPerEvtLaserAnalyzer
0003  *
0004  *  primary author: Julie Malcles - CEA/Saclay
0005  *  author: Gautier Hamel De Monchenault - CEA/Saclay
0006  */
0007
0008 #include "TFile.h"
0009 #include "TTree.h"
0010
0011 #include "EcalPerEvtLaserAnalyzer.h"
0012
0013 #include <sstream>
0014 #include <iomanip>
0015 #include <ctime>
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEEGeom.h"
0020 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEBGeom.h"
0021
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026
0027 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0028 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0029
0030 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h"
0031 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0032
0033 using namespace std;
0034
0035 //========================================================================
0036 EcalPerEvtLaserAnalyzer::EcalPerEvtLaserAnalyzer(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       _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
0053       _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
0054       _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
0055       _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
0056       _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 3)),
0057       _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 7)),
0058       _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0059       _fedid(iConfig.getUntrackedParameter<unsigned int>("fedID", 0)),
0060       _tower(iConfig.getUntrackedParameter<unsigned int>("tower", 1)),
0061       _channel(iConfig.getUntrackedParameter<unsigned int>("channel", 1)),
0062       _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0063       resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0064       refalphabeta_(iConfig.getUntrackedParameter<std::string>("refAlphaBeta")),
0065       nCrys(NCRYSEB),
0066       IsFileCreated(0),
0067       runType(-1),
0068       runNum(0),
0069       dccID(-1),
0070       lightside(2),
0071       doesRefFileExist(0),
0072       ttMat(-1),
0073       peakMat(-1),
0074       peak(-1),
0075       evtMat(-1),
0076       colMat(-1)
0077 //========================================================================
0078
0079 {
0080   if (_ecalPart == "EB") {
0081     ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0082   } else if (_ecalPart == "EE") {
0083     eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0084   }
0085
0086   // Define geometrical constants
0087   //
0088   if (_ecalPart == "EB") {
0089     nCrys = NCRYSEB;
0090   } else {
0091     nCrys = NCRYSEE;
0092   }
0093 }
0094
0095 //========================================================================
0096 EcalPerEvtLaserAnalyzer::~EcalPerEvtLaserAnalyzer() {
0097   //========================================================================
0098
0099   // do anything here that needs to be done at desctruction time
0100   // (e.g. close files, deallocate resources etc.)
0101 }
0102
0103 //========================================================================
0104 void EcalPerEvtLaserAnalyzer::beginJob() {
0105   //========================================================================
0106
0107   // Define temporary files' names
0108
0109   stringstream namefile1;
0110   namefile1 << resdir_ << "/ADCSamples.root";
0111
0112   ADCfile = namefile1.str();
0113
0114   // Create temporary file and trees to save adc samples
0115
0116   ADCFile = new TFile(ADCfile.c_str(), "RECREATE");
0117
0118   stringstream name;
0119   name << "ADCTree";
0120
0121   ADCtrees = new TTree(name.str().c_str(), name.str().c_str());
0122
0123   ADCtrees->Branch("iphi", &phi, "phi/I");
0124   ADCtrees->Branch("ieta", &eta, "eta/I");
0125   ADCtrees->Branch("dccID", &dccID, "dccID/I");
0126   ADCtrees->Branch("event", &event, "event/I");
0127   ADCtrees->Branch("color", &color, "color/I");
0128   ADCtrees->Branch("adc", &adc, "adc[10]/D");
0129   ADCtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
0130   ADCtrees->Branch("peakMatacq", &peak, "peak/D");
0131   ADCtrees->Branch("pn0", &pn0, "pn0/D");
0132   ADCtrees->Branch("pn1", &pn1, "pn1/D");
0133
0134   ADCtrees->SetBranchAddress("ieta", &eta);
0135   ADCtrees->SetBranchAddress("iphi", &phi);
0136   ADCtrees->SetBranchAddress("dccID", &dccID);
0137   ADCtrees->SetBranchAddress("event", &event);
0138   ADCtrees->SetBranchAddress("color", &color);
0139   ADCtrees->SetBranchAddress("adc", adc);
0140   ADCtrees->SetBranchAddress("ttrigMatacq", &ttrig);
0141   ADCtrees->SetBranchAddress("peakMatacq", &peak);
0142   ADCtrees->SetBranchAddress("pn0", &pn0);
0143   ADCtrees->SetBranchAddress("pn1", &pn1);
0144
0145   IsFileCreated = 0;
0146 }
0147
0148 //========================================================================
0149 void EcalPerEvtLaserAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0150   //========================================================================
0151
0152   ++iEvent;
0153
0154   // retrieving DCC header
0155   edm::Handle<EcalRawDataCollection> pDCCHeader;
0156   const EcalRawDataCollection* DCCHeader = nullptr;
0157   e.getByToken(rawDataToken_, pDCCHeader);
0158   if (!pDCCHeader.isValid()) {
0159     edm::LogError("nodata") << "Error! can't get the product  retrieving DCC header" << eventHeaderCollection_.c_str();
0160   } else {
0161     DCCHeader = pDCCHeader.product();
0162   }
0163
0164   // retrieving crystal data from Event
0165   edm::Handle<EBDigiCollection> pEBDigi;
0166   const EBDigiCollection* EBDigi = nullptr;
0167   edm::Handle<EEDigiCollection> pEEDigi;
0168   const EEDigiCollection* EEDigi = nullptr;
0169   if (_ecalPart == "EB") {
0170     e.getByToken(ebDigiToken_, pEBDigi);
0171     if (!pEBDigi.isValid()) {
0172       edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0173     } else {
0174       EBDigi = pEBDigi.product();
0175     }
0176   } else {
0177     e.getByToken(eeDigiToken_, pEEDigi);
0178     if (!pEEDigi.isValid()) {
0179       edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0180     } else {
0181       EEDigi = pEEDigi.product();
0182     }
0183   }
0184
0185   // retrieving crystal PN diodes from Event
0186   edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
0187   const EcalPnDiodeDigiCollection* PNDigi = nullptr;
0188   e.getByToken(pnDiodeDigiToken_, pPNDigi);
0189   if (!pPNDigi.isValid()) {
0190     edm::LogError("nodata") << "Error! can't get the product " << digiPNCollection_.c_str();
0191   } else {
0192     PNDigi = pPNDigi.product();
0193   }
0194
0195   // retrieving electronics mapping
0196   const auto& TheMapping = c.getData(mappingToken_);
0197
0198   // ====================================
0199   // Decode Basic DCCHeader Information
0200   // ====================================
0201
0202   for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0203        ++headerItr) {
0204     // Get run type and run number
0205
0206     int fed = headerItr->fedId();
0207
0208     if (fed != _fedid && _fedid != -999)
0209       continue;
0210
0211     runType = headerItr->getRunType();
0212     runNum = headerItr->getRunNumber();
0213     event = headerItr->getLV1();
0214     dccID = headerItr->getDccInTCCCommand();
0215     fedID = headerItr->fedId();
0216
0217     // take event only if the fed corresponds to the DCC in TCC
0218     if (600 + dccID != fedID)
0219       continue;
0220
0221     // Cut on runType
0222     if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP)
0223       return;
0224
0225     // Define output results files' names
0226
0227     if (IsFileCreated == 0) {
0228       stringstream namefile2;
0229       namefile2 << resdir_ << "/APDAmpl_Run" << runNum << "_" << _fedid << "_" << _tower << "_" << _channel << ".root";
0230       resfile = namefile2.str();
0231
0232       // Get Matacq ttrig
0233
0234       stringstream namefile;
0235       namefile << resdir_ << "/Matacq-Run" << runNum << ".root";
0236
0237       doesRefFileExist = 0;
0238
0239       FILE* test;
0240       test = fopen(namefile.str().c_str(), "r");
0241       if (test)
0242         doesRefFileExist = 1;
0243
0244       if (doesRefFileExist == 1) {
0245         matacqFile = new TFile((namefile.str().c_str()));
0246         matacqTree = (TTree*)matacqFile->Get("MatacqShape");
0247
0248         matacqTree->SetBranchAddress("event", &evtMat);
0249         matacqTree->SetBranchAddress("color", &colMat);
0250         matacqTree->SetBranchAddress("peak", &peakMat);
0251         matacqTree->SetBranchAddress("ttrig", &ttMat);
0252       }
0253
0254       IsFileCreated = 1;
0255     }
0256
0257     // Retrieve laser color and event number
0258
0259     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0260     int color = settings.wavelength;
0261
0262     vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0263     if (iter == colors.end()) {
0264       colors.push_back(color);
0265       edm::LogVerbatim("EcalPerEvtLaserAnalyzer") << " new color found " << color << " " << colors.size();
0266     }
0267   }
0268
0269   // cut on fedID
0270
0271   if (fedID != _fedid && _fedid != -999)
0272     return;
0273
0274   // ======================
0275   // Decode PN Information
0276   // ======================
0277
0278   TPNFit* pnfit = new TPNFit();
0279   pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
0280
0281   double chi2pn = 0;
0282   double ypnrange[50];
0283   double dsum = 0.;
0284   double dsum1 = 0.;
0285   double bl = 0.;
0286   double bl1 = 0.;
0287   double val_max = 0.;
0288   unsigned int samplemax = 0;
0289   unsigned int k;
0290
0291   std::vector<double> allPNAmpl;
0292
0293   for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end();
0294        ++pnItr) {  // Loop on PNs
0295
0296     for (int samId = 0; samId < (*pnItr).size(); samId++) {  // Loop on PN samples
0297       pn[samId] = (*pnItr).sample(samId).adc();
0298     }
0299
0300     for (dsum = 0., k = 0; k < _presamplePN; k++) {
0301       dsum += pn[k];
0302     }
0303     bl = dsum / ((double)_presamplePN);
0304
0305     for (val_max = 0., k = 0; k < _nsamplesPN; k++) {
0306       ypnrange[k] = pn[k] - bl;
0307
0308       if (ypnrange[k] > val_max) {
0309         val_max = ypnrange[k];
0310         samplemax = k;
0311       }
0312     }
0313
0314     chi2pn = pnfit->doFit(samplemax, &ypnrange[0]);
0315
0316     if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
0317       pnAmpl = 0.;
0318     else
0319       pnAmpl = pnfit->getAmpl();
0320
0321     allPNAmpl.push_back(pnAmpl);
0322   }
0323
0324   // ===========
0325   // Get Matacq
0326   // ===========
0327
0328   ttrig = -1.;
0329   peak = -1.;
0330
0331   if (doesRefFileExist == 1) {
0332     // FIXME
0333     if (color == 0)
0334       matacqTree->GetEntry(event - 1);
0335     else if (color == 3)
0336       matacqTree->GetEntry(matacqTree->GetEntries("color==0") + event - 1);
0337     ttrig = ttMat;
0338     peak = peakMat;
0339   }
0340
0341   // ===========================
0342   // Decode EBDigis Information
0343   // ===========================
0344
0345   double yrange[10];
0346   int adcGain = 0;
0347   int side = 0;
0348
0349   if (EBDigi) {
0350     for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
0351          ++digiItr) {  // Loop on crystals
0352
0353       EBDetId id_crystal(digiItr->id());
0354       EBDataFrame df(*digiItr);
0355
0356       int etaG = id_crystal.ieta();  // global
0357       int phiG = id_crystal.iphi();  // global
0358
0359       int etaL;  // local
0360       int phiL;  // local
0361       std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0362
0363       etaL = LocalCoord.first;
0364       phiL = LocalCoord.second;
0365
0366       eta = etaG;
0367       phi = phiG;
0368
0369       side = MEEBGeom::side(etaG, phiG);
0370
0371       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0372
0373       int towerID = elecid_crystal.towerId();
0374       // int channelID=elecid_crystal.channelId()-1;  // FIXME so far for endcap only
0375       int strip = elecid_crystal.stripId();
0376       int xtal = elecid_crystal.xtalId();
0377       int channelID = 5 * (strip - 1) + xtal - 1;  // FIXME
0378
0379       int module = MEEBGeom::lmmod(etaG, phiG);
0380
0381       std::pair<int, int> pnpair = MEEBGeom::pn(module);
0382       unsigned int MyPn0 = pnpair.first;
0383       unsigned int MyPn1 = pnpair.second;
0384
0385       unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0386       assert(channel < nCrys);
0387
0388       double adcmax = 0.0;
0389
0390       if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
0391         continue;
0392       else
0393         channelNumber = channel;
0394
0395       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {  // Loop on adc samples
0396
0397         EcalMGPASample samp_crystal(df.sample(i));
0398         adc[i] = samp_crystal.adc();
0399         adcG[i] = samp_crystal.gainId();
0400
0401         if (i == 0)
0402           adcGain = adcG[i];
0403         if (i > 0)
0404           adcGain = TMath::Max(adcG[i], adcGain);
0405
0406         if (adc[i] > adcmax) {
0407           adcmax = adc[i];
0408         }
0409       }
0410
0411       for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0412         dsum += adc[k];
0413         if (k < _presample - 1)
0414           dsum1 += adc[k];
0415       }
0416
0417       bl = dsum / ((double)_presample);
0418       bl1 = dsum1 / ((double)_presample - 1);
0419
0420       for (val_max = 0., k = 0; k < _nsamples; k++) {
0421         yrange[k] = adc[k] - bl;
0422         if (yrange[k] > val_max) {
0423           val_max = yrange[k];
0424           samplemax = k;
0425         }
0426       }
0427
0428       if (samplemax == 4 || samplemax == 5) {
0429         for (k = 0; k < _nsamples; k++) {
0430           yrange[k] = yrange[k] + bl - bl1;
0431         }
0432       }
0433
0434       for (unsigned int k = 0; k < _nsamples; k++) {
0435         adc[k] = yrange[k];
0436       }
0437
0438       pn0 = allPNAmpl[MyPn0];
0439       pn1 = allPNAmpl[MyPn1];
0440
0441       if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
0442         ADCtrees->Fill();
0443     }
0444
0445   } else {
0446     for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end();
0447          ++digiItr) {  // Loop on crystals
0448
0449       EEDetId id_crystal(digiItr->id());
0450       EEDataFrame df(*digiItr);
0451
0452       phi = id_crystal.ix() - 1;
0453       eta = id_crystal.iy() - 1;
0454       side = 0;  // FIXME
0455
0456       // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
0457
0458       EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0459
0460       int towerID = elecid_crystal.towerId();
0461       int channelID = elecid_crystal.channelId() - 1;
0462
0463       int module = MEEEGeom::lmmod(phi, eta);
0464
0465       std::pair<int, int> pnpair = MEEEGeom::pn(module, _fedid);
0466       unsigned int MyPn0 = pnpair.first;
0467       unsigned int MyPn1 = pnpair.second;
0468
0469       unsigned int channel = MEEEGeom::crystal(phi, eta);
0470       assert(channel < nCrys);
0471
0472       double adcmax = 0.0;
0473
0474       if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
0475         continue;
0476       else
0477         channelNumber = channel;
0478
0479       for (unsigned int i = 0; i < (*digiItr).size(); ++i) {  // Loop on adc samples
0480
0481         EcalMGPASample samp_crystal(df.sample(i));
0482         adc[i] = samp_crystal.adc();
0483         adcG[i] = samp_crystal.gainId();
0484
0485         if (i == 0)
0486           adcGain = adcG[i];
0487         if (i > 0)
0488           adcGain = TMath::Max(adcG[i], adcGain);
0489
0490         if (adc[i] > adcmax) {
0491           adcmax = adc[i];
0492         }
0493       }
0494
0495       for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0496         dsum += adc[k];
0497         if (k < _presample - 1)
0498           dsum1 += adc[k];
0499       }
0500
0501       bl = dsum / ((double)_presample);
0502       bl1 = dsum1 / ((double)_presample - 1);
0503
0504       for (val_max = 0., k = 0; k < _nsamples; k++) {
0505         yrange[k] = adc[k] - bl;
0506         if (yrange[k] > val_max) {
0507           val_max = yrange[k];
0508           samplemax = k;
0509         }
0510       }
0511
0512       if (samplemax == 4 || samplemax == 5) {
0513         for (k = 0; k < _nsamples; k++) {
0514           yrange[k] = yrange[k] + bl - bl1;
0515         }
0516       }
0517
0518       for (unsigned int k = 0; k < _nsamples; k++) {
0519         adc[k] = yrange[k];
0520       }
0521
0522       pn0 = allPNAmpl[MyPn0];
0523       pn1 = allPNAmpl[MyPn1];
0524
0525       if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
0526         ADCtrees->Fill();
0527     }
0528   }
0529
0530 }  // analyze
0531
0532 //========================================================================
0533 void EcalPerEvtLaserAnalyzer::endJob() {
0534   //========================================================================
0535
0536   assert(colors.size() <= nColor);
0537   unsigned int nCol = colors.size();
0538
0539   ADCtrees->Write();
0540
0541   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0542       << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0543   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0544       << "\t+=+       Analyzing laser data: getting per event               +=+";
0545   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0546       << "\t+=+            APD Amplitudes and ADC samples                   +=+";
0547   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0548       << "\t+=+    for fed:" << _fedid << ", tower:" << _tower << ", and channel:" << _channel;
0549
0550   // Define temporary tree to save APD amplitudes
0551
0552   APDFile = new TFile(resfile.c_str(), "RECREATE");
0553
0554   int ieta, iphi, channelID, towerID, flag;
0555   double alpha, beta;
0556
0557   colors.push_back(color);
0558
0559   for (unsigned int i = 0; i < nCol; i++) {
0560     stringstream name1;
0561     name1 << "headerCol" << colors[i];
0562
0563     header[i] = new TTree(name1.str().c_str(), name1.str().c_str());
0564
0565     header[i]->Branch("alpha", &alpha, "alpha/D");
0566     header[i]->Branch("beta", &beta, "beta/D");
0567     header[i]->Branch("iphi", &iphi, "iphi/I");
0568     header[i]->Branch("ieta", &ieta, "ieta/I");
0569     header[i]->Branch("dccID", &dccID, "dccID/I");
0570     header[i]->Branch("towerID", &towerID, "towerID/I");
0571     header[i]->Branch("channelID", &channelID, "channelID/I");
0572
0573     header[i]->SetBranchAddress("alpha", &alpha);
0574     header[i]->SetBranchAddress("beta", &beta);
0575     header[i]->SetBranchAddress("iphi", &iphi);
0576     header[i]->SetBranchAddress("ieta", &ieta);
0577     header[i]->SetBranchAddress("dccID", &dccID);
0578     header[i]->SetBranchAddress("towerID", &towerID);
0579     header[i]->SetBranchAddress("channelID", &channelID);
0580   }
0581
0582   stringstream name2;
0583   name2 << "APDTree";
0584   APDtrees = new TTree(name2.str().c_str(), name2.str().c_str());
0585
0586   //List of branches
0587
0588   APDtrees->Branch("event", &event, "event/I");
0589   APDtrees->Branch("color", &color, "color/I");
0590   APDtrees->Branch("adc", &adc, "adc[10]/D");
0591   APDtrees->Branch("peakMatacq", &peak, "peak/D");
0592   APDtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
0593   APDtrees->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
0594   APDtrees->Branch("apdTime", &apdTime, "apdTime/D");
0595   APDtrees->Branch("flag", &flag, "flag/I");
0596   APDtrees->Branch("pn0", &pn0, "pn0/D");
0597   APDtrees->Branch("pn1", &pn1, "pn1/D");
0598
0599   APDtrees->SetBranchAddress("event", &event);
0600   APDtrees->SetBranchAddress("color", &color);
0601   APDtrees->SetBranchAddress("adc", adc);
0602   APDtrees->SetBranchAddress("peakMatacq", &peak);
0603   APDtrees->SetBranchAddress("ttrigMatacq", &ttrig);
0604   APDtrees->SetBranchAddress("apdAmpl", &apdAmpl);
0605   APDtrees->SetBranchAddress("apdTime", &apdTime);
0606   APDtrees->SetBranchAddress("flag", &flag);
0607   APDtrees->SetBranchAddress("pn0", &pn0);
0608   APDtrees->SetBranchAddress("pn1", &pn1);
0609
0610   // Retrieve alpha and beta for APD amplitudes calculation
0611
0612   TFile* alphaFile = new TFile(refalphabeta_.c_str());
0613   TTree* alphaTree[2];
0614
0615   Double_t alphaRun, betaRun;
0616   int ietaRun, iphiRun, channelIDRun, towerIDRun, dccIDRun, flagRun;
0617
0618   for (unsigned int i = 0; i < nCol; i++) {
0619     stringstream name3;
0620     name3 << "ABCol" << i;
0621     alphaTree[i] = (TTree*)alphaFile->Get(name3.str().c_str());
0622     alphaTree[i]->SetBranchStatus("*", false);
0623     alphaTree[i]->SetBranchStatus("alpha", true);
0624     alphaTree[i]->SetBranchStatus("beta", true);
0625     alphaTree[i]->SetBranchStatus("iphi", true);
0626     alphaTree[i]->SetBranchStatus("ieta", true);
0627     alphaTree[i]->SetBranchStatus("dccID", true);
0628     alphaTree[i]->SetBranchStatus("towerID", true);
0629     alphaTree[i]->SetBranchStatus("channelID", true);
0630     alphaTree[i]->SetBranchStatus("flag", true);
0631
0632     alphaTree[i]->SetBranchAddress("alpha", &alphaRun);
0633     alphaTree[i]->SetBranchAddress("beta", &betaRun);
0634     alphaTree[i]->SetBranchAddress("iphi", &iphiRun);
0635     alphaTree[i]->SetBranchAddress("ieta", &ietaRun);
0636     alphaTree[i]->SetBranchAddress("dccID", &dccIDRun);
0637     alphaTree[i]->SetBranchAddress("towerID", &towerIDRun);
0638     alphaTree[i]->SetBranchAddress("channelID", &channelIDRun);
0639     alphaTree[i]->SetBranchAddress("flag", &flagRun);
0640   }
0641
0642   PulseFitWithFunction* pslsfit = new PulseFitWithFunction();
0643
0644   double chi2;
0645
0646   for (unsigned int icol = 0; icol < nCol; icol++) {
0647     IsThereDataADC[icol] = 1;
0648     stringstream cut;
0649     cut << "color==" << colors.at(icol);
0650     if (ADCtrees->GetEntries(cut.str().c_str()) < 10)
0651       IsThereDataADC[icol] = 0;
0652     IsHeaderFilled[icol] = 0;
0653   }
0654
0655   // Define submodule and channel number inside the submodule (as Patrice)
0656
0657   for (Long64_t jentry = 0; jentry < ADCtrees->GetEntriesFast(); jentry++) {  // Loop on events
0658     ADCtrees->GetEntry(jentry);
0659
0660     int iCry = channelNumber;
0661
0662     // get back color
0663
0664     unsigned int iCol = 0;
0665     for (unsigned int i = 0; i < nCol; i++) {
0666       if (color == colors[i]) {
0667         iCol = i;
0668         i = colors.size();
0669       }
0670     }
0671
0672     alphaTree[iCol]->GetEntry(iCry);
0673
0674     flag = flagRun;
0675     iphi = iphiRun;
0676     ieta = ietaRun;
0677     towerID = towerIDRun;
0678     channelID = channelIDRun;
0679     alpha = alphaRun;
0680     beta = betaRun;
0681
0682     if (IsHeaderFilled[iCol] == 0) {
0683       header[iCol]->Fill();
0684       IsHeaderFilled[iCol] = 1;
0685     }
0686     // Amplitude calculation
0687
0688     apdAmpl = 0;
0689     apdTime = 0;
0690
0691     pslsfit->init(_nsamples, _firstsample, _lastsample, _niter, alphaRun, betaRun);
0692     chi2 = pslsfit->doFit(&adc[0]);
0693
0694     if (chi2 < 0. || chi2 == 102 || chi2 == 101) {
0695       apdAmpl = 0;
0696       apdTime = 0;
0697
0698     } else {
0699       apdAmpl = pslsfit->getAmpl();
0700       apdTime = pslsfit->getTime();
0701     }
0702
0703     APDtrees->Fill();
0704   }
0705
0706   alphaFile->Close();
0707
0708   ADCFile->Close();
0709
0710   APDFile->Write();
0711   APDFile->Close();
0712
0713   // Remove unwanted files
0714
0715   stringstream del;
0716   del << "rm " << ADCfile;
0717   system(del.str().c_str());
0718
0719   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0720       << "\t+=+    .................................................. done  +=+";
0721   edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0722       << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0723 }
0724
0725 DEFINE_FWK_MODULE(EcalPerEvtLaserAnalyzer);