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 "TGraphAsymmErrors.h"
  46 #include "TLine.h"
  47 #include "TCanvas.h"
  48 #include "TLegend.h"
  49 #include "TArrayD.h"
  50 #include "TList.h"
  51 #include "TMinuit.h"
  52 #include "TVirtualFitter.h"
  53 #include "TLegend.h"
  54 #include "TCanvas.h"
  55 #include "TStyle.h"
  56 #include "TLine.h"
  57 #include "TMath.h"
  58 #include "TVirtualFitter.h"
  59 #include "TFitResultPtr.h"
  60 // aliroot includes
  61 #include "AliUnfolding.h"
  62 #include "AliAnaChargedJetResponseMaker.h"
  63 // class includes
  64 #include "AliJetFlowTools.h"
  65 // roo unfold includes (make sure you have these available on your system)
  66 #include "RooUnfold.h"
  67 #include "RooUnfoldResponse.h"
  68 #include "RooUnfoldSvd.h"
  69 #include "RooUnfoldBayes.h"
  70 #include "TSVDUnfold.h"
  71
  72 using namespace std;
  73 //_____________________________________________________________________________
  74 AliJetFlowTools::AliJetFlowTools() :
  75     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
  76     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
  77     fSaveFull           (kTRUE),
  78     fActiveString       (""),
  79     fActiveDir          (0x0),
  80     fInputList          (0x0),
  81     fRefreshInput       (kTRUE),
  82     fOutputFileName     ("UnfoldedSpectra.root"),
  83     fOutputFile         (0x0),
  84     fCentralityArray    (0x0),
  85     fCentralityWeights  (0x0),
  86     fDetectorResponse   (0x0),
  87     fJetFindingEff      (0x0),
  88     fBetaIn             (.1),
  89     fBetaOut            (.1),
  90     fBayesianIterIn     (4),
  91     fBayesianIterOut    (4),
  92     fBayesianSmoothIn   (0.),
  93     fBayesianSmoothOut  (0.),
  94     fAvoidRoundingError (kFALSE),
  95     fUnfoldingAlgorithm (kChi2),
  96     fPrior              (kPriorMeasured),
  97     fPriorUser          (0x0),
  98     fBinsTrue           (0x0),
  99     fBinsRec            (0x0),
 100     fBinsTruePrior      (0x0),
 101     fBinsRecPrior       (0x0),
 102     fSVDRegIn           (5),
 103     fSVDRegOut          (5),
 104     fSVDToy             (kTRUE),
 105     fJetRadius          (0.3),
 106     fEventCount         (-1),
 107     fNormalizeSpectra   (kFALSE),
 108     fSmoothenPrior      (kFALSE),
 109     fFitMin             (60.),
 110     fFitMax             (300.),
 111     fFitStart           (75.),
 112     fSmoothenCounts     (kTRUE),
 113     fTestMode           (kFALSE),
 114     fRawInputProvided   (kFALSE),
 115     fEventPlaneRes      (.63),
 116     fUseDetectorResponse(kTRUE),
 117     fUseDptResponse     (kTRUE),
 118     fTrainPower         (kTRUE),
 119     fDphiUnfolding      (kTRUE),
 120     fDphiDptUnfolding   (kFALSE),
 121     fExLJDpt            (kTRUE),
 122     fSetTreatCorrErrAsUncorrErr(kFALSE),
 123     fTitleFontSize      (-999.),
 124     fRMSSpectrumIn      (0x0),
 125     fRMSSpectrumOut     (0x0),
 126     fRMSRatio           (0x0),
 127     fRMSV2              (0x0),
 128     fDeltaPtDeltaPhi    (0x0),
 129     fJetPtDeltaPhi      (0x0),
 130     fSpectrumIn         (0x0),
 131     fSpectrumOut        (0x0),
 132     fDptInDist          (0x0),
 133     fDptOutDist         (0x0),
 134     fDptIn              (0x0),
 135     fDptOut             (0x0),
 136     fFullResponseIn     (0x0),
 137     fFullResponseOut    (0x0) { // class constructor
 138         // create response maker weight function (tuned to PYTHIA spectrum)
 139         fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
 140         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 141 }
 142 //_____________________________________________________________________________
 143 void AliJetFlowTools::Make() {
 144     // core function of the class
 145     if(fDphiDptUnfolding) {
 146         // to extract the yield as function of Dphi, Dpt - experimental
 147         MakeAU();
 148         return;
 149     }
 150     // 1) rebin the raw output of the jet task to the desired binnings
 151     // 2) calls the unfolding routine
 152     // 3) writes output to file
 153     // can be repeated multiple times with different configurations
 154
 155     // 1) manipulation of input histograms
 156     // check if the input variables are present
 157     if(fRefreshInput) {
 158         if(!PrepareForUnfolding()) {
 159             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
 160             return;
 161         }
 162     }
 163     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
 164     //     parts of the spectrum can end up in over or underflow bins
 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"), fEventPlaneRes));
 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"), fEventPlaneRes));
 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;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
 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", fCentralityArray->At(0)));
 870     if(!fInputList->FindObject(spectrumName.Data())) {
 871         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
 872         return kFALSE;
 873     }
 874
 875     // get the first scaled spectrum
 876     fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
 877     // clone the spectrum on the heap. this is necessary since scale or add change the
 878     // contents of the original histogram
 879     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
 880     fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
 881     printf("Extracted %i wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
 882     // merge subsequent bins (if any)
 883     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 884         spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
 885         printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
 886         fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
 887     }
 888
 889     // in plane spectrum
 890     if(!fDphiUnfolding) {
 891         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
 892         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
 893     } else {
 894         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
 895         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
 896         fSpectrumIn = ProtectHeap(fSpectrumIn);
 897         // out of plane spectrum
 898         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
 899         fSpectrumOut = ProtectHeap(fSpectrumOut);
 900     }
 901     // normalize spectra to event count if requested
 902     if(fNormalizeSpectra) {
 903         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
 904         rho->Scale(fCentralityWeights->At(0));
 905         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 906            rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
 907         }
 908         if(!rho) return 0x0;
 909         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
 910         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
 911         if(fEventCount > 0) {
 912             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
 913             fSpectrumOut->Sumw2();
 914             Double_t pt(0);
 915             Double_t error(0); // lots of issues with the errors here ...
 916             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
 917                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
 918                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
 919                 fSpectrumIn->SetBinContent(1+i, pt);
 920                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
 921                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
 922                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
 923             }
 924             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
 925                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
 926                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
 927                 fSpectrumOut->SetBinContent(1+i, pt);
 928                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
 929                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
 930                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
 931             }
 932         }
 933         if(normalizeToFullSpectrum) fEventCount = -1;
 934     }
 935     // extract the delta pt matrices
 936     TString deltaptName("");
 937     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
 938     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
 939     if(!fDeltaPtDeltaPhi) {
 940         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
 941         printf(" > may be ok, depending no what you want to do < \n");
 942         fRefreshInput = kTRUE;
 943         return kTRUE;
 944     }
 945
 946     // clone the distribution on the heap and if requested merge with other centralities
 947     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
 948     fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
 949     printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
 950     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 951         deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
 952         printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
 953         fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
 954     }
 955
 956     // in plane delta pt distribution
 957     if(!fDphiUnfolding) {
 958         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
 959         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
 960     } else {
 961         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
 962         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
 963         // out of plane delta pt distribution
 964         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
 965         fDptInDist = ProtectHeap(fDptInDist);
 966         fDptOutDist = ProtectHeap(fDptOutDist);
 967         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
 968     }
 969
 970     // create a rec - true smeared response matrix
 971     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
 972     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 973         Bool_t skip = kFALSE;
 974         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 975             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
 976             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 977         }
 978     }
 979     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
 980     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 981         Bool_t skip = kFALSE;
 982         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 983             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
 984             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 985         }
 986     }
 987     fDptIn = new TH2D(*rfIn);
 988     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
 989     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 990     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 991     fDptIn = ProtectHeap(fDptIn);
 992     fDptOut = new TH2D(*rfOut);
 993     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
 994     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 995     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 996     fDptOut = ProtectHeap(fDptOut);
 997
 998     fRefreshInput = kTRUE;     // force cloning of the input
 999     return kTRUE;
