PWGCF/FLOW/Documentation/examples/runFlowTaskExample.C

   1 // Analysis type can be ESD, AOD, MC, ESDMCkineESD, ESDMCkineMC
   2 // in this simple example only ESD is supported
   3 const TString type = "ESD";
   4 // tracks used to determine the Q vector (Global, Tracklet or FMD at the moment)
   5 const TString rptype = "Global";
   6 // Boolean to fill/not fill the QA histograms
   7 Bool_t QA = kTRUE;
   8
   9 // Boolean to use/not use weights for the Q vector
  10 Bool_t WEIGHTS[] = {kFALSE,kFALSE,kFALSE}; //Phi, v'(pt), v'(eta)
  11
  12
  13 void runFlowTaskExample(Int_t nRuns = 2,
  14               Bool_t DATA = kFALSE, const Char_t* dataDir="/Users/snelling/alice_data/Therminator_midcentral", Int_t offset = 0)
  15 {
  16   TStopwatch timer;
  17   timer.Start();
  18
  19   //--------------------------------------
  20   // Load the needed libraries most of them already loaded by aliroot
  21   //--------------------------------------
  22   gSystem->Load("libTree");
  23   gSystem->Load("libGeom");
  24   gSystem->Load("libVMC");
  25   gSystem->Load("libXMLIO");
  26   gSystem->Load("libPhysics");
  27   gSystem->Load("libSTEERBase");
  28   gSystem->Load("libESD");
  29   gSystem->Load("libAOD");
  30   gSystem->Load("libANALYSIS");
  31   gSystem->Load("libANALYSISalice");
  32   gSystem->Load("libCORRFW");
  33   gSystem->Load("libPWGflowBase");
  34   gSystem->Load("libPWGflowTasks");
  35
  36
  37   if (type!="AOD") { TChain* chain = CreateESDChain(dataDir, nRuns, offset);}
  38   else { TChain* chain = CreateAODChain(dataDir, nRuns, offset);}
  39
  40
  41   // CORRFW correction framework cuts
  42   //----------Event cuts----------
  43
  44   TObjArray* mcEventList = new TObjArray(0);
  45   AliCFEventGenCuts* mcEventCuts = new AliCFEventGenCuts("mcEventCuts","MC-level event cuts");
  46   mcEventList->AddLast(mcEventCuts);
  47
  48   TObjArray* recEventList = new TObjArray(0);
  49   AliCFEventRecCuts* recEventCuts = new AliCFEventRecCuts("recEventCuts","rec-level event cuts");
  50   recEventList->AddLast(recEventCuts);
  51
  52   //----------Cuts for RP----------
  53   TObjArray* mcListRP = new TObjArray(0);
  54   AliCFTrackKineCuts *mcKineCutsRP = new AliCFTrackKineCuts("mcKineCutsRP","MC-level kinematic cuts");
  55   mcListRP->AddLast(mcKineCutsRP);
  56   AliCFParticleGenCuts *mcGenCutsRP = new AliCFParticleGenCuts("mcGenCutsRP","MC particle generation cuts for RP");
  57   mcGenCutsRP->SetRequireIsPrimary();
  58   mcListRP->AddLast(mcGenCutsRP);
  59
  60   TObjArray* fPIDCutListRP = new TObjArray(0) ;
  61   AliCFTrackCutPid* cutPidRP = NULL;
  62   fPIDCutListRP->AddLast(cutPidRP);
  63
  64   TObjArray* recListRP = new TObjArray(0) ;
  65   AliCFTrackKineCuts *recKineCutsRP = new AliCFTrackKineCuts("recKineCutsRP","rec-level kine cuts");
  66   recListRP->AddLast(recKineCutsRP);
  67   AliCFTrackQualityCuts *recQualityCutsRP = new AliCFTrackQualityCuts("recQualityCutsRP","rec-level quality cuts");
  68   recQualityCutsRP->SetMinNClusterTPC(70);
  69   recListRP->AddLast(recQualityCutsRP);
  70   AliCFTrackIsPrimaryCuts *recIsPrimaryCutsRP = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsRP","rec-level isPrimary cuts");
  71   recIsPrimaryCutsRP->UseSPDvertex(kTRUE);
  72   recIsPrimaryCutsRP->UseTPCvertex(kTRUE);
  73   recListRP->AddLast(recIsPrimaryCutsRP);
  74
  75   TObjArray* accListRP = new TObjArray(0) ;
  76   AliCFAcceptanceCuts *mcAccCutsRP = new AliCFAcceptanceCuts("mcAccCutsRP","MC acceptance cuts");
  77   mcAccCutsRP->SetMinNHitITS(2);
  78   accListRP->AddLast(mcAccCutsRP);
  79
  80   //----------Cuts for POI----------
  81   TObjArray* mcListPOI = new TObjArray(0);
  82   AliCFTrackKineCuts *mcKineCutsPOI = new AliCFTrackKineCuts("mcKineCutsPOI","MC-level kinematic cuts");
  83   mcListPOI->AddLast(mcKineCutsPOI);
  84   AliCFParticleGenCuts *mcGenCutsPOI = new AliCFParticleGenCuts("mcGenCutsPOI","MC particle generation cuts for POI");
  85   mcGenCutsPOI->SetRequireIsPrimary();
  86   mcListPOI->AddLast(mcGenCutsPOI);
  87
  88   TObjArray* fPIDCutListPOI = new TObjArray(0) ;
  89   AliCFTrackCutPid* cutPidPOI = NULL;
  90   fPIDCutListPOI->AddLast(cutPidPOI);
  91
  92   TObjArray* recListPOI = new TObjArray(0) ;
  93   AliCFTrackKineCuts *recKineCutsPOI = new AliCFTrackKineCuts("recKineCutsPOI","rec-level kine cuts");
  94   recListPOI->AddLast(recKineCutsPOI);
  95   AliCFTrackQualityCuts *recQualityCutsPOI = new AliCFTrackQualityCuts("recQualityCutsPOI","rec-level quality cuts");
  96   recQualityCutsPOI->SetMinNClusterTPC(70);
  97   recListPOI->AddLast(recQualityCutsPOI);
  98   AliCFTrackIsPrimaryCuts *recIsPrimaryCutsPOI = new AliCFTrackIsPrimaryCuts("recIsPrimaryCutsPOI","rec-level isPrimary cuts");
  99   recIsPrimaryCutsPOI->UseSPDvertex(kTRUE);
 100   recIsPrimaryCutsPOI->UseTPCvertex(kTRUE);
 101   recListPOI->AddLast(recIsPrimaryCutsPOI);
 102
 103   TObjArray* accListPOI = new TObjArray(0) ;
 104   AliCFAcceptanceCuts *mcAccCutsPOI = new AliCFAcceptanceCuts("mcAccCutsPOI","MC acceptance cuts");
 105   mcAccCutsPOI->SetMinNHitITS(2);
 106   accListPOI->AddLast(mcAccCutsPOI);
 107
 108   //----------Add Cut Lists to the CF Manager----------
 109   printf("CREATE INTERFACE AND CUTS\n");
 110   AliCFManager* cfmgrRP = new AliCFManager();
 111   cfmgrRP->SetNStepEvent(3);
 112   cfmgrRP->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList);
 113   cfmgrRP->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList);
 114   cfmgrRP->SetNStepParticle(4);
 115   cfmgrRP->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListRP);
 116   cfmgrRP->SetParticleCutsList(AliCFManager::kPartAccCuts,accListRP);
 117   cfmgrRP->SetParticleCutsList(AliCFManager::kPartRecCuts,recListRP);
 118   cfmgrRP->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListRP);
 119
 120   AliCFManager* cfmgrPOI = new AliCFManager();
 121   cfmgrPOI->SetNStepEvent(3);
 122   cfmgrPOI->SetEventCutsList(AliCFManager::kEvtGenCuts,mcEventList);
 123   cfmgrPOI->SetEventCutsList(AliCFManager::kEvtRecCuts,recEventList);
 124   cfmgrPOI->SetNStepParticle(4);
 125   cfmgrPOI->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListPOI);
 126   cfmgrPOI->SetParticleCutsList(AliCFManager::kPartAccCuts,accListPOI);
 127   cfmgrPOI->SetParticleCutsList(AliCFManager::kPartRecCuts,recListPOI);
 128   cfmgrPOI->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutListPOI);
 129
 130
 131   // Make the analysis manager
 132   AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
 133
 134   AliVEventHandler* esdH = new AliESDInputHandler;
 135   mgr->SetInputEventHandler(esdH);
 136   AliMCEventHandler *mc = new AliMCEventHandler();
 137   mgr->SetMCtruthEventHandler(mc);
 138
 139
 140   //  gROOT->LoadMacro("$ALICE_ROOT/OADB/macros/AddTaskPhysicsSelection.C");
 141   //  AliPhysicsSelectionTask* physicsSelTask = AddTaskPhysicsSelection();
 142   //  if (!DATA) {physicsSelTask->GetPhysicsSelection()->SetAnalyzeMC();}
 143
 144   // Add the task which makes the flowevent to the analysis manager
 145   AliAnalysisTaskFlowEvent* taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",rptype,kFALSE,1);
 146   taskFE->SetAnalysisType(type);
 147   //  taskFE->SelectCollisionCandidates();
 148   mgr->AddTask(taskFE);
 149
 150   // Set the cuts for the flowevent
 151   taskFE->SetCFManager1(cfmgrRP);
 152   taskFE->SetCFManager2(cfmgrPOI);
 153
 154   // Instantiate the flow methods and add to the analysis manager
 155   AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane");
 156   mgr->AddTask(taskMCEP);
 157   AliAnalysisTaskQCumulants *taskQC = new AliAnalysisTaskQCumulants("TaskQCumulants",kFALSE);
 158   mgr->AddTask(taskQC);
 159
 160   // Create the output container for the data produced by the task
 161   // Connect to the input and output containers
 162   //===========================================================================
 163   AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
 164
 165   AliAnalysisDataContainer *coutputFE =
 166     mgr->CreateContainer("cobjFlowEventSimple",  AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
 167   mgr->ConnectInput(taskFE,0,cinput1);
 168   mgr->ConnectOutput(taskFE,1,coutputFE);
 169
 170
 171   TString outputMCEP = AliAnalysisManager::GetCommonFileName();
 172   outputMCEP += ":outputMCEPanalysis";
 173   outputMCEP+= type;
 174
 175   AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer("cobjMCEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP);
 176   mgr->ConnectInput(taskMCEP,0,coutputFE);
 177   mgr->ConnectOutput(taskMCEP,1,coutputMCEP);
 178
 179
 180   TString outputQC = AliAnalysisManager::GetCommonFileName();
 181   outputQC += ":outputQCanalysis";
 182   outputQC+= type;
 183
 184   AliAnalysisDataContainer *coutputQC = mgr->CreateContainer("cobjQC", TList::Class(),AliAnalysisManager::kOutputContainer,outputQC);
 185   mgr->ConnectInput(taskQC,0,coutputFE);
 186   mgr->ConnectOutput(taskQC,1,coutputQC);
 187
 188   // Run the analysis
 189   if (!mgr->InitAnalysis()) return;
 190   mgr->PrintStatus();
 191   mgr->StartAnalysis("local",chain);
 192
 193   timer.Stop();
 194   timer.Print();
 195
 196 }
 197
 198 // Helper macros for creating chains
 199 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 200
 201 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 202 {
 203   // creates chain of files in a given directory or file containing a list.
 204   // In case of directory the structure is expected as:
 205   // <aDataDir>/<dir0>/AliESDs.root
 206   // <aDataDir>/<dir1>/AliESDs.root
 207   // ...
 208
 209   if (!aDataDir)
 210     return 0;
 211
 212   Long_t id, size, flags, modtime;
 213   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 214     {
 215       printf("%s not found.\n", aDataDir);
 216       return 0;
 217     }
 218
 219   TChain* chain = new TChain("esdTree");
 220   TChain* chaingAlice = 0;
 221
 222   if (flags & 2)
 223     {
 224       TString execDir(gSystem->pwd());
 225       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 226       TList* dirList            = baseDir->GetListOfFiles();
 227       Int_t nDirs               = dirList->GetEntries();
 228       gSystem->cd(execDir);
 229
 230       Int_t count = 0;
 231
 232       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 233         {
 234           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 235           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 236             continue;
 237
 238           if (offset > 0)
 239             {
 240               --offset;
 241               continue;
 242             }
 243
 244           if (count++ == aRuns)
 245             break;
 246
 247           TString presentDirName(aDataDir);
 248           presentDirName += "/";
 249           presentDirName += presentDir->GetName();
 250           chain->Add(presentDirName + "/AliESDs.root/esdTree");
 251           //  cerr<<presentDirName<<endl;
 252         }
 253
 254     }
 255   else
 256     {
 257       // Open the input stream
 258       ifstream in;
 259       in.open(aDataDir);
 260
 261       Int_t count = 0;
 262
 263       // Read the input list of files and add them to the chain
 264       TString esdfile;
 265       while(in.good()) {
 266         in >> esdfile;
 267         if (!esdfile.Contains("root")) continue; // protection
 268
 269         if (offset > 0)
 270           {
 271             --offset;
 272             continue;
 273           }
 274
 275         if (count++ == aRuns)
 276           break;
 277
 278         // add esd file
 279         chain->Add(esdfile);
 280       }
 281
 282       in.close();
 283     }
 284
 285   return chain;
 286 }
 287
 288
 289 // Helper macros for creating chains
 290 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 291
 292 TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 293 {
 294   // creates chain of files in a given directory or file containing a list.
 295   // In case of directory the structure is expected as:
 296   // <aDataDir>/<dir0>/AliAOD.root
 297   // <aDataDir>/<dir1>/AliAOD.root
 298   // ...
 299
 300   if (!aDataDir)
 301     return 0;
 302
 303   Long_t id, size, flags, modtime;
 304   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 305     {
 306       printf("%s not found.\n", aDataDir);
 307       return 0;
 308     }
 309
 310   TChain* chain = new TChain("aodTree");
 311   TChain* chaingAlice = 0;
 312
 313   if (flags & 2)
 314     {
 315       TString execDir(gSystem->pwd());
 316       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 317       TList* dirList            = baseDir->GetListOfFiles();
 318       Int_t nDirs               = dirList->GetEntries();
 319       gSystem->cd(execDir);
 320
 321       Int_t count = 0;
 322
 323       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 324         {
 325           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 326           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 327             continue;
 328
 329           if (offset > 0)
 330             {
 331               --offset;
 332               continue;
 333             }
 334
 335           if (count++ == aRuns)
 336             break;
 337
 338           TString presentDirName(aDataDir);
 339           presentDirName += "/";
 340           presentDirName += presentDir->GetName();
 341           chain->Add(presentDirName + "/AliAOD.root/aodTree");
 342           // cerr<<presentDirName<<endl;
 343         }
 344
 345     }
 346   else
 347     {
 348       // Open the input stream
 349       ifstream in;
 350       in.open(aDataDir);
 351
 352       Int_t count = 0;
 353
 354       // Read the input list of files and add them to the chain
 355       TString aodfile;
 356       while(in.good()) {
 357         in >> aodfile;
 358         if (!aodfile.Contains("root")) continue; // protection
 359
 360         if (offset > 0)
 361           {
 362             --offset;
 363             continue;
 364           }
 365
 366         if (count++ == aRuns)
 367           break;
 368
 369         // add aod file
 370         chain->Add(aodfile);
 371       }
 372
 373       in.close();
 374     }
 375
 376   return chain;
 377 }
 378