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