PWG2/FLOW/macros/runFlowAnalysisOnTheFly.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  = kTRUE;
  17 //--------------------------------------------------------------------------------------
  18
  19 Bool_t bSameSeed = kFALSE; // use always the same seed for random generators
  20 Bool_t bConstantHarmonics = kFALSE; // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
  21
  22 // Set the event parameters:
  23 Int_t iLoops = 1; // number of times to use each track (to simulate nonflow)
  24
  25 Int_t iMultiplicityOfRP = 500; // multiplicity of RPs
  26 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs
  27 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs)
  28
  29 //......................................................................................
  30 // if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
  31 Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
  32 Double_t dV2RPMax = 0.20; // maximum value of V2(pt) for pt >= 2GeV
  33 //......................................................................................
  34
  35 //......................................................................................
  36 // if you use constant harmonics (bConstantHarmonics = kTRUE):
  37 Double_t dV2RP = 0.05; // elliptic flow of RPs
  38 Double_t dV2SpreadRP = 0.; // elliptic flow spread of RPs
  39
  40 Double_t dV1RP = 0.0; // directed flow of RPs
  41 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
  42
  43 Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
  44 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
  45 //......................................................................................
  46
  47 enum anaModes {mLocal,mLocalSource,mLocalPAR};
  48 // mLocal: Analyze data on your computer using aliroot
  49 // mLocalPAR: Analyze data on your computer using root + PAR files
  50 // mLocalSource: Analyze data on your computer using root + source files
  51
  52 int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=100)
  53 {
  54  TStopwatch timer;
  55  timer.Start();
  56
  57  if (LYZ1 && LYZ2)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
  58  if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  59  if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  60
  61  cout<<endl;
  62  cout<<endl;
  63  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
  64  cout<<endl;
  65
  66  cout<<endl;
  67  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
  68  cout<<endl;
  69  cout<<endl;
  70
  71  LoadLibraries(mode);
  72
  73  // Initialize the seed for random generator
  74  UInt_t sseed = 0;
  75
  76  if(bSameSeed)
  77  {
  78   sseed = 44; // the default constant value for seed for random generators
  79  }
  80
  81  if(!bSameSeed)
  82  {
  83   TTimeStamp dt;
  84   sseed = dt.GetNanoSec()/1000;
  85  }
  86
  87  //---------------------------------------------------------------------------------------
  88  // Initialize the flowevent maker
  89  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
  90  eventMakerOnTheFly->Init();
  91
  92  //---------------------------------------------------------------------------------------
  93  // Initialize all the flow methods:
  94  AliFlowAnalysisWithQCumulants    *qc    = NULL;
  95  AliFlowAnalysisWithCumulants     *gfc   = NULL;
  96  AliFittingQDistribution          *fqd   = NULL;
  97  AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
  98  AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
  99  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 100  AliFlowAnalysisWithScalarProduct *sp    = NULL;
 101  AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;
 102
 103  // MCEP = monte carlo event plane
 104  if (MCEP) {
 105    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 106    mcep->Init();
 107  }
 108
 109   // QC = Q-cumulants
 110  if(QC) {
 111    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 112    qc->Init();
 113  }
 114
 115  // GFC = Generating Function Cumulants
 116  if(GFC) {
 117    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 118    gfc->Init();
 119  }
 120
 121  // FQD = Fitting q-distribution
 122  if(FQD) {
 123    AliFittingQDistribution* fqd = new AliFittingQDistribution();
 124    fqd->Init();
 125  }
 126
 127  // SP = Scalar Product
 128  if(SP) {
 129    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 130    sp->Init();
 131  }
 132
 133  // LYZ1 = Lee-Yang Zeroes first run
 134  if(LYZ1) {
 135    AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
 136    lyz1->SetFirstRun(kTRUE);
 137    lyz1->SetUseSum(kTRUE);
 138    lyz1->Init();
 139  }
 140
 141  // LYZ2 = Lee-Yang Zeroes second run
 142  if(LYZ2) {
 143    AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
 144    // read the input file from the first run
 145    TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
 146    TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
 147    if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) {
 148      cerr << " ERROR: NO First Run file... " << endl ;
 149      break;
 150    }
 151    else {
 152      TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
 153      if (!inputListLYZ2) {cout<<"Input list is NULL pointer!"<<endl; break;}
 154      else {
 155        cout<<"LYZ2 input file/list read..."<<endl;
 156        lyz2->SetFirstRunList(inputListLYZ2);
 157        lyz2->SetFirstRun(kFALSE);
 158        lyz2->SetUseSum(kTRUE);
 159        lyz2->Init();
 160      }
 161    }
 162  }
 163
 164  // LYZEP = Lee-Yang Zeroes event plane
 165  if(LYZEP) {
 166    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 167    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 168    // read the input file from the second lyz run
 169    TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
 170    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 171    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 172      cerr << " ERROR: NO Second Run file... " << endl ;
 173      break;
 174    }
 175    else {
 176      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
 177      if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
 178      else {
 179        cout<<"LYZEP input file/list read..."<<endl;
 180        ep   ->SetSecondRunList(inputListLYZEP);
 181        lyzep->SetSecondRunList(inputListLYZEP);
 182        ep   ->Init();
 183        lyzep->Init();
 184      }
 185    }
 186  }
 187  //---------------------------------------------------------------------------------------
 188
 189  // set the global event parameters:
 190  eventMakerOnTheFly->SetNoOfLoops(iLoops);
 191  eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
 192  eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
 193  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
 194
 195  eventMakerOnTheFly->SetV1RP(dV1RP);
 196  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
 197  eventMakerOnTheFly->SetV4RP(dV4RP);
 198  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
 199
 200  // constant harmonic V2:
 201  if(bConstantHarmonics)
 202  {
 203   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 204   eventMakerOnTheFly->SetV2RP(dV2RP);
 205   eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);
 206  }
 207  // (pt,eta) dependent harmonic V2:
 208  if(!bConstantHarmonics)
 209  {
 210   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 211   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
 212   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);
 213  }
 214
 215  //---------------------------------------------------------------------------------------
 216  // create and analyze events 'on the fly':
 217
 218  for(Int_t i=0;i<nEvts;i++) {
 219    // creating the event with above settings:
 220    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly();
 221
 222    // analyzing the created event 'on the fly':
 223    // do flow analysis for various methods:
 224    if(MCEP) mcep->Make(event);
 225    if(QC) qc->Make(event);
 226    if(GFC) gfc->Make(event);
 227    if(FQD) fqd->Make(event);
 228    if(LYZ1) lyz1->Make(event);
 229    if(LYZ2) lyz2->Make(event);
 230    if(LYZEP) lyzep->Make(event,ep);
 231    if(SP) sp->Make(event);
 232
 233    delete event;
 234  } // end of for(Int_t i=0;i<nEvts;i++)
 235  //---------------------------------------------------------------------------------------
 236
 237
 238
 239  //---------------------------------------------------------------------------------------
 240  // calculating and storing the final results of flow analysis
 241  if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
 242  if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
 243  if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
 244  if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
 245  if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
 246  if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
 247  if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
 248  if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
 249  //---------------------------------------------------------------------------------------
 250
 251
 252
 253  cout<<endl;
 254  cout<<endl;
 255  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 256  cout<<endl;
 257
 258  timer.Stop();
 259  cout << endl;
 260  timer.Print();
 261 }
 262
 263 void SetupPar(char* pararchivename)
 264 {
 265   //Load par files, create analysis libraries
 266   //For testing, if par file already decompressed and modified
 267   //classes then do not decompress.
 268
 269   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 270   TString parpar(Form("%s.par", pararchivename)) ;
 271   if ( gSystem->AccessPathName(parpar.Data()) ) {
 272     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 273     TString processline(Form(".! make %s", parpar.Data())) ;
 274     gROOT->ProcessLine(processline.Data()) ;
 275     gSystem->ChangeDirectory(cdir) ;
 276     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 277     gROOT->ProcessLine(processline.Data()) ;
 278   }
 279   if ( gSystem->AccessPathName(pararchivename) ) {
 280     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 281     gROOT->ProcessLine(processline.Data());
 282   }
 283
 284   TString ocwd = gSystem->WorkingDirectory();
 285   gSystem->ChangeDirectory(pararchivename);
 286
 287   // check for BUILD.sh and execute
 288   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 289     printf("*******************************\n");
 290     printf("*** Building PAR archive    ***\n");
 291     cout<<pararchivename<<endl;
 292     printf("*******************************\n");
 293
 294     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 295       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 296       return -1;
 297     }
 298   }
 299   // check for SETUP.C and execute
 300   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 301     printf("*******************************\n");
 302     printf("*** Setup PAR archive       ***\n");
 303     cout<<pararchivename<<endl;
 304     printf("*******************************\n");
 305     gROOT->Macro("PROOF-INF/SETUP.C");
 306   }
 307
 308   gSystem->ChangeDirectory(ocwd.Data());
 309   printf("Current dir: %s\n", ocwd.Data());
 310 }
 311
 312 void LoadLibraries(const anaModes mode) {
 313
 314   //--------------------------------------
 315   // Load the needed libraries most of them already loaded by aliroot
 316   //--------------------------------------
 317   gSystem->Load("libTree.so");
 318   gSystem->Load("libGeom.so");
 319   gSystem->Load("libVMC.so");
 320   gSystem->Load("libXMLIO.so");
 321   gSystem->Load("libPhysics.so");
 322
 323   //----------------------------------------------------------
 324   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 325   //----------------------------------------------------------
 326   if (mode==mLocal) {
 327     //--------------------------------------------------------
 328     // If you want to use already compiled libraries
 329     // in the aliroot distribution
 330     //--------------------------------------------------------
 331     gSystem->Load("libSTEERBase");
 332     gSystem->Load("libESD");
 333     gSystem->Load("libAOD");
 334     gSystem->Load("libANALYSIS");
 335     gSystem->Load("libANALYSISalice");
 336     gSystem->Load("libCORRFW.so");
 337     cerr<<"libCORRFW.so loaded..."<<endl;
 338     gSystem->Load("libPWG2flowCommon.so");
 339     cerr<<"libPWG2flowCommon.so loaded..."<<endl;
 340     gSystem->Load("libPWG2flowTasks.so");
 341     cerr<<"libPWG2flowTasks.so loaded..."<<endl;
 342   }
 343
 344   else if (mode == mLocalPAR) {
 345     //--------------------------------------------------------
 346     //If you want to use root and par files from aliroot
 347     //--------------------------------------------------------
 348      //If you want to use root and par files from aliroot
 349     //--------------------------------------------------------
 350     SetupPar("STEERBase");
 351     SetupPar("ESD");
 352     SetupPar("AOD");
 353     SetupPar("ANALYSIS");
 354     SetupPar("ANALYSISalice");
 355     SetupPar("PWG2AOD");
 356     SetupPar("CORRFW");
 357     SetupPar("PWG2flowCommon");
 358     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 359     SetupPar("PWG2flowTasks");
 360     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 361   }
 362
 363   //---------------------------------------------------------
 364   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 365   //---------------------------------------------------------
 366   else if (mode==mLocalSource) {
 367
 368     // In root inline compile
 369
 370
 371     // Constants
 372     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 373     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 374     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 375
 376     // Flow event
 377     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 378     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 379     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 380
 381     // Cuts
 382     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 383
 384     // Output histosgrams
 385     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 386     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 387     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 388     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 389
 390     // Functions needed for various methods
 391     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 392     gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
 393     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 394
 395     // Flow Analysis code for various methods
 396     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 397     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 398     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 399     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 400     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 401     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 402     gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");
 403
 404     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 405     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
 406
 407     cout << "finished loading macros!" << endl;
 408
 409   }
 410
 411 }
 412
 413