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