Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-12-23 23:23:31

0001 /* 
0002  *  \class EcalMatacqAnalyzer
0003  *
0004  *  primary author: Gautier Hamel De Monchenault - CEA/Saclay
0005  *  author: Julie Malcles - CEA/Saclay
0006  */
0007 
0008 #include <TFile.h>
0009 #include <TTree.h>
0010 #include <TProfile.h>
0011 #include <TChain.h>
0012 #include <vector>
0013 
0014 #include <CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalMatacqAnalyzer.h>
0015 
0016 #include <sstream>
0017 #include <iostream>
0018 #include <iomanip>
0019 
0020 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0021 #include <FWCore/Utilities/interface/Exception.h>
0022 
0023 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
0024 
0025 #include <FWCore/Framework/interface/Event.h>
0026 #include <FWCore/Framework/interface/MakerMacros.h>
0027 #include <FWCore/ParameterSet/interface/ParameterSet.h>
0028 
0029 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMatacq.h>
0030 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
0031 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMTQ.h>
0032 
0033 //========================================================================
0034 EcalMatacqAnalyzer::EcalMatacqAnalyzer(const edm::ParameterSet& iConfig)
0035     //========================================================================
0036     :
0037 
0038       iEvent(0),
0039 
0040       // framework parameters with default values
0041 
0042       _presample(iConfig.getUntrackedParameter<double>("nPresamples", 6.7)),
0043       _nsamplesaftmax(iConfig.getUntrackedParameter<unsigned int>("nSamplesAftMax", 80)),
0044       _noiseCut(iConfig.getUntrackedParameter<unsigned int>("noiseCut", 7)),
0045       _parabnbefmax(iConfig.getUntrackedParameter<unsigned int>("paraBeforeMax", 8)),
0046       _parabnaftmax(iConfig.getUntrackedParameter<unsigned int>("paraAfterMax", 7)),
0047       _thres(iConfig.getUntrackedParameter<unsigned int>("threshold", 10)),
0048       _lowlev(iConfig.getUntrackedParameter<unsigned int>("lowLevel", 20)),
0049       _highlev(iConfig.getUntrackedParameter<unsigned int>("highLevel", 80)),
0050       _nevlasers(iConfig.getUntrackedParameter<unsigned int>("nEventLaser", 600)),
0051       _timebefmax(iConfig.getUntrackedParameter<unsigned int>("timeBefMax", 100)),
0052       _timeaftmax(iConfig.getUntrackedParameter<unsigned int>("timeAftMax", 250)),
0053       _cutwindow(iConfig.getUntrackedParameter<double>("cutWindow", 0.1)),
0054       _nsamplesshape(iConfig.getUntrackedParameter<unsigned int>("nSamplesShape", 250)),
0055       _presampleshape(iConfig.getUntrackedParameter<unsigned int>("nPresamplesShape", 50)),
0056       _slide(iConfig.getUntrackedParameter<unsigned int>("nSlide", 100)),
0057       _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
0058       _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
0059       resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0060       digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0061       digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0062       eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0063       eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0064       pmatToken_(consumes<EcalMatacqDigiCollection>(edm::InputTag(digiProducer_, digiCollection_))),
0065       dccToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0066       nSides(NSIDES),
0067       lightside(0),
0068       runType(-1),
0069       runNum(0),
0070       event(0),
0071       color(-1),
0072       maxsamp(0),
0073       nsamples(0),
0074       tt(0)
0075 
0076 //========================================================================
0077 {
0078   //now do what ever initialization is needed
0079 }
0080 
0081 //========================================================================
0082 void EcalMatacqAnalyzer::beginJob() {
0083   //========================================================================
0084 
0085   // Define temporary file name
0086 
0087   sampfile = resdir_;
0088   sampfile += "/TmpTreeMatacqAnalyzer.root";
0089 
0090   sampFile = new TFile(sampfile.c_str(), "RECREATE");
0091 
0092   // declaration of the tree to fill
0093 
0094   tree = new TTree("MatacqTree", "MatacqTree");
0095 
0096   //List of branches
0097 
0098   tree->Branch("event", &event, "event/I");
0099   tree->Branch("color", &color, "color/I");
0100   tree->Branch("matacq", &matacq, "matacq[2560]/D");
0101   tree->Branch("nsamples", &nsamples, "nsamples/I");
0102   tree->Branch("maxsamp", &maxsamp, "maxsamp/I");
0103   tree->Branch("tt", &tt, "tt/D");
0104   tree->Branch("lightside", &lightside, "lightside/I");
0105 
0106   tree->SetBranchAddress("event", &event);
0107   tree->SetBranchAddress("color", &color);
0108   tree->SetBranchAddress("matacq", matacq);
0109   tree->SetBranchAddress("nsamples", &nsamples);
0110   tree->SetBranchAddress("maxsamp", &maxsamp);
0111   tree->SetBranchAddress("tt", &tt);
0112   tree->SetBranchAddress("lightside", &lightside);
0113 
0114   // Define output results files' names
0115 
0116   std::stringstream namefile;
0117   namefile << resdir_ << "/MATACQ.root";
0118   outfile = namefile.str();
0119 
0120   // Laser events counter
0121   laserEvents = 0;
0122   isThereMatacq = false;
0123 }
0124 
0125 //========================================================================
0126 void EcalMatacqAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0127   //========================================================================
0128 
0129   ++iEvent;
0130 
0131   if (_debug == 1)
0132     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Entering Analyze -- event= " << iEvent;
0133 
0134   // retrieving MATACQ :
0135   const edm::Handle<EcalMatacqDigiCollection>& pmatacqDigi = e.getHandle(pmatToken_);
0136   const EcalMatacqDigiCollection* matacqDigi = (pmatacqDigi.isValid()) ? pmatacqDigi.product() : nullptr;
0137   if (pmatacqDigi.isValid()) {
0138     if (_debug == 1)
0139       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Matacq Digis Found -- ";
0140   } else {
0141     edm::LogError("EcalMatacqAnalyzzer") << "Error! can't get the product EcalMatacqDigi producer:"
0142                                          << digiProducer_.c_str() << " collection:" << digiCollection_.c_str();
0143     return;
0144   }
0145 
0146   // retrieving DCC header
0147 
0148   const edm::Handle<EcalRawDataCollection>& pDCCHeader = e.getHandle(dccToken_);
0149   const EcalRawDataCollection* DCCHeader = (pDCCHeader.isValid()) ? pDCCHeader.product() : nullptr;
0150   if (!pDCCHeader.isValid()) {
0151     edm::LogError("EcalMatacqAnalyzzer") << "Error! can't get the product EcalRawData producer:"
0152                                          << eventHeaderProducer_.c_str()
0153                                          << " collection:" << eventHeaderCollection_.c_str();
0154     return;
0155   }
0156 
0157   // ====================================
0158   // Decode Basic DCCHeader Information
0159   // ====================================
0160 
0161   if (_debug == 1)
0162     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Before header -- ";
0163 
0164   for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0165        ++headerItr) {
0166     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0167     color = (int)settings.wavelength;
0168     if (color < 0)
0169       return;
0170 
0171     // Get run type and run number
0172 
0173     int fed = headerItr->fedId();
0174 
0175     if (fed != _fedid && _fedid != -999)
0176       continue;
0177 
0178     runType = headerItr->getRunType();
0179     runNum = headerItr->getRunNumber();
0180     event = headerItr->getLV1();
0181 
0182     if (_debug == 1)
0183       edm::LogVerbatim("EcalMatacqAnalyzzer")
0184           << "-- debug test -- runtype:" << runType << " event:" << event << " runNum:" << runNum;
0185 
0186     dccID = headerItr->getDccInTCCCommand();
0187     fedID = headerItr->fedId();
0188     lightside = headerItr->getRtHalf();
0189 
0190     //assert (lightside<2 && lightside>=0);
0191 
0192     if (lightside != 1 && lightside != 0) {
0193       edm::LogVerbatim("EcalMatacqAnalyzzer") << "Unexpected lightside: " << lightside << " for event " << iEvent;
0194       return;
0195     }
0196     if (_debug == 1) {
0197       edm::LogVerbatim("EcalMatacqAnalyzzer")
0198           << "-- debug test -- Inside header before fed cut -- color=" << color << ", dcc=" << dccID
0199           << ", fed=" << fedID << ",  lightside=" << lightside << ", runType=" << runType;
0200     }
0201 
0202     // take event only if the fed corresponds to the DCC in TCC
0203     if (600 + dccID != fedID)
0204       continue;
0205 
0206     if (_debug == 1) {
0207       edm::LogVerbatim("EcalMatacqAnalyzzer")
0208           << "-- debug test -- Inside header after fed cut -- color=" << color << ", dcc=" << dccID << ", fed=" << fedID
0209           << ",  lightside=" << lightside << ", runType=" << runType;
0210     }
0211 
0212     // Cut on runType
0213 
0214     if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
0215         runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
0216       return;
0217 
0218     std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0219     if (iter == colors.end()) {
0220       colors.push_back(color);
0221       edm::LogVerbatim("EcalMatacqAnalyzzer") << " new color found " << color << " " << colors.size();
0222     }
0223   }
0224 
0225   if (_debug == 1)
0226     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Before digis -- Event:" << iEvent;
0227 
0228   // Count laser events
0229   laserEvents++;
0230 
0231   // ===========================
0232   // Decode Matacq Information
0233   // ===========================
0234 
0235   int iCh = 0;
0236   double max = 0;
0237 
0238   for (EcalMatacqDigiCollection::const_iterator it = matacqDigi->begin(); it != matacqDigi->end();
0239        ++it) {  // Loop on matacq channel
0240 
0241     //
0242     const EcalMatacqDigi& digis = *it;
0243 
0244     //if(digis.size()==0 || iCh>=N_channels) continue;
0245     if (_debug == 1) {
0246       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Inside digis -- digi size=" << digis.size();
0247     }
0248 
0249     if (digis.size() == 0)
0250       continue;
0251     else
0252       isThereMatacq = true;
0253 
0254     max = 0;
0255     maxsamp = 0;
0256     nsamples = digis.size();
0257     tt = digis.tTrig();
0258 
0259     for (int i = 0; i < digis.size(); ++i) {  // Loop on matacq samples
0260       matacq[i] = -digis.adcCount(i);
0261       if (matacq[i] > max) {
0262         max = matacq[i];
0263         maxsamp = i;
0264       }
0265     }
0266     if (_debug == 1) {
0267       edm::LogVerbatim("EcalMatacqAnalyzzer")
0268           << "-- debug test -- Inside digis -- nsamples=" << nsamples << ", max=" << max;
0269     }
0270 
0271     iCh++;
0272   }
0273 
0274   if (_debug == 1)
0275     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- After digis -- Event: " << iEvent;
0276   tree->Fill();
0277 
0278 }  // analyze
0279 
0280 //========================================================================
0281 void EcalMatacqAnalyzer::endJob() {
0282   // Don't do anything if there is no events
0283   if (!isThereMatacq) {
0284     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0285     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+     WARNING! NO MATACQ        +=+";
0286     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0287 
0288     // Remove temporary file
0289     FILE* test;
0290     test = fopen(sampfile.c_str(), "r");
0291     if (test) {
0292       std::stringstream del2;
0293       del2 << "rm " << sampfile;
0294       system(del2.str().c_str());
0295     }
0296     return;
0297   }
0298 
0299   assert(colors.size() <= nColor);
0300   unsigned int nCol = colors.size();
0301 
0302   for (unsigned int iCol = 0; iCol < nCol; iCol++) {
0303     for (unsigned int iSide = 0; iSide < nSide; iSide++) {
0304       MTQ[iCol][iSide] = new TMTQ();
0305     }
0306   }
0307 
0308   outFile = new TFile(outfile.c_str(), "RECREATE");
0309 
0310   TProfile* shapeMat = new TProfile("shapeLaser", "shapeLaser", _nsamplesshape, -0.5, double(_nsamplesshape) - 0.5);
0311   TProfile* shapeMatTmp = new TProfile(
0312       "shapeLaserTmp", "shapeLaserTmp", _timeaftmax + _timebefmax, -0.5, double(_timeaftmax + _timebefmax) - 0.5);
0313 
0314   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0315   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+     Analyzing MATACQ data     +=+";
0316 
0317   //
0318   // create output ntuple
0319   //
0320 
0321   mtqShape = new TTree("MatacqShape", "MatacqShape");
0322 
0323   // list of branches
0324   // keep Patrice's notations
0325 
0326   mtqShape->Branch("event", &event, "event/I");
0327   mtqShape->Branch("color", &color, "color/I");
0328   mtqShape->Branch("status", &status, "status/I");
0329   mtqShape->Branch("peak", &peak, "peak/D");
0330   mtqShape->Branch("sigma", &sigma, "sigma/D");
0331   mtqShape->Branch("fit", &fit, "fit/D");
0332   mtqShape->Branch("ampl", &ampl, "ampl/D");
0333   mtqShape->Branch("trise", &trise, "trise/D");
0334   mtqShape->Branch("fwhm", &fwhm, "fwhm/D");
0335   mtqShape->Branch("fw20", &fw20, "fw20/D");
0336   mtqShape->Branch("fw80", &fw80, "fw80/D");
0337   mtqShape->Branch("ped", &ped, "ped/D");
0338   mtqShape->Branch("pedsig", &pedsig, "pedsig/D");
0339   mtqShape->Branch("ttrig", &ttrig, "ttrig/D");
0340   mtqShape->Branch("sliding", &sliding, "sliding/D");
0341 
0342   mtqShape->SetBranchAddress("event", &event);
0343   mtqShape->SetBranchAddress("color", &color);
0344   mtqShape->SetBranchAddress("status", &status);
0345   mtqShape->SetBranchAddress("peak", &peak);
0346   mtqShape->SetBranchAddress("sigma", &sigma);
0347   mtqShape->SetBranchAddress("fit", &fit);
0348   mtqShape->SetBranchAddress("ampl", &ampl);
0349   mtqShape->SetBranchAddress("fwhm", &fwhm);
0350   mtqShape->SetBranchAddress("fw20", &fw20);
0351   mtqShape->SetBranchAddress("fw80", &fw80);
0352   mtqShape->SetBranchAddress("trise", &trise);
0353   mtqShape->SetBranchAddress("ped", &ped);
0354   mtqShape->SetBranchAddress("pedsig", &pedsig);
0355   mtqShape->SetBranchAddress("ttrig", &ttrig);
0356   mtqShape->SetBranchAddress("sliding", &sliding);
0357 
0358   unsigned int endsample;
0359   unsigned int presample;
0360 
0361   // loop over the entries of the tree
0362   TChain* fChain = (TChain*)tree;
0363   Long64_t nentries = fChain->GetEntriesFast();
0364 
0365   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0366     // load the event
0367     Long64_t ientry = fChain->LoadTree(jentry);
0368     if (ientry < 0)
0369       break;
0370     fChain->GetEntry(jentry);
0371 
0372     status = 0;
0373     peak = -1;
0374     sigma = 0;
0375     fit = -1;
0376     ampl = -1;
0377     trise = -1;
0378     ttrig = tt;
0379     fwhm = 0;
0380     fw20 = 0;
0381     fw80 = 0;
0382     ped = 0;
0383     pedsig = 0;
0384     sliding = 0;
0385 
0386     if (_debug == 1)
0387       edm::LogVerbatim("EcalMatacqAnalyzzer")
0388           << "-- debug test -- inside loop 1  -- jentry:" << jentry << " over nentries=" << nentries;
0389 
0390     // create the object for Matacq data analysis
0391 
0392     endsample = maxsamp + _nsamplesaftmax;
0393     presample = int(_presample * nsamples / 100.);
0394     TMatacq* mtq = new TMatacq(nsamples,
0395                                presample,
0396                                endsample,
0397                                _noiseCut,
0398                                _parabnbefmax,
0399                                _parabnaftmax,
0400                                _thres,
0401                                _lowlev,
0402                                _highlev,
0403                                _nevlasers,
0404                                _slide);
0405 
0406     if (_debug == 1)
0407       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 2  -- ";
0408 
0409     // analyse the Matacq data
0410     if (mtq->rawPulseAnalysis(nsamples, &matacq[0]) == 0) {
0411       status = 1;
0412       ped = mtq->getBaseLine();
0413       pedsig = mtq->getsigBaseLine();
0414 
0415       if (_debug == 1)
0416         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 3  -- ped:" << ped;
0417       if (mtq->findPeak() == 0) {
0418         peak = mtq->getTimpeak();
0419         sigma = mtq->getsigTimpeak();
0420       }
0421       if (_debug == 1)
0422         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 4  -- peak:" << peak;
0423       if (mtq->doFit() == 0) {
0424         fit = mtq->getTimax();
0425         ampl = mtq->getAmpl();
0426         fwhm = mtq->getFwhm();
0427         fw20 = mtq->getWidth20();
0428         fw80 = mtq->getWidth80();
0429         sliding = mtq->getSlide();
0430       }
0431       if (_debug == 1)
0432         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 4  -- ampl:" << ampl;
0433       if (mtq->compute_trise() == 0) {
0434         trise = mtq->getTrise();
0435       }
0436       if (_debug == 1)
0437         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 5  -- trise:" << trise;
0438     }
0439 
0440     if (_debug == 1)
0441       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 6  -- status:" << status;
0442 
0443     if (status == 1 && mtq->findPeak() == 0) {
0444       int firstS = int(peak - double(_timebefmax));
0445       int lastS = int(peak + double(_timeaftmax));
0446 
0447       // Fill histo if there are enough samples
0448       if (_debug == 1)
0449         edm::LogVerbatim("EcalMatacqAnalyzzer")
0450             << "-- debug test -- inside loop 7  -- firstS:" << firstS << ", nsamples:" << nsamples;
0451 
0452       if (firstS >= 0 && lastS <= nsamples) {
0453         for (int i = firstS; i < lastS; i++) {
0454           shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0455         }
0456 
0457       } else {  // else extrapolate
0458 
0459         int firstSBis;
0460 
0461         if (firstS < 0) {  // fill first bins with 0
0462 
0463           double thisped;
0464           thisped = (matacq[0] + matacq[1] + matacq[2] + matacq[4] + matacq[5]) / 5.0;
0465 
0466           for (int i = firstS; i < 0; i++) {
0467             shapeMatTmp->Fill(double(i) - firstS, thisped);
0468           }
0469           firstSBis = 0;
0470 
0471         } else {
0472           firstSBis = firstS;
0473         }
0474 
0475         if (lastS > nsamples) {
0476           for (int i = firstSBis; i < int(nsamples); i++) {
0477             shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0478           }
0479 
0480           //extrapolate with expo tail
0481 
0482           double expb = 0.998;
0483           double matacqval = expb * matacq[nsamples - 1];
0484 
0485           for (int i = nsamples; i < lastS; i++) {
0486             shapeMatTmp->Fill(double(i) - firstS, matacqval);
0487             matacqval *= expb;
0488           }
0489 
0490         } else {
0491           for (int i = firstSBis; i < lastS; i++) {
0492             shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0493           }
0494         }
0495       }
0496     }
0497     if (_debug == 1)
0498       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 8";
0499 
0500     // get back color
0501 
0502     int iCol = nCol;
0503     for (unsigned int i = 0; i < nCol; i++) {
0504       if (color == colors[i]) {
0505         iCol = i;
0506         i = nCol;
0507       }
0508     }
0509     if (_debug == 1)
0510       edm::LogVerbatim("EcalMatacqAnalyzzer")
0511           << "-- debug test -- inside loop 8bis color:" << color << " iCol:" << iCol << " nCol:" << nCol;
0512 
0513     // fill TMTQ
0514 
0515     if (status == 1 && mtq->findPeak() == 0 && mtq->doFit() == 0 && mtq->compute_trise() == 0)
0516       MTQ[iCol][lightside]->addEntry(peak, sigma, fit, ampl, trise, fwhm, fw20, fw80, ped, pedsig, sliding);
0517 
0518     // fill the output tree
0519 
0520     if (_debug == 1)
0521       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 9";
0522     mtqShape->Fill();
0523 
0524     // clean up
0525     delete mtq;
0526   }
0527 
0528   if (_debug == 1)
0529     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- after loop ";
0530   sampFile->Close();
0531 
0532   double Peak[6], Sigma[6], Fit[6], Ampl[6], Trise[6], Fwhm[6], Fw20[6], Fw80[6], Ped[6], Pedsig[6], Sliding[6];
0533   int Side;
0534 
0535   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
0536     std::stringstream nametree;
0537     nametree << "MatacqCol" << colors[iColor];
0538     meanTree[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
0539     meanTree[iColor]->Branch("side", &Side, "Side/I");
0540     meanTree[iColor]->Branch("peak", &Peak, "Peak[6]/D");
0541     meanTree[iColor]->Branch("sigma", &Sigma, "Sigma[6]/D");
0542     meanTree[iColor]->Branch("fit", &Fit, "Fit[6]/D");
0543     meanTree[iColor]->Branch("ampl", &Ampl, "Ampl[6]/D");
0544     meanTree[iColor]->Branch("trise", &Trise, "Trise[6]/D");
0545     meanTree[iColor]->Branch("fwhm", &Fwhm, "Fwhm[6]/D");
0546     meanTree[iColor]->Branch("fw20", &Fw20, "Fw20[6]/D");
0547     meanTree[iColor]->Branch("fw80", &Fw80, "Fw80[6]/D");
0548     meanTree[iColor]->Branch("ped", &Ped, "Ped[6]/D");
0549     meanTree[iColor]->Branch("pedsig", &Pedsig, "Pedsig[6]/D");
0550     meanTree[iColor]->Branch("sliding", &Sliding, "Sliding[6]/D");
0551 
0552     meanTree[iColor]->SetBranchAddress("side", &Side);
0553     meanTree[iColor]->SetBranchAddress("peak", Peak);
0554     meanTree[iColor]->SetBranchAddress("sigma", Sigma);
0555     meanTree[iColor]->SetBranchAddress("fit", Fit);
0556     meanTree[iColor]->SetBranchAddress("ampl", Ampl);
0557     meanTree[iColor]->SetBranchAddress("fwhm", Fwhm);
0558     meanTree[iColor]->SetBranchAddress("fw20", Fw20);
0559     meanTree[iColor]->SetBranchAddress("fw80", Fw80);
0560     meanTree[iColor]->SetBranchAddress("trise", Trise);
0561     meanTree[iColor]->SetBranchAddress("ped", Ped);
0562     meanTree[iColor]->SetBranchAddress("pedsig", Pedsig);
0563     meanTree[iColor]->SetBranchAddress("sliding", Sliding);
0564   }
0565 
0566   for (unsigned int iCol = 0; iCol < nCol; iCol++) {
0567     for (unsigned int iSide = 0; iSide < nSides; iSide++) {
0568       Side = iSide;
0569       std::vector<double> val[TMTQ::nOutVar];
0570 
0571       for (int iVar = 0; iVar < TMTQ::nOutVar; iVar++) {
0572         val[iVar] = MTQ[iCol][iSide]->get(iVar);
0573 
0574         for (unsigned int i = 0; i < val[iVar].size(); i++) {
0575           switch (iVar) {
0576             case TMTQ::iPeak:
0577               Peak[i] = val[iVar].at(i);
0578               break;
0579             case TMTQ::iSigma:
0580               Sigma[i] = val[iVar].at(i);
0581               break;
0582             case TMTQ::iFit:
0583               Fit[i] = val[iVar].at(i);
0584               break;
0585             case TMTQ::iAmpl:
0586               Ampl[i] = val[iVar].at(i);
0587               break;
0588             case TMTQ::iFwhm:
0589               Fwhm[i] = val[iVar].at(i);
0590               break;
0591             case TMTQ::iFw20:
0592               Fw20[i] = val[iVar].at(i);
0593               break;
0594             case TMTQ::iFw80:
0595               Fw80[i] = val[iVar].at(i);
0596               break;
0597             case TMTQ::iTrise:
0598               Trise[i] = val[iVar].at(i);
0599               break;
0600             case TMTQ::iPed:
0601               Ped[i] = val[iVar].at(i);
0602               break;
0603             case TMTQ::iPedsig:
0604               Pedsig[i] = val[iVar].at(i);
0605               break;
0606             case TMTQ::iSlide:
0607               Sliding[i] = val[iVar].at(i);
0608               break;
0609           }
0610         }
0611       }
0612       meanTree[iCol]->Fill();
0613       if (_debug == 1)
0614         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside final loop  ";
0615     }
0616   }
0617 
0618   // Calculate maximum with pol 2
0619 
0620   int im = shapeMatTmp->GetMaximumBin();
0621   double q1 = shapeMatTmp->GetBinContent(im - 1);
0622   double q2 = shapeMatTmp->GetBinContent(im);
0623   double q3 = shapeMatTmp->GetBinContent(im + 1);
0624 
0625   double a2 = (q3 + q1) / 2.0 - q2;
0626   double a1 = q2 - q1 + a2 * (1 - 2 * im);
0627   double a0 = q2 - a1 * im - a2 * im * im;
0628 
0629   double am = a0 - a1 * a1 / (4 * a2);
0630 
0631   // Compute pedestal
0632 
0633   double bl = 0;
0634   for (unsigned int i = 1; i < _presampleshape + 1; i++) {
0635     bl += shapeMatTmp->GetBinContent(i);
0636   }
0637   bl /= _presampleshape;
0638 
0639   // Compute and save laser shape
0640 
0641   if (_debug == 1)
0642     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- computing shape  ";
0643 
0644   int firstBin = 0;
0645   double height = 0.0;
0646 
0647   for (unsigned int i = _timebefmax; i > _presampleshape; i--) {
0648     height = shapeMatTmp->GetBinContent(i) - bl;
0649 
0650     if (height < (am - bl) * _cutwindow) {
0651       firstBin = i;
0652       i = _presampleshape;
0653     }
0654   }
0655 
0656   unsigned int lastBin = firstBin + _nsamplesshape;
0657 
0658   for (unsigned int i = firstBin; i < lastBin; i++) {
0659     shapeMat->Fill(i - firstBin, shapeMatTmp->GetBinContent(i) - bl);
0660   }
0661 
0662   mtqShape->Write();
0663   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
0664     meanTree[iColor]->Write();
0665   }
0666   if (_debug == 1)
0667     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- writing  ";
0668   shapeMat->Write();
0669 
0670   // close the output file
0671   outFile->Close();
0672 
0673   // Remove temporary file
0674   FILE* test;
0675   test = fopen(sampfile.c_str(), "r");
0676   if (test) {
0677     std::stringstream del2;
0678     del2 << "rm " << sampfile;
0679     system(del2.str().c_str());
0680   }
0681 
0682   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+    .................... done  +=+";
0683   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0684 }
0685 
0686 DEFINE_FWK_MODULE(EcalMatacqAnalyzer);