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           (kTRUE),
  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     fDphiUnfolding      (kTRUE),
 116     fDphiDptUnfolding   (kFALSE),
 117     fExLJDpt            (kTRUE),
 118     fRMSSpectrumIn      (0x0),
 119     fRMSSpectrumOut     (0x0),
 120     fRMSRatio           (0x0),
 121     fRMSV2              (0x0),
 122     fDeltaPtDeltaPhi    (0x0),
 123     fJetPtDeltaPhi      (0x0),
 124     fSpectrumIn         (0x0),
 125     fSpectrumOut        (0x0),
 126     fDptInDist          (0x0),
 127     fDptOutDist         (0x0),
 128     fDptIn              (0x0),
 129     fDptOut             (0x0),
 130     fFullResponseIn     (0x0),
 131     fFullResponseOut    (0x0) { // class constructor
 132         // create response maker weight function (tuned to PYTHIA spectrum)
 133         fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
 134         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 135 }
 136 //_____________________________________________________________________________
 137 void AliJetFlowTools::Make() {
 138     // core function of the class
 139     if(fDphiDptUnfolding) {
 140         // to extract the yield as function of Dphi, Dpt - experimental
 141         MakeAU();
 142         return;
 143     }
 144     // 1) rebin the raw output of the jet task to the desired binnings
 145     // 2) calls the unfolding routine
 146     // 3) writes output to file
 147     // can be repeated multiple times with different configurations
 148
 149     // 1) manipulation of input histograms
 150     // check if the input variables are present
 151     if(fRefreshInput) {
 152         if(!PrepareForUnfolding()) {
 153             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
 154             return;
 155         }
 156     }
 157     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
 158     //     parts of the spectrum can end up in over or underflow bins
 159 /*
 160     Double_t binsTrue[] = {-50, -45, -40, -35, -30,-25, -20,-15, -10,-5,  0, 5, 10, 15, 20,25, 30,35, 40,45, 50, 55, 60,65, 70, 75,80,85, 90,95, 100};
 161         TArrayD* temparray = new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue);
 162     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, temparray, TString("resized_in_"), kFALSE);
 163     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, temparray,  TString("resized_out_"), kFALSE);
 164     */
 165     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
 166     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec,  TString("resized_out_"), kFALSE);
 167
 168     // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
 169     // the template will be used as a prior for the chi2 unfolding
 170     TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
 171     TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
 172     // get the full response matrix from the dpt and the detector response
 173     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
 174     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
 175     // so that unfolding should return the initial spectrum
 176     if(!fTestMode) {
 177         if(fUseDptResponse && fUseDetectorResponse) {
 178             fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
 179             fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
 180         } else if (fUseDptResponse && !fUseDetectorResponse) {
 181             fFullResponseIn = fDptIn;
 182             fFullResponseOut = fDptOut;
 183         } else if (!fUseDptResponse && fUseDetectorResponse) {
 184             fFullResponseIn = fDetectorResponse;
 185             fFullResponseOut = fDetectorResponse;
 186         } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
 187             printf(" > No response, exiting ! < \n" );
 188             return;
 189         }
 190     } else {
 191         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
 192         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
 193     }
 194     // normalize each slide of the response to one
 195     NormalizeTH2D(fFullResponseIn);
 196     NormalizeTH2D(fFullResponseOut);
 197     // resize to desired binning scheme
 198     TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
 199     TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
 200     // get the kinematic efficiency
 201     TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
 202     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
 203     TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
 204     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
 205     // suppress the errors
 206     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
 207         kinematicEfficiencyIn->SetBinError(1+i, 0.);
 208         kinematicEfficiencyOut->SetBinError(1+i, 0.);
 209     }
 210     TH1D* jetFindingEfficiency(0x0);
 211     if(fJetFindingEff) {
 212         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
 213         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
 214         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
 215     }
 216     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
 217     TH1D* unfoldedJetSpectrumIn(0x0);
 218     TH1D* unfoldedJetSpectrumOut(0x0);
 219     fActiveDir->cd();                   // select active dir
 220     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
 221     dirIn->cd();                        // select inplane subdir
 222     // do the inplane unfolding
 223     unfoldedJetSpectrumIn = UnfoldWrapper(
 224         measuredJetSpectrumIn,
 225         resizedResponseIn,
 226         kinematicEfficiencyIn,
 227         measuredJetSpectrumTrueBinsIn,
 228         TString("in"),
 229         jetFindingEfficiency);
 230     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
 231     resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 232     resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 233     resizedResponseIn = ProtectHeap(resizedResponseIn);
 234     resizedResponseIn->Write();
 235     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
 236     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
 237     kinematicEfficiencyIn->Write();
 238     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 239     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 240     fDetectorResponse->Write();
 241     // optional histograms
 242     if(fSaveFull) {
 243         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
 244         fSpectrumIn->Write();
 245         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
 246         fDptInDist->Write();
 247         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
 248         fDptIn->Write();
 249         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
 250         fFullResponseIn->Write();
 251     }
 252     fActiveDir->cd();
 253     if(fDphiUnfolding) {
 254         TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
 255         dirOut->cd();
 256         // do the out of plane unfolding
 257         unfoldedJetSpectrumOut = UnfoldWrapper(
 258                     measuredJetSpectrumOut,
 259                     resizedResponseOut,
 260                     kinematicEfficiencyOut,
 261                     measuredJetSpectrumTrueBinsOut,
 262                     TString("out"),
 263                     jetFindingEfficiency);
 264         resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
 265         resizedResponseOut->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 266         resizedResponseOut->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 267         resizedResponseOut = ProtectHeap(resizedResponseOut);
 268         resizedResponseOut->Write();
 269         kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
 270         kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
 271         kinematicEfficiencyOut->Write();
 272         fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 273         fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 274         fDetectorResponse->Write();
 275         if(jetFindingEfficiency) jetFindingEfficiency->Write();
 276         // optional histograms
 277         if(fSaveFull) {
 278             fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
 279             fSpectrumOut->Write();
 280             fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
 281             fDptOutDist->Write();
 282             fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
 283             fDptOut->Write();
 284             fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
 285             fFullResponseOut->Write();
 286         }
 287
 288         // write general output histograms to file
 289         fActiveDir->cd();
 290         if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
 291             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
 292             if(ratio) {
 293                 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
 294                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 295                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 296                 ratio = ProtectHeap(ratio);
 297                 ratio->Write();
 298                 // write histo values to RMS files if both routines converged
 299                 // input values are weighted by their uncertainty
 300                 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
 301                     if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
 302                     if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
 303                     if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
 304                }
 305             }
 306             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
 307             if(v2) {
 308                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 309                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 310                 v2->GetYaxis()->SetTitle("v_{2}");
 311                 v2 = ProtectHeap(v2);
 312                 v2->Write();
 313             }
 314         } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
 315             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
 316             if(ratio) {
 317                 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
 318                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 319                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 320                 ratio = ProtectHeap(ratio);
 321                 ratio->Write();
 322             }
 323             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
 324              if(v2) {
 325                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 326                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 327                 v2->GetYaxis()->SetTitle("v_{2}");
 328                 v2 = ProtectHeap(v2);
 329                 v2->Write();
 330             }
 331         }
 332     }   // end of if(fDphiUnfolding)
 333     fDeltaPtDeltaPhi->Write();
 334     unfoldedJetSpectrumIn->Sumw2();
 335     ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
 336     unfoldedJetSpectrumIn->Write();
 337     unfoldedJetSpectrumOut->Sumw2();
 338     ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
 339     unfoldedJetSpectrumOut->Write();
 340     fJetPtDeltaPhi->Write();
 341     // save the current state of the unfolding object
 342     SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
 343     TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
 344     TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
 345     unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
 346     unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
 347     unfoldedJetSpectrumInForSub->Write();
 348
 349 }
 350 //_____________________________________________________________________________
 351 TH1D* AliJetFlowTools::UnfoldWrapper(
 352         const TH1D* measuredJetSpectrum,        // truncated raw jets (same binning as pt rec of response)
 353         const TH2D* resizedResponse,            // response matrix
 354         const TH1D* kinematicEfficiency,        // kinematic efficiency
 355         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
 356         const TString suffix,                   // suffix (in or out of plane)
 357         const TH1D* jetFindingEfficiency)       // jet finding efficiency
 358 {
 359     // wrapper function to call specific unfolding routine
 360     TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
 361     // initialize functon pointer
 362     if(fUnfoldingAlgorithm == kChi2)                    myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
 363     else if(fUnfoldingAlgorithm == kBayesian)           myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
 364     else if(fUnfoldingAlgorithm == kBayesianAli)        myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
 365     else if(fUnfoldingAlgorithm == kSVD)                myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
 366     else if(fUnfoldingAlgorithm == kNone) {
 367         TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
 368         clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
 369         return clone;
 370     }
 371     else return 0x0;
 372     // do the actual unfolding with the selected function
 373     return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
 374 }
 375 //_____________________________________________________________________________
 376 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
 377         const TH1D* measuredJetSpectrum,      // truncated raw jets (same binning as pt rec of response)
 378         const TH2D* resizedResponse,          // response matrix
 379         const TH1D* kinematicEfficiency,      // kinematic efficiency
 380         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
 381         const TString suffix,                 // suffix (in or out of plane)
 382         const TH1D* jetFindingEfficiency)     // jet finding efficiency (optional)
 383 {
 384     // unfold the spectrum using chi2 minimization
 385
 386     // step 0) setup the static members of AliUnfolding
 387     ResetAliUnfolding();                // reset from previous iteration
 388                                         // also deletes and re-creates the global TVirtualFitter
 389     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
 390     if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 391     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 392     if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 393     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 394     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
 395
 396     // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
 397     // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
 398     // denoted by the suffix 'Local'
 399
 400     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
 401     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
 402     // unfolded local will be filled with the result of the unfolding
 403     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
 404
 405     // full response matrix and kinematic efficiency
 406     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
 407     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
 408
 409     // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
 410     TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
 411     // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
 412     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
 413
 414     // step 2) start the unfolding
 415     Int_t status(-1), i(0);
 416     while(status < 0 && i < 100) {
 417         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
 418         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
 419         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
 420         status = AliUnfolding::Unfold(
 421                 resizedResponseLocal,           // response matrix
 422                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
 423                 measuredJetSpectrumLocal,              // measured spectrum
 424                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
 425                 unfoldedLocal);                 // results
 426         // status holds the minuit fit status (where 0 means convergence)
 427         i++;
 428     }
 429     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
 430     if(status == 0 && gMinuit->fISW[1] == 3) {
 431         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
 432         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
 433         if(gMinuit) gMinuit->Command("SET COV");
 434         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
 435         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
 436         pearson->Print();
 437         TH2D *hPearson(new TH2D(*pearson));
 438         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
 439         hPearson = ProtectHeap(hPearson);
 440         hPearson->Write();
 441     } else status = -1;
 442
 443     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
 444     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
 445     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
 446     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
 447     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
 448     if(ratio) {
 449         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
 450         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 451         ratio = ProtectHeap(ratio);
 452         ratio->Write();
 453     }
 454
 455     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
 456     // objects are cloned using 'ProtectHeap()'
 457     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
 458     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 459     measuredJetSpectrumLocal->Write();
 460
 461     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
 462     resizedResponseLocal->Write();
 463
 464     unfoldedLocal = ProtectHeap(unfoldedLocal);
 465     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
 466     unfoldedLocal->Write();
 467
 468     foldedLocal = ProtectHeap(foldedLocal);
 469     foldedLocal->Write();
 470
 471     priorLocal = ProtectHeap(priorLocal);
 472     priorLocal->Write();
 473
 474     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
 475     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));
 476     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
 477     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
 478     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
 479     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
 480     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
 481     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
 482     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
 483     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
 484     fitStatus->Write();
 485     return unfoldedLocal;
 486 }
 487 //_____________________________________________________________________________
 488 TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
 489         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 490         const TH2D* resizedResponse,                   // full response matrix, normalized in slides of pt true
 491         const TH1D* kinematicEfficiency,              // kinematic efficiency
 492         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 493         const TString suffix,                         // suffix (in, out)
 494         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 495 {
 496
 497     TH1D* priorLocal( GetPrior(
 498         measuredJetSpectrum,              // jet pt in pt rec bins
 499         resizedResponse,                  // full response matrix, normalized in slides of pt true
 500         kinematicEfficiency,              // kinematic efficiency
 501         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 502         suffix,                           // suffix (in, out)
 503         jetFindingEfficiency));           // jet finding efficiency (optional)
 504     if(!priorLocal) {
 505         printf(" > couldn't find prior ! < \n");
 506         return 0x0;
 507     } else printf(" 1) retrieved prior \n");
 508
 509     // go back to the 'root' directory of this instance of the SVD unfolding routine
 510     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 511
 512     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 513     // measured jets in pt rec binning
 514     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
 515     // local copie of the response matrix
 516     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
 517     // local copy of response matrix, all true slides normalized to 1
 518     // this response matrix will eventually be used in the re-folding routine
 519     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
 520     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
 521     // kinematic efficiency
 522     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 523     // place holder histos
 524     TH1D *unfoldedLocalSVD(0x0);
 525     TH1D *foldedLocalSVD(0x0);
 526     cout << " 2) setup necessary input " << endl;
 527     // 3) configure routine
 528     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 529     cout << " step 3) configured routine " << endl;
 530
 531     // 4) get transpose matrices
 532     // a) get the transpose of the full response matrix
 533     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
 534     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 535     // normalize it with the prior. this will ensure that high statistics bins will constrain the
 536     // end result most strenuously than bins with limited number of counts
 537     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
 538     cout << " 4) retrieved first transpose matrix " << endl;
 539
 540     // 5) get response for SVD unfolding
 541     RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
 542     cout << " 5) retrieved roo unfold response object " << endl;
 543
 544     // 6) actualy unfolding loop
 545     RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 546     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
 547     // correct the spectrum for the kinematic efficiency
 548     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
 549
 550     // get the pearson coefficients from the covariance matrix
 551     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
 552     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 553     if(pearson) {
 554         TH2D* hPearson(new TH2D(*pearson));
 555         pearson->Print();
 556         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 557         hPearson = ProtectHeap(hPearson);
 558         hPearson->Write();
 559     } else return 0x0;       // return if unfolding didn't converge
 560
 561     // plot singular values and d_i vector
 562     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
 563     TH1* hSVal(svdUnfold->GetSV());
 564     TH1D* hdi(svdUnfold->GetD());
 565     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
 566     hSVal->SetXTitle("singular values");
 567     hSVal->Write();
 568     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
 569     hdi->SetXTitle("|d_{i}^{kreg}|");
 570     hdi->Write();
 571     cout << " plotted singular values and d_i vector " << endl;
 572
 573     // 7) refold the unfolded spectrum
 574     foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
 575     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
 576     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
 577     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 578     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 579     ratio->Write();
 580     cout << " 7) refolded the unfolded spectrum " << endl;
 581
 582     // write the measured, unfolded and re-folded spectra to the output directory
 583     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 584     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 585     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 586     measuredJetSpectrumLocal->Write(); // input spectrum
 587     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 588     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
 589     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
 590     unfoldedLocalSVD->Write();  // unfolded spectrum
 591     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 592     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
 593     foldedLocalSVD->Write();    // re-folded spectrum
 594
 595     // save more general bookkeeeping histograms to the output directory
 596     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
 597     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 598     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 599     responseMatrixLocalTransposePrior->Write();
 600     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
 601     priorLocal->SetXTitle("p_{t} [GeV/c]");
 602     priorLocal = ProtectHeap(priorLocal);
 603     priorLocal->Write();
 604     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
 605     resizedResponseLocalNorm->Write();
 606
 607     // save some info
 608     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));
 609     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 610     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
 611     fitStatus->Write();
 612
 613     return unfoldedLocalSVD;
 614 }
 615 //_____________________________________________________________________________
 616 TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
 617         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 618         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
 619         const TH1D* kinematicEfficiency,              // kinematic efficiency
 620         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 621         const TString suffix,                         // suffix (in, out)
 622         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 623 {
 624     // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
 625     // FIXME careful, not tested yet ! (06122013) FIXME
 626
 627     // step 0) setup the static members of AliUnfolding
 628     ResetAliUnfolding();                // reset from previous iteration
 629                                         // also deletes and re-creates the global TVirtualFitter
 630     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
 631     if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
 632     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
 633     else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
 634     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
 635     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
 636
 637     // 1) get a prior for unfolding and clone all the input histograms
 638     TH1D* priorLocal( GetPrior(
 639         measuredJetSpectrum,              // jet pt in pt rec bins
 640         resizedResponse,                  // full response matrix, normalized in slides of pt true
 641         kinematicEfficiency,              // kinematic efficiency
 642         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 643         suffix,                           // suffix (in, out)
 644         jetFindingEfficiency));           // jet finding efficiency (optional)
 645     if(!priorLocal) {
 646         printf(" > couldn't find prior ! < \n");
 647         return 0x0;
 648     } else printf(" 1) retrieved prior \n");
 649     // switch back to root dir of this unfolding procedure
 650     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 651
 652     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
 653     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
 654     // unfolded local will be filled with the result of the unfolding
 655     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
 656
 657     // full response matrix and kinematic efficiency
 658     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
 659     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
 660
 661     // step 2) start the unfolding
 662     Int_t status(-1), i(0);
 663     while(status < 0 && i < 100) {
 664         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
 665         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
 666         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
 667         status = AliUnfolding::Unfold(
 668                 resizedResponseLocal,           // response matrix
 669                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
 670                 measuredJetSpectrumLocal,              // measured spectrum
 671                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
 672                 unfoldedLocal);                 // results
 673         // status holds the minuit fit status (where 0 means convergence)
 674         i++;
 675     }
 676     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
 677     if(status == 0 && gMinuit->fISW[1] == 3) {
 678         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
 679         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
 680         if(gMinuit) gMinuit->Command("SET COV");
 681         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
 682         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
 683         pearson->Print();
 684         TH2D *hPearson(new TH2D(*pearson));
 685         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
 686         hPearson = ProtectHeap(hPearson);
 687         hPearson->Write();
 688     } else status = -1;
 689
 690     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
 691     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
 692     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
 693     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
 694     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
 695     if(ratio) {
 696         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
 697         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 698         ratio = ProtectHeap(ratio);
 699         ratio->Write();
 700     }
 701
 702     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
 703     // objects are cloned using 'ProtectHeap()'
 704     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
 705     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 706     measuredJetSpectrumLocal->Write();
 707
 708     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
 709     resizedResponseLocal->Write();
 710
 711     unfoldedLocal = ProtectHeap(unfoldedLocal);
 712     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
 713     unfoldedLocal->Write();
 714
 715     foldedLocal = ProtectHeap(foldedLocal);
 716     foldedLocal->Write();
 717
 718     priorLocal = ProtectHeap(priorLocal);
 719     priorLocal->Write();
 720
 721     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
 722     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));
 723     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
 724     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
 725     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
 726     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
 727     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
 728     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
 729     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
 730     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
 731     fitStatus->Write();
 732     return unfoldedLocal;
 733 }
 734 //_____________________________________________________________________________
 735 TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
 736         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 737         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
 738         const TH1D* kinematicEfficiency,              // kinematic efficiency
 739         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 740         const TString suffix,                         // suffix (in, out)
 741         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 742 {
 743     // use bayesian unfolding from the RooUnfold package to unfold jet spectra
 744
 745     // 1) get a prior for unfolding.
 746     TH1D* priorLocal( GetPrior(
 747         measuredJetSpectrum,              // jet pt in pt rec bins
 748         resizedResponse,                  // full response matrix, normalized in slides of pt true
 749         kinematicEfficiency,              // kinematic efficiency
 750         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 751         suffix,                           // suffix (in, out)
 752         jetFindingEfficiency));           // jet finding efficiency (optional)
 753     if(!priorLocal) {
 754         printf(" > couldn't find prior ! < \n");
 755         return 0x0;
 756     } else printf(" 1) retrieved prior \n");
 757     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 758
 759     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 760     // measured jets in pt rec binning
 761     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
 762     // local copie of the response matrix
 763     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
 764     // local copy of response matrix, all true slides normalized to 1
 765     // this response matrix will eventually be used in the re-folding routine
 766     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
 767     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
 768     // kinematic efficiency
 769     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 770     // place holder histos
 771     TH1D *unfoldedLocalBayes(0x0);
 772     TH1D *foldedLocalBayes(0x0);
 773     cout << " 2) setup necessary input " << endl;
 774     // 4) get transpose matrices
 775     // a) get the transpose of the full response matrix
 776     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
 777     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 778     // normalize it with the prior. this will ensure that high statistics bins will constrain the
 779     // end result most strenuously than bins with limited number of counts
 780     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
 781     // 3) get response for Bayesian unfolding
 782     RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
 783
 784     // 4) actualy unfolding loop
 785     RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
 786     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 787     unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
 788     // correct the spectrum for the kinematic efficiency
 789     unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
 790     // get the pearson coefficients from the covariance matrix
 791     TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
 792     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 793     if(pearson) {
 794         TH2D* hPearson(new TH2D(*pearson));
 795         pearson->Print();
 796         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 797         hPearson = ProtectHeap(hPearson);
 798         hPearson->Write();
 799     } else return 0x0;       // return if unfolding didn't converge
 800
 801     // 5) refold the unfolded spectrum
 802     foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
 803     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio  measured / re-folded", kTRUE));
 804     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
 805     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 806     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 807     ratio->Write();
 808     cout << " 7) refolded the unfolded spectrum " << endl;
 809
 810     // write the measured, unfolded and re-folded spectra to the output directory
 811     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 812     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 813     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 814     measuredJetSpectrumLocal->Write(); // input spectrum
 815     unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 816     unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
 817     if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
 818     unfoldedLocalBayes->Write();  // unfolded spectrum
 819     foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 820     foldedLocalBayes = ProtectHeap(foldedLocalBayes);
 821     foldedLocalBayes->Write();    // re-folded spectrum
 822
 823     // save more general bookkeeeping histograms to the output directory
 824     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
 825     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 826     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 827     responseMatrixLocalTransposePrior->Write();
 828     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
 829     priorLocal->SetXTitle("p_{t} [GeV/c]");
 830     priorLocal = ProtectHeap(priorLocal);
 831     priorLocal->Write();
 832     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
 833     resizedResponseLocalNorm->Write();
 834
 835     // save some info
 836     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));
 837     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
 838     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
 839     fitStatus->Write();
 840
 841     return  unfoldedLocalBayes;
 842 }
 843 //_____________________________________________________________________________
 844 Bool_t AliJetFlowTools::PrepareForUnfolding()
 845 {
 846     // prepare for unfolding
 847     if(fRawInputProvided) return kTRUE;
 848     if(!fInputList) {
 849         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
 850         return kFALSE;
 851     }
 852     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
 853     // check if the pt bin for true and rec have been set
 854     if(!fBinsTrue || !fBinsRec) {
 855         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
 856         return kFALSE;
 857     }
 858     if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
 859                           // procedures, these profiles will be nonsensical, user is responsible
 860         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 861         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 862         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 863     }
 864     if(!fTrainPower) {
 865         // clear minuit state to avoid constraining the fit with the results of the previous iteration
 866         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 867     }
 868     // extract the spectra
 869     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
 870     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
 871     if(!fJetPtDeltaPhi) {
 872         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
 873         return kFALSE;
 874     }
 875     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
 876     // in plane spectrum
 877     if(!fDphiUnfolding) {
 878         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
 879         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
 880     } else {
 881         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
 882         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
 883         fSpectrumIn = ProtectHeap(fSpectrumIn);
 884         // out of plane spectrum
 885         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
 886         fSpectrumOut = ProtectHeap(fSpectrumOut);
 887     }
 888     // normalize spectra to event count if requested
 889     if(fNormalizeSpectra) {
 890         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
 891         if(!rho) return 0x0;
 892         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
 893         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
 894         if(fEventCount > 0) {
 895             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
 896             fSpectrumOut->Sumw2();
 897             Double_t pt(0);
 898             Double_t error(0); // lots of issues with the errors here ...
 899             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
 900                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
 901                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
 902                 fSpectrumIn->SetBinContent(1+i, pt);
 903                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
 904                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
 905                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
 906             }
 907             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
 908                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
 909                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
 910                 fSpectrumOut->SetBinContent(1+i, pt);
 911                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
 912                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
 913                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
 914             }
 915         }
 916         if(normalizeToFullSpectrum) fEventCount = -1;
 917     }
 918     // extract the delta pt matrices
 919     TString deltaptName("");
 920     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
 921     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
 922     if(!fDeltaPtDeltaPhi) {
 923         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
 924         printf(" > may be ok, depending no what you want to do < \n");
 925         fRefreshInput = kTRUE;
 926         return kTRUE;
 927     }
 928     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
 929     // in plane delta pt distribution
 930     if(!fDphiUnfolding) {
 931         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
 932         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
 933     } else {
 934         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
 935         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
 936         // out of plane delta pt distribution
 937         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
 938         fDptInDist = ProtectHeap(fDptInDist);
 939         fDptOutDist = ProtectHeap(fDptOutDist);
 940         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
 941     }
 942
 943     // create a rec - true smeared response matrix
 944     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
 945     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 946         Bool_t skip = kFALSE;
 947         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 948             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
 949             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 950         }
 951     }
 952     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
 953     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 954         Bool_t skip = kFALSE;
 955         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 956             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
 957             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 958         }
 959     }
 960     fDptIn = new TH2D(*rfIn);
 961     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
 962     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 963     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 964     fDptIn = ProtectHeap(fDptIn);
 965     fDptOut = new TH2D(*rfOut);
 966     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
 967     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 968     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 969     fDptOut = ProtectHeap(fDptOut);
 970
 971     fRefreshInput = kTRUE;     // force cloning of the input
 972     return kTRUE;
 973 }
 974 //_____________________________________________________________________________
 975 Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
 976     // prepare for unfoldingUA - more robust method to extract input spectra from file
 977     // will replace PrepareForUnfolding eventually (09012014)
 978     if(!fInputList) {
 979         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
 980         return kFALSE;
 981     }
 982     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
 983     // check if the pt bin for true and rec have been set
 984     if(!fBinsTrue || !fBinsRec) {
 985         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
 986         return kFALSE;
 987     }
 988     if(!fTrainPower) {
 989         // clear minuit state to avoid constraining the fit with the results of the previous iteration
 990         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 991     }
 992     // extract the spectra
 993     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
 994     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
 995     if(!fJetPtDeltaPhi) {
 996         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
 997         return kFALSE;
 998     }
 999     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1000     // in plane spectrum
