CMS 3D CMS Logo
EcalMatacqAnalyzer.cc
Go to the documentation of this file.
1 /*
2  * \class EcalMatacqAnalyzer
3  *
4  * primary author: Gautier Hamel De Monchenault - CEA/Saclay
5  * author: Julie Malcles - CEA/Saclay
6  */
7 
8 #include <TFile.h>
9 #include <TTree.h>
10 #include <TProfile.h>
11 #include <TChain.h>
12 #include <vector>
13 
14 #include <CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalMatacqAnalyzer.h>
15 
16 #include <sstream>
17 #include <iostream>
18 #include <iomanip>
19 
20 #include <FWCore/MessageLogger/interface/MessageLogger.h>
21 #include <FWCore/Utilities/interface/Exception.h>
22 
23 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
24 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
25 
26 #include <FWCore/Framework/interface/Event.h>
27 #include <FWCore/Framework/interface/MakerMacros.h>
28 #include <FWCore/ParameterSet/interface/ParameterSet.h>
29 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
30 
31 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMatacq.h>
32 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
33 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMTQ.h>
34 
35 using namespace std;
36 
37 //========================================================================
38 EcalMatacqAnalyzer::EcalMatacqAnalyzer(const edm::ParameterSet& iConfig)
39  //========================================================================
40  :
41 
42  iEvent(0),
43 
44  // framework parameters with default values
45 
46  _presample(iConfig.getUntrackedParameter<double>("nPresamples", 6.7)),
47  _nsamplesaftmax(iConfig.getUntrackedParameter<unsigned int>("nSamplesAftMax", 80)),
48  _noiseCut(iConfig.getUntrackedParameter<unsigned int>("noiseCut", 7)),
49  _parabnbefmax(iConfig.getUntrackedParameter<unsigned int>("paraBeforeMax", 8)),
50  _parabnaftmax(iConfig.getUntrackedParameter<unsigned int>("paraAfterMax", 7)),
51  _thres(iConfig.getUntrackedParameter<unsigned int>("threshold", 10)),
52  _lowlev(iConfig.getUntrackedParameter<unsigned int>("lowLevel", 20)),
53  _highlev(iConfig.getUntrackedParameter<unsigned int>("highLevel", 80)),
54  _nevlasers(iConfig.getUntrackedParameter<unsigned int>("nEventLaser", 600)),
55  _timebefmax(iConfig.getUntrackedParameter<unsigned int>("timeBefMax", 100)),
56  _timeaftmax(iConfig.getUntrackedParameter<unsigned int>("timeAftMax", 250)),
57  _cutwindow(iConfig.getUntrackedParameter<double>("cutWindow", 0.1)),
58  _nsamplesshape(iConfig.getUntrackedParameter<unsigned int>("nSamplesShape", 250)),
59  _presampleshape(iConfig.getUntrackedParameter<unsigned int>("nPresamplesShape", 50)),
60  _slide(iConfig.getUntrackedParameter<unsigned int>("nSlide", 100)),
61  _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
62  _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
63  nSides(NSIDES),
64  lightside(0),
65  runType(-1),
66  runNum(0),
67  event(0),
68  color(-1),
69  maxsamp(0),
70  nsamples(0),
71  tt(0)
72 
73 //========================================================================
74 {
75  //now do what ever initialization is needed
76 
77  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
78 
79  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
80  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
81 
82  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
83  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
84 }
85 
86 //========================================================================
87 EcalMatacqAnalyzer::~EcalMatacqAnalyzer() {
88  //========================================================================
89 
90  // do anything here that needs to be done at desctruction time
91  // (e.g. close files, deallocate resources etc.)
92 }
93 
94 //========================================================================
95 void EcalMatacqAnalyzer::beginJob() {
96  //========================================================================
97 
98  // Define temporary file name
99 
100  sampfile = resdir_;
101  sampfile += "/TmpTreeMatacqAnalyzer.root";
102 
103  sampFile = new TFile(sampfile.c_str(), "RECREATE");
104 
105  // declaration of the tree to fill
106 
107  tree = new TTree("MatacqTree", "MatacqTree");
108 
109  //List of branches
110 
111  tree->Branch("event", &event, "event/I");
112  tree->Branch("color", &color, "color/I");
113  tree->Branch("matacq", &matacq, "matacq[2560]/D");
114  tree->Branch("nsamples", &nsamples, "nsamples/I");
115  tree->Branch("maxsamp", &maxsamp, "maxsamp/I");
116  tree->Branch("tt", &tt, "tt/D");
117  tree->Branch("lightside", &lightside, "lightside/I");
118 
119  tree->SetBranchAddress("event", &event);
120  tree->SetBranchAddress("color", &color);
121  tree->SetBranchAddress("matacq", matacq);
122  tree->SetBranchAddress("nsamples", &nsamples);
123  tree->SetBranchAddress("maxsamp", &maxsamp);
124  tree->SetBranchAddress("tt", &tt);
125  tree->SetBranchAddress("lightside", &lightside);
126 
127  // Define output results files' names
128 
129  std::stringstream namefile;
130  namefile << resdir_ << "/MATACQ.root";
131  outfile = namefile.str();
132 
133  // Laser events counter
134  laserEvents = 0;
135  isThereMatacq = false;
136 }
137 
138 //========================================================================
139 void EcalMatacqAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
140  //========================================================================
141 
142  ++iEvent;
143 
144  if (_debug == 1)
145  std::cout << "-- debug test -- Entering Analyze -- event= " << iEvent << std::endl;
146 
147  // retrieving MATACQ :
148  edm::Handle<EcalMatacqDigiCollection> pmatacqDigi;
149  const EcalMatacqDigiCollection* matacqDigi = nullptr;
150  try {
151  e.getByLabel(digiProducer_, digiCollection_, pmatacqDigi);
152  matacqDigi = pmatacqDigi.product();
153  if (_debug == 1)
154  std::cout << "-- debug test -- Matacq Digis Found -- " << std::endl;
155 
156  } catch (std::exception& ex) {
157  std::cerr << "Error! can't get the product EcalMatacqDigi producer:" << digiProducer_.c_str()
158  << " collection:" << digiCollection_.c_str() << std::endl;
159  if (_debug == 1)
160  std::cout << "-- debug test -- No Matacq Digis Found -- " << std::endl;
161  return;
162  }
163 
164  // retrieving DCC header
165 
166  edm::Handle<EcalRawDataCollection> pDCCHeader;
167  const EcalRawDataCollection* DCCHeader = nullptr;
168  try {
169  e.getByLabel(eventHeaderProducer_, eventHeaderCollection_, pDCCHeader);
170  DCCHeader = pDCCHeader.product();
171  } catch (std::exception& ex) {
172  std::cerr << "Error! can't get the product EcalRawData producer:" << eventHeaderProducer_.c_str()
173  << " collection:" << eventHeaderCollection_.c_str() << std::endl;
174  return;
175  }
176 
177  // ====================================
178  // Decode Basic DCCHeader Information
179  // ====================================
180 
181  if (_debug == 1)
182  std::cout << "-- debug test -- Before header -- " << std::endl;
183 
184  for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
185  ++headerItr) {
186  EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
187  color = (int)settings.wavelength;
188  if (color < 0)
189  return;
190 
191  // Get run type and run number
192 
193  int fed = headerItr->fedId();
194 
195  if (fed != _fedid && _fedid != -999)
196  continue;
197 
198  runType = headerItr->getRunType();
199  runNum = headerItr->getRunNumber();
200  event = headerItr->getLV1();
201 
202  if (_debug == 1)
203  std::cout << "-- debug test -- runtype:" << runType << " event:" << event << " runNum:" << runNum << std::endl;
204 
205  dccID = headerItr->getDccInTCCCommand();
206  fedID = headerItr->fedId();
207  lightside = headerItr->getRtHalf();
208 
209  //assert (lightside<2 && lightside>=0);
210 
211  if (lightside != 1 && lightside != 0) {
212  std::cout << "Unexpected lightside: " << lightside << " for event " << iEvent << std::endl;
213  return;
214  }
215  if (_debug == 1) {
216  std::cout << "-- debug test -- Inside header before fed cut -- color=" << color << ", dcc=" << dccID
217  << ", fed=" << fedID << ", lightside=" << lightside << ", runType=" << runType << std::endl;
218  }
219 
220  // take event only if the fed corresponds to the DCC in TCC
221  if (600 + dccID != fedID)
222  continue;
223 
224  if (_debug == 1) {
225  std::cout << "-- debug test -- Inside header after fed cut -- color=" << color << ", dcc=" << dccID
226  << ", fed=" << fedID << ", lightside=" << lightside << ", runType=" << runType << std::endl;
227  }
228 
229  // Cut on runType
230 
231  if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
232  runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
233  return;
234 
235  std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
236  if (iter == colors.end()) {
237  colors.push_back(color);
238  std::cout << " new color found " << color << " " << colors.size() << std::endl;
239  }
240  }
241 
242  if (_debug == 1)
243  std::cout << "-- debug test -- Before digis -- Event:" << iEvent << std::endl;
244 
245  // Count laser events
246  laserEvents++;
247 
248  // ===========================
249  // Decode Matacq Information
250  // ===========================
251 
252  int iCh = 0;
253  double max = 0;
254 
255  for (EcalMatacqDigiCollection::const_iterator it = matacqDigi->begin(); it != matacqDigi->end();
256  ++it) { // Loop on matacq channel
257 
258  //
259  const EcalMatacqDigi& digis = *it;
260 
261  //if(digis.size()==0 || iCh>=N_channels) continue;
262  if (_debug == 1) {
263  std::cout << "-- debug test -- Inside digis -- digi size=" << digis.size() << std::endl;
264  }
265 
266  if (digis.size() == 0)
267  continue;
268  else
269  isThereMatacq = true;
270 
271  max = 0;
272  maxsamp = 0;
273  nsamples = digis.size();
274  tt = digis.tTrig();
275 
276  for (int i = 0; i < digis.size(); ++i) { // Loop on matacq samples
277  matacq[i] = -digis.adcCount(i);
278  if (matacq[i] > max) {
279  max = matacq[i];
280  maxsamp = i;
281  }
282  }
283  if (_debug == 1) {
284  std::cout << "-- debug test -- Inside digis -- nsamples=" << nsamples << ", max=" << max << std::endl;
285  }
286 
287  iCh++;
288  }
289 
290  if (_debug == 1)
291  std::cout << "-- debug test -- After digis -- Event: " << iEvent << std::endl;
292  tree->Fill();
293 
294 } // analyze
295 
296 //========================================================================
297 void EcalMatacqAnalyzer::endJob() {
298  // Don't do anything if there is no events
299  if (!isThereMatacq) {
300  std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
301  std::cout << "\t+=+ WARNING! NO MATACQ +=+" << std::endl;
302  std::cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
303 
304  // Remove temporary file
305  FILE* test;
306  test = fopen(sampfile.c_str(), "r");
307  if (test) {
308  std::stringstream del2;
309  del2 << "rm " << sampfile;
310  system(del2.str().c_str());
311  }
312  return;
313  }
314 
315  assert(colors.size() <= nColor);
316  unsigned int nCol = colors.size();
317 
318  for (unsigned int iCol = 0; iCol < nCol; iCol++) {
319  for (unsigned int iSide = 0; iSide < nSide; iSide++) {
320  MTQ[iCol][iSide] = new TMTQ();
321  }
322  }
323 
324  outFile = new TFile(outfile.c_str(), "RECREATE");
325 
326  TProfile* shapeMat = new TProfile("shapeLaser", "shapeLaser", _nsamplesshape, -0.5, double(_nsamplesshape) - 0.5);
327  TProfile* shapeMatTmp = new TProfile(
328  "shapeLaserTmp", "shapeLaserTmp", _timeaftmax + _timebefmax, -0.5, double(_timeaftmax + _timebefmax) - 0.5);
329 
330  std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
331  std::cout << "\t+=+ Analyzing MATACQ data +=+" << std::endl;
332 
333  //
334  // create output ntuple
335  //
336 
337  mtqShape = new TTree("MatacqShape", "MatacqShape");
338 
339  // list of branches
340  // keep Patrice's notations
341 
342  mtqShape->Branch("event", &event, "event/I");
343  mtqShape->Branch("color", &color, "color/I");
344  mtqShape->Branch("status", &status, "status/I");
345  mtqShape->Branch("peak", &peak, "peak/D");
346  mtqShape->Branch("sigma", &sigma, "sigma/D");
347  mtqShape->Branch("fit", &fit, "fit/D");
348  mtqShape->Branch("ampl", &ampl, "ampl/D");
349  mtqShape->Branch("trise", &trise, "trise/D");
350  mtqShape->Branch("fwhm", &fwhm, "fwhm/D");
351  mtqShape->Branch("fw20", &fw20, "fw20/D");
352  mtqShape->Branch("fw80", &fw80, "fw80/D");
353  mtqShape->Branch("ped", &ped, "ped/D");
354  mtqShape->Branch("pedsig", &pedsig, "pedsig/D");
355  mtqShape->Branch("ttrig", &ttrig, "ttrig/D");
356  mtqShape->Branch("sliding", &sliding, "sliding/D");
357 
358  mtqShape->SetBranchAddress("event", &event);
359  mtqShape->SetBranchAddress("color", &color);
360  mtqShape->SetBranchAddress("status", &status);
361  mtqShape->SetBranchAddress("peak", &peak);
362  mtqShape->SetBranchAddress("sigma", &sigma);
363  mtqShape->SetBranchAddress("fit", &fit);
364  mtqShape->SetBranchAddress("ampl", &ampl);
365  mtqShape->SetBranchAddress("fwhm", &fwhm);
366  mtqShape->SetBranchAddress("fw20", &fw20);
367  mtqShape->SetBranchAddress("fw80", &fw80);
368  mtqShape->SetBranchAddress("trise", &trise);
369  mtqShape->SetBranchAddress("ped", &ped);
370  mtqShape->SetBranchAddress("pedsig", &pedsig);
371  mtqShape->SetBranchAddress("ttrig", &ttrig);
372  mtqShape->SetBranchAddress("sliding", &sliding);
373 
374  unsigned int endsample;
375  unsigned int presample;
376 
377  // loop over the entries of the tree
378  TChain* fChain = (TChain*)tree;
379  Long64_t nentries = fChain->GetEntriesFast();
380  Long64_t nbytes = 0, nb = 0;
381 
382  for (Long64_t jentry = 0; jentry < nentries; jentry++) {
383  // load the event
384  Long64_t ientry = fChain->LoadTree(jentry);
385  if (ientry < 0)
386  break;
387  nb = fChain->GetEntry(jentry);
388  nbytes += nb;
389 
390  status = 0;
391  peak = -1;
392  sigma = 0;
393  fit = -1;
394  ampl = -1;
395  trise = -1;
396  ttrig = tt;
397  fwhm = 0;
398  fw20 = 0;
399  fw80 = 0;
400  ped = 0;
401  pedsig = 0;
402  sliding = 0;
403 
404  if (_debug == 1)
405  std::cout << "-- debug test -- inside loop 1 -- jentry:" << jentry << " over nentries=" << nentries << std::endl;
406 
407  // create the object for Matacq data analysis
408 
409  endsample = maxsamp + _nsamplesaftmax;
410  presample = int(_presample * nsamples / 100.);
411  TMatacq* mtq = new TMatacq(nsamples,
412  presample,
413  endsample,
414  _noiseCut,
415  _parabnbefmax,
416  _parabnaftmax,
417  _thres,
418  _lowlev,
419  _highlev,
420  _nevlasers,
421  _slide);
422 
423  if (_debug == 1)
424  std::cout << "-- debug test -- inside loop 2 -- " << std::endl;
425 
426  // analyse the Matacq data
427  if (mtq->rawPulseAnalysis(nsamples, &matacq[0]) == 0) {
428  status = 1;
429  ped = mtq->getBaseLine();
430  pedsig = mtq->getsigBaseLine();
431 
432  if (_debug == 1)
433  std::cout << "-- debug test -- inside loop 3 -- ped:" << ped << std::endl;
434  if (mtq->findPeak() == 0) {
435  peak = mtq->getTimpeak();
436  sigma = mtq->getsigTimpeak();
437  }
438  if (_debug == 1)
439  std::cout << "-- debug test -- inside loop 4 -- peak:" << peak << std::endl;
440  if (mtq->doFit() == 0) {
441  fit = mtq->getTimax();
442  ampl = mtq->getAmpl();
443  fwhm = mtq->getFwhm();
444  fw20 = mtq->getWidth20();
445  fw80 = mtq->getWidth80();
446  sliding = mtq->getSlide();
447  }
448  if (_debug == 1)
449  std::cout << "-- debug test -- inside loop 4 -- ampl:" << ampl << std::endl;
450  if (mtq->compute_trise() == 0) {
451  trise = mtq->getTrise();
452  }
453  if (_debug == 1)
454  std::cout << "-- debug test -- inside loop 5 -- trise:" << trise << std::endl;
455  }
456 
457  if (_debug == 1)
458  std::cout << "-- debug test -- inside loop 6 -- status:" << status << std::endl;
459 
460  if (status == 1 && mtq->findPeak() == 0) {
461  int firstS = int(peak - double(_timebefmax));
462  int lastS = int(peak + double(_timeaftmax));
463 
464  // Fill histo if there are enough samples
465  if (_debug == 1)
466  std::cout << "-- debug test -- inside loop 7 -- firstS:" << firstS << ", nsamples:" << nsamples << std::endl;
467 
468  if (firstS >= 0 && lastS <= nsamples) {
469  for (int i = firstS; i < lastS; i++) {
470  shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
471  }
472 
473  } else { // else extrapolate
474 
475  int firstSBis;
476 
477  if (firstS < 0) { // fill first bins with 0
478 
479  double thisped;
480  thisped = (matacq[0] + matacq[1] + matacq[2] + matacq[4] + matacq[5]) / 5.0;
481 
482  for (int i = firstS; i < 0; i++) {
483  shapeMatTmp->Fill(double(i) - firstS, thisped);
484  }
485  firstSBis = 0;
486 
487  } else {
488  firstSBis = firstS;
489  }
490 
491  if (lastS > nsamples) {
492  for (int i = firstSBis; i < int(nsamples); i++) {
493  shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
494  }
495 
496  //extrapolate with expo tail
497 
498  double expb = 0.998;
499  double matacqval = expb * matacq[nsamples - 1];
500 
501  for (int i = nsamples; i < lastS; i++) {
502  shapeMatTmp->Fill(double(i) - firstS, matacqval);
503  matacqval *= expb;
504  }
505 
506  } else {
507  for (int i = firstSBis; i < lastS; i++) {
508  shapeMatTmp->Fill(double(i) - firstS, matacq[i]);
509  }
510  }
511  }
512  }
513  if (_debug == 1)
514  std::cout << "-- debug test -- inside loop 8" << std::endl;
515 
516  // get back color
517 
518  int iCol = nCol;
519  for (unsigned int i = 0; i < nCol; i++) {
520  if (color == colors[i]) {
521  iCol = i;
522  i = nCol;
523  }
524  }
525  if (_debug == 1)
526  std::cout << "-- debug test -- inside loop 8bis color:" << color << " iCol:" << iCol << " nCol:" << nCol
527  << std::endl;
528 
529  // fill TMTQ
530 
531  if (status == 1 && mtq->findPeak() == 0 && mtq->doFit() == 0 && mtq->compute_trise() == 0)
532  MTQ[iCol][lightside]->addEntry(peak, sigma, fit, ampl, trise, fwhm, fw20, fw80, ped, pedsig, sliding);
533 
534  // fill the output tree
535 
536  if (_debug == 1)
537  std::cout << "-- debug test -- inside loop 9" << std::endl;
538  mtqShape->Fill();
539 
540  // clean up
541  delete mtq;
542  }
543 
544  if (_debug == 1)
545  std::cout << "-- debug test -- after loop " << std::endl;
546  sampFile->Close();
547 
548  double Peak[6], Sigma[6], Fit[6], Ampl[6], Trise[6], Fwhm[6], Fw20[6], Fw80[6], Ped[6], Pedsig[6], Sliding[6];
549  int Side;
550 
551  for (unsigned int iColor = 0; iColor < nCol; iColor++) {
552  std::stringstream nametree;
553  nametree << "MatacqCol" << colors[iColor];
554  meanTree[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
555  meanTree[iColor]->Branch("side", &Side, "Side/I");
556  meanTree[iColor]->Branch("peak", &Peak, "Peak[6]/D");
557  meanTree[iColor]->Branch("sigma", &Sigma, "Sigma[6]/D");
558  meanTree[iColor]->Branch("fit", &Fit, "Fit[6]/D");
559  meanTree[iColor]->Branch("ampl", &Ampl, "Ampl[6]/D");
560  meanTree[iColor]->Branch("trise", &Trise, "Trise[6]/D");
561  meanTree[iColor]->Branch("fwhm", &Fwhm, "Fwhm[6]/D");
562  meanTree[iColor]->Branch("fw20", &Fw20, "Fw20[6]/D");
563  meanTree[iColor]->Branch("fw80", &Fw80, "Fw80[6]/D");
564  meanTree[iColor]->Branch("ped", &Ped, "Ped[6]/D");
565  meanTree[iColor]->Branch("pedsig", &Pedsig, "Pedsig[6]/D");
566  meanTree[iColor]->Branch("sliding", &Sliding, "Sliding[6]/D");
567 
568  meanTree[iColor]->SetBranchAddress("side", &Side);
569  meanTree[iColor]->SetBranchAddress("peak", Peak);
570  meanTree[iColor]->SetBranchAddress("sigma", Sigma);
571  meanTree[iColor]->SetBranchAddress("fit", Fit);
572  meanTree[iColor]->SetBranchAddress("ampl", Ampl);
573  meanTree[iColor]->SetBranchAddress("fwhm", Fwhm);
574  meanTree[iColor]->SetBranchAddress("fw20", Fw20);
575  meanTree[iColor]->SetBranchAddress("fw80", Fw80);
576  meanTree[iColor]->SetBranchAddress("trise", Trise);
577  meanTree[iColor]->SetBranchAddress("ped", Ped);
578  meanTree[iColor]->SetBranchAddress("pedsig", Pedsig);
579  meanTree[iColor]->SetBranchAddress("sliding", Sliding);
580  }
581 
582  for (unsigned int iCol = 0; iCol < nCol; iCol++) {
583  for (unsigned int iSide = 0; iSide < nSides; iSide++) {
584  Side = iSide;
585  std::vector<double> val[TMTQ::nOutVar];
586 
587  for (int iVar = 0; iVar < TMTQ::nOutVar; iVar++) {
588  val[iVar] = MTQ[iCol][iSide]->get(iVar);
589 
590  for (unsigned int i = 0; i < val[iVar].size(); i++) {
591  switch (iVar) {
592  case TMTQ::iPeak:
593  Peak[i] = val[iVar].at(i);
594  break;
595  case TMTQ::iSigma:
596  Sigma[i] = val[iVar].at(i);
597  break;
598  case TMTQ::iFit:
599  Fit[i] = val[iVar].at(i);
600  break;
601  case TMTQ::iAmpl:
602  Ampl[i] = val[iVar].at(i);
603  break;
604  case TMTQ::iFwhm:
605  Fwhm[i] = val[iVar].at(i);
606  break;
607  case TMTQ::iFw20:
608  Fw20[i] = val[iVar].at(i);
609  break;
610  case TMTQ::iFw80:
611  Fw80[i] = val[iVar].at(i);
612  break;
613  case TMTQ::iTrise:
614  Trise[i] = val[iVar].at(i);
615  break;
616  case TMTQ::iPed:
617  Ped[i] = val[iVar].at(i);
618  break;
619  case TMTQ::iPedsig:
620  Pedsig[i] = val[iVar].at(i);
621  break;
622  case TMTQ::iSlide:
623  Sliding[i] = val[iVar].at(i);
624  break;
625  }
626  }
627  }
628  meanTree[iCol]->Fill();
629  if (_debug == 1)
630  std::cout << "-- debug test -- inside final loop " << std::endl;
631  }
632  }
633 
634  // Calculate maximum with pol 2
635 
636  int im = shapeMatTmp->GetMaximumBin();
637  double q1 = shapeMatTmp->GetBinContent(im - 1);
638  double q2 = shapeMatTmp->GetBinContent(im);
639  double q3 = shapeMatTmp->GetBinContent(im + 1);
640 
641  double a2 = (q3 + q1) / 2.0 - q2;
642  double a1 = q2 - q1 + a2 * (1 - 2 * im);
643  double a0 = q2 - a1 * im - a2 * im * im;
644 
645  double am = a0 - a1 * a1 / (4 * a2);
646 
647  // Compute pedestal
648 
649  double bl = 0;
650  for (unsigned int i = 1; i < _presampleshape + 1; i++) {
651  bl += shapeMatTmp->GetBinContent(i);
652  }
653  bl /= _presampleshape;
654 
655  // Compute and save laser shape
656 
657  if (_debug == 1)
658  std::cout << "-- debug test -- computing shape " << std::endl;
659 
660  int firstBin = 0;
661  double height = 0.0;
662 
663  for (unsigned int i = _timebefmax; i > _presampleshape; i--) {
664  height = shapeMatTmp->GetBinContent(i) - bl;
665 
666  if (height < (am - bl) * _cutwindow) {
667  firstBin = i;
668  i = _presampleshape;
669  }
670  }
671 
672  unsigned int lastBin = firstBin + _nsamplesshape;
673 
674  for (unsigned int i = firstBin; i < lastBin; i++) {
675  shapeMat->Fill(i - firstBin, shapeMatTmp->GetBinContent(i) - bl);
676  }
677 
678  mtqShape->Write();
679  for (unsigned int iColor = 0; iColor < nCol; iColor++) {
680  meanTree[iColor]->Write();
681  }
682  if (_debug == 1)
683  std::cout << "-- debug test -- writing " << std::endl;
684  shapeMat->Write();
685 
686  // close the output file
687  outFile->Close();
688 
689  // Remove temporary file
690  FILE* test;
691  test = fopen(sampfile.c_str(), "r");
692  if (test) {
693  std::stringstream del2;
694  del2 << "rm " << sampfile;
695  system(del2.str().c_str());
696  }
697 
698  std::cout << "\t+=+ .................... done +=+" << std::endl;
699  std::cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
700 }
701 
702 DEFINE_FWK_MODULE(EcalMatacqAnalyzer);
DeadROC_duringRun.runNum
runNum
Definition: DeadROC_duringRun.py:21
EcalMatacqAnalyzer::color
int color
Definition: EcalMatacqAnalyzer.h:76
EcalMatacqAnalyzer::runType
int runType
Definition: EcalMatacqAnalyzer.h:64
mps_fire.i
i
Definition: mps_fire.py:428
EcalMatacqAnalyzer::nSide
Definition: EcalMatacqAnalyzer.h:24
EcalMatacqDigi::adcCount
float adcCount(const int &i) const
Definition: EcalMatacqDigi.h:53
edm::SortedCollection::const_iterator
std::vector< T >::const_iterator const_iterator
Definition: SortedCollection.h:80
TMTQ::iPeak
Definition: TMTQ.h:10
MessageLogger.h
EcalMatacqAnalyzer::_parabnaftmax
unsigned int _parabnaftmax
Definition: EcalMatacqAnalyzer.h:37
EcalMatacqAnalyzer::beginJob
void beginJob() override
Definition: EcalMatacqAnalyzer.cc:95
colors
Definition: colors.py:1
EcalMatacqAnalyzer::_nsamplesshape
unsigned int _nsamplesshape
Definition: EcalMatacqAnalyzer.h:45
edm::Handle::product
T const * product() const
Definition: Handle.h:70
EcalMatacqDigi
Definition: EcalMatacqDigi.h:16
EcalMatacqAnalyzer::eventHeaderProducer_
std::string eventHeaderProducer_
Definition: EcalMatacqAnalyzer.h:55
groupFilesInBlocks.tt
int tt
Definition: groupFilesInBlocks.py:144
EcalMatacqAnalyzer::digiCollection_
std::string digiCollection_
Definition: EcalMatacqAnalyzer.h:52
TMatacq
Definition: TMatacq.h:6
EcalMatacqAnalyzer::ped
double ped
Definition: EcalMatacqAnalyzer.h:92
TMTQ::nOutVar
Definition: TMTQ.h:10
tree
Definition: tree.py:1
EcalMatacqAnalyzer::iEvent
int iEvent
Definition: EcalMatacqAnalyzer.h:27
gather_cfg.cout
cout
Definition: gather_cfg.py:144
TMatacq::compute_trise
int compute_trise()
Definition: TMatacq.cc:353
EcalDetIdCollections.h
EcalMatacqAnalyzer::tt
double tt
Definition: EcalMatacqAnalyzer.h:80
TMTQ::get
std::vector< double > get(int)
Definition: TMTQ.cc:69
if
if(0==first)
Definition: CAHitNtupletGeneratorKernelsImpl.h:48
cms::cuda::assert
assert(be >=bs)
edm::SortedCollection
Definition: SortedCollection.h:49
EcalMatacqAnalyzer
Definition: EcalMatacqAnalyzer.h:14
EcalMatacqAnalyzer::nSides
unsigned int nSides
Definition: EcalMatacqAnalyzer.h:62
EcalMatacqAnalyzer::_cutwindow
double _cutwindow
Definition: EcalMatacqAnalyzer.h:44
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
EcalMatacqAnalyzer::_timeaftmax
unsigned int _timeaftmax
Definition: EcalMatacqAnalyzer.h:43
TMatacq::doFit
int doFit()
Definition: TMatacq.cc:219
EcalMatacqAnalyzer::outfile
std::string outfile
Definition: EcalMatacqAnalyzer.h:57
TMatacq::getsigBaseLine
double getsigBaseLine()
Definition: TMatacq.h:63
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
edm::Handle
Definition: AssociativeIterator.h:50
EcalMatacqAnalyzer::fwhm
double fwhm
Definition: EcalMatacqAnalyzer.h:92
EcalMatacqAnalyzer::nsamples
int nsamples
Definition: EcalMatacqAnalyzer.h:79
EcalMatacqAnalyzer::_lowlev
unsigned int _lowlev
Definition: EcalMatacqAnalyzer.h:39
EcalMatacqAnalyzer::sliding
double sliding
Definition: EcalMatacqAnalyzer.h:92
TMatacq::getTimpeak
double getTimpeak()
Definition: TMatacq.h:65
EcalMatacqAnalyzer::resdir_
std::string resdir_
Definition: EcalMatacqAnalyzer.h:51
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
EcalMatacqAnalyzer::outFile
TFile * outFile
Definition: EcalMatacqAnalyzer.h:90
TMatacq::rawPulseAnalysis
int rawPulseAnalysis(Int_t, Double_t *)
Definition: TMatacq.cc:57
EcalDCCHeaderBlock::LASER_GAP
Definition: EcalDCCHeaderBlock.h:42
TMTQ
Definition: TMTQ.h:8
TMatacq::getsigTimpeak
double getsigTimpeak()
Definition: TMatacq.h:66
EcalRawDataCollections.h
MakerMacros.h
EcalMatacqAnalyzer::_presampleshape
unsigned int _presampleshape
Definition: EcalMatacqAnalyzer.h:46
test
Definition: SmallWORMDict.h:13
EcalDCCHeaderBlock::LASER_STD
Definition: EcalDCCHeaderBlock.h:26
ctpps_dqm_sourceclient-live_cfg.test
test
Definition: ctpps_dqm_sourceclient-live_cfg.py:7
EcalMatacqAnalyzer::dccID
int dccID
Definition: EcalMatacqAnalyzer.h:66
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
TMTQ::iPed
Definition: TMTQ.h:10
EcalMatacqAnalyzer::sigma
double sigma
Definition: EcalMatacqAnalyzer.h:92
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:262
TMTQ::iTrise
Definition: TMTQ.h:10
Fit
Definition: Fit.h:32
EcalMatacqAnalyzer::ampl
double ampl
Definition: EcalMatacqAnalyzer.h:92
EcalMatacqAnalyzer::status
int status
Definition: EcalMatacqAnalyzer.h:91
EcalMatacqAnalyzer::_noiseCut
unsigned int _noiseCut
Definition: EcalMatacqAnalyzer.h:35
EcalDigiCollections.h
EcalMatacqAnalyzer::_parabnbefmax
unsigned int _parabnbefmax
Definition: EcalMatacqAnalyzer.h:36
EcalMatacqAnalyzer::_nsamplesaftmax
unsigned int _nsamplesaftmax
Definition: EcalMatacqAnalyzer.h:34
EcalMatacqAnalyzer::ttrig
double ttrig
Definition: EcalMatacqAnalyzer.h:92
EcalMatacqAnalyzer::maxsamp
int maxsamp
Definition: EcalMatacqAnalyzer.h:78
EcalMatacqAnalyzer::nColor
Definition: EcalMatacqAnalyzer.h:23
TMTQ::iAmpl
Definition: TMTQ.h:10
NSIDES
#define NSIDES
Definition: EcalLaserAnalyzer.h:44
TMatacq::getTimax
double getTimax()
Definition: TMatacq.h:69
TMatacq::getAmpl
double getAmpl()
Definition: TMatacq.h:68
EcalMatacqAnalyzer.h
q2
double q2[4]
Definition: TauolaWrapper.h:88
EcalMatacqDigi::size
int size() const
Definition: EcalMatacqDigi.h:75
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
EcalDCCHeaderBlock::LASER_POWER_SCAN
Definition: EcalDCCHeaderBlock.h:27
EcalDCCHeaderBlock::EcalDCCEventSettings
Definition: EcalDCCHeaderBlock.h:64
TMatacq::getWidth80
double getWidth80()
Definition: TMatacq.h:74
EcalMatacqAnalyzer::_slide
unsigned int _slide
Definition: EcalMatacqAnalyzer.h:47
cppFunctionSkipper.exception
exception
Definition: cppFunctionSkipper.py:10
TMatacq::getWidth20
double getWidth20()
Definition: TMatacq.h:73
edm::ParameterSet
Definition: ParameterSet.h:47
q1
double q1[4]
Definition: TauolaWrapper.h:87
EcalMatacqAnalyzer::fw80
double fw80
Definition: EcalMatacqAnalyzer.h:92
Event.h
EcalMatacqAnalyzer::laserEvents
int laserEvents
Definition: EcalMatacqAnalyzer.h:70
TMTQ::iFw80
Definition: TMTQ.h:10
EcalMatacqAnalyzer::isThereMatacq
bool isThereMatacq
Definition: EcalMatacqAnalyzer.h:71
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
TMTQ::iFw20
Definition: TMTQ.h:10
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
EcalMatacqAnalyzer::EcalMatacqAnalyzer
EcalMatacqAnalyzer(const edm::ParameterSet &iConfig)
Definition: EcalMatacqAnalyzer.cc:38
TMTQ::iPedsig
Definition: TMTQ.h:10
createfilelist.int
int
Definition: createfilelist.py:10
EcalMatacqAnalyzer::trise
double trise
Definition: EcalMatacqAnalyzer.h:92
iEvent
int iEvent
Definition: GenABIO.cc:224
EcalMatacqAnalyzer::_debug
int _debug
Definition: EcalMatacqAnalyzer.h:49
EcalMatacqAnalyzer::MTQ
TMTQ * MTQ[nColor][nSide]
Definition: EcalMatacqAnalyzer.h:85
edm::EventSetup
Definition: EventSetup.h:58
TMatacq::getSlide
double getSlide()
Definition: TMatacq.h:75
EcalMatacqAnalyzer::analyze
void analyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: EcalMatacqAnalyzer.cc:139
EcalMatacqAnalyzer::_nevlasers
unsigned int _nevlasers
Definition: EcalMatacqAnalyzer.h:41
EcalMatacqAnalyzer::lightside
int lightside
Definition: EcalMatacqAnalyzer.h:63
EcalMatacqAnalyzer::_fedid
int _fedid
Definition: EcalMatacqAnalyzer.h:48
EcalMatacqDigi::tTrig
float tTrig() const
Definition: EcalMatacqDigi.h:105
SurfaceOrientation::Side
Side
Definition: Surface.h:18
EcalMatacqAnalyzer::~EcalMatacqAnalyzer
~EcalMatacqAnalyzer() override
Definition: EcalMatacqAnalyzer.cc:87
TMatacq::findPeak
int findPeak()
Definition: TMatacq.cc:136
EcalMatacqAnalyzer::sampfile
std::string sampfile
Definition: EcalMatacqAnalyzer.h:58
heppy_batch.val
val
Definition: heppy_batch.py:351
std
Definition: JetResolutionObject.h:76
EcalMatacqAnalyzer::_highlev
unsigned int _highlev
Definition: EcalMatacqAnalyzer.h:40
EcalMatacqAnalyzer::fw20
double fw20
Definition: EcalMatacqAnalyzer.h:92
EcalMatacqAnalyzer::matacq
double matacq[2560]
Definition: EcalMatacqAnalyzer.h:77
TMatacq::getFwhm
double getFwhm()
Definition: TMatacq.h:72
EcalMatacqAnalyzer::mtqShape
TTree * mtqShape
Definition: EcalMatacqAnalyzer.h:93
EcalMatacqAnalyzer::runNum
int runNum
Definition: EcalMatacqAnalyzer.h:65
EcalDCCHeaderBlock::EcalDCCEventSettings::wavelength
short wavelength
Definition: EcalDCCHeaderBlock.h:67
EcalMatacqAnalyzer::endJob
void endJob() override
Definition: EcalMatacqAnalyzer.cc:297
EcalMatacqAnalyzer::pedsig
double pedsig
Definition: EcalMatacqAnalyzer.h:92
TMatacq::getBaseLine
double getBaseLine()
Definition: TMatacq.h:62
TMom.h
EcalDCCHeaderBlock::LASER_DELAY_SCAN
Definition: EcalDCCHeaderBlock.h:28
EcalMatacqAnalyzer::_presample
double _presample
Definition: EcalMatacqAnalyzer.h:33
EcalMatacqAnalyzer::_thres
unsigned int _thres
Definition: EcalMatacqAnalyzer.h:38
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
Exception.h
TMatacq.h
TMatacq::getTrise
double getTrise()
Definition: TMatacq.h:71
EcalMatacqAnalyzer::fedID
int fedID
Definition: EcalMatacqAnalyzer.h:67
ParameterSet.h
TMTQ::iFit
Definition: TMTQ.h:10
EcalCondDBWriter_cfi.runType
runType
Definition: EcalCondDBWriter_cfi.py:64
c
auto & c
Definition: CAHitNtupletGeneratorKernelsImpl.h:46
TMTQ::iSlide
Definition: TMTQ.h:10
EcalMatacqAnalyzer::meanTree
TTree * meanTree[nColor]
Definition: EcalMatacqAnalyzer.h:86
event
Definition: event.py:1
edm::Event
Definition: Event.h:73
EcalMatacqAnalyzer::sampFile
TFile * sampFile
Definition: EcalMatacqAnalyzer.h:82
EcalMatacqAnalyzer::digiProducer_
std::string digiProducer_
Definition: EcalMatacqAnalyzer.h:53
EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.cerr
cerr
Definition: EcnaPython_AdcPeg12_S1_10_R170298_1_0_150_Dee0.py:8
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
EcalMatacqAnalyzer::eventHeaderCollection_
std::string eventHeaderCollection_
Definition: EcalMatacqAnalyzer.h:54
fit
Definition: CombinedChiSquaredLikelihood.h:6
TMTQ::iSigma
Definition: TMTQ.h:10
EcalMatacqAnalyzer::_timebefmax
unsigned int _timebefmax
Definition: EcalMatacqAnalyzer.h:42
EcalMatacqAnalyzer::peak
double peak
Definition: EcalMatacqAnalyzer.h:92
TMTQ.h
TMTQ::iFwhm
Definition: TMTQ.h:10
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16