4474ad1f06df10eee33c847d6f7ef45af483fc4c
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.cxx
1 /************************************************************************* 
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * 
3  *                                                                        * 
4  * Author: The ALICE Off-line Project.                                    * 
5  * Contributors are mentioned in the code where appropriate.              * 
6  *                                                                        * 
7  * Permission to use, copy, modify and distribute this software and its   * 
8  * documentation strictly for non-commercial purposes is hereby granted   * 
9  * without fee, provided that the above copyright notice appears in all   * 
10  * copies and that both the copyright notice and this permission notice   * 
11  * appear in the supporting documentation. The authors make no claims     * 
12  * about the suitability of this software for any purpose. It is          * 
13  * provided "as is" without express or implied warranty.                  * 
14  **************************************************************************/
15
16 // Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
17 //         (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
18 //
19 // Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
20 //
21 // The task uses input from two analysis tasks:
22 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
23 //   used to retrieve jet spectra and delta pt distributions
24 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
25 //   used to construct the detector response function
26 // and unfolds jet spectra with respect to the event plane. The user can choose
27 // different alrogithms for unfolding which are available in (ali)root. RooUnfold 
28 // libraries must be present on the system 
29 // ( see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html ).
30 // 
31 // The weak spot of this class is the function PrepareForUnfolding, which will read
32 // output from two output files and expects histograms with certain names and binning. 
33 // Unfolding methods itself are general and should be able to handle any input, therefore one
34 // can forgo the PrepareForUnfolding method, and supply necessary input information via the 
35 // SetRawInput() method
36 //
37 // to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
38
39 // root includes
40 #include "TF1.h"
41 #include "TF2.h"
42 #include "TH1D.h"
43 #include "TH2D.h"
44 #include "TGraph.h"
45 #include "TGraphErrors.h"
46 #include "TGraphAsymmErrors.h"
47 #include "TLine.h"
48 #include "TCanvas.h"
49 #include "TLegend.h"
50 #include "TArrayD.h"
51 #include "TList.h"
52 #include "TMinuit.h"
53 #include "TVirtualFitter.h"
54 #include "TLegend.h"
55 #include "TCanvas.h"
56 #include "TStyle.h"
57 #include "TLine.h"
58 #include "TMath.h"
59 #include "TVirtualFitter.h"
60 #include "TFitResultPtr.h"
61 #include "Minuit2/Minuit2Minimizer.h"
62 #include "Math/Functor.h"
63 // aliroot includes
64 #include "AliUnfolding.h"
65 #include "AliAnaChargedJetResponseMaker.h"
66 // class includes
67 #include "AliJetFlowTools.h"
68 // roo unfold includes (make sure you have these available on your system)
69 #include "RooUnfold.h"
70 #include "RooUnfoldResponse.h"
71 #include "RooUnfoldSvd.h"
72 #include "RooUnfoldBayes.h"
73 #include "TSVDUnfold.h"
74
75 using namespace std;
76 //_____________________________________________________________________________
77 AliJetFlowTools::AliJetFlowTools() :
78     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
79     fRMS                (kTRUE),
80     fSymmRMS            (kTRUE),
81     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
82     fSaveFull           (kTRUE),
83     fActiveString       (""),
84     fActiveDir          (0x0),
85     fInputList          (0x0),
86     fRefreshInput       (kTRUE),
87     fOutputFileName     ("UnfoldedSpectra.root"),
88     fOutputFile         (0x0),
89     fCentralityArray    (0x0),
90     fMergeBinsArray     (0x0),
91     fCentralityWeights  (0x0),
92     fDetectorResponse   (0x0),
93     fJetFindingEff      (0x0),
94     fBetaIn             (.1),
95     fBetaOut            (.1),
96     fBayesianIterIn     (4),
97     fBayesianIterOut    (4),
98     fBayesianSmoothIn   (0.),
99     fBayesianSmoothOut  (0.),
100     fAvoidRoundingError (kFALSE),
101     fUnfoldingAlgorithm (kChi2),
102     fPrior              (kPriorMeasured),
103     fPriorUser          (0x0),
104     fBinsTrue           (0x0),
105     fBinsRec            (0x0),
106     fBinsTruePrior      (0x0),
107     fBinsRecPrior       (0x0),
108     fSVDRegIn           (5),
109     fSVDRegOut          (5),
110     fSVDToy             (kTRUE),
111     fJetRadius          (0.3),
112     fEventCount         (-1),
113     fNormalizeSpectra   (kFALSE),
114     fSmoothenPrior      (kFALSE),
115     fFitMin             (60.),
116     fFitMax             (300.),
117     fFitStart           (75.),
118     fSmoothenCounts     (kTRUE),
119     fTestMode           (kFALSE),
120     fRawInputProvided   (kFALSE),
121     fEventPlaneRes      (.63),
122     fUseDetectorResponse(kTRUE),
123     fUseDptResponse     (kTRUE),
124     fTrainPower         (kTRUE),
125     fDphiUnfolding      (kTRUE),
126     fDphiDptUnfolding   (kFALSE),
127     fExLJDpt            (kTRUE),
128     fTitleFontSize      (-999.),
129     fRMSSpectrumIn      (0x0),
130     fRMSSpectrumOut     (0x0),
131     fRMSRatio           (0x0),
132     fRMSV2              (0x0),
133     fDeltaPtDeltaPhi    (0x0),
134     fJetPtDeltaPhi      (0x0),
135     fSpectrumIn         (0x0),
136     fSpectrumOut        (0x0),
137     fDptInDist          (0x0),
138     fDptOutDist         (0x0),
139     fDptIn              (0x0),
140     fDptOut             (0x0),
141     fFullResponseIn     (0x0),
142     fFullResponseOut    (0x0) { // class constructor
143         // create response maker weight function (tuned to PYTHIA spectrum)
144         fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
145         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
146 }
147 //_____________________________________________________________________________
148 void AliJetFlowTools::Make() {
149     // core function of the class
150     if(fDphiDptUnfolding) {
151         // to extract the yield as function of Dphi, Dpt - experimental
152         MakeAU();
153         return;
154     }
155     // 1) rebin the raw output of the jet task to the desired binnings
156     // 2) calls the unfolding routine
157     // 3) writes output to file
158     // can be repeated multiple times with different configurations
159
160     // 1) manipulation of input histograms
161     // check if the input variables are present
162     if(fRefreshInput) {
163         if(!PrepareForUnfolding()) {
164             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
165             return;
166         }
167     }
168     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
169     //     parts of the spectrum can end up in over or underflow bins
170     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
171     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec,  TString("resized_out_"), kFALSE);
172     
173     // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
174     // the template will be used as a prior for the chi2 unfolding
175     TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
176     TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
177     // get the full response matrix from the dpt and the detector response
178     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
179     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
180     // so that unfolding should return the initial spectrum
181     if(!fTestMode) {
182         if(fUseDptResponse && fUseDetectorResponse) {
183             fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
184             fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
185         } else if (fUseDptResponse && !fUseDetectorResponse) {
186             fFullResponseIn = fDptIn;
187             fFullResponseOut = fDptOut;
188         } else if (!fUseDptResponse && fUseDetectorResponse) {
189             fFullResponseIn = fDetectorResponse;
190             fFullResponseOut = fDetectorResponse;
191         } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
192             printf(" > No response, exiting ! < \n" );
193             return;
194         }
195     } else {
196         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
197         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
198     }
199     // normalize each slide of the response to one
200     NormalizeTH2D(fFullResponseIn);
201     NormalizeTH2D(fFullResponseOut);
202     // resize to desired binning scheme
203     TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
204     TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
205     // get the kinematic efficiency
206     TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
207     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
208     TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
209     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
210     // suppress the errors 
211     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
212         kinematicEfficiencyIn->SetBinError(1+i, 0.);
213         kinematicEfficiencyOut->SetBinError(1+i, 0.);
214     }
215     TH1D* jetFindingEfficiency(0x0);
216     if(fJetFindingEff) {
217         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
218         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
219         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
220     }
221     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
222     TH1D* unfoldedJetSpectrumIn(0x0);
223     TH1D* unfoldedJetSpectrumOut(0x0); 
224     fActiveDir->cd();                   // select active dir
225     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
226     dirIn->cd();                        // select inplane subdir
227     // do the inplane unfolding
228     unfoldedJetSpectrumIn = UnfoldWrapper(
229         measuredJetSpectrumIn,
230         resizedResponseIn,
231         kinematicEfficiencyIn,
232         measuredJetSpectrumTrueBinsIn,
233         TString("in"),
234         jetFindingEfficiency);
235     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
236     resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
237     resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
238     resizedResponseIn = ProtectHeap(resizedResponseIn);
239     resizedResponseIn->Write();
240     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
241     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
242     kinematicEfficiencyIn->Write();
243     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
244     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
245     fDetectorResponse->Write();
246     // optional histograms
247     if(fSaveFull) {
248         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
249         fSpectrumIn->Write();
250         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
251         fDptInDist->Write();
252         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
253         fDptIn->Write();
254         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
255         fFullResponseIn->Write();
256     }
257     fActiveDir->cd();
258     if(fDphiUnfolding) {
259         TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
260         dirOut->cd();
261         // do the out of plane unfolding
262         unfoldedJetSpectrumOut = UnfoldWrapper(
263                     measuredJetSpectrumOut,
264                     resizedResponseOut,
265                     kinematicEfficiencyOut,
266                     measuredJetSpectrumTrueBinsOut,
267                     TString("out"),
268                     jetFindingEfficiency);
269         resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
270         resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
271         resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
272         resizedResponseOut = ProtectHeap(resizedResponseOut);
273         resizedResponseOut->Write();
274         kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
275         kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
276         kinematicEfficiencyOut->Write();
277         fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
278         fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
279         fDetectorResponse->Write();
280         if(jetFindingEfficiency) jetFindingEfficiency->Write();
281         // optional histograms
282         if(fSaveFull) {
283             fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
284             fSpectrumOut->Write();
285             fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
286             fDptOutDist->Write();
287             fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
288             fDptOut->Write();
289             fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
290             fFullResponseOut->Write();
291         }
292
293         // write general output histograms to file
294         fActiveDir->cd();
295         if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
296             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
297             if(ratio) {
298                 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
299                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
300                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
301                 ratio = ProtectHeap(ratio);
302                 ratio->Write();
303                 // write histo values to RMS files if both routines converged
304                 // input values are weighted by their uncertainty
305                 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
306                     if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
307                     if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
308                     if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
309                }
310             }
311             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
312             if(v2) {
313                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
314                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
315                 v2->GetYaxis()->SetTitle("v_{2}");
316                 v2 = ProtectHeap(v2);
317                 v2->Write();
318             }
319         } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
320             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
321             if(ratio) {
322                 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
323                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
324                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
325                 ratio = ProtectHeap(ratio);
326                 ratio->Write();
327             }
328             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
329              if(v2) {
330                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
331                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
332                 v2->GetYaxis()->SetTitle("v_{2}");
333                 v2 = ProtectHeap(v2);
334                 v2->Write();
335             }
336         }
337     }   // end of if(fDphiUnfolding)
338     fDeltaPtDeltaPhi->Write();
339     unfoldedJetSpectrumIn->Sumw2();
340     ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
341     unfoldedJetSpectrumIn->Write();
342     unfoldedJetSpectrumOut->Sumw2();
343     ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
344     unfoldedJetSpectrumOut->Write();
345     fJetPtDeltaPhi->Write();
346     // save the current state of the unfolding object
347     SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
348     TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
349     TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
350     unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
351     unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
352     unfoldedJetSpectrumInForSub->Write();
353
354 }
355 //_____________________________________________________________________________
356 TH1D* AliJetFlowTools::UnfoldWrapper(
357         const TH1D* measuredJetSpectrum,        // truncated raw jets (same binning as pt rec of response) 
358         const TH2D* resizedResponse,            // response matrix
359         const TH1D* kinematicEfficiency,        // kinematic efficiency
360         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
361         const TString suffix,                   // suffix (in or out of plane)
362         const TH1D* jetFindingEfficiency)       // jet finding efficiency
363 {
364     // wrapper function to call specific unfolding routine
365     TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
366     // initialize functon pointer
367     if(fUnfoldingAlgorithm == kChi2)                    myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
368     else if(fUnfoldingAlgorithm == kBayesian)           myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
369     else if(fUnfoldingAlgorithm == kBayesianAli)        myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
370     else if(fUnfoldingAlgorithm == kSVD)                myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
371     else if(fUnfoldingAlgorithm == kNone) {
372         TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
373         clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
374         return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
375     }
376     else return 0x0; 
377     // do the actual unfolding with the selected function
378     return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
379
380 //_____________________________________________________________________________
381 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
382         const TH1D* measuredJetSpectrum,      // truncated raw jets (same binning as pt rec of response) 
383         const TH2D* resizedResponse,          // response matrix
384         const TH1D* kinematicEfficiency,      // kinematic efficiency
385         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
386         const TString suffix,                 // suffix (in or out of plane)
387         const TH1D* jetFindingEfficiency)     // jet finding efficiency (optional)
388 {
389     // unfold the spectrum using chi2 minimization
390
391     // step 0) setup the static members of AliUnfolding
392     ResetAliUnfolding();                // reset from previous iteration
393                                         // also deletes and re-creates the global TVirtualFitter
394     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
395     if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
396     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
397     if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
398     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
399     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
400
401     // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
402     // stay intact. a local copy of a histogram (which only exists in the scope of this function) is 
403     // denoted by the suffix 'Local'
404     
405     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
406     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
407     // unfolded local will be filled with the result of the unfolding
408     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
409
410     // full response matrix and kinematic efficiency
411     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
412     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
413
414     // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
415     TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
416     // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
417     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
418
419     // step 2) start the unfolding
420     Int_t status(-1), i(0);
421     while(status < 0 && i < 100) {
422         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
423         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the 
424         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
425         status = AliUnfolding::Unfold(
426                 resizedResponseLocal,           // response matrix
427                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
428                 measuredJetSpectrumLocal,              // measured spectrum
429                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
430                 unfoldedLocal);                 // results
431         // status holds the minuit fit status (where 0 means convergence)
432         i++;
433     }
434     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
435     TH2D* hPearson(0x0);
436     if(status == 0 && gMinuit->fISW[1] == 3) {
437         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
438         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
439         if(gMinuit) gMinuit->Command("SET COV");
440         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
441         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
442         pearson->Print();
443         hPearson = new TH2D(*pearson);
444         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
445         hPearson = ProtectHeap(hPearson);
446         hPearson->Write();
447         if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
448     } else status = -1; 
449
450     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
451     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
452     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
453     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
454     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
455     if(ratio) {
456         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
457         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
458         ratio = ProtectHeap(ratio);
459         ratio->Write();
460     }
461
462     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap, 
463     // objects are cloned using 'ProtectHeap()'
464     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
465     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
466     measuredJetSpectrumLocal->Write(); 
467
468     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
469     resizedResponseLocal->Write();
470
471     unfoldedLocal = ProtectHeap(unfoldedLocal);
472     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
473     unfoldedLocal->Write();
474
475     foldedLocal = ProtectHeap(foldedLocal);
476     foldedLocal->Write();
477     
478     priorLocal = ProtectHeap(priorLocal);
479     priorLocal->Write();
480
481     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
482     TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 4, -0.5, 3.5));
483     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
484     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
485     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
486     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
487     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
488     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
489     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
490     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
491     fitStatus->Write();
492
493     return unfoldedLocal;
494 }
495 //_____________________________________________________________________________
496 TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
497         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
498         const TH2D* resizedResponse,                   // full response matrix, normalized in slides of pt true
499         const TH1D* kinematicEfficiency,              // kinematic efficiency
500         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
501         const TString suffix,                         // suffix (in, out)
502         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
503 {
504
505     TH1D* priorLocal( GetPrior(
506         measuredJetSpectrum,              // jet pt in pt rec bins 
507         resizedResponse,                  // full response matrix, normalized in slides of pt true
508         kinematicEfficiency,              // kinematic efficiency
509         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
510         suffix,                           // suffix (in, out)
511         jetFindingEfficiency));           // jet finding efficiency (optional)
512     if(!priorLocal) {
513         printf(" > couldn't find prior ! < \n");
514         return 0x0;
515     } else printf(" 1) retrieved prior \n");
516
517     // go back to the 'root' directory of this instance of the SVD unfolding routine
518     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
519
520     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
521     // measured jets in pt rec binning
522     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
523     // local copie of the response matrix
524     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
525     // local copy of response matrix, all true slides normalized to 1 
526     // this response matrix will eventually be used in the re-folding routine
527     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
528     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
529     // kinematic efficiency
530     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
531     // place holder histos
532     TH1D *unfoldedLocalSVD(0x0);
533     TH1D *foldedLocalSVD(0x0);
534     cout << " 2) setup necessary input " << endl;
535     // 3) configure routine
536     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
537     cout << " step 3) configured routine " << endl;
538
539     // 4) get transpose matrices
540     // a) get the transpose of the full response matrix
541     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
542     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
543     // normalize it with the prior. this will ensure that high statistics bins will constrain the
544     // end result most strenuously than bins with limited number of counts
545     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
546     cout << " 4) retrieved first transpose matrix " << endl;
547  
548     // 5) get response for SVD unfolding
549     RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
550     cout << " 5) retrieved roo unfold response object " << endl;
551
552     // 6) actualy unfolding loop
553     RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
554     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
555     // correct the spectrum for the kinematic efficiency
556     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
557
558     // get the pearson coefficients from the covariance matrix
559     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
560     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
561     TH2D* hPearson(0x0);
562     if(pearson) {
563         hPearson = new TH2D(*pearson);
564         pearson->Print();
565         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
566         hPearson = ProtectHeap(hPearson);
567         hPearson->Write();
568         if(fMergeBinsArray) unfoldedLocalSVD = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalSVD, hPearson);
569     } else return 0x0;       // return if unfolding didn't converge
570
571     // plot singular values and d_i vector
572     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
573     TH1* hSVal(svdUnfold->GetSV());
574     TH1D* hdi(svdUnfold->GetD());
575     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
576     hSVal->SetXTitle("singular values");
577     hSVal->Write();
578     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
579     hdi->SetXTitle("|d_{i}^{kreg}|");
580     hdi->Write();
581     cout << " plotted singular values and d_i vector " << endl;
582
583     // 7) refold the unfolded spectrum
584     foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
585     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
586     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
587     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
588     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
589     ratio->Write();
590     cout << " 7) refolded the unfolded spectrum " << endl;
591
592     // write the measured, unfolded and re-folded spectra to the output directory
593     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
594     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
595     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
596     measuredJetSpectrumLocal->Write(); // input spectrum
597     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
598     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
599     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
600     unfoldedLocalSVD->Write();  // unfolded spectrum
601     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
602     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
603     foldedLocalSVD->Write();    // re-folded spectrum
604
605     // save more general bookkeeeping histograms to the output directory
606     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
607     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
608     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
609     responseMatrixLocalTransposePrior->Write();
610     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
611     priorLocal->SetXTitle("p_{t} [GeV/c]");
612     priorLocal = ProtectHeap(priorLocal);
613     priorLocal->Write();
614     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
615     resizedResponseLocalNorm->Write();
616
617     // save some info 
618     TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 1, -0.5, 0.5));
619     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
620     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
621     fitStatus->Write();
622
623     return unfoldedLocalSVD;
624 }
625 //_____________________________________________________________________________
626 TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
627         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
628         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
629         const TH1D* kinematicEfficiency,              // kinematic efficiency
630         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
631         const TString suffix,                         // suffix (in, out)
632         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
633 {
634     // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
635     // FIXME careful, not tested yet ! (06122013) FIXME
636
637     // step 0) setup the static members of AliUnfolding
638     ResetAliUnfolding();                // reset from previous iteration
639                                         // also deletes and re-creates the global TVirtualFitter
640     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
641     if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
642     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
643     else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
644     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
645     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
646
647     // 1) get a prior for unfolding and clone all the input histograms
648     TH1D* priorLocal( GetPrior(
649         measuredJetSpectrum,              // jet pt in pt rec bins 
650         resizedResponse,                  // full response matrix, normalized in slides of pt true
651         kinematicEfficiency,              // kinematic efficiency
652         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
653         suffix,                           // suffix (in, out)
654         jetFindingEfficiency));           // jet finding efficiency (optional)
655     if(!priorLocal) {
656         printf(" > couldn't find prior ! < \n");
657         return 0x0;
658     } else printf(" 1) retrieved prior \n");
659     // switch back to root dir of this unfolding procedure
660     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
661
662     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
663     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
664     // unfolded local will be filled with the result of the unfolding
665     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
666
667     // full response matrix and kinematic efficiency
668     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
669     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
670
671     // step 2) start the unfolding
672     Int_t status(-1), i(0);
673     while(status < 0 && i < 100) {
674         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
675         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the 
676         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
677         status = AliUnfolding::Unfold(
678                 resizedResponseLocal,           // response matrix
679                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
680                 measuredJetSpectrumLocal,              // measured spectrum
681                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
682                 unfoldedLocal);                 // results
683         // status holds the minuit fit status (where 0 means convergence)
684         i++;
685     }
686     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
687     TH2D* hPearson(0x0);
688     if(status == 0 && gMinuit->fISW[1] == 3) {
689         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
690         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
691         if(gMinuit) gMinuit->Command("SET COV");
692         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
693         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
694         pearson->Print();
695         hPearson= new TH2D(*pearson);
696         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
697         hPearson = ProtectHeap(hPearson);
698         hPearson->Write();
699     } else status = -1; 
700     if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
701
702     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
703     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
704     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
705     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
706     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
707     if(ratio) {
708         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
709         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
710         ratio = ProtectHeap(ratio);
711         ratio->Write();
712     }
713
714     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap, 
715     // objects are cloned using 'ProtectHeap()'
716     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
717     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
718     measuredJetSpectrumLocal->Write(); 
719
720     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
721     resizedResponseLocal->Write();
722
723     unfoldedLocal = ProtectHeap(unfoldedLocal);
724     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
725     unfoldedLocal->Write();
726
727     foldedLocal = ProtectHeap(foldedLocal);
728     foldedLocal->Write();
729     
730     priorLocal = ProtectHeap(priorLocal);
731     priorLocal->Write();
732
733     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
734     TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 4, -0.5, 3.5));
735     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
736     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
737     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
738     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
739     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
740     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
741     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
742     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
743     fitStatus->Write();
744
745     return unfoldedLocal;
746 }
747 //_____________________________________________________________________________
748 TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
749         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
750         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
751         const TH1D* kinematicEfficiency,              // kinematic efficiency
752         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
753         const TString suffix,                         // suffix (in, out)
754         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
755 {
756     // use bayesian unfolding from the RooUnfold package to unfold jet spectra
757     
758     // 1) get a prior for unfolding.
759     TH1D* priorLocal( GetPrior(
760         measuredJetSpectrum,              // jet pt in pt rec bins 
761         resizedResponse,                  // full response matrix, normalized in slides of pt true
762         kinematicEfficiency,              // kinematic efficiency
763         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
764         suffix,                           // suffix (in, out)
765         jetFindingEfficiency));           // jet finding efficiency (optional)
766     if(!priorLocal) {
767         printf(" > couldn't find prior ! < \n");
768         return 0x0;
769     } else printf(" 1) retrieved prior \n");
770     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
771
772     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
773     // measured jets in pt rec binning
774     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
775     // local copie of the response matrix
776     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
777     // local copy of response matrix, all true slides normalized to 1 
778     // this response matrix will eventually be used in the re-folding routine
779     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
780     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
781     // kinematic efficiency
782     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
783     // place holder histos
784     TH1D *unfoldedLocalBayes(0x0);
785     TH1D *foldedLocalBayes(0x0);
786     cout << " 2) setup necessary input " << endl;
787     // 4) get transpose matrices
788     // a) get the transpose of the full response matrix
789     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
790     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
791     // normalize it with the prior. this will ensure that high statistics bins will constrain the
792     // end result most strenuously than bins with limited number of counts
793     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
794     // 3) get response for Bayesian unfolding
795     RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
796
797     // 4) actualy unfolding loop
798     RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
799     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
800     unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
801     // correct the spectrum for the kinematic efficiency
802     unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
803     // get the pearson coefficients from the covariance matrix
804     TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
805     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
806     TH2D* hPearson(0x0);
807     if(pearson) {
808         hPearson = new TH2D(*pearson);
809         pearson->Print();
810         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
811         hPearson = ProtectHeap(hPearson);
812         hPearson->Write();
813         if(fMergeBinsArray) unfoldedLocalBayes = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalBayes, hPearson);
814     } else return 0x0;       // return if unfolding didn't converge
815
816     // 5) refold the unfolded spectrum
817     foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
818     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio  measured / re-folded", kTRUE));
819     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
820     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
821     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
822     ratio->Write();
823     cout << " 7) refolded the unfolded spectrum " << endl;
824
825     // write the measured, unfolded and re-folded spectra to the output directory
826     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
827     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
828     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
829     measuredJetSpectrumLocal->Write(); // input spectrum
830     unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
831     unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
832     if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
833     unfoldedLocalBayes->Write();  // unfolded spectrum
834     foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
835     foldedLocalBayes = ProtectHeap(foldedLocalBayes);
836     foldedLocalBayes->Write();    // re-folded spectrum
837
838     // save more general bookkeeeping histograms to the output directory
839     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
840     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
841     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
842     responseMatrixLocalTransposePrior->Write();
843     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
844     priorLocal->SetXTitle("p_{t} [GeV/c]");
845     priorLocal = ProtectHeap(priorLocal);
846     priorLocal->Write();
847     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
848     resizedResponseLocalNorm->Write();
849
850     // save some info 
851     TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 1, -0.5, 0.5));
852     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
853     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
854     fitStatus->Write();
855
856     return  unfoldedLocalBayes; 
857 }
858 //_____________________________________________________________________________
859 Bool_t AliJetFlowTools::PrepareForUnfolding()
860 {
861     // prepare for unfolding
862     if(fRawInputProvided) return kTRUE;
863     if(!fInputList) {
864         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
865         return kFALSE;
866     }
867     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
868     // check if the pt bin for true and rec have been set
869     if(!fBinsTrue || !fBinsRec) {
870         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
871         return kFALSE;
872     }
873     if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
874                           // procedures, these profiles will be nonsensical, user is responsible
875         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
876         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
877         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
878     }
879     if(!fTrainPower) {
880         // clear minuit state to avoid constraining the fit with the results of the previous iteration
881         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
882     }
883     // extract the spectra 
884     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
885     if(!fInputList->FindObject(spectrumName.Data())) {
886         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
887         return kFALSE;
888     }
889
890     // get the first scaled spectrum
891     fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
892     // clone the spectrum on the heap. this is necessary since scale or add change the
893     // contents of the original histogram
894     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
895     fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
896     printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
897     // merge subsequent bins (if any)
898     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
899         spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
900         printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
901         fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
902     }
903
904     // in plane spectrum
905     if(!fDphiUnfolding) {
906         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
907         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
908     } else {
909         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
910         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
911         fSpectrumIn = ProtectHeap(fSpectrumIn);
912         // out of plane spectrum
913         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
914         fSpectrumOut = ProtectHeap(fSpectrumOut);
915     }
916     // normalize spectra to event count if requested
917     if(fNormalizeSpectra) {
918         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
919         rho->Scale(fCentralityWeights->At(0));
920         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
921            rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
922         }
923         if(!rho) return 0x0;
924         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
925         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
926         if(fEventCount > 0) {
927             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
928             fSpectrumOut->Sumw2();
929             Double_t pt(0);            
930             Double_t error(0); // lots of issues with the errors here ...
931             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
932                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
933                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
934                 fSpectrumIn->SetBinContent(1+i, pt);
935                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
936                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
937                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
938             }
939             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
940                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
941                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
942                 fSpectrumOut->SetBinContent(1+i, pt);
943                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
944                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
945                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
946             }
947         }
948         if(normalizeToFullSpectrum) fEventCount = -1;
949     }
950     // extract the delta pt matrices
951     TString deltaptName("");
952     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
953     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
954     if(!fDeltaPtDeltaPhi) {
955         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
956         printf(" > may be ok, depending no what you want to do < \n");
957         fRefreshInput = kTRUE;
958         return kTRUE;
959     }
960
961     // clone the distribution on the heap and if requested merge with other centralities
962     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
963     fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
964     printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
965     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
966         deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
967         printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
968         fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
969     }
970
971     // in plane delta pt distribution
972     if(!fDphiUnfolding) {
973         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
974         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
975     } else {
976         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
977         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
978         // out of plane delta pt distribution
979         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
980         fDptInDist = ProtectHeap(fDptInDist);
981         fDptOutDist = ProtectHeap(fDptOutDist);
982         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
983     }
984
985     // create a rec - true smeared response matrix
986     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
987     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
988         Bool_t skip = kFALSE;
989         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
990             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
991             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
992         }
993     }
994     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
995     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
996         Bool_t skip = kFALSE;
997         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
998             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
999             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1000         }
1001     }
1002     fDptIn = new TH2D(*rfIn);
1003     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1004     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1005     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1006     fDptIn = ProtectHeap(fDptIn);
1007     fDptOut = new TH2D(*rfOut);
1008     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
1009     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1010     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1011     fDptOut = ProtectHeap(fDptOut);
1012     
1013     fRefreshInput = kTRUE;     // force cloning of the input
1014     return kTRUE;
1015 }
1016 //_____________________________________________________________________________
1017 Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
1018     // prepare for unfoldingUA - more robust method to extract input spectra from file
1019     // will replace PrepareForUnfolding eventually (09012014)
1020     if(!fInputList) {
1021         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1022         return kFALSE;
1023     }
1024     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1025     // check if the pt bin for true and rec have been set
1026     if(!fBinsTrue || !fBinsRec) {
1027         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1028         return kFALSE;
1029     }
1030     if(!fTrainPower) {
1031         // clear minuit state to avoid constraining the fit with the results of the previous iteration
1032         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1033     }
1034     // extract the spectra
1035     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
1036     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1037     if(!fJetPtDeltaPhi) {
1038         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1039         return kFALSE;
1040     }
1041     if(fCentralityArray) {
1042         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1043             spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
1044             fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1045         }
1046     }
1047     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1048     // in plane spectrum
1049     fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1050     // extract the delta pt matrices
1051     TString deltaptName("");
1052     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1053     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1054     if(!fDeltaPtDeltaPhi) {
1055         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1056         printf(" > may be ok, depending no what you want to do < \n");
1057         fRefreshInput = kTRUE;
1058         return kTRUE;
1059     }
1060     if(fCentralityArray) {
1061         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1062             deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1063             fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1064         }
1065     }
1066
1067     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1068     // in plane delta pt distribution
1069     fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1070     // create a rec - true smeared response matrix
1071     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1072     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1073         Bool_t skip = kFALSE;
1074         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1075             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1076             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1077         }
1078     }
1079     fDptIn = new TH2D(*rfIn);
1080     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1081     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1082     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1083     fDptIn = ProtectHeap(fDptIn);
1084     
1085     return kTRUE;
1086 }
1087 //_____________________________________________________________________________
1088 TH1D* AliJetFlowTools::GetPrior(
1089         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
1090         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
1091         const TH1D* kinematicEfficiency,              // kinematic efficiency
1092         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
1093         const TString suffix,                         // suffix (in, out)
1094         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
1095 {
1096     // 1) get a prior for unfolding. 
1097     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1098     TH1D* unfolded(0x0);
1099     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1100     dirOut->cd();
1101     switch (fPrior) {    // select the prior for unfolding
1102         case kPriorPythia : {
1103             if(!fPriorUser) {
1104                 printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1105                 return 0x0;
1106             }
1107             // rebin the given prior to the true spectrum (creates a new histo)
1108             return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1109         } break;
1110         case kPriorChi2 : {
1111             TArrayI* placeHolder(0x0);
1112             if(fMergeBinsArray) {
1113                 placeHolder = fMergeBinsArray;
1114                 fMergeBinsArray = 0x0;
1115             }
1116             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
1117                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original binning
1118                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
1119                 TArrayD* tempArrayRec(fBinsRec);
1120                 fBinsRec = fBinsRecPrior;
1121                 // for the prior, do not re-bin the output
1122                 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1123                 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1124                 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1125                 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1126                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1127                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1128                 unfolded= UnfoldSpectrumChi2(
1129                             measuredJetSpectrumChi2,
1130                             resizedResponseChi2,
1131                             kinematicEfficiencyChi2,
1132                             measuredJetSpectrumTrueBinsChi2,    // prior for chi2 unfolding (measured)
1133                             TString(Form("prior_%s", suffix.Data())));
1134                if(!unfolded) {
1135                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1136                     printf("               probably Chi2 unfolding did not converge < \n");
1137                     if(placeHolder) fMergeBinsArray = placeHolder;
1138                     return 0x0;
1139                 }
1140                 fBinsTrue = tempArrayTrue;  // reset bins borders
1141                 fBinsRec = tempArrayRec;
1142                 // if the chi2 prior has a different binning, rebin to the true binning for the  unfolding
1143                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
1144             } else {
1145                 unfolded = UnfoldSpectrumChi2(
1146                             measuredJetSpectrum,
1147                             resizedResponse,
1148                             kinematicEfficiency,
1149                             measuredJetSpectrumTrueBins,        // prior for chi2 unfolding (measured)
1150                             TString(Form("prior_%s", suffix.Data())));
1151                 if(!unfolded) {
1152                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1153                     printf("               probably Chi2 unfolding did not converge < \n");
1154                     if(placeHolder) fMergeBinsArray = placeHolder;
1155                     return 0x0;
1156                 }
1157                 if(placeHolder) fMergeBinsArray = placeHolder;
1158             }
1159             break;
1160
1161         }
1162         case kPriorMeasured : { 
1163             unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
1164         }
1165         default : break;
1166     }
1167     // it can be important that the prior is smooth, this can be achieved by 
1168     // extrapolating the spectrum with a fitted power law when bins are sparsely filed 
1169     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1170     TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1171     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1172     return priorLocal;
1173 }
1174 //_____________________________________________________________________________
1175 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1176     // resize the x-axis of a th1d
1177     if(!histo) {
1178         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1179         return NULL;
1180     } 
1181     // see how many bins we need to copy
1182     TH1D* resized = new TH1D(Form("%s_resized_%s", histo->GetName(), suffix.Data()), Form("%s_resized_%s", histo->GetName(), suffix.Data()), up-low, (double)low, (double)up);
1183     // low is the bin number of the first new bin
1184     Int_t l = histo->GetXaxis()->FindBin(low);
1185     // set the values
1186     Double_t _x(0), _xx(0);
1187     for(Int_t i(0); i < up-low; i++) {
1188         _x = histo->GetBinContent(l+i);
1189         _xx=histo->GetBinError(l+i);
1190         resized->SetBinContent(i+1, _x);
1191         resized->SetBinError(i+1, _xx);
1192     }
1193     return resized;
1194 }
1195 //_____________________________________________________________________________
1196 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1197     // resize the y-axis of a th2d
1198     if(!histo) {
1199         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1200         return NULL;
1201     } 
1202     // see how many bins we need to copy
1203     TH2D* resized = new TH2D(Form("%s_resized_%s", histo->GetName(), suffix.Data()), Form("%s_resized_%s", histo->GetName(), suffix.Data()), x->GetSize()-1, x->GetArray(), y->GetSize()-1, y->GetArray());
1204     // assume only the y-axis has changed
1205     // low is the bin number of the first new bin
1206     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1207     // set the values
1208     Double_t _x(0), _xx(0);
1209     for(Int_t i(0); i < x->GetSize(); i++) {
1210         for(Int_t j(0); j < y->GetSize(); j++) {
1211             _x = histo->GetBinContent(i, low+j);
1212             _xx=histo->GetBinError(i, low+1+j);
1213             resized->SetBinContent(i, j, _x);
1214             resized->SetBinError(i, j, _xx);
1215         }
1216     }
1217     return resized;
1218 }
1219 //_____________________________________________________________________________
1220 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1221     // general method to normalize all vertical slices of a th2 to unity
1222     // i.e. get a probability matrix
1223     if(!histo) {
1224         printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
1225         return NULL;
1226     }
1227     Int_t binsX = histo->GetXaxis()->GetNbins();
1228     Int_t binsY = histo->GetYaxis()->GetNbins();
1229     
1230     // normalize all slices in x
1231     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1232         Double_t weight = 0;
1233         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1234             weight+=histo->GetBinContent(i+1, j+1);
1235         }       // now we know the total weight
1236         for(Int_t j(0); j < binsY; j++) {
1237             if (weight <= 0 ) continue;
1238             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1239             if(noError) histo->SetBinError(  1+i, j+1, 0.);
1240             else histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1241         }
1242     }
1243     return histo;
1244 }
1245 //_____________________________________________________________________________
1246 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1247     // return a TH1D with the supplied histogram rebinned to the supplied bins
1248     // the returned histogram is new, the original is deleted from the heap if kill is true
1249     if(!histo || !bins) {
1250         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1251         return NULL;
1252     }
1253     // create the output histo
1254     TString name = histo->GetName();
1255     name+="_template";
1256     name+=suffix;
1257     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1258     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1259         // loop over the bins of the old histo and fill the new one with its data
1260         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1261     }
1262     if(kill) delete histo;
1263     return rebinned;
1264 }
1265 //_____________________________________________________________________________
1266 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1267     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1268     if(!fResponseMaker || !binsTrue || !binsRec) {
1269         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1270         return 0x0;
1271     }
1272     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1273     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1274 }
1275 //_____________________________________________________________________________
1276 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1277 {
1278     // multiply two matrices
1279     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1280     TH2D* c = (TH2D*)a->Clone("c");
1281     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1282         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1283             Double_t val = 0;
1284             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1285                 Int_t y2 = x1;
1286                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1287             }
1288             c->SetBinContent(x2, y1, val);
1289             c->SetBinError(x2, y1, 0.);
1290         }
1291     }
1292     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1293     return c;
1294 }
1295 //_____________________________________________________________________________
1296 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale) 
1297 {
1298     // normalize a th1d to a certain scale
1299     histo->Sumw2();
1300     Double_t integral = histo->Integral()*scale;
1301     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1302     else if (scale != 1.) histo->Scale(1./scale, "width");
1303     else printf(" > Histogram integral < 0, cannot normalize \n");
1304     return histo;
1305 }
1306 //_____________________________________________________________________________
1307 TH1D* AliJetFlowTools::MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr)
1308 {
1309     // merge a spectrum histogram taking into account the correlation terms
1310     if(!(bins&&spectrum)) {
1311         printf(" > NULL pointer passed as argument in MergeSpectrumBins ! < \n");
1312         return 0x0;
1313     }
1314     Double_t sum(0), error(0), pearson(0);
1315     // take the sum of the bin content
1316     sum += spectrum->GetBinContent(bins->At(0));
1317     sum += spectrum->GetBinContent(bins->At(1));
1318     // quadratically sum the errors
1319     error += TMath::Power(spectrum->GetBinError(bins->At(0)), 2);
1320     error += TMath::Power(spectrum->GetBinError(bins->At(1)), 2);
1321     // add the covariance term
1322     pearson = corr->GetBinContent(bins->At(0), bins->At(1));
1323     if(!corr) {
1324         printf(" > PANIC ! something is wrong with the covariance matrix, assuming full correlation ! < \n ");
1325         pearson = 1;
1326     }
1327     error += 2.*spectrum->GetBinError(bins->At(0))*spectrum->GetBinError(bins->At(1))*pearson;
1328     spectrum->SetBinContent(1, sum);
1329     if(error <= 0) return spectrum;
1330     spectrum->SetBinError(1, TMath::Sqrt(error));
1331     printf(" > sum is %.2f \t error is %.8f < \n", sum, error);
1332     printf(" > REPLACING BIN CONTENT OF FIRST BIN (%i) WITH MERGED BINS (%i) and (%i) < \n", 1, bins->At(0), bins->At(1));
1333     return spectrum;
1334 }
1335 //_____________________________________________________________________________
1336 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix) 
1337 {
1338     // Calculate pearson coefficients from covariance matrix
1339     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1340     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1341     Double_t pearson(0.);
1342     if(nrows==0 && ncols==0) return 0x0;
1343     for(Int_t row = 0; row < nrows; row++) {
1344         for(Int_t col = 0; col<ncols; col++) {
1345         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1346         (*pearsonCoefficients)(row,col) = pearson;
1347         }
1348     }
1349     return pearsonCoefficients;
1350 }
1351 //_____________________________________________________________________________
1352 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1353     // smoothen the spectrum using a user defined function
1354     // returns a clone of the original spectrum if fitting failed
1355     // if kill is kTRUE the input spectrum will be deleted from the heap
1356     // if 'count' is selected, bins are filled with integers (necessary if the 
1357     // histogram is interpreted in a routine which accepts only counts)
1358     if(!spectrum || !function) return 0x0;
1359     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1360     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1361     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1362     TFitResultPtr r = temp->Fit(function, "", "", min, max);
1363     if((int)r == 0) {   // MINUIT status
1364         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1365             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1366                 (counts) ? temp->SetBinContent(i, (int)(function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i))) : temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1367                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1368             }
1369         }
1370     }
1371     if(kill) delete spectrum;
1372     return temp;
1373 }
1374 //_____________________________________________________________________________
1375 void AliJetFlowTools::Style(Bool_t legacy)
1376 {
1377     // set global style for your current aliroot session
1378     if(!gStyle) return;
1379     // legacy style is pleasing to the eye, default is the formal ALICE style
1380     if(legacy) {
1381         gStyle->SetCanvasColor(-1); 
1382         gStyle->SetPadColor(-1); 
1383         gStyle->SetFrameFillColor(-1); 
1384         gStyle->SetHistFillColor(-1); 
1385         gStyle->SetTitleFillColor(-1); 
1386         gStyle->SetFillColor(-1); 
1387         gStyle->SetFillStyle(4000); 
1388         gStyle->SetStatStyle(0); 
1389         gStyle->SetTitleStyle(0); 
1390         gStyle->SetCanvasBorderSize(0); 
1391         gStyle->SetFrameBorderSize(0); 
1392         gStyle->SetLegendBorderSize(0); 
1393         gStyle->SetStatBorderSize(0); 
1394         gStyle->SetTitleBorderSize(0);
1395     } else {
1396         gStyle->Reset("Plain");
1397         gStyle->SetOptTitle(0);
1398         gStyle->SetOptStat(0);
1399         gStyle->SetPalette(1);
1400         gStyle->SetCanvasColor(10);
1401         gStyle->SetCanvasBorderMode(0);
1402         gStyle->SetFrameLineWidth(1);
1403         gStyle->SetFrameFillColor(kWhite);
1404         gStyle->SetPadColor(10);
1405         gStyle->SetPadTickX(1);
1406         gStyle->SetPadTickY(1);
1407         gStyle->SetPadBottomMargin(0.15);
1408         gStyle->SetPadLeftMargin(0.15);
1409         gStyle->SetHistLineWidth(1);
1410         gStyle->SetHistLineColor(kRed);
1411         gStyle->SetFuncWidth(2);
1412         gStyle->SetFuncColor(kGreen);
1413         gStyle->SetLineWidth(2);
1414         gStyle->SetLabelSize(0.045,"xyz");
1415         gStyle->SetLabelOffset(0.01,"y");
1416         gStyle->SetLabelOffset(0.01,"x");
1417         gStyle->SetLabelColor(kBlack,"xyz");
1418         gStyle->SetTitleSize(0.05,"xyz");
1419         gStyle->SetTitleOffset(1.25,"y");
1420         gStyle->SetTitleOffset(1.2,"x");
1421         gStyle->SetTitleFillColor(kWhite);
1422         gStyle->SetTextSizePixels(26);
1423         gStyle->SetTextFont(42);
1424         gStyle->SetLegendBorderSize(0);
1425         gStyle->SetLegendFillColor(kWhite);
1426         gStyle->SetLegendFont(42);
1427     }
1428 }
1429 //_____________________________________________________________________________
1430 void AliJetFlowTools::Style(TCanvas* c, TString style)
1431 {
1432     // set a default style for a canvas
1433     if(!strcmp(style.Data(), "PEARSON")) {
1434         printf(" > style PEARSON canvas < \n");
1435         gStyle->SetOptStat(0);
1436         c->SetGridx();
1437         c->SetGridy();
1438         c->SetTicks();
1439         return;
1440     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1441         printf(" > style SPECTRUM canvas < \n");
1442         gStyle->SetOptStat(0);
1443         c->SetLogy();
1444         c->SetGridx();
1445         c->SetGridy();
1446         c->SetTicks();
1447         return;
1448     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1449 }
1450 //_____________________________________________________________________________
1451 void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
1452 {
1453     // set a default style for a canva
1454     
1455     if(legacy) {
1456         c->SetLeftMargin(.25);
1457         c->SetBottomMargin(.25);
1458     }
1459     else Style();
1460     if(!strcmp(style.Data(), "PEARSON")) {
1461         printf(" > style PEARSON pad < \n");
1462         gStyle->SetOptStat(0);
1463         c->SetGridx();
1464         c->SetGridy();
1465         c->SetTicks();
1466         return;
1467     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1468         printf(" > style SPECTRUM pad < \n");
1469         gStyle->SetOptStat(0);
1470         c->SetLogy();
1471         c->SetGridx();
1472         c->SetGridy();
1473         c->SetTicks();
1474         return;
1475     } else if (!strcmp(style.Data(), "GRID")) {
1476         printf(" > style GRID pad < \n");
1477         gStyle->SetOptStat(0);
1478         c->SetGridx();
1479         c->SetGridy();
1480         c->SetTicks();
1481     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1482 }
1483 //_____________________________________________________________________________
1484 void AliJetFlowTools::Style(TLegend* l)
1485 {
1486     // set a default style for a legend
1487     l->SetFillColor(0);
1488     l->SetBorderSize(0);
1489     if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
1490 }
1491 //_____________________________________________________________________________
1492 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
1493 {
1494     // style a histo
1495     h->GetYaxis()->SetNdivisions(505);
1496     h->SetLineColor(col);
1497     h->SetMarkerColor(col);
1498     h->SetLineWidth(2);
1499     h->SetMarkerSize(1);
1500     if(legacy) {
1501         h->SetTitle("");
1502         h->GetYaxis()->SetLabelSize(0.05);
1503         h->GetXaxis()->SetLabelSize(0.05);
1504         h->GetYaxis()->SetTitleOffset(1.5);
1505         h->GetXaxis()->SetTitleOffset(1.5);
1506         h->GetYaxis()->SetTitleSize(.05);
1507         h->GetXaxis()->SetTitleSize(.05);
1508     } else Style();
1509     switch (type) {
1510         case kInPlaneSpectrum : {
1511             h->SetTitle("IN PLANE");
1512             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1513             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1514         } break;
1515         case kOutPlaneSpectrum : {
1516             h->SetTitle("OUT OF PLANE");
1517             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1518             h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1519        } break;
1520        case kUnfoldedSpectrum : {
1521             h->SetTitle("UNFOLDED");
1522             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1523             h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1524        } break;
1525        case kFoldedSpectrum : {
1526             h->SetTitle("FOLDED");
1527             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1528             h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1529        } break;
1530        case kMeasuredSpectrum : {
1531             h->SetTitle("MEASURED");
1532             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1533             h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{cht} (GeV/#it{c})");
1534        } break;
1535        case kBar : {
1536             h->SetFillColor(col);
1537             h->SetBarWidth(.6);
1538             h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1539             h->SetBarOffset(0.2);
1540        }
1541        case kRatio : {
1542             h->SetMarkerStyle(8);
1543             h->SetMarkerSize(1);
1544        } break;
1545        default : break;
1546     }
1547 }
1548 //_____________________________________________________________________________
1549 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
1550 {
1551     // style a tgraph
1552     h->GetYaxis()->SetNdivisions(505);
1553     h->SetLineColor(col);
1554     h->SetMarkerColor(col);
1555     h->SetLineWidth(2);
1556     h->SetMarkerSize(1);
1557     h->SetTitle("");
1558     h->SetFillColor(kCyan);
1559     if(legacy) {
1560         h->GetYaxis()->SetLabelSize(0.05);
1561         h->GetXaxis()->SetLabelSize(0.05);
1562         h->GetYaxis()->SetTitleOffset(1.6);
1563         h->GetXaxis()->SetTitleOffset(1.6);
1564         h->GetYaxis()->SetTitleSize(.05);
1565         h->GetXaxis()->SetTitleSize(.05);
1566     } else Style();
1567     h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1568     switch (type) {
1569         case kInPlaneSpectrum : {
1570             h->SetTitle("IN PLANE");
1571             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1572         } break;
1573         case kOutPlaneSpectrum : {
1574             h->SetTitle("OUT OF PLANE");
1575             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1576        } break;
1577        case kUnfoldedSpectrum : {
1578             h->SetTitle("UNFOLDED");
1579             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1580        } break;
1581        case kFoldedSpectrum : {
1582             h->SetTitle("FOLDED");
1583             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1584        } break;
1585        case kMeasuredSpectrum : {
1586             h->SetTitle("MEASURED");
1587             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1588        } break;
1589        case kRatio : {
1590 //            h->GetYaxis()->SetTitle("#frac{d#it{N_{in plane}^{jet}}}{d#it{p}_{T}} / #frac{d#it{N_{out of plane}^{jet}}}{d#it{p}_{T}}");
1591             h->GetYaxis()->SetTitle("(d#it{N}^{ch, jet}_{in plane}/(d#it{p}_{T}d#eta))/(d#it{N}^{ch,jet}_{out of plane}/(d#it{p}_{T}d#eta))");
1592        } break;
1593        case kV2 : {
1594 //            h->GetYaxis()->SetTitle("#it{v}_{2} = #frac{1}{#it{R}} #frac{#pi}{4} #frac{#it{N_{in plane}} - #it{N_{out of plane}}}{#it{N_{in plane}} + #it{N_{out of plane}}}");
1595             h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet} \{EP, |#Delta#eta|>0.9 \} ");
1596             h->GetYaxis()->SetRangeUser(-.5, 1.);
1597             h->SetMarkerStyle(8);
1598             h->SetMarkerSize(1);
1599        } break;
1600        default : break;
1601     }
1602 }
1603 //_____________________________________________________________________________
1604 void AliJetFlowTools::GetNominalValues(
1605         TH1D*& ratio,           // pointer reference, output of this function
1606         TGraphErrors*& v2,      // pointer reference, as output of this function
1607         TArrayI* in,
1608         TArrayI* out,
1609         TString inFile,
1610         TString outFile) const
1611 {
1612     // pass clones of the nominal points and nominal v2 values 
1613     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1614     TFile* readMe(new TFile(inFile.Data(), "READ"));                    // open unfolding output read-only
1615     if(readMe->IsZombie()) {
1616         printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1617         return;
1618     }
1619     printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1620     readMe->ls();
1621     TFile* output(new TFile(outFile.Data(), "RECREATE"));   // create new output file
1622     // get some placeholders, have to be initialized but will be deleted
1623     ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1624     TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1625     TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1626     Int_t iIn[] = {in->At(0), in->At(0)};          // first value is the nominal point
1627     Int_t iOut[] = {out->At(0), out->At(0)};
1628
1629     // call the functions and set the relevant pointer references
1630     TH1D* dud(0x0);
1631     DoIntermediateSystematics(
1632             new TArrayI(2, iIn), 
1633             new TArrayI(2, iOut),
1634             dud, dud, dud, dud, dud, dud, 
1635             ratio,              // pointer reference, output of this function 
1636             nominalIn,
1637             nominalOut,
1638             1, 
1639             fBinsTrue->At(0), 
1640             fBinsTrue->At(fBinsTrue->GetSize()-1),
1641             readMe,
1642             "nominal_values");
1643     v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1644
1645     // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1646     ratio->SetDirectory(0);     // disassociate from current gDirectory
1647     readMe->Close();
1648 }
1649 //_____________________________________________________________________________
1650 void AliJetFlowTools::GetCorrelatedUncertainty(
1651         TGraphAsymmErrors*& corrRatio,  // correlated uncertainty function pointer
1652         TGraphAsymmErrors*& corrV2,     // correlated uncertainty function pointer
1653         TArrayI* variationsIn,          // variantions in plnae
1654         TArrayI* variationsOut,         // variantions out of plane
1655         Bool_t sym,                     // treat as symmmetric
1656         TArrayI* variations2ndIn,       // second source of variations
1657         TArrayI* variations2ndOut,      // second source of variations
1658         Bool_t sym2nd,                  // treat as symmetric
1659         TString type,                   // type of uncertaitny
1660         TString type2,
1661         Int_t columns,                  // divide the output canvasses in this many columns
1662         Float_t rangeLow,               // lower pt range
1663         Float_t rangeUp,                // upper pt range
1664         Float_t corr,                   // correlation strength
1665         TString in,                     // input file name (created by this unfolding class)
1666         TString out                     // output file name (which will hold results of the systematic test)
1667         ) const
1668 {
1669     // do full systematics
1670     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1671     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1672     if(readMe->IsZombie()) {
1673         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1674         return;
1675     }
1676     printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1677     readMe->ls();
1678     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1679
1680     // create some null placeholder pointers
1681     TH1D* relativeErrorVariationInLow(0x0);
1682     TH1D* relativeErrorVariationInUp(0x0);
1683     TH1D* relativeErrorVariationOutLow(0x0);
1684     TH1D* relativeErrorVariationOutUp(0x0);
1685     TH1D* relativeError2ndVariationInLow(0x0);
1686     TH1D* relativeError2ndVariationInUp(0x0);
1687     TH1D* relativeError2ndVariationOutLow(0x0);
1688     TH1D* relativeError2ndVariationOutUp(0x0);
1689     TH1D* relativeStatisticalErrorIn(0x0);
1690     TH1D* relativeStatisticalErrorOut(0x0);
1691     // histo for the nominal ratio in / out
1692     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1693     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1694     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1695
1696     // call the functions
1697     if(variationsIn && variationsOut) {
1698         DoIntermediateSystematics(
1699                 variationsIn, 
1700                 variationsOut, 
1701                 relativeErrorVariationInUp,        // pointer reference
1702                 relativeErrorVariationInLow,       // pointer reference
1703                 relativeErrorVariationOutUp,       // pointer reference
1704                 relativeErrorVariationOutLow,      // pointer reference
1705                 relativeStatisticalErrorIn,        // pointer reference
1706                 relativeStatisticalErrorOut,       // pointer reference
1707                 nominal,
1708                 nominalIn,
1709                 nominalOut,
1710                 columns, 
1711                 rangeLow, 
1712                 rangeUp,
1713                 readMe,
1714                 type);
1715         if(relativeErrorVariationInUp) {
1716             // canvas with the error from variation strength
1717             TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1718             relativeErrorVariation->Divide(2);
1719             relativeErrorVariation->cd(1);
1720             Style(gPad, "GRID");
1721             relativeErrorVariationInUp->DrawCopy("b");
1722             relativeErrorVariationInLow->DrawCopy("same b");
1723             Style(AddLegend(gPad));
1724             relativeErrorVariation->cd(2);
1725             Style(gPad, "GRID");
1726             relativeErrorVariationOutUp->DrawCopy("b");
1727             relativeErrorVariationOutLow->DrawCopy("same b");
1728             SavePadToPDF(relativeErrorVariation);
1729             Style(AddLegend(gPad));
1730             relativeErrorVariation->Write();
1731
1732             // now smoothen the diced response error (as it is expected to be flat)
1733             // this is done by fitting a constant to the diced resonse error histo
1734             //
1735             TF1* lin = new TF1("lin", "[0]", rangeLow, rangeUp);
1736             relativeErrorVariationInUp->Fit(lin, "L", "", rangeLow, rangeUp);
1737             if(!gMinuit->fISW[1] == 3) printf(" fit is NOT ok ! " );
1738             for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1739                 relativeErrorVariationInUp->SetBinContent(i+1, lin->GetParameter(0));
1740             }
1741             relativeErrorVariationInLow->Fit(lin, "L", "", rangeLow, rangeUp);
1742             printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1743             for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1744                 relativeErrorVariationInLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1745             }
1746             relativeErrorVariationOutUp->Fit(lin, "L", "", rangeLow, rangeUp);
1747             printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1748             for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1749                 relativeErrorVariationOutUp->SetBinContent(i+1, lin->GetParameter(0));
1750             }
1751             relativeErrorVariationOutLow->Fit(lin, "L", "", rangeLow, rangeUp);
1752             printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1753             for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1754                 relativeErrorVariationOutLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1755             }
1756
1757         
1758         
1759         }
1760     }
1761     // call the functions for a second set of systematic sources
1762     if(variations2ndIn && variations2ndOut) {
1763         DoIntermediateSystematics(
1764                 variations2ndIn, 
1765                 variations2ndOut, 
1766                 relativeError2ndVariationInUp,        // pointer reference
1767                 relativeError2ndVariationInLow,       // pointer reference
1768                 relativeError2ndVariationOutUp,       // pointer reference
1769                 relativeError2ndVariationOutLow,      // pointer reference
1770                 relativeStatisticalErrorIn,           // pointer reference
1771                 relativeStatisticalErrorOut,          // pointer reference
1772                 nominal,
1773                 nominalIn,
1774                 nominalOut,
1775                 columns, 
1776                 rangeLow, 
1777                 rangeUp,
1778                 readMe,
1779                 type2);
1780         if(relativeError2ndVariationInUp) {
1781             // canvas with the error from variation strength
1782             TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
1783             relativeError2ndVariation->Divide(2);
1784             relativeError2ndVariation->cd(1);
1785             Style(gPad, "GRID");
1786             relativeError2ndVariationInUp->DrawCopy("b");
1787             relativeError2ndVariationInLow->DrawCopy("same b");
1788             Style(AddLegend(gPad));
1789             relativeError2ndVariation->cd(2);
1790             Style(gPad, "GRID");
1791             relativeError2ndVariationOutUp->DrawCopy("b");
1792             relativeError2ndVariationOutLow->DrawCopy("same b");
1793             SavePadToPDF(relativeError2ndVariation);
1794             Style(AddLegend(gPad));
1795             relativeError2ndVariation->Write();
1796         }
1797
1798     }
1799
1800     // and the placeholder for the final systematic
1801     TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1802     TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1803     TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1804     TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1805     relativeErrorInUp->SetYTitle("relative uncertainty");
1806     relativeErrorOutUp->SetYTitle("relative uncertainty");
1807     relativeErrorInLow->SetYTitle("relative uncertainty");
1808     relativeErrorOutLow->SetYTitle("relative uncertainty");
1809
1810     // merge the correlated errors (FIXME) trivial for one set 
1811     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1812     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1813     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1814     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1815
1816     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1817         // for the upper bound
1818         if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
1819         if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
1820         if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1821         if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
1822         dInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1823         // for a symmetric first variation
1824         if(sym) dInUp += aInLow*aInLow;
1825         if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1826         dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1827         if(sym) dOutUp += aOutLow*aOutLow;
1828         if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1829         // for the lower bound
1830         if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
1831         if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
1832         if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1833         if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
1834         dInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1835         if(sym) dInLow += aInUp*aInUp;
1836         if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
1837         dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1838         if(sym) dOutLow += aOutUp*aOutUp;
1839         if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1840     }
1841     // project the estimated errors on the nominal ratio
1842     if(nominal) {
1843         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1844         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1845         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1846         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1847         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1848         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1849         Double_t lowErr(0.), upErr(0.);
1850         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1851             // add the in and out of plane errors in quadrature
1852             // propagation is tricky because we should consider two cases:
1853             // [1] the error is uncorrelated, and we propagate upper 1 with lower 2 and vice versa
1854             // [2] the error is correlated - in this case upper 1 should be propagated with upper 2
1855             // as the fluctuations are bound to always to of same sign
1856             if(corr <= 0) {
1857                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
1858                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1859             } else {
1860                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) -2.*corr*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1861                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 2.*corr*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1862             }
1863             ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
1864             ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
1865             // get the bin width (which is the 'error' on x
1866             Double_t binWidth(nominal->GetBinWidth(i+1));
1867             axl[i] = binWidth/2.;
1868             axh[i] = binWidth/2.;
1869             // now get the coordinate for the point
1870             ax[i] = nominal->GetBinCenter(i+1);
1871             ay[i] = nominal->GetBinContent(i+1);
1872         }
1873         // save the nominal ratio
1874         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
1875         corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
1876         nominalError->SetName("correlated uncertainty");
1877         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
1878         nominalCanvas->Divide(2);
1879         nominalCanvas->cd(1);
1880         Style(nominal, kBlack);
1881         Style(nominalError, kCyan, kRatio);
1882         nominalError->SetLineColor(kCyan);
1883         nominalError->SetMarkerColor(kCyan);
1884         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1885         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
1886         nominalError->DrawClone("a2");
1887         nominal->DrawCopy("same E1");
1888         Style(AddLegend(gPad));
1889         Style(gPad, "GRID");
1890         Style(nominalCanvas);
1891         // save nominal v2 and systematic errors
1892         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
1893                     nominalIn,
1894                     nominalOut,
1895                     fEventPlaneRes,
1896                     "v_{2} with shape uncertainty",
1897                     relativeErrorInUp,
1898                     relativeErrorInLow,
1899                     relativeErrorOutUp,
1900                     relativeErrorOutLow,
1901                     corr));
1902         // pass the nominal values to the pointer references
1903         corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
1904         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
1905         nominalCanvas->cd(2);
1906         Style(nominalV2, kBlack);
1907         Style(nominalV2Error, kCyan, kV2);
1908         nominalV2Error->SetLineColor(kCyan);
1909         nominalV2Error->SetMarkerColor(kCyan);
1910         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1911         nominalV2Error->DrawClone("a2");
1912         nominalV2->DrawClone("same E1");
1913         Style(AddLegend(gPad));
1914         Style(gPad, "GRID");
1915         Style(nominalCanvas);
1916         SavePadToPDF(nominalCanvas);
1917         nominalCanvas->Write();
1918     }
1919
1920     TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
1921     relativeError->Divide(2);
1922     relativeError->cd(1);
1923     Style(gPad, "GRID");
1924     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1925     Style(relativeErrorInUp, kBlue, kBar);
1926     Style(relativeErrorInLow, kGreen, kBar);
1927     relativeErrorInUp->DrawCopy("b");
1928     relativeErrorInLow->DrawCopy("same b");
1929     Style(relativeStatisticalErrorIn, kRed);
1930     relativeStatisticalErrorIn->DrawCopy("same");
1931     Style(AddLegend(gPad));
1932     relativeError->cd(2);
1933     Style(gPad, "GRID");
1934     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1935     Style(relativeErrorOutUp, kBlue, kBar);
1936     Style(relativeErrorOutLow, kGreen, kBar);
1937     relativeErrorOutUp->DrawCopy("b");
1938     relativeErrorOutLow->DrawCopy("same b");
1939     Style(relativeStatisticalErrorOut, kRed);
1940     relativeStatisticalErrorOut->DrawCopy("same");
1941     Style(AddLegend(gPad));
1942
1943     // write the buffered file to disk and close the file
1944     SavePadToPDF(relativeError);
1945     relativeError->Write();
1946     output->Write();
1947 //    output->Close();
1948 }
1949 //_____________________________________________________________________________
1950 void AliJetFlowTools::GetShapeUncertainty(
1951         TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
1952         TGraphAsymmErrors*& shapeV2,    // pointer reference to final shape uncertainty of v2
1953         TArrayI* regularizationIn,      // regularization strength systematics, in plane
1954         TArrayI* regularizationOut,     // regularization strenght systematics, out of plane
1955         TArrayI* trueBinIn,             // true bin ranges, in plane
1956         TArrayI* trueBinOut,            // true bin ranges, out of plane
1957         TArrayI* recBinIn,              // rec bin ranges, in plane
1958         TArrayI* recBinOut,             // rec bin ranges, out of plane
1959         TArrayI* methodIn,              // method in
1960         TArrayI* methodOut,             // method out
1961         Int_t columns,                  // divide the output canvasses in this many columns
1962         Float_t rangeLow,               // lower pt range
1963         Float_t rangeUp,                // upper pt range
1964         TString in,                     // input file name (created by this unfolding class)
1965         TString out                     // output file name (which will hold results of the systematic test)
1966         ) const
1967 {
1968     // do full systematics
1969     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1970     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1971     if(readMe->IsZombie()) {
1972         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1973         return;
1974     }
1975     printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1976     readMe->ls();
1977     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1978
1979     // create some null placeholder pointers
1980     TH1D* relativeErrorRegularizationInLow(0x0);
1981     TH1D* relativeErrorRegularizationInUp(0x0);
1982     TH1D* relativeErrorTrueBinInLow(0x0);
1983     TH1D* relativeErrorTrueBinInUp(0x0);
1984     TH1D* relativeErrorRecBinInLow(0x0);
1985     TH1D* relativeErrorRecBinInUp(0x0);
1986     TH1D* relativeErrorMethodInLow(0x0);
1987     TH1D* relativeErrorMethodInUp(0x0);
1988     TH1D* relativeErrorRegularizationOutLow(0x0);
1989     TH1D* relativeErrorRegularizationOutUp(0x0);
1990     TH1D* relativeErrorTrueBinOutLow(0x0);
1991     TH1D* relativeErrorTrueBinOutUp(0x0);
1992     TH1D* relativeErrorRecBinOutLow(0x0);
1993     TH1D* relativeErrorRecBinOutUp(0x0);
1994     TH1D* relativeStatisticalErrorIn(0x0);
1995     TH1D* relativeStatisticalErrorOut(0x0);
1996     TH1D* relativeErrorMethodOutLow(0x0);
1997     TH1D* relativeErrorMethodOutUp(0x0);
1998     // histo for the nominal ratio in / out
1999     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2000     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2001     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2002
2003     // call the functions
2004     if(regularizationIn && regularizationOut) {
2005         DoIntermediateSystematics(
2006                 regularizationIn, 
2007                 regularizationOut, 
2008                 relativeErrorRegularizationInUp,        // pointer reference
2009                 relativeErrorRegularizationInLow,       // pointer reference
2010                 relativeErrorRegularizationOutUp,       // pointer reference
2011                 relativeErrorRegularizationOutLow,      // pointer reference
2012                 relativeStatisticalErrorIn,             // pointer reference
2013                 relativeStatisticalErrorOut,            // pointer reference
2014                 nominal,
2015                 nominalIn,
2016                 nominalOut,
2017                 columns, 
2018                 rangeLow, 
2019                 rangeUp,
2020                 readMe,
2021                 "regularization",
2022                 fRMS);
2023         if(relativeErrorRegularizationInUp) {
2024             // canvas with the error from regularization strength
2025             TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
2026             relativeErrorRegularization->Divide(2);
2027             relativeErrorRegularization->cd(1);
2028             Style(gPad, "GRID");
2029             relativeErrorRegularizationInUp->DrawCopy("b");
2030             relativeErrorRegularizationInLow->DrawCopy("same b");
2031             Style(AddLegend(gPad));
2032             relativeErrorRegularization->cd(2);
2033             Style(gPad, "GRID");
2034             relativeErrorRegularizationOutUp->DrawCopy("b");
2035             relativeErrorRegularizationOutLow->DrawCopy("same b");
2036             SavePadToPDF(relativeErrorRegularization);
2037             Style(AddLegend(gPad));
2038             relativeErrorRegularization->Write();
2039         }
2040     }
2041     if(trueBinIn && trueBinOut) {
2042         DoIntermediateSystematics(
2043                 trueBinIn, 
2044                 trueBinOut, 
2045                 relativeErrorTrueBinInUp,       // pointer reference
2046                 relativeErrorTrueBinInLow,      // pointer reference
2047                 relativeErrorTrueBinOutUp,      // pointer reference
2048                 relativeErrorTrueBinOutLow,     // pointer reference
2049                 relativeStatisticalErrorIn,
2050                 relativeStatisticalErrorOut,
2051                 nominal,
2052                 nominalIn,
2053                 nominalOut,
2054                 columns, 
2055                 rangeLow, 
2056                 rangeUp,
2057                 readMe,
2058                 "trueBin");
2059         if(relativeErrorTrueBinInUp) {
2060             TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
2061             relativeErrorTrueBin->Divide(2);
2062             relativeErrorTrueBin->cd(1);
2063             Style(gPad, "GRID");
2064             relativeErrorTrueBinInUp->DrawCopy("b");
2065             relativeErrorTrueBinInLow->DrawCopy("same b");
2066             Style(AddLegend(gPad));
2067             relativeErrorTrueBin->cd(2);
2068             Style(gPad, "GRID");
2069             relativeErrorTrueBinOutUp->DrawCopy("b");
2070             relativeErrorTrueBinOutLow->DrawCopy("same b");
2071             SavePadToPDF(relativeErrorTrueBin);
2072             Style(AddLegend(gPad));
2073             relativeErrorTrueBin->Write();
2074         }
2075     }
2076     if(recBinIn && recBinOut) {
2077         DoIntermediateSystematics(
2078                 recBinIn, 
2079                 recBinOut, 
2080                 relativeErrorRecBinInUp,       // pointer reference
2081                 relativeErrorRecBinInLow,        // pointer reference
2082                 relativeErrorRecBinOutUp,      // pointer reference
2083                 relativeErrorRecBinOutLow,       // pointer reference
2084                 relativeStatisticalErrorIn,
2085                 relativeStatisticalErrorOut,
2086                 nominal,
2087                 nominalIn,
2088                 nominalOut,
2089                 columns, 
2090                 rangeLow, 
2091                 rangeUp,
2092                 readMe,
2093                 "recBin",
2094                 fRMS);
2095         if(relativeErrorRecBinOutUp) {
2096             // canvas with the error from regularization strength
2097             TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
2098             relativeErrorRecBin->Divide(2);
2099             relativeErrorRecBin->cd(1);
2100             Style(gPad, "GRID");
2101             relativeErrorRecBinInUp->DrawCopy("b");
2102             relativeErrorRecBinInLow->DrawCopy("same b");
2103             Style(AddLegend(gPad));
2104             relativeErrorRecBin->cd(2);
2105             Style(gPad, "GRID");
2106             relativeErrorRecBinOutUp->DrawCopy("b");
2107             relativeErrorRecBinOutLow->DrawCopy("same b");
2108             SavePadToPDF(relativeErrorRecBin);
2109             Style(AddLegend(gPad));
2110             relativeErrorRecBin->Write();
2111         }
2112     }
2113     if(methodIn && methodOut) {
2114         DoIntermediateSystematics(
2115                 methodIn, 
2116                 methodOut, 
2117                 relativeErrorMethodInUp,       // pointer reference
2118                 relativeErrorMethodInLow,      // pointer reference
2119                 relativeErrorMethodOutUp,      // pointer reference
2120                 relativeErrorMethodOutLow,     // pointer reference
2121                 relativeStatisticalErrorIn,
2122                 relativeStatisticalErrorOut,
2123                 nominal,
2124                 nominalIn,
2125                 nominalOut,
2126                 columns, 
2127                 rangeLow, 
2128                 rangeUp,
2129                 readMe,
2130                 "method"
2131                 );
2132         if(relativeErrorMethodInUp) {
2133             TCanvas* relativeErrorMethod(new TCanvas("relativeErrorMethod", "relativeErrorMethod"));
2134             relativeErrorMethod->Divide(2);
2135             relativeErrorMethod->cd(1);
2136             Style(gPad, "GRID");
2137             relativeErrorMethodInUp->DrawCopy("b");
2138             relativeErrorMethodInLow->DrawCopy("same b");
2139             Style(AddLegend(gPad));
2140             relativeErrorMethod->cd(2);
2141             Style(gPad, "GRID");
2142             relativeErrorMethodOutUp->DrawCopy("b");
2143             relativeErrorMethodOutLow->DrawCopy("same b");
2144             SavePadToPDF(relativeErrorMethod);
2145             Style(AddLegend(gPad));
2146             relativeErrorMethod->Write();
2147         }
2148     }
2149
2150     // and the placeholder for the final systematic
2151     TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2152     TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2153     TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2154     TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2155     relativeErrorInUp->SetYTitle("relative uncertainty");
2156     relativeErrorOutUp->SetYTitle("relative uncertainty");
2157     relativeErrorInLow->SetYTitle("relative uncertainty");
2158     relativeErrorOutLow->SetYTitle("relative uncertainty");
2159
2160     // sum of squares for the total systematic error
2161     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.), eInUp(0.);
2162     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.), eOutUp(0.);
2163     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.), eInLow(0.);
2164     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.), eOutLow(0.);
2165
2166     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2167         // for the upper bound
2168         if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
2169         if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
2170         if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
2171         if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
2172         if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
2173         if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
2174         if(relativeErrorMethodInUp) dInUp = relativeErrorMethodInUp->GetBinContent(b+1);
2175         if(relativeErrorMethodOutUp) dOutUp = relativeErrorMethodOutUp->GetBinContent(b+1);
2176         eInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp + dInUp*dInUp;
2177         if(eInUp > 0) relativeErrorInUp->SetBinContent(b+1, 1.*TMath::Sqrt(eInUp));
2178         eOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp + dOutUp*dOutUp;
2179         if(eOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, 1.*TMath::Sqrt(eOutUp));
2180         // for the lower bound
2181         if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
2182         if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
2183         if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
2184         if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
2185         if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
2186         if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
2187         if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2188         if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
2189         if(fSymmRMS) {  // take first category as symmetric
2190             aInLow = aInUp;
2191             aOutLow = aOutUp;
2192             cInLow = cInUp;
2193             cOutLow = cOutUp;   // other sources
2194             if(dInLow < dInUp) dInLow = dInUp;
2195             if(dOutLow < dOutUp) dOutLow = dOutUp;
2196         }
2197
2198         eInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow + dInLow*dInLow;
2199         if(eInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(eInLow));
2200         eOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow + dOutLow*dOutLow;
2201         if(eOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(eOutLow));
2202     }
2203     // project the estimated errors on the nominal ratio
2204     if(nominal) {
2205         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2206         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2207         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2208         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2209         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2210         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2211         Double_t lowErr(0.), upErr(0.);
2212         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2213             // add the in and out of plane errors in quadrature
2214             // take special care here: to propagate the assymetric error, we need to correlate the
2215             // InLow with OutUp (minimum value of ratio) and InUp with OutLow (maximum value of ratio)
2216             lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2217             upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2218             // set the errors 
2219             ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2220             ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2221             // get the bin width (which is the 'error' on x
2222             Double_t binWidth(nominal->GetBinWidth(i+1));
2223             axl[i] = binWidth/2.;
2224             axh[i] = binWidth/2.;
2225             // now get the coordinate for the point
2226             ax[i] = nominal->GetBinCenter(i+1);
2227             ay[i] = nominal->GetBinContent(i+1);
2228         }
2229         // save the nominal ratio
2230         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2231         shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2232         nominalError->SetName("shape uncertainty");
2233         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2234         nominalCanvas->Divide(2);
2235         nominalCanvas->cd(1);
2236         Style(nominal, kBlack);
2237         Style(nominalError, kCyan, kRatio);
2238         nominalError->SetLineColor(kCyan);
2239         nominalError->SetMarkerColor(kCyan);
2240         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2241         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2242         nominalError->DrawClone("a2");
2243         nominal->DrawCopy("same E1");
2244         Style(AddLegend(gPad));
2245         Style(gPad, "GRID");
2246         Style(nominalCanvas);
2247         // save nominal v2 and systematic errors
2248         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2249                     nominalIn,
2250                     nominalOut,
2251                     fEventPlaneRes,
2252                     "v_{2} with shape uncertainty",
2253                     relativeErrorInUp,
2254                     relativeErrorInLow,
2255                     relativeErrorOutUp,
2256                     relativeErrorOutLow));
2257         shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2258         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2259         nominalCanvas->cd(2);
2260         Style(nominalV2, kBlack);
2261         Style(nominalV2Error, kCyan, kV2);
2262         nominalV2Error->SetLineColor(kCyan);
2263         nominalV2Error->SetMarkerColor(kCyan);
2264         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2265         nominalV2Error->DrawClone("a2");
2266         nominalV2->DrawClone("same E1");
2267         Style(AddLegend(gPad));
2268         Style(gPad, "GRID");
2269         Style(nominalCanvas);
2270         SavePadToPDF(nominalCanvas);
2271         nominalCanvas->Write();
2272     }
2273
2274     TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2275     relativeError->Divide(2);
2276     relativeError->cd(1);
2277     Style(gPad, "GRID");
2278     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2279     Style(relativeErrorInUp, kBlue, kBar);
2280     Style(relativeErrorInLow, kGreen, kBar);
2281     relativeErrorInUp->DrawCopy("b");
2282     relativeErrorInLow->DrawCopy("same b");
2283     Style(relativeStatisticalErrorIn, kRed);
2284     relativeStatisticalErrorIn->DrawCopy("same");
2285     Style(AddLegend(gPad));
2286     relativeError->cd(2);
2287     Style(gPad, "GRID");
2288     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2289     Style(relativeErrorOutUp, kBlue, kBar);
2290     Style(relativeErrorOutLow, kGreen, kBar);
2291     relativeErrorOutUp->DrawCopy("b");
2292     relativeErrorOutLow->DrawCopy("same b");
2293     Style(relativeStatisticalErrorOut, kRed);
2294     relativeStatisticalErrorOut->DrawCopy("same");
2295     Style(AddLegend(gPad));
2296
2297     // write the buffered file to disk and close the file
2298     SavePadToPDF(relativeError);
2299     relativeError->Write();
2300     output->Write();
2301 //    output->Close();
2302 }
2303 //_____________________________________________________________________________
2304     void AliJetFlowTools::DoIntermediateSystematics(
2305             TArrayI* variationsIn,                  // variantions in plane
2306             TArrayI* variationsOut,                 // variantions out of plane
2307             TH1D*& relativeErrorInUp,               // pointer reference to minimum relative error histogram in plane
2308             TH1D*& relativeErrorInLow,              // pointer reference to maximum relative error histogram in plane
2309             TH1D*& relativeErrorOutUp,              // pointer reference to minimum relative error histogram out of plane
2310             TH1D*& relativeErrorOutLow,             // pointer reference to maximum relative error histogram out of plane
2311             TH1D*& relativeStatisticalErrorIn,      // relative systematic error on ratio
2312             TH1D*& relativeStatisticalErrorOut,     // relative systematic error on ratio
2313             TH1D*& nominal,                         // clone of the nominal ratio in / out of plane
2314             TH1D*& nominalIn,                       // clone of the nominal in plane yield
2315             TH1D*& nominalOut,                      // clone of the nominal out of plane yield
2316             Int_t columns,                          // divide the output canvasses in this many columns
2317             Float_t rangeLow,                       // lower pt range
2318             Float_t rangeUp,                        // upper pt range
2319             TFile* readMe,                          // input file name (created by this unfolding class)
2320             TString source,                         // source of the variation
2321             Bool_t RMS                              // return RMS of distribution of variations as error
2322             ) const
2323 {
2324    // intermediate systematic check function. first index of supplied array is nominal value
2325    TList* listOfKeys((TList*)readMe->GetListOfKeys());
2326    if(!listOfKeys) {
2327        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2328        return;
2329    }
2330    // check input params
2331    if(variationsIn->GetSize() != variationsOut->GetSize()) {
2332        printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2333        return;
2334    }
2335    TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2336    TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2337    if(!(defRootDirIn && defRootDirOut)) {
2338        printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2339        return;
2340    }
2341    TString defIn(defRootDirIn->GetName());
2342    TString defOut(defRootDirOut->GetName());
2343
2344    // define lines to make the output prettier
2345    TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2346    TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2347    lineLow->SetLineColor(11);
2348    lineUp->SetLineColor(11);
2349    lineLow->SetLineWidth(3);
2350    lineUp->SetLineWidth(3);
2351
2352    // define an output histogram with the maximum relative error from this systematic constribution
2353    // if the option RMS is set to false, sigma is not really a standard deviation but holds the maximum (or minimum) relative value that the data has
2354    // reached in this function call.
2355    // if the option RMS is set to true, sigma holds the RMS value (equal to sigma as the mean is zero for relative errors) of the distribution of variations
2356    // which should correspond to a 68% confidence level
2357    relativeErrorInUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2358    relativeErrorInLow = new TH1D(Form("relative error (low) from  %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2359    relativeErrorOutUp = new TH1D(Form("relative error (up) from  %s", source.Data()), Form("relative error (up)  from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2360    relativeErrorOutLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2361    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2362        if(!RMS) {
2363            relativeErrorInUp->SetBinContent(b+1, 1.);
2364            relativeErrorInUp->SetBinError(b+1, 0.);
2365            relativeErrorOutUp->SetBinContent(b+1, 1.);
2366            relativeErrorOutUp->SetBinError(b+1, .0);
2367            relativeErrorInLow->SetBinContent(b+1, 1.);
2368            relativeErrorInLow->SetBinError(b+1, 0.);
2369            relativeErrorOutLow->SetBinContent(b+1, 1.);
2370            relativeErrorOutLow->SetBinError(b+1, .0);
2371        } else if(RMS) {
2372            relativeErrorInUp->SetBinContent(b+1, 0.);
2373            relativeErrorInUp->SetBinError(b+1, 0.);
2374            relativeErrorOutUp->SetBinContent(b+1, 0.);
2375            relativeErrorOutUp->SetBinError(b+1, 0.);
2376            relativeErrorInLow->SetBinContent(b+1, 0.);
2377            relativeErrorInLow->SetBinError(b+1, 0.);
2378            relativeErrorOutLow->SetBinContent(b+1, 0.);
2379            relativeErrorOutLow->SetBinError(b+1, 0.);
2380        } 
2381    }
2382    Int_t relativeErrorInUpN[100] = {0};
2383    Int_t relativeErrorOutUpN[100] = {0};
2384    Int_t relativeErrorInLowN[100] = {0};
2385    Int_t relativeErrorOutLowN[100] = {0};
2386    
2387    // define an output histogram with the systematic error from this systematic constribution
2388    if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2389        relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "relative statistical error, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2390        relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "relative statistital error, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2391    }
2392
2393    // prepare necessary canvasses
2394    TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2395    TCanvas* canvasOut(0x0);
2396    if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2397    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2398    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2399    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2400    TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data()))); 
2401    TCanvas* canvasSpectraOut(0x0);
2402    if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2403    TCanvas* canvasRatio(0x0);
2404    if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2405    TCanvas* canvasV2(0x0);
2406    if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2407    TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2408    TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2409    TCanvas* canvasMasterOut(0x0);
2410    if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2411    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2412
2413    TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2414    canvasProfiles->Divide(2);
2415    TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2416    TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2417    // get an estimate of the number of outputs and find the default set
2418
2419    Int_t rows = 1;
2420    columns = variationsIn->GetSize()-1;
2421    (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2422    canvasIn->Divide(columns, rows);
2423    if(canvasOut) canvasOut->Divide(columns, rows);
2424    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2425    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2426    canvasSpectraIn->Divide(columns, rows);
2427    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2428    if(canvasRatio) canvasRatio->Divide(columns, rows);
2429    if(canvasV2) canvasV2->Divide(columns, rows);
2430    canvasMasterIn->Divide(columns, rows);
2431    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2432
2433    // prepare a separate set of canvases to hold the nominal points
2434    TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2435    TCanvas* canvasNominalOut(0x0);
2436    if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2437    TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2438    TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2439    if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2440    TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data()))); 
2441    TCanvas* canvasNominalSpectraOut(0x0);
2442    if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()),  Form("Nominal_%s_SpectraOut", source.Data()));
2443    TCanvas* canvasNominalRatio(0x0);
2444    if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2445    TCanvas* canvasNominalV2(0x0);
2446    if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2447    TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2448    TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2449    TCanvas* canvasNominalMasterOut(0x0);
2450    if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2451    (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2452
2453    canvasNominalSpectraIn->Divide(2);
2454    if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2455
2456    canvasNominalMasterIn->Divide(2);
2457    if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2458
2459    // extract the default output 
2460    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2461    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2462    TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2463    TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2464    if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2465    if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2466    printf(" > succesfully extracted default results < \n");
2467
2468    // loop through the directories, only plot the graphs if the deconvolution converged
2469    TDirectoryFile* tempDirIn(0x0); 
2470    TDirectoryFile* tempDirOut(0x0);
2471    TDirectoryFile* tempIn(0x0);
2472    TDirectoryFile* tempOut(0x0);
2473    TH1D* unfoldedSpectrumInForRatio(0x0);
2474    TH1D* unfoldedSpectrumOutForRatio(0x0);
2475    for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2476        tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2477        tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2478        if(!(tempDirIn && tempDirOut)) {
2479            printf(" > DoSystematics: couldn't get a set of variations < \n");
2480            continue;
2481        }
2482        TString dirNameIn(tempDirIn->GetName());
2483        TString dirNameOut(tempDirOut->GetName());
2484        // try to read the in- and out of plane subdirs
2485        tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2486        tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2487        j++;
2488        if(tempIn) { 
2489            // to see if the unfolding converged try to extract the pearson coefficients
2490            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2491            if(pIn) {
2492                printf(" - %s in plane converged \n", dirNameIn.Data());
2493                canvasIn->cd(j);
2494                if(i==0) canvasNominalIn->cd(j);
2495                Style(gPad, "PEARSON");
2496                pIn->DrawCopy("colz");
2497                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2498                if(rIn) {
2499                    Style(rIn);
2500                    printf(" > found RatioRefoldedMeasured < \n");
2501                    canvasRatioMeasuredRefoldedIn->cd(j);
2502                    if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2503                    Style(gPad, "GRID");
2504                    rIn->SetFillColor(kRed);
2505                    rIn->Draw("ap");
2506                }
2507                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2508                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2509                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2510                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2511                if(dvector && avalue && rm && eff) {
2512                    Style(dvector);
2513                    Style(avalue);
2514                    Style(rm);
2515                    Style(eff);
2516                    canvasMISC->cd(1);
2517                    if(i==0) canvasNominalMISC->cd(1);
2518                    Style(gPad, "SPECTRUM");
2519                    dvector->DrawCopy();
2520                    canvasMISC->cd(2);
2521                    if(i==0) canvasNominalMISC->cd(2);
2522                    Style(gPad, "SPECTRUM");
2523                    avalue->DrawCopy();
2524                    canvasMISC->cd(3);
2525                    if(i==0) canvasNominalMISC->cd(3);
2526                    Style(gPad, "PEARSON");
2527                    rm->DrawCopy("colz");
2528                    canvasMISC->cd(4);
2529                    if(i==0) canvasNominalMISC->cd(4);
2530                    Style(gPad, "GRID");
2531                    eff->DrawCopy();
2532                } else if(rm && eff) {
2533                    Style(rm);
2534                    Style(eff);
2535                    canvasMISC->cd(3);
2536                    if(i==0) canvasNominalMISC->cd(3);
2537                    Style(gPad, "PEARSON");
2538                    rm->DrawCopy("colz");
2539                    canvasMISC->cd(4);
2540                    if(i==0) canvasNominalMISC->cd(4);
2541                    Style(gPad, "GRID");
2542                    eff->DrawCopy();
2543                }
2544            }
2545            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2546            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2547            unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2548            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2549            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2550                if(defaultUnfoldedJetSpectrumIn) {
2551                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2552                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2553                    temp->Divide(unfoldedSpectrum);
2554                    // get the absolute relative error
2555                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2556                        if(!RMS) {       // save the maximum deviation that a variation can cause
2557                            // the variation is HIGHER than the nominal point, so the bar goes UP
2558                            if( temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInUp->GetBinContent(b+1)) {
2559                                relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2560                                relativeErrorInUp->SetBinError(b+1, 0.);
2561                            }
2562                            // the variation is LOWER than the nominal point, so the bar goes DOWN
2563                            else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInLow->GetBinContent(b+1)) {
2564                                relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2565                                relativeErrorInLow->SetBinError(b+1, 0.);
2566                            }
2567                        } else if (RMS && !fSymmRMS) { // save info necessary for evaluating the RMS of a distribution of variations
2568                            printf(" oops shouldnt be here \n " );
2569                            if(temp->GetBinContent(b+1) < 1) {
2570                                relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2571                                relativeErrorInUpN[b]++;
2572                            }
2573                            // the variation is LOWER than the nominal point, so the bar goes DOWN
2574                            else if(temp->GetBinContent(b+1) > 1) {
2575                                relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2576                                relativeErrorInLowN[b]++;
2577                            }
2578                        } else if (fSymmRMS) {
2579                            // save symmetric sum of square to get a symmetric rms
2580                            relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2581                            relativeErrorInUpN[b]++;
2582                        }
2583                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2584                    }
2585                    temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2586                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2587                    temp->GetYaxis()->SetTitle("ratio");
2588                    canvasMasterIn->cd(j);
2589                    temp->GetYaxis()->SetRangeUser(0., 2);
2590                    Style(gPad, "GRID");
2591                    temp->DrawCopy();
2592                    canvasNominalMasterIn->cd(1);
2593                    Style(gPad, "GRID");
2594                    if(i > 0 ) {
2595                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2596                        tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2597                        Style(tempSyst, (EColor)(i+2));
2598                        if(i==1) tempSyst->DrawCopy();
2599                        else tempSyst->DrawCopy("same");
2600                    }
2601                }
2602                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2603                canvasSpectraIn->cd(j);
2604                if(i==0) canvasNominalSpectraIn->cd(1);
2605                Style(gPad);
2606                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2607                unfoldedSpectrum->DrawCopy();
2608                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2609                inputSpectrum->DrawCopy("same");
2610                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2611                refoldedSpectrum->DrawCopy("same");
2612                TLegend* l(AddLegend(gPad));
2613                Style(l);
2614                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2615                    Float_t chi(fitStatus->GetBinContent(1));
2616                    Float_t pen(fitStatus->GetBinContent(2));
2617                    Int_t dof((int)fitStatus->GetBinContent(3));
2618                    Float_t beta(fitStatus->GetBinContent(4));
2619                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2620                } else if (fitStatus) { // only available in SVD fit
2621                    Int_t reg((int)fitStatus->GetBinContent(1));
2622                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2623                }
2624                canvasNominalSpectraIn->cd(2);
2625                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2626                tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2627                Style(tempSyst, (EColor)(i+2));
2628                Style(gPad, "SPECTRUM");
2629                if(i==0) tempSyst->DrawCopy();
2630                else tempSyst->DrawCopy("same");
2631            }
2632        }
2633        if(tempOut) {
2634            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2635            if(pOut) {
2636                printf(" - %s out of plane converged \n", dirNameOut.Data());
2637                canvasOut->cd(j);
2638                if(i==0) canvasNominalOut->cd(j);
2639                Style(gPad, "PEARSON");
2640                pOut->DrawCopy("colz");
2641                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2642                if(rOut) {
2643                    Style(rOut);
2644                    printf(" > found RatioRefoldedMeasured < \n");
2645                    canvasRatioMeasuredRefoldedOut->cd(j);
2646                    if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2647                    Style(gPad, "GRID");
2648                    rOut->SetFillColor(kRed);
2649                    rOut->Draw("ap");
2650                }
2651                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2652                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2653                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2654                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2655                if(dvector && avalue && rm && eff) {
2656                    Style(dvector);
2657                    Style(avalue);
2658                    Style(rm);
2659                    Style(eff);
2660                    canvasMISC->cd(5);
2661                    if(i==0) canvasNominalMISC->cd(5);
2662                    Style(gPad, "SPECTRUM");
2663                    dvector->DrawCopy();
2664                    canvasMISC->cd(6);
2665                    if(i==0) canvasNominalMISC->cd(6);
2666                    Style(gPad, "SPECTRUM");
2667                    avalue->DrawCopy();
2668                    canvasMISC->cd(7);
2669                    if(i==0) canvasNominalMISC->cd(7);
2670                    Style(gPad, "PEARSON");
2671                    rm->DrawCopy("colz");
2672                    canvasMISC->cd(8);
2673                    if(i==0) canvasNominalMISC->cd(8);
2674                    Style(gPad, "GRID");
2675                    eff->DrawCopy();
2676                } else if(rm && eff) {
2677                    Style(rm);
2678                    Style(eff);
2679                    canvasMISC->cd(7);
2680                    if(i==0) canvasNominalMISC->cd(7);
2681                    Style(gPad, "PEARSON");
2682                    rm->DrawCopy("colz");
2683                    canvasMISC->cd(8);
2684                    if(i==0) canvasNominalMISC->cd(8);
2685                    Style(gPad, "GRID");
2686                    eff->DrawCopy();
2687                }
2688            }
2689            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2690            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2691            unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2692            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2693            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2694                if(defaultUnfoldedJetSpectrumOut) {
2695                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2696                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2697                    temp->Divide(unfoldedSpectrum);
2698                    // get the absolute relative error 
2699                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2700                        if(!RMS) {
2701                            // check if the error is larger than the current maximum
2702                            if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutUp->GetBinContent(b+1)) {
2703                                relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2704                                relativeErrorOutUp->SetBinError(b+1, 0.);
2705                            }
2706                            // check if the error is smaller than the current minimum
2707                            else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutLow->GetBinContent(b+1)) {
2708                                relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2709                                relativeErrorOutLow->SetBinError(b+1, 0.);
2710                            }
2711                        } else if (RMS && !fSymmRMS) {
2712                            printf(" OOps \n ");
2713                            if(temp->GetBinContent(b+1) < 1) {
2714                                relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2715                                relativeErrorOutUpN[b]++;
2716                            }
2717                            else if(temp->GetBinContent(b+1) > 1) {
2718                                relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2719                                relativeErrorOutLowN[b]++;
2720                            }
2721                        } else if (fSymmRMS) {
2722                            // save symmetric rms value
2723                            relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2724                            relativeErrorOutUpN[b]++;
2725                        }
2726                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2727                    }
2728                    temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2729                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2730                    temp->GetYaxis()->SetTitle("ratio");
2731                    canvasMasterOut->cd(j);
2732                    temp->GetYaxis()->SetRangeUser(0., 2);
2733                    Style(gPad, "GRID");
2734                    temp->DrawCopy();
2735                    canvasNominalMasterOut->cd(1);
2736                    Style(gPad, "GRID");
2737                    if(i > 0 ) {
2738                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2739                        tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2740                        Style(tempSyst, (EColor)(i+2));
2741                        if(i==1) tempSyst->DrawCopy();
2742                        else tempSyst->DrawCopy("same");
2743                    }
2744                }
2745                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2746                canvasSpectraOut->cd(j);
2747                if(i==0) canvasNominalSpectraOut->cd(1);
2748                Style(gPad);
2749                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2750                unfoldedSpectrum->DrawCopy();
2751                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2752                inputSpectrum->DrawCopy("same");
2753                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2754                refoldedSpectrum->DrawCopy("same");
2755                TLegend* l(AddLegend(gPad));
2756                Style(l);
2757                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2758                    Float_t chi(fitStatus->GetBinContent(1));
2759                    Float_t pen(fitStatus->GetBinContent(2));
2760                    Int_t dof((int)fitStatus->GetBinContent(3));
2761                    Float_t beta(fitStatus->GetBinContent(4));
2762                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2763                } else if (fitStatus) { // only available in SVD fit
2764                    Int_t reg((int)fitStatus->GetBinContent(1));
2765                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2766                }
2767                canvasNominalSpectraOut->cd(2);
2768                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2769                tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2770                Style(tempSyst, (EColor)(i+2));
2771                Style(gPad, "SPECTRUM");
2772                if(i==0) tempSyst->DrawCopy();
2773                else tempSyst->DrawCopy("same");
2774            }
2775        }
2776        if(canvasRatio && canvasV2) {
2777            canvasRatio->cd(j);
2778            if(i==0) {
2779                canvasNominalRatio->cd(j);
2780                if(nominal && nominalIn && nominalOut) {
2781                    // if a nominal ratio is requested, delete the placeholder and update the nominal point
2782                    delete nominal;
2783                    delete nominalIn;
2784                    delete nominalOut;
2785                    nominalIn =  (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2786                    nominalOut =  (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2787                    nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2788                    nominal->Divide(nominal, unfoldedSpectrumOutForRatio);       // standard root divide for uncorrelated histos
2789                }
2790            }
2791            TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2792            Double_t _tempx(0), _tempy(0);
2793            if(ratioYield) {
2794                Style(ratioYield);
2795                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2796                    ratioYield->GetPoint(b, _tempx, _tempy);
2797                    ratioProfile->Fill(_tempx, _tempy);
2798                }
2799                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2800                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2801                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2802                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2803                ratioYield->SetFillColor(kRed);
2804                ratioYield->Draw("ap");
2805            }
2806            canvasV2->cd(j);
2807            if(i==0) canvasNominalV2->cd(j);
2808            TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2809            if(ratioV2) {
2810                Style(ratioV2);
2811                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2812                    ratioV2->GetPoint(b, _tempx, _tempy);
2813                    v2Profile->Fill(_tempx, _tempy);
2814
2815                }
2816                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2817                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2818                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2819                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2820                ratioV2->SetFillColor(kRed);
2821                ratioV2->Draw("ap");
2822            }
2823        }
2824        delete unfoldedSpectrumInForRatio;
2825        delete unfoldedSpectrumOutForRatio;
2826    }
2827    // save the canvasses
2828    canvasProfiles->cd(1);
2829    Style(ratioProfile);
2830    ratioProfile->DrawCopy();
2831    canvasProfiles->cd(2);
2832    Style(v2Profile);
2833    v2Profile->DrawCopy();
2834    SavePadToPDF(canvasProfiles);
2835    canvasProfiles->Write(); 
2836    SavePadToPDF(canvasIn);
2837    canvasIn->Write();
2838    if(canvasOut) {
2839        SavePadToPDF(canvasOut);
2840        canvasOut->Write();
2841    }
2842    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2843    canvasRatioMeasuredRefoldedIn->Write();
2844    if(canvasRatioMeasuredRefoldedOut) {
2845        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2846        canvasRatioMeasuredRefoldedOut->Write();
2847    }
2848    SavePadToPDF(canvasSpectraIn);
2849    canvasSpectraIn->Write();
2850    if(canvasSpectraOut) {
2851        SavePadToPDF(canvasSpectraOut);
2852        canvasSpectraOut->Write();
2853        SavePadToPDF(canvasRatio);
2854        canvasRatio->Write();
2855        SavePadToPDF(canvasV2);
2856        canvasV2->Write();
2857    }
2858    SavePadToPDF(canvasMasterIn);
2859    canvasMasterIn->Write();
2860    if(canvasMasterOut) {
2861        SavePadToPDF(canvasMasterOut);
2862        canvasMasterOut->Write();
2863    }
2864    SavePadToPDF(canvasMISC);
2865    canvasMISC->Write();
2866    // save the nomial canvasses
2867    SavePadToPDF(canvasNominalIn);
2868    canvasNominalIn->Write();
2869    if(canvasNominalOut) {
2870        SavePadToPDF(canvasNominalOut);
2871        canvasNominalOut->Write();
2872    }
2873    SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
2874    canvasNominalRatioMeasuredRefoldedIn->Write();
2875    if(canvasNominalRatioMeasuredRefoldedOut) {
2876        SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
2877        canvasNominalRatioMeasuredRefoldedOut->Write();
2878    }
2879    canvasNominalSpectraIn->cd(2);
2880    Style(AddLegend(gPad)); 
2881    SavePadToPDF(canvasNominalSpectraIn);
2882    canvasNominalSpectraIn->Write();
2883    if(canvasNominalSpectraOut) {
2884        canvasNominalSpectraOut->cd(2);
2885        Style(AddLegend(gPad));
2886        SavePadToPDF(canvasNominalSpectraOut);
2887        canvasNominalSpectraOut->Write();
2888        SavePadToPDF(canvasNominalRatio);
2889        canvasNominalRatio->Write();
2890        SavePadToPDF(canvasNominalV2);
2891        canvasNominalV2->Write();
2892    }
2893    canvasNominalMasterIn->cd(1);
2894    Style(AddLegend(gPad));
2895    lineUp->DrawClone("same");
2896    lineLow->DrawClone("same");
2897    SavePadToPDF(canvasNominalMasterIn);
2898    canvasNominalMasterIn->Write();
2899    if(canvasNominalMasterOut) {
2900        canvasNominalMasterOut->cd(1);
2901        Style(AddLegend(gPad));
2902        lineUp->DrawClone("same");
2903        lineLow->DrawClone("same");
2904        SavePadToPDF(canvasNominalMasterOut);
2905        canvasNominalMasterOut->Write();
2906    }
2907    SavePadToPDF(canvasNominalMISC);
2908    canvasNominalMISC->Write();
2909
2910    // save the relative errors
2911    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2912        // to arrive at a min and max from here, combine in up and out low
2913        if(!RMS) {
2914            relativeErrorInUp->SetBinContent(b+1, -1.*(relativeErrorInUp->GetBinContent(b+1)-1));
2915            relativeErrorInUp->SetBinError(b+1, 0.);
2916            relativeErrorOutUp->SetBinContent(b+1, -1.*(relativeErrorOutUp->GetBinContent(b+1)-1));
2917            relativeErrorOutUp->SetBinError(b+1, .0);
2918            relativeErrorInLow->SetBinContent(b+1, -1.*(relativeErrorInLow->GetBinContent(b+1)-1));
2919            relativeErrorInLow->SetBinError(b+1, 0.);
2920            relativeErrorOutLow->SetBinContent(b+1, -1.*(relativeErrorOutLow->GetBinContent(b+1)-1));
2921            relativeErrorOutLow->SetBinError(b+1, .0);
2922        } else if (RMS) {
2923            // these guys are already stored as percentages, so no need to remove the offset of 1
2924            // RMS is defined as sqrt(sum(squared))/N
2925            // min is defined as negative, max is defined as positive
2926            if(!fSymmRMS) {
2927                if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
2928                if(relativeErrorInLowN[b] < 1) relativeErrorInLowN[b] = 1;
2929                if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
2930                if(relativeErrorOutLowN[b] < 1) relativeErrorOutLowN[b] = 1;
2931                relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
2932                relativeErrorInUp->SetBinError(b+1, 0.);
2933                relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
2934                relativeErrorOutUp->SetBinError(b+1, .0);
2935                relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorInLow->GetBinContent(b+1)/relativeErrorInLowN[b]));
2936                relativeErrorInLow->SetBinError(b+1, 0.);
2937                relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorOutLow->GetBinContent(b+1)/relativeErrorOutLowN[b]));
2938                relativeErrorOutLow->SetBinError(b+1, .0);
2939            } else if (fSymmRMS) {
2940                if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
2941                if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
2942                relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
2943                relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
2944            }
2945        }
2946    }
2947    relativeErrorInUp->SetYTitle("relative uncertainty");
2948    relativeErrorOutUp->SetYTitle("relative uncertainty");
2949    relativeErrorInLow->SetYTitle("relative uncertainty");
2950    relativeErrorOutLow->SetYTitle("relative uncertainty");
2951    relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2952    relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2953    relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2954    relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2955
2956    canvasNominalMasterIn->cd(2);
2957    Style(gPad, "GRID");
2958    Style(relativeErrorInUp, kBlue, kBar);
2959    Style(relativeErrorInLow, kGreen, kBar);
2960    relativeErrorInUp->DrawCopy("b");
2961    relativeErrorInLow->DrawCopy("same b");
2962    Style(AddLegend(gPad));
2963    SavePadToPDF(canvasNominalMasterIn);
2964    canvasNominalMasterIn->Write();
2965    canvasNominalMasterOut->cd(2);
2966    Style(gPad, "GRID");
2967    Style(relativeErrorOutUp, kBlue, kBar);
2968    Style(relativeErrorOutLow, kGreen, kBar);
2969    relativeErrorOutUp->DrawCopy("b");
2970    relativeErrorOutLow->DrawCopy("same b");
2971    Style(AddLegend(gPad));
2972    SavePadToPDF(canvasNominalMasterOut);
2973    canvasNominalMasterOut->Write();
2974 }
2975 //_____________________________________________________________________________
2976 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
2977 {
2978    // go through the output file and perform post processing routines
2979    // can either be performed in one go with the unfolding, or at a later stage
2980    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
2981    TFile readMe(in.Data(), "READ");     // open file read-only
2982    if(readMe.IsZombie()) {
2983        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2984        return;
2985    }
2986    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
2987    readMe.ls();
2988    TList* listOfKeys((TList*)readMe.GetListOfKeys());
2989    if(!listOfKeys) {
2990        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2991        return;
2992    }
2993    // prepare necessary canvasses
2994    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
2995    TCanvas* canvasOut(0x0);
2996    if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
2997    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
2998    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2999    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
3000    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn")); 
3001    TCanvas* canvasSpectraOut(0x0);
3002    if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
3003    TCanvas* canvasRatio(0x0);
3004    if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
3005    TCanvas* canvasV2(0x0);
3006    if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
3007    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
3008    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
3009    TCanvas* canvasMasterOut(0x0);
3010    if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
3011    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
3012    TDirectoryFile* defDir(0x0);
3013    TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
3014    canvasProfiles->Divide(2);
3015    TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3016    TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3017    // get an estimate of the number of outputs and find the default set
3018    Int_t cacheMe(0);
3019    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3020        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3021            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3022            cacheMe++;
3023        }
3024    }
3025    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
3026    canvasIn->Divide(columns, rows);
3027    if(canvasOut) canvasOut->Divide(columns, rows);
3028    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
3029    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
3030    canvasSpectraIn->Divide(columns, rows);
3031    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
3032    if(canvasRatio) canvasRatio->Divide(columns, rows);
3033    if(canvasV2) canvasV2->Divide(columns, rows);
3034
3035    canvasMasterIn->Divide(columns, rows);
3036    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
3037    // extract the default output 
3038    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3039    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3040    if(defDir) {
3041        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3042        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3043        if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3044        if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3045        printf(" > succesfully extracted default results < \n");
3046    }
3047
3048    // loop through the directories, only plot the graphs if the deconvolution converged
3049    TDirectoryFile* tempDir(0x0); 
3050    TDirectoryFile* tempIn(0x0);
3051    TDirectoryFile*  tempOut(0x0);
3052    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3053        // read the directory by its name
3054        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3055        if(!tempDir) continue;
3056        TString dirName(tempDir->GetName());
3057        // try to read the in- and out of plane subdirs
3058        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3059        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3060        j++;
3061        if(!tempIn) {    // attempt to get MakeAU output
3062            TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3063            TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
3064            tempPearson->Divide(4,2);
3065            TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
3066            tempRatio->Divide(4,2);
3067            TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
3068            tempResult->Divide(4,2);
3069            for(Int_t q(0); q < 8; q++) {
3070                tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
3071                if(tempIn) {
3072                        // to see if the unfolding converged try to extract the pearson coefficients
3073                        TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3074                        if(pIn) {
3075                        printf(" - %s in plane converged \n", dirName.Data());
3076                            tempPearson->cd(1+q);
3077                             Style(gPad, "PEARSON");
3078                             pIn->DrawCopy("colz");
3079                             TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3080                             if(rIn) {
3081                                 Style(rIn);
3082                                 printf(" > found RatioRefoldedMeasured < \n");
3083                                 tempRatio->cd(q+1);
3084                                 rIn->SetFillColor(kRed);
3085                                 rIn->Draw("ap");
3086                             }
3087                             TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3088                             TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3089                             TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3090                             TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3091                             if(dvector && avalue && rm && eff) {
3092                                 Style(dvector);
3093                                 Style(avalue);
3094                                 Style(rm);
3095                                 Style(eff);
3096                                 canvasMISC->cd(1);
3097                                 Style(gPad, "SPECTRUM");
3098                                 dvector->DrawCopy();
3099                                 canvasMISC->cd(2);
3100                                 Style(gPad, "SPECTRUM");
3101                                 avalue->DrawCopy();
3102                                 canvasMISC->cd(3);
3103                                 Style(gPad, "PEARSON");
3104                                 rm->DrawCopy("colz");
3105                                 canvasMISC->cd(4);
3106                                 Style(gPad, "GRID");
3107                                 eff->DrawCopy();
3108                             } else if(rm && eff) {
3109                                 Style(rm);
3110                                 Style(eff);
3111                                 canvasMISC->cd(3);
3112                                 Style(gPad, "PEARSON");
3113                                 rm->DrawCopy("colz");
3114                                 canvasMISC->cd(4);
3115                                 Style(gPad, "GRID");
3116                                 eff->DrawCopy();
3117                             }
3118                         }
3119                        TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3120                        TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3121                        TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3122                        if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3123                            if(defaultUnfoldedJetSpectrumIn) {
3124                                Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3125                                TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3126                                temp->Divide(unfoldedSpectrum);
3127                                temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3128                                temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3129                                temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3130                                canvasMasterIn->cd(j);
3131                                temp->GetYaxis()->SetRangeUser(0., 2);
3132                                temp->DrawCopy();
3133                            }
3134                            TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3135                            tempResult->cd(q+1);
3136                            Style(gPad);
3137                            Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3138                            unfoldedSpectrum->DrawCopy();
3139                            Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3140                            inputSpectrum->DrawCopy("same");
3141                            Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3142                            refoldedSpectrum->DrawCopy("same");
3143                            TLegend* l(AddLegend(gPad));
3144                            Style(l);
3145                            if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3146                                Float_t chi(fitStatus->GetBinContent(1));
3147                                Float_t pen(fitStatus->GetBinContent(2));
3148                                Int_t dof((int)fitStatus->GetBinContent(3));
3149                                Float_t beta(fitStatus->GetBinContent(4));
3150                                l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3151                            } else if (fitStatus) { // only available in SVD fit
3152                                Int_t reg((int)fitStatus->GetBinContent(1));
3153                                l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3154                            }
3155                        }
3156                }
3157            }
3158        }
3159        if(tempIn) { 
3160            // to see if the unfolding converged try to extract the pearson coefficients
3161            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3162            if(pIn) {
3163                printf(" - %s in plane converged \n", dirName.Data());
3164                canvasIn->cd(j);
3165                Style(gPad, "PEARSON");
3166                pIn->DrawCopy("colz");
3167                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3168                if(rIn) {
3169                    Style(rIn);
3170                    printf(" > found RatioRefoldedMeasured < \n");
3171                    canvasRatioMeasuredRefoldedIn->cd(j);
3172                    rIn->SetFillColor(kRed);
3173                    rIn->Draw("ap");
3174                }
3175                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3176                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3177                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3178                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3179                if(dvector && avalue && rm && eff) {
3180                    Style(dvector);
3181                    Style(avalue);
3182                    Style(rm);
3183                    Style(eff);
3184                    canvasMISC->cd(1);
3185                    Style(gPad, "SPECTRUM");
3186                    dvector->DrawCopy();
3187                    canvasMISC->cd(2);
3188                    Style(gPad, "SPECTRUM");
3189                    avalue->DrawCopy();
3190                    canvasMISC->cd(3);
3191                    Style(gPad, "PEARSON");
3192                    rm->DrawCopy("colz");
3193                    canvasMISC->cd(4);
3194                    Style(gPad, "GRID");
3195                    eff->DrawCopy();
3196                } else if(rm && eff) {
3197                    Style(rm);
3198                    Style(eff);
3199                    canvasMISC->cd(3);
3200                    Style(gPad, "PEARSON");
3201                    rm->DrawCopy("colz");
3202                    canvasMISC->cd(4);
3203                    Style(gPad, "GRID");
3204                    eff->DrawCopy();
3205                }
3206            }
3207            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3208            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3209            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3210            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3211                if(defaultUnfoldedJetSpectrumIn) {
3212                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3213                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3214                    temp->Divide(unfoldedSpectrum);
3215                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3216                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3217                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3218                    canvasMasterIn->cd(j);
3219                    temp->GetYaxis()->SetRangeUser(0., 2);
3220                    temp->DrawCopy();
3221                }
3222                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3223                canvasSpectraIn->cd(j);
3224                Style(gPad);
3225                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3226                unfoldedSpectrum->DrawCopy();
3227                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3228                inputSpectrum->DrawCopy("same");
3229                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3230                refoldedSpectrum->DrawCopy("same");
3231                TLegend* l(AddLegend(gPad));
3232                Style(l);
3233                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3234                    Float_t chi(fitStatus->GetBinContent(1));
3235                    Float_t pen(fitStatus->GetBinContent(2));
3236                    Int_t dof((int)fitStatus->GetBinContent(3));
3237                    Float_t beta(fitStatus->GetBinContent(4));
3238                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3239                } else if (fitStatus) { // only available in SVD fit
3240                    Int_t reg((int)fitStatus->GetBinContent(1));
3241                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3242                }
3243            }
3244        }
3245        if(tempOut) {
3246            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
3247            if(pOut) {
3248                printf(" - %s out of plane converged \n", dirName.Data());
3249                canvasOut->cd(j);
3250                Style(gPad, "PEARSON");
3251                pOut->DrawCopy("colz");
3252                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3253                if(rOut) {
3254                    Style(rOut);
3255                    printf(" > found RatioRefoldedMeasured < \n");
3256                    canvasRatioMeasuredRefoldedOut->cd(j);
3257                    rOut->SetFillColor(kRed);
3258                    rOut->Draw("ap");
3259                }
3260                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
3261                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
3262                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
3263                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
3264                if(dvector && avalue && rm && eff) {
3265                    Style(dvector);
3266                    Style(avalue);
3267                    Style(rm);
3268                    Style(eff);
3269                    canvasMISC->cd(5);
3270                    Style(gPad, "SPECTRUM");
3271                    dvector->DrawCopy();
3272                    canvasMISC->cd(6);
3273                    Style(gPad, "SPECTRUM");
3274                    avalue->DrawCopy();
3275                    canvasMISC->cd(7);
3276                    Style(gPad, "PEARSON");
3277                    rm->DrawCopy("colz");
3278                    canvasMISC->cd(8);
3279                    Style(gPad, "GRID");
3280                    eff->DrawCopy();
3281                } else if(rm && eff) {
3282                    Style(rm);
3283                    Style(eff);
3284                    canvasMISC->cd(7);
3285                    Style(gPad, "PEARSON");
3286                    rm->DrawCopy("colz");
3287                    canvasMISC->cd(8);
3288                    Style(gPad, "GRID");
3289                    eff->DrawCopy();
3290                }
3291            }
3292            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3293            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3294            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3295            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3296                if(defaultUnfoldedJetSpectrumOut) {
3297                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3298                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3299                    temp->Divide(unfoldedSpectrum);
3300                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3301                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3302                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3303                    canvasMasterOut->cd(j);
3304                    temp->GetYaxis()->SetRangeUser(0., 2.);
3305                    temp->DrawCopy();
3306                }
3307                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3308                canvasSpectraOut->cd(j);
3309                Style(gPad);
3310                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3311                unfoldedSpectrum->DrawCopy();
3312                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3313                inputSpectrum->DrawCopy("same");
3314                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3315                refoldedSpectrum->DrawCopy("same");
3316                TLegend* l(AddLegend(gPad));
3317                Style(l);
3318                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3319                    Float_t chi(fitStatus->GetBinContent(1));
3320                    Float_t pen(fitStatus->GetBinContent(2));
3321                    Int_t dof((int)fitStatus->GetBinContent(3));