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 (see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html).
  29 // A test mode is available in which the spectrum is unfolded with a generated unity response
  30 // matrix.
  31 //
  32 // The weak spot of this class is the function PrepareForUnfolding, which will read
  33 // output from two output files and expects histograms with certain names and binning.
  34 // Unfolding methods itself are general and should be able to handle any input, therefore one
  35 // can forgo the PrepareForUnfolding method, and supply necessary input information via the
  36 // SetRawInput() method
  37 //
  38 // to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
  39
  40 // root includes
  41 #include "TF1.h"
  42 #include "TH1D.h"
  43 #include "TH2D.h"
  44 #include "THStack.h"
  45 #include "TGraph.h"
  46 #include "TGraphErrors.h"
  47 #include "TCanvas.h"
  48 #include "TLegend.h"
  49 #include "TArrayD.h"
  50 #include "TList.h"
  51 #include "TMinuit.h"
  52 #include "TVirtualFitter.h"
  53 #include "TLegend.h"
  54 #include "TCanvas.h"
  55 #include "TStyle.h"
  56 #include "TLine.h"
  57 #include "TMath.h"
  58 #include "TVirtualFitter.h"
  59 #include "TFitResultPtr.h"
  60 // aliroot includes
  61 #include "AliUnfolding.h"
  62 #include "AliAnaChargedJetResponseMaker.h"
  63 // class includes
  64 #include "AliJetFlowTools.h"
  65 // roo unfold includes (make sure you have these available on your system)
  66 #include "RooUnfold.h"
  67 #include "RooUnfoldResponse.h"
  68 #include "RooUnfoldSvd.h"
  69 #include "RooUnfoldBayes.h"
  70 #include "TSVDUnfold.h"
  71
  72 using namespace std;
  73 //_____________________________________________________________________________
  74 AliJetFlowTools::AliJetFlowTools() :
  75     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
  76     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
  77     fSaveFull           (kFALSE),
  78     fActiveString       (""),
  79     fActiveDir          (0x0),
  80     fInputList          (0x0),
  81     fRefreshInput       (kTRUE),
  82     fOutputFileName     ("UnfoldedSpectra.root"),
  83     fOutputFile         (0x0),
  84     fCentralityBin      (0),
  85     fDetectorResponse   (0x0),
  86     fJetFindingEff      (0x0),
  87     fBetaIn             (.1),
  88     fBetaOut            (.1),
  89     fBayesianIterIn     (4),
  90     fBayesianIterOut    (4),
  91     fBayesianSmoothIn   (0.),
  92     fBayesianSmoothOut  (0.),
  93     fAvoidRoundingError (kFALSE),
  94     fUnfoldingAlgorithm (kChi2),
  95     fPrior              (kPriorMeasured),
  96     fBinsTrue           (0x0),
  97     fBinsRec            (0x0),
  98     fBinsTruePrior      (0x0),
  99     fBinsRecPrior       (0x0),
 100     fSVDRegIn           (5),
 101     fSVDRegOut          (5),
 102     fSVDToy             (kTRUE),
 103     fJetRadius          (0.3),
 104     fEventCount         (-1),
 105     fNormalizeSpectra   (kFALSE),
 106     fSmoothenPrior      (kFALSE),
 107     fFitMin             (60.),
 108     fFitMax             (300.),
 109     fFitStart           (75.),
 110     fSmoothenCounts     (kTRUE),
 111     fTestMode           (kFALSE),
 112     fNoDphi             (kFALSE),
 113     fRawInputProvided   (kFALSE),
 114     fEventPlaneRes      (.63),
 115     fUseDetectorResponse(kTRUE),
 116     fUseDptResponse     (kTRUE),
 117     fTrainPower         (kTRUE),
 118     fRMSSpectrumIn      (0x0),
 119     fRMSSpectrumOut     (0x0),
 120     fRMSRatio           (0x0),
 121     fRMSV2              (0x0),
 122     fDeltaPtDeltaPhi    (0x0),
 123     fJetPtDeltaPhi      (0x0),
 124     fSpectrumIn         (0x0),
 125     fSpectrumOut        (0x0),
 126     fDptInDist          (0x0),
 127     fDptOutDist         (0x0),
 128     fDptIn              (0x0),
 129     fDptOut             (0x0),
 130     fFullResponseIn     (0x0),
 131     fFullResponseOut    (0x0) { // class constructor
 132     // create response maker weight function (tuned to PYTHIA spectrum)
 133     fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
 134 }
 135 //_____________________________________________________________________________
 136 void AliJetFlowTools::Make() {
 137     // core function of the class
 138     // 1) rebin the raw output of the jet task to the desired binnings
 139     // 2) calls the unfolding routine
 140     // 3) writes output to file
 141     // can be repeated multiple times with different configurations
 142
 143     // 1) manipulation of input histograms
 144     // check if the input variables are present
 145     if(fRefreshInput) {
 146         if(!PrepareForUnfolding()) {
 147             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
 148             return;
 149         }
 150     }
 151     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
 152     //     parts of the spectrum can end up in over or underflow bins
 153     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
 154     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
 155
 156     // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
 157     // the template will be used as a prior for the chi2 unfolding
 158     TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
 159     TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
 160
 161     // get the full response matrix from the dpt and the detector response
 162     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
 163     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
 164     // so that unfolding should return the initial spectrum
 165     if(!fTestMode) {
 166         if(fUseDptResponse && fUseDetectorResponse) {
 167             fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
 168             fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
 169         } else if (fUseDptResponse && !fUseDetectorResponse) {
 170             fFullResponseIn = fDptIn;
 171             fFullResponseOut = fDptOut;
 172         } else if (!fUseDptResponse && fUseDetectorResponse) {
 173             fFullResponseIn = fDetectorResponse;
 174             fFullResponseOut = fDetectorResponse;
 175         } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
 176             printf(" > No response, exiting ! < \n" );
 177             return;
 178         }
 179     } else {
 180         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
 181         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
 182     }
 183     // normalize each slide of the response to one
 184     NormalizeTH2D(fFullResponseIn);
 185     NormalizeTH2D(fFullResponseOut);
 186     // resize to desired binning scheme
 187     TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
 188     TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
 189     // get the kinematic efficiency
 190     TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
 191     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
 192     TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
 193     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
 194     // suppress the errors
 195     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
 196         kinematicEfficiencyIn->SetBinError(1+i, 0.);
 197         kinematicEfficiencyOut->SetBinError(1+i, 0.);
 198     }
 199     TH1D* jetFindingEfficiency(0x0);
 200     if(fJetFindingEff) {
 201         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
 202         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
 203         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
 204     }
 205     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
 206     TH1D* unfoldedJetSpectrumIn(0x0);
 207     TH1D* unfoldedJetSpectrumOut(0x0);
 208     fActiveDir->cd();                   // select active dir
 209     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
 210     dirIn->cd();                        // select inplane subdir
 211     // select the unfolding method
 212     switch (fUnfoldingAlgorithm) {
 213         case kChi2 : {
 214             unfoldedJetSpectrumIn = UnfoldSpectrumChi2(       // do the inplane unfolding
 215                 measuredJetSpectrumIn,
 216                 resizedResponseIn,
 217                 kinematicEfficiencyIn,
 218                 measuredJetSpectrumTrueBinsIn,
 219                 TString("in"),
 220                 jetFindingEfficiency);
 221             printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
 222         } break;
 223         case kBayesian : {
 224             unfoldedJetSpectrumIn = UnfoldSpectrumBayesian(       // do the inplane unfolding
 225                 measuredJetSpectrumIn,
 226                 resizedResponseIn,
 227                 kinematicEfficiencyIn,
 228                 measuredJetSpectrumTrueBinsIn,
 229                 TString("in"),
 230                 jetFindingEfficiency);
 231             printf(" > Spectrum (in plane) unfolded using kBayesian unfolding < \n");
 232         } break;
 233         case kBayesianAli : {
 234             unfoldedJetSpectrumIn = UnfoldSpectrumBayesianAli(       // do the inplane unfolding
 235                 measuredJetSpectrumIn,
 236                 resizedResponseIn,
 237                 kinematicEfficiencyIn,
 238                 measuredJetSpectrumTrueBinsIn,
 239                 TString("in"),
 240                 jetFindingEfficiency);
 241             printf(" > Spectrum (in plane) unfolded using kBayesianAli unfolding < \n");
 242         } break;
 243         case kSVD : {
 244             unfoldedJetSpectrumIn = UnfoldSpectrumSVD(       // do the inplane unfolding
 245                 measuredJetSpectrumIn,
 246                 resizedResponseIn,
 247                 kinematicEfficiencyIn,
 248                 measuredJetSpectrumTrueBinsIn,
 249                 TString("in"),
 250                 jetFindingEfficiency);
 251             printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
 252         } break;
 253         case kNone : {  // do nothing
 254             resizedResponseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
 255             unfoldedJetSpectrumIn = ProtectHeap(measuredJetSpectrumIn, kTRUE, TString("in"));
 256         } break;
 257
 258         default : {
 259             printf(" > Selected unfolding method is not implemented yet ! \n");
 260             return;
 261         }
 262     }
 263     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
 264     resizedResponseIn->SetXTitle("p_{T}^{true} [GeV/c]");
 265     resizedResponseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
 266     resizedResponseIn = ProtectHeap(resizedResponseIn);
 267     resizedResponseIn->Write();
 268     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
 269     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
 270     kinematicEfficiencyIn->Write();
 271     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 272     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 273     fDetectorResponse->Write();
 274     // optional histograms
 275     if(fSaveFull) {
 276         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
 277         fSpectrumIn->Write();
 278         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
 279         fDptInDist->Write();
 280         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
 281         fDptIn->Write();
 282         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
 283         fFullResponseIn->Write();
 284     }
 285     fActiveDir->cd();
 286     TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
 287     dirOut->cd();
 288     switch (fUnfoldingAlgorithm) {
 289         case kChi2 : {
 290             unfoldedJetSpectrumOut = UnfoldSpectrumChi2(
 291                 measuredJetSpectrumOut,
 292                 resizedResponseOut,
 293                 kinematicEfficiencyOut,
 294                 measuredJetSpectrumTrueBinsOut,
 295                 TString("out"),
 296                 jetFindingEfficiency);
 297             printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
 298         } break;
 299         case kBayesian : {
 300             unfoldedJetSpectrumOut = UnfoldSpectrumBayesian(
 301                 measuredJetSpectrumOut,
 302                 resizedResponseOut,
 303                 kinematicEfficiencyOut,
 304                 measuredJetSpectrumTrueBinsOut,
 305                 TString("out"),
 306                 jetFindingEfficiency);
 307             printf(" > Spectrum (out of plane) unfolded using kBayesian < \n");
 308         } break;
 309         case kBayesianAli : {
 310             unfoldedJetSpectrumOut = UnfoldSpectrumBayesianAli(
 311                 measuredJetSpectrumOut,
 312                 resizedResponseOut,
 313                 kinematicEfficiencyOut,
 314                 measuredJetSpectrumTrueBinsOut,
 315                 TString("out"),
 316                 jetFindingEfficiency);
 317             printf(" > Spectrum (out of plane) unfolded using kBayesianAli < \n");
 318         } break;
 319         case kSVD : {
 320             unfoldedJetSpectrumOut = UnfoldSpectrumSVD(
 321                 measuredJetSpectrumOut,
 322                 resizedResponseOut,
 323                 kinematicEfficiencyOut,
 324                 measuredJetSpectrumTrueBinsOut,
 325                 TString("out"),
 326                 jetFindingEfficiency);
 327             printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
 328         } break;
 329         case kNone : {  // do nothing
 330             resizedResponseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
 331             unfoldedJetSpectrumOut = ProtectHeap(measuredJetSpectrumOut, kTRUE, TString("out"));
 332         } break;
 333         default : {
 334             printf(" > Selected unfolding method is not implemented yet ! \n");
 335             return;
 336         }
 337     }
 338     resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
 339     resizedResponseOut->SetXTitle("p_{T}^{true} [GeV/c]");
 340     resizedResponseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
 341     resizedResponseOut = ProtectHeap(resizedResponseOut);
 342     resizedResponseOut->Write();
 343     kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
 344     kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
 345     kinematicEfficiencyOut->Write();
 346     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 347     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 348     fDetectorResponse->Write();
 349     if(jetFindingEfficiency) jetFindingEfficiency->Write();
 350     // optional histograms
 351     if(fSaveFull) {
 352         fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
 353         fSpectrumOut->Write();
 354         fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
 355         fDptOutDist->Write();
 356         fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
 357         fDptOut->Write();
 358         fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
 359         fFullResponseOut->Write();
 360     }
 361
 362     // write general output histograms to file
 363     fActiveDir->cd();
 364     if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
 365         TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
 366         if(ratio) {
 367             ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
 368             ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 369             ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 370             ratio = ProtectHeap(ratio);
 371             ratio->Write();
 372             // write histo values to RMS files if both routines converged
 373             // input values are weighted by their uncertainty
 374             for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
 375                 if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
 376                 if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
 377                 if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
 378            }
 379         }
 380         TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
 381         if(v2) {
 382             v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 383             v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 384             v2->GetYaxis()->SetTitle("v_{2}");
 385             v2 = ProtectHeap(v2);
 386             v2->Write();
 387         }
 388     } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
 389         TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
 390         if(ratio) {
 391             ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
 392             ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 393             ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 394             ratio = ProtectHeap(ratio);
 395             ratio->Write();
 396         }
 397         TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
 398          if(v2) {
 399             v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 400             v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 401             v2->GetYaxis()->SetTitle("v_{2}");
 402             v2 = ProtectHeap(v2);
 403             v2->Write();
 404         }
 405     }
 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) { // 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(fNoDphi) {
 914         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40);
 915         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40);
 916     } else {
 917         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10);
 918         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40));
 919         fSpectrumIn = ProtectHeap(fSpectrumIn);
 920         // out of plane spectrum
 921         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30);
 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(fNoDphi) {
 966         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40);
 967         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40);
 968     } else {
 969         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10);
 970         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40));
 971         // out of plane delta pt distribution
 972         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30);
 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(new TCanvas("PearsonOut", "PearsonOut"));
