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