[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / runBalanceFunctionOnHydro.C

//========================================================//
//Macro that reads the output of the hydro calculations 
//(Therminator) using the format that can be found at
//arXiv:1102.0273
//Author: Panos.Christakoglou@nikhef.nl
//========================================================//

//========================================================//
//Event structure
struct StructEvent {
  UInt_t eventID;//unique identifier of the event (random number)
  UInt_t entries;//number of particles of each event
  UInt_t entriesprev;//set to 0 by default
};
//========================================================//

//========================================================//
//Particle structure
//Primordial particles: fatherid==-1
struct StructParticle {
  Float_t mass;//the mass of the particle (in GeV)
  Float_t t;//the time coordinate of the particle in fm/c
  Float_t x;//the spacial coordinate x in fm/c
  Float_t y;//the spacial coordinate y in fm/c
  Float_t z;//the spacial coordinate z in fm/c
  Float_t e;//the energy in GeV
  Float_t px;//the x-coordinate of the particle's momentum in GeV
  Float_t py;//the y-coordinate of the particle's momentum in GeV
  Float_t pz;//the z-coordinate of the particle's momentum in GeV
  Int_t decayed;//decay flag (1: decayed, 0: no decay)
  Int_t pid;//PDG of particle
  Int_t fatherpid;//PDG of parent
  Int_t rootpid;//root (primordial) particle PDG number
  Int_t eid;//sequence number in the event
  Int_t fathereid;//parent sequence number in the event
  UInt_t eventid;//unique identifier of the event (random number)
};
//========================================================//

//========================================================//
//Balance function analysis variables
Bool_t gRunShuffling=kFALSE;
Bool_t gRunMixing=kFALSE;
Bool_t gRunMixingWithEventPlane=kFALSE;

Double_t gEtaMin = -0.8;
Double_t gEtaMax = 0.8;
Double_t gPtMin = 0.2;
Double_t gPtMax = 20.0;
//========================================================//


void runBalanceFunctionOnHydro(TString aEventDir = "/glusterfs/alice1/alice2/pchrist/Therminator/pA/ChargeBalancing/Centrality1/lhyquid3v-LHCpPb5020s3.1Ti319t0.20Tf150e001/", Int_t aEventFiles = 2) {
  //Macro that reads the themrinator events
  //Author: Panos.Christakoglou@nikhef.nl
  // Time:
  TStopwatch timer;
  timer.Start();

  //========================================================//
  //Load the aliroot libraries
  gSystem->Load("libSTEERBase.so");
  gSystem->Load("libESD.so");
  gSystem->Load("libAOD.so");
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libEventMixing.so");
  gSystem->Load("libCORRFW.so");
  gSystem->Load("libPWGTools.so");
  gSystem->Load("libPWGCFebye.so");
  //========================================================//

  //========================================================//
  //Configure the bf objects
  AliBalancePsi *bf = new AliBalancePsi();
  bf->SetAnalysisLevel("MC");
  bf->SetShuffle(gRunShuffling);
  bf->SetEventClass("EventPlane");
  bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
  bf->InitHistograms();

  AliBalancePsi *bfs = 0x0;
  if(gRunShuffling) {
    bfs = new AliBalancePsi();
    bfs->SetAnalysisLevel("MC");
    bfs->SetEventClass("EventPlane");
    bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
    bfs->InitHistograms();
  }

  AliBalancePsi *bfm = 0x0;
  if(gRunMixing) {
    bfm = new AliBalancePsi();
    bfm->SetAnalysisLevel("MC");
    bfm->SetShuffle(gRunShuffling);
    bfm->SetEventClass("EventPlane");
    bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
    bfm->InitHistograms();
  }
  //========================================================//

  //========================================================//
  //Output objects
  //QA list
  fList = new TList();
  fList->SetName("listQA");
  fList->SetOwner();

  //Balance Function list
  TList *fListBF = new TList();
  fListBF->SetName("listBF");
  fListBF->SetOwner();
  
  //Balance function list: shuffling
  TList *fListBFS = 0x0;
  if(gRunShuffling) {
    fListBFS = new TList();
    fListBFS->SetName("listBFShuffled");
    fListBFS->SetOwner();
    fListBFS->Add(bfs->GetHistNp());
    fListBFS->Add(bfs->GetHistNn());
    fListBFS->Add(bfs->GetHistNpn());
    fListBFS->Add(bfs->GetHistNnn());
    fListBFS->Add(bfs->GetHistNpp());
    fListBFS->Add(bfs->GetHistNnp());
  }

  //Balance function list: event mixing
  TList *fListBFM = 0x0;
  if(gRunMixing) {
    fListBFM = new TList();
    fListBFM->SetName("listBFMixed");
    fListBFM->SetOwner();
    fListBFM->Add(bfm->GetHistNp());
    fListBFM->Add(bfm->GetHistNn());
    fListBFM->Add(bfm->GetHistNpn());
    fListBFM->Add(bfm->GetHistNnn());
    fListBFM->Add(bfm->GetHistNpp());
    fListBFM->Add(bfm->GetHistNnp());
  }
  //========================================================//

  //========================================================//
  //Event Mixing
  if(gRunMixing){
    Int_t trackDepth = 50000;
    Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
    
    // centrality bins
    Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
    Double_t* centbins        = centralityBins;
    Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
    
    // Zvtx bins
    Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
    Double_t* vtxbins     = vertexBins;
    Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
    
    // Event plane angle (Psi) bins
    Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
    Double_t* psibins     = psiBins;
    Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;

    AliEventPoolManager *fPoolMgr = 0x0;
    // run the event mixing also in bins of event plane (statistics!)
    if(gRunMixingEventPlane){
      fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
    }
    else{
      fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
    }
  }
  //========================================================//

  //========================================================//
  //Create the TChain object: events
  TChain *eventChain = new TChain("events");
  StructEvent tStructEvents;
  eventChain->SetBranchAddress("event",&tStructEvents);
  //========================================================//

  //========================================================//
  //Create the TChain object: particles 
  TChain *particleChain = new TChain("particles");
  StructParticle tStructParticles;
  particleChain->SetBranchAddress("particle",&tStructParticles);
  //========================================================//

  //========================================================//
  //Fill the TChain with the files in the directory
  for(Int_t iFile = 1; iFile <= 9; iFile++) {
    TString filename = aEventDir.Data();
    filename += "/event00"; filename += iFile;
    filename += ".root";
    cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
    eventChain->Add(filename.Data());
    particleChain->Add(filename.Data());
  }
  //========================================================//
 
  //========================================================//
  //Histograms
  //Event stats.
  TString gCutName[5] = {"Total","Offline trigger",
                         "Vertex","Analyzed","sel. Centrality"};
  TH2F *fHistEventStats = new TH2F("fHistEventStats",
				   "Event statistics;;Centrality percentile;N_{events}",
				   5,0.5,5.5,220,-5,105);
  for(Int_t i = 1; i <= 5; i++)
    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
  fList->Add(fHistEventStats);

  //Number of accepted particles
  TH2F *fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
  fList->Add(fHistNumberOfAcceptedTracks);
  //========================================================//

  //========================================================//
  //loop over the events
  Int_t nEvents = eventChain->GetEntries();
  cout<<"========================================="<<endl;
  cout<<"Number of events in the chain: "<<nEvents<<endl;
  cout<<"========================================="<<endl;
  Int_t iParticleCounter = 0;
  Int_t nTotalParticles = 0;
  
  //for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
  for(Int_t  // for local changed BF configuration
  //gROOT->LoadMacro("./configBalanceFunctionPsiAnalysis.C");
 iEvent = 0; iEvent < 1; iEvent++) {
    eventChain->GetEntry(iEvent);

    Int_t nParticles = tStructEvents.entries;
    cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<" - Entries(prev.): "<<tStructEvents.entriesprev<<endl;

    //========================================================//
    //loop over particles
    for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
      particleChain->GetEntry(nTotalParticles+iParticle);
      
      Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + TMath::Power(tStructParticles.py,2));
      hPt->Fill(gPt);
      	
      iParticleCounter += 1;
      cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl;
    }//particle loop
    nTotalParticles += nParticles;
  }

  //========================================================//
  //Output file
  TFile *f = TFile::Open("therminator.root","recreate");
  hPdgCode->Write();
  hEta->Write();
  hPt->Write();
  f->Close();  
  //========================================================//

  // Print real and CPU time used for analysis:  
  timer.Stop();
  timer.Print();
}
Commit	Line	Data
ee97f33d	1	//========================================================//
	2	//Macro that reads the output of the hydro calculations
	3	//(Therminator) using the format that can be found at
	4	//arXiv:1102.0273
	5	//Author: Panos.Christakoglou@nikhef.nl
	6	//========================================================//
	7
	8	//========================================================//
	9	//Event structure
	10	struct StructEvent {
	11	UInt_t eventID;//unique identifier of the event (random number)
	12	UInt_t entries;//number of particles of each event
	13	UInt_t entriesprev;//set to 0 by default
	14	};
	15	//========================================================//
	16
	17	//========================================================//
	18	//Particle structure
	19	//Primordial particles: fatherid==-1
	20	struct StructParticle {
	21	Float_t mass;//the mass of the particle (in GeV)
	22	Float_t t;//the time coordinate of the particle in fm/c
	23	Float_t x;//the spacial coordinate x in fm/c
	24	Float_t y;//the spacial coordinate y in fm/c
	25	Float_t z;//the spacial coordinate z in fm/c
	26	Float_t e;//the energy in GeV
	27	Float_t px;//the x-coordinate of the particle's momentum in GeV
	28	Float_t py;//the y-coordinate of the particle's momentum in GeV
	29	Float_t pz;//the z-coordinate of the particle's momentum in GeV
	30	Int_t decayed;//decay flag (1: decayed, 0: no decay)
	31	Int_t pid;//PDG of particle
	32	Int_t fatherpid;//PDG of parent
	33	Int_t rootpid;//root (primordial) particle PDG number
	34	Int_t eid;//sequence number in the event
	35	Int_t fathereid;//parent sequence number in the event
	36	UInt_t eventid;//unique identifier of the event (random number)
	37	};
	38	//========================================================//
	39
	40	//========================================================//
	41	//Balance function analysis variables
	42	Bool_t gRunShuffling=kFALSE;
	43	Bool_t gRunMixing=kFALSE;
	44	Bool_t gRunMixingWithEventPlane=kFALSE;
	45
	46	Double_t gEtaMin = -0.8;
	47	Double_t gEtaMax = 0.8;
	48	Double_t gPtMin = 0.2;
	49	Double_t gPtMax = 20.0;
	50	//========================================================//
	51
	52
	53	void runBalanceFunctionOnHydro(TString aEventDir = "/glusterfs/alice1/alice2/pchrist/Therminator/pA/ChargeBalancing/Centrality1/lhyquid3v-LHCpPb5020s3.1Ti319t0.20Tf150e001/", Int_t aEventFiles = 2) {
	54	//Macro that reads the themrinator events
	55	//Author: Panos.Christakoglou@nikhef.nl
	56	// Time:
	57	TStopwatch timer;
	58	timer.Start();
	59
	60	//========================================================//
	61	//Load the aliroot libraries
	62	gSystem->Load("libSTEERBase.so");
	63	gSystem->Load("libESD.so");
	64	gSystem->Load("libAOD.so");
65	gSystem->Load("libANALYSIS.so");
66	gSystem->Load("libANALYSISalice.so");
67	gSystem->Load("libEventMixing.so");
68	gSystem->Load("libCORRFW.so");
69	gSystem->Load("libPWGTools.so");
70	gSystem->Load("libPWGCFebye.so");
71	//========================================================//
72
73	//========================================================//
74	//Configure the bf objects
75	AliBalancePsi *bf = new AliBalancePsi();
76	bf->SetAnalysisLevel("MC");
77	bf->SetShuffle(gRunShuffling);
78	bf->SetEventClass("EventPlane");
79	bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
80	bf->InitHistograms();
81
82	AliBalancePsi *bfs = 0x0;
83	if(gRunShuffling) {
84	bfs = new AliBalancePsi();
85	bfs->SetAnalysisLevel("MC");
86	bfs->SetEventClass("EventPlane");
87	bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
88	bfs->InitHistograms();
89	}
90
91	AliBalancePsi *bfm = 0x0;
92	if(gRunMixing) {
93	bfm = new AliBalancePsi();
94	bfm->SetAnalysisLevel("MC");
95	bfm->SetShuffle(gRunShuffling);
96	bfm->SetEventClass("EventPlane");
97	bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
98	bfm->InitHistograms();
99	}
100	//========================================================//
101
102	//========================================================//
103	//Output objects
104	//QA list
105	fList = new TList();
106	fList->SetName("listQA");
107	fList->SetOwner();
108
109	//Balance Function list
110	TList *fListBF = new TList();
111	fListBF->SetName("listBF");
112	fListBF->SetOwner();
113
114	//Balance function list: shuffling
115	TList *fListBFS = 0x0;
116	if(gRunShuffling) {
117	fListBFS = new TList();
118	fListBFS->SetName("listBFShuffled");
119	fListBFS->SetOwner();
120	fListBFS->Add(bfs->GetHistNp());
121	fListBFS->Add(bfs->GetHistNn());
122	fListBFS->Add(bfs->GetHistNpn());
123	fListBFS->Add(bfs->GetHistNnn());
124	fListBFS->Add(bfs->GetHistNpp());
125	fListBFS->Add(bfs->GetHistNnp());
126	}
127
128	//Balance function list: event mixing
129	TList *fListBFM = 0x0;
130	if(gRunMixing) {
131	fListBFM = new TList();
132	fListBFM->SetName("listBFMixed");
133	fListBFM->SetOwner();
134	fListBFM->Add(bfm->GetHistNp());
135	fListBFM->Add(bfm->GetHistNn());
136	fListBFM->Add(bfm->GetHistNpn());
137	fListBFM->Add(bfm->GetHistNnn());
138	fListBFM->Add(bfm->GetHistNpp());
139	fListBFM->Add(bfm->GetHistNnp());
140	}
141	//========================================================//
142
143	//========================================================//
144	//Event Mixing
145	if(gRunMixing){
146	Int_t trackDepth = 50000;
147	Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
148
149	// centrality bins
150	Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
151	Double_t* centbins = centralityBins;
152	Int_t nCentralityBins = sizeof(centralityBins) / sizeof(Double_t) - 1;
153
154	// Zvtx bins
155	Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
156	Double_t* vtxbins = vertexBins;
157	Int_t nVertexBins = sizeof(vertexBins) / sizeof(Double_t) - 1;
158
159	// Event plane angle (Psi) bins
160	Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
161	Double_t* psibins = psiBins;
162	Int_t nPsiBins = sizeof(psiBins) / sizeof(Double_t) - 1;
163
164	AliEventPoolManager *fPoolMgr = 0x0;
165	// run the event mixing also in bins of event plane (statistics!)
166	if(gRunMixingEventPlane){
167	fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
168	}
169	else{
170	fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
171	}
172	}
173	//========================================================//
174
175	//========================================================//
176	//Create the TChain object: events
177	TChain *eventChain = new TChain("events");
178	StructEvent tStructEvents;
179	eventChain->SetBranchAddress("event",&tStructEvents);
180	//========================================================//
181
182	//========================================================//
183	//Create the TChain object: particles
184	TChain *particleChain = new TChain("particles");
185	StructParticle tStructParticles;
186	particleChain->SetBranchAddress("particle",&tStructParticles);
187	//========================================================//
188
189	//========================================================//
190	//Fill the TChain with the files in the directory
191	for(Int_t iFile = 1; iFile <= 9; iFile++) {
192	TString filename = aEventDir.Data();
193	filename += "/event00"; filename += iFile;
194	filename += ".root";
195	cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
196	eventChain->Add(filename.Data());
197	particleChain->Add(filename.Data());
198	}
199	//========================================================//
200
201	//========================================================//
202	//Histograms
203	//Event stats.
204	TString gCutName[5] = {"Total","Offline trigger",
205	"Vertex","Analyzed","sel. Centrality"};
206	TH2F *fHistEventStats = new TH2F("fHistEventStats",
207	"Event statistics;;Centrality percentile;N_{events}",
208	5,0.5,5.5,220,-5,105);
209	for(Int_t i = 1; i <= 5; i++)
210	fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
211	fList->Add(fHistEventStats);
212
213	//Number of accepted particles
214	TH2F *fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
215	fList->Add(fHistNumberOfAcceptedTracks);
216	//========================================================//
217
218	//========================================================//
219	//loop over the events
220	Int_t nEvents = eventChain->GetEntries();
221	cout<<"========================================="<<endl;
222	cout<<"Number of events in the chain: "<<nEvents<<endl;
223	cout<<"========================================="<<endl;
224	Int_t iParticleCounter = 0;
225	Int_t nTotalParticles = 0;
226
227	//for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
228	for(Int_t // for local changed BF configuration
229	//gROOT->LoadMacro("./configBalanceFunctionPsiAnalysis.C");
230	iEvent = 0; iEvent < 1; iEvent++) {
231	eventChain->GetEntry(iEvent);
232
233	Int_t nParticles = tStructEvents.entries;
234	cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<" - Entries(prev.): "<<tStructEvents.entriesprev<<endl;
235
236	//========================================================//
237	//loop over particles
238	for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
239	particleChain->GetEntry(nTotalParticles+iParticle);
240
241	Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) + TMath::Power(tStructParticles.py,2));
242	hPt->Fill(gPt);
243
244	iParticleCounter += 1;
245	cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl;
246	}//particle loop
247	nTotalParticles += nParticles;
248	}
249
250	//========================================================//
251	//Output file
252	TFile *f = TFile::Open("therminator.root","recreate");
253	hPdgCode->Write();
254	hEta->Write();
255	hPt->Write();
256	f->Close();
257	//========================================================//
258
259	// Print real and CPU time used for analysis:
260	timer.Stop();
261	timer.Print();
262	}