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

#ifndef AliAnalysisTaskJetFlowMC_H
#define AliAnalysisTaskJetFlowMC_H

// root includes
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h" 
#include "TRandom.h"
// aliroot includes
#include "AliAnalysisTaskSE.h"
// forward declarations
class TList;
class TClonesArray;
class TArrayI;

class AliAnalysisTaskJetFlowMC : public AliAnalysisTaskSE 
{
    public:
        // enumerators
        enum detectorType       {kVZEROA, kVZEROC, kVZEROComb};  // detector that was used
        // constructors, destructor
        AliAnalysisTaskJetFlowMC();
        AliAnalysisTaskJetFlowMC(const char *name);
        virtual ~AliAnalysisTaskJetFlowMC();
        // mirror image of PickTrackMaker
        void    UserCreateOutputObjects();
        TH1F*   BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c = -1, Bool_t append = kTRUE);
        TH2F*   BookTH2F(const char* name, const char* x, const char* y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c = -1, Bool_t append = kTRUE);

        void    UserExec(Option_t *option);
        void    SetDebugMode(Bool_t d)                          { fDebug = d;}
        void    SetTracksInName(const char *name)               { fTracksInName     = name; }
        void    SetTracksOutName(const char *name)              { fTracksOutName    = name; }
        // additional setters
        void    SetCentralityClasses(TArrayI* c)                { fCentralityClasses = c; }
        void    SetReferenceDetector(detectorType type)         { fDetectorType = type; }
        void    SetDifferentialV2(TF1* v2, Int_t c = 0)         { fFuncDiffV2[c] = v2; }
        void    SetDifferentialV3(TF1* v3, Int_t c = 0)         { fFuncDiffV3[c] = v3; }
        void    SetDifferentialV2(TH1* v2, Int_t c = 0)         { fHistDiffV2[c] = v2; }
        void    SetDifferentialV3(TH1* v3, Int_t c = 0)         { fHistDiffV3[c] = v3; }
        void    SetIntegratedV2(TH1* v2)                        { fHistIntV2 = v2; }
        void    SetIntegratedV3(TH1* v3)                        { fHistIntV3 = v3; }
        void    SetTrackSpectrum(TF1* ts)                       { fTrackSpectrum = ts; }
        void    SetRandomizeEta(Bool_t b)                       { fRandomizeEta = b; }
        void    SetMultiplicity(Int_t m)                        { fMult = m; }
        void    SetReuseTracks(Bool_t r)                        { fReuseTracks = r; }
        void    SetSingleFragmentationJetSpectrum(TF1* js)      { fJetSpectrumSF = js; }
        void    SetNoOfSFJets(Int_t n)                          { fNoOfSFJets = n; }
        // additional methods
        void    V2AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;
        void    V3AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;
        void    InjectSingleFragmentationJetSpectrum(Int_t nacc);
        void    CalculateEventPlane();
        // inlined helper calculations
        Double_t        GetTrackPt()       const                { return fTrackSpectrum->GetRandom();}
        Double_t        GetTrackEta()       const               { return gRandom->Uniform(-.9, .9); }
        /* inline */    Double_t GetV2(Double_t pt) const { 
            return (fFuncDiffV2[fCenBin]) ? fFuncDiffV2[fCenBin]->Eval(pt) :
            fHistDiffV2[fCenBin]->GetBinContent(fHistDiffV2[fCenBin]->GetXaxis()->FindBin(pt));
        }
        /* inline */    Double_t GetV3(Double_t pt) const { 
            return (fFuncDiffV3[fCenBin]) ? fFuncDiffV3[fCenBin]->Eval(pt) : 
            fHistDiffV3[fCenBin]->GetBinContent(fHistDiffV3[fCenBin]->GetXaxis()->FindBin(pt));
        }
        /* inline */    void GetFlowFluctuation(Double_t& vn) const {
            vn += TMath::Sqrt(2*(vn*.25)*(vn*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1); 
        }
        /* inline */    Double_t PhaseShift(Double_t x) const {  
            while (x>=TMath::TwoPi())x-=TMath::TwoPi();
            while (x<0.)x+=TMath::TwoPi();
        return x; }
        /* inline */    Double_t PhaseShift(Double_t x, Double_t n) const {
            x = PhaseShift(x);
            if(TMath::Nint(n)==2) while (x>TMath::Pi()) x-=TMath::Pi();
            if(TMath::Nint(n)==3) {
                if(x>2.*TMath::TwoPi()/n) x = TMath::TwoPi() - x;
                if(x>TMath::TwoPi()/n) x = TMath::TwoPi()-(x+TMath::TwoPi()/n);
            }
        return x; }
        /* inline */    void SampleVnFromTF1(Double_t &phi) const {
        phi = (fFuncVn) ? fFuncVn->GetRandom(0., TMath::TwoPi()) : 0; }
        /* inline */    void FillHistogramsOriginalData(Double_t& pt, Double_t& eta, Double_t& phi) const {
            fHistOriginalSpectrum[fCenBin]->Fill(pt);    fHistOriginalEtaPhi[fCenBin]->Fill(eta, phi);
            fHistOriginalDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2));
        }
        /* inline */    void FillHistogramsToyData(Double_t& pt, Double_t& eta, Double_t& phi, Double_t& vn) const {
            fHistToySpectrum[fCenBin]->Fill(pt);         fHistToyEtaPhi[fCenBin]->Fill(eta, phi);
            fHistToyDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2));  fHistToyVn[fCenBin]->Fill(pt, vn);
        }
        void            Terminate(Option_t* option);
        void            PrintInfo() const;
    protected:
        TString         fTracksOutName;         // name of output track array
        TString         fTracksInName;          // name of input track array
        TClonesArray   *fTracksIn;              //!track array in
        TClonesArray   *fTracksOut;             //!track array out
        Bool_t          fReuseTracks;           // use original event as template
        Int_t           fMult;                  // multiplicity of new event
        Int_t           fCenBin;                //! centrality bin
        TArrayI*        fCentralityClasses;     // centrality classes (max 10) 
        TF1*            fFuncVn;                //! vn function
        TList*          fOutputList;            // output list
        TF1*            fTrackSpectrum;         // track pt spectrum
        Bool_t          fRandomizeEta;          // randomize eta
        TF1*            fJetSpectrumSF;         // single fragmentation spectrum of jets
        Int_t           fNoOfSFJets;            // number of single fragmentation jets that will be added
        TF1*            fFuncDiffV2[10];        // differential v2 of charged tracks
        TF1*            fFuncDiffV3[10];        // differential v3 of charged tracks
        TH1*            fHistDiffV2[10];        // differential v2 of charged tracks
        TH1*            fHistDiffV3[10];        // differential v3 of charged tracks
        TH1*            fHistIntV2;             // integrated v2 of charged tracks
        TH1*            fHistIntV3;             // integrated v3 of charged tracks
        Float_t         fFlowFluctuations;      // introduce gaussian flow fluctuations of this magnitude
        Int_t           fMaxNumberOfIterations; // max number of iterations for afterburner
        Double_t        fPsi2;                  //! 2nd order event plane orientation
        Double_t        fPsi3;                  //! 3rd order event plane orientation
        Double_t        fPrecisionPhi;          // afterburner precision
        detectorType    fDetectorType;          // type of detector from which the EP is taken
        // output histograms
        TH1F*           fHistOriginalSpectrum[10];      //! original pt spectrum of accepted tracks
        TH2F*           fHistOriginalEtaPhi[10];        //! original eta phi spectrum of accepted tracks
        TH1F*           fHistToySpectrum[10];           //! spectrum of toy (generated) tracks
        TH2F*           fHistToyEtaPhi[10];             //! eta phi spectrum of toy (generated) tracks
        TH1F*           fHistOriginalDeltaPhi[10];      //! original delta phi spectrum
        TH1F*           fHistToyDeltaPhi[10];           //! delta phi spectrum of toy (generated) tracks
        TH2F*           fHistToyVn[10];                 //! generated differential vn values (should equal the differential spectrum)
        TH1F*           fHistSFJetSpectrum;             //! spectrum of generated sf jets
        TH2F*           fHistSFJetEtaPhi;               //! eta phi of generated sf jets
    private:
        AliAnalysisTaskJetFlowMC(const AliAnalysisTaskJetFlowMC&);            // not implemented
        AliAnalysisTaskJetFlowMC &operator=(const AliAnalysisTaskJetFlowMC&); // not implemented

        ClassDef(AliAnalysisTaskJetFlowMC, 2); // Task to generate toy mc PicoTracks based on real events
};
#endif
Commit	Line	Data
24373b38	1	#ifndef AliAnalysisTaskJetFlowMC_H
	2	#define AliAnalysisTaskJetFlowMC_H
	3
	4	// root includes
	5	#include "TF1.h"
	6	#include "TH1F.h"
	7	#include "TH2F.h"
	8	#include "TRandom.h"
	9	// aliroot includes
	10	#include "AliAnalysisTaskSE.h"
	11	// forward declarations
	12	class TList;
	13	class TClonesArray;
	14	class TArrayI;
	15
	16	class AliAnalysisTaskJetFlowMC : public AliAnalysisTaskSE
	17	{
	18	public:
	19	// enumerators
	20	enum detectorType {kVZEROA, kVZEROC, kVZEROComb}; // detector that was used
	21	// constructors, destructor
	22	AliAnalysisTaskJetFlowMC();
	23	AliAnalysisTaskJetFlowMC(const char *name);
	24	virtual ~AliAnalysisTaskJetFlowMC();
	25	// mirror image of PickTrackMaker
	26	void UserCreateOutputObjects();
	27	TH1F* BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c = -1, Bool_t append = kTRUE);
	28	TH2F* BookTH2F(const char* name, const char* x, const char* y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c = -1, Bool_t append = kTRUE);
	29
	30	void UserExec(Option_t *option);
