[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.h

/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * 
 * See cxx source for full Copyright notice */ 
 /* $Id$ */

/* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
 * see implementation for additional information */

#ifndef AliJetFlowTools_H
#define AliJetFlowTools_H

// root forward declarations
class TF1;
class TH1D;
class TH2D;
class TCanvas;
class TLegend;
class TString;
class TArrayD;
class TGraph;
class TGraphErrors;
class TObjArray;
// aliroot forward declarations
class AliAnaChargedJetResponseMaker;
class AliUnfolding;
// root includes
#include "TMatrixD.h"
#include "TList.h"
#include "TDirectoryFile.h"
#include "TFile.h"
#include "TProfile.h"
#include "TVirtualPad.h"
//_____________________________________________________________________________
class AliJetFlowTools {
    public: 
        AliJetFlowTools();
    protected:
        ~AliJetFlowTools();     // not implemented (deliberately). object ownership is a bit messy in this class
                                // since most (or all) of the objects are owned by the input and output files
    public:
        // enumerators
        enum unfoldingAlgorithm {       // type of unfolding alrogithm
            kChi2,                      // chi^2 unfolding, implemented in AliUnfolding
            kBayesian,                  // Bayesian unfolding, implemented in RooUnfold
            kBayesianAli,               // Bayesian unfolding, implemented in AliUnfolding
            kSVD,                       // SVD unfolding, implemented in RooUnfold
            kNone };                    // no unfolding
        enum prior {                    // prior that is used for unfolding
            kPriorChi2,                 // prior from chi^2 method
            kPriorMeasured,             // use measured spectrum as prior
            kPriorPythia };             // use pythia spectrum as prior
        enum histoType {                // histogram identifier, only used internally
            kInPlaneSpectrum,           // default style for spectrum
            kOutPlaneSpectrum,
            kUnfoldedSpectrum,
            kFoldedSpectrum,
            kMeasuredSpectrum,
            kBar,                       // default style for bar histogram
            kEmpty };                   // default style
        // setters, interface to the class
        void            SetSaveFull(Bool_t b)           {fSaveFull              = b;}
        void            SetInputList(TList* list)       {
            fInputList          = list;
            fRefreshInput       = kTRUE;
        }
        void            SetOutputFileName(TString name) {fOutputFileName        = name;}
        void            CreateOutputList(TString name) {
            // create a new output list and add it to the full output
            if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
            fOutputFile->cd();  // avoid nested dirs
            if(name.EqualTo(fActiveString)) {
                printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
                name+="_2";
            }
            fActiveString = name;
            fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
            fActiveDir->cd();
        }
        void            SetCentralityBin(Int_t bin)             {fCentralityBin         = bin;}
        void            SetDetectorResponse(TH2D* dr)           {fDetectorResponse      = dr;}
        void            SetJetFindingEfficiency(TH1D* e)        {fJetFindingEff         = e;}
        void            SetBinsTrue(TArrayD* bins)              {fBinsTrue              = bins;}
        void            SetBinsRec(TArrayD* bins)               {fBinsRec               = bins;}
        void            SetBinsTruePrior(TArrayD* bins)         {fBinsTruePrior         = bins;}
        void            SetBinsRecPrior(TArrayD* bins)          {fBinsRecPrior          = bins;}
        void            SetSVDReg(Int_t r)                      {fSVDRegIn              = r; fSVDRegOut         = r;}
        void            SetSVDReg(Int_t in, Int_t out)          {fSVDRegIn              = in; fSVDRegOut        = out;}
        void            SetSVDToy(Bool_t b, Float_t r)          {fSVDToy                = b; fJetRadius         = r;}
        void            SetBeta(Double_t b)                     {fBetaIn                = b; fBetaOut           = b;}
        void            SetBeta(Double_t i, Double_t o)         {fBetaIn                = i; fBetaOut           = o;}
        void            SetBayesianIter(Int_t i)                {fBayesianIterIn        = i; fBayesianIterOut    = i;}
        void            SetBayesianIter(Int_t i, Int_t o)       {fBayesianIterIn        = i; fBayesianIterOut    = o;}
        void            SetBayesianSmooth(Float_t s)            {fBayesianSmoothIn      = s; fBayesianSmoothOut = s;}
        void            SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn      = i; fBayesianSmoothOut = o;}
        void            SetAvoidRoundingError(Bool_t r)         {fAvoidRoundingError    = r;}
        void            SetUnfoldingAlgorithm(unfoldingAlgorithm ua)    {fUnfoldingAlgorithm                    = ua;}
        void            SetPrior(prior p)                       {fPrior                 = p;}
        void            SetPrior(prior p, TH1D* spectrum)       {fPrior                 = p; fPriorUser         = spectrum;}
        void            SetNormalizeSpectra(Bool_t b)           {fNormalizeSpectra      = b;}
        void            SetNormalizeSpectra(Int_t e)            { // use to normalize to this no of events
            fEventCount         = e;
            fNormalizeSpectra   = kFALSE;
        }
        void            SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
            fSmoothenPrior      = b;
            fFitMin             = min;
            fFitMax             = max;
            fFitStart           = start;
            fSmoothenCounts     = counts;
        }
        void            SetTestMode(Bool_t t)                   {fTestMode              = t;}
        void            SetEventPlaneResolution(Double_t r)     {fEventPlaneRes         = r;}
        void            SetUseDetectorResponse(Bool_t r)        {fUseDetectorResponse   = r;}
        void            SetUseDptResponse(Bool_t r)             {fUseDptResponse        = r;}
        void            SetTrainPowerFit(Bool_t t)              {fTrainPower            = t;}
        void            SetDphiUnfolding(Bool_t i)              {fDphiUnfolding         = i;}
        void            SetDphiDptUnfolding(Bool_t i)           {fDphiDptUnfolding      = i;}
        void            SetExLJDpt(Bool_t i)                    {fExLJDpt               = i;}
        void            SetWeightFunction(TF1* w)               {fResponseMaker->SetRMMergeWeightFunction(w);}
        void            Make();
        void            MakeAU();       // test function, use with caution (09012014)
        void            Finish() {
            fOutputFile->cd();
            if(fRMSSpectrumIn)  fRMSSpectrumIn->Write();
            if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
            if(fRMSRatio)       fRMSRatio->Write();
            fOutputFile->Close();}
        void            PostProcess(
                TString def,
                Int_t columns = 4,
                Float_t rangeLow = 20,
                Float_t rangeUp = 80,
                TString in = "UnfoldedSpectra.root", 
                TString out = "ProcessedSpectra.root") const;
        void            GetCorrelatedUncertainty(
                TArrayI* variationsIn,
                TArrayI* variationsOut,
                TString type = "",
                Int_t columns = 4,
                Float_t rangeLow = 20,
                Float_t rangeUp = 80,
                TString in = "UnfoldedSpectra.root", 
                TString out = "CorrelatedUncertainty.root") const;
        void            GetShapeUncertainty(
                TArrayI* regularizationIn,
                TArrayI* regularizationOut,
                TArrayI* trueBinIn = 0x0,
                TArrayI* trueBinOut = 0x0,
                TArrayI* recBinIn = 0x0,
                TArrayI* recBinOut = 0x0,
                Int_t columns = 4,
                Float_t rangeLow = 20,
                Float_t rangeUp = 80,
                TString in = "UnfoldedSpectra.root", 
                TString out = "ShapeUncertainty.root") const;
        Bool_t          SetRawInput (
                TH2D* detectorResponse, // detector response matrix
                TH1D* jetPtIn,          // in plane jet spectrum
                TH1D* jetPtOut,         // out of plane jet spectrum
                TH1D* dptIn,            // in plane delta pt distribution
                TH1D* dptOut,           // out of plane delta pt distribution
                Int_t eventCount = 0);  // event count (optional)
        // static const helper functions, mainly histogram manipulation
        static TH1D*    ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
        static TH2D*    ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
        static TH2D*    NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
        static TH1D*    RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
        TH2D*           RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
        static TH2D*    MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
        static TH1D*    NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
        static TGraphErrors*    GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
        static TGraphErrors*    GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = .7, TString name = "");
        static void     WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
        static TH2D*    ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
        static TH2D*    GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
        void            SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
        static TMatrixD*        CalculatePearsonCoefficients(TMatrixD* covmat);
        static TH1D*    SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
        // set style
        void            SetTitleFontSize(Double_t s)    {fTitleFontSize = s;}
        static void     Style();
        static void     Style(TCanvas* c, TString style = "PEARSON");
        static void     Style(TVirtualPad* c, TString style = "SPECTRUM");
        static void     Style(TLegend* l);
        static void     Style(TH1* h, EColor col = kBlue, histoType = kEmpty);
        static void     Style(TGraph* h, EColor col = kBlue, histoType = kEmpty);
        static TLegend* AddLegend(TVirtualPad* p)       {return p->BuildLegend(.565, .663, .882, .883);}
        static void     SavePadToPDF(TVirtualPad* pad)  {pad->SaveAs(Form("%s.pdf", pad->GetName()));}
        // interface to AliUnfolding, not necessary but nice to have all parameters in one place
        static void     SetMinuitStepSize(Float_t s)    {AliUnfolding::SetMinuitStepSize(s);}
        static void     SetMinuitPrecision(Float_t s)   {AliUnfolding::SetMinuitPrecision(s);}
        static void     SetMinuitPrecision(Int_t i)     {AliUnfolding::SetMinuitMaxIterations(i);}
        static void     SetMinuitStrategy(Double_t s)   {AliUnfolding::SetMinuitStrategy(s);}
        static void     SetDebug(Int_t d)               {AliUnfolding::SetDebug(d);}
    private:
        Bool_t          PrepareForUnfolding(); 
        Bool_t          PrepareForUnfolding(Int_t low, Int_t up);
        TH1D*           GetPrior(                       const TH1D* measuredJetSpectrum,
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency);
        TH1D*           UnfoldWrapper(                  const TH1D* measuredJetSpectrum, 
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency = 0x0);
        TH1D*           UnfoldSpectrumChi2(             const TH1D* measuredJetSpectrum, 
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency = 0x0);
        TH1D*           UnfoldSpectrumSVD(              const TH1D* measuredJetSpectrum, 
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency = 0x0);
        TH1D*           UnfoldSpectrumBayesianAli(      const TH1D* measuredJetSpectrum, 
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency = 0x0);
        TH1D*           UnfoldSpectrumBayesian(         const TH1D* measuredJetSpectrum, 
                                                        const TH2D* resizedResponse,
                                                        const TH1D* kinematicEfficiency,
                                                        const TH1D* measuredJetSpectrumTrueBins,
                                                        const TString suffix,
                                                        const TH1D* jetFindingEfficiency = 0x0);
        void            DoIntermediateSystematics(
                TArrayI* variationsIn,
                TArrayI* variationsOut,
                TH1D*& relativeErrorInUp,
                TH1D*& relativeErrorInLow,
                TH1D*& relativeErrorOutUp,
                TH1D*& relativeErrorOutLow,
                TH1D*& relativeSystematicIn,
                TH1D*& relativeSystematicOut,
                Int_t columns,
                Float_t rangeLow,
                Float_t rangeUp,
                TFile* readMe, 
                TString source = "") const;
        static void     ResetAliUnfolding();
        // give object a unique name via the 'protect heap' functions. 
        // may seem redundant, but some internal functions of root (e.g.
        // ProjectionY()) check for existing objects by name and re-use them
        TH1D*           ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
        TH2D*           ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
        TGraphErrors*   ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
        // members, accessible via setters
        AliAnaChargedJetResponseMaker*  fResponseMaker; // utility object
        TF1*                    fPower;                 // smoothening fit
        Bool_t                  fSaveFull;              // save all generated histograms to file
        TString                 fActiveString;          // identifier of active output
        TDirectoryFile*         fActiveDir;             // active directory
        TList*                  fInputList;             // input list
        Bool_t                  fRefreshInput;          // re-read the input (called automatically if input list changes)
        TString                 fOutputFileName;        // output file name
        TFile*                  fOutputFile;            // output file
        Int_t                   fCentralityBin;         // centrality bin
        TH2D*                   fDetectorResponse;      // detector response
        TH1D*                   fJetFindingEff;         // jet finding efficiency
        Double_t                fBetaIn;                // regularization strength, in plane unfolding
        Double_t                fBetaOut;               // regularization strength, out of plane unfoldign
        Int_t                   fBayesianIterIn;        // bayesian regularization parameter, in plane unfolding
        Int_t                   fBayesianIterOut;       // bayesian regularization parameter, out plane unfolding
        Float_t                 fBayesianSmoothIn;      // bayesian smoothening parameter (AliUnfolding)
        Float_t                 fBayesianSmoothOut;     // bayesian smoothening parameter (AliUnfolding)
        Bool_t                  fAvoidRoundingError;    // set dpt to zero for small values far from the diagonal
        unfoldingAlgorithm      fUnfoldingAlgorithm;    // algorithm used for unfolding
        prior                   fPrior;                 // prior for unfolding
        TH1D*                   fPriorUser;             // user supplied prior (e.g. pythia spectrum)
        TArrayD*                fBinsTrue;              // pt true bins
        TArrayD*                fBinsRec;               // pt rec bins
        TArrayD*                fBinsTruePrior;         // holds true bins for the chi2 prior for SVD. setting this is optional
        TArrayD*                fBinsRecPrior;          // holds rec bins for the chi2 prior for SVD. setting this is optional
        Int_t                   fSVDRegIn;              // svd regularization (a good starting point is half of the number of bins)
        Int_t                   fSVDRegOut;             // svd regularization out of plane
        Bool_t                  fSVDToy;                // use toy to estimate coveriance matrix for SVD method
        Float_t                 fJetRadius;             // jet radius (for SVD toy)
        Int_t                   fEventCount;            // number of events
        Bool_t                  fNormalizeSpectra;      // normalize spectra to event count
        Bool_t                  fSmoothenPrior;         // smoothen the tail of the measured spectrum using a powerlaw fit
        Float_t                 fFitMin;                // lower bound of smoothening fit
        Float_t                 fFitMax;                // upper bound of smoothening fit
        Float_t                 fFitStart;              // from this value, use smoothening
        Bool_t                  fSmoothenCounts;        // fill smoothened spectrum with counts
        Bool_t                  fTestMode;              // unfold with unity response for testing
        Bool_t                  fRawInputProvided;      // input histograms provided, not read from file
        Double_t                fEventPlaneRes;         // event plane resolution for current centrality
        Bool_t                  fUseDetectorResponse;   // add detector response to unfolding
        Bool_t                  fUseDptResponse;        // add dpt response to unfolding
        Bool_t                  fTrainPower;            // don't clear the params of fPower for call to Make
                                                        // might give more stable results, but possibly introduces
                                                        // a bias / dependency on previous iterations
        Bool_t                  fDphiUnfolding;         // do the unfolding in in and out of plane orientation
        Bool_t                  fDphiDptUnfolding;      // do the unfolding in dphi and dpt bins (to fit v2)
        Bool_t                  fExLJDpt;               // exclude randon cones with leading jet
        Double_t                fTitleFontSize;         // title font size
        // members, set internally
        TProfile*               fRMSSpectrumIn;         // rms of in plane spectra of converged unfoldings
        TProfile*               fRMSSpectrumOut;        // rms of out of plane spectra of converged unfoldings
        TProfile*               fRMSRatio;              // rms of ratio of converged unfolded spectra
        TProfile*               fRMSV2;                 // rms of v2 of converged unfolded spectra
        TH2D*                   fDeltaPtDeltaPhi;       // delta pt delta phi distribution
        TH2D*                   fJetPtDeltaPhi;         // jet pt delta phi distribution
        TH1D*                   fSpectrumIn;            // in plane jet pt spectrum
        TH1D*                   fSpectrumOut;           // out of plane jet pt spectrum
        TH1D*                   fDptInDist;             // in plane dpt distribution
        TH1D*                   fDptOutDist;            // out of plane dpt distribution
        TH2D*                   fDptIn;                 // in plane dpt matrix
        TH2D*                   fDptOut;                // out plane dpt matrix
        TH2D*                   fFullResponseIn;        // full response matrix, in plane
        TH2D*                   fFullResponseOut;       // full response matrix, out of plane
        // copy and assignment 
        AliJetFlowTools(const AliJetFlowTools&);             // not implemented
        AliJetFlowTools& operator=(const AliJetFlowTools&);  // not implemented
};
#endif
//_____________________________________________________________________________
Commit	Line	Data
dc1455ee	1	/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	2	* See cxx source for full Copyright notice */
	3	/* $Id$ */
	4
