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

// Simple class to generate toy events which can be fed into the jet finder
//
// Author: Redmer Alexander Bertens, Utrecht University, 2013
// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl

// root includes
#include "TChain.h"
#include "TList.h"
#include "TClonesArray.h"
#include "TArrayI.h"
#include "TRandom3.h"
// aliroot includes
#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskJetFlowMC.h"
#include "AliLog.h"
#include "AliPicoTrack.h"

class AliAnalysisTaskJetFlowMC;
using namespace std;

ClassImp(AliAnalysisTaskJetFlowMC)

//_____________________________________________________________________________
AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fTracksOutName("JetFlowMC"),  fTracksInName("PicoTrack"), fTracksIn(0),  fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
    // default constructor for root IO
    for(Int_t i(0); i < 10; i++) {
        fFuncDiffV2[i]                  = 0x0;
        fFuncDiffV3[i]                  = 0x0;
        fHistDiffV2[i]                  = 0x0;
        fHistDiffV3[i]                  = 0x0;
        fHistOriginalSpectrum[i]        = 0x0;
        fHistOriginalEtaPhi[i]          = 0x0;
        fHistToySpectrum[i]             = 0x0;
        fHistToyEtaPhi[i]               = 0x0;
        fHistOriginalDeltaPhi[i]        = 0x0;
        fHistToyDeltaPhi[i]             = 0x0;
        fHistToyVn[i]                   = 0x0;
    }
}
//_____________________________________________________________________________
AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name, Int_t seed) : AliAnalysisTaskSE(name), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
    // constructor
    DefineInput(0, TChain::Class());
    DefineOutput(1, TList::Class());
    for(Int_t i(0); i < 10; i++) {
        fFuncDiffV2[i]                  = 0x0;
        fFuncDiffV3[i]                  = 0x0;
        fHistDiffV2[i]                  = 0x0;
        fHistDiffV3[i]                  = 0x0;
        fHistOriginalSpectrum[i]        = 0x0;
        fHistOriginalEtaPhi[i]          = 0x0;
        fHistToySpectrum[i]             = 0x0;
        fHistToyEtaPhi[i]               = 0x0;
        fHistOriginalDeltaPhi[i]        = 0x0;
        fHistToyDeltaPhi[i]             = 0x0;
        fHistToyVn[i]                   = 0x0;
    }
    // by default construction replace gRandom
    if(seed > -1) {
        if(gRandom) delete gRandom;
        gRandom = new TRandom3(seed);
    }
}
//_____________________________________________________________________________
AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
{
    // desctructor, claim ownership of deleted objects by setting member pointers to NULL
    if(fOutputList) {
      delete fOutputList;
      fOutputList= NULL;
    }
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
{
    // Create my user objects.
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
    fTracksOut = new TClonesArray("AliPicoTrack");
    fTracksOut->SetName(fTracksOutName);
    fOutputList = new TList();
    fOutputList->SetOwner(kTRUE);
    if(!fCentralityClasses) { // classes must be defined at this point
        Int_t c[] = {0, 90};
        fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
    }
    if(fHistIntV2 && !fHistIntV3) {      // define function
        fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
        fFuncVn->SetParameter(0, 1.);        // normalization
        fFuncVn->SetParameter(3, 0.2);       // v2
        fFuncVn->FixParameter(1, 1.);        // constant
        fFuncVn->FixParameter(2, 2.);        // constant
    } else if (!fHistIntV2 && fHistIntV3) {
        fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
        fFuncVn->SetParameter(0, 1.);        // normalization
        fFuncVn->SetParameter(3, 0.2);       // v3
        fFuncVn->FixParameter(1, 1.);        // constant
        fFuncVn->FixParameter(2, 3.);        // constant
    } else if (fHistIntV2 && fHistIntV3) {
        fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
        fFuncVn->SetParameter(0, 1.);       // normalization
        fFuncVn->SetParameter(3, 0.2);      // v2
        fFuncVn->FixParameter(1, 1.);       // constant
        fFuncVn->FixParameter(2, 2.);       // constant
        fFuncVn->FixParameter(5, 3.);       // constant
        fFuncVn->SetParameter(7, 0.2);      // v3
    }
    // add the generator objects that have been added to the task
    if(fTrackSpectrum)  fOutputList->Add(fTrackSpectrum);
    if(fJetSpectrumSF)  fOutputList->Add(fJetSpectrumSF);
    if(fHistIntV2)      fOutputList->Add(fHistIntV2);
    if(fHistIntV3)      fOutputList->Add(fHistIntV3);
    // create the QA histos
    for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
        fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
        fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
        fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
        fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
        fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
        fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
        fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 1000, 0, 200, 200, 0, .2, i);
        // add to outputlist
        if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);
        if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);
        if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);
        if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);
    }
    if(fJetSpectrumSF) {
        fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);
        fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());
    }
    // reset the random seed
    if(gRandom) {
        delete gRandom;
        gRandom = new TRandom3(0);
    }

    fOutputList->Sort();
    PostData(1, fOutputList);
}
//_____________________________________________________________________________
TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
{
    // book a TH1F and connect it to the output container
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    if(!fOutputList) return 0x0;
    TString title(name);
    if(c!=-1) { // format centrality dependent histograms accordingly
        name = Form("%s_%i", name, c);
        title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
    }
    title += Form(";%s;[counts]", x);
    TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
    histogram->Sumw2();
    if(append) fOutputList->Add(histogram);
    return histogram;   
}
//_____________________________________________________________________________
TH2F* AliAnalysisTaskJetFlowMC::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, Bool_t append)
{
    // book a TH2F and connect it to the output container
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    if(!fOutputList) return 0x0;
    TString title(name);
    if(c!=-1) { // format centrality dependent histograms accordingly
        name = Form("%s_%i", name, c);
        title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
    }
    title += Form(";%s;%s", x, y);
    TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
    histogram->Sumw2();
    if(append) fOutputList->Add(histogram);
    return histogram;   
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::UserExec(Option_t *) 
{
    // user exec, called for each event.
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    if(!AliAnalysisManager::GetAnalysisManager()) return;
    // retrieve tracks from input.
    if (!fTracksIn) { 
        fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
        if (!fTracksIn) {
          AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data())); 
          return;
        }
    }
    // get the centrality bin for qa plots
    Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
    fCenBin = -1;
    for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
        if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
            fCenBin = i;
            break; }
    }
    if(fCenBin < 0) return;
    // add tracks to event if not yet there
    fTracksOut->Delete();
    if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);
    fTracksOut->Clear();
    // get the event plane
    CalculateEventPlane();
    Int_t nacc(0);
    if(fReuseTracks) {
        const Int_t Ntracks = fTracksIn->GetEntriesFast();
        for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
            AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));
            if(!track) continue;
            Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(fRandomizeEta ? GetTrackEta() : track->Eta());
            // fill qa histo's before applying any (possible) afterburner
            FillHistogramsOriginalData(pt, eta, phi);
            if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin])        V2AfterBurner(phi, eta, pt);
            else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin])   V3AfterBurner(phi, eta, pt);
            else if(fHistIntV2 || fHistIntV3)                       SampleVnFromTF1(phi);        
            /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1); 
            nacc++;
        }
    } else {
        Double_t pt(0), eta(0), phi(0);
        for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
            pt = GetTrackPt();
            eta = gRandom->Uniform(-.9, .9);
            phi = gRandom->Uniform(0., TMath::TwoPi());
            // fill qa histo's before applying any (possible) afterburner
            FillHistogramsOriginalData(pt, eta, phi);
            if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin])        V2AfterBurner(phi, eta, pt);
            else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin])   V3AfterBurner(phi, eta, pt);
            else if(fHistIntV2 || fHistIntV3)                       SampleVnFromTF1(phi);        
            /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1); 
            nacc++;
        }
    }
    if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
    PostData(1, fOutputList);
    if(fDebug > 0) PrintInfo();
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
{
    // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
    if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    phi = gRandom->Uniform(0, TMath::TwoPi());
    Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);
    if(TMath::AreEqualAbs(v2, 0, 1e-5)) { 
        FillHistogramsToyData(pt, eta, phi, v2);
        return;
    }
    // introduce flow fluctuations (gaussian)
    if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
    for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
        phiprev=phi; //store last value for comparison
        f =  phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
        fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
        phi -= f/fp;
        if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break; 
    }
    FillHistogramsToyData(pt, eta, phi, v2);
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
{
    // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations
    if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    phi = gRandom->Uniform(0, TMath::TwoPi());
    Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);
    if(TMath::AreEqualAbs(v3, 0, 1e-5)) {
        FillHistogramsToyData(pt, eta, phi, v3);
        return;
    }
    // introduce flow fluctuations (gaussian)
    if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
    for (Int_t i=0; i<fMaxNumberOfIterations; i++)  {
        phiprev=phi; //store last value for comparison
        f =  phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));
        fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative
        phi -= f/fp;
        if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
    }
    FillHistogramsToyData(pt, eta, phi, v3);
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
{
    // inject single fragmentation jet spectrum to the tclones array, note that emcal params 
    // equal the barrel kinematics to pass the track and jet cuts later on
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {
        Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());
        /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);
        fHistSFJetSpectrum->Fill(pt);
        fHistSFJetEtaPhi->Fill(eta, phi);
        ++i;
    }
    nacc = 0;
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {
    // grab the event plane orientation from the AliVEvent header
    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);    
    Double_t a(0), b(0), e(0), f(0);
    switch (fDetectorType) {
        case kVZEROA : {
            fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
            fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
            break;
        }
        case kVZEROC : {
            fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
            fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
            break;
                       }
        case kVZEROComb : {
            fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
            fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
            break;
        }
        default : break;
    }
    // if requested, pass the event plane to the phi distribution
    if(fHistIntV2 && fHistIntV3) {
        fFuncVn->SetParameter(4, fPsi2);        fFuncVn->SetParameter(6, fPsi3);
        Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
        Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
        if(fFlowFluctuations > -10) {
            GetFlowFluctuation(v2);             GetFlowFluctuation(v3);
        }
        fFuncVn->SetParameter(3, v2);           fFuncVn->SetParameter(7, v3);
    } else if (fHistIntV2) {
        fFuncVn->SetParameter(4, fPsi2);
        Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
        if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
        fFuncVn->SetParameter(3, v2);
    } else if (fHistIntV3) {
        fFuncVn->SetParameter(4, fPsi3);
        Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
        if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
        fFuncVn->SetParameter(3, v3);
    }
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
{
    // terminate
}
//_____________________________________________________________________________
void AliAnalysisTaskJetFlowMC::PrintInfo() const
{
    // print info
    printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());
    printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());

}
//_____________________________________________________________________________
Commit	Line	Data
5a6bcfc6	1	// Simple class to generate toy events which can be fed into the jet finder
	2	//
	3	// Author: Redmer Alexander Bertens, Utrecht University, 2013
	4	// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
	5
	6	// root includes
	7	#include "TChain.h"
	8	#include "TList.h"
	9	#include "TClonesArray.h"
	10	#include "TArrayI.h"
