CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalTestPulseAnalyzer.cc
Go to the documentation of this file.
1 /*
2  *
3  * \class EcalTestPulseAnalyzer
4  *
5  * primary author: Julie Malcles - CEA/Saclay
6  * author: Gautier Hamel De Monchenault - CEA/Saclay
7  */
8 #include <TFile.h>
9 #include <TTree.h>
10 
11 using namespace std;
12 #include "EcalTestPulseAnalyzer.h"
13 
14 #include <sstream>
15 #include <iostream>
16 #include <iomanip>
17 
18 
19 
20 #include <FWCore/MessageLogger/interface/MessageLogger.h>
21 #include <FWCore/Utilities/interface/Exception.h>
22 
23 #include <FWCore/Framework/interface/Event.h>
24 #include <FWCore/Framework/interface/MakerMacros.h>
25 #include <FWCore/Framework/interface/ESHandle.h>
26 #include <FWCore/ParameterSet/interface/ParameterSet.h>
27 #include <FWCore/Framework/interface/EventSetup.h>
28 
29 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
30 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
31 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
32 
33 #include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
34 #include <Geometry/EcalMapping/interface/EcalMappingRcd.h>
35 #include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
36 
37 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
38 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TSFit.h>
39 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
40 
41 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h>
42 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h>
43 
44 using namespace std;
45 
46 
47 
48 //========================================================================
49 EcalTestPulseAnalyzer::EcalTestPulseAnalyzer(const edm::ParameterSet& iConfig)
50  //========================================================================
51  :
52 iEvent(0),
53 
54 // framework parameters with default values
55 
56 _nsamples( iConfig.getUntrackedParameter< unsigned int >( "nSamples", 10 ) ),
57 _presample( iConfig.getUntrackedParameter< unsigned int >( "nPresamples", 3 ) ),
58 _firstsample( iConfig.getUntrackedParameter< unsigned int >( "firstSample", 1 ) ),
59 _lastsample( iConfig.getUntrackedParameter< unsigned int >( "lastSample", 2 ) ),
60 _samplemin( iConfig.getUntrackedParameter< unsigned int >( "sampleMin", 3 ) ),
61 _samplemax( iConfig.getUntrackedParameter< unsigned int >( "sampleMax", 9 ) ),
62 _nsamplesPN( iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN", 50 ) ),
63 _presamplePN( iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN", 6 ) ),
64 _firstsamplePN( iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN", 7 ) ),
65 _lastsamplePN( iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN", 8 ) ),
66 _niter( iConfig.getUntrackedParameter< unsigned int >( "nIter", 3 ) ),
67 _chi2max( iConfig.getUntrackedParameter< double >( "chi2Max", 10.0 ) ),
68 _timeofmax( iConfig.getUntrackedParameter< double >( "timeOfMax", 4.5 ) ),
69 _ecalPart( iConfig.getUntrackedParameter< std::string >( "ecalPart", "EB" ) ),
70 _fedid( iConfig.getUntrackedParameter< int >( "fedID", -999 ) ),
71 nCrys( NCRYSEB),
72 nTT( NTTEB),
73 nMod( NMODEB),
74 nGainPN( NGAINPN),
75 nGainAPD( NGAINAPD),
76 towerID(-1), channelID(-1),runType(-1), runNum(0), fedID(-1), dccID(-1), side(-1), iZ(1),
77 phi(-1), eta(-1), event(0), apdAmpl(0),apdTime(0),pnAmpl(0),
78 pnID(-1), moduleID(-1), channelIteratorEE(0)
79 
80 
81  //========================================================================
82 
83 {
84 
85  //now do what ever initialization is needed
86 
87  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
88 
89  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
90  digiPNCollection_ = iConfig.getParameter<std::string>("digiPNCollection");
91  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
92 
93  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
94  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
95 
96 
97  // Define geometrical constants
98 
99  if (_ecalPart == "EB") {
100  nCrys = NCRYSEB;
101  nTT = NTTEB;
102  } else {
103  nCrys = NCRYSEE;
104  nTT = NTTEE;
105  }
106 
107  iZ = 1;
108  if( _fedid <= 609 ) iZ=-1;
109 
110  dccMEM = ME::memFromDcc(_fedid);
111  modules = ME::lmmodFromDcc(_fedid);
112  nMod = modules.size();
113 
114 }
115 
116 
117 //========================================================================
118 EcalTestPulseAnalyzer::~EcalTestPulseAnalyzer(){
119  //========================================================================
120 
121  // do anything here that needs to be done at desctruction time
122  // (e.g. close files, deallocate resources etc.)
123 
124 }
125 
126 
127 
128 //========================================================================
129 void EcalTestPulseAnalyzer::beginJob() {
130  //========================================================================
131 
132  // Define temporary file
133 
134  rootfile=resdir_;
135  rootfile+="/TmpTreeTestPulseAnalyzer.root";
136 
137  outFile = new TFile(rootfile.c_str(),"RECREATE");
138 
139  for (unsigned int i=0;i<nCrys;i++){
140 
141  std::stringstream name;
142  name << "Tree" <<i;
143 
144  trees[i]= new TTree(name.str().c_str(),name.str().c_str());
145 
146  //List of branches
147 
148  trees[i]->Branch( "iphi", &phi, "phi/I" );
149  trees[i]->Branch( "ieta", &eta, "eta/I" );
150  trees[i]->Branch( "side", &side, "side/I" );
151  trees[i]->Branch( "dccID", &dccID, "dccID/I" );
152  trees[i]->Branch( "towerID", &towerID, "towerID/I" );
153  trees[i]->Branch( "channelID", &channelID, "channelID/I" );
154  trees[i]->Branch( "event", &event, "event/I" );
155  trees[i]->Branch( "apdGain", &apdGain, "apdGain/I" );
156  trees[i]->Branch( "pnGain", &pnGain, "pnGain/I" );
157  trees[i]->Branch( "apdAmpl", &apdAmpl, "apdAmpl/D" );
158  trees[i]->Branch( "pnAmpl0", &pnAmpl0, "pnAmpl0/D" );
159  trees[i]->Branch( "pnAmpl1", &pnAmpl1, "pnAmpl1/D" );
160 
161  trees[i]->SetBranchAddress( "ieta", &eta );
162  trees[i]->SetBranchAddress( "iphi", &phi );
163  trees[i]->SetBranchAddress( "side", &side );
164  trees[i]->SetBranchAddress( "dccID", &dccID );
165  trees[i]->SetBranchAddress( "towerID", &towerID );
166  trees[i]->SetBranchAddress( "channelID", &channelID );
167  trees[i]->SetBranchAddress( "event", &event );
168  trees[i]->SetBranchAddress( "apdGain", &apdGain );
169  trees[i]->SetBranchAddress( "pnGain", &pnGain );
170  trees[i]->SetBranchAddress( "apdAmpl", &apdAmpl );
171  trees[i]->SetBranchAddress( "pnAmpl0", &pnAmpl0 );
172  trees[i]->SetBranchAddress( "pnAmpl1", &pnAmpl1 );
173 
174  }
175 
176  // Initializations
177 
178  for(unsigned int j=0;j<nCrys;j++){
179  iEta[j]=-1;
180  iPhi[j]=-1;
181  iModule[j]=10;
182  iTowerID[j]=-1;
183  iChannelID[j]=-1;
184  idccID[j]=-1;
185  iside[j]=-1;
186  }
187 
188  for(unsigned int j=0;j<nMod;j++){
189  firstChanMod[j]=0;
190  isFirstChanModFilled[j]=0;
191  }
192 
193  // Define output results file name
194 
195  std::stringstream namefile;
196  namefile << resdir_ <<"/APDPN_TESTPULSE.root";
197  resfile=namefile.str();
198 
199  // TP events counter
200  TPEvents=0;
201 
202 }
203 
204 
205 //========================================================================
206 void EcalTestPulseAnalyzer:: analyze( const edm::Event & e, const edm::EventSetup& c){
207  //========================================================================
208 
209  ++iEvent;
210 
211  // Retrieve DCC header
212  edm::Handle<EcalRawDataCollection> pDCCHeader;
213  const EcalRawDataCollection* DCCHeader=0;
214  try {
215  e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
216  DCCHeader=pDCCHeader.product();
217  }catch ( std::exception& ex ) {
218  std::cerr << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str() << std::endl;
219  }
220 
221 
222  // retrieving crystal EB data from Event
223  edm::Handle<EBDigiCollection> pEBDigi;
224  const EBDigiCollection* EBDigi=0;
225 
226  // retrieving crystal EE data from Event
227  edm::Handle<EEDigiCollection> pEEDigi;
228  const EEDigiCollection* EEDigi=0;
229 
230  if (_ecalPart == "EB") {
231  try {
232  e.getByLabel(digiProducer_,digiCollection_, pEBDigi);
233  EBDigi=pEBDigi.product();
234  }catch ( std::exception& ex ) {
235  std::cerr << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str() << std::endl;
236  }
237  } else {
238  try {
239  e.getByLabel(digiProducer_,digiCollection_, pEEDigi);
240  EEDigi=pEEDigi.product();
241  }catch ( std::exception& ex ) {
242  std::cerr << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str() << std::endl;
243  }
244  }
245 
246 
247  // Retrieve crystal PN diodes from Event
248  edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
249  const EcalPnDiodeDigiCollection* PNDigi=0;
250  try {
251  e.getByLabel(digiProducer_,digiPNCollection_, pPNDigi);
252  PNDigi=pPNDigi.product();
253  }catch ( std::exception& ex ) {
254  std::cerr << "Error! can't get the product " << digiCollection_.c_str() << std::endl;
255  }
256 
257 
258  // retrieving electronics mapping
259  edm::ESHandle< EcalElectronicsMapping > ecalmapping;
260  const EcalElectronicsMapping* TheMapping=0;
261  try{
262  c.get< EcalMappingRcd >().get(ecalmapping);
263  TheMapping = ecalmapping.product();
264  }catch ( std::exception& ex ) {
265  std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
266  }
267 
268 
269  // ====================================
270  // Decode Basic DCCHeader Information
271  // ====================================
272 
273  for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();headerItr != DCCHeader->end();
274  ++headerItr ) {
275 
276  int fed = headerItr->fedId();
277 
278  if(fed!=_fedid && _fedid!=-999) continue;
279 
280  runType=headerItr->getRunType();
281  runNum=headerItr->getRunNumber();
282  event=headerItr->getLV1();
283 
284  dccID=headerItr->getDccInTCCCommand();
285  fedID=headerItr->fedId();
286 
287 
288  if( 600+dccID != fedID ) continue;
289 
290  // Cut on runType
291 
292  if(runType!=EcalDCCHeaderBlock::TESTPULSE_MGPA && runType!=EcalDCCHeaderBlock::TESTPULSE_GAP
293  && runType!=EcalDCCHeaderBlock::TESTPULSE_SCAN_MEM ) return;
294 
295  }
296 
297  // Cut on fedID
298 
299  if(fedID!=_fedid && _fedid!=-999) return;
300 
301  // Count TP events
302  TPEvents++;
303 
304  // ======================
305  // Decode PN Information
306  // ======================
307 
308  TPNFit * pnfit = new TPNFit();
309  pnfit -> init(_nsamplesPN,_firstsamplePN, _lastsamplePN);
310 
311  double chi2pn=0;
312  double ypnrange[50];
313  double dsum=0.;
314  double dsum1=0.;
315  double bl=0.;
316  double val_max=0.;
317  int samplemax=0;
318  unsigned int k;
319  int pnG[50];
320  int pngain=0;
321 
322  std::map <int, std::vector<double> > allPNAmpl;
323  std::map <int, std::vector<int> > allPNGain;
324 
325  for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr ) { // Loop on PNs
326 
327  EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
328 
329  bool isMemRelevant=false;
330  for (unsigned int imem=0;imem<dccMEM.size();imem++){
331  if(pnDetId.iDCCId() == dccMEM[imem]) {
332  isMemRelevant=true;
333  }
334  }
335 
336  // skip mem dcc without relevant data
337  if(!isMemRelevant) continue;
338 
339  for ( int samId=0; samId < (*pnItr).size() ; samId++ ) { // Loop on PN samples
340  pn[samId]=(*pnItr).sample(samId).adc();
341  pnG[samId]=(*pnItr).sample(samId).gainId();
342 
343  if(pnG[samId]!=1) std::cout << "PN gain different from 1 for sample "<<samId<< std::endl;
344  if (samId==0) pngain=pnG[samId];
345  if (samId>0) pngain=TMath::Max(pnG[samId],pngain);
346  }
347 
348  for(dsum=0.,k=0;k<_presamplePN;k++) {
349  dsum+=pn[k];
350  }
351  bl=dsum/((double) _presamplePN);
352 
353 
354  for(val_max=0.,k=0;k<_nsamplesPN;k++) {
355  ypnrange[k]=pn[k] - bl;
356 
357  if(ypnrange[k] > val_max) {
358  val_max= ypnrange[k]; samplemax=k;
359  }
360  }
361 
362  chi2pn = pnfit -> doFit(samplemax,&ypnrange[0]);
363 
364  if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
365  else pnAmpl= pnfit -> getAmpl();
366 
367  allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
368  allPNGain[pnDetId.iDCCId()].push_back(pngain);
369 
370  }
371 
372  // ===========================
373  // Decode EBDigis Information
374  // ===========================
375 
376  TSFit * pstpfit = new TSFit(_nsamples,650);
377  pstpfit -> set_params(_nsamples, _niter, _presample, _samplemin, _samplemax, _timeofmax , _chi2max, _firstsample, _lastsample);
378  pstpfit -> init_errmat(10.);
379 
380  double chi2=0;
381  double yrange[10];
382  int adcgain=0;
383  int adcG[10];
384 
385 
386  if (EBDigi){
387  for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr ) { // Loop on EB crystals
388 
389  EBDetId id_crystal(digiItr->id()) ;
390  EBDataFrame df( *digiItr );
391 
392  int etaG = id_crystal.ieta() ; // global
393  int phiG = id_crystal.iphi() ; // global
394 
395  int etaL ; // local
396  int phiL ; // local
397  std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
398 
399  etaL=LocalCoord.first ;
400  phiL=LocalCoord.second ;
401 
402  eta = etaG;
403  phi = phiG;
404 
405  side=MEEBGeom::side(etaG,phiG);
406 
407  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
408 
409  towerID=elecid_crystal.towerId();
410  int strip=elecid_crystal.stripId();
411  int xtal=elecid_crystal.xtalId();
412  channelID= 5*(strip-1) + xtal-1; // FIXME
413 
414 
415  int module= MEEBGeom::lmmod(etaG, phiG);
416  int iMod = module-1;
417 
418  assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
419 
420 
421  std::pair<int,int> pnpair=MEEBGeom::pn(module);
422  unsigned int MyPn0=pnpair.first;
423  unsigned int MyPn1=pnpair.second;
424 
425  unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
426 
427  if(isFirstChanModFilled[iMod]==0) {
428  firstChanMod[iMod]=channel;
429  isFirstChanModFilled[iMod]=1;
430  }
431 
432  iEta[channel]=eta;
433  iPhi[channel]=phi;
434  iModule[channel]= module ;
435  iTowerID[channel]=towerID;
436  iChannelID[channel]=channelID;
437  idccID[channel]=dccID;
438  iside[channel]=side;
439 
440 
441 
442  // get adc samples
443  //====================
444 
445  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {
446 
447  EcalMGPASample samp_crystal(df.sample(i));
448  adc[i]=samp_crystal.adc() ;
449  adcG[i]=samp_crystal.gainId();
450 
451  if (i==0) adcgain=adcG[i];
452  if (i>0) adcgain=TMath::Max(adcG[i],adcgain);
453  }
454  // Remove pedestal
455  //====================
456  for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
457  dsum+=adc[k];
458  if(k<_presample-1) dsum1+=adc[k];
459  }
460 
461  bl=dsum/((double)_presample);
462 
463  for(val_max=0.,k=0;k<_nsamples;k++) {
464  yrange[k]=adc[k] - bl;
465  if(yrange[k] > val_max) {
466  val_max= yrange[k]; samplemax=k;
467  }
468  }
469 
470  apdGain=adcgain;
471 
472  if(allPNAmpl[dccMEM[0]].size()>MyPn0) pnAmpl0=allPNAmpl[dccMEM[0]][MyPn0];
473  else pnAmpl0=0;
474  if(allPNAmpl[dccMEM[0]].size()>MyPn1) pnAmpl1=allPNAmpl[dccMEM[0]][MyPn1];
475  else pnAmpl1=0;
476 
477  if(allPNGain[dccMEM[0]].size()>MyPn0) pnGain=allPNGain[dccMEM[0]][MyPn0];
478  else pnGain=0;
479 
480  // Perform the fit on apd samples
481  //================================
482 
483  chi2 = pstpfit -> fit_third_degree_polynomial(&yrange[0],ret_data);
484 
485  //Retrieve APD amplitude from fit
486  //================================
487 
488  if( val_max > 100000. || chi2 < 0. || chi2 == 102 ) {
489 
490  apdAmpl=0;
491  apdTime=0;
492 
493  }else{
494 
495  apdAmpl = ret_data[0];
496  apdTime = ret_data[1];
497 
498  }
499 
500  trees[channel]->Fill();
501 
502  }
503 
504  } else {
505 
506  for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr ) { // Loop on EE crystals
507 
508  EEDetId id_crystal(digiItr->id()) ;
509  EEDataFrame df( *digiItr );
510 
511  phi = id_crystal.ix() ;
512  eta = id_crystal.iy() ;
513 
514  int iX = (phi-1)/5+1;
515  int iY = (eta-1)/5+1;
516 
517  side=MEEEGeom::side( iX, iY ,iZ);
518 
519 
520  // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
521 
522  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
523 
524 
525  towerID=elecid_crystal.towerId();
526  channelID=elecid_crystal.channelId()-1;
527 
528  int module=MEEEGeom::lmmod( iX, iY );
529  if( module>=18 && side==1 ) module+=2; // Trick to fix endcap specificity
530  int iMod = module-1;
531 
532  assert( module>=*min_element(modules.begin(),modules.end()) && module<=*max_element(modules.begin(),modules.end()) );
533 
534  std::pair<int,int> pnpair=MEEEGeom::pn( module, _fedid ) ;
535 
536  unsigned int MyPn0=pnpair.first;
537  unsigned int MyPn1=pnpair.second;
538 
539  int hashedIndex=100000*eta+phi;
540 
541  if( channelMapEE.count(hashedIndex) == 0 ){
542  channelMapEE[hashedIndex]=channelIteratorEE;
543  channelIteratorEE++;
544  }
545 
546  unsigned int channel=channelMapEE[hashedIndex];
547 
548  if(isFirstChanModFilled[iMod]==0) {
549  firstChanMod[iMod]=channel;
550  isFirstChanModFilled[iMod]=1;
551  }
552 
553  iEta[channel]=eta;
554  iPhi[channel]=phi;
555  iModule[channel]= module ;
556  iTowerID[channel]=towerID;
557  iChannelID[channel]=channelID;
558  idccID[channel]=dccID;
559  iside[channel]=side;
560 
561  assert (channel < nCrys);
562 
563  // Get adc samples
564  //====================
565 
566  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) {
567 
568  EcalMGPASample samp_crystal(df.sample(i));
569  adc[i]=samp_crystal.adc() ;
570  adcG[i]=samp_crystal.gainId();
571 
572  if (i==0) adcgain=adcG[i];
573  if (i>0) adcgain=TMath::Max(adcG[i],adcgain);
574  }
575 
576 
577  // Remove pedestal
578  //====================
579  for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
580  dsum+=adc[k];
581  if(k<_presample-1) dsum1+=adc[k];
582  }
583 
584  bl=dsum/((double)_presample);
585 
586  for(val_max=0.,k=0;k<_nsamples;k++) {
587  yrange[k]=adc[k] - bl;
588  if(yrange[k] > val_max) {
589  val_max= yrange[k]; samplemax=k;
590  }
591  }
592  apdGain=adcgain;
593 
594  int dccMEMIndex=0;
595  if(side==1) dccMEMIndex+=2; // Trick to fix endcap specificity
596 
597  if(allPNAmpl[dccMEM[dccMEMIndex]].size()>MyPn0) pnAmpl0=allPNAmpl[dccMEM[dccMEMIndex]][MyPn0];
598  else pnAmpl0=0;
599  if(allPNAmpl[dccMEM[dccMEMIndex+1]].size()>MyPn1) pnAmpl1=allPNAmpl[dccMEM[dccMEMIndex+1]][MyPn1];
600  else pnAmpl1=0;
601 
602  if(allPNGain[dccMEM[dccMEMIndex]].size()>MyPn0) pnGain=allPNGain[dccMEM[dccMEMIndex]][MyPn0];
603  else pnGain=0;
604 
605 
606  // Perform the fit on apd samples
607  //=================================
608 
609  chi2 = pstpfit -> fit_third_degree_polynomial(&yrange[0],ret_data);
610 
611  //Retrieve APD amplitude from fit
612  //=================================
613 
614  if( val_max > 100000. || chi2 < 0. || chi2 == 102 ) {
615 
616  apdAmpl=0;
617  apdTime=0;
618 
619  }else{
620 
621  apdAmpl = ret_data[0];
622  apdTime = ret_data[1];
623 
624  }
625 
626  trees[channel]->Fill();
627  }
628 
629  }
630 
631 
632 } // end of analyze
633 
634 
635 //========================================================================
636 void EcalTestPulseAnalyzer::endJob() {
637  //========================================================================
638 
639  // Don't do anything if there is no events
640  if( TPEvents == 0 ){
641 
642  outFile->Close();
643 
644  // Remove temporary file
645 
646  std::stringstream del;
647  del << "rm " <<rootfile;
648  system(del.str().c_str());
649 
650  std::cout << " No TP Events "<< std::endl;
651  return;
652  }
653 
654  std::cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
655  std::cout << "\t+=+ Analyzing test pulse data: getting APD, PN +=+" << std::endl;
656 
657 
658  // Create output ntuples:
659 
660  //std::cout<< "TP Test Name File "<< resfile.c_str() << std::endl;
661 
662  resFile = new TFile(resfile.c_str(),"RECREATE");
663 
664  restrees= new TTree("TPAPD","TPAPD");
665  respntrees= new TTree("TPPN","TPPN");
666 
667  restrees->Branch( "iphi", &iphi, "iphi/I" );
668  restrees->Branch( "ieta", &ieta, "ieta/I" );
669  restrees->Branch( "dccID", &dccID, "dccID/I" );
670  restrees->Branch( "side", &side, "side/I" );
671  restrees->Branch( "towerID", &towerID, "towerID/I" );
672  restrees->Branch( "channelID", &channelID, "channelID/I" );
673  restrees->Branch( "moduleID", &moduleID, "moduleID/I" );
674  restrees->Branch( "flag", &flag, "flag/I" );
675  restrees->Branch( "gain", &gain, "gain/I" );
676  restrees->Branch( "APD", &APD, "APD[6]/D" );
677 
678  respntrees->Branch( "pnID", &pnID, "pnID/I" );
679  respntrees->Branch( "moduleID", &moduleID, "moduleID/I" );
680  respntrees->Branch( "gain", &gain, "gain/I" );
681  respntrees->Branch( "PN", &PN, "PN[6]/D" );
682 
683  restrees->SetBranchAddress( "iphi", &iphi );
684  restrees->SetBranchAddress( "ieta", &ieta );
685  restrees->SetBranchAddress( "dccID", &dccID );
686  restrees->SetBranchAddress( "side", &side );
687  restrees->SetBranchAddress( "towerID", &towerID );
688  restrees->SetBranchAddress( "channelID", &channelID );
689  restrees->SetBranchAddress( "moduleID", &moduleID );
690  restrees->SetBranchAddress( "flag", &flag );
691  restrees->SetBranchAddress( "gain", &gain );
692  restrees->SetBranchAddress( "APD", APD );
693 
694  respntrees->SetBranchAddress( "pnID", &pnID );
695  respntrees->SetBranchAddress( "moduleID", &moduleID );
696  respntrees->SetBranchAddress( "gain", &gain );
697  respntrees->SetBranchAddress( "PN", PN );
698 
699 
700 
701  TMom *APDAnal[1700][10];
702  TMom *PNAnal[9][2][10];
703 
704 
705  for (unsigned int iMod=0;iMod<nMod;iMod++){
706  for (unsigned int ich=0;ich<2;ich++){
707  for (unsigned int ig=0;ig<nGainPN;ig++){
708  PNAnal[iMod][ich][ig]=new TMom();
709  }
710  }
711  }
712 
713  for (unsigned int iCry=0;iCry<nCrys;iCry++){ // Loop on data trees (ie on cristals)
714 
715  for(unsigned int iG=0;iG<nGainAPD;iG++){
716  APDAnal[iCry][iG]=new TMom();
717  }
718 
719 
720 
721  // Define submodule and channel number inside the submodule (as Patrice)
722 
723  unsigned int iMod=iModule[iCry]-1;
724 
725  moduleID=iMod+1;
726  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
727 
728  Long64_t nbytes = 0, nb = 0;
729  for (Long64_t jentry=0; jentry< trees[iCry]->GetEntriesFast();jentry++) {
730  nb = trees[iCry]->GetEntry(jentry); nbytes += nb;
731 
732  // PN Means and RMS
733 
734  if( firstChanMod[iMod]==iCry ){
735  PNAnal[iMod][0][pnGain]->addEntry(pnAmpl0);
736  PNAnal[iMod][1][pnGain]->addEntry(pnAmpl1);
737  }
738 
739  // APD means and RMS
740 
741  APDAnal[iCry][apdGain]->addEntry(apdAmpl);
742 
743  }
744 
745  if (trees[iCry]->GetEntries()<10){
746  flag=-1;
747  for (int j=0;j<6;j++){
748  APD[j]=0.0;
749  }
750  }
751  else flag=1;
752 
753  iphi=iPhi[iCry];
754  ieta=iEta[iCry];
755  dccID=idccID[iCry];
756  side=iside[iCry];
757  towerID=iTowerID[iCry];
758  channelID=iChannelID[iCry];
759 
760  for (unsigned int ig=0;ig<nGainAPD;ig++){
761 
762  APD[0]= APDAnal[iCry][ig]->getMean();
763  APD[1]= APDAnal[iCry][ig]->getRMS();
764  APD[2]= APDAnal[iCry][ig]->getM3();
765  APD[3]= APDAnal[iCry][ig]->getNevt();
766  APD[4]= APDAnal[iCry][ig]->getMin();
767  APD[5]= APDAnal[iCry][ig]->getMax();
768  gain=ig;
769 
770  // Fill APD tree
771 
772  restrees->Fill();
773 
774  }
775  }
776 
777  // Get final results for PN and PN/PN
778 
779  for (unsigned int ig=0;ig<nGainPN;ig++){
780  for (unsigned int iMod=0;iMod<nMod;iMod++){
781  for (int ch=0;ch<2;ch++){
782 
783  pnID=ch;
784  moduleID=iMod;
785  if( moduleID>=20 ) moduleID-=2; // Trick to fix endcap specificity
786 
787  PN[0]= PNAnal[iMod][ch][ig]->getMean();
788  PN[1]= PNAnal[iMod][ch][ig]->getRMS();
789  PN[2]= PNAnal[iMod][ch][ig]->getM3();
790  PN[3]= PNAnal[iMod][ch][ig]->getNevt();
791  PN[4]= PNAnal[iMod][ch][ig]->getMin();
792  PN[5]= PNAnal[iMod][ch][ig]->getMax();
793  gain=ig;
794 
795  // Fill PN tree
796  respntrees->Fill();
797 
798  }
799  }
800  }
801 
802  outFile->Close();
803 
804  // Remove temporary file
805 
806  std::stringstream del;
807  del << "rm " <<rootfile;
808  system(del.str().c_str());
809 
810 
811  // Save final results
812 
813  restrees->Write();
814  respntrees->Write();
815  resFile->Close();
816 
817  std::cout << "\t+=+ ...................................... done +=+" << std::endl;
818  std::cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << std::endl;
819 }
820 
821 
822 DEFINE_FWK_MODULE(EcalTestPulseAnalyzer);
823 
TMom::addEntry
void addEntry(double val)
Definition: TMom.cc:110
MEEBGeom::localCoord
static XYCoord localCoord(int icr)
Definition: MEEBGeom.cc:153
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
EcalTestPulseAnalyzer::channelMapEE
std::map< int, int > channelMapEE
Definition: EcalTestPulseAnalyzer.h:145
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
MEEEGeom::lmmod
static int lmmod(SuperCrysCoord iX, SuperCrysCoord iY)
Definition: MEEEGeom.cc:94
i
int i
Definition: DBlmapReader.cc:9
ME::memFromDcc
static std::vector< ME::DCCid > memFromDcc(ME::DCCid idcc)
Definition: ME.cc:608
NTTEB
#define NTTEB
Definition: EcalPerEvtLaserAnalyzer.h:25
EcalElectronicsId::xtalId
int xtalId() const
get the channel id
Definition: EcalElectronicsId.h:41
edm::DataFrameContainer::const_iterator
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
Definition: DataFrameContainer.h:51
EcalTestPulseAnalyzer::iModule
unsigned int iModule[1700]
Definition: EcalTestPulseAnalyzer.h:158
EcalTestPulseAnalyzer::modules
std::vector< int > modules
Definition: EcalTestPulseAnalyzer.h:147
EcalElectronicsId::stripId
int stripId() const
get the tower id
Definition: EcalElectronicsId.h:39
NGAINPN
#define NGAINPN
Definition: EcalTestPulseAnalyzer.h:32
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
EcalElectronicsId
Ecal readout channel identification [32:20] Unused (so far) [19:13] DCC id [12:6] tower [5:3] strip [...
Definition: EcalElectronicsId.h:20
NTTEE
#define NTTEE
Definition: EcalPerEvtLaserAnalyzer.h:30
NCRYSEE
#define NCRYSEE
Definition: EcalABAnalyzer.h:33
MEEEGeom::side
static int side(SuperCrysCoord iX, SuperCrysCoord iY, int iz)
Definition: MEEEGeom.cc:939
assert
assert(m_qm.get())
init
int init
Definition: HydjetWrapper.h:62
edm::SortedCollection< EcalDCCHeaderBlock >::const_iterator
std::vector< EcalDCCHeaderBlock >::const_iterator const_iterator
Definition: SortedCollection.h:82
EcalElectronicsId::towerId
int towerId() const
get the tower id
Definition: EcalElectronicsId.h:37
edm::DataFrameContainer::begin
const_iterator begin() const
Definition: DataFrameContainer.h:152
NGAINAPD
#define NGAINAPD
Definition: EcalTestPulseAnalyzer.h:33
TMom
Definition: TMom.h:7
EcalTestPulseAnalyzer::EcalTestPulseAnalyzer
EcalTestPulseAnalyzer(const edm::ParameterSet &iConfig)
Definition: EcalTestPulseAnalyzer.cc:49
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
ecal_dqm_sourceclient-live_cfg.cerr
tuple cerr
Definition: ecal_dqm_sourceclient-live_cfg.py:90
MEEBGeom::pn
static std::pair< int, int > pn(int ilmmod)
Definition: MEEBGeom.cc:479
TMom::getM3
double getM3()
Definition: TMom.cc:174
fileinputsource_cfi.runType
tuple runType
Definition: fileinputsource_cfi.py:34
nTT
#define nTT
Definition: TMEGeom.h:6
TMom::getNevt
int getNevt()
Definition: TMom.cc:164
ecalpyutils::hashedIndex
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:42
EcalElectronicsMapping::getElectronicsId
EcalElectronicsId getElectronicsId(const DetId &id) const
Get the electronics id for this det id.
Definition: EcalElectronicsMapping.cc:333
MEEEGeom::pn
static std::pair< int, int > pn(int dee, int ilmod)
Definition: MEEEGeom.cc:615
iEvent
int iEvent
Definition: GenABIO.cc:230
MEEBGeom::electronic_channel
static int electronic_channel(EBLocalCoord ix, EBLocalCoord iy)
Definition: MEEBGeom.cc:341
NCRYSEB
#define NCRYSEB
Definition: EcalABAnalyzer.h:30
NMODEB
#define NMODEB
Definition: EcalLaserAnalyzer.h:37
j
int j
Definition: DBlmapReader.cc:9
event
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
EcalPnDiodeDetId::iDCCId
int iDCCId() const
get the DCCId
Definition: EcalPnDiodeDetId.h:46
MEEBGeom::lmmod
static int lmmod(EBGlobalCoord ieta, EBGlobalCoord iphi)
Definition: MEEBGeom.cc:93
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:413
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:310
Max
T Max(T a, T b)
Definition: MathUtil.h:44
MEEBGeom::side
static int side(EBGlobalCoord ieta, EBGlobalCoord iphi)
Definition: MEEBGeom.cc:114
TMom::getRMS
double getRMS()
Definition: TMom.cc:166
edm::Handle::product
T const * product() const
Definition: Handle.h:81
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
TMom::getMin
double getMin()
Definition: TMom.cc:186
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
phi
Geom::Phi< T > phi() const
Definition: Basic3DVectorLD.h:163
edm::DataFrameContainer::end
const_iterator end() const
Definition: DataFrameContainer.h:156
TMom::getMax
double getMax()
Definition: TMom.cc:187
TPNFit
Definition: TPNFit.h:8
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
ecaldqm::towerID
EcalLogicID towerID(EcalElectronicsId const &)
Definition: LogicIDTranslation.cc:48
EcalTestPulseAnalyzer::analyze
virtual void analyze(const edm::Event &e, const edm::EventSetup &c)
Definition: EcalTestPulseAnalyzer.cc:206
TSFit
Definition: TSFit.h:16
relval_steps.k
list k
Definition: relval_steps.py:1462
EcalTestPulseAnalyzer::isFirstChanModFilled
unsigned int isFirstChanModFilled[9]
Definition: EcalTestPulseAnalyzer.h:162
beam_dqm_sourceclient-live_cfg.chi2
tuple chi2
Definition: beam_dqm_sourceclient-live_cfg.py:335
module
Definition: vlib.h:208
ME::lmmodFromDcc
static std::vector< ME::LMMid > lmmodFromDcc(ME::DCCid idcc)
Definition: ME.cc:623
TMom::getMean
double getMean()
Definition: TMom.cc:147
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
EcalElectronicsId::channelId
int channelId() const
so far for EndCap only :
Definition: EcalElectronicsId.cc:70
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:303
Generated for CMSSW Reference Manual by doxygen 1.8.5