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     fRMS                (kTRUE),
  77     fSymmRMS            (kTRUE),
  78     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
  79     fSaveFull           (kTRUE),
  80     fActiveString       (""),
  81     fActiveDir          (0x0),
  82     fInputList          (0x0),
  83     fRefreshInput       (kTRUE),
  84     fOutputFileName     ("UnfoldedSpectra.root"),
  85     fOutputFile         (0x0),
  86     fCentralityArray    (0x0),
  87     fCentralityWeights  (0x0),
  88     fDetectorResponse   (0x0),
  89     fJetFindingEff      (0x0),
  90     fBetaIn             (.1),
  91     fBetaOut            (.1),
  92     fBayesianIterIn     (4),
  93     fBayesianIterOut    (4),
  94     fBayesianSmoothIn   (0.),
  95     fBayesianSmoothOut  (0.),
  96     fAvoidRoundingError (kFALSE),
  97     fUnfoldingAlgorithm (kChi2),
  98     fPrior              (kPriorMeasured),
  99     fPriorUser          (0x0),
 100     fBinsTrue           (0x0),
 101     fBinsRec            (0x0),
 102     fBinsTruePrior      (0x0),
 103     fBinsRecPrior       (0x0),
 104     fSVDRegIn           (5),
 105     fSVDRegOut          (5),
 106     fSVDToy             (kTRUE),
 107     fJetRadius          (0.3),
 108     fEventCount         (-1),
 109     fNormalizeSpectra   (kFALSE),
 110     fSmoothenPrior      (kFALSE),
 111     fFitMin             (60.),
 112     fFitMax             (300.),
 113     fFitStart           (75.),
 114     fSmoothenCounts     (kTRUE),
 115     fTestMode           (kFALSE),
 116     fRawInputProvided   (kFALSE),
 117     fEventPlaneRes      (.63),
 118     fUseDetectorResponse(kTRUE),
 119     fUseDptResponse     (kTRUE),
 120     fTrainPower         (kTRUE),
 121     fDphiUnfolding      (kTRUE),
 122     fDphiDptUnfolding   (kFALSE),
 123     fExLJDpt            (kTRUE),
 124     fTitleFontSize      (-999.),
 125     fRMSSpectrumIn      (0x0),
 126     fRMSSpectrumOut     (0x0),
 127     fRMSRatio           (0x0),
 128     fRMSV2              (0x0),
 129     fDeltaPtDeltaPhi    (0x0),
 130     fJetPtDeltaPhi      (0x0),
 131     fSpectrumIn         (0x0),
 132     fSpectrumOut        (0x0),
 133     fDptInDist          (0x0),
 134     fDptOutDist         (0x0),
 135     fDptIn              (0x0),
 136     fDptOut             (0x0),
 137     fFullResponseIn     (0x0),
 138     fFullResponseOut    (0x0) { // class constructor
 139         // create response maker weight function (tuned to PYTHIA spectrum)
 140         fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
 141         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 142 }
 143 //_____________________________________________________________________________
 144 void AliJetFlowTools::Make() {
 145     // core function of the class
 146     if(fDphiDptUnfolding) {
 147         // to extract the yield as function of Dphi, Dpt - experimental
 148         MakeAU();
 149         return;
 150     }
 151     // 1) rebin the raw output of the jet task to the desired binnings
 152     // 2) calls the unfolding routine
 153     // 3) writes output to file
 154     // can be repeated multiple times with different configurations
 155
 156     // 1) manipulation of input histograms
 157     // check if the input variables are present
 158     if(fRefreshInput) {
 159         if(!PrepareForUnfolding()) {
 160             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
 161             return;
 162         }
 163     }
 164     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
 165     //     parts of the spectrum can end up in over or underflow bins
 166     TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
 167     TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec,  TString("resized_out_"), kFALSE);
 168
 169     // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
 170     // the template will be used as a prior for the chi2 unfolding
 171     TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
 172     TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
 173     // get the full response matrix from the dpt and the detector response
 174     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
 175     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
 176     // so that unfolding should return the initial spectrum
 177     if(!fTestMode) {
 178         if(fUseDptResponse && fUseDetectorResponse) {
 179             fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
 180             fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
 181         } else if (fUseDptResponse && !fUseDetectorResponse) {
 182             fFullResponseIn = fDptIn;
 183             fFullResponseOut = fDptOut;
 184         } else if (!fUseDptResponse && fUseDetectorResponse) {
 185             fFullResponseIn = fDetectorResponse;
 186             fFullResponseOut = fDetectorResponse;
 187         } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
 188             printf(" > No response, exiting ! < \n" );
 189             return;
 190         }
 191     } else {
 192         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
 193         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
 194     }
 195     // normalize each slide of the response to one
 196     NormalizeTH2D(fFullResponseIn);
 197     NormalizeTH2D(fFullResponseOut);
 198     // resize to desired binning scheme
 199     TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
 200     TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
 201     // get the kinematic efficiency
 202     TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
 203     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
 204     TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
 205     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
 206     // suppress the errors
 207     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
 208         kinematicEfficiencyIn->SetBinError(1+i, 0.);
 209         kinematicEfficiencyOut->SetBinError(1+i, 0.);
 210     }
 211     TH1D* jetFindingEfficiency(0x0);
 212     if(fJetFindingEff) {
 213         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
 214         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
 215         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
 216     }
 217     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
 218     TH1D* unfoldedJetSpectrumIn(0x0);
 219     TH1D* unfoldedJetSpectrumOut(0x0);
 220     fActiveDir->cd();                   // select active dir
 221     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
 222     dirIn->cd();                        // select inplane subdir
 223     // do the inplane unfolding
 224     unfoldedJetSpectrumIn = UnfoldWrapper(
 225         measuredJetSpectrumIn,
 226         resizedResponseIn,
 227         kinematicEfficiencyIn,
 228         measuredJetSpectrumTrueBinsIn,
 229         TString("in"),
 230         jetFindingEfficiency);
 231     resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
 232     resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
 233     resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
 234     resizedResponseIn = ProtectHeap(resizedResponseIn);
 235     resizedResponseIn->Write();
 236     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
 237     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
 238     kinematicEfficiencyIn->Write();
 239     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 240     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 241     fDetectorResponse->Write();
 242     // optional histograms
 243     if(fSaveFull) {
 244         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
 245         fSpectrumIn->Write();
 246         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
 247         fDptInDist->Write();
 248         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
 249         fDptIn->Write();
 250         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
 251         fFullResponseIn->Write();
 252     }
 253     fActiveDir->cd();
 254     if(fDphiUnfolding) {
 255         TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
 256         dirOut->cd();
 257         // do the out of plane unfolding
 258         unfoldedJetSpectrumOut = UnfoldWrapper(
 259                     measuredJetSpectrumOut,
 260                     resizedResponseOut,
 261                     kinematicEfficiencyOut,
 262                     measuredJetSpectrumTrueBinsOut,
 263                     TString("out"),
 264                     jetFindingEfficiency);
 265         resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
 266         resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
 267         resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
 268         resizedResponseOut = ProtectHeap(resizedResponseOut);
 269         resizedResponseOut->Write();
 270         kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
 271         kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
 272         kinematicEfficiencyOut->Write();
 273         fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 274         fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 275         fDetectorResponse->Write();
 276         if(jetFindingEfficiency) jetFindingEfficiency->Write();
 277         // optional histograms
 278         if(fSaveFull) {
 279             fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
 280             fSpectrumOut->Write();
 281             fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
 282             fDptOutDist->Write();
 283             fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
 284             fDptOut->Write();
 285             fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
 286             fFullResponseOut->Write();
 287         }
 288
 289         // write general output histograms to file
 290         fActiveDir->cd();
 291         if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
 292             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
 293             if(ratio) {
 294                 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
 295                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 296                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 297                 ratio = ProtectHeap(ratio);
 298                 ratio->Write();
 299                 // write histo values to RMS files if both routines converged
 300                 // input values are weighted by their uncertainty
 301                 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
 302                     if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
 303                     if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
 304                     if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
 305                }
 306             }
 307             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
 308             if(v2) {
 309                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 310                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 311                 v2->GetYaxis()->SetTitle("v_{2}");
 312                 v2 = ProtectHeap(v2);
 313                 v2->Write();
 314             }
 315         } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
 316             TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
 317             if(ratio) {
 318                 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
 319                 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 320                 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 321                 ratio = ProtectHeap(ratio);
 322                 ratio->Write();
 323             }
 324             TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
 325              if(v2) {
 326                 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 327                 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
 328                 v2->GetYaxis()->SetTitle("v_{2}");
 329                 v2 = ProtectHeap(v2);
 330                 v2->Write();
 331             }
 332         }
 333     }   // end of if(fDphiUnfolding)
 334     fDeltaPtDeltaPhi->Write();
 335     unfoldedJetSpectrumIn->Sumw2();
 336     ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
 337     unfoldedJetSpectrumIn->Write();
 338     unfoldedJetSpectrumOut->Sumw2();
 339     ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
 340     unfoldedJetSpectrumOut->Write();
 341     fJetPtDeltaPhi->Write();
 342     // save the current state of the unfolding object
 343     SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
 344     TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
 345     TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
 346     unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
 347     unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
 348     unfoldedJetSpectrumInForSub->Write();
 349
 350 }
 351 //_____________________________________________________________________________
 352 TH1D* AliJetFlowTools::UnfoldWrapper(
 353         const TH1D* measuredJetSpectrum,        // truncated raw jets (same binning as pt rec of response)
 354         const TH2D* resizedResponse,            // response matrix
 355         const TH1D* kinematicEfficiency,        // kinematic efficiency
 356         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
 357         const TString suffix,                   // suffix (in or out of plane)
 358         const TH1D* jetFindingEfficiency)       // jet finding efficiency
 359 {
 360     // wrapper function to call specific unfolding routine
 361     TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
 362     // initialize functon pointer
 363     if(fUnfoldingAlgorithm == kChi2)                    myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
 364     else if(fUnfoldingAlgorithm == kBayesian)           myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
 365     else if(fUnfoldingAlgorithm == kBayesianAli)        myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
 366     else if(fUnfoldingAlgorithm == kSVD)                myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
 367     else if(fUnfoldingAlgorithm == kNone) {
 368         TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
 369         clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
 370         return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
 371     }
 372     else return 0x0;
 373     // do the actual unfolding with the selected function
 374     return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
 375 }
 376 //_____________________________________________________________________________
 377 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
 378         const TH1D* measuredJetSpectrum,      // truncated raw jets (same binning as pt rec of response)
 379         const TH2D* resizedResponse,          // response matrix
 380         const TH1D* kinematicEfficiency,      // kinematic efficiency
 381         const TH1D* measuredJetSpectrumTrueBins,        // unfolding template: same binning is pt gen of response
 382         const TString suffix,                 // suffix (in or out of plane)
 383         const TH1D* jetFindingEfficiency)     // jet finding efficiency (optional)
 384 {
 385     // unfold the spectrum using chi2 minimization
 386
 387     // step 0) setup the static members of AliUnfolding
 388     ResetAliUnfolding();                // reset from previous iteration
 389                                         // also deletes and re-creates the global TVirtualFitter
 390     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
 391     if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 392     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 393     if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 394     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 395     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
 396
 397     // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
 398     // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
 399     // denoted by the suffix 'Local'
 400
 401     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
 402     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
 403     // unfolded local will be filled with the result of the unfolding
 404     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
 405
 406     // full response matrix and kinematic efficiency
 407     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
 408     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
 409
 410     // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
 411     TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
 412     // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
 413     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
 414
 415     // step 2) start the unfolding
 416     Int_t status(-1), i(0);
 417     while(status < 0 && i < 100) {
 418         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
 419         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
 420         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
 421         status = AliUnfolding::Unfold(
 422                 resizedResponseLocal,           // response matrix
 423                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
 424                 measuredJetSpectrumLocal,              // measured spectrum
 425                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
 426                 unfoldedLocal);                 // results
 427         // status holds the minuit fit status (where 0 means convergence)
 428         i++;
 429     }
 430     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
 431     if(status == 0 && gMinuit->fISW[1] == 3) {
 432         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
 433         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
 434         if(gMinuit) gMinuit->Command("SET COV");
 435         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
 436         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
 437         pearson->Print();
 438         TH2D *hPearson(new TH2D(*pearson));
 439         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
 440         hPearson = ProtectHeap(hPearson);
 441         hPearson->Write();
 442     } else status = -1;
 443
 444     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
 445     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
 446     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
 447     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
 448     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
 449     if(ratio) {
 450         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
 451         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 452         ratio = ProtectHeap(ratio);
 453         ratio->Write();
 454     }
 455
 456     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
 457     // objects are cloned using 'ProtectHeap()'
 458     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
 459     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 460     measuredJetSpectrumLocal->Write();
 461
 462     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
 463     resizedResponseLocal->Write();
 464
 465     unfoldedLocal = ProtectHeap(unfoldedLocal);
 466     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
 467     unfoldedLocal->Write();
 468
 469     foldedLocal = ProtectHeap(foldedLocal);
 470     foldedLocal->Write();
 471
 472     priorLocal = ProtectHeap(priorLocal);
 473     priorLocal->Write();
 474
 475     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
 476     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));
 477     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
 478     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
 479     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
 480     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
 481     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
 482     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
 483     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
 484     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
 485     fitStatus->Write();
 486     return unfoldedLocal;
 487 }
 488 //_____________________________________________________________________________
 489 TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
 490         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 491         const TH2D* resizedResponse,                   // full response matrix, normalized in slides of pt true
 492         const TH1D* kinematicEfficiency,              // kinematic efficiency
 493         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 494         const TString suffix,                         // suffix (in, out)
 495         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 496 {
 497
 498     TH1D* priorLocal( GetPrior(
 499         measuredJetSpectrum,              // jet pt in pt rec bins
 500         resizedResponse,                  // full response matrix, normalized in slides of pt true
 501         kinematicEfficiency,              // kinematic efficiency
 502         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 503         suffix,                           // suffix (in, out)
 504         jetFindingEfficiency));           // jet finding efficiency (optional)
 505     if(!priorLocal) {
 506         printf(" > couldn't find prior ! < \n");
 507         return 0x0;
 508     } else printf(" 1) retrieved prior \n");
 509
 510     // go back to the 'root' directory of this instance of the SVD unfolding routine
 511     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 512
 513     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 514     // measured jets in pt rec binning
 515     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
 516     // local copie of the response matrix
 517     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
 518     // local copy of response matrix, all true slides normalized to 1
 519     // this response matrix will eventually be used in the re-folding routine
 520     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
 521     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
 522     // kinematic efficiency
 523     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 524     // place holder histos
 525     TH1D *unfoldedLocalSVD(0x0);
 526     TH1D *foldedLocalSVD(0x0);
 527     cout << " 2) setup necessary input " << endl;
 528     // 3) configure routine
 529     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 530     cout << " step 3) configured routine " << endl;
 531
 532     // 4) get transpose matrices
 533     // a) get the transpose of the full response matrix
 534     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
 535     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 536     // normalize it with the prior. this will ensure that high statistics bins will constrain the
 537     // end result most strenuously than bins with limited number of counts
 538     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
 539     cout << " 4) retrieved first transpose matrix " << endl;
 540
 541     // 5) get response for SVD unfolding
 542     RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
 543     cout << " 5) retrieved roo unfold response object " << endl;
 544
 545     // 6) actualy unfolding loop
 546     RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 547     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
 548     // correct the spectrum for the kinematic efficiency
 549     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
 550
 551     // get the pearson coefficients from the covariance matrix
 552     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
 553     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 554     if(pearson) {
 555         TH2D* hPearson(new TH2D(*pearson));
 556         pearson->Print();
 557         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 558         hPearson = ProtectHeap(hPearson);
 559         hPearson->Write();
 560     } else return 0x0;       // return if unfolding didn't converge
 561
 562     // plot singular values and d_i vector
 563     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
 564     TH1* hSVal(svdUnfold->GetSV());
 565     TH1D* hdi(svdUnfold->GetD());
 566     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
 567     hSVal->SetXTitle("singular values");
 568     hSVal->Write();
 569     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
 570     hdi->SetXTitle("|d_{i}^{kreg}|");
 571     hdi->Write();
 572     cout << " plotted singular values and d_i vector " << endl;
 573
 574     // 7) refold the unfolded spectrum
 575     foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
 576     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
 577     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
 578     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 579     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 580     ratio->Write();
 581     cout << " 7) refolded the unfolded spectrum " << endl;
 582
 583     // write the measured, unfolded and re-folded spectra to the output directory
 584     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 585     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 586     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 587     measuredJetSpectrumLocal->Write(); // input spectrum
 588     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 589     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
 590     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
 591     unfoldedLocalSVD->Write();  // unfolded spectrum
 592     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 593     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
 594     foldedLocalSVD->Write();    // re-folded spectrum
 595
 596     // save more general bookkeeeping histograms to the output directory
 597     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
 598     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 599     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 600     responseMatrixLocalTransposePrior->Write();
 601     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
 602     priorLocal->SetXTitle("p_{t} [GeV/c]");
 603     priorLocal = ProtectHeap(priorLocal);
 604     priorLocal->Write();
 605     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
 606     resizedResponseLocalNorm->Write();
 607
 608     // save some info
 609     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));
 610     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 611     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
 612     fitStatus->Write();
 613
 614     return unfoldedLocalSVD;
 615 }
 616 //_____________________________________________________________________________
 617 TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
 618         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 619         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
 620         const TH1D* kinematicEfficiency,              // kinematic efficiency
 621         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 622         const TString suffix,                         // suffix (in, out)
 623         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 624 {
 625     // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
 626     // FIXME careful, not tested yet ! (06122013) FIXME
 627
 628     // step 0) setup the static members of AliUnfolding
 629     ResetAliUnfolding();                // reset from previous iteration
 630                                         // also deletes and re-creates the global TVirtualFitter
 631     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
 632     if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
 633     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
 634     else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
 635     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
 636     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
 637
 638     // 1) get a prior for unfolding and clone all the input histograms
 639     TH1D* priorLocal( GetPrior(
 640         measuredJetSpectrum,              // jet pt in pt rec bins
 641         resizedResponse,                  // full response matrix, normalized in slides of pt true
 642         kinematicEfficiency,              // kinematic efficiency
 643         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 644         suffix,                           // suffix (in, out)
 645         jetFindingEfficiency));           // jet finding efficiency (optional)
 646     if(!priorLocal) {
 647         printf(" > couldn't find prior ! < \n");
 648         return 0x0;
 649     } else printf(" 1) retrieved prior \n");
 650     // switch back to root dir of this unfolding procedure
 651     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 652
 653     // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
 654     TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
 655     // unfolded local will be filled with the result of the unfolding
 656     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
 657
 658     // full response matrix and kinematic efficiency
 659     TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
 660     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
 661
 662     // step 2) start the unfolding
 663     Int_t status(-1), i(0);
 664     while(status < 0 && i < 100) {
 665         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
 666         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
 667         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
 668         status = AliUnfolding::Unfold(
 669                 resizedResponseLocal,           // response matrix
 670                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
 671                 measuredJetSpectrumLocal,              // measured spectrum
 672                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
 673                 unfoldedLocal);                 // results
 674         // status holds the minuit fit status (where 0 means convergence)
 675         i++;
 676     }
 677     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
 678     if(status == 0 && gMinuit->fISW[1] == 3) {
 679         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
 680         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
 681         if(gMinuit) gMinuit->Command("SET COV");
 682         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
 683         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
 684         pearson->Print();
 685         TH2D *hPearson(new TH2D(*pearson));
 686         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
 687         hPearson = ProtectHeap(hPearson);
 688         hPearson->Write();
 689     } else status = -1;
 690
 691     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
 692     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
 693     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
 694     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
 695     TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
 696     if(ratio) {
 697         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
 698         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 699         ratio = ProtectHeap(ratio);
 700         ratio->Write();
 701     }
 702
 703     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
 704     // objects are cloned using 'ProtectHeap()'
 705     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
 706     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 707     measuredJetSpectrumLocal->Write();
 708
 709     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
 710     resizedResponseLocal->Write();
 711
 712     unfoldedLocal = ProtectHeap(unfoldedLocal);
 713     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
 714     unfoldedLocal->Write();
 715
 716     foldedLocal = ProtectHeap(foldedLocal);
 717     foldedLocal->Write();
 718
 719     priorLocal = ProtectHeap(priorLocal);
 720     priorLocal->Write();
 721
 722     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
 723     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));
 724     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
 725     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
 726     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
 727     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
 728     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
 729     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
 730     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
 731     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
 732     fitStatus->Write();
 733     return unfoldedLocal;
 734 }
 735 //_____________________________________________________________________________
 736 TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
 737         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
 738         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
 739         const TH1D* kinematicEfficiency,              // kinematic efficiency
 740         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 741         const TString suffix,                         // suffix (in, out)
 742         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 743 {
 744     // use bayesian unfolding from the RooUnfold package to unfold jet spectra
 745
 746     // 1) get a prior for unfolding.
 747     TH1D* priorLocal( GetPrior(
 748         measuredJetSpectrum,              // jet pt in pt rec bins
 749         resizedResponse,                  // full response matrix, normalized in slides of pt true
 750         kinematicEfficiency,              // kinematic efficiency
 751         measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
 752         suffix,                           // suffix (in, out)
 753         jetFindingEfficiency));           // jet finding efficiency (optional)
 754     if(!priorLocal) {
 755         printf(" > couldn't find prior ! < \n");
 756         return 0x0;
 757     } else printf(" 1) retrieved prior \n");
 758     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 759
 760     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 761     // measured jets in pt rec binning
 762     TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
 763     // local copie of the response matrix
 764     TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
 765     // local copy of response matrix, all true slides normalized to 1
 766     // this response matrix will eventually be used in the re-folding routine
 767     TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
 768     resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
 769     // kinematic efficiency
 770     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 771     // place holder histos
 772     TH1D *unfoldedLocalBayes(0x0);
 773     TH1D *foldedLocalBayes(0x0);
 774     cout << " 2) setup necessary input " << endl;
 775     // 4) get transpose matrices
 776     // a) get the transpose of the full response matrix
 777     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
 778     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 779     // normalize it with the prior. this will ensure that high statistics bins will constrain the
 780     // end result most strenuously than bins with limited number of counts
 781     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
 782     // 3) get response for Bayesian unfolding
 783     RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
 784
 785     // 4) actualy unfolding loop
 786     RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
 787     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 788     unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
 789     // correct the spectrum for the kinematic efficiency
 790     unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
 791     // get the pearson coefficients from the covariance matrix
 792     TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
 793     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 794     if(pearson) {
 795         TH2D* hPearson(new TH2D(*pearson));
 796         pearson->Print();
 797         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 798         hPearson = ProtectHeap(hPearson);
 799         hPearson->Write();
 800     } else return 0x0;       // return if unfolding didn't converge
 801
 802     // 5) refold the unfolded spectrum
 803     foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
 804     TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio  measured / re-folded", kTRUE));
 805     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
 806     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 807     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 808     ratio->Write();
 809     cout << " 7) refolded the unfolded spectrum " << endl;
 810
 811     // write the measured, unfolded and re-folded spectra to the output directory
 812     measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 813     measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
 814     measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 815     measuredJetSpectrumLocal->Write(); // input spectrum
 816     unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 817     unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
 818     if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
 819     unfoldedLocalBayes->Write();  // unfolded spectrum
 820     foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 821     foldedLocalBayes = ProtectHeap(foldedLocalBayes);
 822     foldedLocalBayes->Write();    // re-folded spectrum
 823
 824     // save more general bookkeeeping histograms to the output directory
 825     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
 826     responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
 827     responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
 828     responseMatrixLocalTransposePrior->Write();
 829     priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
 830     priorLocal->SetXTitle("p_{t} [GeV/c]");
 831     priorLocal = ProtectHeap(priorLocal);
 832     priorLocal->Write();
 833     resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
 834     resizedResponseLocalNorm->Write();
 835
 836     // save some info
 837     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));
 838     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
 839     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
 840     fitStatus->Write();
 841
 842     return  unfoldedLocalBayes;
 843 }
 844 //_____________________________________________________________________________
 845 Bool_t AliJetFlowTools::PrepareForUnfolding()
 846 {
 847     // prepare for unfolding
 848     if(fRawInputProvided) return kTRUE;
 849     if(!fInputList) {
 850         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
 851         return kFALSE;
 852     }
 853     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
 854     // check if the pt bin for true and rec have been set
 855     if(!fBinsTrue || !fBinsRec) {
 856         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
 857         return kFALSE;
 858     }
 859     if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
 860                           // procedures, these profiles will be nonsensical, user is responsible
 861         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 862         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 863         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 864     }
 865     if(!fTrainPower) {
 866         // clear minuit state to avoid constraining the fit with the results of the previous iteration
 867         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 868     }
 869     // extract the spectra
 870     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
 871     if(!fInputList->FindObject(spectrumName.Data())) {
 872         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
 873         return kFALSE;
 874     }
 875
 876     // get the first scaled spectrum
 877     fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
 878     // clone the spectrum on the heap. this is necessary since scale or add change the
 879     // contents of the original histogram
 880     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
 881     fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
 882     printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
 883     // merge subsequent bins (if any)
 884     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 885         spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
 886         printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
 887         fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
 888     }
 889
 890     // in plane spectrum
 891     if(!fDphiUnfolding) {
 892         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
 893         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
 894     } else {
 895         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
 896         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
 897         fSpectrumIn = ProtectHeap(fSpectrumIn);
 898         // out of plane spectrum
 899         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
 900         fSpectrumOut = ProtectHeap(fSpectrumOut);
 901     }
 902     // normalize spectra to event count if requested
 903     if(fNormalizeSpectra) {
 904         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
 905         rho->Scale(fCentralityWeights->At(0));
 906         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 907            rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
 908         }
 909         if(!rho) return 0x0;
 910         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
 911         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
 912         if(fEventCount > 0) {
 913             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
 914             fSpectrumOut->Sumw2();
 915             Double_t pt(0);
 916             Double_t error(0); // lots of issues with the errors here ...
 917             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
 918                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
 919                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
 920                 fSpectrumIn->SetBinContent(1+i, pt);
 921                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
 922                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
 923                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
 924             }
 925             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
 926                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
 927                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
 928                 fSpectrumOut->SetBinContent(1+i, pt);
 929                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
 930                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
 931                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
 932             }
 933         }
 934         if(normalizeToFullSpectrum) fEventCount = -1;
 935     }
 936     // extract the delta pt matrices
 937     TString deltaptName("");
 938     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
 939     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
 940     if(!fDeltaPtDeltaPhi) {
 941         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
 942         printf(" > may be ok, depending no what you want to do < \n");
 943         fRefreshInput = kTRUE;
 944         return kTRUE;
 945     }
 946
 947     // clone the distribution on the heap and if requested merge with other centralities
 948     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
 949     fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
 950     printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
 951     for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
 952         deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
 953         printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
 954         fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
 955     }
 956
 957     // in plane delta pt distribution
 958     if(!fDphiUnfolding) {
 959         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
 960         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
 961     } else {
 962         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
 963         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
 964         // out of plane delta pt distribution
 965         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
 966         fDptInDist = ProtectHeap(fDptInDist);
 967         fDptOutDist = ProtectHeap(fDptOutDist);
 968         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
 969     }
 970
 971     // create a rec - true smeared response matrix
 972     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
 973     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 974         Bool_t skip = kFALSE;
 975         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 976             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
 977             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 978         }
 979     }
 980     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
 981     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
 982         Bool_t skip = kFALSE;
 983         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
 984             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
 985             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
 986         }
 987     }
 988     fDptIn = new TH2D(*rfIn);
 989     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
 990     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 991     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 992     fDptIn = ProtectHeap(fDptIn);
 993     fDptOut = new TH2D(*rfOut);
 994     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
 995     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
 996     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
 997     fDptOut = ProtectHeap(fDptOut);
 998
 999     fRefreshInput = kTRUE;     // force cloning of the input