20bd6574	11	#include "TRandom3.h"
5a6bcfc6	12	// aliroot includes
	13	#include "AliAODEvent.h"
	14	#include "AliAnalysisManager.h"
	15	#include "AliAnalysisTaskJetFlowMC.h"
	16	#include "AliLog.h"
	17	#include "AliPicoTrack.h"
	18
	19	class AliAnalysisTaskJetFlowMC;
	20	using namespace std;
	21
	22	ClassImp(AliAnalysisTaskJetFlowMC)
	23
	24	//_____________________________________________________________________________
db02017e	25	AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
5a6bcfc6	26	// default constructor for root IO
	27	for(Int_t i(0); i < 10; i++) {
	28	fFuncDiffV2[i] = 0x0;
	29	fFuncDiffV3[i] = 0x0;
	30	fHistDiffV2[i] = 0x0;
	31	fHistDiffV3[i] = 0x0;
	32	fHistOriginalSpectrum[i] = 0x0;
	33	fHistOriginalEtaPhi[i] = 0x0;
	34	fHistToySpectrum[i] = 0x0;
	35	fHistToyEtaPhi[i] = 0x0;
	36	fHistOriginalDeltaPhi[i] = 0x0;
	37	fHistToyDeltaPhi[i] = 0x0;
	38	fHistToyVn[i] = 0x0;
	39	}
	40	}
	41	//_____________________________________________________________________________
