]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliJetFlowTools.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
from Ante Bilandzic:
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.h
CommitLineData
dc1455ee 1/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
3 /* $Id$ */
4
51e6bc5a 5/* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
6 * see implementation for additional information */
7
dc1455ee 8#ifndef AliJetFlowTools_H
9#define AliJetFlowTools_H
10
549b5f40 11// root forward declarations
4292ca60 12class TF1;
dc1455ee 13class TH1D;
14class TH2D;
d7ec324f 15class TCanvas;
dc1455ee 16class TString;
17class TArrayD;
53547ff2 18class TGraph;
dc1455ee 19class TGraphErrors;
20class TObjArray;
4292ca60 21// aliroot forward declarations
22class AliAnaChargedJetResponseMaker;
23class AliUnfolding;
24// root includes
dc1455ee 25#include "TMatrixD.h"
26#include "TList.h"
ad04a83c 27#include "TDirectoryFile.h"
dc1455ee 28#include "TFile.h"
4292ca60 29#include "TProfile.h"
d7ec324f 30#include "TVirtualPad.h"
a39e4b2b 31#include "TPaveText.h"
32#include "TLegend.h"
dc1455ee 33//_____________________________________________________________________________
34class AliJetFlowTools {
35 public:
36 AliJetFlowTools();
37 protected:
53547ff2
RAB		 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
dc1455ee 40 public:
51e6bc5a 41 // enumerators
549b5f40
RAB		 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
18698978 50 kPriorMeasured, // use measured spectrum as prior
51 kPriorPythia }; // use pythia spectrum as prior
549b5f40 52 enum histoType { // histogram identifier, only used internally
18698978 53 kInPlaneSpectrum, // default style for spectrum
53547ff2
RAB		 54 kOutPlaneSpectrum,
55 kUnfoldedSpectrum,
56 kFoldedSpectrum,
57 kMeasuredSpectrum,
18698978 58 kBar, // default style for bar histogram
a39e4b2b 59 kRatio, // default style for ratio
60 kV2, // default style for v2
18698978 61 kEmpty }; // default style
ef12d5a5 62 // setters, interface to the class
4292ca60 63 void SetSaveFull(Bool_t b) {fSaveFull = b;}
64 void SetInputList(TList* list) {
65 fInputList = list;
66 fRefreshInput = kTRUE;
67 }
dc1455ee 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
ad04a83c 71 if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
72 fOutputFile->cd(); // avoid nested dirs
ef12d5a5 73 if(name.EqualTo(fActiveString)) {
74 printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
75 name+="_2";
76 }
dc1455ee 77 fActiveString = name;
ad04a83c 78 fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
79 fActiveDir->cd();
dc1455ee 80 }
549b5f40
RAB		 81 void SetCentralityBin(Int_t bin) {fCentralityBin = bin;}
82 void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;}
83 void SetJetFindingEfficiency(TH1D* e) {fJetFindingEff = e;}
84 void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
85 void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
86 void SetBinsTruePrior(TArrayD* bins) {fBinsTruePrior = bins;}
87 void SetBinsRecPrior(TArrayD* bins) {fBinsRecPrior = bins;}
88 void SetSVDReg(Int_t r) {fSVDRegIn = r; fSVDRegOut = r;}
89 void SetSVDReg(Int_t in, Int_t out) {fSVDRegIn = in; fSVDRegOut = out;}
90 void SetSVDToy(Bool_t b, Float_t r) {fSVDToy = b; fJetRadius = r;}
91 void SetBeta(Double_t b) {fBetaIn = b; fBetaOut = b;}
92 void SetBeta(Double_t i, Double_t o) {fBetaIn = i; fBetaOut = o;}
93 void SetBayesianIter(Int_t i) {fBayesianIterIn = i; fBayesianIterOut = i;}
94 void SetBayesianIter(Int_t i, Int_t o) {fBayesianIterIn = i; fBayesianIterOut = o;}
95 void SetBayesianSmooth(Float_t s) {fBayesianSmoothIn = s; fBayesianSmoothOut = s;}
96 void SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn = i; fBayesianSmoothOut = o;}
97 void SetAvoidRoundingError(Bool_t r) {fAvoidRoundingError = r;}
98 void SetUnfoldingAlgorithm(unfoldingAlgorithm ua) {fUnfoldingAlgorithm = ua;}
99 void SetPrior(prior p) {fPrior = p;}
18698978 100 void SetPrior(prior p, TH1D* spectrum) {fPrior = p; fPriorUser = spectrum;}
549b5f40
RAB		 101 void SetNormalizeSpectra(Bool_t b) {fNormalizeSpectra = b;}
102 void SetNormalizeSpectra(Int_t e) { // use to normalize to this no of events
103 fEventCount = e;
104 fNormalizeSpectra = kFALSE;
4292ca60 105 }
549b5f40
RAB		 106 void SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
107 fSmoothenPrior = b;
4292ca60 108 fFitMin = min;
109 fFitMax = max;
110 fFitStart = start;
549b5f40 111 fSmoothenCounts = counts;
4292ca60 112 }
549b5f40 113 void SetTestMode(Bool_t t) {fTestMode = t;}
ef12d5a5 114 void SetEventPlaneResolution(Double_t r) {fEventPlaneRes = r;}
115 void SetUseDetectorResponse(Bool_t r) {fUseDetectorResponse = r;}
549b5f40
RAB		 116 void SetUseDptResponse(Bool_t r) {fUseDptResponse = r;}
117 void SetTrainPowerFit(Bool_t t) {fTrainPower = t;}
486fb24e 118 void SetDphiUnfolding(Bool_t i) {fDphiUnfolding = i;}
119 void SetDphiDptUnfolding(Bool_t i) {fDphiDptUnfolding = i;}
120 void SetExLJDpt(Bool_t i) {fExLJDpt = i;}
121 void SetWeightFunction(TF1* w) {fResponseMaker->SetRMMergeWeightFunction(w);}
2364cc42 122 void SetTreatCorrErrAsUncorrErr(Bool_t b) {fSetTreatCorrErrAsUncorrErr = b;}
dc1455ee 123 void Make();
486fb24e 124 void MakeAU(); // test function, use with caution (09012014)
4292ca60 125 void Finish() {
126 fOutputFile->cd();
5e11c41c 127 if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
128 if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
129 if(fRMSRatio) fRMSRatio->Write();
4292ca60 130 fOutputFile->Close();}
53547ff2
RAB		 131 void PostProcess(
132 TString def,
486fb24e 133 Int_t columns = 4,
87233f72 134 Float_t rangeLow = 20,
135 Float_t rangeUp = 80,
53547ff2 136 TString in = "UnfoldedSpectra.root",
486fb24e 137 TString out = "ProcessedSpectra.root") const;
a39e4b2b 138 void GetNominalValues(
139 TH1D*& ratio,
140 TGraphErrors*& v2,
141 TArrayI* in,
142 TArrayI* out,
143 TString inFile = "UnfoldedSpectra.root",
144 TString outFile = "Nominal.root") const;
18698978 145 void GetCorrelatedUncertainty(
a39e4b2b 146 TGraphAsymmErrors*& corrRatio,
147 TGraphAsymmErrors*& corrV2,
f3ba6c8e 148 TArrayI* variationsIn,
149 TArrayI* variationsOut,
98d5b614 150 TArrayI* variantions2ndIn,
151 TArrayI* variantions2ndOut,
18698978 152 TString type = "",
f3ba6c8e 153 Int_t columns = 4,
154 Float_t rangeLow = 20,
155 Float_t rangeUp = 80,
156 TString in = "UnfoldedSpectra.root",
18698978 157 TString out = "CorrelatedUncertainty.root") const;
158 void GetShapeUncertainty(
a39e4b2b 159 TGraphAsymmErrors*& shapeRatio,
160 TGraphAsymmErrors*& shapeV2,
18698978 161 TArrayI* regularizationIn,
162 TArrayI* regularizationOut,
163 TArrayI* trueBinIn = 0x0,
164 TArrayI* trueBinOut = 0x0,
165 TArrayI* recBinIn = 0x0,
166 TArrayI* recBinOut = 0x0,
167 Int_t columns = 4,
168 Float_t rangeLow = 20,
169 Float_t rangeUp = 80,
170 TString in = "UnfoldedSpectra.root",
171 TString out = "ShapeUncertainty.root") const;
4292ca60 172 Bool_t SetRawInput (
173 TH2D* detectorResponse, // detector response matrix
174 TH1D* jetPtIn, // in plane jet spectrum
175 TH1D* jetPtOut, // out of plane jet spectrum
176 TH1D* dptIn, // in plane delta pt distribution
177 TH1D* dptOut, // out of plane delta pt distribution
178 Int_t eventCount = 0); // event count (optional)
dc1455ee 179 // static const helper functions, mainly histogram manipulation
180 static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
181 static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
512ced40 182 static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
53547ff2 183 static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
4292ca60 184 TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
d7ec324f 185 static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
dc1455ee 186 static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
4292ca60 187 static TGraphErrors* GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
4e4f12b6 188 static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = .7, TString name = "");
a39e4b2b 189 TGraphAsymmErrors* GetV2WithSystematicErrors(
190 TH1* h1, TH1* h2, Double_t r, TString name,
191 TH1* relativeErrorInUp,
192 TH1* relativeErrorInLow,
193 TH1* relativeErrorOutUp,
194 TH1* relativeErrorOutLow) const;
549b5f40 195 static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
4292ca60 196 static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
ef12d5a5 197 static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
549b5f40 198 void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
d7ec324f 199 static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
549b5f40 200 static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
18698978 201 // set style
202 void SetTitleFontSize(Double_t s) {fTitleFontSize = s;}
203 static void Style();
d7ec324f 204 static void Style(TCanvas* c, TString style = "PEARSON");
205 static void Style(TVirtualPad* c, TString style = "SPECTRUM");
53547ff2
RAB		 206 static void Style(TLegend* l);
207 static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty);
208 static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty);
a39e4b2b 209 static TLegend* AddLegend(TVirtualPad* p, Bool_t style = kFALSE) {
210 if(!style) return p->BuildLegend(.565, .663, .882, .883);
211 else {
212 TLegend* l = AddLegend(p, kFALSE);
213 Style(l);
214 return l;
215 }
216 }
98d5b614 217 static TPaveText* AddTPaveText(TString text, Int_t r = 2) {
a39e4b2b 218 TPaveText* t(new TPaveText(.35, .27, .76, .33,"NDC"));
219// t->SetFillStyle(0);
220 t->SetFillColor(0);
221 t->SetBorderSize(0);
222 t->AddText(0.,0.,text.Data());
98d5b614 223 t->AddText(0., 0., Form("#it{R} = 0.%i #it{k}_{T} charged jets", r));
a39e4b2b 224 t->SetTextColor(kBlack);
225// t->SetTextSize(0.03);
226 t->SetTextFont(42);
227 t->Draw("same");
228 return t;
229 }
512ced40 230 static void SavePadToPDF(TVirtualPad* pad) {pad->SaveAs(Form("%s.pdf", pad->GetName()));}
4292ca60 231 // interface to AliUnfolding, not necessary but nice to have all parameters in one place
232 static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
233 static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
234 static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
235 static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
236 static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
dc1455ee 237 private:
238 Bool_t PrepareForUnfolding();
486fb24e 239 Bool_t PrepareForUnfolding(Int_t low, Int_t up);
549b5f40
RAB		 240 TH1D* GetPrior( const TH1D* measuredJetSpectrum,
241 const TH2D* resizedResponse,
242 const TH1D* kinematicEfficiency,
243 const TH1D* measuredJetSpectrumTrueBins,
244 const TString suffix,
245 const TH1D* jetFindingEfficiency);
486fb24e 246 TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
247 const TH2D* resizedResponse,
248 const TH1D* kinematicEfficiency,
249 const TH1D* measuredJetSpectrumTrueBins,
250 const TString suffix,
251 const TH1D* jetFindingEfficiency = 0x0);
549b5f40
RAB		 252 TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
253 const TH2D* resizedResponse,
254 const TH1D* kinematicEfficiency,
255 const TH1D* measuredJetSpectrumTrueBins,
256 const TString suffix,
257 const TH1D* jetFindingEfficiency = 0x0);
258 TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
259 const TH2D* resizedResponse,
260 const TH1D* kinematicEfficiency,
261 const TH1D* measuredJetSpectrumTrueBins,
262 const TString suffix,
263 const TH1D* jetFindingEfficiency = 0x0);
264 TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
265 const TH2D* resizedResponse,
266 const TH1D* kinematicEfficiency,
267 const TH1D* measuredJetSpectrumTrueBins,
268 const TString suffix,
269 const TH1D* jetFindingEfficiency = 0x0);
270 TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
271 const TH2D* resizedResponse,
272 const TH1D* kinematicEfficiency,
273 const TH1D* measuredJetSpectrumTrueBins,
274 const TString suffix,
275 const TH1D* jetFindingEfficiency = 0x0);
18698978 276 void DoIntermediateSystematics(
277 TArrayI* variationsIn,
278 TArrayI* variationsOut,
279 TH1D*& relativeErrorInUp,
280 TH1D*& relativeErrorInLow,
281 TH1D*& relativeErrorOutUp,
282 TH1D*& relativeErrorOutLow,
283 TH1D*& relativeSystematicIn,
284 TH1D*& relativeSystematicOut,
a39e4b2b 285 TH1D*& nominal,
286 TH1D*& nominalIn,
287 TH1D*& nominalOut,
18698978 288 Int_t columns,
289 Float_t rangeLow,
290 Float_t rangeUp,
291 TFile* readMe,
292 TString source = "") const;
4292ca60 293 static void ResetAliUnfolding();
ef12d5a5 294 // give object a unique name via the 'protect heap' functions.
295 // may seem redundant, but some internal functions of root (e.g.
296 // ProjectionY()) check for existing objects by name and re-use them
f3ba6c8e 297 TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
298 TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
299 TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
dc1455ee 300 // members, accessible via setters
4292ca60 301 AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
302 TF1* fPower; // smoothening fit
303 Bool_t fSaveFull; // save all generated histograms to file
51e6bc5a 304 TString fActiveString; // identifier of active output
305 TDirectoryFile* fActiveDir; // active directory
306 TList* fInputList; // input list
4292ca60 307 Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
51e6bc5a 308 TString fOutputFileName; // output file name
309 TFile* fOutputFile; // output file
310 Int_t fCentralityBin; // centrality bin
51e6bc5a 311 TH2D* fDetectorResponse; // detector response
53547ff2 312 TH1D* fJetFindingEff; // jet finding efficiency
4292ca60 313 Double_t fBetaIn; // regularization strength, in plane unfolding
314 Double_t fBetaOut; // regularization strength, out of plane unfoldign
549b5f40
RAB		 315 Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
316 Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
317 Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
318 Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
51e6bc5a 319 Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
320 unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
321 prior fPrior; // prior for unfolding
18698978 322 TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
51e6bc5a 323 TArrayD* fBinsTrue; // pt true bins
324 TArrayD* fBinsRec; // pt rec bins
ef12d5a5 325 TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
326 TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
4292ca60 327 Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
328 Int_t fSVDRegOut; // svd regularization out of plane
51e6bc5a 329 Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
330 Float_t fJetRadius; // jet radius (for SVD toy)
20abfcc4 331 Int_t fEventCount; // number of events
332 Bool_t fNormalizeSpectra; // normalize spectra to event count
549b5f40 333 Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
4292ca60 334 Float_t fFitMin; // lower bound of smoothening fit
335 Float_t fFitMax; // upper bound of smoothening fit
336 Float_t fFitStart; // from this value, use smoothening
549b5f40 337 Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
ef12d5a5 338 Bool_t fTestMode; // unfold with unity response for testing
4292ca60 339 Bool_t fRawInputProvided; // input histograms provided, not read from file
ef12d5a5 340 Double_t fEventPlaneRes; // event plane resolution for current centrality
341 Bool_t fUseDetectorResponse; // add detector response to unfolding
549b5f40 342 Bool_t fUseDptResponse; // add dpt response to unfolding
ef12d5a5 343 Bool_t fTrainPower; // don't clear the params of fPower for call to Make
486fb24e 344 // might give more stable results, but possibly introduces
ef12d5a5 345 // a bias / dependency on previous iterations
486fb24e 346 Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
347 Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
348 Bool_t fExLJDpt; // exclude randon cones with leading jet
2364cc42 349 Bool_t fSetTreatCorrErrAsUncorrErr; // treat the correlated error as uncorrelated
18698978 350 Double_t fTitleFontSize; // title font size
dc1455ee 351 // members, set internally
4292ca60 352 TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
353 TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
354 TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
ef12d5a5 355 TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
4292ca60 356 TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
357 TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
358 TH1D* fSpectrumIn; // in plane jet pt spectrum
359 TH1D* fSpectrumOut; // out of plane jet pt spectrum
360 TH1D* fDptInDist; // in plane dpt distribution
361 TH1D* fDptOutDist; // out of plane dpt distribution
362 TH2D* fDptIn; // in plane dpt matrix
363 TH2D* fDptOut; // out plane dpt matrix
364 TH2D* fFullResponseIn; // full response matrix, in plane
365 TH2D* fFullResponseOut; // full response matrix, out of plane
dc1455ee 366 // copy and assignment
367 AliJetFlowTools(const AliJetFlowTools&); // not implemented
368 AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
369};
370#endif
371//_____________________________________________________________________________