From 22f32e5d3d14b12ca6c12c201172fd82ac1118e3 Mon Sep 17 00:00:00 2001
From: snelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Tue, 22 Mar 2011 18:27:05 +0000
Subject: [PATCH] added some example cases

---
 .../examples/runFlowOnTheFlyGlauber.C         | 205 ++++++++++++++++++
 ...{runGlauberMCexample2.C => runGlauberMC.C} |  30 ++-
 2 files changed, 223 insertions(+), 12 deletions(-)
 create mode 100644 PWG2/FLOW/Documentation/examples/runFlowOnTheFlyGlauber.C
 rename PWG2/FLOW/Documentation/examples/{runGlauberMCexample2.C => runGlauberMC.C} (60%)

diff --git a/PWG2/FLOW/Documentation/examples/runFlowOnTheFlyGlauber.C b/PWG2/FLOW/Documentation/examples/runFlowOnTheFlyGlauber.C
new file mode 100644
index 00000000000..3dbb5cd9276
--- /dev/null
+++ b/PWG2/FLOW/Documentation/examples/runFlowOnTheFlyGlauber.C
@@ -0,0 +1,205 @@
+#include "TStopwatch.h"
+#include "Riostream.h"
+#include "TFile.h"
+
+int runFlowOnTheFlyGlauber(Int_t nevents=100000, Int_t iseed=1)
+{
+    gSystem->Load("libGeom");
+    gSystem->Load("libVMC");
+    gSystem->Load("libXMLIO");
+    gSystem->Load("libPhysics");
+    gSystem->Load("libPWG2flowCommon");
+    gSystem->Load("libPWG2flowTools");
+    
+    fMyTRandom3 = new TRandom3(iseed);   
+    gRandom->SetSeed(iseed);
+    TFile *outputFile = new TFile("flowOnTheFlyGlauber.root","RECREATE");
+    
+    TStopwatch timer;
+    timer.Start();
+    
+    runOneCentrality(3,nevents);
+    runOneCentrality(4,nevents);
+    runOneCentrality(5,nevents);
+    
+    outputFile->Close();
+    delete outputFile;
+    
+    cout<<endl; cout<<" ---- Fini ---- "<<endl; cout<<endl;
+    timer.Stop(); timer.Print();
+}
+
+void runOneCentrality(Int_t centrality=3, Int_t nEvts=100)
+{
+    
+    //centralities: 0: 0-5
+    //              1: 5-10
+    //              2: 10-20
+    //              3: 20-30
+    //              4: 30-40
+    //              5: 40-50
+    //              6: 50-60
+    //              7: 60-70
+    //              8: 70-80
+    
+    Double_t centPerc[] = {0,5,10,20,30,40,50,60,70,80};
+    Double_t xCross[] = {2.5,7.5,15.,25.,35.,45.,55.,65.,75.};
+    Double_t xCrossErr[9] = {0.};
+    Double_t impact[] = {0., 3.5, 4.95, 6.98, 8.55, 9.88, 11.4, 12.09, 13.06, 13.97};
+    Double_t npart[] = {1000, 356, 305, 220, 155, 105, 67, 40, 21, 10, 4, 0};
+ 
+    Option_t *sysA="Pb"; 
+    Option_t *sysB="Pb";
+    Double_t signn=64; // inelastic nucleon nucleon cross section
+    //const char *fname="GlauberMC_PbPb_ntuple.root"; // name output file
+    //  AliGlauberMC::RunAndSaveNtuple(nevents,sysA,sysB,signn);
+    // run the code to produce an ntuple:
+
+    Double_t mind=0.4;
+    Double_t r=6.62;
+    Double_t a=0.546;
+    
+    AliGlauberMC glauber(sysA,sysB,signn);
+    glauber.SetMinDistance(mind);
+    glauber.Setr(r);
+    glauber.Seta(a);
+    glauber.SetDoPartProduction(kTRUE);
+    glauber.SetBmin(impact[centrality]);
+    glauber.SetBmax(impact[centrality+1]);    
+    glauber.SetdNdEtaType(AliGlauberMC::kNBDSV);
+    glauber.GetdNdEtaParam()[0] = 2.49;    //npp
+    glauber.GetdNdEtaParam()[1] = 1.7;  //ratioSgm2Mu
+    glauber.GetdNdEtaParam()[2] = 0.13; //xhard
+
+    
+    AliFlowAnalysisWithMCEventPlane* mcep2 = new AliFlowAnalysisWithMCEventPlane();
+    mcep2->SetHarmonic(2); // default is v2
+    mcep2->Init();
+    AliFlowAnalysisWithMCEventPlane* mcep3 = new AliFlowAnalysisWithMCEventPlane();
+    mcep3->SetHarmonic(3); // default is v3
+    mcep3->Init();
+    AliFlowAnalysisWithMCEventPlane* mcep4 = new AliFlowAnalysisWithMCEventPlane();
+    mcep4->SetHarmonic(4); // default is v4
+    mcep4->Init();
+    AliFlowAnalysisWithMCEventPlane* mcep5 = new AliFlowAnalysisWithMCEventPlane();
+    mcep5->SetHarmonic(5); // default is v5
+    mcep5->Init();
+    
+    AliFlowAnalysisWithScalarProduct* sp2 = new AliFlowAnalysisWithScalarProduct();
+    sp2->SetHarmonic(2);
+    sp2->Init();
+    AliFlowAnalysisWithScalarProduct* sp3 = new AliFlowAnalysisWithScalarProduct();
+    sp3->SetHarmonic(3);
+    sp3->Init();
+    AliFlowAnalysisWithScalarProduct* sp4 = new AliFlowAnalysisWithScalarProduct();
+    sp4->SetHarmonic(4);
+    sp4->Init();
+    AliFlowAnalysisWithScalarProduct* sp5 = new AliFlowAnalysisWithScalarProduct();
+    sp5->SetHarmonic(5);
+    sp5->Init();
+    
+    AliFlowAnalysisWithQCumulants* qc2 = new AliFlowAnalysisWithQCumulants();
+    qc2->SetHarmonic(2);
+    qc2->Init();
+    AliFlowAnalysisWithQCumulants* qc3 = new AliFlowAnalysisWithQCumulants();
+    qc3->SetHarmonic(3);
+    qc3->Init();
+    AliFlowAnalysisWithQCumulants* qc4 = new AliFlowAnalysisWithQCumulants();
+    qc4->SetHarmonic(4);
+    qc4->Init();
+    AliFlowAnalysisWithQCumulants* qc5 = new AliFlowAnalysisWithQCumulants();
+    qc5->SetHarmonic(5);
+    qc5->Init();
+    
+    TH1F* histv2 = new TH1F(Form("v2 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{2}",100,0,0.5);
+    TH1F* histv3 = new TH1F(Form("v3 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{3}",100,0,0.5);
+    TH1F* histv4 = new TH1F(Form("v4 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{4}",100,0,0.5);
+    TH1F* histv5 = new TH1F(Form("v5 spread %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]),"v_{5}",100,0,0.5);
+    
+    AliFlowTrackSimpleCuts *cutsRP = new AliFlowTrackSimpleCuts();
+    cutsRP->SetPtMin(0.2);
+    cutsRP->SetPtMax(5.0);
+    AliFlowTrackSimpleCuts *cutsPOI = new AliFlowTrackSimpleCuts();
+    cutsPOI->SetPtMin(0.2);
+    cutsPOI->SetPtMax(10.0);
+    
+    Printf("starting the main event loop..");
+    for(Int_t i=0; i<nEvts; i++)
+    {
+        if (!glauber.NextEvent()) continue;
+        Double_t mult = glauber.GetdNdEta();
+        Double_t v2 = 0.2*glauber.GetEpsilon2Part(); histv2->Fill(v2);
+        Double_t v3 = 0.18*glauber.GetEpsilon3Part(); histv3->Fill(v3);
+        Double_t v4 = 0.09*glauber.GetEpsilon4Part(); histv4->Fill(v4);
+        Double_t v5 = 0.03*glauber.GetEpsilon5Part(); histv5->Fill(v5);
+        
+        AliFlowEventSimple* event = new AliFlowEventSimple(mult,AliFlowEventSimple::kGenerate);
+        event->SetMCReactionPlaneAngle(0.);
+        event->AddFlow(0,v2,v3,v4,v5);
+        event->TagRP(cutsRP);
+        event->TagPOI(cutsPOI);
+        
+        // do flow analysis for various methods:
+        mcep2->Make(event);
+        mcep3->Make(event);
+        mcep4->Make(event);
+        mcep5->Make(event);
+        sp2->Make(event);
+        sp3->Make(event);
+        sp4->Make(event);
+        sp5->Make(event);
+        qc2->Make(event);
+        qc3->Make(event);
+        qc4->Make(event);
+        qc5->Make(event);
+        cout <<"Event: " << i+1 << "\r"; cout.flush();
+        delete event;
+    } // end of for(Int_t i=0;i<nEvts;i++)
+    
+    // calculate the final results
+    mcep2->Finish();
+    mcep3->Finish();
+    mcep4->Finish();
+    mcep5->Finish();
+    sp2->Finish();
+    sp3->Finish();
+    sp4->Finish();
+    sp5->Finish();
+    qc2->Finish();
+    qc3->Finish();
+    qc4->Finish();
+    qc5->Finish();
+    
+    // open a new file which will hold the final results of all methods:
+    const Int_t nMethods = 12;
+    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"};
+    TDirectoryFile *dirFileFinal[nMethods] = {NULL};
+    TString fileName[nMethods];
+    for(Int_t i=0; i<nMethods; i++)
+    {
+        // form a file name for each method:
+        fileName[i]+=method[i].Data();
+        fileName[i]+=Form(" %.0f-%.0f",centPerc[centrality],centPerc[centrality+1]);
+        dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
+    }
+    
+    // store the final results
+    mcep2->WriteHistograms(dirFileFinal[0]);
+    mcep3->WriteHistograms(dirFileFinal[1]);
+    mcep4->WriteHistograms(dirFileFinal[2]);
+    mcep5->WriteHistograms(dirFileFinal[3]);
+    sp2->WriteHistograms(dirFileFinal[4]);
+    sp3->WriteHistograms(dirFileFinal[5]);
+    sp4->WriteHistograms(dirFileFinal[6]);
+    sp5->WriteHistograms(dirFileFinal[7]);
+    qc2->WriteHistograms(dirFileFinal[8]);
+    qc3->WriteHistograms(dirFileFinal[9]);
+    qc4->WriteHistograms(dirFileFinal[10]);
+    qc5->WriteHistograms(dirFileFinal[11]);
+    
+    histv2->Write();
+    histv3->Write();
+    histv4->Write();
+    histv5->Write();
+}
+
diff --git a/PWG2/FLOW/Documentation/examples/runGlauberMCexample2.C b/PWG2/FLOW/Documentation/examples/runGlauberMC.C
similarity index 60%
rename from PWG2/FLOW/Documentation/examples/runGlauberMCexample2.C
rename to PWG2/FLOW/Documentation/examples/runGlauberMC.C
index 861fa41080b..42675929dc5 100644
--- a/PWG2/FLOW/Documentation/examples/runGlauberMCexample2.C
+++ b/PWG2/FLOW/Documentation/examples/runGlauberMC.C
@@ -1,17 +1,17 @@
+void runGlauberMC()
 {
   //load libraries
- gSystem->Load("libPWG2flowTools");
- //  gSystem->SetBuildDir("/tmp");
- //  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberNucleon.cxx+");
- //  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberNucleus.cxx+");
- //  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberMC.cxx+");
+  gSystem->Load("libVMC");
+  gSystem->Load("libPhysics");
+  gSystem->Load("libTree");
+  gSystem->Load("libPWG2flowTools");
 
   //set the random seed from current time
   TTimeStamp time;
-  Int_t seed = time->GetSec();
+  Int_t seed = time.GetSec();
   gRandom->SetSeed(seed);
 
-  Int_t nevents = 1000; // number of events to simulate 
+  Int_t nevents = 1000000; // number of events to simulate 
   // supported systems are e.g. "p", "d", "Si", "Au", "Pb", "U" 
   Option_t *sysA="Pb"; 
   Option_t *sysB="Pb";
@@ -30,14 +30,20 @@
   mcg.SetMinDistance(mind);
   mcg.Setr(r);
   mcg.Seta(a);
-  mcg.SetDoPartProduction(kFALSE);
-  //mcg.SetDoPartProduction(kTRUE);
+  mcg.SetDoPartProduction(kTRUE);
+  
+  //////////////////
+  mcg.SetdNdEtaType(AliGlauberMC::kNBDSV);
+  mcg.GetdNdEtaParam()[0] = 2.49;    //npp
+  mcg.GetdNdEtaParam()[1] = 1.7;  //ratioSgm2Mu
+  mcg.GetdNdEtaParam()[2] = 0.13; //xhard
+  //////////////////
+
   mcg.Run(nevents);
 
-  TNtuple  *nt=mcg.GetNtuple();
+  TNtuple  *nt = mcg.GetNtuple();
   TFile out(fname,"recreate",fname,9);
   if(nt) nt->Write();
-  printf("total cross section with a nucleon-nucleon cross section \t%f is \t%f",signn,mcg.GetTotXSect());
+  printf("total cross section with a nucleon-nucleon cross section %.4f is %.4f\n\n",signn,mcg.GetTotXSect());
   out.Close();
-
 }
-- 
2.39.3