CbmRoot
PairAnalysisHelper.cxx
Go to the documentation of this file.
1 // //
3 // helper functions wrapped in a namespace.
4 //
5 //
6 // Authors:
7 // * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
8 // Julian Book <Julian.Book@cern.ch>
9 //
10 // //
12 
13 
14 #include <TDatabasePDG.h>
15 #include <TError.h>
16 #include <TF1.h>
17 #include <TFormula.h>
18 #include <TGraph.h>
19 #include <TH1.h>
20 #include <TH2.h>
21 #include <THashList.h>
22 #include <TMCProcess.h>
23 #include <TMath.h>
24 #include <TObjArray.h>
25 #include <TObjString.h>
26 #include <TPad.h>
27 #include <TParticlePDG.h>
28 #include <TProfile.h>
29 #include <TRandom.h>
30 #include <TVectorD.h>
31 
32 
33 #include "CbmModuleList.h"
34 #include "PairAnalysisHelper.h"
35 #include "PairAnalysisStyler.h"
36 #include "PairAnalysisVarManager.h"
37 
38 //_____________________________________________________________________________
39 TVectorD*
40 PairAnalysisHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) {
41  //
42  // Make logarithmic binning
43  // the user has to delete the array afterwards!!!
44  //
45 
46  //check limits
47  if (xmin < 1e-20 || xmax < 1e-20) {
48  Error("PairAnalysisHelper::MakeLogBinning",
49  "For Log binning xmin and xmax must be > 1e-20. Using linear binning "
50  "instead!");
51  return PairAnalysisHelper::MakeLinBinning(nbinsX, xmin, xmax);
52  }
53  if (xmax < xmin) {
54  Double_t tmp = xmin;
55  xmin = xmax;
56  xmax = tmp;
57  }
58  TVectorD* binLim = new TVectorD(nbinsX + 1);
59  Double_t first = xmin;
60  Double_t last = xmax;
61  Double_t expMax = TMath::Log(last / first);
62  for (Int_t i = 0; i < nbinsX + 1; ++i) {
63  (*binLim)[i] = first * TMath::Exp(expMax / nbinsX * (Double_t) i);
64  }
65  return binLim;
66 }
67 
68 //_____________________________________________________________________________
69 TVectorD*
70 PairAnalysisHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) {
71  //
72  // Make linear binning
73  // the user has to delete the array afterwards!!!
74  //
75  if (xmax < xmin) {
76  Double_t tmp = xmin;
77  xmin = xmax;
78  xmax = tmp;
79  }
80  TVectorD* binLim = new TVectorD(nbinsX + 1);
81  Double_t first = xmin;
82  Double_t last = xmax;
83  Double_t binWidth = (last - first) / nbinsX;
84  for (Int_t i = 0; i < nbinsX + 1; ++i) {
85  (*binLim)[i] = first + binWidth * (Double_t) i;
86  }
87  return binLim;
88 }
89 
90 //_____________________________________________________________________________
91 TVectorD* PairAnalysisHelper::MakeArbitraryBinning(const char* bins) {
92  //
93  // Make arbitrary binning, bins separated by a ','
94  //
95  TString limits(bins);
96  if (limits.IsNull()) {
97  Error("PairAnalysisHelper::MakeArbitraryBinning",
98  "Bin Limit string is empty, cannot add the variable");
99  return 0x0;
100  }
101 
102  TObjArray* arr = limits.Tokenize(",");
103  Int_t nLimits = arr->GetEntries();
104  if (nLimits < 2) {
105  Error("PairAnalysisHelper::MakeArbitraryBinning",
106  "Need at leas 2 bin limits, cannot add the variable");
107  delete arr;
108  return 0x0;
109  }
110 
111  TVectorD* binLimits = new TVectorD(nLimits);
112  for (Int_t iLim = 0; iLim < nLimits; ++iLim) {
113  (*binLimits)[iLim] =
114  (static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
115  }
116 
117  delete arr;
118  return binLimits;
119 }
120 
121 //_____________________________________________________________________________
122 TVectorD* PairAnalysisHelper::MakeGausBinning(Int_t nbinsX,
123  Double_t mean,
124  Double_t sigma) {
125  //
126  // Make gaussian binning
127  // the user has to delete the array afterwards!!!
128  //
129 
130  TVectorD* binLim = new TVectorD(nbinsX + 1);
131 
132  TF1 g("g", "gaus", mean - 5 * sigma, mean + 5 * sigma);
133  g.SetParameters(1, mean, sigma);
134  Double_t sum = g.Integral(mean - 5 * sigma, mean + 5 * sigma);
135  //printf("full integral gaussian: %f \n",sum);
136 
137  TF1 g2("g2", "gaus(0)/[3]", mean - 5 * sigma, mean + 5 * sigma);
138  g2.SetParameters(1, mean, sigma, sum);
139 
140  // Double_t *params=g2.GetParameters();
141 
142  Double_t epsilon = sigma / 10000; // step size
143  Double_t xt = mean - 5 * sigma; // bin limit
144  Double_t pint = 0.0; // previous integral
145 
146  Int_t bin = 0; // current entry
147  Double_t lim = epsilon; // requested integral values (start,..., end values)
148 
149  // calculate intergral until you found all limits
150  while ((xt += epsilon) <= mean + 5 * sigma) {
151 
152  // current integral
153  // Double_t cint = g2.Integral(mean-5*sigma,xt,params,epsilon); //fast, but NOT root 6 (Integral(min,max,epsilon)) compatible
154  Double_t cint =
155  g2.Integral(mean - 5 * sigma, xt); //,params,epsilon); //slow
156  // printf(" integral to %f: %f , search limit: %f \n",xt,cint,lim);
157 
159  if (cint >= lim && pint < lim) {
160  // printf(" %d bin found for %f with requested integral %f, actual integral is %f \n",bin,xt,lim,cint);
162  (*binLim)[bin] = xt;
163 
165  bin++;
166  lim = bin * (1. / nbinsX);
168  if (bin == nbinsX) lim = 1. - epsilon;
169  }
170 
172  pint = cint;
173  }
174 
175  // binLim->Print();
176  return binLim;
177 }
178 
179 //_____________________________________________________________________________
180 TVectorD* PairAnalysisHelper::CombineBinning(TVectorD* low, TVectorD* high) {
181  //
182  // Make a new combined binning of "low" and "high"
183  // the user has to delete the returned array afterwards!!!
184  //
185 
186  // fill final vector
187  Int_t lastEqualsFirst =
188  (Int_t)(TMath::Abs((*low)[low->GetNrows() - 1] - (*high)[0]) < 1.e-15);
189  TVectorD* binLim =
190  new TVectorD(low->GetNrows() + high->GetNrows() - lastEqualsFirst);
191  for (Int_t i = 0; i < low->GetNrows(); i++)
192  (*binLim)[i] = (*low)[i];
193  for (Int_t i = 0; i < high->GetNrows() - lastEqualsFirst; i++)
194  (*binLim)[low->GetNrows() + i] = (*high)[i + lastEqualsFirst];
195 
196  // clean up
197  // delete low; delete high;
198  // binLim->Print();
199  return binLim;
200 }
201 
202 //_____________________________________________________________________________
203 TArrayD* PairAnalysisHelper::MakeStatBinLimits(TH1* h, Double_t stat) {
204  //
205  // get bin limits for stat. error less than 'stat'
206  //
207  if (!h || stat > 1.) return 0x0;
208 
209  Double_t cont = 0.0;
210  Double_t err = 0.0;
211  Double_t from = h->GetBinLowEdge(1);
212  Double_t to = 0.0;
213 
214  TArrayD* vBins = new TArrayD(1 + 1);
215  vBins->AddAt(from, 0);
216  Int_t vEle = 1;
217 
218  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
219 
220  if (h->GetBinContent(i) == 0.0 && h->GetBinError(i) == 0.) continue;
221 
222  to = h->GetBinLowEdge(i + 1);
223  cont += h->GetBinContent(i);
224  err += (h->GetBinError(i) * h->GetBinError(i));
225  vBins->AddAt(to, vEle);
226 
227  // Printf("cont %f/%f=%f err %f(%f) sum of -> rel err %f%% (current: %f%%)",
228  // h->GetBinContent(i),h->Integral(),h->GetBinContent(i)/h->Integral(), h->GetBinError(i), TMath::Sqrt(h->GetBinContent(i)),h->GetBinError(i)/h->GetBinContent(i)*100, TMath::Sqrt(err)/cont*100);
229 
230  // check for new bin
231  if (
232  TMath::Sqrt(err) / cont <= stat
233  // || (h->GetBinContent(i)/h->Integral()) < 0.01
234  // || (h->GetBinContent(i)==0.0 && h->GetBinError(i)==0. && vEle==1)
235  ) {
236  // Printf("bin from %f to %f with err %f%%",from,to,TMath::Sqrt(err)/cont*100);
237  err = 0.0;
238  cont = 0.0;
239  vEle++;
240  from = to;
241  vBins->Set(vEle + 1);
242  }
243  }
244 
245  vBins->AddAt(h->GetXaxis()->GetXmax(), vBins->GetSize() - 1);
246 
247  // for(Int_t i=0;i<vBins->GetSize();i++) Printf("%d %f",i,vBins->At(i));
248  return vBins;
249 }
250 
251 //_____________________________________________________________________________
252 TVectorD* PairAnalysisHelper::MakePdgBinning() {
253  //
254  // Make arbitrary binning using defined PDG codes
255  //
256 
257  // array of pdgcodes stored in TDatabasePDG
258  TDatabasePDG* pdg = new TDatabasePDG(); //::Instance();
259  pdg->ReadPDGTable();
260  TIter next(pdg->ParticleList());
261  TGraph gr;
262  TParticlePDG* p;
263  Int_t i = 0;
264  while ((p = (TParticlePDG*) next())) {
265  if (TMath::Abs(p->PdgCode()) < 1e+6) {
266  // printf("%s -> %d \n",p->GetName(),p->PdgCode());
267  gr.SetPoint(i++, p->PdgCode(), 1.);
268  }
269  }
270  gr.Sort();
271  TVectorD* vec = new TVectorD(gr.GetN(), gr.GetX());
272  // vec->Print();
273  delete pdg;
274  return vec;
275 }
276 
277 //_____________________________________________________________________________
278 Double_t PairAnalysisHelper::EvalFormula(TFormula* form,
279  const Double_t* values) {
280  //
281  // evaluate return value for formula with current parameter values
282  //
283  Double_t params[10] = {0.};
284  //put parameter values into array using variables stored as the parameters
285  for (Int_t ip = 0; ip < form->GetNpar(); ip++)
286  params[ip] = values[(UInt_t) form->GetParameter(ip)];
287  return (form->EvalPar(0x0, params));
288 }
289 
290 //_____________________________________________________________________________
291 TString PairAnalysisHelper::GetFormulaTitle(TFormula* form) {
292  //
293  // evaluate the formula in a readable way (for axis etc)
294  //
295  // TODO: add units, switch to TMathText, get ride of obsolete parentheses
296  TString tform(form->GetExpFormula());
297  // TMathText
298  // tform.ReplaceAll("*","\\cdot"); // multiplication sign
299  // tform.ReplaceAll("TMath::","\\"); // get ride of TMath::
300  // TLatex
301  tform.ReplaceAll("*", "#upoint"); // multiplication sign
302  // tform.ReplaceAll("Sqrt","sqrt"); // sqrt sign
303  tform.ReplaceAll("TMath::", ""); // get ride of TMath::
304  tform.ReplaceAll("()", ""); // remove function parenthesis
305  tform.ToLower(); // lower cases (e.g. Cos -> cos)
306  // replace parameter variables with proper labels
307  for (Int_t j = 0; j < form->GetNpar(); j++)
308  tform.ReplaceAll(
309  Form("[%d]", j),
310  PairAnalysisVarManager::GetValueLabel((UInt_t) form->GetParameter(j)));
311  return (tform);
312 }
313 
314 //_____________________________________________________________________________
315 TString PairAnalysisHelper::GetFormulaName(TFormula* form) {
316  //
317  // build formula key with parameter names
318  //
319  TString tform("f(");
320  for (Int_t j = 0; j < form->GetNpar(); j++) {
321  tform +=
322  PairAnalysisVarManager::GetValueName((UInt_t) form->GetParameter(j));
323  if (j != form->GetNpar() - 1) tform += ",";
324  }
325  tform += ")";
326  return (tform);
327 }
328 
329 //_____________________________________________________________________________
330 TFormula* PairAnalysisHelper::GetFormula(const char* name,
331  const char* formula) {
332  //
333  // build a TFormula object
334  //
335  TString check(formula);
336  if (check.IsNull()) return 0x0;
337  TFormula* form = new TFormula(name, formula);
338  // compile function
339  if (form->Compile()) return 0x0;
340  //set parameter/variable identifier
341  for (Int_t i = 0; i < form->GetNpar(); i++) {
342  form->SetParName(
343  i, PairAnalysisVarManager::GetValueName(form->GetParameter(i)));
344  // fUsedVars->SetBitNumber((Int_t)form->GetParameter(i),kTRUE);
345  }
346  return form;
347 }
348 
349 //_____________________________________________________________________________
350 void PairAnalysisHelper::SetPDGBinLabels(TH1* hist, Bool_t clean) {
351  //
352  // build formula key with parameter names
353  //
354 
355  TH1* hLbl = PairAnalysisStyler::GetFirstHistogram();
356  if (!hLbl) hLbl = hist;
357 
358  TDatabasePDG* pdg = TDatabasePDG::Instance();
359  TAxis* xaxis = hLbl->GetXaxis();
360  for (Int_t i = 1; i < hist->GetNbinsX() + 1; i++) {
361  // printf("bin %d: low edge: %.0f --> %s \n",i,xaxis->GetBinLowEdge(i), pdg->GetParticle(xaxis->GetBinLowEdge(i))->GetName());
362  if (clean && !hist->GetBinContent(i)) continue;
363  Int_t pdgCode = (Int_t) xaxis->GetBinLowEdge(i);
364  TParticlePDG* p = pdg->GetParticle(pdgCode);
365  // xaxis->SetBinLabel(i,(p?p->GetName():""));
366  xaxis->SetBinLabel(i, (p ? GetPDGlabel(pdgCode) : ""));
367  }
368 
369  if (hLbl->GetDimension() == 1) return;
370 
371  TAxis* yaxis = hLbl->GetYaxis();
372  TString keyY = yaxis->GetName();
373  if (keyY.Contains("pdg", TString::kIgnoreCase)) {
374  for (Int_t i = 1; i < hist->GetNbinsY() + 1; i++) {
375  // printf("bin %d: low edge: %.0f --> %s \n",i,yaxis->GetBinLowEdge(i), pdg->GetParticle(yaxis->GetBinLowEdge(i))->GetName());
376  if (clean && !hist->GetBinContent(i)) continue;
377  Int_t pdgCode = (Int_t) yaxis->GetBinLowEdge(i);
378  TParticlePDG* p = pdg->GetParticle(pdgCode);
379  // yaxis->SetBinLabel(i,(p?p->GetName():""));
380  yaxis->SetBinLabel(i, (p ? GetPDGlabel(pdgCode) : ""));
381  }
382  }
383 }
384 
385 //_____________________________________________________________________________
386 TString PairAnalysisHelper::GetPDGlabel(Int_t pdg) {
387  //
388  // return the label in latex format corresponding to pdg code
389  //
390 
391  TString name = TDatabasePDG::Instance()->GetParticle(pdg)->GetName();
392  name.ReplaceAll("dd_1_bar", "primary");
393  name.ReplaceAll("proton", "p");
394  // correct greek letters
395  /*
396  if(name.Contains("delta",TString::kIgnoreCase) ||
397  name.Contains("gamma",TString::kIgnoreCase) ||
398  name.Contains("sigma",TString::kIgnoreCase) ||
399  name.Contains("xi",TString::kIgnoreCase) ||
400  name.Contains("lambda",TString::kIgnoreCase) ||
401  name.Contains("omega",TString::kIgnoreCase) ||
402  name.Contains("eta",TString::kIgnoreCase) ||
403  name.Contains("tau",TString::kIgnoreCase) ||
404  name.Contains("phi",TString::kIgnoreCase) ||
405  name.Contains("eta",TString::kIgnoreCase) ||
406  name.Contains("upsilon",TString::kIgnoreCase) ||
407  name.Contains("nu",TString::kIgnoreCase) ||
408  name.Contains("pi",TString::kIgnoreCase) ||
409  name.Contains("rho",TString::kIgnoreCase) ) name.Prepend("#");
410  */
411  name.ReplaceAll("delta", "#delta");
412  name.ReplaceAll("gamma", "#gamma");
413  name.ReplaceAll("psi", "#psi");
414  name.ReplaceAll("sigma", "#sigma");
415  name.ReplaceAll("xi", "#xi");
416  name.ReplaceAll("lambda", "#lambda");
417  name.ReplaceAll("omega", "#omega");
418  name.ReplaceAll("eta", "#eta");
419  name.ReplaceAll("tau", "#tau");
420  name.ReplaceAll("phi", "#phi");
421  name.ReplaceAll("upsilon", "#upsilon");
422  name.ReplaceAll("nu", "#nu");
423  name.ReplaceAll("mu", "#mu");
424  name.ReplaceAll("pi", "#pi");
425  name.ReplaceAll("rho", "#rho");
426 
427  // correct anti particles
428  if (name.Contains("anti")) {
429  name.ReplaceAll("anti", "#bar{");
430  name.Append("}");
431  } else if (name.Contains("_bar")) {
432  name.ReplaceAll("_bar", "}");
433  name.Prepend("#bar{");
434  }
435  // correct indices
436  name.ReplaceAll("+", "^{+}");
437  name.ReplaceAll("-", "^{-}");
438  name.ReplaceAll("0", "^{0}");
439  name.ReplaceAll("_s", "_{s}");
440  name.ReplaceAll("_c", "_{c}");
441  name.ReplaceAll("_b", "_{b}");
442  name.ReplaceAll("_1", "_{1}");
443  // specials
444  name.ReplaceAll("/psi", "/#psi");
445  // Printf(" %d = %s",pdg,name.Data());
446  return name;
447 }
448 
449 //_____________________________________________________________________________
450 void PairAnalysisHelper::SetGEANTBinLabels(TH1* hist) {
451  //
452  // build formula key with parameter names
453  //
454  TAxis* xaxis = hist->GetXaxis();
455  for (Int_t i = 1; i < hist->GetNbinsX() + 1; i++) {
456  xaxis->SetBinLabel(i, TMCProcessName[i - 1]);
457  }
458  xaxis->LabelsOption("v"); // vertical labels
459  if (gPad) {
460  xaxis->SetTitleOffset(3.6);
461  gPad->SetTopMargin(0.01);
462  gPad->SetBottomMargin(0.4);
463  }
464 }
465 
466 
467 //_____________________________________________________________________________
468 TString PairAnalysisHelper::GetDetName(ECbmModuleId det) {
469  //
470  // get detector name
471  //
472  if (det == static_cast<ECbmModuleId>(9)) return ("");
473  TString name = CbmModuleList::GetModuleNameCaps(det);
474 
475  // TString name = CbmModuleList::GetModuleName(det);
476  return (name);
477 }
478 
479 //_____________________________________________________________________________
480 Double_t PairAnalysisHelper::GetContentMinimum(TH1* h, Bool_t inclErr) {
481  //
482  // get minimum bin content of histogram (having entries)
483  //
484  Int_t bin, binx, biny, binz;
485  Int_t xfirst = h->GetXaxis()->GetFirst();
486  Int_t xlast = h->GetXaxis()->GetLast();
487  Int_t yfirst = h->GetYaxis()->GetFirst();
488  Int_t ylast = h->GetYaxis()->GetLast();
489  Int_t zfirst = h->GetZaxis()->GetFirst();
490  Int_t zlast = h->GetZaxis()->GetLast();
491  Double_t minimum = FLT_MAX, value = 0., error = 0.;
492  for (binz = zfirst; binz <= zlast; binz++) {
493  for (biny = yfirst; biny <= ylast; biny++) {
494  for (binx = xfirst; binx <= xlast; binx++) {
495  bin = h->GetBin(binx, biny, binz);
496  value = h->GetBinContent(bin);
497  error = h->GetBinError(bin);
498  // Printf(" \t hist%s bin%d value%f error%f \n",h->GetTitle(),bin,value,error);
499  if (gPad->GetLogy() && (value - error) <= 0.) continue;
500  if (error > value * 0.9) continue;
501  if (inclErr) value -= h->GetBinError(bin);
502  if (value < minimum
503  && TMath::Abs(h->GetBinError(bin) - 1.e-15) > 1.e-15) {
504  minimum = value;
505  }
506  }
507  }
508  }
509  // Printf(" RETURN VALUE: hist%s %f \n",h->GetTitle(),minimum);
510  return minimum;
511 }
512 
513 Double_t PairAnalysisHelper::GetContentMaximum(TH1* h, Bool_t inclErr) {
514  //
515  // get maximum bin content+error of histogram (having entries)
516  //
517  Int_t bin, binx, biny, binz;
518  Int_t xfirst = h->GetXaxis()->GetFirst();
519  Int_t xlast = h->GetXaxis()->GetLast();
520  Int_t yfirst = h->GetYaxis()->GetFirst();
521  Int_t ylast = h->GetYaxis()->GetLast();
522  Int_t zfirst = h->GetZaxis()->GetFirst();
523  Int_t zlast = h->GetZaxis()->GetLast();
524  Double_t maximum = -1. * FLT_MAX, value = 0., error = 0.;
525  for (binz = zfirst; binz <= zlast; binz++) {
526  for (biny = yfirst; biny <= ylast; biny++) {
527  for (binx = xfirst; binx <= xlast; binx++) {
528  bin = h->GetBin(binx, biny, binz);
529  value = h->GetBinContent(bin);
530  error = h->GetBinError(bin);
531  if (inclErr) value += h->GetBinError(bin);
532  if (value > maximum && TMath::Abs(error - 1.e-15) > 1.e-15) {
533  maximum = value;
534  }
535  }
536  }
537  }
538  return maximum;
539 }
540 
541 Double_t PairAnalysisHelper::GetQuantile(TH1* h1, Double_t p /*=0.5*/) {
542  //
543  // calculates the quantile of the bin contents, p=0.5 -> Median
544  // useful functionallity for plotting 2D distibutions with some extreme outliers
545  //
546  if (p < 0.0 || p > 1.) return -1.;
547  Int_t nbinsX = h1->GetNbinsX();
548  Int_t nbinsY = h1->GetNbinsY();
549  Int_t nbinsZ = h1->GetNbinsZ();
550  Int_t nbins = (nbinsX * (nbinsY ? nbinsY : 1) * (nbinsZ ? nbinsZ : 1));
551  Int_t xbin = -1;
552  Int_t ybin = -1;
553  Int_t zbin = -1;
554  Double_t val[nbins];
555  Int_t idx[nbins];
556  Int_t nfilled = 0;
557  for (Int_t i = 1; i <= nbins; i++) {
558  h1->GetBinXYZ(i, xbin, ybin, zbin);
559  if (xbin < h1->GetXaxis()->GetFirst() || xbin > h1->GetXaxis()->GetLast())
560  continue;
561  if (ybin < h1->GetYaxis()->GetFirst() || ybin > h1->GetYaxis()->GetLast())
562  continue;
563  Double_t con = h1->GetBinContent(i);
564  Double_t err = h1->GetBinError(i);
565  if (err != 0.0) {
566  // printf("bin%d %.f+-%.f \n",i,con,err);
567  val[nfilled] =
568  con + (h1->GetDimension() < 2 ? err : 0.0); // w or w/o err?
569  nfilled++;
570  }
571  }
572  if (nfilled == 0) return -1.;
573  TMath::Sort(nfilled, val, idx, kFALSE); // kFALSE=increasing numbers
574  Int_t pos = (Int_t)((Double_t) nfilled * p);
575  //for(int i=0; i<nfilled; i++) cout << i << " " << idx[i] << " " << val[idx[i]] << endl;
576  //printf("nfilled %d quantile %f pos %d: %f \n",nfilled,p,pos,val[idx[pos]]);
577  return val[idx[pos]];
578 }
579 
580 void PairAnalysisHelper::NormalizeSlicesY(TH2* h) {
581  //
582  // normalize slices along Y in case of a 2D histogram
583  //
584  TH2* hsum = (TH2*) h->Clone("SliceInts");
585  hsum->Reset("CE");
586  for (Int_t ix = 0; ix <= hsum->GetNbinsX() + 1; ix++) {
587  Double_t ysum = h->Integral(ix, ix);
588  for (Int_t iy = 0; iy <= hsum->GetNbinsY() + 1; iy++) {
589  hsum->SetBinContent(ix, iy, ysum);
590  }
591  }
593  h->Divide(hsum);
594  delete hsum;
595 }
596 
597 void PairAnalysisHelper::CumulateSlicesX(TH2* h, Bool_t reverse, Bool_t norm) {
598  //
599  // caluclate cumulative sum of bins (normalized to one)
600  // for slices along X in case of a 2D histogram
601  //
603  Int_t xstart = (reverse ? h->GetNbinsX() + 1 : 0 - 1);
604  Int_t xend = (reverse ? 0 : h->GetNbinsX() + 1 + 1);
605  Int_t xincr = (reverse ? -1 : +1);
606 
607  Double_t integral = 1.;
608  for (Int_t iy = 0; iy <= h->GetNbinsY() + 1; iy++) {
609  if (norm) integral = h->Integral(0, h->GetNbinsX() + 1, iy, iy);
610 
611  Double_t cumInt = 0;
612  for (Int_t ix = xstart; ix != xend; ix += xincr) {
613  cumInt += h->GetBinContent(ix, iy);
614  h->SetBinContent(ix, iy, cumInt / integral);
615  }
616  }
617 }
618 
619 void PairAnalysisHelper::Cumulate(TH1* h, Bool_t reverse, Bool_t norm) {
620  //
621  // caluclate cumulative sum of bins (normalized to one)
622  //
624  Int_t xstart = (reverse ? h->GetNbinsX() + 1 : 0 - 1);
625  Int_t xend = (reverse ? 0 : h->GetNbinsX() + 1 + 1);
626  Int_t xincr = (reverse ? -1 : +1);
627 
628  Double_t integral = 1.;
629  if (norm) integral = h->Integral(0, h->GetNbinsX() + 1);
630 
631  Double_t cumInt = 0;
632  for (Int_t ix = xstart; ix != xend; ix += xincr) {
633  cumInt += h->GetBinContent(ix);
634  h->SetBinContent(ix, cumInt / integral);
635  }
636 }
637 
638 
639 TObject* PairAnalysisHelper::FindObjectByTitle(TObjArray* arrhist,
640  TString ref) {
641  //
642  // shortcut to find a certain pair type object in array
643  //
644  for (Int_t i = 0; i < arrhist->GetEntriesFast(); i++) {
645  if (!ref.CompareTo(arrhist->UncheckedAt(i)->GetTitle())) {
646  return arrhist->UncheckedAt(i);
647  }
648  }
649  return 0x0;
650 
651  // return ( arrhist->FindObject(Form("Pair.%s",PairAnalysis::PairClassName(type))) );
652  // TString ref=Form("Pair.%s",PairAnalysis::PairClassName(type));
653 }
654 
655 
656 //_____________________________________________________________________________
657 Int_t PairAnalysisHelper::GetPrecision(Double_t value) {
658  //
659  // computes the precision of a double
660  // usefull for axis ranges etc
661 
662  Bool_t bfnd = kFALSE;
663  Int_t precision = 0;
664  value *= 1e6;
665  while (!bfnd) {
666  // Printf(" value %f precision %d bfnd %d",TMath::Abs(value*TMath::Power(10,precision)), precision, bfnd);
667  bfnd =
668  ((TMath::Abs(value * TMath::Power(10, precision)) / 1e6
669  - TMath::Floor(TMath::Abs(value * TMath::Power(10, precision)) / 1e6))
670  != 0.0
671  ? kFALSE
672  : kTRUE);
673  if (!bfnd) precision++;
674  }
675 
676  // Printf("precision for %f found to be %d", value, precision);
677  return precision;
678 }
PairAnalysisHelper::EvalFormula
Double_t EvalFormula(TFormula *form, const Double_t *values)
Definition: PairAnalysisHelper.cxx:278
PairAnalysisHelper.h
PairAnalysisHelper::SetGEANTBinLabels
void SetGEANTBinLabels(TH1 *hist)
Definition: PairAnalysisHelper.cxx:450
PairAnalysisHelper::GetDetName
TString GetDetName(ECbmModuleId det)
Definition: PairAnalysisHelper.cxx:468
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
gr
Char_t * gr
Definition: CbmEvDisTracks.cxx:289
PairAnalysisHelper::NormalizeSlicesY
void NormalizeSlicesY(TH2 *h)
Definition: PairAnalysisHelper.cxx:580
ECbmModuleId
ECbmModuleId
Definition: CbmDefs.h:33
PairAnalysisHelper::MakeStatBinLimits
TArrayD * MakeStatBinLimits(TH1 *h, Double_t stat)
Definition: PairAnalysisHelper.cxx:203
PairAnalysisVarManager::GetValueName
static const char * GetValueName(Int_t i)
Definition: PairAnalysisVarManager.h:377
PairAnalysisHelper::MakeLinBinning
TVectorD * MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
Definition: PairAnalysisHelper.cxx:70
PairAnalysisHelper::GetPrecision
Int_t GetPrecision(Double_t value)
Definition: PairAnalysisHelper.cxx:657
PairAnalysisHelper::MakeLogBinning
TVectorD * MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
Definition: PairAnalysisHelper.cxx:40
PairAnalysisHelper::GetFormulaName
TString GetFormulaName(TFormula *form)
Definition: PairAnalysisHelper.cxx:315
PairAnalysisHelper::GetFormula
TFormula * GetFormula(const char *name, const char *formula)
Definition: PairAnalysisHelper.cxx:330
CbmModuleList::GetModuleNameCaps
static TString GetModuleNameCaps(ECbmModuleId moduleId)
Definition: CbmModuleList.cxx:77
PairAnalysisVarManager.h
h
Data class with information on a STS local track.
PairAnalysisHelper::CombineBinning
TVectorD * CombineBinning(TVectorD *low, TVectorD *high)
Definition: PairAnalysisHelper.cxx:180
PairAnalysisVarManager::GetValueLabel
static const char * GetValueLabel(Int_t i)
Definition: PairAnalysisVarManager.h:380
PairAnalysisHelper::MakeGausBinning
TVectorD * MakeGausBinning(Int_t nbinsX, Double_t mean, Double_t sigma)
Definition: PairAnalysisHelper.cxx:122
PairAnalysisHelper::CumulateSlicesX
void CumulateSlicesX(TH2 *h, Bool_t reverse=kFALSE, Bool_t norm=kFALSE)
Definition: PairAnalysisHelper.cxx:597
CbmModuleList.h
first
bool first
Definition: LKFMinuit.cxx:143
PairAnalysisHelper::FindObjectByTitle
TObject * FindObjectByTitle(TObjArray *arrhist, TString ref)
Definition: PairAnalysisHelper.cxx:639
PairAnalysisHelper::MakePdgBinning
TVectorD * MakePdgBinning()
Definition: PairAnalysisHelper.cxx:252
PairAnalysisStyler.h
PairAnalysisHelper::GetPDGlabel
TString GetPDGlabel(Int_t pdg)
Definition: PairAnalysisHelper.cxx:386
PairAnalysisStyler::GetFirstHistogram
TH1 * GetFirstHistogram()
Definition: PairAnalysisStyler.cxx:605
PairAnalysisHelper::MakeArbitraryBinning
TVectorD * MakeArbitraryBinning(const char *bins)
Definition: PairAnalysisHelper.cxx:91
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
PairAnalysisHelper::Cumulate
void Cumulate(TH1 *h, Bool_t reverse=kFALSE, Bool_t norm=kFALSE)
Definition: PairAnalysisHelper.cxx:619
PairAnalysisHelper::GetContentMaximum
Double_t GetContentMaximum(TH1 *h, Bool_t inclErr=kTRUE)
Definition: PairAnalysisHelper.cxx:513
PairAnalysisHelper::SetPDGBinLabels
void SetPDGBinLabels(TH1 *hist, Bool_t clean)
Definition: PairAnalysisHelper.cxx:350
PairAnalysisHelper::GetContentMinimum
Double_t GetContentMinimum(TH1 *h, Bool_t inclErr=kTRUE)
Definition: PairAnalysisHelper.cxx:480
PairAnalysisHelper::GetQuantile
Double_t GetQuantile(TH1 *h1, Double_t p=0.5)
Definition: PairAnalysisHelper.cxx:541
PairAnalysisHelper::GetFormulaTitle
TString GetFormulaTitle(TFormula *form)
Definition: PairAnalysisHelper.cxx:291