1 // from CreateESDChain.C - instead of #include "CreateESDChain.C"
2 TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 20, Int_t offset = 0) ;
3 void LookupWrite(TChain* chain, const char* target) ;
7 //Flow analysis method can be:
9 // LYZ1 = Lee Yang Zeroes first run
10 // LYZ2 = Lee Yang Zeroes second run
11 // LYZEP = Lee Yang Zeroes Event Plane
13 const TString method = "SP";
15 //Type of analysis can be:
16 // ESD, AOD, MC, ESDMC0, ESDMC1
17 const TString type = "ESDMC1";
23 const Double_t ptmin1 = 0.0;
24 const Double_t ptmax1 = 1000.0;
25 const Double_t ymin1 = -2.;
26 const Double_t ymax1 = 2.;
27 const Int_t mintrackrefsTPC1 = 2;
28 const Int_t mintrackrefsITS1 = 3;
29 const Int_t charge1 = 1;
30 const Int_t PDG1 = 211;
31 const Int_t minclustersTPC1 = 50;
32 const Int_t maxnsigmatovertex1 = 3;
34 //for differential flow
35 const Double_t ptmin2 = 0.0;
36 const Double_t ptmax2 = 1000.0;
37 const Double_t ymin2 = -2.;
38 const Double_t ymax2 = 2.;
39 const Int_t mintrackrefsTPC2 = 2;
40 const Int_t mintrackrefsITS2 = 3;
41 const Int_t charge2 = 1;
42 const Int_t PDG2 = 321;
43 const Int_t minclustersTPC2 = 50;
44 const Int_t maxnsigmatovertex2 = 3;
48 void runAliAnalysisTaskFlow(Int_t nRuns = 10, const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX", Int_t offset = 0)
54 // include path (to find the .h files when compiling)
55 gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
56 gSystem->AddIncludePath("-I$ROOTSYS/include") ;
58 // load needed libraries
59 gSystem->Load("libTree.so");
60 gSystem->Load("libESD.so");
61 cerr<<"libESD loaded..."<<endl;
62 gSystem->Load("libANALYSIS.so");
63 cerr<<"libANALYSIS.so loaded..."<<endl;
64 gSystem->Load("libANALYSISRL.so");
65 cerr<<"libANALYSISRL.so loaded..."<<endl;
66 gSystem->Load("libCORRFW.so");
67 cerr<<"libCORRFW.so loaded..."<<endl;
68 gSystem->Load("libPWG2flow.so");
69 cerr<<"libPWG2flow.so loaded..."<<endl;
71 // create the TChain. CreateESDChain() is defined in CreateESDChain.C
72 TChain* chain = CreateESDChain(dataDir, nRuns, offset);
73 cout<<"chain ("<<chain<<")"<<endl;
75 //____________________________________________//
76 //Create cuts using correction framework
78 //############# cuts on MC
79 AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
80 mcKineCuts1->SetPtRange(ptmin1,ptmax1);
81 mcKineCuts1->SetRapidityRange(ymin1,ymax1);
82 mcKineCuts1->SetChargeMC(charge1);
84 AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
85 mcKineCuts2->SetPtRange(ptmin2,ptmax2);
86 mcKineCuts2->SetRapidityRange(ymin2,ymax2);
87 mcKineCuts2->SetChargeMC(charge2);
89 AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
90 mcGenCuts1->SetRequireIsPrimary();
91 mcGenCuts1->SetRequirePdgCode(PDG1);
93 AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
94 mcGenCuts2->SetRequireIsPrimary();
95 mcGenCuts2->SetRequirePdgCode(PDG2);
97 //############# Acceptance Cuts ????????
98 AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts","MC acceptance cuts");
99 mcAccCuts->SetMinNHitITS(mintrackrefsITS1);
100 mcAccCuts->SetMinNHitTPC(mintrackrefsTPC1);
102 //############# Rec-Level kinematic cuts
103 AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
104 recKineCuts1->SetPtRange(ptmin1,ptmax1);
105 recKineCuts1->SetRapidityRange(ymin1,ymax1);
106 recKineCuts1->SetChargeRec(charge1);
108 AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
109 recKineCuts2->SetPtRange(ptmin2,ptmax2);
110 recKineCuts2->SetRapidityRange(ymin2,ymax2);
111 recKineCuts2->SetChargeRec(charge2);
113 AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
114 recQualityCuts->SetMinNClusterTPC(minclustersTPC1);
115 recQualityCuts->SetRequireITSRefit(kTRUE);
117 AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
118 recIsPrimaryCuts->SetMaxNSigmaToVertex(maxnsigmatovertex1);
120 AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID") ;
121 AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
122 int n_species = AliPID::kSPECIES ;
123 Double_t* prior = new Double_t[n_species];
125 prior[0] = 0.0244519 ;
126 prior[1] = 0.0143988 ;
127 prior[2] = 0.805747 ;
128 prior[3] = 0.0928785 ;
129 prior[4] = 0.0625243 ;
131 cutPID1->SetPriors(prior);
132 cutPID1->SetProbabilityCut(0.0);
133 cutPID1->SetDetectors("TPC TOF");
134 switch(TMath::Abs(PDG1)) {
135 case 11 : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
136 case 13 : cutPID1->SetParticleType(AliPID::kMuon , kTRUE); break;
137 case 211 : cutPID1->SetParticleType(AliPID::kPion , kTRUE); break;
138 case 321 : cutPID1->SetParticleType(AliPID::kKaon , kTRUE); break;
139 case 2212 : cutPID1->SetParticleType(AliPID::kProton , kTRUE); break;
140 default : printf("UNDEFINED PID\n"); break;
143 cutPID2->SetPriors(prior);
144 cutPID2->SetProbabilityCut(0.0);
145 cutPID2->SetDetectors("TPC TOF");
146 switch(TMath::Abs(PDG2)) {
147 case 11 : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
148 case 13 : cutPID2->SetParticleType(AliPID::kMuon , kTRUE); break;
149 case 211 : cutPID2->SetParticleType(AliPID::kPion , kTRUE); break;
150 case 321 : cutPID2->SetParticleType(AliPID::kKaon , kTRUE); break;
151 case 2212 : cutPID2->SetParticleType(AliPID::kProton , kTRUE); break;
152 default : printf("UNDEFINED PID\n"); break;
156 printf("CREATE MC KINE CUTS\n");
157 TObjArray* mcList1 = new TObjArray(0);
158 mcList1->AddLast(mcKineCuts1);
159 mcList1->AddLast(mcGenCuts1);
161 TObjArray* mcList2 = new TObjArray(0);
162 mcList2->AddLast(mcKineCuts2);
163 mcList2->AddLast(mcGenCuts2);
165 printf("CREATE ACCEPTANCE CUTS\n");
166 TObjArray* accList = new TObjArray(0) ;
167 accList->AddLast(mcAccCuts);
169 printf("CREATE RECONSTRUCTION CUTS\n");
170 TObjArray* recList1 = new TObjArray(0) ;
171 recList1->AddLast(recKineCuts1);
172 recList1->AddLast(recQualityCuts);
173 recList1->AddLast(recIsPrimaryCuts);
175 TObjArray* recList2 = new TObjArray(0) ;
176 recList2->AddLast(recKineCuts2);
177 recList2->AddLast(recQualityCuts);
178 recList2->AddLast(recIsPrimaryCuts);
180 printf("CREATE PID CUTS\n");
181 TObjArray* fPIDCutList1 = new TObjArray(0) ;
182 fPIDCutList1->AddLast(cutPID1);
184 TObjArray* fPIDCutList2 = new TObjArray(0) ;
185 fPIDCutList2->AddLast(cutPID2);
187 printf("CREATE INTERFACE AND CUTS\n");
188 AliCFManager* cfmgr1 = new AliCFManager();
189 cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
190 //cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
191 cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);
192 cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);
194 AliCFManager* cfmgr2 = new AliCFManager();
195 cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
196 //cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
197 cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
198 cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
201 if (method == "LYZ2"){
202 // read the input files from the first run
203 TString firstRunFileName = "outputLYZ1analysis" ;
204 firstRunFileName += type;
205 firstRunFileName += "_firstrun.root";
206 cout<<"The input file is "<<firstRunFileName.Data()<<endl;
207 TFile* fFirstRunFile = new TFile(firstRunFileName.Data(),"READ");
208 if(!fFirstRunFile || fFirstRunFile->IsZombie()) { cerr << " ERROR: NO First Run file... " << endl ; }
209 cout<<"input files read..."<<endl;
213 //____________________________________________//
214 // Make the analysis manager
215 AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
218 AliVEventHandler* esdH = new AliESDInputHandler;
219 mgr->SetInputEventHandler(esdH); }
222 AliVEventHandler* aodH = new AliAODInputHandler;
223 mgr->SetInputEventHandler(aodH); }
225 if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
226 AliVEventHandler* esdH = new AliESDInputHandler;
227 mgr->SetInputEventHandler(esdH);
229 AliMCEventHandler *mc = new AliMCEventHandler();
230 mgr->SetMCtruthEventHandler(mc); }
232 //____________________________________________//
235 AliAnalysisTaskScalarProduct *task1 = new AliAnalysisTaskScalarProduct("TaskScalarProduct");
236 task1->SetAnalysisType(type);
237 task1->SetCFManager1(cfmgr1);
238 task1->SetCFManager2(cfmgr2);
241 else if (method == "LYZ1"){
242 AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE);
243 task1->SetAnalysisType(type);
244 task1->SetFirstRunLYZ(kTRUE);
245 task1->SetUseSumLYZ(kTRUE);
246 task1->SetCFManager1(cfmgr1);
247 task1->SetCFManager2(cfmgr2);
250 else if (method == "LYZ2"){
251 AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE);
252 task1->SetAnalysisType(type);
253 task1->SetFirstRunLYZ(kFALSE);
254 task1->SetUseSumLYZ(kTRUE);
255 task1->SetCFManager1(cfmgr1);
256 task1->SetCFManager2(cfmgr2);
259 else if (method == "LYZEP"){
260 AliAnalysisTaskLYZEventPlane *task1 = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane");
261 task1->SetAnalysisType(type);
262 //task1->SetCFManager1(cfmgr1);
263 //task1->SetCFManager2(cfmgr2);
266 else if (method == "CUM"){
267 AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
268 task1->SetAnalysisType(type);
269 //task1->SetCFManager1(cfmgr1);
270 //task1->SetCFManager2(cfmgr2);
275 // Create containers for input/output
276 AliAnalysisDataContainer *cinput1 =
277 mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
279 if (method == "LYZ2"){ AliAnalysisDataContainer *cinput2 =
280 mgr->CreateContainer("cobj2",TList::Class(),AliAnalysisManager::kInputContainer); }
283 TString outputName = "output";
285 outputName+= "analysis";
287 if (method == "LYZ1") outputName+= "_firstrun";
288 if (method == "LYZ2") outputName+= "_secondrun";
289 outputName+= ".root";
290 AliAnalysisDataContainer *coutput1 =
291 mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer,outputName);
293 //____________________________________________//
294 mgr->ConnectInput(task1,0,cinput1);
295 if (method == "LYZ2") { mgr->ConnectInput(task1,1,cinput2); }
296 mgr->ConnectOutput(task1,0,coutput1);
298 if (method == "LYZ2"){ cinput2->SetData(fFirstRunFile);}
300 if (!mgr->InitAnalysis()) return;
302 mgr->StartAnalysis("local",chain);
309 // Helper macros for creating chains
310 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
312 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
314 // creates chain of files in a given directory or file containing a list.
315 // In case of directory the structure is expected as:
316 // <aDataDir>/<dir0>/AliESDs.root
317 // <aDataDir>/<dir1>/AliESDs.root
323 Long_t id, size, flags, modtime;
324 if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
326 printf("%s not found.\n", aDataDir);
330 TChain* chain = new TChain("esdTree");
331 TChain* chaingAlice = 0;
335 TString execDir(gSystem->pwd());
336 TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
337 TList* dirList = baseDir->GetListOfFiles();
338 Int_t nDirs = dirList->GetEntries();
339 gSystem->cd(execDir);
343 for (Int_t iDir=0; iDir<nDirs; ++iDir)
345 TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
346 if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
355 if (count++ == aRuns)
358 TString presentDirName(aDataDir);
359 presentDirName += "/";
360 presentDirName += presentDir->GetName();
361 chain->Add(presentDirName + "/AliESDs.root/esdTree");
362 cerr<<presentDirName<<endl;
368 // Open the input stream
374 // Read the input list of files and add them to the chain
378 if (!esdfile.Contains("root")) continue; // protection
386 if (count++ == aRuns)
399 void LookupWrite(TChain* chain, const char* target)
401 // looks up the chain and writes the remaining files to the text file target
405 TObjArray* list = chain->GetListOfFiles();
406 TIterator* iter = list->MakeIterator();
410 outfile.open(target);
412 while ((obj = iter->Next()))
413 outfile << obj->GetTitle() << "#AliESDs.root" << endl;