PWG2/FLOW/macros/runFlowAnalysisOnTheFly.C

   1 #include "TStopwatch.h"
   2 #include "TObjArray"
   3 #include "Riostream.h"
   4 #include "TFile.h"
   5
   6 //--------------------------------------------------------------------------------------
   7 // RUN SETTINGS
   8 // flow analysis method can be: (set to kTRUE or kFALSE)
   9 Bool_t SP       = kTRUE;
  10 Bool_t LYZ1SUM  = kFALSE;
  11 Bool_t LYZ1PROD = kFALSE;
  12 Bool_t LYZ2SUM  = kFALSE;
  13 Bool_t LYZ2PROD = kFALSE;
  14 Bool_t LYZEP    = kFALSE;
  15 Bool_t GFC      = kFALSE;
  16 Bool_t QC       = kTRUE;
  17 Bool_t FQD      = kFALSE;
  18 Bool_t MCEP     = kTRUE;
  19 //--------------------------------------------------------------------------------------
  20
  21 // Weights
  22 // use weights for Q-vector:
  23 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
  24 Bool_t usePtWeights  = kFALSE; // pt weights
  25 Bool_t useEtaWeights = kFALSE; // eta weights
  26
  27 // Run same flow analysis method but with different settings/aims
  28 // You will have to label each setting/aim with your own label (see examples bellow):
  29 Bool_t GFC_Additional_Analysis = kFALSE;
  30 Bool_t QC_Additional_Analysis  = kFALSE;
  31 Bool_t FQD_Additional_Analysis = kFALSE;
  32
  33 // Define the range for eta subevents
  34 Double_t minA = -0.9;
  35 Double_t maxA = -0.01;
  36 Double_t minB = 0.01;
  37 Double_t maxB = 0.9;
  38
  39 // Parameters for the simulation of events 'on the fly':
  40 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs)
  41
  42 //===SEED========================================================
  43 // use always the same seed for random generators.
  44 // usage of same seed (kTRUE) is relevant in two cases:
  45 // 1.) If you want to use LYZ method to calcualte differential flow;
  46 // 2.) If you want to use phi weights for GFC, QC and FQD
  47 Bool_t bSameSeed = kTRUE;
  48
  49
  50 //===NONFLOW=============================================
  51 // number of times to use each track (to simulate nonflow)
  52 Int_t iLoops = 1;
  53
  54
  55 //===FLOW HARMONICS===============================================================
  56 // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
  57 Bool_t bConstantHarmonics = kTRUE;
  58
  59 // if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
  60 Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
  61 Double_t dV2RPMax = 0.20; // maximum value of V2(pt) for pt >= 2GeV
  62
  63 // 1.) if kTRUE  = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
  64 //                 mean = dV2RP and spread = dV2SpreadRP
  65 // 2.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from
  66 //                 interval [dMinV2RP,dMaxV2RP]
  67 // 3.) for a fixed elliptic flow use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP
  68 Bool_t bV2DistrOfRPsIsGauss = kTRUE;
  69
  70 // if you use constant harmonics (bConstantHarmonics = kTRUE) (i.e. no pt dependence):
  71 Double_t dV2RP = 0.05;       // elliptic flow of RPs (if sampled from Gaussian)
  72 Double_t dV2SpreadRP = 0.0;  // elliptic flow spread of RPs (if sampled from Gaussian)
  73 Double_t dMinV2RP = 0.02;    // minimal elliptic flow of RPs (if sampled uniformly)
  74 Double_t dMaxV2RP = 0.08;    // maximal elliptic flow of RPs (if sampled uniformly)
  75
  76 Double_t dV1RP = 0.0; // directed flow of RPs
  77 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
  78
  79 Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
  80 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
  81
  82
  83 //===MULTIPLICITY===============================================================
  84 // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
  85 //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
  86 // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from
  87 //                 interval [iMinMultOfRP,iMaxMultOfRP]
  88 // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
  89 Bool_t bMultDistrOfRPsIsGauss = kTRUE;
  90
  91 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
  92 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
  93 Int_t iMinMultOfRP = 50;             // minimal multiplicity of RPs (if sampled uniformly)
  94 Int_t iMaxMultOfRP = 500;             // maximal multiplicity of RPs (if sampled uniformly)
  95
  96
  97
  98 //===DETECTOR ACCEPTANCE===============================================================
  99
 100 // 1.) if kTRUE = detectors has uniform azimuthal acceptance
 101 // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors.
 102 //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles
 103 //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
 104 //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2
 105 //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is
 106 //                 simulated.
 107 Bool_t uniformAcceptance = kTRUE;
 108
 109 // settings for non-uniform acceptance:
 110 // Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
 111
 112 // 1st non-uniform sector:
 113 Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
 114 Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
 115 Double_t p1 = 0.6;     // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
 116                         // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
 117
 118 // 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):
 119 Double_t phimin2 = 200.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
 120 Double_t phimax2 = 210.0; // second non-uniform sector ends at this azimuth
 121 Double_t p2 = 0.01;
 122
 123 //=====================================================================================================
 124
 125 enum anaModes {mLocal,mLocalSource,mLocalPAR};
 126 // mLocal: Analyze data on your computer using aliroot
 127 // mLocalPAR: Analyze data on your computer using root + PAR files
 128 // mLocalSource: Analyze data on your computer using root + source files
 129
 130 int runFlowAnalysisOnTheFly( Int_t nEvts=44, Int_t mode=mLocal)
 131 {
 132  TStopwatch timer;
 133  timer.Start();
 134
 135  if (LYZ1SUM && LYZ2SUM)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 136  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 137  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 138  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 139
 140  if(!uniformAcceptance && phimin1 > phimax1)
 141  {
 142   cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
 143   break;
 144  }
 145
 146  if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
 147  {
 148   cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
 149   break;
 150  }
 151
 152  if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) ||
 153     (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
 154  {
 155   cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
 156   break;
 157  }
 158
 159  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
 160  {
 161   cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
 162   break;
 163  }
 164
 165  cout<<endl;
 166  cout<<endl;
 167  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
 168  cout<<endl;
 169
 170  cout<<endl;
 171  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
 172  cout<<endl;
 173  cout<<endl;
 174
 175  LoadLibraries(mode);
 176
 177  // Initialize the seed for random generator
 178  UInt_t sseed = 0;
 179
 180  if(bSameSeed)
 181  {
 182   sseed = 44; // the default constant value for seed for random generators
 183  }
 184
 185  if(!bSameSeed)
 186  {
 187   TTimeStamp dt;
 188   sseed = dt.GetNanoSec()/1000;
 189  }
 190
 191  cout<<endl;
 192  cout<<"Seed for the random generators is "<<sseed<<endl;
 193  cout<<endl;
 194
 195  //---------------------------------------------------------------------------------------
 196  // If the weights are used:
 197  TFile *fileWithWeights = NULL;
 198  TList *listWithWeights = NULL;
 199
 200  if(usePhiWeights||usePtWeights||useEtaWeights) {
 201    fileWithWeights = TFile::Open("weights.root","READ");
 202    if(fileWithWeights) {
 203      listWithWeights = (TList*)fileWithWeights->Get("weights");
 204    }
 205    else
 206      {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 207       break;
 208      }
 209  }
 210
 211  //---------------------------------------------------------------------------------------
 212  // Initialize the flowevent maker
 213  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
 214  eventMakerOnTheFly->Init();
 215  //set the range for eta subevents
 216  eventMakerOnTheFly->SetSubeventEtaRange(minA,maxA,minB,maxB);
 217
 218  //---------------------------------------------------------------------------------------
 219  // Initialize all the flow methods for default analysis:
 220  AliFlowAnalysisWithQCumulants *qc = NULL;
 221  AliFlowAnalysisWithCumulants *gfc = NULL;
 222  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
 223  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
 224  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
 225  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
 226  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
 227  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 228  AliFlowAnalysisWithScalarProduct *sp = NULL;
 229  AliFlowAnalysisWithMCEventPlane *mcep = NULL;
 230
 231  // MCEP = monte carlo event plane
 232  if (MCEP) {
 233    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 234    // mcep->SetHarmonic(2); // default is v2
 235    mcep->Init();
 236  }
 237
 238  // QC = Q-cumulants
 239  if(QC) {
 240    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 241    if(listWithWeights) qc->SetWeightsList(listWithWeights);
 242    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 243    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 244    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 245    // qc->SetHarmonic(2); // default is v2
 246    // qc->SetApplyCorrectionForNUA(kTRUE); // default
 247    // qc->SetCalculate2DFlow(kFALSE); // default
 248    // qc->SetMultiplicityWeight("combinations"); // default
 249    // qc->SetMultiplicityWeight("unit");
 250    // qc->SetMultiplicityWeight("multiplicity");
 251    qc->Init();
 252  }
 253
 254  // GFC = Generating Function Cumulants
 255  if(GFC) {
 256    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 257    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 258    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 259    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 260    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 261    gfc->Init();
 262  }
 263
 264  // FQD = Fitting q-distribution
 265  if(FQD) {
 266    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 267    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 268    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 269    fqd->Init();
 270  }
 271
 272  // SP = Scalar Product
 273  if(SP) {
 274    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 275    if(listWithWeights) sp->SetWeightsList(listWithWeights);
 276    sp->SetUsePhiWeights(usePhiWeights);
 277    sp->Init();
 278  }
 279
 280  // LYZ1 = Lee-Yang Zeroes first run
 281  if(LYZ1SUM) {
 282    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 283    lyz1sum->SetFirstRun(kTRUE);
 284    lyz1sum->SetUseSum(kTRUE);
 285    lyz1sum->Init();
 286  }
 287  if(LYZ1PROD) {
 288    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 289    lyz1prod->SetFirstRun(kTRUE);
 290    lyz1prod->SetUseSum(kFALSE);
 291    lyz1prod->Init();
 292  }
 293  // LYZ2 = Lee-Yang Zeroes second run
 294  if(LYZ2SUM) {
 295    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 296    // read the input file from the first run
 297    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 298    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 299    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 300      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 301      break;
 302    }
 303    else {
 304      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 305      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
 306      else {
 307        cout<<"LYZ2SUM input file/list read..."<<endl;
 308        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 309        lyz2sum->SetFirstRun(kFALSE);
 310        lyz2sum->SetUseSum(kTRUE);
 311        lyz2sum->Init();
 312      }
 313    }
 314  }
 315  if(LYZ2PROD) {
 316    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 317    // read the input file from the first run
 318    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 319    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 320    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 321      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 322      break;
 323    }
 324    else {
 325      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 326      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
 327      else {
 328        cout<<"LYZ2PROD input file/list read..."<<endl;
 329        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 330        lyz2prod->SetFirstRun(kFALSE);
 331        lyz2prod->SetUseSum(kFALSE);
 332        lyz2prod->Init();
 333      }
 334    }
 335  }
 336
 337  // LYZEP = Lee-Yang Zeroes event plane
 338  if(LYZEP) {
 339    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 340    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 341    // read the input file from the second lyz run
 342    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 343    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 344    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 345      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 346      break;
 347    }
 348    else {
 349      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 350      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 351      else {
 352        cout<<"LYZEP input file/list read..."<<endl;
 353        ep   ->SetSecondRunList(inputListLYZEP);
 354        lyzep->SetSecondRunList(inputListLYZEP);
 355        ep   ->Init();
 356        lyzep->Init();
 357      }
 358    }
 359  }
 360  //---------------------------------------------------------------------------------------
 361
 362  //---------------------------------------------------------------------------------------
 363  // Initialize all the flow methods for additional analysis with different settings/aims
 364  // Label each setting/aim with different label !!!!
 365
 366  // GFC:
 367  TString gfcDefaultName = "outputGFCanalysis";
 368  // 1.) GFC analysis for elliptic flow with r0 = 1.5:
 369  AliFlowAnalysisWithCumulants *gfc_1;
 370  TString gfcAnalysisLabels_1 = "_r0_1.5"; // all histograms and output file name will have this label
 371  TString gfcOutputFileName_1;
 372  gfcOutputFileName_1 = gfcDefaultName.Data();
 373  gfcOutputFileName_1 += gfcAnalysisLabels_1.Data();
 374  gfcOutputFileName_1 += ".root";
 375  if(GFC_Additional_Analysis)
 376  {
 377   gfc_1 = new AliFlowAnalysisWithCumulants();
 378   //gfc_1->SetAnalysisLabel(gfcAnalysisLabels_1.Data());
 379   if(listWithWeights) gfc_1->SetWeightsList(listWithWeights);
 380   if(usePhiWeights) gfc_1->SetUsePhiWeights(usePhiWeights);
 381   if(usePtWeights) gfc_1->SetUsePtWeights(usePtWeights);
 382   if(useEtaWeights) gfc_1->SetUseEtaWeights(useEtaWeights);
 383   gfc_1->Init();
 384  }
 385
 386  // QC:
 387  TString qcDefaultName = "outputQCanalysis";
 388  // 1.) QC analysis for directed flow:
 389  AliFlowAnalysisWithQCumulants *qc_1;
 390  TString qcAnalysisLabels_1 = "_v1"; // all histograms and output file name will have this label
 391  TString qcOutputFileName_1;
 392  qcOutputFileName_1 = qcDefaultName.Data();
 393  qcOutputFileName_1 += qcAnalysisLabels_1.Data();
 394  qcOutputFileName_1 += ".root";
 395  if(QC_Additional_Analysis)
 396  {
 397   qc_1 = new AliFlowAnalysisWithQCumulants();
 398   //qc_1->SetAnalysisLabel(qcAnalysisLabels_1->Data());
 399   if(listWithWeights) qc_1->SetWeightsList(listWithWeights);
 400   if(usePhiWeights) qc_1->SetUsePhiWeights(usePhiWeights);
 401   qc_1->Init();
 402  }
 403
 404  // FQD:
 405  TString fqdDefaultName = "outputFQDanalysis";
 406  // 1.) FQD fitting with fixed sigma:
 407  AliFlowAnalysisWithFittingQDistribution *fqd_1;
 408  TString fqdAnalysisLabels_1 = "_fixedSigma"; // all histograms and output file name will have this label
 409  TString fqdOutputFileName_1;
 410  fqdOutputFileName_1 = fqdDefaultName.Data();
 411  fqdOutputFileName_1 += fqdAnalysisLabels_1.Data();
 412  fqdOutputFileName_1 += ".root";
 413  if(FQD_Additional_Analysis)
 414  {
 415   fqd_1 = new AliFlowAnalysisWithFittingQDistribution();
 416   //fqd_1->SetAnalysisLabel(fqdAnalysisLabels_1->Data());
 417   if(listWithWeights) fqd_1->SetWeightsList(listWithWeights);
 418   if(usePhiWeights) fqd_1->SetUsePhiWeights(usePhiWeights);
 419   fqd_1->Init();
 420  }
 421  //---------------------------------------------------------------------------------------
 422
 423  // set the global event parameters:
 424  eventMakerOnTheFly->SetNoOfLoops(iLoops);
 425
 426  if(bMultDistrOfRPsIsGauss)
 427  {
 428   eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 429   eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
 430   eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
 431  } else
 432    {
 433     eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 434     eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
 435     eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP);
 436    }
 437
 438  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
 439
 440  eventMakerOnTheFly->SetV1RP(dV1RP);
 441  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
 442  eventMakerOnTheFly->SetV4RP(dV4RP);
 443  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
 444
 445  // constant harmonic V2:
 446  if(bConstantHarmonics)
 447  {
 448   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 449   if(bV2DistrOfRPsIsGauss)
 450   {
 451    eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 452    eventMakerOnTheFly->SetV2RP(dV2RP);
 453    eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);
 454   } else
 455     {
 456      eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 457      eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
 458      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);
 459     }
 460  }
 461
 462  // (pt,eta) dependent harmonic V2:
 463  if(!bConstantHarmonics)
 464  {
 465   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 466   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
 467   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);
 468  }
 469
 470  // non-uniform acceptance:
 471  if(!uniformAcceptance)
 472  {
 473   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
 474   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
 475   eventMakerOnTheFly->SetFirstSectorProbability(p1);
 476   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
 477   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
 478   eventMakerOnTheFly->SetSecondSectorProbability(p2);
 479  }
 480
 481  //---------------------------------------------------------------------------------------
 482  // create and analyze events 'on the fly':
 483
 484  for(Int_t i=0;i<nEvts;i++) {
 485    // creating the event with above settings:
 486    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly();
 487
 488    // analyzing the created event 'on the fly':
 489    // do flow analysis for various methods:
 490    if(MCEP)    mcep->Make(event);
 491    if(QC)      qc->Make(event);
 492    if(GFC)     gfc->Make(event);
 493    if(FQD)     fqd->Make(event);
 494    if(LYZ1SUM) lyz1sum->Make(event);
 495    if(LYZ1PROD)lyz1prod->Make(event);
 496    if(LYZ2SUM) lyz2sum->Make(event);
 497    if(LYZ2PROD)lyz2prod->Make(event);
 498    if(LYZEP)   lyzep->Make(event,ep);
 499    if(SP)      sp->Make(event);
 500
 501    if(GFC_Additional_Analysis)
 502    {
 503     // r0 = 1.5:
 504     gfc_1->Make(event);
 505    }
 506    if(QC_Additional_Analysis)
 507    {
 508     // v1:
 509     qc_1->Make(event);
 510    }
 511    if(FQD_Additional_Analysis)
 512    {
 513     // fixed sigma:
 514     fqd_1->Make(event);
 515    }
 516
 517    delete event;
 518  } // end of for(Int_t i=0;i<nEvts;i++)
 519  //---------------------------------------------------------------------------------------
 520
 521
 522
 523  //---------------------------------------------------------------------------------------
 524  // create a new file which will hold the final results of all methods:
 525  TString outputFileName = "AnalysisResults.root";
 526  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 527  // create a new file for each method wich will hold list with final results:
 528  const Int_t nMethods = 10;
 529  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP"};
 530  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 531  TString fileName[nMethods];
 532  for(Int_t i=0;i<nMethods;i++)
 533  {
 534   // form a file name for each method:
 535   fileName[i]+="output";
 536   fileName[i]+=method[i].Data();
 537   fileName[i]+="analysis";
 538   dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 539  }
 540
 541  // calculating and storing the final results of default flow analysis:
 542  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
 543  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
 544  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
 545  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
 546  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
 547  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
 548  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
 549  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
 550  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
 551  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
 552  //---------------------------------------------------------------------------------------
 553
 554  outputFile->Close();
 555  delete outputFile;
 556
 557  cout<<endl;
 558  cout<<endl;
 559  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 560  cout<<endl;
 561
 562  timer.Stop();
 563  cout << endl;
 564  timer.Print();
 565 }
 566
 567 void SetupPar(char* pararchivename)
 568 {
 569   //Load par files, create analysis libraries
 570   //For testing, if par file already decompressed and modified
 571   //classes then do not decompress.
 572
 573   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 574   TString parpar(Form("%s.par", pararchivename)) ;
 575   if ( gSystem->AccessPathName(parpar.Data()) ) {
 576     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 577     TString processline(Form(".! make %s", parpar.Data())) ;
 578     gROOT->ProcessLine(processline.Data()) ;
 579     gSystem->ChangeDirectory(cdir) ;
 580     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 581     gROOT->ProcessLine(processline.Data()) ;
 582   }
 583   if ( gSystem->AccessPathName(pararchivename) ) {
 584     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 585     gROOT->ProcessLine(processline.Data());
 586   }
 587
 588   TString ocwd = gSystem->WorkingDirectory();
 589   gSystem->ChangeDirectory(pararchivename);
 590
 591   // check for BUILD.sh and execute
 592   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 593     printf("*******************************\n");
 594     printf("*** Building PAR archive    ***\n");
 595     cout<<pararchivename<<endl;
 596     printf("*******************************\n");
 597
 598     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 599       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 600       return -1;
 601     }
 602   }
 603   // check for SETUP.C and execute
 604   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 605     printf("*******************************\n");
 606     printf("*** Setup PAR archive       ***\n");
 607     cout<<pararchivename<<endl;
 608     printf("*******************************\n");
 609     gROOT->Macro("PROOF-INF/SETUP.C");
 610   }
 611
 612   gSystem->ChangeDirectory(ocwd.Data());
 613   printf("Current dir: %s\n", ocwd.Data());
 614 }
 615
 616 void LoadLibraries(const anaModes mode) {
 617
 618   //--------------------------------------
 619   // Load the needed libraries most of them already loaded by aliroot
 620   //--------------------------------------
 621   //gSystem->Load("libTree");
 622   gSystem->Load("libGeom");
 623   gSystem->Load("libVMC");
 624   gSystem->Load("libXMLIO");
 625   gSystem->Load("libPhysics");
 626
 627   //----------------------------------------------------------
 628   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 629   //----------------------------------------------------------
 630   if (mode==mLocal) {
 631     //--------------------------------------------------------
 632     // If you want to use already compiled libraries
 633     // in the aliroot distribution
 634     //--------------------------------------------------------
 635     gSystem->Load("libSTEERBase");
 636     gSystem->Load("libESD");
 637     gSystem->Load("libAOD");
 638     gSystem->Load("libANALYSIS");
 639     gSystem->Load("libANALYSISalice");
 640     gSystem->Load("libCORRFW");
 641     cerr<<"libCORRFW loaded..."<<endl;
 642     gSystem->Load("libPWG2flowCommon");
 643     cerr<<"libPWG2flowCommon loaded..."<<endl;
 644     gSystem->Load("libPWG2flowTasks");
 645     cerr<<"libPWG2flowTasks loaded..."<<endl;
 646   }
 647
 648   else if (mode == mLocalPAR) {
 649     //--------------------------------------------------------
 650     //If you want to use root and par files from aliroot
 651     //--------------------------------------------------------
 652      //If you want to use root and par files from aliroot
 653     //--------------------------------------------------------
 654     SetupPar("STEERBase");
 655     SetupPar("ESD");
 656     SetupPar("AOD");
 657     SetupPar("ANALYSIS");
 658     SetupPar("ANALYSISalice");
 659     SetupPar("PWG2AOD");
 660     SetupPar("CORRFW");
 661     SetupPar("PWG2flowCommon");
 662     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 663     SetupPar("PWG2flowTasks");
 664     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 665   }
 666
 667   //---------------------------------------------------------
 668   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 669   //---------------------------------------------------------
 670   else if (mode==mLocalSource) {
 671
 672     // In root inline compile
 673
 674
 675     // Constants
 676     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 677     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 678     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 679
 680     // Flow event
 681     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 682     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 683     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 684
 685     // Cuts
 686     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 687
 688     // Output histosgrams
 689     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 690     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 691     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 692     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 693
 694     // Functions needed for various methods
 695     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 696     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 697
 698     // Flow Analysis code for various methods
 699     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 700     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 701     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 702     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 703     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 704     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 705     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
 706
 707     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 708     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
 709
 710     cout << "finished loading macros!" << endl;
 711
 712   }
 713
 714 }
 715
 716