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