CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalSimpleTBAnalyzer.cc
Go to the documentation of this file.
1 
8 //
9 // $Id: EcalSimpleTBAnalyzer.cc,v 1.8 2007/12/21 16:10:55 ferriff Exp $
10 //
11 //
12 
20 
21 
22 //#include<fstream>
23 
24 #include "TFile.h"
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TF1.h"
28 
29 #include <iostream>
30 #include <string>
31 #include <stdexcept>
32 //
33 // constants, enums and typedefs
34 //
35 
36 //
37 // static data member definitions
38 //
39 
40 //
41 // constructors and destructor
42 //
43 
44 
45 //========================================================================
47 //========================================================================
48 {
49  //now do what ever initialization is needed
50  rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile","ecalSimpleTBanalysis.root");
51  digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
52  digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
53  hitCollection_ = iConfig.getParameter<std::string>("hitCollection");
54  hitProducer_ = iConfig.getParameter<std::string>("hitProducer");
55  hodoRecInfoCollection_ = iConfig.getParameter<std::string>("hodoRecInfoCollection");
56  hodoRecInfoProducer_ = iConfig.getParameter<std::string>("hodoRecInfoProducer");
57  tdcRecInfoCollection_ = iConfig.getParameter<std::string>("tdcRecInfoCollection");
58  tdcRecInfoProducer_ = iConfig.getParameter<std::string>("tdcRecInfoProducer");
59  eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
60  eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
61 
62 
63  std::cout << "EcalSimpleTBAnalyzer: fetching hitCollection: " << hitCollection_.c_str()
64  << " produced by " << hitProducer_.c_str() << std::endl;
65 
66 }
67 
68 
69 //========================================================================
71 //========================================================================
72 {
73  // do anything here that needs to be done at desctruction time
74  // (e.g. close files, deallocate resources etc.)
75  // Amplitude vs TDC offset
76 // if (h_ampltdc)
77 // delete h_ampltdc;
78 
79 // // Reconstructed energies
80 // delete h_e1x1;
81 // delete h_e3x3;
82 // delete h_e5x5;
83 
84 // delete h_bprofx;
85 // delete h_bprofy;
86 
87 // delete h_qualx;
88 // delete h_qualy;
89 
90 // delete h_slopex;
91 // delete h_slopey;
92 
93 // delete h_mapx;
94 // delete h_mapy;
95 
96 }
97 
98 //========================================================================
99 void
101 //========================================================================
102 
103  // Amplitude vs TDC offset
104  h_ampltdc = new TH2F("h_ampltdc","Max Amplitude vs TDC offset", 100,0.,1.,1000, 0., 4000.);
105 
106  // Reconstructed energies
107  h_tableIsMoving = new TH1F("h_tableIsMoving","TableIsMoving", 100000, 0., 100000.);
108 
109  h_e1x1 = new TH1F("h_e1x1","E1x1 energy", 1000, 0., 4000.);
110  h_e3x3 = new TH1F("h_e3x3","E3x3 energy", 1000, 0., 4000.);
111  h_e5x5 = new TH1F("h_e5x5","E5x5 energy", 1000, 0., 4000.);
112 
113  h_e1x1_center = new TH1F("h_e1x1_center","E1x1 energy", 1000, 0., 4000.);
114  h_e3x3_center = new TH1F("h_e3x3_center","E3x3 energy", 1000, 0., 4000.);
115  h_e5x5_center = new TH1F("h_e5x5_center","E5x5 energy", 1000, 0., 4000.);
116 
117  h_e1e9 = new TH1F("h_e1e9","E1/E9 ratio", 600, 0., 1.2);
118  h_e1e25 = new TH1F("h_e1e25","E1/E25 ratio", 600, 0., 1.2);
119  h_e9e25 = new TH1F("h_e9e25","E9/E25 ratio", 600, 0., 1.2);
120 
121  h_bprofx = new TH1F("h_bprofx","Beam Profile X",100,-20.,20.);
122  h_bprofy = new TH1F("h_bprofy","Beam Profile Y",100,-20.,20.);
123 
124  h_qualx = new TH1F("h_qualx","Beam Quality X",5000,0.,5.);
125  h_qualy = new TH1F("h_qualy","Beam Quality X",5000,0.,5.);
126 
127  h_slopex = new TH1F("h_slopex","Beam Slope X",500, -5e-4 , 5e-4 );
128  h_slopey = new TH1F("h_slopey","Beam Slope Y",500, -5e-4 , 5e-4 );
129 
130  char hname[50];
131  char htitle[50];
132  for (unsigned int icry=0;icry<25;icry++)
133  {
134  sprintf(hname,"h_mapx_%d",icry);
135  sprintf(htitle,"Max Amplitude vs X %d",icry);
136  h_mapx[icry] = new TH2F(hname,htitle,80,-20,20,1000,0.,4000.);
137  sprintf(hname,"h_mapy_%d",icry);
138  sprintf(htitle,"Max Amplitude vs Y %d",icry);
139  h_mapy[icry] = new TH2F(hname,htitle,80,-20,20,1000,0.,4000.);
140  }
141 
142  h_e1e9_mapx = new TH2F("h_e1e9_mapx","E1/E9 vs X",80,-20,20,600,0.,1.2);
143  h_e1e9_mapy = new TH2F("h_e1e9_mapy","E1/E9 vs Y",80,-20,20,600,0.,1.2);
144 
145  h_e1e25_mapx = new TH2F("h_e1e25_mapx","E1/E25 vs X",80,-20,20,600,0.,1.2);
146  h_e1e25_mapy = new TH2F("h_e1e25_mapy","E1/E25 vs Y",80,-20,20,600,0.,1.2);
147 
148  h_e9e25_mapx = new TH2F("h_e9e25_mapx","E9/E25 vs X",80,-20,20,600,0.,1.2);
149  h_e9e25_mapy = new TH2F("h_e9e25_mapy","E9/E25 vs Y",80,-20,20,600,0.,1.2);
150 
151  h_Shape_ = new TH2F("h_Shape_","Xtal in Beam Shape",250,0,10,350,0,3500);
152 
153 }
154 
155 //========================================================================
156 void
158 //========================================================================
159 
160  TFile f(rootfile_.c_str(),"RECREATE");
161 
162  // Amplitude vs TDC offset
163  h_ampltdc->Write();
164 
165  // Reconstructed energies
166  h_e1x1->Write();
167  h_e3x3->Write();
168  h_e5x5->Write();
169 
170  h_e1x1_center->Write();
171  h_e3x3_center->Write();
172  h_e5x5_center->Write();
173 
174  h_e1e9->Write();
175  h_e1e25->Write();
176  h_e9e25->Write();
177 
178  h_bprofx->Write();
179  h_bprofy->Write();
180 
181  h_qualx->Write();
182  h_qualy->Write();
183 
184  h_slopex->Write();
185  h_slopey->Write();
186 
187  h_Shape_->Write();
188 
189  for (unsigned int icry=0;icry<25;icry++)
190  {
191  h_mapx[icry]->Write();
192  h_mapy[icry]->Write();
193  }
194 
195  h_e1e9_mapx->Write();
196  h_e1e9_mapy->Write();
197 
198  h_e1e25_mapx->Write();
199  h_e1e25_mapy->Write();
200 
201  h_e9e25_mapx->Write();
202  h_e9e25_mapy->Write();
203 
204  h_tableIsMoving->Write();
205 
206  f.Close();
207 }
208 
209 //
210 // member functions
211 //
212 
213 //========================================================================
214 void
216 //========================================================================
217 
218  using namespace edm;
219  using namespace cms;
220 
221 
222 
224  const EBDigiCollection* digis=0;
225  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
226  iEvent.getByLabel( digiProducer_, digiCollection_,pdigis);
227  if ( pdigis.isValid() ) {
228  digis = pdigis.product(); // get a ptr to the product
229  //iEvent.getByLabel( hitProducer_, phits);
230  } else {
231  edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << digiCollection_;
232  }
233 
234  // fetch the digis and compute signal amplitude
236  const EBUncalibratedRecHitCollection* hits=0;
237  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
238  iEvent.getByLabel( hitProducer_, hitCollection_,phits);
239  if (phits.isValid()) {
240  hits = phits.product(); // get a ptr to the product
241  //iEvent.getByLabel( hitProducer_, phits);
242  } else {
243  edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hitCollection_;
244  }
245 
247  const EcalTBHodoscopeRecInfo* recHodo=0;
248  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
250  if ( pHodo.isValid() ) {
251  recHodo = pHodo.product(); // get a ptr to the product
252  } else {
253  edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hodoRecInfoCollection_;
254  }
255 
257  const EcalTBTDCRecInfo* recTDC=0;
258  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
260  if ( pTDC.isValid() ) {
261  recTDC = pTDC.product(); // get a ptr to the product
262  } else {
263  edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << tdcRecInfoCollection_;
264  }
265 
266  Handle<EcalTBEventHeader> pEventHeader;
267  const EcalTBEventHeader* evtHeader=0;
268  //std::cout << "EcalSimpleTBAnalyzer::analyze getting product with label: " << digiProducer_.c_str()<< " prodname: " << digiCollection_.c_str() << endl;
269  iEvent.getByLabel( eventHeaderProducer_ , pEventHeader );
270  if ( pEventHeader.isValid() ) {
271  evtHeader = pEventHeader.product(); // get a ptr to the product
272  } else {
273  edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << eventHeaderProducer_;
274  }
275 
276  if (!hits)
277  return;
278 
279  if (!recTDC)
280  return;
281 
282  if (!recHodo)
283  return;
284 
285  if (!evtHeader)
286  return;
287 
288  if (hits->size() == 0)
289  return;
290 
291  if (evtHeader->tableIsMoving())
292  h_tableIsMoving->Fill(evtHeader->eventNumber());
293 
294  // Crystal hit by beam
295  if (xtalInBeam_.null())
296  {
297  xtalInBeam_ = EBDetId(1,evtHeader->crystalInBeam(),EBDetId::SMCRYSTALMODE);
298  std::cout<< "Xtal In Beam is " << xtalInBeam_.ic() << std::endl;
299  }
300  else if (xtalInBeam_ != EBDetId(1,evtHeader->crystalInBeam(),EBDetId::SMCRYSTALMODE))
301  return;
302 
303  if (evtHeader->tableIsMoving())
304  return;
305 
306 
307 // EBDetId maxHitId(0);
308 // float maxHit= -999999.;
309 
310 // for(EBUncalibratedRecHitCollection::const_iterator ithit = hits->begin(); ithit != hits->end(); ++ithit)
311 // {
312 // if (ithit->amplitude()>=maxHit)
313 // {
314 // maxHit=ithit->amplitude();
315 // maxHitId=ithit->id();
316 // }
317 
318 // }
319 
320 // if (maxHitId==EBDetId(0))
321 // return;
322 
323 // EBDetId maxHitId(1,704,EBDetId::SMCRYSTALMODE);
324 
325  //Find EBDetId in a 5x5 Matrix (to be substituted by the Selector code)
326  // Something like
327  // EBFixedWindowSelector<EcalUncalibratedRecHit> Simple5x5Matrix(hits,maxHitId,5,5);
328  // std::vector<EcalUncalibratedRecHit> Energies5x5 = Simple5x5Matrix.getHits();
329 
330 
331  EBDetId Xtals5x5[25];
332  for (unsigned int icry=0;icry<25;icry++)
333  {
334  unsigned int row = icry / 5;
335  unsigned int column= icry %5;
336  int ieta=xtalInBeam_.ieta()+column-2;
337  int iphi=xtalInBeam_.iphi()+row-2;
338  EBDetId tempId(ieta, iphi,EBDetId::ETAPHIMODE);
339  if (tempId.ism()==1)
340  Xtals5x5[icry]=tempId;
341  else
342  Xtals5x5[icry]=EBDetId(0);
344  }
345 
346 
347 
348  bool gain_switch = false;
349  double samples_save[10]; for(int i=0; i < 10; ++i) samples_save[i]=0.0;
350  double gain_save[10]; for(int i=0; i < 10; ++i) gain_save[i]=0.0;
351 
352  // find the rechit corresponding digi and the max sample
353  EBDigiCollection::const_iterator myDg = digis->find(xtalInBeam_);
354  int sMax = -1;
355  double eMax = 0.;
356  if (myDg != digis->end())
357  {
358  EBDataFrame myDigi = (*myDg);
359  for (int sample = 0; sample < myDigi.size(); ++sample)
360  {
361  double analogSample = myDigi.sample(sample).adc();
362  double gainSample = myDigi.sample(sample).gainId();
363  samples_save[sample] = analogSample;
364  gain_save[sample] = gainSample;
365  // std::cout << analogSample << " ";
366  if ( eMax < analogSample )
367  {
368  eMax = analogSample;
369  sMax = sample;
370  }
371  if(gainSample != 1) gain_switch = true;
372  }
373  // std::cout << std::endl;
374  }
375 
376  for(int i =0; i < 10; ++i) {
377  h_Shape_->Fill(double(i)+recTDC->offset(),samples_save[i]);
378  }
379 
380  double amplitude[25];
381 
382  double amplitude3x3=0;
383  double amplitude5x5=0;
384 
385  for (unsigned int icry=0;icry<25;icry++)
386  {
387  if (!Xtals5x5[icry].null())
388  {
389  amplitude[icry]=(hits->find(Xtals5x5[icry]))->amplitude();
390  amplitude5x5 += amplitude[icry];
391  // Is in 3x3?
392  if ( icry == 6 || icry == 7 || icry == 8 ||
393  icry == 11 || icry == 12 || icry ==13 ||
394  icry == 16 || icry == 17 || icry ==18 )
395  {
396  amplitude3x3+=amplitude[icry];
397  }
398  }
399  }
400 
401 
402  h_e1x1->Fill(amplitude[12]);
403  h_e3x3->Fill(amplitude3x3);
404  h_e5x5->Fill(amplitude5x5);
405 
406  h_e1e9->Fill(amplitude[12]/amplitude3x3);
407  h_e1e25->Fill(amplitude[12]/amplitude5x5);
408  h_e9e25->Fill(amplitude3x3/amplitude5x5);
409 
410  if (recTDC)
411  h_ampltdc->Fill(recTDC->offset(),amplitude[12]);
412 
413  if (recHodo)
414  {
415  float x=recHodo->posX();
416  float y=recHodo->posY();
417  float xslope=recHodo->slopeX();
418  float yslope=recHodo->slopeY();
419  float xqual=recHodo->qualX();
420  float yqual=recHodo->qualY();
421 
422  //Filling beam profiles
423  h_bprofx->Fill(x);
424  h_bprofy->Fill(y);
425  h_qualx->Fill(xqual);
426  h_qualy->Fill(yqual);
427  h_slopex->Fill(xslope);
428  h_slopey->Fill(yslope);
429 
430  //Fill central events
431 
432 
433  if ( fabs(x + 2.5) < 2.5 && fabs(y + 0.5) < 2.5)
434  {
435  h_e1x1_center->Fill(amplitude[12]);
436  h_e3x3_center->Fill(amplitude3x3);
437  h_e5x5_center->Fill(amplitude5x5);
438 
439  h_e1e9->Fill(amplitude[12]/amplitude3x3);
440  h_e1e25->Fill(amplitude[12]/amplitude5x5);
441  h_e9e25->Fill(amplitude3x3/amplitude5x5);
442  }
443 
444  for (unsigned int icry=0;icry<25;icry++)
445  {
446  h_mapx[icry]->Fill(x,amplitude[icry]);
447  h_mapy[icry]->Fill(y,amplitude[icry]);
448  }
449 
450  h_e1e9_mapx->Fill(x,amplitude[12]/amplitude3x3);
451  h_e1e9_mapy->Fill(y,amplitude[12]/amplitude3x3);
452 
453  h_e1e25_mapx->Fill(x,amplitude[12]/amplitude5x5);
454  h_e1e25_mapy->Fill(y,amplitude[12]/amplitude5x5);
455 
456  h_e9e25_mapx->Fill(x,amplitude3x3/amplitude5x5);
457  h_e9e25_mapy->Fill(y,amplitude3x3/amplitude5x5);
458  }
459 
460 }
461 
462 
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
std::string hodoRecInfoCollection_
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
std::string eventHeaderCollection_
virtual void analyze(const edm::Event &, const edm::EventSetup &)
EcalMGPASample sample(int i) const
Definition: EcalDataFrame.h:28
int gainId() const
get the gainId (2 bits)
int ism() const
get the ECAL/SM id
Definition: EBDetId.cc:79
int size() const
Definition: EcalDataFrame.h:25
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
int iEvent
Definition: GenABIO.cc:243
double f[11][100]
int ieta() const
get the crystal ieta
Definition: EBDetId.h:44
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
static const int ETAPHIMODE
Definition: EBDetId.h:145
int ic() const
get ECAL/crystal number inside SM
Definition: EBDetId.cc:94
EcalSimpleTBAnalyzer(const edm::ParameterSet &)
bool null() const
is this a null id ?
Definition: DetId.h:47
tuple cout
Definition: gather_cfg.py:41
Definition: DDAxes.h:10
static const int SMCRYSTALMODE
Definition: EBDetId.h:146
int adc() const
get the ADC sample (12 bits)