CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalPerEvtLaserAnalyzer.cc
Go to the documentation of this file.
1 /*
2  * \class EcalPerEvtLaserAnalyzer
3  *
4  * primary author: Julie Malcles - CEA/Saclay
5  * author: Gautier Hamel De Monchenault - CEA/Saclay
6  */
7 
8 #include <TFile.h>
9 #include <TTree.h>
10 
11 #include "EcalPerEvtLaserAnalyzer.h"
12 
13 #include <sstream>
14 #include <iostream>
15 #include <iomanip>
16 #include <ctime>
17 
18 #include <FWCore/MessageLogger/interface/MessageLogger.h>
19 #include <FWCore/Utilities/interface/Exception.h>
20 
21 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEEEGeom.h>
22 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/MEEBGeom.h>
23 
24 #include <FWCore/Framework/interface/Event.h>
25 #include <FWCore/Framework/interface/MakerMacros.h>
26 #include <FWCore/ParameterSet/interface/ParameterSet.h>
27 #include <FWCore/Framework/interface/EventSetup.h>
28 #include <FWCore/Framework/interface/ESHandle.h>
29 
30 #include <Geometry/EcalMapping/interface/EcalElectronicsMapping.h>
31 #include <Geometry/EcalMapping/interface/EcalMappingRcd.h>
32 
33 #include <DataFormats/EcalDigi/interface/EcalDigiCollections.h>
34 #include <DataFormats/EcalDetId/interface/EcalElectronicsId.h>
35 #include <DataFormats/EcalDetId/interface/EcalDetIdCollections.h>
36 #include <DataFormats/EcalRawData/interface/EcalRawDataCollections.h>
37 
38 #include <vector>
39 
40 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h>
41 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h>
42 
43 using namespace std;
44 
45 
46 //========================================================================
47 EcalPerEvtLaserAnalyzer::EcalPerEvtLaserAnalyzer(const edm::ParameterSet& iConfig)
48  //========================================================================
49  :
50 iEvent(0),
51 
52 // framework parameters with default values
53 _nsamples( iConfig.getUntrackedParameter< unsigned int >( "nSamples", 10 ) ),
54 _presample( iConfig.getUntrackedParameter< unsigned int >( "nPresamples", 3 ) ),
55 _firstsample( iConfig.getUntrackedParameter< unsigned int >( "firstSample", 1 ) ),
56 _lastsample( iConfig.getUntrackedParameter< unsigned int >( "lastSample", 2 ) ),
57 _samplemin( iConfig.getUntrackedParameter< unsigned int >( "sampleMin", 3 ) ),
58 _samplemax( iConfig.getUntrackedParameter< unsigned int >( "sampleMax", 9 ) ),
59 _nsamplesPN( iConfig.getUntrackedParameter< unsigned int >( "nSamplesPN", 50 ) ),
60 _presamplePN( iConfig.getUntrackedParameter< unsigned int >( "nPresamplesPN", 6 ) ),
61 _firstsamplePN( iConfig.getUntrackedParameter< unsigned int >( "firstSamplePN", 7 ) ),
62 _lastsamplePN( iConfig.getUntrackedParameter< unsigned int >( "lastSamplePN", 8 ) ),
63 _timingcutlow( iConfig.getUntrackedParameter< unsigned int >( "timingCutLow", 3 ) ),
64 _timingcuthigh( iConfig.getUntrackedParameter< unsigned int >( "timingCutHigh", 7 ) ),
65 _niter( iConfig.getUntrackedParameter< unsigned int >( "nIter", 3 ) ),
66 _fedid( iConfig.getUntrackedParameter< unsigned int >( "fedID", 0 ) ),
67 _tower( iConfig.getUntrackedParameter< unsigned int >( "tower", 1 ) ),
68 _channel( iConfig.getUntrackedParameter< unsigned int >( "channel", 1 ) ),
69 _ecalPart( iConfig.getUntrackedParameter< std::string >( "ecalPart", "EB" ) ),
70 nCrys( NCRYSEB),
71 nTT( NTTEB),
72 nSides( NSIDES),
73 IsFileCreated(0), nCh(0),
74 runType(-1), runNum(0), dccID(-1), lightside(2), doesRefFileExist(0),
75 ttMat(-1), peakMat(-1), peak(-1), evtMat(-1), colMat(-1)
76  //========================================================================
77 
78 {
79 
80  //now do what ever initialization is needed
81 
82  resdir_ = iConfig.getUntrackedParameter<std::string>("resDir");
83  refalphabeta_ = iConfig.getUntrackedParameter<std::string>("refAlphaBeta");
84 
85  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
86  digiPNCollection_ = iConfig.getParameter<std::string>("digiPNCollection");
87  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
88 
89  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
90  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
91 
92  // Define geometrical constants
93 
94  if (_ecalPart == "EB") {
95  nCrys = NCRYSEB;
96  nTT = NTTEB;
97  } else {
98  nCrys = NCRYSEE;
99  nTT = NTTEE;
100  }
101 
102 }
103 
104 //========================================================================
105 EcalPerEvtLaserAnalyzer::~EcalPerEvtLaserAnalyzer(){
106  //========================================================================
107 
108 
109  // do anything here that needs to be done at desctruction time
110  // (e.g. close files, deallocate resources etc.)
111 
112 }
113 
114 
115 
116 //========================================================================
117 void EcalPerEvtLaserAnalyzer::beginJob() {
118  //========================================================================
119 
120 
121  // Define temporary files' names
122 
123  stringstream namefile1;
124  namefile1 <<resdir_<< "/ADCSamples.root" ;
125 
126  ADCfile=namefile1.str();
127 
128  // Create temporary file and trees to save adc samples
129 
130  ADCFile = new TFile(ADCfile.c_str(),"RECREATE");
131 
132 
133  stringstream name;
134  name <<"ADCTree";
135 
136  ADCtrees= new TTree(name.str().c_str(),name.str().c_str());
137 
138  ADCtrees->Branch( "iphi", &phi, "phi/I" );
139  ADCtrees->Branch( "ieta", &eta, "eta/I" );
140  ADCtrees->Branch( "dccID", &dccID, "dccID/I" );
141  ADCtrees->Branch( "event", &event, "event/I" );
142  ADCtrees->Branch( "color", &color, "color/I" );
143  ADCtrees->Branch( "adc", &adc , "adc[10]/D" );
144  ADCtrees->Branch( "ttrigMatacq", &ttrig, "ttrig/D" );
145  ADCtrees->Branch( "peakMatacq", &peak, "peak/D" );
146  ADCtrees->Branch( "pn0", &pn0, "pn0/D" );
147  ADCtrees->Branch( "pn1", &pn1, "pn1/D" );
148 
149  ADCtrees->SetBranchAddress( "ieta", &eta );
150  ADCtrees->SetBranchAddress( "iphi", &phi );
151  ADCtrees->SetBranchAddress( "dccID", &dccID );
152  ADCtrees->SetBranchAddress( "event", &event );
153  ADCtrees->SetBranchAddress( "color", &color );
154  ADCtrees->SetBranchAddress( "adc", adc );
155  ADCtrees->SetBranchAddress( "ttrigMatacq", &ttrig );
156  ADCtrees->SetBranchAddress( "peakMatacq", &peak );
157  ADCtrees->SetBranchAddress( "pn0", &pn0 );
158  ADCtrees->SetBranchAddress( "pn1", &pn1 );
159 
160  IsFileCreated=0;
161  nCh=nCrys;
162 
163 }
164 
165 
166 //========================================================================
167 void EcalPerEvtLaserAnalyzer:: analyze( const edm::Event & e, const edm::EventSetup& c){
168  //========================================================================
169 
170  ++iEvent;
171 
172  // retrieving DCC header
173  edm::Handle<EcalRawDataCollection> pDCCHeader;
174  const EcalRawDataCollection* DCCHeader=0;
175  try {
176  e.getByLabel(eventHeaderProducer_,eventHeaderCollection_, pDCCHeader);
177  DCCHeader=pDCCHeader.product();
178  }catch ( std::exception& ex ) {
179  std::cerr << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str() << std::endl;
180  }
181 
182  // retrieving crystal data from Event
183  edm::Handle<EBDigiCollection> pEBDigi;
184  const EBDigiCollection* EBDigi=0;
185  edm::Handle<EEDigiCollection> pEEDigi;
186  const EEDigiCollection* EEDigi=0;
187 
188  if (_ecalPart == "EB") {
189  try {
190  e.getByLabel(digiProducer_,digiCollection_, pEBDigi);
191  EBDigi=pEBDigi.product();
192  }catch ( std::exception& ex ) {
193  std::cerr << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str() << std::endl;
194  }
195  } else {
196  try {
197  e.getByLabel(digiProducer_,digiCollection_, pEEDigi);
198  EEDigi=pEEDigi.product();
199  }catch ( std::exception& ex ) {
200  std::cerr << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str() << std::endl;
201  }
202  }
203 
204 
205  // retrieving crystal PN diodes from Event
206  edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
207  const EcalPnDiodeDigiCollection* PNDigi=0;
208  try {
209  e.getByLabel(digiProducer_,digiPNCollection_, pPNDigi);
210  PNDigi=pPNDigi.product();
211  }catch ( std::exception& ex ) {
212  std::cerr << "Error! can't get the product " << digiPNCollection_.c_str() << std::endl;
213  }
214 
215  // retrieving electronics mapping
216  edm::ESHandle< EcalElectronicsMapping > ecalmapping;
217  const EcalElectronicsMapping* TheMapping=0;
218  try{
219  c.get< EcalMappingRcd >().get(ecalmapping);
220  TheMapping = ecalmapping.product();
221  }catch ( std::exception& ex ) {
222  std::cerr << "Error! can't get the product EcalMappingRcd"<< std::endl;
223  }
224 
225  // ====================================
226  // Decode Basic DCCHeader Information
227  // ====================================
228 
229  for ( EcalRawDataCollection::const_iterator headerItr= DCCHeader->begin();headerItr != DCCHeader->end();
230  ++headerItr ) {
231 
232  // Get run type and run number
233 
234  int fed = headerItr->fedId();
235 
236  if(fed!=_fedid && _fedid!=-999) continue;
237 
238  runType=headerItr->getRunType();
239  runNum=headerItr->getRunNumber();
240  event=headerItr->getLV1();
241  dccID=headerItr->getDccInTCCCommand();
242  fedID=headerItr->fedId();
243 
244  // take event only if the fed corresponds to the DCC in TCC
245  if( 600+dccID != fedID ) continue;
246 
247  // Cut on runType
248  if(runType!=EcalDCCHeaderBlock::LASER_STD && runType!=EcalDCCHeaderBlock::LASER_GAP) return;
249 
250  // Define output results files' names
251 
252  if (IsFileCreated==0){
253 
254  stringstream namefile2;
255  namefile2 << resdir_ <<"/APDAmpl_Run"<<runNum<<"_"<<_fedid <<"_"<<_tower<<"_"<<_channel<<".root";
256  resfile=namefile2.str();
257 
258  // Get Matacq ttrig
259 
260  stringstream namefile;
261  namefile << resdir_ <<"/Matacq-Run"<<runNum<<".root";
262 
263  doesRefFileExist=0;
264 
265  FILE *test;
266  test = fopen(namefile.str().c_str(),"r");
267  if (test) doesRefFileExist=1;
268 
269  if(doesRefFileExist==1){
270  matacqFile = new TFile((namefile.str().c_str()));
271  matacqTree=(TTree*)matacqFile->Get("MatacqShape");
272 
273  matacqTree->SetBranchAddress( "event", &evtMat );
274  matacqTree->SetBranchAddress( "color", &colMat );
275  matacqTree->SetBranchAddress( "peak", &peakMat );
276  matacqTree->SetBranchAddress( "ttrig", &ttMat );
277  }
278 
279  IsFileCreated=1;
280 
281  }
282 
283  // Retrieve laser color and event number
284 
285  EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
286  int color = settings.wavelength;
287 
288  vector<int>::iterator iter = find( colors.begin(), colors.end(), color );
289  if( iter==colors.end() ){
290  colors.push_back( color );
291  cout <<" new color found "<< color<<" "<< colors.size()<< endl;
292  }
293 
294  }
295 
296  // cut on fedID
297 
298  if(fedID!=_fedid && _fedid!=-999) return;
299 
300 
301  // ======================
302  // Decode PN Information
303  // ======================
304 
305  TPNFit * pnfit = new TPNFit();
306  pnfit -> init(_nsamplesPN,_firstsamplePN, _lastsamplePN);
307 
308  double chi2pn=0;
309  double ypnrange[50];
310  double dsum=0.;
311  double dsum1=0.;
312  double bl=0.;
313  double bl1=0.;
314  double val_max=0.;
315  unsigned int samplemax=0;
316  unsigned int k;
317 
318  std::vector<double> allPNAmpl;
319 
320  for ( EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr ) { // Loop on PNs
321 
322 
323  for ( int samId=0; samId < (*pnItr).size() ; samId++ ) { // Loop on PN samples
324  pn[samId]=(*pnItr).sample(samId).adc();
325  }
326 
327 
328  for(dsum=0.,k=0;k<_presamplePN;k++) {
329  dsum+=pn[k];
330  }
331  bl=dsum/((double) _presamplePN);
332 
333  for(val_max=0.,k=0;k<_nsamplesPN;k++) {
334  ypnrange[k]=pn[k] - bl;
335 
336  if(ypnrange[k] > val_max) {
337  val_max= ypnrange[k]; samplemax=k;
338  }
339  }
340 
341  chi2pn = pnfit -> doFit(samplemax,&ypnrange[0]);
342 
343  if(chi2pn == 101 || chi2pn == 102 || chi2pn == 103) pnAmpl=0.;
344  else pnAmpl= pnfit -> getAmpl();
345 
346  allPNAmpl.push_back(pnAmpl);
347 
348  }
349 
350  // ===========
351  // Get Matacq
352  // ===========
353 
354  ttrig=-1.;
355  peak=-1.;
356 
357  if(doesRefFileExist==1){
358  // FIXME
359  if (color==0) matacqTree->GetEntry(event-1);
360  else if(color==3) matacqTree->GetEntry(matacqTree->GetEntries("color==0")+event-1);
361  ttrig=ttMat;
362  peak=peakMat;
363 
364  }
365 
366 
367  // ===========================
368  // Decode EBDigis Information
369  // ===========================
370 
371  double yrange[10];
372  int adcGain=0;
373  int side=0;
374 
375  if (EBDigi){
376 
377  for ( EBDigiCollection::const_iterator digiItr= EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr ) { // Loop on crystals
378 
379  EBDetId id_crystal(digiItr->id()) ;
380  EBDataFrame df( *digiItr );
381 
382  int etaG = id_crystal.ieta() ; // global
383  int phiG = id_crystal.iphi() ; // global
384 
385  int etaL ; // local
386  int phiL ; // local
387  std::pair<int, int> LocalCoord=MEEBGeom::localCoord( etaG , phiG );
388 
389  etaL=LocalCoord.first ;
390  phiL=LocalCoord.second ;
391 
392  eta = etaG;
393  phi = phiG;
394 
395  side=MEEBGeom::side(etaG,phiG);
396 
397  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
398 
399  int towerID=elecid_crystal.towerId();
400  // int channelID=elecid_crystal.channelId()-1; // FIXME so far for endcap only
401  int strip=elecid_crystal.stripId();
402  int xtal=elecid_crystal.xtalId();
403  int channelID= 5*(strip-1) + xtal-1; // FIXME
404 
405  int module= MEEBGeom::lmmod(etaG, phiG);
406 
407  std::pair<int,int> pnpair=MEEBGeom::pn(module);
408  unsigned int MyPn0=pnpair.first;
409  unsigned int MyPn1=pnpair.second;
410 
411  unsigned int channel=MEEBGeom::electronic_channel( etaL, phiL );
412  assert(channel < nCrys);
413 
414  double adcmax=0.0;
415 
416  if(towerID!=int(_tower) || channelID!=int(_channel) || dccID!=int(_fedid-600)) continue;
417  else channelNumber=channel;
418 
419  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { // Loop on adc samples
420 
421  EcalMGPASample samp_crystal(df.sample(i));
422  adc[i]=samp_crystal.adc() ;
423  adcG[i]=samp_crystal.gainId();
424 
425  if (i==0) adcGain=adcG[i];
426  if (i>0) adcGain=TMath::Max(adcG[i],adcGain);
427 
428  if (adc[i]>adcmax) {
429  adcmax=adc[i];
430  }
431  }
432 
433  for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
434  dsum+=adc[k];
435  if(k<_presample-1) dsum1+=adc[k];
436  }
437 
438  bl=dsum/((double)_presample);
439  bl1=dsum1/((double)_presample-1);
440 
441  for(val_max=0.,k=0;k<_nsamples;k++) {
442  yrange[k]=adc[k] - bl;
443  if(yrange[k] > val_max) {
444  val_max= yrange[k]; samplemax=k;
445  }
446  }
447 
448  if(samplemax==4 || samplemax==5) {
449  val_max=val_max+bl-bl1;
450  for(k=0;k<_nsamples;k++) {
451  yrange[k]=yrange[k]+bl-bl1;
452  }
453  }
454 
455  for(unsigned int k=0;k<_nsamples;k++) {
456  adc[k]=yrange[k];
457  }
458 
459  pn0=allPNAmpl[MyPn0];
460  pn1=allPNAmpl[MyPn1];
461 
462  if(samplemax>=_timingcutlow && samplemax<=_timingcuthigh && lightside==side) ADCtrees->Fill();
463 
464  }
465 
466  } else {
467 
468  for ( EEDigiCollection::const_iterator digiItr= EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr ) { // Loop on crystals
469 
470  EEDetId id_crystal(digiItr->id()) ;
471  EEDataFrame df( *digiItr );
472 
473  phi = id_crystal.ix()-1 ;
474  eta = id_crystal.iy()-1 ;
475  side=0; // FIXME
476 
477  // Recover the TT id and the electronic crystal numbering from EcalElectronicsMapping
478 
479  EcalElectronicsId elecid_crystal = TheMapping->getElectronicsId(id_crystal);
480 
481  int towerID=elecid_crystal.towerId();
482  int channelID=elecid_crystal.channelId()-1;
483 
484  int module=MEEEGeom::lmmod( phi, eta );
485 
486  std::pair<int,int> pnpair=MEEEGeom::pn( module, _fedid ) ;
487  unsigned int MyPn0=pnpair.first;
488  unsigned int MyPn1=pnpair.second;
489 
490  unsigned int channel=MEEEGeom::crystal( phi, eta );
491  assert(channel < nCrys);
492 
493  double adcmax=0.0;
494 
495  if(towerID!=int(_tower) || channelID!=int(_channel) || dccID!=int(_fedid-600)) continue;
496  else channelNumber=channel;
497 
498  for (unsigned int i=0; i< (*digiItr).size() ; ++i ) { // Loop on adc samples
499 
500  EcalMGPASample samp_crystal(df.sample(i));
501  adc[i]=samp_crystal.adc() ;
502  adcG[i]=samp_crystal.gainId();
503 
504  if (i==0) adcGain=adcG[i];
505  if (i>0) adcGain=TMath::Max(adcG[i],adcGain);
506 
507  if (adc[i]>adcmax) {
508  adcmax=adc[i];
509  }
510  }
511 
512  for(dsum=0.,dsum1=0.,k=0;k<_presample;k++) {
513  dsum+=adc[k];
514  if(k<_presample-1) dsum1+=adc[k];
515  }
516 
517  bl=dsum/((double)_presample);
518  bl1=dsum1/((double)_presample-1);
519 
520  for(val_max=0.,k=0;k<_nsamples;k++) {
521  yrange[k]=adc[k] - bl;
522  if(yrange[k] > val_max) {
523  val_max= yrange[k]; samplemax=k;
524  }
525  }
526 
527  if(samplemax==4 || samplemax==5) {
528  val_max=val_max+bl-bl1;
529  for(k=0;k<_nsamples;k++) {
530  yrange[k]=yrange[k]+bl-bl1;
531  }
532  }
533 
534  for(unsigned int k=0;k<_nsamples;k++) {
535  adc[k]=yrange[k];
536  }
537 
538  pn0=allPNAmpl[MyPn0];
539  pn1=allPNAmpl[MyPn1];
540 
541  if(samplemax>=_timingcutlow && samplemax<=_timingcuthigh && lightside==side) ADCtrees->Fill();
542 
543  }
544  }
545 
546 }// analyze
547 
548 
549 //========================================================================
550 void EcalPerEvtLaserAnalyzer::endJob() {
551 //========================================================================
552 
553 
554  assert( colors.size()<= nColor );
555  unsigned int nCol=colors.size();
556 
557 
558  ADCtrees->Write();
559 
560 
561  cout << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
562  cout << "\t+=+ Analyzing laser data: getting per event +=+" << endl;
563  cout << "\t+=+ APD Amplitudes and ADC samples +=+"<< endl;
564  cout << "\t+=+ for fed:" <<_fedid<<", tower:"<<_tower<<", and channel:" << _channel << endl;
565 
566 
567  // Define temporary tree to save APD amplitudes
568 
569  APDFile = new TFile(resfile.c_str(),"RECREATE");
570 
571  int ieta, iphi, channelID, towerID, flag;
572  double alpha, beta;
573 
574 
575  colors.push_back( color );
576 
577  for (unsigned int i=0;i<nCol;i++){
578 
579  stringstream name1;
580  name1<<"headerCol"<<colors[i];
581 
582  header[i] = new TTree(name1.str().c_str(),name1.str().c_str());
583 
584  header[i]->Branch( "alpha", &alpha, "alpha/D" );
585  header[i]->Branch( "beta", &beta, "beta/D" );
586  header[i]->Branch( "iphi", &iphi, "iphi/I" );
587  header[i]->Branch( "ieta", &ieta, "ieta/I" );
588  header[i]->Branch( "dccID", &dccID, "dccID/I" );
589  header[i]->Branch( "towerID", &towerID, "towerID/I" );
590  header[i]->Branch( "channelID", &channelID, "channelID/I" );
591 
592  header[i]->SetBranchAddress( "alpha", &alpha );
593  header[i]->SetBranchAddress( "beta", &beta );
594  header[i]->SetBranchAddress( "iphi", &iphi );
595  header[i]->SetBranchAddress( "ieta", &ieta );
596  header[i]->SetBranchAddress( "dccID", &dccID );
597  header[i]->SetBranchAddress( "towerID", &towerID );
598  header[i]->SetBranchAddress( "channelID", &channelID );
599 
600  }
601 
602  stringstream name2;
603  name2<<"APDTree";
604  APDtrees= new TTree(name2.str().c_str(),name2.str().c_str());
605 
606  //List of branches
607 
608  APDtrees->Branch( "event", &event, "event/I" );
609  APDtrees->Branch( "color", &color, "color/I" );
610  APDtrees->Branch( "adc", &adc , "adc[10]/D" );
611  APDtrees->Branch( "peakMatacq", &peak , "peak/D" );
612  APDtrees->Branch( "ttrigMatacq", &ttrig , "ttrig/D" );
613  APDtrees->Branch( "apdAmpl", &apdAmpl, "apdAmpl/D" );
614  APDtrees->Branch( "apdTime", &apdTime, "apdTime/D" );
615  APDtrees->Branch( "flag", &flag, "flag/I" );
616  APDtrees->Branch( "pn0", &pn0, "pn0/D" );
617  APDtrees->Branch( "pn1", &pn1, "pn1/D" );
618 
619  APDtrees->SetBranchAddress( "event", &event );
620  APDtrees->SetBranchAddress( "color", &color );
621  APDtrees->SetBranchAddress( "adc", adc );
622  APDtrees->SetBranchAddress( "peakMatacq", &peak );
623  APDtrees->SetBranchAddress( "ttrigMatacq", &ttrig );
624  APDtrees->SetBranchAddress( "apdAmpl", &apdAmpl );
625  APDtrees->SetBranchAddress( "apdTime", &apdTime );
626  APDtrees->SetBranchAddress( "flag", &flag );
627  APDtrees->SetBranchAddress( "pn0", &pn0 );
628  APDtrees->SetBranchAddress( "pn1", &pn1 );
629 
630  // Retrieve alpha and beta for APD amplitudes calculation
631 
632 
633  TFile *alphaFile = new TFile(refalphabeta_.c_str());
634  TTree *alphaTree[2];
635 
636  Double_t alphaRun, betaRun;
637  int ietaRun, iphiRun, channelIDRun, towerIDRun, dccIDRun, flagRun;
638 
639  for( unsigned int i=0;i<nCol;i++){
640 
641  stringstream name3;
642  name3<<"ABCol"<<i;
643  alphaTree[i]=(TTree*)alphaFile->Get(name3.str().c_str());
644  alphaTree[i]->SetBranchStatus( "*", 0 );
645  alphaTree[i]->SetBranchStatus( "alpha", 1 );
646  alphaTree[i]->SetBranchStatus( "beta", 1 );
647  alphaTree[i]->SetBranchStatus( "iphi", 1 );
648  alphaTree[i]->SetBranchStatus( "ieta", 1 );
649  alphaTree[i]->SetBranchStatus( "dccID", 1 );
650  alphaTree[i]->SetBranchStatus( "towerID", 1 );
651  alphaTree[i]->SetBranchStatus( "channelID", 1 );
652  alphaTree[i]->SetBranchStatus( "flag", 1 );
653 
654  alphaTree[i]->SetBranchAddress( "alpha", &alphaRun );
655  alphaTree[i]->SetBranchAddress( "beta", &betaRun );
656  alphaTree[i]->SetBranchAddress( "iphi", &iphiRun );
657  alphaTree[i]->SetBranchAddress( "ieta", &ietaRun );
658  alphaTree[i]->SetBranchAddress( "dccID", &dccIDRun );
659  alphaTree[i]->SetBranchAddress( "towerID", &towerIDRun );
660  alphaTree[i]->SetBranchAddress( "channelID", &channelIDRun );
661  alphaTree[i]->SetBranchAddress( "flag", &flagRun );
662 
663  }
664 
665  PulseFitWithFunction * pslsfit = new PulseFitWithFunction();
666 
667  double chi2;
668 
669  for (unsigned int icol=0;icol<nCol;icol++){
670 
671  IsThereDataADC[icol]=1;
672  stringstream cut;
673  cut <<"color=="<<colors.at(icol);
674  if(ADCtrees->GetEntries(cut.str().c_str())<10) IsThereDataADC[icol]=0;
675  IsHeaderFilled[icol]=0;
676 
677  }
678 
679  // Define submodule and channel number inside the submodule (as Patrice)
680 
681 
682  Long64_t nbytes = 0, nb = 0;
683  for (Long64_t jentry=0; jentry< ADCtrees->GetEntriesFast();jentry++) { // Loop on events
684  nb = ADCtrees->GetEntry(jentry); nbytes += nb;
685 
686  int iCry=channelNumber;
687 
688  // get back color
689 
690  unsigned int iCol=0;
691  for(unsigned int i=0;i<nCol;i++){
692  if(color==colors[i]) {
693  iCol=i;
694  i=colors.size();
695  }
696  }
697 
698  alphaTree[iCol]->GetEntry(iCry);
699 
700 
701  flag=flagRun;
702  iphi=iphiRun;
703  ieta=ietaRun;
704  towerID=towerIDRun;
705  channelID=channelIDRun;
706  alpha=alphaRun;
707  beta=betaRun;
708 
709  if (IsHeaderFilled[iCol]==0){
710  header[iCol]->Fill();
711  IsHeaderFilled[iCol]=1;
712 
713  }
714  // Amplitude calculation
715 
716  apdAmpl=0;
717  apdTime=0;
718 
719  pslsfit -> init(_nsamples,_firstsample,_lastsample,_niter, alphaRun, betaRun);
720  chi2 = pslsfit -> doFit(&adc[0]);
721 
722  if( chi2 < 0. || chi2 == 102 || chi2 == 101 ) {
723 
724  apdAmpl=0;
725  apdTime=0;
726 
727  }else{
728 
729  apdAmpl = pslsfit -> getAmpl();
730  apdTime = pslsfit -> getTime();
731 
732  }
733 
734  APDtrees->Fill();
735 
736  }
737 
738  alphaFile->Close();
739 
740  ADCFile->Close();
741 
742  APDFile->Write();
743  APDFile->Close();
744 
745 
746 
747 
748  // Remove unwanted files
749 
750  stringstream del;
751  del << "rm " <<ADCfile;
752  system(del.str().c_str());
753 
754  cout << "\t+=+ .................................................. done +=+" << endl;
755  cout << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
756 }
757 
758 
759 
760 DEFINE_FWK_MODULE(EcalPerEvtLaserAnalyzer);
761 
beta
const double beta
Definition: CosmicMuonParameters.h:69
MEEBGeom::localCoord
static XYCoord localCoord(int icr)
Definition: MEEBGeom.cc:153
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
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
alpha
float alpha
Definition: AMPTWrapper.h:95
MEEEGeom::crystal
static int crystal(CrysCoord ix, CrysCoord iy)
Definition: MEEEGeom.cc:398
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
NSIDES
#define NSIDES
Definition: EcalLaserAnalyzer.h:44
EcalElectronicsId::stripId
int stripId() const
get the tower id
Definition: EcalElectronicsId.h:39
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
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
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
MEEBGeom::pn
static std::pair< int, int > pn(int ilmmod)
Definition: MEEBGeom.cc:479
nTT
#define nTT
Definition: TMEGeom.h:6
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
EcalPerEvtLaserAnalyzer::analyze
virtual void analyze(const edm::Event &e, const edm::EventSetup &c)
Definition: EcalPerEvtLaserAnalyzer.cc:167
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:405
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
environment_file_cff.runType
tuple runType
Definition: environment_file_cff.py:30
EcalPerEvtLaserAnalyzer::EcalPerEvtLaserAnalyzer
EcalPerEvtLaserAnalyzer(const edm::ParameterSet &iConfig)
Definition: EcalPerEvtLaserAnalyzer.cc:47
edm::Handle::product
T const * product() const
Definition: Handle.h:81
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
run_regression.test
tuple test
Definition: run_regression.py:386
dtNoiseDBValidation_cfg.cerr
tuple cerr
Definition: dtNoiseDBValidation_cfg.py:22
edm::DataFrameContainer::end
const_iterator end() const
Definition: DataFrameContainer.h:156
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
relval_steps.k
list k
Definition: relval_steps.py:1469
module
Definition: vlib.h:208
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
GOODCOLL_filter_cfg.cut
tuple cut
Definition: GOODCOLL_filter_cfg.py:65
Generated for CMSSW Reference Manual by doxygen 1.8.5