Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57: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 
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   double max = 0;
0236 
0237   for (EcalMatacqDigiCollection::const_iterator it = matacqDigi->begin(); it != matacqDigi->end();
0238        ++it) {  // Loop on matacq channel
0239 
0240     //
0241     const EcalMatacqDigi& digis = *it;
0242 
0243     if (_debug == 1) {
0244       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- Inside digis -- digi size=" << digis.size();
0245     }
0246 
0247     if (digis.size() == 0)
0248       continue;
0249     else
0250       isThereMatacq = true;
0251 
0252     max = 0;
0253     maxsamp = 0;
0254     nsamples = digis.size();
0255     tt = digis.tTrig();
0256 
0257     for (int i = 0; i < digis.size(); ++i) {  // Loop on matacq samples
0258       matacq[i] = -digis.adcCount(i);
0259       if (matacq[i] > max) {
0260         max = matacq[i];
0261         maxsamp = i;
0262       }
0263     }
0264     if (_debug == 1) {
0265       edm::LogVerbatim("EcalMatacqAnalyzzer")
0266           << "-- debug test -- Inside digis -- nsamples=" << nsamples << ", max=" << max;
0267     }
0268   }
0269 
0270   if (_debug == 1)
0271     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- After digis -- Event: " << iEvent;
0272   tree->Fill();
0273 
0274 }  // analyze
0275 
0276 //========================================================================
0277 void EcalMatacqAnalyzer::endJob() {
0278   // Don't do anything if there is no events
0279   if (!isThereMatacq) {
0280     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0281     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+     WARNING! NO MATACQ        +=+";
0282     edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0283 
0284     // Remove temporary file
0285     FILE* test;
0286     test = fopen(sampfile.c_str(), "r");
0287     if (test) {
0288       std::stringstream del2;
0289       del2 << "rm " << sampfile;
0290       system(del2.str().c_str());
0291     }
0292     return;
0293   }
0294 
0295   assert(colors.size() <= nColor);
0296   unsigned int nCol = colors.size();
0297 
0298   for (unsigned int iCol = 0; iCol < nCol; iCol++) {
0299     for (unsigned int iSide = 0; iSide < nSide; iSide++) {
0300       MTQ[iCol][iSide] = new TMTQ();
0301     }
0302   }
0303 
0304   outFile = new TFile(outfile.c_str(), "RECREATE");
0305 
0306   TProfile* shapeMat = new TProfile("shapeLaser", "shapeLaser", _nsamplesshape, -0.5, double(_nsamplesshape) - 0.5);
0307   TProfile* shapeMatTmp = new TProfile(
0308       "shapeLaserTmp", "shapeLaserTmp", _timeaftmax + _timebefmax, -0.5, double(_timeaftmax + _timebefmax) - 0.5);
0309 
0310   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0311   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+     Analyzing MATACQ data     +=+";
0312 
0313   //
0314   // create output ntuple
0315   //
0316 
0317   mtqShape = new TTree("MatacqShape", "MatacqShape");
0318 
0319   // list of branches
0320   // keep Patrice's notations
0321 
0322   mtqShape->Branch("event", &event, "event/I");
0323   mtqShape->Branch("color", &color, "color/I");
0324   mtqShape->Branch("status", &status, "status/I");
0325   mtqShape->Branch("peak", &peak, "peak/D");
0326   mtqShape->Branch("sigma", &sigma, "sigma/D");
0327   mtqShape->Branch("fit", &fit, "fit/D");
0328   mtqShape->Branch("ampl", &ampl, "ampl/D");
0329   mtqShape->Branch("trise", &trise, "trise/D");
0330   mtqShape->Branch("fwhm", &fwhm, "fwhm/D");
0331   mtqShape->Branch("fw20", &fw20, "fw20/D");
0332   mtqShape->Branch("fw80", &fw80, "fw80/D");
0333   mtqShape->Branch("ped", &ped, "ped/D");
0334   mtqShape->Branch("pedsig", &pedsig, "pedsig/D");
0335   mtqShape->Branch("ttrig", &ttrig, "ttrig/D");
0336   mtqShape->Branch("sliding", &sliding, "sliding/D");
0337 
0338   mtqShape->SetBranchAddress("event", &event);
0339   mtqShape->SetBranchAddress("color", &color);
0340   mtqShape->SetBranchAddress("status", &status);
0341   mtqShape->SetBranchAddress("peak", &peak);
0342   mtqShape->SetBranchAddress("sigma", &sigma);
0343   mtqShape->SetBranchAddress("fit", &fit);
0344   mtqShape->SetBranchAddress("ampl", &ampl);
0345   mtqShape->SetBranchAddress("fwhm", &fwhm);
0346   mtqShape->SetBranchAddress("fw20", &fw20);
0347   mtqShape->SetBranchAddress("fw80", &fw80);
0348   mtqShape->SetBranchAddress("trise", &trise);
0349   mtqShape->SetBranchAddress("ped", &ped);
0350   mtqShape->SetBranchAddress("pedsig", &pedsig);
0351   mtqShape->SetBranchAddress("ttrig", &ttrig);
0352   mtqShape->SetBranchAddress("sliding", &sliding);
0353 
0354   unsigned int endsample;
0355   unsigned int presample;
0356 
0357   // loop over the entries of the tree
0358   TChain* fChain = (TChain*)tree;
0359   Long64_t nentries = fChain->GetEntriesFast();
0360 
0361   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0362     // load the event
0363     Long64_t ientry = fChain->LoadTree(jentry);
0364     if (ientry < 0)
0365       break;
0366     fChain->GetEntry(jentry);
0367 
0368     status = 0;
0369     peak = -1;
0370     sigma = 0;
0371     fit = -1;
0372     ampl = -1;
0373     trise = -1;
0374     ttrig = tt;
0375     fwhm = 0;
0376     fw20 = 0;
0377     fw80 = 0;
0378     ped = 0;
0379     pedsig = 0;
0380     sliding = 0;
0381 
0382     if (_debug == 1)
0383       edm::LogVerbatim("EcalMatacqAnalyzzer")
0384           << "-- debug test -- inside loop 1  -- jentry:" << jentry << " over nentries=" << nentries;
0385 
0386     // create the object for Matacq data analysis
0387 
0388     endsample = maxsamp + _nsamplesaftmax;
0389     presample = int(_presample * nsamples / 100.);
0390     TMatacq* mtq = new TMatacq(nsamples,
0391                                presample,
0392                                endsample,
0393                                _noiseCut,
0394                                _parabnbefmax,
0395                                _parabnaftmax,
0396                                _thres,
0397                                _lowlev,
0398                                _highlev,
0399                                _nevlasers,
0400                                _slide);
0401 
0402     if (_debug == 1)
0403       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 2  -- ";
0404 
0405     // analyse the Matacq data
0406     if (mtq->rawPulseAnalysis(nsamples, &matacq[0]) == 0) {
0407       status = 1;
0408       ped = mtq->getBaseLine();
0409       pedsig = mtq->getsigBaseLine();
0410 
0411       if (_debug == 1)
0412         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 3  -- ped:" << ped;
0413       if (mtq->findPeak() == 0) {
0414         peak = mtq->getTimpeak();
0415         sigma = mtq->getsigTimpeak();
0416       }
0417       if (_debug == 1)
0418         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 4  -- peak:" << peak;
0419       if (mtq->doFit() == 0) {
0420         fit = mtq->getTimax();
0421         ampl = mtq->getAmpl();
0422         fwhm = mtq->getFwhm();
0423         fw20 = mtq->getWidth20();
0424         fw80 = mtq->getWidth80();
0425         sliding = mtq->getSlide();
0426       }
0427       if (_debug == 1)
0428         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 4  -- ampl:" << ampl;
0429       if (mtq->compute_trise() == 0) {
0430         trise = mtq->getTrise();
0431       }
0432       if (_debug == 1)
0433         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 5  -- trise:" << trise;
0434     }
0435 
0436     if (_debug == 1)
0437       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 6  -- status:" << status;
0438 
0439     if (status == 1 && mtq->findPeak() == 0) {
0440       int firstS = int(peak - double(_timebefmax));
0441       int lastS = int(peak + double(_timeaftmax));
0442 
0443       // Fill histo if there are enough samples
0444       if (_debug == 1)
0445         edm::LogVerbatim("EcalMatacqAnalyzzer")
0446             << "-- debug test -- inside loop 7  -- firstS:" << firstS << ", nsamples:" << nsamples;
0447 
0448       if (firstS >= 0 && lastS <= nsamples) {
0449         for (int i = firstS; i < lastS; i++) {
0450           shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0451         }
0452 
0453       } else {  // else extrapolate
0454 
0455         int firstSBis;
0456 
0457         if (firstS < 0) {  // fill first bins with 0
0458 
0459           double thisped;
0460           thisped = (matacq[0] + matacq[1] + matacq[2] + matacq[4] + matacq[5]) / 5.0;
0461 
0462           for (int i = firstS; i < 0; i++) {
0463             shapeMatTmp->Fill(double(i) - firstS, thisped);
0464           }
0465           firstSBis = 0;
0466 
0467         } else {
0468           firstSBis = firstS;
0469         }
0470 
0471         if (lastS > nsamples) {
0472           for (int i = firstSBis; i < int(nsamples); i++) {
0473             shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0474           }
0475 
0476           //extrapolate with expo tail
0477 
0478           double expb = 0.998;
0479           double matacqval = expb * matacq[nsamples - 1];
0480 
0481           for (int i = nsamples; i < lastS; i++) {
0482             shapeMatTmp->Fill(double(i) - firstS, matacqval);
0483             matacqval *= expb;
0484           }
0485 
0486         } else {
0487           for (int i = firstSBis; i < lastS; i++) {
0488             shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
0489           }
0490         }
0491       }
0492     }
0493     if (_debug == 1)
0494       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 8";
0495 
0496     // get back color
0497 
0498     int iCol = nCol;
0499     for (unsigned int i = 0; i < nCol; i++) {
0500       if (color == colors[i]) {
0501         iCol = i;
0502         i = nCol;
0503       }
0504     }
0505     if (_debug == 1)
0506       edm::LogVerbatim("EcalMatacqAnalyzzer")
0507           << "-- debug test -- inside loop 8bis color:" << color << " iCol:" << iCol << " nCol:" << nCol;
0508 
0509     // fill TMTQ
0510 
0511     if (status == 1 && mtq->findPeak() == 0 && mtq->doFit() == 0 && mtq->compute_trise() == 0)
0512       MTQ[iCol][lightside]->addEntry(peak, sigma, fit, ampl, trise, fwhm, fw20, fw80, ped, pedsig, sliding);
0513 
0514     // fill the output tree
0515 
0516     if (_debug == 1)
0517       edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside loop 9";
0518     mtqShape->Fill();
0519 
0520     // clean up
0521     delete mtq;
0522   }
0523 
0524   if (_debug == 1)
0525     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- after loop ";
0526   sampFile->Close();
0527 
0528   double Peak[6], Sigma[6], Fit[6], Ampl[6], Trise[6], Fwhm[6], Fw20[6], Fw80[6], Ped[6], Pedsig[6], Sliding[6];
0529   int Side;
0530 
0531   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
0532     std::stringstream nametree;
0533     nametree << "MatacqCol" << colors[iColor];
0534     meanTree[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
0535     meanTree[iColor]->Branch("side", &Side, "Side/I");
0536     meanTree[iColor]->Branch("peak", &Peak, "Peak[6]/D");
0537     meanTree[iColor]->Branch("sigma", &Sigma, "Sigma[6]/D");
0538     meanTree[iColor]->Branch("fit", &Fit, "Fit[6]/D");
0539     meanTree[iColor]->Branch("ampl", &Ampl, "Ampl[6]/D");
0540     meanTree[iColor]->Branch("trise", &Trise, "Trise[6]/D");
0541     meanTree[iColor]->Branch("fwhm", &Fwhm, "Fwhm[6]/D");
0542     meanTree[iColor]->Branch("fw20", &Fw20, "Fw20[6]/D");
0543     meanTree[iColor]->Branch("fw80", &Fw80, "Fw80[6]/D");
0544     meanTree[iColor]->Branch("ped", &Ped, "Ped[6]/D");
0545     meanTree[iColor]->Branch("pedsig", &Pedsig, "Pedsig[6]/D");
0546     meanTree[iColor]->Branch("sliding", &Sliding, "Sliding[6]/D");
0547 
0548     meanTree[iColor]->SetBranchAddress("side", &Side);
0549     meanTree[iColor]->SetBranchAddress("peak", Peak);
0550     meanTree[iColor]->SetBranchAddress("sigma", Sigma);
0551     meanTree[iColor]->SetBranchAddress("fit", Fit);
0552     meanTree[iColor]->SetBranchAddress("ampl", Ampl);
0553     meanTree[iColor]->SetBranchAddress("fwhm", Fwhm);
0554     meanTree[iColor]->SetBranchAddress("fw20", Fw20);
0555     meanTree[iColor]->SetBranchAddress("fw80", Fw80);
0556     meanTree[iColor]->SetBranchAddress("trise", Trise);
0557     meanTree[iColor]->SetBranchAddress("ped", Ped);
0558     meanTree[iColor]->SetBranchAddress("pedsig", Pedsig);
0559     meanTree[iColor]->SetBranchAddress("sliding", Sliding);
0560   }
0561 
0562   for (unsigned int iCol = 0; iCol < nCol; iCol++) {
0563     for (unsigned int iSide = 0; iSide < nSides; iSide++) {
0564       Side = iSide;
0565       std::vector<double> val[TMTQ::nOutVar];
0566 
0567       for (int iVar = 0; iVar < TMTQ::nOutVar; iVar++) {
0568         val[iVar] = MTQ[iCol][iSide]->get(iVar);
0569 
0570         for (unsigned int i = 0; i < val[iVar].size(); i++) {
0571           switch (iVar) {
0572             case TMTQ::iPeak:
0573               Peak[i] = val[iVar].at(i);
0574               break;
0575             case TMTQ::iSigma:
0576               Sigma[i] = val[iVar].at(i);
0577               break;
0578             case TMTQ::iFit:
0579               Fit[i] = val[iVar].at(i);
0580               break;
0581             case TMTQ::iAmpl:
0582               Ampl[i] = val[iVar].at(i);
0583               break;
0584             case TMTQ::iFwhm:
0585               Fwhm[i] = val[iVar].at(i);
0586               break;
0587             case TMTQ::iFw20:
0588               Fw20[i] = val[iVar].at(i);
0589               break;
0590             case TMTQ::iFw80:
0591               Fw80[i] = val[iVar].at(i);
0592               break;
0593             case TMTQ::iTrise:
0594               Trise[i] = val[iVar].at(i);
0595               break;
0596             case TMTQ::iPed:
0597               Ped[i] = val[iVar].at(i);
0598               break;
0599             case TMTQ::iPedsig:
0600               Pedsig[i] = val[iVar].at(i);
0601               break;
0602             case TMTQ::iSlide:
0603               Sliding[i] = val[iVar].at(i);
0604               break;
0605           }
0606         }
0607       }
0608       meanTree[iCol]->Fill();
0609       if (_debug == 1)
0610         edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- inside final loop  ";
0611     }
0612   }
0613 
0614   // Calculate maximum with pol 2
0615 
0616   int im = shapeMatTmp->GetMaximumBin();
0617   double q1 = shapeMatTmp->GetBinContent(im - 1);
0618   double q2 = shapeMatTmp->GetBinContent(im);
0619   double q3 = shapeMatTmp->GetBinContent(im + 1);
0620 
0621   double a2 = (q3 + q1) / 2.0 - q2;
0622   double a1 = q2 - q1 + a2 * (1 - 2 * im);
0623   double a0 = q2 - a1 * im - a2 * im * im;
0624 
0625   double am = a0 - a1 * a1 / (4 * a2);
0626 
0627   // Compute pedestal
0628 
0629   double bl = 0;
0630   for (unsigned int i = 1; i < _presampleshape + 1; i++) {
0631     bl += shapeMatTmp->GetBinContent(i);
0632   }
0633   bl /= _presampleshape;
0634 
0635   // Compute and save laser shape
0636 
0637   if (_debug == 1)
0638     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- computing shape  ";
0639 
0640   int firstBin = 0;
0641   double height = 0.0;
0642 
0643   for (unsigned int i = _timebefmax; i > _presampleshape; i--) {
0644     height = shapeMatTmp->GetBinContent(i) - bl;
0645 
0646     if (height < (am - bl) * _cutwindow) {
0647       firstBin = i;
0648       i = _presampleshape;
0649     }
0650   }
0651 
0652   unsigned int lastBin = firstBin + _nsamplesshape;
0653 
0654   for (unsigned int i = firstBin; i < lastBin; i++) {
0655     shapeMat->Fill(i - firstBin, shapeMatTmp->GetBinContent(i) - bl);
0656   }
0657 
0658   mtqShape->Write();
0659   for (unsigned int iColor = 0; iColor < nCol; iColor++) {
0660     meanTree[iColor]->Write();
0661   }
0662   if (_debug == 1)
0663     edm::LogVerbatim("EcalMatacqAnalyzzer") << "-- debug test -- writing  ";
0664   shapeMat->Write();
0665 
0666   // close the output file
0667   outFile->Close();
0668 
0669   // Remove temporary file
0670   FILE* test;
0671   test = fopen(sampfile.c_str(), "r");
0672   if (test) {
0673     std::stringstream del2;
0674     del2 << "rm " << sampfile;
0675     system(del2.str().c_str());
0676   }
0677 
0678   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+    .................... done  +=+";
0679   edm::LogVerbatim("EcalMatacqAnalyzzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0680 }
0681 
0682 DEFINE_FWK_MODULE(EcalMatacqAnalyzer);
[ 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!