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