PWGCF/FLOW/macros/runFlowAnalysisOnTheFly.C

   1 // Settings for the simulation of events 'on the fly':
   2 //  a) Determine how many events you want to create;
   3 //  b) Set random or same seed for random generator;
   4 //  c) Determine multiplicites of events;
   5 //  d) Parametrize the phi distribution;
   6 //   d1) Enable/disable uniform event-wise fluctuations of v2;
   7 //   d2) Enable/diable pt dependence of v2;
   8 //  e) Parametrize the pt distribution;
   9 //  f) Determine how many times each sampled particle will be taken (simulating nonflow);
  10 //  g) Configure detector's:
  11 //   g1) acceptance;
  12 //   g2) efficiency;
  13 //  h) Decide which flow analysis methods you will use;
  14 //  i) Define simple cuts for Reference Particle (RP) selection;
  15 //  j) Define simple cuts for Particle of Interest (POI) selection;
  16 //  k) Define the ranges for two subevents separated with eta gap (needed only for SP method);
  17 //  l) Enable/disable usage of particle weights.
  18
  19 // a) Determine how many events you want to create:
  20 Int_t iNevts = 1000; // total statistics
  21
  22 // b) Set random or same seed for random generator:
  23 Bool_t bSameSeed = kFALSE; // if kTRUE, the created events are the same when re-doing flow analysis 'on the fly'
  24
  25 // c) Determine multiplicites of events:
  26 //    Remark 1: Multiplicity M for each event is sampled uniformly from interval iMinMult <= M < iMaxMult;
  27 //    Remark 2: For constant M of e.g. 500 for each event, set iMinMult = 500 and iMaxMult = 501.
  28 Int_t iMinMult = 500; // uniformly sampled multiplicity is >= iMinMult
  29 Int_t iMaxMult = 501; // uniformly sampled multiplicity is < iMaxMult
  30
  31 // d) Parametrize the phi distribution:
  32 //    Remark 1: Hardwired is Fourier-like distribution f(phi) = (1/2pi)(1+sum_{n=1}^{4} 2v_n cos[n(phi-rp)]),
  33 //              where reaction plane (rp) is sampled uniformly for each event from interval [0,2pi]
  34 Double_t dV1 = 0.0; // constant harmonic v1
  35 Double_t dV2 = 0.05; // constant harmonic v2
  36 Double_t dV3 = 0.0; // constant harmonic v3
  37 Double_t dV4 = 0.0; // constant harmonic v4
  38 Double_t dV5 = 0.0; // constant harmonic v5
  39 Double_t dV6 = 0.0; // constant harmonic v6
  40 //    Remark 2: By default all harmonics are constant for each event and for each particle. However, for v2
  41 //              the uniform event-wise fluctuations or pt dependence can be enabled:
  42 //  d1) Enable/disable uniform event-wise fluctuations of v2:
  43 Bool_t bUniformFluctuationsV2 = kFALSE; // enable uniform event-wise flow fluctuations (set than also dMinV2 and dMaxV2 bellow)
  44 Double_t dMinV2 = 0.04; // lower boundary on v2, when bUniformFluctuationsV2 = kTRUE
  45 Double_t dMaxV2 = 0.06; // upper boundary on v2, when bUniformFluctuationsV2 = kTRUE
  46 //  d2) Enable/disable pt dependence of v2:
  47 Bool_t bPtDependentV2 = kFALSE; // enable pt dependence of v2 (set then also dV2vsPtMax and dV2vsPtCutOff bellow)
  48 Double_t dV2vsPtCutOff = 2.0; // up to pt = dV2vsPtCutOff v2 is growing linearly as a function of pt
  49 Double_t dV2vsPtMax = 0.20; // for pt >= dV2vsPtCutOff, v2(pt) = dV2vsPtMax
  50
  51 // e) Parametrize the pt distribution:
  52 //    Remark: Hardwired is Boltzmann distribution f(pt) = pt*exp[-sqrt(dMass^2+pt^2)/dT]
  53 Double_t dMass = 0.13957; // mass in GeV/c^2 (e.g. m_{pions} = 0.13957)
  54 Double_t dTemperature = 0.44; // "temperature" in GeV/c (increase this parameter to get more high pt particles)
  55
  56 // f) Determine how many times each sampled particle will be taken in the analysis (simulating nonflow):
  57 Int_t nTimes = 1; // e.g. for nTimes = 2, strong 2-particle nonflow correlations are introduced
  58
  59 // g1) Configure detector's acceptance:
  60 Bool_t uniformAcceptance = kTRUE; // if kTRUE: detectors has uniform azimuthal acceptance.
  61                                   // if kFALSE: you will simulate detector with non-uniform acceptance in one or
  62                                   // two sectors. For each sector you specify phiMin, phiMax and probability p.
  63                                   // Then all particles emitted in direction phiMin < phi < phiMax will be taken
  64                                   // with probability p. If p = 0, that sector is completely blocked. Set bellow
  65                                   // phiMin1, phiMax1, p1 for the first sector and phiMin2, phiMax2, p2 for the second
  66                                   // sector. If you set phiMin2 = phiMax2 = p2 = 0, only first non-uniform sector is
  67                                   // simulated.
  68 // 1st non-uniform sector:
  69 Double_t phiMin1 = 60; // first non-uniform sector starts at this azimuth (in degrees)
  70 Double_t phiMax1 = 120; // first non-uniform sector ends at this azimuth (in degrees)
  71 Double_t p1 = 0.5; // probablitity that particles emitted in [phiMin1,phiMax1] are taken
  72 // 2nd non-uniform sector:
  73 Double_t phiMin2 = 0.; // first non-uniform sector starts at this azimuth (in degrees)
  74 Double_t phiMax2 = 0.; // first non-uniform sector ends at this azimuth (in degrees)
  75 Double_t p2 = 0.; // probablitity that particles emitted in [phiMin2,phiMax2] are taken
  76
  77 // g2) Configure detector's efficiency:
  78 Bool_t uniformEfficiency = kTRUE; // if kTRUE: detectors has uniform pT efficiency
  79                                   // if kFALSE: you will simulate detector with non-uniform pT efficiency.
  80                                   // Then all particles emitted in ptMin <= pt < ptMax will be taken
  81                                   // with probability p, to be specified in lines just below.
  82 Double_t ptMin = 0.8; // non-uniform efficiency vs pT starts at pT = fPtMin
  83 Double_t ptMax = 1.2; // non-uniform efficiency vs pT ends at pT = fPtMax
  84 Double_t p = 0.5; // probablitity that particles emitted in [ptMin,ptMax> are taken
  85
  86 // h) Decide which flow analysis methods you will use:
  87 Bool_t MCEP     = kTRUE; // Monte Carlo Event Plane
  88 Bool_t SP       = kTRUE; // Scalar Product (a.k.a 'flow analysis with eta gaps')
  89 Bool_t GFC      = kTRUE; // Generating Function Cumulants
  90 Bool_t QC       = kTRUE; // Q-cumulants
  91 Bool_t FQD      = kTRUE; // Fitted q-distribution
  92 Bool_t LYZ1SUM  = kTRUE; // Lee-Yang Zero (sum generating function), first pass over the data
  93 Bool_t LYZ1PROD = kTRUE; // Lee-Yang Zero (product generating function), first pass over the data
  94 Bool_t LYZ2SUM  = kFALSE; // Lee-Yang Zero (sum generating function), second pass over the data
  95 Bool_t LYZ2PROD = kFALSE; // Lee-Yang Zero (product generating function), second pass over the data
  96 Bool_t LYZEP    = kFALSE; // Lee-Yang Zero Event Plane
  97 Bool_t MH       = kFALSE; // Mixed Harmonics (used for strong parity violation studies)
  98 Bool_t NL       = kFALSE; // Nested Loops (neeed for debugging, only for developers)
  99 Bool_t MPC      = kTRUE; // Multi-particle correlations (NEW!)
 100
 101 // i) Define simple cuts for Reference Particle (RP) selection:
 102 Double_t ptMinRP = 0.0; // in GeV
 103 Double_t ptMaxRP = 10.0; // in GeV
 104 Double_t etaMinRP = -1.;
 105 Double_t etaMaxRP = 1.;
 106 Double_t phiMinRP = 0.0; // in degrees
 107 Double_t phiMaxRP = 360.0; // in degrees
 108 Bool_t bUseChargeRP = kFALSE; // if kFALSE, RPs with both sign of charges are taken
 109 Int_t chargeRP = 1; // +1 or -1
 110
 111 // j) Define simple cuts for Particle of Interest (POI) selection:
 112 Double_t ptMinPOI = 0.0; // in GeV
 113 Double_t ptMaxPOI = 10.0; // in GeV
 114 Double_t etaMinPOI = -1.; //
 115 Double_t etaMaxPOI = 1.;
 116 Double_t phiMinPOI = 0.0; // in degrees
 117 Double_t phiMaxPOI = 360.0; // in degrees
 118 Bool_t bUseChargePOI = kFALSE; // if kFALSE, POIs with both sign of charges are taken
 119 Int_t chargePOI = -1; // +1 or -1
 120
 121 // k) Define the ranges for two subevents separated with eta gap (needed only for SP method):
 122 Double_t etaMinA = -0.8; // minimum eta of subevent A
 123 Double_t etaMaxA = -0.5; // maximum eta of subevent A
 124 Double_t etaMinB = 0.5; // minimum eta of subevent B
 125 Double_t etaMaxB = 0.8; // maximum eta of subevent B
 126
 127 // l) Enable/disable usage of particle weights:
 128 Bool_t usePhiWeights = kFALSE; // phi weights
 129 Bool_t usePtWeights  = kFALSE; // pt weights
 130 Bool_t useEtaWeights = kFALSE; // eta weights
 131
 132 enum anaModes {mLocal,mLocalSource,mLocalPAR};
 133 // mLocal: Analyze data on your computer using aliroot
 134 // mLocalPAR: Analyze data on your computer using root + PAR files
 135 // mLocalSource: Analyze data on your computer using root + source files
 136
 137 #include "TStopwatch.h"
 138 #include "TObjArray"
 139 #include "Riostream.h"
 140 #include "TFile.h"
 141
 142 int runFlowAnalysisOnTheFly(Int_t mode=mLocal)
 143 {
 144  // Beging analysis 'on the fly'.
 145
 146  // a) Formal necessities....;
 147  // b) Initialize the flow event maker 'on the fly';
 148  // c) If enabled, access particle weights from external file;
 149  // d) Configure the flow analysis methods;
 150  // e) Simple cuts for RPs;
 151  // f) Simple cuts for POIs;
 152  // g) Create and analyse events 'on the fly';
 153  // h) Create the output file and directory structure for the final results of all methods;
 154  // i) Calculate and store the final results of all methods.
 155
 156  // a) Formal necessities....:
 157  CheckUserSettings();
 158  WelcomeMessage();
 159  TStopwatch timer;
 160  timer.Start();
 161  LoadLibraries(mode);
 162
 163  // b) Initialize the flow event maker 'on the fly':
 164  UInt_t uiSeed = 0; // if uiSeed is 0, the seed is determined uniquely in space and time via TUUID
 165  if(bSameSeed){uiSeed = 44;}
 166  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(uiSeed);
 167  eventMakerOnTheFly->SetMinMult(iMinMult);
 168  eventMakerOnTheFly->SetMaxMult(iMaxMult);
 169  eventMakerOnTheFly->SetMass(dMass);
 170  eventMakerOnTheFly->SetTemperature(dTemperature);
 171  eventMakerOnTheFly->SetV1(dV1);
 172  eventMakerOnTheFly->SetV2(dV2);
 173  eventMakerOnTheFly->SetV3(dV3);
 174  eventMakerOnTheFly->SetV4(dV4);
 175  eventMakerOnTheFly->SetV5(dV5);
 176  eventMakerOnTheFly->SetV6(dV6);
 177  if(bUniformFluctuationsV2)
 178  {
 179   eventMakerOnTheFly->SetUniformFluctuationsV2(bUniformFluctuationsV2);
 180   eventMakerOnTheFly->SetMinV2(dMinV2);
 181   eventMakerOnTheFly->SetMaxV2(dMaxV2);
 182  }
 183  if(bPtDependentV2)
 184  {
 185   eventMakerOnTheFly->SetPtDependentV2(bPtDependentV2);
 186   eventMakerOnTheFly->SetV2vsPtCutOff(dV2vsPtCutOff);
 187   eventMakerOnTheFly->SetV2vsPtMax(dV2vsPtMax);
 188  }
 189  eventMakerOnTheFly->SetSubeventEtaRange(etaMinA,etaMaxA,etaMinB,etaMaxB);
 190  eventMakerOnTheFly->SetNTimes(nTimes);
 191  if(!uniformAcceptance)
 192  {
 193   eventMakerOnTheFly->SetUniformAcceptance(kFALSE);
 194   eventMakerOnTheFly->SetFirstSectorPhiMin(phiMin1);
 195   eventMakerOnTheFly->SetFirstSectorPhiMax(phiMax1);
 196   eventMakerOnTheFly->SetFirstSectorProbability(p1);
 197   eventMakerOnTheFly->SetSecondSectorPhiMin(phiMin2);
 198   eventMakerOnTheFly->SetSecondSectorPhiMax(phiMax2);
 199   eventMakerOnTheFly->SetSecondSectorProbability(p2);
 200  }
 201  if(!uniformEfficiency)
 202  {
 203   eventMakerOnTheFly->SetUniformEfficiency(kFALSE);
 204   eventMakerOnTheFly->SetPtMin(ptMin);
 205   eventMakerOnTheFly->SetPtMax(ptMax);
 206   eventMakerOnTheFly->SetPtProbability(p);
 207  }
 208  eventMakerOnTheFly->Init();
 209
 210  // c) If enabled, access particle weights from external file:
 211  TFile *fileWithWeights = NULL;
 212  TList *listWithWeights = NULL;
 213  if(usePhiWeights||usePtWeights||useEtaWeights)
 214  {
 215   fileWithWeights = TFile::Open("weights.root","READ");
 216   if(fileWithWeights)
 217   {
 218    listWithWeights = (TList*)fileWithWeights->Get("weights");
 219   }
 220   else
 221   {
 222    cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 223    break;
 224   }
 225  } // end of if(usePhiWeights||usePtWeights||useEtaWeights)
 226
 227  // d) Configure the flow analysis methods:
 228  AliFlowAnalysisWithQCumulants *qc = NULL;
 229  AliFlowAnalysisWithCumulants *gfc = NULL;
 230  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
 231  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
 232  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
 233  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
 234  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
 235  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 236  AliFlowAnalysisWithScalarProduct *sp = NULL;
 237  AliFlowAnalysisWithMixedHarmonics *mh = NULL;
 238  AliFlowAnalysisWithNestedLoops *nl = NULL;
 239  AliFlowAnalysisWithMCEventPlane *mcep = NULL;
 240  AliFlowAnalysisWithMultiparticleCorrelations *mpc = NULL;
 241  // MCEP = monte carlo event plane
 242  if(MCEP)
 243  {
 244   //AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 245   mcep = new AliFlowAnalysisWithMCEventPlane();
 246   mcep->SetHarmonic(2); // default is v2
 247   mcep->Init();
 248  } // end of if(MCEP)
 249  // SP = Scalar Product
 250  if(SP)
 251  {
 252   sp = new AliFlowAnalysisWithScalarProduct();
 253   if(listWithWeights){sp->SetWeightsList(listWithWeights);}
 254   sp->SetUsePhiWeights(usePhiWeights);
 255   sp->SetHarmonic(2);
 256   sp->SetApplyCorrectionForNUA(kFALSE);
 257   sp->Init();
 258  } // end of if(SP)
 259  // QC = Q-cumulants
 260  if(QC)
 261  {
 262   qc = new AliFlowAnalysisWithQCumulants();
 263   if(listWithWeights){qc->SetWeightsList(listWithWeights);}
 264   if(usePhiWeights){qc->SetUsePhiWeights(usePhiWeights);}
 265   if(usePtWeights){qc->SetUsePtWeights(usePtWeights);}
 266   if(useEtaWeights){qc->SetUseEtaWeights(useEtaWeights);}
 267   qc->SetHarmonic(2);
 268   qc->SetCalculateDiffFlow(kTRUE);
 269   qc->SetCalculate2DDiffFlow(kFALSE); // vs (pt,eta)
 270   qc->SetApplyCorrectionForNUA(kFALSE);
 271   qc->SetFillMultipleControlHistograms(kFALSE);
 272   qc->SetMultiplicityWeight("combinations"); // default (other supported options are "unit" and "multiplicity")
 273   qc->SetCalculateCumulantsVsM(kFALSE);
 274   qc->SetCalculateAllCorrelationsVsM(kFALSE); // calculate all correlations in mixed harmonics "vs M"
 275   qc->SetnBinsMult(10000);
 276   qc->SetMinMult(0);
 277   qc->SetMaxMult(10000);
 278   qc->SetBookOnlyBasicCCH(kFALSE); // book only basic common control histograms
 279   qc->SetCalculateDiffFlowVsEta(kTRUE); // if you set kFALSE only differential flow vs pt is calculated
 280   qc->SetCalculateMixedHarmonics(kFALSE); // calculate all multi-partice mixed-harmonics correlators
 281   qc->Init();
 282  } // end of if(QC)
 283  // GFC = Generating Function Cumulants
 284  if(GFC)
 285  {
 286   gfc = new AliFlowAnalysisWithCumulants();
 287   if(listWithWeights){gfc->SetWeightsList(listWithWeights);}
 288   if(usePhiWeights){gfc->SetUsePhiWeights(usePhiWeights);}
 289   if(usePtWeights){gfc->SetUsePtWeights(usePtWeights);}
 290   if(useEtaWeights){gfc->SetUseEtaWeights(useEtaWeights);}
 291   // calculation vs multiplicity:
 292   gfc->SetCalculateVsMultiplicity(kFALSE);
 293   gfc->SetnBinsMult(10000);
 294   gfc->SetMinMult(0);
 295   gfc->SetMaxMult(10000);
 296   gfc->SetHarmonic(2);
 297   // tuning of interpolating parameters:
 298   gfc->SetTuneParameters(kFALSE);
 299   Double_t r0[10] = {1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7}; // up to 10 values allowed
 300   for(Int_t r=0;r<10;r++){gfc->SetTuningR0(r0[r],r);}
 301   gfc->Init();
 302  } // end of if(GFC)
 303  // FQD = Fitting q-distribution
 304  if(FQD)
 305  {
 306   fqd = new AliFlowAnalysisWithFittingQDistribution();
 307   if(listWithWeights){fqd->SetWeightsList(listWithWeights);}
 308   if(usePhiWeights){fqd->SetUsePhiWeights(usePhiWeights);}
 309   fqd->SetHarmonic(2);
 310   fqd->Init();
 311  } // end of if(FQD)
 312  // LYZ1 = Lee-Yang Zeroes first run
 313  if(LYZ1SUM)
 314  {
 315   lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 316   lyz1sum->SetFirstRun(kTRUE);
 317   lyz1sum->SetUseSum(kTRUE);
 318   lyz1sum->Init();
 319  } // end of if(LYZ1SUM)
 320  if(LYZ1PROD)
 321  {
 322   lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 323   lyz1prod->SetFirstRun(kTRUE);
 324   lyz1prod->SetUseSum(kFALSE);
 325   lyz1prod->Init();
 326  } // end of if(LYZ1PROD)
 327  // LYZ2 = Lee-Yang Zeroes second run
 328  if(LYZ2SUM)
 329  {
 330   lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 331   // read the input file from the first run
 332   TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 333   TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 334   if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie())
 335   {
 336    cerr<<" ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first."<<endl;
 337    break;
 338   } else
 339     {
 340      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 341      if(!inputListLYZ2SUM){cout<<"Input list LYZ2SUM is NULL pointer!"<<endl;break;}
 342      else
 343      {
 344       cout<<"LYZ2SUM input file/list read..."<<endl;
 345       lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 346       lyz2sum->SetFirstRun(kFALSE);
 347       lyz2sum->SetUseSum(kTRUE);
 348       lyz2sum->Init();
 349      }
 350    }
 351  } // end of if(LYZ2SUM)
 352  if(LYZ2PROD)
 353  {
 354   lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 355   // read the input file from the first run
 356   TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 357   TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 358   if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie())
 359   {
 360    cerr<<" ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first."<<endl;
 361    break;
 362   } else
 363     {
 364      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 365      if(!inputListLYZ2PROD){cout<<"Input list LYZ2PROD is NULL pointer!"<<endl;break;}
 366      else
 367      {
 368       cout<<"LYZ2PROD input file/list read..."<<endl;
 369       lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 370       lyz2prod->SetFirstRun(kFALSE);
 371       lyz2prod->SetUseSum(kFALSE);
 372       lyz2prod->Init();
 373      }
 374    }
 375  } // end of if(LYZ2PROD)
 376  // LYZEP = Lee-Yang Zeroes event plane
 377  if(LYZEP)
 378  {
 379   AliFlowLYZEventPlane *ep = new AliFlowLYZEventPlane() ;
 380   AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 381    // read the input file from the second lyz run
 382    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 383    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 384    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 385      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 386      break;
 387    }
 388    else {
 389      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 390      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 391      else {
 392        cout<<"LYZEP input file/list read..."<<endl;
 393        ep   ->SetSecondRunList(inputListLYZEP);
 394        lyzep->SetSecondRunList(inputListLYZEP);
 395        ep   ->Init();
 396        lyzep->Init();
 397      }
 398    }
 399  }
 400  // Mixed Harmonics:
 401  if(MH)
 402  {
 403   mh = new AliFlowAnalysisWithMixedHarmonics();
 404   mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
 405   mh->SetMinMultiplicity(100);
 406   mh->SetNoOfMultipicityBins(5);
 407   mh->SetMultipicityBinWidth(200);
 408   mh->Init();
 409  } // end of if(MH)
 410  // NL = Nested Loops:
 411  if(NL)
 412  {
 413   nl = new AliFlowAnalysisWithNestedLoops();
 414   nl->Init();
 415  } // end of if(NL)
 416  // MPC = Multi-particle correlations
 417  if(MPC)
 418  {
 419   mpc = new AliFlowAnalysisWithMultiparticleCorrelations();
 420   // Weights:
 421   //TFile *file = TFile::Open("weights.root","READ");
 422   //TH1D *weightsHist = (TH1D*) file->FindObjectAny("Phi weights");
 423   //mpc->SetPhiWeightsHist(weightsHist);
 424   // Control histograms:
 425   mpc->SetFillControlHistograms(kTRUE);
 426   mpc->SetFillKinematicsHist(kTRUE);
 427   // Q-vector:
 428   mpc->SetCalculateQvector(kTRUE);
 429   // Correlations:
 430   mpc->SetCalculateCorrelations(kTRUE);
 431   mpc->SetSkipZeroHarmonics(kTRUE);
 432   mpc->SetCalculateIsotropic(kTRUE);
 433   // Standard candles:
 434   mpc->SetCalculateStandardCandles(kTRUE);
 435   // Init:
 436   mpc->Init();
 437  } // end of if(MPC)
 438
 439  // e) Simple cuts for RPs:
 440  AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
 441  cutsRP->SetPtMax(ptMaxRP);
 442  cutsRP->SetPtMin(ptMinRP);
 443  cutsRP->SetEtaMax(etaMaxRP);
 444  cutsRP->SetEtaMin(etaMinRP);
 445  cutsRP->SetPhiMax(phiMaxRP*TMath::Pi()/180.);
 446  cutsRP->SetPhiMin(phiMinRP*TMath::Pi()/180.);
 447  if(bUseChargeRP){cutsRP->SetCharge(chargeRP);}
 448
 449  // f) Simple cuts for POIs:
 450  AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
 451  cutsPOI->SetPtMax(ptMaxPOI);
 452  cutsPOI->SetPtMin(ptMinPOI);
 453  cutsPOI->SetEtaMax(etaMaxPOI);
 454  cutsPOI->SetEtaMin(etaMinPOI);
 455  cutsPOI->SetPhiMax(phiMaxPOI*TMath::Pi()/180.);
 456  cutsPOI->SetPhiMin(phiMinPOI*TMath::Pi()/180.);
 457  if(bUseChargePOI){cutsPOI->SetCharge(chargePOI);}
 458
 459  // g) Create and analyse events 'on the fly':
 460  for(Int_t i=0;i<iNevts;i++)
 461  {
 462   // Creating the event 'on the fly':
 463   AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly(cutsRP,cutsPOI);
 464   // Passing the created event to flow analysis methods:
 465   if(MCEP){mcep->Make(event);}
 466   if(QC){qc->Make(event);}
 467   if(GFC){gfc->Make(event);}
 468   if(FQD){fqd->Make(event);}
 469   if(LYZ1SUM){lyz1sum->Make(event);}
 470   if(LYZ1PROD){lyz1prod->Make(event);}
 471   if(LYZ2SUM){lyz2sum->Make(event);}
 472   if(LYZ2PROD){lyz2prod->Make(event);}
 473   if(LYZEP){lyzep->Make(event,ep);}
 474   if(SP){sp->Make(event);}
 475   if(MH){mh->Make(event);}
 476   if(NL){nl->Make(event);}
 477   if(MPC){mpc->Make(event);}
 478   delete event;
 479  } // end of for(Int_t i=0;i<iNevts;i++)
 480
 481  // h) Create the output file and directory structure for the final results of all methods:
 482  TString outputFileName = "AnalysisResults.root";
 483  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 484  const Int_t nMethods = 13;
 485  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MPC"};
 486  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 487  TString fileName[nMethods];
 488  for(Int_t i=0;i<nMethods;i++)
 489  {
 490   fileName[i]+="output";
 491   fileName[i]+=method[i].Data();
 492   fileName[i]+="analysis";
 493   dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 494  }
 495
 496  // i) Calculate and store the final results of all methods:
 497  if(MCEP){mcep->Finish();mcep->WriteHistograms(dirFileFinal[0]);}
 498  if(SP){sp->Finish();sp->WriteHistograms(dirFileFinal[1]);}
 499  if(GFC){gfc->Finish();gfc->WriteHistograms(dirFileFinal[2]);}
 500  if(QC){qc->Finish();qc->WriteHistograms(dirFileFinal[3]);}
 501  if(FQD){fqd->Finish();fqd->WriteHistograms(dirFileFinal[4]);}
 502  if(LYZ1SUM){lyz1sum->Finish();lyz1sum->WriteHistograms(dirFileFinal[5]);}
 503  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
 504  if(LYZ2SUM){lyz2sum->Finish();lyz2sum->WriteHistograms(dirFileFinal[7]);}
 505  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
 506  if(LYZEP){lyzep->Finish();lyzep->WriteHistograms(dirFileFinal[9]);}
 507  if(MH){mh->Finish();mh->WriteHistograms(dirFileFinal[10]);}
 508  if(NL){nl->Finish();nl->WriteHistograms(dirFileFinal[11]);}
 509  if(MPC){mpc->Finish();mpc->WriteHistograms(dirFileFinal[12]);}
 510
 511  outputFile->Close();
 512  delete outputFile;
 513
 514  cout<<endl;
 515  cout<<endl;
 516  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 517  cout<<endl;
 518
 519  timer.Stop();
 520  cout << endl;
 521  timer.Print();
 522
 523 } // end of int runFlowAnalysisOnTheFly(Int_t mode=mLocal)
 524
 525 void SetupPar(char* pararchivename)
 526 {
 527   //Load par files, create analysis libraries
 528   //For testing, if par file already decompressed and modified
 529   //classes then do not decompress.
 530
 531   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 532   TString parpar(Form("%s.par", pararchivename)) ;
 533   if ( gSystem->AccessPathName(parpar.Data()) ) {
 534     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 535     TString processline(Form(".! make %s", parpar.Data())) ;
 536     gROOT->ProcessLine(processline.Data()) ;
 537     gSystem->ChangeDirectory(cdir) ;
 538     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 539     gROOT->ProcessLine(processline.Data()) ;
 540   }
 541   if ( gSystem->AccessPathName(pararchivename) ) {
 542     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 543     gROOT->ProcessLine(processline.Data());
 544   }
 545
 546   TString ocwd = gSystem->WorkingDirectory();
 547   gSystem->ChangeDirectory(pararchivename);
 548
 549   // check for BUILD.sh and execute
 550   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 551     printf("*******************************\n");
 552     printf("*** Building PAR archive    ***\n");
 553     cout<<pararchivename<<endl;
 554     printf("*******************************\n");
 555
 556     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 557       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 558       return -1;
 559     }
 560   }
 561   // check for SETUP.C and execute
 562   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 563     printf("*******************************\n");
 564     printf("*** Setup PAR archive       ***\n");
 565     cout<<pararchivename<<endl;
 566     printf("*******************************\n");
 567     gROOT->Macro("PROOF-INF/SETUP.C");
 568   }
 569
 570   gSystem->ChangeDirectory(ocwd.Data());
 571   printf("Current dir: %s\n", ocwd.Data());
 572 }
 573
 574 void CheckUserSettings()
 575 {
 576  // Check if user settings make sense before taking off.
 577
 578  if(iNevts <= 0)
 579  {
 580   printf("\n WARNING: nEvts <= 0 !!!! Please check your settings before taking off.\n\n");
 581   exit(0);
 582  }
 583  if(iMinMult < 0.)
 584  {
 585   printf("\n WARNING: iMinMult < 0 !!!! Please check your settings before taking off.\n\n");
 586   exit(0);
 587  }
 588  if(iMaxMult <= 0.)
 589  {
 590   printf("\n WARNING: iMaxMult <= 0 !!!! Please check your settings before taking off.\n\n");
 591   exit(0);
 592  }
 593  if(iMinMult >= iMaxMult)
 594  {
 595   printf("\n WARNING: iMinMult >= iMaxMult !!!! Please check your settings before taking off.\n\n");
 596   exit(0);
 597  }
 598  if(dMass < 0.)
 599  {
 600   printf("\n WARNING: dMass < 0 !!!! Please check your settings before taking off.\n\n");
 601   exit(0);
 602  }
 603  if(dTemperature <= 1e-44)
 604  {
 605   printf("\n WARNING: dTemperature <= 0 !!!! Please check your settings before taking off.\n\n");
 606   exit(0);
 607  }
 608  if(TMath::Abs(dV1) > 0.5)
 609  {
 610   printf("\n WARNING: |dV1| > 0.5 !!!! Please check your settings before taking off.\n\n");
 611   exit(0);
 612  }
 613  if(TMath::Abs(dV2) > 0.5)
 614  {
 615   printf("\n WARNING: |dV2| > 0.5 !!!! Please check your settings before taking off.\n\n");
 616   exit(0);
 617  }
 618  if(TMath::Abs(dV3) > 0.5)
 619  {
 620   printf("\n WARNING: |dV3| > 0.5 !!!! Please check your settings before taking off.\n\n");
 621   exit(0);
 622  }
 623  if(TMath::Abs(dV4) > 0.5)
 624  {
 625   printf("\n WARNING: |dV4| > 0.5 !!!! Please check your settings before taking off.\n\n");
 626   exit(0);
 627  }
 628  if(TMath::Abs(dV5) > 0.5)
 629  {
 630   printf("\n WARNING: |dV5| > 0.5 !!!! Please check your settings before taking off.\n\n");
 631   exit(0);
 632  }
 633  if(TMath::Abs(dV6) > 0.5)
 634  {
 635   printf("\n WARNING: |dV6| > 0.5 !!!! Please check your settings before taking off.\n\n");
 636   exit(0);
 637  }
 638  if(LYZ1SUM && LYZ2SUM)
 639  {
 640   cout<<" WARNING: You cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl;
 641   exit(0);
 642  }
 643  if(LYZ1PROD && LYZ2PROD)
 644  {
 645   cout<<" WARNING: You cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl;
 646   exit(0);
 647  }
 648  if(LYZ2SUM && LYZEP)
 649  {
 650   cout<<" WARNING: You cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl;
 651   exit(0);
 652  }
 653  if(LYZ1SUM && LYZEP)
 654  {
 655   cout<<" WARNING: You cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl;
 656   exit(0);
 657  }
 658
 659  if(!uniformAcceptance && phiMin1 > phiMax1)
 660  {
 661   cout<<" WARNING: You must have phiMin1 < phiMax1 !!!!"<<endl;
 662   exit(0);
 663  }
 664  if(!uniformAcceptance && !((TMath::Abs(phiMin2) < 1.e-44) && (TMath::Abs(phiMax2) < 1.e-44) && (TMath::Abs(p2) < 1.e-44))
 665     && (phiMin2 < phiMax1 || phiMin2 > phiMax2))
 666  {
 667   cout<<" WARNING: You must have phiMin2 > phiMax1 and phiMin2 < phiMax2 !!!!"<<endl;
 668   exit(0);
 669  }
 670  if((phiMin1 < 0 || phiMin1 > 360) || (phiMax1 < 0 || phiMax1 > 360) ||
 671     (phiMin2 < 0 || phiMin2 > 360) || (phiMax2 < 0 || phiMax2 > 360) )
 672  {
 673   cout<<" WARNING: You must take azimuthal angles from interval [0,360] !!!!"<<endl;
 674   exit(0);
 675  }
 676  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
 677  {
 678   cout<<" WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
 679   exit(0);
 680  }
 681  if(bPtDependentV2 && bUniformFluctuationsV2)
 682  {
 683   cout<<" WARNING: Uniform fluctuations not supported for pt denependent v2 !!!!"<<endl;
 684   exit(0);
 685  }
 686
 687 } // end of void CheckUserSettings()
 688
 689 void WelcomeMessage()
 690 {
 691  // Welcome.
 692
 693  cout<<endl;
 694  cout<<endl;
 695  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
 696  cout<<endl;
 697
 698  gSystem->Sleep(1544);
 699
 700  cout<<endl;
 701  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
 702  cout<<endl;
 703  cout<<endl;
 704
 705  gSystem->Sleep(1544);
 706
 707 } // end of void WelcomeMessage()
 708
 709 void LoadLibraries(const anaModes mode) {
 710
 711   //--------------------------------------
 712   // Load the needed libraries most of them already loaded by aliroot
 713   //--------------------------------------
 714   //gSystem->Load("libTree");
 715   gSystem->Load("libGeom");
 716   gSystem->Load("libVMC");
 717   gSystem->Load("libXMLIO");
 718   gSystem->Load("libPhysics");
 719
 720   //----------------------------------------------------------
 721   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 722   //----------------------------------------------------------
 723   if (mode==mLocal) {
 724     //--------------------------------------------------------
 725     // If you want to use already compiled libraries
 726     // in the aliroot distribution
 727     //--------------------------------------------------------
 728     gSystem->Load("libSTEERBase");
 729     gSystem->Load("libESD");
 730     gSystem->Load("libAOD");
 731     gSystem->Load("libANALYSIS");
 732     gSystem->Load("libANALYSISalice");
 733     gSystem->Load("libCORRFW");
 734     cerr<<"libCORRFW loaded..."<<endl;
 735     gSystem->Load("libPWGflowBase");
 736     cerr<<"libPWGflowBase loaded..."<<endl;
 737     gSystem->Load("libPWGflowTasks");
 738     cerr<<"libPWGflowTasks loaded..."<<endl;
 739   }
 740
 741   else if (mode == mLocalPAR) {
 742     //--------------------------------------------------------
 743     //If you want to use root and par files from aliroot
 744     //--------------------------------------------------------
 745      //If you want to use root and par files from aliroot
 746     //--------------------------------------------------------
 747     SetupPar("STEERBase");
 748     SetupPar("ESD");
 749     SetupPar("AOD");
 750     SetupPar("ANALYSIS");
 751     SetupPar("ANALYSISalice");
 752     SetupPar("PWG2AOD");
 753     SetupPar("CORRFW");
 754     SetupPar("PWGflowBase");
 755     cerr<<"PWGflowBase.par loaded..."<<endl;
 756     SetupPar("PWGflowTasks");
 757     cerr<<"PWGflowTasks.par loaded..."<<endl;
 758   }
 759
 760   //---------------------------------------------------------
 761   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 762   //---------------------------------------------------------
 763   else if (mode==mLocalSource) {
 764
 765     // In root inline compile
 766
 767     // Constants
 768     gROOT->LoadMacro("Base/AliFlowCommonConstants.cxx+");
 769     gROOT->LoadMacro("Base/AliFlowLYZConstants.cxx+");
 770
 771     // Flow event
 772     gROOT->LoadMacro("Base/AliFlowVector.cxx+");
 773     gROOT->LoadMacro("Base/AliFlowTrackSimple.cxx+");
 774     gROOT->LoadMacro("Base/AliFlowTrackSimpleCuts.cxx+");
 775     gROOT->LoadMacro("Base/AliFlowEventSimple.cxx+");
 776
 777     // Output histosgrams
 778     gROOT->LoadMacro("Base/AliFlowCommonHist.cxx+");
 779     gROOT->LoadMacro("Base/AliFlowCommonHistResults.cxx+");
 780     gROOT->LoadMacro("Base/AliFlowLYZHist1.cxx+");
 781     gROOT->LoadMacro("Base/AliFlowLYZHist2.cxx+");
 782
 783     // Functions needed for various methods
 784     gROOT->LoadMacro("Base/AliCumulantsFunctions.cxx+");
 785     gROOT->LoadMacro("Base/AliFlowLYZEventPlane.cxx+");
 786
 787     // Flow Analysis code for various methods
 788     gROOT->LoadMacro("Base/AliFlowAnalysisWithMCEventPlane.cxx+");
 789     gROOT->LoadMacro("Base/AliFlowAnalysisWithScalarProduct.cxx+");
 790     gROOT->LoadMacro("Base/AliFlowAnalysisWithLYZEventPlane.cxx+");
 791     gROOT->LoadMacro("Base/AliFlowAnalysisWithLeeYangZeros.cxx+");
 792     gROOT->LoadMacro("Base/AliFlowAnalysisWithCumulants.cxx+");
 793     gROOT->LoadMacro("Base/AliFlowAnalysisWithQCumulants.cxx+");
 794     gROOT->LoadMacro("Base/AliFlowAnalysisWithFittingQDistribution.cxx+");
 795     gROOT->LoadMacro("Base/AliFlowAnalysisWithMixedHarmonics.cxx+");
 796     gROOT->LoadMacro("Base/AliFlowAnalysisWithNestedLoops.cxx+");
 797
 798     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 799     gROOT->LoadMacro("Base/AliFlowEventSimpleMakerOnTheFly.cxx+");
 800
 801     cout << "finished loading macros!" << endl;
 802
 803   }
 804
 805 }
 806
 807