51e6bc5a	5	/* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
	6	* see implementation for additional information */
	7
dc1455ee	8	#ifndef AliJetFlowTools_H
	9	#define AliJetFlowTools_H
	10
549b5f40	11	// root forward declarations
4292ca60	12	class TF1;
dc1455ee	13	class TH1D;
dc1455ee	14	class TH2D;
d7ec324f	15	class TCanvas;
d7ec324f	16	class TLegend;
dc1455ee	17	class TString;
dc1455ee	18	class TArrayD;
53547ff2	19	class TGraph;
dc1455ee	20	class TGraphErrors;
dc1455ee	21	class TObjArray;
4292ca60	22	// aliroot forward declarations
	23	class AliAnaChargedJetResponseMaker;
	24	class AliUnfolding;
	25	// root includes
dc1455ee	26	#include "TMatrixD.h"
dc1455ee	27	#include "TList.h"
ad04a83c	28	#include "TDirectoryFile.h"
dc1455ee	29	#include "TFile.h"
4292ca60	30	#include "TProfile.h"
d7ec324f	31	#include "TVirtualPad.h"
dc1455ee	32	//_____________________________________________________________________________
	33	class AliJetFlowTools {
	34	public:
	35	AliJetFlowTools();
	36	protected:
53547ff2 RAB	37	~AliJetFlowTools(); // not implemented (deliberately). object ownership is a bit messy in this class
53547ff2 RAB	38	// since most (or all) of the objects are owned by the input and output files
dc1455ee	39	public:
51e6bc5a	40	// enumerators
549b5f40 RAB	41	enum unfoldingAlgorithm { // type of unfolding alrogithm
	42	kChi2, // chi^2 unfolding, implemented in AliUnfolding
	43	kBayesian, // Bayesian unfolding, implemented in RooUnfold
	44	kBayesianAli, // Bayesian unfolding, implemented in AliUnfolding
	45	kSVD, // SVD unfolding, implemented in RooUnfold
	46	kNone }; // no unfolding
	47	enum prior { // prior that is used for unfolding
	48	kPriorChi2, // prior from chi^2 method
18698978	49	kPriorMeasured, // use measured spectrum as prior
18698978	50	kPriorPythia }; // use pythia spectrum as prior
549b5f40	51	enum histoType { // histogram identifier, only used internally
18698978	52	kInPlaneSpectrum, // default style for spectrum
53547ff2 RAB	53	kOutPlaneSpectrum,
	54	kUnfoldedSpectrum,
	55	kFoldedSpectrum,
	56	kMeasuredSpectrum,
18698978	57	kBar, // default style for bar histogram
18698978	58	kEmpty }; // default style
ef12d5a5	59	// setters, interface to the class
4292ca60	60	void SetSaveFull(Bool_t b) {fSaveFull = b;}
	61	void SetInputList(TList* list) {
	62	fInputList = list;
	63	fRefreshInput = kTRUE;
	64	}
dc1455ee	65	void SetOutputFileName(TString name) {fOutputFileName = name;}
	66	void CreateOutputList(TString name) {
	67	// create a new output list and add it to the full output
ad04a83c	68	if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
ad04a83c	69	fOutputFile->cd(); // avoid nested dirs
ef12d5a5	70	if(name.EqualTo(fActiveString)) {
	71	printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
	72	name+="_2";
	73	}
dc1455ee	74	fActiveString = name;
ad04a83c	75	fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
ad04a83c	76	fActiveDir->cd();
dc1455ee	77	}
549b5f40 RAB	78	void SetCentralityBin(Int_t bin) {fCentralityBin = bin;}
	79	void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;}
	80	void SetJetFindingEfficiency(TH1D* e) {fJetFindingEff = e;}
	81	void SetBinsTrue(TArrayD* bins) {fBinsTrue = bins;}
	82	void SetBinsRec(TArrayD* bins) {fBinsRec = bins;}
	83	void SetBinsTruePrior(TArrayD* bins) {fBinsTruePrior = bins;}
	84	void SetBinsRecPrior(TArrayD* bins) {fBinsRecPrior = bins;}
	85	void SetSVDReg(Int_t r) {fSVDRegIn = r; fSVDRegOut = r;}
	86	void SetSVDReg(Int_t in, Int_t out) {fSVDRegIn = in; fSVDRegOut = out;}
	87	void SetSVDToy(Bool_t b, Float_t r) {fSVDToy = b; fJetRadius = r;}
	88	void SetBeta(Double_t b) {fBetaIn = b; fBetaOut = b;}
	89	void SetBeta(Double_t i, Double_t o) {fBetaIn = i; fBetaOut = o;}
	90	void SetBayesianIter(Int_t i) {fBayesianIterIn = i; fBayesianIterOut = i;}
	91	void SetBayesianIter(Int_t i, Int_t o) {fBayesianIterIn = i; fBayesianIterOut = o;}
	92	void SetBayesianSmooth(Float_t s) {fBayesianSmoothIn = s; fBayesianSmoothOut = s;}
	93	void SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn = i; fBayesianSmoothOut = o;}
	94	void SetAvoidRoundingError(Bool_t r) {fAvoidRoundingError = r;}
	95	void SetUnfoldingAlgorithm(unfoldingAlgorithm ua) {fUnfoldingAlgorithm = ua;}
	96	void SetPrior(prior p) {fPrior = p;}
