1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
5 /* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
6 * see implementation for additional information */
8 #ifndef AliJetFlowTools_H
9 #define AliJetFlowTools_H
11 // root forward declarations
20 // aliroot forward declarations
21 class AliAnaChargedJetResponseMaker;
27 #include "TDirectoryFile.h"
30 #include "TVirtualPad.h"
31 #include "TPaveText.h"
36 //_____________________________________________________________________________
37 class AliJetFlowTools {
41 ~AliJetFlowTools(); // not implemented (deliberately). object ownership is a bit messy in this class
42 // since most (or all) of the objects are owned by the input and output files
45 enum unfoldingAlgorithm { // type of unfolding alrogithm
46 kChi2, // chi^2 unfolding, implemented in AliUnfolding
47 kBayesian, // Bayesian unfolding, implemented in RooUnfold
48 kBayesianAli, // Bayesian unfolding, implemented in AliUnfolding
49 kSVD, // SVD unfolding, implemented in RooUnfold
50 kNone }; // no unfolding
51 enum prior { // prior that is used for unfolding
52 kPriorChi2, // prior from chi^2 method
53 kPriorMeasured, // use measured spectrum as prior
54 kPriorPythia, // use pythia spectrum as prior
55 kPriorTF1 }; // use properly binned TF1 as prior
56 enum histoType { // histogram identifier, only used internally
57 kInPlaneSpectrum, // default style for spectrum
62 kBar, // default style for bar histogram
63 kRatio, // default style for ratio
64 kV2, // default style for v2
65 kDeltaPhi, // default for delta phi
66 kEmpty }; // default style
67 // setters, interface to the class
68 void SetSaveFull(Bool_t b) {fSaveFull = b;}
69 void SetInputList(TList* list) {
71 fRefreshInput = kTRUE;
73 void SetOutputFileName(TString name) {fOutputFileName = name;}
74 void CreateOutputList(TString name) {
75 // create a new output list and add it to the full output
76 if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
77 fOutputFile->cd(); // avoid nested dirs
78 if(name.EqualTo(fActiveString)) {
79 printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
83 fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
86 void SetCentralityBin(Int_t bin) {
87 // in case of one centraltiy
88 fCentralityArray = new TArrayI(1);
89 fCentralityArray->AddAt(bin, 0);
90 // for one centrality there's no need for weights
91 fCentralityWeights = new TArrayD(1);
92 fCentralityWeights->AddAt(1., 0);
94 void SetMergeSpectrumBins(TArrayI* a) {fMergeBinsArray = a;}
95 void SetCentralityBin(TArrayI* bins) {
96 fCentralityArray = bins;
98 void SetCentralityWeight(TArrayD* weights) {
99 fCentralityWeights = weights;
100 if(!fCentralityArray) printf(" > Warning: centrality weights set, but bins are not defined! \n");
102 void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;}
103 void SetJetFindingEfficiency(TH1D* e) {fJetFindingEff = e;}
104 void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
105 void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
106 void SetBinsTruePrior(TArrayD* bins) {fBinsTruePrior = bins;}
107 void SetBinsRecPrior(TArrayD* bins) {fBinsRecPrior = bins;}
108 void SetSVDReg(Int_t r) {fSVDRegIn = r; fSVDRegOut = r;}
109 void SetSVDReg(Int_t in, Int_t out) {fSVDRegIn = in; fSVDRegOut = out;}
110 void SetSVDToy(Bool_t b, Float_t r) {fSVDToy = b; fJetRadius = r;}
111 void SetBeta(Double_t b) {fBetaIn = b; fBetaOut = b;}
112 void SetBeta(Double_t i, Double_t o) {fBetaIn = i; fBetaOut = o;}
113 void SetBayesianIter(Int_t i) {fBayesianIterIn = i; fBayesianIterOut = i;}
114 void SetBayesianIter(Int_t i, Int_t o) {fBayesianIterIn = i; fBayesianIterOut = o;}
115 void SetBayesianSmooth(Float_t s) {fBayesianSmoothIn = s; fBayesianSmoothOut = s;}
116 void SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn = i; fBayesianSmoothOut = o;}
117 void SetAvoidRoundingError(Bool_t r) {fAvoidRoundingError = r;}
118 void SetUnfoldingAlgorithm(unfoldingAlgorithm ua) {fUnfoldingAlgorithm = ua;}
119 void SetPrior(prior p) {fPrior = p;}
120 void SetPrior(prior p, TF1* function, TArrayD* bins) {
122 // set prior to value supplied in TF1
123 fPriorUser = new TH1D("prior_user", "prior_user", bins->GetSize()-1, bins->GetArray());
124 // loop over bins and fill the histo from the tf1
125 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));
127 void SetPrior(prior p, TH1D* spectrum) {fPrior = p; fPriorUser = spectrum;}
128 void SetNormalizeSpectra(Bool_t b) {fNormalizeSpectra = b;}
129 void SetNormalizeSpectra(Int_t e) { // use to normalize to this no of events
131 fNormalizeSpectra = kFALSE;
133 void SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
138 fSmoothenCounts = counts;
140 void SetTestMode(Bool_t t) {fTestMode = t;}
141 void SetEventPlaneResolution(Double_t r) {fEventPlaneRes = r;}
142 void SetUseDetectorResponse(Bool_t r) {fUseDetectorResponse = r;}
143 void SetUseDptResponse(Bool_t r) {fUseDptResponse = r;}
144 void SetTrainPowerFit(Bool_t t) {fTrainPower = t;}
145 void SetDphiUnfolding(Bool_t i) {fDphiUnfolding = i;}
146 void SetDphiDptUnfolding(Bool_t i) {fDphiDptUnfolding = i;}
147 void SetExLJDpt(Bool_t i) {fExLJDpt = i;}
148 void SetWeightFunction(TF1* w) {fResponseMaker->SetRMMergeWeightFunction(w);}
149 void SetRMS(Bool_t r) {fRMS = r;}
150 void SetSymmRMS(Bool_t r) {fSymmRMS = r; fRMS = r;}
151 void SetRho0(Bool_t r) {fRho0 = r;}
152 void SetBootstrap(Bool_t b, Bool_t r = kTRUE) {
153 // note that setting this option will not lead to true resampling
154 // but rather to randomly drawing a new distribution from a pdf
155 // of the measured distribution
157 // by default fully randomize randomizer from system time
160 gRandom = new TRandom3(0);
164 void MakeAU(); // test function, use with caution (09012014)
167 if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
168 if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
169 if(fRMSRatio) fRMSRatio->Write();
170 fOutputFile->Close();}
174 Float_t rangeLow = 20,
175 Float_t rangeUp = 80,
176 TString in = "UnfoldedSpectra.root",
177 TString out = "ProcessedSpectra.root") const;
178 void BootstrapSpectra(
180 TString in = "UnfoldedSpectra.root",
181 TString out = "BootstrapSpectra.root") const;
182 void GetNominalValues(
187 TString inFile = "UnfoldedSpectra.root",
188 TString outFile = "Nominal.root") const;
189 void GetCorrelatedUncertainty(
190 TGraphAsymmErrors*& corrRatio,
191 TGraphAsymmErrors*& corrV2,
192 TArrayI* variationsIn,
193 TArrayI* variationsOut,
195 TArrayI* variantions2ndIn,
196 TArrayI* variantions2ndOut,
201 Float_t rangeLow = 20,
202 Float_t rangeUp = 80,
204 TString in = "UnfoldedSpectra.root",
205 TString out = "CorrelatedUncertainty.root") const;
206 void GetShapeUncertainty(
207 TGraphAsymmErrors*& shapeRatio,
208 TGraphAsymmErrors*& shapeV2,
209 TArrayI* regularizationIn,
210 TArrayI* regularizationOut,
211 TArrayI* trueBinIn = 0x0,
212 TArrayI* trueBinOut = 0x0,
213 TArrayI* recBinIn = 0x0,
214 TArrayI* recBinOut = 0x0,
215 TArrayI* methodIn = 0x0,
216 TArrayI* methodOut = 0x0,
218 Float_t rangeLow = 20,
219 Float_t rangeUp = 80,
220 TString in = "UnfoldedSpectra.root",
221 TString out = "ShapeUncertainty.root") const;
223 TH2D* detectorResponse, // detector response matrix
224 TH1D* jetPtIn, // in plane jet spectrum
225 TH1D* jetPtOut, // out of plane jet spectrum
226 TH1D* dptIn, // in plane delta pt distribution
227 TH1D* dptOut, // out of plane delta pt distribution
228 Int_t eventCount = 0); // event count (optional)
229 // static const helper functions, mainly histogram manipulation
230 static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
231 static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
232 static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
233 static TH1* Bootstrap(TH1* hist, Bool_t kill = kTRUE);
234 static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
235 TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
236 static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
237 static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
238 static TH1D* MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr);
239 static TGraphErrors* GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
240 static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
241 void ReplaceBins(TArrayI* array, TGraphAsymmErrors* graph);
242 void ReplaceBins(TArrayI* array, TGraphErrors* graph);
243 TGraphAsymmErrors* GetV2WithSystematicErrors(
244 TH1* h1, TH1* h2, Double_t r, TString name,
245 TH1* relativeErrorInUp,
246 TH1* relativeErrorInLow,
247 TH1* relativeErrorOutUp,
248 TH1* relativeErrorOutLow,
249 Float_t rho = 0.) const;
250 static void GetSignificance(
251 TGraphErrors* n, // points with stat error
252 TGraphAsymmErrors* shape, // points with shape error
253 TGraphAsymmErrors* corr, // corr with stat error
254 Int_t low, // pt lower level
255 Int_t up // pt upper level
257 static void MinimizeChi22d();
258 static Double_t PhenixChi22d(const Double_t *xx );
259 static Double_t ConstructFunction2d(Double_t *x, Double_t *par);
260 static TF2* ReturnFunction2d();
261 static void MinimizeChi2nd();
262 static Double_t PhenixChi2nd(const Double_t *xx );
263 static Double_t ConstructFunctionnd(Double_t *x, Double_t *par);
264 static void MinimizeChi2();
265 static TF2* ReturnFunctionnd();
266 static Double_t PhenixChi2(const Double_t *xx );
267 static Double_t ConstructFunction(Double_t *x, Double_t *par);
268 static TF1* ReturnFunction();
269 static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
270 static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
271 static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
272 void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
273 static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
274 static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
276 void SetTitleFontSize(Double_t s) {fTitleFontSize = s;}
277 static void Style(Bool_t legacy = kFALSE);
278 static void Style(TCanvas* c, TString style = "PEARSON");
279 static void Style(TVirtualPad* c, TString style = "SPECTRUM", Bool_t legacy = kFALSE);
280 static void Style(TLegend* l);
281 static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
282 static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty, Bool_t legacy = kFALSE);
283 static TLegend* AddLegend(TVirtualPad* p, Bool_t style = kFALSE) {
284 if(!style) return p->BuildLegend(.565, .663, .882, .883);
286 TLegend* l = AddLegend(p, kFALSE);
291 static TPaveText* AddTPaveText(
292 // this text appears under the logo
299 TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
302 t->AddText(0.,0.,text.Data());
303 t->AddText(0., 0., Form("#it{R} = 0.%i anti-#it{k}_{T}, |#eta_{jet}|<%.1f", r, .9-r/10.));
304 t->SetTextColor(kBlack);
306 t->SetTextSize(gStyle->GetTextSize()*.8);
310 static TPaveText* AddText(
316 Double_t d = .6211) {
317 TPaveText* t(new TPaveText(a, b, c, d, "NDC"));
320 t->AddText(0.,0.,text.Data());
321 t->SetTextColor(col);
323 t->SetTextSize(gStyle->GetTextSize()*.8);
327 static TLatex* AddLogo(Int_t logo, Double_t xmin = .59, Double_t ymax = .81) {
328 if(logo == 0) return AddTLatex(xmin, ymax, "ALICE");
329 else if (logo == 1) return AddTLatex(xmin, ymax, "ALICE Preliminary");
330 else if (logo == 2) return AddTLatex(xmin, ymax, "ALICE Simulation");
333 static TLatex* AddSystem() {
334 return AddTLatex(0.55, 87, "Pb-Pb #sqrt{#it{s}}}_{NN} = 2.76 TeV");
336 static TLatex* AddTLatex(Double_t xmin, Double_t ymax, TString string) {
337 TLatex* tex = new TLatex(xmin, ymax, string.Data());
339 tex->SetTextFont(42);
344 static void SavePadToPDF(TVirtualPad* pad) {return;/*pad->SaveAs(Form("%s.pdf", pad->GetName()));*/}
345 // interface to AliUnfolding, not necessary but nice to have all parameters in one place
346 static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
347 static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
348 static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
349 static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
350 static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
352 Bool_t PrepareForUnfolding();
353 Bool_t PrepareForUnfolding(Int_t low, Int_t up);
354 TH1D* GetPrior( const TH1D* measuredJetSpectrum,
355 const TH2D* resizedResponse,
356 const TH1D* kinematicEfficiency,
357 const TH1D* measuredJetSpectrumTrueBins,
358 const TString suffix,
359 const TH1D* jetFindingEfficiency);
360 TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
361 const TH2D* resizedResponse,
362 const TH1D* kinematicEfficiency,
363 const TH1D* measuredJetSpectrumTrueBins,
364 const TString suffix,
365 const TH1D* jetFindingEfficiency = 0x0);
366 TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
367 const TH2D* resizedResponse,
368 const TH1D* kinematicEfficiency,
369 const TH1D* measuredJetSpectrumTrueBins,
370 const TString suffix,
371 const TH1D* jetFindingEfficiency = 0x0);
372 TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
373 const TH2D* resizedResponse,
374 const TH1D* kinematicEfficiency,
375 const TH1D* measuredJetSpectrumTrueBins,
376 const TString suffix,
377 const TH1D* jetFindingEfficiency = 0x0);
378 TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
379 const TH2D* resizedResponse,
380 const TH1D* kinematicEfficiency,
381 const TH1D* measuredJetSpectrumTrueBins,
382 const TString suffix,
383 const TH1D* jetFindingEfficiency = 0x0);
384 TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
385 const TH2D* resizedResponse,
386 const TH1D* kinematicEfficiency,
387 const TH1D* measuredJetSpectrumTrueBins,
388 const TString suffix,
389 const TH1D* jetFindingEfficiency = 0x0);
390 void DoIntermediateSystematics(
391 TArrayI* variationsIn,
392 TArrayI* variationsOut,
393 TH1D*& relativeErrorInUp,
394 TH1D*& relativeErrorInLow,
395 TH1D*& relativeErrorOutUp,
396 TH1D*& relativeErrorOutLow,
397 TH1D*& relativeSystematicIn,
398 TH1D*& relativeSystematicOut,
407 Bool_t RMS = kFALSE) const;
408 static void ResetAliUnfolding();
409 // give object a unique name via the 'protect heap' functions.
410 // may seem redundant, but some internal functions of root (e.g.
411 // ProjectionY()) check for existing objects by name and re-use them
412 TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
413 TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
414 TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
415 // members, accessible via setters
416 AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
417 Bool_t fRMS; // systematic method
418 Bool_t fSymmRMS; // symmetric systematic method
419 Bool_t fRho0; // use the result obtained with the 'classic' fixed rho
420 Bool_t fBootstrap; // use bootstrap resampling of input data
421 TF1* fPower; // smoothening fit
422 Bool_t fSaveFull; // save all generated histograms to file
423 TString fActiveString; // identifier of active output
424 TDirectoryFile* fActiveDir; // active directory
425 TList* fInputList; // input list
426 Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
427 TString fOutputFileName; // output file name
428 TFile* fOutputFile; // output file
429 TArrayI* fCentralityArray; // array of centrality bins that are merged
430 TArrayI* fMergeBinsArray; // array of pt bins that are merged
431 TArrayD* fCentralityWeights; // array of centrality weights
432 TH2D* fDetectorResponse; // detector response
433 TH1D* fJetFindingEff; // jet finding efficiency
434 Double_t fBetaIn; // regularization strength, in plane unfolding
435 Double_t fBetaOut; // regularization strength, out of plane unfoldign
436 Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
437 Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
438 Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
439 Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
440 Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
441 unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
442 prior fPrior; // prior for unfolding
443 TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
444 TArrayD* fBinsTrue; // pt true bins
445 TArrayD* fBinsRec; // pt rec bins
446 TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
447 TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
448 Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
449 Int_t fSVDRegOut; // svd regularization out of plane
450 Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
451 Float_t fJetRadius; // jet radius (for SVD toy)
452 Int_t fEventCount; // number of events
453 Bool_t fNormalizeSpectra; // normalize spectra to event count
454 Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
455 Float_t fFitMin; // lower bound of smoothening fit
456 Float_t fFitMax; // upper bound of smoothening fit
457 Float_t fFitStart; // from this value, use smoothening
458 Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
459 Bool_t fTestMode; // unfold with unity response for testing
460 Bool_t fRawInputProvided; // input histograms provided, not read from file
461 Double_t fEventPlaneRes; // event plane resolution for current centrality
462 Bool_t fUseDetectorResponse; // add detector response to unfolding
463 Bool_t fUseDptResponse; // add dpt response to unfolding
464 Bool_t fTrainPower; // don't clear the params of fPower for call to Make
465 // might give more stable results, but possibly introduces
466 // a bias / dependency on previous iterations
467 Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
468 Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
469 Bool_t fExLJDpt; // exclude randon cones with leading jet
470 Double_t fTitleFontSize; // title font size
471 // members, set internally
472 TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
473 TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
474 TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
475 TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
476 TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
477 TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
478 TH1D* fSpectrumIn; // in plane jet pt spectrum
479 TH1D* fSpectrumOut; // out of plane jet pt spectrum
480 TH1D* fDptInDist; // in plane dpt distribution
481 TH1D* fDptOutDist; // out of plane dpt distribution
482 TH2D* fDptIn; // in plane dpt matrix
483 TH2D* fDptOut; // out plane dpt matrix
484 TH2D* fFullResponseIn; // full response matrix, in plane
485 TH2D* fFullResponseOut; // full response matrix, out of plane
487 // copy and assignment
488 AliJetFlowTools(const AliJetFlowTools&); // not implemented
489 AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
492 //_____________________________________________________________________________