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