18698978	97	void SetPrior(prior p, TH1D* spectrum) {fPrior = p; fPriorUser = spectrum;}
549b5f40 RAB	98	void SetNormalizeSpectra(Bool_t b) {fNormalizeSpectra = b;}
	99	void SetNormalizeSpectra(Int_t e) { // use to normalize to this no of events
	100	fEventCount = e;
	101	fNormalizeSpectra = kFALSE;
4292ca60	102	}
549b5f40 RAB	103	void SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
549b5f40 RAB	104	fSmoothenPrior = b;
4292ca60	105	fFitMin = min;
	106	fFitMax = max;
	107	fFitStart = start;
549b5f40	108	fSmoothenCounts = counts;
4292ca60	109	}
549b5f40	110	void SetTestMode(Bool_t t) {fTestMode = t;}
ef12d5a5	111	void SetEventPlaneResolution(Double_t r) {fEventPlaneRes = r;}
ef12d5a5	112	void SetUseDetectorResponse(Bool_t r) {fUseDetectorResponse = r;}
549b5f40 RAB	113	void SetUseDptResponse(Bool_t r) {fUseDptResponse = r;}
549b5f40 RAB	114	void SetTrainPowerFit(Bool_t t) {fTrainPower = t;}
486fb24e	115	void SetDphiUnfolding(Bool_t i) {fDphiUnfolding = i;}
	116	void SetDphiDptUnfolding(Bool_t i) {fDphiDptUnfolding = i;}
	117	void SetExLJDpt(Bool_t i) {fExLJDpt = i;}
	118	void SetWeightFunction(TF1* w) {fResponseMaker->SetRMMergeWeightFunction(w);}
