PWGCF/FLOW/macros/runFlowAnalysis.C

   1 #include <Riostream.h>
   2 #include <TStopwatch.h>
   3 #include <TObjArray.h>
   4 #include <TFile.h>
   5 #include <TRandom3.h>
   6 #include <TTimeStamp.h>
   7
   8 #include <Base/AliFlowCommonConstants.h>
   9 #include <Base/AliFlowLYZConstants.h>
  10 #include <Base/AliFlowCumuConstants.h>
  11 #include <Base/AliFlowVector.h>
  12 #include <Base/AliFlowTrackSimple.h>
  13 #include <Base/AliFlowEvent.h>
  14 #include <Base/AliFlowEventSimple.h>
  15 #include <Base/AliFlowTrackSimpleCuts.h>
  16 #include <Base/AliFlowCommonHist.h>
  17 #include <Base/AliFlowCommonHistResults.h>
  18 #include <Base/AliFlowLYZHist1.h>
  19 #include <Base/AliFlowLYZHist2.h>
  20 #include <Base/AliCumulantsFunctions.h>
  21 #include <Base/AliFlowLYZEventPlane.h>
  22 #include <Base/AliFlowAnalysisWithMCEventPlane.h>
  23 #include <Base/AliFlowAnalysisWithScalarProduct.h>
  24 #include <Base/AliFlowAnalysisWithLYZEventPlane.h>
  25 #include <Base/AliFlowAnalysisWithLeeYangZeros.h>
  26 #include <Base/AliFlowAnalysisWithCumulants.h>
  27 #include <Base/AliFlowAnalysisWithQCumulants.h>
  28 #include <Base/AliFlowAnalysisWithFittingQDistribution.h>
  29 #include <Base/AliFlowAnalysisWithMixedHarmonics.h>
  30 #include <Base/AliFlowAnalysisWithNestedLoops.h>
  31
  32 //--------------------------------------------------------------------------------------
  33 // Run flow analysis on local data with custom FlowEvent maker
  34 // RUN SETTINGS
  35 //flow analysis method can be: (set to kTRUE or kFALSE)
  36 Bool_t SP       = kTRUE;
  37 Bool_t LYZ1SUM  = kTRUE;
  38 Bool_t LYZ1PROD = kTRUE;
  39 Bool_t LYZ2SUM  = kFALSE;
  40 Bool_t LYZ2PROD = kFALSE;
  41 Bool_t LYZEP    = kFALSE;
  42 Bool_t GFC      = kTRUE;
  43 Bool_t QC       = kTRUE;
  44 Bool_t FQD      = kTRUE;
  45 Bool_t MH       = kTRUE;
  46 Bool_t NL       = kFALSE;
  47 Bool_t MCEP     = kFALSE; //does not work yet 24/12/08
  48 //--------------------------------------------------------------------------------------
  49
  50 // Weights
  51 // Use weights for Q vector
  52 Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
  53 Bool_t usePtWeights  = kFALSE; //v'(pt) (differential flow in pt)
  54 Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
  55
  56 //--------------------------------------------------------------------------------------
  57 // CUT SETTINGS
  58 //integrated selection
  59 Double_t ptMaxInt  = 10.;
  60 Double_t ptMinInt  = 0.;
  61 Double_t etaMaxInt = 1.;
  62 Double_t etaMinInt = -1.;
  63 Double_t phiMaxInt = 7.5;
  64 Double_t phiMinInt = 0.;
  65 Int_t PIDInt       = 211;
  66
  67 //differential selection
  68 Double_t ptMaxDiff  = 10.;
  69 Double_t ptMinDiff  = 0.;
  70 Double_t etaMaxDiff = 1.;
  71 Double_t etaMinDiff = -1.;
  72 Double_t phiMaxDiff = 7.5;
  73 Double_t phiMinDiff = 0.;
  74 Int_t PIDDiff       = 211;
  75 /*
  76 //--------------------------------------------------------------------------------------
  77 // FLOW SETTINGS (R.Rietkerk)
  78 Int_t nLoops=1;                 // Number of times to use the same particle (nonflow).
  79 Double_t xEllipticFlowValue=0.1;// Add Elliptic Flow. Must be in range [0,1].
  80 Int_t nMultiplicityOfEvent=500; // Set Average Multiplicity.
  81 Double_t xSigmaFlow=0.00;       // Add Elliptic Flow. Must be in range [0,1].
  82 Int_t nSigmaMult=50;            // Set Average Multiplicity.
  83 //--------------------------------------------------------------------------------------
  84 */
  85
  86 enum anaModes {mLocal,mLocalSource,mLocalPAR,};
  87 //mLocal: Analyze data on your computer using aliroot
  88 //mLocalPAR: Analyze data on your computer using root + PAR files
  89 //mLocalSource: Analyze data on your computer using root + source files
  90
  91 void LoadLibraries(const anaModes mode);
  92
  93 Int_t offset = 0;
  94
  95 int runFlowAnalysis(const anaModes mode=mLocal, Int_t aRuns = 100, const char*
  96                     dir="/data/alice1/kolk/KineOnly3/")
  97                     //              dir="/Users/snelling/alice_data/KineOnly3/")
  98                     //              dir="/Users/snelling/alice_data/stoomboot/5b/")
  99 {
 100   TStopwatch timer;
 101   timer.Start();
 102
 103   if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
 104   if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
 105   if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
 106   if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
 107
 108
 109   cout<<endl;
 110   cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
 111   cout<<endl;
 112
 113   LoadLibraries(mode);
 114
 115   TRandom3 random3Temp; //init for manual settings (R.Rietkerk)
 116   TTimeStamp dt;
 117   Int_t sseed = dt.GetNanoSec()/1000;
 118   random3Temp.SetSeed(sseed);
 119
 120   if (mode == mLocal || mode == mLocalPAR) {
 121   }
 122   else if (mode == mLocalSource) {
 123   }
 124   else{
 125     cout << "No supported running mode selected!" << endl;
 126     break;
 127   }
 128
 129
 130   //------------------------------------------------------------------------
 131   //cuts:
 132   AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
 133   cutsInt->SetPtMax(ptMaxInt);
 134   cutsInt->SetPtMin(ptMinInt);
 135   cutsInt->SetEtaMax(etaMaxInt);
 136   cutsInt->SetEtaMin(etaMinInt);
 137   cutsInt->SetPhiMax(phiMaxInt);
 138   cutsInt->SetPhiMin(phiMinInt);
 139   cutsInt->SetPID(PIDInt);
 140
 141   AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
 142   cutsDiff->SetPtMax(ptMaxDiff);
 143   cutsDiff->SetPtMin(ptMinDiff);
 144   cutsDiff->SetEtaMax(etaMaxDiff);
 145   cutsDiff->SetEtaMin(etaMinDiff);
 146   cutsDiff->SetPhiMax(phiMaxDiff);
 147   cutsDiff->SetPhiMin(phiMinDiff);
 148   cutsDiff->SetPID(PIDDiff);
 149
 150   //if the weights are used:
 151   TFile *fileWithWeights = NULL;
 152   TList *listWithWeights = NULL;
 153
 154   if(usePhiWeights||usePtWeights||useEtaWeights) {
 155     fileWithWeights = TFile::Open("weights.root","READ");
 156     if(fileWithWeights) {
 157       listWithWeights = (TList*)fileWithWeights->Get("weights");
 158     }
 159     else
 160       {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 161         break;
 162       }
 163   }
 164
 165   //flow methods:
 166   AliFlowAnalysisWithQCumulants *qc = NULL;
 167   AliFlowAnalysisWithCumulants *gfc = NULL;
 168   AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
 169   AliFlowAnalysisWithLeeYangZeros *lyz1sum = NULL;
 170   AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
 171   AliFlowAnalysisWithLeeYangZeros *lyz2sum = NULL;
 172   AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
 173   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 174   AliFlowAnalysisWithScalarProduct *sp = NULL;
 175   AliFlowAnalysisWithMCEventPlane *mcep = NULL;
 176   AliFlowAnalysisWithMixedHarmonics *mh = NULL;
 177   AliFlowAnalysisWithNestedLoops *nl = NULL;
 178
 179   //MCEP = monte carlo event plane
 180   if (MCEP) {
 181     AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 182     mcep->Init();
 183   }
 184
 185   //QC = Q-cumulants
 186   if(QC) {
 187     AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 188     if(listWithWeights) qc->SetWeightsList(listWithWeights);
 189     if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 190     if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 191     if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 192     qc->Init();
 193   }
 194
 195   //GFC = Generating Function Cumulants
 196   if(GFC) {
 197     AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 198     if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 199     if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 200     if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 201     if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 202     gfc->Init();
 203   }
 204
 205   //FQD = Fitting q-distribution
 206   if(FQD) {
 207     AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 208     if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 209     if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 210     fqd->Init();
 211   }
 212
 213   //SP = Scalar Product
 214   if(SP) {
 215     AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 216     if(usePhiWeights) sp->SetUsePhiWeights(usePhiWeights);
 217     sp->Init();
 218   }
 219
 220   //LYZ1 = Lee-Yang Zeroes first run
 221   if(LYZ1SUM) {
 222     AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 223     lyz1sum->SetFirstRun(kTRUE);
 224     lyz1sum->SetUseSum(kTRUE);
 225     lyz1sum->Init();
 226   }
 227   if(LYZ1PROD) {
 228     AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 229     lyz1prod->SetFirstRun(kTRUE);
 230     lyz1prod->SetUseSum(kFALSE);
 231     lyz1prod->Init();
 232   }
 233   //LYZ2 = Lee-Yang Zeroes second run
 234   if(LYZ2SUM) {
 235     AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 236     // read the input file from the first run
 237     TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 238     TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 239     if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 240       cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 241       break;
 242     }
 243     else {
 244       TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 245       if (!inputListLYZ2SUM) {cout<<"SUM Input list is NULL pointer!"<<endl; break;}
 246       else {
 247         cout<<"LYZ2SUM input file/list read..."<<endl;
 248         lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 249         lyz2sum->SetFirstRun(kFALSE);
 250         lyz2sum->SetUseSum(kTRUE);
 251         lyz2sum->Init();
 252       }
 253     }
 254   }
 255   if(LYZ2PROD) {
 256     AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 257     // read the input file from the first run
 258     TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 259     TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 260     if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 261       cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 262       break;
 263     }
 264     else {
 265       TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 266       if (!inputListLYZ2PROD) {cout<<"PROD Input list is NULL pointer!"<<endl; break;}
 267       else {
 268         cout<<"LYZ2PROD input file/list read..."<<endl;
 269         lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 270         lyz2prod->SetFirstRun(kFALSE);
 271         lyz2prod->SetUseSum(kTRUE);
 272         lyz2prod->Init();
 273       }
 274     }
 275   }
 276  //LYZEP = Lee-Yang Zeroes event plane
 277   if(LYZEP) {
 278     AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 279     AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 280     // read the input file from the second lyz run
 281     TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 282     TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 283     if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 284       cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 285       break;
 286     }
 287     else {
 288       TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 289       if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
 290       else {
 291         cout<<"LYZEP input file/list read..."<<endl;
 292         ep   ->SetSecondRunList(inputListLYZEP);
 293         lyzep->SetSecondRunList(inputListLYZEP);
 294         ep   ->Init();
 295         lyzep->Init();
 296       }
 297     }
 298   }
 299   // MH = Mixed Harmonics:
 300   if(MH) {
 301     AliFlowAnalysisWithMixedHarmonics* mh = new AliFlowAnalysisWithMixedHarmonics();
 302     if(listWithWeights) mh->SetWeightsList(listWithWeights);
 303     //if(usePhiWeights) mh->SetUsePhiWeights(usePhiWeights); // to be improved (enabled)
 304     //if(usePtWeights) mh->SetUsePtWeights(usePtWeights); // to be improved (enabled)
 305     //if(useEtaWeights) mh->SetUseEtaWeights(useEtaWeights); // to be improved (enabled)
 306     mh->Init();
 307   }
 308   // NL = Nested Loops:
 309   if(NL) {
 310     AliFlowAnalysisWithNestedLoops* nl = new AliFlowAnalysisWithNestedLoops();
 311     if(listWithWeights) nl->SetWeightsList(listWithWeights);
 312     //if(usePhiWeights) nl->SetUsePhiWeights(usePhiWeights); // to be improved (enabled)
 313     //if(usePtWeights) nl->SetUsePtWeights(usePtWeights); // to be improved (enabled)
 314     //if(useEtaWeights) nl->SetUseEtaWeights(useEtaWeights); // to be improved (enabled)
 315     nl->Init();
 316   }
 317
 318   //------------------------------------------------------------------------
 319
 320
 321   //standard code to read files in directory
 322   Int_t fCount = 0;
 323   TString execDir(gSystem->pwd());
 324   TString targetDir(dir);
 325   TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
 326   TList* dirList           = baseDir->GetListOfFiles();
 327   if (!dirList) {
 328     cout << endl << "No input files in: " << targetDir.Data() << endl;
 329     break;
 330   }
 331   Int_t nDirs              = dirList->GetEntries();
 332   cout<<endl;
 333   cout<<"Int_t nDirs = "<<nDirs<<endl;
 334   gSystem->cd(execDir);
 335
 336   for(Int_t iDir=0;iDir<nDirs;++iDir)
 337     {
 338       TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
 339       if(!presentDir || !presentDir->IsDirectory() ||
 340          strcmp(presentDir->GetName(), ".") == 0 ||
 341          strcmp(presentDir->GetName(), "..") == 0)
 342         {
 343           cout << endl;
 344           cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
 345             " - Skipping ... " << endl;
 346           continue ;
 347         }
 348
 349       if(offset > 0) { --offset ; continue ; }
 350       if((aRuns > 0) && (fCount >= aRuns)) { break ; }
 351
 352       TString presentDirName(dir); // aDataDir
 353       presentDirName += presentDir->GetName();
 354       presentDirName += "/";
 355       //cerr<<" presentDirName = "<<presentDirName<<endl;
 356
 357       TString fileName = presentDirName;
 358       fileName += "galice.root";
 359       Long_t *id, *size, *flags, *modtime;
 360       if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
 361         {
 362           cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
 363                << endl;
 364           continue;
 365         }
 366       //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
 367
 368       //loop (simulations in the present dir)
 369       TSystemDirectory* evtsDir = new TSystemDirectory(".",
 370                                                        presentDirName.Data());
 371       TList* fileList       = evtsDir->GetListOfFiles();
 372       Int_t nFiles                  = fileList->GetEntries();
 373       //cout<<" Int_t nFiles = "<<nFiles<<endl;
 374       gSystem->cd(execDir);
 375       for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
 376         {
 377           TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
 378           TString presentFileName(presentDirName);
 379           presentFileName += presentFile->GetName();
 380
 381           if(!(presentFileName.Contains("Kinematics") &&
 382                presentFileName.Contains("root"))) { continue ; }
 383
 384           //cout << " found: " << presentFileName.Data() << endl;
 385
 386           TFile* kineFile = new TFile(presentFileName.Data(), "READ");
 387           // kineFile->ls();
 388           Int_t nEvts = kineFile->GetNkeys() ;
 389           //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
 390           TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
 391           TTree* kTree;
 392           TIter next(kineEventsList);
 393           TKey* key;
 394
 395           Double_t xRPAngle;
 396           //loop over the events
 397           while( key=(TKey *)next() )
 398             {
 399               TDirectory* tDir = (TDirectory*)key->ReadObj();
 400               if(!tDir) break;
 401
 402               TString evtDir(tDir->GetName());
 403               //cout << "  . . found: " << tDir->GetName() << endl;
 404
 405               kTree = (TTree *)tDir->Get("TreeK");
 406               if(!kTree) break;
 407
 408               Int_t nPart = kTree->GetEntries();
 409               //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl;
 410
 411               //-----------------------------------------------------------
 412               //fill and save the flow event
 413
 414               /*
 415               Int_t nNewMultOfEvent = random3Temp.Gaus(nMultiplicityOfEvent,nSigmaMult);
 416               cout << "new multiplicity: " << nNewMultOfEvent << endl;
 417               Double_t xNewFlowValue = random3Temp.Gaus(xEllipticFlowValue,xSigmaFlow);
 418               if ( (fCount % 100) == 0) {
 419                 cout << "new multiplicity: " << nNewMultOfEvent << endl;
 420                 cout << "new flow value: " << xNewFlowValue << endl;
 421               }
 422               */
 423
 424               AliFlowEventSimple *fEvent = new AliFlowEventSimple(kTree, cutsInt, cutsDiff);
 425               //xRPAngle=random3Temp.Uniform(0.0,TMath::TwoPi());
 426         //fEvent->SetMCReactionPlaneAngle(xRPAngle);
 427         //fEvent->SetV2(xNewFlowValue);
 428         //fEvent->CloneTracks(nLoops);
 429
 430               // do flow analysis for various methods
 431               if(MCEP)    mcep->Make(fEvent);
 432               if(QC)      qc->Make(fEvent);
 433               if(GFC)     gfc->Make(fEvent);
 434               if(FQD)     fqd->Make(fEvent);
 435               if(LYZ1SUM) lyz1sum->Make(fEvent);
 436               if(LYZ1PROD)lyz1prod->Make(fEvent);
 437               if(LYZ2SUM) lyz2sum->Make(fEvent);
 438               if(LYZ2PROD)lyz2prod->Make(fEvent);
 439               if(LYZEP)   lyzep->Make(fEvent,ep);
 440               if(SP)      sp->Make(fEvent);
 441               if(MH)      mh->Make(fEvent);
 442               if(NL)      nl->Make(fEvent);
 443               //-----------------------------------------------------------
 444               fCount++;
 445               //cout << "# " << fCount << " events processed" << endl;
 446               delete kTree;
 447               delete fEvent;
 448             }
 449           delete kineFile ;
 450         }
 451       delete evtsDir ;
 452     }
 453
 454  //---------------------------------------------------------------------------------------
 455  // create a new file which will hold the final results of all methods:
 456  TString outputFileName = "AnalysisResults.root";
 457  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 458  // create a new file for each method wich will hold list with final results:
 459  const Int_t nMethods = 12;
 460  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL"};
 461  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 462  TString fileNameMethod[nMethods];
 463  for(Int_t i=0;i<nMethods;i++)
 464  {
 465   // form a file name for each method:
 466   fileNameMethod[i]+="output";
 467   fileNameMethod[i]+=method[i].Data();
 468   fileNameMethod[i]+="analysis";
 469   dirFileFinal[i] = new TDirectoryFile(fileNameMethod[i].Data(),fileNameMethod[i].Data());
 470  }
 471
 472  // calculating and storing the final results of default flow analysis:
 473  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
 474  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
 475  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
 476  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
 477  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
 478  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
 479  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
 480  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
 481  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
 482  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
 483  if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
 484  if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
 485  //---------------------------------------------------------------------------------------
 486
 487  outputFile->Close();
 488  delete outputFile;
 489
 490   cout << endl;
 491   cout << " Fini ... " << endl;
 492   cout << endl;
 493
 494   timer.Stop();
 495   cout << endl;
 496   timer.Print();
 497 }
 498
 499 void SetupPar(char* pararchivename)
 500 {
 501   //Load par files, create analysis libraries
 502   //For testing, if par file already decompressed and modified
 503   //classes then do not decompress.
 504
 505   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 506   TString parpar(Form("%s.par", pararchivename)) ;
 507   if ( gSystem->AccessPathName(parpar.Data()) ) {
 508     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 509     TString processline(Form(".! make %s", parpar.Data())) ;
 510     gROOT->ProcessLine(processline.Data()) ;
 511     gSystem->ChangeDirectory(cdir) ;
 512     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 513     gROOT->ProcessLine(processline.Data()) ;
 514   }
 515   if ( gSystem->AccessPathName(pararchivename) ) {
 516     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 517     gROOT->ProcessLine(processline.Data());
 518   }
 519
 520   TString ocwd = gSystem->WorkingDirectory();
 521   gSystem->ChangeDirectory(pararchivename);
 522
 523   // check for BUILD.sh and execute
 524   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 525     printf("*******************************\n");
 526     printf("*** Building PAR archive    ***\n");
 527     cout<<pararchivename<<endl;
 528     printf("*******************************\n");
 529
 530     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 531       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 532       return -1;
 533     }
 534   }
 535   // check for SETUP.C and execute
 536   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 537     printf("*******************************\n");
 538     printf("*** Setup PAR archive       ***\n");
 539     cout<<pararchivename<<endl;
 540     printf("*******************************\n");
 541     gROOT->Macro("PROOF-INF/SETUP.C");
 542   }
 543
 544   gSystem->ChangeDirectory(ocwd.Data());
 545   printf("Current dir: %s\n", ocwd.Data());
 546 }
 547
 548 void LoadLibraries(const anaModes mode) {
 549
 550   //--------------------------------------
 551   // Load the needed libraries most of them already loaded by aliroot
 552   //--------------------------------------
 553   //gSystem->Load("libTree");
 554   gSystem->Load("libGeom");
 555   gSystem->Load("libVMC");
 556   gSystem->Load("libXMLIO");
 557   gSystem->Load("libPhysics");
 558
 559   //----------------------------------------------------------
 560   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 561   //----------------------------------------------------------
 562   if (mode==mLocal) {
 563     //--------------------------------------------------------
 564     // If you want to use already compiled libraries
 565     // in the aliroot distribution
 566     //--------------------------------------------------------
 567     gSystem->Load("libSTEERBase");
 568     gSystem->Load("libESD");
 569     gSystem->Load("libAOD");
 570     gSystem->Load("libANALYSIS");
 571     gSystem->Load("libANALYSISalice");
 572     gSystem->Load("libCORRFW");
 573     cerr<<"libCORRFW loaded..."<<endl;
 574     gSystem->Load("libPWGflowBase");
 575     cerr<<"libPWGflowBase loaded..."<<endl;
 576     gSystem->Load("libPWGflowTasks");
 577     cerr<<"libPWGflowTasks loaded..."<<endl;
 578   }
 579
 580   else if (mode == mLocalPAR) {
 581     //--------------------------------------------------------
 582     //If you want to use root and par files from aliroot
 583     //--------------------------------------------------------
 584      //If you want to use root and par files from aliroot
 585     //--------------------------------------------------------
 586     SetupPar("STEERBase");
 587     SetupPar("ESD");
 588     SetupPar("AOD");
 589     SetupPar("ANALYSIS");
 590     SetupPar("ANALYSISalice");
 591     SetupPar("PWG2AOD");
 592     SetupPar("CORRFW");
 593     SetupPar("PWGflowBase");
 594     cerr<<"PWGflowBase.par loaded..."<<endl;
 595     SetupPar("PWGflowTasks");
 596     cerr<<"PWGflowTasks.par loaded..."<<endl;
 597   }
 598
 599   //---------------------------------------------------------
 600   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 601   //---------------------------------------------------------
 602   else if (mode==mLocalSource) {
 603
 604     // In root inline compile
 605
 606
 607     // Constants
 608     gROOT->LoadMacro("Base/AliFlowCommonConstants.cxx+");
 609     gROOT->LoadMacro("Base/AliFlowLYZConstants.cxx+");
 610
 611     // Flow event
 612     gROOT->LoadMacro("Base/AliFlowVector.cxx+");
 613     gROOT->LoadMacro("Base/AliFlowTrackSimple.cxx+");
 614     gROOT->LoadMacro("Base/AliFlowTrackSimpleCuts.cxx+");
 615     gROOT->LoadMacro("Base/AliFlowEventSimple.cxx+");
 616
 617     // Output histosgrams
 618     gROOT->LoadMacro("Base/AliFlowCommonHist.cxx+");
 619     gROOT->LoadMacro("Base/AliFlowCommonHistResults.cxx+");
 620     gROOT->LoadMacro("Base/AliFlowLYZHist1.cxx+");
 621     gROOT->LoadMacro("Base/AliFlowLYZHist2.cxx+");
 622
 623     // Functions needed for various methods
 624     gROOT->LoadMacro("Base/AliCumulantsFunctions.cxx+");
 625     gROOT->LoadMacro("Base/AliFlowLYZEventPlane.cxx+");
 626
 627     // Flow Analysis code for various methods
 628     gROOT->LoadMacro("Base/AliFlowAnalysisWithMCEventPlane.cxx+");
 629     gROOT->LoadMacro("Base/AliFlowAnalysisWithScalarProduct.cxx+");
 630     gROOT->LoadMacro("Base/AliFlowAnalysisWithLYZEventPlane.cxx+");
 631     gROOT->LoadMacro("Base/AliFlowAnalysisWithLeeYangZeros.cxx+");
 632     gROOT->LoadMacro("Base/AliFlowAnalysisWithCumulants.cxx+");
 633     gROOT->LoadMacro("Base/AliFlowAnalysisWithQCumulants.cxx+");
 634     gROOT->LoadMacro("Base/AliFlowAnalysisWithFittingQDistribution.cxx+");
 635     gROOT->LoadMacro("Base/AliFlowAnalysisWithMixedHarmonics.cxx+");
 636     gROOT->LoadMacro("Base/AliFlowAnalysisWithNestedLoops.cxx+");
 637
 638     cout << "finished loading macros!" << endl;
 639
 640   }
 641
 642 }
 643
 644