#include "TTree.h"
#include "TFile.h" //needed as include
#include "TList.h"
+#include "TH2F.h"
#include "TRandom3.h"
#include "TTimeStamp.h"
// ALICE Analysis Framework
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskSE.h"
+#include "AliESDtrackCuts.h"
// ESD interface
#include "AliESDEvent.h"
#include "AliGenGeVSimEventHeader.h"
#include "AliGenEposEventHeader.h"
+// Interface to Load short life particles
+#include "TObjArray.h"
+#include "AliFlowCandidateTrack.h"
+
// Interface to make the Flow Event Simple used in the flow analysis methods
#include "AliFlowEvent.h"
-
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventCuts.h"
+#include "AliFlowCommonConstants.h"
#include "AliAnalysisTaskFlowEvent.h"
#include "AliLog.h"
AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
AliAnalysisTaskSE(),
// fOutputFile(NULL),
- fAnalysisType("ESD"),
+ fAnalysisType("AUTOMATIC"),
+ fRPType("Global"),
fCFManager1(NULL),
fCFManager2(NULL),
- fQAInt(NULL),
- fQADiff(NULL),
+ fCutsEvent(NULL),
+ fCutsRP(NULL),
+ fCutsPOI(NULL),
+ fQAList(NULL),
fMinMult(0),
fMaxMult(10000000),
fMinA(-1.0),
fMaxA(-0.01),
fMinB(0.01),
fMaxB(1.0),
- fQA(kFALSE),
- fMCReactionPlaneAngle(0.),
- fCount(0),
- fNoOfLoops(1),
- fEllipticFlowValue(0.),
- fSigmaEllipticFlowValue(0.),
- fMultiplicityOfEvent(1000000000),
- fSigmaMultiplicityOfEvent(0),
- fMyTRandom3(NULL),
- fbAfterburnerOn(kFALSE)
+ fQAon(kFALSE),
+ fLoadCandidates(kFALSE),
+ fNbinsMult(10000),
+ fNbinsPt(100),
+ fNbinsPhi(100),
+ fNbinsEta(200),
+ fNbinsQ(500),
+ fMultMin(0.),
+ fMultMax(10000.),
+ fPtMin(0.),
+ fPtMax(10.),
+ fPhiMin(0.),
+ fPhiMax(TMath::TwoPi()),
+ fEtaMin(-5.),
+ fEtaMax(5.),
+ fQMin(0.),
+ fQMax(3.),
+ fExcludedEtaMin(0.),
+ fExcludedEtaMax(0.),
+ fExcludedPhiMin(0.),
+ fExcludedPhiMax(0.),
+ fAfterburnerOn(kFALSE),
+ fNonFlowNumberOfTrackClones(0),
+ fV1(0.),
+ fV2(0.),
+ fV3(0.),
+ fV4(0.),
+ fV5(0.),
+ fFlowEvent(NULL),
+ fMyTRandom3(NULL)
{
// Constructor
cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
}
//________________________________________________________________________
-AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed) :
+AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
AliAnalysisTaskSE(name),
// fOutputFile(NULL),
- fAnalysisType("ESD"),
+ fAnalysisType("AUTOMATIC"),
+ fRPType(RPtype),
fCFManager1(NULL),
fCFManager2(NULL),
- fQAInt(NULL),
- fQADiff(NULL),
+ fCutsEvent(NULL),
+ fCutsRP(NULL),
+ fCutsPOI(NULL),
+ fQAList(NULL),
fMinMult(0),
fMaxMult(10000000),
fMinA(-1.0),
fMaxA(-0.01),
fMinB(0.01),
fMaxB(1.0),
- fQA(on),
- fMCReactionPlaneAngle(0.),
- fCount(0),
- fNoOfLoops(1),
- fEllipticFlowValue(0.),
- fSigmaEllipticFlowValue(0.),
- fMultiplicityOfEvent(1000000000),
- fSigmaMultiplicityOfEvent(0),
- fMyTRandom3(NULL),
- fbAfterburnerOn(kFALSE)
+ fQAon(on),
+ fLoadCandidates(bCandidates),
+ fNbinsMult(10000),
+ fNbinsPt(100),
+ fNbinsPhi(100),
+ fNbinsEta(200),
+ fNbinsQ(500),
+ fMultMin(0.),
+ fMultMax(10000.),
+ fPtMin(0.),
+ fPtMax(10.),
+ fPhiMin(0.),
+ fPhiMax(TMath::TwoPi()),
+ fEtaMin(-5.),
+ fEtaMax(5.),
+ fQMin(0.),
+ fQMax(3.),
+ fExcludedEtaMin(0.),
+ fExcludedEtaMax(0.),
+ fExcludedPhiMin(0.),
+ fExcludedPhiMax(0.),
+ fAfterburnerOn(kFALSE),
+ fNonFlowNumberOfTrackClones(0),
+ fV1(0.),
+ fV2(0.),
+ fV3(0.),
+ fV4(0.),
+ fV5(0.),
+ fFlowEvent(NULL),
+ fMyTRandom3(NULL)
{
// Constructor
cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
fMyTRandom3 = new TRandom3(iseed);
gRandom->SetSeed(fMyTRandom3->Integer(65539));
+ int availableINslot=1;
+ //FMD input slot
+ if (strcmp(RPtype,"FMD")==0) {
+ DefineInput(availableINslot++, TList::Class());
+ }
+ //Candidates input slot
+ if( fLoadCandidates )
+ DefineInput(availableINslot, TObjArray::Class());
// Define output slots here
// Define here the flow event output
DefineOutput(1, AliFlowEventSimple::Class());
- if(on)
- {
- DefineOutput(2, TList::Class());
- DefineOutput(3, TList::Class());
- }
+ DefineOutput(2, TList::Class());
+
// and for testing open an output file
// fOutputFile = new TFile("FlowEvents.root","RECREATE");
//
// Destructor
//
- if (fMyTRandom3) delete fMyTRandom3;
+ delete fMyTRandom3;
+ delete fFlowEvent;
+ delete fCutsEvent;
+ delete fCutsRP;
+ delete fCutsPOI;
+ delete fQAList;
// objects in the output list are deleted
// by the TSelector dtor (I hope)
}
+//________________________________________________________________________
+void AliAnalysisTaskFlowEvent::NotifyRun()
+{
+ //at the beginning of each new run
+ AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
+ if (!fESD) return;
+
+ Int_t run = fESD->GetRunNumber();
+ AliInfo(Form("Stariting run #%i",run));
+}
+
//________________________________________________________________________
void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
{
// Called at every worker node to initialize
cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
- if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC"))
+ if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
{
- AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD and MC are allowed.");
+ AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
exit(1);
}
+
+ //set the common constants
+ AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+ cc->SetNbinsMult(fNbinsMult);
+ cc->SetNbinsPt(fNbinsPt);
+ cc->SetNbinsPhi(fNbinsPhi);
+ cc->SetNbinsEta(fNbinsEta);
+ cc->SetNbinsQ(fNbinsQ);
+ cc->SetMultMin(fMultMin);
+ cc->SetMultMax(fMultMax);
+ cc->SetPtMin(fPtMin);
+ cc->SetPtMax(fPtMax);
+ cc->SetPhiMin(fPhiMin);
+ cc->SetPhiMax(fPhiMax);
+ cc->SetEtaMin(fEtaMin);
+ cc->SetEtaMax(fEtaMax);
+ cc->SetQMin(fQMin);
+ cc->SetQMax(fQMax);
+
+ fFlowEvent = new AliFlowEvent(3000);
+
+ if (fQAon)
+ {
+ fQAList=new TList();
+ fQAList->SetName(Form("%s QA",GetName()));
+ if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA());
+ if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA());
+ if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA());
+ PostData(2,fQAList);
+ }
}
//________________________________________________________________________
{
// Main loop
// Called for each event
- AliFlowEvent* flowEvent = NULL;
+ //delete fFlowEvent;
AliMCEvent* mcEvent = MCEvent(); // from TaskSE
AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
+ AliMultiplicity* myTracklets = NULL;
+ AliESDPmdTrack* pmdtracks = NULL;//pmd
+
+ int availableINslot=1;
+
+ if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
+ {
+ AliError("cuts not set");
+ return;
+ }
+
+ //DEFAULT - automatically takes care of everything
+ if (fAnalysisType == "AUTOMATIC")
+ {
+ //check event cuts
+ if (!fCutsEvent->IsSelected(InputEvent())) return;
+
+ //first attach all possible information to the cuts
+ fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
+ fCutsPOI->SetEvent( InputEvent(), MCEvent() );
+
+ //then make the event
+ fFlowEvent->Fill( fCutsRP, fCutsPOI );
+ //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
+
+ if (myESD)
+ fFlowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent()));
+ if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+ }
// Make the FlowEvent for MC input
- if (fAnalysisType == "MC")
+ else if (fAnalysisType == "MC")
{
// Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
// This handler can return the current MC event
AliWarning("Event does not pass multiplicity cuts"); return;
}
//make event
- flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
+ fFlowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
}
+
// Make the FlowEvent for ESD input
else if (fAnalysisType == "ESD")
{
}
//check the offline trigger (check if the event has the correct trigger)
- AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
+ AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
//check multiplicity
if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
}
//make event
- flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
+ if (fRPType == "Global") {
+ fFlowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
+ }
+ else if (fRPType == "TPCOnly") {
+ fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
+ }
+ else if (fRPType == "TPCHybrid") {
+ fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
+ }
+ else if (fRPType == "Tracklet"){
+ fFlowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
+ }
+ else if (fRPType == "FMD"){
+ TList* FMDdata = dynamic_cast<TList*>(GetInputData(availableINslot++));
+ if(!FMDdata) {
+ cout<<" No FMDdata "<<endl;
+ return;
+ }
+ TH2F* histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
+ if (!histFMD) {
+ cout<< "No histFMD"<<endl;
+ return;
+ }
+ fFlowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
+ }
+ else if (fRPType == "PMD"){
+ fFlowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
+ }
+ else return;
+
+ // if monte carlo event get reaction plane from monte carlo (depends on generator)
+ if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+ //set reference multiplicity, TODO: maybe move it to the constructor?
+ fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
}
+
// Make the FlowEvent for ESD input combined with MC info
else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
{
//make event
if (fAnalysisType == "ESDMCkineESD")
{
- flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
+ fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
}
else if (fAnalysisType == "ESDMCkineMC")
{
- flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
+ fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
}
+ if (!fFlowEvent) return;
+ // if monte carlo event get reaction plane from monte carlo (depends on generator)
+ if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+ //set reference multiplicity, TODO: maybe move it to the constructor?
+ fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
}
// Make the FlowEventSimple for AOD input
else if (fAnalysisType == "AOD")
AliError("ERROR: AOD not available");
return;
}
- AliInfo(Form("There are %d tracks in this event", myAOD->GetNumberOfTracks()));
- flowEvent = new AliFlowEvent(myAOD);
+ AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
+ fFlowEvent = new AliFlowEvent(myAOD);
}
+ //inject candidates
+ if(fLoadCandidates) {
+ TObjArray* Candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
+ //if(Candidates->GetEntriesFast())
+ //printf("I received %d candidates\n",Candidates->GetEntriesFast());
+ if (Candidates)
+ {
+ for(int iCand=0; iCand!=Candidates->GetEntriesFast(); ++iCand ) {
+ AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(Candidates->At(iCand));
+ if (!cand) continue;
+ cand->SetForPOISelection(kTRUE);
+ cand->SetForRPSelection(kFALSE);
+ //printf(" Ⱶ Checking at candidate %d with %d daughters\n",iCand,cand->GetNDaughters());
+ for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
+ //printf(" Ⱶ Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
+ for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs ) {
+ AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
+ if (!iRP) continue;
+ if( !iRP->InRPSelection() )
+ continue;
+ if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
+ //printf(" was in RP set");
+ cand->SetDaughter( iDau, iRP );
+ //temporarily untagging all daugters
+ iRP->SetForRPSelection(kFALSE);
+ }
+ }
+ //printf("\n");
+ }
+ fFlowEvent->AddTrack(cand);
+ }
+ }
+ }
+
+ if (!fFlowEvent) return; //shuts up coverity
+
//check final event cuts
- Int_t mult = flowEvent->NumberOfTracks();
- if (mult<fMinMult && mult>fMaxMult)
+ Int_t mult = fFlowEvent->NumberOfTracks();
+ AliInfo(Form("FlowEvent has %i tracks",mult));
+ if (mult<fMinMult || mult>fMaxMult)
{
AliWarning("FlowEvent cut on multiplicity"); return;
}
- //tag subevents
- flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
-
- ////TODO afterburner
- //if (fbAfterburnerOn && fMyTRandom3) {
- // // set the new value of the values using a after burner
- // cout << "settings for afterburner in TaskFlowEvent.cxx:" << endl;
- // cout << "fCount = " << fCount << endl;
- // cout << "fNoOfLoops = " << fNoOfLoops << endl;
- // cout << "fEllipticFlowValue = " << fEllipticFlowValue << endl;
- // cout << "fSigmaEllipticFlowValue = " << fSigmaEllipticFlowValue << endl;
- // cout << "fMultiplicityOflowEvent = " << fMultiplicityOflowEvent << endl;
- // cout << "fSigmaMultiplicityOflowEvent = " << fSigmaMultiplicityOflowEvent << endl;
- // Double_t xRPangle;
- // if (!mcEvent) { xRPangle = TMath::TwoPi()*(fMyTRandom3->Rndm()); }
- // else { xRPangle = fRP; }
- // Double_t xNewFlowValue = fMyTRandom3->Gaus(fEllipticFlowValue,fSigmaEllipticFlowValue);
- // Int_t nNewMultOflowEvent = TMath::Nint(fMyTRandom3->Gaus(fMultiplicityOflowEvent,fSigmaMultiplicityOflowEvent));
- // cout << "xRPangle = " << xRPangle << endl;
- // cout << "xNewFlowValue = " << xNewFlowValue << endl;
- // cout << "nNewMultOflowEvent = " << nNewMultOflowEvent << endl;
- // cout << "settings for after burner" << endl;
- // flowEventMaker->SetMCReactionPlaneAngle(xRPangle);
- // flowEventMaker->SetNoOfLoops(fNoOfLoops);
- // flowEventMaker->SetEllipticFlowValue(xNewFlowValue);
- // flowEventMaker->SetMultiplicityOflowEvent(nNewMultOflowEvent);
- // //end settings afterburner
- //}
-
- // if monte carlo event get reaction plane from monte carlo (depends on generator)
- if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
+ //define dead zone
+ fFlowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
- //fListHistos->Print();
- // fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
- PostData(1,flowEvent);
- if (fQA)
+
+ //////////////////////////////////////////////////////////////////////////////
+ ///////////////////////////AFTERBURNER
+ if (fAfterburnerOn)
{
- PostData(2,fQAInt);
- PostData(3,fQADiff);
+ //if reaction plane not set from elsewhere randomize it before adding flow
+ if (!fFlowEvent->IsSetMCReactionPlaneAngle())
+ fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
+
+ fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
+ fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
}
+ //////////////////////////////////////////////////////////////////////////////
+
+ //tag subEvents
+ fFlowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
+
+ //fListHistos->Print();
+ //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
+ PostData(1,fFlowEvent);
}
//________________________________________________________________________
// Called once at the end of the query -- do not call in case of CAF
}
-