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