PWG2/FLOW/macros/runFlowTask.C

   1 enum anaModes {mLocal,mLocalPAR,mPROOF,mGRID};
   2 //mLocal: Analyze locally files in your computer using aliroot
   3 //mLocalPAR: Analyze locally files in your computer using root + PAR files
   4 //mPROOF: Analyze CAF files with PROOF
   5
   6 // RUN SETTINGS
   7
   8 // Flow analysis method can be:(set to kTRUE or kFALSE)
   9 Bool_t MCEP     = kTRUE;  // correlation with Monte Carlo reaction plane
  10 Bool_t SP       = kTRUE;  // scalar product method (similar to eventplane method)
  11 Bool_t GFC      = kTRUE;  // cumulants based on generating function
  12 Bool_t QC       = kTRUE;  // cumulants using Q vectors
  13 Bool_t FQD      = kTRUE;  // fit of the distribution of the Q vector (only integrated v)
  14 Bool_t LYZ1SUM  = kTRUE;  // Lee Yang Zeroes using sum generating function (integrated v)
  15 Bool_t LYZ1PROD = kTRUE;  // Lee Yang Zeroes using product generating function (integrated v)
  16 Bool_t LYZ2SUM  = kFALSE; // Lee Yang Zeroes using sum generating function (second pass differential v)
  17 Bool_t LYZ2PROD = kFALSE; // Lee Yang Zeroes using product generating function (second pass differential v)
  18 Bool_t LYZEP    = kFALSE; // Lee Yang Zeroes Event plane using sum generating function (gives eventplane + weight)
  19 Bool_t MH       = kTRUE;  // azimuthal correlators in mixed harmonics
  20 Bool_t NL       = kTRUE;  // nested loops (for instance distribution of phi1-phi2 for all distinct pairs)
  21
  22 Bool_t METHODS[] = {SP,LYZ1SUM,LYZ1PROD,LYZ2SUM,LYZ2PROD,LYZEP,GFC,QC,FQD,MCEP,MH,NL};
  23
  24 // Analysis type can be ESD, AOD, MC, ESDMCkineESD, ESDMCkineMC
  25 const TString type = "ESD";
  26
  27 // Boolean to fill/not fill the QA histograms
  28 Bool_t QA = kTRUE;
  29
  30 // Boolean to use/not use weights for the Q vector
  31 Bool_t WEIGHTS[] = {kFALSE,kFALSE,kFALSE}; //Phi, v'(pt), v'(eta)
  32
  33
  34 //void runFlowTask(Int_t mode=mLocal, Int_t nRuns = 2,
  35 //Bool_t DATA = kTRUE, const Char_t* dataDir="/data/alice2/kolk/PP/data/LHC09d/104892/test", Int_t offset = 0)
  36 //              Bool_t DATA = kFALSE, const Char_t* dataDir="/data/alice2/kolk/PP/LHC09d10/104873", Int_t offset = 0)
  37
  38 void runFlowTask(Int_t mode = mPROOF, Int_t nRuns = 50000000,
  39                  //Bool_t DATA = kFALSE, const Char_t* dataDir="/PWG2/akisiel/Therminator_midcentral_ESD", Int_t offset=0)
  40                  //Bool_t DATA = kFALSE, const Char_t* dataDir="/PWG2/akisiel/LHC10d6_0.9TeV_EPOS_12502X", Int_t offset=0)
  41                  //Bool_t DATA = kFALSE, const Char_t* dataDir="/alice/sim/LHC10d2_117048", Int_t offset=0) //phojet 7 TeV
  42                  //Bool_t DATA = kTRUE, const Char_t* dataDir="/alice/data/LHC09d_000104792_p6", Int_t offset=0) //data 0.9 TeV
  43                  Bool_t DATA = kFALSE, const Char_t* dataDir="/PWG4/morsch/HIJING_CENT_4EV", Int_t offset=0) //hijing Pb Pb pilot
  44
  45 //void runFlowTask(Int_t mode = mGRID, Bool_t DATA = kTRUE)
  46 {
  47   TStopwatch timer;
  48   timer.Start();
  49
  50   CrossCheckUserSettings(DATA);
  51
  52   LoadLibraries(mode);
  53
  54   if (mode == mGRID) {
  55     // Create and configure the alien handler plugin
  56     gROOT->LoadMacro("CreateAlienHandler.C");
  57     AliAnalysisGrid *alienHandler = CreateAlienHandler();
  58     if (!alienHandler) return;
  59   }
  60
  61   if (mode==mLocal || mode == mLocalPAR) {
  62     if (type!="AOD") { TChain* chain = CreateESDChain(dataDir, nRuns, offset);}
  63     else { TChain* chain = CreateAODChain(dataDir, nRuns, offset);}
  64   }
  65   //____________________________________________//
  66   // Make the analysis manager
  67   AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
  68
  69   if (mode == mGRID) {
  70     // Connect plug-in to the analysis manager
  71     mgr->SetGridHandler(alienHandler);
  72   }
  73
  74   if (type == "ESD"){
  75     AliVEventHandler* esdH = new AliESDInputHandler;
  76     mgr->SetInputEventHandler(esdH);
  77     if (MCEP) {
  78       AliMCEventHandler *mc = new AliMCEventHandler();
  79       mgr->SetMCtruthEventHandler(mc);
  80     }
  81   }
  82
  83   if (type == "AOD"){
  84     AliVEventHandler* aodH = new AliAODInputHandler;
  85     mgr->SetInputEventHandler(aodH);
  86     if (MCEP) {
  87       AliMCEventHandler *mc = new AliMCEventHandler();
  88       mgr->SetMCtruthEventHandler(mc);
  89     }
  90   }
  91
  92   if (type == "MC" || type == "ESDMCkineESD" || type == "ESDMCkineMC"){
  93     AliVEventHandler* esdH = new AliESDInputHandler;
  94     mgr->SetInputEventHandler(esdH);
  95
  96     AliMCEventHandler *mc = new AliMCEventHandler();
  97     mgr->SetMCtruthEventHandler(mc);
  98   }
  99
 100   //____________________________________________//
 101   // Load the analysis task
 102   gROOT->LoadMacro("AddTaskFlow.C");
 103   AliAnalysisTaskFlowEvent* taskFE = AddTaskFlow(type,METHODS,QA,WEIGHTS);
 104
 105
 106   // Task to check the offline trigger
 107   if (mode == mLocal || mode == mGRID) {
 108     gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); }
 109   else if (mode == mPROOF || mode == mLocalPAR) {
 110     gROOT->LoadMacro("AddTaskPhysicsSelection.C"); }
 111   AliPhysicsSelectionTask* physicsSelTask = AddTaskPhysicsSelection();
 112   if (!DATA) {physicsSelTask->GetPhysicsSelection()->SetAnalyzeMC();}
 113
 114
 115   // Enable debug printouts
 116   mgr->SetDebugLevel(2);
 117
 118
 119   //____________________________________________//
 120   // Run the analysis
 121   if (!mgr->InitAnalysis()) return;
 122   mgr->PrintStatus();
 123
 124   if (mode == mLocal || mode == mLocalPAR) {
 125     mgr->StartAnalysis("local",chain);
 126   }
 127   else if (mode == mPROOF) {
 128     mgr->StartAnalysis("proof",dataDir,nRuns,offset);
 129   }
 130   else if (mode == mGRID) {
 131     mgr->StartAnalysis("grid");
 132   }
 133
 134   timer.Stop();
 135   timer.Print();
 136
 137 }
 138
 139 void CrossCheckUserSettings(Bool_t bData)
 140 {
 141  // Check in this method if the user settings make sense.
 142
 143  if(MCEP==kTRUE && bData==kTRUE)
 144  {
 145   cout<<endl;
 146   cout<<"WARNING: In real datasets there is no Monte Carlo information available !!!!"<<endl;
 147   cout<<"         Set for real data analysis DATA = kTRUE and MCEP = kFALSE in the macro."<<endl;
 148   cout<<endl;
 149   exit(0);
 150  }
 151
 152 } // end of void CrossCheckUserSettings()
 153
 154 void LoadLibraries(const anaModes mode) {
 155
 156   //--------------------------------------
 157   // Load the needed libraries most of them already loaded by aliroot
 158   //--------------------------------------
 159   gSystem->Load("libTree");
 160   gSystem->Load("libGeom");
 161   gSystem->Load("libVMC");
 162   gSystem->Load("libXMLIO");
 163   gSystem->Load("libPhysics");
 164
 165   //----------------------------------------------------------
 166   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 167   //----------------------------------------------------------
 168   if (mode==mLocal || mode==mGRID) {
 169     //--------------------------------------------------------
 170     // If you want to use already compiled libraries
 171     // in the aliroot distribution
 172     //--------------------------------------------------------
 173     gSystem->Load("libSTEERBase");
 174     gSystem->Load("libESD");
 175     gSystem->Load("libAOD");
 176     gSystem->Load("libANALYSIS");
 177     gSystem->Load("libANALYSISalice");
 178     gSystem->Load("libCORRFW");
 179     gSystem->Load("libPWG2forward");
 180     if (mode==mLocal) {
 181       gSystem->Load("libPWG2flowCommon");
 182       cerr<<"libPWG2flowCommon loaded..."<<endl;
 183       gSystem->Load("libPWG2flowTasks");
 184       cerr<<"libPWG2flowTasks loaded..."<<endl;
 185     }
 186     if (mode==mGRID) {
 187       SetupPar("PWG2flowCommon");
 188       cerr<<"PWG2flowCommon.par loaded..."<<endl;
 189       SetupPar("PWG2flowTasks");
 190       cerr<<"PWG2flowTasks.par loaded..."<<endl;
 191     }
 192   }
 193
 194   else if (mode == mLocalPAR) {
 195     //--------------------------------------------------------
 196     //If you want to use root and par files from aliroot
 197     //--------------------------------------------------------
 198     SetupPar("STEERBase");
 199     SetupPar("ESD");
 200     SetupPar("AOD");
 201     SetupPar("ANALYSIS");
 202     SetupPar("ANALYSISalice");
 203         SetupPar("CORRFW");
 204     SetupPar("PWG2flowCommon");
 205     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 206     SetupPar("PWG2flowTasks");
 207     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 208   }
 209
 210   //---------------------------------------------------------
 211   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 212   //---------------------------------------------------------
 213   else if (mode==mPROOF) {
 214     //
 215     //gEnv->SetValue("XSec.GSI.DelegProxy","2");
 216     //  set to debug root versus if needed
 217     //TProof::Mgr("alicecaf")->SetROOTVersion("v5-24-00a_dbg");
 218     //TProof::Mgr("alicecaf")->SetROOTVersion("v5-24-00a");
 219     //TProof::Reset("proof://snelling@alicecaf.cern.ch");
 220     // Connect to proof
 221     printf("*** Connect to PROOF ***\n");
 222     gEnv->SetValue("XSec.GSI.DelegProxy","2");
 223     // Put appropriate username here
 224     //TProof::Open("abilandz@alicecaf.cern.ch");
 225     //TProof::Open("nkolk@alicecaf.cern.ch");
 226     //TProof::Open("snelling@localhost");
 227     TProof::Open("alice-caf.cern.ch");
 228     //TProof::Open("skaf.saske.sk");
 229     //TProof::Open("prf000-iep-grid.saske.sk");
 230     //Info("runSKAF.C","Loading libs on proof (may take while, around 1 min) ...");
 231     // list the data available
 232     //gProof->ShowDataSets("/*/*");
 233     //gProof->ShowDataSets("/alice/sim/"); //for MC Data
 234     //gProof->ShowDataSets("/alice/data/"); //for REAL Data
 235
 236     // Clear the Packages
 237
 238     gProof->ClearPackage("STEERBase.par");
 239     gProof->ClearPackage("ESD.par");
 240     gProof->ClearPackage("AOD.par");
 241     gProof->ClearPackage("ANALYSIS.par");
 242     gProof->ClearPackage("ANALYSISalice.par");
 243     gProof->ClearPackage("CORRFW.par");
 244
 245     gProof->ClearPackage("PWG2flowCommon");
 246     gProof->ClearPackage("PWG2flowTasks");
 247
 248     // Upload the Packages
 249     gProof->UploadPackage("STEERBase.par");
 250     gProof->UploadPackage("ESD.par");
 251     gProof->UploadPackage("AOD.par");
 252
 253     gProof->UploadPackage("ANALYSIS.par");
 254     gProof->UploadPackage("ANALYSISalice.par");
 255     gProof->UploadPackage("CORRFW.par");
 256     gProof->UploadPackage("PWG2flowCommon.par");
 257     gProof->UploadPackage("PWG2flowTasks.par");
 258
 259     // Enable the Packages
 260     // The global package
 261     //gProof->EnablePackage("aliroot_v4-19-05-AN",kTRUE);
 262     // Or separate
 263
 264     gProof->EnablePackage("STEERBase");
 265     gProof->EnablePackage("ESD");
 266     gProof->EnablePackage("AOD");
 267
 268     // Always needed
 269     gProof->EnablePackage("ANALYSIS");
 270     gProof->EnablePackage("ANALYSISalice");
 271     gProof->EnablePackage("CORRFW");
 272     gProof->EnablePackage("PWG2flowCommon");
 273     gProof->EnablePackage("PWG2flowTasks");
 274
 275     // Show enables Packages
 276     gProof->ShowEnabledPackages();
 277   }
 278
 279 }
 280
 281 void SetupPar(char* pararchivename) {
 282   //Load par files, create analysis libraries
 283   //For testing, if par file already decompressed and modified
 284   //classes then do not decompress.
 285
 286   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 287   TString parpar(Form("%s.par", pararchivename)) ;
 288   if ( gSystem->AccessPathName(parpar.Data()) ) {
 289     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 290     TString processline(Form(".! make %s", parpar.Data())) ;
 291     gROOT->ProcessLine(processline.Data()) ;
 292     gSystem->ChangeDirectory(cdir) ;
 293     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 294     gROOT->ProcessLine(processline.Data()) ;
 295   }
 296   if ( gSystem->AccessPathName(pararchivename) ) {
 297     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 298     gROOT->ProcessLine(processline.Data());
 299   }
 300
 301   TString ocwd = gSystem->WorkingDirectory();
 302   gSystem->ChangeDirectory(pararchivename);
 303
 304   // check for BUILD.sh and execute
 305   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 306     printf("*******************************\n");
 307     printf("*** Building PAR archive    ***\n");
 308     cout<<pararchivename<<endl;
 309     printf("*******************************\n");
 310     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 311       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 312       return -1;
 313     }
 314   }
 315   // check for SETUP.C and execute
 316   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 317     printf("*******************************\n");
 318     printf("*** Setup PAR archive       ***\n");
 319     cout<<pararchivename<<endl;
 320     printf("*******************************\n");
 321     gROOT->Macro("PROOF-INF/SETUP.C");
 322   }
 323
 324   gSystem->ChangeDirectory(ocwd.Data());
 325   printf("Current dir: %s\n", ocwd.Data());
 326 }
 327
 328
 329 // Helper macros for creating chains
 330 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 331
 332 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 333 {
 334   // creates chain of files in a given directory or file containing a list.
 335   // In case of directory the structure is expected as:
 336   // <aDataDir>/<dir0>/AliESDs.root
 337   // <aDataDir>/<dir1>/AliESDs.root
 338   // ...
 339
 340   if (!aDataDir)
 341     return 0;
 342
 343   Long_t id, size, flags, modtime;
 344   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 345     {
 346       printf("%s not found.\n", aDataDir);
 347       return 0;
 348     }
 349
 350   TChain* chain = new TChain("esdTree");
 351   TChain* chaingAlice = 0;
 352
 353   if (flags & 2)
 354     {
 355       TString execDir(gSystem->pwd());
 356       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 357       TList* dirList            = baseDir->GetListOfFiles();
 358       Int_t nDirs               = dirList->GetEntries();
 359       gSystem->cd(execDir);
 360
 361       Int_t count = 0;
 362
 363       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 364         {
 365           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 366           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 367             continue;
 368
 369           if (offset > 0)
 370             {
 371               --offset;
 372               continue;
 373             }
 374
 375           if (count++ == aRuns)
 376             break;
 377
 378           TString presentDirName(aDataDir);
 379           presentDirName += "/";
 380           presentDirName += presentDir->GetName();
 381           chain->Add(presentDirName + "/AliESDs.root/esdTree");
 382           //  cerr<<presentDirName<<endl;
 383         }
 384
 385     }
 386   else
 387     {
 388       // Open the input stream
 389       ifstream in;
 390       in.open(aDataDir);
 391
 392       Int_t count = 0;
 393
 394       // Read the input list of files and add them to the chain
 395       TString esdfile;
 396       while(in.good()) {
 397         in >> esdfile;
 398         if (!esdfile.Contains("root")) continue; // protection
 399
 400         if (offset > 0)
 401           {
 402             --offset;
 403             continue;
 404           }
 405
 406         if (count++ == aRuns)
 407           break;
 408
 409         // add esd file
 410         chain->Add(esdfile);
 411       }
 412
 413       in.close();
 414     }
 415
 416   return chain;
 417 }
 418
 419
 420 // Helper macros for creating chains
 421 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 422
 423 TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 424 {
 425   // creates chain of files in a given directory or file containing a list.
 426   // In case of directory the structure is expected as:
 427   // <aDataDir>/<dir0>/AliAOD.root
 428   // <aDataDir>/<dir1>/AliAOD.root
 429   // ...
 430
 431   if (!aDataDir)
 432     return 0;
 433
 434   Long_t id, size, flags, modtime;
 435   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 436     {
 437       printf("%s not found.\n", aDataDir);
 438       return 0;
 439     }
 440
 441   TChain* chain = new TChain("aodTree");
 442   TChain* chaingAlice = 0;
 443
 444   if (flags & 2)
 445     {
 446       TString execDir(gSystem->pwd());
 447       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 448       TList* dirList            = baseDir->GetListOfFiles();
 449       Int_t nDirs               = dirList->GetEntries();
 450       gSystem->cd(execDir);
 451
 452       Int_t count = 0;
 453
 454       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 455         {
 456           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 457           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 458             continue;
 459
 460           if (offset > 0)
 461             {
 462               --offset;
 463               continue;
 464             }
 465
 466           if (count++ == aRuns)
 467             break;
 468
 469           TString presentDirName(aDataDir);
 470           presentDirName += "/";
 471           presentDirName += presentDir->GetName();
 472           chain->Add(presentDirName + "/AliAOD.root/aodTree");
 473           // cerr<<presentDirName<<endl;
 474         }
 475
 476     }
 477   else
 478     {
 479       // Open the input stream
 480       ifstream in;
 481       in.open(aDataDir);
 482
 483       Int_t count = 0;
 484
 485       // Read the input list of files and add them to the chain
 486       TString aodfile;
 487       while(in.good()) {
 488         in >> aodfile;
 489         if (!aodfile.Contains("root")) continue; // protection
 490
 491         if (offset > 0)
 492           {
 493             --offset;
 494             continue;
 495           }
 496
 497         if (count++ == aRuns)
 498           break;
 499
 500         // add aod file
 501         chain->Add(aodfile);
 502       }
 503
 504       in.close();
 505     }
 506
 507   return chain;
 508 }
 509