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