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