Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:46:49

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