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