CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 <math.h>
13 #include <time.h>
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 {
27  init(alpha0, beta0, width0, chi20);
28 }
29 
30 TShapeAnalysis::TShapeAnalysis(TTree *tAB, double alpha0, double beta0, double width0, double chi20)
31 {
32  init(tAB, alpha0, beta0, width0, chi20);
33 }
34 
35 // Destructor
36 TShapeAnalysis::~TShapeAnalysis()
37 {
38 }
39 
40 void TShapeAnalysis::init(double alpha0, double beta0, double width0, double chi20)
41 {
42  tABinit=NULL;
43  nchsel=fNchsel;
44  for(int cris=0;cris<fNchsel;cris++){
45 
46  index[cris]=-1;
47  npass[cris]=0;
48  npassok[cris]=0.;
49 
50  alpha_val[cris]=alpha0;
51  beta_val[cris]=beta0;
52  width_val[cris]=width0;
53  chi2_val[cris]=chi20;
54  flag_val[cris]=0;
55 
56  alpha_init[cris]=alpha0;
57  beta_init[cris]=beta0;
58  width_init[cris]=width0;
59  chi2_init[cris]=chi20;
60  flag_init[cris]=0;
61 
62  phi_init[cris]=0;
63  eta_init[cris]=0;
64  side_init[cris]=0;
65  dcc_init[cris]=0;
66  tower_init[cris]=0;
67  ch_init[cris]=0;
68  assignChannel(cris,cris);
69 
70  }
71 }
72 
73 void TShapeAnalysis::init(TTree *tAB, double alpha0, double beta0, double width0, double chi20)
74 {
75  init( alpha0, beta0, width0, chi20 );
76 
77  tABinit=tAB->CloneTree();
78 
79  // Declaration of leaf types
80  Int_t sidei;
81  Int_t iphii;
82  Int_t ietai;
83  Int_t dccIDi;
84  Int_t towerIDi;
85  Int_t channelIDi;
86  Double_t alphai;
87  Double_t betai;
88  Double_t widthi;
89  Double_t chi2i;
90  Int_t flagi;
91 
92  // List of branches
93  TBranch *b_iphi;
94  TBranch *b_ieta;
95  TBranch *b_side;
96  TBranch *b_dccID;
97  TBranch *b_towerID;
98  TBranch *b_channelID;
99  TBranch *b_alpha;
100  TBranch *b_beta;
101  TBranch *b_width;
102  TBranch *b_chi2;
103  TBranch *b_flag;
104 
105  if(tABinit){
106 
107  tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
108  tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
109  tABinit->SetBranchAddress("side", &sidei, &b_side);
110  tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
111  tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
112  tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
113  tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
114  tABinit->SetBranchAddress("beta", &betai, &b_beta);
115  tABinit->SetBranchAddress("width", &widthi, &b_width);
116  tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
117  tABinit->SetBranchAddress("flag", &flagi, &b_flag);
118 
119  nchsel=tABinit->GetEntries();
120  assert(nchsel<=fNchsel);
121 
122  for(int cris=0;cris<nchsel;cris++){
123 
124  tABinit->GetEntry(cris);
125 
126  // std::cout<< "Loop 1 "<< cris<<" "<<alphai<< std::endl;
127 
128  putalphaVal(cris,alphai);
129  putchi2Val(cris,chi2i);
130  putbetaVal(cris,betai);
131  putwidthVal(cris,widthi);
132  putflagVal(cris,flagi);
133 
134  putalphaInit(cris,alphai);
135  putchi2Init(cris,chi2i);
136  putbetaInit(cris,betai);
137  putwidthInit(cris,widthi);
138  putflagInit(cris,flagi);
139  putetaInit(cris,ietai);
140  putphiInit(cris,iphii);
141 
142  }
143  }
144 }
145 
146 void TShapeAnalysis::set_const(int ns, int ns1, int ns2, int ps, int nevtmax, double noise_val, double chi2_cut)
147 {
148  nsamplecristal= ns;
149  presample=ps;
150  sampbmax= ns1;
151  sampamax= ns2;
152  nevt= nevtmax;
153  noise= noise_val;
154  chi2cut=chi2_cut;
155 }
156 
157 
158 void TShapeAnalysis::set_presample(int ps)
159 {
160  presample=ps;
161 }
162 void TShapeAnalysis::set_nch(int nch){
163 
164  assert (nch<=fNchsel);
165  if(tABinit) assert(nch==nchsel);
166  nchsel=nch;
167 
168 }
169 void TShapeAnalysis::assignChannel(int n, int ch)
170 {
171  if(n >= nchsel)
172  printf(" number of channels exceed maximum allowed\n");
173 
174  index[n]=ch;
175 }
176 
177 void TShapeAnalysis::putDateStart(long int timecur)
178 {
179  timestart=timecur;
180 }
181 
182 void TShapeAnalysis::putDateStop(long int timecur)
183 {
184  timestop=timecur;
185 }
186 
187 void TShapeAnalysis::getDateStart()
188 {
189  time_t t,timecur;
190  timecur= time(&t);
191  timestart= ((long int) timecur);
192 }
193 
194 void TShapeAnalysis::getDateStop()
195 {
196  time_t t,timecur;
197  timecur= time(&t);
198  timestop= ((long int) timecur);
199 }
200 
201 void TShapeAnalysis::putAllVals(int ch, double* sampl, int ieta, int iphi, int dcc, int side, int tower, int chid)
202 {
203  dcc_init[ch]=dcc;
204  tower_init[ch]=side;
205  ch_init[ch]=chid;
206  side_init[ch]=side;
207  eta_init[ch]=ieta;
208  phi_init[ch]=iphi;
209  putAllVals(ch, sampl, ieta, iphi);
210 
211 }
212 
213 void TShapeAnalysis::putAllVals(int ch, double* sampl, int ieta, int iphi)
214 {
215 
216  int i,k;
217  int n=-1;
218  for(i=0;i<nchsel;i++)
219  if(index[i] == ch) n=i;
220 
221  if(n >= 0) {
222  if(npass[n] < nevt) {
223 
224  for(k=0;k<nsamplecristal;k++) {
225  rawsglu[n][npass[n]][k] = sampl[k];
226  }
227 
228  npass[n]++;
229  }
230  } else {
231  printf("no index found for ch=%d\n",ch);
232  }
233 }
234 
235 void TShapeAnalysis::computeShape(string namefile, TTree *tAB)
236 {
237 
238  double tm_atmax[200];
239  double parout[3];
240 
241  double chi2_all=0.;
242 
243  double **dbi ;
244  dbi = new double *[200];
245  for(int k=0;k<200;k++) dbi[k] = new double[2];
246 
247  double **signalu ;
248  signalu = new double *[200];
249  for(int k=0 ;k<200;k++) signalu[k] = new double[10];
250 
251  TFParams *pjf = new TFParams() ;
252 
253  for(int i=0;i<nchsel;i++) {
254 
255  if(index[i] >= 0) {
256 
257  if(npass[i] <= 10) {
258 
259  putalphaVal(i, alpha_init[i]);
260  putbetaVal(i, beta_init[i]);
261  putwidthVal(i,width_init[i]);
262  putchi2Val(i,chi2_init[i]);
263  putflagVal(i,0);
264 
265  } else {
266 
267  pjf->set_const(nsamplecristal, sampbmax, sampamax, alpha_init[i], beta_init[i], npass[i]);
268 
269  for(int pass=0;pass<npass[i];pass++){
270 
271  double ped=0;
272  for(int k=0;k<presample;k++) {
273  ped+= rawsglu[i][pass][k];
274  }
275  ped/=double(presample);
276 
277  for(int k=0;k<nsamplecristal;k++) {
278  signalu[pass][k]= rawsglu[i][pass][k]-ped;
279  }
280  }
281 
282  int debug=0;
283  chi2_all= pjf->fitpj(signalu,&parout[0],dbi,noise, debug);
284 
285  if(parout[0]>=0.0 && parout[1]>=0.0 && chi2_all<=chi2cut && chi2_all>0.0){
286 
287  putalphaVal(i,parout[0]);
288  putbetaVal(i,parout[1]);
289  putchi2Val(i,chi2_all);
290  putflagVal(i,1);
291 
292  }else{
293 
294  putalphaVal(i,alpha_init[i]);
295  putbetaVal(i,beta_init[i]);
296  putwidthVal(i,width_init[i]);
297  putchi2Val(i,chi2_init[i]);
298  putflagVal(i,0);
299 
300  }
301 
302  for(int kj=0;kj<npass[i];kj++) { // last event wrong here
303  tm_atmax[kj]= dbi[kj][1];
304  }
305  computetmaxVal(i,&tm_atmax[0]);
306 
307  }
308  }
309  }
310 
311  if(tAB) tABinit=tAB->CloneTree();
312 
313  // Declaration of leaf types
314  Int_t sidei;
315  Int_t iphii;
316  Int_t ietai;
317  Int_t dccIDi;
318  Int_t towerIDi;
319  Int_t channelIDi;
320  Double_t alphai;
321  Double_t betai;
322  Double_t widthi;
323  Double_t chi2i;
324  Int_t flagi;
325 
326  // List of branches
327  TBranch *b_iphi;
328  TBranch *b_ieta;
329  TBranch *b_side;
330  TBranch *b_dccID;
331  TBranch *b_towerID;
332  TBranch *b_channelID;
333  TBranch *b_alpha;
334  TBranch *b_beta;
335  TBranch *b_width;
336  TBranch *b_chi2;
337  TBranch *b_flag;
338 
339 
340  if(tABinit){
341  tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
342  tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
343  tABinit->SetBranchAddress("side", &sidei, &b_side);
344  tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
345  tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
346  tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
347  tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
348  tABinit->SetBranchAddress("beta", &betai, &b_beta);
349  tABinit->SetBranchAddress("width", &widthi, &b_width);
350  tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
351  tABinit->SetBranchAddress("flag", &flagi, &b_flag);
352  }
353 
354  TFile *fABout = new TFile(namefile.c_str(), "RECREATE");
355  tABout=new TTree("ABCol0","ABCol0");
356 
357  // Declaration of leaf types
358  Int_t side;
359  Int_t iphi;
360  Int_t ieta;
361  Int_t dccID;
362  Int_t towerID;
363  Int_t channelID;
364  Double_t alpha;
365  Double_t beta;
366  Double_t width;
367  Double_t chi2;
368  Int_t flag;
369 
370  tABout->Branch( "iphi", &iphi, "iphi/I" );
371  tABout->Branch( "ieta", &ieta, "ieta/I" );
372  tABout->Branch( "side", &side, "side/I" );
373  tABout->Branch( "dccID", &dccID, "dccID/I" );
374  tABout->Branch( "towerID", &towerID, "towerID/I" );
375  tABout->Branch( "channelID", &channelID, "channelID/I" );
376  tABout->Branch( "alpha", &alpha, "alpha/D" );
377  tABout->Branch( "beta", &beta, "beta/D" );
378  tABout->Branch( "width", &width, "width/D" );
379  tABout->Branch( "chi2", &chi2, "chi2/D" );
380  tABout->Branch( "flag", &flag, "flag/I" );
381 
382  tABout->SetBranchAddress( "ieta", &ieta );
383  tABout->SetBranchAddress( "iphi", &iphi );
384  tABout->SetBranchAddress( "side", &side );
385  tABout->SetBranchAddress( "dccID", &dccID );
386  tABout->SetBranchAddress( "towerID", &towerID );
387  tABout->SetBranchAddress( "channelID", &channelID );
388  tABout->SetBranchAddress( "alpha", &alpha );
389  tABout->SetBranchAddress( "beta", &beta );
390  tABout->SetBranchAddress( "width", &width );
391  tABout->SetBranchAddress( "chi2", &chi2 );
392  tABout->SetBranchAddress( "flag", &flag );
393 
394  for(int i=0;i<nchsel;i++) {
395 
396  if(tABinit){
397 
398  tABinit->GetEntry(i);
399  iphi=iphii;
400  ieta=ietai;
401  side=sidei;
402  dccID=dccIDi;
403  towerID=towerIDi;
404  channelID=channelIDi;
405 
406  }else{
407 
408  iphi=phi_init[i];
409  ieta=eta_init[i];
410  side=side_init[i];
411  dccID=dcc_init[i];
412  towerID=tower_init[i];
413  channelID=ch_init[i];
414 
415  }
416 
417  alpha=alpha_val[i];
418  beta=beta_val[i];
419  width=width_val[i];
420  chi2=chi2_val[i];
421  flag=flag_val[i];
422 
423  tABout->Fill();
424  }
425 
426 
427  tABout->Write();
428  fABout->Close();
429 
430  delete pjf;
431 }
432 
433 void TShapeAnalysis::computetmaxVal(int i, double* tm_val)
434 {
435  double tm_mean=0; //double tm_sig=0;
436 
437  double tm=0.; double sigtm=0.;
438  for(int k=0;k<npass[i]-1;k++) {
439  if(1. < tm_val[k] && tm_val[k] < 10.) {
440  npassok[i]++;
441  tm+= tm_val[k];
442  sigtm+= tm_val[k]*tm_val[k];
443  }
444  }
445  if(npassok[i] <= 0) {
446  tm_mean=0.; //tm_sig=0.;
447  } else {
448  for(int k=0;k<npass[i]-1;k++) {
449  if(1. < tm_val[k] && tm_val[k] < 10.) {
450  double ss= (sigtm/npassok[i]-tm/npassok[i]*tm/npassok[i]);
451  if(ss < 0.) ss=0.;
452  //tm_sig=sqrt(ss);
453  tm_mean= tm/npassok[i];
454  }
455  }
456  }
457  //printf("npassok[%d]=%f tm_mean=%f tm_sig=%f\n",i,npassok[i],tm_mean,tm_sig);
458  putwidthVal(i,tm_mean);
459 
460 }
461 
462 void TShapeAnalysis::putalphaVal(int n, double val)
463 {
464  alpha_val[n]= val;
465 }
466 
467 void TShapeAnalysis::putchi2Val(int n, double val)
468 {
469  chi2_val[n]= val;
470 }
471 void TShapeAnalysis::putbetaVal(int n, double val)
472 {
473  beta_val[n]= val;
474 }
475 
476 void TShapeAnalysis::putwidthVal(int n, double val)
477 {
478  width_val[n]= val;
479 }
480 
481 void TShapeAnalysis::putflagVal(int n, int val)
482 {
483  flag_val[n]= val;
484 }
485 
486 void TShapeAnalysis::putalphaInit(int n, double val)
487 {
488  alpha_init[n]= val;
489 }
490 
491 void TShapeAnalysis::putchi2Init(int n, double val)
492 {
493  chi2_init[n]= val;
494 }
495 void TShapeAnalysis::putbetaInit(int n, double val)
496 {
497  beta_init[n]= val;
498 }
499 
500 void TShapeAnalysis::putwidthInit(int n, double val)
501 {
502  width_init[n]= val;
503 }
504 
505 void TShapeAnalysis::putetaInit(int n, int val )
506 {
507  eta_init[n]= val;
508 }
509 
510 void TShapeAnalysis::putphiInit(int n, int val)
511 {
512  phi_init[n]= val;
513 }
514 
515 void TShapeAnalysis::putflagInit(int n, int val)
516 {
517  flag_init[n]= val;
518 }
519 std::vector<double> TShapeAnalysis::getVals(int n)
520 {
521 
522  std::vector<double> v;
523 
524  v.push_back(alpha_val[n]);
525  v.push_back(beta_val[n]);
526  v.push_back(width_val[n]);
527  v.push_back(chi2_val[n]);
528  v.push_back(flag_val[n]);
529 
530  return v;
531 }
532 std::vector<double> TShapeAnalysis::getInitVals(int n)
533 {
534 
535  std::vector<double> v;
536 
537  v.push_back(alpha_init[n]);
538  v.push_back(beta_init[n]);
539  v.push_back(width_init[n]);
540  v.push_back(chi2_init[n]);
541  v.push_back(flag_init[n]);
542 
543  return v;
544 }
545 
546 void TShapeAnalysis::printshapeData(int gRunNumber)
547 {
548  FILE *fd;
549  int nev;
550  sprintf(filename,"runABW%d.pedestal",gRunNumber);
551  fd = fopen(filename, "w");
552  if(fd == NULL) printf("Error while opening file : %s\n",filename);
553 
554  for(int i=0; i<nchsel;i++) {
555  if(index[i] >= 0) {
556  nev= (int) npassok[i];
557  double trise= alpha_val[i]*beta_val[i];
558  fprintf( fd, "%d %d 1 %ld %ld %f %f %f %f\n",
559  index[i],nev,timestart,timestop,alpha_val[i],beta_val[i],trise,width_val[i]);
560  }
561  }
562  int iret=fclose(fd);
563  printf(" Closing file : %d\n",iret);
564 
565 }
566 
beta
const double beta
Definition: CosmicMuonParameters.h:69
TShapeAnalysis::putbetaInit
void putbetaInit(int, double)
Definition: TShapeAnalysis.cc:495
TShapeAnalysis::putchi2Val
void putchi2Val(int, double)
Definition: TShapeAnalysis.cc:467
i
int i
Definition: DBlmapReader.cc:9
alpha
float alpha
Definition: AMPTWrapper.h:95
TShapeAnalysis::putwidthVal
void putwidthVal(int, double)
Definition: TShapeAnalysis.cc:476
HLT_25ns14e33_v1_cff.noise
tuple noise
Definition: HLT_25ns14e33_v1_cff.py:2166
TShapeAnalysis::putflagVal
void putflagVal(int, int)
Definition: TShapeAnalysis.cc:481
create_public_lumi_plots.width
int width
Definition: create_public_lumi_plots.py:1008
TShapeAnalysis::~TShapeAnalysis
virtual ~TShapeAnalysis()
Definition: TShapeAnalysis.cc:36
TShapeAnalysis::putDateStart
void putDateStart(long int)
Definition: TShapeAnalysis.cc:177
assert
assert(m_qm.get())
init
int init
Definition: HydjetWrapper.h:62
TShapeAnalysis::init
void init(double, double, double, double)
Definition: TShapeAnalysis.cc:40
NULL
#define NULL
Definition: scimark2.h:8
TShapeAnalysis::computeShape
void computeShape(std::string namefile, TTree *)
Definition: TShapeAnalysis.cc:235
TShapeAnalysis::putAllVals
void putAllVals(int, double *, int, int)
Definition: TShapeAnalysis.cc:213
TShapeAnalysis::putbetaVal
void putbetaVal(int, double)
Definition: TShapeAnalysis.cc:471
fNchsel
#define fNchsel
Definition: TShapeAnalysis.h:8
TShapeAnalysis::putDateStop
void putDateStop(long int)
Definition: TShapeAnalysis.cc:182
TShapeAnalysis::printshapeData
void printshapeData(int)
Definition: TShapeAnalysis.cc:546
TShapeAnalysis::TShapeAnalysis
TShapeAnalysis(double, double, double, double)
Definition: TShapeAnalysis.cc:25
TShapeAnalysis::assignChannel
void assignChannel(int, int)
Definition: TShapeAnalysis.cc:169
TShapeAnalysis::putetaInit
void putetaInit(int, int)
Definition: TShapeAnalysis.cc:505
TFParams::set_const
void set_const(int, int, int, double, double, int)
Definition: TFParams.cc:552
TShapeAnalysis::putflagInit
void putflagInit(int, int)
Definition: TShapeAnalysis.cc:515
TShapeAnalysis::getInitVals
std::vector< double > getInitVals(int)
Definition: TShapeAnalysis.cc:532
TShapeAnalysis::putalphaVal
void putalphaVal(int, double)
Definition: TShapeAnalysis.cc:462
debug
#define debug
Definition: HDRShower.cc:19
TShapeAnalysis::putchi2Init
void putchi2Init(int, double)
Definition: TShapeAnalysis.cc:491
TShapeAnalysis::putalphaInit
void putalphaInit(int, double)
Definition: TShapeAnalysis.cc:486
TShapeAnalysis::putwidthInit
void putwidthInit(int, double)
Definition: TShapeAnalysis.cc:500
TShapeAnalysis::getVals
std::vector< double > getVals(int)
Definition: TShapeAnalysis.cc:519
TShapeAnalysis::set_const
void set_const(int, int, int, int, int, double, double)
Definition: TShapeAnalysis.cc:146
TShapeAnalysis::set_presample
void set_presample(int)
Definition: TShapeAnalysis.cc:158
lut2db_cfg.filename
tuple filename
Definition: lut2db_cfg.py:20
TShapeAnalysis::computetmaxVal
void computetmaxVal(int, double *)
Definition: TShapeAnalysis.cc:433
ecaldqm::towerID
EcalLogicID towerID(EcalElectronicsId const &)
Definition: LogicIDTranslation.cc:48
relval_steps.k
list k
Definition: relval_steps.py:1469
TFParams::fitpj
double fitpj(double **, double *, double **, double noise_val, int debug)
Definition: TFParams.cc:38
TShapeAnalysis::putphiInit
void putphiInit(int, int)
Definition: TShapeAnalysis.cc:510
Generated for CMSSW Reference Manual by doxygen 1.8.5