5a35c43f	42	AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name, Int_t seed) : AliAnalysisTaskSE(name), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
5a6bcfc6	43	// constructor
	44	DefineInput(0, TChain::Class());
	45	DefineOutput(1, TList::Class());
	46	for(Int_t i(0); i < 10; i++) {
	47	fFuncDiffV2[i] = 0x0;
	48	fFuncDiffV3[i] = 0x0;
	49	fHistDiffV2[i] = 0x0;
	50	fHistDiffV3[i] = 0x0;
	51	fHistOriginalSpectrum[i] = 0x0;
	52	fHistOriginalEtaPhi[i] = 0x0;
	53	fHistToySpectrum[i] = 0x0;
	54	fHistToyEtaPhi[i] = 0x0;
	55	fHistOriginalDeltaPhi[i] = 0x0;
	56	fHistToyDeltaPhi[i] = 0x0;
	57	fHistToyVn[i] = 0x0;
	58	}
5a35c43f	59	// by default construction replace gRandom
	60	if(seed > -1) {
	61	if(gRandom) delete gRandom;
	62	gRandom = new TRandom3(seed);
	63	}
5a6bcfc6	64	}
	65	//_____________________________________________________________________________
	66	AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
	67	{
	68	// desctructor, claim ownership of deleted objects by setting member pointers to NULL
	69	if(fOutputList) {
	70	delete fOutputList;
	71	fOutputList= NULL;
	72	}
	73	}
	74	//_____________________________________________________________________________
	75	void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
	76	{
	77	// Create my user objects.
	78	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	79	fTracksOut = new TClonesArray("AliPicoTrack");
	80	fTracksOut->SetName(fTracksOutName);
	81	fOutputList = new TList();
	82	fOutputList->SetOwner(kTRUE);
	83	if(!fCentralityClasses) { // classes must be defined at this point
	84	Int_t c[] = {0, 90};
	85	fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
	86	}
	87	if(fHistIntV2 && !fHistIntV3) { // define function
	88	fFuncVn = new TF1("kV2", "[0]([1]+[2][3]TMath::Cos([2](x-[4])))", 0, TMath::TwoPi());
	89	fFuncVn->SetParameter(0, 1.); // normalization
	90	fFuncVn->SetParameter(3, 0.2); // v2
	91	fFuncVn->FixParameter(1, 1.); // constant
	92	fFuncVn->FixParameter(2, 2.); // constant
	93	} else if (!fHistIntV2 && fHistIntV3) {
	94	fFuncVn = new TF1("kV3", "[0]([1]+[2][3]TMath::Cos([2](x-[4])))", 0, TMath::TwoPi());
	95	fFuncVn->SetParameter(0, 1.); // normalization
	96	fFuncVn->SetParameter(3, 0.2); // v3
	97	fFuncVn->FixParameter(1, 1.); // constant
	98	fFuncVn->FixParameter(2, 3.); // constant
	99	} else if (fHistIntV2 && fHistIntV3) {
	100	fFuncVn = new TF1("kCombined", "[0]([1]+[2]([3]TMath::Cos([2](x-[4]))+[7]TMath::Cos([5](x-[6]))))", 0, TMath::TwoPi());
	101	fFuncVn->SetParameter(0, 1.); // normalization
	102	fFuncVn->SetParameter(3, 0.2); // v2
	103	fFuncVn->FixParameter(1, 1.); // constant
	104	fFuncVn->FixParameter(2, 2.); // constant
	105	fFuncVn->FixParameter(5, 3.); // constant
	106	fFuncVn->SetParameter(7, 0.2); // v3
	107	}
	108	// add the generator objects that have been added to the task
	109	if(fTrackSpectrum) fOutputList->Add(fTrackSpectrum);
	110	if(fJetSpectrumSF) fOutputList->Add(fJetSpectrumSF);
	111	if(fHistIntV2) fOutputList->Add(fHistIntV2);
	112	if(fHistIntV3) fOutputList->Add(fHistIntV3);
	113	// create the QA histos
	114	for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
d0ef8ed2	115	fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
5a6bcfc6	116	fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
d0ef8ed2	117	fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
5a6bcfc6	118	fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
	119	fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
	120	fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
3880af5d	121	fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 1000, 0, 200, 200, 0, .2, i);
5a6bcfc6	122	// add to outputlist
	123	if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);
	124	if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);
	125	if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);
	126	if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);
	127	}
	128	if(fJetSpectrumSF) {
	129	fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);
	130	fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());
	131	}
