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