1001     fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1002     // extract the delta pt matrices
1003     TString deltaptName("");
1004     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
1005     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1006     if(!fDeltaPtDeltaPhi) {
1007         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1008         printf(" > may be ok, depending no what you want to do < \n");
1009         fRefreshInput = kTRUE;
1010         return kTRUE;
1011     }
1012     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1013     // in plane delta pt distribution
1014     fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1015     // create a rec - true smeared response matrix
1016     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1017     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1018         Bool_t skip = kFALSE;
1019         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1020             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1021             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1022         }
1023     }
1024     fDptIn = new TH2D(*rfIn);
1025     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1026     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1027     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1028     fDptIn = ProtectHeap(fDptIn);
1029
1030     return kTRUE;
1031 }
1032 //_____________________________________________________________________________
1033 TH1D* AliJetFlowTools::GetPrior(
1034         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
1035         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
1036         const TH1D* kinematicEfficiency,              // kinematic efficiency
1037         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
1038         const TString suffix,                         // suffix (in, out)
1039         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
1040 {
1041     // 1) get a prior for unfolding.
1042     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1043     TH1D* unfolded(0x0);
1044     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1045     dirOut->cd();
1046     switch (fPrior) {    // select the prior for unfolding
1047         case kPriorChi2 : {
1048             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
1049                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original binning
1050                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
1051                 TArrayD* tempArrayRec(fBinsRec);
1052                 fBinsRec = fBinsRecPrior;
1053                 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1054                 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1055                 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1056                 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1057                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1058                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1059                 unfolded= UnfoldSpectrumChi2(
1060                             measuredJetSpectrumChi2,
1061                             resizedResponseChi2,
1062                             kinematicEfficiencyChi2,
1063                             measuredJetSpectrumTrueBinsChi2,    // prior for chi2 unfolding (measured)
1064                             TString(Form("prior_%s", suffix.Data())));
1065                if(!unfolded) {
1066                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1067                     printf("               probably Chi2 unfolding did not converge < \n");
1068                     return 0x0;
1069                 }
1070                 fBinsTrue = tempArrayTrue;  // reset bins borders
1071                 fBinsRec = tempArrayRec;
1072                 // if the chi2 prior has a different binning, rebin to the true binning for the  unfolding
1073                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
1074             } else {
1075                 unfolded = UnfoldSpectrumChi2(
1076                             measuredJetSpectrum,
1077                             resizedResponse,
1078                             kinematicEfficiency,
1079                             measuredJetSpectrumTrueBins,        // prior for chi2 unfolding (measured)
1080                             TString(Form("prior_%s", suffix.Data())));
1081                 if(!unfolded) {
1082                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1083                     printf("               probably Chi2 unfolding did not converge < \n");
1084                     return 0x0;
1085                 }
1086             }
1087             break;
1088         }
1089         case kPriorMeasured : {
1090             unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
1091         }
1092         default : break;
1093     }
1094     // it can be important that the prior is smooth, this can be achieved by
1095     // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1096     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1097     TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1098     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1099     return priorLocal;
1100 }
1101 //_____________________________________________________________________________
1102 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1103     // resize the x-axis of a th1d
1104     if(!histo) {
1105         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1106         return NULL;
1107     }
1108     // see how many bins we need to copy
1109     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);
1110     // low is the bin number of the first new bin
1111     Int_t l = histo->GetXaxis()->FindBin(low);
1112     // set the values
1113     Double_t _x(0), _xx(0);
1114     for(Int_t i(0); i < up-low; i++) {
1115         _x = histo->GetBinContent(l+i);
1116         _xx=histo->GetBinError(l+i);
1117         resized->SetBinContent(i+1, _x);
1118         resized->SetBinError(i+1, _xx);
1119     }
1120     return resized;
1121 }
1122 //_____________________________________________________________________________
1123 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1124     // resize the y-axis of a th2d
1125     if(!histo) {
1126         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1127         return NULL;
1128     }
1129     // see how many bins we need to copy
1130     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());
1131     // assume only the y-axis has changed
1132     // low is the bin number of the first new bin
1133     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1134     // set the values
1135     Double_t _x(0), _xx(0);
1136     for(Int_t i(0); i < x->GetSize(); i++) {
1137         for(Int_t j(0); j < y->GetSize(); j++) {
1138             _x = histo->GetBinContent(i, low+j);
1139             _xx=histo->GetBinError(i, low+1+j);
1140             resized->SetBinContent(i, j, _x);
1141             resized->SetBinError(i, j, _xx);
1142         }
1143     }
1144     return resized;
1145 }
1146 //_____________________________________________________________________________
1147 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1148     // general method to normalize all vertical slices of a th2 to unity
1149     // i.e. get a probability matrix
1150     if(!histo) {
1151         printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1152         return NULL;
1153     }
1154     Int_t binsX = histo->GetXaxis()->GetNbins();
1155     Int_t binsY = histo->GetYaxis()->GetNbins();
1156
1157     // normalize all slices in x
1158     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1159         Double_t weight = 0;
1160         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1161             weight+=histo->GetBinContent(i+1, j+1);
1162         }       // now we know the total weight
1163         for(Int_t j(0); j < binsY; j++) {
1164             if (weight <= 0 ) continue;
1165             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1166             if(noError) histo->SetBinError(  1+i, j+1, 0.);
1167             else histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1168         }
1169     }
1170     return histo;
1171 }
1172 //_____________________________________________________________________________
1173 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1174     // return a TH1D with the supplied histogram rebinned to the supplied bins
1175     // the returned histogram is new, the original is deleted from the heap if kill is true
1176     if(!histo || !bins) {
1177         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1178         return NULL;
1179     }
1180     // create the output histo
1181     TString name = histo->GetName();
1182     name+="_template";
1183     name+=suffix;
1184     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1185     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1186         // loop over the bins of the old histo and fill the new one with its data
1187         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1188     }
1189     if(kill) delete histo;
1190     return rebinned;
1191 }
1192 //_____________________________________________________________________________
1193 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1194     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1195     if(!fResponseMaker || !binsTrue || !binsRec) {
1196         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1197         return 0x0;
1198     }
1199     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1200     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1201 }
1202 //_____________________________________________________________________________
1203 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1204 {
1205     // multiply two matrices
1206     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1207     TH2D* c = (TH2D*)a->Clone("c");
1208     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1209         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1210             Double_t val = 0;
1211             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1212                 Int_t y2 = x1;
1213                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1214             }
1215             c->SetBinContent(x2, y1, val);
1216             c->SetBinError(x2, y1, 0.);
1217         }
1218     }
1219     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1220     return c;
1221 }
1222 //_____________________________________________________________________________
1223 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1224 {
1225     // normalize a th1d to a certain scale
1226     histo->Sumw2();
1227     Double_t integral = histo->Integral()*scale;
1228     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1229     else if (scale != 1.) histo->Scale(1./scale, "width");
1230     else printf(" > Histogram integral < 0, cannot normalize \n");
1231     return histo;
1232 }
1233 //_____________________________________________________________________________
1234 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1235 {
1236     // Calculate pearson coefficients from covariance matrix
1237     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1238     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1239     Double_t pearson(0.);
1240     if(nrows==0 && ncols==0) return 0x0;
1241     for(Int_t row = 0; row < nrows; row++) {
1242         for(Int_t col = 0; col<ncols; col++) {
1243         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1244         (*pearsonCoefficients)(row,col) = pearson;
1245         }
1246     }
1247     return pearsonCoefficients;
1248 }
1249 //_____________________________________________________________________________
1250 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1251     // smoothen the spectrum using a user defined function
1252     // returns a clone of the original spectrum if fitting failed
1253     // if kill is kTRUE the input spectrum will be deleted from the heap
1254     // if 'count' is selected, bins are filled with integers (necessary if the
1255     // histogram is interpreted in a routine which accepts only counts)
1256     if(!spectrum || !function) return 0x0;
1257     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1258     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1259     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1260     TFitResultPtr r = temp->Fit(function, "", "", min, max);
1261     if((int)r == 0) {   // MINUIT status
1262         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1263             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1264                 (counts) ? temp->SetBinContent(i, (int)(function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i))) : temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1265                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1266             }
1267         }
1268     }
1269     if(kill) delete spectrum;
1270     return temp;
1271 }
1272 //_____________________________________________________________________________
1273 void AliJetFlowTools::Style(TCanvas* c, TString style)
1274 {
1275     // set a default style for a canvas
1276     if(!strcmp(style.Data(), "PEARSON")) {
1277         printf(" > style PEARSON canvas < \n");
1278         gStyle->SetOptStat(0);
1279         c->SetGridx();
1280         c->SetGridy();
1281         c->SetTicks();
1282         return;
1283     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1284         printf(" > style SPECTRUM canvas < \n");
1285         gStyle->SetOptStat(0);
1286         c->SetLogy();
1287         c->SetGridx();
1288         c->SetGridy();
1289         c->SetTicks();
1290         return;
1291     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1292 }
1293 //_____________________________________________________________________________
1294 void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1295 {
1296     // set a default style for a canvas
1297     if(!strcmp(style.Data(), "PEARSON")) {
1298         printf(" > style PEARSON pad < \n");
1299         gStyle->SetOptStat(0);
1300         c->SetGridx();
1301         c->SetGridy();
1302         c->SetTicks();
1303         return;
1304     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1305         printf(" > style SPECTRUM pad < \n");
1306         gStyle->SetOptStat(0);
1307         c->SetLogy();
1308         c->SetGridx();
1309         c->SetGridy();
1310         c->SetTicks();
1311         return;
1312     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1313 }
1314 //_____________________________________________________________________________
1315 void AliJetFlowTools::Style(TLegend* l)
1316 {
1317     // set a default style for a legend
1318     l->SetTextSize(.06);
1319     l->SetFillColor(0);
1320 }
1321 //_____________________________________________________________________________
1322 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1323 {
1324     // style a histo
1325     h->SetLineColor(col);
1326     h->SetMarkerColor(col);
1327     h->SetLineWidth(2.);
1328     h->SetMarkerSize(2.);
1329     h->SetTitleSize(0.06);
1330     h->GetYaxis()->SetLabelSize(0.06);
1331     h->GetXaxis()->SetLabelSize(0.06);
1332     h->GetYaxis()->SetTitleSize(0.06);
1333     h->GetXaxis()->SetTitleSize(0.06);
1334     h->GetYaxis()->SetTitleOffset(1.);
1335     h->GetXaxis()->SetTitleOffset(.9);
1336     switch (type) {
1337         case kInPlaneSpectrum : {
1338             h->SetTitle("IN PLANE");
1339             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1340             h->GetYaxis()->SetTitle("counts");
1341         } break;
1342         case kOutPlaneSpectrum : {
1343             h->SetTitle("OUT OF PLANE");
1344             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1345             h->GetYaxis()->SetTitle("counts");
1346        } break;
1347        case kUnfoldedSpectrum : {
1348             h->SetTitle("UNFOLDED");
1349             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1350             h->GetYaxis()->SetTitle("counts");
1351        } break;
1352        case kFoldedSpectrum : {
1353             h->SetTitle("FOLDED");
1354             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1355             h->GetYaxis()->SetTitle("counts");
1356        } break;
1357        case kMeasuredSpectrum : {
1358             h->SetTitle("MEASURED");
1359             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1360             h->GetYaxis()->SetTitle("counts");
1361        } break;
1362        default : break;
1363     }
1364 }
1365 //_____________________________________________________________________________
1366 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1367 {
1368     // style a histo
1369     h->SetLineColor(col);
1370     h->SetMarkerColor(col);
1371     h->SetLineWidth(2.);
1372     h->SetMarkerSize(2.);
1373     h->GetYaxis()->SetLabelSize(0.06);
1374     h->GetXaxis()->SetLabelSize(0.06);
1375     h->GetYaxis()->SetTitleSize(0.06);
1376     h->GetXaxis()->SetTitleSize(0.06);
1377     h->GetYaxis()->SetTitleOffset(1.);
1378     h->GetXaxis()->SetTitleOffset(.9);
1379     switch (type) {
1380         case kInPlaneSpectrum : {
1381             h->SetTitle("IN PLANE");
1382             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1383             h->GetYaxis()->SetTitle("counts");
1384         } break;
1385         case kOutPlaneSpectrum : {
1386             h->SetTitle("OUT OF PLANE");
1387             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1388             h->GetYaxis()->SetTitle("counts");
1389        } break;
1390        case kUnfoldedSpectrum : {
1391             h->SetTitle("UNFOLDED");
1392             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1393             h->GetYaxis()->SetTitle("counts");
1394        } break;
1395        case kFoldedSpectrum : {
1396             h->SetTitle("FOLDED");
1397             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1398             h->GetYaxis()->SetTitle("counts");
1399        } break;
1400        case kMeasuredSpectrum : {
1401             h->SetTitle("MEASURED");
1402             h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1403             h->GetYaxis()->SetTitle("counts");
1404        } break;
1405        default : break;
1406     }
1407 }
1408 //_____________________________________________________________________________
1409 void AliJetFlowTools::DoSystematics(Int_t nominalIn, Int_t nominalOut, TArrayI* variationsIn, TArrayI* variationsOut, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
1410 {
1411    // do systematics. structure of this function is similar to that of PostProcess, however
1412    // in this case a nomial index can be supplied as well as an array of systematic checks
1413    //
1414    // FIXME make separate canvas for nominal value
1415    //
1416    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
1417    TFile readMe(in.Data(), "READ");     // open file read-only
1418    if(readMe.IsZombie()) {
1419        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1420        return;
1421    }
1422    printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1423    readMe.ls();
1424    TList* listOfKeys((TList*)readMe.GetListOfKeys());
1425    if(!listOfKeys) {
1426        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1427        return;
1428    }
1429    // check input params
1430    if(variationsIn->GetSize() != variationsOut->GetSize()) {
1431        printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
1432        return;
1433    }
1434    TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(nominalIn)->GetName())));
1435    TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(nominalOut)->GetName())));
1436    if(!(defRootDirIn && defRootDirOut)) {
1437        printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
1438        return;
1439    }
1440    TString defIn(defRootDirIn->GetName());
1441    TString defOut(defRootDirOut->GetName());
1442    // prepare necessary canvasses
1443    TCanvas* canvasIn(new TCanvas("SYSTPearsonIn", "SYSTPearsonIn"));
1444    TCanvas* canvasOut(0x0);
1445    if(fDphiUnfolding) canvasOut = new TCanvas("SYSTPearsonOut", "SYSTPearsonOut");
1446    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("SYSTRefoldedIn", "SYSTRefoldedIn"));
1447    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
1448    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("SYSTRefoldedOut", "SYSTRefoldedOut");
1449    TCanvas* canvasSpectraIn(new TCanvas("SYSTSpectraIn", "SYSTSpectraIn"));
1450    TCanvas* canvasSpectraOut(0x0);
1451    if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SYSTSpectraOut", "SYSTSpectraOut");
1452    TCanvas* canvasRatio(0x0);
1453    if(fDphiUnfolding) canvasRatio = new TCanvas("SYSTRatio", "SYSTRatio");
1454    TCanvas* canvasV2(0x0);
1455    if(fDphiUnfolding) canvasV2 = new TCanvas("SYSTV2", "SYSTV2");
1456    TCanvas* canvasMISC(new TCanvas("SYSTMISC", "SYSTMISC"));
1457    TCanvas* canvasMasterIn(new TCanvas("SYSTdefaultIn", "SYSTdefaultIn"));
1458    TCanvas* canvasMasterOut(0x0);
1459    if(fDphiUnfolding) canvasMasterOut = new TCanvas("SYSTdefaultOut", "SYSTdefaultOut");
1460    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
1461
1462    TCanvas* canvasProfiles(new TCanvas("SYSTcanvasProfiles", "SYSTcanvasProfiles"));
1463    canvasProfiles->Divide(2);
1464    TProfile* ratioProfile(new TProfile("SYSTratioProfile", "SYSTratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1465    TProfile* v2Profile(new TProfile("SYSTv2Profile", "SYSTv2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1466    // get an estimate of the number of outputs and find the default set
1467    Int_t rows(TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
1468    canvasIn->Divide(columns, rows);
1469    if(canvasOut) canvasOut->Divide(columns, rows);
1470    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1471    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1472    canvasSpectraIn->Divide(columns, rows);
1473    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
1474    if(canvasRatio) canvasRatio->Divide(columns, rows);
1475    if(canvasV2) canvasV2->Divide(columns, rows);
1476    canvasMasterIn->Divide(columns, rows);
1477    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
1478
1479    // prepare a separate set of canvases to hold the nominal points
1480    TCanvas* canvasNominalIn(new TCanvas("NominalPearsonIn", "NominalPearsonIn"));
1481    TCanvas* canvasNominalOut(0x0);
1482    if(fDphiUnfolding) canvasNominalOut = new TCanvas("NominalPearsonOut", "NominalPearsonOut");
1483    TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas("NominalRefoldedIn", "NominalRefoldedIn"));
1484    TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
1485    if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas("NominalRefoldedOut", "NominalRefoldedOut");
1486    TCanvas* canvasNominalSpectraIn(new TCanvas("NominalSpectraIn", "NominalSpectraIn"));
1487    TCanvas* canvasNominalSpectraOut(0x0);
1488    if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas("NominalSpectraOut", "NominalSpectraOut");
1489    TCanvas* canvasNominalRatio(0x0);
1490    if(fDphiUnfolding) canvasNominalRatio = new TCanvas("NominalRatio", "NominalRatio");
1491    TCanvas* canvasNominalV2(0x0);
1492    if(fDphiUnfolding) canvasNominalV2 = new TCanvas("NominalV2", "NominalV2");
1493    TCanvas* canvasNominalMISC(new TCanvas("NominalMISC", "NominalMISC"));
1494    TCanvas* canvasNominalMasterIn(new TCanvas("NominaldefaultIn", "NominaldefaultIn"));
1495    TCanvas* canvasNominalMasterOut(0x0);
1496    if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas("NominaldefaultOut", "NominaldefaultOut");
1497    (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
1498
1499    columns = 1;         rows = 1;
1500    canvasNominalIn->Divide(columns, rows);
1501    if(canvasNominalOut) canvasNominalOut->Divide(columns, rows);
1502    canvasNominalRatioMeasuredRefoldedIn->Divide(columns, rows);
1503    if(canvasNominalRatioMeasuredRefoldedOut) canvasNominalRatioMeasuredRefoldedOut->Divide(columns, rows);
1504    canvasNominalSpectraIn->Divide(columns, rows);
1505    if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(columns, rows);
1506    if(canvasNominalRatio) canvasNominalRatio->Divide(columns, rows);
1507    if(canvasNominalV2) canvasNominalV2->Divide(columns, rows);
1508
1509    canvasNominalMasterIn->Divide(columns, rows);
1510    if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(columns, rows);
1511
1512    // extract the default output
1513    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
1514    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
1515    TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
1516    TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
1517    if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
1518    if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
1519    printf(" > succesfully extracted default results < \n");
1520
1521    // loop through the directories, only plot the graphs if the deconvolution converged
1522    TDirectoryFile* tempDirIn(0x0);
1523    TDirectoryFile* tempDirOut(0x0);
1524    TDirectoryFile* tempIn(0x0);
1525    TDirectoryFile* tempOut(0x0);
1526    TH1D* unfoldedSpectrumInForRatio(0x0);
1527    TH1D* unfoldedSpectrumOutForRatio(0x0);
1528    for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
1529        tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(variationsIn->At(i))->GetName())));
1530        tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(variationsOut->At(i))->GetName())));
1531        if(!(tempDirIn && tempDirOut)) {
1532            printf(" > DoSystematics: couldn't get a set of variations < \n");
1533            continue;
1534        }
1535        TString dirNameIn(tempDirIn->GetName());
1536        TString dirNameOut(tempDirOut->GetName());
1537        // try to read the in- and out of plane subdirs
1538        tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
1539        tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
1540        j++;
1541        if(tempIn) {
1542            // to see if the unfolding converged try to extract the pearson coefficients
1543            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
1544            if(pIn) {
1545                printf(" - %s in plane converged \n", dirNameIn.Data());
1546                canvasIn->cd(j);
1547                if(i==0) canvasNominalIn->cd(j);
1548                Style(gPad, "PEARSON");
1549                pIn->DrawCopy("colz");
1550                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
1551                if(rIn) {
1552                    Style(rIn);
1553                    printf(" > found RatioRefoldedMeasured < \n");
1554                    canvasRatioMeasuredRefoldedIn->cd(j);
1555                    if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
1556                    rIn->SetFillColor(kRed);
1557                    rIn->Draw("ap");
1558                }
1559                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1560                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1561                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
1562                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
1563                if(dvector && avalue && rm && eff) {
1564                    Style(dvector);
1565                    Style(avalue);
1566                    Style(rm);
1567                    Style(eff);
1568                    canvasMISC->cd(1);
1569                    if(i==0) canvasNominalMISC->cd(1);
1570                    Style(gPad, "SPECTRUM");
1571                    dvector->DrawCopy();
1572                    canvasMISC->cd(2);
1573                    if(i==0) canvasNominalMISC->cd(2);
1574                    Style(gPad, "SPECTRUM");
1575                    avalue->DrawCopy();
1576                    canvasMISC->cd(3);
1577                    if(i==0) canvasNominalMISC->cd(3);
1578                    Style(gPad, "PEARSON");
1579                    rm->DrawCopy("colz");
1580                    canvasMISC->cd(4);
1581                    if(i==0) canvasNominalMISC->cd(4);
1582                    eff->DrawCopy();
1583                } else if(rm && eff) {
1584                    Style(rm);
1585                    Style(eff);
1586                    canvasMISC->cd(3);
1587                    if(i==0) canvasNominalMISC->cd(3);
1588                    Style(gPad, "PEARSON");
1589                    rm->DrawCopy("colz");
1590                    canvasMISC->cd(4);
1591                    if(i==0) canvasNominalMISC->cd(4);
1592                    eff->DrawCopy();
1593                }
1594            }
1595            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
1596            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
1597            unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
1598            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
1599            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1600                if(defaultUnfoldedJetSpectrumIn) {
1601                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
1602                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
1603                    temp->Divide(unfoldedSpectrum);
1604                    temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
1605                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1606                    temp->GetYaxis()->SetTitle("ratio");
1607                    canvasMasterIn->cd(j);
1608                    if(i==0) canvasNominalMasterIn->cd(j);
1609                    temp->GetYaxis()->SetRangeUser(0., 2);
1610                    temp->DrawCopy();
1611                }
1612                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
1613                canvasSpectraIn->cd(j);
1614                if(i==0) canvasNominalSpectraIn->cd(j);
1615                Style(gPad);
1616                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1617                unfoldedSpectrum->DrawCopy();
1618                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1619                inputSpectrum->DrawCopy("same");
1620                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1621                refoldedSpectrum->DrawCopy("same");
1622                TLegend* l(AddLegend(gPad));
1623                Style(l);
1624                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1625                    Float_t chi(fitStatus->GetBinContent(1));
1626                    Float_t pen(fitStatus->GetBinContent(2));
1627                    Int_t dof((int)fitStatus->GetBinContent(3));
1628                    Float_t beta(fitStatus->GetBinContent(4));
1629                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1630                } else if (fitStatus) { // only available in SVD fit
1631                    Int_t reg((int)fitStatus->GetBinContent(1));
1632                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1633                }
1634            }
1635        }
1636        if(tempOut) {
1637            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
1638            if(pOut) {
1639                printf(" - %s out of plane converged \n", dirNameOut.Data());
1640                canvasOut->cd(j);
1641                if(i==0) canvasNominalOut->cd(j);
1642                Style(gPad, "PEARSON");
1643                pOut->DrawCopy("colz");
1644                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
1645                if(rOut) {
1646                    Style(rOut);
1647                    printf(" > found RatioRefoldedMeasured < \n");
1648                    canvasRatioMeasuredRefoldedOut->cd(j);
1649                    if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
1650                    rOut->SetFillColor(kRed);
1651                    rOut->Draw("ap");
1652                }
1653                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1654                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1655                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
1656                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
1657                if(dvector && avalue && rm && eff) {
1658                    Style(dvector);
1659                    Style(avalue);
1660                    Style(rm);
1661                    Style(eff);
1662                    canvasMISC->cd(5);
1663                    if(i==0) canvasNominalMISC->cd(5);
1664                    Style(gPad, "SPECTRUM");
1665                    dvector->DrawCopy();
1666                    canvasMISC->cd(6);
1667                    if(i==0) canvasNominalMISC->cd(6);
1668                    Style(gPad, "SPECTRUM");
1669                    avalue->DrawCopy();
1670                    canvasMISC->cd(7);
1671                    if(i==0) canvasNominalMISC->cd(7);
1672                    Style(gPad, "PEARSON");
1673                    rm->DrawCopy("colz");
1674                    canvasMISC->cd(8);
1675                    if(i==0) canvasNominalMISC->cd(8);
1676                    eff->DrawCopy();
1677                } else if(rm && eff) {
1678                    Style(rm);
1679                    Style(eff);
1680                    canvasMISC->cd(7);
1681                    if(i==0) canvasNominalMISC->cd(7);
1682                    Style(gPad, "PEARSON");
1683                    rm->DrawCopy("colz");
1684                    canvasMISC->cd(8);
1685                    if(i==0) canvasNominalMISC->cd(8);
1686                    eff->DrawCopy();
1687                }
1688            }
1689            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
1690            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
1691            unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
1692            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
1693            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1694                if(defaultUnfoldedJetSpectrumOut) {
1695                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
1696                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
1697                    temp->Divide(unfoldedSpectrum);
1698                    temp->SetTitle(Form("ratio nominal [%s] / [%s]", defOut.Data(), dirNameOut.Data()));
1699                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1700                    temp->GetYaxis()->SetTitle("ratio");
1701                    canvasMasterOut->cd(j);
1702                    if(i==0) canvasNominalMasterOut->cd(j);
1703                    temp->GetYaxis()->SetRangeUser(0., 2.);
1704                    temp->DrawCopy();
1705                }
1706                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
1707                canvasSpectraOut->cd(j);
1708                if(i==0) canvasNominalSpectraOut->cd(j);
1709                Style(gPad);
1710                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1711                unfoldedSpectrum->DrawCopy();
1712                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1713                inputSpectrum->DrawCopy("same");
1714                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1715                refoldedSpectrum->DrawCopy("same");
1716                TLegend* l(AddLegend(gPad));
1717                Style(l);
1718                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1719                    Float_t chi(fitStatus->GetBinContent(1));
1720                    Float_t pen(fitStatus->GetBinContent(2));
1721                    Int_t dof((int)fitStatus->GetBinContent(3));
1722                    Float_t beta(fitStatus->GetBinContent(4));
1723                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1724                } else if (fitStatus) { // only available in SVD fit
1725                    Int_t reg((int)fitStatus->GetBinContent(1));
1726                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1727                }
1728            }
1729        }
1730        if(canvasRatio && canvasV2) {
1731            canvasRatio->cd(j);
1732            if(i==0) canvasNominalRatio->cd(j);
1733            TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
1734            Double_t _tempx(0), _tempy(0);
1735            if(ratioYield) {
1736                Style(ratioYield);
1737                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
1738                    ratioYield->GetPoint(b+1, _tempx, _tempy);
1739                    ratioProfile->Fill(_tempx, _tempy);
1740                }
1741                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
1742                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1743                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
1744                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1745                ratioYield->SetFillColor(kRed);
1746                ratioYield->Draw("ap");
1747            }
1748            canvasV2->cd(j);
1749            if(i==0) canvasNominalV2->cd(j);
1750            TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, .7, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
1751            if(ratioV2) {
1752                Style(ratioV2);
1753                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
1754                    ratioV2->GetPoint(b+1, _tempx, _tempy);
1755                    v2Profile->Fill(_tempx, _tempy);
1756
1757                }
1758                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
1759                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1760                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
1761                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1762                ratioV2->SetFillColor(kRed);
1763                ratioV2->Draw("ap");
1764            }
1765        }
1766        delete unfoldedSpectrumInForRatio;
1767        delete unfoldedSpectrumOutForRatio;
1768    }
1769    TFile output(out.Data(), "RECREATE");
1770    // save the canvasses
1771    canvasProfiles->cd(1);
1772    Style(ratioProfile);
1773    ratioProfile->DrawCopy();
1774    canvasProfiles->cd(2);
1775    Style(v2Profile);
1776    v2Profile->DrawCopy();
1777    SavePadToPDF(canvasProfiles);
1778    canvasProfiles->Write();
1779    SavePadToPDF(canvasIn);
1780    canvasIn->Write();
1781    if(canvasOut) {
1782        SavePadToPDF(canvasOut);
1783        canvasOut->Write();
1784    }
1785    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
1786    canvasRatioMeasuredRefoldedIn->Write();
1787    if(canvasRatioMeasuredRefoldedOut) {
1788        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
1789        canvasRatioMeasuredRefoldedOut->Write();
1790    }
1791    SavePadToPDF(canvasSpectraIn);
1792    canvasSpectraIn->Write();
1793    if(canvasSpectraOut) {
1794        SavePadToPDF(canvasSpectraOut);
1795        canvasSpectraOut->Write();
1796        SavePadToPDF(canvasRatio);
1797        canvasRatio->Write();
1798        SavePadToPDF(canvasV2);
1799        canvasV2->Write();
1800    }
1801    SavePadToPDF(canvasMasterIn);
1802    canvasMasterIn->Write();
1803    if(canvasMasterOut) {
1804        SavePadToPDF(canvasMasterOut);
1805        canvasMasterOut->Write();
1806    }
1807    SavePadToPDF(canvasMISC);
1808    canvasMISC->Write();
1809    // save the nomial canvasses
1810    SavePadToPDF(canvasNominalIn);
1811    canvasNominalIn->Write();
1812    if(canvasNominalOut) {
1813        SavePadToPDF(canvasNominalOut);
1814        canvasNominalOut->Write();
1815    }
1816    SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
1817    canvasNominalRatioMeasuredRefoldedIn->Write();
1818    if(canvasNominalRatioMeasuredRefoldedOut) {
1819        SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
1820        canvasNominalRatioMeasuredRefoldedOut->Write();
1821    }
1822    SavePadToPDF(canvasNominalSpectraIn);
1823    canvasNominalSpectraIn->Write();
1824    if(canvasNominalSpectraOut) {
1825        SavePadToPDF(canvasNominalSpectraOut);
1826        canvasNominalSpectraOut->Write();
1827        SavePadToPDF(canvasNominalRatio);
1828        canvasNominalRatio->Write();
1829        SavePadToPDF(canvasNominalV2);
1830        canvasNominalV2->Write();
1831    }
1832    SavePadToPDF(canvasNominalMasterIn);
1833    canvasNominalMasterIn->Write();
1834    if(canvasNominalMasterOut) {
1835        SavePadToPDF(canvasNominalMasterOut);
1836        canvasNominalMasterOut->Write();
1837    }
1838    SavePadToPDF(canvasNominalMISC);
1839    canvasNominalMISC->Write();
1840
1841    output.Write();
1842    output.Close();
1843 }
1844 //_____________________________________________________________________________
1845 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
1846 {
1847    // go through the output file and perform post processing routines
1848    // can either be performed in one go with the unfolding, or at a later stage
1849    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
1850    TFile readMe(in.Data(), "READ");     // open file read-only
1851    if(readMe.IsZombie()) {
1852        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1853        return;
1854    }
1855    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
1856    readMe.ls();
1857    TList* listOfKeys((TList*)readMe.GetListOfKeys());
1858    if(!listOfKeys) {
1859        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1860        return;
1861    }
1862    // prepare necessary canvasses
1863    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
1864    TCanvas* canvasOut(0x0);
1865    if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
1866    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
1867    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
1868    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
1869    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
1870    TCanvas* canvasSpectraOut(0x0);
1871    if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
1872    TCanvas* canvasRatio(0x0);
1873    if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
1874    TCanvas* canvasV2(0x0);
1875    if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
1876    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1877    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
1878    TCanvas* canvasMasterOut(0x0);
1879    if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
1880    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
1881    TDirectoryFile* defDir(0x0);
1882    TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
1883    canvasProfiles->Divide(2);
1884    TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1885    TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1886    // get an estimate of the number of outputs and find the default set
1887    Int_t cacheMe(0);
1888    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1889        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1890            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1891            cacheMe++;
1892        }
1893    }
1894    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
1895    canvasIn->Divide(columns, rows);
1896    if(canvasOut) canvasOut->Divide(columns, rows);
1897    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1898    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1899    canvasSpectraIn->Divide(columns, rows);
1900    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
1901    if(canvasRatio) canvasRatio->Divide(columns, rows);
1902    if(canvasV2) canvasV2->Divide(columns, rows);
1903
1904    canvasMasterIn->Divide(columns, rows);
1905    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
1906    // extract the default output
1907    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
1908    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
1909    if(defDir) {
1910        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1911        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
1912        if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1913        if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
1914        printf(" > succesfully extracted default results < \n");
1915    }
1916
1917    // loop through the directories, only plot the graphs if the deconvolution converged
1918    TDirectoryFile* tempDir(0x0);
1919    TDirectoryFile* tempIn(0x0);
1920    TDirectoryFile*  tempOut(0x0);
1921    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
1922        // read the directory by its name
1923        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1924        if(!tempDir) continue;
1925        TString dirName(tempDir->GetName());
1926        // try to read the in- and out of plane subdirs
1927        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1928        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1929        j++;
1930        if(!tempIn) {    // attempt to get MakeAU output
1931            TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
1932            TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
1933            tempPearson->Divide(4,2);
1934            TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
1935            tempRatio->Divide(4,2);
1936            TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
1937            tempResult->Divide(4,2);
1938            for(Int_t q(0); q < 8; q++) {
1939                tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
1940                if(tempIn) {
1941                     // to see if the unfolding converged try to extract the pearson coefficients
1942                     TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1943                     if(pIn) {
1944                        printf(" - %s in plane converged \n", dirName.Data());
1945                         tempPearson->cd(1+q);
1946                         Style(gPad, "PEARSON");
1947                         pIn->DrawCopy("colz");
1948                         TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1949                         if(rIn) {
1950                             Style(rIn);
1951                             printf(" > found RatioRefoldedMeasured < \n");
1952                             tempRatio->cd(q+1);
1953                             rIn->SetFillColor(kRed);
1954                             rIn->Draw("ap");
1955                         }
1956                         TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1957                         TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1958                         TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1959                         TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1960                         if(dvector && avalue && rm && eff) {
1961                             Style(dvector);
1962                             Style(avalue);
1963                             Style(rm);
1964                             Style(eff);
1965                             canvasMISC->cd(1);
1966                             Style(gPad, "SPECTRUM");
1967                             dvector->DrawCopy();
1968                             canvasMISC->cd(2);
1969                             Style(gPad, "SPECTRUM");
1970                             avalue->DrawCopy();
1971                             canvasMISC->cd(3);
1972                             Style(gPad, "PEARSON");
1973                             rm->DrawCopy("colz");
1974                             canvasMISC->cd(4);
1975                             eff->DrawCopy();
1976                         } else if(rm && eff) {
1977                             Style(rm);
1978                             Style(eff);
1979                             canvasMISC->cd(3);
1980                             Style(gPad, "PEARSON");
1981                             rm->DrawCopy("colz");
1982                             canvasMISC->cd(4);
1983                             eff->DrawCopy();
1984                         }
1985                     }
1986                     TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1987                     TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1988                     TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1989                     if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1990                         if(defaultUnfoldedJetSpectrumIn) {
1991                             Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
1992                             TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
1993                             temp->Divide(unfoldedSpectrum);
1994                             temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1995                             temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1996                             temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1997                             canvasMasterIn->cd(j);
1998                             temp->GetYaxis()->SetRangeUser(0., 2);
1999                             temp->DrawCopy();
2000                         }
2001                         TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2002                         tempResult->cd(q+1);
2003                         Style(gPad);
2004                         Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2005                         unfoldedSpectrum->DrawCopy();
2006                         Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2007                         inputSpectrum->DrawCopy("same");
2008                         Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2009                         refoldedSpectrum->DrawCopy("same");
2010                         TLegend* l(AddLegend(gPad));
2011                         Style(l);
2012                         if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2013                             Float_t chi(fitStatus->GetBinContent(1));
2014                             Float_t pen(fitStatus->GetBinContent(2));
2015                             Int_t dof((int)fitStatus->GetBinContent(3));
2016                             Float_t beta(fitStatus->GetBinContent(4));
2017                             l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2018                         } else if (fitStatus) { // only available in SVD fit
2019                             Int_t reg((int)fitStatus->GetBinContent(1));
2020                             l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2021                         }
2022                     }
2023                }
2024            }
2025        }
2026        if(tempIn) {
2027            // to see if the unfolding converged try to extract the pearson coefficients
2028            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2029            if(pIn) {
2030                printf(" - %s in plane converged \n", dirName.Data());
2031                canvasIn->cd(j);
2032                Style(gPad, "PEARSON");
2033                pIn->DrawCopy("colz");
2034                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2035                if(rIn) {
2036                    Style(rIn);
2037                    printf(" > found RatioRefoldedMeasured < \n");
2038                    canvasRatioMeasuredRefoldedIn->cd(j);
2039                    rIn->SetFillColor(kRed);
2040                    rIn->Draw("ap");
2041                }
2042                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2043                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2044                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2045                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2046                if(dvector && avalue && rm && eff) {
2047                    Style(dvector);
2048                    Style(avalue);
2049                    Style(rm);
2050                    Style(eff);
2051                    canvasMISC->cd(1);
2052                    Style(gPad, "SPECTRUM");
2053                    dvector->DrawCopy();
2054                    canvasMISC->cd(2);
2055                    Style(gPad, "SPECTRUM");
2056                    avalue->DrawCopy();
2057                    canvasMISC->cd(3);
2058                    Style(gPad, "PEARSON");
2059                    rm->DrawCopy("colz");
2060                    canvasMISC->cd(4);
2061                    eff->DrawCopy();
2062                } else if(rm && eff) {
2063                    Style(rm);
2064                    Style(eff);
2065                    canvasMISC->cd(3);
2066                    Style(gPad, "PEARSON");
2067                    rm->DrawCopy("colz");
2068                    canvasMISC->cd(4);
2069                    eff->DrawCopy();
2070                }
2071            }
2072            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2073            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2074            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2075            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2076                if(defaultUnfoldedJetSpectrumIn) {
2077                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2078                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
2079                    temp->Divide(unfoldedSpectrum);
2080                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
2081                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2082                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2083                    canvasMasterIn->cd(j);
2084                    temp->GetYaxis()->SetRangeUser(0., 2);
2085                    temp->DrawCopy();
2086                }
2087                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2088                canvasSpectraIn->cd(j);
2089                Style(gPad);
2090                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2091                unfoldedSpectrum->DrawCopy();
2092                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2093                inputSpectrum->DrawCopy("same");
2094                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2095                refoldedSpectrum->DrawCopy("same");
2096                TLegend* l(AddLegend(gPad));
2097                Style(l);
2098                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2099                    Float_t chi(fitStatus->GetBinContent(1));
2100                    Float_t pen(fitStatus->GetBinContent(2));
2101                    Int_t dof((int)fitStatus->GetBinContent(3));
2102                    Float_t beta(fitStatus->GetBinContent(4));
2103                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2104                } else if (fitStatus) { // only available in SVD fit
2105                    Int_t reg((int)fitStatus->GetBinContent(1));
2106                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2107                }
2108            }
2109        }
2110        if(tempOut) {
2111            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
2112            if(pOut) {
2113                printf(" - %s out of plane converged \n", dirName.Data());
2114                canvasOut->cd(j);
2115                Style(gPad, "PEARSON");
2116                pOut->DrawCopy("colz");
2117                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2118                if(rOut) {
2119                    Style(rOut);
2120                    printf(" > found RatioRefoldedMeasured < \n");
2121                    canvasRatioMeasuredRefoldedOut->cd(j);
2122                    rOut->SetFillColor(kRed);
2123                    rOut->Draw("ap");
2124                }
2125                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2126                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2127                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
2128                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
2129                if(dvector && avalue && rm && eff) {
2130                    Style(dvector);
2131                    Style(avalue);
2132                    Style(rm);
2133                    Style(eff);
2134                    canvasMISC->cd(5);
2135                    Style(gPad, "SPECTRUM");
2136                    dvector->DrawCopy();
2137                    canvasMISC->cd(6);
2138                    Style(gPad, "SPECTRUM");
2139                    avalue->DrawCopy();
2140                    canvasMISC->cd(7);
2141                    Style(gPad, "PEARSON");
2142                    rm->DrawCopy("colz");
2143                    canvasMISC->cd(8);
2144                    eff->DrawCopy();
2145                } else if(rm && eff) {
2146                    Style(rm);
2147                    Style(eff);
2148                    canvasMISC->cd(7);
2149                    Style(gPad, "PEARSON");
2150                    rm->DrawCopy("colz");
2151                    canvasMISC->cd(8);
2152                    eff->DrawCopy();
2153                }
2154            }
2155            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
2156            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
2157            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
2158            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2159                if(defaultUnfoldedJetSpectrumOut) {
2160                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2161                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
2162                    temp->Divide(unfoldedSpectrum);
2163                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
2164                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2165                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2166                    canvasMasterOut->cd(j);
2167                    temp->GetYaxis()->SetRangeUser(0., 2.);
2168                    temp->DrawCopy();
2169                }
2170                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
2171                canvasSpectraOut->cd(j);
2172                Style(gPad);
2173                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2174                unfoldedSpectrum->DrawCopy();
2175                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2176                inputSpectrum->DrawCopy("same");
2177                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2178                refoldedSpectrum->DrawCopy("same");
2179                TLegend* l(AddLegend(gPad));
2180                Style(l);
2181                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2182                    Float_t chi(fitStatus->GetBinContent(1));
2183                    Float_t pen(fitStatus->GetBinContent(2));
2184                    Int_t dof((int)fitStatus->GetBinContent(3));
2185                    Float_t beta(fitStatus->GetBinContent(4));
2186                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2187                } else if (fitStatus) { // only available in SVD fit
2188                    Int_t reg((int)fitStatus->GetBinContent(1));
2189                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2190                }
2191            }
2192        }
2193        if(canvasRatio && canvasV2) {
2194            canvasRatio->cd(j);
2195            TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
2196            Double_t _tempx(0), _tempy(0);
2197            if(ratioYield) {
2198                Style(ratioYield);
2199                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2200                    ratioYield->GetPoint(b+1, _tempx, _tempy);
2201                    ratioProfile->Fill(_tempx, _tempy);
2202                }
2203                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2204                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2205                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2206                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2207                ratioYield->SetFillColor(kRed);
2208                ratioYield->Draw("ap");
2209            }
2210            canvasV2->cd(j);
2211            TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
2212            if(ratioV2) {
2213                Style(ratioV2);
2214                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2215                    ratioV2->GetPoint(b+1, _tempx, _tempy);
2216                    v2Profile->Fill(_tempx, _tempy);
2217
2218                }
2219                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2220                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2221                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2222                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2223                ratioV2->SetFillColor(kRed);
2224                ratioV2->Draw("ap");
2225            }
2226        }
2227    }
2228    TFile output(out.Data(), "RECREATE");
2229    canvasProfiles->cd(1);
2230    Style(ratioProfile);
2231    ratioProfile->DrawCopy();
2232    canvasProfiles->cd(2);
2233    Style(v2Profile);
2234    v2Profile->DrawCopy();
2235    SavePadToPDF(canvasProfiles);
2236    canvasProfiles->Write();
2237    SavePadToPDF(canvasIn);
2238    canvasIn->Write();
2239    if(canvasOut) {
2240        SavePadToPDF(canvasOut);
2241        canvasOut->Write();
2242    }
2243    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2244    canvasRatioMeasuredRefoldedIn->Write();
2245    if(canvasRatioMeasuredRefoldedOut) {
2246        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2247        canvasRatioMeasuredRefoldedOut->Write();
2248    }
2249    SavePadToPDF(canvasSpectraIn);
2250    canvasSpectraIn->Write();
2251    if(canvasSpectraOut) {
2252        SavePadToPDF(canvasSpectraOut);
2253        canvasSpectraOut->Write();
2254        SavePadToPDF(canvasRatio);
2255        canvasRatio->Write();
2256        SavePadToPDF(canvasV2);
2257        canvasV2->Write();
2258    }
2259    SavePadToPDF(canvasMasterIn);
2260    canvasMasterIn->Write();
2261    if(canvasMasterOut) {
2262        SavePadToPDF(canvasMasterOut);
2263        canvasMasterOut->Write();
2264    }
2265    SavePadToPDF(canvasMISC);
2266    canvasMISC->Write();
2267    output.Write();
2268    output.Close();
2269 }
2270 //_____________________________________________________________________________
2271 Bool_t AliJetFlowTools::SetRawInput (
2272         TH2D* detectorResponse,  // detector response matrix
2273         TH1D* jetPtIn,           // in plane jet spectrum
2274         TH1D* jetPtOut,          // out of plane jet spectrum
2275         TH1D* dptIn,             // in plane delta pt distribution
2276         TH1D* dptOut,            // out of plane delta pt distribution
2277         Int_t eventCount) {
2278     // set input histograms manually
2279     fDetectorResponse   = detectorResponse;
2280     fSpectrumIn         = jetPtIn;
2281     fSpectrumOut        = jetPtOut;
2282     fDptInDist          = dptIn;
2283     fDptOutDist         = dptOut;
2284     fRawInputProvided   = kTRUE;
2285     // check if all data is provided
2286     if(!fDetectorResponse) {
2287         printf(" fDetectorResponse not found \n ");
2288         return kFALSE;
2289     }
2290     // check if the pt bin for true and rec have been set
2291     if(!fBinsTrue || !fBinsRec) {
2292         printf(" No true or rec bins set, please set binning ! \n");
2293         return kFALSE;
2294     }
2295     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
2296                           // procedures, these profiles will be nonsensical, user is responsible
2297         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2298         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2299         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2300     }
2301     // normalize spectra to event count if requested
2302     if(fNormalizeSpectra) {
2303         fEventCount = eventCount;
2304         if(fEventCount > 0) {
2305             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
2306             fSpectrumOut->Sumw2();
2307             fSpectrumIn->Scale(1./((double)fEventCount));
2308             fSpectrumOut->Scale(1./((double)fEventCount));
2309         }
2310     }
2311     if(!fNormalizeSpectra && fEventCount > 0) {
2312         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
2313         fSpectrumOut->Sumw2();
2314         fSpectrumIn->Scale(1./((double)fEventCount));
2315         fSpectrumOut->Scale(1./((double)fEventCount));
2316     }
2317     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
2318     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
2319     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
2320     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
2321     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
2322     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
2323     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
2324     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
2325
2326     return kTRUE;
2327 }
2328 //_____________________________________________________________________________
2329 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
2330 {
2331     // return ratio of h1 / h2
2332     // histograms can have different binning. errors are propagated as uncorrelated
2333     if(!(h1 && h2) ) {
2334         printf(" GetRatio called with NULL argument(s) \n ");
2335         return 0x0;
2336     }
2337     Int_t j(0);
2338     TGraphErrors *gr = new TGraphErrors();
2339     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2340     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
2341     TH1* dud((TH1*)h1->Clone("dud"));
2342     dud->Sumw2();
2343     h1->Sumw2();
2344     h2->Sumw2();
2345     if(!dud->Divide(h1, h2)) {
2346         printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
2347         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2348             binCent = h1->GetXaxis()->GetBinCenter(i);
2349             if(xmax > 0. && binCent > xmax) continue;
2350             j = h2->FindBin(binCent);
2351             binWidth = h1->GetXaxis()->GetBinWidth(i);
2352             if(h2->GetBinContent(j) > 0.) {
2353                 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
2354             /* original propagation of uncertainty changed 08012014
2355             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
2356             Double_t B = 0.;
2357             if(h2->GetBinError(j)>0.) {
2358                 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
2359                 error2 = A*A + B*B;
2360             } else error2 = A*A;        */
2361                 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
2362                 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
2363                 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
2364                 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
2365                 gr->SetPoint(gr->GetN(), binCent, ratio);
2366                 gr->SetPointError(gr->GetN()-1, 0.5*binWidth, error2);
2367             }
2368         }
2369     } else {
2370         printf( " > using ROOT to divide two histograms < \n");
2371         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2372             binCent = dud->GetXaxis()->GetBinCenter(i);
2373             if(xmax > 0. && binCent > xmax) continue;
2374             binWidth = dud->GetXaxis()->GetBinWidth(i);
2375             gr->SetPoint(gr->GetN(),binCent,dud->GetBinContent(i));
2376             gr->SetPointError(gr->GetN()-1,0.5*binWidth,dud->GetBinError(i));
2377         }
2378     }
2379
2380     if(appendFit) {
2381         TF1* fit(new TF1("lin", "pol0", 10, 100));
2382         gr->Fit(fit);
2383     }
2384     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
2385     if(dud) delete dud;
2386     return gr;
2387 }
2388 //_____________________________________________________________________________
2389 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
2390 {
2391     // get v2 from difference of in plane, out of plane yield
2392     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
2393     // different binning is allowed but will mean that the error
2394     // propagation is unreliable
2395     // r is the event plane resolution for the chosen centrality
2396     if(!(h1 && h2) ) {
2397         printf(" GetV2 called with NULL argument(s) \n ");
2398         return 0x0;
2399     }
2400     Int_t j(0);
2401     TGraphErrors *gr = new TGraphErrors();
2402     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2403     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
2404     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
2405     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2406         binCent = h1->GetXaxis()->GetBinCenter(i);
2407         j = h2->FindBin(binCent);
2408         binWidth = h1->GetXaxis()->GetBinWidth(i);
2409         if(h2->GetBinContent(j) > 0.) {
2410             in = h1->GetBinContent(i);
2411             ein = h1->GetBinError(i);
2412             out = h2->GetBinContent(j);
2413             eout = h2->GetBinError(j);
2414             ratio = pre*((in-out)/(in+out));
2415             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);
2416             if(error2 > 0) error2 = TMath::Sqrt(error2);
2417             gr->SetPoint(gr->GetN(),binCent,ratio);
2418             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
2419         }
2420     }
2421     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
2422     return gr;
2423 }
2424 //_____________________________________________________________________________
2425 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
2426     // write object, if a unique identifier is given the object is cloned
2427     // and the clone is saved. setting kill to true will delete the original obect from the heap
2428     if(!object) {
2429         printf(" > WriteObject:: called with NULL arguments \n ");
2430         return;
2431     } else if(!strcmp("", suffix.Data())) object->Write();
2432     else {
2433         TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
2434         newObject->Write();
2435     }
2436     if(kill) delete object;
2437 }
2438 //_____________________________________________________________________________
2439 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
2440     // construt a delta pt response matrix from supplied dpt distribution
2441     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
2442     // do a weighted rebinning to a (coarser) dpt distribution
2443     // be careful with the binning of the dpt response: it should be equal to that
2444     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
2445     //
2446     // the response matrix will be square and have the same binning
2447     // (min, max and granularity) of the input histogram
2448     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
2449     Double_t _bins[bins+1];             // prepare array with bin borders
2450     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
2451     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
2452     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
2453     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
2454         Bool_t skip = kFALSE;
2455         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
2456             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
2457             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
2458         }
2459     }
2460     return res;
2461 }
2462 //_____________________________________________________________________________
2463 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
2464     if(!binsTrue || !binsRec) {
2465         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
2466         return 0x0;
2467     }
2468     TString name(Form("unityResponse_%s", suffix.Data()));
2469     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
2470     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
2471         for(Int_t j(0); j < binsRec->GetSize(); j++) {
2472             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
2473         }
2474     }
2475     return unity;
2476 }
2477 //_____________________________________________________________________________
2478 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
2479     // save configuration parameters to histogram
2480     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
2481     summary->SetBinContent(1, fBetaIn);
2482     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
2483     summary->SetBinContent(2, fBetaOut);
2484     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
2485     summary->SetBinContent(3, fCentralityBin);
2486     summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
2487     summary->SetBinContent(4, (int)convergedIn);
2488     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
2489     summary->SetBinContent(5, (int)convergedOut);
2490     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
2491     summary->SetBinContent(6, (int)fAvoidRoundingError);
2492     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
2493     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
2494     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
2495     summary->SetBinContent(8, (int)fPrior);
2496     summary->GetXaxis()->SetBinLabel(8, "fPrior");
2497     summary->SetBinContent(9, fSVDRegIn);
2498     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
2499     summary->SetBinContent(10, fSVDRegOut);
2500     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
2501     summary->SetBinContent(11, (int)fSVDToy);
2502     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
2503     summary->SetBinContent(12, fJetRadius);
2504     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
2505     summary->SetBinContent(13, (int)fNormalizeSpectra);
2506     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
2507     summary->SetBinContent(14, (int)fSmoothenPrior);
2508     summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
2509     summary->SetBinContent(15, (int)fTestMode);
2510     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
2511     summary->SetBinContent(16, (int)fUseDetectorResponse);
2512     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
2513     summary->SetBinContent(17, fBayesianIterIn);
2514     summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
2515     summary->SetBinContent(18, fBayesianIterOut);
2516     summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
2517     summary->SetBinContent(19, fBayesianSmoothIn);
2518     summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
2519     summary->SetBinContent(20, fBayesianSmoothOut);
2520     summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
2521 }
2522 //_____________________________________________________________________________
2523 void AliJetFlowTools::ResetAliUnfolding() {
2524      // ugly function: reset all unfolding parameters
2525      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
2526      if(fitter) {
2527          printf(" > Found fitter, will delete it < \n");
2528          delete fitter;
2529      }
2530      if(gMinuit) {
2531          printf(" > Found gMinuit, will re-create it < \n");
2532          delete gMinuit;
2533          gMinuit = new TMinuit;
2534      }
2535      AliUnfolding::fgCorrelationMatrix = 0;
2536      AliUnfolding::fgCorrelationMatrixSquared = 0;
2537      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
2538      AliUnfolding::fgCurrentESDVector = 0;
2539      AliUnfolding::fgEntropyAPriori = 0;
2540      AliUnfolding::fgEfficiency = 0;
2541      AliUnfolding::fgUnfoldedAxis = 0;
2542      AliUnfolding::fgMeasuredAxis = 0;
2543      AliUnfolding::fgFitFunction = 0;
2544      AliUnfolding::fgMaxInput  = -1;
2545      AliUnfolding::fgMaxParams = -1;
2546      AliUnfolding::fgOverflowBinLimit = -1;
2547      AliUnfolding::fgRegularizationWeight = 10000;
2548      AliUnfolding::fgSkipBinsBegin = 0;
2549      AliUnfolding::fgMinuitStepSize = 0.1;
2550      AliUnfolding::fgMinuitPrecision = 1e-6;
2551      AliUnfolding::fgMinuitMaxIterations = 1000000;
2552      AliUnfolding::fgMinuitStrategy = 1.;
2553      AliUnfolding::fgMinimumInitialValue = kFALSE;
2554      AliUnfolding::fgMinimumInitialValueFix = -1;
2555      AliUnfolding::fgNormalizeInput = kFALSE;
2556      AliUnfolding::fgNotFoundEvents = 0;
2557      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
2558      AliUnfolding::fgBayesianSmoothing  = 1;
2559      AliUnfolding::fgBayesianIterations = 10;
2560      AliUnfolding::fgDebug = kFALSE;
2561      AliUnfolding::fgCallCount = 0;
2562      AliUnfolding::fgPowern = 5;
2563      AliUnfolding::fChi2FromFit = 0.;
2564      AliUnfolding::fPenaltyVal  = 0.;
2565      AliUnfolding::fAvgResidual = 0.;
2566      AliUnfolding::fgPrintChi2Details = 0;
2567      AliUnfolding::fgCanvas = 0;
2568      AliUnfolding::fghUnfolded = 0;
2569      AliUnfolding::fghCorrelation = 0;
2570      AliUnfolding::fghEfficiency = 0;
2571      AliUnfolding::fghMeasured = 0;
2572      AliUnfolding::SetMinuitStepSize(1.);
2573      AliUnfolding::SetMinuitPrecision(1e-6);
2574      AliUnfolding::SetMinuitMaxIterations(100000);
2575      AliUnfolding::SetMinuitStrategy(2.);
2576      AliUnfolding::SetDebug(1);
2577 }
2578 //_____________________________________________________________________________
2579 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
2580     // protect heap by adding unique qualifier to name
2581     if(!protect) return 0x0;
2582     TH1D* p = (TH1D*)protect->Clone();
2583     TString tempString(fActiveString);
2584     tempString+=suffix;
2585     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2586     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2587     if(kill) delete protect;
2588     return p;
2589 }
2590 //_____________________________________________________________________________
2591 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
2592     // protect heap by adding unique qualifier to name
2593     if(!protect) return 0x0;
2594     TH2D* p = (TH2D*)protect->Clone();
2595     TString tempString(fActiveString);
2596     tempString+=suffix;
2597     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2598     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2599     if(kill) delete protect;
2600     return p;
2601 }
2602 //_____________________________________________________________________________
2603 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
2604     // protect heap by adding unique qualifier to name
2605     if(!protect) return 0x0;
2606     TGraphErrors* p = (TGraphErrors*)protect->Clone();
2607     TString tempString(fActiveString);
2608     tempString+=suffix;
2609     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2610     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
2611     if(kill) delete protect;
2612     return p;
2613 }
2614 //_____________________________________________________________________________
2615 void AliJetFlowTools::MakeAU() {
2616     // === azimuthal unfolding ===
2617     //
2618     // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
2619     // in transverse momentum and azimuthal correlation space to extract v2 params
2620     // settings are equal to the ones used for 'Make()'
2621     //
2622     // basic steps that are followed:
2623     // 1) rebin the raw output of the jet task to the desired binnings
2624     // 2) calls the unfolding routine
2625     // 3) writes output to file
2626     // can be repeated multiple times with different configurations
2627
2628     Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
2629     Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
2630     TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
2631     TH1D* dPtdPhi[6];
2632     for(Int_t i(0); i < 6; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
2633
2634     for(Int_t i(0); i < 8; i++) {
2635         // 1) manipulation of input histograms
2636         // check if the input variables are present
2637         if(!PrepareForUnfolding(low[i], up[i])) return;
2638         // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
2639         //     parts of the spectrum can end up in over or underflow bins
2640         TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
2641
2642         // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
2643         // the template will be used as a prior for the chi2 unfolding
2644         TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
2645
2646         // get the full response matrix from the dpt and the detector response
2647         fDetectorResponse = NormalizeTH2D(fDetectorResponse);
2648         // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
2649         // so that unfolding should return the initial spectrum
2650         if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
2651         else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
2652         else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
2653         else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
2654         // normalize each slide of the response to one
2655         NormalizeTH2D(fFullResponseIn);
2656         // resize to desired binning scheme
2657         TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
2658         // get the kinematic efficiency
2659         TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
2660         kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
2661         // suppress the errors
2662         for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
2663         TH1D* jetFindingEfficiency(0x0);
2664         if(fJetFindingEff) {
2665             jetFindingEfficiency = ProtectHeap(fJetFindingEff);
2666             jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
2667             jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
2668         }
2669         // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
2670         TH1D* unfoldedJetSpectrumIn(0x0);
2671         fActiveDir->cd();                   // select active dir
2672         TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
2673         dirIn->cd();                        // select inplane subdir
2674         // select the unfolding method
2675         unfoldedJetSpectrumIn = UnfoldWrapper(
2676             measuredJetSpectrumIn,
2677             resizedResponseIn,
2678             kinematicEfficiencyIn,
2679             measuredJetSpectrumTrueBinsIn,
2680             TString("dPtdPhiUnfolding"),
2681             jetFindingEfficiency);
2682         if(i==5) {
2683             resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
2684             resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
2685             resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
2686             resizedResponseIn = ProtectHeap(resizedResponseIn);
2687             resizedResponseIn->Write();
2688             kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
2689             kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
2690             kinematicEfficiencyIn->Write();
2691             fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
2692             fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
2693             fDetectorResponse->Write();
2694             // optional histograms
2695             if(fSaveFull) {
2696                 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
2697                 fSpectrumIn->Write();
2698                 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
2699                 fDptInDist->Write();
2700                 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
2701                 fDptIn->Write();
2702                 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
2703                 fFullResponseIn->Write();
2704             }
2705         }
2706         fActiveDir->cd();
2707         fDeltaPtDeltaPhi->Write();
2708         fJetPtDeltaPhi->Write();
2709
2710         TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
2711         Double_t integralError(0);
2712         for(Int_t j(0); j < 6; j++) {
2713             // get the integrated
2714            Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
2715            dPtdPhi[j]->SetBinContent(i+1, integral);
2716            dPtdPhi[j]->SetBinError(i+1, integralError);
2717         }
2718         dud->Write();
2719         // save the current state of the unfolding object
2720         SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
2721     }
2722     TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
2723     TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2724     for(Int_t i(0); i < 6; i++) {
2725         dPtdPhi[i]->Fit(fourier, "VI");
2726         v2->SetBinContent(1+i, fourier->GetParameter(1));
2727         v2->SetBinError(1+i, fourier->GetParError(1));
2728         dPtdPhi[i]->Write();
2729     }
2730     v2->Write();
2731 }
2732 //_____________________________________________________________________________