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 MCEP     = kTRUE;
  10 Bool_t SP       = kTRUE;
  11 Bool_t GFC      = kTRUE;
  12 Bool_t QC       = kTRUE;
  13 Bool_t FQD      = kTRUE;
  14 Bool_t LYZ1SUM  = kTRUE;
  15 Bool_t LYZ1PROD = kTRUE;
  16 Bool_t LYZ2SUM  = kFALSE;
  17 Bool_t LYZ2PROD = kFALSE;
  18 Bool_t LYZEP    = kFALSE;
  19 Bool_t MH       = kTRUE; // mixed harmonics
  20 Bool_t NL       = kFALSE; // nested loops
  21 Bool_t MCEP_AH  = kFALSE; // MCEP in another harmonic
  22 //--------------------------------------------------------------------------------------
  23
  24 // Weights
  25 // use weights for Q-vector:
  26 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
  27 Bool_t usePtWeights  = kFALSE; // pt weights
  28 Bool_t useEtaWeights = kFALSE; // eta weights
  29
  30 // Define the range for eta subevents
  31 Double_t minA = -0.9;
  32 Double_t maxA = -0.01;
  33 Double_t minB = 0.01;
  34 Double_t maxB = 0.9;
  35
  36 // Define simple cuts for RP selection:
  37 Double_t ptMinRP = 0.0; // in GeV
  38 Double_t ptMaxRP = 10.0; // in GeV
  39 Double_t etaMinRP  = -1.;
  40 Double_t etaMaxRP  = 1.;
  41 Double_t phiMinRP  = 0.0; // in degrees
  42 Double_t phiMaxRP  = 360.0; // in degrees
  43 Int_t pidRP = -1; // to be improved (supported eventually)
  44
  45 // Define simple cuts for POI selection:
  46 Double_t ptMinPOI = 0.0; // in GeV
  47 Double_t ptMaxPOI = 10.0; // in GeV
  48 Double_t etaMinPOI  = -1.;
  49 Double_t etaMaxPOI  = 1.;
  50 Double_t phiMinPOI  = 0.0; // in degrees
  51 Double_t phiMaxPOI  = 360.0; // in degrees
  52 Int_t pidPOI = -1; // to be improved (supported eventually)
  53
  54 // Parameters for the simulation of events 'on the fly':
  55 //===SEED========================================================
  56 // use always the same seed for random generators.
  57 // usage of same seed (kTRUE) is relevant in two cases:
  58 // 1.) If you want to use LYZ method to calcualte differential flow;
  59 // 2.) If you want to use phi weights for GFC, QC and FQD
  60 Bool_t bSameSeed = kFALSE;
  61
  62 //===NONFLOW=============================================
  63 // number of times to use each track (to simulate nonflow)
  64 Int_t iLoops = 1;
  65 Double_t phiRange = 0.0; // If the original track with azimuthal angle phi is splitted,
  66                          // the azimuthal angle of splitted track is sampled uniformly
  67                          // from (phi-phiRange, phi+phiRange). If phiRange = 0.0, the
  68                          // azimuthal angle of splitted track is the same as azimuthal
  69                          // angle of original track. (Remark: phiRange's unit is degree.)
  70 Double_t ptRange = 0.0; // If the original track with momentum pt is splitted,
  71                         // the momentum of splitted track is sampled uniformly
  72                         // from (pt-ptRange, pt+ptRange). If ptRange = 0.0, the
  73                         // momentum of splitted track is the same as momentum
  74                         // of original track. (Remark: ptRange's unit is GeV.)
  75 Double_t etaRange = 0.0; // If the original track with pseudorapidity eta is splitted,
  76                          // the pseudorapidity of splitted track is sampled uniformly
  77                          // from (eta-etaRange, eta+etaRange). If etaRange = 0.0, the
  78                          // pseudorapidity of splitted track will be the same as the
  79                          // pseudorapidity of original track.
  80 // in addition one can simulate nonflow only in a certain detector's sector
  81 // ranging from nonflowSectorMin to nonflowSectorMax. Outside this sector the
  82 // tracks will not be splitted. Use the following two settings only if iLoops>1:
  83 Double_t nonflowSectorMin = 0.0; // detector's sector in which tracks will be splitted starts at this angle (in degrees)
  84 Double_t nonflowSectorMax = 360.0; // detector's sector in which tracks will be splitted ends at this angle (in degrees)
  85 // to study nonflow fluctuations event-by-event one can create the jet-like structure in the following way:
  86 Bool_t bCreateJets = kFALSE; // set kTRUE if you want to study nonflow fluctuations via jet creation
  87 Double_t jetProbability = 0.2; // probability than event will have a jet
  88 Double_t jetTracksFraction = 0.1; // percentage of total tracks in an event belonging to the jet
  89 Double_t jetCone = 20; // angular size of jet cone (in degrees)
  90
  91 //===GLAUBER MODEL================================================================
  92 Bool_t bUseGlauberModel = kFALSE; // 1.) if kTRUE = multiplicity and constant flow harmonics are
  93                                   //                determined event-by-event from Glauber model
  94                                   //                (pt and/or eta dependence of flow harmonics not supported;
  95                                   //                settings for multiplicity and flow harmonics bellow are irrelevant)
  96                                   // 2.) if kFALSE = multiplicity and flow harmonics are determined
  97                                   //                 independently event-by-event with bellow settings
  98                                   //                 (pt and/or eta dependence of flow harmonics supported)
  99
 100 //===FLOW HARMONICS===============================================================
 101 // harmonics V1, V2, V4... are constants or functions of pt and eta:
 102 Bool_t bPtDependentHarmonicV1 = kFALSE;
 103 Bool_t bEtaDependentHarmonicV1 = kFALSE;
 104 Bool_t bPtDependentHarmonicV2 = kFALSE;
 105 Bool_t bEtaDependentHarmonicV2 = kFALSE;
 106 Bool_t bPtDependentHarmonicV4 = kFALSE;
 107 Bool_t bEtaDependentHarmonicV4 = kFALSE;
 108 // 1.) if you use constant harmonics (bPtDependentHarmonic* = kFALSE, bEtaDependentHarmonic* = kFALSE)
 109 //     you can additionally select if V2 will be sampled e-b-e from Gaussian or from uniform distribution
 110 //     (constant harmonics V1 and V4 will be always sampled from Gaussian):
 111 Bool_t bConstantV2IsSampledFromGauss = kTRUE;
 112  //  1a.) if kTRUE = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
 113  //                  mean = dV2RP and spread = dV2SpreadRP (analogously for V1 and V4).
 114  //  1b.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from
 115  //                   interval [dMinV2RP,dMaxV2RP] (not supported for V1 and V4).
 116  //  1c.) for a fixed elliptic flow e-b-e use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP
 117  // V2:
 118  Double_t dV2RP = 0.05;       // elliptic flow of RPs (if sampled from Gaussian)
 119  Double_t dV2SpreadRP = 0.0;  // elliptic flow spread of RPs (if sampled from Gaussian)
 120  Double_t dMinV2RP = 0.04;    // minimal elliptic flow of RPs (if sampled uniformly)
 121  Double_t dMaxV2RP = 0.06;    // maximal elliptic flow of RPs (if sampled uniformly)
 122  // V1:
 123  Double_t dV1RP = 0.0; // directed flow of RPs
 124  Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
 125  // V4:
 126  Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
 127  Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
 128 // 2.) if you use (pt,eta) dependent harmonic V1 (bPtDependentHarmonicV1 = kTRUE or bEtaDependentHarmonicV1 = kTRUE):
 129 //  2a.) V1(pt) is linear up to pt = dV1PtCutOff and for pt > dV1PtCutOff it is constant, V1(pt) = dV1vsPtEtaMax:
 130 //  2b.) V1(eta) is determined from formula V1(eta) = -eta
 131  Double_t dV1vsPtEtaMax = 0.10; // maximum value of V1(pt,eta) (for pt >= dV1PtCutOff and at eta = 0)
 132  Double_t dV1PtCutOff = 2.0; // up to this pt value V1(pt) is growing linearly versus pt
 133 // 3.) if you use (pt,eta) dependent harmonic V2 (bPtDependentHarmonicV2 = kTRUE or bEtaDependentHarmonicV2 = kTRUE):
 134 //  2a.) V2(pt) is linear up to pt = dV2PtCutOff and for pt > dV2PtCutOff it is constant, V2(pt) = dV2vsPtEtaMax:
 135 //  2b.) V2(eta) is Gaussian centered at midrapidity (eta=0) with V2(0) = dV2vsPtEtaMax and spread = dV2vsEtaSpread:
 136  Double_t dV2vsPtEtaMax = 0.20; // maximum value of V2(pt,eta) (for pt >= dV2PtCutOff and at eta = 0)
 137  Double_t dV2PtCutOff = 2.0; // up to this pt value V2(pt) is growing linearly versus pt
 138  Double_t dV2vsEtaSpread = 0.75; // V2(eta) is Gaussian centered at midrapidity (eta=0) with spread = dV2vsEtaSpread
 139 // 4.) if you use (pt,eta) dependent harmonic V4 (bPtDependentHarmonicV4 = kTRUE or bEtaDependentHarmonicV4 = kTRUE):
 140 //     V4(pt,eta) is determined as V4(pt,eta) = pow(V2(pt,eta),2.)
 141
 142 //===MULTIPLICITY===============================================================
 143 // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
 144 //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
 145 // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from
 146 //                 interval [iMinMultOfRP,iMaxMultOfRP]
 147 // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
 148 Bool_t bMultDistrOfRPsIsGauss = kTRUE;
 149
 150 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
 151 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
 152 Int_t iMinMultOfRP = 400;              // minimal multiplicity of RPs(if sampled uniformly)
 153 Int_t iMaxMultOfRP = 600;             // maximal multiplicity of RPs (if sampled uniformly)
 154
 155 //===DETECTOR ACCEPTANCE===============================================================
 156
 157 // 1.) if kTRUE = detectors has uniform azimuthal acceptance
 158 // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors.
 159 //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles
 160 //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
 161 //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2
 162 //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is
 163 //                 simulated.
 164 Bool_t uniformAcceptance = kTRUE;
 165
 166 // settings for non-uniform acceptance:
 167 // Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
 168
 169 // 1st non-uniform sector:
 170 Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
 171 Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
 172 Double_t p1 = 0.5;     // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
 173                         // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
 174
 175 // 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):
 176 Double_t phimin2 = 0.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
 177 Double_t phimax2 = 0.0; // second non-uniform sector ends at this azimuth
 178 Double_t p2 = 0.0;
 179
 180
 181 //===MODIFYING Pt SPECTRA===============================================================
 182 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs)
 183
 184 enum anaModes {mLocal,mLocalSource,mLocalPAR};
 185 // mLocal: Analyze data on your computer using aliroot
 186 // mLocalPAR: Analyze data on your computer using root + PAR files
 187 // mLocalSource: Analyze data on your computer using root + source files
 188
 189 int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
 190 {
 191  CheckUserSettings();
 192  TStopwatch timer;
 193  timer.Start();
 194
 195  cout<<endl;
 196  cout<<endl;
 197  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
 198  cout<<endl;
 199
 200  cout<<endl;
 201  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
 202  cout<<endl;
 203  cout<<endl;
 204
 205  LoadLibraries(mode);
 206
 207  // Initialize the seed for random generator
 208  UInt_t sseed = 0;
 209
 210  if(bSameSeed)
 211  {
 212   sseed = 44; // the default constant value for seed for random generators
 213  }
 214
 215  if(!bSameSeed)
 216  {
 217   TTimeStamp dt;
 218   sseed = dt.GetNanoSec()/1000;
 219  }
 220
 221  cout<<endl;
 222  cout<<"Seed for the random generators is "<<sseed<<endl;
 223  cout<<endl;
 224
 225  //---------------------------------------------------------------------------------------
 226  // If the weights are used:
 227  TFile *fileWithWeights = NULL;
 228  TList *listWithWeights = NULL;
 229
 230  if(usePhiWeights||usePtWeights||useEtaWeights) {
 231    fileWithWeights = TFile::Open("weights.root","READ");
 232    if(fileWithWeights) {
 233      listWithWeights = (TList*)fileWithWeights->Get("weights");
 234    }
 235    else
 236      {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 237       break;
 238      }
 239  }
 240
 241  //---------------------------------------------------------------------------------------
 242  // Initialize the flowevent maker
 243  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
 244  eventMakerOnTheFly->Init();
 245  //set the range for eta subevents
 246  eventMakerOnTheFly->SetSubeventEtaRange(minA,maxA,minB,maxB);
 247
 248  //---------------------------------------------------------------------------------------
 249  // Initialize all the flow methods for default analysis:
 250  AliFlowAnalysisWithQCumulants *qc = NULL;
 251  AliFlowAnalysisWithCumulants *gfc = NULL;
 252  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
 253  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
 254  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
 255  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
 256  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
 257  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 258  AliFlowAnalysisWithScalarProduct *sp = NULL;
 259  AliFlowAnalysisWithMixedHarmonics *mh = NULL;
 260  AliFlowAnalysisWithNestedLoops *nl = NULL;
 261  AliFlowAnalysisWithMCEventPlane *mcep = NULL;
 262  AliFlowAnalysisWithMCEventPlane *mcep_AH = NULL;
 263
 264  // MCEP = monte carlo event plane
 265  if (MCEP) {
 266    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 267    //mcep->SetHarmonic(2); // default is v2
 268    //mcep->SetFlowOfResonances(kTRUE);
 269    mcep->Init();
 270  }
 271
 272  // MCEP = monte carlo event plane in another harmonic:
 273  if(MCEP_AH)
 274  {
 275   AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
 276   mcep_ah->SetHarmonic(1);
 277   mcep_ah->Init();
 278  }
 279
 280  // Mixed harmonics:
 281  if(MH)
 282  {
 283   AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
 284   mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
 285   mh->SetMinMultiplicity(100);
 286   mh->SetNoOfMultipicityBins(5);
 287   mh->SetMultipicityBinWidth(200);
 288   mh->Init();
 289  }
 290
 291  // NL = nested loops:
 292  if(NL) {
 293    AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
 294    nl->Init();
 295  }
 296
 297  // QC = Q-cumulants
 298  if(QC) {
 299    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 300    if(listWithWeights) qc->SetWeightsList(listWithWeights);
 301    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 302    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 303    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 304    // qc->SetHarmonic(2); // default is v2
 305    // qc->SetApplyCorrectionForNUA(kTRUE); // default
 306    // qc->SetCalculate2DFlow(kFALSE); // default
 307    // qc->SetMultiplicityWeight("combinations"); // default
 308    // qc->SetMultiplicityWeight("unit");
 309    // qc->SetMultiplicityWeight("multiplicity");
 310    qc->SetnBinsMult(10000);
 311    qc->SetMinMult(0);
 312    qc->SetMaxMult(10000);
 313    qc->Init();
 314  }
 315
 316  // GFC = Generating Function Cumulants
 317  if(GFC) {
 318    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 319    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 320    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 321    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 322    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 323    // calculation vs multiplicity:
 324    gfc->SetCalculateVsMultiplicity(kFALSE);
 325    gfc->SetnBinsMult(10000);
 326    gfc->SetMinMult(0);
 327    gfc->SetMaxMult(10000);
 328    // tuning of interpolating parameters:
 329    gfc->SetTuneParameters(kFALSE);
 330    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
 331    for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
 332    gfc->Init();
 333  }
 334
 335  // FQD = Fitting q-distribution
 336  if(FQD) {
 337    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 338    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 339    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 340    fqd->Init();
 341  }
 342
 343  // SP = Scalar Product
 344  if(SP) {
 345    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 346    if(listWithWeights) sp->SetWeightsList(listWithWeights);
 347    sp->SetUsePhiWeights(usePhiWeights);
 348    sp->Init();
 349  }
 350
 351  // LYZ1 = Lee-Yang Zeroes first run
 352  if(LYZ1SUM) {
 353    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 354    lyz1sum->SetFirstRun(kTRUE);
 355    lyz1sum->SetUseSum(kTRUE);
 356    lyz1sum->Init();
 357  }
 358  if(LYZ1PROD) {
 359    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 360    lyz1prod->SetFirstRun(kTRUE);
 361    lyz1prod->SetUseSum(kFALSE);
 362    lyz1prod->Init();
 363  }
 364  // LYZ2 = Lee-Yang Zeroes second run
 365  if(LYZ2SUM) {
 366    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 367    // read the input file from the first run
 368    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 369    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 370    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 371      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 372      break;
 373    }
 374    else {
 375      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 376      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
 377      else {
 378        cout<<"LYZ2SUM input file/list read..."<<endl;
 379        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 380        lyz2sum->SetFirstRun(kFALSE);
 381        lyz2sum->SetUseSum(kTRUE);
 382        lyz2sum->Init();
 383      }
 384    }
 385  }
 386  if(LYZ2PROD) {
 387    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 388    // read the input file from the first run
 389    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 390    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 391    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 392      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 393      break;
 394    }
 395    else {
 396      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 397      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
 398      else {
 399        cout<<"LYZ2PROD input file/list read..."<<endl;
 400        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 401        lyz2prod->SetFirstRun(kFALSE);
 402        lyz2prod->SetUseSum(kFALSE);
 403        lyz2prod->Init();
 404      }
 405    }
 406  }
 407
 408  // LYZEP = Lee-Yang Zeroes event plane
 409  if(LYZEP) {
 410    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 411    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 412    // read the input file from the second lyz run
 413    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 414    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 415    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 416      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 417      break;
 418    }
 419    else {
 420      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 421      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 422      else {
 423        cout<<"LYZEP input file/list read..."<<endl;
 424        ep   ->SetSecondRunList(inputListLYZEP);
 425        lyzep->SetSecondRunList(inputListLYZEP);
 426        ep   ->Init();
 427        lyzep->Init();
 428      }
 429    }
 430  }
 431  //---------------------------------------------------------------------------------------
 432
 433  // set the global event parameters:
 434  eventMakerOnTheFly->SetNoOfLoops(iLoops);
 435  eventMakerOnTheFly->SetPhiRange(phiRange*TMath::Pi()/180.);
 436  eventMakerOnTheFly->SetPtRange(ptRange);
 437  eventMakerOnTheFly->SetEtaRange(etaRange);
 438  eventMakerOnTheFly->SetNonflowSectorMin(nonflowSectorMin*TMath::Pi()/180.);
 439  eventMakerOnTheFly->SetNonflowSectorMax(nonflowSectorMax*TMath::Pi()/180.);
 440  eventMakerOnTheFly->SetUseGlauberModel(bUseGlauberModel);
 441  eventMakerOnTheFly->SetCreateJets(bCreateJets);
 442  eventMakerOnTheFly->SetJetProbability(jetProbability);
 443  eventMakerOnTheFly->SetJetTracksFraction(jetTracksFraction);
 444  eventMakerOnTheFly->SetJetCone(jetCone);
 445
 446  if(!bUseGlauberModel)
 447  {
 448   if(bMultDistrOfRPsIsGauss)
 449   {
 450    eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 451    eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
 452    eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
 453   } else
 454     {
 455      eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 456      eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
 457      eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP);
 458     }
 459  } // end of if(!bUseGlauberModel)
 460
 461  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
 462  eventMakerOnTheFly->SetPtDependentHarmonicV1(bPtDependentHarmonicV1);
 463  eventMakerOnTheFly->SetEtaDependentHarmonicV1(bEtaDependentHarmonicV1);
 464  eventMakerOnTheFly->SetPtDependentHarmonicV2(bPtDependentHarmonicV2);
 465  eventMakerOnTheFly->SetEtaDependentHarmonicV2(bEtaDependentHarmonicV2);
 466  eventMakerOnTheFly->SetPtDependentHarmonicV4(bPtDependentHarmonicV4);
 467  eventMakerOnTheFly->SetEtaDependentHarmonicV4(bEtaDependentHarmonicV4);
 468  // V1:
 469  if(!(bPtDependentHarmonicV1 || bEtaDependentHarmonicV1))
 470  {
 471   if(!bUseGlauberModel)
 472   {
 473    eventMakerOnTheFly->SetV1RP(dV1RP);
 474    eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
 475   }
 476  } else // (pt,eta) dependent V1
 477    {
 478     eventMakerOnTheFly->SetV1vsPtEtaMax(dV1vsPtEtaMax);
 479     if(bPtDependentHarmonicV1)
 480     {
 481      eventMakerOnTheFly->SetV1PtCutOff(dV1PtCutOff);
 482     }
 483    }
 484  // V2:
 485  if(!(bPtDependentHarmonicV2 || bEtaDependentHarmonicV2)) // constant V2
 486  {
 487   if(!bUseGlauberModel)
 488   {
 489    eventMakerOnTheFly->SetConstantV2IsSampledFromGauss(bConstantV2IsSampledFromGauss);
 490    if(bConstantV2IsSampledFromGauss) // Gauss
 491    {
 492     eventMakerOnTheFly->SetV2RP(dV2RP);
 493     eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);
 494    } else // uniform
 495      {
 496       eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
 497       eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);
 498      }
 499   } // end of if(!bUseGlauberModel)
 500  } else // (pt,eta) dependent V2
 501    {
 502     eventMakerOnTheFly->SetV2vsPtEtaMax(dV2vsPtEtaMax);
 503     if(bPtDependentHarmonicV2)
 504     {
 505      eventMakerOnTheFly->SetV2PtCutOff(dV2PtCutOff);
 506     }
 507     if(bEtaDependentHarmonicV2)
 508     {
 509      eventMakerOnTheFly->SetV2vsEtaSpread(dV2vsEtaSpread);
 510     }
 511    }
 512  // V4:
 513  if(!(bPtDependentHarmonicV4 || bEtaDependentHarmonicV4))
 514  {
 515   if(!bUseGlauberModel)
 516   {
 517    eventMakerOnTheFly->SetV4RP(dV4RP);
 518    eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
 519   }
 520  }
 521
 522  // non-uniform acceptance:
 523  if(!uniformAcceptance)
 524  {
 525   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
 526   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
 527   eventMakerOnTheFly->SetFirstSectorProbability(p1);
 528   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
 529   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
 530   eventMakerOnTheFly->SetSecondSectorProbability(p2);
 531  }
 532
 533  // simple cuts for RPs and POIs:
 534  AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
 535  cutsRP->SetPtMax(ptMaxRP);
 536  cutsRP->SetPtMin(ptMinRP);
 537  cutsRP->SetEtaMax(etaMaxRP);
 538  cutsRP->SetEtaMin(etaMinRP);
 539  cutsRP->SetPhiMax(phiMaxRP*TMath::Pi()/180.);
 540  cutsRP->SetPhiMin(phiMinRP*TMath::Pi()/180.);
 541  cutsRP->SetPID(pidRP);
 542
 543  AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
 544  cutsPOI->SetPtMax(ptMaxPOI);
 545  cutsPOI->SetPtMin(ptMinPOI);
 546  cutsPOI->SetEtaMax(etaMaxPOI);
 547  cutsPOI->SetEtaMin(etaMinPOI);
 548  cutsPOI->SetPhiMax(phiMaxPOI*TMath::Pi()/180.);
 549  cutsPOI->SetPhiMin(phiMinPOI*TMath::Pi()/180.);
 550  cutsPOI->SetPID(pidPOI);
 551
 552  //---------------------------------------------------------------------------------------
 553  // create and analyze events 'on the fly':
 554
 555  for(Int_t i=0;i<nEvts;i++) {
 556    // creating the event with above settings:
 557    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly(cutsRP,cutsPOI);
 558
 559    // analyzing the created event 'on the fly':
 560    // do flow analysis for various methods:
 561    if(MCEP)    mcep->Make(event);
 562    if(QC)      qc->Make(event);
 563    if(GFC)     gfc->Make(event);
 564    if(FQD)     fqd->Make(event);
 565    if(LYZ1SUM) lyz1sum->Make(event);
 566    if(LYZ1PROD)lyz1prod->Make(event);
 567    if(LYZ2SUM) lyz2sum->Make(event);
 568    if(LYZ2PROD)lyz2prod->Make(event);
 569    if(LYZEP)   lyzep->Make(event,ep);
 570    if(SP)      sp->Make(event);
 571    if(MH)      mh->Make(event);
 572    if(NL)      nl->Make(event);
 573    if(MCEP_AH) mcep_ah->Make(event);
 574
 575    delete event;
 576  } // end of for(Int_t i=0;i<nEvts;i++)
 577  //---------------------------------------------------------------------------------------
 578
 579  //---------------------------------------------------------------------------------------
 580  // create a new file which will hold the final results of all methods:
 581  TString outputFileName = "AnalysisResults.root";
 582  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 583  // create a new file for each method wich will hold list with final results:
 584  const Int_t nMethods = 13;
 585  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
 586  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 587  TString fileName[nMethods];
 588  for(Int_t i=0;i<nMethods;i++)
 589  {
 590   // form a file name for each method:
 591   fileName[i]+="output";
 592   fileName[i]+=method[i].Data();
 593   fileName[i]+="analysis";
 594   dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 595  }
 596
 597  // calculating and storing the final results of default flow analysis:
 598  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
 599  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
 600  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
 601  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
 602  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
 603  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
 604  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
 605  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
 606  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
 607  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
 608  if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
 609  if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
 610  if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
 611  //---------------------------------------------------------------------------------------
 612
 613  outputFile->Close();
 614  delete outputFile;
 615
 616  cout<<endl;
 617  cout<<endl;
 618  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 619  cout<<endl;
 620
 621  timer.Stop();
 622  cout << endl;
 623  timer.Print();
 624 }
 625
 626 void SetupPar(char* pararchivename)
 627 {
 628   //Load par files, create analysis libraries
 629   //For testing, if par file already decompressed and modified
 630   //classes then do not decompress.
 631
 632   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 633   TString parpar(Form("%s.par", pararchivename)) ;
 634   if ( gSystem->AccessPathName(parpar.Data()) ) {
 635     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 636     TString processline(Form(".! make %s", parpar.Data())) ;
 637     gROOT->ProcessLine(processline.Data()) ;
 638     gSystem->ChangeDirectory(cdir) ;
 639     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 640     gROOT->ProcessLine(processline.Data()) ;
 641   }
 642   if ( gSystem->AccessPathName(pararchivename) ) {
 643     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 644     gROOT->ProcessLine(processline.Data());
 645   }
 646
 647   TString ocwd = gSystem->WorkingDirectory();
 648   gSystem->ChangeDirectory(pararchivename);
 649
 650   // check for BUILD.sh and execute
 651   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 652     printf("*******************************\n");
 653     printf("*** Building PAR archive    ***\n");
 654     cout<<pararchivename<<endl;
 655     printf("*******************************\n");
 656
 657     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 658       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 659       return -1;
 660     }
 661   }
 662   // check for SETUP.C and execute
 663   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 664     printf("*******************************\n");
 665     printf("*** Setup PAR archive       ***\n");
 666     cout<<pararchivename<<endl;
 667     printf("*******************************\n");
 668     gROOT->Macro("PROOF-INF/SETUP.C");
 669   }
 670
 671   gSystem->ChangeDirectory(ocwd.Data());
 672   printf("Current dir: %s\n", ocwd.Data());
 673 }
 674
 675 void CheckUserSettings()
 676 {
 677  // Check if user settings make sense before taking off:
 678
 679  if (LYZ1SUM && LYZ2SUM)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 680  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 681  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 682  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 683
 684  if(!uniformAcceptance && phimin1 > phimax1)
 685  {
 686   cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
 687   break;
 688  }
 689
 690  if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
 691  {
 692   cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
 693   break;
 694  }
 695
 696  if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) ||
 697     (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
 698  {
 699   cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
 700   break;
 701  }
 702
 703  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
 704  {
 705   cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
 706   break;
 707  }
 708
 709  if(bUseGlauberModel)
 710  {
 711   if(bPtDependentHarmonicV1||bEtaDependentHarmonicV1||
 712      bPtDependentHarmonicV2||bEtaDependentHarmonicV2||
 713      bPtDependentHarmonicV4||bEtaDependentHarmonicV4)
 714   {
 715    cout<<"WARNING: When using Glauber model pt and/or eta dependence of flow harmonics is not supported !!!!"<<endl;
 716    cout<<"         Set all booleans bPtDependentHarmonic* to kFALSE in the macro."<<endl;
 717    exit(0);
 718   }
 719  }
 720
 721  if(bPtDependentHarmonicV4 && !bPtDependentHarmonicV2))
 722  {
 723   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
 724   cout<<"         You must also set bPtDependentHarmonicV2 = kTRUE in the macro."<<endl;
 725   exit(0);
 726  }
 727
 728  if(bEtaDependentHarmonicV4 && !bEtaDependentHarmonicV2))
 729  {
 730   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
 731   cout<<"         You must also set bEtaDependentHarmonicV2 = kTRUE in the macro."<<endl;
 732   exit(0);
 733  }
 734
 735  if(bPtDependentHarmonicV4 && (bEtaDependentHarmonicV4 != bEtaDependentHarmonicV2))
 736  {
 737   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
 738   cout<<"         You must also have bEtaDependentHarmonicV2 = bEtaDependentHarmonicV4 in the macro."<<endl;
 739   exit(0);
 740  }
 741
 742  if(bEtaDependentHarmonicV4 && (bPtDependentHarmonicV4 != bPtDependentHarmonicV2))
 743  {
 744   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
 745   cout<<"         You must also have bPtDependentHarmonicV2 = bPtDependentHarmonicV4 in the macro."<<endl;
 746   exit(0);
 747  }
 748
 749  if(!uniformAcceptance && bCreateJets)
 750  {
 751   cout<<"WARNING: Jet creation is supported only for uniform acceptance !!!!"<<endl;
 752   exit(0);
 753  }
 754
 755 } // end of void CheckUserSettings()
 756
 757 void LoadLibraries(const anaModes mode) {
 758
 759   //--------------------------------------
 760   // Load the needed libraries most of them already loaded by aliroot
 761   //--------------------------------------
 762   //gSystem->Load("libTree");
 763   gSystem->Load("libGeom");
 764   gSystem->Load("libVMC");
 765   gSystem->Load("libXMLIO");
 766   gSystem->Load("libPhysics");
 767
 768   //----------------------------------------------------------
 769   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 770   //----------------------------------------------------------
 771   if (mode==mLocal) {
 772     //--------------------------------------------------------
 773     // If you want to use already compiled libraries
 774     // in the aliroot distribution
 775     //--------------------------------------------------------
 776     gSystem->Load("libSTEERBase");
 777     gSystem->Load("libESD");
 778     gSystem->Load("libAOD");
 779     gSystem->Load("libANALYSIS");
 780     gSystem->Load("libANALYSISalice");
 781     gSystem->Load("libCORRFW");
 782     cerr<<"libCORRFW loaded..."<<endl;
 783     gSystem->Load("libPWG2flowCommon");
 784     cerr<<"libPWG2flowCommon loaded..."<<endl;
 785     gSystem->Load("libPWG2flowTasks");
 786     cerr<<"libPWG2flowTasks loaded..."<<endl;
 787   }
 788
 789   else if (mode == mLocalPAR) {
 790     //--------------------------------------------------------
 791     //If you want to use root and par files from aliroot
 792     //--------------------------------------------------------
 793      //If you want to use root and par files from aliroot
 794     //--------------------------------------------------------
 795     SetupPar("STEERBase");
 796     SetupPar("ESD");
 797     SetupPar("AOD");
 798     SetupPar("ANALYSIS");
 799     SetupPar("ANALYSISalice");
 800     SetupPar("PWG2AOD");
 801     SetupPar("CORRFW");
 802     SetupPar("PWG2flowCommon");
 803     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 804     SetupPar("PWG2flowTasks");
 805     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 806   }
 807
 808   //---------------------------------------------------------
 809   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 810   //---------------------------------------------------------
 811   else if (mode==mLocalSource) {
 812
 813     // In root inline compile
 814
 815     // Constants
 816     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 817     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 818
 819     // Flow event
 820     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 821     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 822     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 823     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 824
 825     // Output histosgrams
 826     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 827     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 828     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 829     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 830
 831     // Functions needed for various methods
 832     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 833     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 834
 835     // Flow Analysis code for various methods
 836     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 837     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 838     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 839     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 840     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 841     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 842     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
 843     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx+");
 844     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithNestedLoops.cxx+");
 845
 846     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 847     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
 848
 849     cout << "finished loading macros!" << endl;
 850
 851   }
 852
 853 }
 854
 855