7ffa833f	31	void SetDebugMode(Bool_t d) { fDebug = d;}
24373b38	32	void SetTracksInName(const char *name) { fTracksInName = name; }
	33	void SetTracksOutName(const char *name) { fTracksOutName = name; }
	34	// additional setters
	35	void SetCentralityClasses(TArrayI* c) { fCentralityClasses = c; }
	36	void SetReferenceDetector(detectorType type) { fDetectorType = type; }
	37	void SetDifferentialV2(TF1* v2, Int_t c = 0) { fFuncDiffV2[c] = v2; }
	38	void SetDifferentialV3(TF1* v3, Int_t c = 0) { fFuncDiffV3[c] = v3; }
	39	void SetDifferentialV2(TH1* v2, Int_t c = 0) { fHistDiffV2[c] = v2; }
	40	void SetDifferentialV3(TH1* v3, Int_t c = 0) { fHistDiffV3[c] = v3; }
	41	void SetIntegratedV2(TH1* v2) { fHistIntV2 = v2; }
	42	void SetIntegratedV3(TH1* v3) { fHistIntV3 = v3; }
	43	void SetTrackSpectrum(TF1* ts) { fTrackSpectrum = ts; }
7ffa833f	44	void SetRandomizeEta(Bool_t b) { fRandomizeEta = b; }
db02017e	45	void SetMultiplicity(Int_t m) { fMult = m; }
db02017e	46	void SetReuseTracks(Bool_t r) { fReuseTracks = r; }
24373b38	47	void SetSingleFragmentationJetSpectrum(TF1* js) { fJetSpectrumSF = js; }
	48	void SetNoOfSFJets(Int_t n) { fNoOfSFJets = n; }
	49	// additional methods
	50	void V2AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;
	51	void V3AfterBurner(Double_t& phi, Double_t& eta, Double_t& pt) const;
	52	void InjectSingleFragmentationJetSpectrum(Int_t nacc);
	53	void CalculateEventPlane();
	54	// inlined helper calculations
	55	Double_t GetTrackPt() const { return fTrackSpectrum->GetRandom();}