20bd6574	132	// reset the random seed
	133	if(gRandom) {
	134	delete gRandom;
	135	gRandom = new TRandom3(0);
	136	}
	137
5a6bcfc6	138	fOutputList->Sort();
	139	PostData(1, fOutputList);
	140	}
	141	//_____________________________________________________________________________
	142	TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
	143	{
	144	// book a TH1F and connect it to the output container
	145	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	146	if(!fOutputList) return 0x0;
	147	TString title(name);
	148	if(c!=-1) { // format centrality dependent histograms accordingly
	149	name = Form("%s_%i", name, c);
	150	title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
	151	}
	152	title += Form(";%s;[counts]", x);
	153	TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
	154	histogram->Sumw2();
	155	if(append) fOutputList->Add(histogram);
	156	return histogram;
	157	}
	158	//_____________________________________________________________________________
	159	TH2F* AliAnalysisTaskJetFlowMC::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, Bool_t append)
	160	{
	161	// book a TH2F and connect it to the output container
	162	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	163	if(!fOutputList) return 0x0;
	164	TString title(name);
	165	if(c!=-1) { // format centrality dependent histograms accordingly
	166	name = Form("%s_%i", name, c);
	167	title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
	168	}
	169	title += Form(";%s;%s", x, y);
	170	TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
	171	histogram->Sumw2();
	172	if(append) fOutputList->Add(histogram);
	173	return histogram;
	174	}
	175	//_____________________________________________________________________________
	176	void AliAnalysisTaskJetFlowMC::UserExec(Option_t *)
	177	{
	178	// user exec, called for each event.
	179	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	180	if(!AliAnalysisManager::GetAnalysisManager()) return;
	181	// retrieve tracks from input.
	182	if (!fTracksIn) {
	183	fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
	184	if (!fTracksIn) {
	185	AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data()));
	186	return;
	187	}
	188	}
	189	// get the centrality bin for qa plots
	190	Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
	191	fCenBin = -1;
	192	for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
	193	if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
	194	fCenBin = i;
	195	break; }
	196	}
	197	if(fCenBin < 0) return;
	198	// add tracks to event if not yet there
	199	fTracksOut->Delete();
	200	if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);
	201	fTracksOut->Clear();
202	// get the event plane
203	CalculateEventPlane();
5a6bcfc6	204	Int_t nacc(0);
db02017e	205	if(fReuseTracks) {
	206	const Int_t Ntracks = fTracksIn->GetEntriesFast();
	207	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
	208	AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));
	209	if(!track) continue;
	210	Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(fRandomizeEta ? GetTrackEta() : track->Eta());
	211	// fill qa histo's before applying any (possible) afterburner
	212	FillHistogramsOriginalData(pt, eta, phi);
	213	if(fHistDiffV2[fCenBin] \|\| fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
	214	else if(fHistDiffV3[fCenBin] \|\| fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
	215	else if(fHistIntV2 \|\| fHistIntV3) SampleVnFromTF1(phi);
	216	/AliPicoTrack picotrack =/ new ((fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1);
	217	nacc++;
	218	}
	219	} else {
	220	Double_t pt(0), eta(0), phi(0);
	221	for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
	222	pt = GetTrackPt();
	223	eta = gRandom->Uniform(-.9, .9);
	224	phi = gRandom->Uniform(0., TMath::TwoPi());
	225	// fill qa histo's before applying any (possible) afterburner
	226	FillHistogramsOriginalData(pt, eta, phi);
	227	if(fHistDiffV2[fCenBin] \|\| fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
	228	else if(fHistDiffV3[fCenBin] \|\| fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
	229	else if(fHistIntV2 \|\| fHistIntV3) SampleVnFromTF1(phi);
	230	/AliPicoTrack picotrack =/ new ((fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1);
	231	nacc++;
	232	}
5a6bcfc6	233	}
	234	if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
	235	PostData(1, fOutputList);
	236	if(fDebug > 0) PrintInfo();
	237	}
	238	//_____________________________________________________________________________
	239	void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
	240	{
	241	// similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
	242	if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	243	phi = gRandom->Uniform(0, TMath::TwoPi());
	244	Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);
	245	if(TMath::AreEqualAbs(v2, 0, 1e-5)) {
	246	FillHistogramsToyData(pt, eta, phi, v2);
	247	return;
	248	}
	249	// introduce flow fluctuations (gaussian)
	250	if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
	251	for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
	252	phiprev=phi; //store last value for comparison
	253	f = phi-phi0+v2TMath::Sin(2.(phi-fPsi2));
	254	fp = 1.0+2.0v2TMath::Cos(2.*(phi-fPsi2)); //first derivative
	255	phi -= f/fp;
	256	if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
	257	}
	258	FillHistogramsToyData(pt, eta, phi, v2);
	259	}
	260	//_____________________________________________________________________________
	261	void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
	262	{
	263	// similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations
	264	if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	265	phi = gRandom->Uniform(0, TMath::TwoPi());
	266	Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);
	267	if(TMath::AreEqualAbs(v3, 0, 1e-5)) {
	268	FillHistogramsToyData(pt, eta, phi, v3);
	269	return;
	270	}
	271	// introduce flow fluctuations (gaussian)
	272	if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
	273	for (Int_t i=0; i<fMaxNumberOfIterations; i++) {
	274	phiprev=phi; //store last value for comparison
	275	f = phi-phi0+2./3.v3TMath::Sin(3.*(phi-fPsi3));
	276	fp = 1.0+2.0v3TMath::Cos(3.*(phi-fPsi3)); //first derivative
	277	phi -= f/fp;
	278	if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
	279	}
	280	FillHistogramsToyData(pt, eta, phi, v3);
	281	}
	282	//_____________________________________________________________________________
	283	void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
	284	{
	285	// inject single fragmentation jet spectrum to the tclones array, note that emcal params
	286	// equal the barrel kinematics to pass the track and jet cuts later on
	287	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
	288	for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {
	289	Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());
	290	/AliPicoTrack picotrack =/ new ((fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);
	291	fHistSFJetSpectrum->Fill(pt);
	292	fHistSFJetEtaPhi->Fill(eta, phi);
	293	++i;
	294	}
	295	nacc = 0;
	296	}
297	//_____________________________________________________________________________
298	void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {
299	// grab the event plane orientation from the AliVEvent header
300	if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
301	Double_t a(0), b(0), e(0), f(0);
302	switch (fDetectorType) {
303	case kVZEROA : {
304	fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
305	fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
306	break;
307	}
308	case kVZEROC : {
309	fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
310	fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
311	break;
312	}
313	case kVZEROComb : {
314	fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
315	fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
316	break;
317	}
318	default : break;
319	}
320	// if requested, pass the event plane to the phi distribution
321	if(fHistIntV2 && fHistIntV3) {
322	fFuncVn->SetParameter(4, fPsi2); fFuncVn->SetParameter(6, fPsi3);
323	Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
324	Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
325	if(fFlowFluctuations > -10) {
326	GetFlowFluctuation(v2); GetFlowFluctuation(v3);
327	}
328	fFuncVn->SetParameter(3, v2); fFuncVn->SetParameter(7, v3);
329	} else if (fHistIntV2) {
330	fFuncVn->SetParameter(4, fPsi2);
331	Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
332	if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
333	fFuncVn->SetParameter(3, v2);
334	} else if (fHistIntV3) {
335	fFuncVn->SetParameter(4, fPsi3);
336	Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
337	if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
338	fFuncVn->SetParameter(3, v3);
339	}
340	}
341	//_____________________________________________________________________________
342	void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
343	{
344	// terminate
345	}
346	//_____________________________________________________________________________
347	void AliAnalysisTaskJetFlowMC::PrintInfo() const
348	{
349	// print info
350	printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());
351	printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());
352
353	}
354	//_____________________________________________________________________________