CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ESTimingTask.cc
Go to the documentation of this file.
1 #include <memory>
2 #include <fstream>
3 #include <iostream>
4 
19 
20 #include "TStyle.h"
21 #include "TH2F.h"
22 #include "TMath.h"
23 #include "TGraph.h"
24 
25 using namespace cms;
26 using namespace edm;
27 using namespace std;
28 
29 // fit function
30 double fitf(double *x, double *par) {
31 
32  double wc = par[2];
33  double n = par[3]; // n-1 (in fact)
34  double v1 = pow(wc/n*(x[0]-par[1]), n);
35  double v2 = TMath::Exp(n-wc*(x[0]-par[1]));
36  double v = par[0]*v1*v2;
37 
38  if (x[0] < par[1]) v = 0;
39 
40  return v;
41 }
42 
44 
45  rechitlabel_ = ps.getParameter<InputTag>("RecHitLabel");
46  digilabel_ = ps.getParameter<InputTag>("DigiLabel");
47  prefixME_ = ps.getUntrackedParameter<string>("prefixME", "EcalPreshower");
48 
49  dqmStore_ = Service<DQMStore>().operator->();
50  eCount_ = 0;
51 
52  fit_ = new TF1("fitShape", fitf, -200, 200, 4);
53  fit_->SetParameters(50, 10, 0, 0);
54 
55  //Histogram init
56  for (int i = 0; i < 2; ++i)
57  for (int j = 0; j < 2; ++j)
58  hTiming_[i][j] = 0;
59 
60  dqmStore_->setCurrentFolder(prefixME_ + "/ESTimingTask");
61 
62  //Booking Histograms
63  //Notice: Change ESRenderPlugin under DQM/RenderPlugins/src if you change this histogram name.
64  char histo[200];
65  for (int i=0 ; i<2; ++i)
66  for (int j=0 ; j<2; ++j) {
67  int iz = (i==0)? 1:-1;
68  sprintf(histo, "ES Timing Z %d P %d", iz, j+1);
69  hTiming_[i][j] = dqmStore_->book1D(histo, histo, 81, -20.5, 20.5);
70  hTiming_[i][j]->setAxisTitle("ES Timing (ns)", 1);
71  }
72 
73  sprintf(histo, "ES 2D Timing");
74  h2DTiming_ = dqmStore_->book2D(histo, histo, 81, -20.5, 20.5, 81, -20.5, 20.5);
75  h2DTiming_->setAxisTitle("ES- Timing (ns)", 1);
76  h2DTiming_->setAxisTitle("ES+ Timing (ns)", 2);
77 
78  htESP_ = new TH1F("htESP", "Timing ES+", 81, -20.5, 20.5);
79  htESM_ = new TH1F("htESM", "Timing ES-", 81, -20.5, 20.5);
80 }
81 
83  delete htESP_;
84  delete htESM_;
85 }
86 
88 }
89 
91 }
92 
93 void ESTimingTask::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
94 
95  set(iSetup);
96 
97  runNum_ = e.id().run();
98  eCount_++;
99 
100  htESP_->Reset();
101  htESM_->Reset();
102 
103  //Digis
104  int zside, plane, ix, iy, is;
105  double adc[3];
106  // double para[10];
107  //double tx[3] = {-5., 20., 45.};
109  if ( e.getByLabel(digilabel_, digis) ) {
110 
111  for (ESDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
112 
113  ESDataFrame dataframe = (*digiItr);
114  ESDetId id = dataframe.id();
115 
116  zside = id.zside();
117  plane = id.plane();
118  ix = id.six();
119  iy = id.siy();
120  is = id.strip();
121 
122  //if (zside==1 && plane==1 && ix==15 && iy==6) continue;
123  if (zside==1 && plane==1 && ix==7 && iy==28) continue;
124  if (zside==1 && plane==1 && ix==24 && iy==9 && is==21) continue;
125  if (zside==-1 && plane==2 && ix==35 && iy==17 && is==23) continue;
126 
127  int i = (zside==1)? 0:1;
128  int j = plane-1;
129 
130  for (int k=0; k<dataframe.size(); ++k)
131  adc[k] = dataframe.sample(k).adc();
132 
133  double status = 0;
134  if (adc[1] < 200) status = 1;
135  if (fabs(adc[0]) > 10) status = 1;
136  if (adc[1] < 0 || adc[2] < 0) status = 1;
137  if (adc[0] > adc[1] || adc[0] > adc[2]) status = 1;
138  if (adc[2] > adc[1]) status = 1;
139 
140  if (int(status) == 0) {
141 
142  double A1 = adc[1];
143  double A2 = adc[2];
144  double DeltaT = 25.;
145  double aaa = (A2 > 0 && A1 > 0) ? log(A2/A1)/n_ : 20.; // if A1=0, t0=20
146  double bbb = wc_/n_*DeltaT;
147  double ccc= exp(aaa+bbb);
148 
149  double t0 = (2.-ccc)/(1.-ccc) * DeltaT - 5;
150  hTiming_[i][j]->Fill(t0);
151  //cout<<"t0 : "<<t0<<endl;
152  /*
153  TGraph *gr = new TGraph(3, tx, adc);
154  fit_->SetParameters(50, 10, wc_, n_);
155  fit_->FixParameter(2, wc_);
156  fit_->FixParameter(3, n_);
157  fit_->Print();
158  gr->Fit("fitShape", "MQ");
159  fit_->GetParameters(para);
160  delete gr;
161  //hTiming_[i][j]->Fill(para[1]);
162  */
163  //cout<<"ADC : "<<zside<<" "<<plane<<" "<<ix<<" "<<iy<<" "<<is<<" "<<adc[0]<<" "<<adc[1]<<" "<<adc[2]<<" "<<para[1]<<" "<<wc_<<" "<<n_<<endl;
164 
165  if (zside == 1) htESP_->Fill(t0);
166  else if (zside == -1) htESM_->Fill(t0);
167  }
168 
169  }
170  } else {
171  LogWarning("ESTimingTask") << digilabel_ << " not available";
172  }
173 
174  if (htESP_->GetEntries() > 0 && htESM_->GetEntries() > 0)
175  h2DTiming_->Fill(htESM_->GetMean(), htESP_->GetMean());
176 
177 }
178 
180 
181 
182  es.get<ESGainRcd>().get(esgain_);
183  const ESGain *gain = esgain_.product();
184 
185  int ESGain = (int) gain->getESGain();
186 
187  if (ESGain == 1) { // LG
188  wc_ = 0.0837264;
189  n_ = 2.016;
190  } else { // HG
191  wc_ = 0.07291;
192  n_ = 1.798;
193  }
194 
195  //cout<<"gain : "<<ESGain<<endl;
196  //cout<<wc_<<" "<<n_<<endl;
197 
198 }
199 
200 //define this as a plug-in
int adc(sample_type sample)
get the ADC sample (12 bits)
RunNumber_t run() const
Definition: EventID.h:42
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
virtual void beginJob(void)
Definition: ESTimingTask.cc:87
double fitf(double *x, double *par)
Definition: ESTimingTask.cc:30
DEFINE_FWK_MODULE(TauMET)
Definition: ESGain.h:5
const ESDetId & id() const
Definition: ESDataFrame.h:21
void set(const edm::EventSetup &es)
int size() const
Definition: ESDataFrame.h:23
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: ESTimingTask.cc:93
virtual ~ESTimingTask()
Definition: ESTimingTask.cc:82
const ESSample & sample(int i) const
Definition: ESDataFrame.h:26
int j
Definition: DBlmapReader.cc:9
int zside() const
Definition: ESDetId.h:33
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
int k[5][pyjets_maxn]
float getESGain() const
Definition: ESGain.h:11
const T & get() const
Definition: EventSetup.h:55
edm::EventID id() const
Definition: EventBase.h:56
virtual void endJob(void)
Definition: ESTimingTask.cc:90
int adc() const
get the ADC sample (singed 16 bits)
Definition: ESSample.h:18
Definition: DDAxes.h:10
tuple status
Definition: ntuplemaker.py:245
ESTimingTask(const edm::ParameterSet &ps)
Definition: ESTimingTask.cc:43
mathSSE::Vec4< T > v
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void set(const std::string &name, int value)
set the flag, with a run-time name