+// Simple class to generate toy events which can be fed into the jet finder\r
+//\r
+// Author: Redmer Alexander Bertens, Utrecht University, 2013\r
+// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl\r
+\r
+// root includes\r
+#include "TChain.h"\r
+#include "TList.h"\r
+#include "TClonesArray.h"\r
+#include "TArrayI.h"\r
+// aliroot includes\r
+#include "AliAODEvent.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAnalysisTaskJetFlowMC.h"\r
+#include "AliLog.h"\r
+#include "AliPicoTrack.h"\r
+\r
+class AliAnalysisTaskJetFlowMC;\r
+using namespace std;\r
+\r
+ClassImp(AliAnalysisTaskJetFlowMC)\r
+\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), 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) {\r
+ // default constructor for root IO\r
+ for(Int_t i(0); i < 10; i++) {\r
+ fFuncDiffV2[i] = 0x0;\r
+ fFuncDiffV3[i] = 0x0;\r
+ fHistDiffV2[i] = 0x0;\r
+ fHistDiffV3[i] = 0x0;\r
+ fHistOriginalSpectrum[i] = 0x0;\r
+ fHistOriginalEtaPhi[i] = 0x0;\r
+ fHistToySpectrum[i] = 0x0;\r
+ fHistToyEtaPhi[i] = 0x0;\r
+ fHistOriginalDeltaPhi[i] = 0x0;\r
+ fHistToyDeltaPhi[i] = 0x0;\r
+ fHistToyVn[i] = 0x0;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name) : AliAnalysisTaskSE(name), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), 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) {\r
+ // constructor\r
+ DefineInput(0, TChain::Class());\r
+ DefineOutput(1, TList::Class());\r
+ for(Int_t i(0); i < 10; i++) {\r
+ fFuncDiffV2[i] = 0x0;\r
+ fFuncDiffV3[i] = 0x0;\r
+ fHistDiffV2[i] = 0x0;\r
+ fHistDiffV3[i] = 0x0;\r
+ fHistOriginalSpectrum[i] = 0x0;\r
+ fHistOriginalEtaPhi[i] = 0x0;\r
+ fHistToySpectrum[i] = 0x0;\r
+ fHistToyEtaPhi[i] = 0x0;\r
+ fHistOriginalDeltaPhi[i] = 0x0;\r
+ fHistToyDeltaPhi[i] = 0x0;\r
+ fHistToyVn[i] = 0x0;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()\r
+{\r
+ // desctructor, claim ownership of deleted objects by setting member pointers to NULL\r
+ if(fOutputList) {\r
+ delete fOutputList;\r
+ fOutputList= NULL;\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()\r
+{\r
+ // Create my user objects.\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);\r
+ fTracksOut = new TClonesArray("AliPicoTrack");\r
+ fTracksOut->SetName(fTracksOutName);\r
+ fOutputList = new TList();\r
+ fOutputList->SetOwner(kTRUE);\r
+ if(!fCentralityClasses) { // classes must be defined at this point\r
+ Int_t c[] = {0, 90};\r
+ fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);\r
+ }\r
+ if(fHistIntV2 && !fHistIntV3) { // define function\r
+ fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v2\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 2.); // constant\r
+ } else if (!fHistIntV2 && fHistIntV3) {\r
+ fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v3\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 3.); // constant\r
+ } else if (fHistIntV2 && fHistIntV3) {\r
+ fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());\r
+ fFuncVn->SetParameter(0, 1.); // normalization\r
+ fFuncVn->SetParameter(3, 0.2); // v2\r
+ fFuncVn->FixParameter(1, 1.); // constant\r
+ fFuncVn->FixParameter(2, 2.); // constant\r
+ fFuncVn->FixParameter(5, 3.); // constant\r
+ fFuncVn->SetParameter(7, 0.2); // v3\r
+ }\r
+ // add the generator objects that have been added to the task\r
+ if(fTrackSpectrum) fOutputList->Add(fTrackSpectrum);\r
+ if(fJetSpectrumSF) fOutputList->Add(fJetSpectrumSF);\r
+ if(fHistIntV2) fOutputList->Add(fHistIntV2);\r
+ if(fHistIntV3) fOutputList->Add(fHistIntV3);\r
+ // create the QA histos\r
+ for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {\r
+ fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 100, 0, 100, i);\r
+ fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);\r
+ fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 100, 0, 100, i);\r
+ fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);\r
+ fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);\r
+ fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);\r
+ fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 100, 0, 20, 80, 0, .8, i);\r
+ // add to outputlist\r
+ if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);\r
+ if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);\r
+ if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);\r
+ if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);\r
+ }\r
+ if(fJetSpectrumSF) {\r
+ fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);\r
+ fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());\r
+ }\r
+ fOutputList->Sort();\r
+ PostData(1, fOutputList);\r
+}\r
+//_____________________________________________________________________________\r
+TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)\r
+{\r
+ // book a TH1F and connect it to the output container\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!fOutputList) return 0x0;\r
+ TString title(name);\r
+ if(c!=-1) { // format centrality dependent histograms accordingly\r
+ name = Form("%s_%i", name, c);\r
+ title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));\r
+ }\r
+ title += Form(";%s;[counts]", x);\r
+ TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);\r
+ histogram->Sumw2();\r
+ if(append) fOutputList->Add(histogram);\r
+ return histogram; \r
+}\r
+//_____________________________________________________________________________\r
+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)\r
+{\r
+ // book a TH2F and connect it to the output container\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!fOutputList) return 0x0;\r
+ TString title(name);\r
+ if(c!=-1) { // format centrality dependent histograms accordingly\r
+ name = Form("%s_%i", name, c);\r
+ title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));\r
+ }\r
+ title += Form(";%s;%s", x, y);\r
+ TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);\r
+ histogram->Sumw2();\r
+ if(append) fOutputList->Add(histogram);\r
+ return histogram; \r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::UserExec(Option_t *) \r
+{\r
+ // user exec, called for each event.\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); if(!AliAnalysisManager::GetAnalysisManager()) return;\r
+ // retrieve tracks from input.\r
+ if (!fTracksIn) { \r
+ fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));\r
+ if (!fTracksIn) {\r
+ AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data())); \r
+ return;\r
+ }\r
+ }\r
+ // get the centrality bin for qa plots\r
+ Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));\r
+ fCenBin = -1;\r
+ for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {\r
+ if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {\r
+ fCenBin = i;\r
+ break; }\r
+ }\r
+ if(fCenBin < 0) return;\r
+ // add tracks to event if not yet there\r
+ fTracksOut->Delete();\r
+ if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);\r
+ fTracksOut->Clear();\r
+ // get the event plane\r
+ CalculateEventPlane();\r
+ const Int_t Ntracks = fTracksIn->GetEntriesFast();\r
+ Int_t nacc(0);\r
+ for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {\r
+ AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));\r
+ if(!track) continue;\r
+ Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(track->Eta());\r
+ // fill qa histo's before applying any (possible) afterburner\r
+ FillHistogramsOriginalData(pt, eta, phi);\r
+ if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);\r
+ else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);\r
+ else if(fHistDiffV2 || fHistDiffV3) SampleVnFromTF1(phi); \r
+ /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1); \r
+ nacc++;\r
+ }\r
+ if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);\r
+ PostData(1, fOutputList);\r
+ if(fDebug > 0) PrintInfo();\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const\r
+{\r
+ // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations\r
+ if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); phi = PhaseShift(gRandom->Uniform(0, TMath::TwoPi()) + fPsi2);\r
+ Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);\r
+ if(TMath::AreEqualAbs(v2, 0, 1e-5)) { \r
+ FillHistogramsToyData(pt, eta, phi, v2);\r
+ return;\r
+ }\r
+ // introduce flow fluctuations (gaussian)\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v2);\r
+ for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {\r
+ phiprev=phi; //store last value for comparison\r
+ f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));\r
+ fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative\r
+ phi -= f/fp;\r
+ if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break; \r
+ }\r
+ FillHistogramsToyData(pt, eta, phi, v2);\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const\r
+{\r
+ // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations\r
+ if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); phi = PhaseShift(gRandom->Uniform(0, TMath::TwoPi()) + fPsi3);\r
+ Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);\r
+ if(TMath::AreEqualAbs(v3, 0, 1e-5)) {\r
+ FillHistogramsToyData(pt, eta, phi, v3);\r
+ return;\r
+ }\r
+ // introduce flow fluctuations (gaussian)\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v3);\r
+ for (Int_t i=0; i<fMaxNumberOfIterations; i++) {\r
+ phiprev=phi; //store last value for comparison\r
+ f = phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));\r
+ fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative\r
+ phi -= f/fp;\r
+ if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;\r
+ }\r
+ FillHistogramsToyData(pt, eta, phi, v3);\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)\r
+{\r
+ // inject single fragmentation jet spectrum to the tclones array, note that emcal params \r
+ // equal the barrel kinematics to pass the track and jet cuts later on\r
+ if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__); for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {\r
+ Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());\r
+ /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);\r
+ fHistSFJetSpectrum->Fill(pt);\r
+ fHistSFJetEtaPhi->Fill(eta, phi);\r
+ ++i;\r
+ }\r
+ nacc = 0;\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {\r
+ // grab the event plane orientation from the AliVEvent header\r
+ 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);\r
+ switch (fDetectorType) {\r
+ case kVZEROA : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);\r
+ break;\r
+ }\r
+ case kVZEROC : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);\r
+ break;\r
+ }\r
+ case kVZEROComb : {\r
+ fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;\r
+ fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);\r
+ break;\r
+ }\r
+ default : break;\r
+ }\r
+ // if requested, pass the event plane to the phi distribution\r
+ if(fHistIntV2 && fHistIntV3) {\r
+ fFuncVn->SetParameter(4, fPsi2); fFuncVn->SetParameter(6, fPsi3);\r
+ Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) {\r
+ GetFlowFluctuation(v2); GetFlowFluctuation(v3);\r
+ }\r
+ fFuncVn->SetParameter(3, v2); fFuncVn->SetParameter(7, v3);\r
+ } else if (fHistIntV2) {\r
+ fFuncVn->SetParameter(4, fPsi2);\r
+ Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v2);\r
+ fFuncVn->SetParameter(3, v2);\r
+ } else if (fHistIntV3) {\r
+ fFuncVn->SetParameter(4, fPsi3);\r
+ Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));\r
+ if(fFlowFluctuations > -10) GetFlowFluctuation(v3);\r
+ fFuncVn->SetParameter(3, v3);\r
+ }\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)\r
+{\r
+ // terminate\r
+}\r
+//_____________________________________________________________________________\r
+void AliAnalysisTaskJetFlowMC::PrintInfo() const\r
+{\r
+ // print info\r
+ printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());\r
+ printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());\r
+\r
+}\r
+//_____________________________________________________________________________\r