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