1000     return kTRUE;
1001 }
1002 //_____________________________________________________________________________
1003 Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
1004     // prepare for unfoldingUA - more robust method to extract input spectra from file
1005     // will replace PrepareForUnfolding eventually (09012014)
1006     if(!fInputList) {
1007         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1008         return kFALSE;
1009     }
1010     if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1011     // check if the pt bin for true and rec have been set
1012     if(!fBinsTrue || !fBinsRec) {
1013         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1014         return kFALSE;
1015     }
1016     if(!fTrainPower) {
1017         // clear minuit state to avoid constraining the fit with the results of the previous iteration
1018         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1019     }
1020     // extract the spectra
1021     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
1022     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1023     if(!fJetPtDeltaPhi) {
1024         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1025         return kFALSE;
1026     }
1027     if(fCentralityArray) {
1028         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1029             spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
1030             fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1031         }
1032     }
1033     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1034     // in plane spectrum
1035     fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1036     // extract the delta pt matrices
1037     TString deltaptName("");
1038     deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1039     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1040     if(!fDeltaPtDeltaPhi) {
1041         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1042         printf(" > may be ok, depending no what you want to do < \n");
1043         fRefreshInput = kTRUE;
1044         return kTRUE;
1045     }
1046     if(fCentralityArray) {
1047         for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1048             deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1049             fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1050         }
1051     }
1052
1053     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1054     // in plane delta pt distribution
1055     fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1056     // create a rec - true smeared response matrix
1057     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1058     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1059         Bool_t skip = kFALSE;
1060         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1061             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1062             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1063         }
1064     }
1065     fDptIn = new TH2D(*rfIn);
1066     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1067     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1068     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1069     fDptIn = ProtectHeap(fDptIn);
1070
1071     return kTRUE;
1072 }
1073 //_____________________________________________________________________________
1074 TH1D* AliJetFlowTools::GetPrior(
1075         const TH1D* measuredJetSpectrum,              // jet pt in pt rec bins
1076         const TH2D* resizedResponse,                  // full response matrix, normalized in slides of pt true
1077         const TH1D* kinematicEfficiency,              // kinematic efficiency
1078         const TH1D* measuredJetSpectrumTrueBins,      // jet pt in pt true bins, also the prior when measured is chosen as prior
1079         const TString suffix,                         // suffix (in, out)
1080         const TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
1081 {
1082     // 1) get a prior for unfolding.
1083     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1084     TH1D* unfolded(0x0);
1085     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1086     dirOut->cd();
1087     switch (fPrior) {    // select the prior for unfolding
1088         case kPriorPythia : {
1089             if(!fPriorUser) {
1090                 printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1091                 return 0x0;
1092             }
1093             // rebin the given prior to the true spectrum (creates a new histo)
1094             return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1095         } break;
1096         case kPriorChi2 : {
1097             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
1098                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original binning
1099                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
1100                 TArrayD* tempArrayRec(fBinsRec);
1101                 fBinsRec = fBinsRecPrior;
1102                 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1103                 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1104                 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1105                 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1106                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1107                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1108                 unfolded= UnfoldSpectrumChi2(
1109                             measuredJetSpectrumChi2,
1110                             resizedResponseChi2,
1111                             kinematicEfficiencyChi2,
1112                             measuredJetSpectrumTrueBinsChi2,    // prior for chi2 unfolding (measured)
1113                             TString(Form("prior_%s", suffix.Data())));
1114                if(!unfolded) {
1115                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1116                     printf("               probably Chi2 unfolding did not converge < \n");
1117                     return 0x0;
1118                 }
1119                 fBinsTrue = tempArrayTrue;  // reset bins borders
1120                 fBinsRec = tempArrayRec;
1121                 // if the chi2 prior has a different binning, rebin to the true binning for the  unfolding
1122                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
1123             } else {
1124                 unfolded = UnfoldSpectrumChi2(
1125                             measuredJetSpectrum,
1126                             resizedResponse,
1127                             kinematicEfficiency,
1128                             measuredJetSpectrumTrueBins,        // prior for chi2 unfolding (measured)
1129                             TString(Form("prior_%s", suffix.Data())));
1130                 if(!unfolded) {
1131                     printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1132                     printf("               probably Chi2 unfolding did not converge < \n");
1133                     return 0x0;
1134                 }
1135             }
1136             break;
1137         }
1138         case kPriorMeasured : {
1139             unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
1140         }
1141         default : break;
1142     }
1143     // it can be important that the prior is smooth, this can be achieved by
1144     // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1145     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1146     TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1147     if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1148     return priorLocal;
1149 }
1150 //_____________________________________________________________________________
1151 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1152     // resize the x-axis of a th1d
1153     if(!histo) {
1154         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1155         return NULL;
1156     }
1157     // see how many bins we need to copy
1158     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);
1159     // low is the bin number of the first new bin
1160     Int_t l = histo->GetXaxis()->FindBin(low);
1161     // set the values
1162     Double_t _x(0), _xx(0);
1163     for(Int_t i(0); i < up-low; i++) {
1164         _x = histo->GetBinContent(l+i);
1165         _xx=histo->GetBinError(l+i);
1166         resized->SetBinContent(i+1, _x);
1167         resized->SetBinError(i+1, _xx);
1168     }
1169     return resized;
1170 }
1171 //_____________________________________________________________________________
1172 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1173     // resize the y-axis of a th2d
1174     if(!histo) {
1175         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1176         return NULL;
1177     }
1178     // see how many bins we need to copy
1179     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());
1180     // assume only the y-axis has changed
1181     // low is the bin number of the first new bin
1182     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1183     // set the values
1184     Double_t _x(0), _xx(0);
1185     for(Int_t i(0); i < x->GetSize(); i++) {
1186         for(Int_t j(0); j < y->GetSize(); j++) {
1187             _x = histo->GetBinContent(i, low+j);
1188             _xx=histo->GetBinError(i, low+1+j);
1189             resized->SetBinContent(i, j, _x);
1190             resized->SetBinError(i, j, _xx);
1191         }
1192     }
1193     return resized;
1194 }
1195 //_____________________________________________________________________________
1196 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1197     // general method to normalize all vertical slices of a th2 to unity
1198     // i.e. get a probability matrix
1199     if(!histo) {
1200         printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
1201         return NULL;
1202     }
1203     Int_t binsX = histo->GetXaxis()->GetNbins();
1204     Int_t binsY = histo->GetYaxis()->GetNbins();
1205
1206     // normalize all slices in x
1207     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1208         Double_t weight = 0;
1209         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1210             weight+=histo->GetBinContent(i+1, j+1);
1211         }       // now we know the total weight
1212         for(Int_t j(0); j < binsY; j++) {
1213             if (weight <= 0 ) continue;
1214             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1215             if(noError) histo->SetBinError(  1+i, j+1, 0.);
1216             else histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1217         }
1218     }
1219     return histo;
1220 }
1221 //_____________________________________________________________________________
1222 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1223     // return a TH1D with the supplied histogram rebinned to the supplied bins
1224     // the returned histogram is new, the original is deleted from the heap if kill is true
1225     if(!histo || !bins) {
1226         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1227         return NULL;
1228     }
1229     // create the output histo
1230     TString name = histo->GetName();
1231     name+="_template";
1232     name+=suffix;
1233     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1234     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1235         // loop over the bins of the old histo and fill the new one with its data
1236         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1237     }
1238     if(kill) delete histo;
1239     return rebinned;
1240 }
1241 //_____________________________________________________________________________
1242 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1243     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1244     if(!fResponseMaker || !binsTrue || !binsRec) {
1245         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1246         return 0x0;
1247     }
1248     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1249     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1250 }
1251 //_____________________________________________________________________________
1252 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1253 {
1254     // multiply two matrices
1255     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1256     TH2D* c = (TH2D*)a->Clone("c");
1257     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1258         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1259             Double_t val = 0;
1260             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1261                 Int_t y2 = x1;
1262                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1263             }
1264             c->SetBinContent(x2, y1, val);
1265             c->SetBinError(x2, y1, 0.);
1266         }
1267     }
1268     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1269     return c;
1270 }
1271 //_____________________________________________________________________________
1272 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1273 {
1274     // normalize a th1d to a certain scale
1275     histo->Sumw2();
1276     Double_t integral = histo->Integral()*scale;
1277     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1278     else if (scale != 1.) histo->Scale(1./scale, "width");
1279     else printf(" > Histogram integral < 0, cannot normalize \n");
1280     return histo;
1281 }
1282 //_____________________________________________________________________________
1283 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1284 {
1285     // Calculate pearson coefficients from covariance matrix
1286     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1287     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1288     Double_t pearson(0.);
1289     if(nrows==0 && ncols==0) return 0x0;
1290     for(Int_t row = 0; row < nrows; row++) {
1291         for(Int_t col = 0; col<ncols; col++) {
1292         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1293         (*pearsonCoefficients)(row,col) = pearson;
1294         }
1295     }
1296     return pearsonCoefficients;
1297 }
1298 //_____________________________________________________________________________
1299 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1300     // smoothen the spectrum using a user defined function
1301     // returns a clone of the original spectrum if fitting failed
1302     // if kill is kTRUE the input spectrum will be deleted from the heap
1303     // if 'count' is selected, bins are filled with integers (necessary if the
1304     // histogram is interpreted in a routine which accepts only counts)
1305     if(!spectrum || !function) return 0x0;
1306     if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1307     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1308     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1309     TFitResultPtr r = temp->Fit(function, "", "", min, max);
1310     if((int)r == 0) {   // MINUIT status
1311         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1312             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1313                 (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));
1314                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1315             }
1316         }
1317     }
1318     if(kill) delete spectrum;
1319     return temp;
1320 }
1321 //_____________________________________________________________________________
1322 void AliJetFlowTools::Style(Bool_t legacy)
1323 {
1324     // set global style for your current aliroot session
1325     if(!gStyle) return;
1326     // legacy style is pleasing to the eye, default is the formal ALICE style
1327     if(legacy) {
1328         gStyle->SetCanvasColor(-1);
1329         gStyle->SetPadColor(-1);
1330         gStyle->SetFrameFillColor(-1);
1331         gStyle->SetHistFillColor(-1);
1332         gStyle->SetTitleFillColor(-1);
1333         gStyle->SetFillColor(-1);
1334         gStyle->SetFillStyle(4000);
1335         gStyle->SetStatStyle(0);
1336         gStyle->SetTitleStyle(0);
1337         gStyle->SetCanvasBorderSize(0);
1338         gStyle->SetFrameBorderSize(0);
1339         gStyle->SetLegendBorderSize(0);
1340         gStyle->SetStatBorderSize(0);
1341         gStyle->SetTitleBorderSize(0);
1342     } else {
1343         gStyle->Reset("Plain");
1344         gStyle->SetOptTitle(0);
1345         gStyle->SetOptStat(0);
1346         gStyle->SetPalette(1);
1347         gStyle->SetCanvasColor(10);
1348         gStyle->SetCanvasBorderMode(0);
1349         gStyle->SetFrameLineWidth(1);
1350         gStyle->SetFrameFillColor(kWhite);
1351         gStyle->SetPadColor(10);
1352         gStyle->SetPadTickX(1);
1353         gStyle->SetPadTickY(1);
1354         gStyle->SetPadBottomMargin(0.15);
1355         gStyle->SetPadLeftMargin(0.15);
1356         gStyle->SetHistLineWidth(1);
1357         gStyle->SetHistLineColor(kRed);
1358         gStyle->SetFuncWidth(2);
1359         gStyle->SetFuncColor(kGreen);
1360         gStyle->SetLineWidth(2);
1361         gStyle->SetLabelSize(0.045,"xyz");
1362         gStyle->SetLabelOffset(0.01,"y");
1363         gStyle->SetLabelOffset(0.01,"x");
1364         gStyle->SetLabelColor(kBlack,"xyz");
1365         gStyle->SetTitleSize(0.05,"xyz");
1366         gStyle->SetTitleOffset(1.25,"y");
1367         gStyle->SetTitleOffset(1.2,"x");
1368         gStyle->SetTitleFillColor(kWhite);
1369         gStyle->SetTextSizePixels(26);
1370         gStyle->SetTextFont(42);
1371         gStyle->SetLegendBorderSize(0);
1372         gStyle->SetLegendFillColor(kWhite);
1373         gStyle->SetLegendFont(42);
1374     }
1375 }
1376 //_____________________________________________________________________________
1377 void AliJetFlowTools::Style(TCanvas* c, TString style)
1378 {
1379     // set a default style for a canvas
1380     if(!strcmp(style.Data(), "PEARSON")) {
1381         printf(" > style PEARSON canvas < \n");
1382         gStyle->SetOptStat(0);
1383         c->SetGridx();
1384         c->SetGridy();
1385         c->SetTicks();
1386         return;
1387     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1388         printf(" > style SPECTRUM canvas < \n");
1389         gStyle->SetOptStat(0);
1390         c->SetLogy();
1391         c->SetGridx();
1392         c->SetGridy();
1393         c->SetTicks();
1394         return;
1395     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1396 }
1397 //_____________________________________________________________________________
1398 void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
1399 {
1400     // set a default style for a canva
1401
1402     if(legacy) {
1403         c->SetLeftMargin(.25);
1404         c->SetBottomMargin(.25);
1405     }
1406     else Style();
1407     if(!strcmp(style.Data(), "PEARSON")) {
1408         printf(" > style PEARSON pad < \n");
1409         gStyle->SetOptStat(0);
1410         c->SetGridx();
1411         c->SetGridy();
1412         c->SetTicks();
1413         return;
1414     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1415         printf(" > style SPECTRUM pad < \n");
1416         gStyle->SetOptStat(0);
1417         c->SetLogy();
1418         c->SetGridx();
1419         c->SetGridy();
1420         c->SetTicks();
1421         return;
1422     } else if (!strcmp(style.Data(), "GRID")) {
1423         printf(" > style GRID pad < \n");
1424         gStyle->SetOptStat(0);
1425         c->SetGridx();
1426         c->SetGridy();
1427         c->SetTicks();
1428     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1429 }
1430 //_____________________________________________________________________________
1431 void AliJetFlowTools::Style(TLegend* l)
1432 {
1433     // set a default style for a legend
1434     l->SetFillColor(0);
1435     l->SetBorderSize(0);
1436     if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
1437 }
1438 //_____________________________________________________________________________
1439 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
1440 {
1441     // style a histo
1442     h->SetLineColor(col);
1443     h->SetMarkerColor(col);
1444     h->SetLineWidth(2);
1445     h->SetMarkerSize(1);
1446     if(legacy) {
1447         h->SetTitle("");
1448         h->GetYaxis()->SetLabelSize(0.05);
1449         h->GetXaxis()->SetLabelSize(0.05);
1450         h->GetYaxis()->SetTitleOffset(1.5);
1451         h->GetXaxis()->SetTitleOffset(1.5);
1452         h->GetYaxis()->SetTitleSize(.05);
1453         h->GetXaxis()->SetTitleSize(.05);
1454     } else Style();
1455     switch (type) {
1456         case kInPlaneSpectrum : {
1457             h->SetTitle("IN PLANE");
1458             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1459             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1460         } break;
1461         case kOutPlaneSpectrum : {
1462             h->SetTitle("OUT OF PLANE");
1463             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1464             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1465        } break;
1466        case kUnfoldedSpectrum : {
1467             h->SetTitle("UNFOLDED");
1468             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1469             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1470        } break;
1471        case kFoldedSpectrum : {
1472             h->SetTitle("FOLDED");
1473             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1474             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1475        } break;
1476        case kMeasuredSpectrum : {
1477             h->SetTitle("MEASURED");
1478             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1479             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1480        } break;
1481        case kBar : {
1482             h->SetFillColor(col);
1483             h->SetBarWidth(.6);
1484             h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1485             h->SetBarOffset(0.2);
1486        }
1487        default : break;
1488     }
1489 }
1490 //_____________________________________________________________________________
1491 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
1492 {
1493     // style a tgraph
1494     h->SetLineColor(col);
1495     h->SetMarkerColor(col);
1496     h->SetLineWidth(2);
1497     h->SetMarkerSize(1);
1498     h->SetTitle("");
1499     h->SetFillColor(kCyan);
1500     if(legacy) {
1501         h->GetYaxis()->SetLabelSize(0.05);
1502         h->GetXaxis()->SetLabelSize(0.05);
1503         h->GetYaxis()->SetTitleOffset(1.6);
1504         h->GetXaxis()->SetTitleOffset(1.6);
1505         h->GetYaxis()->SetTitleSize(.05);
1506         h->GetXaxis()->SetTitleSize(.05);
1507     } else Style();
1508     h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1509     switch (type) {
1510         case kInPlaneSpectrum : {
1511             h->SetTitle("IN PLANE");
1512             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1513         } break;
1514         case kOutPlaneSpectrum : {
1515             h->SetTitle("OUT OF PLANE");
1516             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1517        } break;
1518        case kUnfoldedSpectrum : {
1519             h->SetTitle("UNFOLDED");
1520             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1521        } break;
1522        case kFoldedSpectrum : {
1523             h->SetTitle("FOLDED");
1524             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1525        } break;
1526        case kMeasuredSpectrum : {
1527             h->SetTitle("MEASURED");
1528             h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1529        } break;
1530        case kRatio : {
1531 //            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}}");
1532             h->GetYaxis()->SetTitle("(d#it{N}^{ch, jet}_{in plane}/(d#it{p}_{T}d#eta))/(d#it{N}^{ch,jet}_{out of plane}/(d#it{p}_{T}d#eta))");
1533        } break;
1534        case kV2 : {
1535 //            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}}}");
1536             h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet}");
1537             h->GetYaxis()->SetRangeUser(-.5, 1.);
1538        } break;
1539        default : break;
1540     }
1541 }
1542 //_____________________________________________________________________________
1543 void AliJetFlowTools::GetNominalValues(
1544         TH1D*& ratio,           // pointer reference, output of this function
1545         TGraphErrors*& v2,      // pointer reference, as output of this function
1546         TArrayI* in,
1547         TArrayI* out,
1548         TString inFile,
1549         TString outFile) const
1550 {
1551     // pass clones of the nominal points and nominal v2 values
1552     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1553     TFile* readMe(new TFile(inFile.Data(), "READ"));                    // open unfolding output read-only
1554     if(readMe->IsZombie()) {
1555         printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1556         return;
1557     }
1558     printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1559     readMe->ls();
1560     TFile* output(new TFile(outFile.Data(), "RECREATE"));   // create new output file
1561     // get some placeholders, have to be initialized but will be deleted
1562     ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1563     TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1564     TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1565     Int_t iIn[] = {in->At(0), in->At(0)};          // first value is the nominal point
1566     Int_t iOut[] = {out->At(0), out->At(0)};
1567
1568     // call the functions and set the relevant pointer references
1569     TH1D* dud(0x0);
1570     DoIntermediateSystematics(
1571             new TArrayI(2, iIn),
1572             new TArrayI(2, iOut),
1573             dud, dud, dud, dud, dud, dud,
1574             ratio,              // pointer reference, output of this function
1575             nominalIn,
1576             nominalOut,
1577             1,
1578             fBinsTrue->At(0),
1579             fBinsTrue->At(fBinsTrue->GetSize()-1),
1580             readMe,
1581             "nominal_values");
1582     v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1583
1584     // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1585     ratio->SetDirectory(0);     // disassociate from current gDirectory
1586     readMe->Close();
1587 }
1588 //_____________________________________________________________________________
1589 void AliJetFlowTools::GetCorrelatedUncertainty(
1590         TGraphAsymmErrors*& corrRatio,  // correlated uncertainty function pointer
1591         TGraphAsymmErrors*& corrV2,     // correlated uncertainty function pointer
1592         TArrayI* variationsIn,          // variantions in plnae
1593         TArrayI* variationsOut,         // variantions out of plane
1594         Bool_t sym,                     // treat as symmmetric
1595         TArrayI* variations2ndIn,       // second source of variations
1596         TArrayI* variations2ndOut,      // second source of variations
1597         Bool_t sym2nd,                  // treat as symmetric
1598         TString type,                   // type of uncertaitny
1599         TString type2,
1600         Int_t columns,                  // divide the output canvasses in this many columns
1601         Float_t rangeLow,               // lower pt range
1602         Float_t rangeUp,                // upper pt range
1603         Float_t corr,                   // correlation strength
1604         TString in,                     // input file name (created by this unfolding class)
1605         TString out                     // output file name (which will hold results of the systematic test)
1606         ) const
1607 {
1608     // do full systematics
1609     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1610     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1611     if(readMe->IsZombie()) {
1612         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1613         return;
1614     }
1615     printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1616     readMe->ls();
1617     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1618
1619     // create some null placeholder pointers
1620     TH1D* relativeErrorVariationInLow(0x0);
1621     TH1D* relativeErrorVariationInUp(0x0);
1622     TH1D* relativeErrorVariationOutLow(0x0);
1623     TH1D* relativeErrorVariationOutUp(0x0);
1624     TH1D* relativeError2ndVariationInLow(0x0);
1625     TH1D* relativeError2ndVariationInUp(0x0);
1626     TH1D* relativeError2ndVariationOutLow(0x0);
1627     TH1D* relativeError2ndVariationOutUp(0x0);
1628     TH1D* relativeStatisticalErrorIn(0x0);
1629     TH1D* relativeStatisticalErrorOut(0x0);
1630     // histo for the nominal ratio in / out
1631     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1632     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1633     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1634
1635     // call the functions
1636     if(variationsIn && variationsOut) {
1637         DoIntermediateSystematics(
1638                 variationsIn,
1639                 variationsOut,
1640                 relativeErrorVariationInUp,        // pointer reference
1641                 relativeErrorVariationInLow,       // pointer reference
1642                 relativeErrorVariationOutUp,       // pointer reference
1643                 relativeErrorVariationOutLow,      // pointer reference
1644                 relativeStatisticalErrorIn,        // pointer reference
1645                 relativeStatisticalErrorOut,       // pointer reference
1646                 nominal,
1647                 nominalIn,
1648                 nominalOut,
1649                 columns,
1650                 rangeLow,
1651                 rangeUp,
1652                 readMe,
1653                 type);
1654         if(relativeErrorVariationInUp) {
1655             // canvas with the error from variation strength
1656             TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1657             relativeErrorVariation->Divide(2);
1658             relativeErrorVariation->cd(1);
1659             Style(gPad, "GRID");
1660             relativeErrorVariationInUp->DrawCopy("b");
1661             relativeErrorVariationInLow->DrawCopy("same b");
1662             Style(AddLegend(gPad));
1663             relativeErrorVariation->cd(2);
1664             Style(gPad, "GRID");
1665             relativeErrorVariationOutUp->DrawCopy("b");
1666             relativeErrorVariationOutLow->DrawCopy("same b");
1667             SavePadToPDF(relativeErrorVariation);
1668             Style(AddLegend(gPad));
1669             relativeErrorVariation->Write();
1670         }
1671     }
1672     // call the functions for a second set of systematic sources
1673     if(variations2ndIn && variations2ndOut) {
1674         DoIntermediateSystematics(
1675                 variations2ndIn,
1676                 variations2ndOut,
1677                 relativeError2ndVariationInUp,        // pointer reference
1678                 relativeError2ndVariationInLow,       // pointer reference
1679                 relativeError2ndVariationOutUp,       // pointer reference
1680                 relativeError2ndVariationOutLow,      // pointer reference
1681                 relativeStatisticalErrorIn,           // pointer reference
1682                 relativeStatisticalErrorOut,          // pointer reference
1683                 nominal,
1684                 nominalIn,
1685                 nominalOut,
1686                 columns,
1687                 rangeLow,
1688                 rangeUp,
1689                 readMe,
1690                 type2);
1691         if(relativeError2ndVariationInUp) {
1692             // canvas with the error from variation strength
1693             TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
1694             relativeError2ndVariation->Divide(2);
1695             relativeError2ndVariation->cd(1);
1696             Style(gPad, "GRID");
1697             relativeError2ndVariationInUp->DrawCopy("b");
1698             relativeError2ndVariationInLow->DrawCopy("same b");
1699             Style(AddLegend(gPad));
1700             relativeError2ndVariation->cd(2);
1701             Style(gPad, "GRID");
1702             relativeError2ndVariationOutUp->DrawCopy("b");
1703             relativeError2ndVariationOutLow->DrawCopy("same b");
1704             SavePadToPDF(relativeError2ndVariation);
1705             Style(AddLegend(gPad));
1706             relativeError2ndVariation->Write();
1707         }
1708
1709     }
1710
1711     // and the placeholder for the final systematic
1712     TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1713     TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1714     TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1715     TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1716     relativeErrorInUp->SetYTitle("relative uncertainty");
1717     relativeErrorOutUp->SetYTitle("relative uncertainty");
1718     relativeErrorInLow->SetYTitle("relative uncertainty");
1719     relativeErrorOutLow->SetYTitle("relative uncertainty");
1720
1721     // merge the correlated errors (FIXME) trivial for one set
1722     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1723     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1724     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1725     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1726
1727     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1728         // for the upper bound
1729         if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
1730         if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
1731         if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1732         if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
1733         dInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1734         // for a symmetric first variation
1735         if(sym) dInUp += aInLow*aInLow;
1736         if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1737         dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1738         if(sym) dOutUp += aOutLow*aOutLow;
1739         if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1740         // for the lower bound
1741         if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
1742         if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
1743         if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1744         if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
1745         dInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1746         if(sym) dInLow += aInUp*aInUp;
1747         if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
1748         dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1749         if(sym) dOutLow += aOutUp*aOutUp;
1750         if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1751     }
1752     // project the estimated errors on the nominal ratio
1753     if(nominal) {
1754         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1755         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1756         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1757         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1758         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1759         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1760         Double_t lowErr(0.), upErr(0.);
1761         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1762             // add the in and out of plane errors in quadrature
1763             // propagation is tricky because we should consider two cases:
1764             // [1] the error is uncorrelated, and we propagate upper 1 with lower 2 and vice versa
1765             // [2] the error is correlated - in this case upper 1 should be propagated with upper 2
1766             // as the fluctuations are bound to always to of same sign
1767             if(corr <= 0) {
1768                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
1769                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1770             } else {
1771                 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) -2.*corr*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1772                 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 2.*corr*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1773             }
1774             ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
1775             ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
1776             // get the bin width (which is the 'error' on x
1777             Double_t binWidth(nominal->GetBinWidth(i+1));
1778             axl[i] = binWidth/2.;
1779             axh[i] = binWidth/2.;
1780             // now get the coordinate for the point
1781             ax[i] = nominal->GetBinCenter(i+1);
1782             ay[i] = nominal->GetBinContent(i+1);
1783         }
1784         // save the nominal ratio
1785         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
1786         corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
1787         nominalError->SetName("correlated uncertainty");
1788         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
1789         nominalCanvas->Divide(2);
1790         nominalCanvas->cd(1);
1791         Style(nominal, kBlack);
1792         Style(nominalError, kCyan, kRatio);
1793         nominalError->SetLineColor(kCyan);
1794         nominalError->SetMarkerColor(kCyan);
1795         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1796         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
1797         nominalError->DrawClone("a2");
1798         nominal->DrawCopy("same E1");
1799         Style(AddLegend(gPad));
1800         Style(gPad, "GRID");
1801         Style(nominalCanvas);
1802         // save nominal v2 and systematic errors
1803         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
1804                     nominalIn,
1805                     nominalOut,
1806                     fEventPlaneRes,
1807                     "v_{2} with shape uncertainty",
1808                     relativeErrorInUp,
1809                     relativeErrorInLow,
1810                     relativeErrorOutUp,
1811                     relativeErrorOutLow,
1812                     corr));
1813         // pass the nominal values to the pointer references
1814         corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
1815         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
1816         nominalCanvas->cd(2);
1817         Style(nominalV2, kBlack);
1818         Style(nominalV2Error, kCyan, kV2);
1819         nominalV2Error->SetLineColor(kCyan);
1820         nominalV2Error->SetMarkerColor(kCyan);
1821         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1822         nominalV2Error->DrawClone("a2");
1823         nominalV2->DrawClone("same E1");
1824         Style(AddLegend(gPad));
1825         Style(gPad, "GRID");
1826         Style(nominalCanvas);
1827         SavePadToPDF(nominalCanvas);
1828         nominalCanvas->Write();
1829     }
1830
1831     TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
1832     relativeError->Divide(2);
1833     relativeError->cd(1);
1834     Style(gPad, "GRID");
1835     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1836     Style(relativeErrorInUp, kBlue, kBar);
1837     Style(relativeErrorInLow, kGreen, kBar);
1838     relativeErrorInUp->DrawCopy("b");
1839     relativeErrorInLow->DrawCopy("same b");
1840     Style(relativeStatisticalErrorIn, kRed);
1841     relativeStatisticalErrorIn->DrawCopy("same");
1842     Style(AddLegend(gPad));
1843     relativeError->cd(2);
1844     Style(gPad, "GRID");
1845     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1846     Style(relativeErrorOutUp, kBlue, kBar);
1847     Style(relativeErrorOutLow, kGreen, kBar);
1848     relativeErrorOutUp->DrawCopy("b");
1849     relativeErrorOutLow->DrawCopy("same b");
1850     Style(relativeStatisticalErrorOut, kRed);
1851     relativeStatisticalErrorOut->DrawCopy("same");
1852     Style(AddLegend(gPad));
1853
1854     // write the buffered file to disk and close the file
1855     SavePadToPDF(relativeError);
1856     relativeError->Write();
1857     output->Write();
1858 //    output->Close();
1859 }
1860 //_____________________________________________________________________________
1861 void AliJetFlowTools::GetShapeUncertainty(
1862         TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
1863         TGraphAsymmErrors*& shapeV2,    // pointer reference to final shape uncertainty of v2
1864         TArrayI* regularizationIn,      // regularization strength systematics, in plane
1865         TArrayI* regularizationOut,     // regularization strenght systematics, out of plane
1866         TArrayI* trueBinIn,             // true bin ranges, in plane
1867         TArrayI* trueBinOut,            // true bin ranges, out of plane
1868         TArrayI* recBinIn,              // rec bin ranges, in plane
1869         TArrayI* recBinOut,             // rec bin ranges, out of plane
1870         TArrayI* methodIn,              // method in
1871         TArrayI* methodOut,             // method out
1872         Int_t columns,                  // divide the output canvasses in this many columns
1873         Float_t rangeLow,               // lower pt range
1874         Float_t rangeUp,                // upper pt range
1875         TString in,                     // input file name (created by this unfolding class)
1876         TString out                     // output file name (which will hold results of the systematic test)
1877         ) const
1878 {
1879     // do full systematics
1880     if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();   // if for some weird reason the unfolding output is still mutable
1881     TFile* readMe(new TFile(in.Data(), "READ"));                        // open unfolding output read-only
1882     if(readMe->IsZombie()) {
1883         printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1884         return;
1885     }
1886     printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1887     readMe->ls();
1888     TFile* output(new TFile(out.Data(), "RECREATE"));   // create new output file
1889
1890     // create some null placeholder pointers
1891     TH1D* relativeErrorRegularizationInLow(0x0);
1892     TH1D* relativeErrorRegularizationInUp(0x0);
1893     TH1D* relativeErrorTrueBinInLow(0x0);
1894     TH1D* relativeErrorTrueBinInUp(0x0);
1895     TH1D* relativeErrorRecBinInLow(0x0);
1896     TH1D* relativeErrorRecBinInUp(0x0);
1897     TH1D* relativeErrorMethodInLow(0x0);
1898     TH1D* relativeErrorMethodInUp(0x0);
1899     TH1D* relativeErrorRegularizationOutLow(0x0);
1900     TH1D* relativeErrorRegularizationOutUp(0x0);
1901     TH1D* relativeErrorTrueBinOutLow(0x0);
1902     TH1D* relativeErrorTrueBinOutUp(0x0);
1903     TH1D* relativeErrorRecBinOutLow(0x0);
1904     TH1D* relativeErrorRecBinOutUp(0x0);
1905     TH1D* relativeStatisticalErrorIn(0x0);
1906     TH1D* relativeStatisticalErrorOut(0x0);
1907     TH1D* relativeErrorMethodOutLow(0x0);
1908     TH1D* relativeErrorMethodOutUp(0x0);
1909     // histo for the nominal ratio in / out
1910     TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1911     TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1912     TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1913
1914     // call the functions
1915     if(regularizationIn && regularizationOut) {
1916         DoIntermediateSystematics(
1917                 regularizationIn,
1918                 regularizationOut,
1919                 relativeErrorRegularizationInUp,        // pointer reference
1920                 relativeErrorRegularizationInLow,       // pointer reference
1921                 relativeErrorRegularizationOutUp,       // pointer reference
1922                 relativeErrorRegularizationOutLow,      // pointer reference
1923                 relativeStatisticalErrorIn,             // pointer reference
1924                 relativeStatisticalErrorOut,            // pointer reference
1925                 nominal,
1926                 nominalIn,
1927                 nominalOut,
1928                 columns,
1929                 rangeLow,
1930                 rangeUp,
1931                 readMe,
1932                 "regularization",
1933                 fRMS);
1934         if(relativeErrorRegularizationInUp) {
1935             // canvas with the error from regularization strength
1936             TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
1937             relativeErrorRegularization->Divide(2);
1938             relativeErrorRegularization->cd(1);
1939             Style(gPad, "GRID");
1940             relativeErrorRegularizationInUp->DrawCopy("b");
1941             relativeErrorRegularizationInLow->DrawCopy("same b");
1942             Style(AddLegend(gPad));
1943             relativeErrorRegularization->cd(2);
1944             Style(gPad, "GRID");
1945             relativeErrorRegularizationOutUp->DrawCopy("b");
1946             relativeErrorRegularizationOutLow->DrawCopy("same b");
1947             SavePadToPDF(relativeErrorRegularization);
1948             Style(AddLegend(gPad));
1949             relativeErrorRegularization->Write();
1950         }
1951     }
1952     if(trueBinIn && trueBinOut) {
1953         DoIntermediateSystematics(
1954                 trueBinIn,
1955                 trueBinOut,
1956                 relativeErrorTrueBinInUp,       // pointer reference
1957                 relativeErrorTrueBinInLow,      // pointer reference
1958                 relativeErrorTrueBinOutUp,      // pointer reference
1959                 relativeErrorTrueBinOutLow,     // pointer reference
1960                 relativeStatisticalErrorIn,
1961                 relativeStatisticalErrorOut,
1962                 nominal,
1963                 nominalIn,
1964                 nominalOut,
1965                 columns,
1966                 rangeLow,
1967                 rangeUp,
1968                 readMe,
1969                 "trueBin");
1970         if(relativeErrorTrueBinInUp) {
1971             TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
1972             relativeErrorTrueBin->Divide(2);
1973             relativeErrorTrueBin->cd(1);
1974             Style(gPad, "GRID");
1975             relativeErrorTrueBinInUp->DrawCopy("b");
1976             relativeErrorTrueBinInLow->DrawCopy("same b");
1977             Style(AddLegend(gPad));
1978             relativeErrorTrueBin->cd(2);
1979             Style(gPad, "GRID");
1980             relativeErrorTrueBinOutUp->DrawCopy("b");
1981             relativeErrorTrueBinOutLow->DrawCopy("same b");
1982             SavePadToPDF(relativeErrorTrueBin);
1983             Style(AddLegend(gPad));
1984             relativeErrorTrueBin->Write();
1985         }
1986     }
1987     if(recBinIn && recBinOut) {
1988         DoIntermediateSystematics(
1989                 recBinIn,
1990                 recBinOut,
1991                 relativeErrorRecBinInUp,       // pointer reference
1992                 relativeErrorRecBinInLow,        // pointer reference
1993                 relativeErrorRecBinOutUp,      // pointer reference
1994                 relativeErrorRecBinOutLow,       // pointer reference
1995                 relativeStatisticalErrorIn,
1996                 relativeStatisticalErrorOut,
1997                 nominal,
1998                 nominalIn,
1999                 nominalOut,
2000                 columns,
2001                 rangeLow,
2002                 rangeUp,
2003                 readMe,
2004                 "recBin");
2005         if(relativeErrorRecBinOutUp) {
2006             // canvas with the error from regularization strength
2007             TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
2008             relativeErrorRecBin->Divide(2);
2009             relativeErrorRecBin->cd(1);
2010             Style(gPad, "GRID");
2011             relativeErrorRecBinInUp->DrawCopy("b");
2012             relativeErrorRecBinInLow->DrawCopy("same b");
2013             Style(AddLegend(gPad));
2014             relativeErrorRecBin->cd(2);
2015             Style(gPad, "GRID");
2016             relativeErrorRecBinOutUp->DrawCopy("b");
2017             relativeErrorRecBinOutLow->DrawCopy("same b");
2018             SavePadToPDF(relativeErrorRecBin);
2019             Style(AddLegend(gPad));
2020             relativeErrorRecBin->Write();
2021         }
2022     }
2023     if(methodIn && methodOut) {
2024         DoIntermediateSystematics(
2025                 methodIn,
2026                 methodOut,
2027                 relativeErrorMethodInUp,       // pointer reference
2028                 relativeErrorMethodInLow,      // pointer reference
2029                 relativeErrorMethodOutUp,      // pointer reference
2030                 relativeErrorMethodOutLow,     // pointer reference
2031                 relativeStatisticalErrorIn,
2032                 relativeStatisticalErrorOut,
2033                 nominal,
2034                 nominalIn,
2035                 nominalOut,
2036                 columns,
2037                 rangeLow,
2038                 rangeUp,
2039                 readMe,
2040                 "method"
2041                 );
2042         if(relativeErrorMethodInUp) {
2043             TCanvas* relativeErrorMethod(new TCanvas("relativeErrorMethod", "relativeErrorMethod"));
2044             relativeErrorMethod->Divide(2);
2045             relativeErrorMethod->cd(1);
2046             Style(gPad, "GRID");
2047             relativeErrorMethodInUp->DrawCopy("b");
2048             relativeErrorMethodInLow->DrawCopy("same b");
2049             Style(AddLegend(gPad));
2050             relativeErrorMethod->cd(2);
2051             Style(gPad, "GRID");
2052             relativeErrorMethodOutUp->DrawCopy("b");
2053             relativeErrorMethodOutLow->DrawCopy("same b");
2054             SavePadToPDF(relativeErrorMethod);
2055             Style(AddLegend(gPad));
2056             relativeErrorMethod->Write();
2057         }
2058     }
2059
2060     // and the placeholder for the final systematic
2061     TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2062     TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2063     TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2064     TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2065     relativeErrorInUp->SetYTitle("relative uncertainty");
2066     relativeErrorOutUp->SetYTitle("relative uncertainty");
2067     relativeErrorInLow->SetYTitle("relative uncertainty");
2068     relativeErrorOutLow->SetYTitle("relative uncertainty");
2069
2070     // sum of squares for the total systematic error
2071     Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.), eInUp(0.);
2072     Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.), eOutUp(0.);
2073     Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.), eInLow(0.);
2074     Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.), eOutLow(0.);
2075
2076     for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2077         // for the upper bound
2078         if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
2079         if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
2080         if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
2081         if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
2082         if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
2083         if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
2084         if(relativeErrorMethodInUp) dInUp = relativeErrorMethodInUp->GetBinContent(b+1);
2085         if(relativeErrorMethodOutUp) dOutUp = relativeErrorMethodOutUp->GetBinContent(b+1);
2086         eInUp  = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp + dInUp*dInUp;
2087         if(eInUp > 0) relativeErrorInUp->SetBinContent(b+1, 1.*TMath::Sqrt(eInUp));
2088         eOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp + dOutUp*dOutUp;
2089         if(eOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, 1.*TMath::Sqrt(eOutUp));
2090         // for the lower bound
2091         if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
2092         if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
2093         if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
2094         if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
2095         if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
2096         if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
2097         if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2098         if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
2099         if(fSymmRMS) {  // take first category as symmetric
2100 //            aInLow = aInUp*1.5;
2101 //            aOutLow = aOutUp*1.5;
2102             aInLow = aInUp;
2103             aOutLow = aOutUp;
2104             if(dInLow < dInUp) dInLow = dInUp;
2105             if(dOutLow < dOutUp) dOutLow = dOutUp;
2106         }
2107
2108         eInLow  = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow + dInLow*dInLow;
2109         if(eInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(eInLow));
2110         eOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow + dOutLow*dOutLow;
2111         if(eOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(eOutLow));
2112     }
2113     // project the estimated errors on the nominal ratio
2114     if(nominal) {
2115         Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2116         Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2117         Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2118         Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2119         Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2120         Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2121         Double_t lowErr(0.), upErr(0.);
2122         for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2123             // add the in and out of plane errors in quadrature
2124             // take special care here: to propagate the assymetric error, we need to correlate the
2125             // InLow with OutUp (minimum value of ratio) and InUp with OutLow (maximum value of ratio)
2126             lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2127             upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2128             // set the errors
2129             ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2130             ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2131             // get the bin width (which is the 'error' on x
2132             Double_t binWidth(nominal->GetBinWidth(i+1));
2133             axl[i] = binWidth/2.;
2134             axh[i] = binWidth/2.;
2135             // now get the coordinate for the point
2136             ax[i] = nominal->GetBinCenter(i+1);
2137             ay[i] = nominal->GetBinContent(i+1);
2138         }
2139         // save the nominal ratio
2140         TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2141         shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2142         nominalError->SetName("shape uncertainty");
2143         TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2144         nominalCanvas->Divide(2);
2145         nominalCanvas->cd(1);
2146         Style(nominal, kBlack);
2147         Style(nominalError, kCyan, kRatio);
2148         nominalError->SetLineColor(kCyan);
2149         nominalError->SetMarkerColor(kCyan);
2150         nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2151         nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2152         nominalError->DrawClone("a2");
2153         nominal->DrawCopy("same E1");
2154         Style(AddLegend(gPad));
2155         Style(gPad, "GRID");
2156         Style(nominalCanvas);
2157         // save nominal v2 and systematic errors
2158         TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2159                     nominalIn,
2160                     nominalOut,
2161                     fEventPlaneRes,
2162                     "v_{2} with shape uncertainty",
2163                     relativeErrorInUp,
2164                     relativeErrorInLow,
2165                     relativeErrorOutUp,
2166                     relativeErrorOutLow));
2167         shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2168         TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2169         nominalCanvas->cd(2);
2170         Style(nominalV2, kBlack);
2171         Style(nominalV2Error, kCyan, kV2);
2172         nominalV2Error->SetLineColor(kCyan);
2173         nominalV2Error->SetMarkerColor(kCyan);
2174         nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2175         nominalV2Error->DrawClone("a2");
2176         nominalV2->DrawClone("same E1");
2177         Style(AddLegend(gPad));
2178         Style(gPad, "GRID");
2179         Style(nominalCanvas);
2180         SavePadToPDF(nominalCanvas);
2181         nominalCanvas->Write();
2182     }
2183
2184     TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2185     relativeError->Divide(2);
2186     relativeError->cd(1);
2187     Style(gPad, "GRID");
2188     relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2189     Style(relativeErrorInUp, kBlue, kBar);
2190     Style(relativeErrorInLow, kGreen, kBar);
2191     relativeErrorInUp->DrawCopy("b");
2192     relativeErrorInLow->DrawCopy("same b");
2193     Style(relativeStatisticalErrorIn, kRed);
2194     relativeStatisticalErrorIn->DrawCopy("same");
2195     Style(AddLegend(gPad));
2196     relativeError->cd(2);
2197     Style(gPad, "GRID");
2198     relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2199     Style(relativeErrorOutUp, kBlue, kBar);
2200     Style(relativeErrorOutLow, kGreen, kBar);
2201     relativeErrorOutUp->DrawCopy("b");
2202     relativeErrorOutLow->DrawCopy("same b");
2203     Style(relativeStatisticalErrorOut, kRed);
2204     relativeStatisticalErrorOut->DrawCopy("same");
2205     Style(AddLegend(gPad));
2206
2207     // write the buffered file to disk and close the file
2208     SavePadToPDF(relativeError);
2209     relativeError->Write();
2210     output->Write();
2211 //    output->Close();
2212 }
2213 //_____________________________________________________________________________
2214     void AliJetFlowTools::DoIntermediateSystematics(
2215             TArrayI* variationsIn,                  // variantions in plane
2216             TArrayI* variationsOut,                 // variantions out of plane
2217             TH1D*& relativeErrorInUp,               // pointer reference to minimum relative error histogram in plane
2218             TH1D*& relativeErrorInLow,              // pointer reference to maximum relative error histogram in plane
2219             TH1D*& relativeErrorOutUp,              // pointer reference to minimum relative error histogram out of plane
2220             TH1D*& relativeErrorOutLow,             // pointer reference to maximum relative error histogram out of plane
2221             TH1D*& relativeStatisticalErrorIn,      // relative systematic error on ratio
2222             TH1D*& relativeStatisticalErrorOut,     // relative systematic error on ratio
2223             TH1D*& nominal,                         // clone of the nominal ratio in / out of plane
2224             TH1D*& nominalIn,                       // clone of the nominal in plane yield
2225             TH1D*& nominalOut,                      // clone of the nominal out of plane yield
2226             Int_t columns,                          // divide the output canvasses in this many columns
2227             Float_t rangeLow,                       // lower pt range
2228             Float_t rangeUp,                        // upper pt range
2229             TFile* readMe,                          // input file name (created by this unfolding class)
2230             TString source,                         // source of the variation
2231             Bool_t RMS                              // return RMS of distribution of variations as error
2232             ) const
2233 {
2234    // intermediate systematic check function. first index of supplied array is nominal value
2235    TList* listOfKeys((TList*)readMe->GetListOfKeys());
2236    if(!listOfKeys) {
2237        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2238        return;
2239    }
2240    // check input params
2241    if(variationsIn->GetSize() != variationsOut->GetSize()) {
2242        printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2243        return;
2244    }
2245    TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2246    TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2247    if(!(defRootDirIn && defRootDirOut)) {
2248        printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2249        return;
2250    }
2251    TString defIn(defRootDirIn->GetName());
2252    TString defOut(defRootDirOut->GetName());
2253
2254    // define lines to make the output prettier
2255    TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2256    TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2257    lineLow->SetLineColor(11);
2258    lineUp->SetLineColor(11);
2259    lineLow->SetLineWidth(3);
2260    lineUp->SetLineWidth(3);
2261
2262    // define an output histogram with the maximum relative error from this systematic constribution
2263    // if the option RMS is set to false, sigma is not really a standard deviation but holds the maximum (or minimum) relative value that the data has
2264    // reached in this function call.
2265    // if the option RMS is set to true, sigma holds the RMS value (equal to sigma as the mean is zero for relative errors) of the distribution of variations
2266    // which should correspond to a 68% confidence level
2267    relativeErrorInUp = new TH1D(Form("max #sigma/|x| from %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2268    relativeErrorInLow = new TH1D(Form("min #sigma/|x| from  %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2269    relativeErrorOutUp = new TH1D(Form("max #sigma/|x| from  %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2270    relativeErrorOutLow = new TH1D(Form("min #sigma/|x| from %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2271    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2272        if(!RMS) {
2273            relativeErrorInUp->SetBinContent(b+1, 1.);
2274            relativeErrorInUp->SetBinError(b+1, 0.);
2275            relativeErrorOutUp->SetBinContent(b+1, 1.);
2276            relativeErrorOutUp->SetBinError(b+1, .0);
2277            relativeErrorInLow->SetBinContent(b+1, 1.);
2278            relativeErrorInLow->SetBinError(b+1, 0.);
2279            relativeErrorOutLow->SetBinContent(b+1, 1.);
2280            relativeErrorOutLow->SetBinError(b+1, .0);
2281        } else if(RMS) {
2282            relativeErrorInUp->SetBinContent(b+1, 0.);
2283            relativeErrorInUp->SetBinError(b+1, 0.);
2284            relativeErrorOutUp->SetBinContent(b+1, 0.);
2285            relativeErrorOutUp->SetBinError(b+1, 0.);
2286            relativeErrorInLow->SetBinContent(b+1, 0.);
2287            relativeErrorInLow->SetBinError(b+1, 0.);
2288            relativeErrorOutLow->SetBinContent(b+1, 0.);
2289            relativeErrorOutLow->SetBinError(b+1, 0.);
2290        }
2291    }
2292    Int_t relativeErrorInUpN[100] = {0};
2293    Int_t relativeErrorOutUpN[100] = {0};
2294    Int_t relativeErrorInLowN[100] = {0};
2295    Int_t relativeErrorOutLowN[100] = {0};
2296
2297    // define an output histogram with the systematic error from this systematic constribution
2298    if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2299        relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "#sigma/|x|, statistical, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2300        relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "#sigma/|x|, statistical, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2301    }
2302
2303    // prepare necessary canvasses
2304    TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2305    TCanvas* canvasOut(0x0);
2306    if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2307    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2308    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2309    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2310    TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
2311    TCanvas* canvasSpectraOut(0x0);
2312    if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2313    TCanvas* canvasRatio(0x0);
2314    if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2315    TCanvas* canvasV2(0x0);
2316    if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2317    TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2318    TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2319    TCanvas* canvasMasterOut(0x0);
2320    if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2321    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2322
2323    TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2324    canvasProfiles->Divide(2);
2325    TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2326    TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2327    // get an estimate of the number of outputs and find the default set
2328
2329    Int_t rows = 1;
2330    columns = variationsIn->GetSize()-1;
2331    (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2332    canvasIn->Divide(columns, rows);
2333    if(canvasOut) canvasOut->Divide(columns, rows);
2334    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2335    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2336    canvasSpectraIn->Divide(columns, rows);
2337    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2338    if(canvasRatio) canvasRatio->Divide(columns, rows);
2339    if(canvasV2) canvasV2->Divide(columns, rows);
2340    canvasMasterIn->Divide(columns, rows);
2341    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2342
2343    // prepare a separate set of canvases to hold the nominal points
2344    TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2345    TCanvas* canvasNominalOut(0x0);
2346    if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2347    TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2348    TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2349    if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2350    TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
2351    TCanvas* canvasNominalSpectraOut(0x0);
2352    if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()),  Form("Nominal_%s_SpectraOut", source.Data()));
2353    TCanvas* canvasNominalRatio(0x0);
2354    if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2355    TCanvas* canvasNominalV2(0x0);
2356    if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2357    TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2358    TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2359    TCanvas* canvasNominalMasterOut(0x0);
2360    if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2361    (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2362
2363    canvasNominalSpectraIn->Divide(2);
2364    if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2365
2366    canvasNominalMasterIn->Divide(2);
2367    if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2368
2369    // extract the default output
2370    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2371    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2372    TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2373    TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2374    if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2375    if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2376    printf(" > succesfully extracted default results < \n");
2377
2378    // loop through the directories, only plot the graphs if the deconvolution converged
2379    TDirectoryFile* tempDirIn(0x0);
2380    TDirectoryFile* tempDirOut(0x0);
2381    TDirectoryFile* tempIn(0x0);
2382    TDirectoryFile* tempOut(0x0);
2383    TH1D* unfoldedSpectrumInForRatio(0x0);
2384    TH1D* unfoldedSpectrumOutForRatio(0x0);
2385    for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2386        tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2387        tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2388        if(!(tempDirIn && tempDirOut)) {
2389            printf(" > DoSystematics: couldn't get a set of variations < \n");
2390            continue;
2391        }
2392        TString dirNameIn(tempDirIn->GetName());
2393        TString dirNameOut(tempDirOut->GetName());
2394        // try to read the in- and out of plane subdirs
2395        tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2396        tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2397        j++;
2398        if(tempIn) {
2399            // to see if the unfolding converged try to extract the pearson coefficients
2400            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2401            if(pIn) {
2402                printf(" - %s in plane converged \n", dirNameIn.Data());
2403                canvasIn->cd(j);
2404                if(i==0) canvasNominalIn->cd(j);
2405                Style(gPad, "PEARSON");
2406                pIn->DrawCopy("colz");
2407                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2408                if(rIn) {
2409                    Style(rIn);
2410                    printf(" > found RatioRefoldedMeasured < \n");
2411                    canvasRatioMeasuredRefoldedIn->cd(j);
2412                    if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2413                    Style(gPad, "GRID");
2414                    rIn->SetFillColor(kRed);
2415                    rIn->Draw("ap");
2416                }
2417                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2418                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2419                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2420                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2421                if(dvector && avalue && rm && eff) {
2422                    Style(dvector);
2423                    Style(avalue);
2424                    Style(rm);
2425                    Style(eff);
2426                    canvasMISC->cd(1);
2427                    if(i==0) canvasNominalMISC->cd(1);
2428                    Style(gPad, "SPECTRUM");
2429                    dvector->DrawCopy();
2430                    canvasMISC->cd(2);
2431                    if(i==0) canvasNominalMISC->cd(2);
2432                    Style(gPad, "SPECTRUM");
2433                    avalue->DrawCopy();
2434                    canvasMISC->cd(3);
2435                    if(i==0) canvasNominalMISC->cd(3);
2436                    Style(gPad, "PEARSON");
2437                    rm->DrawCopy("colz");
2438                    canvasMISC->cd(4);
2439                    if(i==0) canvasNominalMISC->cd(4);
2440                    Style(gPad, "GRID");
2441                    eff->DrawCopy();
2442                } else if(rm && eff) {
2443                    Style(rm);
2444                    Style(eff);
2445                    canvasMISC->cd(3);
2446                    if(i==0) canvasNominalMISC->cd(3);
2447                    Style(gPad, "PEARSON");
2448                    rm->DrawCopy("colz");
2449                    canvasMISC->cd(4);
2450                    if(i==0) canvasNominalMISC->cd(4);
2451                    Style(gPad, "GRID");
2452                    eff->DrawCopy();
2453                }
2454            }
2455            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2456            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2457            unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2458            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2459            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2460                if(defaultUnfoldedJetSpectrumIn) {
2461                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2462                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2463                    temp->Divide(unfoldedSpectrum);
2464                    // get the absolute relative error
2465                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2466                        if(!RMS) {       // save the maximum deviation that a variation can cause
2467                            // the variation is HIGHER than the nominal point, so the bar goes UP
2468                            if( temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInUp->GetBinContent(b+1)) {
2469                                relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2470                                relativeErrorInUp->SetBinError(b+1, 0.);
2471                            }
2472                            // the variation is LOWER than the nominal point, so the bar goes DOWN
2473                            else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInLow->GetBinContent(b+1)) {
2474                                relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2475                                relativeErrorInLow->SetBinError(b+1, 0.);
2476                            }
2477                        } else if (RMS && !fSymmRMS) { // save info necessary for evaluating the RMS of a distribution of variations
2478                            printf(" oops shouldnt be here \n " );
2479                            if(temp->GetBinContent(b+1) < 1) {
2480                                relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2481                                relativeErrorInUpN[b]++;
2482                            }
2483                            // the variation is LOWER than the nominal point, so the bar goes DOWN
2484                            else if(temp->GetBinContent(b+1) > 1) {
2485                                relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2486                                relativeErrorInLowN[b]++;
2487                            }
2488                        } else if (fSymmRMS) {
2489                            // save symmetric sum of square to get a symmetric rms
2490                            relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2491                            relativeErrorInUpN[b]++;
2492                        }
2493                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2494                    }
2495                    temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2496                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2497                    temp->GetYaxis()->SetTitle("ratio");
2498                    canvasMasterIn->cd(j);
2499                    temp->GetYaxis()->SetRangeUser(0., 2);
2500                    Style(gPad, "GRID");
2501                    temp->DrawCopy();
2502                    canvasNominalMasterIn->cd(1);
2503                    Style(gPad, "GRID");
2504                    if(i > 0 ) {
2505                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2506                        tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2507                        Style(tempSyst, (EColor)(i+2));
2508                        if(i==1) tempSyst->DrawCopy();
2509                        else tempSyst->DrawCopy("same");
2510                    }
2511                }
2512                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2513                canvasSpectraIn->cd(j);
2514                if(i==0) canvasNominalSpectraIn->cd(1);
2515                Style(gPad);
2516                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2517                unfoldedSpectrum->DrawCopy();
2518                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2519                inputSpectrum->DrawCopy("same");
2520                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2521                refoldedSpectrum->DrawCopy("same");
2522                TLegend* l(AddLegend(gPad));
2523                Style(l);
2524                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2525                    Float_t chi(fitStatus->GetBinContent(1));
2526                    Float_t pen(fitStatus->GetBinContent(2));
2527                    Int_t dof((int)fitStatus->GetBinContent(3));
2528                    Float_t beta(fitStatus->GetBinContent(4));
2529                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2530                } else if (fitStatus) { // only available in SVD fit
2531                    Int_t reg((int)fitStatus->GetBinContent(1));
2532                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2533                }
2534                canvasNominalSpectraIn->cd(2);
2535                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2536                tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2537                Style(tempSyst, (EColor)(i+2));
2538                Style(gPad, "SPECTRUM");
2539                if(i==0) tempSyst->DrawCopy();
2540                else tempSyst->DrawCopy("same");
2541            }
2542        }
2543        if(tempOut) {
2544            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2545            if(pOut) {
2546                printf(" - %s out of plane converged \n", dirNameOut.Data());
2547                canvasOut->cd(j);
2548                if(i==0) canvasNominalOut->cd(j);
2549                Style(gPad, "PEARSON");
2550                pOut->DrawCopy("colz");
2551                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2552                if(rOut) {
2553                    Style(rOut);
2554                    printf(" > found RatioRefoldedMeasured < \n");
2555                    canvasRatioMeasuredRefoldedOut->cd(j);
2556                    if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2557                    Style(gPad, "GRID");
2558                    rOut->SetFillColor(kRed);
2559                    rOut->Draw("ap");
2560                }
2561                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2562                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2563                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2564                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2565                if(dvector && avalue && rm && eff) {
2566                    Style(dvector);
2567                    Style(avalue);
2568                    Style(rm);
2569                    Style(eff);
2570                    canvasMISC->cd(5);
2571                    if(i==0) canvasNominalMISC->cd(5);
2572                    Style(gPad, "SPECTRUM");
2573                    dvector->DrawCopy();
2574                    canvasMISC->cd(6);
2575                    if(i==0) canvasNominalMISC->cd(6);
2576                    Style(gPad, "SPECTRUM");
2577                    avalue->DrawCopy();
2578                    canvasMISC->cd(7);
2579                    if(i==0) canvasNominalMISC->cd(7);
2580                    Style(gPad, "PEARSON");
2581                    rm->DrawCopy("colz");
2582                    canvasMISC->cd(8);
2583                    if(i==0) canvasNominalMISC->cd(8);
2584                    Style(gPad, "GRID");
2585                    eff->DrawCopy();
2586                } else if(rm && eff) {
2587                    Style(rm);
2588                    Style(eff);
2589                    canvasMISC->cd(7);
2590                    if(i==0) canvasNominalMISC->cd(7);
2591                    Style(gPad, "PEARSON");
2592                    rm->DrawCopy("colz");
2593                    canvasMISC->cd(8);
2594                    if(i==0) canvasNominalMISC->cd(8);
2595                    Style(gPad, "GRID");
2596                    eff->DrawCopy();
2597                }
2598            }
2599            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2600            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2601            unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2602            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2603            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2604                if(defaultUnfoldedJetSpectrumOut) {
2605                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2606                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2607                    temp->Divide(unfoldedSpectrum);
2608                    // get the absolute relative error
2609                    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2610                        if(!RMS) {
2611                            // check if the error is larger than the current maximum
2612                            if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutUp->GetBinContent(b+1)) {
2613                                relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2614                                relativeErrorOutUp->SetBinError(b+1, 0.);
2615                            }
2616                            // check if the error is smaller than the current minimum
2617                            else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutLow->GetBinContent(b+1)) {
2618                                relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2619                                relativeErrorOutLow->SetBinError(b+1, 0.);
2620                            }
2621                        } else if (RMS && !fSymmRMS) {
2622                            printf(" OOps \n ");
2623                            if(temp->GetBinContent(b+1) < 1) {
2624                                relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2625                                relativeErrorOutUpN[b]++;
2626                            }
2627                            else if(temp->GetBinContent(b+1) > 1) {
2628                                relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2629                                relativeErrorOutLowN[b]++;
2630                            }
2631                        } else if (fSymmRMS) {
2632                            // save symmetric rms value
2633                            relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2634                            relativeErrorOutUpN[b]++;
2635                        }
2636                        if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2637                    }
2638                    temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2639                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2640                    temp->GetYaxis()->SetTitle("ratio");
2641                    canvasMasterOut->cd(j);
2642                    temp->GetYaxis()->SetRangeUser(0., 2);
2643                    Style(gPad, "GRID");
2644                    temp->DrawCopy();
2645                    canvasNominalMasterOut->cd(1);
2646                    Style(gPad, "GRID");
2647                    if(i > 0 ) {
2648                        TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2649                        tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2650                        Style(tempSyst, (EColor)(i+2));
2651                        if(i==1) tempSyst->DrawCopy();
2652                        else tempSyst->DrawCopy("same");
2653                    }
2654                }
2655                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2656                canvasSpectraOut->cd(j);
2657                if(i==0) canvasNominalSpectraOut->cd(1);
2658                Style(gPad);
2659                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2660                unfoldedSpectrum->DrawCopy();
2661                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2662                inputSpectrum->DrawCopy("same");
2663                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2664                refoldedSpectrum->DrawCopy("same");
2665                TLegend* l(AddLegend(gPad));
2666                Style(l);
2667                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2668                    Float_t chi(fitStatus->GetBinContent(1));
2669                    Float_t pen(fitStatus->GetBinContent(2));
2670                    Int_t dof((int)fitStatus->GetBinContent(3));
2671                    Float_t beta(fitStatus->GetBinContent(4));
2672                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2673                } else if (fitStatus) { // only available in SVD fit
2674                    Int_t reg((int)fitStatus->GetBinContent(1));
2675                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2676                }
2677                canvasNominalSpectraOut->cd(2);
2678                TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2679                tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2680                Style(tempSyst, (EColor)(i+2));
2681                Style(gPad, "SPECTRUM");
2682                if(i==0) tempSyst->DrawCopy();
2683                else tempSyst->DrawCopy("same");
2684            }
2685        }
2686        if(canvasRatio && canvasV2) {
2687            canvasRatio->cd(j);
2688            if(i==0) {
2689                canvasNominalRatio->cd(j);
2690                if(nominal && nominalIn && nominalOut) {
2691                    // if a nominal ratio is requested, delete the placeholder and update the nominal point
2692                    delete nominal;
2693                    delete nominalIn;
2694                    delete nominalOut;
2695                    nominalIn =  (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2696                    nominalOut =  (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2697                    nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2698                    nominal->Divide(nominal, unfoldedSpectrumOutForRatio);       // standard root divide for uncorrelated histos
2699                }
2700            }
2701            TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2702            Double_t _tempx(0), _tempy(0);
2703            if(ratioYield) {
2704                Style(ratioYield);
2705                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2706                    ratioYield->GetPoint(b, _tempx, _tempy);
2707                    ratioProfile->Fill(_tempx, _tempy);
2708                }
2709                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2710                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2711                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2712                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2713                ratioYield->SetFillColor(kRed);
2714                ratioYield->Draw("ap");
2715            }
2716            canvasV2->cd(j);
2717            if(i==0) canvasNominalV2->cd(j);
2718            TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2719            if(ratioV2) {
2720                Style(ratioV2);
2721                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2722                    ratioV2->GetPoint(b, _tempx, _tempy);
2723                    v2Profile->Fill(_tempx, _tempy);
2724
2725                }
2726                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2727                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2728                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2729                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2730                ratioV2->SetFillColor(kRed);
2731                ratioV2->Draw("ap");
2732            }
2733        }
2734        delete unfoldedSpectrumInForRatio;
2735        delete unfoldedSpectrumOutForRatio;
2736    }
2737    // save the canvasses
2738    canvasProfiles->cd(1);
2739    Style(ratioProfile);
2740    ratioProfile->DrawCopy();
2741    canvasProfiles->cd(2);
2742    Style(v2Profile);
2743    v2Profile->DrawCopy();
2744    SavePadToPDF(canvasProfiles);
2745    canvasProfiles->Write();
2746    SavePadToPDF(canvasIn);
2747    canvasIn->Write();
2748    if(canvasOut) {
2749        SavePadToPDF(canvasOut);
2750        canvasOut->Write();
2751    }
2752    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2753    canvasRatioMeasuredRefoldedIn->Write();
2754    if(canvasRatioMeasuredRefoldedOut) {
2755        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2756        canvasRatioMeasuredRefoldedOut->Write();
2757    }
2758    SavePadToPDF(canvasSpectraIn);
2759    canvasSpectraIn->Write();
2760    if(canvasSpectraOut) {
2761        SavePadToPDF(canvasSpectraOut);
2762        canvasSpectraOut->Write();
2763        SavePadToPDF(canvasRatio);
2764        canvasRatio->Write();
2765        SavePadToPDF(canvasV2);
2766        canvasV2->Write();
2767    }
2768    SavePadToPDF(canvasMasterIn);
2769    canvasMasterIn->Write();
2770    if(canvasMasterOut) {
2771        SavePadToPDF(canvasMasterOut);
2772        canvasMasterOut->Write();
2773    }
2774    SavePadToPDF(canvasMISC);
2775    canvasMISC->Write();
2776    // save the nomial canvasses
2777    SavePadToPDF(canvasNominalIn);
2778    canvasNominalIn->Write();
2779    if(canvasNominalOut) {
2780        SavePadToPDF(canvasNominalOut);
2781        canvasNominalOut->Write();
2782    }
2783    SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
2784    canvasNominalRatioMeasuredRefoldedIn->Write();
2785    if(canvasNominalRatioMeasuredRefoldedOut) {
2786        SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
2787        canvasNominalRatioMeasuredRefoldedOut->Write();
2788    }
2789    canvasNominalSpectraIn->cd(2);
2790    Style(AddLegend(gPad));
2791    SavePadToPDF(canvasNominalSpectraIn);
2792    canvasNominalSpectraIn->Write();
2793    if(canvasNominalSpectraOut) {
2794        canvasNominalSpectraOut->cd(2);
2795        Style(AddLegend(gPad));
2796        SavePadToPDF(canvasNominalSpectraOut);
2797        canvasNominalSpectraOut->Write();
2798        SavePadToPDF(canvasNominalRatio);
2799        canvasNominalRatio->Write();
2800        SavePadToPDF(canvasNominalV2);
2801        canvasNominalV2->Write();
2802    }
2803    canvasNominalMasterIn->cd(1);
2804    Style(AddLegend(gPad));
2805    lineUp->DrawClone("same");
2806    lineLow->DrawClone("same");
2807    SavePadToPDF(canvasNominalMasterIn);
2808    canvasNominalMasterIn->Write();
2809    if(canvasNominalMasterOut) {
2810        canvasNominalMasterOut->cd(1);
2811        Style(AddLegend(gPad));
2812        lineUp->DrawClone("same");
2813        lineLow->DrawClone("same");
2814        SavePadToPDF(canvasNominalMasterOut);
2815        canvasNominalMasterOut->Write();
2816    }
2817    SavePadToPDF(canvasNominalMISC);
2818    canvasNominalMISC->Write();
2819
2820    // save the relative errors
2821    for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2822        // to arrive at a min and max from here, combine in up and out low
2823        if(!RMS) {
2824            relativeErrorInUp->SetBinContent(b+1, -1.*(relativeErrorInUp->GetBinContent(b+1)-1));
2825            relativeErrorInUp->SetBinError(b+1, 0.);
2826            relativeErrorOutUp->SetBinContent(b+1, -1.*(relativeErrorOutUp->GetBinContent(b+1)-1));
2827            relativeErrorOutUp->SetBinError(b+1, .0);
2828            relativeErrorInLow->SetBinContent(b+1, -1.*(relativeErrorInLow->GetBinContent(b+1)-1));
2829            relativeErrorInLow->SetBinError(b+1, 0.);
2830            relativeErrorOutLow->SetBinContent(b+1, -1.*(relativeErrorOutLow->GetBinContent(b+1)-1));
2831            relativeErrorOutLow->SetBinError(b+1, .0);
2832        } else if (RMS) {
2833            // these guys are already stored as percentages, so no need to remove the offset of 1
2834            // RMS is defined as sqrt(sum(squared))/N
2835            // min is defined as negative, max is defined as positive
2836            if(!fSymmRMS) {
2837                if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
2838                if(relativeErrorInLowN[b] < 1) relativeErrorInLowN[b] = 1;
2839                if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
2840                if(relativeErrorOutLowN[b] < 1) relativeErrorOutLowN[b] = 1;
2841                relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
2842                relativeErrorInUp->SetBinError(b+1, 0.);
2843                relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
2844                relativeErrorOutUp->SetBinError(b+1, .0);
2845                relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorInLow->GetBinContent(b+1)/relativeErrorInLowN[b]));
2846                relativeErrorInLow->SetBinError(b+1, 0.);
2847                relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorOutLow->GetBinContent(b+1)/relativeErrorOutLowN[b]));
2848                relativeErrorOutLow->SetBinError(b+1, .0);
2849            } else if (fSymmRMS) {
2850                if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
2851                if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
2852                relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
2853                relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
2854            }
2855        }
2856    }
2857    relativeErrorInUp->SetYTitle("relative uncertainty");
2858    relativeErrorOutUp->SetYTitle("relative uncertainty");
2859    relativeErrorInLow->SetYTitle("relative uncertainty");
2860    relativeErrorOutLow->SetYTitle("relative uncertainty");
2861    relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2862    relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2863    relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2864    relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2865
2866    canvasNominalMasterIn->cd(2);
2867    Style(gPad, "GRID");
2868    Style(relativeErrorInUp, kBlue, kBar);
2869    Style(relativeErrorInLow, kGreen, kBar);
2870    relativeErrorInUp->DrawCopy("b");
2871    relativeErrorInLow->DrawCopy("same b");
2872    Style(AddLegend(gPad));
2873    SavePadToPDF(canvasNominalMasterIn);
2874    canvasNominalMasterIn->Write();
2875    canvasNominalMasterOut->cd(2);
2876    Style(gPad, "GRID");
2877    Style(relativeErrorOutUp, kBlue, kBar);
2878    Style(relativeErrorOutLow, kGreen, kBar);
2879    relativeErrorOutUp->DrawCopy("b");
2880    relativeErrorOutLow->DrawCopy("same b");
2881    Style(AddLegend(gPad));
2882    SavePadToPDF(canvasNominalMasterOut);
2883    canvasNominalMasterOut->Write();
2884 }
2885 //_____________________________________________________________________________
2886 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
2887 {
2888    // go through the output file and perform post processing routines
2889    // can either be performed in one go with the unfolding, or at a later stage
2890    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
2891    TFile readMe(in.Data(), "READ");     // open file read-only
2892    if(readMe.IsZombie()) {
2893        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2894        return;
2895    }
2896    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
2897    readMe.ls();
2898    TList* listOfKeys((TList*)readMe.GetListOfKeys());
2899    if(!listOfKeys) {
2900        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2901        return;
2902    }
2903    // prepare necessary canvasses
2904    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
2905    TCanvas* canvasOut(0x0);
2906    if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
2907    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
2908    TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2909    if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
2910    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
2911    TCanvas* canvasSpectraOut(0x0);
2912    if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
2913    TCanvas* canvasRatio(0x0);
2914    if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
2915    TCanvas* canvasV2(0x0);
2916    if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
2917    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
2918    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
2919    TCanvas* canvasMasterOut(0x0);
2920    if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
2921    (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2922    TDirectoryFile* defDir(0x0);
2923    TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
2924    canvasProfiles->Divide(2);
2925    TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2926    TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2927    // get an estimate of the number of outputs and find the default set
2928    Int_t cacheMe(0);
2929    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
2930        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
2931            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2932            cacheMe++;
2933        }
2934    }
2935    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
2936    canvasIn->Divide(columns, rows);
2937    if(canvasOut) canvasOut->Divide(columns, rows);
2938    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2939    if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2940    canvasSpectraIn->Divide(columns, rows);
2941    if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2942    if(canvasRatio) canvasRatio->Divide(columns, rows);
2943    if(canvasV2) canvasV2->Divide(columns, rows);
2944
2945    canvasMasterIn->Divide(columns, rows);
2946    if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2947    // extract the default output
2948    TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2949    TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2950    if(defDir) {
2951        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
2952        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
2953        if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
2954        if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
2955        printf(" > succesfully extracted default results < \n");
2956    }
2957
2958    // loop through the directories, only plot the graphs if the deconvolution converged
2959    TDirectoryFile* tempDir(0x0);
2960    TDirectoryFile* tempIn(0x0);
2961    TDirectoryFile*  tempOut(0x0);
2962    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
2963        // read the directory by its name
2964        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2965        if(!tempDir) continue;
2966        TString dirName(tempDir->GetName());
2967        // try to read the in- and out of plane subdirs
2968        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
2969        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
2970        j++;
2971        if(!tempIn) {    // attempt to get MakeAU output
2972            TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
2973            TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
2974            tempPearson->Divide(4,2);
2975            TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
2976            tempRatio->Divide(4,2);
2977            TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
2978            tempResult->Divide(4,2);
2979            for(Int_t q(0); q < 8; q++) {
2980                tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
2981                if(tempIn) {
2982                        // to see if the unfolding converged try to extract the pearson coefficients
2983                        TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2984                        if(pIn) {
2985                        printf(" - %s in plane converged \n", dirName.Data());
2986                            tempPearson->cd(1+q);
2987                             Style(gPad, "PEARSON");
2988                             pIn->DrawCopy("colz");
2989                             TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2990                             if(rIn) {
2991                                 Style(rIn);
2992                                 printf(" > found RatioRefoldedMeasured < \n");
2993                                 tempRatio->cd(q+1);
2994                                 rIn->SetFillColor(kRed);
2995                                 rIn->Draw("ap");
2996                             }
2997                             TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2998                             TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2999                             TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3000                             TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3001                             if(dvector && avalue && rm && eff) {
3002                                 Style(dvector);
3003                                 Style(avalue);
3004                                 Style(rm);
3005                                 Style(eff);
3006                                 canvasMISC->cd(1);
3007                                 Style(gPad, "SPECTRUM");
3008                                 dvector->DrawCopy();
3009                                 canvasMISC->cd(2);
3010                                 Style(gPad, "SPECTRUM");
3011                                 avalue->DrawCopy();
3012                                 canvasMISC->cd(3);
3013                                 Style(gPad, "PEARSON");
3014                                 rm->DrawCopy("colz");
3015                                 canvasMISC->cd(4);
3016                                 Style(gPad, "GRID");
3017                                 eff->DrawCopy();
3018                             } else if(rm && eff) {
3019                                 Style(rm);
3020                                 Style(eff);
3021                                 canvasMISC->cd(3);
3022                                 Style(gPad, "PEARSON");
3023                                 rm->DrawCopy("colz");
3024                                 canvasMISC->cd(4);
3025                                 Style(gPad, "GRID");
3026                                 eff->DrawCopy();
3027                             }
3028                         }
3029                        TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3030                        TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3031                        TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3032                        if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3033                            if(defaultUnfoldedJetSpectrumIn) {
3034                                Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3035                                TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3036                                temp->Divide(unfoldedSpectrum);
3037                                temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3038                                temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3039                                temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3040                                canvasMasterIn->cd(j);
3041                                temp->GetYaxis()->SetRangeUser(0., 2);
3042                                temp->DrawCopy();
3043                            }
3044                            TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3045                            tempResult->cd(q+1);
3046                            Style(gPad);
3047                            Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3048                            unfoldedSpectrum->DrawCopy();
3049                            Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3050                            inputSpectrum->DrawCopy("same");
3051                            Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3052                            refoldedSpectrum->DrawCopy("same");
3053                            TLegend* l(AddLegend(gPad));
3054                            Style(l);
3055                            if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3056                                Float_t chi(fitStatus->GetBinContent(1));
3057                                Float_t pen(fitStatus->GetBinContent(2));
3058                                Int_t dof((int)fitStatus->GetBinContent(3));
3059                                Float_t beta(fitStatus->GetBinContent(4));
3060                                l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3061                            } else if (fitStatus) { // only available in SVD fit
3062                                Int_t reg((int)fitStatus->GetBinContent(1));
3063                                l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3064                            }
3065                        }
3066                }
3067            }
3068        }
3069        if(tempIn) {
3070            // to see if the unfolding converged try to extract the pearson coefficients
3071            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3072            if(pIn) {
3073                printf(" - %s in plane converged \n", dirName.Data());
3074                canvasIn->cd(j);
3075                Style(gPad, "PEARSON");
3076                pIn->DrawCopy("colz");
3077                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3078                if(rIn) {
3079                    Style(rIn);
3080                    printf(" > found RatioRefoldedMeasured < \n");
3081                    canvasRatioMeasuredRefoldedIn->cd(j);
3082                    rIn->SetFillColor(kRed);
3083                    rIn->Draw("ap");
3084                }
3085                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3086                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3087                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3088                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3089                if(dvector && avalue && rm && eff) {
3090                    Style(dvector);
3091                    Style(avalue);
3092                    Style(rm);
3093                    Style(eff);
3094                    canvasMISC->cd(1);
3095                    Style(gPad, "SPECTRUM");
3096                    dvector->DrawCopy();
3097                    canvasMISC->cd(2);
3098                    Style(gPad, "SPECTRUM");
3099                    avalue->DrawCopy();
3100                    canvasMISC->cd(3);
3101                    Style(gPad, "PEARSON");
3102                    rm->DrawCopy("colz");
3103                    canvasMISC->cd(4);
3104                    Style(gPad, "GRID");
3105                    eff->DrawCopy();
3106                } else if(rm && eff) {
3107                    Style(rm);
3108                    Style(eff);
3109                    canvasMISC->cd(3);
3110                    Style(gPad, "PEARSON");
3111                    rm->DrawCopy("colz");
3112                    canvasMISC->cd(4);
3113                    Style(gPad, "GRID");
3114                    eff->DrawCopy();
3115                }
3116            }
3117            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3118            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3119            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3120            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3121                if(defaultUnfoldedJetSpectrumIn) {
3122                    Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3123                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3124                    temp->Divide(unfoldedSpectrum);
3125                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3126                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3127                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3128                    canvasMasterIn->cd(j);
3129                    temp->GetYaxis()->SetRangeUser(0., 2);
3130                    temp->DrawCopy();
3131                }
3132                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3133                canvasSpectraIn->cd(j);
3134                Style(gPad);
3135                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3136                unfoldedSpectrum->DrawCopy();
3137                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3138                inputSpectrum->DrawCopy("same");
3139                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3140                refoldedSpectrum->DrawCopy("same");
3141                TLegend* l(AddLegend(gPad));
3142                Style(l);
3143                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3144                    Float_t chi(fitStatus->GetBinContent(1));
3145                    Float_t pen(fitStatus->GetBinContent(2));
3146                    Int_t dof((int)fitStatus->GetBinContent(3));
3147                    Float_t beta(fitStatus->GetBinContent(4));
3148                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3149                } else if (fitStatus) { // only available in SVD fit
3150                    Int_t reg((int)fitStatus->GetBinContent(1));
3151                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3152                }
3153            }
3154        }
3155        if(tempOut) {
3156            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
3157            if(pOut) {
3158                printf(" - %s out of plane converged \n", dirName.Data());
3159                canvasOut->cd(j);
3160                Style(gPad, "PEARSON");
3161                pOut->DrawCopy("colz");
3162                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3163                if(rOut) {
3164                    Style(rOut);
3165                    printf(" > found RatioRefoldedMeasured < \n");
3166                    canvasRatioMeasuredRefoldedOut->cd(j);
3167                    rOut->SetFillColor(kRed);
3168                    rOut->Draw("ap");
3169                }
3170                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
3171                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
3172                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
3173                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
3174                if(dvector && avalue && rm && eff) {
3175                    Style(dvector);
3176                    Style(avalue);
3177                    Style(rm);
3178                    Style(eff);
3179                    canvasMISC->cd(5);
3180                    Style(gPad, "SPECTRUM");
3181                    dvector->DrawCopy();
3182                    canvasMISC->cd(6);
3183                    Style(gPad, "SPECTRUM");
3184                    avalue->DrawCopy();
3185                    canvasMISC->cd(7);
3186                    Style(gPad, "PEARSON");
3187                    rm->DrawCopy("colz");
3188                    canvasMISC->cd(8);
3189                    Style(gPad, "GRID");
3190                    eff->DrawCopy();
3191                } else if(rm && eff) {
3192                    Style(rm);
3193                    Style(eff);
3194                    canvasMISC->cd(7);
3195                    Style(gPad, "PEARSON");
3196                    rm->DrawCopy("colz");
3197                    canvasMISC->cd(8);
3198                    Style(gPad, "GRID");
3199                    eff->DrawCopy();
3200                }
3201            }
3202            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3203            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3204            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3205            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3206                if(defaultUnfoldedJetSpectrumOut) {
3207                    Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3208                    TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3209                    temp->Divide(unfoldedSpectrum);
3210                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3211                    temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3212                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3213                    canvasMasterOut->cd(j);
3214                    temp->GetYaxis()->SetRangeUser(0., 2.);
3215                    temp->DrawCopy();
3216                }
3217                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3218                canvasSpectraOut->cd(j);
3219                Style(gPad);
3220                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3221                unfoldedSpectrum->DrawCopy();
3222                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3223                inputSpectrum->DrawCopy("same");
3224                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3225                refoldedSpectrum->DrawCopy("same");
3226                TLegend* l(AddLegend(gPad));
3227                Style(l);
3228                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3229                    Float_t chi(fitStatus->GetBinContent(1));
3230                    Float_t pen(fitStatus->GetBinContent(2));
3231                    Int_t dof((int)fitStatus->GetBinContent(3));
3232                    Float_t beta(fitStatus->GetBinContent(4));
3233                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3234                } else if (fitStatus) { // only available in SVD fit
3235                    Int_t reg((int)fitStatus->GetBinContent(1));
3236                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3237                }
3238            }
3239        }
3240        if(canvasRatio && canvasV2) {
3241            canvasRatio->cd(j);
3242            TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
3243            Double_t _tempx(0), _tempy(0);
3244            if(ratioYield) {
3245                Style(ratioYield);
3246                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3247                    ratioYield->GetPoint(b, _tempx, _tempy);
3248                    ratioProfile->Fill(_tempx, _tempy);
3249                }
3250                ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3251                ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3252                ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3253                ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3254                ratioYield->SetFillColor(kRed);
3255                ratioYield->Draw("ap");
3256            }
3257            canvasV2->cd(j);
3258            TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3259            if(ratioV2) {
3260                Style(ratioV2);
3261                for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3262                    ratioV2->GetPoint(b, _tempx, _tempy);
3263                    v2Profile->Fill(_tempx, _tempy);
3264
3265                }
3266                v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3267                v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3268                ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3269                ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3270                ratioV2->SetFillColor(kRed);
3271                ratioV2->Draw("ap");
3272            }
3273        }
3274    }
3275    TFile output(out.Data(), "RECREATE");
3276    canvasProfiles->cd(1);
3277    Style(ratioProfile);
3278    ratioProfile->DrawCopy();
3279    canvasProfiles->cd(2);
3280    Style(v2Profile);
3281    v2Profile->DrawCopy();
3282    SavePadToPDF(canvasProfiles);
3283    canvasProfiles->Write();
3284    SavePadToPDF(canvasIn);
3285    canvasIn->Write();
3286    if(canvasOut) {
3287        SavePadToPDF(canvasOut);
3288        canvasOut->Write();
3289    }
3290    SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3291    canvasRatioMeasuredRefoldedIn->Write();
3292    if(canvasRatioMeasuredRefoldedOut) {
3293        SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3294        canvasRatioMeasuredRefoldedOut->Write();
3295    }
3296    SavePadToPDF(canvasSpectraIn);
3297    canvasSpectraIn->Write();
3298    if(canvasSpectraOut) {
3299        SavePadToPDF(canvasSpectraOut);
3300        canvasSpectraOut->Write();
3301        SavePadToPDF(canvasRatio);
3302        canvasRatio->Write();
3303        SavePadToPDF(canvasV2);
3304        canvasV2->Write();
3305    }
3306    SavePadToPDF(canvasMasterIn);
3307    canvasMasterIn->Write();
3308    if(canvasMasterOut) {
3309        SavePadToPDF(canvasMasterOut);
3310        canvasMasterOut->Write();
3311    }
3312    SavePadToPDF(canvasMISC);
3313    canvasMISC->Write();
3314    output.Write();
3315    output.Close();
3316 }
3317 //_____________________________________________________________________________
3318 Bool_t AliJetFlowTools::SetRawInput (
3319         TH2D* detectorResponse,  // detector response matrix
3320         TH1D* jetPtIn,           // in plane jet spectrum
3321         TH1D* jetPtOut,          // out of plane jet spectrum
3322         TH1D* dptIn,             // in plane delta pt distribution
3323         TH1D* dptOut,            // out of plane delta pt distribution
3324         Int_t eventCount) {
3325     // set input histograms manually
3326     fDetectorResponse   = detectorResponse;
3327     fSpectrumIn         = jetPtIn;
3328     fSpectrumOut        = jetPtOut;
3329     fDptInDist          = dptIn;
3330     fDptOutDist         = dptOut;
3331     fRawInputProvided   = kTRUE;
3332     // check if all data is provided
3333     if(!fDetectorResponse) {
3334         printf(" fDetectorResponse not found \n ");
3335         return kFALSE;
3336     }
3337     // check if the pt bin for true and rec have been set
3338     if(!fBinsTrue || !fBinsRec) {
3339         printf(" No true or rec bins set, please set binning ! \n");
3340         return kFALSE;
3341     }
3342     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
3343                         // procedures, these profiles will be nonsensical, user is responsible
3344         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3345         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3346         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3347     }
3348     // normalize spectra to event count if requested
3349     if(fNormalizeSpectra) {
3350         fEventCount = eventCount;
3351         if(fEventCount > 0) {
3352             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
3353             fSpectrumOut->Sumw2();
3354             fSpectrumIn->Scale(1./((double)fEventCount));
3355             fSpectrumOut->Scale(1./((double)fEventCount));
3356         }
3357     }
3358     if(!fNormalizeSpectra && fEventCount > 0) {
3359         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
3360         fSpectrumOut->Sumw2();
3361         fSpectrumIn->Scale(1./((double)fEventCount));
3362         fSpectrumOut->Scale(1./((double)fEventCount));
3363     }
3364     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
3365     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
3366     fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3367     fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3368     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
3369     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
3370     fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3371     fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3372
3373     return kTRUE;
3374 }
3375 //_____________________________________________________________________________
3376 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
3377 {
3378     // return ratio of h1 / h2
3379     // histograms can have different binning. errors are propagated as uncorrelated
3380     if(!(h1 && h2) ) {
3381         printf(" GetRatio called with NULL argument(s) \n ");
3382         return 0x0;
3383     }
3384     Int_t j(0);
3385     TGraphErrors *gr = new TGraphErrors();
3386     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3387     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3388     TH1* dud((TH1*)h1->Clone("dud"));
3389     dud->Sumw2();
3390     h1->Sumw2();
3391     h2->Sumw2();
3392     if(!dud->Divide(h1, h2)) {
3393         printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
3394         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3395             binCent = h1->GetXaxis()->GetBinCenter(i);
3396             if(xmax > 0. && binCent > xmax) continue;
3397             j = h2->FindBin(binCent);
3398             binWidth = h1->GetXaxis()->GetBinWidth(i);
3399             if(h2->GetBinContent(j) > 0.) {
3400                 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
3401                 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
3402                 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
3403                 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
3404                 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
3405                 gr->SetPoint(i-1, binCent, ratio);
3406                 gr->SetPointError(i-1, 0.5*binWidth, error2);
3407             }
3408         }
3409     } else {
3410         printf( " > using ROOT to divide two histograms < \n");
3411         for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3412             binCent = dud->GetXaxis()->GetBinCenter(i);
3413             if(xmax > 0. && binCent > xmax) continue;
3414             binWidth = dud->GetXaxis()->GetBinWidth(i);
3415             gr->SetPoint(i-1,binCent,dud->GetBinContent(i));
3416             gr->SetPointError(i-1, 0.5*binWidth,dud->GetBinError(i));
3417         }
3418     }
3419
3420     if(appendFit) {
3421         TF1* fit(new TF1("lin", "pol0", 10, 100));
3422         gr->Fit(fit);
3423     }
3424     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3425     if(dud) delete dud;
3426     return gr;
3427 }
3428 //_____________________________________________________________________________
3429 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
3430 {
3431     // get v2 from difference of in plane, out of plane yield
3432     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
3433     // different binning is allowed but will mean that the error
3434     // propagation is unreliable
3435     // r is the event plane resolution for the chosen centrality
3436     if(!(h1 && h2) ) {
3437         printf(" GetV2 called with NULL argument(s) \n ");
3438         return 0x0;
3439     }
3440     TGraphErrors *gr = new TGraphErrors();
3441     gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3442     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3443     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
3444
3445     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3446         binCent = h1->GetXaxis()->GetBinCenter(i);
3447         binWidth = h1->GetXaxis()->GetBinWidth(i);
3448         if(h2->GetBinContent(i) > 0.) {
3449             in = h1->GetBinContent(i);
3450             ein = h1->GetBinError(i);
3451             out = h2->GetBinContent(i);
3452             eout = h2->GetBinError(i);
3453             ratio = pre*((in-out)/(in+out));
3454             error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
3455             error2 = error2*pre*pre;
3456             if(error2 > 0) error2 = TMath::Sqrt(error2);
3457             gr->SetPoint(i-1,binCent,ratio);
3458             gr->SetPointError(i-1,0.5*binWidth,error2);
3459         }
3460     }
3461     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3462     return gr;
3463 }
3464 //_____________________________________________________________________________
3465 TGraphAsymmErrors* AliJetFlowTools::GetV2WithSystematicErrors(
3466         TH1* h1, TH1* h2, Double_t r, TString name,
3467         TH1* relativeErrorInUp,
3468         TH1* relativeErrorInLow,
3469         TH1* relativeErrorOutUp,
3470         TH1* relativeErrorOutLow,
3471         Float_t rho) const
3472 {
3473     // get v2 with asymmetric systematic error
3474     // note that this is ONLY the systematic error, no statistical error!
3475     // rho is the pearson correlation coefficient
3476     TGraphErrors* tempV2(GetV2(h1, h2, r, name));
3477     Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
3478     Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
3479     Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
3480     Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
3481     Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
3482     Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
3483     Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
3484     // loop through the bins and do error propagation
3485     for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3486         // extract the absolute errors
3487         in = h1->GetBinContent(i+1);
3488         einUp = TMath::Abs(in*relativeErrorInUp->GetBinContent(i+1));
3489         einLow = TMath::Abs(in*relativeErrorInLow->GetBinContent(1+i));
3490         out = h2->GetBinContent(i+1);
3491         eoutUp = TMath::Abs(out*relativeErrorOutUp->GetBinContent(1+i));
3492         eoutLow = TMath::Abs(out*relativeErrorOutLow->GetBinContent(1+i));
3493         // get the error squared
3494         if(rho <= 0) {
3495             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)))*eoutLow*eoutLow);
3496             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)))*eoutUp*eoutUp);
3497         } else {
3498             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);
3499             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);
3500         }
3501         if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
3502         if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
3503         // set the errors
3504         ayh[i] = error2Up;
3505         ayl[i] = error2Low;
3506         // get the bin width (which is the 'error' on x)
3507         Double_t binWidth(h1->GetBinWidth(i+1));
3508         axl[i] = binWidth/2.;
3509         axh[i] = binWidth/2.;
3510         // now get the coordinate for the poin
3511         tempV2->GetPoint(i, ax[i], ay[i]);
3512     }
3513     // save the nominal ratio
3514     TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
3515     nominalError->SetName("v_{2} with shape uncertainty");
3516     // do some memory management
3517     delete tempV2;
3518     delete[] ax;
3519     delete[] ay;
3520     delete[] axh;
3521     delete[] axl;
3522     delete[] ayh;
3523     delete[] ayl;
3524
3525     return nominalError;
3526 }
3527 //_____________________________________________________________________________
3528 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
3529     // write object, if a unique identifier is given the object is cloned
3530     // and the clone is saved. setting kill to true will delete the original obect from the heap
3531     if(!object) {
3532         printf(" > WriteObject:: called with NULL arguments \n ");
3533         return;
3534     } else if(!strcmp("", suffix.Data())) object->Write();
3535     else {
3536         TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
3537         newObject->Write();
3538     }
3539     if(kill) delete object;
3540 }
3541 //_____________________________________________________________________________
3542 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
3543     // construt a delta pt response matrix from supplied dpt distribution
3544     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
3545     // do a weighted rebinning to a (coarser) dpt distribution
3546     // be careful with the binning of the dpt response: it should be equal to that
3547     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
3548     //
3549     // the response matrix will be square and have the same binning
3550     // (min, max and granularity) of the input histogram
3551     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
3552     Double_t _bins[bins+1];             // prepare array with bin borders
3553     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
3554     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
3555     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
3556     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
3557         Bool_t skip = kFALSE;
3558         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
3559             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
3560             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
3561         }
3562     }
3563     return res;
3564 }
3565 //_____________________________________________________________________________
3566 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
3567     if(!binsTrue || !binsRec) {
3568         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
3569         return 0x0;
3570     }
3571     TString name(Form("unityResponse_%s", suffix.Data()));
3572     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
3573     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
3574         for(Int_t j(0); j < binsRec->GetSize(); j++) {
3575             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
3576         }
3577     }
3578     return unity;
3579 }
3580 //_____________________________________________________________________________
3581 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
3582     // save configuration parameters to histogram
3583     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
3584     summary->SetBinContent(1, fBetaIn);
3585     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
3586     summary->SetBinContent(2, fBetaOut);
3587     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
3588     summary->SetBinContent(3, fCentralityArray->At(0));
3589     summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]");
3590     summary->SetBinContent(4, (int)convergedIn);
3591     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
3592     summary->SetBinContent(5, (int)convergedOut);
3593     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
3594     summary->SetBinContent(6, (int)fAvoidRoundingError);
3595     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
3596     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
3597     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
3598     summary->SetBinContent(8, (int)fPrior);
3599     summary->GetXaxis()->SetBinLabel(8, "fPrior");
3600     summary->SetBinContent(9, fSVDRegIn);
3601     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
3602     summary->SetBinContent(10, fSVDRegOut);
3603     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
3604     summary->SetBinContent(11, (int)fSVDToy);
3605     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
3606     summary->SetBinContent(12, fJetRadius);
3607     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
3608     summary->SetBinContent(13, (int)fNormalizeSpectra);
3609     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
3610     summary->SetBinContent(14, (int)fSmoothenPrior);
3611     summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
3612     summary->SetBinContent(15, (int)fTestMode);
3613     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
3614     summary->SetBinContent(16, (int)fUseDetectorResponse);
3615     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
3616     summary->SetBinContent(17, fBayesianIterIn);
3617     summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
3618     summary->SetBinContent(18, fBayesianIterOut);
3619     summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
3620     summary->SetBinContent(19, fBayesianSmoothIn);
3621     summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
3622     summary->SetBinContent(20, fBayesianSmoothOut);
3623     summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
3624 }
3625 //_____________________________________________________________________________
3626 void AliJetFlowTools::ResetAliUnfolding() {
3627      // ugly function: reset all unfolding parameters
3628      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
3629      if(fitter) {
3630          printf(" > Found fitter, will delete it < \n");
3631          delete fitter;
3632      }
3633      if(gMinuit) {
3634          printf(" > Found gMinuit, will re-create it < \n");
3635          delete gMinuit;
3636          gMinuit = new TMinuit;
3637      }
3638      AliUnfolding::fgCorrelationMatrix = 0;
3639      AliUnfolding::fgCorrelationMatrixSquared = 0;
3640      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
3641      AliUnfolding::fgCurrentESDVector = 0;
3642      AliUnfolding::fgEntropyAPriori = 0;
3643      AliUnfolding::fgEfficiency = 0;
3644      AliUnfolding::fgUnfoldedAxis = 0;
3645      AliUnfolding::fgMeasuredAxis = 0;
3646      AliUnfolding::fgFitFunction = 0;
3647      AliUnfolding::fgMaxInput  = -1;
3648      AliUnfolding::fgMaxParams = -1;
3649      AliUnfolding::fgOverflowBinLimit = -1;
3650      AliUnfolding::fgRegularizationWeight = 10000;
3651      AliUnfolding::fgSkipBinsBegin = 0;
3652      AliUnfolding::fgMinuitStepSize = 0.1;
3653      AliUnfolding::fgMinuitPrecision = 1e-6;
3654      AliUnfolding::fgMinuitMaxIterations = 1000000;
3655      AliUnfolding::fgMinuitStrategy = 1.;
3656      AliUnfolding::fgMinimumInitialValue = kFALSE;
3657      AliUnfolding::fgMinimumInitialValueFix = -1;
3658      AliUnfolding::fgNormalizeInput = kFALSE;
3659      AliUnfolding::fgNotFoundEvents = 0;
3660      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
3661      AliUnfolding::fgBayesianSmoothing  = 1;
3662      AliUnfolding::fgBayesianIterations = 10;
3663      AliUnfolding::fgDebug = kFALSE;
3664      AliUnfolding::fgCallCount = 0;
3665      AliUnfolding::fgPowern = 5;
3666      AliUnfolding::fChi2FromFit = 0.;
3667      AliUnfolding::fPenaltyVal  = 0.;
3668      AliUnfolding::fAvgResidual = 0.;
3669      AliUnfolding::fgPrintChi2Details = 0;
3670      AliUnfolding::fgCanvas = 0;
3671      AliUnfolding::fghUnfolded = 0;
3672      AliUnfolding::fghCorrelation = 0;
3673      AliUnfolding::fghEfficiency = 0;
3674      AliUnfolding::fghMeasured = 0;
3675      AliUnfolding::SetMinuitStepSize(1.);
3676      AliUnfolding::SetMinuitPrecision(1e-6);
3677      AliUnfolding::SetMinuitMaxIterations(100000);
3678      AliUnfolding::SetMinuitStrategy(2.);
3679      AliUnfolding::SetDebug(1);
3680 }
3681 //_____________________________________________________________________________
3682 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
3683     // protect heap by adding unique qualifier to name
3684     if(!protect) return 0x0;
3685     TH1D* p = (TH1D*)protect->Clone();
3686     TString tempString(fActiveString);
3687     tempString+=suffix;
3688     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3689     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3690     if(kill) delete protect;
3691     return p;
3692 }
3693 //_____________________________________________________________________________
3694 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
3695     // protect heap by adding unique qualifier to name
3696     if(!protect) return 0x0;
3697     TH2D* p = (TH2D*)protect->Clone();
3698     TString tempString(fActiveString);
3699     tempString+=suffix;
3700     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3701     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3702     if(kill) delete protect;
3703     return p;
3704 }
3705 //_____________________________________________________________________________
3706 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
3707     // protect heap by adding unique qualifier to name
3708     if(!protect) return 0x0;
3709     TGraphErrors* p = (TGraphErrors*)protect->Clone();
3710     TString tempString(fActiveString);
3711     tempString+=suffix;
3712     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3713     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
3714     if(kill) delete protect;
3715     return p;
3716 }
3717 //_____________________________________________________________________________
3718 void AliJetFlowTools::MakeAU() {
3719     // === azimuthal unfolding ===
3720     //
3721     // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
3722     // in transverse momentum and azimuthal correlation space to extract v2 params
3723     // settings are equal to the ones used for 'Make()'
3724     //
3725     // basic steps that are followed:
3726     // 1) rebin the raw output of the jet task to the desired binnings
3727     // 2) calls the unfolding routine
3728     // 3) writes output to file
3729     // can be repeated multiple times with different configurations
3730
3731     Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
3732     Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
3733     TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3734     TH1D* dPtdPhi[8];
3735     for(Int_t i(0); i < 8; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
3736
3737     for(Int_t i(0); i < 8; i++) {
3738         // 1) manipulation of input histograms
3739         // check if the input variables are present
3740         if(!PrepareForUnfolding(low[i], up[i])) return;
3741         // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
3742         //     parts of the spectrum can end up in over or underflow bins
3743         TH1D* measuredJetSpectrumIn  = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
3744
3745         // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
3746         // the template will be used as a prior for the chi2 unfolding
3747         TH1D* measuredJetSpectrumTrueBinsIn  = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
3748
3749         // get the full response matrix from the dpt and the detector response
3750         fDetectorResponse = NormalizeTH2D(fDetectorResponse);
3751         // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
3752         // so that unfolding should return the initial spectrum
3753         if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
3754         else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
3755         else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
3756         else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
3757         // normalize each slide of the response to one
3758         NormalizeTH2D(fFullResponseIn);
3759         // resize to desired binning scheme
3760         TH2D* resizedResponseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
3761         // get the kinematic efficiency
3762         TH1D* kinematicEfficiencyIn  = resizedResponseIn->ProjectionX();
3763         kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
3764         // suppress the errors
3765         for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
3766         TH1D* jetFindingEfficiency(0x0);
3767         if(fJetFindingEff) {
3768             jetFindingEfficiency = ProtectHeap(fJetFindingEff);
3769             jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
3770             jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
3771         }
3772         // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
3773         TH1D* unfoldedJetSpectrumIn(0x0);
3774         fActiveDir->cd();                   // select active dir
3775         TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
3776         dirIn->cd();                        // select inplane subdir
3777         // select the unfolding method
3778         unfoldedJetSpectrumIn = UnfoldWrapper(
3779             measuredJetSpectrumIn,
3780             resizedResponseIn,
3781             kinematicEfficiencyIn,
3782             measuredJetSpectrumTrueBinsIn,
3783             TString("dPtdPhiUnfolding"),
3784             jetFindingEfficiency);
3785         if(i==5) {
3786             resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
3787             resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
3788             resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
3789             resizedResponseIn = ProtectHeap(resizedResponseIn);
3790             resizedResponseIn->Write();
3791             kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
3792             kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
3793             kinematicEfficiencyIn->Write();
3794             fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
3795             fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
3796             fDetectorResponse->Write();
3797             // optional histograms
3798             if(fSaveFull) {
3799                 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
3800                 fSpectrumIn->Write();
3801                 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
3802                 fDptInDist->Write();
3803                 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
3804                 fDptIn->Write();
3805                 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
3806                 fFullResponseIn->Write();
3807             }
3808         }
3809         fActiveDir->cd();
3810         fDeltaPtDeltaPhi->Write();
3811         fJetPtDeltaPhi->Write();
3812
3813         TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
3814         Double_t integralError(0);
3815         for(Int_t j(0); j < 6; j++) {
3816             // get the integrated
3817             Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
3818             dPtdPhi[j]->SetBinContent(i+1, integral);
3819             dPtdPhi[j]->SetBinError(i+1, integralError);
3820         }
3821         dud->Write();
3822         // save the current state of the unfolding object
3823         SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
3824     }
3825     TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
3826     TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3827     for(Int_t i(0); i < 6; i++) {
3828         dPtdPhi[i]->Fit(fourier, "VI");
3829         v2->SetBinContent(1+i, fourier->GetParameter(1));
3830         v2->SetBinError(1+i, fourier->GetParError(1));
3831         dPtdPhi[i]->Write();
3832     }
3833     v2->Write();
3834 }
3835 //_____________________________________________________________________________