CMS 3D CMS Logo
TShapeAnalysis.cc
Go to the documentation of this file.
1 /*
2  * \class TShapeAnalysis
3  *
4  * original author: Patrice Verrecchia
5  * modified by Julie Malcles - CEA/Saclay
6  */
7 
8 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h>
9 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TFParams.h>
10 
11 #include <iostream>
12 #include <cmath>
13 #include <ctime>
14 #include <cassert>
15 
16 #include <TFile.h>
17 #include <TTree.h>
18 #include <TBranch.h>
19 
20 //ClassImp(TShapeAnalysis)
21 
22 using namespace std;
23 
24 // Constructor...
25 TShapeAnalysis::TShapeAnalysis(double alpha0, double beta0, double width0, double chi20) {
26  init(alpha0, beta0, width0, chi20);
27 }
28 
29 TShapeAnalysis::TShapeAnalysis(TTree *tAB, double alpha0, double beta0, double width0, double chi20) {
30  init(tAB, alpha0, beta0, width0, chi20);
31 }
32 
33 // Destructor
34 TShapeAnalysis::~TShapeAnalysis() {}
35 
36 void TShapeAnalysis::init(double alpha0, double beta0, double width0, double chi20) {
37  tABinit = nullptr;
38  nchsel = fNchsel;
39  for (int cris = 0; cris < fNchsel; cris++) {
40  index[cris] = -1;
41  npass[cris] = 0;
42  npassok[cris] = 0.;
43 
44  alpha_val[cris] = alpha0;
45  beta_val[cris] = beta0;
46  width_val[cris] = width0;
47  chi2_val[cris] = chi20;
48  flag_val[cris] = 0;
49 
50  alpha_init[cris] = alpha0;
51  beta_init[cris] = beta0;
52  width_init[cris] = width0;
53  chi2_init[cris] = chi20;
54  flag_init[cris] = 0;
55 
56  phi_init[cris] = 0;
57  eta_init[cris] = 0;
58  side_init[cris] = 0;
59  dcc_init[cris] = 0;
60  tower_init[cris] = 0;
61  ch_init[cris] = 0;
62  assignChannel(cris, cris);
63  }
64 }
65 
66 void TShapeAnalysis::init(TTree *tAB, double alpha0, double beta0, double width0, double chi20) {
67  init(alpha0, beta0, width0, chi20);
68 
69  tABinit = tAB->CloneTree();
70 
71  // Declaration of leaf types
72  Int_t sidei;
73  Int_t iphii;
74  Int_t ietai;
75  Int_t dccIDi;
76  Int_t towerIDi;
77  Int_t channelIDi;
78  Double_t alphai;
79  Double_t betai;
80  Double_t widthi;
81  Double_t chi2i;
82  Int_t flagi;
83 
84  // List of branches
85  TBranch *b_iphi;
86  TBranch *b_ieta;
87  TBranch *b_side;
88  TBranch *b_dccID;
89  TBranch *b_towerID;
90  TBranch *b_channelID;
91  TBranch *b_alpha;
92  TBranch *b_beta;
93  TBranch *b_width;
94  TBranch *b_chi2;
95  TBranch *b_flag;
96 
97  if (tABinit) {
98  tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
99  tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
100  tABinit->SetBranchAddress("side", &sidei, &b_side);
101  tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
102  tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
103  tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
104  tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
105  tABinit->SetBranchAddress("beta", &betai, &b_beta);
106  tABinit->SetBranchAddress("width", &widthi, &b_width);
107  tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
108  tABinit->SetBranchAddress("flag", &flagi, &b_flag);
109 
110  nchsel = tABinit->GetEntries();
111  assert(nchsel <= fNchsel);
112 
113  for (int cris = 0; cris < nchsel; cris++) {
114  tABinit->GetEntry(cris);
115 
116  // std::cout<< "Loop 1 "<< cris<<" "<<alphai<< std::endl;
117 
118  putalphaVal(cris, alphai);
119  putchi2Val(cris, chi2i);
120  putbetaVal(cris, betai);
121  putwidthVal(cris, widthi);
122  putflagVal(cris, flagi);
123 
124  putalphaInit(cris, alphai);
125  putchi2Init(cris, chi2i);
126  putbetaInit(cris, betai);
127  putwidthInit(cris, widthi);
128  putflagInit(cris, flagi);
129  putetaInit(cris, ietai);
130  putphiInit(cris, iphii);
131  }
132  }
133 }
134 
135 void TShapeAnalysis::set_const(int ns, int ns1, int ns2, int ps, int nevtmax, double noise_val, double chi2_cut) {
136  nsamplecristal = ns;
137  presample = ps;
138  sampbmax = ns1;
139  sampamax = ns2;
140  nevt = nevtmax;
141  noise = noise_val;
142  chi2cut = chi2_cut;
143 }
144 
145 void TShapeAnalysis::set_presample(int ps) { presample = ps; }
146 void TShapeAnalysis::set_nch(int nch) {
147  assert(nch <= fNchsel);
148  if (tABinit)
149  assert(nch == nchsel);
150  nchsel = nch;
151 }
152 void TShapeAnalysis::assignChannel(int n, int ch) {
153  if (n >= nchsel)
154  printf(" number of channels exceed maximum allowed\n");
155 
156  index[n] = ch;
157 }
158 
159 void TShapeAnalysis::putDateStart(long int timecur) { timestart = timecur; }
160 
161 void TShapeAnalysis::putDateStop(long int timecur) { timestop = timecur; }
162 
163 void TShapeAnalysis::getDateStart() {
164  time_t t, timecur;
165  timecur = time(&t);
166  timestart = ((long int)timecur);
167 }
168 
169 void TShapeAnalysis::getDateStop() {
170  time_t t, timecur;
171  timecur = time(&t);
172  timestop = ((long int)timecur);
173 }
174 
175 void TShapeAnalysis::putAllVals(int ch, double *sampl, int ieta, int iphi, int dcc, int side, int tower, int chid) {
176  dcc_init[ch] = dcc;
177  tower_init[ch] = side;
178  ch_init[ch] = chid;
179  side_init[ch] = side;
180  eta_init[ch] = ieta;
181  phi_init[ch] = iphi;
182  putAllVals(ch, sampl, ieta, iphi);
183 }
184 
185 void TShapeAnalysis::putAllVals(int ch, double *sampl, int ieta, int iphi) {
186  int i, k;
187  int n = -1;
188  for (i = 0; i < nchsel; i++)
189  if (index[i] == ch)
190  n = i;
191 
192  if (n >= 0) {
193  if (npass[n] < nevt) {
194  for (k = 0; k < nsamplecristal; k++) {
195  rawsglu[n][npass[n]][k] = sampl[k];
196  }
197 
198  npass[n]++;
199  }
200  } else {
201  printf("no index found for ch=%d\n", ch);
202  }
203 }
204 
205 void TShapeAnalysis::computeShape(string namefile, TTree *tAB) {
206  double tm_atmax[200];
207  double parout[3];
208 
209  double chi2_all = 0.;
210 
211  double **dbi;
212  dbi = new double *[200];
213  for (int k = 0; k < 200; k++)
214  dbi[k] = new double[2];
215 
216  double **signalu;
217  signalu = new double *[200];
218  for (int k = 0; k < 200; k++)
219  signalu[k] = new double[10];
220 
221  TFParams *pjf = new TFParams();
222 
223  for (int i = 0; i < nchsel; i++) {
224  if (index[i] >= 0) {
225  if (npass[i] <= 10) {
226  putalphaVal(i, alpha_init[i]);
227  putbetaVal(i, beta_init[i]);
228  putwidthVal(i, width_init[i]);
229  putchi2Val(i, chi2_init[i]);
230  putflagVal(i, 0);
231 
232  } else {
233  pjf->set_const(nsamplecristal, sampbmax, sampamax, alpha_init[i], beta_init[i], npass[i]);
234 
235  for (int pass = 0; pass < npass[i]; pass++) {
236  double ped = 0;
237  for (int k = 0; k < presample; k++) {
238  ped += rawsglu[i][pass][k];
239  }
240  ped /= double(presample);
241 
242  for (int k = 0; k < nsamplecristal; k++) {
243  signalu[pass][k] = rawsglu[i][pass][k] - ped;
244  }
245  }
246 
247  int debug = 0;
248  chi2_all = pjf->fitpj(signalu, &parout[0], dbi, noise, debug);
249 
250  if (parout[0] >= 0.0 && parout[1] >= 0.0 && chi2_all <= chi2cut && chi2_all > 0.0) {
251  putalphaVal(i, parout[0]);
252  putbetaVal(i, parout[1]);
253  putchi2Val(i, chi2_all);
254  putflagVal(i, 1);
255 
256  } else {
257  putalphaVal(i, alpha_init[i]);
258  putbetaVal(i, beta_init[i]);
259  putwidthVal(i, width_init[i]);
260  putchi2Val(i, chi2_init[i]);
261  putflagVal(i, 0);
262  }
263 
264  for (int kj = 0; kj < npass[i]; kj++) { // last event wrong here
265  tm_atmax[kj] = dbi[kj][1];
266  }
267  computetmaxVal(i, &tm_atmax[0]);
268  }
269  }
270  }
271 
272  if (tAB)
273  tABinit = tAB->CloneTree();
274 
275  // Declaration of leaf types
276  Int_t sidei;
277  Int_t iphii;
278  Int_t ietai;
279  Int_t dccIDi;
280  Int_t towerIDi;
281  Int_t channelIDi;
282  Double_t alphai;
283  Double_t betai;
284  Double_t widthi;
285  Double_t chi2i;
286  Int_t flagi;
287 
288  // List of branches
289  TBranch *b_iphi;
290  TBranch *b_ieta;
291  TBranch *b_side;
292  TBranch *b_dccID;
293  TBranch *b_towerID;
294  TBranch *b_channelID;
295  TBranch *b_alpha;
296  TBranch *b_beta;
297  TBranch *b_width;
298  TBranch *b_chi2;
299  TBranch *b_flag;
300 
301  if (tABinit) {
302  tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
303  tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
304  tABinit->SetBranchAddress("side", &sidei, &b_side);
305  tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
306  tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
307  tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
308  tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
309  tABinit->SetBranchAddress("beta", &betai, &b_beta);
310  tABinit->SetBranchAddress("width", &widthi, &b_width);
311  tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
312  tABinit->SetBranchAddress("flag", &flagi, &b_flag);
313  }
314 
315  TFile *fABout = new TFile(namefile.c_str(), "RECREATE");
316  tABout = new TTree("ABCol0", "ABCol0");
317 
318  // Declaration of leaf types
319  Int_t side;
320  Int_t iphi;
321  Int_t ieta;
322  Int_t dccID;
323  Int_t towerID;
324  Int_t channelID;
325  Double_t alpha;
326  Double_t beta;
327  Double_t width;
328  Double_t chi2;
329  Int_t flag;
330 
331  tABout->Branch("iphi", &iphi, "iphi/I");
332  tABout->Branch("ieta", &ieta, "ieta/I");
333  tABout->Branch("side", &side, "side/I");
334  tABout->Branch("dccID", &dccID, "dccID/I");
335  tABout->Branch("towerID", &towerID, "towerID/I");
336  tABout->Branch("channelID", &channelID, "channelID/I");
337  tABout->Branch("alpha", &alpha, "alpha/D");
338  tABout->Branch("beta", &beta, "beta/D");
339  tABout->Branch("width", &width, "width/D");
340  tABout->Branch("chi2", &chi2, "chi2/D");
341  tABout->Branch("flag", &flag, "flag/I");
342 
343  tABout->SetBranchAddress("ieta", &ieta);
344  tABout->SetBranchAddress("iphi", &iphi);
345  tABout->SetBranchAddress("side", &side);
346  tABout->SetBranchAddress("dccID", &dccID);
347  tABout->SetBranchAddress("towerID", &towerID);
348  tABout->SetBranchAddress("channelID", &channelID);
349  tABout->SetBranchAddress("alpha", &alpha);
350  tABout->SetBranchAddress("beta", &beta);
351  tABout->SetBranchAddress("width", &width);
352  tABout->SetBranchAddress("chi2", &chi2);
353  tABout->SetBranchAddress("flag", &flag);
354 
355  for (int i = 0; i < nchsel; i++) {
356  if (tABinit) {
357  tABinit->GetEntry(i);
358  iphi = iphii;
359  ieta = ietai;
360  side = sidei;
361  dccID = dccIDi;
362  towerID = towerIDi;
363  channelID = channelIDi;
364 
365  } else {
366  iphi = phi_init[i];
367  ieta = eta_init[i];
368  side = side_init[i];
369  dccID = dcc_init[i];
370  towerID = tower_init[i];
371  channelID = ch_init[i];
372  }
373 
374  alpha = alpha_val[i];
375  beta = beta_val[i];
376  width = width_val[i];
377  chi2 = chi2_val[i];
378  flag = flag_val[i];
379 
380  tABout->Fill();
381  }
382 
383  tABout->Write();
384  fABout->Close();
385 
386  delete pjf;
387 }
388 
389 void TShapeAnalysis::computetmaxVal(int i, double *tm_val) {
390  double tm_mean = 0; //double tm_sig=0;
391 
392  double tm = 0.;
393  double sigtm = 0.;
394  for (int k = 0; k < npass[i] - 1; k++) {
395  if (1. < tm_val[k] && tm_val[k] < 10.) {
396  npassok[i]++;
397  tm += tm_val[k];
398  sigtm += tm_val[k] * tm_val[k];
399  }
400  }
401  if (npassok[i] <= 0) {
402  tm_mean = 0.; //tm_sig=0.;
403  } else {
404  for (int k = 0; k < npass[i] - 1; k++) {
405  if (1. < tm_val[k] && tm_val[k] < 10.) {
406  double ss = (sigtm / npassok[i] - tm / npassok[i] * tm / npassok[i]);
407  if (ss < 0.)
408  ss = 0.;
409  //tm_sig=sqrt(ss);
410  tm_mean = tm / npassok[i];
411  }
412  }
413  }
414  //printf("npassok[%d]=%f tm_mean=%f tm_sig=%f\n",i,npassok[i],tm_mean,tm_sig);
415  putwidthVal(i, tm_mean);
416 }
417 
418 void TShapeAnalysis::putalphaVal(int n, double val) { alpha_val[n] = val; }
419 
420 void TShapeAnalysis::putchi2Val(int n, double val) { chi2_val[n] = val; }
421 void TShapeAnalysis::putbetaVal(int n, double val) { beta_val[n] = val; }
422 
423 void TShapeAnalysis::putwidthVal(int n, double val) { width_val[n] = val; }
424 
425 void TShapeAnalysis::putflagVal(int n, int val) { flag_val[n] = val; }
426 
427 void TShapeAnalysis::putalphaInit(int n, double val) { alpha_init[n] = val; }
428 
429 void TShapeAnalysis::putchi2Init(int n, double val) { chi2_init[n] = val; }
430 void TShapeAnalysis::putbetaInit(int n, double val) { beta_init[n] = val; }
431 
432 void TShapeAnalysis::putwidthInit(int n, double val) { width_init[n] = val; }
433 
434 void TShapeAnalysis::putetaInit(int n, int val) { eta_init[n] = val; }
435 
436 void TShapeAnalysis::putphiInit(int n, int val) { phi_init[n] = val; }
437 
438 void TShapeAnalysis::putflagInit(int n, int val) { flag_init[n] = val; }
439 std::vector<double> TShapeAnalysis::getVals(int n) {
440  std::vector<double> v;
441 
442  v.push_back(alpha_val[n]);
443  v.push_back(beta_val[n]);
444  v.push_back(width_val[n]);
445  v.push_back(chi2_val[n]);
446  v.push_back(flag_val[n]);
447 
448  return v;
449 }
450 std::vector<double> TShapeAnalysis::getInitVals(int n) {
451  std::vector<double> v;
452 
453  v.push_back(alpha_init[n]);
454  v.push_back(beta_init[n]);
455  v.push_back(width_init[n]);
456  v.push_back(chi2_init[n]);
457  v.push_back(flag_init[n]);
458 
459  return v;
460 }
461 
462 void TShapeAnalysis::printshapeData(int gRunNumber) {
463  FILE *fd;
464  int nev;
465  sprintf(filename, "runABW%d.pedestal", gRunNumber);
466  fd = fopen(filename, "w");
467  if (fd == nullptr)
468  printf("Error while opening file : %s\n", filename);
469 
470  for (int i = 0; i < nchsel; i++) {
471  if (index[i] >= 0) {
472  nev = (int)npassok[i];
473  double trise = alpha_val[i] * beta_val[i];
474  fprintf(fd,
475  "%d %d 1 %ld %ld %f %f %f %f\n",
476  index[i],
477  nev,
478  timestart,
479  timestop,
480  alpha_val[i],
481  beta_val[i],
482  trise,
483  width_val[i]);
484  }
485  }
486  int iret = fclose(fd);
487  printf(" Closing file : %d\n", iret);
488 }
TShapeAnalysis::putflagInit
void putflagInit(int, int)
Definition: TShapeAnalysis.cc:438
ApeEstimator_cff.width
width
Definition: ApeEstimator_cff.py:24
init
int init
Definition: HydjetWrapper.h:64
mps_fire.i
i
Definition: mps_fire.py:428
TShapeAnalysis::getDateStart
void getDateStart()
Definition: TShapeAnalysis.cc:163
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
TShapeAnalysis::putalphaVal
void putalphaVal(int, double)
Definition: TShapeAnalysis.cc:418
TShapeAnalysis::putbetaVal
void putbetaVal(int, double)
Definition: TShapeAnalysis.cc:421
zMuMuMuonUserData.alpha
alpha
zGenParticlesMatch = cms.InputTag(""),
Definition: zMuMuMuonUserData.py:9
trackAssociatorByChi2_cfi.chi2cut
chi2cut
Definition: trackAssociatorByChi2_cfi.py:4
zMuMuMuonUserData.beta
beta
Definition: zMuMuMuonUserData.py:10
TShapeAnalysis::putphiInit
void putphiInit(int, int)
Definition: TShapeAnalysis.cc:436
TShapeAnalysis::putAllVals
void putAllVals(int, double *, int, int)
Definition: TShapeAnalysis.cc:185
fNchsel
#define fNchsel
Definition: TShapeAnalysis.h:8
TShapeAnalysis::putDateStop
void putDateStop(long int)
Definition: TShapeAnalysis.cc:161
hgcalTowerProducer_cfi.tower
tower
Definition: hgcalTowerProducer_cfi.py:4
TShapeAnalysis::putDateStart
void putDateStart(long int)
Definition: TShapeAnalysis.cc:159
cms::cuda::assert
assert(be >=bs)
TShapeAnalysis::assignChannel
void assignChannel(int, int)
Definition: TShapeAnalysis.cc:152
TShapeAnalysis::putflagVal
void putflagVal(int, int)
Definition: TShapeAnalysis.cc:425
findQualityFiles.v
v
Definition: findQualityFiles.py:179
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
TShapeAnalysis::putbetaInit
void putbetaInit(int, double)
Definition: TShapeAnalysis.cc:430
TShapeAnalysis::getVals
std::vector< double > getVals(int)
Definition: TShapeAnalysis.cc:439
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
contentValuesCheck.ss
ss
Definition: contentValuesCheck.py:33
ztee.fd
fd
Definition: ztee.py:136
debug
#define debug
Definition: HDRShower.cc:19
TShapeAnalysis::computetmaxVal
void computetmaxVal(int, double *)
Definition: TShapeAnalysis.cc:389
TShapeAnalysis::putetaInit
void putetaInit(int, int)
Definition: TShapeAnalysis.cc:434
TShapeAnalysis::printshapeData
void printshapeData(int)
Definition: TShapeAnalysis.cc:462
corrVsCorr.filename
filename
Definition: corrVsCorr.py:123
dqmdumpme.k
k
Definition: dqmdumpme.py:60
TShapeAnalysis::getInitVals
std::vector< double > getInitVals(int)
Definition: TShapeAnalysis.cc:450
TShapeAnalysis::~TShapeAnalysis
~TShapeAnalysis() override
Definition: TShapeAnalysis.cc:34
TShapeAnalysis.h
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
TShapeAnalysis::set_presample
void set_presample(int)
Definition: TShapeAnalysis.cc:145
TShapeAnalysis::putalphaInit
void putalphaInit(int, double)
Definition: TShapeAnalysis.cc:427
TShapeAnalysis::putchi2Val
void putchi2Val(int, double)
Definition: TShapeAnalysis.cc:420
hgcalDigitizer_cfi.noise
noise
Definition: hgcalDigitizer_cfi.py:155
TShapeAnalysis::computeShape
void computeShape(std::string namefile, TTree *)
Definition: TShapeAnalysis.cc:205
TShapeAnalysis::putwidthVal
void putwidthVal(int, double)
Definition: TShapeAnalysis.cc:423
createfilelist.int
int
Definition: createfilelist.py:10
TShapeAnalysis::putwidthInit
void putwidthInit(int, double)
Definition: TShapeAnalysis.cc:432
nevt
int nevt
Definition: ReggeGribovPartonMCHadronizer.h:66
heppy_batch.val
val
Definition: heppy_batch.py:351
std
Definition: JetResolutionObject.h:76
TShapeAnalysis::putchi2Init
void putchi2Init(int, double)
Definition: TShapeAnalysis.cc:429
TFParams::set_const
void set_const(int, int, int, double, double, int)
Definition: TFParams.cc:509
TShapeAnalysis::set_nch
void set_nch(int)
Definition: TShapeAnalysis.cc:146
TFParams
Definition: TFParams.h:45
ecaldqm::towerID
EcalLogicID towerID(EcalElectronicsId const &)
Definition: LogicIDTranslation.cc:37
TShapeAnalysis::set_const
void set_const(int, int, int, int, int, double, double)
Definition: TShapeAnalysis.cc:135
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
TFParams.h
ntuplemaker.time
time
Definition: ntuplemaker.py:310
TShapeAnalysis::init
void init(double, double, double, double)
Definition: TShapeAnalysis.cc:36
TFParams::fitpj
double fitpj(double **, double *, double **, double noise_val, int debug)
Definition: TFParams.cc:34
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
TShapeAnalysis::TShapeAnalysis
TShapeAnalysis(double, double, double, double)
Definition: TShapeAnalysis.cc:25
TShapeAnalysis::getDateStop
void getDateStop()
Definition: TShapeAnalysis.cc:169
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:116
Generated for CMSSW Reference Manual by doxygen 1.8.16