dc1455ee	119	void Make();
486fb24e	120	void MakeAU(); // test function, use with caution (09012014)
4292ca60	121	void Finish() {
4292ca60	122	fOutputFile->cd();
5e11c41c	123	if(fRMSSpectrumIn) fRMSSpectrumIn->Write();
	124	if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
	125	if(fRMSRatio) fRMSRatio->Write();
4292ca60	126	fOutputFile->Close();}
53547ff2 RAB	127	void PostProcess(
53547ff2 RAB	128	TString def,
486fb24e	129	Int_t columns = 4,
87233f72	130	Float_t rangeLow = 20,
87233f72	131	Float_t rangeUp = 80,
53547ff2	132	TString in = "UnfoldedSpectra.root",
486fb24e	133	TString out = "ProcessedSpectra.root") const;
18698978	134	void GetCorrelatedUncertainty(
f3ba6c8e	135	TArrayI* variationsIn,
f3ba6c8e	136	TArrayI* variationsOut,
18698978	137	TString type = "",
f3ba6c8e	138	Int_t columns = 4,
	139	Float_t rangeLow = 20,
	140	Float_t rangeUp = 80,
	141	TString in = "UnfoldedSpectra.root",
18698978	142	TString out = "CorrelatedUncertainty.root") const;
	143	void GetShapeUncertainty(
	144	TArrayI* regularizationIn,
	145	TArrayI* regularizationOut,
	146	TArrayI* trueBinIn = 0x0,
	147	TArrayI* trueBinOut = 0x0,
	148	TArrayI* recBinIn = 0x0,
	149	TArrayI* recBinOut = 0x0,
	150	Int_t columns = 4,
	151	Float_t rangeLow = 20,
	152	Float_t rangeUp = 80,
	153	TString in = "UnfoldedSpectra.root",
	154	TString out = "ShapeUncertainty.root") const;
