.so cleanup: more gSystem->Load()
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / jetFlowTools.C
CommitLineData
dc1455ee 1void jetFlowTools() {
2 // load and compile the libraries
3 Load();
4
a39e4b2b 5
6 const Int_t SVDbestValueIn = 5;
7 const Int_t SVDbestValueOut = 4;
8 const Double_t bestBetaIn = 1.25;
9 const Double_t bestBetaOut = 1.25;
10
11 Bool_t runUnfolding = 0;
12 Bool_t doSystematics = (!runUnfolding);
13
14 // read detector response from output of matching taks
dc1455ee 15 // AliAnalysisTaskJetMatching
16 TString drInputName = "response.root";
17 printf("- Reading file %s ... \n", drInputName.Data());
18 TFile drInput(drInputName.Data()); // detector response input matrix
19 if(drInput.IsZombie()) {
20 printf(" > read error ! < \n");
21 return;
22 }
23 TH2D* detres = (TH2D*)drInput.Get("detector_response");
24 if(!detres) {
25 printf(" > failed to extract detector respose < \n");
26 return;
27 } else printf(" > Found detector response < \n");
28
0396b3c6 29 // get a TList from the AliAnalysisJetV2 task
dc1455ee 30 TFile f("AnalysisResults.root");
31 if(f.IsZombie()) {
32 printf(" > read error ! < \n");
33 return;
34 }
a39e4b2b 35 TList* l = (TList*)f.Get("RhoVnMod_R02_kCombined_Jet_AKTChargedR020_PicoTracks_pT0150_pt_scheme_TpcRho_ExLJ_TPC_PWGJE");
dc1455ee 36 if(!l) {
37 printf(" > failed to find output list ! \n");
38 return;
39 }
ad04a83c 40 // create an instance of the Tools class
dc1455ee 41 AliJetFlowTools* tools = new AliJetFlowTools();
ef12d5a5 42
a39e4b2b 43 // set some common variables
44 tools->SetCentralityBin(1);
45 tools->SetDetectorResponse(detres);
ef12d5a5 46 // set the true (unfolded) bins
a39e4b2b 47 Double_t binsTrue[] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 170};
ef12d5a5 48 tools->SetBinsTrue(new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue));
ef12d5a5 49 // set the measured (folded) bins
a39e4b2b 50 Double_t binsRec[] = {20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
ef12d5a5 51 tools->SetBinsRec(new TArrayD(sizeof(binsRec)/sizeof(binsRec[0]), binsRec));
a39e4b2b 52 // connect input
51e6bc5a 53 tools->SetInputList(l);
a39e4b2b 54
55 if(runUnfolding) {
56 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kNone);
57 tools->CreateOutputList(TString("DONOTHING"));
58 tools->Make();
59 // optional: smoothen the spectrum
60 tools->SetSaveFull(kTRUE);
61 tools->SetSmoothenPrior(kFALSE, 50, 250, 70, kFALSE);
62 // optional: normalize the spectrum
63 tools->SetUseDetectorResponse(kTRUE);
64 // optional: save a lot of raw output
65 tools->SetExLJDpt(kTRUE);
66 // do some /chi2 unfolding
67 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kChi2);
68 // first step: fishnet, see what good unfolding regularizations are
69 Double_t beta[] = {.001, .01 .1, .25, 1.25};
70 Int_t kReg[] = {9, 8, 7, 6, 5, 4, 3, 5};
71 // for out
72 Double_t betaOut[] = {.001, .01, .1, .25, 1.25};
73 Int_t kRegOut[] = {9, 8, 7, 6, 5, 4, 3, 4};
74 for(Int_t b(0); b < sizeof(beta)/sizeof(beta[0]); b++) {
75 tools->CreateOutputList(TString(Form("#beta = %.4f", beta[b])));
76 tools->SetBeta(beta[b], betaOut[b]);
77 tools->Make();
78 }
79 tools->SetPrior(AliJetFlowTools::kPriorChi2);
80 for(Int_t j(0); j < sizeof(kReg)/sizeof(kReg[0]); j++) {
81 // do some SVD unfolding
82 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kSVD);
83 tools->CreateOutputList(TString(Form("SVD kReg %i", kReg[j])));
84 tools->SetSVDReg(kReg[j]);
85 tools->Make();
86 }
87 // after fishnet:
88 // here we change the pt binning, using optimal svd and beta values
89 tools->SetBeta(bestBetaIn, bestBetaOut);
90 tools->SetSVDReg(SVDbestValueIn, SVDbestValueOut);
91 // ###### change the true (unfolded) binning ########
92 Double_t binsTrue2[] = {5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 170};
93 tools->SetBinsTrue(new TArrayD(sizeof(binsTrue2)/sizeof(binsTrue2[0]), binsTrue2));
94 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kSVD);
95 tools->CreateOutputList(TString("true bin removed"));
96 tools->Make();
97 // revert the true bin settings to their default ones
98 tools->SetBinsTrue(new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue));
99 // ####### change the measured binning
100 // remove a bin at low pt
101 Double_t binsRech[] = {25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
102 tools->SetBinsRec(new TArrayD(sizeof(binsRech)/sizeof(binsRech[0]), binsRech));
103 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kSVD);
104 tools->CreateOutputList(TString("measured bin removed"));
105 tools->Make();
106 // add a bin at low pt
107 Double_t binsRecl[] = {15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
108 tools->SetBinsRec(new TArrayD(sizeof(binsRecl)/sizeof(binsRecl[0]), binsRecl));
109 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kSVD);
110 tools->CreateOutputList(TString("measured bin added"));
111 tools->Make();
112
113 // revert to the original binning
114 tools->SetBinsTrue(new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue));
115 tools->SetBinsRec(new TArrayD(sizeof(binsRec)/sizeof(binsRec[0]), binsRec));
5e11c41c 116
a39e4b2b 117 // unfold using different tracking efficiency
118 TString drInputName95 = "/home/rbertens/Documents/CERN/jet-flow/RESPONSE/R02/95_pct_efficiency/5gev_leading_track_bias/response.root";
119 TFile drInput95(drInputName95.Data());
120 if(drInput95.IsZombie()) return;
121 TH2D* detres95 = (TH2D*)drInput95.Get("detector_response");
122 tools->SetDetectorResponse(detres95);
123 tools->CreateOutputList(TString("diced response"));
ef12d5a5 124 tools->Make();
5e11c41c 125
a39e4b2b 126 // switch back to the original detector response
127 tools->SetDetectorResponse(detres);
128
129
130 // now do the svd unfolding with a pythia spectrum as a prior
131 tools->SetUnfoldingAlgorithm(AliJetFlowTools::kSVD);
132 gROOT->LoadMacro("/home/rbertens/Documents/CERN/jet-flow/RESPONSE/pythia.C");
133 tools->SetPrior(AliJetFlowTools::kPriorPythia, pythia());
134 tools->CreateOutputList(TString("pythia prior"));
135 tools->Make();
136
137
138 // ########### systematics are done ################
139 // write the output to file
140 tools->Finish();
141 // do some postprocessing
142 tools->PostProcess(TString("SVD kReg 4"), 4, 20, 100);
143 } // end of run unfolding
144
145 if(doSystematics) {
146 /**
147 * evaluate systematics */
148 // first element of array should point to the nominal value
149 const Int_t reg[] = {9, 4, 8, 10, 16};
150 TArrayI* regArray = new TArrayI(sizeof(reg)/sizeof(reg[0]), reg);
151 const Int_t rec[] = {9, 12};
152 TArrayI* recArray = new TArrayI(sizeof(rec)/sizeof(rec[0]), rec);
153 const Int_t tru[] = {9, 13, 14};
154 TArrayI* truArray = new TArrayI(sizeof(tru)/sizeof(tru[0]), tru);
155
156
157 // place holder pointers. these will be assigned by using pointer references in the relevant functions
158 //
159 // for the nominal points
160 TH1D* nominalRatio (0x0);
161 TGraphErrors* nominalV2 (0x0);
162
163 // for the shape uncertainty
164 TGraphAsymmErrors* shapeRatio (0x0);
165 TGraphAsymmErrors* shapeV2 (0x0);
166
167 // for the correlated uncertainty
168 TGraphAsymmErrors* corrRatio (0x0);
169 TGraphAsymmErrors* corrV2 (0x0);
170 // get the actual values
171
172 tools->GetNominalValues(
173 nominalRatio,
174 nominalV2,
175 regArray, // doesn't matter which array is passed, as long as first element points to nominal value
176 regArray);
177
178 tools->GetShapeUncertainty(
179 shapeRatio,
180 shapeV2,
181 regArray, // systematics from regularization
182 regArray,
183 recArray, // from true spectrum variation
184 recArray,
185 truArray, // from rec spectrum variation
186 truArray,
187 4,
188 20,
189 100);
190
191 const Int_t cor[] = {9, 15};
192 TArrayI* corArray = new TArrayI(sizeof(cor)/sizeof(cor[0]), cor);
193
194 tools->GetCorrelatedUncertainty(
195 corrRatio,
196 corrV2,
197 corArray, // correlated systematics
198 corArray,
199 "diced respose", // name of systematic source
200 4,
201 20,
202 100);
203
204
205
206
207
208 using namespace AliJetFlowTools;
209 Double_t rangeLow(20.);
210 Double_t rangeHigh(90.);
211
212 TFile FinalResults = TFile("FinalResults.root", "RECREATE");
213 // combine the final results and write them to a file
214 TCanvas* full = new TCanvas("full", "full");
215 full->Divide(2);
216 full->cd(1);
217 AliJetFlowTools::Style(gPad, "PEARSON");
218 // shape uncertianty, full boxes
219 Style(shapeRatio, kYellow, AliJetFlowTools::kRatio);
220 shapeRatio->SetTitle("shape uncertainty");
221 shapeRatio->GetXaxis()->SetRangeUser(rangeLow, rangeHigh);
222 shapeRatio->GetYaxis()->SetRangeUser(.7, 2.2);
223 shapeRatio->DrawClone("a2");
224
225 // correlated uncertainty, open boxes
226 Style(corrRatio, kGray, AliJetFlowTools::kRatio);
227 corrRatio->SetTitle("correlated uncertainty");
228 corrRatio->SetFillStyle(0);
229 corrRatio->SetFillColor(kWhite);
230 corrRatio->DrawClone("5");
231
232 // ratio itself
233 Style(nominalRatio, kBlack, AliJetFlowTools::kRatio);
234 nominalRatio->DrawCopy("same E1");
235 nominalRatio->SetTitle("in / out of plane jet yield");
236
237 AddTPaveText("0-10 \% cc Pb-Pb #sqrt{s_{NN}} = 2.76 TeV");
238 AliJetFlowTools::AddLegend(gPad, kTRUE);
239 full->cd(2);
240 AliJetFlowTools::Style(gPad, "PEARSON");
241
242 // shape uncertainto on v2
243 Style(shapeV2, kYellow, AliJetFlowTools::kV2);
244 shapeV2->SetTitle("shape uncertainty");
245 shapeV2->GetXaxis()->SetRangeUser(rangeLow, rangeHigh);
246 shapeV2->GetYaxis()->SetRangeUser(-.5, 1.);
247 shapeV2->DrawClone("a2");
248
249 // correlated uncertainty
250 Style(corrV2, kGray, AliJetFlowTools::kV2);
251 corrV2->SetFillColor(kWhite);
252 corrV2->SetLineStyle(0);
253 corrV2->SetFillStyle(0);
254 corrV2->SetTitle("correlated uncertainty");
255 corrV2->DrawClone("5");
256
257 // v2 itself
258 Style(nominalV2, kBlack, AliJetFlowTools::kV2);
259 nominalV2->SetTitle("jet #it{v}_{2}");
260 nominalV2->SetFillColor(kWhite);
261 nominalV2->DrawClone("same E1");
262
263 AddTPaveText("0-10 \% cc Pb-Pb #sqrt{s_{NN}} = 2.76 TeV");
264 AliJetFlowTools::AddLegend(gPad, kTRUE);
265 gStyle->SetTitleStyle(0);
266 gStyle->SetStatStyle(0);
267
268 full->Write();
269 FinalResults.Close();
270 }
dc1455ee 271}
272
273//_____________________________________________________________________________
274void Load() {
4070f709 275 gSystem->Load("libTree");
276 gSystem->Load("libGeom");
277 gSystem->Load("libVMC");
dc1455ee 278 gSystem->Load("libPhysics");
279
4070f709 280 gSystem->Load("libSTEERBase");
281 gSystem->Load("libESD");
282 gSystem->Load("libAOD");
283 gSystem->Load("libANALYSIS");
284 gSystem->Load("libANALYSISalice");
285
286 gSystem->Load("libEMCALUtils");
287 gSystem->Load("libPHOSUtils");
288 gSystem->Load("libCGAL");
289 gSystem->Load("libfastjet");
290 gSystem->Load("libsiscone");
291 gSystem->Load("libSISConePlugin");
292
293 gSystem->Load("libCORRFW");
294 gSystem->Load("libPWGTools");
295 gSystem->Load("libJETAN");
296 gSystem->Load("libFASTJETAN");
297 gSystem->Load("libPWGJE");
dc1455ee 298
299 // include paths, necessary for compilation
300 gSystem->AddIncludePath("-Wno-deprecated");
301 gSystem->AddIncludePath("-I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/EMCAL");
302 gSystem->AddIncludePath("-I$ALICE_ROOT/PWGDQ/dielectron -I$ALICE_ROOT/PWGHF/hfe");
303 gSystem->AddIncludePath("-I$ALICE_ROOT/JETAN -I$ALICE_ROOT/JETAN/fastjet");
304
51e6bc5a 305 // attempt to load RooUnfold libraries,
b0635849 306 gSystem->Load("/home/rbertens/Documents/CERN/alice/BUILDS/ROOUNFOLD/RooUnfold-1.1.1/libRooUnfold");
a39e4b2b 307 gSystem->AddIncludePath("-I/home/rbertens/Documents/CERN/alice/BUILDS/ROOUNFOLD/RooUnfold-1.1.1/src/");
308 // compile unfolding class (only if there are local changes or the .o is not found)
5e11c41c 309 gROOT->LoadMacro("$ALICE_ROOT/PWG/FLOW/Tasks/AliJetFlowTools.cxx+");
dc1455ee 310}
311//_____________________________________________________________________________