379048b7ee20fe113185db8835eb64c382a51f8e
[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / jetFlowTools.C
1 void jetFlowTools() {
2     // load and compile the libraries
3     Load();
4
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
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
29     // get a TList from the AliAnalysisJetV2 task
30     TFile f("AnalysisResults.root");
31     if(f.IsZombie()) {
32         printf(" > read error ! < \n");
33         return;
34     }
35     TList* l = (TList*)f.Get("RhoVnMod_R02_kCombined_Jet_AKTChargedR020_PicoTracks_pT0150_pt_scheme_TpcRho_ExLJ_TPC_PWGJE");
36     if(!l) {
37         printf(" > failed to find output list ! \n");
38         return;
39     }
40     // create an instance of the Tools class
41     AliJetFlowTools* tools = new AliJetFlowTools();
42
43     // set some common variables
44     tools->SetCentralityBin(1);
45     tools->SetDetectorResponse(detres);
46     // set the true (unfolded) bins
47     Double_t binsTrue[] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 170};
48     tools->SetBinsTrue(new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue));
49     // set the measured (folded) bins
50     Double_t binsRec[] = {20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90};
51     tools->SetBinsRec(new TArrayD(sizeof(binsRec)/sizeof(binsRec[0]), binsRec));
52     // connect input
53     tools->SetInputList(l);
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));
116  
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"));
124         tools->Make();
125
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     }    
271 }
272
273 //_____________________________________________________________________________
274 void Load() {
275     gSystem->Load("libTree");
276     gSystem->Load("libGeom");
277     gSystem->Load("libVMC");
278     gSystem->Load("libPhysics");
279
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");
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
305     // attempt to load RooUnfold libraries, 
306     gSystem->Load("/home/rbertens/Documents/CERN/alice/BUILDS/ROOUNFOLD/RooUnfold-1.1.1/libRooUnfold.so");
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)
309     gROOT->LoadMacro("$ALICE_ROOT/PWG/FLOW/Tasks/AliJetFlowTools.cxx+");
310 }
311 //_____________________________________________________________________________