1000 }
1001 //_____________________________________________________________________________
1002 Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
1003     // prepare for unfoldingUA - more robust method to extract input spectra from file
1004     // will replace PrepareForUnfolding eventually (09012014)
1005     if(!fInputList) {
1006         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1007         return kFALSE;
1008     }
1009     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1010     // check if the pt bin for true and rec have been set
1011     if(!fBinsTrue || !fBinsRec) {
1012         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1013         return kFALSE;
1014     }
1015     if(!fTrainPower) {
1016         // clear minuit state to avoid constraining the fit with the results of the previous iteration
1017         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1018     }
1019     // extract the spectra
1020     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
1021     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1022     if(!fJetPtDeltaPhi) {
1023         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1024         return kFALSE;
1025     }
1026     if(fCentralityArray) {
1027         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1028             spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
1029             fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1030         }
1031     }
1032     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1033     // in plane spectrum
1034     fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1035     // extract the delta pt matrices
1036     TString deltaptName("");
1037     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1038     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1039     if(!fDeltaPtDeltaPhi) {
1040         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1041         printf(" > may be ok, depending no what you want to do < \n");
1042         fRefreshInput = kTRUE;
1043         return kTRUE;
1044     }
1045     if(fCentralityArray) {
1046         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1047             deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1048             fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1049         }
1050     }
1051
1052     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1053     // in plane delta pt distribution
1054     fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1055     // create a rec - true smeared response matrix
1056     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1057     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1058         Bool_t skip = kFALSE;
1059         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1060             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1061             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1062         }
1063     }
1064     fDptIn = new TH2D(*rfIn);
1065     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1066     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1067     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1068     fDptIn = ProtectHeap(fDptIn);
1069
1070     return kTRUE;
1071 }
1072 //_____________________________________________________________________________
1073 TH1D* AliJetFlowTools::GetPrior(
1074         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
1075         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
1076         const TH1D* kinematicEfficiency,              // kinematic efficiency
1077         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
1078         const TString suffix,                         // suffix (in, out)
1079         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
1080 {
1081     // 1) get a prior for unfolding.
1082     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1083     TH1D* unfolded(0x0);
1084     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1085     dirOut->cd();
1086     switch (fPrior) {    // select the prior for unfolding
1087         case kPriorPythia : {
1088             if(!fPriorUser) {
1089                 printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1090                 return 0x0;
1091             }
1092             // rebin the given prior to the true spectrum (creates a new histo)
1093             return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1094         } break;
1095         case kPriorChi2 : {
1096             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
1097                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original binning
1098                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
1099                 TArrayD* tempArrayRec(fBinsRec);
1100                 fBinsRec = fBinsRecPrior;
1101                 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1102                 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1103                 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1104                 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1105                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1106                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1107                 unfolded= UnfoldSpectrumChi2(
1108                             measuredJetSpectrumChi2,
1109                             resizedResponseChi2,
1110                             kinematicEfficiencyChi2,
1111                             measuredJetSpectrumTrueBinsChi2,    // prior for chi2 unfolding (measured)
1112                             TString(Form("prior_%s", suffix.Data())));
1113                if(!unfolded) {
1114                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1115                     printf("               probably Chi2 unfolding did not converge < \n");
1116                     return 0x0;
1117                 }
1118                 fBinsTrue = tempArrayTrue;  // reset bins borders
1119                 fBinsRec = tempArrayRec;
1120                 // if the chi2 prior has a different binning, rebin to the true binning for the  unfolding
1121                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
1122             } else {
1123                 unfolded = UnfoldSpectrumChi2(
1124                             measuredJetSpectrum,
1125                             resizedResponse,
1126                             kinematicEfficiency,
1127                             measuredJetSpectrumTrueBins,        // prior for chi2 unfolding (measured)
1128                             TString(Form("prior_%s", suffix.Data())));
1129                 if(!unfolded) {
1130                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1131                     printf("               probably Chi2 unfolding did not converge < \n");
1132                     return 0x0;
1133                 }
1134             }
1135             break;
1136         }
1137         case kPriorMeasured : {
1138             unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
1139         }
1140         default : break;
1141     }
1142     // it can be important that the prior is smooth, this can be achieved by
1143     // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1144     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1145     TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1146     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1147     return priorLocal;
1148 }
1149 //_____________________________________________________________________________
1150 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1151     // resize the x-axis of a th1d
1152     if(!histo) {
1153         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1154         return NULL;
1155     }
1156     // see how many bins we need to copy
1157     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);
1158     // low is the bin number of the first new bin
1159     Int_t l = histo->GetXaxis()->FindBin(low);
1160     // set the values
1161     Double_t _x(0), _xx(0);
1162     for(Int_t i(0); i < up-low; i++) {
1163         _x = histo->GetBinContent(l+i);
1164         _xx=histo->GetBinError(l+i);
1165         resized->SetBinContent(i+1, _x);
1166         resized->SetBinError(i+1, _xx);
1167     }
1168     return resized;
1169 }
1170 //_____________________________________________________________________________
1171 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1172     // resize the y-axis of a th2d
1173     if(!histo) {
1174         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1175         return NULL;
1176     }
1177     // see how many bins we need to copy
1178     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());
1179     // assume only the y-axis has changed
1180     // low is the bin number of the first new bin
1181     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1182     // set the values
1183     Double_t _x(0), _xx(0);
1184     for(Int_t i(0); i < x->GetSize(); i++) {
1185         for(Int_t j(0); j < y->GetSize(); j++) {
1186             _x = histo->GetBinContent(i, low+j);
1187             _xx=histo->GetBinError(i, low+1+j);
1188             resized->SetBinContent(i, j, _x);
1189             resized->SetBinError(i, j, _xx);
1190         }
1191     }
1192     return resized;
1193 }
1194 //_____________________________________________________________________________
1195 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1196     // general method to normalize all vertical slices of a th2 to unity
1197     // i.e. get a probability matrix
1198     if(!histo) {
1199         printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1200         return NULL;
1201     }
1202     Int_t binsX = histo->GetXaxis()->GetNbins();
1203     Int_t binsY = histo->GetYaxis()->GetNbins();
1204
1205     // normalize all slices in x
1206     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1207         Double_t weight = 0;
1208         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1209             weight+=histo->GetBinContent(i+1, j+1);
1210         }       // now we know the total weight
1211         for(Int_t j(0); j < binsY; j++) {
1212             if (weight <= 0 ) continue;
1213             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1214             if(noError) histo->SetBinError(  1+i, j+1, 0.);
1215             else histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1216         }
1217     }
1218     return histo;
1219 }
1220 //_____________________________________________________________________________
1221 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1222     // return a TH1D with the supplied histogram rebinned to the supplied bins
1223     // the returned histogram is new, the original is deleted from the heap if kill is true
1224     if(!histo || !bins) {
1225         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1226         return NULL;
1227     }
1228     // create the output histo
1229     TString name = histo->GetName();
1230     name+="_template";
1231     name+=suffix;
1232     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1233     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1234         // loop over the bins of the old histo and fill the new one with its data
1235         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1236     }
1237     if(kill) delete histo;
1238     return rebinned;
1239 }
1240 //_____________________________________________________________________________
1241 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1242     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1243     if(!fResponseMaker || !binsTrue || !binsRec) {
1244         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1245         return 0x0;
1246     }
1247     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1248     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1249 }
1250 //_____________________________________________________________________________
1251 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1252 {
1253     // multiply two matrices
1254     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1255     TH2D* c = (TH2D*)a->Clone("c");
1256     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1257         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1258             Double_t val = 0;
1259             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1260                 Int_t y2 = x1;
1261                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1262             }
1263             c->SetBinContent(x2, y1, val);
1264             c->SetBinError(x2, y1, 0.);
1265         }
1266     }
1267     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1268     return c;
1269 }
1270 //_____________________________________________________________________________
1271 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1272 {
1273     // normalize a th1d to a certain scale
1274     histo->Sumw2();
1275     Double_t integral = histo->Integral()*scale;
1276     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1277     else if (scale != 1.) histo->Scale(1./scale, "width");
1278     else printf(" > Histogram integral < 0, cannot normalize \n");
1279     return histo;
1280 }
1281 //_____________________________________________________________________________
1282 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1283 {
1284     // Calculate pearson coefficients from covariance matrix
1285     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1286     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1287     Double_t pearson(0.);
1288     if(nrows==0 && ncols==0) return 0x0;
1289     for(Int_t row = 0; row < nrows; row++) {
1290         for(Int_t col = 0; col<ncols; col++) {
1291         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1292         (*pearsonCoefficients)(row,col) = pearson;
1293         }
1294     }
1295     return pearsonCoefficients;
1296 }
1297 //_____________________________________________________________________________
1298 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1299     // smoothen the spectrum using a user defined function
1300     // returns a clone of the original spectrum if fitting failed
1301     // if kill is kTRUE the input spectrum will be deleted from the heap
1302     // if 'count' is selected, bins are filled with integers (necessary if the
1303     // histogram is interpreted in a routine which accepts only counts)
1304     if(!spectrum || !function) return 0x0;
1305     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1306     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1307     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1308     TFitResultPtr r = temp->Fit(function, "", "", min, max);
1309     if((int)r == 0) {   // MINUIT status
1310         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1311             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1312                 (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));
1313                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1314             }
1315         }
1316     }
1317     if(kill) delete spectrum;
1318     return temp;
1319 }
1320 //_____________________________________________________________________________
1321 void AliJetFlowTools::Style()
1322 {
1323     // set global style for your current aliroot session
1324     if(!gStyle) return;
1325     gStyle->SetCanvasColor(-1);
1326     gStyle->SetPadColor(-1);
1327     gStyle->SetFrameFillColor(-1);
1328     gStyle->SetHistFillColor(-1);
1329     gStyle->SetTitleFillColor(-1);
1330     gStyle->SetFillColor(-1);
1331     gStyle->SetFillStyle(4000);
1332     gStyle->SetStatStyle(0);
1333     gStyle->SetTitleStyle(0);
1334     gStyle->SetCanvasBorderSize(0);
1335     gStyle->SetFrameBorderSize(0);
1336     gStyle->SetLegendBorderSize(0);
1337     gStyle->SetStatBorderSize(0);
1338     gStyle->SetTitleBorderSize(0);
1339 }
1340 //_____________________________________________________________________________
1341 void AliJetFlowTools::Style(TCanvas* c, TString style)
1342 {
1343     // set a default style for a canvas
1344     if(!strcmp(style.Data(), "PEARSON")) {
1345         printf(" > style PEARSON canvas < \n");
1346         gStyle->SetOptStat(0);
1347         c->SetGridx();
1348         c->SetGridy();
1349         c->SetTicks();
1350         return;
1351     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1352         printf(" > style SPECTRUM canvas < \n");
1353         gStyle->SetOptStat(0);
1354         c->SetLogy();
1355         c->SetGridx();
1356         c->SetGridy();
1357         c->SetTicks();
1358         return;
1359     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1360 }
1361 //_____________________________________________________________________________
1362 void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1363 {
1364     // set a default style for a canvas
1365     c->SetLeftMargin(.25);
1366     c->SetBottomMargin(.25);
1367     if(!strcmp(style.Data(), "PEARSON")) {
1368         printf(" > style PEARSON pad < \n");
1369         gStyle->SetOptStat(0);
1370         c->SetGridx();
1371         c->SetGridy();
1372         c->SetTicks();
1373         return;
1374     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1375         printf(" > style SPECTRUM pad < \n");
1376         gStyle->SetOptStat(0);
1377         c->SetLogy();
1378         c->SetGridx();
1379         c->SetGridy();
1380         c->SetTicks();
1381         return;
1382     } else if (!strcmp(style.Data(), "GRID")) {
1383         printf(" > style GRID pad < \n");
1384         gStyle->SetOptStat(0);
1385         c->SetGridx();
1386         c->SetGridy();
1387         c->SetTicks();
1388     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1389 }
1390 //_____________________________________________________________________________
1391 void AliJetFlowTools::Style(TLegend* l)
1392 {
1393     // set a default style for a legend
1394 //    l->SetTextSize(.06);
1395     l->SetFillColor(0);
1396 //    l->SetFillStyle(4050);
1397     l->SetBorderSize(0);
1398 }
1399 //_____________________________________________________________________________
1400 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1401 {
1402     // style a histo
1403     h->SetLineColor(col);
1404     h->SetMarkerColor(col);
1405     h->SetLineWidth(2.);
1406     h->SetMarkerSize(1.);
1407     h->SetTitle("");
1408     h->GetYaxis()->SetLabelSize(0.05);
1409     h->GetXaxis()->SetLabelSize(0.05);
1410     h->GetYaxis()->SetTitleOffset(1.5);
1411     h->GetXaxis()->SetTitleOffset(1.5);
1412     h->GetYaxis()->SetTitleSize(.05);
1413     h->GetXaxis()->SetTitleSize(.05);
1414     switch (type) {
1415         case kInPlaneSpectrum : {
1416             h->SetTitle("IN PLANE");
1417             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1418             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1419         } break;
1420         case kOutPlaneSpectrum : {
1421             h->SetTitle("OUT OF PLANE");
1422             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1423             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1424        } break;
1425        case kUnfoldedSpectrum : {
1426             h->SetTitle("UNFOLDED");
1427             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1428             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1429        } break;
1430        case kFoldedSpectrum : {
1431             h->SetTitle("FOLDED");
1432             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1433             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1434        } break;
1435        case kMeasuredSpectrum : {
1436             h->SetTitle("MEASURED");
1437             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1438             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1439        } break;
1440        case kBar : {
1441             h->SetFillColor(col);
1442             h->SetBarWidth(.6);
1443             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1444             h->SetBarOffset(0.2);
1445        }
1446        default : break;
1447     }
1448 }
1449 //_____________________________________________________________________________
1450 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1451 {
1452     // style a histo
1453     h->SetLineColor(col);
1454     h->SetMarkerColor(col);
1455     h->SetLineWidth(2.);
1456     h->SetMarkerSize(1.);
1457     h->SetTitle("");
1458     h->SetFillColor(kYellow);
1459     h->GetYaxis()->SetLabelSize(0.05);
1460     h->GetXaxis()->SetLabelSize(0.05);
1461     h->GetYaxis()->SetTitleOffset(1.6);
1462     h->GetXaxis()->SetTitleOffset(1.6);
1463     h->GetYaxis()->SetTitleSize(.05);
1464     h->GetXaxis()->SetTitleSize(.05);
1465     h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]");
1466     switch (type) {
1467         case kInPlaneSpectrum : {
1468             h->SetTitle("IN PLANE");
1469             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1470         } break;
1471         case kOutPlaneSpectrum : {
1472             h->SetTitle("OUT OF PLANE");
1473             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1474        } break;
1475        case kUnfoldedSpectrum : {
1476             h->SetTitle("UNFOLDED");
1477             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1478        } break;
1479        case kFoldedSpectrum : {
1480             h->SetTitle("FOLDED");
1481             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1482        } break;
1483        case kMeasuredSpectrum : {
1484             h->SetTitle("MEASURED");
1485             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1486        } break;
1487        case kRatio : {
1488             h->GetYaxis()->SetTitle("#frac{d#it{N_{in plane}^{jet}}}{d#it{p}_{T}} / #frac{d#it{N_{out of plane}^{jet}}}{d#it{p}_{T}}");
1489        } break;
1490        case kV2 : {
1491             h->GetYaxis()->SetTitle("#it{v}_{2} = #frac{1}{#it{R}} #frac{#pi}{4} #frac{#it{N_{in plane}} - #it{N_{out of plane}}}{#it{N_{in plane}} + #it{N_{out of plane}}}");
1492             h->GetYaxis()->SetRangeUser(-.5, 1.);
1493        } break;
1494        default : break;
1495     }
1496 }
1497 //_____________________________________________________________________________
1498 void AliJetFlowTools::GetNominalValues(
1499         TH1D*& ratio,           // pointer reference, output of this function
1500         TGraphErrors*& v2,      // pointer reference, as output of this function
1501         TArrayI* in,
1502         TArrayI* out,
1503         TString inFile,
1504         TString outFile) const
1505 {
1506     // pass clones of the nominal points and nominal v2 values
1507     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1508     TFile* readMe(new TFile(inFile.Data(), "READ"));                    // open unfolding output read-only
1509     if(readMe->IsZombie()) {
1510         printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1511         return;
1512     }
1513     printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1514     readMe->ls();
1515     TFile* output(new TFile(outFile.Data(), "RECREATE"));   // create new output file
1516     // get some placeholders, have to be initialized but will be deleted
1517     ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1518     TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1519     TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1520     Int_t iIn[] = {in->At(0), in->At(0)};          // first value is the nominal point
1521     Int_t iOut[] = {out->At(0), out->At(0)};
1522
1523     // call the functions and set the relevant pointer references
1524     TH1D* dud(0x0);
1525     DoIntermediateSystematics(
1526             new TArrayI(2, iIn),
1527             new TArrayI(2, iOut),
1528             dud, dud, dud, dud, dud, dud,
1529             ratio,              // pointer reference, output of this function
1530             nominalIn,
1531             nominalOut,
1532             1,
1533             fBinsTrue->At(0),
1534             fBinsTrue->At(fBinsTrue->GetSize()),
1535             readMe,
1536             "nominal_values");
1537     v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1538
1539     // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1540     ratio->SetDirectory(0);     // disassociate from current gDirectory
1541     readMe->Close();
1542 }
1543 //_____________________________________________________________________________
1544 void AliJetFlowTools::GetCorrelatedUncertainty(
1545         TGraphAsymmErrors*& corrRatio,  // correlated uncertainty function pointer
1546         TGraphAsymmErrors*& corrV2,     // correlated uncertainty function pointer
1547         TArrayI* variationsIn,          // variantions in plnae
1548         TArrayI* variationsOut,         // variantions out of plane
1549         TArrayI* variations2ndIn,      // second source of variations
1550         TArrayI* variations2ndOut,     // second source of variations
1551         TString type,                   // type of uncertaitny
1552         Int_t columns,                  // divide the output canvasses in this many columns
1553         Float_t rangeLow,               // lower pt range
1554         Float_t rangeUp,                // upper pt range
1555         TString in,                     // input file name (created by this unfolding class)
1556         TString out                     // output file name (which will hold results of the systematic test)
1557         ) const
1558 {
1559     // do full systematics
1560     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1561     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1562     if(readMe->IsZombie()) {
1563         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1564         return;
1565     }
1566     printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1567     readMe->ls();
1568     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1569
1570     // create some null placeholder pointers
1571     TH1D* relativeErrorVariationInLow(0x0);
1572     TH1D* relativeErrorVariationInUp(0x0);
1573     TH1D* relativeErrorVariationOutLow(0x0);
1574     TH1D* relativeErrorVariationOutUp(0x0);
1575     TH1D* relativeError2ndVariationInLow(0x0);
1576     TH1D* relativeError2ndVariationInUp(0x0);
1577     TH1D* relativeError2ndVariationOutLow(0x0);
1578     TH1D* relativeError2ndVariationOutUp(0x0);
1579     TH1D* relativeStatisticalErrorIn(0x0);
1580     TH1D* relativeStatisticalErrorOut(0x0);
1581     // histo for the nominal ratio in / out
1582     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1583     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1584     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1585
1586     // call the functions
1587     if(variationsIn && variationsOut) {
1588         DoIntermediateSystematics(
1589                 variationsIn,
1590                 variationsOut,
1591                 relativeErrorVariationInUp,        // pointer reference
1592                 relativeErrorVariationInLow,       // pointer reference
1593                 relativeErrorVariationOutUp,       // pointer reference
1594                 relativeErrorVariationOutLow,      // pointer reference
1595                 relativeStatisticalErrorIn,        // pointer reference
1596                 relativeStatisticalErrorOut,       // pointer reference
1597                 nominal,
1598                 nominalIn,
1599                 nominalOut,
1600                 columns,
1601                 rangeLow,
1602                 rangeUp,
1603                 readMe,
1604                 type);
1605         if(relativeErrorVariationInUp) {
1606             // canvas with the error from variation strength
1607             TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1608             relativeErrorVariation->Divide(2);
1609             relativeErrorVariation->cd(1);
1610             Style(gPad, "GRID");
1611             relativeErrorVariationInUp->DrawCopy("b");
1612             relativeErrorVariationInLow->DrawCopy("same b");
1613             Style(AddLegend(gPad));
1614             relativeErrorVariation->cd(2);
1615             Style(gPad, "GRID");
1616             relativeErrorVariationOutUp->DrawCopy("b");
1617             relativeErrorVariationOutLow->DrawCopy("same b");
1618             SavePadToPDF(relativeErrorVariation);
1619             Style(AddLegend(gPad));
1620             relativeErrorVariation->Write();
1621         }
1622
1623     }
1624     // call the functions for a second set of systematic sources
1625     if(variations2ndIn && variations2ndOut) {
1626         DoIntermediateSystematics(
1627                 variations2ndIn,
1628                 variations2ndOut,
1629                 relativeError2ndVariationInUp,        // pointer reference
1630                 relativeError2ndVariationInLow,       // pointer reference
1631                 relativeError2ndVariationOutUp,       // pointer reference
1632                 relativeError2ndVariationOutLow,      // pointer reference
1633                 relativeStatisticalErrorIn,           // pointer reference
1634                 relativeStatisticalErrorOut,          // pointer reference
1635                 nominal,
1636                 nominalIn,
1637                 nominalOut,
1638                 columns,
1639                 rangeLow,
1640                 rangeUp,
1641                 readMe,
1642                 type);
1643         if(relativeError2ndVariationInUp) {
1644             // canvas with the error from variation strength
1645             TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type.Data()), Form("relativeError2nd_%s", type.Data())));
1646             relativeError2ndVariation->Divide(2);
1647             relativeError2ndVariation->cd(1);
1648             Style(gPad, "GRID");
1649             relativeError2ndVariationInUp->DrawCopy("b");
1650             relativeError2ndVariationInLow->DrawCopy("same b");
1651             Style(AddLegend(gPad));
1652             relativeError2ndVariation->cd(2);
1653             Style(gPad, "GRID");
1654             relativeError2ndVariationOutUp->DrawCopy("b");
1655             relativeError2ndVariationOutLow->DrawCopy("same b");
1656             SavePadToPDF(relativeError2ndVariation);
1657             Style(AddLegend(gPad));
1658             relativeError2ndVariation->Write();
1659         }
1660
1661     }
1662
1663     // and the placeholder for the final systematic
1664     TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1665     TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1666     TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1667     TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1668     relativeErrorInUp->SetYTitle("relative uncertainty");
1669     relativeErrorOutUp->SetYTitle("relative uncertainty");
1670     relativeErrorInLow->SetYTitle("relative uncertainty");
1671     relativeErrorOutLow->SetYTitle("relative uncertainty");
1672
1673     // merge the correlated errors (FIXME) trivial for one set
1674     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1675     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1676     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1677     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1678
1679     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1680         // for the upper bound
1681         if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
1682         if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
1683         if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1684         if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
1685         dInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1686         if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1687         dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1688         if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1689         // for the lower bound
1690         if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
1691         if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
1692         if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1693         if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
1694         dInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1695         if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(dInLow));
1696         dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1697         if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1698     }
1699     // project the estimated errors on the nominal ratio
1700     if(nominal) {
1701         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1702         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1703         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1704         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1705         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1706         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1707         Double_t lowErr(0.), upErr(0.);
1708         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1709             // add the in and out of plane errors in quadrature
1710             if(fSetTreatCorrErrAsUncorrErr) {
1711                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1);
1712                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1713             } else {
1714                 // the in and out of plane correlated errors will be fully correlated, so take the correlation coefficient equal to 1
1715                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) - 1.*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1716                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 1.*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1717             }
1718             // set the errors
1719             ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
1720             ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
1721             // get the bin width (which is the 'error' on x
1722             Double_t binWidth(nominal->GetBinWidth(i+1));
1723             axl[i] = binWidth/2.;
1724             axh[i] = binWidth/2.;
1725             // now get the coordinate for the point
1726             ax[i] = nominal->GetBinCenter(i+1);
1727             ay[i] = nominal->GetBinContent(i+1);
1728         }
1729         // save the nominal ratio
1730         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
1731         corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
1732         nominalError->SetName("correlated uncertainty");
1733         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
1734         nominalCanvas->Divide(2);
1735         nominalCanvas->cd(1);
1736         Style(nominal, kBlack);
1737         Style(nominalError, kYellow, kRatio);
1738         nominalError->SetLineColor(kYellow);
1739         nominalError->SetMarkerColor(kYellow);
1740         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1741         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
1742         nominalError->DrawClone("a2");
1743         nominal->DrawCopy("same E1");
1744         Style(AddLegend(gPad));
1745         Style(gPad, "GRID");
1746         Style(nominalCanvas);
1747         // save nominal v2 and systematic errors
1748         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
1749                     nominalIn,
1750                     nominalOut,
1751                     fEventPlaneRes,
1752                     "v_{2} with correlated uncertainty",
1753                     relativeErrorInUp,
1754                     relativeErrorInLow,
1755                     relativeErrorOutUp,
1756                     relativeErrorOutLow,
1757                     .5));
1758         // pass the nominal values to the pointer references
1759         corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
1760         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
1761         nominalCanvas->cd(2);
1762         Style(nominalV2, kBlack);
1763         Style(nominalV2Error, kYellow, kV2);
1764         nominalV2Error->SetLineColor(kYellow);
1765         nominalV2Error->SetMarkerColor(kYellow);
1766         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1767         nominalV2Error->DrawClone("a2");
1768         nominalV2->DrawClone("same E1");
1769         Style(AddLegend(gPad));
1770         Style(gPad, "GRID");
1771         Style(nominalCanvas);
1772         SavePadToPDF(nominalCanvas);
1773         nominalCanvas->Write();
1774     }
1775
1776     TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
1777     relativeError->Divide(2);
1778     relativeError->cd(1);
1779     Style(gPad, "GRID");
1780     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1781     Style(relativeErrorInUp, kBlue, kBar);
1782     Style(relativeErrorInLow, kGreen, kBar);
1783     relativeErrorInUp->DrawCopy("b");
1784     relativeErrorInLow->DrawCopy("same b");
1785     Style(relativeStatisticalErrorIn, kRed);
1786     relativeStatisticalErrorIn->DrawCopy("same");
1787     Style(AddLegend(gPad));
1788     relativeError->cd(2);
1789     Style(gPad, "GRID");
1790     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1791     Style(relativeErrorOutUp, kBlue, kBar);
1792     Style(relativeErrorOutLow, kGreen, kBar);
1793     relativeErrorOutUp->DrawCopy("b");
1794     relativeErrorOutLow->DrawCopy("same b");
1795     Style(relativeStatisticalErrorOut, kRed);
1796     relativeStatisticalErrorOut->DrawCopy("same");
1797     Style(AddLegend(gPad));
1798
1799     // write the buffered file to disk and close the file
1800     SavePadToPDF(relativeError);
1801     relativeError->Write();
1802     output->Write();
1803 //    output->Close();
1804 }
1805 //_____________________________________________________________________________
1806 void AliJetFlowTools::GetShapeUncertainty(
1807         TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
1808         TGraphAsymmErrors*& shapeV2,    // pointer reference to final shape uncertainty of v2
1809         TArrayI* regularizationIn,      // regularization strength systematics, in plane
1810         TArrayI* regularizationOut,     // regularization strenght systematics, out of plane
1811         TArrayI* trueBinIn,             // true bin ranges, in plane
1812         TArrayI* trueBinOut,            // true bin ranges, out of plane
1813         TArrayI* recBinIn,              // rec bin ranges, in plane
1814         TArrayI* recBinOut,             // rec bin ranges, out of plane
1815         Int_t columns,                  // divide the output canvasses in this many columns
1816         Float_t rangeLow,               // lower pt range
1817         Float_t rangeUp,                // upper pt range
1818         TString in,                     // input file name (created by this unfolding class)
1819         TString out                     // output file name (which will hold results of the systematic test)
1820         ) const
1821 {
1822     // do full systematics
1823     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1824     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1825     if(readMe->IsZombie()) {
1826         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1827         return;
1828     }
1829     printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1830     readMe->ls();
1831     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1832
1833     // create some null placeholder pointers
1834     TH1D* relativeErrorRegularizationInLow(0x0);
1835     TH1D* relativeErrorRegularizationInUp(0x0);
1836     TH1D* relativeErrorTrueBinInLow(0x0);
1837     TH1D* relativeErrorTrueBinInUp(0x0);
1838     TH1D* relativeErrorRecBinInLow(0x0);
1839     TH1D* relativeErrorRecBinInUp(0x0);
1840     TH1D* relativeErrorRegularizationOutLow(0x0);
1841     TH1D* relativeErrorRegularizationOutUp(0x0);
1842     TH1D* relativeErrorTrueBinOutLow(0x0);
1843     TH1D* relativeErrorTrueBinOutUp(0x0);
1844     TH1D* relativeErrorRecBinOutLow(0x0);
1845     TH1D* relativeErrorRecBinOutUp(0x0);
1846     TH1D* relativeStatisticalErrorIn(0x0);
1847     TH1D* relativeStatisticalErrorOut(0x0);
1848     // histo for the nominal ratio in / out
1849     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1850     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1851     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1852
1853     // call the functions
1854     if(regularizationIn && regularizationOut) {
1855         DoIntermediateSystematics(
1856                 regularizationIn,
1857                 regularizationOut,
1858                 relativeErrorRegularizationInUp,        // pointer reference
1859                 relativeErrorRegularizationInLow,       // pointer reference
1860                 relativeErrorRegularizationOutUp,       // pointer reference
1861                 relativeErrorRegularizationOutLow,      // pointer reference
1862                 relativeStatisticalErrorIn,             // pointer reference
1863                 relativeStatisticalErrorOut,            // pointer reference
1864                 nominal,
1865                 nominalIn,
1866                 nominalOut,
1867                 columns,
1868                 rangeLow,
1869                 rangeUp,
1870                 readMe,
1871                 "regularization");
1872         if(relativeErrorRegularizationInUp) {
1873             // canvas with the error from regularization strength
1874             TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
1875             relativeErrorRegularization->Divide(2);
1876             relativeErrorRegularization->cd(1);
1877             Style(gPad, "GRID");
1878             relativeErrorRegularizationInUp->DrawCopy("b");
1879             relativeErrorRegularizationInLow->DrawCopy("same b");
1880             Style(AddLegend(gPad));
1881             relativeErrorRegularization->cd(2);
1882             Style(gPad, "GRID");
1883             relativeErrorRegularizationOutUp->DrawCopy("b");
1884             relativeErrorRegularizationOutLow->DrawCopy("same b");
1885             SavePadToPDF(relativeErrorRegularization);
1886             Style(AddLegend(gPad));
1887             relativeErrorRegularization->Write();
1888         }
1889     }
1890     if(trueBinIn && trueBinOut) {
1891         DoIntermediateSystematics(
1892                 trueBinIn,
1893                 trueBinOut,
1894                 relativeErrorTrueBinInUp,       // pointer reference
1895                 relativeErrorTrueBinInLow,      // pointer reference
1896                 relativeErrorTrueBinOutUp,      // pointer reference
1897                 relativeErrorTrueBinOutLow,     // pointer reference
1898                 relativeStatisticalErrorIn,
1899                 relativeStatisticalErrorOut,
1900                 nominal,
1901                 nominalIn,
1902                 nominalOut,
1903                 columns,
1904                 rangeLow,
1905                 rangeUp,
1906                 readMe,
1907                 "trueBin");
1908         if(relativeErrorTrueBinInUp) {
1909             TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
1910             relativeErrorTrueBin->Divide(2);
1911             relativeErrorTrueBin->cd(1);
1912             Style(gPad, "GRID");
1913             relativeErrorTrueBinInUp->DrawCopy("b");
1914             relativeErrorTrueBinInLow->DrawCopy("same b");
1915             Style(AddLegend(gPad));
1916             relativeErrorTrueBin->cd(2);
1917             Style(gPad, "GRID");
1918             relativeErrorTrueBinOutUp->DrawCopy("b");
1919             relativeErrorTrueBinOutLow->DrawCopy("same b");
1920             SavePadToPDF(relativeErrorTrueBin);
1921             Style(AddLegend(gPad));
1922             relativeErrorTrueBin->Write();
1923         }
1924     }
1925     if(recBinIn && recBinOut) {
1926         DoIntermediateSystematics(
1927                 recBinIn,
1928                 recBinOut,
1929                 relativeErrorRecBinInLow,       // pointer reference
1930                 relativeErrorRecBinInUp,        // pointer reference
1931                 relativeErrorRecBinOutLow,      // pointer reference
1932                 relativeErrorRecBinOutUp,       // pointer reference
1933                 relativeStatisticalErrorIn,
1934                 relativeStatisticalErrorOut,
1935                 nominal,
1936                 nominalIn,
1937                 nominalOut,
1938                 columns,
1939                 rangeLow,
1940                 rangeUp,
1941                 readMe,
1942                 "recBin");
1943         if(relativeErrorRecBinOutUp) {
1944             // canvas with the error from regularization strength
1945             TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
1946             relativeErrorRecBin->Divide(2);
1947             relativeErrorRecBin->cd(1);
1948             Style(gPad, "GRID");
1949             relativeErrorRecBinInUp->DrawCopy("b");
1950             relativeErrorRecBinInLow->DrawCopy("same b");
1951             Style(AddLegend(gPad));
1952             relativeErrorRecBin->cd(2);
1953             Style(gPad, "GRID");
1954             relativeErrorRecBinOutUp->DrawCopy("b");
1955             relativeErrorRecBinOutLow->DrawCopy("same b");
1956             SavePadToPDF(relativeErrorRecBin);
1957             Style(AddLegend(gPad));
1958             relativeErrorRecBin->Write();
1959         }
1960     }
1961
1962     // and the placeholder for the final systematic
1963     TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1964     TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1965     TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1966     TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1967     relativeErrorInUp->SetYTitle("relative uncertainty");
1968     relativeErrorOutUp->SetYTitle("relative uncertainty");
1969     relativeErrorInLow->SetYTitle("relative uncertainty");
1970     relativeErrorOutLow->SetYTitle("relative uncertainty");
1971
1972     // sum of squares for the total systematic error
1973     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1974     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1975     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1976     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1977
1978     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1979         // for the upper bound
1980         if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
1981         if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
1982         if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
1983         if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
1984         if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
1985         if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
1986         dInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1987         if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1988         dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1989         if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1990         // for the lower bound
1991         if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
1992         if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
1993         if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
1994         if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
1995         if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
1996         if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
1997         dInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1998         if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(dInLow));
1999         dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
2000         if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
2001     }
2002     // project the estimated errors on the nominal ratio
2003     if(nominal) {
2004         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2005         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2006         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2007         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2008         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2009         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2010         Double_t lowErr(0.), upErr(0.);
2011         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2012             // add the in and out of plane errors in quadrature
2013             lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1);
2014             upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
2015             // set the errors
2016             ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2017             ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2018             // get the bin width (which is the 'error' on x
2019             Double_t binWidth(nominal->GetBinWidth(i+1));
2020             axl[i] = binWidth/2.;
2021             axh[i] = binWidth/2.;
2022             // now get the coordinate for the point
2023             ax[i] = nominal->GetBinCenter(i+1);
2024             ay[i] = nominal->GetBinContent(i+1);
2025         }
2026         // save the nominal ratio
2027         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2028         shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2029         nominalError->SetName("shape uncertainty");
2030         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2031         nominalCanvas->Divide(2);
2032         nominalCanvas->cd(1);
2033         Style(nominal, kBlack);
2034         Style(nominalError, kYellow, kRatio);
2035         nominalError->SetLineColor(kYellow);
2036         nominalError->SetMarkerColor(kYellow);
2037         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2038         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2039         nominalError->DrawClone("a2");
2040         nominal->DrawCopy("same E1");
2041         Style(AddLegend(gPad));
2042         Style(gPad, "GRID");
2043         Style(nominalCanvas);
2044         // save nominal v2 and systematic errors
2045         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2046                     nominalIn,
2047                     nominalOut,
2048                     fEventPlaneRes,
2049                     "v_{2} with shape uncertainty",
2050                     relativeErrorInUp,
2051                     relativeErrorInLow,
2052                     relativeErrorOutUp,
2053                     relativeErrorOutLow));
2054         shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2055         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2056         nominalCanvas->cd(2);
2057         Style(nominalV2, kBlack);
2058         Style(nominalV2Error, kYellow, kV2);
2059         nominalV2Error->SetLineColor(kYellow);
2060         nominalV2Error->SetMarkerColor(kYellow);
2061         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2062         nominalV2Error->DrawClone("a2");
2063         nominalV2->DrawClone("same E1");
2064         Style(AddLegend(gPad));
2065         Style(gPad, "GRID");
2066         Style(nominalCanvas);
2067         SavePadToPDF(nominalCanvas);
2068         nominalCanvas->Write();
2069     }
2070
2071     TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2072     relativeError->Divide(2);
2073     relativeError->cd(1);
2074     Style(gPad, "GRID");
2075     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2076     Style(relativeErrorInUp, kBlue, kBar);
2077     Style(relativeErrorInLow, kGreen, kBar);
2078     relativeErrorInUp->DrawCopy("b");
2079     relativeErrorInLow->DrawCopy("same b");
2080     Style(relativeStatisticalErrorIn, kRed);
2081     relativeStatisticalErrorIn->DrawCopy("same");
2082     Style(AddLegend(gPad));
2083     relativeError->cd(2);
2084     Style(gPad, "GRID");
2085     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2086     Style(relativeErrorOutUp, kBlue, kBar);
2087     Style(relativeErrorOutLow, kGreen, kBar);
2088     relativeErrorOutUp->DrawCopy("b");
2089     relativeErrorOutLow->DrawCopy("same b");
2090     Style(relativeStatisticalErrorOut, kRed);
2091     relativeStatisticalErrorOut->DrawCopy("same");
2092     Style(AddLegend(gPad));
2093
2094     // write the buffered file to disk and close the file
2095     SavePadToPDF(relativeError);
2096     relativeError->Write();
2097     output->Write();
2098 //    output->Close();
2099 }
2100 //_____________________________________________________________________________
2101     void AliJetFlowTools::DoIntermediateSystematics(
2102             TArrayI* variationsIn,                  // variantions in plane
2103             TArrayI* variationsOut,                 // variantions out of plane
2104             TH1D*& relativeErrorInUp,               // pointer reference to minimum relative error histogram in plane
2105             TH1D*& relativeErrorInLow,              // pointer reference to maximum relative error histogram in plane
2106             TH1D*& relativeErrorOutUp,              // pointer reference to minimum relative error histogram out of plane
2107             TH1D*& relativeErrorOutLow,             // pointer reference to maximum relative error histogram out of plane
2108             TH1D*& relativeStatisticalErrorIn,      // relative systematic error on ratio
2109             TH1D*& relativeStatisticalErrorOut,     // relative systematic error on ratio
2110             TH1D*& nominal,                         // clone of the nominal ratio in / out of plane
2111             TH1D*& nominalIn,                       // clone of the nominal in plane yield
2112             TH1D*& nominalOut,                      // clone of the nominal out of plane yield
2113             Int_t columns,                          // divide the output canvasses in this many columns
2114             Float_t rangeLow,                       // lower pt range
2115             Float_t rangeUp,                        // upper pt range
2116             TFile* readMe,                          // input file name (created by this unfolding class)
2117             TString source                          // source of the variation
2118             ) const
2119 {
2120    // intermediate systematic check function. first index of supplied array is nominal value
2121    //
2122    TList* listOfKeys((TList*)readMe->GetListOfKeys());
2123    if(!listOfKeys) {
2124        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2125        return;
2126    }
2127    // check input params
2128    if(variationsIn->GetSize() != variationsOut->GetSize()) {
2129        printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2130        return;
2131    }
2132    TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2133    TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2134    if(!(defRootDirIn && defRootDirOut)) {
2135        printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2136        return;
2137    }
2138    TString defIn(defRootDirIn->GetName());
2139    TString defOut(defRootDirOut->GetName());
2140
2141    // define lines to make the output prettier
2142    TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2143    TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2144    lineLow->SetLineColor(11);
2145    lineUp->SetLineColor(11);
2146    lineLow->SetLineWidth(3);
2147    lineUp->SetLineWidth(3);
2148
2149    // define an output histogram with the maximum relative error from this systematic constribution
2150    relativeErrorInUp = new TH1D(Form("max #sigma/|x| from %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2151    relativeErrorInLow = new TH1D(Form("min #sigma/|x| from  %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2152    relativeErrorOutUp = new TH1D(Form("max #sigma/|x| from  %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2153    relativeErrorOutLow = new TH1D(Form("min #sigma/|x| from %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2154    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2155        relativeErrorInUp->SetBinContent(b+1, 1.);
2156        relativeErrorInUp->SetBinError(b+1, 0.);
2157        relativeErrorOutUp->SetBinContent(b+1, 1.);
2158        relativeErrorOutUp->SetBinError(b+1, .0);
2159        relativeErrorInLow->SetBinContent(b+1, 1.);
2160        relativeErrorInLow->SetBinError(b+1, 0.);
2161        relativeErrorOutLow->SetBinContent(b+1, 1.);
2162        relativeErrorOutLow->SetBinError(b+1, .0);
2163    }
2164    // define an output histogram with the systematic error from this systematic constribution
2165    if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2166        relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "#sigma/|x|, statistical, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2167        relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "#sigma/|x|, statistical, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2168    }
2169
2170    // prepare necessary canvasses
2171    TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2172    TCanvas* canvasOut(0x0);
2173    if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2174    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2175    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2176    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2177    TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
2178    TCanvas* canvasSpectraOut(0x0);
2179    if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2180    TCanvas* canvasRatio(0x0);
2181    if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2182    TCanvas* canvasV2(0x0);
2183    if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2184    TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2185    TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2186    TCanvas* canvasMasterOut(0x0);
2187    if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2188    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2189
2190    TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2191    canvasProfiles->Divide(2);
2192    TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2193    TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2194    // get an estimate of the number of outputs and find the default set
2195
2196    Int_t rows = 1;
2197    columns = variationsIn->GetSize()-1;
2198    (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2199    canvasIn->Divide(columns, rows);
2200    if(canvasOut) canvasOut->Divide(columns, rows);
2201    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2202    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2203    canvasSpectraIn->Divide(columns, rows);
2204    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2205    if(canvasRatio) canvasRatio->Divide(columns, rows);
2206    if(canvasV2) canvasV2->Divide(columns, rows);
2207    canvasMasterIn->Divide(columns, rows);
2208    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2209
2210    // prepare a separate set of canvases to hold the nominal points
2211    TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2212    TCanvas* canvasNominalOut(0x0);
2213    if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2214    TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2215    TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2216    if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2217    TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
2218    TCanvas* canvasNominalSpectraOut(0x0);
2219    if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()),  Form("Nominal_%s_SpectraOut", source.Data()));
2220    TCanvas* canvasNominalRatio(0x0);
2221    if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2222    TCanvas* canvasNominalV2(0x0);
2223    if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2224    TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2225    TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2226    TCanvas* canvasNominalMasterOut(0x0);
2227    if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2228    (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2229
2230    canvasNominalSpectraIn->Divide(2);
2231    if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2232
2233    canvasNominalMasterIn->Divide(2);
2234    if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2235
2236    // extract the default output
2237    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2238    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2239    TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2240    TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2241    if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2242    if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2243    printf(" > succesfully extracted default results < \n");
2244
2245    // loop through the directories, only plot the graphs if the deconvolution converged
2246    TDirectoryFile* tempDirIn(0x0);
2247    TDirectoryFile* tempDirOut(0x0);
2248    TDirectoryFile* tempIn(0x0);
2249    TDirectoryFile* tempOut(0x0);
2250    TH1D* unfoldedSpectrumInForRatio(0x0);
2251    TH1D* unfoldedSpectrumOutForRatio(0x0);
2252    for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2253        tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2254        tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2255        if(!(tempDirIn && tempDirOut)) {
2256            printf(" > DoSystematics: couldn't get a set of variations < \n");
2257            continue;
2258        }
2259        TString dirNameIn(tempDirIn->GetName());
2260        TString dirNameOut(tempDirOut->GetName());
2261        // try to read the in- and out of plane subdirs
2262        tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2263        tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2264        j++;
2265        if(tempIn) {
2266            // to see if the unfolding converged try to extract the pearson coefficients
2267            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2268            if(pIn) {
2269                printf(" - %s in plane converged \n", dirNameIn.Data());
2270                canvasIn->cd(j);
2271                if(i==0) canvasNominalIn->cd(j);
2272                Style(gPad, "PEARSON");
2273                pIn->DrawCopy("colz");
2274                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2275                if(rIn) {
2276                    Style(rIn);
2277                    printf(" > found RatioRefoldedMeasured < \n");
2278                    canvasRatioMeasuredRefoldedIn->cd(j);
2279                    if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2280                    Style(gPad, "GRID");
2281                    rIn->SetFillColor(kRed);
2282                    rIn->Draw("ap");
2283                }
2284                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2285                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2286                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2287                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2288                if(dvector && avalue && rm && eff) {
2289                    Style(dvector);
2290                    Style(avalue);
2291                    Style(rm);
2292                    Style(eff);
2293                    canvasMISC->cd(1);
2294                    if(i==0) canvasNominalMISC->cd(1);
2295                    Style(gPad, "SPECTRUM");
2296                    dvector->DrawCopy();
2297                    canvasMISC->cd(2);
2298                    if(i==0) canvasNominalMISC->cd(2);
2299                    Style(gPad, "SPECTRUM");
2300                    avalue->DrawCopy();
2301                    canvasMISC->cd(3);
2302                    if(i==0) canvasNominalMISC->cd(3);
2303                    Style(gPad, "PEARSON");
2304                    rm->DrawCopy("colz");
2305                    canvasMISC->cd(4);
2306                    if(i==0) canvasNominalMISC->cd(4);
2307                    Style(gPad, "GRID");
2308                    eff->DrawCopy();
2309                } else if(rm && eff) {
2310                    Style(rm);
2311                    Style(eff);
2312                    canvasMISC->cd(3);
2313                    if(i==0) canvasNominalMISC->cd(3);
2314                    Style(gPad, "PEARSON");
2315                    rm->DrawCopy("colz");
2316                    canvasMISC->cd(4);
2317                    if(i==0) canvasNominalMISC->cd(4);
2318                    Style(gPad, "GRID");
2319                    eff->DrawCopy();
2320                }
2321            }
2322            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2323            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2324            unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2325            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2326            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2327                if(defaultUnfoldedJetSpectrumIn) {
2328                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2329                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2330                    temp->Divide(unfoldedSpectrum);
2331                    // get the absolute relative error
2332                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2333                        // check if the error is larger than the current maximum
2334                        if( temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInUp->GetBinContent(b+1)) {
2335                            relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2336                            relativeErrorInUp->SetBinError(b+1, 0.);
2337                        }
2338                        // check if the error is smaller than the current minimum
2339                        else if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInLow->GetBinContent(b+1)) {
2340                            relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2341                            relativeErrorInLow->SetBinError(b+1, 0.);
2342                        }
2343                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2344                    }
2345                    temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2346                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2347                    temp->GetYaxis()->SetTitle("ratio");
2348                    canvasMasterIn->cd(j);
2349                    temp->GetYaxis()->SetRangeUser(0., 2);
2350                    Style(gPad, "GRID");
2351                    temp->DrawCopy();
2352                    canvasNominalMasterIn->cd(1);
2353                    Style(gPad, "GRID");
2354                    if(i > 0 ) {
2355                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2356                        tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2357                        Style(tempSyst, (EColor)(i+2));
2358                        if(i==1) tempSyst->DrawCopy();
2359                        else tempSyst->DrawCopy("same");
2360                    }
2361                }
2362                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2363                canvasSpectraIn->cd(j);
2364                if(i==0) canvasNominalSpectraIn->cd(1);
2365                Style(gPad);
2366                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2367                unfoldedSpectrum->DrawCopy();
2368                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2369                inputSpectrum->DrawCopy("same");
2370                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2371                refoldedSpectrum->DrawCopy("same");
2372                TLegend* l(AddLegend(gPad));
2373                Style(l);
2374                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2375                    Float_t chi(fitStatus->GetBinContent(1));
2376                    Float_t pen(fitStatus->GetBinContent(2));
2377                    Int_t dof((int)fitStatus->GetBinContent(3));
2378                    Float_t beta(fitStatus->GetBinContent(4));
2379                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2380                } else if (fitStatus) { // only available in SVD fit
2381                    Int_t reg((int)fitStatus->GetBinContent(1));
2382                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2383                }
2384                canvasNominalSpectraIn->cd(2);
2385                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2386                tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2387                Style(tempSyst, (EColor)(i+2));
2388                Style(gPad, "SPECTRUM");
2389                if(i==0) tempSyst->DrawCopy();
2390                else tempSyst->DrawCopy("same");
2391            }
2392        }
2393        if(tempOut) {
2394            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2395            if(pOut) {
2396                printf(" - %s out of plane converged \n", dirNameOut.Data());
2397                canvasOut->cd(j);
2398                if(i==0) canvasNominalOut->cd(j);
2399                Style(gPad, "PEARSON");
2400                pOut->DrawCopy("colz");
2401                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2402                if(rOut) {
2403                    Style(rOut);
2404                    printf(" > found RatioRefoldedMeasured < \n");
2405                    canvasRatioMeasuredRefoldedOut->cd(j);
2406                    if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2407                    Style(gPad, "GRID");
2408                    rOut->SetFillColor(kRed);
2409                    rOut->Draw("ap");
2410                }
2411                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2412                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2413                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2414                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2415                if(dvector && avalue && rm && eff) {
2416                    Style(dvector);
2417                    Style(avalue);
2418                    Style(rm);
2419                    Style(eff);
2420                    canvasMISC->cd(5);
2421                    if(i==0) canvasNominalMISC->cd(5);
2422                    Style(gPad, "SPECTRUM");
2423                    dvector->DrawCopy();
2424                    canvasMISC->cd(6);
2425                    if(i==0) canvasNominalMISC->cd(6);
2426                    Style(gPad, "SPECTRUM");
2427                    avalue->DrawCopy();
2428                    canvasMISC->cd(7);
2429                    if(i==0) canvasNominalMISC->cd(7);
2430                    Style(gPad, "PEARSON");
2431                    rm->DrawCopy("colz");
2432                    canvasMISC->cd(8);
2433                    if(i==0) canvasNominalMISC->cd(8);
2434                    Style(gPad, "GRID");
2435                    eff->DrawCopy();
2436                } else if(rm && eff) {
2437                    Style(rm);
2438                    Style(eff);
2439                    canvasMISC->cd(7);
2440                    if(i==0) canvasNominalMISC->cd(7);
2441                    Style(gPad, "PEARSON");
2442                    rm->DrawCopy("colz");
2443                    canvasMISC->cd(8);
2444                    if(i==0) canvasNominalMISC->cd(8);
2445                    Style(gPad, "GRID");
2446                    eff->DrawCopy();
2447                }
2448            }
2449            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2450            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2451            unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2452            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2453            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2454                if(defaultUnfoldedJetSpectrumOut) {
2455                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2456                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2457                    temp->Divide(unfoldedSpectrum);
2458                    // get the absolute relative error
2459                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2460                        // check if the error is larger than the current maximum
2461                        if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutUp->GetBinContent(b+1)) {
2462                            relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2463                            relativeErrorOutUp->SetBinError(b+1, 0.);
2464                        }
2465                        // check if the error is smaller than the current minimum
2466                        else if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutLow->GetBinContent(b+1)) {
2467                            relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2468                            relativeErrorOutLow->SetBinError(b+1, 0.);
2469                        }
2470                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2471                     }
2472                    temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2473                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2474                    temp->GetYaxis()->SetTitle("ratio");
2475                    canvasMasterOut->cd(j);
2476                    temp->GetYaxis()->SetRangeUser(0., 2);
2477                    Style(gPad, "GRID");
2478                    temp->DrawCopy();
2479                    canvasNominalMasterOut->cd(1);
2480                    Style(gPad, "GRID");
2481                    if(i > 0 ) {
2482                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2483                        tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2484                        Style(tempSyst, (EColor)(i+2));
2485                        if(i==1) tempSyst->DrawCopy();
2486                        else tempSyst->DrawCopy("same");
2487                    }
2488                }
2489                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2490                canvasSpectraOut->cd(j);
2491                if(i==0) canvasNominalSpectraOut->cd(1);
2492                Style(gPad);
2493                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2494                unfoldedSpectrum->DrawCopy();
2495                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2496                inputSpectrum->DrawCopy("same");
2497                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2498                refoldedSpectrum->DrawCopy("same");
2499                TLegend* l(AddLegend(gPad));
2500                Style(l);
2501                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2502                    Float_t chi(fitStatus->GetBinContent(1));
2503                    Float_t pen(fitStatus->GetBinContent(2));
2504                    Int_t dof((int)fitStatus->GetBinContent(3));
2505                    Float_t beta(fitStatus->GetBinContent(4));
2506                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2507                } else if (fitStatus) { // only available in SVD fit
2508                    Int_t reg((int)fitStatus->GetBinContent(1));
2509                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2510                }
2511                canvasNominalSpectraOut->cd(2);
2512                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2513                tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2514                Style(tempSyst, (EColor)(i+2));
2515                Style(gPad, "SPECTRUM");
2516                if(i==0) tempSyst->DrawCopy();
2517                else tempSyst->DrawCopy("same");
2518            }
2519        }
2520        if(canvasRatio && canvasV2) {
2521            canvasRatio->cd(j);
2522            if(i==0) {
2523                canvasNominalRatio->cd(j);
2524                if(nominal && nominalIn && nominalOut) {
2525                    // if a nominal ratio is requested, delete the placeholder and update the nominal point
2526                    delete nominal;
2527                    delete nominalIn;
2528                    delete nominalOut;
2529                    nominalIn =  (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2530                    nominalOut =  (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2531                    nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2532                    nominal->Divide(nominal, unfoldedSpectrumOutForRatio);       // standard root divide for uncorrelated histos
2533                }
2534            }
2535            TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2536            Double_t _tempx(0), _tempy(0);
2537            if(ratioYield) {
2538                Style(ratioYield);
2539                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2540                    ratioYield->GetPoint(b, _tempx, _tempy);
2541                    ratioProfile->Fill(_tempx, _tempy);
2542                }
2543                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2544                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2545                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2546                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2547                ratioYield->SetFillColor(kRed);
2548                ratioYield->Draw("ap");
2549            }
2550            canvasV2->cd(j);
2551            if(i==0) canvasNominalV2->cd(j);
2552            TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2553            if(ratioV2) {
2554                Style(ratioV2);
2555                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2556                    ratioV2->GetPoint(b, _tempx, _tempy);
2557                    v2Profile->Fill(_tempx, _tempy);
2558
2559                }
2560                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2561                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2562                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2563                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2564                ratioV2->SetFillColor(kRed);
2565                ratioV2->Draw("ap");
2566            }
2567        }
2568        delete unfoldedSpectrumInForRatio;
2569        delete unfoldedSpectrumOutForRatio;
2570    }
2571    // save the canvasses
2572    canvasProfiles->cd(1);
2573    Style(ratioProfile);
2574    ratioProfile->DrawCopy();
2575    canvasProfiles->cd(2);
2576    Style(v2Profile);
2577    v2Profile->DrawCopy();
2578    SavePadToPDF(canvasProfiles);
2579    canvasProfiles->Write();
2580    SavePadToPDF(canvasIn);
2581    canvasIn->Write();
2582    if(canvasOut) {
2583        SavePadToPDF(canvasOut);
2584        canvasOut->Write();
2585    }
2586    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2587    canvasRatioMeasuredRefoldedIn->Write();
2588    if(canvasRatioMeasuredRefoldedOut) {
2589        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2590        canvasRatioMeasuredRefoldedOut->Write();
2591    }
2592    SavePadToPDF(canvasSpectraIn);
2593    canvasSpectraIn->Write();
2594    if(canvasSpectraOut) {
2595        SavePadToPDF(canvasSpectraOut);
2596        canvasSpectraOut->Write();
2597        SavePadToPDF(canvasRatio);
2598        canvasRatio->Write();
2599        SavePadToPDF(canvasV2);
2600        canvasV2->Write();
2601    }
2602    SavePadToPDF(canvasMasterIn);
2603    canvasMasterIn->Write();
2604    if(canvasMasterOut) {
2605        SavePadToPDF(canvasMasterOut);
2606        canvasMasterOut->Write();
2607    }
2608    SavePadToPDF(canvasMISC);
2609    canvasMISC->Write();
2610    // save the nomial canvasses
2611    SavePadToPDF(canvasNominalIn);
2612    canvasNominalIn->Write();
2613    if(canvasNominalOut) {
2614        SavePadToPDF(canvasNominalOut);
2615        canvasNominalOut->Write();
2616    }
2617    SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
2618    canvasNominalRatioMeasuredRefoldedIn->Write();
2619    if(canvasNominalRatioMeasuredRefoldedOut) {
2620        SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
2621        canvasNominalRatioMeasuredRefoldedOut->Write();
2622    }
2623    canvasNominalSpectraIn->cd(2);
2624    Style(AddLegend(gPad));
2625    SavePadToPDF(canvasNominalSpectraIn);
2626    canvasNominalSpectraIn->Write();
2627    if(canvasNominalSpectraOut) {
2628        canvasNominalSpectraOut->cd(2);
2629        Style(AddLegend(gPad));
2630        SavePadToPDF(canvasNominalSpectraOut);
2631        canvasNominalSpectraOut->Write();
2632        SavePadToPDF(canvasNominalRatio);
2633        canvasNominalRatio->Write();
2634        SavePadToPDF(canvasNominalV2);
2635        canvasNominalV2->Write();
2636    }
2637    canvasNominalMasterIn->cd(1);
2638    Style(AddLegend(gPad));
2639    lineUp->DrawClone("same");
2640    lineLow->DrawClone("same");
2641    SavePadToPDF(canvasNominalMasterIn);
2642    canvasNominalMasterIn->Write();
2643    if(canvasNominalMasterOut) {
2644        canvasNominalMasterOut->cd(1);
2645        Style(AddLegend(gPad));
2646        lineUp->DrawClone("same");
2647        lineLow->DrawClone("same");
2648        SavePadToPDF(canvasNominalMasterOut);
2649        canvasNominalMasterOut->Write();
2650    }
2651    SavePadToPDF(canvasNominalMISC);
2652    canvasNominalMISC->Write();
2653
2654    // save the relative errors
2655    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2656        relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)-1);
2657        relativeErrorInUp->SetBinError(b+1, 0.);
2658        relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)-1);
2659        relativeErrorOutUp->SetBinError(b+1, .0);
2660        relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)-1);
2661        relativeErrorInLow->SetBinError(b+1, 0.);
2662        relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)-1);
2663        relativeErrorOutLow->SetBinError(b+1, .0);
2664    }
2665    relativeErrorInUp->SetYTitle("relative uncertainty");
2666    relativeErrorOutUp->SetYTitle("relative uncertainty");
2667    relativeErrorInLow->SetYTitle("relative uncertainty");
2668    relativeErrorOutLow->SetYTitle("relative uncertainty");
2669    relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2670    relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2671    relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2672    relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2673
2674    canvasNominalMasterIn->cd(2);
2675    Style(gPad, "GRID");
2676    Style(relativeErrorInUp, kBlue, kBar);
2677    Style(relativeErrorInLow, kGreen, kBar);
2678    relativeErrorInUp->DrawCopy("b");
2679    relativeErrorInLow->DrawCopy("same b");
2680    Style(AddLegend(gPad));
2681    SavePadToPDF(canvasNominalMasterIn);
2682    canvasNominalMasterIn->Write();
2683    canvasNominalMasterOut->cd(2);
2684    Style(gPad, "GRID");
2685    Style(relativeErrorOutUp, kBlue, kBar);
2686    Style(relativeErrorOutLow, kGreen, kBar);
2687    relativeErrorOutUp->DrawCopy("b");
2688    relativeErrorOutLow->DrawCopy("same b");
2689    Style(AddLegend(gPad));
2690    SavePadToPDF(canvasNominalMasterOut);
2691    canvasNominalMasterOut->Write();
2692 }
2693 //_____________________________________________________________________________
2694 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
2695 {
2696    // go through the output file and perform post processing routines
2697    // can either be performed in one go with the unfolding, or at a later stage
2698    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
2699    TFile readMe(in.Data(), "READ");     // open file read-only
2700    if(readMe.IsZombie()) {
2701        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2702        return;
2703    }
2704    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
2705    readMe.ls();
2706    TList* listOfKeys((TList*)readMe.GetListOfKeys());
2707    if(!listOfKeys) {
2708        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2709        return;
2710    }
2711    // prepare necessary canvasses
2712    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
2713    TCanvas* canvasOut(0x0);
2714    if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
2715    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
2716    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2717    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
2718    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
2719    TCanvas* canvasSpectraOut(0x0);
2720    if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
2721    TCanvas* canvasRatio(0x0);
2722    if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
2723    TCanvas* canvasV2(0x0);
2724    if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
2725    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
2726    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
2727    TCanvas* canvasMasterOut(0x0);
2728    if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
2729    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2730    TDirectoryFile* defDir(0x0);
2731    TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
2732    canvasProfiles->Divide(2);
2733    TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2734    TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2735    // get an estimate of the number of outputs and find the default set
2736    Int_t cacheMe(0);
2737    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
2738        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
2739            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2740            cacheMe++;
2741        }
2742    }
2743    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
2744    canvasIn->Divide(columns, rows);
2745    if(canvasOut) canvasOut->Divide(columns, rows);
2746    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2747    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2748    canvasSpectraIn->Divide(columns, rows);
2749    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2750    if(canvasRatio) canvasRatio->Divide(columns, rows);
2751    if(canvasV2) canvasV2->Divide(columns, rows);
2752
2753    canvasMasterIn->Divide(columns, rows);
2754    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2755    // extract the default output
2756    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2757    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2758    if(defDir) {
2759        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
2760        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
2761        if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
2762        if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
2763        printf(" > succesfully extracted default results < \n");
2764    }
2765
2766    // loop through the directories, only plot the graphs if the deconvolution converged
2767    TDirectoryFile* tempDir(0x0);
2768    TDirectoryFile* tempIn(0x0);
2769    TDirectoryFile*  tempOut(0x0);
2770    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
2771        // read the directory by its name
2772        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2773        if(!tempDir) continue;
2774        TString dirName(tempDir->GetName());
2775        // try to read the in- and out of plane subdirs
2776        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
2777        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
2778        j++;
2779        if(!tempIn) {    // attempt to get MakeAU output
2780            TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
2781            TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
2782            tempPearson->Divide(4,2);
2783            TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
2784            tempRatio->Divide(4,2);
2785            TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
2786            tempResult->Divide(4,2);
2787            for(Int_t q(0); q < 8; q++) {
2788                tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
2789                if(tempIn) {
2790                        // to see if the unfolding converged try to extract the pearson coefficients
2791                        TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2792                        if(pIn) {
2793                        printf(" - %s in plane converged \n", dirName.Data());
2794                            tempPearson->cd(1+q);
2795                             Style(gPad, "PEARSON");
2796                             pIn->DrawCopy("colz");
2797                             TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2798                             if(rIn) {
2799                                 Style(rIn);
2800                                 printf(" > found RatioRefoldedMeasured < \n");
2801                                 tempRatio->cd(q+1);
2802                                 rIn->SetFillColor(kRed);
2803                                 rIn->Draw("ap");
2804                             }
2805                             TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2806                             TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2807                             TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2808                             TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2809                             if(dvector && avalue && rm && eff) {
2810                                 Style(dvector);
2811                                 Style(avalue);
2812                                 Style(rm);
2813                                 Style(eff);
2814                                 canvasMISC->cd(1);
2815                                 Style(gPad, "SPECTRUM");
2816                                 dvector->DrawCopy();
2817                                 canvasMISC->cd(2);
2818                                 Style(gPad, "SPECTRUM");
2819                                 avalue->DrawCopy();
2820                                 canvasMISC->cd(3);
2821                                 Style(gPad, "PEARSON");
2822                                 rm->DrawCopy("colz");
2823                                 canvasMISC->cd(4);
2824                                 Style(gPad, "GRID");
2825                                 eff->DrawCopy();
2826                             } else if(rm && eff) {
2827                                 Style(rm);
2828                                 Style(eff);
2829                                 canvasMISC->cd(3);
2830                                 Style(gPad, "PEARSON");
2831                                 rm->DrawCopy("colz");
2832                                 canvasMISC->cd(4);
2833                                 Style(gPad, "GRID");
2834                                 eff->DrawCopy();
2835                             }
2836                         }
2837                        TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2838                        TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2839                        TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2840                        if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2841                            if(defaultUnfoldedJetSpectrumIn) {
2842                                Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2843                                TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
2844                                temp->Divide(unfoldedSpectrum);
2845                                temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
2846                                temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2847                                temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2848                                canvasMasterIn->cd(j);
2849                                temp->GetYaxis()->SetRangeUser(0., 2);
2850                                temp->DrawCopy();
2851                            }
2852                            TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2853                            tempResult->cd(q+1);
2854                            Style(gPad);
2855                            Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2856                            unfoldedSpectrum->DrawCopy();
2857                            Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2858                            inputSpectrum->DrawCopy("same");
2859                            Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2860                            refoldedSpectrum->DrawCopy("same");
2861                            TLegend* l(AddLegend(gPad));
2862                            Style(l);
2863                            if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2864                                Float_t chi(fitStatus->GetBinContent(1));
2865                                Float_t pen(fitStatus->GetBinContent(2));
2866                                Int_t dof((int)fitStatus->GetBinContent(3));
2867                                Float_t beta(fitStatus->GetBinContent(4));
2868                                l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2869                            } else if (fitStatus) { // only available in SVD fit
2870                                Int_t reg((int)fitStatus->GetBinContent(1));
2871                                l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2872                            }
2873                        }
2874                }
2875            }
2876        }
2877        if(tempIn) {
2878            // to see if the unfolding converged try to extract the pearson coefficients
2879            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2880            if(pIn) {
2881                printf(" - %s in plane converged \n", dirName.Data());
2882                canvasIn->cd(j);
2883                Style(gPad, "PEARSON");
2884                pIn->DrawCopy("colz");
2885                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2886                if(rIn) {
2887                    Style(rIn);
2888                    printf(" > found RatioRefoldedMeasured < \n");
2889                    canvasRatioMeasuredRefoldedIn->cd(j);
2890                    rIn->SetFillColor(kRed);
2891                    rIn->Draw("ap");
2892                }
2893                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2894                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2895                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2896                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2897                if(dvector && avalue && rm && eff) {
2898                    Style(dvector);
2899                    Style(avalue);
2900                    Style(rm);
2901                    Style(eff);
2902                    canvasMISC->cd(1);
2903                    Style(gPad, "SPECTRUM");
2904                    dvector->DrawCopy();
2905                    canvasMISC->cd(2);
2906                    Style(gPad, "SPECTRUM");
2907                    avalue->DrawCopy();
2908                    canvasMISC->cd(3);
2909                    Style(gPad, "PEARSON");
2910                    rm->DrawCopy("colz");
2911                    canvasMISC->cd(4);
2912                    Style(gPad, "GRID");
2913                    eff->DrawCopy();
2914                } else if(rm && eff) {
2915                    Style(rm);
2916                    Style(eff);
2917                    canvasMISC->cd(3);
2918                    Style(gPad, "PEARSON");
2919                    rm->DrawCopy("colz");
2920                    canvasMISC->cd(4);
2921                    Style(gPad, "GRID");
2922                    eff->DrawCopy();
2923                }
2924            }
2925            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2926            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2927            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2928            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2929                if(defaultUnfoldedJetSpectrumIn) {
2930                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2931                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
2932                    temp->Divide(unfoldedSpectrum);
2933                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
2934                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2935                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2936                    canvasMasterIn->cd(j);
2937                    temp->GetYaxis()->SetRangeUser(0., 2);
2938                    temp->DrawCopy();
2939                }
2940                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2941                canvasSpectraIn->cd(j);
2942                Style(gPad);
2943                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2944                unfoldedSpectrum->DrawCopy();
2945                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2946                inputSpectrum->DrawCopy("same");
2947                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2948                refoldedSpectrum->DrawCopy("same");
2949                TLegend* l(AddLegend(gPad));
2950                Style(l);
2951                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2952                    Float_t chi(fitStatus->GetBinContent(1));
2953                    Float_t pen(fitStatus->GetBinContent(2));
2954                    Int_t dof((int)fitStatus->GetBinContent(3));
2955                    Float_t beta(fitStatus->GetBinContent(4));
2956                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2957                } else if (fitStatus) { // only available in SVD fit
2958                    Int_t reg((int)fitStatus->GetBinContent(1));
2959                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2960                }
2961            }
2962        }
2963        if(tempOut) {
2964            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
2965            if(pOut) {
2966                printf(" - %s out of plane converged \n", dirName.Data());
2967                canvasOut->cd(j);
2968                Style(gPad, "PEARSON");
2969                pOut->DrawCopy("colz");
2970                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2971                if(rOut) {
2972                    Style(rOut);
2973                    printf(" > found RatioRefoldedMeasured < \n");
2974                    canvasRatioMeasuredRefoldedOut->cd(j);
2975                    rOut->SetFillColor(kRed);
2976                    rOut->Draw("ap");
2977                }
2978                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2979                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2980                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
2981                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
2982                if(dvector && avalue && rm && eff) {
2983                    Style(dvector);
2984                    Style(avalue);
2985                    Style(rm);
2986                    Style(eff);
2987                    canvasMISC->cd(5);
2988                    Style(gPad, "SPECTRUM");
2989                    dvector->DrawCopy();
2990                    canvasMISC->cd(6);
2991                    Style(gPad, "SPECTRUM");
2992                    avalue->DrawCopy();
2993                    canvasMISC->cd(7);
2994                    Style(gPad, "PEARSON");
2995                    rm->DrawCopy("colz");
2996                    canvasMISC->cd(8);
2997                    Style(gPad, "GRID");
2998                    eff->DrawCopy();
2999                } else if(rm && eff) {
3000                    Style(rm);
3001                    Style(eff);
3002                    canvasMISC->cd(7);
3003                    Style(gPad, "PEARSON");
3004                    rm->DrawCopy("colz");
3005                    canvasMISC->cd(8);
3006                    Style(gPad, "GRID");
3007                    eff->DrawCopy();
3008                }
3009            }
3010            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3011            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3012            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3013            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3014                if(defaultUnfoldedJetSpectrumOut) {
3015                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3016                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3017                    temp->Divide(unfoldedSpectrum);
3018                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3019                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3020                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3021                    canvasMasterOut->cd(j);
3022                    temp->GetYaxis()->SetRangeUser(0., 2.);
3023                    temp->DrawCopy();
3024                }
3025                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3026                canvasSpectraOut->cd(j);
3027                Style(gPad);
3028                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3029                unfoldedSpectrum->DrawCopy();
3030                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3031                inputSpectrum->DrawCopy("same");
3032                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3033                refoldedSpectrum->DrawCopy("same");
3034                TLegend* l(AddLegend(gPad));
3035                Style(l);
3036                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3037                    Float_t chi(fitStatus->GetBinContent(1));
3038                    Float_t pen(fitStatus->GetBinContent(2));
3039                    Int_t dof((int)fitStatus->GetBinContent(3));
3040                    Float_t beta(fitStatus->GetBinContent(4));
3041                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3042                } else if (fitStatus) { // only available in SVD fit
3043                    Int_t reg((int)fitStatus->GetBinContent(1));
3044                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3045                }
3046            }
3047        }
3048        if(canvasRatio && canvasV2) {
3049            canvasRatio->cd(j);
3050            TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
3051            Double_t _tempx(0), _tempy(0);
3052            if(ratioYield) {
3053                Style(ratioYield);
3054                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3055                    ratioYield->GetPoint(b, _tempx, _tempy);
3056                    ratioProfile->Fill(_tempx, _tempy);
3057                }
3058                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3059                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3060                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3061                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3062                ratioYield->SetFillColor(kRed);
3063                ratioYield->Draw("ap");
3064            }
3065            canvasV2->cd(j);
3066            TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3067            if(ratioV2) {
3068                Style(ratioV2);
3069                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3070                    ratioV2->GetPoint(b, _tempx, _tempy);
3071                    v2Profile->Fill(_tempx, _tempy);
3072
3073                }
3074                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3075                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3076                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3077                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3078                ratioV2->SetFillColor(kRed);
3079                ratioV2->Draw("ap");
3080            }
3081        }
3082    }
3083    TFile output(out.Data(), "RECREATE");
3084    canvasProfiles->cd(1);
3085    Style(ratioProfile);
3086    ratioProfile->DrawCopy();
3087    canvasProfiles->cd(2);
3088    Style(v2Profile);
3089    v2Profile->DrawCopy();
3090    SavePadToPDF(canvasProfiles);
3091    canvasProfiles->Write();
3092    SavePadToPDF(canvasIn);
3093    canvasIn->Write();
3094    if(canvasOut) {
3095        SavePadToPDF(canvasOut);
3096        canvasOut->Write();
3097    }
3098    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3099    canvasRatioMeasuredRefoldedIn->Write();
3100    if(canvasRatioMeasuredRefoldedOut) {
3101        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3102        canvasRatioMeasuredRefoldedOut->Write();
3103    }
3104    SavePadToPDF(canvasSpectraIn);
3105    canvasSpectraIn->Write();
3106    if(canvasSpectraOut) {
3107        SavePadToPDF(canvasSpectraOut);
3108        canvasSpectraOut->Write();
3109        SavePadToPDF(canvasRatio);
3110        canvasRatio->Write();
3111        SavePadToPDF(canvasV2);
3112        canvasV2->Write();
3113    }
3114    SavePadToPDF(canvasMasterIn);
3115    canvasMasterIn->Write();
3116    if(canvasMasterOut) {
3117        SavePadToPDF(canvasMasterOut);
3118        canvasMasterOut->Write();
3119    }
3120    SavePadToPDF(canvasMISC);
3121    canvasMISC->Write();
3122    output.Write();
3123    output.Close();
3124 }
3125 //_____________________________________________________________________________
3126 Bool_t AliJetFlowTools::SetRawInput (
3127         TH2D* detectorResponse,  // detector response matrix
3128         TH1D* jetPtIn,           // in plane jet spectrum
3129         TH1D* jetPtOut,          // out of plane jet spectrum
3130         TH1D* dptIn,             // in plane delta pt distribution
3131         TH1D* dptOut,            // out of plane delta pt distribution
3132         Int_t eventCount) {
3133     // set input histograms manually
3134     fDetectorResponse   = detectorResponse;
3135     fSpectrumIn         = jetPtIn;
3136     fSpectrumOut        = jetPtOut;
3137     fDptInDist          = dptIn;
3138     fDptOutDist         = dptOut;
3139     fRawInputProvided   = kTRUE;
3140     // check if all data is provided
3141     if(!fDetectorResponse) {
3142         printf(" fDetectorResponse not found \n ");
3143         return kFALSE;
3144     }
3145     // check if the pt bin for true and rec have been set
3146     if(!fBinsTrue || !fBinsRec) {
3147         printf(" No true or rec bins set, please set binning ! \n");
3148         return kFALSE;
3149     }
3150     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
3151                         // procedures, these profiles will be nonsensical, user is responsible
3152         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3153         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3154         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3155     }
3156     // normalize spectra to event count if requested
3157     if(fNormalizeSpectra) {
3158         fEventCount = eventCount;
3159         if(fEventCount > 0) {
3160             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
3161             fSpectrumOut->Sumw2();
3162             fSpectrumIn->Scale(1./((double)fEventCount));
3163             fSpectrumOut->Scale(1./((double)fEventCount));
3164         }
3165     }
3166     if(!fNormalizeSpectra && fEventCount > 0) {
3167         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
3168         fSpectrumOut->Sumw2();
3169         fSpectrumIn->Scale(1./((double)fEventCount));
3170         fSpectrumOut->Scale(1./((double)fEventCount));
3171     }
3172     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
3173     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
3174     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3175     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3176     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
3177     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
3178     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3179     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3180
3181     return kTRUE;
3182 }
3183 //_____________________________________________________________________________
3184 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
3185 {
3186     // return ratio of h1 / h2
3187     // histograms can have different binning. errors are propagated as uncorrelated
3188     if(!(h1 && h2) ) {
3189         printf(" GetRatio called with NULL argument(s) \n ");
3190         return 0x0;
3191     }
3192     Int_t j(0);
3193     TGraphErrors *gr = new TGraphErrors();
3194     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3195     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3196     TH1* dud((TH1*)h1->Clone("dud"));
3197     dud->Sumw2();
3198     h1->Sumw2();
3199     h2->Sumw2();
3200     if(!dud->Divide(h1, h2)) {
3201         printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
3202         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3203             binCent = h1->GetXaxis()->GetBinCenter(i);
3204             if(xmax > 0. && binCent > xmax) continue;
3205             j = h2->FindBin(binCent);
3206             binWidth = h1->GetXaxis()->GetBinWidth(i);
3207             if(h2->GetBinContent(j) > 0.) {
3208                 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
3209             /* original propagation of uncertainty changed 08012014
3210             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
3211             Double_t B = 0.;
3212             if(h2->GetBinError(j)>0.) {
3213                 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
3214                 error2 = A*A + B*B;
3215             } else error2 = A*A;        */
3216                 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
3217                 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
3218                 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
3219                 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
3220                 gr->SetPoint(i-1, binCent, ratio);
3221                 gr->SetPointError(i-1, 0.5*binWidth, error2);
3222             }
3223         }
3224     } else {
3225         printf( " > using ROOT to divide two histograms < \n");
3226         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3227             binCent = dud->GetXaxis()->GetBinCenter(i);
3228             if(xmax > 0. && binCent > xmax) continue;
3229             binWidth = dud->GetXaxis()->GetBinWidth(i);
3230             gr->SetPoint(i-1,binCent,dud->GetBinContent(i));
3231             gr->SetPointError(i-1, 0.5*binWidth,dud->GetBinError(i));
3232         }
3233     }
3234
3235     if(appendFit) {
3236         TF1* fit(new TF1("lin", "pol0", 10, 100));
3237         gr->Fit(fit);
3238     }
3239     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3240     if(dud) delete dud;
3241     return gr;
3242 }
3243 //_____________________________________________________________________________
3244 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
3245 {
3246     // get v2 from difference of in plane, out of plane yield
3247     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
3248     // different binning is allowed but will mean that the error
3249     // propagation is unreliable
3250     // r is the event plane resolution for the chosen centrality
3251     if(!(h1 && h2) ) {
3252         printf(" GetV2 called with NULL argument(s) \n ");
3253         return 0x0;
3254     }
3255     TGraphErrors *gr = new TGraphErrors();
3256     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3257     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3258     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
3259
3260     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3261         binCent = h1->GetXaxis()->GetBinCenter(i);
3262         binWidth = h1->GetXaxis()->GetBinWidth(i);
3263         if(h2->GetBinContent(i) > 0.) {
3264             in = h1->GetBinContent(i);
3265             ein = h1->GetBinError(i);
3266             out = h2->GetBinContent(i);
3267             eout = h2->GetBinError(i);
3268             ratio = pre*((in-out)/(in+out));
3269             error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
3270             error2 = error2*pre*pre;
3271             if(error2 > 0) error2 = TMath::Sqrt(error2);
3272             gr->SetPoint(i-1,binCent,ratio);
3273             gr->SetPointError(i-1,0.5*binWidth,error2);
3274         }
3275     }
3276     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3277     return gr;
3278 }
3279 //_____________________________________________________________________________
3280 TGraphAsymmErrors* AliJetFlowTools::GetV2WithSystematicErrors(
3281         TH1* h1, TH1* h2, Double_t r, TString name,
3282         TH1* relativeErrorInUp,
3283         TH1* relativeErrorInLow,
3284         TH1* relativeErrorOutUp,
3285         TH1* relativeErrorOutLow,
3286         Float_t rho) const
3287 {
3288     // get v2 with asymmetric systematic error
3289     // note that this is ONLY the systematic error, no statistical error!
3290     // rho is the pearson correlation coefficient
3291     TGraphErrors* tempV2(GetV2(h1, h2, r, name));
3292     Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
3293     Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
3294     Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
3295     Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
3296     Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
3297     Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
3298     Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
3299     // loop through the bins and do error propagation
3300     for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3301         // extract the absolute errors
3302         in = h1->GetBinContent(i+1);
3303         einUp = in*relativeErrorInUp->GetBinContent(i+1);
3304         einLow = in*relativeErrorInLow->GetBinContent(1+i);
3305         out = h2->GetBinContent(i+1);
3306         eoutUp = out*relativeErrorOutUp->GetBinContent(1+i);
3307         eoutLow = out*relativeErrorOutLow->GetBinContent(1+i);
3308         // get the error squared
3309         error2Up = TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einUp*eoutUp);
3310         error2Low =TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einLow*eoutLow);
3311         if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
3312         if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
3313         // set the errors
3314         ayh[i] = error2Up;
3315         ayl[i] = error2Low;
3316         // get the bin width (which is the 'error' on x)
3317         Double_t binWidth(h1->GetBinWidth(i+1));
3318         axl[i] = binWidth/2.;
3319         axh[i] = binWidth/2.;
3320         // now get the coordinate for the poin
3321         tempV2->GetPoint(i, ax[i], ay[i]);
3322     }
3323     // save the nominal ratio
3324     TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
3325     nominalError->SetName("v_{2} with shape uncertainty");
3326     // do some memory management
3327     delete tempV2;
3328     delete[] ax;
3329     delete[] ay;
3330     delete[] axh;
3331     delete[] axl;
3332     delete[] ayh;
3333     delete[] ayl;
3334
3335     return nominalError;
3336 }
3337 //_____________________________________________________________________________
3338 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
3339     // write object, if a unique identifier is given the object is cloned
3340     // and the clone is saved. setting kill to true will delete the original obect from the heap
3341     if(!object) {
3342         printf(" > WriteObject:: called with NULL arguments \n ");
3343         return;
3344     } else if(!strcmp("", suffix.Data())) object->Write();
3345     else {
3346         TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
3347         newObject->Write();
3348     }
3349     if(kill) delete object;
3350 }
3351 //_____________________________________________________________________________
3352 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
3353     // construt a delta pt response matrix from supplied dpt distribution
3354     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
3355     // do a weighted rebinning to a (coarser) dpt distribution
3356     // be careful with the binning of the dpt response: it should be equal to that
3357     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
3358     //
3359     // the response matrix will be square and have the same binning
3360     // (min, max and granularity) of the input histogram
3361     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
3362     Double_t _bins[bins+1];             // prepare array with bin borders
3363     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
3364     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
3365     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
3366     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
3367         Bool_t skip = kFALSE;
3368         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
3369             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
3370             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
3371         }
3372     }
3373     return res;
3374 }
3375 //_____________________________________________________________________________
3376 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
3377     if(!binsTrue || !binsRec) {
3378         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
3379         return 0x0;
3380     }
3381     TString name(Form("unityResponse_%s", suffix.Data()));
3382     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
3383     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
3384         for(Int_t j(0); j < binsRec->GetSize(); j++) {
3385             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
3386         }
3387     }
3388     return unity;
3389 }
3390 //_____________________________________________________________________________
3391 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
3392     // save configuration parameters to histogram
3393     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
3394     summary->SetBinContent(1, fBetaIn);
3395     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
3396     summary->SetBinContent(2, fBetaOut);
3397     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
3398     summary->SetBinContent(3, fCentralityArray->At(0));
3399     summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]");
3400     summary->SetBinContent(4, (int)convergedIn);
3401     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
3402     summary->SetBinContent(5, (int)convergedOut);
3403     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
3404     summary->SetBinContent(6, (int)fAvoidRoundingError);
3405     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
3406     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
3407     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
3408     summary->SetBinContent(8, (int)fPrior);
3409     summary->GetXaxis()->SetBinLabel(8, "fPrior");
3410     summary->SetBinContent(9, fSVDRegIn);
3411     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
3412     summary->SetBinContent(10, fSVDRegOut);
3413     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
3414     summary->SetBinContent(11, (int)fSVDToy);
3415     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
3416     summary->SetBinContent(12, fJetRadius);
3417     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
3418     summary->SetBinContent(13, (int)fNormalizeSpectra);
3419     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
3420     summary->SetBinContent(14, (int)fSmoothenPrior);
3421     summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
3422     summary->SetBinContent(15, (int)fTestMode);
3423     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
3424     summary->SetBinContent(16, (int)fUseDetectorResponse);
3425     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
3426     summary->SetBinContent(17, fBayesianIterIn);
3427     summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
3428     summary->SetBinContent(18, fBayesianIterOut);
3429     summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
3430     summary->SetBinContent(19, fBayesianSmoothIn);
3431     summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
3432     summary->SetBinContent(20, fBayesianSmoothOut);
3433     summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
3434 }
3435 //_____________________________________________________________________________
3436 void AliJetFlowTools::ResetAliUnfolding() {
3437      // ugly function: reset all unfolding parameters
3438      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
3439      if(fitter) {
3440          printf(" > Found fitter, will delete it < \n");
3441          delete fitter;
3442      }
3443      if(gMinuit) {
3444          printf(" > Found gMinuit, will re-create it < \n");
3445          delete gMinuit;
3446          gMinuit = new TMinuit;
3447      }
3448      AliUnfolding::fgCorrelationMatrix = 0;
3449      AliUnfolding::fgCorrelationMatrixSquared = 0;
3450      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
3451      AliUnfolding::fgCurrentESDVector = 0;
3452      AliUnfolding::fgEntropyAPriori = 0;
3453      AliUnfolding::fgEfficiency = 0;
3454      AliUnfolding::fgUnfoldedAxis = 0;
3455      AliUnfolding::fgMeasuredAxis = 0;
3456      AliUnfolding::fgFitFunction = 0;
3457      AliUnfolding::fgMaxInput  = -1;
3458      AliUnfolding::fgMaxParams = -1;
3459      AliUnfolding::fgOverflowBinLimit = -1;
3460      AliUnfolding::fgRegularizationWeight = 10000;
3461      AliUnfolding::fgSkipBinsBegin = 0;
3462      AliUnfolding::fgMinuitStepSize = 0.1;
3463      AliUnfolding::fgMinuitPrecision = 1e-6;
3464      AliUnfolding::fgMinuitMaxIterations = 1000000;
3465      AliUnfolding::fgMinuitStrategy = 1.;
3466      AliUnfolding::fgMinimumInitialValue = kFALSE;
3467      AliUnfolding::fgMinimumInitialValueFix = -1;
3468      AliUnfolding::fgNormalizeInput = kFALSE;
3469      AliUnfolding::fgNotFoundEvents = 0;
3470      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
3471      AliUnfolding::fgBayesianSmoothing  = 1;
3472      AliUnfolding::fgBayesianIterations = 10;
3473      AliUnfolding::fgDebug = kFALSE;
3474      AliUnfolding::fgCallCount = 0;
3475      AliUnfolding::fgPowern = 5;
3476      AliUnfolding::fChi2FromFit = 0.;
3477      AliUnfolding::fPenaltyVal  = 0.;
3478      AliUnfolding::fAvgResidual = 0.;
3479      AliUnfolding::fgPrintChi2Details = 0;
3480      AliUnfolding::fgCanvas = 0;
3481      AliUnfolding::fghUnfolded = 0;
3482      AliUnfolding::fghCorrelation = 0;
3483      AliUnfolding::fghEfficiency = 0;
3484      AliUnfolding::fghMeasured = 0;
3485      AliUnfolding::SetMinuitStepSize(1.);
3486      AliUnfolding::SetMinuitPrecision(1e-6);
3487      AliUnfolding::SetMinuitMaxIterations(100000);
3488      AliUnfolding::SetMinuitStrategy(2.);
3489      AliUnfolding::SetDebug(1);
3490 }
3491 //_____________________________________________________________________________
3492 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
3493     // protect heap by adding unique qualifier to name
3494     if(!protect) return 0x0;
3495     TH1D* p = (TH1D*)protect->Clone();
3496     TString tempString(fActiveString);
3497     tempString+=suffix;
3498     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3499     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3500     if(kill) delete protect;
3501     return p;
3502 }
3503 //_____________________________________________________________________________
3504 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
3505     // protect heap by adding unique qualifier to name
3506     if(!protect) return 0x0;
3507     TH2D* p = (TH2D*)protect->Clone();
3508     TString tempString(fActiveString);
3509     tempString+=suffix;
3510     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3511     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3512     if(kill) delete protect;
3513     return p;
3514 }
3515 //_____________________________________________________________________________
3516 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
3517     // protect heap by adding unique qualifier to name
3518     if(!protect) return 0x0;
3519     TGraphErrors* p = (TGraphErrors*)protect->Clone();
3520     TString tempString(fActiveString);
3521     tempString+=suffix;
3522     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3523     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3524     if(kill) delete protect;
3525     return p;
3526 }
3527 //_____________________________________________________________________________
3528 void AliJetFlowTools::MakeAU() {
3529     // === azimuthal unfolding ===
3530     //
3531     // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
3532     // in transverse momentum and azimuthal correlation space to extract v2 params
3533     // settings are equal to the ones used for 'Make()'
3534     //
3535     // basic steps that are followed:
3536     // 1) rebin the raw output of the jet task to the desired binnings
3537     // 2) calls the unfolding routine
3538     // 3) writes output to file
3539     // can be repeated multiple times with different configurations
3540
3541     Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
3542     Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
3543     TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3544     TH1D* dPtdPhi[8];
3545     for(Int_t i(0); i < 8; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
3546
3547     for(Int_t i(0); i < 8; i++) {
3548         // 1) manipulation of input histograms
3549         // check if the input variables are present
3550         if(!PrepareForUnfolding(low[i], up[i])) return;
3551         // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
3552         //     parts of the spectrum can end up in over or underflow bins
3553         TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
3554
3555         // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
3556         // the template will be used as a prior for the chi2 unfolding
3557         TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
3558
3559         // get the full response matrix from the dpt and the detector response
3560         fDetectorResponse = NormalizeTH2D(fDetectorResponse);
3561         // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
3562         // so that unfolding should return the initial spectrum
3563         if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
3564         else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
3565         else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
3566         else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
3567         // normalize each slide of the response to one
3568         NormalizeTH2D(fFullResponseIn);
3569         // resize to desired binning scheme
3570         TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
3571         // get the kinematic efficiency
3572         TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
3573         kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
3574         // suppress the errors
3575         for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
3576         TH1D* jetFindingEfficiency(0x0);
3577         if(fJetFindingEff) {
3578             jetFindingEfficiency = ProtectHeap(fJetFindingEff);
3579             jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
3580             jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
3581         }
3582         // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
3583         TH1D* unfoldedJetSpectrumIn(0x0);
3584         fActiveDir->cd();                   // select active dir
3585         TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
3586         dirIn->cd();                        // select inplane subdir
3587         // select the unfolding method
3588         unfoldedJetSpectrumIn = UnfoldWrapper(
3589             measuredJetSpectrumIn,
3590             resizedResponseIn,
3591             kinematicEfficiencyIn,
3592             measuredJetSpectrumTrueBinsIn,
3593             TString("dPtdPhiUnfolding"),
3594             jetFindingEfficiency);
3595         if(i==5) {
3596             resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
3597             resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
3598             resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
3599             resizedResponseIn = ProtectHeap(resizedResponseIn);
3600             resizedResponseIn->Write();
3601             kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
3602             kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
3603             kinematicEfficiencyIn->Write();
3604             fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
3605             fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
3606             fDetectorResponse->Write();
3607             // optional histograms
3608             if(fSaveFull) {
3609                 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
3610                 fSpectrumIn->Write();
3611                 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
3612                 fDptInDist->Write();
3613                 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
3614                 fDptIn->Write();
3615                 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
3616                 fFullResponseIn->Write();
3617             }
3618         }
3619         fActiveDir->cd();
3620         fDeltaPtDeltaPhi->Write();
3621         fJetPtDeltaPhi->Write();
3622
3623         TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
3624         Double_t integralError(0);
3625         for(Int_t j(0); j < 6; j++) {
3626             // get the integrated
3627             Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
3628             dPtdPhi[j]->SetBinContent(i+1, integral);
3629             dPtdPhi[j]->SetBinError(i+1, integralError);
3630         }
3631         dud->Write();
3632         // save the current state of the unfolding object
3633         SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
3634     }
3635     TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
3636     TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3637     for(Int_t i(0); i < 6; i++) {
3638         dPtdPhi[i]->Fit(fourier, "VI");
3639         v2->SetBinContent(1+i, fourier->GetParameter(1));
3640         v2->SetBinError(1+i, fourier->GetParError(1));
3641         dPtdPhi[i]->Write();
3642     }
3643     v2->Write();
3644 }
3645 //_____________________________________________________________________________