PWG2/FLOW/macros/runAliAnalysisTaskFlow.C

   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) ;
   4
   5 //SETTING THE ANALYSIS
   6
   7 //Flow analysis method can be:
   8 // SP    = Scalar Product
   9 // LYZ1  = Lee Yang Zeroes first run
  10 // LYZ2  = Lee Yang Zeroes second run
  11 // LYZEP = Lee Yang Zeroes Event Plane
  12 // CUM   = Cumulants
  13 const TString method = "SP";
  14
  15 //Type of analysis can be:
  16 // ESD, AOD, MC, ESDMC0, ESDMC1
  17 const TString type = "ESDMC1";
  18
  19
  20 //SETTING THE CUTS
  21
  22 //for integrated flow
  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;
  33
  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;
  45
  46
  47
  48 void runAliAnalysisTaskFlow(Int_t nRuns = 10, const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX", Int_t offset = 0)
  49
  50 {
  51   TStopwatch timer;
  52   timer.Start();
  53
  54   // include path (to find the .h files when compiling)
  55   gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
  56   gSystem->AddIncludePath("-I$ROOTSYS/include") ;
  57
  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;
  70
  71   // create the TChain. CreateESDChain() is defined in CreateESDChain.C
  72   TChain* chain = CreateESDChain(dataDir, nRuns, offset);
  73   cout<<"chain ("<<chain<<")"<<endl;
  74
  75   //____________________________________________//
  76   //Create cuts using correction framework
  77
  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);
  83
  84   AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
  85   mcKineCuts2->SetPtRange(ptmin2,ptmax2);
  86   mcKineCuts2->SetRapidityRange(ymin2,ymax2);
  87   mcKineCuts2->SetChargeMC(charge2);
  88
  89   AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
  90   mcGenCuts1->SetRequireIsPrimary();
  91   mcGenCuts1->SetRequirePdgCode(PDG1);
  92
  93   AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
  94   mcGenCuts2->SetRequireIsPrimary();
  95   mcGenCuts2->SetRequirePdgCode(PDG2);
  96
  97   //############# Acceptance Cuts  ????????
  98   AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts","MC acceptance cuts");
  99   mcAccCuts->SetMinNHitITS(mintrackrefsITS1);
 100   mcAccCuts->SetMinNHitTPC(mintrackrefsTPC1);
 101
 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);
 107
 108   AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
 109   recKineCuts2->SetPtRange(ptmin2,ptmax2);
 110   recKineCuts2->SetRapidityRange(ymin2,ymax2);
 111   recKineCuts2->SetChargeRec(charge2);
 112
 113   AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
 114   recQualityCuts->SetMinNClusterTPC(minclustersTPC1);
 115   recQualityCuts->SetRequireITSRefit(kTRUE);
 116
 117   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
 118   recIsPrimaryCuts->SetMaxNSigmaToVertex(maxnsigmatovertex1);
 119
 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];
 124
 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 ;
 130
 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;
 141   }
 142
 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;
 153   }
 154
 155
 156   printf("CREATE MC KINE CUTS\n");
 157   TObjArray* mcList1 = new TObjArray(0);
 158   mcList1->AddLast(mcKineCuts1);
 159   mcList1->AddLast(mcGenCuts1);
 160
 161   TObjArray* mcList2 = new TObjArray(0);
 162   mcList2->AddLast(mcKineCuts2);
 163   mcList2->AddLast(mcGenCuts2);
 164
 165   printf("CREATE ACCEPTANCE CUTS\n");
 166   TObjArray* accList = new TObjArray(0) ;
 167   accList->AddLast(mcAccCuts);
 168
 169   printf("CREATE RECONSTRUCTION CUTS\n");
 170   TObjArray* recList1 = new TObjArray(0) ;
 171   recList1->AddLast(recKineCuts1);
 172   recList1->AddLast(recQualityCuts);
 173   recList1->AddLast(recIsPrimaryCuts);
 174
 175   TObjArray* recList2 = new TObjArray(0) ;
 176   recList2->AddLast(recKineCuts2);
 177   recList2->AddLast(recQualityCuts);
 178   recList2->AddLast(recIsPrimaryCuts);
 179
 180   printf("CREATE PID CUTS\n");
 181   TObjArray* fPIDCutList1 = new TObjArray(0) ;
 182   fPIDCutList1->AddLast(cutPID1);
 183
 184   TObjArray* fPIDCutList2 = new TObjArray(0) ;
 185   fPIDCutList2->AddLast(cutPID2);
 186
 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);
 193
 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);
 199
 200
 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;
 210   }
 211
 212
 213   //____________________________________________//
 214   // Make the analysis manager
 215   AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
 216
 217   if (type == "ESD"){
 218     AliVEventHandler* esdH = new AliESDInputHandler;
 219     mgr->SetInputEventHandler(esdH); }
 220
 221   if (type == "AOD"){
 222     AliVEventHandler* aodH = new AliAODInputHandler;
 223     mgr->SetInputEventHandler(aodH); }
 224
 225   if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
 226     AliVEventHandler* esdH = new AliESDInputHandler;
 227     mgr->SetInputEventHandler(esdH);
 228
 229     AliMCEventHandler *mc = new AliMCEventHandler();
 230     mgr->SetMCtruthEventHandler(mc); }
 231
 232   //____________________________________________//
 233   // 1st MC EP task
 234   if (method == "SP"){
 235     AliAnalysisTaskScalarProduct *task1 = new AliAnalysisTaskScalarProduct("TaskScalarProduct");
 236     task1->SetAnalysisType(type);
 237     task1->SetCFManager1(cfmgr1);
 238     task1->SetCFManager2(cfmgr2);
 239     mgr->AddTask(task1);
 240   }
 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);
 248     mgr->AddTask(task1);
 249   }
 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);
 257     mgr->AddTask(task1);
 258   }
 259   else if (method == "LYZEP"){
 260     AliAnalysisTaskLYZEventPlane *task1 = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane");
 261     task1->SetAnalysisType(type);
 262     //task1->SetCFManager1(cfmgr1);
 263     //task1->SetCFManager2(cfmgr2);
 264     mgr->AddTask(task1);
 265   }
 266   else if (method == "CUM"){
 267     AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
 268     task1->SetAnalysisType(type);
 269     //task1->SetCFManager1(cfmgr1);
 270     //task1->SetCFManager2(cfmgr2);
 271     mgr->AddTask(task1);
 272   }
 273
 274
 275   // Create containers for input/output
 276   AliAnalysisDataContainer *cinput1 =
 277     mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
 278
 279   if (method == "LYZ2"){ AliAnalysisDataContainer *cinput2 =
 280                     mgr->CreateContainer("cobj2",TList::Class(),AliAnalysisManager::kInputContainer); }
 281
 282
 283   TString outputName = "output";
 284   outputName+= method;
 285   outputName+= "analysis";
 286   outputName+= type;
 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);
 292
 293   //____________________________________________//
 294   mgr->ConnectInput(task1,0,cinput1);
 295   if (method == "LYZ2") { mgr->ConnectInput(task1,1,cinput2); }
 296   mgr->ConnectOutput(task1,0,coutput1);
 297
 298   if (method == "LYZ2"){ cinput2->SetData(fFirstRunFile);}
 299
 300   if (!mgr->InitAnalysis()) return;
 301   mgr->PrintStatus();
 302   mgr->StartAnalysis("local",chain);
 303
 304   timer.Stop();
 305   timer.Print();
 306 }
 307
 308
 309 // Helper macros for creating chains
 310 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 311
 312 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 313 {
 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
 318   // ...
 319
 320   if (!aDataDir)
 321     return 0;
 322
 323   Long_t id, size, flags, modtime;
 324   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 325     {
 326       printf("%s not found.\n", aDataDir);
 327       return 0;
 328     }
 329
 330   TChain* chain = new TChain("esdTree");
 331   TChain* chaingAlice = 0;
 332
 333   if (flags & 2)
 334     {
 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);
 340
 341       Int_t count = 0;
 342
 343       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 344         {
 345           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 346           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 347             continue;
 348
 349           if (offset > 0)
 350             {
 351               --offset;
 352               continue;
 353             }
 354
 355           if (count++ == aRuns)
 356             break;
 357
 358           TString presentDirName(aDataDir);
 359           presentDirName += "/";
 360           presentDirName += presentDir->GetName();
 361           chain->Add(presentDirName + "/AliESDs.root/esdTree");
 362           cerr<<presentDirName<<endl;
 363         }
 364
 365     }
 366   else
 367     {
 368       // Open the input stream
 369       ifstream in;
 370       in.open(aDataDir);
 371
 372       Int_t count = 0;
 373
 374       // Read the input list of files and add them to the chain
 375       TString esdfile;
 376       while(in.good()) {
 377         in >> esdfile;
 378         if (!esdfile.Contains("root")) continue; // protection
 379
 380         if (offset > 0)
 381           {
 382             --offset;
 383             continue;
 384           }
 385
 386         if (count++ == aRuns)
 387           break;
 388
 389         // add esd file
 390         chain->Add(esdfile);
 391       }
 392
 393       in.close();
 394     }
 395
 396   return chain;
 397 }
 398
 399 void LookupWrite(TChain* chain, const char* target)
 400 {
 401   // looks up the chain and writes the remaining files to the text file target
 402
 403   chain->Lookup();
 404
 405   TObjArray* list = chain->GetListOfFiles();
 406   TIterator* iter = list->MakeIterator();
 407   TObject* obj = 0;
 408
 409   ofstream outfile;
 410   outfile.open(target);
 411
 412   while ((obj = iter->Next()))
 413     outfile << obj->GetTitle() << "#AliESDs.root" << endl;
 414
 415   outfile.close();
 416
 417   delete iter;
 418 }