AliPhysics  e59a9ba (e59a9ba)
 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
74  void SetOffsetStart(Int_t g) {gOffsetStart = g;}
75  void SetOffsetStop(Int_t g) {gOffsetStop = g;}
76  void SetReductionFactor(Float_t g) {gReductionFactor = g;}
78  void SetPwrtTo(Float_t p) {gPwrtTo = p;}
79  void SetPwrtToArray(TArrayD* a) {
80  gPwrtToArray = a; // set the array
81  gPwrtTo = -9999; // tells instance to not use this value
82  }
83  void SetPwrtToStatArray(TArrayD* a) {gPwrtToStatArray = a;}
84  void SetPivot(Float_t p) {fPivot = p;}
85  void SetConstantUE(Bool_t ue) {fConstantUE = ue;}
86  void SetSubdueError(Bool_t b) {fSubdueError = b;}
87  void SetSaveFull(Bool_t b) {fSaveFull = b;}
88  void SetInputList(TList* list) {
89  fInputList = list;
90  fRefreshInput = kTRUE;
91  }
92  void SetOutputFileName(TString name) {fOutputFileName = name;}
93  void CreateOutputList(TString name) {
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  }
106  void SetCentralityBin(Int_t bin) {
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  }
114  void SetMergeSpectrumBins(TArrayI* a) {fMergeBinsArray = a;}
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  }
130  void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;}
132  void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
133  void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
134  void SetBinsTruePrior(TArrayD* bins) {fBinsTruePrior = 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;}
142  void SetBayesianIter(Int_t i, Int_t o) {fBayesianIterIn = i; fBayesianIterOut = o;}
144  void SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn = i; fBayesianSmoothOut = 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;}
169  void SetEventPlaneResolution(Double_t r) {fEventPlaneRes = r;}
171  void SetUseDptResponse(Bool_t r) {fUseDptResponse = r;}
172  void SetTrainPowerFit(Bool_t t) {fTrainPower = t;}
173  void SetDphiUnfolding(Bool_t i) {fDphiUnfolding = i;}
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  // main function. buffers about 5mb per call!
192  void Make(TH1* customIn = 0x0, TH1* customOut = 0x0);
193  void MakeAU(); // test function, use with caution (09012014)
194  void Finish() {
195  fOutputFile->cd();
196  if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
197  if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
198  if(fRMSRatio) fRMSRatio->Write();
199  fOutputFile->Close();}
200  void PostProcess(
201  TString def,
202  Int_t columns = 4,
203  Float_t rangeLow = 20,
204  Float_t rangeUp = 80,
205  TString in = "UnfoldedSpectra.root",
206  TString out = "ProcessedSpectra.root") const;
207  void BootstrapSpectra(
208  TString def,
209  TString in = "UnfoldedSpectra.root",
210  TString out = "BootstrapSpectra.root") const;
211  void GetNominalValues(
212  TH1D*& ratio,
213  TGraphErrors*& v2,
214  TArrayI* in,
215  TArrayI* out,
216  TString inFile = "UnfoldedSpectra.root",
217  TString outFile = "Nominal.root") const;
219  TGraphAsymmErrors*& corrRatio,
220  TGraphAsymmErrors*& corrV2,
221  TArrayI* variationsIn,
222  TArrayI* variationsOut,
223  Bool_t sym,
224  TArrayI* variantions2ndIn,
225  TArrayI* variantions2ndOut,
226  Bool_t sym2nd,
227  TString type = "",
228  TString type2 = "",
229  Int_t columns = 4,
230  Float_t rangeLow = 20,
231  Float_t rangeUp = 80,
232  Float_t corr = .5,
233  TString in = "UnfoldedSpectra.root",
234  TString out = "CorrelatedUncertainty.root") const;
235  void GetShapeUncertainty(
236  TGraphAsymmErrors*& shapeRatio,
237  TGraphAsymmErrors*& shapeV2,
238  TArrayI* regularizationIn,
239  TArrayI* regularizationOut,
240  TArrayI* recBinIn = 0x0,
241  TArrayI* recBinOut = 0x0,
242  TArrayI* methodIn = 0x0,
243  TArrayI* methodOut = 0x0,
244  Int_t columns = 4,
245  Float_t rangeLow = 20,
246  Float_t rangeUp = 80,
247  Float_t corr = .0,
248  TString in = "UnfoldedSpectra.root",
249  TString out = "ShapeUncertainty.root",
250  Bool_t regularizationOnV2 = kTRUE // get uncertainty on yields separately or v2 directly
251  ) const;
252  Bool_t SetRawInput (
253  TH2D* detectorResponse, // detector response matrix
254  TH1D* jetPtIn, // in plane jet spectrum
255  TH1D* jetPtOut, // out of plane jet spectrum
256  TH1D* dptIn, // in plane delta pt distribution
257  TH1D* dptOut, // out of plane delta pt distribution
258  Int_t eventCount = 0); // event count (optional)
261  // static const helper functions, mainly histogram manipulation
262  static void RemoveSign(Double_t& d) {d = TMath::Abs(d);}
263  static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
264  static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
265  static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
266  static TH1* Bootstrap(TH1* hist, Bool_t kill = kTRUE);
267  static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
268  TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
269  static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
270  static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
271  static TH1D* MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr);
272  static TGraphErrors* GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
273  static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
274  static TH1D* GetV2Histo(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
275  static TH1F* ConvertGraphToHistogram(TGraphErrors* g);
276  static TGraphAsymmErrors* AddHistoErrorsToGraphErrors(TGraphAsymmErrors* g, TH1D* h);
277  static Double_t GetRMSOfTH1(TH1* h, Double_t a, Double_t b);
278  static TF1* GetErrorFromFit(TH1* h1, TH1* h2, Double_t a, Double_t b,
279  Float_t pivot = 50., Bool_t subdueError = kFALSE,
280  TString str = "", Bool_t setContent = kTRUE);
281  void ReplaceBins(TArrayI* array, TGraphAsymmErrors* graph);
282  void ReplaceBins(TArrayI* array, TGraphErrors* graph);
283  TGraphAsymmErrors* GetV2WithSystematicErrors(
284  TH1* h1, TH1* h2, Double_t r, TString name,
285  TH1* relativeErrorInUp,
286  TH1* relativeErrorInLow,
287  TH1* relativeErrorOutUp,
288  TH1* relativeErrorOutLow,
289  Float_t rho = 0.) const;
290  static void GetSignificance(
291  TGraphErrors* n, // points with stat error
292  TGraphAsymmErrors* shape, // points with shape error
293  TGraphAsymmErrors* corr, // corr with stat error
294  Int_t low, // pt lower level
295  Int_t up // pt upper level
296  );
297  static void MinimizeChi2nd();
298  static Double_t PhenixChi2nd(const Double_t *xx );
299  static Double_t ConstructFunctionnd(Double_t *x, Double_t *par);
300  static TF2* ReturnFunctionnd(Double_t &p);
301  static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
302  static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
303  static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
304  void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
305  static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
306  static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
307  // set style
308  void SetTitleFontSize(Double_t s) {fTitleFontSize = s;}
309  static void Style(Bool_t legacy = kFALSE);
310  static void Style(TCanvas* c, TString style = "PEARSON");
311  static void Style(TVirtualPad* c, TString style = "SPECTRUM", Bool_t legacy = kFALSE);
312  static void Style(TLegend* l);
313  static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
314  static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
315  static TLegend* AddLegend(TVirtualPad* p, Bool_t style = kFALSE) {
316  if(!style) return p->BuildLegend(.565, .663, .882, .883);
317  else {
318  TLegend* l = AddLegend(p, kFALSE);
319  Style(l);
320  return l;
321  }
322  }
323  static TPaveText* AddTPaveText(
324  // this text appears under the logo
325  TString text,
326  Int_t r = 2,
327  Double_t a = .587,
328  Double_t b = .695,
329  Double_t c = .872,
330  Double_t d = .801) {
331  TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
332  t->SetFillColor(0);
333  t->SetBorderSize(0);
334  t->AddText(0.,0.,text.Data());
335  t->AddText(0., 0., Form("#it{R} = 0.%i anti-#it{k}_{T}, |#eta_{jet}|<%.1f", r, .9-r/10.));
336  t->SetTextColor(kBlack);
337  t->SetTextFont(42);
338  t->SetTextSize(gStyle->GetTextSize()*.8);
339  t->Draw("same");
340  return t;
341  }
342  static TPaveText* AddText(
343  TString text,
344  EColor col,
345  Double_t a = .2098,
346  Double_t b = .5601,
347  Double_t c = .613,
348  Double_t d = .6211) {
349  TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
350  t->SetFillColor(0);
351  t->SetBorderSize(0);
352  t->AddText(0.,0.,text.Data());
353  t->SetTextColor(col);
354  t->SetTextFont(42);
355  t->SetTextSize(gStyle->GetTextSize()*.8);
356  t->Draw("same");
357  return t;
358  }
359  static TLatex* AddLogo(Int_t logo, Double_t xmin = .59, Double_t ymax = .81) {
360  if(logo == 0) return AddTLatex(xmin, ymax, "ALICE");
361  else if (logo == 1) return AddTLatex(xmin, ymax, "ALICE Preliminary");
362  else if (logo == 2) return AddTLatex(xmin, ymax, "ALICE Simulation");
363  else if (logo == 3) return AddTLatex(xmin, ymax, "work in progress");
364  return 0x0;
365  }
366  static TLatex* AddSystem() {
367  return AddTLatex(0.55, 87, "Pb-Pb #sqrt{#it{s}}}_{NN} = 2.76 TeV");
368  }
369  static TLatex* AddTLatex(Double_t xmin, Double_t ymax, TString string) {
370 TLatex* tex = new TLatex(xmin, ymax, string.Data());
371  tex->SetNDC();
372  tex->SetTextFont(42);
373  tex->Draw("same");
374  return tex;
375  }
376 
377  static void SavePadToPDF(TVirtualPad* pad) {
378  if(pad) return;/*pad->SaveAs(Form("%s.pdf", pad->GetName()));*/
379  else return;}
380  // interface to AliUnfolding, not necessary but nice to have all parameters in one place
381  static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
382  static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
383  static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
384  static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
385  static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
386  private:
387  Bool_t PrepareForUnfolding(TH1* customIn = 0x0, TH1* customOut = 0x0);
388  Bool_t PrepareForUnfolding(Int_t low, Int_t up);
389  TH1D* GetPrior( const TH1D* measuredJetSpectrum,
390  const TH2D* resizedResponse,
391  const TH1D* kinematicEfficiency,
392  const TH1D* measuredJetSpectrumTrueBins,
393  const TString suffix,
394  const TH1D* jetFindingEfficiency);
395  TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
396  const TH2D* resizedResponse,
397  const TH1D* kinematicEfficiency,
398  const TH1D* measuredJetSpectrumTrueBins,
399  const TString suffix,
400  const TH1D* jetFindingEfficiency = 0x0);
401  TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
402  const TH2D* resizedResponse,
403  const TH1D* kinematicEfficiency,
404  const TH1D* measuredJetSpectrumTrueBins,
405  const TString suffix,
406  const TH1D* jetFindingEfficiency = 0x0);
407  TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
408  const TH2D* resizedResponse,
409  const TH1D* kinematicEfficiency,
410  const TH1D* measuredJetSpectrumTrueBins,
411  const TString suffix,
412  const TH1D* jetFindingEfficiency = 0x0);
413  TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
414  const TH2D* resizedResponse,
415  const TH1D* kinematicEfficiency,
416  const TH1D* measuredJetSpectrumTrueBins,
417  const TString suffix,
418  const TH1D* jetFindingEfficiency = 0x0);
419  TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
420  const TH2D* resizedResponse,
421  const TH1D* kinematicEfficiency,
422  const TH1D* measuredJetSpectrumTrueBins,
423  const TString suffix,
424  const TH1D* jetFindingEfficiency = 0x0);
425  TH1D* FoldSpectrum( const TH1D* measuredJetSpectrum,
426  const TH2D* resizedResponse,
427  const TH1D* kinematicEfficiency,
428  const TH1D* measuredJetSpectrumTrueBins,
429  const TString suffix,
430  const TH1D* jetFindingEfficiency = 0x0);
431  void SystematicsWrapper(
432  TArrayI* variationsIn,
433  TArrayI* variationsOut,
434  TH1D*& relativeErrorInUp,
435  TH1D*& relativeErrorInLow,
436  TH1D*& relativeErrorOutUp,
437  TH1D*& relativeErrorOutLow,
438  TH1D*& relativeSystematicIn,
439  TH1D*& relativeSystematicOut,
440  TH1D*& nominal,
441  TH1D*& nominalIn,
442  TH1D*& nominalOut,
443  Int_t columns,
444  Float_t rangeLow,
445  Float_t rangeUp,
446  TFile* readMe,
447  TString source = "",
448  Bool_t RMS = kFALSE,
449  Bool_t onRatio = kTRUE) const;
451  TArrayI* variationsIn,
452  TArrayI* variationsOut,
453  TH1D*& relativeErrorInUp,
454  TH1D*& relativeErrorInLow,
455  TH1D*& relativeErrorOutUp,
456  TH1D*& relativeErrorOutLow,
457  TH1D*& relativeSystematicIn,
458  TH1D*& relativeSystematicOut,
459  TH1D*& nominal,
460  TH1D*& nominalIn,
461  TH1D*& nominalOut,
462  Int_t columns,
463  Float_t rangeLow,
464  Float_t rangeUp,
465  TFile* readMe,
466  TString source = "",
467  Bool_t RMS = kFALSE) const;
469  TArrayI* variationsIn,
470  TArrayI* variationsOut,
471  TH1D*& relativeErrorInUp,
472  TH1D*& relativeErrorInLow,
473  TH1D*& relativeErrorOutUp,
474  TH1D*& relativeErrorOutLow,
475  TH1D*& relativeSystematicIn,
476  TH1D*& relativeSystematicOut,
477  TH1D*& nominal,
478  TH1D*& nominalIn,
479  TH1D*& nominalOut,
480  Int_t columns,
481  Float_t rangeLow,
482  Float_t rangeUp,
483  TFile* readMe,
484  TString source = "",
485  Bool_t RMS = kFALSE) const;
486 
487  static void ResetAliUnfolding();
488  static void SquelchWarning() {
489  printf(" >> I squelched a warning, jay, I'm contributing ! << \n");
490  return;
491  }
492  // give object a unique name via the 'protect heap' functions.
493  // may seem redundant, but some internal functions of root (e.g.
494  // ProjectionY()) check for existing objects by name and re-use them
495  TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
496  TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
497  TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
498  // members, accessible via setters
499  Int_t fListPrefix; // prefix for list readability
500  AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
501  Bool_t fRMS; // systematic method
502  Bool_t fSymmRMS; // symmetric systematic method
503  Bool_t fConstantUE; // assign a constant unfolding error
504  Bool_t fRho0; // use the result obtained with the 'classic' fixed rho
505  Bool_t fBootstrap; // use bootstrap resampling of input data
506  TF1* fPower; // smoothening fit
507  Bool_t fSaveFull; // save all generated histograms to file
508  TString fActiveString; // identifier of active output
509  TDirectoryFile* fActiveDir; // active directory
510  TList* fInputList; // input list
511  Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
512  TString fOutputFileName; // output file name
513  TFile* fOutputFile; // output file
514  TArrayI* fCentralityArray; // array of centrality bins that are merged
515  TArrayI* fMergeBinsArray; // array of pt bins that are merged
516  TArrayD* fCentralityWeights; // array of centrality weights
517  TList* fMergeWithList; // additional input list
518  Int_t fMergeWithCen; // centrality bin of additional input list
519  Float_t fMergeWithWeight; // weight of additional input list
520  TH2D* fDetectorResponse; // detector response
521  TH1D* fJetFindingEff; // jet finding efficiency
522  Double_t fBetaIn; // regularization strength, in plane unfolding
523  Double_t fBetaOut; // regularization strength, out of plane unfoldign
524  Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
525  Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
526  Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
527  Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
528  Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
529  unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
530  prior fPrior; // prior for unfolding
531  TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
532  TArrayD* fBinsTrue; // pt true bins
533  TArrayD* fBinsRec; // pt rec bins
534  TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
535  TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
536  Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
537  Int_t fSVDRegOut; // svd regularization out of plane
538  Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
539  Float_t fJetRadius; // jet radius (for SVD toy)
540  Int_t fEventCount; // number of events
541  Bool_t fNormalizeSpectra; // normalize spectra to event count
542  Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
543  Float_t fFitMin; // lower bound of smoothening fit
544  Float_t fFitMax; // upper bound of smoothening fit
545  Float_t fFitStart; // from this value, use smoothening
546  Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
547  Bool_t fTestMode; // unfold with unity response for testing
548  Bool_t fRawInputProvided; // input histograms provided, not read from file
549  Double_t fEventPlaneRes; // event plane resolution for current centrality
550  Bool_t fUseDetectorResponse; // add detector response to unfolding
551  Bool_t fUseDptResponse; // add dpt response to unfolding
552  Bool_t fTrainPower; // don't clear the params of fPower for call to Make
553  // might give more stable results, but possibly introduces
554  // a bias / dependency on previous iterations
555  Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
556  Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
557  Bool_t fExLJDpt; // exclude randon cones with leading jet
558  Double_t fTitleFontSize; // title font size
559  // members, set internally
560  TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
561  TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
562  TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
563  TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
564  TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
565  TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
566  TH1D* fSpectrumIn; // in plane jet pt spectrum
567  TH1D* fSpectrumOut; // out of plane jet pt spectrum
568  TH1D* fDptInDist; // in plane dpt distribution
569  TH1D* fDptOutDist; // out of plane dpt distribution
570  TH2D* fDptIn; // in plane dpt matrix
571  TH2D* fDptOut; // out plane dpt matrix
572  TH2D* fFullResponseIn; // full response matrix, in plane
573  TH2D* fFullResponseOut; // full response matrix, out of plane
574  Float_t fPivot; // pivot of step-function
575  Bool_t fSubdueError; // no error > pivot
576  TH1* fUnfoldedSpectrumIn; // unfolded spectrum in plane
577  TH1* fUnfoldedSpectrumOut; // unfolded spectrum out of plane
578 
579  static TArrayD* gV2; // internal use only, do not touch these
580  static TArrayD* gStat; // internal use only, do not touch these
581  static TArrayD* gShape; // internal use only, do not touch these
582  static TArrayD* gCorr; // internal use only, do not touch these
583 
584  static Int_t gOffsetStart; // see initialization below
585  static Int_t gOffsetStop; // see initialization below
586  static Float_t gReductionFactor; // multiply shape uncertainty by this factor
587  static Float_t gReductionFactorCorr; // multiply corr uncertainty by this factor
588  static Float_t gPwrtTo; // p-value will be evaluated wrt y = gPwrtTo
589  static TArrayD* gPwrtToArray; // p-value will be evaluated wrt y = gPwrtToArray[i]
590  static TArrayD* gPwrtToStatArray; // stat uncertainty of previous array
591 
592  // copy and assignment
593  AliJetFlowTools(const AliJetFlowTools&); // not implemented
594  AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
595 
596 };
597 // initialize the static members
598 TArrayD* AliJetFlowTools::gV2 = new TArrayD(7); // DO NOT TOUCH - these arrays are filled by the
599 TArrayD* AliJetFlowTools::gStat = new TArrayD(7); // 'GetSignificance' function and
600 TArrayD* AliJetFlowTools::gShape = new TArrayD(7); // then used in the chi2 minimization routine
601 TArrayD* AliJetFlowTools::gCorr = new TArrayD(7); // to calculate the significance of the results
602 Int_t AliJetFlowTools::gOffsetStart = 0; // start chi2 fit from this bin w.r.t. the binning supplied in the 'GetCorr/GetShape' functions
603 Int_t AliJetFlowTools::gOffsetStop = 0; // stop chi2 fit at this bin w.r.t. the binning supplied in the 'GetCorr/GetShape' functions
604 Float_t AliJetFlowTools::gReductionFactor = 1.; // multiply shape uncertainty by this factor
605 Float_t AliJetFlowTools::gReductionFactorCorr = 1.; // multiply corr uncertainty by this factor
606 Float_t AliJetFlowTools::gPwrtTo = 0.; // p-value will be evaluated wrt y = gPwrtTo
607 TArrayD* AliJetFlowTools::gPwrtToArray = new TArrayD(7);// array of values w.r.t. p value is evaluated
608 TArrayD* AliJetFlowTools::gPwrtToStatArray = new TArrayD(7);// array of stat errors on values w.r.t. p values is evaluated
609 #endif
610 //_____________________________________________________________________________
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)
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
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)
static void Style(Bool_t legacy=kFALSE)
Double_t fTitleFontSize
static void SetMinuitStrategy(Double_t s)
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)
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 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)
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)
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)