CMS 3D CMS Logo

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