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: 2012/02/09 10:07:36 $
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 
844  for (unsigned int iCry=0;iCry<nCrys;iCry++){
845  for (unsigned int icol=0;icol<nCol;icol++){
846 
847  // Declare APD stuff
848  //===================
849 
850  APDFirstAnal[iCry][icol]=new TAPD();
851  IsThereDataADC[iCry][icol]=1;
852  stringstream cut;
853  cut <<"color=="<<colors.at(icol);
854  if(ADCtrees[iCry]->GetEntries(cut.str().c_str())<10) IsThereDataADC[iCry][icol]=0;
855 
856  }
857 
858  unsigned int iMod=iModule[iCry]-1;
859  assert(iMod<=nMod);
860 
862 
863  // Loop on events
864  //================
865 
866  Long64_t nbytes = 0, nb = 0;
867  for (Long64_t jentry=0; jentry< ADCtrees[iCry]->GetEntriesFast();jentry++) { // Loop on events
868  nb = ADCtrees[iCry]->GetEntry(jentry); nbytes += nb;
869 
870  flagfit=1;
871  apdAmpl=0.0;
872  apdTime=0.0;
873  ieta=eta;
874  iphi=phi;
875 
876  // Get back color
877 
878  unsigned int iCol=0;
879  for(unsigned int i=0;i<nCol;i++){
880  if(color==colors[i]) {
881  iCol=i;
882  i=colors.size();
883  }
884  }
885 
886  // Amplitude calculation
887 
890 
891 
892  if(isSPRFine && APDPulse->isPulseOK()) {
893 
894  psfit -> doFit(&adcNoPed[0]);
895  apdAmpl = psfit -> getAmpl();
896  apdTime = psfit -> getTime();
897 
898  }else{
899 
900  apdAmpl=0;
901  apdTime=0;
902  flagfit=0;
903 
904  }
905 
906  if (_debug>=1) cout <<"-- debug test -- endJob -- apdAmpl:"<<apdAmpl
907  <<" apdTime:"<< apdTime<< endl;
908 
909  double pnmean;
910  if (pn0<10 && pn1>10) {
911  pnmean=pn1;
912  }else if (pn1<10 && pn0>10){
913  pnmean=pn0;
914  }else pnmean=0.5*(pn0+pn1);
915 
916  if (_debug>=1) cout <<"-- debug test -- endJob -- pnMean:"<<pnmean << endl;
917 
918  // Fill PN stuff
919  //===============
920 
921  if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){
922  for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
923  PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
924  }
925  }
926 
927  // Fill APD stuff
928  //================
929 
930  if(apdAmpl!=0.0) APDFirstAnal[iCry][iCol]->addEntry(apdAmpl,pnmean, pn0, pn1, apdTime);
931  if (_debug>=1) cout <<"-- debug test -- endJob -- filling APDTree"<< endl;
932 
933  APDtrees[iCry]->Fill();
934 
935  // Fill reference trees
936  //=====================
937 
938  if( apdRefMap[0][iMod+1]==iCry || apdRefMap[1][iMod+1]==iCry) {
939 
940  apdAmplA=0.0;
941  apdAmplB=0.0;
942  eventref=event;
943  colorref=color;
944 
945  for(unsigned int ir=0;ir<nRefChan;ir++){
946 
947  if (_debug>=1) cout <<"-- debug test -- ir:" << ir <<" tt:"<< towerID<<" refmap:"<<apdRefMap[ir][iMod+1]<< " iCry:"<<iCry<<endl;
948 
949  if(apdRefMap[ir][iMod+1]==iCry) {
950  if (_debug>=1) cout <<"-- debug test -- cut passed " <<endl;
951  if (ir==0) apdAmplA=apdAmpl;
952  else if(ir==1) apdAmplB=apdAmpl;
953  if (_debug>=1) cout <<"-- debug test -- apdAmplA=" <<apdAmplA<<endl;
954  if (_debug>=1) cout <<"-- debug test -- apdAmplB=" <<apdAmplB<<endl;
955  if (_debug>=1) cout <<"-- debug test -- color=" <<color<<", event:"<< event<<", ir:" << ir <<" tt-1:"<< towerID-1<< endl;
956 
957  RefAPDtrees[ir][iMod+1]->Fill();
958 
959  if (_debug>=1) cout <<"-- debug test -- tree filled"<< event<<endl;
960  }
961  }
962  }
963  }
964  }
965 
966  delete psfit;
967 
968  ADCFile->Close();
969 
970  if (_debug==1) cout <<"-- debug test -- endJob -- after apdAmpl Loop"<< endl;
971 
972  // Remove temporary file
973  //=======================
974 
975  stringstream del;
976  del << "rm " <<ADCfile;
977  system(del.str().c_str());
978 
979 
980  // Create output trees
981  //=====================
982 
983  resFile = new TFile(resfile.c_str(),"RECREATE");
984 
985  for (unsigned int iColor=0;iColor<nCol;iColor++){
986 
987  stringstream nametree;
988  nametree <<"APDCol"<<colors.at(iColor);
989  stringstream nametree2;
990  nametree2 <<"PNCol"<<colors.at(iColor);
991 
992  restrees[iColor]= new TTree(nametree.str().c_str(),nametree.str().c_str());
993  respntrees[iColor]= new TTree(nametree2.str().c_str(),nametree2.str().c_str());
994 
995  restrees[iColor]->Branch( "iphi", &iphi, "iphi/I" );
996  restrees[iColor]->Branch( "ieta", &ieta, "ieta/I" );
997  restrees[iColor]->Branch( "side", &side, "side/I" );
998  restrees[iColor]->Branch( "dccID", &dccID, "dccID/I" );
999  restrees[iColor]->Branch( "moduleID", &moduleID, "moduleID/I" );
1000  restrees[iColor]->Branch( "towerID", &towerID, "towerID/I" );
1001  restrees[iColor]->Branch( "channelID", &channelID, "channelID/I" );
1002  restrees[iColor]->Branch( "APD", &APD, "APD[6]/D" );
1003  restrees[iColor]->Branch( "Time", &Time, "Time[6]/D" );
1004  restrees[iColor]->Branch( "APDoPN", &APDoPN, "APDoPN[6]/D" );
1005  restrees[iColor]->Branch( "APDoPNA", &APDoPNA, "APDoPNA[6]/D" );
1006  restrees[iColor]->Branch( "APDoPNB", &APDoPNB, "APDoPNB[6]/D" );
1007  restrees[iColor]->Branch( "APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D" );
1008  restrees[iColor]->Branch( "APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D" );
1009  restrees[iColor]->Branch( "ShapeCor", &ShapeCor, "ShapeCor/D" );
1010  restrees[iColor]->Branch( "flag", &flag, "flag/I" );
1011 
1012 
1013  respntrees[iColor]->Branch( "moduleID", &moduleID, "moduleID/I" );
1014  respntrees[iColor]->Branch( "pnID", &pnID, "pnID/I" );
1015  respntrees[iColor]->Branch( "PN", &PN, "PN[6]/D" );
1016  respntrees[iColor]->Branch( "PNoPN", &PNoPN, "PNoPN[6]/D" );
1017  respntrees[iColor]->Branch( "PNoPNA", &PNoPNA, "PNoPNA[6]/D" );
1018  respntrees[iColor]->Branch( "PNoPNB", &PNoPNB, "PNoPNB[6]/D" );
1019 
1020  restrees[iColor]->SetBranchAddress( "iphi", &iphi );
1021  restrees[iColor]->SetBranchAddress( "ieta", &ieta );
1022  restrees[iColor]->SetBranchAddress( "dccID", &dccID );
1023  restrees[iColor]->SetBranchAddress( "moduleID", &moduleID );
1024  restrees[iColor]->SetBranchAddress( "towerID", &towerID );
1025  restrees[iColor]->SetBranchAddress( "channelID", &channelID );
1026  restrees[iColor]->SetBranchAddress( "APD", APD );
1027  restrees[iColor]->SetBranchAddress( "Time", Time );
1028  restrees[iColor]->SetBranchAddress( "APDoPN", APDoPN );
1029  restrees[iColor]->SetBranchAddress( "APDoPNA", APDoPNA );
1030  restrees[iColor]->SetBranchAddress( "APDoPNB", APDoPNB );
1031  restrees[iColor]->SetBranchAddress( "APDoAPDA", APDoAPDA );
1032  restrees[iColor]->SetBranchAddress( "APDoAPDB", APDoAPDB );
1033  restrees[iColor]->SetBranchAddress( "ShapeCor", &ShapeCor );
1034  restrees[iColor]->SetBranchAddress( "flag", &flag );
1035 
1036  respntrees[iColor]->SetBranchAddress( "moduleID", &moduleID );
1037  respntrees[iColor]->SetBranchAddress( "pnID", &pnID );
1038  respntrees[iColor]->SetBranchAddress( "PN", PN );
1039  respntrees[iColor]->SetBranchAddress( "PNoPN", PNoPN );
1040  respntrees[iColor]->SetBranchAddress( "PNoPNA", PNoPNA );
1041  respntrees[iColor]->SetBranchAddress( "PNoPNB", PNoPNB );
1042 
1043  }
1044 
1045 
1046  // Set Cuts for PN stuff
1047  //=======================
1048 
1049  for (unsigned int iM=0;iM<nMod;iM++){
1050  unsigned int iMod=modules[iM]-1;
1051 
1052  for (unsigned int ich=0;ich<nPNPerMod;ich++){
1053  for (unsigned int icol=0;icol<nCol;icol++){
1054  PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),PNFirstAnal[iMod][ich][icol]->getPN().at(1));
1055  }
1056  }
1057  }
1058 
1059  // Build ref trees indexes
1060  //========================
1061 
1062  for(unsigned int imod=0;imod<nMod;imod++){
1063  int jmod=modules[imod];
1064  if( RefAPDtrees[0][jmod]->GetEntries()!=0 && RefAPDtrees[1][jmod]->GetEntries()!=0 ){
1065  RefAPDtrees[0][jmod]->BuildIndex("eventref");
1066  RefAPDtrees[1][jmod]->BuildIndex("eventref");
1067  }
1068  }
1069 
1070 
1071  // Final loop on crystals
1072  //=======================
1073 
1074  for (unsigned int iCry=0;iCry<nCrys;iCry++){
1075 
1076  unsigned int iMod=iModule[iCry]-1;
1077 
1078  // Set cuts on APD stuff
1079  //=======================
1080 
1081  for(unsigned int iCol=0;iCol<nCol;iCol++){
1082 
1083  std::vector<double> lowcut;
1084  std::vector<double> highcut;
1085  double cutMin;
1086  double cutMax;
1087 
1088  cutMin=APDFirstAnal[iCry][iCol]->getAPD().at(0)-2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
1089  if(cutMin<0) cutMin=0;
1090  cutMax=APDFirstAnal[iCry][iCol]->getAPD().at(0)+2.0*APDFirstAnal[iCry][iCol]->getAPD().at(1);
1091 
1092  lowcut.push_back(cutMin);
1093  highcut.push_back(cutMax);
1094 
1095  cutMin=APDFirstAnal[iCry][iCol]->getTime().at(0)-2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
1096  cutMax=APDFirstAnal[iCry][iCol]->getTime().at(0)+2.0*APDFirstAnal[iCry][iCol]->getTime().at(1);
1097  lowcut.push_back(cutMin);
1098  highcut.push_back(cutMax);
1099 
1100 
1101  APDAnal[iCry][iCol]=new TAPD();
1102  APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),APDFirstAnal[iCry][iCol]->getAPD().at(1));
1103  APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
1104  APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
1105  APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
1106  APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),APDFirstAnal[iCry][iCol]->getTime().at(1));
1107 
1108 
1109  APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut,highcut);
1110  APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut,highcut);
1111 
1112  }
1113 
1114  // Final loop on events
1115  //=======================
1116 
1117  Long64_t nbytes = 0, nb = 0;
1118  for (Long64_t jentry=0; jentry< APDtrees[iCry]->GetEntriesFast();jentry++) {
1119  nb = APDtrees[iCry]->GetEntry(jentry); nbytes += nb;
1120 
1121  double pnmean;
1122  if (pn0<10 && pn1>10) {
1123  pnmean=pn1;
1124  }else if (pn1<10 && pn0>10){
1125  pnmean=pn0;
1126  }else pnmean=0.5*(pn0+pn1);
1127 
1128  // Get back color
1129  //===============
1130 
1131  unsigned int iCol=0;
1132  for(unsigned int i=0;i<nCol;i++){
1133  if(color==colors[i]) {
1134  iCol=i;
1135  i=colors.size();
1136  }
1137  }
1138 
1139  // Fill PN stuff
1140  //===============
1141 
1142  if( firstChanMod[iMod]==iCry && IsThereDataADC[iCry][iCol]==1 ){
1143  for (unsigned int ichan=0;ichan<nPNPerMod;ichan++){
1144  PNAnal[iMod][ichan][iCol]->addEntry(pnmean,pn0,pn1);
1145  }
1146  }
1147 
1148  // Get ref amplitudes
1149  //===================
1150 
1151  if (_debug>=1) cout <<"-- debug test -- LastLoop event:"<<event<<" apdAmpl:"<< apdAmpl<< endl;
1152  apdAmplA = 0.0;
1153  apdAmplB = 0.0;
1154 
1155  for (unsigned int iRef=0;iRef<nRefChan;iRef++){
1156  RefAPDtrees[iRef][iMod+1]->GetEntryWithIndex(event);
1157  }
1158 
1159  if (_debug==1 ) cout <<"-- debug test -- LastLoop apdAmplA:"<<apdAmplA<< " apdAmplB:"<< apdAmplB<<", event:"<< event<<", eventref:"<< eventref<< endl;
1160 
1161  // Fill APD stuff
1162  //===============
1163 
1164  APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
1165  }
1166 
1167  moduleID=iMod+1;
1168 
1169  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
1170 
1171  // Get final results for APD
1172  //===========================
1173 
1174  for(unsigned int iColor=0;iColor<nCol;iColor++){
1175 
1176 
1177  std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
1178  std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
1179  std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
1180  std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
1181  std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
1182  std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
1183  std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
1184 
1185 
1186  for(unsigned int i=0;i<apdvec.size();i++){
1187 
1188  APD[i]=apdvec.at(i);
1189  APDoPN[i]=apdpnvec.at(i);
1190  APDoPNA[i]=apdpn0vec.at(i);
1191  APDoPNB[i]=apdpn1vec.at(i);
1192  APDoAPDA[i]=apdapd0vec.at(i);
1193  APDoAPDB[i]=apdapd1vec.at(i);
1194  Time[i]=timevec.at(i);
1196 
1197  }
1198 
1199  // Fill APD results trees
1200  //========================
1201 
1202  iphi=iPhi[iCry];
1203  ieta=iEta[iCry];
1204  dccID=idccID[iCry];
1205  towerID=iTowerID[iCry];
1206  side=iside[iCry];
1207  channelID=iChannelID[iCry];
1208 
1209  if( !wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor]==0 ) flag=1;
1210  else flag=0;
1211 
1212  if (_debug>=1) cout <<"-- debug test -- endJob -- APD[0]"<< APD[0]<<" APDoPN[0] "<<APDoPN[0]<<" APDoAPDA[0] "<<APDoAPDA[0]<< " flag "<< flag<< endl;
1213  restrees[iColor]->Fill();
1214 
1215  }
1216 
1217  }
1218 
1219  // Get final results for PN
1220  //==========================
1221 
1222  for (unsigned int iM=0;iM<nMod;iM++){
1223  unsigned int iMod=modules[iM]-1;
1224 
1225  side=iside[firstChanMod[iMod]];
1226 
1227  for (unsigned int ch=0;ch<nPNPerMod;ch++){
1228 
1229  pnID=ch;
1230  moduleID=iMod+1;
1231 
1232  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
1233 
1234  for(unsigned int iColor=0;iColor<nCol;iColor++){
1235 
1236  std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
1237  std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
1238  std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
1239  std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
1240 
1241  for(unsigned int i=0;i<pnvec.size();i++){
1242 
1243  PN[i]=pnvec.at(i);
1244  PNoPN[i]=pnopnvec.at(i);
1245  PNoPNA[i]=pnopn0vec.at(i);
1246  PNoPNB[i]=pnopn1vec.at(i);
1247  }
1248 
1249  if (_debug>=1) cout <<"-- debug test -- endJob -- filling pn results'tree: PN[0]:"<<PN[0]<<" iModule:" << iMod<<" iColor:"<<iColor<<" ch:"<< ch<< endl;
1250 
1251  // Fill PN results trees
1252  //========================
1253 
1254  respntrees[iColor]->Fill();
1255  }
1256  }
1257  }
1258 
1259  // Remove temporary files
1260  //========================
1261 
1262  if(!_saveallevents){
1263 
1264  APDFile->Close();
1265  stringstream del2;
1266  del2 << "rm " <<APDfile;
1267  system(del2.str().c_str());
1268 
1269  }else {
1270 
1271  APDFile->cd();
1272  APDtrees[0]->Write();
1273  APDFile->Close();
1274  resFile->cd();
1275  }
1276 
1277 
1278  // Save results
1279  //===============
1280 
1281  for (unsigned int i=0;i<nCol;i++){
1282  restrees[i]->Write();
1283  respntrees[i]->Write();
1284  }
1285 
1286  resFile->Close();
1287 
1288  cout << "\t+=+ .................................................. done +=+" << endl;
1289  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
1290 }
1291 
1292 
1293 
1294 //========================================================================
1296 //========================================================================
1297 
1298 
1299  // Get Pulse From Matacq Analysis:
1300  //================================
1301 
1302  bool IsMatacqOK=false;
1303 
1304  int doesMatFileExist=0;
1305  int doesMatShapeExist=0;
1306  FILE *test2;
1307  TProfile *laserShape=0;
1308  test2 = fopen(matfile.c_str(),"r");
1309  if (test2) doesMatFileExist=1;
1310 
1311  TFile *MatShapeFile;
1312  if (doesMatFileExist==1){
1313  MatShapeFile = new TFile(matfile.c_str());
1314  laserShape= (TProfile*) MatShapeFile->Get("shapeLaser");
1315  if(laserShape){
1316  doesMatShapeExist=1;
1317  double y=laserShape->Integral("w");
1318  if(y!=0)laserShape->Scale(1.0/y);
1319  }
1320  }else{
1321 
1322  cout <<" ERROR! Matacq shape file not found !"<< endl;
1323 
1324  }
1325  if (doesMatShapeExist) IsMatacqOK=true;
1326 
1327  // Get SPR from the average elec shape in SM6:
1328  //============================================
1329 
1330  int doesElecFileExist=0;
1331  FILE *test;
1332  test = fopen(elecfile_.c_str(),"r");
1333  if (test) doesElecFileExist=1;
1334 
1335  TFile *ElecShapesFile;
1336  TH1D* elecShape=0 ;
1337 
1338  if (doesElecFileExist==1){
1339  ElecShapesFile = new TFile(elecfile_.c_str());
1340  stringstream name;
1341  name << "MeanElecShape";
1342  elecShape=(TH1D*) ElecShapesFile->Get(name.str().c_str());
1343  if(elecShape && doesMatShapeExist==1){
1344  double x=elecShape->GetMaximum();
1345  if (x!=0) elecShape->Scale(1.0/x);
1346  isSPRFine=true;
1347  }else{
1348  isSPRFine=false;
1349  }
1350 
1351  }else{
1352 
1353  cout <<" ERROR! Elec shape file not found !"<< endl;
1354 
1355  }
1356 
1357 
1358  if(IsMatacqOK){
1359 
1360  ShapeFile = new TFile(shapefile.c_str(),"RECREATE");
1361 
1362  unsigned int nBins=int(laserShape->GetEntries());
1363  assert( nSamplesShapes == nBins);
1364  double elec_jj, laser_iiMinusjj;
1365  double sum_jj;
1366 
1367  if(isSPRFine==true){
1368 
1369  unsigned int nBins2=int(elecShape->GetNbinsX());
1370 
1371  if(nBins2<nBins){
1372  cout<< "EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins"<< std::endl;
1373  isSPRFine=false;
1374  }
1375  assert( nSamplesShapes == nBins2 );
1376 
1377  stringstream name;
1378  name << "PulseShape";
1379 
1380  PulseShape=new TProfile(name.str().c_str(),name.str().c_str(),nBins,-0.5,double(nBins)-0.5);
1381 
1382  // shift shapes to have max close to the real APD max
1383 
1384  for(int ii=0;ii<50;ii++){
1385  shapes[ii]=0.0;
1386  PulseShape->Fill(ii,0.0);
1387  }
1388 
1389  for(unsigned int ii=0;ii<nBins-50;ii++){
1390  sum_jj=0.0;
1391  for(unsigned int jj=0;jj<ii;jj++){
1392  elec_jj=elecShape->GetBinContent(jj+1);
1393  laser_iiMinusjj=laserShape->GetBinContent(ii-jj+1);
1394  sum_jj+=elec_jj*laser_iiMinusjj;
1395  }
1396  PulseShape->Fill(ii+50,sum_jj);
1397  shapes[ii+50]=sum_jj;
1398  }
1399 
1400  double scale= PulseShape->GetMaximum();
1402 
1403  if(scale!=0) {
1404  PulseShape->Scale(1.0/scale);
1405  for(unsigned int ii=0;ii<nBins;ii++){
1406  shapesVec.push_back(shapes[ii]/scale);
1407  }
1408  }
1409 
1410  if(_saveshapes) PulseShape->Write();
1411  }
1412  }
1413  ShapeFile->Close();
1414 
1415  if(!_saveshapes) {
1416 
1417  stringstream del;
1418  del << "rm " <<shapefile;
1419  system(del.str().c_str());
1420 
1421  }
1422 
1423  return IsMatacqOK;
1424 }
1425 
1426 
1427 void EcalLaserAnalyzer2::setGeomEB(int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr){
1428 
1429  side=MEEBGeom::side(etaG,phiG);
1430 
1431  assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
1432 
1433  eta = etaG;
1434  phi = phiG;
1435  channelID=5*(strip-1) + xtal-1;
1436  towerID=tower;
1437 
1438  vector<int> apdRefChan=ME::apdRefChannels(module, lmr);
1439  for (unsigned int iref=0;iref<nRefChan;iref++){
1440  if(channelID==apdRefChan[iref] && towerID==apdRefTT
1441  && apdRefMap[iref].count(module)==0){
1442  apdRefMap[iref][module]=channel;
1443  }
1444  }
1445 
1446  if(isFirstChanModFilled[module-1]==0) {
1447  firstChanMod[module-1]=channel;
1448  isFirstChanModFilled[module-1]=1;
1449  }
1450 
1451  iEta[channel]=eta;
1452  iPhi[channel]=phi;
1453  iModule[channel]= module ;
1454  iTowerID[channel]=towerID;
1455  iChannelID[channel]=channelID;
1456  idccID[channel]=dccID;
1457  iside[channel]=side;
1458 
1459 }
1460 
1461 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){
1462 
1463  side=MEEEGeom::side(iX, iY, iZ);
1464 
1465  assert( module>=*min_element(modules.begin(),modules.end())
1466  && module<=*max_element(modules.begin(),modules.end()) );
1467 
1468  eta = etaG;
1469  phi = phiG;
1470  channelID=ch;
1471  towerID=tower;
1472 
1473  vector<int> apdRefChan=ME::apdRefChannels(module, lmr);
1474  for (unsigned int iref=0;iref<nRefChan;iref++){
1475  if(channelID==apdRefChan[iref] && towerID==apdRefTT
1476  && apdRefMap[iref].count(module)==0){
1477  apdRefMap[iref][module]=channel;
1478  }
1479  }
1480 
1481  if(isFirstChanModFilled[module-1]==0) {
1482  firstChanMod[module-1]=channel;
1483  isFirstChanModFilled[module-1]=1;
1484  }
1485 
1486  iEta[channel]=eta;
1487  iPhi[channel]=phi;
1488  iModule[channel]= module ;
1489  iTowerID[channel]=towerID;
1490  iChannelID[channel]=channelID;
1491  idccID[channel]=dccID;
1492  iside[channel]=side;
1493 
1494 }
1496 
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)
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:356
JetCorrectorParameters corr
Definition: classes.h:9
bool wasTimingOK[NCRYSEB]
std::vector< double > getAPD()
Definition: TAPD.cc:236
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: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:121
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