PWG/FLOW/Tasks/AliJetFlowTools.cxx

   1 /*************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
  17 //         (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
  18 //
  19 // Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
  20 //
  21 // The task uses input from two analysis tasks:
  22 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
  23 //   used to retrieve jet spectra and delta pt distributions
  24 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
  25 //   used to construct the detector response function
  26 // and unfolds jet spectra with respect to the event plane. The user can choose
  27 // different alrogithms for unfolding which are available in (ali)root. RooUnfold
  28 // libraries must be present on the system (see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html).
  29 // A test mode is available in which the spectrum is unfolded with a generated unity response
  30 // matrix.
  31 //
  32 // The weak spot of this class is the function PrepareForUnfolding, which will read
  33 // output from two output files and expects histograms with certain names and binning.
  34 // Unfolding methods itself are general and should be able to handle any input, therefore one
  35 // can forgo the PrepareForUnfolding method, and supply necessary input information via the
  36 // SetRawInput() method
  37 //
  38 // to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
  39
  40 // root includes
  41 #include "TF1.h"
  42 #include "TH1D.h"
  43 #include "TH2D.h"
  44 #include "THStack.h"
  45 #include "TGraph.h"
  46 #include "TGraphErrors.h"
  47 #include "TCanvas.h"
  48 #include "TLegend.h"
  49 #include "TArrayD.h"
  50 #include "TList.h"
  51 #include "TMinuit.h"
  52 #include "TVirtualFitter.h"
  53 #include "TLegend.h"
  54 #include "TCanvas.h"
  55 #include "TStyle.h"
  56 #include "TLine.h"
  57 #include "TMath.h"
  58 #include "TVirtualFitter.h"
  59 #include "TFitResultPtr.h"
  60 // aliroot includes
  61 #include "AliUnfolding.h"
  62 #include "AliAnaChargedJetResponseMaker.h"
  63 // class includes
  64 #include "AliJetFlowTools.h"
  65 // roo unfold includes (make sure you have these available on your system)
  66 #include "RooUnfold.h"
  67 #include "RooUnfoldResponse.h"
  68 #include "RooUnfoldSvd.h"
  69 #include "TSVDUnfold.h"
  70
  71 using namespace std;
  72 //_____________________________________________________________________________
  73 AliJetFlowTools::AliJetFlowTools() :
  74     fResponseMaker      (new AliAnaChargedJetResponseMaker()),
  75     fPower              (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
  76     fSaveFull           (kFALSE),
  77     fActiveString       (""),
  78     fActiveDir          (0x0),
  79     fInputList          (0x0),
  80     fRefreshInput       (kTRUE),
  81     fOutputFileName     ("UnfoldedSpectra.root"),
  82     fOutputFile         (0x0),
  83     fCentralityBin      (0),
  84     fDetectorResponse   (0x0),
  85     fJetFindingEff      (0x0),
  86     fBetaIn             (.1),
  87     fBetaOut            (.1),
  88     fAvoidRoundingError (kFALSE),
  89     fUnfoldingAlgorithm (kChi2),
  90     fPrior              (kPriorMeasured),
  91     fBinsTrue           (0x0),
  92     fBinsRec            (0x0),
  93     fBinsTruePrior      (0x0),
  94     fBinsRecPrior       (0x0),
  95     fSVDRegIn           (5),
  96     fSVDRegOut          (5),
  97     fSVDToy             (kTRUE),
  98     fJetRadius          (0.3),
  99     fEventCount         (-1),
 100     fNormalizeSpectra   (kTRUE),
 101     fSmoothenSpectrum   (kTRUE),
 102     fFitMin             (60.),
 103     fFitMax             (105.),
 104     fFitStart           (75.),
 105     fTestMode           (kFALSE),
 106     fNoDphi             (kFALSE),
 107     fRawInputProvided   (kFALSE),
 108     fEventPlaneRes      (.63),
 109     fUseDetectorResponse(kTRUE),
 110     fTrainPower         (kTRUE),
 111     fRMSSpectrumIn      (0x0),
 112     fRMSSpectrumOut     (0x0),
 113     fRMSRatio           (0x0),
 114     fRMSV2              (0x0),
 115     fDeltaPtDeltaPhi    (0x0),
 116     fJetPtDeltaPhi      (0x0),
 117     fSpectrumIn         (0x0),
 118     fSpectrumOut        (0x0),
 119     fDptInDist          (0x0),
 120     fDptOutDist         (0x0),
 121     fDptIn              (0x0),
 122     fDptOut             (0x0),
 123     fFullResponseIn     (0x0),
 124     fFullResponseOut    (0x0),
 125     fUnfoldedIn         (0x0),
 126     fUnfoldedOut        (0x0) { // class constructor
 127     // create response maker weight function
 128     fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
 129 }
 130 //_____________________________________________________________________________
 131 void AliJetFlowTools::Make() {
 132     // core function of the class
 133     // 1) rebin the raw output of the jet task to the desired binnings
 134     // 2) calls the unfolding routine
 135     // 3) writes output to file
 136     // can be repeated multiple times with different configurations
 137
 138     // 1) manipulation of input histograms
 139     // check if the input variables are present
 140     if(fRefreshInput) {
 141         if(!PrepareForUnfolding()) {
 142             printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
 143             return;
 144         }
 145     }
 146     // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
 147     //     parts of the spectrum can end up in over or underflow bins
 148     TH1D* resizedJetPtIn  = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
 149     TH1D* resizedJetPtOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
 150
 151     // 1b) get the unfolding template
 152     // the template will be used as a prior for the chi2 unfolding
 153     // it holds thie rec spectrum, but is rebinned to the gen binning scheme
 154     TH1D* unfoldingTemplateIn  = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
 155     TH1D* unfoldingTemplateOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
 156
 157     // get the full response matrix from the dpt and the detector response
 158     fDetectorResponse = NormalizeTH2D(fDetectorResponse);
 159     // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
 160     // so that unfolding should return the initial spectrum
 161     if(!fTestMode) {
 162         fFullResponseIn  = (fUseDetectorResponse) ? MatrixMultiplication(fDptIn, fDetectorResponse) : fDptIn;
 163         fFullResponseOut = (fUseDetectorResponse) ? MatrixMultiplication(fDptOut, fDetectorResponse) : fDptOut;
 164     } else {
 165         fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
 166         fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
 167     }
 168     // normalize each slide of the response to one
 169     NormalizeTH2D(fFullResponseIn);
 170     NormalizeTH2D(fFullResponseOut);
 171     // resize to desired binning scheme
 172     TH2D* resizedResonseIn  = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
 173     TH2D* resizedResonseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
 174     // get the kinematic efficiency
 175     TH1D* kinematicEfficiencyIn  = resizedResonseIn->ProjectionX();
 176     kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
 177     TH1D* kinematicEfficiencyOut = resizedResonseOut->ProjectionX();
 178     kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
 179     // suppress the errors
 180     for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
 181         kinematicEfficiencyIn->SetBinError(1+i, 0.);
 182         kinematicEfficiencyOut->SetBinError(1+i, 0.);
 183     }
 184     TH1D* jetFindingEfficiency(0x0);
 185     if(fJetFindingEff) {
 186         jetFindingEfficiency = ProtectHeap(fJetFindingEff);
 187         jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
 188         jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
 189     }
 190     // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
 191     Bool_t convergedIn(kFALSE), convergedOut(kFALSE);
 192
 193     fActiveDir->cd();                   // select active dir
 194     TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
 195     dirIn->cd();                        // select inplane subdir
 196     // select the unfolding method
 197     switch (fUnfoldingAlgorithm) {
 198         case kChi2 : {
 199             convergedIn = UnfoldSpectrumChi2(       // do the inplane unfolding
 200                 resizedJetPtIn,
 201                 resizedResonseIn,
 202                 kinematicEfficiencyIn,
 203                 unfoldingTemplateIn,
 204                 fUnfoldedIn,
 205                 TString("in"),
 206                 jetFindingEfficiency);
 207             printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
 208         } break;
 209         case kSVD : {
 210             convergedIn = UnfoldSpectrumSVD(       // do the inplane unfolding
 211                 resizedJetPtIn,
 212                 resizedResonseIn,
 213                 kinematicEfficiencyIn,
 214                 unfoldingTemplateIn,
 215                 fUnfoldedIn,
 216                 TString("in"),
 217                 jetFindingEfficiency);
 218             printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
 219         } break;
 220         case kSVDlegacy : {
 221             convergedIn = UnfoldSpectrumSVDlegacy(       // do the inplane unfolding
 222                 resizedJetPtIn,
 223                 resizedResonseIn,
 224                 kinematicEfficiencyIn,
 225                 unfoldingTemplateIn,
 226                 fUnfoldedIn,
 227                 TString("in"),
 228                 jetFindingEfficiency);
 229             printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
 230         } break;
 231         case kNone : {  // do nothing, just rebin and optionally smooothen the spectrum
 232             resizedResonseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
 233             if(fSmoothenSpectrum) resizedJetPtIn = SmoothenSpectrum(resizedJetPtIn, fPower, fFitMin, fFitMin, fFitStart);
 234             fUnfoldedIn = ProtectHeap(resizedJetPtIn, kTRUE, TString("in"));
 235             convergedIn = kTRUE;
 236         } break;
 237
 238         default : {
 239             printf(" > Selected unfolding method is not implemented yet ! \n");
 240             return;
 241         }
 242     }
 243     resizedResonseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
 244     resizedResonseIn->SetXTitle("p_{T}^{true} [GeV/c]");
 245     resizedResonseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
 246     resizedResonseIn = ProtectHeap(resizedResonseIn);
 247     resizedResonseIn->Write();
 248     kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
 249     kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
 250     kinematicEfficiencyIn->Write();
 251     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 252     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 253     fDetectorResponse->Write();
 254     // optional histograms
 255     if(fSaveFull) {
 256         fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
 257         fSpectrumIn->Write();
 258         fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
 259         fDptInDist->Write();
 260         fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
 261         fDptIn->Write();
 262         fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
 263         fFullResponseIn->Write();
 264     }
 265     fActiveDir->cd();
 266     TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
 267     dirOut->cd();
 268     switch (fUnfoldingAlgorithm) {
 269         case kChi2 : {
 270             convergedOut = UnfoldSpectrumChi2(
 271                 resizedJetPtOut,
 272                 resizedResonseOut,
 273                 kinematicEfficiencyOut,
 274                 unfoldingTemplateOut,
 275                 fUnfoldedOut,
 276                 TString("out"),
 277                 jetFindingEfficiency);
 278             printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
 279         } break;
 280         case kSVD : {
 281             convergedOut = UnfoldSpectrumSVD(
 282                 resizedJetPtOut,
 283                 resizedResonseOut,
 284                 kinematicEfficiencyOut,
 285                 unfoldingTemplateOut,
 286                 fUnfoldedOut,
 287                 TString("out"),
 288                 jetFindingEfficiency);
 289             printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
 290         } break;
 291         case kSVDlegacy : {
 292             convergedOut = UnfoldSpectrumSVDlegacy(
 293                 resizedJetPtOut,
 294                 resizedResonseOut,
 295                 kinematicEfficiencyOut,
 296                 unfoldingTemplateOut,
 297                 fUnfoldedOut,
 298                 TString("out"),
 299                 jetFindingEfficiency);
 300             printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
 301         } break;
 302         case kNone : {  // do nothing, just rebin and optionally smooothen the spectrum
 303             resizedResonseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
 304             if(fSmoothenSpectrum) resizedJetPtOut = SmoothenSpectrum(resizedJetPtOut, fPower, fFitMin, fFitMin, fFitStart);
 305             fUnfoldedOut = ProtectHeap(resizedJetPtOut, kTRUE, TString("out"));
 306             convergedOut = kTRUE;
 307         } break;
 308         default : {
 309             printf(" > Selected unfolding method is not implemented yet ! \n");
 310             return;
 311         }
 312     }
 313     resizedResonseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
 314     resizedResonseOut->SetXTitle("p_{T}^{true} [GeV/c]");
 315     resizedResonseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
 316     resizedResonseOut = ProtectHeap(resizedResonseOut);
 317     resizedResonseOut->Write();
 318     kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
 319     kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
 320     kinematicEfficiencyOut->Write();
 321     fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
 322     fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
 323     fDetectorResponse->Write();
 324     if(jetFindingEfficiency) jetFindingEfficiency->Write();
 325     // optional histograms
 326     if(fSaveFull) {
 327         fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
 328         fSpectrumOut->Write();
 329         fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
 330         fDptOutDist->Write();
 331         fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
 332         fDptOut->Write();
 333         fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
 334         fFullResponseOut->Write();
 335     }
 336
 337     // write general output histograms to file
 338     fActiveDir->cd();
 339     if(convergedIn && convergedOut && fUnfoldedIn && fUnfoldedOut) {
 340         TGraphErrors* ratio(GetRatio((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_in"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_out")));
 341         if(ratio) {
 342             ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
 343             ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 344             ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 345             ratio = ProtectHeap(ratio);
 346             ratio->Write();
 347             // write histo values to RMS files if both routines converged
 348             // input values are weighted by their uncertainty
 349             for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
 350                 if(fUnfoldedIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), fUnfoldedIn->GetBinContent(i+1), 1./TMath::Power(fUnfoldedIn->GetBinError(i+1), 2.));
 351                 if(fUnfoldedOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), fUnfoldedOut->GetBinContent(i+1), 1./TMath::Power(fUnfoldedOut->GetBinError(i+1), 2.));
 352                 if(fUnfoldedOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), fUnfoldedIn->GetBinContent(i+1) / fUnfoldedOut->GetBinContent(i+1));
 353            }
 354         }
 355         TGraphErrors* v2(GetV2((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_inv2"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_outv2")));
 356         if(v2) {
 357             v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 358             v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 359             v2->GetYaxis()->SetTitle("v_{2}");
 360             v2 = ProtectHeap(v2);
 361             v2->Write();
 362         }
 363     } else if (fUnfoldedOut && fUnfoldedIn) {
 364         TGraphErrors* ratio(GetRatio((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_in"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
 365         if(ratio) {
 366             ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
 367             ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 368             ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
 369             ratio = ProtectHeap(ratio);
 370             ratio->Write();
 371         }
 372         TGraphErrors* v2(GetV2((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_inv2"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_outv2")));
 373          if(v2) {
 374             v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
 375             v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 376             v2->GetYaxis()->SetTitle("v_{2}");
 377             v2 = ProtectHeap(v2);
 378             v2->Write();
 379         }
 380     }
 381     fDeltaPtDeltaPhi->Write();
 382     fJetPtDeltaPhi->Write();
 383     SaveConfiguration(convergedIn, convergedOut);
 384 }
 385 //_____________________________________________________________________________
 386 Bool_t AliJetFlowTools::UnfoldSpectrumChi2(
 387         TH1D* resizedJetPt,             // truncated raw jets (same binning as pt rec of response)
 388         TH2D* resizedResonse,           // response matrix
 389         TH1D* kinematicEfficiency,      // kinematic efficiency
 390         TH1D* unfoldingTemplate,        // unfolding template: same binning is pt gen of response
 391         TH1D *&unfolded,                // will point to the unfolded spectrum
 392         TString suffix,                 // suffix (in or out of plane)
 393         TH1D* jetFindingEfficiency)     // jet finding efficiency (optional)
 394 {
 395     // unfold the spectrum using chi2 minimization
 396
 397     // step 0) setup the static members of AliUnfolding
 398     ResetAliUnfolding();                // reset from previous iteration
 399                                         // also deletes and re-creates the global TVirtualFitter
 400     AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
 401     if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 402     else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 403     if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
 404     else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
 405     AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
 406
 407     // step 1) clone all input histograms.
 408
 409     // resizedJetPtLocal holds the spectrum that needs to be unfolded
 410     TH1D *resizedJetPtLocal = (TH1D*)resizedJetPt->Clone(Form("resizedJetPtLocal_%s", suffix.Data()));
 411     if(fSmoothenSpectrum) resizedJetPtLocal = SmoothenSpectrum(resizedJetPtLocal, fPower, fFitMin, fFitMax, fFitStart);
 412     // unfolded local will be filled with the result of the unfolding
 413     TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
 414
 415     // full response matrix and kinematic efficiency
 416     TH2D* resizedResponseLocal = (TH2D*)resizedResonse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
 417     TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
 418
 419     // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
 420     TH1D *priorLocal = (TH1D*)unfoldingTemplate->Clone(Form("priorLocal_%s", suffix.Data()));
 421     if(fSmoothenSpectrum) priorLocal = SmoothenSpectrum(priorLocal, fPower, fFitMin, fFitMax, fFitStart);
 422
 423     // step 2) start the unfolding
 424     Int_t status(-1), i(0);
 425     while(status < 0 && i < 100) {
 426         // i > 0 means that the first iteration didn't converge. in that case, the result of the first
 427         // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
 428         if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
 429         status = AliUnfolding::Unfold(
 430                 resizedResponseLocal,           // response matrix
 431                 kinematicEfficiencyLocal,       // efficiency applied on the unfolded spectrum (can be NULL)
 432                 resizedJetPtLocal,              // measured spectrum
 433                 priorLocal,                     // initial conditions (set NULL to use measured spectrum)
 434                 unfoldedLocal);                 // results
 435         // status holds the minuit fit status (where 0 means convergence)
 436         i++;
 437     }
 438     // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
 439     if(status == 0 && gMinuit->fISW[1] == 3) {
 440         // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
 441         TVirtualFitter *fitter(TVirtualFitter::GetFitter());
 442         if(gMinuit) gMinuit->Command("SET COV");
 443         TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
 444         TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
 445         pearson->Print();
 446         TH2D *hPearson(new TH2D(*pearson));
 447         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
 448         hPearson = ProtectHeap(hPearson);
 449         hPearson->Write();
 450     } else status = -1;
 451
 452     // step 3) refold the unfolded spectrum and save the ratio measured / refolded
 453     TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
 454     foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
 455     unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
 456     TGraphErrors* ratio(GetRatio(foldedLocal, resizedJetPtLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
 457     if(ratio) {
 458         ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
 459         ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 460         ratio = ProtectHeap(ratio);
 461         ratio->Write();
 462     }
 463
 464     // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
 465     // objects are cloned using 'ProtectHeap()'
 466     resizedJetPtLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
 467     resizedJetPtLocal = ProtectHeap(resizedJetPtLocal);
 468     resizedJetPtLocal->Write();
 469
 470     resizedResponseLocal = ProtectHeap(resizedResponseLocal);
 471     resizedResponseLocal->Write();
 472
 473     unfoldedLocal = ProtectHeap(unfoldedLocal);
 474     if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
 475     unfoldedLocal->Write();
 476
 477     foldedLocal = ProtectHeap(foldedLocal);
 478     foldedLocal->Write();
 479
 480     priorLocal = ProtectHeap(priorLocal);
 481     priorLocal->Write();
 482
 483     // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
 484     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));
 485     fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
 486     fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
 487     fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
 488     fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
 489     fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
 490     fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
 491     fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
 492     fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
 493     fitStatus->Write();
 494     unfolded = unfoldedLocal;
 495     return (status == 0) ? kTRUE : kFALSE;
 496 }
 497 //_____________________________________________________________________________
 498 Bool_t AliJetFlowTools::UnfoldSpectrumSVDlegacy(
 499         TH1D* resizedJetPt,                     // jet pt in pt rec bins
 500         TH2D* resizedResonse,                   // full response matrix, normalized in slides of pt true
 501         TH1D* kinematicEfficiency,              // kinematic efficiency
 502         TH1D* unfoldingTemplate,                // jet pt in pt true bins, also the prior when measured is chosen as prior
 503         TH1D *&unfolded,                        // will point to result. temporarily holds prior when chi2 is chosen as prior
 504         TString suffix,                         // suffix (in, out)
 505         TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 506 {
 507     // use SVD (singular value decomposition) method to unfold spectra
 508
 509     // 1) get a prior for unfolding.
 510     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
 511     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
 512     dirOut->cd();
 513     switch (fPrior) {    // select the prior for unfolding
 514         case kPriorChi2 : {
 515             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
 516                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original (SVD) binning
 517                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
 518                 TArrayD* tempArrayRec(fBinsRec);
 519                 fBinsRec = fBinsRecPrior;
 520                 TH1D* resizedJetPtChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
 521                 TH1D* unfoldingTemplateChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
 522                 TH2D* resizedResonseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
 523                 TH1D* kinematicEfficiencyChi2(resizedResonseChi2->ProjectionX());
 524                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
 525                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
 526                 if(! UnfoldSpectrumChi2(
 527                             resizedJetPtChi2,
 528                             resizedResonseChi2,
 529                             kinematicEfficiencyChi2,
 530                             unfoldingTemplateChi2,  // prior for chi2 unfolding (measured)
 531                             unfolded,               // will hold the result from chi2 (is prior for SVD)
 532                             TString(Form("prior_%s", suffix.Data()))) ) {
 533                     printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
 534                     printf("               probably Chi2 unfolding did not converge < \n");
 535                     return kFALSE;
 536                 }
 537                 fBinsTrue = tempArrayTrue;  // reset bins borders
 538                 fBinsRec = tempArrayRec;
 539                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
 540             } else if(! UnfoldSpectrumChi2(
 541                         resizedJetPt,
 542                         resizedResonse,
 543                         kinematicEfficiency,
 544                         unfoldingTemplate,      // prior for chi2 unfolding (measured)
 545                         unfolded,               // will hold the result from chi2 (is prior for SVD)
 546                         TString(Form("prior_%s", suffix.Data()))) ) {
 547                 printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
 548                 printf("               probably Chi2 unfolding did not converge < \n");
 549                 return kFALSE;
 550             }
 551             if(!unfolded) {
 552                 printf(" > UnfoldSVD:: panic, Chi2 unfolding converged but the prior is NULL ! < " );
 553                 return kFALSE;
 554             }
 555             break;
 556         }
 557         case kPriorMeasured : {
 558             unfolded = (TH1D*)unfoldingTemplate->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
 559             if(fSmoothenSpectrum) {     // optionally smoothen the measured prior
 560                 unfolded->Sumw2();
 561                 TFitResultPtr r = unfolded->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
 562                 if((int)r == 0) {
 563                     for(Int_t i(1); i < unfolded->GetNbinsX() + 1; i++) {
 564                         if(unfolded->GetBinCenter(i) > fFitStart) {     // from this pt value use extrapolation
 565                             unfolded->SetBinContent(i,fPower->Integral(unfolded->GetXaxis()->GetBinLowEdge(i),unfolded->GetXaxis()->GetBinUpEdge(i))/unfolded->GetXaxis()->GetBinWidth(i));
 566                         }
 567                     }
 568                 }else printf(" > PANIC, SMOOTHENING FAILED < \n");
 569             }
 570         }
 571         default : break;
 572     }
 573     // note: true and measured spectrum must have same binning for SVD unfolding
 574     // a sane starting point for regularization is nbins / 2 (but the user has to set this ! )
 575     if(unfoldingTemplate->GetXaxis()->GetNbins() != resizedJetPt->GetXaxis()->GetNbins()) {
 576         printf(" > UnfoldSpectrumSVD:: PANIC, true and measured spectrum must have same numer of bins ! < \n ");
 577     }
 578     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 579     cout << " 1) retrieved prior " << endl;
 580
 581     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 582     // prior
 583     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
 584     TH1D *unfoldedLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
 585     // raw jets in pt rec binning
 586     TH1D *cachedRawJetLocal((TH1D*)resizedJetPt->Clone(Form("jets_%s", suffix.Data())));
 587     // raw jets in pt true binning
 588     TH1D *cachedRawJetLocalCoarse((TH1D*)unfoldingTemplate->Clone(Form("unfoldingTemplate_%s", suffix.Data())));
 589     // copy of raw jets in pt true binning
 590     TH1D *cachedRawJetLocalCoarseOrig((TH1D*)cachedRawJetLocalCoarse->Clone(Form("cachedRawJetLocalCoarseOrig_%s", suffix.Data())));
 591     // local copies response matrix
 592     TH2D *cachedResponseLocal((TH2D*)resizedResonse->Clone(Form("cachedResponseLocal_%s", suffix.Data())));
 593     // local copy of response matrix, all true slides normalized to 1 (correction for the efficiency)
 594     TH2D *cachedResponseLocalNorm((TH2D*)resizedResonse->Clone(Form("cachedResponseLocalNorm_%s", suffix.Data())));
 595     cachedResponseLocalNorm = NormalizeTH2D(cachedResponseLocalNorm);
 596     // kinematic efficiency
 597     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 598     // place holder histos
 599     TH1D *unfoldedLocalSVD(0x0);
 600     TH1D *foldedLocalSVD(0x0);
 601     cout << " 2) setup necessary input " << endl;
 602     // 3) configure routine
 603     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 604     // prior: use fit for where the histogram is sparsely filled
 605     if(fSmoothenSpectrum) cachedRawJetLocalCoarse = SmoothenSpectrum(cachedRawJetLocalCoarse, fPower, fFitMin, fFitMax, fFitStart);
 606     if(fSmoothenSpectrum) cachedRawJetLocal = SmoothenSpectrum(cachedRawJetLocal, fPower, fFitMin, fFitMax, fFitStart);
 607     if(fSmoothenSpectrum) unfoldedLocal = SmoothenSpectrum(unfoldedLocal, fPower, fFitMin, fFitMax, fFitStart);
 608     cout << " step 3) configured routine " << endl;
 609
 610     // 4) get transpose matrices
 611     // a) get the transpose matrix for the prior
 612     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(cachedResponseLocal));
 613     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 614     // normalize it with the prior
 615     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, unfoldedLocal);
 616     cout << " 4a) retrieved first transpose matrix " << endl;
 617     // b) prior norm
 618     TH2* responseMatrixLocalTransposePriorNorm(fResponseMaker->GetTransposeResponsMatrix(cachedResponseLocalNorm));
 619     responseMatrixLocalTransposePriorNorm->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePriorNorm->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePriorNorm->GetName(), suffix.Data()));
 620     // normalize with the prior
 621     responseMatrixLocalTransposePriorNorm = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePriorNorm, unfoldedLocal);
 622     cout << " 4b) retrieved second transpose matrix " << endl;
 623
 624     // 5) get response for SVD unfolding
 625     RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
 626
 627     // change to inplane dir
 628     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 629
 630     cout << " 5) retrieved roo unfold response object " << endl;
 631     // 6) actualy unfolding loop
 632     RooUnfoldSvd unfoldSVD(&responseSVD, cachedRawJetLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 633     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
 634     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
 635     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 636     cout << " Pearson coeffs" << endl;
 637     // create the unfolding qa plots
 638     cout << " 6) unfolded spectrum " << endl;
 639     if(pearson) {
 640         TH2D* hPearson = new TH2D(*pearson);
 641         pearson->Print();
 642         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 643         hPearson = ProtectHeap(hPearson);
 644         hPearson->Write();
 645     } else return kFALSE;       // return if unfolding didn't converge
 646     // correct for the efficiency
 647     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
 648
 649     // plot singular values and d_i vector
 650     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
 651     TH1* hSVal(svdUnfold->GetSV());
 652     TH1D* hdi(svdUnfold->GetD());
 653     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
 654     hSVal->SetXTitle("singular values");
 655     hSVal->Write();
 656     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
 657     hdi->SetXTitle("|d_{i}^{kreg}|");
 658     hdi->Write();
 659     cout << " plotted singular values and d_i vector " << endl;
 660
 661     // 7) refold the unfolded spectrum
 662     foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, cachedResponseLocalNorm,kinematicEfficiencyLocal);
 663     TGraphErrors* ratio(GetRatio(cachedRawJetLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
 664     ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
 665     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 666     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 667     ratio->Write();
 668     cout << " 7) refolded the unfolded spectrum " << endl;
 669
 670     // write to output
 671     cachedRawJetLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 672     cachedRawJetLocal = ProtectHeap(cachedRawJetLocal);
 673     cachedRawJetLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 674     cachedRawJetLocal->Write(); // input spectrum
 675     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 676     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
 677     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
 678     unfoldedLocalSVD->Write();  // unfolded spectrum
 679     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 680     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
 681     foldedLocalSVD->Write();    // re-folded spectrum
 682
 683     // switch back to active root directory
 684     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 685     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix");
 686     responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
 687     responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
 688     responseMatrixLocalTransposePrior->Write();
 689     responseMatrixLocalTransposePriorNorm->SetNameTitle("TransposeResponseMatrixNorm", "Transpose of response matrix normalized with prior");
 690     responseMatrixLocalTransposePriorNorm->SetXTitle("p_{T}^{true} [GeV/c]");
 691     responseMatrixLocalTransposePriorNorm->SetYTitle("p_{T}^{rec} [GeV/c]");
 692     responseMatrixLocalTransposePriorNorm->Write();
 693     cachedRawJetLocal->SetNameTitle("PriorOriginal", "Prior, original");
 694     cachedRawJetLocal->SetXTitle("p_{t} [GeV/c]");
 695     cachedRawJetLocalCoarse->SetNameTitle("PriorSmoothened", "Prior, smoothened");
 696     cachedRawJetLocalCoarse->SetXTitle("p_{t} [GeV/c]");
 697     cachedRawJetLocalCoarse->Write();
 698     cachedRawJetLocalCoarseOrig->SetNameTitle("Prior", "Prior");
 699     cachedRawJetLocalCoarseOrig->SetXTitle("p_{t} [GeV/c]");
 700     cachedRawJetLocalCoarseOrig->Write();
 701     unfolded = unfoldedLocalSVD;
 702     cachedResponseLocalNorm = ProtectHeap(cachedResponseLocalNorm);
 703     cachedResponseLocalNorm->Write();
 704
 705     // save some info
 706     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));
 707     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 708     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
 709     fitStatus->Write();
 710
 711     return (unfoldedLocalSVD) ? kTRUE : kFALSE;
 712 }
 713 //_____________________________________________________________________________
 714 Bool_t AliJetFlowTools::UnfoldSpectrumSVD(
 715         TH1D* resizedJetPt,                     // jet pt in pt rec bins
 716         TH2D* resizedResonse,                   // full response matrix, normalized in slides of pt true
 717         TH1D* kinematicEfficiency,              // kinematic efficiency
 718         TH1D* unfoldingTemplate,                // jet pt in pt true bins, also the prior when measured is chosen as prior
 719         TH1D *&unfolded,                        // will point to result. temporarily holds prior when chi2 is chosen as prior
 720         TString suffix,                         // suffix (in, out)
 721         TH1D* jetFindingEfficiency)             // jet finding efficiency (optional)
 722 {
 723     // use SVD (singular value decomposition) method to unfold spectra
 724
 725     // 1) get a prior for unfolding.
 726     // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
 727     TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
 728     dirOut->cd();
 729     switch (fPrior) {    // select the prior for unfolding
 730         case kPriorChi2 : {
 731             if(fBinsTruePrior && fBinsRecPrior) {       // if set, use different binning for the prior
 732                 TArrayD* tempArrayTrue(fBinsTrue);      // temporarily cache the original (SVD) binning
 733                 fBinsTrue = fBinsTruePrior;             // switch binning schemes (will be used in UnfoldSpectrumChi2())
 734                 TArrayD* tempArrayRec(fBinsRec);
 735                 fBinsRec = fBinsRecPrior;
 736                 TH1D* resizedJetPtChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
 737                 TH1D* unfoldingTemplateChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
 738                 TH2D* resizedResonseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
 739                 TH1D* kinematicEfficiencyChi2(resizedResonseChi2->ProjectionX());
 740                 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
 741                 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
 742                 if(! UnfoldSpectrumChi2(
 743                             resizedJetPtChi2,
 744                             resizedResonseChi2,
 745                             kinematicEfficiencyChi2,
 746                             unfoldingTemplateChi2,  // prior for chi2 unfolding (measured)
 747                             unfolded,               // will hold the result from chi2 (is prior for SVD)
 748                             TString(Form("prior_%s", suffix.Data()))) ) {
 749                     printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
 750                     printf("               probably Chi2 unfolding did not converge < \n");
 751                     return kFALSE;
 752                 }
 753                 fBinsTrue = tempArrayTrue;  // reset bins borders
 754                 fBinsRec = tempArrayRec;
 755                 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data())));     // rebin unfolded
 756             } else if(! UnfoldSpectrumChi2(
 757                         resizedJetPt,
 758                         resizedResonse,
 759                         kinematicEfficiency,
 760                         unfoldingTemplate,      // prior for chi2 unfolding (measured)
 761                         unfolded,               // will hold the result from chi2 (is prior for SVD)
 762                         TString(Form("prior_%s", suffix.Data()))) ) {
 763                 printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
 764                 printf("               probably Chi2 unfolding did not converge < \n");
 765                 return kFALSE;
 766             }
 767             if(!unfolded) {
 768                 printf(" > UnfoldSVD:: panic, Chi2 unfolding converged but the prior is NULL ! < " );
 769                 return kFALSE;
 770             }
 771             break;
 772         }
 773         case kPriorMeasured : {
 774             unfolded = (TH1D*)unfoldingTemplate->Clone(Form("kPriorMeasured_%s", suffix.Data()));       // copy template to unfolded to use as prior
 775             if(fSmoothenSpectrum) unfolded = SmoothenSpectrum(unfolded, fPower, fFitMin, fFitMax, fFitStart);
 776         }
 777         default : break;
 778     }
 779     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 780     cout << " 1) retrieved prior " << endl;
 781
 782     // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
 783     // prior
 784     if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
 785     TH1D *unfoldedLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
 786     // raw jets in pt rec binning
 787     TH1D *cachedRawJetLocal((TH1D*)resizedJetPt->Clone(Form("jets_%s", suffix.Data())));
 788     // raw jets in pt true binning
 789     TH1D *cachedRawJetLocalCoarse((TH1D*)unfoldingTemplate->Clone(Form("unfoldingTemplate_%s", suffix.Data())));
 790     // copy of raw jets in pt true binning
 791     TH1D *cachedRawJetLocalCoarseOrig((TH1D*)cachedRawJetLocalCoarse->Clone(Form("cachedRawJetLocalCoarseOrig_%s", suffix.Data())));
 792     // local copies response matrix
 793     TH2D *cachedResponseLocal((TH2D*)resizedResonse->Clone(Form("cachedResponseLocal_%s", suffix.Data())));
 794     // kinematic efficiency
 795     TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
 796     // place holder histos
 797     TH1D *unfoldedLocalSVD(0x0);
 798     TH1D *foldedLocalSVD(0x0);
 799     cout << " 2) setup necessary input " << endl;
 800     // 3) configure routine
 801     RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
 802     // prior: use fit for where the histogram is sparsely filled
 803     if(fSmoothenSpectrum) cachedRawJetLocalCoarse = SmoothenSpectrum(cachedRawJetLocalCoarse, fPower, fFitMin, fFitMax, fFitStart);
 804     if(fSmoothenSpectrum) cachedRawJetLocal = SmoothenSpectrum(cachedRawJetLocal, fPower, fFitMin, fFitMax, fFitStart);
 805     if(fSmoothenSpectrum) unfoldedLocal = SmoothenSpectrum(unfoldedLocal, fPower, fFitMin, fFitMax, fFitStart);
 806     cout << " 3) configured routine " << endl;
 807
 808     // 4) get transpose matrices, where the y-axis corresponds to the true binning
 809     // and the x-axis to the measured binning
 810     TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(cachedResponseLocal));
 811     responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
 812     // normalize the transpose matrix with the prior in the y-direction (truth)
 813     TH1D* tempUnfoldedLocal = static_cast<TH1D*>(unfoldedLocal->Clone("temp"));
 814     tempUnfoldedLocal->Multiply(kinematicEfficiency);
 815     responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, tempUnfoldedLocal);
 816     delete tempUnfoldedLocal;
 817
 818     // get the jet spectrum response matrix in the form of a RooUnfoldResponse object
 819     RooUnfoldResponse responseSVD(0, unfoldedLocal, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
 820
 821     // change to inplane dir
 822     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 823     cout << " 5) retrieved roo unfold response object " << endl;
 824
 825     RooUnfoldSvd unfoldSVD(&responseSVD, cachedRawJetLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 826     unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
 827
 828     TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
 829     TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
 830     cout << " Pearson coeffs" << endl;
 831     // create the unfolding qa plots
 832     cout << " 6) unfolded spectrum " << endl;
 833     if(pearson) {
 834         TH2D* hPearson = new TH2D(*pearson);
 835         pearson->Print();
 836         hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
 837         hPearson = ProtectHeap(hPearson);
 838         hPearson->Write();
 839     } else return kFALSE;       // return if unfolding didn't converge
 840     // correct for the efficiency
 841     unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
 842     // plot singular values and d_i vector
 843     TSVDUnfold* svdUnfold(unfoldSVD.Impl());
 844     TH1* hSVal(svdUnfold->GetSV());
 845     TH1D* hdi(svdUnfold->GetD());
 846     hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
 847     hSVal->SetXTitle("singular values");
 848     hSVal->Write();
 849     hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
 850     hdi->SetXTitle("|d_{i}^{kreg}|");
 851     hdi->Write();
 852     cout << " plotted singular values and d_i vector " << endl;
 853
 854     // 7 refold the unfolded spectrum with the RooUnfold object
 855     TH1D* unfolded_eff = static_cast<TH1D*>(unfoldedLocalSVD->Clone("unfolded_eff"));
 856     unfolded_eff->Multiply(kinematicEfficiencyLocal);
 857     RooUnfoldResponse rooRefold(0, 0, responseMatrixLocalTransposePrior, Form("rooRefold_%s", suffix.Data()), Form("rooRefold_%s", suffix.Data()));
 858     foldedLocalSVD = static_cast<TH1D*>(rooRefold.ApplyToTruth(unfolded_eff, "refolded"));
 859     delete unfolded_eff;
 860     TGraphErrors* ratio(GetRatio(cachedRawJetLocal, foldedLocalSVD, "ratio  measured / re-folded", kTRUE));
 861     ratio->SetName(Form("RatioRefoldedMeasured_%s", fActiveString.Data()));
 862     ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
 863     ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
 864     ratio->Write();
 865     cout << " 7) refolded the unfolded spectrum " << endl;
 866
 867     // write to output
 868     cachedRawJetLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
 869     cachedRawJetLocal = ProtectHeap(cachedRawJetLocal);
 870     cachedRawJetLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
 871     cachedRawJetLocal->Write(); // input spectrum
 872     unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
 873     unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
 874     if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
 875     unfoldedLocalSVD->Write();  // unfolded spectrum
 876     foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
 877     foldedLocalSVD = ProtectHeap(foldedLocalSVD);
 878     foldedLocalSVD->Write();    // re-folded spectrum
 879    // switch back to active root directory
 880     (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
 881     responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix");
 882     responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
 883     responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
 884     responseMatrixLocalTransposePrior->Write();
 885     cachedRawJetLocal->SetNameTitle("PriorOriginal", "Prior, original");
 886     cachedRawJetLocal->SetXTitle("p_{t} [GeV/c]");
 887     cachedRawJetLocalCoarse->SetNameTitle("PriorSmoothened", "Prior, smoothened");
 888     cachedRawJetLocalCoarse->SetXTitle("p_{t} [GeV/c]");
 889     cachedRawJetLocalCoarse->Write();
 890     cachedRawJetLocalCoarseOrig->SetNameTitle("Prior", "Prior");
 891     cachedRawJetLocalCoarseOrig->SetXTitle("p_{t} [GeV/c]");
 892     cachedRawJetLocalCoarseOrig->Write();
 893     unfolded = unfoldedLocalSVD;
 894
 895     // save some info
 896     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));
 897     fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
 898     fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
 899     fitStatus->Write();
 900
 901     return (unfoldedLocalSVD) ? kTRUE : kFALSE;
 902 }
 903 //_____________________________________________________________________________
 904 Bool_t AliJetFlowTools::PrepareForUnfolding()
 905 {
 906     // prepare for unfolding
 907     if(fRawInputProvided) return kTRUE;
 908     if(!fInputList) {
 909         printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
 910         return kFALSE;
 911     }
 912     if(!fDetectorResponse) {
 913         printf(" AliJetFlowTools::PrepareForUnfolding() fDetectorResponse not found \n - Set detector response using AliJetFlowTools::SetDetectorResponse() \n ");
 914         return kFALSE;
 915     }
 916     // check if the pt bin for true and rec have been set
 917     if(!fBinsTrue || !fBinsRec) {
 918         printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
 919         return kFALSE;
 920     }
 921     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
 922                           // procedures, these profiles will be nonsensical, user is responsible
 923         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 924         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 925         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
 926     }
 927     if(!fTrainPower) {
 928         for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
 929     }
 930     // extract the spectra
 931     TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
 932     fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
 933     if(!fJetPtDeltaPhi) {
 934         printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
 935         return kFALSE;
 936     }
 937     fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
 938     // in plane spectrum
 939     if(fNoDphi) {
 940         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40);
 941         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40);
 942     } else {
 943         fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10);
 944         fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40));
 945         fSpectrumIn = ProtectHeap(fSpectrumIn);
 946         // out of plane spectrum
 947         fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30);
 948         fSpectrumOut = ProtectHeap(fSpectrumOut);
 949     }
 950     // normalize spectra to event count if requested
 951     if(fNormalizeSpectra) {
 952         TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
 953         if(!rho) return 0x0;
 954         Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
 955         if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
 956         if(fEventCount > 0) {
 957             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
 958             fSpectrumOut->Sumw2();
 959             Double_t pt(0);
 960             Double_t error(0); // lots of issues with the errors here ...
 961             for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
 962                 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount;       // normalized count
 963                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
 964                 fSpectrumIn->SetBinContent(1+i, pt);
 965                 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
 966                 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
 967                 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
 968             }
 969             for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
 970                 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount;       // normalized count
 971                 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
 972                 fSpectrumOut->SetBinContent(1+i, pt);
 973                 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
 974                 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
 975                 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
 976             }
 977         }
 978         if(normalizeToFullSpectrum) fEventCount = -1;
 979     }
 980     // extract the delta pt matrices
 981     TString deltaptName(Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin));
 982     fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
 983     if(!fDeltaPtDeltaPhi) {
 984         printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
 985     }
 986     fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
 987     // in plane delta pt distribution
 988     if(fNoDphi) {
 989         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40);
 990         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40);
 991     } else {
 992         fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10);
 993         fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40));
 994         // out of plane delta pt distribution
 995         fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30);
 996         fDptInDist = ProtectHeap(fDptInDist);
 997         fDptOutDist = ProtectHeap(fDptOutDist);
 998         // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
 999     }
1000
1001     // create a rec - true smeared response matrix
1002     TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1003     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1004         Bool_t skip = kFALSE;
1005         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1006             (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1007             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1008         }
1009     }
1010     TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
1011     for(Int_t j(-50); j < 250; j++) {   // loop on pt true slices j
1012         Bool_t skip = kFALSE;
1013         for(Int_t k(-50); k < 250; k++) {       // loop on pt gen slices k
1014             (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
1015             if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1016         }
1017     }
1018     fDptIn = new TH2D(*rfIn);
1019     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1020     fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1021     fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1022     fDptIn = ProtectHeap(fDptIn);
1023     fDptOut = new TH2D(*rfOut);
1024     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1025     fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1026     fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1027     fDptOut = ProtectHeap(fDptOut);
1028
1029     fRefreshInput = kTRUE;     // force cloning of the input
1030     return kTRUE;
1031 }
1032 //_____________________________________________________________________________
1033 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1034     // resize the x-axis of a th1d
1035     if(!histo) {
1036         printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1037         return NULL;
1038     }
1039     // see how many bins we need to copy
1040     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);
1041     // low is the bin number of the first new bin
1042     Int_t l = histo->GetXaxis()->FindBin(low);
1043     // set the values
1044     Double_t _x(0), _xx(0);
1045     for(Int_t i(0); i < up-low; i++) {
1046         _x = histo->GetBinContent(l+i);
1047         _xx=histo->GetBinError(l+i);
1048         resized->SetBinContent(i+1, _x);
1049         resized->SetBinError(i+1, _xx);
1050     }
1051     return resized;
1052 }
1053 //_____________________________________________________________________________
1054 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1055     // resize the y-axis of a th2d
1056     if(!histo) {
1057         printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1058         return NULL;
1059     }
1060     // see how many bins we need to copy
1061     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());
1062     // assume only the y-axis has changed
1063     // low is the bin number of the first new bin
1064     Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1065     // set the values
1066     Double_t _x(0), _xx(0);
1067     for(Int_t i(0); i < x->GetSize(); i++) {
1068         for(Int_t j(0); j < y->GetSize(); j++) {
1069             _x = histo->GetBinContent(i, low+j);
1070             _xx=histo->GetBinError(i, low+1+j);
1071             resized->SetBinContent(i, j, _x);
1072             resized->SetBinError(i, j, _xx);
1073         }
1074     }
1075     return resized;
1076 }
1077 //_____________________________________________________________________________
1078 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo) {
1079     // general method to normalize all vertical slices of a th2 to unity
1080     // i.e. get a probability matrix
1081     if(!histo) {
1082         printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1083         return NULL;
1084     }
1085     Int_t binsX = histo->GetXaxis()->GetNbins();
1086     Int_t binsY = histo->GetYaxis()->GetNbins();
1087
1088     // normalize all slices in x
1089     for(Int_t i(0); i < binsX; i++) {   // for each vertical slice
1090         Double_t weight = 0;
1091         for(Int_t j(0); j < binsY; j++) {       // loop over all the horizontal components
1092             weight+=histo->GetBinContent(i+1, j+1);
1093         }       // now we know the total weight
1094         for(Int_t j(0); j < binsY; j++) {
1095             if (weight <= 0 ) continue;
1096             histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1097             histo->SetBinError(  1+i, j+1, histo->GetBinError(  1+i, j+1)/weight);
1098         }
1099     }
1100     return histo;
1101 }
1102 //_____________________________________________________________________________
1103 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1104     // return a TH1D with the supplied histogram rebinned to the supplied bins
1105     // the returned histogram is new, the original is deleted from the heap if kill is true
1106     if(!histo || !bins) {
1107         printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1108         return NULL;
1109     }
1110     // create the output histo
1111     TString name = histo->GetName();
1112     name+="_template";
1113     name+=suffix;
1114     TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1115     for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1116         // loop over the bins of the old histo and fill the new one with its data
1117         rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1118     }
1119     if(kill) delete histo;
1120     return rebinned;
1121 }
1122 //_____________________________________________________________________________
1123 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1124     // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1125     if(!fResponseMaker || !binsTrue || !binsRec) {
1126         printf(" > RebinTH2D:: function called with NULL arguments < \n");
1127         return 0x0;
1128     }
1129     TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1130     return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1131 }
1132 //_____________________________________________________________________________
1133 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1134 {
1135     // multiply two matrices
1136     if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1137     TH2D* c = (TH2D*)a->Clone("c");
1138     for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1139         for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1140             Double_t val = 0;
1141             for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1142                 Int_t y2 = x1;
1143                 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1144             }
1145             c->SetBinContent(x2, y1, val);
1146         }
1147     }
1148     if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1149     return c;
1150 }
1151 //_____________________________________________________________________________
1152 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1153 {
1154     // normalize a th1d to a certain scale
1155     histo->Sumw2();
1156     Double_t integral = histo->Integral()*scale;
1157     if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1158     else if (scale != 1.) histo->Scale(1./scale, "width");
1159     else printf(" > Histogram integral < 0, cannot normalize \n");
1160     return histo;
1161 }
1162 //_____________________________________________________________________________
1163 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1164 {
1165     // Calculate pearson coefficients from covariance matrix
1166     TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1167     Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1168     Double_t pearson(0.);
1169     if(nrows==0 && ncols==0) return 0x0;
1170     for(Int_t row = 0; row < nrows; row++) {
1171         for(Int_t col = 0; col<ncols; col++) {
1172         if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1173         (*pearsonCoefficients)(row,col) = pearson;
1174         }
1175     }
1176     return pearsonCoefficients;
1177 }
1178 //_____________________________________________________________________________
1179 TH1D* AliJetFlowTools::SmoothenSpectrum(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1180     // smoothen the spectrum using a user defined function
1181     // returns a clone of the original spectrum if fitting failed
1182     // if kill is kTRUE the input spectrum will be deleted from the heap
1183     // if 'count' is selected, bins are filled with integers (necessary if the
1184     // histogram is interpreted in a routine which accepts only counts)
1185     TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1186     temp->Sumw2();      // if already called on the original, this will give off a warning but do nothing
1187     TFitResultPtr r = temp->Fit(function, "QWILS", "", min, max);
1188     if((int)r == 0) {   // MINUIT status
1189         for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1190             if(temp->GetBinCenter(i) > start) {     // from this pt value use extrapolation
1191                 if(counts) temp->SetBinContent(i, (int)function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1192                 else temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1193                 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1194             }
1195         }
1196     }
1197     if(kill) delete spectrum;
1198     return temp;
1199 }
1200 //_____________________________________________________________________________
1201 void AliJetFlowTools::Style(TCanvas* c, TString style)
1202 {
1203     // set a default style for a canvas
1204     if(!strcmp(style.Data(), "PEARSON")) {
1205         printf(" > style PEARSON canvas < \n");
1206         gStyle->SetOptStat(0);
1207         c->SetGridx();
1208         c->SetGridy();
1209         c->SetTicks();
1210         return;
1211     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1212         printf(" > style SPECTRUM canvas < \n");
1213         gStyle->SetOptStat(0);
1214         c->SetLogy();
1215         c->SetGridx();
1216         c->SetGridy();
1217         c->SetTicks();
1218         return;
1219     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1220 }
1221 //_____________________________________________________________________________
1222 void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1223 {
1224     // set a default style for a canvas
1225     if(!strcmp(style.Data(), "PEARSON")) {
1226         printf(" > style PEARSON pad < \n");
1227         gStyle->SetOptStat(0);
1228         c->SetGridx();
1229         c->SetGridy();
1230         c->SetTicks();
1231         return;
1232     } else if(!strcmp(style.Data(), "SPECTRUM")) {
1233         printf(" > style SPECTRUM pad < \n");
1234         gStyle->SetOptStat(0);
1235         c->SetLogy();
1236         c->SetGridx();
1237         c->SetGridy();
1238         c->SetTicks();
1239         return;
1240     } else printf(" > Style called with unknown option %s \n    returning < \n", style.Data());
1241 }
1242 //_____________________________________________________________________________
1243 void AliJetFlowTools::Style(TLegend* l)
1244 {
1245     // set a default style for a legend
1246     l->SetTextSize(.06);
1247     l->SetFillColor(0);
1248 }
1249 //_____________________________________________________________________________
1250 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1251 {
1252     // style a histo
1253     h->SetLineColor(col);
1254     h->SetMarkerColor(col);
1255     h->SetLineWidth(2.);
1256     h->SetMarkerSize(2.);
1257     h->SetTitleSize(0.06);
1258     h->GetYaxis()->SetLabelSize(0.06);
1259     h->GetXaxis()->SetLabelSize(0.06);
1260     h->GetYaxis()->SetTitleSize(0.06);
1261     h->GetXaxis()->SetTitleSize(0.06);
1262     h->GetYaxis()->SetTitleOffset(1.);
1263     h->GetXaxis()->SetTitleOffset(.9);
1264     switch (type) {
1265         case kInPlaneSpectrum : {
1266             h->SetTitle("IN PLANE");
1267             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1268             h->GetYaxis()->SetTitle("counts");
1269         } break;
1270         case kOutPlaneSpectrum : {
1271             h->SetTitle("OUT OF PLANE");
1272             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1273             h->GetYaxis()->SetTitle("counts");
1274        } break;
1275        case kUnfoldedSpectrum : {
1276             h->SetTitle("UNFOLDED");
1277             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1278             h->GetYaxis()->SetTitle("counts");
1279        } break;
1280        case kFoldedSpectrum : {
1281             h->SetTitle("FOLDED");
1282             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1283             h->GetYaxis()->SetTitle("counts");
1284        } break;
1285        case kMeasuredSpectrum : {
1286             h->SetTitle("MEASURED");
1287             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1288             h->GetYaxis()->SetTitle("counts");
1289        } break;
1290        default : break;
1291     }
1292 }
1293 //_____________________________________________________________________________
1294 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1295 {
1296     // style a histo
1297     h->SetLineColor(col);
1298     h->SetMarkerColor(col);
1299     h->SetLineWidth(2.);
1300     h->SetMarkerSize(2.);
1301     h->GetYaxis()->SetLabelSize(0.06);
1302     h->GetXaxis()->SetLabelSize(0.06);
1303     h->GetYaxis()->SetTitleSize(0.06);
1304     h->GetXaxis()->SetTitleSize(0.06);
1305     h->GetYaxis()->SetTitleOffset(1.);
1306     h->GetXaxis()->SetTitleOffset(.9);
1307     switch (type) {
1308         case kInPlaneSpectrum : {
1309             h->SetTitle("IN PLANE");
1310             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1311             h->GetYaxis()->SetTitle("counts");
1312         } break;
1313         case kOutPlaneSpectrum : {
1314             h->SetTitle("OUT OF PLANE");
1315             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1316             h->GetYaxis()->SetTitle("counts");
1317        } break;
1318        case kUnfoldedSpectrum : {
1319             h->SetTitle("UNFOLDED");
1320             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1321             h->GetYaxis()->SetTitle("counts");
1322        } break;
1323        case kFoldedSpectrum : {
1324             h->SetTitle("FOLDED");
1325             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1326             h->GetYaxis()->SetTitle("counts");
1327        } break;
1328        case kMeasuredSpectrum : {
1329             h->SetTitle("MEASURED");
1330             h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1331             h->GetYaxis()->SetTitle("counts");
1332        } break;
1333        default : break;
1334     }
1335 }
1336 //_____________________________________________________________________________
1337 void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t columns)
1338 {
1339    // go through the output file and perform post processing routines
1340    // can either be performed in one go with the unfolding, or at a later stage
1341    if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
1342    TFile readMe(in.Data(), "READ");     // open file read-only
1343    if(readMe.IsZombie()) {
1344        printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1345        return;
1346    }
1347    printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
1348    readMe.ls();
1349    TList* listOfKeys((TList*)readMe.GetListOfKeys());
1350    if(!listOfKeys) {
1351        printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1352        return;
1353    }
1354    // prepare necessary canvasses
1355    TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
1356    TCanvas* canvasOut(new TCanvas("PearsonOut", "PearsonOut"));
1357    TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
1358    TCanvas* canvasRatioMeasuredRefoldedOut(new TCanvas("RefoldedOut", "RefoldedOut"));
1359    TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
1360    TCanvas* canvasSpectraOut(new TCanvas("SpectraOut", "SpectraOut"));
1361    TCanvas* canvasRatio(new TCanvas("Ratio", "Ratio"));
1362    TCanvas* canvasV2(new TCanvas("V2", "V2"));
1363    TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1364    TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
1365    TCanvas* canvasMasterOut(new TCanvas("defaultOut", "defaultOut"));
1366    canvasMISC->Divide(4, 2);
1367    TDirectoryFile* defDir(0x0);
1368
1369    // get an estimate of the number of outputs and find the default set
1370    Int_t cacheMe(0);
1371    for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1372        if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1373            if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1374            cacheMe++;
1375        }
1376    }
1377    Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%4)>0));
1378    canvasIn->Divide(columns, rows);
1379    canvasOut->Divide(columns, rows);
1380    canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1381    canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1382    canvasSpectraIn->Divide(columns, rows);
1383    canvasSpectraOut->Divide(columns, rows);
1384    canvasRatio->Divide(columns, rows);
1385    canvasV2->Divide(columns, rows);
1386
1387    canvasMasterIn->Divide(columns, rows);
1388    canvasMasterOut->Divide(columns, rows);
1389    // extract the default output
1390    TH1D* defUnfoldedIn(0x0);
1391    TH1D* defUnfoldedOut(0x0);
1392    THStack stackIn("StackRatioIn","StackRatioIn");
1393    THStack stackOut("StackRatioOut", "StackRatioOut");
1394    if(defDir) {
1395        TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1396        TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
1397        if(defDirIn) defUnfoldedIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1398        if(defUnfoldedIn) stackIn.Add(defUnfoldedIn);
1399        if(defDirOut) defUnfoldedOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
1400        if(defUnfoldedOut) stackOut.Add(defUnfoldedOut);
1401        printf(" > succesfully extracted default results < \n");
1402    }
1403
1404    // loop through the directories, only plot the graphs if the deconvolution converged
1405    TDirectoryFile* tempDir(0x0);
1406    TDirectoryFile* tempIn(0x0);
1407    TDirectoryFile*  tempOut(0x0);
1408    for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
1409        tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1410        if(!tempDir) continue;
1411        TString dirName(tempDir->GetName());
1412        tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1413        tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1414        j++;
1415        if(tempIn) {
1416            // to see if the unfolding converged try to extract the pearson coefficients
1417            TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1418            if(pIn) {
1419                printf(" - %s in plane converged \n", dirName.Data());
1420                canvasIn->cd(j);
1421                Style(gPad, "PEARSON");
1422                pIn->DrawCopy("colz");
1423                TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1424                if(rIn) {
1425                    Style(rIn);
1426                    printf(" > found RatioRefoldedMeasured < \n");
1427                    canvasRatioMeasuredRefoldedIn->cd(j);
1428                    rIn->Draw("ac");
1429                }
1430                TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1431                TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1432                TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1433                TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1434                if(dvector && avalue && rm && eff) {
1435                    Style(dvector);
1436                    Style(avalue);
1437                    Style(rm);
1438                    Style(eff);
1439                    canvasMISC->cd(1);
1440                    Style(gPad, "SPECTRUM");
1441                    dvector->DrawCopy();
1442                    canvasMISC->cd(2);
1443                    Style(gPad, "SPECTRUM");
1444                    avalue->DrawCopy();
1445                    canvasMISC->cd(3);
1446                    Style(gPad, "PEARSON");
1447                    rm->DrawCopy("colz");
1448                    canvasMISC->cd(4);
1449                    eff->DrawCopy();
1450                } else if(rm && eff) {
1451                    Style(rm);
1452                    Style(eff);
1453                    canvasMISC->cd(3);
1454                    Style(gPad, "PEARSON");
1455                    rm->DrawCopy("colz");
1456                    canvasMISC->cd(4);
1457                    eff->DrawCopy();
1458                }
1459            }
1460            TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1461            TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1462            TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1463            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1464                if(defUnfoldedIn) {
1465                    TH1D* temp((TH1D*)defUnfoldedIn->Clone(Form("defUnfoldedIn_%s", dirName.Data())));
1466                    temp->Divide(unfoldedSpectrum);
1467                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1468                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1469                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1470                    canvasMasterIn->cd(j);
1471                    temp->GetXaxis()->SetRangeUser(0., 2);
1472                    temp->DrawCopy();
1473                }
1474                TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
1475                canvasSpectraIn->cd(j);
1476                Style(gPad);
1477                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1478                unfoldedSpectrum->DrawCopy();
1479                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1480                inputSpectrum->DrawCopy("same");
1481                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1482                refoldedSpectrum->DrawCopy("same");
1483                TLegend* l(AddLegend(gPad));
1484                Style(l);
1485                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1486                    Float_t chi(fitStatus->GetBinContent(1));
1487                    Float_t pen(fitStatus->GetBinContent(2));
1488                    Int_t dof((int)fitStatus->GetBinContent(3));
1489                    Float_t beta(fitStatus->GetBinContent(4));
1490                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1491                } else if (fitStatus) { // only available in SVD fit
1492                    Int_t reg((int)fitStatus->GetBinContent(1));
1493                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1494                }
1495            }
1496        }
1497        if(tempOut) {
1498            TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
1499            if(pOut) {
1500                printf(" - %s out of plane converged \n", dirName.Data());
1501                canvasOut->cd(j);
1502                Style(gPad, "PEARSON");
1503                pOut->DrawCopy("colz");
1504                TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1505                if(rOut) {
1506                    Style(rOut);
1507                    printf(" > found RatioRefoldedMeasured < \n");
1508                    canvasRatioMeasuredRefoldedOut->cd(j);
1509                    rOut->Draw("ac");
1510                }
1511                TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1512                TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1513                TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
1514                TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
1515                if(dvector && avalue && rm && eff) {
1516                    Style(dvector);
1517                    Style(avalue);
1518                    Style(rm);
1519                    Style(eff);
1520                    canvasMISC->cd(5);
1521                    Style(gPad, "SPECTRUM");
1522                    dvector->DrawCopy();
1523                    canvasMISC->cd(6);
1524                    Style(gPad, "SPECTRUM");
1525                    avalue->DrawCopy();
1526                    canvasMISC->cd(7);
1527                    Style(gPad, "PEARSON");
1528                    rm->DrawCopy("colz");
1529                    canvasMISC->cd(8);
1530                    eff->DrawCopy();
1531                } else if(rm && eff) {
1532                    Style(rm);
1533                    Style(eff);
1534                    canvasMISC->cd(3);
1535                    Style(gPad, "PEARSON");
1536                    rm->DrawCopy("colz");
1537                    canvasMISC->cd(4);
1538                    eff->DrawCopy();
1539                }
1540            }
1541            TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
1542            TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
1543            TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
1544            if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1545                if(defUnfoldedOut) {
1546                    TH1D* temp((TH1D*)defUnfoldedOut->Clone(Form("defUnfoldedOut_%s", dirName.Data())));
1547                    temp->Divide(unfoldedSpectrum);
1548                    temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1549                    temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1550                    temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1551                    canvasMasterOut->cd(j);
1552                    temp->GetXaxis()->SetRangeUser(0., 2.);
1553                    temp->DrawCopy();
1554                }
1555                TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
1556                canvasSpectraOut->cd(j);
1557                Style(gPad);
1558                Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1559                unfoldedSpectrum->DrawCopy();
1560                Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1561                inputSpectrum->DrawCopy("same");
1562                Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1563                refoldedSpectrum->DrawCopy("same");
1564                TLegend* l(AddLegend(gPad));
1565                Style(l);
1566                if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1567                    Float_t chi(fitStatus->GetBinContent(1));
1568                    Float_t pen(fitStatus->GetBinContent(2));
1569                    Int_t dof((int)fitStatus->GetBinContent(3));
1570                    Float_t beta(fitStatus->GetBinContent(4));
1571                    l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1572                } else if (fitStatus) { // only available in SVD fit
1573                    Int_t reg((int)fitStatus->GetBinContent(1));
1574                    l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1575                }
1576            }
1577        }
1578        canvasRatio->cd(j);
1579        TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
1580        if(ratioYield) {
1581            Style(ratioYield);
1582            ratioYield->Draw("ac");
1583        }
1584        canvasV2->cd(j);
1585        TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
1586        if(ratioV2) {
1587            Style(ratioV2);
1588            ratioV2->Draw("ac");
1589        }
1590    }
1591    TFile output(out.Data(), "RECREATE");
1592    canvasIn->Write();
1593    canvasOut->Write();
1594    canvasRatioMeasuredRefoldedIn->Write();
1595    canvasRatioMeasuredRefoldedOut->Write();
1596    canvasSpectraIn->Write();
1597    canvasSpectraOut->Write();
1598    canvasRatio->Write();
1599    canvasV2->Write();
1600    canvasMasterIn->Write();
1601    canvasMasterOut->Write();
1602    canvasMISC->Write();
1603    output.Write();
1604    output.Close();
1605 }
1606 //_____________________________________________________________________________
1607 Bool_t AliJetFlowTools::SetRawInput (
1608         TH2D* detectorResponse,  // detector response matrix
1609         TH1D* jetPtIn,           // in plane jet spectrum
1610         TH1D* jetPtOut,          // out of plane jet spectrum
1611         TH1D* dptIn,             // in plane delta pt distribution
1612         TH1D* dptOut,            // out of plane delta pt distribution
1613         Int_t eventCount) {
1614     // set input histograms manually
1615     fDetectorResponse   = detectorResponse;
1616     fSpectrumIn         = jetPtIn;
1617     fSpectrumOut        = jetPtOut;
1618     fDptInDist          = dptIn;
1619     fDptOutDist         = dptOut;
1620     fRawInputProvided   = kTRUE;
1621     // check if all data is provided
1622     if(!fDetectorResponse) {
1623         printf(" fDetectorResponse not found \n ");
1624         return kFALSE;
1625     }
1626     // check if the pt bin for true and rec have been set
1627     if(!fBinsTrue || !fBinsRec) {
1628         printf(" No true or rec bins set, please set binning ! \n");
1629         return kFALSE;
1630     }
1631     if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
1632                           // procedures, these profiles will be nonsensical, user is responsible
1633         fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1634         fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1635         fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1636     }
1637     // normalize spectra to event count if requested
1638     if(fNormalizeSpectra) {
1639         fEventCount = eventCount;
1640         if(fEventCount > 0) {
1641             fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1642             fSpectrumOut->Sumw2();
1643             fSpectrumIn->Scale(1./((double)fEventCount));
1644             fSpectrumOut->Scale(1./((double)fEventCount));
1645         }
1646     }
1647     if(!fNormalizeSpectra && fEventCount > 0) {
1648         fSpectrumIn->Sumw2();       // necessary for correct error propagation of scale
1649         fSpectrumOut->Sumw2();
1650         fSpectrumIn->Scale(1./((double)fEventCount));
1651         fSpectrumOut->Scale(1./((double)fEventCount));
1652     }
1653     fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
1654     fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1655     fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1656     fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1657     fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
1658     fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1659     fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1660     fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1661
1662     return kTRUE;
1663 }
1664 //_____________________________________________________________________________
1665 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
1666 {
1667     // return ratio of h1 / h2
1668     // histograms can have different binning. errors are propagated as uncorrelated
1669     if(!(h1 && h2) ) {
1670         printf(" GetRatio called with NULL argument(s) \n ");
1671         return 0x0;
1672     }
1673     Int_t j(0);
1674     TGraphErrors *gr = new TGraphErrors();
1675     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1676     Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1677     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1678         binCent = h1->GetXaxis()->GetBinCenter(i);
1679         if(xmax > 0. && binCent > xmax) continue;
1680         j = h2->FindBin(binCent);
1681         binWidth = h1->GetXaxis()->GetBinWidth(i);
1682         if(h2->GetBinContent(j) > 0.) {
1683             ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
1684             Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
1685             Double_t B = 0.;
1686             if(h2->GetBinError(j)>0.) {
1687                 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
1688                 error2 = A*A + B*B;
1689             } else error2 = A*A;
1690             if(error2 > 0 ) error2 = TMath::Sqrt(error2);
1691             gr->SetPoint(gr->GetN(),binCent,ratio);
1692             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1693         }
1694     }
1695     if(appendFit) {
1696         TF1* fit(new TF1("lin", "pol0", 10, 100));
1697         gr->Fit(fit);
1698     }
1699     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1700     return gr;
1701 }
1702 //_____________________________________________________________________________
1703 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
1704 {
1705     // get v2 from difference of in plane, out of plane yield
1706     // h1 must hold the in-plane yield, h2 holds the out of plane  yield
1707     // different binning is allowed
1708     // r is the event plane resolution for the chosen centrality
1709     if(!(h1 && h2) ) {
1710         printf(" GetV2 called with NULL argument(s) \n ");
1711         return 0x0;
1712     }
1713     Int_t j(0);
1714     TGraphErrors *gr = new TGraphErrors();
1715     gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1716     Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1717     Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
1718     for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1719         binCent = h1->GetXaxis()->GetBinCenter(i);
1720         j = h2->FindBin(binCent);
1721         binWidth = h1->GetXaxis()->GetBinWidth(i);
1722         if(h2->GetBinContent(j) > 0.) {
1723             in = h1->GetBinContent(i);
1724             ein = h1->GetBinError(i);
1725             out = h2->GetBinContent(j);
1726             eout = h2->GetBinError(j);
1727             ratio = pre*((in-out)/(in+out));
1728             error2 = (r*4.)/(TMath::Pi())*((out*out/(TMath::Power(in+out, 4)))*ein*ein+(in*in/(TMath::Power(in+out, 4)))*eout*eout);
1729             if(error2 > 0) error2 = TMath::Sqrt(error2);
1730             gr->SetPoint(gr->GetN(),binCent,ratio);
1731             gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1732         }
1733     }
1734     if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1735     return gr;
1736 }
1737 //_____________________________________________________________________________
1738 void AliJetFlowTools::WriteObject(TObject* object) {
1739     // write object with unique identifier to active TDirectoryFile
1740     if(!object) {
1741         printf(" > WriteObject:: called with NULL arguments \n ");
1742         return;
1743     } else object->Write();
1744     // FIXME to be implememnted
1745 }
1746 //_____________________________________________________________________________
1747 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
1748     // construt a delta pt response matrix from supplied dpt distribution
1749     // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
1750     // do a weighted rebinning to a (coarser) dpt distribution
1751     // be careful with the binning of the dpt response: it should be equal to that
1752     // of the response matrix, otherwise dpt and response matrices cannot be multiplied
1753     //
1754     // the response matrix will be square and have the same binning
1755     // (min, max and granularity) of the input histogram
1756     Int_t bins(dpt->GetXaxis()->GetNbins());        // number of bins, will also be no of rows, columns
1757     Double_t _bins[bins+1];             // prepare array with bin borders
1758     for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
1759     _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1);    // get upper edge
1760     TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
1761     for(Int_t j(0); j < bins+1; j++) {   // loop on pt true slices j
1762         Bool_t skip = kFALSE;
1763         for(Int_t k(0); k < bins+1; k++) {       // loop on pt gen slices k
1764             (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
1765             if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
1766         }
1767     }
1768     return res;
1769 }
1770 //_____________________________________________________________________________
1771 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1772     if(!binsTrue || !binsRec) {
1773         printf(" > GetUnityResponse:: function called with NULL arguments < \n");
1774         return 0x0;
1775     }
1776     TString name(Form("unityResponse_%s", suffix.Data()));
1777     TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
1778     for(Int_t i(0); i < binsTrue->GetSize(); i++) {
1779         for(Int_t j(0); j < binsRec->GetSize(); j++) {
1780             if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
1781         }
1782     }
1783     return unity;
1784 }
1785 //_____________________________________________________________________________
1786 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) {
1787     // save configuration parameters to histogram
1788     TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 16, -.5, 16.5);
1789     summary->SetBinContent(1, fBetaIn);
1790     summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
1791     summary->SetBinContent(2, fBetaOut);
1792     summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
1793     summary->SetBinContent(3, fCentralityBin);
1794     summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
1795     summary->SetBinContent(4, (int)convergedIn);
1796     summary->GetXaxis()->SetBinLabel(4, "convergedIn");
1797     summary->SetBinContent(5, (int)convergedOut);
1798     summary->GetXaxis()->SetBinLabel(5, "convergedOut");
1799     summary->SetBinContent(6, (int)fAvoidRoundingError);
1800     summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
1801     summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
1802     summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
1803     summary->SetBinContent(8, (int)fPrior);
1804     summary->GetXaxis()->SetBinLabel(8, "fPrior");
1805     summary->SetBinContent(9, fSVDRegIn);
1806     summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
1807     summary->SetBinContent(10, fSVDRegOut);
1808     summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
1809     summary->SetBinContent(11, (int)fSVDToy);
1810     summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
1811     summary->SetBinContent(12, fJetRadius);
1812     summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
1813     summary->SetBinContent(13, (int)fNormalizeSpectra);
1814     summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
1815     summary->SetBinContent(14, (int)fSmoothenSpectrum);
1816     summary->GetXaxis()->SetBinLabel(14, "fSmoothenSpectrum");
1817     summary->SetBinContent(15, (int)fTestMode);
1818     summary->GetXaxis()->SetBinLabel(15, "fTestMode");
1819     summary->SetBinContent(16, (int)fUseDetectorResponse);
1820     summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
1821     summary->Write();
1822 }
1823 //_____________________________________________________________________________
1824 void AliJetFlowTools::ResetAliUnfolding() {
1825      // ugly function: reset all unfolding parameters
1826      TVirtualFitter* fitter(TVirtualFitter::GetFitter());
1827      if(fitter) {
1828          printf(" > Found fitter, will delete it < \n");
1829          delete fitter;
1830      }
1831      if(gMinuit) {
1832          printf(" > Found gMinuit, will re-create it < \n");
1833          delete gMinuit;
1834          gMinuit = new TMinuit;
1835      }
1836      AliUnfolding::fgCorrelationMatrix = 0;
1837      AliUnfolding::fgCorrelationMatrixSquared = 0;
1838      AliUnfolding::fgCorrelationCovarianceMatrix = 0;
1839      AliUnfolding::fgCurrentESDVector = 0;
1840      AliUnfolding::fgEntropyAPriori = 0;
1841      AliUnfolding::fgEfficiency = 0;
1842      AliUnfolding::fgUnfoldedAxis = 0;
1843      AliUnfolding::fgMeasuredAxis = 0;
1844      AliUnfolding::fgFitFunction = 0;
1845      AliUnfolding::fgMaxInput  = -1;
1846      AliUnfolding::fgMaxParams = -1;
1847      AliUnfolding::fgOverflowBinLimit = -1;
1848      AliUnfolding::fgRegularizationWeight = 10000;
1849      AliUnfolding::fgSkipBinsBegin = 0;
1850      AliUnfolding::fgMinuitStepSize = 0.1;
1851      AliUnfolding::fgMinuitPrecision = 1e-6;
1852      AliUnfolding::fgMinuitMaxIterations = 1000000;
1853      AliUnfolding::fgMinuitStrategy = 1.;
1854      AliUnfolding::fgMinimumInitialValue = kFALSE;
1855      AliUnfolding::fgMinimumInitialValueFix = -1;
1856      AliUnfolding::fgNormalizeInput = kFALSE;
1857      AliUnfolding::fgNotFoundEvents = 0;
1858      AliUnfolding::fgSkipBin0InChi2 = kFALSE;
1859      AliUnfolding::fgBayesianSmoothing  = 1;
1860      AliUnfolding::fgBayesianIterations = 10;
1861      AliUnfolding::fgDebug = kFALSE;
1862      AliUnfolding::fgCallCount = 0;
1863      AliUnfolding::fgPowern = 5;
1864      AliUnfolding::fChi2FromFit = 0.;
1865      AliUnfolding::fPenaltyVal  = 0.;
1866      AliUnfolding::fAvgResidual = 0.;
1867      AliUnfolding::fgPrintChi2Details = 0;
1868      AliUnfolding::fgCanvas = 0;
1869      AliUnfolding::fghUnfolded = 0;
1870      AliUnfolding::fghCorrelation = 0;
1871      AliUnfolding::fghEfficiency = 0;
1872      AliUnfolding::fghMeasured = 0;
1873      AliUnfolding::SetMinuitStepSize(1.);
1874      AliUnfolding::SetMinuitPrecision(1e-6);
1875      AliUnfolding::SetMinuitMaxIterations(100000);
1876      AliUnfolding::SetMinuitStrategy(2.);
1877      AliUnfolding::SetDebug(1);
1878 }
1879 //_____________________________________________________________________________
1880 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) {
1881     // protect heap by adding unique qualifier to name
1882     if(!protect) return 0x0;
1883     TH1D* p = (TH1D*)protect->Clone();
1884     TString tempString(fActiveString);
1885     tempString+=suffix;
1886     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1887     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1888     if(kill) delete protect;
1889     return p;
1890 }
1891 //_____________________________________________________________________________
1892 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) {
1893     // protect heap by adding unique qualifier to name
1894     if(!protect) return 0x0;
1895     TH2D* p = (TH2D*)protect->Clone();
1896     TString tempString(fActiveString);
1897     tempString+=suffix;
1898     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1899     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1900     if(kill) delete protect;
1901     return p;
1902 }
1903 //_____________________________________________________________________________
1904 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) {
1905     // protect heap by adding unique qualifier to name
1906     if(!protect) return 0x0;
1907     TGraphErrors* p = (TGraphErrors*)protect->Clone();
1908     TString tempString(fActiveString);
1909     tempString+=suffix;
1910     p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1911     p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
1912     if(kill) delete protect;
1913     return p;
1914 }
1915 //_____________________________________________________________________________
1916