--- /dev/null
+// Analysis type can be ESD, AOD, MC, ESDMCkineESD, ESDMCkineMC
+// in this simple example only ESD is supported
+const TString type = "ESD";
+// tracks used to determine the Q vector (Global, Tracklet or FMD at the moment)
+const TString rptype = "Global";
+// Boolean to fill/not fill the QA histograms
+Bool_t QA = kTRUE;
+
+// Boolean to use/not use weights for the Q vector
+Bool_t WEIGHTS[] = {kFALSE,kFALSE,kFALSE}; //Phi, v'(pt), v'(eta)
+
+
+void runFlowTaskExample(Int_t nRuns = 2,
+ Bool_t DATA = kFALSE, const Char_t* dataDir="/Users/snelling/alice_data/Therminator_midcentral", Int_t offset = 0)
+{
+ TStopwatch timer;
+ timer.Start();
+
+ //--------------------------------------
+ // Load the needed libraries most of them already loaded by aliroot
+ //--------------------------------------
+ gSystem->Load("libTree");
+ gSystem->Load("libGeom");
+ gSystem->Load("libVMC");
+ gSystem->Load("libXMLIO");
+ gSystem->Load("libPhysics");
+ gSystem->Load("libSTEERBase");
+ gSystem->Load("libESD");
+ gSystem->Load("libAOD");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libCORRFW");
+ gSystem->Load("libPWG2flowCommon");
+ gSystem->Load("libPWG2flowTasks");
+
+
+ if (type!="AOD") { TChain* chain = CreateESDChain(dataDir, nRuns, offset);}
+ else { TChain* chain = CreateAODChain(dataDir, nRuns, offset);}
+
+
+ // CORRFW correction framework cuts
+ //----------Event cuts----------
+
+ TObjArray* mcEventList = new TObjArray(0);
+ AliCFEventGenCuts* mcEventCuts = new AliCFEventGenCuts("mcEventCuts","MC-level event cuts");
+ mcEventList->AddLast(mcEventCuts);
+
+ TObjArray* recEventList = new TObjArray(0);
+ AliCFEventRecCuts* recEventCuts = new AliCFEventRecCuts("recEventCuts","rec-level event cuts");
+ recEventList->AddLast(recEventCuts);
+
+ //----------Cuts for RP----------
+ TObjArray* mcListRP = new TObjArray(0);
+ AliCFTrackKineCuts *mcKineCutsRP = new AliCFTrackKineCuts("mcKineCutsRP","MC-level kinematic cuts");
+ mcListRP->AddLast(mcKineCutsRP);
+ AliCFParticleGenCuts *mcGenCutsRP = new AliCFParticleGenCuts("mcGenCutsRP","MC particle generation cuts for RP");
+ mcGenCutsRP->SetRequireIsPrimary();
+ mcListRP->AddLast(mcGenCutsRP);
+
+ TObjArray* fPIDCutListRP = new TObjArray(0) ;
+ AliCFTrackCutPid* cutPidRP = NULL;
+ fPIDCutListRP->AddLast(cutPidRP);
+
+ TObjArray* recListRP = new TObjArray(0) ;
+ AliCFTrackKineCuts *recKineCutsRP = new AliCFTrackKineCuts("recKineCutsRP","rec-level kine cuts");
+ recListRP->AddLast(recKineCutsRP);
+ AliCFTrackQualityCuts *recQualityCutsRP = new AliCFTrackQualityCuts("recQualityCutsRP","rec-level quality cuts");
+ recQualityCutsRP->SetMinNClusterTPC(70);
+ recListRP->AddLast(recQualityCutsRP);
+ AliCFTrackIsPrimaryCuts *recIsPrimaryCutsRP = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsRP","rec-level isPrimary cuts");
+ recIsPrimaryCutsRP->UseSPDvertex(kTRUE);
+ recIsPrimaryCutsRP->UseTPCvertex(kTRUE);
+ recListRP->AddLast(recIsPrimaryCutsRP);
+
+ TObjArray* accListRP = new TObjArray(0) ;
+ AliCFAcceptanceCuts *mcAccCutsRP = new AliCFAcceptanceCuts("mcAccCutsRP","MC acceptance cuts");
+ mcAccCutsRP->SetMinNHitITS(2);
+ accListRP->AddLast(mcAccCutsRP);
+
+ //----------Cuts for POI----------
+ TObjArray* mcListPOI = new TObjArray(0);
+ AliCFTrackKineCuts *mcKineCutsPOI = new AliCFTrackKineCuts("mcKineCutsPOI","MC-level kinematic cuts");
+ mcListPOI->AddLast(mcKineCutsPOI);
+ AliCFParticleGenCuts *mcGenCutsPOI = new AliCFParticleGenCuts("mcGenCutsPOI","MC particle generation cuts for POI");
+ mcGenCutsPOI->SetRequireIsPrimary();
+ mcListPOI->AddLast(mcGenCutsPOI);
+
+ TObjArray* fPIDCutListPOI = new TObjArray(0) ;
+ AliCFTrackCutPid* cutPidPOI = NULL;
+ fPIDCutListPOI->AddLast(cutPidPOI);
+
+ TObjArray* recListPOI = new TObjArray(0) ;
+ AliCFTrackKineCuts *recKineCutsPOI = new AliCFTrackKineCuts("recKineCutsPOI","rec-level kine cuts");
+ recListPOI->AddLast(recKineCutsPOI);
+ AliCFTrackQualityCuts *recQualityCutsPOI = new AliCFTrackQualityCuts("recQualityCutsPOI","rec-level quality cuts");
+ recQualityCutsPOI->SetMinNClusterTPC(70);
+ recListPOI->AddLast(recQualityCutsPOI);
+ AliCFTrackIsPrimaryCuts *recIsPrimaryCutsPOI = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsPOI","rec-level isPrimary cuts");
+ recIsPrimaryCutsPOI->UseSPDvertex(kTRUE);
+ recIsPrimaryCutsPOI->UseTPCvertex(kTRUE);
+ recListPOI->AddLast(recIsPrimaryCutsPOI);
+
+ TObjArray* accListPOI = new TObjArray(0) ;
+ AliCFAcceptanceCuts *mcAccCutsPOI = new AliCFAcceptanceCuts("mcAccCutsPOI","MC acceptance cuts");
+ mcAccCutsPOI->SetMinNHitITS(2);
+ accListPOI->AddLast(mcAccCutsPOI);
+
+ //----------Add Cut Lists to the CF Manager----------
+ printf("CREATE INTERFACE AND CUTS\n");
+ AliCFManager* cfmgrRP = new AliCFManager();
+ cfmgrRP->SetNStepEvent(3);
+ cfmgrRP->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList);
+ cfmgrRP->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList);
+ cfmgrRP->SetNStepParticle(4);
+ cfmgrRP->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListRP);
+ cfmgrRP->SetParticleCutsList(AliCFManager::kPartAccCuts,accListRP);
+ cfmgrRP->SetParticleCutsList(AliCFManager::kPartRecCuts,recListRP);
+ cfmgrRP->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListRP);
+
+ AliCFManager* cfmgrPOI = new AliCFManager();
+ cfmgrPOI->SetNStepEvent(3);
+ cfmgrPOI->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList);
+ cfmgrPOI->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList);
+ cfmgrPOI->SetNStepParticle(4);
+ cfmgrPOI->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListPOI);
+ cfmgrPOI->SetParticleCutsList(AliCFManager::kPartAccCuts,accListPOI);
+ cfmgrPOI->SetParticleCutsList(AliCFManager::kPartRecCuts,recListPOI);
+ cfmgrPOI->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListPOI);
+
+
+ // Make the analysis manager
+ AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
+
+ AliVEventHandler* esdH = new AliESDInputHandler;
+ mgr->SetInputEventHandler(esdH);
+ AliMCEventHandler *mc = new AliMCEventHandler();
+ mgr->SetMCtruthEventHandler(mc);
+
+
+ // gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ // AliPhysicsSelectionTask* physicsSelTask = AddTaskPhysicsSelection();
+ // if (!DATA) {physicsSelTask->GetPhysicsSelection()->SetAnalyzeMC();}
+
+ // Add the task which makes the flowevent to the analysis manager
+ AliAnalysisTaskFlowEvent* taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",rptype,kFALSE,1);
+ taskFE->SetAnalysisType(type);
+ // taskFE->SelectCollisionCandidates();
+ mgr->AddTask(taskFE);
+
+ // Set the cuts for the flowevent
+ taskFE->SetCFManager1(cfmgrRP);
+ taskFE->SetCFManager2(cfmgrPOI);
+
+ // Instantiate the flow methods and add to the analysis manager
+ AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane");
+ mgr->AddTask(taskMCEP);
+ AliAnalysisTaskQCumulants *taskQC = new AliAnalysisTaskQCumulants("TaskQCumulants",kFALSE);
+ mgr->AddTask(taskQC);
+
+ // Create the output container for the data produced by the task
+ // Connect to the input and output containers
+ //===========================================================================
+ AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+
+ AliAnalysisDataContainer *coutputFE =
+ mgr->CreateContainer("cobjFlowEventSimple", AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
+ mgr->ConnectInput(taskFE,0,cinput1);
+ mgr->ConnectOutput(taskFE,1,coutputFE);
+
+
+ TString outputMCEP = AliAnalysisManager::GetCommonFileName();
+ outputMCEP += ":outputMCEPanalysis";
+ outputMCEP+= type;
+
+ AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer("cobjMCEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP);
+ mgr->ConnectInput(taskMCEP,0,coutputFE);
+ mgr->ConnectOutput(taskMCEP,1,coutputMCEP);
+
+
+ TString outputQC = AliAnalysisManager::GetCommonFileName();
+ outputQC += ":outputQCanalysis";
+ outputQC+= type;
+
+ AliAnalysisDataContainer *coutputQC = mgr->CreateContainer("cobjQC", TList::Class(),AliAnalysisManager::kOutputContainer,outputQC);
+ mgr->ConnectInput(taskQC,0,coutputFE);
+ mgr->ConnectOutput(taskQC,1,coutputQC);
+
+ // Run the analysis
+ if (!mgr->InitAnalysis()) return;
+ mgr->PrintStatus();
+ mgr->StartAnalysis("local",chain);
+
+ timer.Stop();
+ timer.Print();
+
+}
+
+// Helper macros for creating chains
+// from: CreateESDChain.C,v 1.10 jgrosseo Exp
+
+TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+ // creates chain of files in a given directory or file containing a list.
+ // In case of directory the structure is expected as:
+ // <aDataDir>/<dir0>/AliESDs.root
+ // <aDataDir>/<dir1>/AliESDs.root
+ // ...
+
+ if (!aDataDir)
+ return 0;
+
+ Long_t id, size, flags, modtime;
+ if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+ {
+ printf("%s not found.\n", aDataDir);
+ return 0;
+ }
+
+ TChain* chain = new TChain("esdTree");
+ TChain* chaingAlice = 0;
+
+ if (flags & 2)
+ {
+ TString execDir(gSystem->pwd());
+ TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+ TList* dirList = baseDir->GetListOfFiles();
+ Int_t nDirs = dirList->GetEntries();
+ gSystem->cd(execDir);
+
+ Int_t count = 0;
+
+ for (Int_t iDir=0; iDir<nDirs; ++iDir)
+ {
+ TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+ if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+ continue;
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ TString presentDirName(aDataDir);
+ presentDirName += "/";
+ presentDirName += presentDir->GetName();
+ chain->Add(presentDirName + "/AliESDs.root/esdTree");
+ // cerr<<presentDirName<<endl;
+ }
+
+ }
+ else
+ {
+ // Open the input stream
+ ifstream in;
+ in.open(aDataDir);
+
+ Int_t count = 0;
+
+ // Read the input list of files and add them to the chain
+ TString esdfile;
+ while(in.good()) {
+ in >> esdfile;
+ if (!esdfile.Contains("root")) continue; // protection
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ // add esd file
+ chain->Add(esdfile);
+ }
+
+ in.close();
+ }
+
+ return chain;
+}
+
+
+// Helper macros for creating chains
+// from: CreateESDChain.C,v 1.10 jgrosseo Exp
+
+TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+ // creates chain of files in a given directory or file containing a list.
+ // In case of directory the structure is expected as:
+ // <aDataDir>/<dir0>/AliAOD.root
+ // <aDataDir>/<dir1>/AliAOD.root
+ // ...
+
+ if (!aDataDir)
+ return 0;
+
+ Long_t id, size, flags, modtime;
+ if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+ {
+ printf("%s not found.\n", aDataDir);
+ return 0;
+ }
+
+ TChain* chain = new TChain("aodTree");
+ TChain* chaingAlice = 0;
+
+ if (flags & 2)
+ {
+ TString execDir(gSystem->pwd());
+ TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+ TList* dirList = baseDir->GetListOfFiles();
+ Int_t nDirs = dirList->GetEntries();
+ gSystem->cd(execDir);
+
+ Int_t count = 0;
+
+ for (Int_t iDir=0; iDir<nDirs; ++iDir)
+ {
+ TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+ if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+ continue;
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ TString presentDirName(aDataDir);
+ presentDirName += "/";
+ presentDirName += presentDir->GetName();
+ chain->Add(presentDirName + "/AliAOD.root/aodTree");
+ // cerr<<presentDirName<<endl;
+ }
+
+ }
+ else
+ {
+ // Open the input stream
+ ifstream in;
+ in.open(aDataDir);
+
+ Int_t count = 0;
+
+ // Read the input list of files and add them to the chain
+ TString aodfile;
+ while(in.good()) {
+ in >> aodfile;
+ if (!aodfile.Contains("root")) continue; // protection
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ // add aod file
+ chain->Add(aodfile);
+ }
+
+ in.close();
+ }
+
+ return chain;
+}
+
--- /dev/null
+void runStarFlowExample( Int_t maxNumberOfEvents = 1000, const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
+{
+ gSystem->Load("libTree.so");
+ gSystem->Load("libVMC.so");
+ gSystem->Load("libPhysics.so");
+ gSystem->Load("libPWG2flowCommon");
+
+ TStopwatch timer;
+ timer.Start();
+
+ //define reference particles
+ AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
+ rpCuts->SetPtMin(0.05);
+ rpCuts->SetPtMax(10.);
+
+ //define particles of interest
+ AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
+ poiCuts->SetPtMin(0.05);
+ poiCuts->SetPtMax(10.);
+
+ //define event cuts
+ AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
+ starEventCuts-> SetCentralityIDMax(3);
+ starEventCuts-> SetCentralityIDMin(3);
+
+ AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
+ AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants(); //mcep->SetHarmonic(2); // default is v2
+
+ mcep->Init(); qc->Init();
+
+ Int_t i=0;
+ AliStarEventReader starReader(inputDataFiles) ;
+ while ( starReader.GetNextEvent() ) // Get next event
+ {
+ AliStarEvent* starEvent = starReader.GetEvent();
+ if ( !starEventCuts->PassesCuts(starEvent) ) continue; // Test if the event is good
+ AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts); // make a flow event from a star event
+ // do flow analysis for various methods
+
+ mcep->Make(flowEvent);
+ qc->Make(flowEvent);
+
+ delete flowEvent;
+
+ i++;
+ cout <<"Event: " << i << "\r"; cout.flush();
+ if (i>maxNumberOfEvents) break;
+ }
+
+ //---------------------------------------------------------------------------------------
+ // create a new file which will hold the final results of all methods:
+ TString outputFileName = "AnalysisResults.root";
+ TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
+ // create a new file for each method wich will hold list with final results:
+ const Int_t nMethods = 2;
+ TString method[nMethods] = {"MCEP","QC"};
+ TDirectoryFile *dirFileFinal[nMethods] = {NULL};
+ TString fileName[nMethods];
+ for(Int_t i=0;i<nMethods;i++)
+ {
+ // form a file name for each method:
+ fileName[i]+="output";
+ fileName[i]+=method[i].Data();
+ fileName[i]+="analysis";
+ dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
+ }
+
+ // calculating and storing the final results of default flow analysis:
+ mcep->Finish(); mcep->WriteHistograms(dirFileFinal[0]);
+ qc->Finish(); qc->WriteHistograms(dirFileFinal[1]);
+ //---------------------------------------------------------------------------------------
+
+ outputFile->Close();
+ delete outputFile;
+
+ timer.Stop();
+ cout << endl;
+ timer.Print();
+
+ delete rpCuts;
+ delete poiCuts;
+ delete starEventCuts;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff