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
21 // aliroot forward declarations
22 class AliAnaChargedJetResponseMaker;
27 #include "TDirectoryFile.h"
30 #include "TVirtualPad.h"
31 #include "TPaveText.h"
33 //_____________________________________________________________________________
34 class AliJetFlowTools {
38 ~AliJetFlowTools(); // not implemented (deliberately). object ownership is a bit messy in this class
39 // since most (or all) of the objects are owned by the input and output files
42 enum unfoldingAlgorithm { // type of unfolding alrogithm
43 kChi2, // chi^2 unfolding, implemented in AliUnfolding
44 kBayesian, // Bayesian unfolding, implemented in RooUnfold
45 kBayesianAli, // Bayesian unfolding, implemented in AliUnfolding
46 kSVD, // SVD unfolding, implemented in RooUnfold
47 kNone }; // no unfolding
48 enum prior { // prior that is used for unfolding
49 kPriorChi2, // prior from chi^2 method
50 kPriorMeasured, // use measured spectrum as prior
51 kPriorPythia }; // use pythia spectrum as prior
52 enum histoType { // histogram identifier, only used internally
53 kInPlaneSpectrum, // default style for spectrum
58 kBar, // default style for bar histogram
59 kRatio, // default style for ratio
60 kV2, // default style for v2
61 kEmpty }; // default style
62 // setters, interface to the class
63 void SetSaveFull(Bool_t b) {fSaveFull = b;}
64 void SetInputList(TList* list) {
66 fRefreshInput = kTRUE;
68 void SetOutputFileName(TString name) {fOutputFileName = name;}
69 void CreateOutputList(TString name) {
70 // create a new output list and add it to the full output
71 if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
72 fOutputFile->cd(); // avoid nested dirs
73 if(name.EqualTo(fActiveString)) {
74 printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
78 fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
81 void SetCentralityBin(Int_t bin) {
82 // in case of one centraltiy
83 fCentralityArray = new TArrayI(1);
84 fCentralityArray->AddAt(bin, 0);
85 // for one centrality there's no need for weights
86 fCentralityWeights = new TArrayD(1);
87 fCentralityWeights->AddAt(1., 0);
89 void SetCentralityBin(TArrayI* bins) {
90 fCentralityArray = bins;
92 void SetCentralityWeight(TArrayD* weights) {
93 fCentralityWeights = weights;
94 if(!fCentralityArray) printf(" > Warning: centrality weights set, but bins are not defined! \n");
96 void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;}
97 void SetJetFindingEfficiency(TH1D* e) {fJetFindingEff = e;}
98 void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
99 void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
100 void SetBinsTruePrior(TArrayD* bins) {fBinsTruePrior = bins;}
101 void SetBinsRecPrior(TArrayD* bins) {fBinsRecPrior = bins;}
102 void SetSVDReg(Int_t r) {fSVDRegIn = r; fSVDRegOut = r;}
103 void SetSVDReg(Int_t in, Int_t out) {fSVDRegIn = in; fSVDRegOut = out;}
104 void SetSVDToy(Bool_t b, Float_t r) {fSVDToy = b; fJetRadius = r;}
105 void SetBeta(Double_t b) {fBetaIn = b; fBetaOut = b;}
106 void SetBeta(Double_t i, Double_t o) {fBetaIn = i; fBetaOut = o;}
107 void SetBayesianIter(Int_t i) {fBayesianIterIn = i; fBayesianIterOut = i;}
108 void SetBayesianIter(Int_t i, Int_t o) {fBayesianIterIn = i; fBayesianIterOut = o;}
109 void SetBayesianSmooth(Float_t s) {fBayesianSmoothIn = s; fBayesianSmoothOut = s;}
110 void SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn = i; fBayesianSmoothOut = o;}
111 void SetAvoidRoundingError(Bool_t r) {fAvoidRoundingError = r;}
112 void SetUnfoldingAlgorithm(unfoldingAlgorithm ua) {fUnfoldingAlgorithm = ua;}
113 void SetPrior(prior p) {fPrior = p;}
114 void SetPrior(prior p, TH1D* spectrum) {fPrior = p; fPriorUser = spectrum;}
115 void SetNormalizeSpectra(Bool_t b) {fNormalizeSpectra = b;}
116 void SetNormalizeSpectra(Int_t e) { // use to normalize to this no of events
118 fNormalizeSpectra = kFALSE;
120 void SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
125 fSmoothenCounts = counts;
127 void SetTestMode(Bool_t t) {fTestMode = t;}
128 void SetEventPlaneResolution(Double_t r) {fEventPlaneRes = r;}
129 void SetUseDetectorResponse(Bool_t r) {fUseDetectorResponse = r;}
130 void SetUseDptResponse(Bool_t r) {fUseDptResponse = r;}
131 void SetTrainPowerFit(Bool_t t) {fTrainPower = t;}
132 void SetDphiUnfolding(Bool_t i) {fDphiUnfolding = i;}
133 void SetDphiDptUnfolding(Bool_t i) {fDphiDptUnfolding = i;}
134 void SetExLJDpt(Bool_t i) {fExLJDpt = i;}
135 void SetWeightFunction(TF1* w) {fResponseMaker->SetRMMergeWeightFunction(w);}
136 void SetTreatCorrErrAsUncorrErr(Bool_t b) {fSetTreatCorrErrAsUncorrErr = b;}
138 void MakeAU(); // test function, use with caution (09012014)
141 if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
142 if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
143 if(fRMSRatio) fRMSRatio->Write();
144 fOutputFile->Close();}
148 Float_t rangeLow = 20,
149 Float_t rangeUp = 80,
150 TString in = "UnfoldedSpectra.root",
151 TString out = "ProcessedSpectra.root") const;
152 void GetNominalValues(
157 TString inFile = "UnfoldedSpectra.root",
158 TString outFile = "Nominal.root") const;
159 void GetCorrelatedUncertainty(
160 TGraphAsymmErrors*& corrRatio,
161 TGraphAsymmErrors*& corrV2,
162 TArrayI* variationsIn,
163 TArrayI* variationsOut,
164 TArrayI* variantions2ndIn,
165 TArrayI* variantions2ndOut,
168 Float_t rangeLow = 20,
169 Float_t rangeUp = 80,
170 TString in = "UnfoldedSpectra.root",
171 TString out = "CorrelatedUncertainty.root") const;
172 void GetShapeUncertainty(
173 TGraphAsymmErrors*& shapeRatio,
174 TGraphAsymmErrors*& shapeV2,
175 TArrayI* regularizationIn,
176 TArrayI* regularizationOut,
177 TArrayI* trueBinIn = 0x0,
178 TArrayI* trueBinOut = 0x0,
179 TArrayI* recBinIn = 0x0,
180 TArrayI* recBinOut = 0x0,
182 Float_t rangeLow = 20,
183 Float_t rangeUp = 80,
184 TString in = "UnfoldedSpectra.root",
185 TString out = "ShapeUncertainty.root") const;
187 TH2D* detectorResponse, // detector response matrix
188 TH1D* jetPtIn, // in plane jet spectrum
189 TH1D* jetPtOut, // out of plane jet spectrum
190 TH1D* dptIn, // in plane delta pt distribution
191 TH1D* dptOut, // out of plane delta pt distribution
192 Int_t eventCount = 0); // event count (optional)
193 // static const helper functions, mainly histogram manipulation
194 static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
195 static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
196 static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
197 static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
198 TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
199 static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
200 static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
201 static TGraphErrors* GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
202 static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = 0., TString name = "");
203 TGraphAsymmErrors* GetV2WithSystematicErrors(
204 TH1* h1, TH1* h2, Double_t r, TString name,
205 TH1* relativeErrorInUp,
206 TH1* relativeErrorInLow,
207 TH1* relativeErrorOutUp,
208 TH1* relativeErrorOutLow,
209 Float_t rho = 0.) const;
210 static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
211 static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
212 static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
213 void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
214 static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
215 static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
217 void SetTitleFontSize(Double_t s) {fTitleFontSize = s;}
219 static void Style(TCanvas* c, TString style = "PEARSON");
220 static void Style(TVirtualPad* c, TString style = "SPECTRUM");
221 static void Style(TLegend* l);
222 static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty);
223 static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty);
224 static TLegend* AddLegend(TVirtualPad* p, Bool_t style = kFALSE) {
225 if(!style) return p->BuildLegend(.565, .663, .882, .883);
227 TLegend* l = AddLegend(p, kFALSE);
232 static TPaveText* AddTPaveText(TString text, Int_t r = 2) {
233 TPaveText* t(new TPaveText(.35, .27, .76, .33,"NDC"));
234 // t->SetFillStyle(0);
237 t->AddText(0.,0.,text.Data());
238 t->AddText(0., 0., Form("#it{R} = 0.%i #it{k}_{T} charged jets", r));
239 t->SetTextColor(kBlack);
240 // t->SetTextSize(0.03);
245 static void SavePadToPDF(TVirtualPad* pad) {pad->SaveAs(Form("%s.pdf", pad->GetName()));}
246 // interface to AliUnfolding, not necessary but nice to have all parameters in one place
247 static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
248 static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
249 static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
250 static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
251 static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
253 Bool_t PrepareForUnfolding();
254 Bool_t PrepareForUnfolding(Int_t low, Int_t up);
255 TH1D* GetPrior( const TH1D* measuredJetSpectrum,
256 const TH2D* resizedResponse,
257 const TH1D* kinematicEfficiency,
258 const TH1D* measuredJetSpectrumTrueBins,
259 const TString suffix,
260 const TH1D* jetFindingEfficiency);
261 TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
262 const TH2D* resizedResponse,
263 const TH1D* kinematicEfficiency,
264 const TH1D* measuredJetSpectrumTrueBins,
265 const TString suffix,
266 const TH1D* jetFindingEfficiency = 0x0);
267 TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
268 const TH2D* resizedResponse,
269 const TH1D* kinematicEfficiency,
270 const TH1D* measuredJetSpectrumTrueBins,
271 const TString suffix,
272 const TH1D* jetFindingEfficiency = 0x0);
273 TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
274 const TH2D* resizedResponse,
275 const TH1D* kinematicEfficiency,
276 const TH1D* measuredJetSpectrumTrueBins,
277 const TString suffix,
278 const TH1D* jetFindingEfficiency = 0x0);
279 TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
280 const TH2D* resizedResponse,
281 const TH1D* kinematicEfficiency,
282 const TH1D* measuredJetSpectrumTrueBins,
283 const TString suffix,
284 const TH1D* jetFindingEfficiency = 0x0);
285 TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
286 const TH2D* resizedResponse,
287 const TH1D* kinematicEfficiency,
288 const TH1D* measuredJetSpectrumTrueBins,
289 const TString suffix,
290 const TH1D* jetFindingEfficiency = 0x0);
291 void DoIntermediateSystematics(
292 TArrayI* variationsIn,
293 TArrayI* variationsOut,
294 TH1D*& relativeErrorInUp,
295 TH1D*& relativeErrorInLow,
296 TH1D*& relativeErrorOutUp,
297 TH1D*& relativeErrorOutLow,
298 TH1D*& relativeSystematicIn,
299 TH1D*& relativeSystematicOut,
307 TString source = "") const;
308 static void ResetAliUnfolding();
309 // give object a unique name via the 'protect heap' functions.
310 // may seem redundant, but some internal functions of root (e.g.
311 // ProjectionY()) check for existing objects by name and re-use them
312 TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
313 TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
314 TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
315 // members, accessible via setters
316 AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
317 TF1* fPower; // smoothening fit
318 Bool_t fSaveFull; // save all generated histograms to file
319 TString fActiveString; // identifier of active output
320 TDirectoryFile* fActiveDir; // active directory
321 TList* fInputList; // input list
322 Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
323 TString fOutputFileName; // output file name
324 TFile* fOutputFile; // output file
325 TArrayI* fCentralityArray; // array of bins that are merged
326 TArrayD* fCentralityWeights; // array of centrality weights
327 TH2D* fDetectorResponse; // detector response
328 TH1D* fJetFindingEff; // jet finding efficiency
329 Double_t fBetaIn; // regularization strength, in plane unfolding
330 Double_t fBetaOut; // regularization strength, out of plane unfoldign
331 Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
332 Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
333 Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
334 Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
335 Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
336 unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
337 prior fPrior; // prior for unfolding
338 TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
339 TArrayD* fBinsTrue; // pt true bins
340 TArrayD* fBinsRec; // pt rec bins
341 TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
342 TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
343 Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
344 Int_t fSVDRegOut; // svd regularization out of plane
345 Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
346 Float_t fJetRadius; // jet radius (for SVD toy)
347 Int_t fEventCount; // number of events
348 Bool_t fNormalizeSpectra; // normalize spectra to event count
349 Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
350 Float_t fFitMin; // lower bound of smoothening fit
351 Float_t fFitMax; // upper bound of smoothening fit
352 Float_t fFitStart; // from this value, use smoothening
353 Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
354 Bool_t fTestMode; // unfold with unity response for testing
355 Bool_t fRawInputProvided; // input histograms provided, not read from file
356 Double_t fEventPlaneRes; // event plane resolution for current centrality
357 Bool_t fUseDetectorResponse; // add detector response to unfolding
358 Bool_t fUseDptResponse; // add dpt response to unfolding
359 Bool_t fTrainPower; // don't clear the params of fPower for call to Make
360 // might give more stable results, but possibly introduces
361 // a bias / dependency on previous iterations
362 Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
363 Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
364 Bool_t fExLJDpt; // exclude randon cones with leading jet
365 Bool_t fSetTreatCorrErrAsUncorrErr; // treat the correlated error as uncorrelated
366 Double_t fTitleFontSize; // title font size
367 // members, set internally
368 TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
369 TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
370 TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
371 TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
372 TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
373 TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
374 TH1D* fSpectrumIn; // in plane jet pt spectrum
375 TH1D* fSpectrumOut; // out of plane jet pt spectrum
376 TH1D* fDptInDist; // in plane dpt distribution
377 TH1D* fDptOutDist; // out of plane dpt distribution
378 TH2D* fDptIn; // in plane dpt matrix
379 TH2D* fDptOut; // out plane dpt matrix
380 TH2D* fFullResponseIn; // full response matrix, in plane
381 TH2D* fFullResponseOut; // full response matrix, out of plane
382 // copy and assignment
383 AliJetFlowTools(const AliJetFlowTools&); // not implemented
384 AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
387 //_____________________________________________________________________________