1407    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
1408    TCanvas* canvasRatioMeasuredRefoldedOut(new TCanvas("RefoldedOut", "RefoldedOut"));
1409    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
1410    TCanvas* canvasSpectraOut(new TCanvas("SpectraOut", "SpectraOut"));
1411    TCanvas* canvasRatio(new TCanvas("Ratio", "Ratio"));
1412    TCanvas* canvasV2(new TCanvas("V2", "V2"));
1413    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1414    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
1415    TCanvas* canvasMasterOut(new TCanvas("defaultOut", "defaultOut"));
1416    canvasMISC->Divide(4, 2);
1417    TDirectoryFile* defDir(0x0);
1418
1419    // get an estimate of the number of outputs and find the default set
1420    Int_t cacheMe(0);
1421    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1422        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1423            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1424            cacheMe++;
1425        }
1426    }
1427    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%4)>0));
1428    canvasIn->Divide(columns, rows);
1429    canvasOut->Divide(columns, rows);
1430    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1431    canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1432    canvasSpectraIn->Divide(columns, rows);
1433    canvasSpectraOut->Divide(columns, rows);
1434    canvasRatio->Divide(columns, rows);
1435    canvasV2->Divide(columns, rows);
1436
1437    canvasMasterIn->Divide(columns, rows);
1438    canvasMasterOut->Divide(columns, rows);
1439    // extract the default output
1440    TH1D* deunfoldedJetSpectrumIn(0x0);
1441    TH1D* deunfoldedJetSpectrumOut(0x0);
1442    THStack stackIn("StackRatioIn","StackRatioIn");
1443    THStack stackOut("StackRatioOut", "StackRatioOut");
1444    if(defDir) {
1445        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1446        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
1447        if(defDirIn) deunfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1448        if(deunfoldedJetSpectrumIn) stackIn.Add(deunfoldedJetSpectrumIn);
1449        if(defDirOut) deunfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
1450        if(deunfoldedJetSpectrumOut) stackOut.Add(deunfoldedJetSpectrumOut);
1451        printf(" > succesfully extracted default results < \n");
1452    }
1453
1454    // loop through the directories, only plot the graphs if the deconvolution converged
1455    TDirectoryFile* tempDir(0x0);
1456    TDirectoryFile* tempIn(0x0);
1457    TDirectoryFile*  tempOut(0x0);
1458    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
1459        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1460        if(!tempDir) continue;
1461        TString dirName(tempDir->GetName());
1462        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1463        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1464        j++;
1465        if(tempIn) {
1466            // to see if the unfolding converged try to extract the pearson coefficients
1467            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1468            if(pIn) {
1469                printf(" - %s in plane converged \n", dirName.Data());
1470                canvasIn->cd(j);
1471                Style(gPad, "PEARSON");
1472                pIn->DrawCopy("colz");
1473                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1474                if(rIn) {
1475                    Style(rIn);
1476                    printf(" > found RatioRefoldedMeasured < \n");
1477                    canvasRatioMeasuredRefoldedIn->cd(j);
1478                    rIn->Draw("ac");
1479                }
1480                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1481                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1482                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1483                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1484                if(dvector && avalue && rm && eff) {
1485                    Style(dvector);
1486                    Style(avalue);
1487                    Style(rm);
1488                    Style(eff);
1489                    canvasMISC->cd(1);
1490                    Style(gPad, "SPECTRUM");
1491                    dvector->DrawCopy();
1492                    canvasMISC->cd(2);
1493                    Style(gPad, "SPECTRUM");
1494                    avalue->DrawCopy();
1495                    canvasMISC->cd(3);
1496                    Style(gPad, "PEARSON");
1497                    rm->DrawCopy("colz");
1498                    canvasMISC->cd(4);
1499                    eff->DrawCopy();
1500                } else if(rm && eff) {
1501                    Style(rm);
1502                    Style(eff);
1503                    canvasMISC->cd(3);
1504                    Style(gPad, "PEARSON");
1505                    rm->DrawCopy("colz");
1506                    canvasMISC->cd(4);
1507                    eff->DrawCopy();
1508                }
1509            }
1510            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1511            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1512            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1513            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1514                if(deunfoldedJetSpectrumIn) {
1515                    TH1D* temp((TH1D*)deunfoldedJetSpectrumIn->Clone(Form("deunfoldedJetSpectrumIn_%s", dirName.Data())));
1516                    temp->Divide(unfoldedSpectrum);
1517                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1518                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1519                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1520                    canvasMasterIn->cd(j);
1521                    temp->GetXaxis()->SetRangeUser(0., 2);
1522                    temp->DrawCopy();
1523                }
1524                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
1525                canvasSpectraIn->cd(j);
1526                Style(gPad);
1527                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1528                unfoldedSpectrum->DrawCopy();
1529                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1530                inputSpectrum->DrawCopy("same");
1531                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1532                refoldedSpectrum->DrawCopy("same");
1533                TLegend* l(AddLegend(gPad));
1534                Style(l);
1535                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1536                    Float_t chi(fitStatus->GetBinContent(1));
1537                    Float_t pen(fitStatus->GetBinContent(2));
1538                    Int_t dof((int)fitStatus->GetBinContent(3));
1539                    Float_t beta(fitStatus->GetBinContent(4));
1540                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1541                } else if (fitStatus) { // only available in SVD fit
1542                    Int_t reg((int)fitStatus->GetBinContent(1));
1543                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1544                }
1545            }
1546        }
1547        if(tempOut) {
1548            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
1549            if(pOut) {
1550                printf(" - %s out of plane converged \n", dirName.Data());
1551                canvasOut->cd(j);
1552                Style(gPad, "PEARSON");
1553                pOut->DrawCopy("colz");
1554                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1555                if(rOut) {
1556                    Style(rOut);
1557                    printf(" > found RatioRefoldedMeasured < \n");
1558                    canvasRatioMeasuredRefoldedOut->cd(j);
1559                    rOut->Draw("ac");
1560                }
1561                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1562                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1563                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
1564                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
1565                if(dvector && avalue && rm && eff) {
1566                    Style(dvector);
1567                    Style(avalue);
1568                    Style(rm);
1569                    Style(eff);
1570                    canvasMISC->cd(5);
1571                    Style(gPad, "SPECTRUM");
1572                    dvector->DrawCopy();
1573                    canvasMISC->cd(6);
1574                    Style(gPad, "SPECTRUM");
1575                    avalue->DrawCopy();
1576                    canvasMISC->cd(7);
1577                    Style(gPad, "PEARSON");
1578                    rm->DrawCopy("colz");
1579                    canvasMISC->cd(8);
1580                    eff->DrawCopy();
1581                } else if(rm && eff) {
1582                    Style(rm);
1583                    Style(eff);
1584                    canvasMISC->cd(7);
1585                    Style(gPad, "PEARSON");
1586                    rm->DrawCopy("colz");
1587                    canvasMISC->cd(8);
1588                    eff->DrawCopy();
1589                }
1590            }
1591            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
1592            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
1593            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
1594            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1595                if(deunfoldedJetSpectrumOut) {
1596                    TH1D* temp((TH1D*)deunfoldedJetSpectrumOut->Clone(Form("deunfoldedJetSpectrumOut_%s", dirName.Data())));
1597                    temp->Divide(unfoldedSpectrum);
1598                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1599                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1600                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1601                    canvasMasterOut->cd(j);
1602                    temp->GetXaxis()->SetRangeUser(0., 2.);
1603                    temp->DrawCopy();
1604                }
1605                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
1606                canvasSpectraOut->cd(j);
1607                Style(gPad);
1608                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1609                unfoldedSpectrum->DrawCopy();
1610                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1611                inputSpectrum->DrawCopy("same");
1612                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1613                refoldedSpectrum->DrawCopy("same");
1614                TLegend* l(AddLegend(gPad));
1615                Style(l);
1616                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1617                    Float_t chi(fitStatus->GetBinContent(1));
1618                    Float_t pen(fitStatus->GetBinContent(2));
1619                    Int_t dof((int)fitStatus->GetBinContent(3));
1620                    Float_t beta(fitStatus->GetBinContent(4));
1621                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1622                } else if (fitStatus) { // only available in SVD fit
1623                    Int_t reg((int)fitStatus->GetBinContent(1));
1624                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1625                }
1626            }
1627        }
1628        canvasRatio->cd(j);
1629        TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
1630        if(ratioYield) {
1631            Style(ratioYield);
1632            ratioYield->GetYaxis()->SetRangeUser(-1., 3.);
1633            ratioYield->Draw("ac");
1634        }
1635        canvasV2->cd(j);
1636        TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
1637        if(ratioV2) {
1638            Style(ratioV2);
1639            ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
1640            ratioV2->Draw("ac");
1641        }
1642    }
1643    TFile output(out.Data(), "RECREATE");
1644    SavePadToPDF(canvasIn);
1645    canvasIn->Write();
1646    SavePadToPDF(canvasOut);
1647    canvasOut->Write();
1648    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
1649    canvasRatioMeasuredRefoldedIn->Write();
1650    SavePadToPDF(canvasRatioMeasuredRefoldedOut);
1651    canvasRatioMeasuredRefoldedOut->Write();
1652    SavePadToPDF(canvasSpectraIn);
1653    canvasSpectraIn->Write();
1654    SavePadToPDF(canvasSpectraOut);
1655    canvasSpectraOut->Write();
1656    SavePadToPDF(canvasRatio);
1657    canvasRatio->Write();
1658    SavePadToPDF(canvasV2);
1659    canvasV2->Write();
1660    SavePadToPDF(canvasMasterIn);
1661    canvasMasterIn->Write();
1662    SavePadToPDF(canvasMasterIn);
1663    canvasMasterOut->Write();
1664    SavePadToPDF(canvasMISC);
1665    canvasMISC->Write();
1666    output.Write();
1667    output.Close();
1668 }
1669 //_____________________________________________________________________________
1670 Bool_t AliJetFlowTools::SetRawInput (
1671         TH2D* detectorResponse,  // detector response matrix
1672         TH1D* jetPtIn,           // in plane jet spectrum
1673         TH1D* jetPtOut,          // out of plane jet spectrum
1674         TH1D* dptIn,             // in plane delta pt distribution
1675         TH1D* dptOut,            // out of plane delta pt distribution
1676         Int_t eventCount) {
1677     // set input histograms manually
1678     fDetectorResponse   = detectorResponse;
1679     fSpectrumIn         = jetPtIn;
1680     fSpectrumOut        = jetPtOut;
1681     fDptInDist          = dptIn;
1682     fDptOutDist         = dptOut;
1683     fRawInputProvided   = kTRUE;
1684     // check if all data is provided
1685     if(!fDetectorResponse) {
1686         printf(" fDetectorResponse not found \n ");
1687         return kFALSE;
1688     }
1689     // check if the pt bin for true and rec have been set
1690     if(!fBinsTrue || !fBinsRec) {
1691         printf(" No true or rec bins set, please set binning ! \n");
1692         return kFALSE;
1693     }
1694     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
1695                           // procedures, these profiles will be nonsensical, user is responsible
1696         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1697         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1698         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1699     }
1700     // normalize spectra to event count if requested
1701     if(fNormalizeSpectra) {
1702         fEventCount = eventCount;
1703         if(fEventCount > 0) {
1704             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1705             fSpectrumOut->Sumw2();
1706             fSpectrumIn->Scale(1./((double)fEventCount));
1707             fSpectrumOut->Scale(1./((double)fEventCount));
1708         }
1709     }
1710     if(!fNormalizeSpectra && fEventCount > 0) {
1711         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1712         fSpectrumOut->Sumw2();
1713         fSpectrumIn->Scale(1./((double)fEventCount));
1714         fSpectrumOut->Scale(1./((double)fEventCount));
1715     }
1716     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
1717     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1718     fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1719     fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1720     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
1721     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1722     fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1723     fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1724
1725     return kTRUE;
1726 }
1727 //_____________________________________________________________________________
1728 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
1729 {
1730     // return ratio of h1 / h2
1731     // histograms can have different binning. errors are propagated as uncorrelated
1732     if(!(h1 && h2) ) {
1733         printf(" GetRatio called with NULL argument(s) \n ");
1734         return 0x0;
1735     }
1736     Int_t j(0);
1737     TGraphErrors *gr = new TGraphErrors();
1738     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1739     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1740     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1741         binCent = h1->GetXaxis()->GetBinCenter(i);
1742         if(xmax > 0. && binCent > xmax) continue;
1743         j = h2->FindBin(binCent);
1744         binWidth = h1->GetXaxis()->GetBinWidth(i);
1745         if(h2->GetBinContent(j) > 0.) {
1746             ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
1747             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
1748             Double_t B = 0.;
1749             if(h2->GetBinError(j)>0.) {
1750                 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
1751                 error2 = A*A + B*B;
1752             } else error2 = A*A;
1753             if(error2 > 0 ) error2 = TMath::Sqrt(error2);
1754             gr->SetPoint(gr->GetN(),binCent,ratio);
1755             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1756         }
1757     }
1758     if(appendFit) {
1759         TF1* fit(new TF1("lin", "pol0", 10, 100));
1760         gr->Fit(fit);
1761     }
1762     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1763     return gr;
1764 }
1765 //_____________________________________________________________________________
1766 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
1767 {
1768     // get v2 from difference of in plane, out of plane yield
1769     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
1770     // different binning is allowed but will mean that the error
1771     // propagation is unreliable
1772     // r is the event plane resolution for the chosen centrality
1773     if(!(h1 && h2) ) {
1774         printf(" GetV2 called with NULL argument(s) \n ");
1775         return 0x0;
1776     }
1777     Int_t j(0);
1778     TGraphErrors *gr = new TGraphErrors();
1779     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1780     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1781     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
1782     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1783         binCent = h1->GetXaxis()->GetBinCenter(i);
1784         j = h2->FindBin(binCent);
1785         binWidth = h1->GetXaxis()->GetBinWidth(i);
1786         if(h2->GetBinContent(j) > 0.) {
1787             in = h1->GetBinContent(i);
1788             ein = h1->GetBinError(i);
1789             out = h2->GetBinContent(j);
1790             eout = h2->GetBinError(j);
1791             ratio = pre*((in-out)/(in+out));
1792             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);
1793             if(error2 > 0) error2 = TMath::Sqrt(error2);
1794             gr->SetPoint(gr->GetN(),binCent,ratio);
1795             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1796         }
1797     }
1798     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1799     return gr;
1800 }
1801 //_____________________________________________________________________________
1802 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
1803     // write object, if a unique identifier is given the object is cloned
1804     // and the clone is saved. setting kill to true will delete the original obect from the heap
1805     if(!object) {
1806         printf(" > WriteObject:: called with NULL arguments \n ");
1807         return;
1808     } else if(!strcmp("", suffix.Data())) object->Write();
1809     else {
1810         TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
1811         newObject->Write();
1812     }
1813     if(kill) delete object;
1814 }
1815 //_____________________________________________________________________________
1816 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
1817     // construt a delta pt response matrix from supplied dpt distribution
1818     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
1819     // do a weighted rebinning to a (coarser) dpt distribution
1820     // be careful with the binning of the dpt response: it should be equal to that
1821     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
1822     //
1823     // the response matrix will be square and have the same binning
1824     // (min, max and granularity) of the input histogram
1825     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
1826     Double_t _bins[bins+1];             // prepare array with bin borders
1827     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
1828     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
1829     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
1830     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
1831         Bool_t skip = kFALSE;
1832         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
1833             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
1834             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
1835         }
1836     }
1837     return res;
1838 }
1839 //_____________________________________________________________________________
1840 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1841     if(!binsTrue || !binsRec) {
1842         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
1843         return 0x0;
1844     }
1845     TString name(Form("unityResponse_%s", suffix.Data()));
1846     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
1847     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
1848         for(Int_t j(0); j < binsRec->GetSize(); j++) {
1849             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
1850         }
1851     }
1852     return unity;
1853 }
1854 //_____________________________________________________________________________
1855 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
1856     // save configuration parameters to histogram
1857     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
1858     summary->SetBinContent(1, fBetaIn);
1859     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
1860     summary->SetBinContent(2, fBetaOut);
1861     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
1862     summary->SetBinContent(3, fCentralityBin);
1863     summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
1864     summary->SetBinContent(4, (int)convergedIn);
1865     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
1866     summary->SetBinContent(5, (int)convergedOut);
1867     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
1868     summary->SetBinContent(6, (int)fAvoidRoundingError);
1869     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
1870     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
1871     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
1872     summary->SetBinContent(8, (int)fPrior);
1873     summary->GetXaxis()->SetBinLabel(8, "fPrior");
1874     summary->SetBinContent(9, fSVDRegIn);
1875     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
1876     summary->SetBinContent(10, fSVDRegOut);
1877     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
1878     summary->SetBinContent(11, (int)fSVDToy);
1879     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
1880     summary->SetBinContent(12, fJetRadius);
1881     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
1882     summary->SetBinContent(13, (int)fNormalizeSpectra);
1883     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
1884     summary->SetBinContent(14, (int)fSmoothenPrior);
1885     summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
1886     summary->SetBinContent(15, (int)fTestMode);
1887     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
1888     summary->SetBinContent(16, (int)fUseDetectorResponse);
1889     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
1890     summary->SetBinContent(17, fBayesianIterIn);
1891     summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
1892     summary->SetBinContent(18, fBayesianIterOut);
1893     summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
1894     summary->SetBinContent(19, fBayesianSmoothIn);
1895     summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
1896     summary->SetBinContent(20, fBayesianSmoothOut);
1897     summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
1898 }
1899 //_____________________________________________________________________________
1900 void AliJetFlowTools::ResetAliUnfolding() {
1901      // ugly function: reset all unfolding parameters
1902      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
1903      if(fitter) {
1904          printf(" > Found fitter, will delete it < \n");
1905          delete fitter;
1906      }
1907      if(gMinuit) {
1908          printf(" > Found gMinuit, will re-create it < \n");
1909          delete gMinuit;
1910          gMinuit = new TMinuit;
1911      }
1912      AliUnfolding::fgCorrelationMatrix = 0;
1913      AliUnfolding::fgCorrelationMatrixSquared = 0;
1914      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
1915      AliUnfolding::fgCurrentESDVector = 0;
1916      AliUnfolding::fgEntropyAPriori = 0;
1917      AliUnfolding::fgEfficiency = 0;
1918      AliUnfolding::fgUnfoldedAxis = 0;
1919      AliUnfolding::fgMeasuredAxis = 0;
1920      AliUnfolding::fgFitFunction = 0;
1921      AliUnfolding::fgMaxInput  = -1;
1922      AliUnfolding::fgMaxParams = -1;
1923      AliUnfolding::fgOverflowBinLimit = -1;
1924      AliUnfolding::fgRegularizationWeight = 10000;
1925      AliUnfolding::fgSkipBinsBegin = 0;
1926      AliUnfolding::fgMinuitStepSize = 0.1;
1927      AliUnfolding::fgMinuitPrecision = 1e-6;
1928      AliUnfolding::fgMinuitMaxIterations = 1000000;
1929      AliUnfolding::fgMinuitStrategy = 1.;
1930      AliUnfolding::fgMinimumInitialValue = kFALSE;
1931      AliUnfolding::fgMinimumInitialValueFix = -1;
1932      AliUnfolding::fgNormalizeInput = kFALSE;
1933      AliUnfolding::fgNotFoundEvents = 0;
1934      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
1935      AliUnfolding::fgBayesianSmoothing  = 1;
1936      AliUnfolding::fgBayesianIterations = 10;
1937      AliUnfolding::fgDebug = kFALSE;
1938      AliUnfolding::fgCallCount = 0;
1939      AliUnfolding::fgPowern = 5;
1940      AliUnfolding::fChi2FromFit = 0.;
1941      AliUnfolding::fPenaltyVal  = 0.;
1942      AliUnfolding::fAvgResidual = 0.;
1943      AliUnfolding::fgPrintChi2Details = 0;
1944      AliUnfolding::fgCanvas = 0;
1945      AliUnfolding::fghUnfolded = 0;
1946      AliUnfolding::fghCorrelation = 0;
1947      AliUnfolding::fghEfficiency = 0;
1948      AliUnfolding::fghMeasured = 0;
1949      AliUnfolding::SetMinuitStepSize(1.);
1950      AliUnfolding::SetMinuitPrecision(1e-6);
1951      AliUnfolding::SetMinuitMaxIterations(100000);
1952      AliUnfolding::SetMinuitStrategy(2.);
1953      AliUnfolding::SetDebug(1);
1954 }
1955 //_____________________________________________________________________________
1956 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) {
1957     // protect heap by adding unique qualifier to name
1958     if(!protect) return 0x0;
1959     TH1D* p = (TH1D*)protect->Clone();
1960     TString tempString(fActiveString);
1961     tempString+=suffix;
1962     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1963     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1964     if(kill) delete protect;
1965     return p;
1966 }
1967 //_____________________________________________________________________________
1968 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) {
1969     // protect heap by adding unique qualifier to name
1970     if(!protect) return 0x0;
1971     TH2D* p = (TH2D*)protect->Clone();
1972     TString tempString(fActiveString);
1973     tempString+=suffix;
1974     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1975     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1976     if(kill) delete protect;
1977     return p;
1978 }
1979 //_____________________________________________________________________________
1980 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) {
1981     // protect heap by adding unique qualifier to name
1982     if(!protect) return 0x0;
1983     TGraphErrors* p = (TGraphErrors*)protect->Clone();
1984     TString tempString(fActiveString);
1985     tempString+=suffix;
1986     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1987     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1988     if(kill) delete protect;
1989     return p;
1990 }
1991 //_____________________________________________________________________________