PWG/FLOW/Tasks/AliJetFlowTools.cxx

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