CMS 3D CMS Logo

EcalLaserAnalyzer2.cc
Go to the documentation of this file.
1 /*
2  * \class EcalLaserAnalyzer2
3  * primary author: Julie Malcles - CEA/Saclay
4  * author: Gautier Hamel De Monchenault - CEA/Saclay
5  */
6 
7 #include <TAxis.h>
8 #include <TH1.h>
9 #include <TProfile.h>
10 #include <TTree.h>
11 #include <TChain.h>
12 #include <TFile.h>
13 #include <TMath.h>
14 
16 
17 #include <sstream>
18 #include <iostream>
19 #include <iomanip>
20 #include <fstream>
21 
24 
27 
30 
38 
39 #include <vector>
40 
52 
53 
54 using namespace std;
55 
56 
57 //========================================================================
59  //========================================================================
60  :
61 iEvent(0),
62 
63 
64 // framework parameters with default values
65 _nsamples( iConfig.getUntrackedParameter< unsigned int >( "nSamples", 10 ) ),
66 _presample( iConfig.getUntrackedParameter< unsigned int >( "nPresamples", 2 ) ),
67 _firstsample( iConfig.getUntrackedParameter< unsigned int >( "firstSample", 1 ) ),
68 _lastsample( iConfig.getUntrackedParameter< unsigned int >( "lastSample", 2 ) ),
69 _nsamplesPN( iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN", 50 ) ),
70 _presamplePN( iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN", 6 ) ),
71 _firstsamplePN( iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN", 7 ) ),
72 _lastsamplePN( iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN", 8 ) ),
73 _timingcutlow( iConfig.getUntrackedParameter< unsigned int >( "timingCutLow", 2 ) ),
74 _timingcuthigh( iConfig.getUntrackedParameter< unsigned int >( "timingCutHigh", 9 ) ),
75 _timingquallow( iConfig.getUntrackedParameter< unsigned int >( "timingQualLow", 3 ) ),
76 _timingqualhigh( iConfig.getUntrackedParameter< unsigned int >( "timingQualHigh", 8 ) ),
77 _ratiomincutlow( iConfig.getUntrackedParameter< double >( "ratioMinCutLow", 0.4 ) ),
78 _ratiomincuthigh( iConfig.getUntrackedParameter< double >( "ratioMinCutHigh", 0.95 ) ),
79 _ratiomaxcutlow( iConfig.getUntrackedParameter< double >( "ratioMaxCutLow", 0.8 ) ),
80 _presamplecut( iConfig.getUntrackedParameter< double >( "presampleCut", 5.0 ) ),
81 _niter( iConfig.getUntrackedParameter< unsigned int >( "nIter", 5 ) ),
82 _noise( iConfig.getUntrackedParameter< double >( "noise", 2.0 ) ),
83 _ecalPart( iConfig.getUntrackedParameter< std::string >( "ecalPart", "EB" ) ),
84 _saveshapes( iConfig.getUntrackedParameter< bool >( "saveShapes", true ) ),
85 _docorpn( iConfig.getUntrackedParameter< bool >( "doCorPN", false ) ),
86 _fedid( iConfig.getUntrackedParameter< int >( "fedID", -999 ) ),
87 _saveallevents( iConfig.getUntrackedParameter< bool >( "saveAllEvents", false ) ),
88 _qualpercent( iConfig.getUntrackedParameter< double >( "qualPercent", 0.2 ) ),
89 _debug( iConfig.getUntrackedParameter< int >( "debug", 0 ) ),
90 nCrys( NCRYSEB),
91 nPNPerMod( NPNPERMOD),
92 nMod( NMODEB),
93 nSides( NSIDES),
94 nSamplesShapes( NSAMPSHAPES),
95 IsMatacqOK(false),
96 runType(-1), runNum(0),towerID(-1), channelID(-1), fedID(-1), dccID(-1), side(2), lightside(2),
97 iZ(1), phi(-1), eta(-1), event(0), color(0),pn0(0), pn1(0), apdAmpl(0),apdTime(0),pnAmpl(0),
98 pnID(-1), moduleID(-1), flag(0), channelIteratorEE(0), ShapeCor(0)
99 
100 
101  //========================================================================
102 
103 {
104 
105 
106  // Initialization from cfg file
107 
108  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
109  elecfile_ = iConfig.getUntrackedParameter<std::string>("elecFile");
110  pncorfile_ = iConfig.getUntrackedParameter<std::string>("pnCorFile");
111 
112  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
113  digiPNCollection_ = iConfig.getParameter<std::string>("digiPNCollection");
114  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
115 
116  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
117  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
118 
119 
120  // Geometrical constants initialization
121  if (_ecalPart == "EB") {
122  nCrys = NCRYSEB;
123  } else {
124  nCrys = NCRYSEE;
125  }
126  iZ = 1;
127  if(_fedid <= 609 )
128  iZ = -1;
130  nMod = modules.size();
131  nRefChan = NREFCHAN;
132 
133  for(unsigned int j=0;j<nCrys;j++){
134  iEta[j]=-1;
135  iPhi[j]=-1;
136  iModule[j]=10;
137  iTowerID[j]=-1;
138  iChannelID[j]=-1;
139  idccID[j]=-1;
140  iside[j]=-1;
141  wasTimingOK[j]=true;
142  wasGainOK[j]=true;
143  }
144 
145  for(unsigned int j=0;j<nMod;j++){
146  int ii= modules[j];
147  firstChanMod[ii-1]=0;
148  isFirstChanModFilled[ii-1]=0;
149  }
150 
151  // Quality check flags
152 
153  isGainOK=true;
154  isTimingOK=true;
155 
156  // PN linearity corrector
157 
158  pnCorrector = new TPNCor(pncorfile_);
159 
160 
161  // Objects dealing with pulses
162 
167 
168  // Object dealing with MEM numbering
169 
170  Mem = new TMem(_fedid);
171 
172  // Objects needed for npresample calculation
173 
174  Delta01=new TMom();
175  Delta12=new TMom();
176 
177 }
178 
179 //========================================================================
181  //========================================================================
182 
183 
184  // do anything here that needs to be done at desctruction time
185  // (e.g. close files, deallocate resources etc.)
186 
187 }
188 
189 
190 
191 //========================================================================
193  //========================================================================
194 
195 
196  // Create temporary files and trees to save adc samples
197  //======================================================
198 
200  ADCfile+="/APDSamplesLaser.root";
201 
203  APDfile+="/APDPNLaserAllEvents.root";
204 
205  ADCFile = new TFile(ADCfile.c_str(),"RECREATE");
206 
207  for (unsigned int i=0;i<nCrys;i++){
208 
209  stringstream name;
210  name << "ADCTree" <<i+1;
211  ADCtrees[i]= new TTree(name.str().c_str(),name.str().c_str());
212 
213  ADCtrees[i]->Branch( "iphi", &phi, "phi/I" );
214  ADCtrees[i]->Branch( "ieta", &eta, "eta/I" );
215  ADCtrees[i]->Branch( "towerID", &towerID, "towerID/I" );
216  ADCtrees[i]->Branch( "channelID", &channelID, "channelID/I" );
217  ADCtrees[i]->Branch( "dccID", &dccID, "dccID/I" );
218  ADCtrees[i]->Branch( "side", &side, "side/I" );
219  ADCtrees[i]->Branch( "event", &event, "event/I" );
220  ADCtrees[i]->Branch( "color", &color, "color/I" );
221  ADCtrees[i]->Branch( "adc", &adc , "adc[10]/D" );
222  ADCtrees[i]->Branch( "pn0", &pn0 , "pn0/D" );
223  ADCtrees[i]->Branch( "pn1", &pn1 , "pn1/D" );
224 
225 
226  ADCtrees[i]->SetBranchAddress( "ieta", &eta );
227  ADCtrees[i]->SetBranchAddress( "iphi", &phi );
228  ADCtrees[i]->SetBranchAddress( "towerID", &towerID );
229  ADCtrees[i]->SetBranchAddress( "channelID", &channelID );
230  ADCtrees[i]->SetBranchAddress( "dccID", &dccID );
231  ADCtrees[i]->SetBranchAddress( "side", &side );
232  ADCtrees[i]->SetBranchAddress( "event", &event );
233  ADCtrees[i]->SetBranchAddress( "color", &color );
234  ADCtrees[i]->SetBranchAddress( "adc", adc );
235  ADCtrees[i]->SetBranchAddress( "pn0", &pn0 );
236  ADCtrees[i]->SetBranchAddress( "pn1", &pn1 );
237 
238  }
239 
240  // Define output results filenames
241  //==================================
242  stringstream namefile1;
243  namefile1 << resdir_ <<"/SHAPE_LASER.root";
244  shapefile=namefile1.str();
245 
246  stringstream namefile2;
247  namefile2 << resdir_ <<"/APDPN_LASER.root";
248  resfile=namefile2.str();
249 
250  stringstream namefile3;
251  namefile3 << resdir_ <<"/MATACQ.root";
252 
253  matfile=namefile3.str();
254 
255  // Get Pulse Shapes
256  //==================
257 
259  if(!IsMatacqOK){
260  cout <<" ERROR! No matacq shape available: analysis aborted !"<< endl;
261  return;
262  }
263 
264  // Laser events counter
265 
266  laserEvents=0;
267 
268 
269 }
270 
271 
272 //========================================================================
274  //========================================================================
275 
276  ++iEvent;
277 
278 
279  // retrieving DCC header
281  const EcalRawDataCollection* DCCHeader=nullptr;
282  try {
284  DCCHeader=pDCCHeader.product();
285  }catch ( std::exception& ex ) {
286  std::cerr << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str() << std::endl;
287  }
288 
289  //retrieving crystal data from Event
291  const EBDigiCollection* EBDigi=nullptr;
293  const EEDigiCollection* EEDigi=nullptr;
294 
295  if (_ecalPart == "EB") {
296  try {
298  EBDigi=pEBDigi.product();
299  }catch ( std::exception& ex ) {
300  std::cerr << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str() << std::endl;
301  }
302  } else if (_ecalPart == "EE") {
303  try {
305  EEDigi=pEEDigi.product();
306  }catch ( std::exception& ex ) {
307  std::cerr << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str() << std::endl;
308  }
309  } else {
310  cout <<" Wrong ecalPart in cfg file " << endl;
311  return;
312  }
313 
314 
315  // retrieving crystal PN diodes from Event
316 
318  const EcalPnDiodeDigiCollection* PNDigi=nullptr;
319  try {
321  PNDigi=pPNDigi.product();
322  }catch ( std::exception& ex ) {
323  std::cerr << "Error! can't get the product " << digiPNCollection_.c_str() << std::endl;
324  }
325 
326  // retrieving electronics mapping
328  const EcalElectronicsMapping* TheMapping=nullptr;
329  try{
330  c.get< EcalMappingRcd >().get(ecalmapping);
331  TheMapping = ecalmapping.product();
332  }catch ( std::exception& ex ) {
333  std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
334  }
335 
336  // ====================================
337  // Decode Basic DCCHeader Information
338  // ====================================
339 
340 
341  for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();headerItr != DCCHeader->end();
342  ++headerItr ) {
343 
344  // Get run type and run number
345 
346  int fed = headerItr->fedId();
347  if(fed!=_fedid && _fedid!=-999) continue;
348 
349  runType=headerItr->getRunType();
350  runNum=headerItr->getRunNumber();
351  event=headerItr->getLV1();
352 
353  dccID=headerItr->getDccInTCCCommand();
354  fedID=headerItr->fedId();
355  lightside=headerItr->getRtHalf();
356 
357  // Check fed corresponds to the DCC in TCC
358 
359  if( 600+dccID != fedID ) continue;
360 
361  // Cut on runType
362 
367 
368  // Retrieve laser color and event number
369 
370  EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
371  color = settings.wavelength;
372  if( color<0 ) return;
373 
374  vector<int>::iterator iter = find( colors.begin(), colors.end(), color );
375  if( iter==colors.end() ){
376  colors.push_back( color );
377  cout <<" new color found "<< color<<" "<< colors.size()<< endl;
378  }
379  }
380 
381  // Check Matacq shape exists
382 
383  if(!IsMatacqOK) return;
384 
385 
386  // Cut on fedID
387 
388  if(fedID!=_fedid && _fedid!=-999) return;
389 
390  // Count laser events
391 
392  laserEvents++;
393 
394 
395  // ======================
396  // Decode PN Information
397  // ======================
398 
399  TPNFit * pnfit = new TPNFit();
401 
402  double chi2pn=0;
403  unsigned int samplemax=0;
404  int pnGain=0;
405 
406  map <int, vector<double> > allPNAmpl;
407  map <int, vector<double> > allPNGain;
408 
409 
410  // Loop on PNs digis
411 
412  for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr ) {
413 
414  EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
415 
416  if (_debug==1) cout <<"-- debug test -- Inside PNDigi - pnID=" <<
417  pnDetId.iPnId()<<", dccID="<< pnDetId.iDCCId()<< endl;
418 
419  // Skip MEM DCC without relevant data
420 
421  bool isMemRelevant=Mem->isMemRelevant(pnDetId.iDCCId());
422  if(!isMemRelevant) continue;
423 
424  // Loop on PN samples
425 
426  for ( int samId=0; samId < (*pnItr).size() ; samId++ ) {
427  pn[samId]=(*pnItr).sample(samId).adc();
428  pnG[samId]=(*pnItr).sample(samId).gainId();
429  if (samId==0) pnGain=pnG[samId];
430  if (samId>0) pnGain=int(TMath::Max(pnG[samId],pnGain));
431  }
432 
433  if(pnGain!=1) cout << "PN gain different from 1"<< endl;
434 
435  // Calculate amplitude from pulse
436 
437  PNPulse->setPulse(pn);
439  samplemax=PNPulse->getMaxSample();
440  chi2pn = pnfit -> doFit(samplemax,&pnNoPed[0]);
441  if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
442  else pnAmpl= pnfit -> getAmpl();
443 
444  // Apply linearity correction
445 
446  double corr=1.0;
447  if( _docorpn ) corr=pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
448  pnAmpl*=corr;
449 
450  // Fill PN ampl vector
451 
452  allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
453 
454  if (_debug==1) cout <<"-- debug -- Inside PNDigi - PNampl=" <<
455  pnAmpl<<", PNgain="<< pnGain<<endl;
456 
457  }
458 
459 
460 
461  // ===========================
462  // Decode EBDigis Information
463  // ===========================
464 
465  int adcGain=0;
466 
467  if (EBDigi){
468 
469  // Loop on crystals
470  //===================
471 
472  for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr ) { // Loop on crystals
473 
474 
475  // Retrieve geometry
476  //===================
477 
478  EBDetId id_crystal(digiItr->id()) ;
479  EBDataFrame df( *digiItr );
480  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
481 
482  int etaG = id_crystal.ieta() ; // global
483  int phiG = id_crystal.iphi() ; // global
484 
485  std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
486 
487  int etaL=LocalCoord.first ; // local
488  int phiL=LocalCoord.second ;// local
489 
490  int strip=elecid_crystal.stripId();
491  int xtal=elecid_crystal.xtalId();
492 
493  int module= MEEBGeom::lmmod(etaG, phiG);
494  int tower=elecid_crystal.towerId();
495 
496  int apdRefTT=MEEBGeom::apdRefTower(module);
497 
498  std::pair<int,int> pnpair=MEEBGeom::pn(module);
499  unsigned int MyPn0=pnpair.first;
500  unsigned int MyPn1=pnpair.second;
501 
502  int lmr=MEEBGeom::lmr( etaG,phiG );
503  unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
504  assert( channel < nCrys );
505 
506  setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);
507 
508  if (_debug==1) cout << "-- debug -- Inside EBDigi - towerID:"<< towerID<<
509  " channelID:" <<channelID<<" module:"<< module<<
510  " modules:"<<modules.size()<< endl;
511 
512  // APD Pulse
513  //===========
514 
515  // Loop on adc samples
516 
517  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {
518 
519  EcalMGPASample samp_crystal(df.sample(i));
520  adc[i]=samp_crystal.adc() ;
521  adcG[i]=samp_crystal.gainId();
522  adc[i]*=adcG[i];
523  if (i==0) adcGain=adcG[i];
524  if (i>0) adcGain=TMath::Max(adcG[i],adcGain);
525  }
526 
528 
529 
530  // Quality checks
531  //================
532 
533  if(adcGain!=1) nEvtBadGain[channel]++;
534  if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
535  nEvtTot[channel]++;
536 
537  // Associate PN ampl
538  //===================
539 
540  int mem0=Mem->Mem(lmr,0);
541  int mem1=Mem->Mem(lmr,1);
542 
543  if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
544  else pn0=0;
545  if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
546  else pn1=0;
547 
548  // Fill if Pulse is fine
549  //=======================
550 
551  if( APDPulse->isPulseOK() && lightside==side){
552 
553  ADCtrees[channel]->Fill();
554 
557 
558  }
559  }
560  } else if (EEDigi) {
561 
562  // Loop on crystals
563  //===================
564 
565  for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin();
566  digiItr != EEDigi->end(); ++digiItr ) {
567 
568  // Retrieve geometry
569  //===================
570 
571  EEDetId id_crystal(digiItr->id()) ;
572  EEDataFrame df( *digiItr );
573  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
574 
575  int etaG = id_crystal.iy() ;
576  int phiG = id_crystal.ix() ;
577 
578  int iX = (phiG-1)/5+1;
579  int iY = (etaG-1)/5+1;
580 
581  int tower=elecid_crystal.towerId();
582  int ch=elecid_crystal.channelId()-1;
583 
584  int module=MEEEGeom::lmmod( iX, iY );
585  if( module>=18 && side==1 ) module+=2;
586  int lmr=MEEEGeom::lmr( iX, iY ,iZ);
587  int dee=MEEEGeom::dee(lmr);
588  int apdRefTT=MEEEGeom::apdRefTower(lmr, module);
589 
590  std::pair<int,int> pnpair=MEEEGeom::pn( dee, module ) ;
591  unsigned int MyPn0=pnpair.first;
592  unsigned int MyPn1=pnpair.second;
593 
594  int hashedIndex=100000*eta+phi;
595  if( channelMapEE.count(hashedIndex) == 0 ){
598  }
599  unsigned int channel=channelMapEE[hashedIndex];
600  assert ( channel < nCrys );
601 
602  setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);
603 
604 
605  if (_debug==1) cout << "-- debug -- Inside EEDigi - towerID:"<< towerID<<
606  " channelID:" <<channelID<<" module:"<< module<<
607  " modules:"<<modules.size()<< endl;
608 
609  // APD Pulse
610  //===========
611 
612  if( (*digiItr).size()>10) cout <<"SAMPLES SIZE > 10!" << (*digiItr).size()<< endl;
613 
614  // Loop on adc samples
615 
616  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {
617 
618  EcalMGPASample samp_crystal(df.sample(i));
619  adc[i]=samp_crystal.adc() ;
620  adcG[i]=samp_crystal.gainId();
621  adc[i]*=adcG[i];
622 
623  if (i==0) adcGain=adcG[i];
624  if (i>0) adcGain=TMath::Max(adcG[i],adcGain);
625  }
626 
628 
629  // Quality checks
630  //================
631 
632  if(adcGain!=1) nEvtBadGain[channel]++;
633  if(!APDPulse->isTimingQualOK()) nEvtBadTiming[channel]++;
634  nEvtTot[channel]++;
635 
636 
637  // Associate PN ampl
638  //===================
639 
640  int mem0=Mem->Mem(lmr,0);
641  int mem1=Mem->Mem(lmr,1);
642 
643  if(allPNAmpl[mem0].size()>MyPn0) pn0=allPNAmpl[mem0][MyPn0];
644  else pn0=0;
645  if(allPNAmpl[mem1].size()>MyPn1) pn1=allPNAmpl[mem1][MyPn1];
646  else pn1=0;
647 
648  // Fill if Pulse is fine
649  //=======================
650 
651  if( APDPulse->isPulseOK() && lightside==side){
652  ADCtrees[channel]->Fill();
653 
656 
657  }
658  }
659  }
660 }
661 // analyze
662 
663 
664 //========================================================================
666 //========================================================================
667 
668  if(!IsMatacqOK){
669 
670  cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
671  cout << "\t+=+ WARNING! NO MATACQ +=+" << endl;
672  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
673  return;
674 
675  }
676 
677  // Adjust channel numbers for EE
678  //===============================
679 
680  if( _ecalPart == "EE" ) {
681  nCrys=channelMapEE.size();
682  }
683 
684  // Set presamples number
685  //======================
686 
687  double delta01=Delta01->getMean();
688  double delta12=Delta12->getMean();
689  if(delta12>_presamplecut) {
690  _presample=2;
691  if(delta01>_presamplecut) _presample=1;
692  }
693 
695 
696  // Don't do anything if there is no events
697  //=========================================
698 
699  if( laserEvents == 0 ){
700 
701  ADCFile->Close();
702  stringstream del;
703  del << "rm " <<ADCfile;
704  system(del.str().c_str());
705  cout << " No Laser Events "<< endl;
706  return;
707  }
708 
709  // Set quality flags for gains and timing
710  //=========================================
711 
712  double BadGainEvtPercentage=0.0;
713  double BadTimingEvtPercentage=0.0;
714 
715  int nChanBadGain=0;
716  int nChanBadTiming=0;
717 
718  for (unsigned int i=0;i<nCrys;i++){
719  if(nEvtTot[i]!=0){
720  BadGainEvtPercentage=double(nEvtBadGain[i])/double(nEvtTot[i]);
721  BadTimingEvtPercentage=double(nEvtBadTiming[i])/double(nEvtTot[i]);
722  }
723  if(BadGainEvtPercentage>_qualpercent) {
724  wasGainOK[i]=false;
725  nChanBadGain++;
726  }
727  if(BadTimingEvtPercentage>_qualpercent){
728  wasTimingOK[i]=false;
729  nChanBadTiming++;
730  }
731  }
732 
733  double BadGainChanPercentage=double(nChanBadGain)/double(nCrys);
734  double BadTimingChanPercentage=double(nChanBadTiming)/double(nCrys);
735 
736  if(BadGainChanPercentage>_qualpercent) isGainOK = false;
737  if(BadTimingChanPercentage>_qualpercent) isTimingOK = false;
738 
739  // Analyze adc samples to get amplitudes
740  //=======================================
741 
742  cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
743  cout << "\t+=+ Analyzing laser data: getting APD, PN, APD/PN, PN/PN +=+" << endl;
744 
745  if( !isGainOK )
746  cout << "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1 +=+" << endl;
747  if( !isTimingOK )
748  cout << "\t+=+ ............................ WARNING! TIMING WAS BAD +=+" << endl;
749 
750  APDFile = new TFile(APDfile.c_str(),"RECREATE");
751 
752  int ieta, iphi;
753 
754  int flagfit;
755  for (unsigned int i=0;i<nCrys;i++){
756 
757  stringstream name;
758  name << "APDTree" <<i+1;
759 
760  APDtrees[i]= new TTree(name.str().c_str(),name.str().c_str());
761 
762 
763  //List of branches
764 
765  APDtrees[i]->Branch( "event", &event, "event/I" );
766  APDtrees[i]->Branch( "color", &color, "color/I" );
767  APDtrees[i]->Branch( "iphi", &iphi, "iphi/I" );
768  APDtrees[i]->Branch( "ieta", &ieta, "ieta/I" );
769  APDtrees[i]->Branch( "side", &side, "side/I" );
770  APDtrees[i]->Branch( "dccID", &dccID, "dccID/I" );
771  APDtrees[i]->Branch( "towerID", &towerID, "towerID/I" );
772  APDtrees[i]->Branch( "channelID", &channelID, "channelID/I" );
773  APDtrees[i]->Branch( "apdAmpl", &apdAmpl, "apdAmpl/D" );
774  APDtrees[i]->Branch( "apdTime", &apdTime, "apdTime/D" );
775  if(_saveallevents) APDtrees[i]->Branch( "adc", &adc ,"adc[10]/D" );
776  APDtrees[i]->Branch( "flagfit", &flagfit, "flagfit/I" );
777  APDtrees[i]->Branch( "pn0", &pn0, "pn0/D" );
778  APDtrees[i]->Branch( "pn1", &pn1, "pn1/D" );
779 
780  APDtrees[i]->SetBranchAddress( "event", &event );
781  APDtrees[i]->SetBranchAddress( "color", &color );
782  APDtrees[i]->SetBranchAddress( "iphi", &iphi );
783  APDtrees[i]->SetBranchAddress( "ieta", &ieta );
784  APDtrees[i]->SetBranchAddress( "side", &side );
785  APDtrees[i]->SetBranchAddress( "dccID", &dccID );
786  APDtrees[i]->SetBranchAddress( "towerID", &towerID );
787  APDtrees[i]->SetBranchAddress( "channelID", &channelID );
788  APDtrees[i]->SetBranchAddress( "apdAmpl", &apdAmpl );
789  APDtrees[i]->SetBranchAddress( "apdTime", &apdTime );
790  if(_saveallevents)APDtrees[i]->SetBranchAddress( "adc", adc );
791  APDtrees[i]->SetBranchAddress( "flagfit", &flagfit );
792  APDtrees[i]->SetBranchAddress( "pn0", &pn0 );
793  APDtrees[i]->SetBranchAddress( "pn1", &pn1 );
794 
795 
796  }
797 
798  for (unsigned int iref=0;iref<nRefChan;iref++){
799  for (unsigned int imod=0;imod<nMod;imod++){
800 
801  int jmod=modules[imod];
802 
803  stringstream nameref;
804  nameref << "refAPDTree" <<imod<<"_"<<iref;
805 
806  RefAPDtrees[iref][jmod]= new TTree(nameref.str().c_str(),nameref.str().c_str());
807 
808  RefAPDtrees[iref][jmod]->Branch( "eventref", &eventref, "eventref/I" );
809  RefAPDtrees[iref][jmod]->Branch( "colorref", &colorref, "colorref/I" );
810  if(iref==0) RefAPDtrees[iref][jmod]->Branch( "apdAmplA", &apdAmplA, "apdAmplA/D" );
811  if(iref==1) RefAPDtrees[iref][jmod]->Branch( "apdAmplB", &apdAmplB, "apdAmplB/D" );
812 
813  RefAPDtrees[iref][jmod]->SetBranchAddress( "eventref", &eventref );
814  RefAPDtrees[iref][jmod]->SetBranchAddress( "colorref", &colorref );
815  if(iref==0)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplA", &apdAmplA );
816  if(iref==1)RefAPDtrees[iref][jmod]->SetBranchAddress( "apdAmplB", &apdAmplB );
817 
818  }
819  }
820 
821 
822  assert( colors.size()<= nColor );
823  unsigned int nCol=colors.size();
824 
825  // Declare PN stuff
826  //===================
827 
828  for (unsigned int iM=0;iM<nMod;iM++){
829  unsigned int iMod=modules[iM]-1;
830  for (unsigned int ich=0;ich<nPNPerMod;ich++){
831  for (unsigned int icol=0;icol<nCol;icol++){
832  PNFirstAnal[iMod][ich][icol]=new TPN(ich);
833  PNAnal[iMod][ich][icol]=new TPN(ich);
834  }
835  }
836  }
837 
838  // Declare function for APD ampl fit
839  //===================================
840 
841  PulseFitWithShape* psfit = new PulseFitWithShape();
842 
843  for (unsigned int iCry=0;iCry<nCrys;iCry++){
844  for (unsigned int icol=0;icol<nCol;icol++){
845 
846  // Declare APD stuff
847  //===================
848 
849  APDFirstAnal[iCry][icol]=new TAPD();
850  IsThereDataADC[iCry][icol]=1;
851  stringstream cut;
852  cut <<"color=="<<colors.at(icol);
853  if(ADCtrees[iCry]->GetEntries(cut.str().c_str())<10) IsThereDataADC[iCry][icol]=0;
854 
855  }
856 
857  unsigned int iMod=iModule[iCry]-1;
858  assert(iMod<=nMod);
859 
861 
862  // Loop on events
863  //================
864 
865  Long64_t nbytes = 0, nb = 0;
866  for (Long64_t jentry=0; jentry< ADCtrees[iCry]->GetEntriesFast();jentry++) { // Loop on events
867  nb = ADCtrees[iCry]->GetEntry(jentry); nbytes += nb;
868 
869  flagfit=1;
870  apdAmpl=0.0;
871  apdTime=0.0;
872  ieta=eta;
873  iphi=phi;
874 
875  // Get back color
876 
877  unsigned int iCol=0;
878  for(unsigned int i=0;i<nCol;i++){
879  if(color==colors[i]) {
880  iCol=i;
881  i=colors.size();
882  }
883  }
884 
885  // Amplitude calculation
886 
889 
890 
891  if(isSPRFine && APDPulse->isPulseOK()) {
892 
893  psfit -> doFit(&adcNoPed[0]);
894  apdAmpl = psfit -> getAmpl();
895  apdTime = psfit -> getTime();
896 
897  }else{
898 
899  apdAmpl=0;
900  apdTime=0;
901  flagfit=0;
902 
903  }
904 
905  if (_debug>=1) cout <<"-- debug test -- endJob -- apdAmpl:"<<apdAmpl
906  <<" apdTime:"<< apdTime<< endl;
907 
908  double pnmean;
909  if (pn0<10 && pn1>10) {
910  pnmean=pn1;
911  }else if (pn1<10 && pn0>10){
912  pnmean=pn0;
913  }else pnmean=0.5*(pn0+pn1);
914 
915  if (_debug>=1) cout <<"-- debug test -- endJob -- pnMean:"<<pnmean << endl;
916 
917  // Fill PN stuff
918  //===============
919 
920  if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){
921  for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
922  PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
923  }
924  }
925 
926  // Fill APD stuff
927  //================
928 
929  if(apdAmpl!=0.0) APDFirstAnal[iCry][iCol]->addEntry(apdAmpl,pnmean, pn0, pn1, apdTime);
930  if (_debug>=1) cout <<"-- debug test -- endJob -- filling APDTree"<< endl;
931 
932  APDtrees[iCry]->Fill();
933 
934  // Fill reference trees
935  //=====================
936 
937  if( apdRefMap[0][iMod+1]==iCry || apdRefMap[1][iMod+1]==iCry) {
938 
939  apdAmplA=0.0;
940  apdAmplB=0.0;
941  eventref=event;
942  colorref=color;
943 
944  for(unsigned int ir=0;ir<nRefChan;ir++){
945 
946  if (_debug>=1) cout <<"-- debug test -- ir:" << ir <<" tt:"<< towerID<<" refmap:"<<apdRefMap[ir][iMod+1]<< " iCry:"<<iCry<<endl;
947 
948  if(apdRefMap[ir][iMod+1]==iCry) {
949  if (_debug>=1) cout <<"-- debug test -- cut passed " <<endl;
950  if (ir==0) apdAmplA=apdAmpl;
951  else if(ir==1) apdAmplB=apdAmpl;
952  if (_debug>=1) cout <<"-- debug test -- apdAmplA=" <<apdAmplA<<endl;
953  if (_debug>=1) cout <<"-- debug test -- apdAmplB=" <<apdAmplB<<endl;
954  if (_debug>=1) cout <<"-- debug test -- color=" <<color<<", event:"<< event<<", ir:" << ir <<" tt-1:"<< towerID-1<< endl;
955 
956  RefAPDtrees[ir][iMod+1]->Fill();
957 
958  if (_debug>=1) cout <<"-- debug test -- tree filled"<< event<<endl;
959  }
960  }
961  }
962  }
963  }
964 
965  delete psfit;
966 
967  ADCFile->Close();
968 
969  if (_debug==1) cout <<"-- debug test -- endJob -- after apdAmpl Loop"<< endl;
970 
971  // Remove temporary file
972  //=======================
973 
974  stringstream del;
975  del << "rm " <<ADCfile;
976  system(del.str().c_str());
977 
978 
979  // Create output trees
980  //=====================
981 
982  resFile = new TFile(resfile.c_str(),"RECREATE");
983 
984  for (unsigned int iColor=0;iColor<nCol;iColor++){
985 
986  stringstream nametree;
987  nametree <<"APDCol"<<colors.at(iColor);
988  stringstream nametree2;
989  nametree2 <<"PNCol"<<colors.at(iColor);
990 
991  restrees[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
992  respntrees[iColor]= new TTree(nametree2.str().c_str(),nametree2.str().c_str());
993 
994  restrees[iColor]->Branch( "iphi", &iphi, "iphi/I" );
995  restrees[iColor]->Branch( "ieta", &ieta, "ieta/I" );
996  restrees[iColor]->Branch( "side", &side, "side/I" );
997  restrees[iColor]->Branch( "dccID", &dccID, "dccID/I" );
998  restrees[iColor]->Branch( "moduleID", &moduleID, "moduleID/I" );
999  restrees[iColor]->Branch( "towerID", &towerID, "towerID/I" );
1000  restrees[iColor]->Branch( "channelID", &channelID, "channelID/I" );
1001  restrees[iColor]->Branch( "APD", &APD, "APD[6]/D" );
1002  restrees[iColor]->Branch( "Time", &Time, "Time[6]/D" );
1003  restrees[iColor]->Branch( "APDoPN", &APDoPN, "APDoPN[6]/D" );
1004  restrees[iColor]->Branch( "APDoPNA", &APDoPNA, "APDoPNA[6]/D" );
1005  restrees[iColor]->Branch( "APDoPNB", &APDoPNB, "APDoPNB[6]/D" );
1006  restrees[iColor]->Branch( "APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D" );
1007  restrees[iColor]->Branch( "APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D" );
1008  restrees[iColor]->Branch( "ShapeCor", &ShapeCor, "ShapeCor/D" );
1009  restrees[iColor]->Branch( "flag", &flag, "flag/I" );
1010 
1011 
1012  respntrees[iColor]->Branch( "moduleID", &moduleID, "moduleID/I" );
1013  respntrees[iColor]->Branch( "pnID", &pnID, "pnID/I" );
1014  respntrees[iColor]->Branch( "PN", &PN, "PN[6]/D" );
1015  respntrees[iColor]->Branch( "PNoPN", &PNoPN, "PNoPN[6]/D" );
1016  respntrees[iColor]->Branch( "PNoPNA", &PNoPNA, "PNoPNA[6]/D" );
1017  respntrees[iColor]->Branch( "PNoPNB", &PNoPNB, "PNoPNB[6]/D" );
1018 
1019  restrees[iColor]->SetBranchAddress( "iphi", &iphi );
1020  restrees[iColor]->SetBranchAddress( "ieta", &ieta );
1021  restrees[iColor]->SetBranchAddress( "dccID", &dccID );
1022  restrees[iColor]->SetBranchAddress( "moduleID", &moduleID );
1023  restrees[iColor]->SetBranchAddress( "towerID", &towerID );
1024  restrees[iColor]->SetBranchAddress( "channelID", &channelID );
1025  restrees[iColor]->SetBranchAddress( "APD", APD );
1026  restrees[iColor]->SetBranchAddress( "Time", Time );
1027  restrees[iColor]->SetBranchAddress( "APDoPN", APDoPN );
1028  restrees[iColor]->SetBranchAddress( "APDoPNA", APDoPNA );
1029  restrees[iColor]->SetBranchAddress( "APDoPNB", APDoPNB );
1030  restrees[iColor]->SetBranchAddress( "APDoAPDA", APDoAPDA );
1031  restrees[iColor]->SetBranchAddress( "APDoAPDB", APDoAPDB );
1032  restrees[iColor]->SetBranchAddress( "ShapeCor", &ShapeCor );
1033  restrees[iColor]->SetBranchAddress( "flag", &flag );
1034 
1035  respntrees[iColor]->SetBranchAddress( "moduleID", &moduleID );
1036  respntrees[iColor]->SetBranchAddress( "pnID", &pnID );
1037  respntrees[iColor]->SetBranchAddress( "PN", PN );
1038  respntrees[iColor]->SetBranchAddress( "PNoPN", PNoPN );
1039  respntrees[iColor]->SetBranchAddress( "PNoPNA", PNoPNA );
1040  respntrees[iColor]->SetBranchAddress( "PNoPNB", PNoPNB );
1041 
1042  }
1043 
1044 
1045  // Set Cuts for PN stuff
1046  //=======================
1047 
1048  for (unsigned int iM=0;iM<nMod;iM++){
1049  unsigned int iMod=modules[iM]-1;
1050 
1051  for (unsigned int ich=0;ich<nPNPerMod;ich++){
1052  for (unsigned int icol=0;icol<nCol;icol++){
1053  PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),PNFirstAnal[iMod][ich][icol]->getPN().at(1));
1054  }
1055  }
1056  }
1057 
1058  // Build ref trees indexes
1059  //========================
1060 
1061  for(unsigned int imod=0;imod<nMod;imod++){
1062  int jmod=modules[imod];
1063  if( RefAPDtrees[0][jmod]->GetEntries()!=0 && RefAPDtrees[1][jmod]->GetEntries()!=0 ){
1064  RefAPDtrees[0][jmod]->BuildIndex("eventref");
1065  RefAPDtrees[1][jmod]->BuildIndex("eventref");
1066  }
1067  }
1068 
1069 
1070  // Final loop on crystals
1071  //=======================
1072 
1073  for (unsigned int iCry=0;iCry<nCrys;iCry++){
1074 
1075  unsigned int iMod=iModule[iCry]-1;
1076 
1077  // Set cuts on APD stuff
1078  //=======================
1079 
1080  for(unsigned int iCol=0;iCol<nCol;iCol++){
1081 
1082  std::vector<double> lowcut;
1083  std::vector<double> highcut;
1084  double cutMin;
1085  double cutMax;
1086 
1087  cutMin=APDFirstAnal[iCry][iCol]->getAPD().at(0)-2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
1088  if(cutMin<0) cutMin=0;
1089  cutMax=APDFirstAnal[iCry][iCol]->getAPD().at(0)+2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
1090 
1091  lowcut.push_back(cutMin);
1092  highcut.push_back(cutMax);
1093 
1094  cutMin=APDFirstAnal[iCry][iCol]->getTime().at(0)-2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
1095  cutMax=APDFirstAnal[iCry][iCol]->getTime().at(0)+2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
1096  lowcut.push_back(cutMin);
1097  highcut.push_back(cutMax);
1098 
1099 
1100  APDAnal[iCry][iCol]=new TAPD();
1101  APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),APDFirstAnal[iCry][iCol]->getAPD().at(1));
1102  APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
1103  APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
1104  APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
1105  APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),APDFirstAnal[iCry][iCol]->getTime().at(1));
1106 
1107 
1108  APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut,highcut);
1109  APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut,highcut);
1110 
1111  }
1112 
1113  // Final loop on events
1114  //=======================
1115 
1116  Long64_t nbytes = 0, nb = 0;
1117  for (Long64_t jentry=0; jentry< APDtrees[iCry]->GetEntriesFast();jentry++) {
1118  nb = APDtrees[iCry]->GetEntry(jentry); nbytes += nb;
1119 
1120  double pnmean;
1121  if (pn0<10 && pn1>10) {
1122  pnmean=pn1;
1123  }else if (pn1<10 && pn0>10){
1124  pnmean=pn0;
1125  }else pnmean=0.5*(pn0+pn1);
1126 
1127  // Get back color
1128  //===============
1129 
1130  unsigned int iCol=0;
1131  for(unsigned int i=0;i<nCol;i++){
1132  if(color==colors[i]) {
1133  iCol=i;
1134  i=colors.size();
1135  }
1136  }
1137 
1138  // Fill PN stuff
1139  //===============
1140 
1141  if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){
1142  for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
1143  PNAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
1144  }
1145  }
1146 
1147  // Get ref amplitudes
1148  //===================
1149 
1150  if (_debug>=1) cout <<"-- debug test -- LastLoop event:"<<event<<" apdAmpl:"<< apdAmpl<< endl;
1151  apdAmplA = 0.0;
1152  apdAmplB = 0.0;
1153 
1154  for (unsigned int iRef=0;iRef<nRefChan;iRef++){
1155  RefAPDtrees[iRef][iMod+1]->GetEntryWithIndex(event);
1156  }
1157 
1158  if (_debug==1 ) cout <<"-- debug test -- LastLoop apdAmplA:"<<apdAmplA<< " apdAmplB:"<< apdAmplB<<", event:"<< event<<", eventref:"<< eventref<< endl;
1159 
1160  // Fill APD stuff
1161  //===============
1162 
1163  APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
1164  }
1165 
1166  moduleID=iMod+1;
1167 
1168  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
1169 
1170  // Get final results for APD
1171  //===========================
1172 
1173  for(unsigned int iColor=0;iColor<nCol;iColor++){
1174 
1175 
1176  std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
1177  std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
1178  std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
1179  std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
1180  std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
1181  std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
1182  std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
1183 
1184 
1185  for(unsigned int i=0;i<apdvec.size();i++){
1186 
1187  APD[i]=apdvec.at(i);
1188  APDoPN[i]=apdpnvec.at(i);
1189  APDoPNA[i]=apdpn0vec.at(i);
1190  APDoPNB[i]=apdpn1vec.at(i);
1191  APDoAPDA[i]=apdapd0vec.at(i);
1192  APDoAPDB[i]=apdapd1vec.at(i);
1193  Time[i]=timevec.at(i);
1195 
1196  }
1197 
1198  // Fill APD results trees
1199  //========================
1200 
1201  iphi=iPhi[iCry];
1202  ieta=iEta[iCry];
1203  dccID=idccID[iCry];
1204  towerID=iTowerID[iCry];
1205  side=iside[iCry];
1206  channelID=iChannelID[iCry];
1207 
1208  if( !wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor]==0 ) flag=1;
1209  else flag=0;
1210 
1211  if (_debug>=1) cout <<"-- debug test -- endJob -- APD[0]"<< APD[0]<<" APDoPN[0] "<<APDoPN[0]<<" APDoAPDA[0] "<<APDoAPDA[0]<< " flag "<< flag<< endl;
1212  restrees[iColor]->Fill();
1213 
1214  }
1215 
1216  }
1217 
1218  // Get final results for PN
1219  //==========================
1220 
1221  for (unsigned int iM=0;iM<nMod;iM++){
1222  unsigned int iMod=modules[iM]-1;
1223 
1224  side=iside[firstChanMod[iMod]];
1225 
1226  for (unsigned int ch=0;ch<nPNPerMod;ch++){
1227 
1228  pnID=ch;
1229  moduleID=iMod+1;
1230 
1231  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
1232 
1233  for(unsigned int iColor=0;iColor<nCol;iColor++){
1234 
1235  std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
1236  std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
1237  std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
1238  std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
1239 
1240  for(unsigned int i=0;i<pnvec.size();i++){
1241 
1242  PN[i]=pnvec.at(i);
1243  PNoPN[i]=pnopnvec.at(i);
1244  PNoPNA[i]=pnopn0vec.at(i);
1245  PNoPNB[i]=pnopn1vec.at(i);
1246  }
1247 
1248  if (_debug>=1) cout <<"-- debug test -- endJob -- filling pn results'tree: PN[0]:"<<PN[0]<<" iModule:" << iMod<<" iColor:"<<iColor<<" ch:"<< ch<< endl;
1249 
1250  // Fill PN results trees
1251  //========================
1252 
1253  respntrees[iColor]->Fill();
1254  }
1255  }
1256  }
1257 
1258  // Remove temporary files
1259  //========================
1260 
1261  if(!_saveallevents){
1262 
1263  APDFile->Close();
1264  stringstream del2;
1265  del2 << "rm " <<APDfile;
1266  system(del2.str().c_str());
1267 
1268  }else {
1269 
1270  APDFile->cd();
1271  APDtrees[0]->Write();
1272  APDFile->Close();
1273  resFile->cd();
1274  }
1275 
1276 
1277  // Save results
1278  //===============
1279 
1280  for (unsigned int i=0;i<nCol;i++){
1281  restrees[i]->Write();
1282  respntrees[i]->Write();
1283  }
1284 
1285  resFile->Close();
1286 
1287  cout << "\t+=+ .................................................. done +=+" << endl;
1288  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
1289 }
1290 
1291 
1292 
1293 //========================================================================
1295 //========================================================================
1296 
1297 
1298  // Get Pulse From Matacq Analysis:
1299  //================================
1300 
1301  bool IsMatacqOK=false;
1302 
1303  int doesMatFileExist=0;
1304  int doesMatShapeExist=0;
1305  FILE *test2;
1306  TProfile *laserShape=nullptr;
1307  test2 = fopen(matfile.c_str(),"r");
1308  if (test2) doesMatFileExist=1;
1309 
1310  TFile *MatShapeFile;
1311  if (doesMatFileExist==1){
1312  MatShapeFile = new TFile(matfile.c_str());
1313  laserShape= (TProfile*) MatShapeFile->Get("shapeLaser");
1314  if(laserShape){
1315  doesMatShapeExist=1;
1316  double y=laserShape->Integral("w");
1317  if(y!=0)laserShape->Scale(1.0/y);
1318  }
1319  }else{
1320 
1321  cout <<" ERROR! Matacq shape file not found !"<< endl;
1322 
1323  }
1324  if (doesMatShapeExist) IsMatacqOK=true;
1325 
1326  // Get SPR from the average elec shape in SM6:
1327  //============================================
1328 
1329  int doesElecFileExist=0;
1330  FILE *test;
1331  test = fopen(elecfile_.c_str(),"r");
1332  if (test) doesElecFileExist=1;
1333 
1334  TFile *ElecShapesFile;
1335  TH1D* elecShape=nullptr ;
1336 
1337  if (doesElecFileExist==1){
1338  ElecShapesFile = new TFile(elecfile_.c_str());
1339  stringstream name;
1340  name << "MeanElecShape";
1341  elecShape=(TH1D*) ElecShapesFile->Get(name.str().c_str());
1342  if(elecShape && doesMatShapeExist==1){
1343  double x=elecShape->GetMaximum();
1344  if (x!=0) elecShape->Scale(1.0/x);
1345  isSPRFine=true;
1346  }else{
1347  isSPRFine=false;
1348  }
1349 
1350  }else{
1351 
1352  cout <<" ERROR! Elec shape file not found !"<< endl;
1353 
1354  }
1355 
1356 
1357  if(IsMatacqOK){
1358 
1359  ShapeFile = new TFile(shapefile.c_str(),"RECREATE");
1360 
1361  unsigned int nBins=int(laserShape->GetEntries());
1362  assert( nSamplesShapes == nBins);
1363  double elec_jj, laser_iiMinusjj;
1364  double sum_jj;
1365 
1366  if(isSPRFine==true){
1367 
1368  unsigned int nBins2=int(elecShape->GetNbinsX());
1369 
1370  if(nBins2<nBins){
1371  cout<< "EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins"<< std::endl;
1372  isSPRFine=false;
1373  }
1374  assert( nSamplesShapes == nBins2 );
1375 
1376  stringstream name;
1377  name << "PulseShape";
1378 
1379  PulseShape=new TProfile(name.str().c_str(),name.str().c_str(),nBins,-0.5,double(nBins)-0.5);
1380 
1381  // shift shapes to have max close to the real APD max
1382 
1383  for(int ii=0;ii<50;ii++){
1384  shapes[ii]=0.0;
1385  PulseShape->Fill(ii,0.0);
1386  }
1387 
1388  for(unsigned int ii=0;ii<nBins-50;ii++){
1389  sum_jj=0.0;
1390  for(unsigned int jj=0;jj<ii;jj++){
1391  elec_jj=elecShape->GetBinContent(jj+1);
1392  laser_iiMinusjj=laserShape->GetBinContent(ii-jj+1);
1393  sum_jj+=elec_jj*laser_iiMinusjj;
1394  }
1395  PulseShape->Fill(ii+50,sum_jj);
1396  shapes[ii+50]=sum_jj;
1397  }
1398 
1399  double scale= PulseShape->GetMaximum();
1401 
1402  if(scale!=0) {
1403  PulseShape->Scale(1.0/scale);
1404  for(unsigned int ii=0;ii<nBins;ii++){
1405  shapesVec.push_back(shapes[ii]/scale);
1406  }
1407  }
1408 
1409  if(_saveshapes) PulseShape->Write();
1410  }
1411  }
1412  ShapeFile->Close();
1413 
1414  if(!_saveshapes) {
1415 
1416  stringstream del;
1417  del << "rm " <<shapefile;
1418  system(del.str().c_str());
1419 
1420  }
1421 
1422  return IsMatacqOK;
1423 }
1424 
1425 
1426 void EcalLaserAnalyzer2::setGeomEB(int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr){
1427 
1428  side=MEEBGeom::side(etaG,phiG);
1429 
1430  assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
1431 
1432  eta = etaG;
1433  phi = phiG;
1434  channelID=5*(strip-1) + xtal-1;
1435  towerID=tower;
1436 
1437  vector<int> apdRefChan=ME::apdRefChannels(module, lmr);
1438  for (unsigned int iref=0;iref<nRefChan;iref++){
1439  if(channelID==apdRefChan[iref] && towerID==apdRefTT
1440  && apdRefMap[iref].count(module)==0){
1441  apdRefMap[iref][module]=channel;
1442  }
1443  }
1444 
1445  if(isFirstChanModFilled[module-1]==0) {
1446  firstChanMod[module-1]=channel;
1447  isFirstChanModFilled[module-1]=1;
1448  }
1449 
1450  iEta[channel]=eta;
1451  iPhi[channel]=phi;
1452  iModule[channel]= module ;
1453  iTowerID[channel]=towerID;
1454  iChannelID[channel]=channelID;
1455  idccID[channel]=dccID;
1456  iside[channel]=side;
1457 
1458 }
1459 
1460 void EcalLaserAnalyzer2::setGeomEE(int etaG, int phiG,int iX, int iY, int iZ, int module, int tower, int ch , int apdRefTT, int channel, int lmr){
1461 
1462  side=MEEEGeom::side(iX, iY, iZ);
1463 
1464  assert( module>=*min_element(modules.begin(),modules.end())
1465  && module<=*max_element(modules.begin(),modules.end()) );
1466 
1467  eta = etaG;
1468  phi = phiG;
1469  channelID=ch;
1470  towerID=tower;
1471 
1472  vector<int> apdRefChan=ME::apdRefChannels(module, lmr);
1473  for (unsigned int iref=0;iref<nRefChan;iref++){
1474  if(channelID==apdRefChan[iref] && towerID==apdRefTT
1475  && apdRefMap[iref].count(module)==0){
1476  apdRefMap[iref][module]=channel;
1477  }
1478  }
1479 
1480  if(isFirstChanModFilled[module-1]==0) {
1481  firstChanMod[module-1]=channel;
1482  isFirstChanModFilled[module-1]=1;
1483  }
1484 
1485  iEta[channel]=eta;
1486  iPhi[channel]=phi;
1487  iModule[channel]= module ;
1488  iTowerID[channel]=towerID;
1489  iChannelID[channel]=channelID;
1490  idccID[channel]=dccID;
1491  iside[channel]=side;
1492 
1493 }
1495 
size
Write out results.
void addEntry(double val)
Definition: TMom.cc:110
static XYCoord localCoord(int icr)
Definition: MEEBGeom.cc:153
T getParameter(std::string const &) const
TPN * PNAnal[9][2][nColor]
T getUntrackedParameter(std::string const &, T const &) const
static int lmmod(SuperCrysCoord iX, SuperCrysCoord iY)
Definition: MEEEGeom.cc:94
void analyze(const edm::Event &e, const edm::EventSetup &c) override
Definition: TPN.h:8
void setAPDCut(double, double)
Definition: TAPD.cc:147
std::map< int, int > channelMapEE
int xtalId() const
get the channel id
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
std::vector< double > getPN()
Definition: TPN.cc:96
#define NSIDES
Definition: TAPD.h:8
static int apdRefTower(int ilmr, int ilmmod)
Definition: MEEEGeom.cc:1006
void addEntry(double, double, double)
Definition: TPN.cc:40
void beginJob() override
int stripId() const
get the tower id
double getDelta(int, int)
Definition: TAPDPulse.cc:95
TAPD * APDAnal[1700][nColor]
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< double > getPNoPN1()
Definition: TPN.cc:99
Ecal readout channel identification [32:20] Unused (so far) [19:13] DCC id [12:6] tower [5:3] strip [...
#define NCRYSEE
bool setPulse(double *)
Definition: TAPDPulse.cc:64
TAPD * APDFirstAnal[1700][nColor]
static int side(SuperCrysCoord iX, SuperCrysCoord iY, int iz)
Definition: MEEEGeom.cc:939
TPN * PNFirstAnal[9][2][nColor]
void setPNCut(double, double)
Definition: TPN.cc:72
int init
Definition: HydjetWrapper.h:67
void addEntry(double, double, double, double, double, double, double)
Definition: TAPD.cc:53
unsigned int _timingcuthigh
double getPNCorrectionFactor(double val0, int gain)
Definition: TPNCor.cc:72
std::vector< EcalDCCHeaderBlock >::const_iterator const_iterator
int towerId() const
get the tower id
double * getAdcWithoutPedestal()
Definition: TPNPulse.cc:99
const_iterator begin() const
unsigned int isFirstChanModFilled[21]
void setGeomEE(int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr)
static int lmr(EBGlobalCoord ieta, EBGlobalCoord iphi)
Definition: MEEBGeom.cc:121
Definition: TMom.h:7
unsigned int _timingquallow
#define nullptr
std::string digiPNCollection_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
bool setPulse(double *)
Definition: TPNPulse.cc:55
std::vector< double > getAPDoPN0()
Definition: TAPD.cc:237
void setTimeCut(double, double)
Definition: TAPD.cc:151
unsigned int _presamplePN
static std::pair< int, int > pn(int ilmmod)
Definition: MEEBGeom.cc:479
void set2DAPDoAPD1Cut(const std::vector< double > &, const std::vector< double > &)
Definition: TAPD.cc:192
void setAPDoPN0Cut(double, double)
Definition: TAPD.cc:149
std::vector< double > getPNoPN0()
Definition: TPN.cc:98
std::vector< double > getAPDoPN1()
Definition: TAPD.cc:238
int iPnId() const
get the PnId
Definition: TMem.h:7
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:42
EcalElectronicsId getElectronicsId(const DetId &id) const
Get the electronics id for this det id.
static std::pair< int, int > pn(int dee, int ilmod)
Definition: MEEEGeom.cc:615
int Mem(int, int)
Definition: TMem.cc:52
int iEvent
Definition: GenABIO.cc:230
std::vector< double > getAPDoAPD0()
Definition: TAPD.cc:240
void set2DAPDoAPD0Cut(const std::vector< double > &, const std::vector< double > &)
Definition: TAPD.cc:183
static int electronic_channel(EBLocalCoord ix, EBLocalCoord iy)
Definition: MEEBGeom.cc:341
#define NCRYSEB
std::string eventHeaderCollection_
unsigned int _nsamplesPN
#define NMODEB
double * getAdcWithoutPedestal()
Definition: TAPDPulse.cc:209
EcalLaserAnalyzer2(const edm::ParameterSet &iConfig)
static int apdRefTower(int ilmmod)
Definition: MEEBGeom.cc:520
#define NREFCHAN
unsigned int _lastsample
int iDCCId() const
get the DCCId
static int lmmod(EBGlobalCoord ieta, EBGlobalCoord iphi)
Definition: MEEBGeom.cc:93
std::vector< double > getAPDoAPD1()
Definition: TAPD.cc:245
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:475
JetCorrectorParameters corr
Definition: classes.h:5
std::vector< double > getAPD()
Definition: TAPD.cc:235
ii
Definition: cuy.py:588
const_iterator end() const
T Max(T a, T b)
Definition: MathUtil.h:44
unsigned int firstChanMod[21]
#define NSAMPSHAPES
void setGeomEB(int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr)
static int side(EBGlobalCoord ieta, EBGlobalCoord iphi)
Definition: MEEBGeom.cc:114
bool isTimingQualOK()
Definition: TAPDPulse.cc:123
T const * product() const
Definition: Handle.h:81
unsigned int nSamplesShapes
int IsThereDataADC[1700][nColor]
const T & get() const
Definition: EventSetup.h:58
void endJob() override
TTree * respntrees[nColor]
TTree * RefAPDtrees[2][21]
#define NPNPERMOD
unsigned int _firstsamplePN
const_iterator end() const
std::vector< double > shapesVec
void setPresamples(int)
Definition: TAPDPulse.cc:222
bool isPulseOK()
Definition: TAPDPulse.cc:139
Definition: colors.py:1
TTree * restrees[nColor]
std::vector< double > getPNoPN()
Definition: TPN.cc:97
Definition: TPNFit.h:8
static int dee(SuperCrysCoord iX, SuperCrysCoord iY, int iz)
Definition: MEEEGeom.cc:284
bool isMemRelevant(int)
Definition: TMem.cc:40
unsigned int _firstsample
int getMaxSample()
Definition: TPNPulse.cc:81
unsigned int _timingcutlow
unsigned int iModule[1700]
EcalLogicID towerID(EcalElectronicsId const &)
std::string eventHeaderProducer_
std::vector< double > getAPDoPN()
Definition: TAPD.cc:236
Definition: TPNCor.h:7
std::map< int, unsigned int > apdRefMap[2]
std::vector< double > getTime()
Definition: TAPD.cc:239
T const * product() const
Definition: ESHandle.h:86
Definition: vlib.h:208
unsigned int _timingqualhigh
static std::vector< int > apdRefChannels(ME::LMMid ilmmod, ME::LMRid ilmr)
Definition: ME.cc:588
unsigned int _lastsamplePN
static std::vector< ME::LMMid > lmmodFromDcc(ME::DCCid idcc)
Definition: ME.cc:623
double getMean()
Definition: TMom.cc:147
static int lmr(SuperCrysCoord iX, SuperCrysCoord iY, int iz)
Definition: MEEEGeom.cc:250
int channelId() const
so far for EndCap only :
const_iterator begin() const
void setAPDoPNCut(double, double)
Definition: TAPD.cc:148
Definition: event.py:1
void setAPDoPN1Cut(double, double)
Definition: TAPD.cc:150