CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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_.c_str());
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=0;
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=0;
293  const EEDigiCollection* EEDigi=0;
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=0;
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=0;
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=0;
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=0 ;
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 
void addEntry(double val)
Definition: TMom.cc:110
static XYCoord localCoord(int icr)
Definition: MEEBGeom.cc:153
T getParameter(std::string const &) const
TTree * RefAPDtrees[NREFCHAN][NMODEE]
T getUntrackedParameter(std::string const &, T const &) const
static int lmmod(SuperCrysCoord iX, SuperCrysCoord iY)
Definition: MEEEGeom.cc:94
int i
Definition: DBlmapReader.cc:9
virtual void analyze(const edm::Event &e, const edm::EventSetup &c)
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
int stripId() const
get the tower id
double getDelta(int, int)
Definition: TAPDPulse.cc:95
#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 [...
TAPD * APDAnal[NCRYSEB][nColor]
#define NCRYSEE
bool setPulse(double *)
Definition: TAPDPulse.cc:64
static int side(SuperCrysCoord iX, SuperCrysCoord iY, int iz)
Definition: MEEEGeom.cc:939
void setPNCut(double, double)
Definition: TPN.cc:72
int init
Definition: HydjetWrapper.h:62
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
TTree * ADCtrees[NCRYSEB]
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
T eta() const
Definition: TMom.h:7
unsigned int _timingquallow
std::string digiPNCollection_
int ii
Definition: cuy.py:588
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
bool setPulse(double *)
Definition: TPNPulse.cc:55
std::vector< double > getAPDoPN0()
Definition: TAPD.cc:237
void setTimeCut(double, double)
Definition: TAPD.cc:151
TPN * PNAnal[NMODEB][NPNPERMOD][nColor]
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
double shapes[NSAMPSHAPES]
#define NCRYSEB
std::string eventHeaderCollection_
unsigned int _nsamplesPN
#define NMODEB
int IsThereDataADC[NCRYSEB][nColor]
double * getAdcWithoutPedestal()
Definition: TAPDPulse.cc:209
EcalLaserAnalyzer2(const edm::ParameterSet &iConfig)
int j
Definition: DBlmapReader.cc:9
static int apdRefTower(int ilmmod)
Definition: MEEBGeom.cc:520
#define NREFCHAN
unsigned int _lastsample
std::vector< int > modules
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
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:390
JetCorrectorParameters corr
Definition: classes.h:5
bool wasTimingOK[NCRYSEB]
std::vector< double > getAPD()
Definition: TAPD.cc:235
const_iterator end() const
TTree * APDtrees[NCRYSEB]
#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
unsigned int firstChanMod[NMODEE]
TAPD * APDFirstAnal[NCRYSEB][nColor]
bool isTimingQualOK()
Definition: TAPDPulse.cc:123
unsigned int nSamplesShapes
unsigned int iModule[NCRYSEB]
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
unsigned int isFirstChanModFilled[NMODEE]
TTree * respntrees[nColor]
#define NPNPERMOD
T const * product() const
Definition: Handle.h:81
unsigned int _firstsamplePN
const_iterator end() const
std::vector< double > shapesVec
void setPresamples(int)
Definition: TAPDPulse.cc:222
bool isPulseOK()
Definition: TAPDPulse.cc:139
TTree * restrees[nColor]
std::vector< double > getPNoPN()
Definition: TPN.cc:97
static double getTime()
Definition: TPNFit.h:8
std::vector< int > colors
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
tuple cout
Definition: gather_cfg.py:121
EcalLogicID towerID(EcalElectronicsId const &)
std::string eventHeaderProducer_
std::vector< double > getAPDoPN()
Definition: TAPD.cc:236
volatile std::atomic< bool > shutdown_flag false
Definition: DDAxes.h:10
Definition: TPNCor.h:7
TPN * PNFirstAnal[NMODEB][NPNPERMOD][nColor]
std::map< int, unsigned int > apdRefMap[2]
std::vector< double > getTime()
Definition: TAPD.cc:239
int nEvtBadTiming[NCRYSEB]
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
tuple size
Write out results.
int channelId() const
so far for EndCap only :
const_iterator begin() const
void setAPDoPNCut(double, double)
Definition: TAPD.cc:148
void setAPDoPN1Cut(double, double)
Definition: TAPD.cc:150
list at
Definition: asciidump.py:428
Definition: DDAxes.h:10