[u/mrichter/AliRoot.git] / PWGCF / FLOW / Documentation / examples / runFlowOnTheFlyGlauber.C

#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("libPWGflowBase");
    gSystem->Load("libPWGGlauber");
    
    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();
}
Commit	Line	Data
22f32e5d	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");
2e311896	11	gSystem->Load("libPWGflowBase");
2e311896	12	gSystem->Load("libPWGGlauber");
22f32e5d	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