4292ca60	155	Bool_t SetRawInput (
	156	TH2D* detectorResponse, // detector response matrix
	157	TH1D* jetPtIn, // in plane jet spectrum
	158	TH1D* jetPtOut, // out of plane jet spectrum
	159	TH1D* dptIn, // in plane delta pt distribution
	160	TH1D* dptOut, // out of plane delta pt distribution
	161	Int_t eventCount = 0); // event count (optional)
dc1455ee	162	// static const helper functions, mainly histogram manipulation
	163	static TH1D* ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
	164	static TH2D* ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
512ced40	165	static TH2D* NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
53547ff2	166	static TH1D* RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
4292ca60	167	TH2D* RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
d7ec324f	168	static TH2D* MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
dc1455ee	169	static TH1D* NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
4292ca60	170	static TGraphErrors* GetRatio(TH1 h1 = 0x0, TH1 h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
4e4f12b6	171	static TGraphErrors* GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = .7, TString name = "");
549b5f40	172	static void WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
4292ca60	173	static TH2D* ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
ef12d5a5	174	static TH2D* GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
549b5f40	175	void SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
d7ec324f	176	static TMatrixD* CalculatePearsonCoefficients(TMatrixD* covmat);
549b5f40	177	static TH1D* SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
18698978	178	// set style
	179	void SetTitleFontSize(Double_t s) {fTitleFontSize = s;}
	180	static void Style();
