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