AliPhysics  cda3415 (cda3415)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliJetFlowTools.h
Go to the documentation of this file.
1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2  * See cxx source for full Copyright notice */
3  /* $Id$ */
4 
5 /* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
6  * see implementation for additional information */
7 
8 #ifndef AliJetFlowTools_H
9 #define AliJetFlowTools_H
10 
11 // root forward declarations
12 class TH1;
13 class TF2;
14 class TH2D;
15 class TCanvas;
16 class TString;
17 class TArrayD;
18 class TGraph;
19 class TGraphErrors;
20 class TObjArray;
21 // aliroot forward declarations
22 class AliAnaChargedJetResponseMaker;
23 class AliUnfolding;
24 // root includes
25 #include "TRandom3.h"
26 #include "TMatrixD.h"
27 #include "TList.h"
28 #include "TDirectoryFile.h"
29 #include "TFile.h"
30 #include "TProfile.h"
31 #include "TVirtualPad.h"
32 #include "TPaveText.h"
33 #include "TLegend.h"
34 #include "TLatex.h"
35 #include "TF1.h"
36 #include "TH1D.h"
37 #include "TMath.h"
38 // define the following variable to build with debug flags
39 // #define ALIJETFLOWTOOLS_DEBUG_FLAG
40 
41 //_____________________________________________________________________________
43  public:
45  protected:
46  ~AliJetFlowTools(); // not implemented (deliberately). object ownership is a bit messy in this class
47  // since most (or all) of the objects are owned by the input and output files
48  public:
49  // enumerators
50  enum unfoldingAlgorithm { // type of unfolding alrogithm
51  kChi2, // chi^2 unfolding, implemented in AliUnfolding
52  kBayesian, // Bayesian unfolding, implemented in RooUnfold
53  kBayesianAli, // Bayesian unfolding, implemented in AliUnfolding
54  kSVD, // SVD unfolding, implemented in RooUnfold
55  kFold, // instead of unfolding, fold the input with the response
56  kNone }; // no unfolding
57  enum prior { // prior that is used for unfolding
58  kPriorChi2, // prior from chi^2 method
59  kPriorMeasured, // use measured spectrum as prior
60  kPriorPythia, // use pythia spectrum as prior
61  kPriorTF1 }; // use properly binned TF1 as prior
62  enum histoType { // histogram identifier, only used internally
63  kInPlaneSpectrum, // default style for spectrum
68  kBar, // default style for bar histogram
69  kRatio, // default style for ratio
70  kV2, // default style for v2
71  kDeltaPhi, // default for delta phi
72  kEmpty }; // default style
73  // setters, interface to the class
78  void SetPwrtTo(Float_t p) {gPwrtTo = p;}
80  gPwrtToArray = a; // set the array
81  gPwrtTo = -9999; // tells instance to not use this value
82  }
84  void SetPivot(Float_t p) {fPivot = p;}
85  void SetConstantUE(Bool_t ue) {fConstantUE = ue;}
87  void SetSaveFull(Bool_t b) {fSaveFull = b;}
89  fInputList = list;
90  fRefreshInput = kTRUE;
91  }
94  // create a new output list and add it to the full output
95  fListPrefix++;
96  if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
97  fOutputFile->cd(); // avoid nested dirs
98  if(name.EqualTo(fActiveString)) {
99  printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
100  name+="_2";
101  }
102  fActiveString = Form("%i__%s", fListPrefix, name.Data());
103  fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
104  fActiveDir->cd();
105  }
107  // in case of one centraltiy
108  fCentralityArray = new TArrayI(1);
109  fCentralityArray->AddAt(bin, 0);
110  // for one centrality there's no need for weights
111  fCentralityWeights = new TArrayD(1);
112  fCentralityWeights->AddAt(1., 0);
113  }
115  void SetCentralityBin(TArrayI* bins) {
116  fCentralityArray = bins;
117  }
118  void SetCentralityWeight(TArrayD* weights) {
119  fCentralityWeights = weights;
120  if(!fCentralityArray) printf(" > Warning: centrality weights set, but bins are not defined! \n");
121  }
123  TList* l,
124  Int_t c,
125  Float_t weight) {
126  fMergeWithList = l;
127  fMergeWithCen = c;
128  fMergeWithWeight = weight;
129  }
132  void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
133  void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
135  void SetBinsRecPrior(TArrayD* bins) {fBinsRecPrior = bins;}
136  void SetSVDReg(Int_t r) {fSVDRegIn = r; fSVDRegOut = r;}
137  void SetSVDReg(Int_t in, Int_t out) {fSVDRegIn = in; fSVDRegOut = out;}
138  void SetSVDToy(Bool_t b, Float_t r) {fSVDToy = b; fJetRadius = r;}
139  void SetBeta(Double_t b) {fBetaIn = b; fBetaOut = b;}
140  void SetBeta(Double_t i, Double_t o) {fBetaIn = i; fBetaOut = o;}
147  void SetPrior(prior p) {fPrior = p;}
148  void SetPrior(prior p, TF1* function, TArrayD* bins) {
149  fPrior = p;
150  // set prior to value supplied in TF1
151  fPriorUser = new TH1D("prior_user", "prior_user", bins->GetSize()-1, bins->GetArray());
152  // loop over bins and fill the histo from the tf1
153  for(Int_t i(0); i < fPriorUser->GetNbinsX() + 1; i++) fPriorUser->SetBinContent(i, function->Integral(fPriorUser->GetXaxis()->GetBinLowEdge(i), fPriorUser->GetXaxis()->GetBinUpEdge(i))/fPriorUser->GetXaxis()->GetBinWidth(i));
154  }
155  void SetPrior(prior p, TH1D* spectrum) {fPrior = p; fPriorUser = spectrum;}
157  void SetNormalizeSpectra(Int_t e) { // use to normalize to this no of events
158  fEventCount = e;
159  fNormalizeSpectra = kFALSE;
160  }
161  void SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
162  fSmoothenPrior = b;
163  fFitMin = min;
164  fFitMax = max;
165  fFitStart = start;
166  fSmoothenCounts = counts;
167  }
168  void SetTestMode(Bool_t t) {fTestMode = t;}
175  void SetExLJDpt(Bool_t i) {fExLJDpt = i;}
176  void SetWeightFunction(TF1* w) {fResponseMaker->SetRMMergeWeightFunction(w);}
177  void SetRMS(Bool_t r) {fRMS = r;}
178  void SetSymmRMS(Bool_t r) {fSymmRMS = r;}
179  void SetRho0(Bool_t r) {fRho0 = r;}
180  void SetBootstrap(Bool_t b, Bool_t r = kTRUE) {
181  // note that setting this option will not lead to true resampling
182  // but rather to randomly drawing a new distribution from a pdf
183  // of the measured distribution
184  fBootstrap = b;
185  // by default fully randomize randomizer from system time
186  if(r) {
187  delete gRandom;
188  gRandom = new TRandom3(0);
189  }
190  }
191  void SetHarmonic(Int_t n) {fHarmonic = n;}
192  // main function. buffers about 5mb per call!
193  void Make(TH1* customIn = 0x0, TH1* customOut = 0x0);
194  void MakeAU(); // test function, use with caution (09012014)
195  void Finish() {
196  fOutputFile->cd();
197  if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
198  if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
199  if(fRMSRatio) fRMSRatio->Write();
200  fOutputFile->Close();}
201  void PostProcess(
202  TString def,
203  Int_t columns = 4,
204  Float_t rangeLow = 20,
205  Float_t rangeUp = 80,
206  TString in = "UnfoldedSpectra.root",
207  TString out = "ProcessedSpectra.root") const;
208  void BootstrapSpectra(
209  TString def,
210  TString in = "UnfoldedSpectra.root",
211  TString out = "BootstrapSpectra.root") const;
212  void GetNominalValues(
213  TH1D*& ratio,
214  TGraphErrors*& v2,
215  TArrayI* in,
216  TArrayI* out,
217  TString inFile = "UnfoldedSpectra.root",
218  TString outFile = "Nominal.root") const;
220  TGraphAsymmErrors*& corrRatio,
221  TGraphAsymmErrors*& corrV2,
222  TArrayI* variationsIn,
223  TArrayI* variationsOut,
224  Bool_t sym,
225  TArrayI* variantions2ndIn,
226  TArrayI* variantions2ndOut,
227  Bool_t sym2nd,
228  TString type = "",
229  TString type2 = "",
230  Int_t columns = 4,
231  Float_t rangeLow = 20,
232  Float_t rangeUp = 80,
233  Float_t corr = .5,
234  TString in = "UnfoldedSpectra.root",
235  TString out = "CorrelatedUncertainty.root") const;
236  void GetShapeUncertainty(
237  TGraphAsymmErrors*& shapeRatio,
238  TGraphAsymmErrors*& shapeV2,
239  TArrayI* regularizationIn,
240  TArrayI* regularizationOut,
241  TArrayI* recBinIn = 0x0,
242  TArrayI* recBinOut = 0x0,
243  TArrayI* methodIn = 0x0,
244  TArrayI* methodOut = 0x0,
245  Int_t columns = 4,
246  Float_t rangeLow = 20,
247  Float_t rangeUp = 80,
248  Float_t corr = .0,
249  TString in = "UnfoldedSpectra.root",
250  TString out = "ShapeUncertainty.root",
251  Bool_t regularizationOnV2 = kTRUE // get uncertainty on yields separately or v2 directly
252  ) const;
254  TH2D* detectorResponse, // detector response matrix
255  TH1D* jetPtIn, // in plane jet spectrum
256  TH1D* jetPtOut, // out of plane jet spectrum
257  TH1D* dptIn, // in plane delta pt distribution
258  TH1D* dptOut, // out of plane delta pt distribution
259  Int_t eventCount = 0); // event count (optional)
262  // static const helper functions, mainly histogram manipulation
263  static void RemoveSign(Double_t& d) {d = TMath::Abs(d);}
264  static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
265  static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
266  static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
267  static TH1* Bootstrap(TH1* hist, Bool_t kill = kTRUE);
268  static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
269  TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
270  static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
271  static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
272  static TH1D* MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr);
273  static TGraphErrors* GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
274  static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
275  static TH1D* GetV2Histo(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
276  static TH1F* ConvertGraphToHistogram(TGraphErrors* g);
278  static Double_t GetRMSOfTH1(TH1* h, Double_t a, Double_t b);
279  static TF1* GetErrorFromFit(TH1* h1, TH1* h2, Double_t a, Double_t b,
280  Float_t pivot = 50., Bool_t subdueError = kFALSE,
281  TString str = "", Bool_t setContent = kTRUE);
282  void ReplaceBins(TArrayI* array, TGraphAsymmErrors* graph);
283  void ReplaceBins(TArrayI* array, TGraphErrors* graph);
285  TH1* h1, TH1* h2, Double_t r, TString name,
286  TH1* relativeErrorInUp,
287  TH1* relativeErrorInLow,
288  TH1* relativeErrorOutUp,
289  TH1* relativeErrorOutLow,
290  Float_t rho = 0.) const;
291  static void GetSignificance(
292  TGraphErrors* n, // points with stat error
293  TGraphAsymmErrors* shape, // points with shape error
294  TGraphAsymmErrors* corr, // corr with stat error
295  Int_t low, // pt lower level
296  Int_t up // pt upper level
297  );
298  static void MinimizeChi2nd();
299  static Double_t PhenixChi2nd(const Double_t *xx );
301  static TF2* ReturnFunctionnd(Double_t &p);
302  static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
303  static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
304  static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
305  void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
306  static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
307  static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
308  // set style
310  static void Style(Bool_t legacy = kFALSE);
311  static void Style(TCanvas* c, TString style = "PEARSON");
312  static void Style(TVirtualPad* c, TString style = "SPECTRUM", Bool_t legacy = kFALSE);
313  static void Style(TLegend* l);
314  static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
315  static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
316  static TLegend* AddLegend(TVirtualPad* p, Bool_t style = kFALSE) {
317  if(!style) return p->BuildLegend(.565, .663, .882, .883);
318  else {
319  TLegend* l = AddLegend(p, kFALSE);
320  Style(l);
321  return l;
322  }
323  }
324  static TPaveText* AddTPaveText(
325  // this text appears under the logo
326  TString text,
327  Int_t r = 2,
328  Double_t a = .587,
329  Double_t b = .695,
330  Double_t c = .872,
331  Double_t d = .801) {
332  TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
333  t->SetFillColor(0);
334  t->SetBorderSize(0);
335  t->AddText(0.,0.,text.Data());
336  t->AddText(0., 0., Form("#it{R} = 0.%i anti-#it{k}_{T}, |#eta_{jet}|<%.1f", r, .9-r/10.));
337  t->SetTextColor(kBlack);
338  t->SetTextFont(42);
339  t->SetTextSize(gStyle->GetTextSize()*.8);
340  t->Draw("same");
341  return t;
342  }
343  static TPaveText* AddText(
344  TString text,
345  EColor col,
346  Double_t a = .2098,
347  Double_t b = .5601,
348  Double_t c = .613,
349  Double_t d = .6211) {
350  TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
351  t->SetFillColor(0);
352  t->SetBorderSize(0);
353  t->AddText(0.,0.,text.Data());
354  t->SetTextColor(col);
355  t->SetTextFont(42);
356  t->SetTextSize(gStyle->GetTextSize()*.8);
357  t->Draw("same");
358  return t;
359  }
360  static TLatex* AddLogo(Int_t logo, Double_t xmin = .59, Double_t ymax = .81) {
361  if(logo == 0) return AddTLatex(xmin, ymax, "ALICE");
362  else if (logo == 1) return AddTLatex(xmin, ymax, "ALICE Preliminary");
363  else if (logo == 2) return AddTLatex(xmin, ymax, "ALICE Simulation");
364  else if (logo == 3) return AddTLatex(xmin, ymax, "work in progress");
365  return 0x0;
366  }
367  static TLatex* AddSystem() {
368  return AddTLatex(0.55, 87, "Pb-Pb #sqrt{#it{s}}}_{NN} = 2.76 TeV");
369  }
370  static TLatex* AddTLatex(Double_t xmin, Double_t ymax, TString string) {
371 TLatex* tex = new TLatex(xmin, ymax, string.Data());
372  tex->SetNDC();
373  tex->SetTextFont(42);
374  tex->Draw("same");
375  return tex;
376  }
377 
378  static void SavePadToPDF(TVirtualPad* pad) {
379  if(pad) return;/*pad->SaveAs(Form("%s.pdf", pad->GetName()));*/
380  else return;}
381  // interface to AliUnfolding, not necessary but nice to have all parameters in one place
382  static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
383  static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
384  static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
385  static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
386  static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
387  private:
388  Bool_t PrepareForUnfolding(TH1* customIn = 0x0, TH1* customOut = 0x0);
390  TH1D* GetPrior( const TH1D* measuredJetSpectrum,
391  const TH2D* resizedResponse,
392  const TH1D* kinematicEfficiency,
393  const TH1D* measuredJetSpectrumTrueBins,
394  const TString suffix,
395  const TH1D* jetFindingEfficiency);
396  TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
397  const TH2D* resizedResponse,
398  const TH1D* kinematicEfficiency,
399  const TH1D* measuredJetSpectrumTrueBins,
400  const TString suffix,
401  const TH1D* jetFindingEfficiency = 0x0);
402  TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
403  const TH2D* resizedResponse,
404  const TH1D* kinematicEfficiency,
405  const TH1D* measuredJetSpectrumTrueBins,
406  const TString suffix,
407  const TH1D* jetFindingEfficiency = 0x0);
408  TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
409  const TH2D* resizedResponse,
410  const TH1D* kinematicEfficiency,
411  const TH1D* measuredJetSpectrumTrueBins,
412  const TString suffix,
413  const TH1D* jetFindingEfficiency = 0x0);
414  TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
415  const TH2D* resizedResponse,
416  const TH1D* kinematicEfficiency,
417  const TH1D* measuredJetSpectrumTrueBins,
418  const TString suffix,
419  const TH1D* jetFindingEfficiency = 0x0);
420  TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
421  const TH2D* resizedResponse,
422  const TH1D* kinematicEfficiency,
423  const TH1D* measuredJetSpectrumTrueBins,
424  const TString suffix,
425  const TH1D* jetFindingEfficiency = 0x0);
426  TH1D* FoldSpectrum( const TH1D* measuredJetSpectrum,
427  const TH2D* resizedResponse,
428  const TH1D* kinematicEfficiency,
429  const TH1D* measuredJetSpectrumTrueBins,
430  const TString suffix,
431  const TH1D* jetFindingEfficiency = 0x0);
432  void SystematicsWrapper(
433  TArrayI* variationsIn,
434  TArrayI* variationsOut,
435  TH1D*& relativeErrorInUp,
436  TH1D*& relativeErrorInLow,
437  TH1D*& relativeErrorOutUp,
438  TH1D*& relativeErrorOutLow,
439  TH1D*& relativeSystematicIn,
440  TH1D*& relativeSystematicOut,
441  TH1D*& nominal,
442  TH1D*& nominalIn,
443  TH1D*& nominalOut,
444  Int_t columns,
445  Float_t rangeLow,
446  Float_t rangeUp,
447  TFile* readMe,
448  TString source = "",
449  Bool_t RMS = kFALSE,
450  Bool_t onRatio = kTRUE) const;
452  TArrayI* variationsIn,
453  TArrayI* variationsOut,
454  TH1D*& relativeErrorInUp,
455  TH1D*& relativeErrorInLow,
456  TH1D*& relativeErrorOutUp,
457  TH1D*& relativeErrorOutLow,
458  TH1D*& relativeSystematicIn,
459  TH1D*& relativeSystematicOut,
460  TH1D*& nominal,
461  TH1D*& nominalIn,
462  TH1D*& nominalOut,
463  Int_t columns,
464  Float_t rangeLow,
465  Float_t rangeUp,
466  TFile* readMe,
467  TString source = "",
468  Bool_t RMS = kFALSE) const;
470  TArrayI* variationsIn,
471  TArrayI* variationsOut,
472  TH1D*& relativeErrorInUp,
473  TH1D*& relativeErrorInLow,
474  TH1D*& relativeErrorOutUp,
475  TH1D*& relativeErrorOutLow,
476  TH1D*& relativeSystematicIn,
477  TH1D*& relativeSystematicOut,
478  TH1D*& nominal,
479  TH1D*& nominalIn,
480  TH1D*& nominalOut,
481  Int_t columns,
482  Float_t rangeLow,
483  Float_t rangeUp,
484  TFile* readMe,
485  TString source = "",
486  Bool_t RMS = kFALSE) const;
487 
488  static void ResetAliUnfolding();
489  static void SquelchWarning() {
490  printf(" >> I squelched a warning, jay, I'm contributing ! << \n");
491  return;
492  }
493  // give object a unique name via the 'protect heap' functions.
494  // may seem redundant, but some internal functions of root (e.g.
495  // ProjectionY()) check for existing objects by name and re-use them
496  TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
497  TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
498  TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
499  // members, accessible via setters
500  Int_t fListPrefix; // prefix for list readability
501  AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
502  Bool_t fRMS; // systematic method
503  Bool_t fSymmRMS; // symmetric systematic method
504  Bool_t fConstantUE; // assign a constant unfolding error
505  Bool_t fRho0; // use the result obtained with the 'classic' fixed rho
506  Bool_t fBootstrap; // use bootstrap resampling of input data
507  TF1* fPower; // smoothening fit
508  Bool_t fSaveFull; // save all generated histograms to file
509  TString fActiveString; // identifier of active output
510  TDirectoryFile* fActiveDir; // active directory
511  TList* fInputList; // input list
512  Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
513  TString fOutputFileName; // output file name
514  TFile* fOutputFile; // output file
515  TArrayI* fCentralityArray; // array of centrality bins that are merged
516  TArrayI* fMergeBinsArray; // array of pt bins that are merged
517  TArrayD* fCentralityWeights; // array of centrality weights
518  TList* fMergeWithList; // additional input list
519  Int_t fMergeWithCen; // centrality bin of additional input list
520  Float_t fMergeWithWeight; // weight of additional input list
521  TH2D* fDetectorResponse; // detector response
522  TH1D* fJetFindingEff; // jet finding efficiency
523  Double_t fBetaIn; // regularization strength, in plane unfolding
524  Double_t fBetaOut; // regularization strength, out of plane unfoldign
525  Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
526  Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
527  Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
528  Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
529  Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
530  unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
531  prior fPrior; // prior for unfolding
532  TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
533  TArrayD* fBinsTrue; // pt true bins
534  TArrayD* fBinsRec; // pt rec bins
535  TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
536  TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
537  Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
538  Int_t fSVDRegOut; // svd regularization out of plane
539  Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
540  Float_t fJetRadius; // jet radius (for SVD toy)
541  Int_t fEventCount; // number of events
542  Bool_t fNormalizeSpectra; // normalize spectra to event count
543  Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
544  Float_t fFitMin; // lower bound of smoothening fit
545  Float_t fFitMax; // upper bound of smoothening fit
546  Float_t fFitStart; // from this value, use smoothening
547  Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
548  Bool_t fTestMode; // unfold with unity response for testing
549  Bool_t fRawInputProvided; // input histograms provided, not read from file
550  Double_t fEventPlaneRes; // event plane resolution for current centrality
551  Bool_t fUseDetectorResponse; // add detector response to unfolding
552  Bool_t fUseDptResponse; // add dpt response to unfolding
553  Bool_t fTrainPower; // don't clear the params of fPower for call to Make
554  // might give more stable results, but possibly introduces
555  // a bias / dependency on previous iterations
556  Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
557  Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
558  Bool_t fExLJDpt; // exclude randon cones with leading jet
559  Double_t fTitleFontSize; // title font size
560  // members, set internally
561  TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
562  TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
563  TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
564  TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
565  TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
566  TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
567  TH1D* fSpectrumIn; // in plane jet pt spectrum
568  TH1D* fSpectrumOut; // out of plane jet pt spectrum
569  TH1D* fDptInDist; // in plane dpt distribution
570  TH1D* fDptOutDist; // out of plane dpt distribution
571  TH2D* fDptIn; // in plane dpt matrix
572  TH2D* fDptOut; // out plane dpt matrix
573  TH2D* fFullResponseIn; // full response matrix, in plane
574  TH2D* fFullResponseOut; // full response matrix, out of plane
575  Float_t fPivot; // pivot of step-function
576  Bool_t fSubdueError; // no error > pivot
577  TH1* fUnfoldedSpectrumIn; // unfolded spectrum in plane
578  TH1* fUnfoldedSpectrumOut; // unfolded spectrum out of plane
579  Int_t fHarmonic; // vn harmonic
580 
581  static TArrayD* gV2; // internal use only, do not touch these
582  static TArrayD* gStat; // internal use only, do not touch these
583  static TArrayD* gShape; // internal use only, do not touch these
584  static TArrayD* gCorr; // internal use only, do not touch these
585 
586  static Int_t gOffsetStart; // see initialization below
587  static Int_t gOffsetStop; // see initialization below
588  static Float_t gReductionFactor; // multiply shape uncertainty by this factor
589  static Float_t gReductionFactorCorr; // multiply corr uncertainty by this factor
590  static Float_t gPwrtTo; // p-value will be evaluated wrt y = gPwrtTo
591  static TArrayD* gPwrtToArray; // p-value will be evaluated wrt y = gPwrtToArray[i]
592  static TArrayD* gPwrtToStatArray; // stat uncertainty of previous array
593 
594  // copy and assignment
595  AliJetFlowTools(const AliJetFlowTools&); // not implemented
596  AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
597 
598 };
599 // initialize the static members
600 TArrayD* AliJetFlowTools::gV2 = new TArrayD(7); // DO NOT TOUCH - these arrays are filled by the
601 TArrayD* AliJetFlowTools::gStat = new TArrayD(7); // 'GetSignificance' function and
602 TArrayD* AliJetFlowTools::gShape = new TArrayD(7); // then used in the chi2 minimization routine
603 TArrayD* AliJetFlowTools::gCorr = new TArrayD(7); // to calculate the significance of the results
604 Int_t AliJetFlowTools::gOffsetStart = 0; // start chi2 fit from this bin w.r.t. the binning supplied in the 'GetCorr/GetShape' functions
605 Int_t AliJetFlowTools::gOffsetStop = 0; // stop chi2 fit at this bin w.r.t. the binning supplied in the 'GetCorr/GetShape' functions
606 Float_t AliJetFlowTools::gReductionFactor = 1.; // multiply shape uncertainty by this factor
607 Float_t AliJetFlowTools::gReductionFactorCorr = 1.; // multiply corr uncertainty by this factor
608 Float_t AliJetFlowTools::gPwrtTo = 0.; // p-value will be evaluated wrt y = gPwrtTo
609 TArrayD* AliJetFlowTools::gPwrtToArray = new TArrayD(7);// array of values w.r.t. p value is evaluated
610 TArrayD* AliJetFlowTools::gPwrtToStatArray = new TArrayD(7);// array of stat errors on values w.r.t. p values is evaluated
611 #endif
612 //_____________________________________________________________________________
TH1D * UnfoldSpectrumSVD(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
void SetNormalizeSpectra(Bool_t b)
return jsonbuilder str().c_str()
void SetOutputFileName(TString name)
void SetPivot(Float_t p)
void SetBayesianSmooth(Float_t s)
void SetBinsTrue(TArrayD *bins)
void SetPrior(prior p, TF1 *function, TArrayD *bins)
static TH1D * SmoothenPrior(TH1D *spectrum, TF1 *function, Double_t min, Double_t max, Double_t start, Bool_t kill=kTRUE, Bool_t counts=kTRUE)
const Double_t ymax
Definition: AddTaskCFDStar.C:7
static void WriteObject(TObject *object, TString suffix="", Bool_t kill=kTRUE)
void SetCentralityBin(TArrayI *bins)
double Double_t
Definition: External.C:58
static TH2D * ResizeYaxisTH2D(TH2D *histo, TArrayD *x, TArrayD *y, TString suffix="")
TH1D * UnfoldSpectrumBayesian(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
void SetMergeWith(TList *l, Int_t c, Float_t weight)
Bool_t fUseDetectorResponse
void SetBeta(Double_t i, Double_t o)
static void SetMinuitPrecision(Int_t i)
static TH2D * NormalizeTH2D(TH2D *histo, Bool_t noError=kTRUE)
static Float_t gPwrtTo
void SetSVDReg(Int_t r)
void SetJetFindingEfficiency(TH1D *e)
TArrayD * fBinsRecPrior
void SetReductionFactorCorr(Float_t g)
static Int_t gOffsetStart
static TLatex * AddLogo(Int_t logo, Double_t xmin=.59, Double_t ymax=.81)
static TLatex * AddTLatex(Double_t xmin, Double_t ymax, TString string)
void SetWeightFunction(TF1 *w)
TDirectoryFile * fActiveDir
void SetBinsRecPrior(TArrayD *bins)
static Double_t ConstructFunctionnd(Double_t *x, Double_t *par)
TList * list
TH1D * UnfoldSpectrumBayesianAli(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
static TGraphAsymmErrors * AddHistoErrorsToGraphErrors(TGraphAsymmErrors *g, TH1D *h)
static Double_t PhenixChi2nd(const Double_t *xx)
void SetBinsTruePrior(TArrayD *bins)
AliAnaChargedJetResponseMaker * fResponseMaker
void GetNominalValues(TH1D *&ratio, TGraphErrors *&v2, TArrayI *in, TArrayI *out, TString inFile="UnfoldedSpectra.root", TString outFile="Nominal.root") const
TCanvas * c
Definition: TestFitELoss.C:172
Float_t fBayesianSmoothOut
static TArrayD * gStat
void SetExLJDpt(Bool_t i)
TProfile * fRMSSpectrumIn
static TH1D * MergeSpectrumBins(TArrayI *bins, TH1D *spectrum, TH2D *corr)
TArrayD * fBinsTruePrior
static TH1D * GetV2Histo(TH1 *h1=0x0, TH1 *h2=0x0, Double_t r=0., TString name="")
void SetDphiDptUnfolding(Bool_t i)
static Float_t gReductionFactor
static Float_t gReductionFactorCorr
void BootstrapSpectra(TString def, TString in="UnfoldedSpectra.root", TString out="BootstrapSpectra.root") const
void ReplaceBins(TArrayI *array, TGraphAsymmErrors *graph)
void DoIntermediateSystematicsOnV2(TArrayI *variationsIn, TArrayI *variationsOut, TH1D *&relativeErrorInUp, TH1D *&relativeErrorInLow, TH1D *&relativeErrorOutUp, TH1D *&relativeErrorOutLow, TH1D *&relativeSystematicIn, TH1D *&relativeSystematicOut, TH1D *&nominal, TH1D *&nominalIn, TH1D *&nominalOut, Int_t columns, Float_t rangeLow, Float_t rangeUp, TFile *readMe, TString source="", Bool_t RMS=kFALSE) const
TRandom * gRandom
static TH2D * MatrixMultiplication(TH2D *a, TH2D *b, TString name="CombinedResponse")
static TPaveText * AddTPaveText(TString text, Int_t r=2, Double_t a=.587, Double_t b=.695, Double_t c=.872, Double_t d=.801)
void SetDphiUnfolding(Bool_t i)
static TH2D * ConstructDPtResponseFromTH1D(TH1D *dpt, Bool_t AvoidRoundingError)
Float_t fMergeWithWeight
void SetAvoidRoundingError(Bool_t r)
void SetSaveFull(Bool_t b)
Bool_t PrepareForUnfolding(TH1 *customIn=0x0, TH1 *customOut=0x0)
void PostProcess(TString def, Int_t columns=4, Float_t rangeLow=20, Float_t rangeUp=80, TString in="UnfoldedSpectra.root", TString out="ProcessedSpectra.root") const
static Int_t gOffsetStop
void SetOffsetStart(Int_t g)
static TArrayD * gCorr
static void GetSignificance(TGraphErrors *n, TGraphAsymmErrors *shape, TGraphAsymmErrors *corr, Int_t low, Int_t up)
void SetBeta(Double_t b)
int Int_t
Definition: External.C:63
static void Style(Bool_t legacy=kFALSE)
Double_t fTitleFontSize
static void SetMinuitStrategy(Double_t s)
float Float_t
Definition: External.C:68
void Make(TH1 *customIn=0x0, TH1 *customOut=0x0)
void SetMergeSpectrumBins(TArrayI *a)
static void SavePadToPDF(TVirtualPad *pad)
void GetShapeUncertainty(TGraphAsymmErrors *&shapeRatio, TGraphAsymmErrors *&shapeV2, TArrayI *regularizationIn, TArrayI *regularizationOut, TArrayI *recBinIn=0x0, TArrayI *recBinOut=0x0, TArrayI *methodIn=0x0, TArrayI *methodOut=0x0, Int_t columns=4, Float_t rangeLow=20, Float_t rangeUp=80, Float_t corr=.0, TString in="UnfoldedSpectra.root", TString out="ShapeUncertainty.root", Bool_t regularizationOnV2=kTRUE) const
static void SetMinuitStepSize(Float_t s)
void SetBinsRec(TArrayD *bins)
static TH1D * RebinTH1D(TH1D *histo, TArrayD *bins, TString suffix="", Bool_t kill=kTRUE)
void SetPwrtTo(Float_t p)
void SetPwrtToStatArray(TArrayD *a)
static TH1F * ConvertGraphToHistogram(TGraphErrors *g)
void SetUnfoldingAlgorithm(unfoldingAlgorithm ua)
Definition: External.C:228
Definition: External.C:212
void SetTrainPowerFit(Bool_t t)
void SetPwrtToArray(TArrayD *a)
void SetSubdueError(Bool_t b)
void SetUseDptResponse(Bool_t r)
void GetCorrelatedUncertainty(TGraphAsymmErrors *&corrRatio, TGraphAsymmErrors *&corrV2, TArrayI *variationsIn, TArrayI *variationsOut, Bool_t sym, TArrayI *variantions2ndIn, TArrayI *variantions2ndOut, Bool_t sym2nd, TString type="", TString type2="", Int_t columns=4, Float_t rangeLow=20, Float_t rangeUp=80, Float_t corr=.5, TString in="UnfoldedSpectra.root", TString out="CorrelatedUncertainty.root") const
unfoldingAlgorithm fUnfoldingAlgorithm
void SetConstantUE(Bool_t ue)
static TArrayD * gShape
void SetDetectorResponse(TH2D *dr)
static TGraphErrors * GetRatio(TH1 *h1=0x0, TH1 *h2=0x0, TString name="", Bool_t appendFit=kFALSE, Int_t xmax=-1)
Bool_t fAvoidRoundingError
void SetBayesianIter(Int_t i, Int_t o)
TH1D * GetPrior(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency)
static TF1 * GetErrorFromFit(TH1 *h1, TH1 *h2, Double_t a, Double_t b, Float_t pivot=50., Bool_t subdueError=kFALSE, TString str="", Bool_t setContent=kTRUE)
void DoIntermediateSystematics(TArrayI *variationsIn, TArrayI *variationsOut, TH1D *&relativeErrorInUp, TH1D *&relativeErrorInLow, TH1D *&relativeErrorOutUp, TH1D *&relativeErrorOutLow, TH1D *&relativeSystematicIn, TH1D *&relativeSystematicOut, TH1D *&nominal, TH1D *&nominalIn, TH1D *&nominalOut, Int_t columns, Float_t rangeLow, Float_t rangeUp, TFile *readMe, TString source="", Bool_t RMS=kFALSE) const
static TLegend * AddLegend(TVirtualPad *p, Bool_t style=kFALSE)
TH1D * UnfoldSpectrumChi2(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
static TH1D * ResizeXaxisTH1D(TH1D *histo, Int_t low, Int_t up, TString suffix="")
TH1D * UnfoldWrapper(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
static TGraphErrors * GetV2(TH1 *h1=0x0, TH1 *h2=0x0, Double_t r=0., TString name="")
Float_t fBayesianSmoothIn
TH1D * ProtectHeap(TH1D *protect, Bool_t kill=kTRUE, TString suffix="") const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void SetSymmRMS(Bool_t r)
static void SetDebug(Int_t d)
void SetCentralityWeight(TArrayD *weights)
void SetInputList(TList *list)
TGraphAsymmErrors * GetV2WithSystematicErrors(TH1 *h1, TH1 *h2, Double_t r, TString name, TH1 *relativeErrorInUp, TH1 *relativeErrorInLow, TH1 *relativeErrorOutUp, TH1 *relativeErrorOutLow, Float_t rho=0.) const
void CreateOutputList(TString name)
static TArrayD * gPwrtToStatArray
void SetBayesianSmooth(Float_t i, Float_t o)
TArrayI * fMergeBinsArray
void SetSVDToy(Bool_t b, Float_t r)
void SetRho0(Bool_t r)
static TH2D * GetUnityResponse(TArrayD *binsTrue, TArrayD *binsRec, TString suffix="")
static void RemoveSign(Double_t &d)
void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const
TH2D * RebinTH2D(TH2D *histo, TArrayD *binsTrue, TArrayD *binsRec, TString suffix="")
static TMatrixD * CalculatePearsonCoefficients(TMatrixD *covmat)
static TLatex * AddSystem()
static TH1 * Bootstrap(TH1 *hist, Bool_t kill=kTRUE)
TArrayI * fCentralityArray
void SetRMS(Bool_t r)
void SetTitleFontSize(Double_t s)
static TArrayD * gV2
static TF2 * ReturnFunctionnd(Double_t &p)
void SetTestMode(Bool_t t)
void SetHarmonic(Int_t n)
void SetBootstrap(Bool_t b, Bool_t r=kTRUE)
void SetNormalizeSpectra(Int_t e)
AliJetFlowTools & operator=(const AliJetFlowTools &)
Double_t fEventPlaneRes
void SetUseDetectorResponse(Bool_t r)
bool Bool_t
Definition: External.C:53
TProfile * fRMSSpectrumOut
void SetPrior(prior p, TH1D *spectrum)
Bool_t SetRawInput(TH2D *detectorResponse, TH1D *jetPtIn, TH1D *jetPtOut, TH1D *dptIn, TH1D *dptOut, Int_t eventCount=0)
TH1D * FoldSpectrum(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
TH1 * GetUnfoldedSpectrumIn() const
static void ResetAliUnfolding()
TH1 * GetUnfoldedSpectrumOut() const
void SystematicsWrapper(TArrayI *variationsIn, TArrayI *variationsOut, TH1D *&relativeErrorInUp, TH1D *&relativeErrorInLow, TH1D *&relativeErrorOutUp, TH1D *&relativeErrorOutLow, TH1D *&relativeSystematicIn, TH1D *&relativeSystematicOut, TH1D *&nominal, TH1D *&nominalIn, TH1D *&nominalOut, Int_t columns, Float_t rangeLow, Float_t rangeUp, TFile *readMe, TString source="", Bool_t RMS=kFALSE, Bool_t onRatio=kTRUE) const
void SetReductionFactor(Float_t g)
void SetSmoothenPrior(Bool_t b, Float_t min=50., Float_t max=100., Float_t start=75., Bool_t counts=kTRUE)
static TPaveText * AddText(TString text, EColor col, Double_t a=.2098, Double_t b=.5601, Double_t c=.613, Double_t d=.6211)
TProfile * fRMSRatio
static void SquelchWarning()
static void MinimizeChi2nd()
void SetBayesianIter(Int_t i)
Definition: External.C:196
void SetOffsetStop(Int_t g)
static TH1D * NormalizeTH1D(TH1D *histo, Double_t scale=1.)
TArrayD * fCentralityWeights
static void SetMinuitPrecision(Float_t s)
void SetSVDReg(Int_t in, Int_t out)
void SetPrior(prior p)
void SetEventPlaneResolution(Double_t r)
void SetCentralityBin(Int_t bin)
static TArrayD * gPwrtToArray
static Double_t GetRMSOfTH1(TH1 *h, Double_t a, Double_t b)