7ffa833f	56	Double_t GetTrackEta() const { return gRandom->Uniform(-.9, .9); }
24373b38	57	/* inline */ Double_t GetV2(Double_t pt) const {
	58	return (fFuncDiffV2[fCenBin]) ? fFuncDiffV2[fCenBin]->Eval(pt) :
	59	fHistDiffV2[fCenBin]->GetBinContent(fHistDiffV2[fCenBin]->GetXaxis()->FindBin(pt));
	60	}
	61	/* inline */ Double_t GetV3(Double_t pt) const {
	62	return (fFuncDiffV3[fCenBin]) ? fFuncDiffV3[fCenBin]->Eval(pt) :
	63	fHistDiffV3[fCenBin]->GetBinContent(fHistDiffV3[fCenBin]->GetXaxis()->FindBin(pt));
	64	}
	65	/* inline */ void GetFlowFluctuation(Double_t& vn) const {
	66	vn += TMath::Sqrt(2(vn.25)(vn.25))TMath::ErfInverse(2(gRandom->Uniform(0, fFlowFluctuations))-1);
	67	}
	68	/* inline */ Double_t PhaseShift(Double_t x) const {
	69	while (x>=TMath::TwoPi())x-=TMath::TwoPi();
	70	while (x<0.)x+=TMath::TwoPi();
	71	return x; }
	72	/* inline */ Double_t PhaseShift(Double_t x, Double_t n) const {
	73	x = PhaseShift(x);
	74	if(TMath::Nint(n)==2) while (x>TMath::Pi()) x-=TMath::Pi();
	75	if(TMath::Nint(n)==3) {
	76	if(x>2.*TMath::TwoPi()/n) x = TMath::TwoPi() - x;
	77	if(x>TMath::TwoPi()/n) x = TMath::TwoPi()-(x+TMath::TwoPi()/n);
	78	}
	79	return x; }
	80	/* inline */ void SampleVnFromTF1(Double_t &phi) const {
	81	phi = (fFuncVn) ? fFuncVn->GetRandom(0., TMath::TwoPi()) : 0; }
	82	/* inline */ void FillHistogramsOriginalData(Double_t& pt, Double_t& eta, Double_t& phi) const {
	83	fHistOriginalSpectrum[fCenBin]->Fill(pt); fHistOriginalEtaPhi[fCenBin]->Fill(eta, phi);
	84	fHistOriginalDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2));
	85	}
	86	/* inline */ void FillHistogramsToyData(Double_t& pt, Double_t& eta, Double_t& phi, Double_t& vn) const {
	87	fHistToySpectrum[fCenBin]->Fill(pt); fHistToyEtaPhi[fCenBin]->Fill(eta, phi);
	88	fHistToyDeltaPhi[fCenBin]->Fill(PhaseShift(phi-fPsi2, 2)); fHistToyVn[fCenBin]->Fill(pt, vn);
	89	}
	90	void Terminate(Option_t* option);
	91	void PrintInfo() const;
	92	protected:
	93	TString fTracksOutName; // name of output track array
	94	TString fTracksInName; // name of input track array
	95	TClonesArray *fTracksIn; //!track array in
	96	TClonesArray *fTracksOut; //!track array out
