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