AliJetFlowTools.cxx - use ROOT error propagation by default for histogram operations...
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.cxx
1 /************************************************************************* 
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * 
3  *                                                                        * 
4  * Author: The ALICE Off-line Project.                                    * 
5  * Contributors are mentioned in the code where appropriate.              * 
6  *                                                                        * 
7  * Permission to use, copy, modify and distribute this software and its   * 
8  * documentation strictly for non-commercial purposes is hereby granted   * 
9  * without fee, provided that the above copyright notice appears in all   * 
10  * copies and that both the copyright notice and this permission notice   * 
11  * appear in the supporting documentation. The authors make no claims     * 
12  * about the suitability of this software for any purpose. It is          * 
13  * provided "as is" without express or implied warranty.                  * 
14  **************************************************************************/
15
16 // Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
17 //         (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
18 //
19 // Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
20 //
21 // The task uses input from two analysis tasks:
22 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
23 //   used to retrieve jet spectra and delta pt distributions
24 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
25 //   used to construct the detector response function
26 // and unfolds jet spectra with respect to the event plane. The user can choose
27 // different alrogithms for unfolding which are available in (ali)root. RooUnfold 
28 // libraries must be present on the system 
29 // ( see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html ).
30 // 
31 // The weak spot of this class is the function PrepareForUnfolding, which will read
32 // output from two output files and expects histograms with certain names and binning. 
33 // Unfolding methods itself are general and should be able to handle any input, therefore one
34 // can forgo the PrepareForUnfolding method, and supply necessary input information via the 
35 // SetRawInput() method
36 //
37 // to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
38
39 // root includes
40 #include "TF1.h"
41 #include "TH1D.h"
42 #include "TH2D.h"
43 #include "TGraph.h"
44 #include "TGraphErrors.h"
45 #include "TCanvas.h"
46 #include "TLegend.h"
47 #include "TArrayD.h"
48 #include "TList.h"
49 #include "TMinuit.h"
50 #include "TVirtualFitter.h"
51 #include "TLegend.h"
52 #include "TCanvas.h"
53 #include "TStyle.h"
54 #include "TLine.h"
55 #include "TMath.h"
56 #include "TVirtualFitter.h"
57 #include "TFitResultPtr.h"
58 // aliroot includes
59 #include "AliUnfolding.h"
60 #include "AliAnaChargedJetResponseMaker.h"
61 // class includes
62 #include "AliJetFlowTools.h"
63 // roo unfold includes (make sure you have these available on your system)
64 #include "RooUnfold.h"
65 #include "RooUnfoldResponse.h"
66 #include "RooUnfoldSvd.h"
67 #include "RooUnfoldBayes.h"
68 #include "TSVDUnfold.h"
69
70 using namespace std;
71 //_____________________________________________________________________________
72 AliJetFlowTools::AliJetFlowTools() :
73     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
74     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
75     fSaveFull           (kFALSE),
76     fActiveString       (""),
77     fActiveDir          (0x0),
78     fInputList          (0x0),
79     fRefreshInput       (kTRUE),
80     fOutputFileName     ("UnfoldedSpectra.root"),
81     fOutputFile         (0x0),
82     fCentralityBin      (0),
83     fDetectorResponse   (0x0),
84     fJetFindingEff      (0x0),
85     fBetaIn             (.1),
86     fBetaOut            (.1),
87     fBayesianIterIn     (4),
88     fBayesianIterOut    (4),
89     fBayesianSmoothIn   (0.),
90     fBayesianSmoothOut  (0.),
91     fAvoidRoundingError (kFALSE),
92     fUnfoldingAlgorithm (kChi2),
93     fPrior              (kPriorMeasured),
94     fBinsTrue           (0x0),
95     fBinsRec            (0x0),
96     fBinsTruePrior      (0x0),
97     fBinsRecPrior       (0x0),
98     fSVDRegIn           (5),
99     fSVDRegOut          (5),
100     fSVDToy             (kTRUE),
101     fJetRadius          (0.3),
102     fEventCount         (-1),
103     fNormalizeSpectra   (kFALSE),
104     fSmoothenPrior      (kFALSE),
105     fFitMin             (60.),
106     fFitMax             (300.),
107     fFitStart           (75.),
108     fSmoothenCounts     (kTRUE),
109     fTestMode           (kFALSE),
110     fRawInputProvided   (kFALSE),
111     fEventPlaneRes      (.63),
112     fUseDetectorResponse(kTRUE),
113     fUseDptResponse     (kTRUE),
114     fTrainPower         (kTRUE),
115     fInOutUnfolding     (kTRUE),
116     fRMSSpectrumIn      (0x0),
117     fRMSSpectrumOut     (0x0),
118     fRMSRatio           (0x0),
119     fRMSV2              (0x0),
120     fDeltaPtDeltaPhi    (0x0),
121     fJetPtDeltaPhi      (0x0),
122     fSpectrumIn         (0x0),
123     fSpectrumOut        (0x0),
124     fDptInDist          (0x0),
125     fDptOutDist         (0x0),
126     fDptIn              (0x0),
127     fDptOut             (0x0),
128     fFullResponseIn     (0x0),
129     fFullResponseOut    (0x0) { // class constructor
130     // create response maker weight function (tuned to PYTHIA spectrum)
131     fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
132 }
133 //_____________________________________________________________________________
134 void AliJetFlowTools::Make() {
135     // core function of the class
136     // 1) rebin the raw output of the jet task to the desired binnings
137     // 2) calls the unfolding routine
138     // 3) writes output to file
139     // can be repeated multiple times with different configurations
140
141     // 1) manipulation of input histograms
142     // check if the input variables are present
143     if(fRefreshInput) {
144         if(!PrepareForUnfolding()) {
145             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
146             return;
147         }
148     }
149     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
150     //     parts of the spectrum can end up in over or underflow bins
151     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
152     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
153
154     // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
155     // the template will be used as a prior for the chi2 unfolding
156     TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
157     TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
158
159     // get the full response matrix from the dpt and the detector response
160     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
161     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
162     // so that unfolding should return the initial spectrum
163     if(!fTestMode) {
164         if(fUseDptResponse && fUseDetectorResponse) {
165             fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
166             fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
167         } else if (fUseDptResponse && !fUseDetectorResponse) {
168             fFullResponseIn = fDptIn;
169             fFullResponseOut = fDptOut;
170         } else if (!fUseDptResponse && fUseDetectorResponse) {
171             fFullResponseIn = fDetectorResponse;
172             fFullResponseOut = fDetectorResponse;
173         } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
174             printf(" > No response, exiting ! < \n" );
175             return;
176         }
177     } else {
178         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
179         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
180     }
181     // normalize each slide of the response to one
182     NormalizeTH2D(fFullResponseIn);
183     NormalizeTH2D(fFullResponseOut);
184     // resize to desired binning scheme
185     TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
186     TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
187     // get the kinematic efficiency
188     TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
189     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
190     TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
191     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
192     // suppress the errors 
193     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
194         kinematicEfficiencyIn->SetBinError(1+i, 0.);
195         kinematicEfficiencyOut->SetBinError(1+i, 0.);
196     }
197     TH1D* jetFindingEfficiency(0x0);
198     if(fJetFindingEff) {
199         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
200         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
201         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
202     }
203     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
204     TH1D* unfoldedJetSpectrumIn(0x0);
205     TH1D* unfoldedJetSpectrumOut(0x0); 
206     fActiveDir->cd();                   // select active dir
207     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
208     dirIn->cd();                        // select inplane subdir
209     // select the unfolding method
210     switch (fUnfoldingAlgorithm) {
211         case kChi2 : {
212             unfoldedJetSpectrumIn = UnfoldSpectrumChi2(       // do the inplane unfolding
213                 measuredJetSpectrumIn,
214                 resizedResponseIn,
215                 kinematicEfficiencyIn,
216                 measuredJetSpectrumTrueBinsIn,
217                 TString("in"),
218                 jetFindingEfficiency);
219             printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
220         } break;
221         case kBayesian : {
222             unfoldedJetSpectrumIn = UnfoldSpectrumBayesian(       // do the inplane unfolding
223                 measuredJetSpectrumIn,
224                 resizedResponseIn,
225                 kinematicEfficiencyIn,
226                 measuredJetSpectrumTrueBinsIn,
227                 TString("in"),
228                 jetFindingEfficiency);
229             printf(" > Spectrum (in plane) unfolded using kBayesian unfolding < \n");
230         } break;
231         case kBayesianAli : {
232             unfoldedJetSpectrumIn = UnfoldSpectrumBayesianAli(       // do the inplane unfolding
233                 measuredJetSpectrumIn,
234                 resizedResponseIn,
235                 kinematicEfficiencyIn,
236                 measuredJetSpectrumTrueBinsIn,
237                 TString("in"),
238                 jetFindingEfficiency);
239             printf(" > Spectrum (in plane) unfolded using kBayesianAli unfolding < \n");
240         } break;
241         case kSVD : {
242             unfoldedJetSpectrumIn = UnfoldSpectrumSVD(       // do the inplane unfolding
243                 measuredJetSpectrumIn,
244                 resizedResponseIn,
245                 kinematicEfficiencyIn,
246                 measuredJetSpectrumTrueBinsIn,
247                 TString("in"),
248                 jetFindingEfficiency);
249             printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
250         } break;
251         case kNone : {  // do nothing 
252             resizedResponseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
253             unfoldedJetSpectrumIn = ProtectHeap(measuredJetSpectrumIn, kTRUE, TString("in"));
254         } break;
255         
256         default : {
257             printf(" > Selected unfolding method is not implemented yet ! \n");
258             return;
259         }
260     }
261     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
262     resizedResponseIn->SetXTitle("p_{T}^{true} [GeV/c]");
263     resizedResponseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
264     resizedResponseIn = ProtectHeap(resizedResponseIn);
265     resizedResponseIn->Write();
266     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
267     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
268     kinematicEfficiencyIn->Write();
269     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
270     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
271     fDetectorResponse->Write();
272     // optional histograms
273     if(fSaveFull) {
274         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
275         fSpectrumIn->Write();
276         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
277         fDptInDist->Write();
278         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
279         fDptIn->Write();
280         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
281         fFullResponseIn->Write();
282     }
283     fActiveDir->cd();
284     if(fInOutUnfolding) {
285         TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
286         dirOut->cd();
287         switch (fUnfoldingAlgorithm) {
288             case kChi2 : {
289                 unfoldedJetSpectrumOut = UnfoldSpectrumChi2(
290                     measuredJetSpectrumOut,
291                     resizedResponseOut,
292                     kinematicEfficiencyOut,
293                     measuredJetSpectrumTrueBinsOut,
294                     TString("out"),
295                     jetFindingEfficiency);
296                 printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
297             } break;
298             case kBayesian : {
299                 unfoldedJetSpectrumOut = UnfoldSpectrumBayesian(
300                     measuredJetSpectrumOut,
301                     resizedResponseOut,
302                     kinematicEfficiencyOut,
303                     measuredJetSpectrumTrueBinsOut,
304                     TString("out"),
305                     jetFindingEfficiency);
306                 printf(" > Spectrum (out of plane) unfolded using kBayesian < \n");
307             } break;
308             case kBayesianAli : {
309                 unfoldedJetSpectrumOut = UnfoldSpectrumBayesianAli(
310                     measuredJetSpectrumOut,
311                     resizedResponseOut,
312                     kinematicEfficiencyOut,
313                     measuredJetSpectrumTrueBinsOut,
314                     TString("out"),
315                     jetFindingEfficiency);
316                 printf(" > Spectrum (out of plane) unfolded using kBayesianAli < \n");
317             } break;
318             case kSVD : {
319                 unfoldedJetSpectrumOut = UnfoldSpectrumSVD(
320                     measuredJetSpectrumOut,
321                     resizedResponseOut,
322                     kinematicEfficiencyOut,
323                     measuredJetSpectrumTrueBinsOut,
324                     TString("out"),
325                     jetFindingEfficiency);
326                 printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
327             } break;
328             case kNone : {  // do nothing
329                 resizedResponseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
330                 unfoldedJetSpectrumOut = ProtectHeap(measuredJetSpectrumOut, kTRUE, TString("out"));
331             } break;
332             default : {
333                 printf(" > Selected unfolding method is not implemented yet ! \n");
334                 return;
335             }
336         }
337         resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
338         resizedResponseOut->SetXTitle("p_{T}^{true} [GeV/c]");
339         resizedResponseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
340         resizedResponseOut = ProtectHeap(resizedResponseOut);
341         resizedResponseOut->Write();
342         kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
343         kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
344         kinematicEfficiencyOut->Write();
345         fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
346         fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
347         fDetectorResponse->Write();
348         if(jetFindingEfficiency) jetFindingEfficiency->Write();
349         // optional histograms
350         if(fSaveFull) {
351             fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
352             fSpectrumOut->Write();
353             fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
354             fDptOutDist->Write();
355             fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
356             fDptOut->Write();
357             fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
358             fFullResponseOut->Write();
359         }
360
361         // write general output histograms to file
362         fActiveDir->cd();
363         if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
364             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
365             if(ratio) {
366                 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
367                 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
368                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
369                 ratio = ProtectHeap(ratio);
370                 ratio->Write();
371                 // write histo values to RMS files if both routines converged
372                 // input values are weighted by their uncertainty
373                 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
374                     if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
375                     if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
376                     if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
377                }
378             }
379             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
380             if(v2) {
381                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
382                 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
383                 v2->GetYaxis()->SetTitle("v_{2}");
384                 v2 = ProtectHeap(v2);
385                 v2->Write();
386             }
387         } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
388             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
389             if(ratio) {
390                 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
391                 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
392                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
393                 ratio = ProtectHeap(ratio);
394                 ratio->Write();
395             }
396             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
397              if(v2) {
398                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
399                 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
400                 v2->GetYaxis()->SetTitle("v_{2}");
401                 v2 = ProtectHeap(v2);
402                 v2->Write();
403             }
404         }
405     }   // end of if(fInOutUnfolding)
406     fDeltaPtDeltaPhi->Write();
407     fJetPtDeltaPhi->Write();
408     // save the current state of the unfolding object
409     SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
410 }
411 //_____________________________________________________________________________
412 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
413         const TH1D* measuredJetSpectrum,      // truncated raw jets (same binning as pt rec of response) 
414         const TH2D* resizedResponse,           // response matrix
415         const TH1D* kinematicEfficiency,      // kinematic efficiency
416         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
417         const TString suffix,                 // suffix (in or out of plane)
418         const TH1D* jetFindingEfficiency)     // jet finding efficiency (optional)
419 {
420     // unfold the spectrum using chi2 minimization
421
422     // step 0) setup the static members of AliUnfolding
423     ResetAliUnfolding();                // reset from previous iteration
424                                         // also deletes and re-creates the global TVirtualFitter
425     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
426     if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
427     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
428     if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
429     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
430     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
431
432     // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
433     // stay intact. a local copy of a histogram (which only exists in the scope of this function) is 
434     // denoted by the suffix 'Local'
435     
436     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
437     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
438     // unfolded local will be filled with the result of the unfolding
439     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
440
441     // full response matrix and kinematic efficiency
442     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
443     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
444
445     // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
446     TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
447     // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
448     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
449
450     // step 2) start the unfolding
451     Int_t status(-1), i(0);
452     while(status < 0 && i < 100) {
453         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
454         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the 
455         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
456         status = AliUnfolding::Unfold(
457                 resizedResponseLocal,           // response matrix
458                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
459                 measuredJetSpectrumLocal,              // measured spectrum
460                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
461                 unfoldedLocal);                 // results
462         // status holds the minuit fit status (where 0 means convergence)
463         i++;
464     }
465     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
466     if(status == 0 && gMinuit->fISW[1] == 3) {
467         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
468         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
469         if(gMinuit) gMinuit->Command("SET COV");
470         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
471         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
472         pearson->Print();
473         TH2D *hPearson(new TH2D(*pearson));
474         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
475         hPearson = ProtectHeap(hPearson);
476         hPearson->Write();
477     } else status = -1; 
478
479     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
480     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
481     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
482     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
483     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
484     if(ratio) {
485         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
486         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
487         ratio = ProtectHeap(ratio);
488         ratio->Write();
489     }
490
491     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap, 
492     // objects are cloned using 'ProtectHeap()'
493     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
494     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
495     measuredJetSpectrumLocal->Write(); 
496
497     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
498     resizedResponseLocal->Write();
499
500     unfoldedLocal = ProtectHeap(unfoldedLocal);
501     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
502     unfoldedLocal->Write();
503
504     foldedLocal = ProtectHeap(foldedLocal);
505     foldedLocal->Write();
506     
507     priorLocal = ProtectHeap(priorLocal);
508     priorLocal->Write();
509
510     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
511     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));
512     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
513     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
514     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
515     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
516     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
517     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
518     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
519     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
520     fitStatus->Write();
521     return unfoldedLocal;
522 }
523 //_____________________________________________________________________________
524 TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
525         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
526         const TH2D* resizedResponse,                   // full response matrix, normalized in slides of pt true
527         const TH1D* kinematicEfficiency,              // kinematic efficiency
528         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
529         const TString suffix,                         // suffix (in, out)
530         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
531 {
532
533     TH1D* priorLocal( GetPrior(
534         measuredJetSpectrum,              // jet pt in pt rec bins 
535         resizedResponse,                  // full response matrix, normalized in slides of pt true
536         kinematicEfficiency,              // kinematic efficiency
537         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
538         suffix,                           // suffix (in, out)
539         jetFindingEfficiency));           // jet finding efficiency (optional)
540     if(!priorLocal) {
541         printf(" > couldn't find prior ! < \n");
542         return 0x0;
543     } else printf(" 1) retrieved prior \n");
544
545     // go back to the 'root' directory of this instance of the SVD unfolding routine
546     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
547
548     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549     // measured jets in pt rec binning
550     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
551     // local copie of the response matrix
552     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
553     // local copy of response matrix, all true slides normalized to 1 
554     // this response matrix will eventually be used in the re-folding routine
555     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
556     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
557     // kinematic efficiency
558     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
559     // place holder histos
560     TH1D *unfoldedLocalSVD(0x0);
561     TH1D *foldedLocalSVD(0x0);
562     cout << " 2) setup necessary input " << endl;
563     // 3) configure routine
564     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
565     cout << " step 3) configured routine " << endl;
566
567     // 4) get transpose matrices
568     // a) get the transpose of the full response matrix
569     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
570     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
571     // normalize it with the prior. this will ensure that high statistics bins will constrain the
572     // end result most strenuously than bins with limited number of counts
573     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
574     cout << " 4) retrieved first transpose matrix " << endl;
575  
576     // 5) get response for SVD unfolding
577     RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
578     cout << " 5) retrieved roo unfold response object " << endl;
579
580     // 6) actualy unfolding loop
581     RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
582     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
583     // correct the spectrum for the kinematic efficiency
584     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
585
586     // get the pearson coefficients from the covariance matrix
587     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
588     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
589     if(pearson) {
590         TH2D* hPearson(new TH2D(*pearson));
591         pearson->Print();
592         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
593         hPearson = ProtectHeap(hPearson);
594         hPearson->Write();
595     } else return 0x0;       // return if unfolding didn't converge
596
597     // plot singular values and d_i vector
598     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
599     TH1* hSVal(svdUnfold->GetSV());
600     TH1D* hdi(svdUnfold->GetD());
601     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
602     hSVal->SetXTitle("singular values");
603     hSVal->Write();
604     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
605     hdi->SetXTitle("|d_{i}^{kreg}|");
606     hdi->Write();
607     cout << " plotted singular values and d_i vector " << endl;
608
609     // 7) refold the unfolded spectrum
610     foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
611     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
612     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
613     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
614     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
615     ratio->Write();
616     cout << " 7) refolded the unfolded spectrum " << endl;
617
618     // write the measured, unfolded and re-folded spectra to the output directory
619     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
620     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
621     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
622     measuredJetSpectrumLocal->Write(); // input spectrum
623     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
624     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
625     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
626     unfoldedLocalSVD->Write();  // unfolded spectrum
627     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
628     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
629     foldedLocalSVD->Write();    // re-folded spectrum
630
631     // save more general bookkeeeping histograms to the output directory
632     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
633     responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
634     responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
635     responseMatrixLocalTransposePrior->Write();
636     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
637     priorLocal->SetXTitle("p_{t} [GeV/c]");
638     priorLocal = ProtectHeap(priorLocal);
639     priorLocal->Write();
640     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
641     resizedResponseLocalNorm->Write();
642
643     // save some info 
644     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));
645     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
646     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
647     fitStatus->Write();
648
649     return unfoldedLocalSVD;
650 }
651 //_____________________________________________________________________________
652 TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
653         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
654         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
655         const TH1D* kinematicEfficiency,              // kinematic efficiency
656         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
657         const TString suffix,                         // suffix (in, out)
658         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
659 {
660     // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
661     // FIXME careful, not tested yet ! (06122013) FIXME
662
663     // step 0) setup the static members of AliUnfolding
664     ResetAliUnfolding();                // reset from previous iteration
665                                         // also deletes and re-creates the global TVirtualFitter
666     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
667     if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
668     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
669     else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
670     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
671     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
672
673     // 1) get a prior for unfolding and clone all the input histograms
674     TH1D* priorLocal( GetPrior(
675         measuredJetSpectrum,              // jet pt in pt rec bins 
676         resizedResponse,                  // full response matrix, normalized in slides of pt true
677         kinematicEfficiency,              // kinematic efficiency
678         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
679         suffix,                           // suffix (in, out)
680         jetFindingEfficiency));           // jet finding efficiency (optional)
681     if(!priorLocal) {
682         printf(" > couldn't find prior ! < \n");
683         return 0x0;
684     } else printf(" 1) retrieved prior \n");
685     // switch back to root dir of this unfolding procedure
686     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
687
688     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
689     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
690     // unfolded local will be filled with the result of the unfolding
691     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
692
693     // full response matrix and kinematic efficiency
694     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
695     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
696
697     // step 2) start the unfolding
698     Int_t status(-1), i(0);
699     while(status < 0 && i < 100) {
700         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
701         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the 
702         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
703         status = AliUnfolding::Unfold(
704                 resizedResponseLocal,           // response matrix
705                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
706                 measuredJetSpectrumLocal,              // measured spectrum
707                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
708                 unfoldedLocal);                 // results
709         // status holds the minuit fit status (where 0 means convergence)
710         i++;
711     }
712     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
713     if(status == 0 && gMinuit->fISW[1] == 3) {
714         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
715         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
716         if(gMinuit) gMinuit->Command("SET COV");
717         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
718         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
719         pearson->Print();
720         TH2D *hPearson(new TH2D(*pearson));
721         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
722         hPearson = ProtectHeap(hPearson);
723         hPearson->Write();
724     } else status = -1; 
725
726     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
727     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
728     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
729     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
730     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
731     if(ratio) {
732         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
733         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
734         ratio = ProtectHeap(ratio);
735         ratio->Write();
736     }
737
738     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap, 
739     // objects are cloned using 'ProtectHeap()'
740     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
741     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
742     measuredJetSpectrumLocal->Write(); 
743
744     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
745     resizedResponseLocal->Write();
746
747     unfoldedLocal = ProtectHeap(unfoldedLocal);
748     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
749     unfoldedLocal->Write();
750
751     foldedLocal = ProtectHeap(foldedLocal);
752     foldedLocal->Write();
753     
754     priorLocal = ProtectHeap(priorLocal);
755     priorLocal->Write();
756
757     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
758     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));
759     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
760     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
761     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
762     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
763     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
764     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
765     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
766     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
767     fitStatus->Write();
768     return unfoldedLocal;
769 }
770 //_____________________________________________________________________________
771 TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
772         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
773         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
774         const TH1D* kinematicEfficiency,              // kinematic efficiency
775         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
776         const TString suffix,                         // suffix (in, out)
777         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
778 {
779     // use bayesian unfolding from the RooUnfold package to unfold jet spectra
780     
781     // 1) get a prior for unfolding.
782     TH1D* priorLocal( GetPrior(
783         measuredJetSpectrum,              // jet pt in pt rec bins 
784         resizedResponse,                  // full response matrix, normalized in slides of pt true
785         kinematicEfficiency,              // kinematic efficiency
786         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
787         suffix,                           // suffix (in, out)
788         jetFindingEfficiency));           // jet finding efficiency (optional)
789     if(!priorLocal) {
790         printf(" > couldn't find prior ! < \n");
791         return 0x0;
792     } else printf(" 1) retrieved prior \n");
793     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
794
795     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
796     // measured jets in pt rec binning
797     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
798     // local copie of the response matrix
799     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
800     // local copy of response matrix, all true slides normalized to 1 
801     // this response matrix will eventually be used in the re-folding routine
802     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
803     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
804     // kinematic efficiency
805     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
806     // place holder histos
807     TH1D *unfoldedLocalBayes(0x0);
808     TH1D *foldedLocalBayes(0x0);
809     cout << " 2) setup necessary input " << endl;
810     // 4) get transpose matrices
811     // a) get the transpose of the full response matrix
812     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
813     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
814     // normalize it with the prior. this will ensure that high statistics bins will constrain the
815     // end result most strenuously than bins with limited number of counts
816     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
817     // 3) get response for Bayesian unfolding
818     RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
819
820     // 4) actualy unfolding loop
821     RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
822     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
823     unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
824     // correct the spectrum for the kinematic efficiency
825     unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
826     // get the pearson coefficients from the covariance matrix
827     TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
828     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
829     if(pearson) {
830         TH2D* hPearson(new TH2D(*pearson));
831         pearson->Print();
832         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
833         hPearson = ProtectHeap(hPearson);
834         hPearson->Write();
835     } else return 0x0;       // return if unfolding didn't converge
836
837     // 5) refold the unfolded spectrum
838     foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
839     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio  measured / re-folded", kTRUE));
840     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
841     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
842     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
843     ratio->Write();
844     cout << " 7) refolded the unfolded spectrum " << endl;
845
846     // write the measured, unfolded and re-folded spectra to the output directory
847     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
848     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
849     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
850     measuredJetSpectrumLocal->Write(); // input spectrum
851     unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
852     unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
853     if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
854     unfoldedLocalBayes->Write();  // unfolded spectrum
855     foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
856     foldedLocalBayes = ProtectHeap(foldedLocalBayes);
857     foldedLocalBayes->Write();    // re-folded spectrum
858
859     // save more general bookkeeeping histograms to the output directory
860     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
861     responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
862     responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
863     responseMatrixLocalTransposePrior->Write();
864     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
865     priorLocal->SetXTitle("p_{t} [GeV/c]");
866     priorLocal = ProtectHeap(priorLocal);
867     priorLocal->Write();
868     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
869     resizedResponseLocalNorm->Write();
870
871     // save some info 
872     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));
873     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
874     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
875     fitStatus->Write();
876
877     return  unfoldedLocalBayes; 
878 }
879 //_____________________________________________________________________________
880 Bool_t AliJetFlowTools::PrepareForUnfolding()
881 {
882     // prepare for unfolding
883     if(fRawInputProvided) return kTRUE;
884     if(!fInputList) {
885         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
886         return kFALSE;
887     }
888     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
889     // check if the pt bin for true and rec have been set
890     if(!fBinsTrue || !fBinsRec) {
891         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
892         return kFALSE;
893     }
894     if(!fRMSSpectrumIn && fInOutUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
895                           // procedures, these profiles will be nonsensical, user is responsible
896         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
897         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
898         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
899     }
900     if(!fTrainPower) {
901         // clear minuit state to avoid constraining the fit with the results of the previous iteration
902         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
903     }
904     // extract the spectra
905     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
906     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
907     if(!fJetPtDeltaPhi) {
908         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
909         return kFALSE;
910     }
911     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
912     // in plane spectrum
913     if(!fInOutUnfolding) {
914         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
915         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
916     } else {
917         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
918         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
919         fSpectrumIn = ProtectHeap(fSpectrumIn);
920         // out of plane spectrum
921         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
922         fSpectrumOut = ProtectHeap(fSpectrumOut);
923     }
924     // normalize spectra to event count if requested
925     if(fNormalizeSpectra) {
926         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
927         if(!rho) return 0x0;
928         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
929         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
930         if(fEventCount > 0) {
931             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
932             fSpectrumOut->Sumw2();
933             Double_t pt(0);            
934             Double_t error(0); // lots of issues with the errors here ...
935             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
936                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
937                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
938                 fSpectrumIn->SetBinContent(1+i, pt);
939                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
940                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
941                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
942             }
943             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
944                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
945                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
946                 fSpectrumOut->SetBinContent(1+i, pt);
947                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
948                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
949                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
950             }
951         }
952         if(normalizeToFullSpectrum) fEventCount = -1;
953     }
954     // extract the delta pt matrices
955     TString deltaptName(Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin));
956     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
957     if(!fDeltaPtDeltaPhi) {
958         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
959         printf(" > may be ok, depending no what you want to do < \n");
960         fRefreshInput = kTRUE;
961         return kTRUE;
962     }
963     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
964     // in plane delta pt distribution
965     if(!fInOutUnfolding) {
966         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
967         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
968     } else {
969         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
970         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
971         // out of plane delta pt distribution
972         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
973         fDptInDist = ProtectHeap(fDptInDist);
974         fDptOutDist = ProtectHeap(fDptOutDist);
975         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
976     }
977
978     // create a rec - true smeared response matrix
979     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
980     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
981         Bool_t skip = kFALSE;
982         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
983             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
984             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
985         }
986     }
987     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
988     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
989         Bool_t skip = kFALSE;
990         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
991             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
992             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
993         }
994     }
995     fDptIn = new TH2D(*rfIn);
996     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
997     fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
998     fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
999     fDptIn = ProtectHeap(fDptIn);
1000     fDptOut = new TH2D(*rfOut);
1001     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1002     fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1003     fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1004     fDptOut = ProtectHeap(fDptOut);
1005     
1006     fRefreshInput = kTRUE;     // force cloning of the input
1007     return kTRUE;
1008 }
1009 //_____________________________________________________________________________
1010 TH1D* AliJetFlowTools::GetPrior(
1011         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins 
1012         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
1013         const TH1D* kinematicEfficiency,              // kinematic efficiency
1014         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
1015         const TString suffix,                         // suffix (in, out)
1016         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
1017 {
1018     // 1) get a prior for unfolding. 
1019     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1020     TH1D* unfolded(0x0);
1021     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1022     dirOut->cd();
1023     switch (fPrior) {    // select the prior for unfolding
1024         case kPriorChi2 : {
1025             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
1026                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original binning
1027                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
1028                 TArrayD* tempArrayRec(fBinsRec);
1029                 fBinsRec = fBinsRecPrior;
1030                 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1031                 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1032                 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1033                 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1034                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1035                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1036                 unfolded= UnfoldSpectrumChi2(
1037                             measuredJetSpectrumChi2,
1038                             resizedResponseChi2,
1039                             kinematicEfficiencyChi2,
1040                             measuredJetSpectrumTrueBinsChi2,    // prior for chi2 unfolding (measured)
1041                             TString(Form("prior_%s", suffix.Data())));
1042                if(!unfolded) {
1043                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1044                     printf("               probably Chi2 unfolding did not converge < \n");
1045                     return 0x0;
1046                 }
1047                 fBinsTrue = tempArrayTrue;  // reset bins borders
1048                 fBinsRec = tempArrayRec;
1049                 // if the chi2 prior has a different binning, rebin to the true binning for the  unfolding
1050                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
1051             } else {
1052                 unfolded = UnfoldSpectrumChi2(
1053                             measuredJetSpectrum,
1054                             resizedResponse,
1055                             kinematicEfficiency,
1056                             measuredJetSpectrumTrueBins,        // prior for chi2 unfolding (measured)
1057                             TString(Form("prior_%s", suffix.Data())));
1058                 if(!unfolded) {
1059                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1060                     printf("               probably Chi2 unfolding did not converge < \n");
1061                     return 0x0;
1062                 }
1063             }
1064             break;
1065         }
1066         case kPriorMeasured : { 
1067             unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
1068         }
1069         default : break;
1070     }
1071     // it can be important that the prior is smooth, this can be achieved by 
1072     // extrapolating the spectrum with a fitted power law when bins are sparsely filed 
1073     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1074     TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1075     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1076     return priorLocal;
1077 }
1078 //_____________________________________________________________________________
1079 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1080     // resize the x-axis of a th1d
1081     if(!histo) {
1082         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1083         return NULL;
1084     } 
1085     // see how many bins we need to copy
1086     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);
1087     // low is the bin number of the first new bin
1088     Int_t l = histo->GetXaxis()->FindBin(low);
1089     // set the values
1090     Double_t _x(0), _xx(0);
1091     for(Int_t i(0); i < up-low; i++) {
1092         _x = histo->GetBinContent(l+i);
1093         _xx=histo->GetBinError(l+i);
1094         resized->SetBinContent(i+1, _x);
1095         resized->SetBinError(i+1, _xx);
1096     }
1097     return resized;
1098 }
1099 //_____________________________________________________________________________
1100 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1101     // resize the y-axis of a th2d
1102     if(!histo) {
1103         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1104         return NULL;
1105     } 
1106     // see how many bins we need to copy
1107     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());
1108     // assume only the y-axis has changed
1109     // low is the bin number of the first new bin
1110     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1111     // set the values
1112     Double_t _x(0), _xx(0);
1113     for(Int_t i(0); i < x->GetSize(); i++) {
1114         for(Int_t j(0); j < y->GetSize(); j++) {
1115             _x = histo->GetBinContent(i, low+j);
1116             _xx=histo->GetBinError(i, low+1+j);
1117             resized->SetBinContent(i, j, _x);
1118             resized->SetBinError(i, j, _xx);
1119         }
1120     }
1121     return resized;
1122 }
1123 //_____________________________________________________________________________
1124 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1125     // general method to normalize all vertical slices of a th2 to unity
1126     // i.e. get a probability matrix
1127     if(!histo) {
1128         printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1129         return NULL;
1130     }
1131     Int_t binsX = histo->GetXaxis()->GetNbins();
1132     Int_t binsY = histo->GetYaxis()->GetNbins();
1133     
1134     // normalize all slices in x
1135     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1136         Double_t weight = 0;
1137         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1138             weight+=histo->GetBinContent(i+1, j+1);
1139         }       // now we know the total weight
1140         for(Int_t j(0); j < binsY; j++) {
1141             if (weight <= 0 ) continue;
1142             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1143             if(noError) histo->SetBinError(  1+i, j+1, 0.);
1144             else histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1145         }
1146     }
1147     return histo;
1148 }
1149 //_____________________________________________________________________________
1150 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1151     // return a TH1D with the supplied histogram rebinned to the supplied bins
1152     // the returned histogram is new, the original is deleted from the heap if kill is true
1153     if(!histo || !bins) {
1154         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1155         return NULL;
1156     }
1157     // create the output histo
1158     TString name = histo->GetName();
1159     name+="_template";
1160     name+=suffix;
1161     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1162     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1163         // loop over the bins of the old histo and fill the new one with its data
1164         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1165     }
1166     if(kill) delete histo;
1167     return rebinned;
1168 }
1169 //_____________________________________________________________________________
1170 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1171     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1172     if(!fResponseMaker || !binsTrue || !binsRec) {
1173         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1174         return 0x0;
1175     }
1176     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1177     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1178 }
1179 //_____________________________________________________________________________
1180 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1181 {
1182     // multiply two matrices
1183     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1184     TH2D* c = (TH2D*)a->Clone("c");
1185     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1186         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1187             Double_t val = 0;
1188             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1189                 Int_t y2 = x1;
1190                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1191             }
1192             c->SetBinContent(x2, y1, val);
1193             c->SetBinError(x2, y1, 0.);
1194         }
1195     }
1196     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1197     return c;
1198 }
1199 //_____________________________________________________________________________
1200 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale) 
1201 {
1202     // normalize a th1d to a certain scale
1203     histo->Sumw2();
1204     Double_t integral = histo->Integral()*scale;
1205     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1206     else if (scale != 1.) histo->Scale(1./scale, "width");
1207     else printf(" > Histogram integral < 0, cannot normalize \n");
1208     return histo;
1209 }
1210 //_____________________________________________________________________________
1211 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix) 
1212 {
1213     // Calculate pearson coefficients from covariance matrix
1214     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1215     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1216     Double_t pearson(0.);
1217     if(nrows==0 && ncols==0) return 0x0;
1218     for(Int_t row = 0; row < nrows; row++) {
1219         for(Int_t col = 0; col<ncols; col++) {
1220         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1221         (*pearsonCoefficients)(row,col) = pearson;
1222         }
1223     }
1224     return pearsonCoefficients;
1225 }
1226 //_____________________________________________________________________________
1227 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1228     // smoothen the spectrum using a user defined function
1229     // returns a clone of the original spectrum if fitting failed
1230     // if kill is kTRUE the input spectrum will be deleted from the heap
1231     // if 'count' is selected, bins are filled with integers (necessary if the 
1232     // histogram is interpreted in a routine which accepts only counts)
1233     if(!spectrum || !function) return 0x0;
1234     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1235     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1236     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1237     TFitResultPtr r = temp->Fit(function, "WLIS", "", min, max);
1238     if((int)r == 0) {   // MINUIT status
1239         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1240             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1241                 if(counts) temp->SetBinContent(i, (int)function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1242                 else temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1243                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1244             }
1245         }
1246     }
1247     if(kill) delete spectrum;
1248     return temp;
1249 }
1250 //_____________________________________________________________________________
1251 void AliJetFlowTools::Style(TCanvas* c, TString style)
1252 {
1253     // set a default style for a canvas
1254     if(!strcmp(style.Data(), "PEARSON")) {
1255         printf(" > style PEARSON canvas < \n");
1256         gStyle->SetOptStat(0);
1257         c->SetGridx();
1258         c->SetGridy();
1259         c->SetTicks();
1260         return;
1261     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1262         printf(" > style SPECTRUM canvas < \n");
1263         gStyle->SetOptStat(0);
1264         c->SetLogy();
1265         c->SetGridx();
1266         c->SetGridy();
1267         c->SetTicks();
1268         return;
1269     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1270 }
1271 //_____________________________________________________________________________
1272 void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1273 {
1274     // set a default style for a canvas
1275     if(!strcmp(style.Data(), "PEARSON")) {
1276         printf(" > style PEARSON pad < \n");
1277         gStyle->SetOptStat(0);
1278         c->SetGridx();
1279         c->SetGridy();
1280         c->SetTicks();
1281         return;
1282     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1283         printf(" > style SPECTRUM pad < \n");
1284         gStyle->SetOptStat(0);
1285         c->SetLogy();
1286         c->SetGridx();
1287         c->SetGridy();
1288         c->SetTicks();
1289         return;
1290     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1291 }
1292 //_____________________________________________________________________________
1293 void AliJetFlowTools::Style(TLegend* l)
1294 {
1295     // set a default style for a legend
1296     l->SetTextSize(.06);
1297     l->SetFillColor(0);
1298 }
1299 //_____________________________________________________________________________
1300 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1301 {
1302     // style a histo
1303     h->SetLineColor(col);
1304     h->SetMarkerColor(col);
1305     h->SetLineWidth(2.);
1306     h->SetMarkerSize(2.);
1307     h->SetTitleSize(0.06);
1308     h->GetYaxis()->SetLabelSize(0.06);
1309     h->GetXaxis()->SetLabelSize(0.06);
1310     h->GetYaxis()->SetTitleSize(0.06);
1311     h->GetXaxis()->SetTitleSize(0.06);
1312     h->GetYaxis()->SetTitleOffset(1.);
1313     h->GetXaxis()->SetTitleOffset(.9);
1314     switch (type) {
1315         case kInPlaneSpectrum : {
1316             h->SetTitle("IN PLANE");
1317             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1318             h->GetYaxis()->SetTitle("counts");
1319         } break;
1320         case kOutPlaneSpectrum : {
1321             h->SetTitle("OUT OF PLANE");
1322             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1323             h->GetYaxis()->SetTitle("counts");
1324        } break;
1325        case kUnfoldedSpectrum : {
1326             h->SetTitle("UNFOLDED");
1327             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1328             h->GetYaxis()->SetTitle("counts");
1329        } break;
1330        case kFoldedSpectrum : {
1331             h->SetTitle("FOLDED");
1332             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1333             h->GetYaxis()->SetTitle("counts");
1334        } break;
1335        case kMeasuredSpectrum : {
1336             h->SetTitle("MEASURED");
1337             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1338             h->GetYaxis()->SetTitle("counts");
1339        } break;
1340        default : break;
1341     }
1342 }
1343 //_____________________________________________________________________________
1344 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1345 {
1346     // style a histo
1347     h->SetLineColor(col);
1348     h->SetMarkerColor(col);
1349     h->SetLineWidth(2.);
1350     h->SetMarkerSize(2.);
1351     h->GetYaxis()->SetLabelSize(0.06);
1352     h->GetXaxis()->SetLabelSize(0.06);
1353     h->GetYaxis()->SetTitleSize(0.06);
1354     h->GetXaxis()->SetTitleSize(0.06);
1355     h->GetYaxis()->SetTitleOffset(1.);
1356     h->GetXaxis()->SetTitleOffset(.9);
1357     switch (type) {
1358         case kInPlaneSpectrum : {
1359             h->SetTitle("IN PLANE");
1360             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1361             h->GetYaxis()->SetTitle("counts");
1362         } break;
1363         case kOutPlaneSpectrum : {
1364             h->SetTitle("OUT OF PLANE");
1365             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1366             h->GetYaxis()->SetTitle("counts");
1367        } break;
1368        case kUnfoldedSpectrum : {
1369             h->SetTitle("UNFOLDED");
1370             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1371             h->GetYaxis()->SetTitle("counts");
1372        } break;
1373        case kFoldedSpectrum : {
1374             h->SetTitle("FOLDED");
1375             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1376             h->GetYaxis()->SetTitle("counts");
1377        } break;
1378        case kMeasuredSpectrum : {
1379             h->SetTitle("MEASURED");
1380             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1381             h->GetYaxis()->SetTitle("counts");
1382        } break;
1383        default : break;
1384     }
1385 }
1386 //_____________________________________________________________________________
1387 void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t columns) const
1388 {
1389    // go through the output file and perform post processing routines
1390    // can either be performed in one go with the unfolding, or at a later stage
1391    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
1392    TFile readMe(in.Data(), "READ");     // open file read-only
1393    if(readMe.IsZombie()) {
1394        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1395        return;
1396    }
1397    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
1398    readMe.ls();
1399    TList* listOfKeys((TList*)readMe.GetListOfKeys());
1400    if(!listOfKeys) {
1401        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1402        return;
1403    }
1404    // prepare necessary canvasses
1405    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
1406    TCanvas* canvasOut(0x0);
1407    if(fInOutUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
1408    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
1409    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
1410    if(fInOutUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
1411    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn")); 
1412    TCanvas* canvasSpectraOut(0x0);
1413    if(fInOutUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
1414    TCanvas* canvasRatio(0x0);
1415    if(fInOutUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
1416    TCanvas* canvasV2(0x0);
1417    if(fInOutUnfolding) canvasV2 = new TCanvas("V2", "V2");
1418    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1419    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
1420    TCanvas* canvasMasterOut(0x0);
1421    if(fInOutUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
1422    (fInOutUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
1423    TDirectoryFile* defDir(0x0);
1424    
1425    // get an estimate of the number of outputs and find the default set
1426    Int_t cacheMe(0);
1427    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1428        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1429            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1430            cacheMe++;
1431        }
1432    }
1433    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%4)>0));
1434    canvasIn->Divide(columns, rows);
1435    if(canvasOut) canvasOut->Divide(columns, rows);
1436    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1437    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1438    canvasSpectraIn->Divide(columns, rows);
1439    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
1440    if(canvasRatio) canvasRatio->Divide(columns, rows);
1441    if(canvasV2) canvasV2->Divide(columns, rows);
1442
1443    canvasMasterIn->Divide(columns, rows);
1444    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
1445    // extract the default output 
1446    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
1447    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
1448    if(defDir) {
1449        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1450        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
1451        if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1452        if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
1453        printf(" > succesfully extracted default results < \n");
1454    }
1455  
1456    // loop through the directories, only plot the graphs if the deconvolution converged
1457    TDirectoryFile* tempDir(0x0); 
1458    TDirectoryFile* tempIn(0x0);
1459    TDirectoryFile*  tempOut(0x0);
1460    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
1461        // read the directory by its name
1462        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1463        if(!tempDir) continue;
1464        TString dirName(tempDir->GetName());
1465        // try to read the in- and out of plane subdirs
1466        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1467        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1468        j++;
1469        if(tempIn) { 
1470            // to see if the unfolding converged try to extract the pearson coefficients
1471            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1472            if(pIn) {
1473                printf(" - %s in plane converged \n", dirName.Data());
1474                canvasIn->cd(j);
1475                Style(gPad, "PEARSON");
1476                pIn->DrawCopy("colz");
1477                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1478                if(rIn) {
1479                    Style(rIn);
1480                    printf(" > found RatioRefoldedMeasured < \n");
1481                    canvasRatioMeasuredRefoldedIn->cd(j);
1482                    rIn->SetFillColor(kRed);
1483                    rIn->Draw("ap");
1484                }
1485                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1486                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1487                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1488                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1489                if(dvector && avalue && rm && eff) {
1490                    Style(dvector);
1491                    Style(avalue);
1492                    Style(rm);
1493                    Style(eff);
1494                    canvasMISC->cd(1);
1495                    Style(gPad, "SPECTRUM");
1496                    dvector->DrawCopy();
1497                    canvasMISC->cd(2);
1498                    Style(gPad, "SPECTRUM");
1499                    avalue->DrawCopy();
1500                    canvasMISC->cd(3);
1501                    Style(gPad, "PEARSON");
1502                    rm->DrawCopy("colz");
1503                    canvasMISC->cd(4);
1504                    eff->DrawCopy();
1505                } else if(rm && eff) {
1506                    Style(rm);
1507                    Style(eff);
1508                    canvasMISC->cd(3);
1509                    Style(gPad, "PEARSON");
1510                    rm->DrawCopy("colz");
1511                    canvasMISC->cd(4);
1512                    eff->DrawCopy();
1513                }
1514            }
1515            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1516            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1517            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1518            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1519                if(defaultUnfoldedJetSpectrumIn) {
1520                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
1521                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
1522                    temp->Divide(unfoldedSpectrum);
1523                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1524                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1525                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1526                    canvasMasterIn->cd(j);
1527                    temp->GetYaxis()->SetRangeUser(0., 2);
1528                    temp->DrawCopy();
1529                }
1530                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
1531                canvasSpectraIn->cd(j);
1532                Style(gPad);
1533                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1534                unfoldedSpectrum->DrawCopy();
1535                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1536                inputSpectrum->DrawCopy("same");
1537                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1538                refoldedSpectrum->DrawCopy("same");
1539                TLegend* l(AddLegend(gPad));
1540                Style(l);
1541                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1542                    Float_t chi(fitStatus->GetBinContent(1));
1543                    Float_t pen(fitStatus->GetBinContent(2));
1544                    Int_t dof((int)fitStatus->GetBinContent(3));
1545                    Float_t beta(fitStatus->GetBinContent(4));
1546                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1547                } else if (fitStatus) { // only available in SVD fit
1548                    Int_t reg((int)fitStatus->GetBinContent(1));
1549                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1550                }
1551            }
1552        }
1553        if(tempOut) {
1554            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
1555            if(pOut) {
1556                printf(" - %s out of plane converged \n", dirName.Data());
1557                canvasOut->cd(j);
1558                Style(gPad, "PEARSON");
1559                pOut->DrawCopy("colz");
1560                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1561                if(rOut) {
1562                    Style(rOut);
1563                    printf(" > found RatioRefoldedMeasured < \n");
1564                    canvasRatioMeasuredRefoldedOut->cd(j);
1565                    rOut->SetFillColor(kRed);
1566                    rOut->Draw("ap");
1567                }
1568                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1569                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1570                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
1571                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
1572                if(dvector && avalue && rm && eff) {
1573                    Style(dvector);
1574                    Style(avalue);
1575                    Style(rm);
1576                    Style(eff);
1577                    canvasMISC->cd(5);
1578                    Style(gPad, "SPECTRUM");
1579                    dvector->DrawCopy();
1580                    canvasMISC->cd(6);
1581                    Style(gPad, "SPECTRUM");
1582                    avalue->DrawCopy();
1583                    canvasMISC->cd(7);
1584                    Style(gPad, "PEARSON");
1585                    rm->DrawCopy("colz");
1586                    canvasMISC->cd(8);
1587                    eff->DrawCopy();
1588                } else if(rm && eff) {
1589                    Style(rm);
1590                    Style(eff);
1591                    canvasMISC->cd(7);
1592                    Style(gPad, "PEARSON");
1593                    rm->DrawCopy("colz");
1594                    canvasMISC->cd(8);
1595                    eff->DrawCopy();
1596                }
1597            }
1598            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
1599            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
1600            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
1601            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1602                if(defaultUnfoldedJetSpectrumOut) {
1603                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
1604                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
1605                    temp->Divide(unfoldedSpectrum);
1606                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1607                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1608                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1609                    canvasMasterOut->cd(j);
1610                    temp->GetYaxis()->SetRangeUser(0., 2.);
1611                    temp->DrawCopy();
1612                }
1613                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
1614                canvasSpectraOut->cd(j);
1615                Style(gPad);
1616                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1617                unfoldedSpectrum->DrawCopy();
1618                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1619                inputSpectrum->DrawCopy("same");
1620                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1621                refoldedSpectrum->DrawCopy("same");
1622                TLegend* l(AddLegend(gPad));
1623                Style(l);
1624                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1625                    Float_t chi(fitStatus->GetBinContent(1));
1626                    Float_t pen(fitStatus->GetBinContent(2));
1627                    Int_t dof((int)fitStatus->GetBinContent(3));
1628                    Float_t beta(fitStatus->GetBinContent(4));
1629                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1630                } else if (fitStatus) { // only available in SVD fit
1631                    Int_t reg((int)fitStatus->GetBinContent(1));
1632                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1633                }
1634            }
1635        }
1636        if(canvasRatio && canvasV2) {
1637            canvasRatio->cd(j);
1638            TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
1639            if(ratioYield) {
1640                Style(ratioYield);
1641                ratioYield->GetYaxis()->SetRangeUser(-1., 3.);
1642                ratioYield->SetFillColor(kRed);
1643                ratioYield->Draw("ap");
1644            }
1645            canvasV2->cd(j);
1646            TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
1647            if(ratioV2) {
1648                Style(ratioV2);
1649                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
1650                ratioV2->SetFillColor(kRed);
1651                ratioV2->Draw("ap");
1652            }
1653        }
1654    }
1655    TFile output(out.Data(), "RECREATE");
1656    SavePadToPDF(canvasIn);
1657    canvasIn->Write();
1658    if(canvasOut) {
1659        SavePadToPDF(canvasOut);
1660        canvasOut->Write();
1661    }
1662    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
1663    canvasRatioMeasuredRefoldedIn->Write();
1664    if(canvasRatioMeasuredRefoldedOut) {
1665        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
1666        canvasRatioMeasuredRefoldedOut->Write();
1667    }
1668    SavePadToPDF(canvasSpectraIn);
1669    canvasSpectraIn->Write();
1670    if(canvasSpectraOut) {
1671        SavePadToPDF(canvasSpectraOut);
1672        canvasSpectraOut->Write();
1673        SavePadToPDF(canvasRatio);
1674        canvasRatio->Write();
1675        SavePadToPDF(canvasV2);
1676        canvasV2->Write();
1677    }
1678    SavePadToPDF(canvasMasterIn);
1679    canvasMasterIn->Write();
1680    if(canvasMasterOut) {
1681        SavePadToPDF(canvasMasterOut);
1682        canvasMasterOut->Write();
1683    }
1684    SavePadToPDF(canvasMISC);
1685    canvasMISC->Write();
1686    output.Write();
1687    output.Close();
1688 }
1689 //_____________________________________________________________________________
1690 Bool_t AliJetFlowTools::SetRawInput (
1691         TH2D* detectorResponse,  // detector response matrix
1692         TH1D* jetPtIn,           // in plane jet spectrum
1693         TH1D* jetPtOut,          // out of plane jet spectrum
1694         TH1D* dptIn,             // in plane delta pt distribution
1695         TH1D* dptOut,            // out of plane delta pt distribution
1696         Int_t eventCount) {
1697     // set input histograms manually
1698     fDetectorResponse   = detectorResponse;
1699     fSpectrumIn         = jetPtIn;
1700     fSpectrumOut        = jetPtOut;
1701     fDptInDist          = dptIn;
1702     fDptOutDist         = dptOut;
1703     fRawInputProvided   = kTRUE;
1704     // check if all data is provided
1705     if(!fDetectorResponse) {
1706         printf(" fDetectorResponse not found \n ");
1707         return kFALSE;
1708     }
1709     // check if the pt bin for true and rec have been set
1710     if(!fBinsTrue || !fBinsRec) {
1711         printf(" No true or rec bins set, please set binning ! \n");
1712         return kFALSE;
1713     }
1714     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
1715                           // procedures, these profiles will be nonsensical, user is responsible
1716         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1717         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1718         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1719     }
1720     // normalize spectra to event count if requested
1721     if(fNormalizeSpectra) {
1722         fEventCount = eventCount;
1723         if(fEventCount > 0) {
1724             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1725             fSpectrumOut->Sumw2();
1726             fSpectrumIn->Scale(1./((double)fEventCount));
1727             fSpectrumOut->Scale(1./((double)fEventCount));
1728         }
1729     }
1730     if(!fNormalizeSpectra && fEventCount > 0) {
1731         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1732         fSpectrumOut->Sumw2();
1733         fSpectrumIn->Scale(1./((double)fEventCount));
1734         fSpectrumOut->Scale(1./((double)fEventCount));
1735     }
1736     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
1737     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1738     fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1739     fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1740     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
1741     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1742     fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1743     fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1744     
1745     return kTRUE;
1746 }
1747 //_____________________________________________________________________________
1748 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax) 
1749 {
1750     // return ratio of h1 / h2
1751     // histograms can have different binning. errors are propagated as uncorrelated
1752     if(!(h1 && h2) ) {
1753         printf(" GetRatio called with NULL argument(s) \n ");
1754         return 0x0;
1755     }
1756     Int_t j(0);
1757     TGraphErrors *gr = new TGraphErrors();
1758     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1759     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1760     TH1* dud((TH1*)h1->Clone("dud"));
1761     if(!dud->Divide(h1, h2)) {
1762         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1763             binCent = h1->GetXaxis()->GetBinCenter(i);
1764             if(xmax > 0. && binCent > xmax) continue;
1765             j = h2->FindBin(binCent);
1766             binWidth = h1->GetXaxis()->GetBinWidth(i);
1767             if(h2->GetBinContent(j) > 0.) {
1768                 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
1769             /* original propagation of uncertainty
1770             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
1771             Double_t B = 0.;
1772             if(h2->GetBinError(j)>0.) {
1773                 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
1774                 error2 = A*A + B*B;
1775             } else error2 = A*A;        */
1776                 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
1777                 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
1778                 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
1779                 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
1780                 gr->SetPoint(gr->GetN(),binCent,ratio);
1781                 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1782             }
1783         }
1784     } else {
1785         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1786             binCent = dud->GetXaxis()->GetBinCenter(i);
1787             if(xmax > 0. && binCent > xmax) continue;
1788             binWidth = dud->GetXaxis()->GetBinWidth(i);
1789             gr->SetPoint(gr->GetN(),binCent,dud->GetBinContent(i));
1790             gr->SetPointError(gr->GetN()-1,0.5*binWidth,dud->GetBinError(i));
1791         }
1792     }
1793
1794     if(appendFit) {
1795         TF1* fit(new TF1("lin", "pol0", 10, 100));
1796         gr->Fit(fit);
1797     }
1798     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1799     if(dud) delete dud;
1800     return gr;
1801 }
1802 //_____________________________________________________________________________
1803 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name) 
1804 {
1805     // get v2 from difference of in plane, out of plane yield
1806     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
1807     // different binning is allowed but will mean that the error
1808     // propagation is unreliable
1809     // r is the event plane resolution for the chosen centrality
1810     if(!(h1 && h2) ) {
1811         printf(" GetV2 called with NULL argument(s) \n ");
1812         return 0x0;
1813     }
1814     Int_t j(0);
1815     TGraphErrors *gr = new TGraphErrors();
1816     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1817     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1818     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
1819     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1820         binCent = h1->GetXaxis()->GetBinCenter(i);
1821         j = h2->FindBin(binCent);
1822         binWidth = h1->GetXaxis()->GetBinWidth(i);
1823         if(h2->GetBinContent(j) > 0.) {
1824             in = h1->GetBinContent(i);
1825             ein = h1->GetBinError(i);
1826             out = h2->GetBinContent(j);
1827             eout = h2->GetBinError(j);
1828             ratio = pre*((in-out)/(in+out));
1829             error2 =((r*4.)/(TMath::Pi()))*((4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout);
1830             if(error2 > 0) error2 = TMath::Sqrt(error2);
1831             gr->SetPoint(gr->GetN(),binCent,ratio);
1832             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1833         }
1834     }
1835     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1836     return gr;
1837 }
1838 //_____________________________________________________________________________
1839 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
1840     // write object, if a unique identifier is given the object is cloned
1841     // and the clone is saved. setting kill to true will delete the original obect from the heap
1842     if(!object) {
1843         printf(" > WriteObject:: called with NULL arguments \n ");
1844         return;
1845     } else if(!strcmp("", suffix.Data())) object->Write();
1846     else {
1847         TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
1848         newObject->Write();
1849     }
1850     if(kill) delete object;
1851 }
1852 //_____________________________________________________________________________
1853 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
1854     // construt a delta pt response matrix from supplied dpt distribution
1855     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to 
1856     // do a weighted rebinning to a (coarser) dpt distribution
1857     // be careful with the binning of the dpt response: it should be equal to that
1858     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
1859     //
1860     // the response matrix will be square and have the same binning
1861     // (min, max and granularity) of the input histogram
1862     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
1863     Double_t _bins[bins+1];             // prepare array with bin borders
1864     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
1865     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
1866     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
1867     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
1868         Bool_t skip = kFALSE;
1869         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
1870             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
1871             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
1872         }
1873     }
1874     return res;
1875 }
1876 //_____________________________________________________________________________
1877 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1878     if(!binsTrue || !binsRec) {
1879         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
1880         return 0x0;
1881     }
1882     TString name(Form("unityResponse_%s", suffix.Data()));
1883     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
1884     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
1885         for(Int_t j(0); j < binsRec->GetSize(); j++) {
1886             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
1887         }
1888     }
1889     return unity;
1890 }
1891 //_____________________________________________________________________________
1892 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
1893     // save configuration parameters to histogram
1894     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
1895     summary->SetBinContent(1, fBetaIn);
1896     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
1897     summary->SetBinContent(2, fBetaOut);
1898     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
1899     summary->SetBinContent(3, fCentralityBin);
1900     summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
1901     summary->SetBinContent(4, (int)convergedIn);
1902     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
1903     summary->SetBinContent(5, (int)convergedOut);
1904     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
1905     summary->SetBinContent(6, (int)fAvoidRoundingError);
1906     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
1907     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
1908     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
1909     summary->SetBinContent(8, (int)fPrior);
1910     summary->GetXaxis()->SetBinLabel(8, "fPrior");
1911     summary->SetBinContent(9, fSVDRegIn);
1912     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
1913     summary->SetBinContent(10, fSVDRegOut);
1914     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
1915     summary->SetBinContent(11, (int)fSVDToy);
1916     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
1917     summary->SetBinContent(12, fJetRadius);
1918     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
1919     summary->SetBinContent(13, (int)fNormalizeSpectra);
1920     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
1921     summary->SetBinContent(14, (int)fSmoothenPrior);
1922     summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
1923     summary->SetBinContent(15, (int)fTestMode);
1924     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
1925     summary->SetBinContent(16, (int)fUseDetectorResponse);
1926     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
1927     summary->SetBinContent(17, fBayesianIterIn);
1928     summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
1929     summary->SetBinContent(18, fBayesianIterOut);
1930     summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
1931     summary->SetBinContent(19, fBayesianSmoothIn);
1932     summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
1933     summary->SetBinContent(20, fBayesianSmoothOut);
1934     summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
1935 }
1936 //_____________________________________________________________________________
1937 void AliJetFlowTools::ResetAliUnfolding() {
1938      // ugly function: reset all unfolding parameters 
1939      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
1940      if(fitter) {
1941          printf(" > Found fitter, will delete it < \n");
1942          delete fitter;
1943      }
1944      if(gMinuit) {
1945          printf(" > Found gMinuit, will re-create it < \n");
1946          delete gMinuit;
1947          gMinuit = new TMinuit;
1948      }
1949      AliUnfolding::fgCorrelationMatrix = 0;
1950      AliUnfolding::fgCorrelationMatrixSquared = 0;
1951      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
1952      AliUnfolding::fgCurrentESDVector = 0;
1953      AliUnfolding::fgEntropyAPriori = 0;
1954      AliUnfolding::fgEfficiency = 0;
1955      AliUnfolding::fgUnfoldedAxis = 0;
1956      AliUnfolding::fgMeasuredAxis = 0;
1957      AliUnfolding::fgFitFunction = 0;
1958      AliUnfolding::fgMaxInput  = -1;
1959      AliUnfolding::fgMaxParams = -1;
1960      AliUnfolding::fgOverflowBinLimit = -1;
1961      AliUnfolding::fgRegularizationWeight = 10000;
1962      AliUnfolding::fgSkipBinsBegin = 0;
1963      AliUnfolding::fgMinuitStepSize = 0.1;
1964      AliUnfolding::fgMinuitPrecision = 1e-6;
1965      AliUnfolding::fgMinuitMaxIterations = 1000000;
1966      AliUnfolding::fgMinuitStrategy = 1.;
1967      AliUnfolding::fgMinimumInitialValue = kFALSE;
1968      AliUnfolding::fgMinimumInitialValueFix = -1;
1969      AliUnfolding::fgNormalizeInput = kFALSE;
1970      AliUnfolding::fgNotFoundEvents = 0;
1971      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
1972      AliUnfolding::fgBayesianSmoothing  = 1;
1973      AliUnfolding::fgBayesianIterations = 10;
1974      AliUnfolding::fgDebug = kFALSE;
1975      AliUnfolding::fgCallCount = 0;
1976      AliUnfolding::fgPowern = 5;
1977      AliUnfolding::fChi2FromFit = 0.;
1978      AliUnfolding::fPenaltyVal  = 0.;
1979      AliUnfolding::fAvgResidual = 0.;
1980      AliUnfolding::fgPrintChi2Details = 0;
1981      AliUnfolding::fgCanvas = 0;
1982      AliUnfolding::fghUnfolded = 0;     
1983      AliUnfolding::fghCorrelation = 0;  
1984      AliUnfolding::fghEfficiency = 0;
1985      AliUnfolding::fghMeasured = 0;   
1986      AliUnfolding::SetMinuitStepSize(1.);
1987      AliUnfolding::SetMinuitPrecision(1e-6);
1988      AliUnfolding::SetMinuitMaxIterations(100000);
1989      AliUnfolding::SetMinuitStrategy(2.);
1990      AliUnfolding::SetDebug(1);
1991 }
1992 //_____________________________________________________________________________
1993 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) {
1994     // protect heap by adding unique qualifier to name
1995     if(!protect) return 0x0;
1996     TH1D* p = (TH1D*)protect->Clone();
1997     TString tempString(fActiveString);
1998     tempString+=suffix;
1999     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2000     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2001     if(kill) delete protect;
2002     return p;
2003 }
2004 //_____________________________________________________________________________
2005 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) {
2006     // protect heap by adding unique qualifier to name
2007     if(!protect) return 0x0;
2008     TH2D* p = (TH2D*)protect->Clone();
2009     TString tempString(fActiveString);
2010     tempString+=suffix;
2011     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2012     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2013     if(kill) delete protect;
2014     return p;
2015 }
2016 //_____________________________________________________________________________
2017 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) {
2018     // protect heap by adding unique qualifier to name
2019     if(!protect) return 0x0;
2020     TGraphErrors* p = (TGraphErrors*)protect->Clone();
2021     TString tempString(fActiveString);
2022     tempString+=suffix;
2023     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2024     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2025     if(kill) delete protect;
2026     return p;
2027 }
2028 //_____________________________________________________________________________