PWG2/FLOW/macros/runFlowAnalysis.C

   1 #include "TStopwatch.h"
   2 #include "TObjArray"
   3 #include "Riostream.h"
   4 #include "TFile.h"
   5
   6 //--------------------------------------------------------------------------------------
   7 // RUN SETTINGS
   8 //flow analysis method can be: (set to kTRUE or kFALSE)
   9 Bool_t SP    = kTRUE;
  10 Bool_t LYZ1  = kTRUE;
  11 Bool_t LYZ2  = kFALSE;
  12 Bool_t LYZEP = kFALSE;
  13 Bool_t GFC   = kTRUE;
  14 Bool_t QC    = kTRUE;
  15 Bool_t FQD   = kTRUE;
  16 Bool_t MCEP  = kFALSE; //does not work yet 24/12/08
  17 //--------------------------------------------------------------------------------------
  18
  19 // Weights
  20 // Use weights for Q vector
  21 Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
  22 Bool_t usePtWeights  = kFALSE; //v'(pt) (differential flow in pt)
  23 Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
  24
  25 //--------------------------------------------------------------------------------------
  26 // CUT SETTINGS
  27 //integrated selection
  28 Double_t ptMaxInt  = 10.;
  29 Double_t ptMinInt  = 0.;
  30 Double_t etaMaxInt = 1.;
  31 Double_t etaMinInt = -1.;
  32 Double_t phiMaxInt = 7.5;
  33 Double_t phiMinInt = 0.;
  34 Int_t PIDInt       = 211;
  35
  36 //differential selection
  37 Double_t ptMaxDiff  = 10.;
  38 Double_t ptMinDiff  = 0.;
  39 Double_t etaMaxDiff = 1.;
  40 Double_t etaMinDiff = -1.;
  41 Double_t phiMaxDiff = 7.5;
  42 Double_t phiMinDiff = 0.;
  43 Int_t PIDDiff       = 211;
  44 //--------------------------------------------------------------------------------------
  45
  46 enum anaModes {mLocal,mLocalSource,mLocalPAR,};
  47 //mLocal: Analyze data on your computer using aliroot
  48 //mLocalPAR: Analyze data on your computer using root + PAR files
  49 //mLocalSource: Analyze data on your computer using root + source files
  50
  51
  52
  53 Int_t offset = 0;
  54
  55 int runFlowAnalysis(Int_t mode=mLocal, Int_t aRuns = 100, const char*
  56                     //                    dir="/data/alice1/kolk/KineOnly3/")
  57                     dir="/Users/snelling/alice_data/KineOnly3/")
  58 {
  59   TStopwatch timer;
  60   timer.Start();
  61
  62   if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
  63   if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  64   if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  65
  66
  67   cout<<endl;
  68   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
  69   cout<<endl;
  70
  71   LoadLibraries(mode);
  72
  73   if (mode == mLocal || mode == mLocalPAR) {
  74     // AliFlow event in aliroot or with pars
  75     AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
  76   }
  77   else if (mode == mLocalSource) {
  78     // flow event in source mode
  79     FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
  80   }
  81   else{
  82     cout << "No supported running mode selected!" << endl;
  83     break;
  84   }
  85
  86
  87   //------------------------------------------------------------------------
  88   //cuts:
  89   AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
  90   cutsInt->SetPtMax(ptMaxInt);
  91   cutsInt->SetPtMin(ptMinInt);
  92   cutsInt->SetEtaMax(etaMaxInt);
  93   cutsInt->SetEtaMin(etaMinInt);
  94   cutsInt->SetPhiMax(phiMaxInt);
  95   cutsInt->SetPhiMin(phiMinInt);
  96   cutsInt->SetPID(PIDInt);
  97
  98   AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
  99   cutsDiff->SetPtMax(ptMaxDiff);
 100   cutsDiff->SetPtMin(ptMinDiff);
 101   cutsDiff->SetEtaMax(etaMaxDiff);
 102   cutsDiff->SetEtaMin(etaMinDiff);
 103   cutsDiff->SetPhiMax(phiMaxDiff);
 104   cutsDiff->SetPhiMin(phiMinDiff);
 105   cutsDiff->SetPID(PIDDiff);
 106
 107   //if the weights are used:
 108   TFile *fileWithWeights = NULL;
 109   TList *listWithWeights = NULL;
 110
 111   if(usePhiWeights||usePtWeights||useEtaWeights) {
 112     fileWithWeights = TFile::Open("weights.root","READ");
 113     if(fileWithWeights) {
 114       listWithWeights = (TList*)fileWithWeights->Get("weights");
 115     }
 116     else
 117       {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 118         break;
 119       }
 120   }
 121
 122   //flow methods:
 123   AliFlowAnalysisWithQCumulants    *qc    = NULL;
 124   AliFlowAnalysisWithCumulants     *gfc   = NULL;
 125   AliFittingQDistribution          *fqd   = NULL;
 126   AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
 127   AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
 128   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 129   AliFlowAnalysisWithScalarProduct *sp    = NULL;
 130   AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;
 131
 132   //MCEP = monte carlo event plane
 133   if (MCEP) {
 134     AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 135     mcep->Init();
 136   }
 137
 138   //QC = Q-cumulants
 139   if(QC) {
 140     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 141     qc->Init();
 142     if(listWithWeights) qc->SetWeightsList(listWithWeights);
 143     if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 144     if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 145     if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 146   }
 147
 148   //GFC = Generating Function Cumulants
 149   if(GFC) {
 150     AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 151     gfc->Init();
 152     if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 153     if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 154     if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 155     if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 156   }
 157
 158   //FQD = Fitting q-distribution
 159   if(FQD) {
 160     AliFittingQDistribution* fqd = new AliFittingQDistribution();
 161     fqd->Init();
 162     if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 163     if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 164   }
 165
 166   //SP = Scalar Product
 167   if(SP) {
 168     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 169     sp->Init();
 170   }
 171
 172   //LYZ1 = Lee-Yang Zeroes first run
 173   if(LYZ1) {
 174     AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
 175     lyz1->SetFirstRun(kTRUE);
 176     lyz1->SetUseSum(kTRUE);
 177     lyz1->Init();
 178   }
 179
 180   //LYZ2 = Lee-Yang Zeroes second run
 181   if(LYZ2) {
 182     AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
 183     // read the input file from the first run
 184     TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
 185     TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
 186     if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
 187       cerr << " ERROR: NO First Run file... " << endl ;
 188       break;
 189     }
 190     else {
 191       TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
 192       if (!inputListLYZ2) {cout<<"Input list is NULL pointer!"<<endl; break;}
 193       else {
 194         cout<<"LYZ2 input file/list read..."<<endl;
 195         lyz2->SetFirstRunList(inputListLYZ2);
 196         lyz2->SetFirstRun(kFALSE);
 197         lyz2->SetUseSum(kTRUE);
 198         lyz2->Init();
 199       }
 200     }
 201   }
 202
 203  //LYZEP = Lee-Yang Zeroes event plane
 204   if(LYZEP) {
 205     AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 206     AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 207     // read the input file from the second lyz run
 208     TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
 209     TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 210     if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 211       cerr << " ERROR: NO Second Run file... " << endl ;
 212       break;
 213     }
 214     else {
 215       TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
 216       if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
 217       else {
 218         cout<<"LYZEP input file/list read..."<<endl;
 219         ep   ->SetSecondRunList(inputListLYZEP);
 220         lyzep->SetSecondRunList(inputListLYZEP);
 221         ep   ->Init();
 222         lyzep->Init();
 223       }
 224     }
 225   }
 226   //------------------------------------------------------------------------
 227
 228
 229   //standard code
 230   Int_t fCount = 0;
 231   TString execDir(gSystem->pwd());
 232   TString targetDir(dir);
 233   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
 234   TList* dirList           = baseDir->GetListOfFiles();
 235   if (!dirList) {
 236     cout << endl << "No input files in: " << targetDir.Data() << endl;
 237     break;
 238   }
 239   Int_t nDirs              = dirList->GetEntries();
 240   cout<<endl;
 241   cout<<"Int_t nDirs = "<<nDirs<<endl;
 242   gSystem->cd(execDir);
 243
 244   for(Int_t iDir=0;iDir<nDirs;++iDir)
 245     {
 246       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
 247       if(!presentDir || !presentDir->IsDirectory() ||
 248          strcmp(presentDir->GetName(), ".") == 0 ||
 249          strcmp(presentDir->GetName(), "..") == 0)
 250         {
 251           cout << endl;
 252           cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
 253             " - Skipping ... " << endl;
 254           continue ;
 255         }
 256
 257       if(offset > 0) { --offset ; continue ; }
 258       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
 259
 260       TString presentDirName(dir); // aDataDir
 261       presentDirName += presentDir->GetName();
 262       presentDirName += "/";
 263       //cerr<<" presentDirName = "<<presentDirName<<endl;
 264
 265       TString fileName = presentDirName;
 266       fileName += "galice.root";
 267       Long_t *id, *size, *flags, *modtime;
 268       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
 269         {
 270           cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
 271                << endl;
 272           continue;
 273         }
 274       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
 275
 276       //loop (simulations in the present dir)
 277       TSystemDirectory* evtsDir = new TSystemDirectory(".",
 278                                                        presentDirName.Data());
 279       TList* fileList       = evtsDir->GetListOfFiles();
 280       Int_t nFiles                  = fileList->GetEntries();
 281       //cout<<" Int_t nFiles = "<<nFiles<<endl;
 282       gSystem->cd(execDir);
 283       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
 284         {
 285           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
 286           TString presentFileName(presentDirName);
 287           presentFileName += presentFile->GetName();
 288
 289           if(!(presentFileName.Contains("Kinematics") &&
 290                presentFileName.Contains("root"))) { continue ; }
 291
 292           //cout << " found: " << presentFileName.Data() << endl;
 293
 294           TFile* kineFile = new TFile(presentFileName.Data(), "READ");
 295           // kineFile->ls();
 296           Int_t nEvts = kineFile->GetNkeys() ;
 297           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
 298           TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
 299           TTree* kTree;
 300           TIter next(kineEventsList);
 301           TKey* key;
 302
 303           //loop over the events
 304           while( key=(TKey *)next() )
 305             {
 306               TDirectory* tDir = (TDirectory*)key->ReadObj();
 307               if(!tDir) break;
 308
 309               TString evtDir(tDir->GetName());
 310               //cout << "  . . found: " << tDir->GetName() << endl;
 311
 312               kTree = (TTree *)tDir->Get("TreeK");
 313               if(!kTree) break;
 314
 315               Int_t nPart = kTree->GetEntries();
 316               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl;
 317
 318               //-----------------------------------------------------------
 319               //fill and save the flow event
 320               AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
 321
 322               if(MCEP) mcep->Make(fEvent);
 323               if(QC) qc->Make(fEvent);
 324               if(GFC) gfc->Make(fEvent);
 325               if(FQD) fqd->Make(fEvent);
 326               if(LYZ1) lyz1->Make(fEvent);
 327               if(LYZ2) lyz2->Make(fEvent);
 328               if(LYZEP) lyzep->Make(fEvent,ep);
 329               if(SP) sp->Make(fEvent);
 330               //-----------------------------------------------------------
 331               fCount++;
 332               //cout << "# " << fCount << " events processed" << endl;
 333               delete kTree;
 334               delete fEvent;
 335             }
 336           delete kineFile ;
 337         }
 338       delete evtsDir ;
 339     }
 340
 341   //--------------------------------------------------------------
 342   //calculating and storing the final results of flow analysis
 343   if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
 344   if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
 345   if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
 346   if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
 347   if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
 348   if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
 349   if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
 350   if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
 351
 352   //--------------------------------------------------------------
 353
 354   cout << endl;
 355   cout << " Fini ... " << endl;
 356   cout << endl;
 357
 358   timer.Stop();
 359   cout << endl;
 360   timer.Print();
 361 }
 362
 363 void SetupPar(char* pararchivename)
 364 {
 365   //Load par files, create analysis libraries
 366   //For testing, if par file already decompressed and modified
 367   //classes then do not decompress.
 368
 369   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 370   TString parpar(Form("%s.par", pararchivename)) ;
 371   if ( gSystem->AccessPathName(parpar.Data()) ) {
 372     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 373     TString processline(Form(".! make %s", parpar.Data())) ;
 374     gROOT->ProcessLine(processline.Data()) ;
 375     gSystem->ChangeDirectory(cdir) ;
 376     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 377     gROOT->ProcessLine(processline.Data()) ;
 378   }
 379   if ( gSystem->AccessPathName(pararchivename) ) {
 380     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 381     gROOT->ProcessLine(processline.Data());
 382   }
 383
 384   TString ocwd = gSystem->WorkingDirectory();
 385   gSystem->ChangeDirectory(pararchivename);
 386
 387   // check for BUILD.sh and execute
 388   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 389     printf("*******************************\n");
 390     printf("*** Building PAR archive    ***\n");
 391     cout<<pararchivename<<endl;
 392     printf("*******************************\n");
 393
 394     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 395       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 396       return -1;
 397     }
 398   }
 399   // check for SETUP.C and execute
 400   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 401     printf("*******************************\n");
 402     printf("*** Setup PAR archive       ***\n");
 403     cout<<pararchivename<<endl;
 404     printf("*******************************\n");
 405     gROOT->Macro("PROOF-INF/SETUP.C");
 406   }
 407
 408   gSystem->ChangeDirectory(ocwd.Data());
 409   printf("Current dir: %s\n", ocwd.Data());
 410 }
 411
 412 void LoadLibraries(const anaModes mode) {
 413
 414   //--------------------------------------
 415   // Load the needed libraries most of them already loaded by aliroot
 416   //--------------------------------------
 417   gSystem->Load("libTree.so");
 418   gSystem->Load("libGeom.so");
 419   gSystem->Load("libVMC.so");
 420   gSystem->Load("libXMLIO.so");
 421   gSystem->Load("libPhysics.so");
 422
 423   //----------------------------------------------------------
 424   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 425   //----------------------------------------------------------
 426   if (mode==mLocal) {
 427     //--------------------------------------------------------
 428     // If you want to use already compiled libraries
 429     // in the aliroot distribution
 430     //--------------------------------------------------------
 431     gSystem->Load("libSTEERBase");
 432     gSystem->Load("libESD");
 433     gSystem->Load("libAOD");
 434     gSystem->Load("libANALYSIS");
 435     gSystem->Load("libANALYSISalice");
 436     gSystem->Load("libCORRFW.so");
 437     cerr<<"libCORRFW.so loaded..."<<endl;
 438     gSystem->Load("libPWG2flowCommon.so");
 439     cerr<<"libPWG2flowCommon.so loaded..."<<endl;
 440     gSystem->Load("libPWG2flowTasks.so");
 441     cerr<<"libPWG2flowTasks.so loaded..."<<endl;
 442   }
 443
 444   else if (mode == mLocalPAR) {
 445     //--------------------------------------------------------
 446     //If you want to use root and par files from aliroot
 447     //--------------------------------------------------------
 448      //If you want to use root and par files from aliroot
 449     //--------------------------------------------------------
 450     SetupPar("STEERBase");
 451     SetupPar("ESD");
 452     SetupPar("AOD");
 453     SetupPar("ANALYSIS");
 454     SetupPar("ANALYSISalice");
 455     SetupPar("PWG2AOD");
 456     SetupPar("CORRFW");
 457     SetupPar("PWG2flowCommon");
 458     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 459     SetupPar("PWG2flowTasks");
 460     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 461   }
 462
 463   //---------------------------------------------------------
 464   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 465   //---------------------------------------------------------
 466   else if (mode==mLocalSource) {
 467
 468     // In root inline compile
 469
 470
 471     // Constants
 472     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 473     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 474     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 475
 476     // Flow event
 477     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 478     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 479     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 480
 481     // Cuts
 482     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 483
 484     // Output histosgrams
 485     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 486     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 487     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 488     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 489
 490     // Functions needed for various methods
 491     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 492     gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
 493     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 494
 495     // Flow Analysis code for various methods
 496     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 497     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 498     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 499     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 500     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 501     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 502     gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");
 503
 504     // Class to fill the FlowEvent without aliroot dependence
 505     // can be found in the directory FlowEventMakers
 506     gROOT->LoadMacro("FlowEventMakers/FlowEventSimpleMaker.cxx+");
 507
 508     cout << "finished loading macros!" << endl;
 509
 510   }
 511
 512 }
 513
 514