ATLAS Offline Software
Calibrator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /********************************************************************
6 
7  NAME: Calibrator.cxx
8 PACKAGE: TRTCalibTools
9 
10 AUTHORS: Johan Lundquist
11 CREATED: 27-03-2009
12 
13 PURPOSE: Class for calibrating a TRT sub-level
14 
15 ********************************************************************/
16 
17 #include "Calibrator.h"
18 #include "CxxUtils/checker_macros.h"
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <TH1D.h>
22 #include <TDirectory.h>
23 #include <TF1.h>
24 #include <TStyle.h>
25 #include <TNtuple.h>
26 #include <TGraphErrors.h>
27 #include <iostream>
28 #include <iomanip>
29 #include <fstream>
30 #include <cmath>
31 
32 
33 caldata::caldata(){
34  res = -999;
35  resMean = -999;
36  reserr = -999;
37  tres = -999;
38  tresMean = -999;
39  t0 = -999;
40  t0err = -999;
41  reft0 = -999;
42  t0off = -999;
43  rtt0 = -999;
44  nhits = -999;
45  x = -999;
46  y = -999;
47  z = -999;
48  oldt02 = -999;
49  sumt0 = -999;
50  sumx = -999;
51  sumy = -999;
52  sumz = -999;
53  calflag = -999;
54  rtflag = -999;
55  t0flag = -999;
56  det = -999;
57  lay = -999;
58  mod = -999;
59  brd = -999;
60  chp = -999;
61  stl = -999;
62  stw = -999;
63  sid = -999;
64  ntres = -999;
65  nrt = -999;
66  nres = -999;
67  t0fittype = -999;
68  rtgraph = nullptr;
69 }
70 
71 caldata::~caldata(){
72 }
73 
74 caldata::caldata(bool makehist, int nbinst, int nbinsr){
75  res = -999;
76  resMean = -999;
77  reserr = -999;
78  tres = -999;
79  tresMean = -999;
80  t0 = -999;
81  t0err = -999;
82  reft0 = -999;
83  t0off = -999;
84  rtt0 = -999;
85  nhits = -999;
86  x = -999;
87  y = -999;
88  z = -999;
89  oldt02 = -999;
90  sumt0 = -999;
91  sumx = -999;
92  sumy = -999;
93  sumz = -999;
94  calflag = -999;
95  rtflag = -999;
96  t0flag = -999;
97  det = -999;
98  lay = -999;
99  mod = -999;
100  brd = -999;
101  chp = -999;
102  stl = -999;
103  stw = -999;
104  sid = -999;
105  ntres = -999;
106  nrt = -999;
107  nres = -999;
108  t0fittype= -999;
109  rtgraph = nullptr;
110 
111  if (makehist) {
112  m_treshist.resize (100);
113  reshist.resize (100);
114  rthist.resize (nbinsr*nbinst+200);
115  }
116 }
117 
118 RtGraph::RtGraph(TH2F* rtHist, int binvar, const char* binlabel, bool pflag, TDirectory* dir){
119 
120  int fitresult;
121  npoints = binvar==0 ? rtHist->GetNbinsX() : rtHist->GetNbinsY() ;
122 
123 
124  m_mean = -20.0;
125  m_t = -20.0;
126  m_d = -20.0;
127  m_et = -20.0;
128  m_ed = -20.0;
129  // m_ff = nullptr ;
130 
131 
132  TH1D** hslizes = new TH1D*[npoints];
133  m_btype = new bintype[npoints];
134  m_tv = new float[npoints];
135  m_dv = new float[npoints];
136  m_etv = new float[npoints];
137  m_edv = new float[npoints];
138  m_rightsig = new float[npoints];
139  m_leftsig = new float[npoints];
140  m_leftval = new float[npoints];
141  m_rightval = new float[npoints];
142  m_maxbin = new int[npoints];
143  m_maxval = new float[npoints];
144 
145  m_ipoint=0;
146  mintime=999;
147  m_mindistance=0;
148 
149  if (pflag){
150  TDirectory* binhistdir = dir->mkdir("binhist");
151  binhistdir->cd();
152  std::cout << " Calibrator::RtGraph: " << " Directory " << std::string(binhistdir->GetPath()) << " binlabel " << std::string(binlabel) << std::endl;
153  }
154 
155  std::cout << " Calibrator::RtGraph: " << " number of points " << npoints << " binvar " << binvar << std::endl;
156  // check and classify the bin-histograms
157  for (int i=0;i<npoints;i++) {
158 
159  m_chname = std::string(Form("%s slize_%i",binlabel,i));
160  m_chtit = binvar==0 ? std::string(Form("bin %i : %4.2f < %s < %4.2f [ns]",i,rtHist->GetXaxis()->GetBinLowEdge(i+1),binlabel,rtHist->GetXaxis()->GetBinUpEdge(i+1))) :
161  std::string(Form("bin %i : %4.2f < %s < %4.2f [mm]",i,rtHist->GetYaxis()->GetBinLowEdge(i+1),binlabel,rtHist->GetYaxis()->GetBinUpEdge(i+1))) ;
162  hslizes[i] = binvar==0 ? rtHist->ProjectionY(m_chname.data(),i+1,i+1) :
163  rtHist->ProjectionX(m_chname.data(),i+1,i+1);
164  hslizes[i]->SetTitle(m_chtit.data());
165  if (!pflag){
166  hslizes[i]->SetDirectory(nullptr);
167  }
168 
169  m_maxbin[i]=1; m_maxval[i]=0;
170  for (int j=1;j<=hslizes[i]->GetNbinsX();j++) {
171  float val=hslizes[i]->GetBinContent(j);
172  if (val>m_maxval[i]){
173  m_maxval[i]=val;
174  m_maxbin[i]=j;
175  m_leftval[i] = hslizes[i]->GetBinContent(j-1);
176  m_rightval[i] = hslizes[i]->GetBinContent(j+1);
177  }
178  }
179 
180  if(m_maxval[i]==0) m_maxval[i]=1.0e-4;
181  m_rightsig[i] = m_rightval[i]==0 ? 100 : std::sqrt(m_maxval[i])/m_maxval[i] + std::sqrt(m_rightval[i])/m_rightval[i];
182  m_leftsig[i] = m_leftval[i]==0 ? 100 : std::sqrt(m_maxval[i])/m_maxval[i] + std::sqrt(m_leftval[i])/m_leftval[i];
183 
184  m_btype[i]=EMPTY;
185  if(true){
186  m_btype[i]=GOOD;
187  m_btype[i-1]=GOOD;
188  }
189  else{
190  if (hslizes[i]->GetEntries()>0){
191  if (i<(float)hslizes[i]->GetNbinsX()/2) m_btype[i]=LOW;
192  else m_btype[i]=HIGH;
193  }
194  }
195  int rtEntries = rtHist->GetEntries();
196  if(rtEntries==0) rtEntries=1;
197  printf("%s ... %8f %8i %4i %8f %8f %4i\n", m_chtit.data(), (float)hslizes[i]->GetEntries()/(float)rtEntries, (int)hslizes[i]->GetEntries(), m_maxbin[i], m_leftsig[i], m_rightsig[i], m_btype[i]);
198  }
199 
200  float frmin=0,frmax=0;
201  for (int i=0;i<npoints;i++) {
202 
203  // prepare initial fit parameters
204  m_mean=hslizes[i]->GetBinCenter(m_maxbin[i]);
205  if (m_btype[i]==LOW) {frmin = hslizes[i]->GetBinCenter(1); frmax = hslizes[i]->GetBinCenter(m_maxbin[i]+4);}
206  if (m_btype[i]==GOOD) {frmin = hslizes[i]->GetBinCenter(m_maxbin[i]-3); frmax = hslizes[i]->GetBinCenter(m_maxbin[i]+3);}
207  if (m_btype[i]==HIGH) {frmin = hslizes[i]->GetBinCenter(m_maxbin[i]-3); frmax = hslizes[i]->GetBinCenter(m_maxbin[i]+4);}
208 
209  std::unique_ptr<TF1> ff(new TF1("dtfitfunc","gaus"));
210 
211  ff->SetRange(frmin,frmax);
212 
213  m_t = binvar==0 ? rtHist->GetXaxis()->GetBinCenter(i+1) : rtHist->GetYaxis()->GetBinCenter(i+1);
214  m_et = binvar==0 ? rtHist->GetXaxis()->GetBinWidth(1)/std::sqrt(12.0) : rtHist->GetYaxis()->GetBinWidth(1)/std::sqrt(12.0);
215  m_d = 0;
216  m_ed = 0;
217 
218  if (m_btype[i]==LOW){
219  m_d = 0;
220  m_ed = 0;
221  }
222  if (m_btype[i]==HIGH){
223  m_d = 2;
224  m_ed = 0;
225  }
226  if (m_btype[i]==EMPTY){
227  m_d = 0;
228  m_ed = 0;
229  }
230  if (m_btype[i]==GOOD && hslizes[i]->GetEntries()>3){
231  if (m_btype[i-1]==GOOD){
232  fitresult=hslizes[i]->Fit("dtfitfunc","QR");
233  ff->SetRange(m_mean-1.0*ff->GetParameter(2), m_mean+1.0*ff->GetParameter(2));
234  fitresult=hslizes[i]->Fit("dtfitfunc","QR");
235  m_d = ff->GetParameter(1);
236  m_ed = ff->GetParError(1);
237  std::cout << fitresult << " " << m_ed << std::endl;
238  }
239  else{
240  m_d = 0;
241  m_ed = 0;
242  }
243  }
244 
245 
246  if (m_btype[i]==GOOD && m_ed!=0){
247  m_tv[m_ipoint]=m_t;
248  m_dv[m_ipoint]=m_d;
249  m_etv[m_ipoint]=m_et;
250  m_edv[m_ipoint]=m_ed;
251  m_ipoint++;
252  if (binvar==0) {
253  rval.push_back(m_d);
254  tval.push_back(m_t);
255  }
256  else {
257  rval.push_back(m_t);
258  tval.push_back(m_d);
259  }
260  }
261 
262 
263  if(m_d>=m_mindistance && m_t<mintime) mintime = m_t;
264 
265  }
266 
267  rtgr = binvar==0 ? new TGraphErrors(m_ipoint,m_tv,m_dv,m_etv,m_edv) : new TGraphErrors(m_ipoint,m_dv,m_tv,m_edv,m_etv);
268  rtgr->SetTitle((std::string("binning in ") + binlabel).data());
269  trgr = binvar==0 ? new TGraphErrors(m_ipoint,m_dv,m_tv,m_edv,m_etv) : new TGraphErrors(m_ipoint,m_tv,m_dv,m_etv,m_edv);
270  trgr->SetTitle((std::string("binning in ") + binlabel).data());
271 
272  rtgr->SetMarkerStyle(1) ;
273  rtgr->SetMarkerColor(2) ;
274  rtgr->SetLineColor(2) ;
275  rtgr->SetName("rtgraph") ;
276  trgr->SetMarkerStyle(1) ;
277  trgr->SetMarkerColor(2) ;
278  trgr->SetLineColor(2) ;
279  trgr->SetName("trgraph") ;
280 
281  dir->cd();
282  delete [] hslizes;
283 }
284 
285 
286 
287 RtGraph::~RtGraph(){
288 
289  delete [] m_btype ;
290  delete [] m_tv ;
291  delete [] m_dv ;
292  delete [] m_etv ;
293  delete [] m_edv ;
294  delete [] m_rightsig ;
295  delete [] m_leftsig ;
296  delete [] m_leftval ;
297  delete [] m_rightval ;
298  delete [] m_maxbin ;
299  delete [] m_maxval ;
300 }
301 
302 
303 double pol3deg(double *x, double *par) {
304  double r = x[0];
305  double t = par[0]+r*(par[1]+(r*(par[2]+r*par[3])));
306  return t;
307 }
308 
309 double pol3deg2(double *x, double *par) {
310  double t = x[0];
311  double r = par[1]*(t-par[0]) + par[2]*(t-par[0])*(t-par[0]) + par[3]*(t-par[0])*(t-par[0])*(t-par[0]);
312  return r;
313 }
314 
315 
316 double trrel_dines(double *x, double *par){
317  double rmin0 = par[0];
318  double rho = par[1];
319  double v = par[2];
320  double t_const = par[3];
321  double r = x[0];
322  double t = t_const + 2*rho/v*std::asin(std::sqrt(rmin0*rmin0*(1-0.25*r*r)+r*r)/(2*rho));
323  return t;
324 }
325 
326 double rtrel_dines(double *x, double*par){
327  double rmin0 = par[0];
328  double rho = par[1];
329  double v = par[2];
330  double t_const = par[3];
331  double t = x[0];
332  double r_squared = (4*rho*rho*std::sin(v*(t-t_const)/(2*rho))-rmin0*rmin0)/(1-0.25*rmin0*rmin0);
333  double r=r_squared>0 ? std::sqrt(r_squared) : 0.0;
334  return r;
335 }
336 
337 
338 
339 Calibrator::Calibrator()
340  : dort(false),
341  dot0(false),
342  dores(false),
343  nort(false),
344  not0(false),
345  usebref(false),
346  bequiet(false),
347  printlog(false),
348  printt0(false),
349  printrt(false),
350  usep0(false),
351  floatp3(false),
352  useshortstraws(true),
353 
354  m_name ("None"),
355  m_rtbinning ("None"),
356  m_minrtstat (-10),
357  m_mint0stat (-10),
358  m_t0shift (-100.),
359 
360  m_nbinsr(100),
361  m_nbinst(55),
362  m_nbinstres(100),
363  m_nbinsres(100),
364  m_minr(0),
365  m_maxr(2),
366  m_mint(-5),
367  m_maxt(50),
368  m_mintres(-10),
369  m_maxtres(10),
370  m_minres(-0.6),
371  m_maxres(0.6),
372 
373  m_ntreshits(0),
374  m_nreshits(0),
375  m_nrthits(0),
376  m_nhits(0),
377 
378  m_isdines ( false)
379 {
380  selection.insert(-3);
381  level =-10;
382 }
383 
384 Calibrator::Calibrator(int lev, const std::string& nme, int mint0, int minrt, const std::string& rtr, const std::string& rtb, float t0sft)
385  : dort(false),
386  dot0(false),
387  dores(false),
388  nort(false),
389  not0(false),
390  usebref(false),
391  bequiet(false),
392  printlog(false),
393  printt0(false),
394  printrt(false),
395  usep0(false),
396  floatp3(false),
397  useshortstraws(true),
398 
399  m_name(nme),
400  m_rtbinning(rtb),
401  m_minrtstat(minrt),
402  m_mint0stat(mint0),
403  m_t0shift(t0sft),
404 
405  m_nbinsr(100),
406  m_nbinst(55),
407  m_nbinstres(100),
408  m_nbinsres(100),
409  m_minr(0),
410  m_maxr(2),
411  m_mint(-5),
412  m_maxt(50),
413  m_mintres(-10),
414  m_maxtres(10),
415  m_minres(-0.6),
416  m_maxres(0.6),
417 
418  m_ntreshits(0),
419  m_nreshits(0),
420  m_nrthits(0),
421  m_nhits(0),
422 
423  m_isdines ( rtr.find("dines")!=std::string::npos)
424 {
425  level=lev;
426  selection.insert(-3);
427 }
428 
429 Calibrator::~Calibrator(){
430 }
431 
432 int Calibrator::Simple1dHist(float min, float max, int nbins, float value){
433  if ( (value<min) || (value>max) ) return -1;
434  if ( max-min == 0 ) return -1;
435  int binno=(int)(nbins*((value-min)/(max-min)));
436  return binno;
437 }
438 
439 int Calibrator::Simple2dHist(float minx, float maxx, int nbinsx, float miny, float maxy, int nbinsy, float valuex, float valuey){
440  if ( (valuex<minx) || (valuex>maxx) || (valuey<miny) || (valuey>maxy) ) return -1;
441  if(maxx-minx == 0 || maxy-miny==0) return -1;
442  int binnox=(int)(nbinsx*((valuex-minx)/(maxx-minx)));
443  int binnoy=(int)(nbinsy*((valuey-miny)/(maxy-miny)));
444  return binnoy*nbinsx+binnox;
445 }
446 
447 float Calibrator::AccumulativeMean(float n, float oldmean, float newvalue){
448  if((int)n==0) return oldmean;
449  return oldmean*((n-1)/n)+newvalue/n;
450 }
451 
452 bool Calibrator::HasKey(const std::string &key) const {
453 
454  return data.find(key) != data.end();
455 }
456 
457 bool Calibrator::CheckSelection(int level){
458  if (selection.find(level)!=selection.end() || selection.find(-3)!=selection.end()) return true;
459  else return false;
460 }
461 
462 bool Calibrator::Skip(){
463  if (selection.find(-4)!=selection.end()) return true;
464  else return false;
465 }
466 
467 std::string Calibrator::PrintInfo(){
468  std::string yn[2]={"NO","YES"};
469  std::string info = std::string(Form("CONFIGURATION %-16s: dort=%-3s, dot0=%-3s, dores=%-3s, selection=",m_name.data(),yn[dort].data(),yn[dot0].data(),yn[dores].data()));
470  for (int isel : selection) {
471  if (isel==-3) info += std::string("all");
472  else if (isel==-4) info += std::string("none");
473  else info += std::string(Form("%2i,",isel));
474  }
475  return info;
476 }
477 
478 std::string Calibrator::PrintStat(){
479  int units = (int)data.size();
480  if(units==0) units++;
481  std::string info = std::string(Form("STATISTICS %16s: nunits=%8i, nhits=%9i, hits/unit=%11.1f", m_name.data(), units, m_nhits, (float)m_nhits/(float)units ));
482  return info;
483 }
484 
485 std::string Calibrator::GetOptString() const {
486  std::string optstr="";
487  if (dort) optstr += 'R';
488  if (dot0) optstr += 'T';
489  if (printlog) optstr +='P';
490  if (printt0) optstr += 'F';
491  if (printrt) optstr += 'G';
492  if (bequiet) optstr += 'Q';
493  if (usebref) optstr += 'B';
494  if (useshortstraws) optstr += 'S';
495  if (usep0 && dort) optstr += '0';
496  if (floatp3 && dort) optstr += '3';
497  if (!(dort || dot0 || usebref || bequiet || printlog || printt0 || printrt)) optstr += 'N';
498 
499  return optstr;
500 }
501 
502 int Calibrator::GetNDirs(){
503  return data.size();
504 }
505 
506 std::string Calibrator::GetSelString(){
507  std::string selstr="";
508  if (selection.find(-3)!=selection.end()) selstr="*";
509  else if (selection.find(-4)!=selection.end()) selstr="-";
510  else for (int isel : selection) selstr += std::string(Form("%i",isel));
511  return selstr;
512 }
513 
514 std::string Calibrator::GetBinnedRt(const std::string& key){
515 
516  int nbins=data[key].rtgraph->tval.size();
517  double maxt=data[key].rtgraph->tval[nbins-1];
518  double mint=data[key].rtgraph->tval[0];
519 
520  std::string brt = "";
521  brt += Form("%f %f %i ",mint,maxt,nbins);
522 
523  for (int ip=0; ip<nbins;ip++){
524  brt += Form("%f ",data[key].rtgraph->rval[ip]);
525  }
526  return brt;
527 }
528 
529 int Calibrator::UpdateOldConstants(){
530 
531  std::string line;
532  std::ifstream oldconstfile("calib_constants_in.txt");
533  std::string key;
534  float t0;
535 
536  if (oldconstfile.is_open())
537  {
538  while (!oldconstfile.eof() )
539  {
540  oldconstfile >> key >> t0;
541  if (data.find(key) != data.end()){
542  std::cout << "UPDATED OLD T0: " << key << " " << data[key].oldt02 << " -> " << t0 << std::endl;
543  data[key].oldt02=t0;
544  }
545  }
546  oldconstfile.close();
547  return 0;
548  }
549  else {
550  std::cout << "NO OLD T0S FILE FOUND. USING AVERAGE VALUES!" << std::endl;
551  return -1;
552  }
553 }
554 
555 float Calibrator::FitRt ATLAS_NOT_THREAD_SAFE (const std::string& key, const std::string& opt, TH2F* rtHist, TDirectory* dir){ // Global gStyle is used.
556 
557  float rtpars[4];
558 
559  //create r-m_t and m_t-r graphs
560  RtGraph* rtg = m_rtbinning.find('t')==std::string::npos ? new RtGraph(rtHist,1,std::string("abs(rtrack)").data(),!bequiet,dir) : new RtGraph(rtHist,0,std::string("t-t0").data(),!bequiet,dir);
561 
562  TF1 dtfitfunc("dtfitfunc","gaus(0)");
563 
564  gStyle->SetOptFit(1111);
565  gStyle->SetPalette(1);
566 
567  // select type of R-m_t relation
568  TF1 *trfunc,*rtfunc,*oldrtfunc;
569  if (m_isdines){
570  trfunc = new TF1("trfunc",trrel_dines,0,3,4);
571  rtfunc = new TF1("rtfunc",rtrel_dines,rtg->mintime,200,4);
572  oldrtfunc = new TF1("oldrtfunc",rtrel_dines,rtg->mintime,200,4);
573  trfunc->SetParameters(0.3, 1.0, 0.05, -3.);
574  double rmin0Limits[2]={0.0,2.0};
575  double rhoLimits[2]={0.0,10.};
576  double vLimits[2]={0.0,1.0};
577  trfunc->SetParLimits(0,rmin0Limits[0],rmin0Limits[1]);
578  trfunc->SetParLimits(1,rhoLimits[0],rhoLimits[1]);
579  trfunc->SetParLimits(2,vLimits[0],vLimits[1]);
580  rtg->trgr->SetTitle("Dines' R-t sqrt(..)");
581  }
582  else {
583  trfunc = new TF1("trfunc",pol3deg,0,3,4);
584  rtfunc = new TF1("rtfunc",pol3deg,rtg->mintime,200,4);
585  oldrtfunc = new TF1("oldrtfunc",pol3deg,rtg->mintime,200,4);
586  trfunc->SetParameters(0.0, 0.0, 0.0, 0.0);
587  rtg->trgr->SetTitle("Polynomial R-t");
588  }
589 
590  // m_t-r relation
591  trfunc->SetRange(0,2);
592  if (opt.find('3')==std::string::npos) trfunc->FixParameter(3,0);
593  trfunc->SetLineColor(4) ;
594  rtg->trgr->Fit(trfunc,"QR"); // always fit the m_t-r relation
595  if (!bequiet) rtg->trgr->Write();
596 
597 
598  // r-m_t relation
599  rtfunc->SetRange(0,45);
600  if (opt.find('3')==std::string::npos) rtfunc->FixParameter(3,0);
601  rtfunc->SetLineColor(4) ;
602  if (m_isdines) {
603  rtfunc->SetParameters(trfunc->GetParameter(0),trfunc->GetParameter(1),trfunc->GetParameter(2),trfunc->GetParameter(3));
604 
605  }
606  else { // else do a fit
607  rtfunc->SetParameters(0.000000e+00, 6.269950e-02, -3.370054e-04, -1.244642e-07);
608  if(rtg->rtgr != nullptr) {
609  rtg->rtgr->Fit(rtfunc,"QR"); //fit Rt first time
610  for (int ipnt=0; ipnt<rtg->rtgr->GetN(); ipnt++){ //calculate m_t-error
611  float deriv = rtfunc->Derivative(rtg->rtgr->GetX()[ipnt]);
612  if(fabs(deriv)<1.0e-24) deriv=0.05;
613  rtg->rtgr->SetPointError(ipnt , rtg->rtgr->GetEY()[ipnt]/deriv , rtg->rtgr->GetEY()[ipnt]);
614  }
615  rtg->rtgr->Fit(rtfunc,"R"); //fit again
616  }
617 
618  }
619  if (!bequiet and (rtg->rtgr != nullptr)) rtg->rtgr->Write();
620 
621  // old r-m_t relation
622  oldrtfunc->SetRange(0,45);
623  oldrtfunc->SetLineColor(1) ;
624  oldrtfunc->SetLineStyle(2) ;
625  oldrtfunc->SetLineWidth(1) ;
626  oldrtfunc->SetParameters(data[key].oldrtpar[0],data[key].oldrtpar[1],data[key].oldrtpar[2],data[key].oldrtpar[3]);
627 
628  if (!bequiet) oldrtfunc->Write();
629 
630  rtpars[0] = rtfunc->GetParameter(0);
631  rtpars[1] = rtfunc->GetParameter(1);
632  rtpars[2] = rtfunc->GetParameter(2);
633  rtpars[3] = rtfunc->GetParameter(3);
634 
635  //get the t0 from the tr relation
636  float tdrift = 20;
637 
638  if (!m_isdines) {
639  float rdrift = 0.0;
640  int ntries = 0;
641  float precision = 0.0001;
642  int maxtries = 500;
643  float drdt;
644  float driftradius = rtpars[0]+tdrift*(rtpars[1]+tdrift*(rtpars[2]));
645  float residual = std::abs(rdrift) - driftradius;
646  while (std::abs(residual) > precision) {
647 
648  drdt = rtpars[1]+tdrift*(2*rtpars[2]);
649  if(fabs(drdt)<1.0e-24) drdt=0.05;
650  tdrift = tdrift + residual/drdt;
651 
652  driftradius = rtpars[0]+tdrift*(rtpars[1]+tdrift*(rtpars[2]));
653  residual = rdrift - driftradius;
654 
655  ntries = ntries + 1;
656  if (ntries>maxtries){
657  break;
658  }
659  }
660  }
661 
662  if (opt.find('0')==std::string::npos) {
663  if (m_isdines){
664  data[key].rtpar[0] = rtpars[0];
665  data[key].rtpar[1] = rtpars[1];
666  data[key].rtpar[2] = rtpars[2];
667  data[key].rtpar[3] = 0;
668  } else {
669  data[key].rtpar[0] = 0;
670  data[key].rtpar[1] = rtpars[1] + 2*tdrift*rtpars[2] + 3*tdrift*tdrift*rtpars[3];
671  data[key].rtpar[2] = rtpars[2] + 3*tdrift*rtpars[3];
672  data[key].rtpar[3] = rtpars[3];
673  }
674  }
675  else {
676  data[key].rtpar[0] = rtpars[0];
677  data[key].rtpar[1] = rtpars[1];
678  data[key].rtpar[2] = rtpars[2];
679  data[key].rtpar[3] = rtpars[3];
680  }
681 
682  data[key].rtgraph=rtg;
683 
684  delete rtfunc;
685  delete oldrtfunc;
686  delete trfunc;
687 
688  return 0.0;
689 }
690 
691 
692 float Calibrator::FitTimeResidual ATLAS_NOT_THREAD_SAFE (const std::string& key, TH1F* tresHist){ // Global gStyle is used.
693 
694  float mean = tresHist->GetMean();
695  float rms = tresHist->GetRMS();
696  float val,maxval=0;
697  int maxbin=1;
698  for (int j=1; j<=tresHist->GetNbinsX()+1;j++){
699  val=tresHist->GetBinContent(j);
700  if (val>maxval){
701  maxval=val;
702  maxbin=j;
703  }
704  }
705 
706 
707  data[key].tresMean = 0;
708  float fitmean=tresHist->GetBinCenter(maxbin);
709  if (std::abs(fitmean)>5.0){
710  data[key].tresMean = mean;
711  data[key].t0err = rms;
712  data[key].t0fittype = 1;
713  return mean;
714  }
715 
716  float fitrms=rms;
717  if (fitrms>5.0) fitrms=5.0;
718 
719  tresHist->Fit("gaus","QR","",mean-rms,mean+rms);
720  if (tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2) < m_mintres ||
721  tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2) > m_maxtres) {
722  data[key].tresMean = mean;
723  data[key].t0err = rms;
724  data[key].t0fittype = 1;
725  return mean;
726  }
727 
728  tresHist->Fit("gaus","QR","",tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2),tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2));
729  if (tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2) < m_mintres ||
730  tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2) > m_maxtres) {
731  data[key].tres = tresHist->GetFunction("gaus")->GetParameter(2);
732  data[key].tresMean = tresHist->GetFunction("gaus")->GetParameter(1);
733  data[key].t0err = tresHist->GetFunction("gaus")->GetParError(1);
734  data[key].t0fittype = 2;
735  return tresHist->GetFunction("gaus")->GetParameter(1);
736  }
737 
738  tresHist->Fit("gaus","QR","",tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2),tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2));
739  if (tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2) < m_mintres ||
740  tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2) > m_maxtres) {
741  data[key].tres = tresHist->GetFunction("gaus")->GetParameter(2);
742  data[key].tresMean = tresHist->GetFunction("gaus")->GetParameter(1);
743  data[key].t0err = tresHist->GetFunction("gaus")->GetParError(1);
744  data[key].t0fittype = 2;
745  return tresHist->GetFunction("gaus")->GetParameter(1);
746  }
747 
748  tresHist->Fit("gaus","QR","",tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2),tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2));
749  if (tresHist->GetFunction("gaus")->GetParameter(1) - tresHist->GetFunction("gaus")->GetParameter(2) < m_mintres ||
750  tresHist->GetFunction("gaus")->GetParameter(1) + tresHist->GetFunction("gaus")->GetParameter(2) > m_maxtres) {
751  data[key].tresMean = mean;
752  data[key].t0err = rms;
753  data[key].t0fittype = 1;
754  return mean;
755  }
756 
757  data[key].tres = tresHist->GetFunction("gaus")->GetParameter(2);
758  data[key].tresMean = tresHist->GetFunction("gaus")->GetParameter(1);
759  data[key].t0err = tresHist->GetFunction("gaus")->GetParError(1);
760  data[key].t0fittype = 2;
761 
762  // protection against wrong fits:
763  if ( std::abs(tresHist->GetFunction("gaus")->GetParameter(2))>10 || std::abs(tresHist->GetFunction("gaus")->GetParameter(1))>10 ) {
764 
765  data[key].tres = tresHist->GetRMS();
766  data[key].tresMean = tresHist->GetMean();
767  data[key].t0err = tresHist->GetRMS();
768  data[key].t0fittype = 6;
769  return tresHist->GetMean();
770  }
771 
772 
773  gStyle->SetOptFit(1111);
774 
775  return tresHist->GetFunction("gaus")->GetParameter(1);
776  }
777 
778 
779 float Calibrator::FitResidual ATLAS_NOT_THREAD_SAFE (const std::string& key, TH1F* resHist){ // Global gStyle is used.
780 
781  float mean = resHist->GetMean();
782 
783  resHist->Fit("gaus","QRI","",mean-0.3,mean+0.3);
784  resHist->Fit("gaus","QRI","",resHist->GetFunction("gaus")->GetParameter(1) - 1.5*resHist->GetFunction("gaus")->GetParameter(2),resHist->GetFunction("gaus")->GetParameter(1) + 1.5*resHist->GetFunction("gaus")->GetParameter(2));
785  resHist->Fit("gaus","QRI","",resHist->GetFunction("gaus")->GetParameter(1) - 1.5*resHist->GetFunction("gaus")->GetParameter(2),resHist->GetFunction("gaus")->GetParameter(1) + 1.5*resHist->GetFunction("gaus")->GetParameter(2));
786  resHist->Fit("gaus","QRI","",resHist->GetFunction("gaus")->GetParameter(1) - 1.5*resHist->GetFunction("gaus")->GetParameter(2),resHist->GetFunction("gaus")->GetParameter(1) + 1.5*resHist->GetFunction("gaus")->GetParameter(2));
787 
788  gStyle->SetOptFit(1111);
789 
790  data[key].res=resHist->GetFunction("gaus")->GetParameter(2);
791  data[key].resMean=resHist->GetFunction("gaus")->GetParameter(1);
792  data[key].reserr=resHist->GetFunction("gaus")->GetParError(2);
793 
794  return resHist->GetFunction("gaus")->GetParameter(2);
795 
796 }
797 
798 TDirectory* Calibrator::Calibrate ATLAS_NOT_THREAD_SAFE (TDirectory* dir, std::string key, const std::string& opt, caldata * caldata_above){ // Thread unsafe FitResidual, FitRt, FitTimeResidual are used.
799 
800  //set some bool flags
801  bool calrt=opt.find('R')!=std::string::npos;
802  bool calt0=opt.find('T')!=std::string::npos;
803  bool donothing=opt.find('N')!=std::string::npos;
804  bool prnt=opt.find('P')!=std::string::npos;
805  bool useref=opt.find('B')!=std::string::npos;
806  bool useshortstw=opt.find('S')!=std::string::npos;
807 
808  if (donothing) return dir;
809 
810  data[key].calflag=true;
811 
812  bool enough_t0 = data[key].ntres>=m_mint0stat;
813  bool enough_rt = data[key].nrt>=m_minrtstat;
814 
815  if ((enough_rt && calrt) || (enough_t0 && calt0)) {
816  m_hdirs[key] = dir->mkdir(Form("%s%s",m_name.data(),key.data()));
817  m_hdirs[key]->cd();
818  if(level<2) {
819  std::cout << " Calibrator::Calibrate Parent Directory: " << std::string(dir->GetPath()) << std::endl;
820  std::cout << " Calibrator::Calibrate key " << key << " name " << m_name << " entries for rt " << data[key].nrt << " min req " << m_minrtstat << std::endl;
821  std::cout << " Calibrator::Calibrate key " << key << " entries for t0 " << data[key].ntres << " min req " << m_mint0stat << std::endl;
822  std::cout << " Calibrator::Calibrate Directory: " << std::string(m_hdirs[key]->GetPath()) << std::endl;
823  }
824  }
825  else m_hdirs[key]=dir;
826 
827  //Fit also the residual if an rt or t0 calibration is made
828  if ((enough_rt && calrt) || (enough_t0 && calt0)) {
829 
830  m_resHists[key] = new TH1F("residual","residual",m_nbinsres,m_minres,m_maxres);
831  int ncont=0;
832  for (int i=0;i<100;i++) {
833  m_resHists[key]->SetBinContent(i+1,data[key].reshist[i]);
834  ncont+=data[key].reshist[i];
835  }
836  m_resHists[key]->SetEntries((int)data[key].nres);
837  if (bequiet) m_resHists[key]->SetDirectory(nullptr);
838  if(ncont==0) {
839  std::cout << " No bin content " << std::endl;
840  } else {
841  FitResidual(key,m_resHists[key]);
842  }
843 
844  }
845 
846 
847  if (prnt) printf("Calibrator::calibrate %8s %-14s: \n",m_name.data(),key.data());
848 
849  //Calibrate r-m_t
850  if (nort){
851  //use old data
852  //std::cout << " use old data " << std::endl;
853  data[key].rtflag=true;
854  data[key].rtt0=caldata_above->rtt0;
855  data[key].rtgraph=caldata_above->rtgraph;
856  for (int i=0;i<4;i++) data[key].rtpar[i]=data[key].oldrtpar[i];
857  if (prnt) printf("RT << %7i (%8.1e) %8.1e %8.1e %8.1e, %3.2f : ", data[key].nrt,data[key].rtpar[0],data[key].rtpar[1],data[key].rtpar[2],data[key].rtpar[3], data[key].rtt0);
858  }
859  else{
860  //do fit
861  if ((calrt && enough_rt) || level==0) {
862  data[key].rtflag=true;
863  m_rtHists[key] = new TH2F("rt-relation","rt-relation",m_nbinst,m_mint,m_maxt,m_nbinsr,m_minr,m_maxr); //make a root histogram
864  for (int i=0;i<m_nbinst;i++) {
865  for (int j=0;j<m_nbinsr;j++) {
866  m_rtHists[key]->SetBinContent(i+1,j+1,data[key].rthist[i+m_nbinst*j]);
867  }
868  }
869  m_rtHists[key]->SetEntries(data[key].nrt);
870  if(level<2) std::cout << " Calibrator::calibrate do rt fit for key " << key << " with number of entries " << m_rtHists[key]->GetEntries() << std::endl;
871  if (bequiet) {
872  m_rtHists[key]->SetDirectory(nullptr);
873  if(level<2) std::cout << " remove this level from calout.root file " << std::endl;
874  }
875  data[key].rtt0=FitRt(key,opt,m_rtHists[key],m_hdirs[key]); //do the fit
876  if (prnt) printf("RT %7i (%8.1e) %8.1e %8.1e %8.1e, %3.2f : ", data[key].nrt, data[key].rtpar[0], data[key].rtpar[1], data[key].rtpar[2], data[key]. rtpar[3], data[key].rtt0);
877  }
878  else{
879  //use data from level above
880  data[key].rtgraph=caldata_above->rtgraph;
881  data[key].rtt0=caldata_above->rtt0;
882  for (int i=0;i<4;i++) data[key].rtpar[i]=caldata_above->rtpar[i];
883  if (prnt) printf("RT /\\ %7i (%8.1e) %8.1e %8.1e %8.1e, %3.2f : ", data[key].nrt,data[key].rtpar[0],data[key].rtpar[1],data[key].rtpar[2],data[key].rtpar[3], data[key].rtt0);
884  }
885  }
886 
887 
888  //Calibrate t0
889  if (not0){
890  //use old data
891  //std::cout << " use old data for T0 " std::endl;
892  data[key].t0flag=true;
893  data[key].t0=data[key].oldt02 + data[key].rtt0;
894  data[key].t0err=0;
895  data[key].t0off=0;
896  data[key].t0fittype = 5;
897  if (prnt) printf("T0 << %7i %05.2f%+05.2f%+05.2f=%05.2f \n", data[key].ntres, data[key].oldt02, 0.0, data[key].rtt0, data[key].t0);
898  }
899  else{
900  if (useref && level==5){
901  //use chip reference values
902  //std::cout << " use chip reference for T0 " << std::endl;
903  data[key].t0flag=true;
904  data[key].t0=caldata_above->t0 + data[key].reft0 + data[key].rtt0;
905  data[key].t0err=caldata_above->t0err;
906  data[key].t0off=data[key].t0-caldata_above->t0;
907  data[key].t0fittype = 3;
908  if (prnt) printf("T0 ** %7i %05.2f%+05.2f%+05.2f=%05.2f \n", data[key].ntres, caldata_above->t0, data[key].reft0, data[key].rtt0, data[key].t0);
909  }
910  else {
911  //do fit
912  if ((calt0 && enough_t0) || level==0) {
913  data[key].t0flag=true;
914  m_tresHists[key] = new TH1F("timeresidual","time residual",m_nbinstres,m_mintres,m_maxtres); //make a root histogram
915 
916  for (int i=0;i<100;i++) {
917  m_tresHists[key]->SetBinContent(i+1,data[key].m_treshist[i]);
918  }
919  m_tresHists[key]->SetEntries(data[key].ntres);
920  if(level<2) std::cout << " Calibrator::calibrate do t0 fit for key " << key << " Numb. entries " << m_rtHists[key]->GetEntries() << std::endl;
921  if (bequiet) {
922  if(level<2) std::cout << " Calibrator::calibrate: remove this level from calibout.root " << std::endl;
923  m_tresHists[key]->SetDirectory(nullptr);
924  }
925 
926  data[key].t0=data[key].oldt02 + FitTimeResidual(key,m_tresHists[key]) + data[key].rtt0 + m_t0shift; //do the fit and modify t0
927  data[key].t0off=data[key].t0-caldata_above->t0; //calculate t0 offset from level above
928  if (data[key].t0<0) data[key].t0=0;
929 
930  if (prnt) printf("T0 %7i %05.2f%+05.2f%+05.2f%+05.2f=%05.2f \n", data[key].ntres, data[key].oldt02, data[key].t0-data[key].oldt02-data[key].rtt0, data[key].rtt0, m_t0shift, data[key].t0);
931 
932 
933  }
934  //use data from level above
935  else {
936  //std::cout << " use data from the level above " << std::endl;
937  //TEMP FIX to not destroy right T0s
938  if (data[key].oldt02 + (caldata_above->t0 - caldata_above->oldt02) >0) data[key].t0=data[key].oldt02 + (caldata_above->t0 - caldata_above->oldt02);
939  else data[key].t0= 0;
940  //TEMP FIX to not destroy right T0s
941 
942  //add the short straw correction here. In this way, the shift is only done when contants at STRAW level come from level above.
943  if ((level == 6 && useshortstw) && std::abs(data[key].det)<2 && (data[key].lay==0 && data[key].stl<9) ) data[key].t0=caldata_above->t0-0.75;
944  data[key].t0err=caldata_above->t0err;
945  data[key].t0off=data[key].t0-caldata_above->t0 + data[key].rtt0;
946  data[key].t0fittype = 4;
947  if (prnt) printf("T0 /\\ %7i %05.2f%+05.2f%+05.2f=%05.2f \n", data[key].ntres, caldata_above->t0, caldata_above->oldt02, data[key].oldt02 , data[key].t0);
948  }
949  }
950  }
951 
952 
953  if (prnt && !bequiet) std::cout << " H";
954 
955  if (prnt) std::cout << std::endl;
956 
957  return m_hdirs[key];
958 
959 }
960 
961 
962 int Calibrator::AddHit(const std::string& key, const databundle & d, int* binhist, bool makehist){
963 
964  int tresbin=Simple1dHist(m_mintres,m_maxtres,m_nbinstres,d.tres);
965  int resbin=Simple1dHist(m_minres,m_maxres,m_nbinsres,d.res);
966  int rtbin=Simple2dHist(m_mint,m_maxt,m_nbinst,m_minr,m_maxr,m_nbinsr,d.t,d.r);
967  int npop;
968  int ibin;
969  int value;
970 
971  //if the it is the first hit or histogram
972  if (data.find(key) == data.end()){
973 
974  //create a new histogram object
975  caldata* hist=new caldata(makehist,m_nbinst,m_nbinsr);
976 
977  //out of memory?
978  if (hist == nullptr){
979  std::cout << "OUT OF MEMORY!" << std::endl;
980  return -1;
981  }
982 
983  hist->ntres=0;
984  hist->nres=0;
985  hist->nrt=0;
986  hist->sumt0=0;
987 
988  //zero histogram bins
989  for (int i=0;i<m_nbinsr*m_nbinst+200;i++) {
990  if (makehist) hist->rthist[i]=0;
991  if (i<100) {
992  if (makehist) hist->m_treshist[i]=0;
993  if (makehist) hist->reshist[i]=0;
994  }
995  }
996 
997  if (binhist==nullptr){ //if it is a hit
998 
999  if (tresbin>=0){
1000  if (makehist) hist->m_treshist[tresbin]=d.weight; //add value to bin
1001  hist->ntres=1; //set bin content to 1
1002  m_nhits=1;
1003  }
1004  if (resbin>=0){
1005  if (makehist) hist->reshist[resbin]=1;
1006  hist->nres=1;
1007  }
1008  if (rtbin>=0){
1009  if (makehist) hist->rthist[rtbin]=1;
1010  hist->nrt=1;
1011  }
1012 
1013  hist->sumt0=d.t0;
1014  }
1015  else { //if it is a histogram
1016 
1017  npop=binhist[0]; //the number of populated (non-zero) bins
1018  if(level<2) std::cout << " Calibrator::AddHit called for key " << key << " add histogram with "<< npop << " non-zero bins " << std::endl;
1019 
1020  for(int ipop=2;ipop<2*npop+2;ipop=ipop+2) { //loop over the data
1021 
1022  ibin=binhist[ipop]; //bin number
1023  value=binhist[ipop+1]; //value
1024 
1025  if (ibin<100) { //timeresidual histogram
1026  m_ntreshits+=value;
1027  hist->ntres+=value;
1028  m_nhits+=value;
1029  if (makehist) hist->m_treshist[ibin]=(float)value;
1030  hist->sumt0+=d.t0*value;
1031  }
1032  else if (ibin>=100 && ibin<200) { //residual histogram
1033  m_nreshits+=value;
1034  hist->nres+=value;
1035  if (makehist) hist->reshist[ibin-100]=(float)value;
1036  }
1037  else { //rt histogram
1038  m_nrthits+=value;
1039  hist->nrt+=value;
1040  if (makehist) hist->rthist[ibin-200]=(float)value;
1041  }
1042 
1043  }
1044  }
1045 
1046  //initialize parameters
1047  if (level>0) hist->det=d.det; else hist->det=-3; //
1048  if (level>1) hist->lay=d.lay; else hist->lay=-3;
1049  if (level>2) hist->mod=d.mod; else hist->mod=-3;
1050  if (level>3) hist->brd=d.brd; else hist->brd=-3;
1051  if (level>4) hist->chp=d.chp; else hist->chp=-3;
1052  if (level>5) hist->stw=d.stw; else hist->stw=-3;
1053  if (level>2) hist->stl=d.stl; else hist->stl=-3;
1054  if (level>5) hist->sid=d.sid; else hist->sid=-3;
1055  hist->res=0;
1056  hist->resMean=0;
1057  hist->reserr=0;
1058  hist->tres=0;
1059  hist->tresMean=0;
1060  hist->t0=0;
1061  if (level==5) hist->reft0=d.rt0; else hist->reft0=0;
1062  hist->t0off=0;
1063  hist->t0err=0;
1064  hist->nhits=1;
1065  hist->t0fittype=0;
1066 
1067  hist->x=d.x; //straw x coordinate
1068  hist->y=d.y; //straw y coordinate
1069  hist->z=d.z; //straw z coordinate
1070 
1071  hist->oldt02=d.t0; //straw old t0 value
1072  hist->oldrtpar[0]=d.rtpar[0];
1073  hist->oldrtpar[1]=d.rtpar[1];
1074  hist->oldrtpar[2]=d.rtpar[2];
1075  hist->oldrtpar[3]=d.rtpar[3];
1076 
1077  hist->sumx=0;
1078  hist->sumy=0;
1079  hist->sumz=0;
1080 
1081  hist->rtflag=false;
1082  hist->t0flag=false;
1083  hist->calflag=false;
1084 
1085 
1086 
1087  for (unsigned int i =0; i < 4; i++){
1088  hist->rtpar[i]=-10.0;
1089  }
1090 
1091  data[key]=*hist; //save the histogram in the map
1092 
1093 
1094 
1095  data[key].oldt02 = AccumulativeMean(data[key].nhits, data[key].oldt02, d.t0); //update old t0 m_mean value
1096 
1097  data[key].nhits++; // increment m_nhits
1098 
1099  if(level<2) std::cout << " First hit for det " << data[key].det << " lay " << data[key].lay << " t0 " << data[key].t0 << " upd t0 " << data[key].oldt02 << " nhits: " << data[key].nhits << std::endl;
1100 
1101  delete hist;
1102 
1103  return 1;
1104 }
1105  else { //if not the first hit or histogram
1106 
1107  //increment histogram bins
1108  if (binhist==nullptr){
1109 
1110  if (tresbin>=0){
1111  if (makehist) data[key].m_treshist[tresbin]=data[key].m_treshist[tresbin]+d.weight;
1112  data[key].ntres++;
1113  m_nhits=1;
1114  }
1115  if (resbin>=0){
1116  if (makehist) data[key].reshist[resbin]++;
1117  data[key].nres++;
1118  }
1119  if (rtbin>=0){
1120  if (makehist) data[key].rthist[rtbin]++;
1121  data[key].nrt++;
1122  }
1123 
1124  data[key].sumt0+=d.t0;
1125 
1126  }
1127 
1128  else {
1129 
1130  npop=binhist[0];
1131 
1132  for(int ipop=2;ipop<2*npop+2;ipop=ipop+2) {
1133 
1134  ibin=binhist[ipop];
1135  value=binhist[ipop+1];
1136 
1137  if (ibin<100) {
1138  m_ntreshits+=value;
1139  m_nhits+=value;
1140  data[key].ntres+=value;
1141  if (makehist) data[key].m_treshist[ibin]+=(float)value;
1142  data[key].sumt0+=d.t0*value;
1143  }
1144  else if (ibin>=100 && ibin<200) {
1145  m_nreshits+=value;
1146  data[key].nres+=value;
1147  if (makehist) data[key].reshist[ibin-100]+=(float)value;
1148  }
1149  else {
1150  m_nrthits+=value;
1151  data[key].nrt+=value;
1152  if (makehist) data[key].rthist[ibin-200]+=(float)value;
1153  }
1154 
1155  }
1156 
1157  }
1158 
1159 
1160  data[key].sumx+=d.x;
1161  data[key].sumy+=d.y;
1162  data[key].sumz+=d.z;
1163  data[key].oldt02 = AccumulativeMean(data[key].nhits, data[key].oldt02, d.t0);
1164  data[key].oldrtpar[0] = AccumulativeMean(data[key].nhits, data[key].oldrtpar[0], d.rtpar[0]);
1165  data[key].oldrtpar[1] = AccumulativeMean(data[key].nhits, data[key].oldrtpar[1], d.rtpar[1]);
1166  data[key].oldrtpar[2] = AccumulativeMean(data[key].nhits, data[key].oldrtpar[2], d.rtpar[2]);
1167  data[key].oldrtpar[3] = AccumulativeMean(data[key].nhits, data[key].oldrtpar[3], d.rtpar[3]);
1168  data[key].x = AccumulativeMean(data[key].nhits, data[key].x, d.x);
1169  data[key].y = AccumulativeMean(data[key].nhits, data[key].y, d.y);
1170  data[key].z = AccumulativeMean(data[key].nhits, data[key].z, d.z);
1171 
1172  data[key].nhits++; //increment hit counts
1173 
1174  return 0;
1175  }
1176 
1177 
1178 }
1179 
1180 void Calibrator::WriteStat(TDirectory* dir){
1181  dir->cd();
1182  TNtuple* stattup = new TNtuple(Form("%stuple",m_name.data()),"statistics","det:lay:mod:brd:chp:sid:stl:stw:rtflag:t0flag:t0:oldt0:rt0:dt0:t0offset:ftype:nhits:nt0:nrt:res:resMean:dres:tres:tresMean:x:y:z");
1183  for(std::map<std::string,caldata>::iterator ihist=data.begin(); ihist!=data.end(); ++ihist){
1184  if ((ihist->second).calflag) {
1185  //std::cout << " Writing TNtuple for name: " << m_name << " det: " << (ihist->second).det << " lay " << (ihist->second).lay << " mod: " << (ihist->second).mod << " brd " << (ihist->second).brd << " res width " << (ihist->second).res << " tres mean " << (ihist->second).tresMean << std::endl;
1186  float const ntvar[27]={
1187  float((ihist->second).det),
1188  float((ihist->second).lay),
1189  float((ihist->second).mod),
1190  float((ihist->second).brd),
1191  float((ihist->second).chp),
1192  float((ihist->second).sid),
1193  float((ihist->second).stl),
1194  float((ihist->second).stw),
1195  float((int)(ihist->second).rtflag),
1196  float((int)(ihist->second).t0flag),
1197  (ihist->second).t0,
1198  float((ihist->second).oldt02),
1199  //oldt0(ihist->first),
1200  (ihist->second).reft0,
1201  (ihist->second).t0err,
1202  (ihist->second).t0off,
1203  float((ihist->second).t0fittype),
1204  (ihist->second).nhits,
1205  (float)(ihist->second).ntres,
1206  (float)(ihist->second).nrt,
1207  (float)(ihist->second).res,
1208  (float)(ihist->second).resMean,
1209  (ihist->second).reserr,
1210  (ihist->second).tres,
1211  (ihist->second).tresMean,
1212  (ihist->second).x,
1213  (ihist->second).y,
1214  (ihist->second).z
1215  };
1216  stattup->Fill(ntvar);
1217  }
1218 
1219  }
1220 
1221  stattup->Write();
1222 
1223 }
1224 
1225 void Calibrator::DumpConstants(){
1226 
1227  std::ofstream dumpfile( "calib_constants_out.txt", std::ios::out | std::ios::app);
1228 
1229  for (const std::pair<const std::string, caldata>& p : data) {
1230  dumpfile << p.first << " " << p.second.t0 << " " << std::endl;
1231 
1232  }
1233 
1234  dumpfile.close();
1235 
1236 }
grepfile.info
info
Definition: grepfile.py:38
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Calibrator::m_name
std::string m_name
The name of the Calibrator instance.
Definition: Calibrator.h:359
beamspotman.r
def r
Definition: beamspotman.py:676
Calibrator::printt0
bool printt0
if true a t0 entry in the calibration output file is prined for each sub-module in this sub-level
Definition: Calibrator.h:349
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
Calibrator::dort
bool dort
if true an rt fit is made, if false the value from the level above is used
Definition: Calibrator.h:341
Calibrator::m_maxres
float m_maxres
upper limit of 1D residual histogram
Definition: Calibrator.h:375
Calibrator::m_nhits
int m_nhits
...
Definition: Calibrator.h:384
pol3deg
double pol3deg(double *x, double *par)
Definition: Calibrator.cxx:303
checkFileSG.line
line
Definition: checkFileSG.py:75
TH1F::GetBinContent
double GetBinContent(int) const
Definition: rootspy.cxx:326
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
Calibrator::m_nbinsr
int m_nbinsr
number of r-bins in the 2D rt histogram
Definition: Calibrator.h:364
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
max
#define max(a, b)
Definition: cfImp.cxx:41
caldata::rthist
std::vector< float > rthist
the 2D rt histogram (20x32 bins)
Definition: Calibrator.h:147
caldata::t0fittype
int t0fittype
the type of time residual fit that was made
Definition: Calibrator.h:120
Calibrator::printlog
bool printlog
if true a log entry is prined for each sub-modile in this sub-level
Definition: Calibrator.h:348
caldata::nhits
float nhits
the number of straws in the sub-module
Definition: Calibrator.h:133
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
caldata::ntres
int ntres
number of time residual histogram entries
Definition: Calibrator.h:117
caldata::brd
int brd
board
Definition: Calibrator.h:112
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
caldata::sid
int sid
straw ID
Definition: Calibrator.h:116
TH2F
Definition: rootspy.cxx:420
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
Calibrator::m_maxtres
float m_maxtres
upper limit of 1D time residual histogram
Definition: Calibrator.h:373
Calibrator::m_mintres
float m_mintres
lower limit of 1D time residual histogram
Definition: Calibrator.h:372
hist_file_dump.d
d
Definition: hist_file_dump.py:137
RtGraph::~RtGraph
~RtGraph()
the destructor
Definition: Calibrator.cxx:287
Calibrator::UpdateOldConstants
int UpdateOldConstants()
...
Definition: Calibrator.cxx:529
plotmaker.hist
hist
Definition: plotmaker.py:148
Calibrator::HasKey
bool HasKey(const std::string &) const
...
Definition: Calibrator.cxx:452
pol3deg2
double pol3deg2(double *x, double *par)
Definition: Calibrator.cxx:309
RtGraph::HIGH
@ HIGH
Definition: Calibrator.h:51
ALFA_EventTPCnv_Dict::t0
std::vector< ALFA_RawData_p1 > t0
Definition: ALFA_EventTPCnvDict.h:42
Calibrator::m_maxt
float m_maxt
upper limit of t in 2D rt histogram
Definition: Calibrator.h:371
TH1D
Definition: rootspy.cxx:342
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
caldata
A structure to contain data associated with the calibration of a certain sub-module.
Definition: Calibrator.h:104
Calibrator::CheckSelection
bool CheckSelection(int)
Checks if a given sub-module is selected for calibration.
Definition: Calibrator.cxx:457
Calibrator::printrt
bool printrt
if true an rt entry in the calibration output file is prined for each sub-module in this sub-level
Definition: Calibrator.h:350
caldata::x
float x
sub-module x position (average of all straws in the module)
Definition: Calibrator.h:134
Calibrator::m_ntreshits
int m_ntreshits
...
Definition: Calibrator.h:381
athena.value
value
Definition: athena.py:122
RtGraph::GOOD
@ GOOD
Definition: Calibrator.h:51
caldata::det
int det
detector (barrel or end-cap)
Definition: Calibrator.h:109
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Calibrator::dot0
bool dot0
if true a time residual fit is made, if false the value from the level above is used
Definition: Calibrator.h:342
RtGraph::m_tv
float * m_tv
Definition: Calibrator.h:66
caldata::resMean
float resMean
the residual mean
Definition: Calibrator.h:122
caldata::oldt02
double oldt02
the old sub-module t0 (average of t0 for all straws in the module)
Definition: Calibrator.h:137
caldata::sumt0
double sumt0
...
Definition: Calibrator.h:138
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
caldata::lay
int lay
layer
Definition: Calibrator.h:110
x
#define x
Calibrator::m_nrthits
int m_nrthits
...
Definition: Calibrator.h:383
Calibrator::m_nbinst
int m_nbinst
number of t-bins in the 2D rt histogram
Definition: Calibrator.h:365
Calibrator::usebref
bool usebref
if true chip reference t0 values (offset from board mean) are used
Definition: Calibrator.h:346
Calibrator::floatp3
bool floatp3
if true the 3rd order coeficcient of the rt fit function is not fixed to 0
Definition: Calibrator.h:352
RtGraph::tval
std::vector< double > tval
the r values
Definition: Calibrator.h:43
python.iconfTool.models.loaders.level
level
Definition: loaders.py:20
caldata::y
float y
sub-module y position (average of all straws in the module)
Definition: Calibrator.h:135
caldata::rtflag
bool rtflag
flag indicating if an R-t calibration has been made
Definition: Calibrator.h:143
Calibrator::Skip
bool Skip()
...
Definition: Calibrator.cxx:462
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
Calibrator::PrintInfo
std::string PrintInfo()
Definition: Calibrator.cxx:467
TH1D::GetBinContent
double GetBinContent(int) const
Definition: rootspy.cxx:347
Calibrator::dores
bool dores
if true a residual fit is made
Definition: Calibrator.h:343
caldata::sumz
double sumz
...
Definition: Calibrator.h:141
Calibrator::usep0
bool usep0
if true the 0th order coeficcient of the rt fit function is not set to 0 in the calibration output fi...
Definition: Calibrator.h:351
caldata::t0off
float t0off
the t0 offset from the level below
Definition: Calibrator.h:129
RtGraph::EMPTY
@ EMPTY
Definition: Calibrator.h:51
caldata::stl
int stl
straw-layer
Definition: Calibrator.h:114
caldata::caldata
caldata()
Definition: Calibrator.cxx:33
caldata::t0flag
bool t0flag
flag indicating if a t0 calibration has been made
Definition: Calibrator.h:144
Calibrator::m_minr
float m_minr
lower limit of r in 2D rt histogram
Definition: Calibrator.h:368
RtGraph::m_rightval
float * m_rightval
Definition: Calibrator.h:57
Calibrator::GetOptString
std::string GetOptString() const
Creates a string summarizing what is being done at this sub-level.
Definition: Calibrator.cxx:485
lumiFormat.i
int i
Definition: lumiFormat.py:92
z
#define z
D3PDSizeSummary.ff
ff
Definition: D3PDSizeSummary.py:305
caldata::mod
int mod
phi module
Definition: Calibrator.h:111
beamspotman.n
n
Definition: beamspotman.py:731
caldata::t0
float t0
the new t0
Definition: Calibrator.h:126
caldata::rtgraph
RtGraph * rtgraph
the rt graph
Definition: Calibrator.h:148
caldata::rtpar
std::array< float, 4 > rtpar
the new rt-parameters
Definition: Calibrator.h:132
caldata::nrt
int nrt
number of rt histogram entries
Definition: Calibrator.h:118
WritePulseShapeToCool.det
det
Definition: WritePulseShapeToCool.py:204
Calibrator::DumpConstants
void DumpConstants()
Definition: Calibrator.cxx:1225
RtGraph::m_btype
bintype * m_btype
Definition: Calibrator.h:52
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:14
RtGraph::m_mindistance
float m_mindistance
Definition: Calibrator.h:59
Calibrator::Calibrator
Calibrator()
Definition: Calibrator.cxx:339
RtGraph::m_leftval
float * m_leftval
Definition: Calibrator.h:56
RtGraph::m_chname
std::string m_chname
the histogram name
Definition: Calibrator.h:49
Calibrator::m_mint
float m_mint
lower limit of t in 2D rt histogram
Definition: Calibrator.h:370
RtGraph::m_d
float m_d
Definition: Calibrator.h:63
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
ATLAS_NOT_THREAD_SAFE
float Calibrator::FitRt ATLAS_NOT_THREAD_SAFE(const std::string &key, const std::string &opt, TH2F *rtHist, TDirectory *dir)
Definition: Calibrator.cxx:555
RtGraph::m_rightsig
float * m_rightsig
Definition: Calibrator.h:53
RtGraph::m_leftsig
float * m_leftsig
Definition: Calibrator.h:54
caldata::sumx
double sumx
...
Definition: Calibrator.h:139
trrel_dines
double trrel_dines(double *x, double *par)
Definition: Calibrator.cxx:316
Calibrator::Simple2dHist
int Simple2dHist(float, float, int, float, float, int, float, float)
A 2D histograming function.
Definition: Calibrator.cxx:439
Calibrator::selection
std::set< int > selection
a set containing the sub-modules to be calibrated
Definition: Calibrator.h:354
RtGraph::m_et
float m_et
Definition: Calibrator.h:64
beamspotman.dir
string dir
Definition: beamspotman.py:623
caldata::t0err
float t0err
the new to error
Definition: Calibrator.h:127
min
#define min(a, b)
Definition: cfImp.cxx:40
caldata::tres
float tres
the time residual
Definition: Calibrator.h:124
caldata::tresMean
float tresMean
the time residual mean
Definition: Calibrator.h:125
RtGraph::m_etv
float * m_etv
Definition: Calibrator.h:68
Calibrator::WriteStat
void WriteStat(TDirectory *)
Creates an ntuple with entries containing data associated with the sub-modules in a sub level.
Definition: Calibrator.cxx:1180
Calibrator::Simple1dHist
int Simple1dHist(float, float, int, float)
A 1D histograming function.
Definition: Calibrator.cxx:432
fitman.fitresult
fitresult
Definition: fitman.py:590
caldata::chp
int chp
chip
Definition: Calibrator.h:113
caldata::sumy
double sumy
...
Definition: Calibrator.h:140
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
pmontree.opt
opt
Definition: pmontree.py:16
Calibrator::useshortstraws
bool useshortstraws
if true a shift of -0.75 ns is applied for straws in layer ==0, strawlayer < 9, and when doing calibr...
Definition: Calibrator.h:353
Calibrator::m_maxr
float m_maxr
upper limit of r in 2D rt histogram
Definition: Calibrator.h:369
caldata::reserr
float reserr
the residual error
Definition: Calibrator.h:123
Calibrator::PrintStat
std::string PrintStat()
Prints some sub-level statistics.
Definition: Calibrator.cxx:478
databundle
A structure to contain hit data.
Definition: Calibrator.h:77
RtGraph::RtGraph
RtGraph(TH2F *, int, const char *, bool, TDirectory *)
the constructor
Definition: Calibrator.cxx:118
RtGraph::m_chtit
std::string m_chtit
the histogram title
Definition: Calibrator.h:50
RtGraph::m_mean
float m_mean
Definition: Calibrator.h:58
perfmonmt-refit.units
string units
Definition: perfmonmt-refit.py:77
RtGraph::mintime
float mintime
the minimum t-value
Definition: Calibrator.h:47
caldata::calflag
bool calflag
flag indicating if any calibration has been made
Definition: Calibrator.h:142
Calibrator::AccumulativeMean
float AccumulativeMean(float, float, float)
Increments a single bin in all three histograms in a sub-module.
Definition: Calibrator.cxx:447
RtGraph::LOW
@ LOW
Definition: Calibrator.h:51
RtGraph::trgr
TGraphErrors * trgr
the t(r) graph
Definition: Calibrator.h:44
python.PyAthena.v
v
Definition: PyAthena.py:157
Calibrator::AddHit
int AddHit(const std::string &, const databundle &, int *, bool)
Adds hits to a sub-module either in the form of a single straw hit or a histogram containing all the ...
Definition: Calibrator.cxx:962
RtGraph::m_maxval
float * m_maxval
Definition: Calibrator.h:55
Calibrator.h
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
Calibrator::GetSelString
std::string GetSelString()
Creates a string summarizing which modules are being calibrated at this sub-level.
Definition: Calibrator.cxx:506
y
#define y
Calibrator::m_minres
float m_minres
lower limit of 1D residual histogram
Definition: Calibrator.h:374
TH1F
Definition: rootspy.cxx:320
RtGraph::m_dv
float * m_dv
Definition: Calibrator.h:67
caldata::m_treshist
std::vector< float > m_treshist
the 1D time residual histogram (100 bins)
Definition: Calibrator.h:145
RtGraph::rtgr
TGraphErrors * rtgr
the r(t) graph
Definition: Calibrator.h:45
caldata::nres
int nres
number of residual histogram entries
Definition: Calibrator.h:119
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
generate::GetEntries
double GetEntries(TH1D *h, int ilow, int ihi)
Definition: rmsFrac.cxx:20
RtGraph::m_edv
float * m_edv
Definition: Calibrator.h:69
caldata::stw
int stw
straw number (within the strawlayer)
Definition: Calibrator.h:115
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:569
Calibrator::~Calibrator
~Calibrator()
Destructor.
Definition: Calibrator.cxx:429
caldata::reshist
std::vector< float > reshist
the 1D residual histogram (100 bins)
Definition: Calibrator.h:146
beamspotnt.rms
rms
Definition: bin/beamspotnt.py:1266
Calibrator::m_nbinsres
int m_nbinsres
number of bins in the 1D residual histogram
Definition: Calibrator.h:367
caldata::rtt0
float rtt0
the t0 fron the R-t fit
Definition: Calibrator.h:130
RtGraph::m_ipoint
int m_ipoint
Definition: Calibrator.h:61
RtGraph::rval
std::vector< double > rval
array of the histograms for all bins
Definition: Calibrator.h:42
caldata::~caldata
~caldata()
Definition: Calibrator.cxx:71
caldata::reft0
float reft0
the reference t0 (offset from board mean)
Definition: Calibrator.h:128
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
RtGraph
A class for generating a r-t and t-r graphs by binning the 2D histograms in Calibrator::rtHists in r ...
Definition: Calibrator.h:34
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
RtGraph::m_ed
float m_ed
Definition: Calibrator.h:65
Calibrator::bequiet
bool bequiet
if true no histograms are written to the root file
Definition: Calibrator.h:347
Calibrator::m_nbinstres
int m_nbinstres
number of bins in the 1D time residual histogram
Definition: Calibrator.h:366
checker_macros.h
Define macros for attributes used to control the static checker.
dumpfile
std::ofstream dumpfile("dumpfile.log")
Calibrator::level
int level
the sub-level of the Calibrator instance
Definition: Calibrator.h:355
RtGraph::npoints
int npoints
the number of graph points
Definition: Calibrator.h:46
RtGraph::m_maxbin
int * m_maxbin
Definition: Calibrator.h:60
rtrel_dines
double rtrel_dines(double *x, double *par)
Definition: Calibrator.cxx:326
Calibrator::m_nreshits
int m_nreshits
...
Definition: Calibrator.h:382
caldata::z
float z
sub-module z position (average of all straws in the module)
Definition: Calibrator.h:136
readCCLHist.float
float
Definition: readCCLHist.py:83
Calibrator::GetBinnedRt
std::string GetBinnedRt(const std::string &)
Definition: Calibrator.cxx:514
RtGraph::bintype
bintype
Definition: Calibrator.h:51
Calibrator::data
std::map< std::string, caldata > data
A map between the sub-module identifier string and the calibration data structure (caldata)
Definition: Calibrator.h:339
fitman.rho
rho
Definition: fitman.py:532
RtGraph::m_t
float m_t
Definition: Calibrator.h:62
Calibrator::GetNDirs
int GetNDirs()
...
Definition: Calibrator.cxx:502
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37
Generated on Tue Jul 2 2024 21:11:23 for ATLAS Offline Software by   doxygen 1.8.18