d7ec324f	181	static void Style(TCanvas* c, TString style = "PEARSON");
d7ec324f	182	static void Style(TVirtualPad* c, TString style = "SPECTRUM");
53547ff2 RAB	183	static void Style(TLegend* l);
	184	static void Style(TH1* h, EColor col = kBlue, histoType = kEmpty);
	185	static void Style(TGraph* h, EColor col = kBlue, histoType = kEmpty);
18698978	186	static TLegend* AddLegend(TVirtualPad* p) {return p->BuildLegend(.565, .663, .882, .883);}
512ced40	187	static void SavePadToPDF(TVirtualPad* pad) {pad->SaveAs(Form("%s.pdf", pad->GetName()));}
4292ca60	188	// interface to AliUnfolding, not necessary but nice to have all parameters in one place
	189	static void SetMinuitStepSize(Float_t s) {AliUnfolding::SetMinuitStepSize(s);}
	190	static void SetMinuitPrecision(Float_t s) {AliUnfolding::SetMinuitPrecision(s);}
	191	static void SetMinuitPrecision(Int_t i) {AliUnfolding::SetMinuitMaxIterations(i);}
	192	static void SetMinuitStrategy(Double_t s) {AliUnfolding::SetMinuitStrategy(s);}
	193	static void SetDebug(Int_t d) {AliUnfolding::SetDebug(d);}