db02017e	97	Bool_t fReuseTracks; // use original event as template
db02017e	98	Int_t fMult; // multiplicity of new event
7ffa833f	99	Int_t fCenBin; //! centrality bin
24373b38	100	TArrayI* fCentralityClasses; // centrality classes (max 10)
	101	TF1* fFuncVn; //! vn function
	102	TList* fOutputList; // output list
	103	TF1* fTrackSpectrum; // track pt spectrum
7ffa833f	104	Bool_t fRandomizeEta; // randomize eta
24373b38	105	TF1* fJetSpectrumSF; // single fragmentation spectrum of jets
	106	Int_t fNoOfSFJets; // number of single fragmentation jets that will be added
	107	TF1* fFuncDiffV2[10]; // differential v2 of charged tracks
	108	TF1* fFuncDiffV3[10]; // differential v3 of charged tracks
	109	TH1* fHistDiffV2[10]; // differential v2 of charged tracks
	110	TH1* fHistDiffV3[10]; // differential v3 of charged tracks
	111	TH1* fHistIntV2; // integrated v2 of charged tracks
	112	TH1* fHistIntV3; // integrated v3 of charged tracks
	113	Float_t fFlowFluctuations; // introduce gaussian flow fluctuations of this magnitude
	114	Int_t fMaxNumberOfIterations; // max number of iterations for afterburner
	115	Double_t fPsi2; //! 2nd order event plane orientation
	116	Double_t fPsi3; //! 3rd order event plane orientation
	117	Double_t fPrecisionPhi; // afterburner precision
	118	detectorType fDetectorType; // type of detector from which the EP is taken
	119	// output histograms
	120	TH1F* fHistOriginalSpectrum[10]; //! original pt spectrum of accepted tracks
	121	TH2F* fHistOriginalEtaPhi[10]; //! original eta phi spectrum of accepted tracks
	122	TH1F* fHistToySpectrum[10]; //! spectrum of toy (generated) tracks
	123	TH2F* fHistToyEtaPhi[10]; //! eta phi spectrum of toy (generated) tracks
	124	TH1F* fHistOriginalDeltaPhi[10]; //! original delta phi spectrum
	125	TH1F* fHistToyDeltaPhi[10]; //! delta phi spectrum of toy (generated) tracks
	126	TH2F* fHistToyVn[10]; //! generated differential vn values (should equal the differential spectrum)
	127	TH1F* fHistSFJetSpectrum; //! spectrum of generated sf jets
	128	TH2F* fHistSFJetEtaPhi; //! eta phi of generated sf jets
	129	private:
	130	AliAnalysisTaskJetFlowMC(const AliAnalysisTaskJetFlowMC&); // not implemented
	131	AliAnalysisTaskJetFlowMC &operator=(const AliAnalysisTaskJetFlowMC&); // not implemented
	132
db02017e	133	ClassDef(AliAnalysisTaskJetFlowMC, 2); // Task to generate toy mc PicoTracks based on real events
24373b38	134	};
24373b38	135	#endif