PWG2/FLOW/Documentation/examples/runFlowOnTheFlyGlauber.C

   1 #include "TStopwatch.h"
   2 #include "Riostream.h"
   3 #include "TFile.h"
   4
   5 int runFlowOnTheFlyGlauber(Int_t nevents=100000, Int_t iseed=1)
   6 {
   7     gSystem->Load("libGeom");
   8     gSystem->Load("libVMC");
   9     gSystem->Load("libXMLIO");
  10     gSystem->Load("libPhysics");
  11     gSystem->Load("libPWG2flowCommon");
  12     gSystem->Load("libPWG2flowTools");
  13
  14     fMyTRandom3 = new TRandom3(iseed);
  15     gRandom->SetSeed(iseed);
  16     TFile *outputFile = new TFile("flowOnTheFlyGlauber.root","RECREATE");
  17
  18     TStopwatch timer;
  19     timer.Start();
  20
  21     runOneCentrality(3,nevents);
  22     runOneCentrality(4,nevents);
  23     runOneCentrality(5,nevents);
  24
  25     outputFile->Close();
  26     delete outputFile;
  27
  28     cout<<endl; cout<<" ---- Fini ---- "<<endl; cout<<endl;
  29     timer.Stop(); timer.Print();
  30 }
  31
  32 void runOneCentrality(Int_t centrality=3, Int_t nEvts=100)
  33 {
  34
  35     //centralities: 0: 0-5
  36     //              1: 5-10
  37     //              2: 10-20
  38     //              3: 20-30
  39     //              4: 30-40
  40     //              5: 40-50
  41     //              6: 50-60
  42     //              7: 60-70
  43     //              8: 70-80
  44
  45     Double_t centPerc[] = {0,5,10,20,30,40,50,60,70,80};
  46     Double_t xCross[] = {2.5,7.5,15.,25.,35.,45.,55.,65.,75.};
  47     Double_t xCrossErr[9] = {0.};
  48     Double_t impact[] = {0., 3.5, 4.95, 6.98, 8.55, 9.88, 11.4, 12.09, 13.06, 13.97};
  49     Double_t npart[] = {1000, 356, 305, 220, 155, 105, 67, 40, 21, 10, 4, 0};
  50
  51     Option_t *sysA="Pb";
  52     Option_t *sysB="Pb";
  53     Double_t signn=64; // inelastic nucleon nucleon cross section
  54     //const char *fname="GlauberMC_PbPb_ntuple.root"; // name output file
  55     //  AliGlauberMC::RunAndSaveNtuple(nevents,sysA,sysB,signn);
  56     // run the code to produce an ntuple:
  57
  58     Double_t mind=0.4;
  59     Double_t r=6.62;
  60     Double_t a=0.546;
  61
  62     AliGlauberMC glauber(sysA,sysB,signn);
  63     glauber.SetMinDistance(mind);
  64     glauber.Setr(r);
  65     glauber.Seta(a);
  66     glauber.SetDoPartProduction(kTRUE);
  67     glauber.SetBmin(impact[centrality]);
  68     glauber.SetBmax(impact[centrality+1]);
  69     glauber.SetdNdEtaType(AliGlauberMC::kNBDSV);
  70     glauber.GetdNdEtaParam()[0] = 2.49;    //npp
  71     glauber.GetdNdEtaParam()[1] = 1.7;  //ratioSgm2Mu
  72     glauber.GetdNdEtaParam()[2] = 0.13; //xhard
  73
  74
  75     AliFlowAnalysisWithMCEventPlane* mcep2 = new AliFlowAnalysisWithMCEventPlane();
  76     mcep2->SetHarmonic(2); // default is v2
  77     mcep2->Init();
  78     AliFlowAnalysisWithMCEventPlane* mcep3 = new AliFlowAnalysisWithMCEventPlane();
  79     mcep3->SetHarmonic(3); // default is v3
  80     mcep3->Init();
  81     AliFlowAnalysisWithMCEventPlane* mcep4 = new AliFlowAnalysisWithMCEventPlane();
  82     mcep4->SetHarmonic(4); // default is v4
  83     mcep4->Init();
  84     AliFlowAnalysisWithMCEventPlane* mcep5 = new AliFlowAnalysisWithMCEventPlane();
  85     mcep5->SetHarmonic(5); // default is v5
  86     mcep5->Init();
  87
  88     AliFlowAnalysisWithScalarProduct* sp2 = new AliFlowAnalysisWithScalarProduct();
  89     sp2->SetHarmonic(2);
  90     sp2->Init();
  91     AliFlowAnalysisWithScalarProduct* sp3 = new AliFlowAnalysisWithScalarProduct();
  92     sp3->SetHarmonic(3);
  93     sp3->Init();
  94     AliFlowAnalysisWithScalarProduct* sp4 = new AliFlowAnalysisWithScalarProduct();
  95     sp4->SetHarmonic(4);
  96     sp4->Init();
  97     AliFlowAnalysisWithScalarProduct* sp5 = new AliFlowAnalysisWithScalarProduct();
  98     sp5->SetHarmonic(5);
  99     sp5->Init();
 100
 101     AliFlowAnalysisWithQCumulants* qc2 = new AliFlowAnalysisWithQCumulants();
 102     qc2->SetHarmonic(2);
 103     qc2->Init();
 104     AliFlowAnalysisWithQCumulants* qc3 = new AliFlowAnalysisWithQCumulants();
 105     qc3->SetHarmonic(3);
 106     qc3->Init();
 107     AliFlowAnalysisWithQCumulants* qc4 = new AliFlowAnalysisWithQCumulants();
 108     qc4->SetHarmonic(4);
 109     qc4->Init();
 110     AliFlowAnalysisWithQCumulants* qc5 = new AliFlowAnalysisWithQCumulants();
 111     qc5->SetHarmonic(5);
 112     qc5->Init();
 113
 114     TH1F* histv2 = new TH1F(Form("v2 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{2}",100,0,0.5);
 115     TH1F* histv3 = new TH1F(Form("v3 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{3}",100,0,0.5);
 116     TH1F* histv4 = new TH1F(Form("v4 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{4}",100,0,0.5);
 117     TH1F* histv5 = new TH1F(Form("v5 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{5}",100,0,0.5);
 118
 119     AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
 120     cutsRP->SetPtMin(0.2);
 121     cutsRP->SetPtMax(5.0);
 122     AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
 123     cutsPOI->SetPtMin(0.2);
 124     cutsPOI->SetPtMax(10.0);
 125
 126     Printf("starting the main event loop..");
 127     for(Int_t i=0; i<nEvts; i++)
 128     {
 129         if (!glauber.NextEvent()) continue;
 130         Double_t mult = glauber.GetdNdEta();
 131         Double_t v2 = 0.2*glauber.GetEpsilon2Part(); histv2->Fill(v2);
 132         Double_t v3 = 0.18*glauber.GetEpsilon3Part(); histv3->Fill(v3);
 133         Double_t v4 = 0.09*glauber.GetEpsilon4Part(); histv4->Fill(v4);
 134         Double_t v5 = 0.03*glauber.GetEpsilon5Part(); histv5->Fill(v5);
 135
 136         AliFlowEventSimple* event = new AliFlowEventSimple(mult,AliFlowEventSimple::kGenerate);
 137         event->SetMCReactionPlaneAngle(0.);
 138         event->AddFlow(0,v2,v3,v4,v5);
 139         event->TagRP(cutsRP);
 140         event->TagPOI(cutsPOI);
 141
 142         // do flow analysis for various methods:
 143         mcep2->Make(event);
 144         mcep3->Make(event);
 145         mcep4->Make(event);
 146         mcep5->Make(event);
 147         sp2->Make(event);
 148         sp3->Make(event);
 149         sp4->Make(event);
 150         sp5->Make(event);
 151         qc2->Make(event);
 152         qc3->Make(event);
 153         qc4->Make(event);
 154         qc5->Make(event);
 155         cout <<"Event: " << i+1 << "\r"; cout.flush();
 156         delete event;
 157     } // end of for(Int_t i=0;i<nEvts;i++)
 158
 159     // calculate the final results
 160     mcep2->Finish();
 161     mcep3->Finish();
 162     mcep4->Finish();
 163     mcep5->Finish();
 164     sp2->Finish();
 165     sp3->Finish();
 166     sp4->Finish();
 167     sp5->Finish();
 168     qc2->Finish();
 169     qc3->Finish();
 170     qc4->Finish();
 171     qc5->Finish();
 172
 173     // open a new file which will hold the final results of all methods:
 174     const Int_t nMethods = 12;
 175     TString method[nMethods] = {"v2 MCEP","v3 MCEP","v4 MCEP","v5 MCEP","v2 SP","v3 SP","v4 SP","v5 SP","v2 QC","v3 QC","v4 QC","v5 QC"};
 176     TDirectoryFile *dirFileFinal[nMethods] = {NULL};
 177     TString fileName[nMethods];
 178     for(Int_t i=0; i<nMethods; i++)
 179     {
 180         // form a file name for each method:
 181         fileName[i]+=method[i].Data();
 182         fileName[i]+=Form(" %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]);
 183         dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
 184     }
 185
 186     // store the final results
 187     mcep2->WriteHistograms(dirFileFinal[0]);
 188     mcep3->WriteHistograms(dirFileFinal[1]);
 189     mcep4->WriteHistograms(dirFileFinal[2]);
 190     mcep5->WriteHistograms(dirFileFinal[3]);
 191     sp2->WriteHistograms(dirFileFinal[4]);
 192     sp3->WriteHistograms(dirFileFinal[5]);
 193     sp4->WriteHistograms(dirFileFinal[6]);
 194     sp5->WriteHistograms(dirFileFinal[7]);
 195     qc2->WriteHistograms(dirFileFinal[8]);
 196     qc3->WriteHistograms(dirFileFinal[9]);
 197     qc4->WriteHistograms(dirFileFinal[10]);
 198     qc5->WriteHistograms(dirFileFinal[11]);
 199
 200     histv2->Write();
 201     histv3->Write();
 202     histv4->Write();
 203     histv5->Write();
 204 }
 205