dc1455ee	194	private:
dc1455ee	195	Bool_t PrepareForUnfolding();
486fb24e	196	Bool_t PrepareForUnfolding(Int_t low, Int_t up);
549b5f40 RAB	197	TH1D* GetPrior( const TH1D* measuredJetSpectrum,
	198	const TH2D* resizedResponse,
	199	const TH1D* kinematicEfficiency,
	200	const TH1D* measuredJetSpectrumTrueBins,
	201	const TString suffix,
	202	const TH1D* jetFindingEfficiency);
486fb24e	203	TH1D* UnfoldWrapper( const TH1D* measuredJetSpectrum,
	204	const TH2D* resizedResponse,
	205	const TH1D* kinematicEfficiency,
	206	const TH1D* measuredJetSpectrumTrueBins,
	207	const TString suffix,
	208	const TH1D* jetFindingEfficiency = 0x0);
549b5f40 RAB	209	TH1D* UnfoldSpectrumChi2( const TH1D* measuredJetSpectrum,
	210	const TH2D* resizedResponse,
	211	const TH1D* kinematicEfficiency,
	212	const TH1D* measuredJetSpectrumTrueBins,
	213	const TString suffix,
	214	const TH1D* jetFindingEfficiency = 0x0);
	215	TH1D* UnfoldSpectrumSVD( const TH1D* measuredJetSpectrum,
	216	const TH2D* resizedResponse,
	217	const TH1D* kinematicEfficiency,
	218	const TH1D* measuredJetSpectrumTrueBins,
	219	const TString suffix,
	220	const TH1D* jetFindingEfficiency = 0x0);
	221	TH1D* UnfoldSpectrumBayesianAli( const TH1D* measuredJetSpectrum,
	222	const TH2D* resizedResponse,
	223	const TH1D* kinematicEfficiency,
	224	const TH1D* measuredJetSpectrumTrueBins,
	225	const TString suffix,
	226	const TH1D* jetFindingEfficiency = 0x0);
	227	TH1D* UnfoldSpectrumBayesian( const TH1D* measuredJetSpectrum,
	228	const TH2D* resizedResponse,
	229	const TH1D* kinematicEfficiency,
	230	const TH1D* measuredJetSpectrumTrueBins,
	231	const TString suffix,
	232	const TH1D* jetFindingEfficiency = 0x0);
18698978	233	void DoIntermediateSystematics(
	234	TArrayI* variationsIn,
	235	TArrayI* variationsOut,
	236	TH1D*& relativeErrorInUp,
	237	TH1D*& relativeErrorInLow,
	238	TH1D*& relativeErrorOutUp,
	239	TH1D*& relativeErrorOutLow,
	240	TH1D*& relativeSystematicIn,
	241	TH1D*& relativeSystematicOut,
	242	Int_t columns,
	243	Float_t rangeLow,
	244	Float_t rangeUp,
	245	TFile* readMe,
	246	TString source = "") const;
4292ca60	247	static void ResetAliUnfolding();
ef12d5a5	248	// give object a unique name via the 'protect heap' functions.
	249	// may seem redundant, but some internal functions of root (e.g.
	250	// ProjectionY()) check for existing objects by name and re-use them
f3ba6c8e	251	TH1D* ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
	252	TH2D* ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "") const;
	253	TGraphErrors* ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "") const;
dc1455ee	254	// members, accessible via setters
4292ca60	255	AliAnaChargedJetResponseMaker* fResponseMaker; // utility object
	256	TF1* fPower; // smoothening fit
	257	Bool_t fSaveFull; // save all generated histograms to file
51e6bc5a	258	TString fActiveString; // identifier of active output
	259	TDirectoryFile* fActiveDir; // active directory
	260	TList* fInputList; // input list
4292ca60	261	Bool_t fRefreshInput; // re-read the input (called automatically if input list changes)
51e6bc5a	262	TString fOutputFileName; // output file name
	263	TFile* fOutputFile; // output file
	264	Int_t fCentralityBin; // centrality bin
