PWG2/FLOW/macros/runStarFlowAnalysis.C

   1 //example script on what to do with the star events
   2 //run e.g. like this:
   3 //                    root runStarFlowAnalysis.C
   4 // flow analysis method can be: (set to kTRUE or kFALSE)
   5 Bool_t MCEP     = kTRUE;
   6 Bool_t SP       = kTRUE;
   7 Bool_t GFC      = kTRUE;
   8 Bool_t QC       = kTRUE;
   9 Bool_t FQD      = kTRUE;
  10 Bool_t LYZ1SUM  = kTRUE;
  11 Bool_t LYZ1PROD = kTRUE;
  12 Bool_t LYZ2SUM  = kFALSE;
  13 Bool_t LYZ2PROD = kFALSE;
  14 Bool_t LYZEP    = kFALSE;
  15 Bool_t MH       = kFALSE; // mixed harmonics
  16 Bool_t NL       = kFALSE; // nested loops
  17 Bool_t MCEP_AH  = kFALSE; // MCEP in another harmonic
  18 //--------------------------------------------------------------------------------------
  19
  20 // Weights
  21 // use weights for Q-vector:
  22 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
  23 Bool_t usePtWeights  = kFALSE; // pt weights
  24 Bool_t useEtaWeights = kFALSE; // eta weights
  25
  26 // Define the range for eta subevents
  27 Double_t minA = -0.9;
  28 Double_t maxA = -0.01;
  29 Double_t minB = 0.01;
  30 Double_t maxB = 0.9;
  31
  32
  33 void  runStarFlowAnalysis(const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
  34 {
  35   gSystem->Load("libTree.so");
  36   gSystem->Load("libVMC.so");
  37   gSystem->Load("libPhysics.so");
  38   gSystem->Load("libPWG2flowCommon");
  39
  40   Int_t maxNumberOfEvents = 1000000;
  41
  42   TStopwatch timer;
  43   timer.Start();
  44
  45   if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
  46   if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
  47   if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
  48   if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
  49
  50
  51
  52   //define reference particles
  53   AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
  54   rpCuts->SetPtMin(0.05);
  55   rpCuts->SetPtMax(10.);
  56
  57   //define particles of interest
  58   AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
  59   poiCuts->SetPtMin(0.05);
  60   poiCuts->SetPtMax(10.);
  61
  62   //define event cuts
  63   AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
  64   starEventCuts-> SetCentralityIDMax(3);
  65   starEventCuts-> SetCentralityIDMin(3);
  66
  67   //if the weights are used:
  68   TFile *fileWithWeights = NULL;
  69   TList *listWithWeights = NULL;
  70
  71   if(usePhiWeights||usePtWeights||useEtaWeights) {
  72     fileWithWeights = TFile::Open("weights.root","READ");
  73     if(fileWithWeights) {
  74       listWithWeights = (TList*)fileWithWeights->Get("weights");
  75     }
  76     else
  77       {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
  78         break;
  79       }
  80   }
  81
  82   //---------------------------------------------------------------------------------------
  83   // Initialize all the flow methods for default analysis:
  84   AliFlowAnalysisWithQCumulants *qc = NULL;
  85   AliFlowAnalysisWithCumulants *gfc = NULL;
  86   AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
  87   AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
  88   AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
  89   AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
  90   AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
  91   AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
  92   AliFlowAnalysisWithScalarProduct *sp = NULL;
  93   AliFlowAnalysisWithMixedHarmonics *mh = NULL;
  94   AliFlowAnalysisWithNestedLoops *nl = NULL;
  95   AliFlowAnalysisWithMCEventPlane *mcep = NULL;
  96   AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;
  97
  98 // MCEP = monte carlo event plane
  99  if (MCEP) {
 100    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 101    //mcep->SetHarmonic(2); // default is v2
 102    //mcep->SetFlowOfResonances(kTRUE);
 103    mcep->Init();
 104  }
 105
 106  // MCEP = monte carlo event plane in another harmonic:
 107  if(MCEP_AH)
 108  {
 109   AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
 110   mcep_ah->SetHarmonic(1);
 111   mcep_ah->Init();
 112  }
 113
 114  // Mixed harmonics:
 115  if(MH)
 116  {
 117   AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
 118   mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
 119   mh->SetMinMultiplicity(100);
 120   mh->SetNoOfMultipicityBins(5);
 121   mh->SetMultipicityBinWidth(200);
 122   mh->Init();
 123  }
 124
 125  // NL = nested loops:
 126  if(NL) {
 127    AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
 128    nl->Init();
 129  }
 130
 131  // QC = Q-cumulants
 132  if(QC) {
 133    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 134    if(listWithWeights) qc->SetWeightsList(listWithWeights);
 135    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 136    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 137    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 138    // qc->SetHarmonic(2); // default is v2
 139    // qc->SetApplyCorrectionForNUA(kTRUE); // default
 140    // qc->SetCalculate2DFlow(kFALSE); // default
 141    // qc->SetMultiplicityWeight("combinations"); // default
 142    // qc->SetMultiplicityWeight("unit");
 143    // qc->SetMultiplicityWeight("multiplicity");
 144    qc->SetnBinsMult(10000);
 145    qc->SetMinMult(0);
 146    qc->SetMaxMult(10000);
 147    qc->Init();
 148  }
 149
 150  // GFC = Generating Function Cumulants
 151  if(GFC) {
 152    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 153    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 154    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 155    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 156    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 157    // calculation vs multiplicity:
 158    gfc->SetCalculateVsMultiplicity(kFALSE);
 159    gfc->SetnBinsMult(10000);
 160    gfc->SetMinMult(0);
 161    gfc->SetMaxMult(10000);
 162    // tuning of interpolating parameters:
 163    gfc->SetTuneParameters(kFALSE);
 164    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
 165    for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
 166    gfc->Init();
 167  }
 168
 169  // FQD = Fitting q-distribution
 170  if(FQD) {
 171    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 172    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 173    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 174    fqd->Init();
 175  }
 176
 177  // SP = Scalar Product
 178  if(SP) {
 179    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 180    if(listWithWeights) sp->SetWeightsList(listWithWeights);
 181    sp->SetUsePhiWeights(usePhiWeights);
 182    sp->Init();
 183  }
 184
 185  // LYZ1 = Lee-Yang Zeroes first run
 186  if(LYZ1SUM) {
 187    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 188    lyz1sum->SetFirstRun(kTRUE);
 189    lyz1sum->SetUseSum(kTRUE);
 190    lyz1sum->Init();
 191  }
 192  if(LYZ1PROD) {
 193    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 194    lyz1prod->SetFirstRun(kTRUE);
 195    lyz1prod->SetUseSum(kFALSE);
 196    lyz1prod->Init();
 197  }
 198  // LYZ2 = Lee-Yang Zeroes second run
 199  if(LYZ2SUM) {
 200    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 201    // read the input file from the first run
 202    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 203    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 204    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 205      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 206      break;
 207    }
 208    else {
 209      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 210      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
 211      else {
 212        cout<<"LYZ2SUM input file/list read..."<<endl;
 213        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 214        lyz2sum->SetFirstRun(kFALSE);
 215        lyz2sum->SetUseSum(kTRUE);
 216        lyz2sum->Init();
 217      }
 218    }
 219  }
 220  if(LYZ2PROD) {
 221    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 222    // read the input file from the first run
 223    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 224    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 225    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 226      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 227      break;
 228    }
 229    else {
 230      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 231      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
 232      else {
 233        cout<<"LYZ2PROD input file/list read..."<<endl;
 234        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 235        lyz2prod->SetFirstRun(kFALSE);
 236        lyz2prod->SetUseSum(kFALSE);
 237        lyz2prod->Init();
 238      }
 239    }
 240  }
 241
 242  // LYZEP = Lee-Yang Zeroes event plane
 243  if(LYZEP) {
 244    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 245    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 246    // read the input file from the second lyz run
 247    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 248    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 249    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 250      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 251      break;
 252    }
 253    else {
 254      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 255      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 256      else {
 257        cout<<"LYZEP input file/list read..."<<endl;
 258        ep   ->SetSecondRunList(inputListLYZEP);
 259        lyzep->SetSecondRunList(inputListLYZEP);
 260        ep   ->Init();
 261        lyzep->Init();
 262      }
 263    }
 264  }
 265  //---------------------------------------------------------------------------------------
 266
 267
 268   Int_t i=0;
 269   AliStarEventReader starReader(inputDataFiles) ;
 270   while ( starReader.GetNextEvent() )                                // Get next event
 271   {
 272     AliStarEvent* starEvent = starReader.GetEvent();
 273     if ( !starEventCuts->PassesCuts(starEvent) ) continue;              // Test if the event is good
 274
 275     AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts);  // make a flow event from a star event (aka "the magic")
 276     flowEvent->TagSubeventsInEta(minA, maxA, minB, maxB );
 277
 278     /////analysis here////////////////
 279
 280     // do flow analysis for various methods
 281     if(MCEP)    mcep->Make(flowEvent);
 282     if(QC)      qc->Make(flowEvent);
 283     if(GFC)     gfc->Make(flowEvent);
 284     if(FQD)     fqd->Make(flowEvent);
 285     if(LYZ1SUM) lyz1sum->Make(flowEvent);
 286     if(LYZ1PROD)lyz1prod->Make(flowEvent);
 287     if(LYZ2SUM) lyz2sum->Make(flowEvent);
 288     if(LYZ2PROD)lyz2prod->Make(flowEvent);
 289     if(LYZEP)   lyzep->Make(flowEvent,ep);
 290     if(SP)      sp->Make(flowEvent);
 291     if(MH)      mh->Make(flowEvent);
 292     if(NL)      nl->Make(flowEvent);
 293     if(MCEP_AH) mcep_ah->Make(flowEvent);
 294
 295     //////////////////////////////////
 296
 297     flowEvent->Print();
 298
 299     delete flowEvent;
 300
 301     i++;
 302     if (i>maxNumberOfEvents) break;
 303   }
 304
 305   //---------------------------------------------------------------------------------------
 306   // create a new file which will hold the final results of all methods:
 307   TString outputFileName = "AnalysisResults.root";
 308   TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
 309   // create a new file for each method wich will hold list with final results:
 310   const Int_t nMethods = 13;
 311   TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
 312   TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 313   TString fileName[nMethods];
 314   for(Int_t i=0;i<nMethods;i++)
 315     {
 316       // form a file name for each method:
 317       fileName[i]+="output";
 318       fileName[i]+=method[i].Data();
 319       fileName[i]+="analysis";
 320       dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 321     }
 322
 323   // calculating and storing the final results of default flow analysis:
 324   if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
 325   if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
 326   if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
 327   if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
 328   if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
 329   if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
 330   if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
 331   if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
 332   if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
 333   if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
 334   if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
 335   if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
 336   if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
 337   //---------------------------------------------------------------------------------------
 338
 339   outputFile->Close();
 340   delete outputFile;
 341
 342   cout<<endl;
 343   cout<<endl;
 344   cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 345   cout<<endl;
 346
 347   timer.Stop();
 348   cout << endl;
 349   timer.Print();
 350
 351   delete rpCuts;
 352   delete poiCuts;
 353   delete starEventCuts;
 354 }