51e6bc5a	265	TH2D* fDetectorResponse; // detector response
53547ff2	266	TH1D* fJetFindingEff; // jet finding efficiency
4292ca60	267	Double_t fBetaIn; // regularization strength, in plane unfolding
4292ca60	268	Double_t fBetaOut; // regularization strength, out of plane unfoldign
549b5f40 RAB	269	Int_t fBayesianIterIn; // bayesian regularization parameter, in plane unfolding
	270	Int_t fBayesianIterOut; // bayesian regularization parameter, out plane unfolding
	271	Float_t fBayesianSmoothIn; // bayesian smoothening parameter (AliUnfolding)
	272	Float_t fBayesianSmoothOut; // bayesian smoothening parameter (AliUnfolding)
51e6bc5a	273	Bool_t fAvoidRoundingError; // set dpt to zero for small values far from the diagonal
	274	unfoldingAlgorithm fUnfoldingAlgorithm; // algorithm used for unfolding
	275	prior fPrior; // prior for unfolding
18698978	276	TH1D* fPriorUser; // user supplied prior (e.g. pythia spectrum)
51e6bc5a	277	TArrayD* fBinsTrue; // pt true bins
51e6bc5a	278	TArrayD* fBinsRec; // pt rec bins
ef12d5a5	279	TArrayD* fBinsTruePrior; // holds true bins for the chi2 prior for SVD. setting this is optional
ef12d5a5	280	TArrayD* fBinsRecPrior; // holds rec bins for the chi2 prior for SVD. setting this is optional
4292ca60	281	Int_t fSVDRegIn; // svd regularization (a good starting point is half of the number of bins)
4292ca60	282	Int_t fSVDRegOut; // svd regularization out of plane
51e6bc5a	283	Bool_t fSVDToy; // use toy to estimate coveriance matrix for SVD method
51e6bc5a	284	Float_t fJetRadius; // jet radius (for SVD toy)
20abfcc4	285	Int_t fEventCount; // number of events
20abfcc4	286	Bool_t fNormalizeSpectra; // normalize spectra to event count
549b5f40	287	Bool_t fSmoothenPrior; // smoothen the tail of the measured spectrum using a powerlaw fit
4292ca60	288	Float_t fFitMin; // lower bound of smoothening fit
	289	Float_t fFitMax; // upper bound of smoothening fit
	290	Float_t fFitStart; // from this value, use smoothening
549b5f40	291	Bool_t fSmoothenCounts; // fill smoothened spectrum with counts
ef12d5a5	292	Bool_t fTestMode; // unfold with unity response for testing
4292ca60	293	Bool_t fRawInputProvided; // input histograms provided, not read from file
ef12d5a5	294	Double_t fEventPlaneRes; // event plane resolution for current centrality
ef12d5a5	295	Bool_t fUseDetectorResponse; // add detector response to unfolding
549b5f40	296	Bool_t fUseDptResponse; // add dpt response to unfolding
ef12d5a5	297	Bool_t fTrainPower; // don't clear the params of fPower for call to Make
486fb24e	298	// might give more stable results, but possibly introduces
ef12d5a5	299	// a bias / dependency on previous iterations
486fb24e	300	Bool_t fDphiUnfolding; // do the unfolding in in and out of plane orientation
	301	Bool_t fDphiDptUnfolding; // do the unfolding in dphi and dpt bins (to fit v2)
	302	Bool_t fExLJDpt; // exclude randon cones with leading jet
18698978	303	Double_t fTitleFontSize; // title font size
dc1455ee	304	// members, set internally
4292ca60	305	TProfile* fRMSSpectrumIn; // rms of in plane spectra of converged unfoldings
	306	TProfile* fRMSSpectrumOut; // rms of out of plane spectra of converged unfoldings
	307	TProfile* fRMSRatio; // rms of ratio of converged unfolded spectra
ef12d5a5	308	TProfile* fRMSV2; // rms of v2 of converged unfolded spectra
4292ca60	309	TH2D* fDeltaPtDeltaPhi; // delta pt delta phi distribution
	310	TH2D* fJetPtDeltaPhi; // jet pt delta phi distribution
	311	TH1D* fSpectrumIn; // in plane jet pt spectrum
	312	TH1D* fSpectrumOut; // out of plane jet pt spectrum
	313	TH1D* fDptInDist; // in plane dpt distribution
	314	TH1D* fDptOutDist; // out of plane dpt distribution
	315	TH2D* fDptIn; // in plane dpt matrix
	316	TH2D* fDptOut; // out plane dpt matrix
	317	TH2D* fFullResponseIn; // full response matrix, in plane
	318	TH2D* fFullResponseOut; // full response matrix, out of plane
dc1455ee	319	// copy and assignment
	320	AliJetFlowTools(const AliJetFlowTools&); // not implemented
	321	AliJetFlowTools& operator=(const AliJetFlowTools&); // not implemented
	322	};
	323	#endif
	324	//_____________________________________________________________________________