[u/mrichter/AliRoot.git] / PWGPP / EVCHAR / GlauberSNM / plotGlauberCenVarsZPA.C

///==========================================================================
///
///    macro to plot centrality bin values 
///==========================================================================
///
#include <cstdlib>

const Int_t nbins;

double centPercent[]={5.,10.,20.,40.,60.,80.,99.9999999}; 
//double centPercent[]={10.,20.,40.,60.,80.,99.9999999}; 


nbins = sizeof(centPercent)/sizeof(double);

TArrayI* binUp = new TArrayI(nbins);
TArrayF* Multbin = new TArrayF(nbins);

void plotGlauberCenVarsZPA()    
{
 TGraphErrors *gNpart=new TGraphErrors(0);
 gNpart->SetName("gNpart"); 
 TGraphErrors *gNcoll=new TGraphErrors(0);
 gNcoll->SetName("gNcoll"); 
 TGraphErrors *gtAA=new TGraphErrors(0);
 gtAA->SetName("gtAA"); 


  TFile *f=new TFile("ZPA_ntuple_188359.root");
  
  TNtuple *nt=(TNtuple*)f->Get("gnt");
  
  Float_t B;
  Float_t Npart;
  Float_t Ncoll;
  Float_t tAA;
  Float_t ntot;
  
  nt->SetBranchAddress("B",&B);
  nt->SetBranchAddress("Npart",&Npart);
  nt->SetBranchAddress("Ncoll",&Ncoll);
  nt->SetBranchAddress("tAA",&tAA);
  nt->SetBranchAddress("ntot",&ntot);


  Int_t totev=nt->GetEntries();
 
 
  TH1F*hB=new TH1F("hB","hB",28000,0.,14.);
  TH1F*hNpart=new TH1F("hNpart","hNpart",420,-0.5,419.5);
  TH1F*hNcoll=new TH1F("hNcoll","hNcoll",3500,0,35);
  TH1F*htAA=new TH1F("htAA","htAA",100000,-0.5,39.5);
  TH1F*hntot=new TH1F("hntot","hntot",100000,-0.5,39.5);
 
  
  for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
    nt->GetEvent(iEvent); 
    
    Float_t nu = (Float_t) ntot; 
    nu = gRandom->Gaus(nu,0.5); 

    if(nu>0) {
  
    hNpart->Fill(Npart);
    hNcoll->Fill(Ncoll);
    hB->Fill(B);
    htAA->Fill(tAA);
    hntot->Fill(nu);
    }
  } 

   new TCanvas();
   hNcoll->Draw();

   //---------------------------------------------------
   getCentrality(hntot);
   //---------------------------------------------------


 TH1F* hnpartcutb[nbins];
 char histtitp[100];
 for (Int_t i=0; i<binUp->GetSize(); i++) {
   sprintf(histtitp,"npartcutb%d",i);
   hnpartcutb[i] = new TH1F(histtitp,histtitp,10000,-0.5,9999.5);
   hnpartcutb[i]->SetLineWidth(1);
   hnpartcutb[i]->SetStats(1);

 }

 TH1F* hncollcutb[nbins];
 char histtitc[100];
 for (Int_t i=0; i<binUp->GetSize(); i++) {
   sprintf(histtitc,"ncollcutb%d",i);
   hncollcutb[i] = new TH1F(histtitc,histtitc,10000,-0.5,9999.5);
   hncollcutb[i]->SetLineWidth(1);
   hncollcutb[i]->SetStats(1);

 }

 TH1F* htaacutb[nbins];
 char histtitt[100];
 for (Int_t i=0; i<binUp->GetSize(); i++) {
   sprintf(histtitt,"taacutb%d",i);
   htaacutb[i] = new TH1F(histtitt,histtitt,100000,-0.5,39.5);
   htaacutb[i]->SetLineWidth(1);
   htaacutb[i]->SetStats(1);

 }

 TH1F* hbcutb[nbins];
 char histtitt[100];
 for (Int_t i=0; i<binUp->GetSize(); i++) {
   sprintf(histtitt,"bcutb%d",i);
   hbcutb[i] = new TH1F(histtitt,histtitt,100000,-0.5,39.5);
   hbcutb[i]->SetLineWidth(1);
   hbcutb[i]->SetStats(1);

 }


 for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
   nt->GetEvent(iEvent);
   for (int ibin=0; ibin<nbins; ibin++) {
     //     if (B<Multbin->At(ibin)) {

     Float_t nu = (Float_t) Ncoll; 
     nu = nu+1.*gRandom->Rndm();
  
     if (nu>Multbin->At(ibin)) {
       hnpartcutb[ibin]->Fill(Npart); 
       hncollcutb[ibin]->Fill(Ncoll); 
       htaacutb[ibin]->Fill(tAA); 
       hbcutb[ibin]->Fill(B); 
       break;
     } 
   } 
 } 
 

 // superimpose cut histograms
// new TCanvas();
// hnpartcutb[0]->Draw("");
 for (Int_t icentr=0; icentr<nbins;icentr++) {
   hnpartcutb[icentr]->SetLineColor(icentr);
   hnpartcutb[icentr]->Draw("same");
 } 

 for (Int_t icentr=0; icentr<nbins;icentr++) {
   new TCanvas();
   hnpartcutb[icentr]->SetLineColor(icentr);
   hnpartcutb[icentr]->SetStats(1); 
   hnpartcutb[icentr]->Draw("");
   gNpart->SetPoint(icentr,Float_t(icentr),hnpartcutb[icentr]->GetMean());
   gNpart->SetPointError(icentr,0,hnpartcutb[icentr]->GetRMS());
 }
 cout << endl;

 for (Int_t icentr=0; icentr<nbins;icentr++) {
   new TCanvas();
   hncollcutb[icentr]->SetLineColor(icentr);
   hncollcutb[icentr]->SetStats(1); 
   hncollcutb[icentr]->Draw("");
   gNcoll->SetPoint(icentr,Float_t(icentr),hncollcutb[icentr]->GetMean());
   gNcoll->SetPointError(icentr,0,hncollcutb[icentr]->GetRMS());
 }
 cout << endl;

 for (Int_t icentr=0; icentr<nbins;icentr++) {
   new TCanvas();
   htaacutb[icentr]->SetLineColor(icentr);
   htaacutb[icentr]->SetStats(1); 
   htaacutb[icentr]->Draw("");
   gtAA->SetPoint(icentr,Float_t(icentr),htaacutb[icentr]->GetMean());
   gtAA->SetPointError(icentr,0,htaacutb[icentr]->GetRMS());
 }

 for (Int_t icentr=0; icentr<nbins;icentr++) 
   cout<<icentr<<" | "<<setprecision(3)<<centPercent[icentr]<<" | "
       <<Multbin->At(icentr-1)<<" | "<<Multbin->At(icentr)<<" | "
       <<hnpartcutb[icentr]->GetMean()<<" | " <<hnpartcutb[icentr]->GetRMS()<< " | " 
       <<hncollcutb[icentr]->GetMean()<<" | " <<hncollcutb[icentr]->GetRMS()<< " | " 
       <<hbcutb[icentr]->GetMean()    <<" | " <<hbcutb[icentr]->GetRMS()    << " | "
       <<htaacutb[icentr]->GetMean()  <<" | " <<htaacutb[icentr]->GetRMS()  << " | "<<endl;


 //TString suffixhisto=Form("/home/atoia/GlauberNtuple/62/GlauberMC_PbPb_histoB_sigma%d_mind%d_r%d_a%d_%d.root",mysignn,mymind,myr,mya,nbins);
 //TString suffixhisto=Form("/home/atoia/GlauberNtuple/GlauberMC_PbPb_histoB_sigma%d_mind%d_r%d_a%d_%d.root",mysignn,mymind,myr,mya,nbins);
 //const Char_t* filehistoname=suffixhisto.Data();
 //TFile*filefinal=new TFile(filehistoname,"recreate");
 // --FOR TEST --
 //TFile*filefinal=new TFile("GlauberInfo.root","recreate");
 
 // gNpart->Write();
 // gNcoll->Write();
 // gtAA->Write();
 
 // hNpart->Write();
 // hNcoll->Write();
 // hB    ->Write();
 // htAA  ->Write();
  
 // for (int i=0; i<nbins; i++) {
 //   hnpartcutb[i]->Write();
 //   hncollcutb[i]->Write();
 //   htaacutb[i]  ->Write();
 // }
 
}


void getCentrality(TH1 *histNch, Float_t ff=1.0)
// histNch - histo of multiplicity distribution (better with binsize=1)x
// ff fraction of accepted events. All losses are assumed to occur in most
// peripheral bin
{

 //double sum= histNch->GetEntries() - histNch->GetBinContent(1);
 double sum= histNch->Integral(); 
 int nbinsx=histNch->GetNbinsX();
 double frac=0.;
 int ic=0;
 //for (int ib=1;ib<=nbinsx;ib++){
 for (int ib=nbinsx;ib>0;ib--){
   frac += histNch->GetBinContent(ib)/sum*100.*ff;
   if(frac > centPercent[ic]){
     binUp->SetAt(ib,ic);
     Multbin->SetAt(histNch->GetBinCenter(ib),ic);
     cout<<" centrality="<<centPercent[ic]<<"   ncoll <="<< histNch->GetBinCenter(ib) <<endl;
     //cout<<" centrality="<<centPercent[ic]<<" impact parameter <="<< histNch->GetBinCenter(ib) <<endl;
     ic++;
   }
   if(ic==nbins) break;
 }
 
 printf(" \n float binUp[%i] = {",nbins);  
 // cout <<" \n float multCent[nbins] = {";

 for (int ic=nbins-1; ic>-1; ic--){
   cout<< binUp->At(ic);
   if (ic!=0) cout<<", ";
 }
 cout<<"};\n"<<endl;


 printf(" \n float multCent[%i] = {",nbins);  
 // cout <<" \n float multCent[nbins] = {";

 for (int ic=nbins-1; ic>-1; ic--){
   cout<< histNch->GetBinCenter(binUp->At(ic));
   if (ic!=0) cout<<", ";
 }
 cout<<"};\n"<<endl;
}
Commit	Line	Data
ae7e261d	1	///==========================================================================
	2	///
	3	/// macro to plot centrality bin values
	4	///==========================================================================
	5	///
	6	#include <cstdlib>
	7
	8	const Int_t nbins;
	9
	10	double centPercent[]={5.,10.,20.,40.,60.,80.,99.9999999};
	11	//double centPercent[]={10.,20.,40.,60.,80.,99.9999999};
	12
	13
	14	nbins = sizeof(centPercent)/sizeof(double);
	15
	16	TArrayI* binUp = new TArrayI(nbins);
	17	TArrayF* Multbin = new TArrayF(nbins);
	18
	19	void plotGlauberCenVarsZPA()
	20	{
	21	TGraphErrors *gNpart=new TGraphErrors(0);
	22	gNpart->SetName("gNpart");
	23	TGraphErrors *gNcoll=new TGraphErrors(0);
	24	gNcoll->SetName("gNcoll");
	25	TGraphErrors *gtAA=new TGraphErrors(0);
	26	gtAA->SetName("gtAA");
	27
	28
	29	TFile *f=new TFile("ZPA_ntuple_188359.root");
	30
	31	TNtuple nt=(TNtuple)f->Get("gnt");
	32
	33	Float_t B;
	34	Float_t Npart;
	35	Float_t Ncoll;
	36	Float_t tAA;
	37	Float_t ntot;
	38
	39	nt->SetBranchAddress("B",&B);
	40	nt->SetBranchAddress("Npart",&Npart);
	41	nt->SetBranchAddress("Ncoll",&Ncoll);
	42	nt->SetBranchAddress("tAA",&tAA);
	43	nt->SetBranchAddress("ntot",&ntot);
	44
	45
	46	Int_t totev=nt->GetEntries();
	47
	48
	49	TH1F*hB=new TH1F("hB","hB",28000,0.,14.);
	50	TH1F*hNpart=new TH1F("hNpart","hNpart",420,-0.5,419.5);
	51	TH1F*hNcoll=new TH1F("hNcoll","hNcoll",3500,0,35);
	52	TH1F*htAA=new TH1F("htAA","htAA",100000,-0.5,39.5);
	53	TH1F*hntot=new TH1F("hntot","hntot",100000,-0.5,39.5);
	54
	55
	56	for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
	57	nt->GetEvent(iEvent);
	58
	59	Float_t nu = (Float_t) ntot;
	60	nu = gRandom->Gaus(nu,0.5);
	61
	62	if(nu>0) {
	63
	64	hNpart->Fill(Npart);
65	hNcoll->Fill(Ncoll);
66	hB->Fill(B);
67	htAA->Fill(tAA);
68	hntot->Fill(nu);
69	}
70	}
71
72	new TCanvas();
73	hNcoll->Draw();
74
75	//---------------------------------------------------
76	getCentrality(hntot);
77	//---------------------------------------------------
78
79
80	TH1F* hnpartcutb[nbins];
81	char histtitp[100];
82	for (Int_t i=0; i<binUp->GetSize(); i++) {
83	sprintf(histtitp,"npartcutb%d",i);
84	hnpartcutb[i] = new TH1F(histtitp,histtitp,10000,-0.5,9999.5);
85	hnpartcutb[i]->SetLineWidth(1);
86	hnpartcutb[i]->SetStats(1);
87
88	}
89
90	TH1F* hncollcutb[nbins];
91	char histtitc[100];
92	for (Int_t i=0; i<binUp->GetSize(); i++) {
93	sprintf(histtitc,"ncollcutb%d",i);
94	hncollcutb[i] = new TH1F(histtitc,histtitc,10000,-0.5,9999.5);
95	hncollcutb[i]->SetLineWidth(1);
96	hncollcutb[i]->SetStats(1);
97
98	}
99
100	TH1F* htaacutb[nbins];
101	char histtitt[100];
102	for (Int_t i=0; i<binUp->GetSize(); i++) {
103	sprintf(histtitt,"taacutb%d",i);
104	htaacutb[i] = new TH1F(histtitt,histtitt,100000,-0.5,39.5);
105	htaacutb[i]->SetLineWidth(1);
106	htaacutb[i]->SetStats(1);
107
108	}
109
110	TH1F* hbcutb[nbins];
111	char histtitt[100];
112	for (Int_t i=0; i<binUp->GetSize(); i++) {
113	sprintf(histtitt,"bcutb%d",i);
114	hbcutb[i] = new TH1F(histtitt,histtitt,100000,-0.5,39.5);
115	hbcutb[i]->SetLineWidth(1);
116	hbcutb[i]->SetStats(1);
117
118	}
119
120
121
122	for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
123	nt->GetEvent(iEvent);
124	for (int ibin=0; ibin<nbins; ibin++) {
125	// if (B<Multbin->At(ibin)) {
126
127	Float_t nu = (Float_t) Ncoll;
128	nu = nu+1.*gRandom->Rndm();
129
130	if (nu>Multbin->At(ibin)) {
131	hnpartcutb[ibin]->Fill(Npart);
132	hncollcutb[ibin]->Fill(Ncoll);
133	htaacutb[ibin]->Fill(tAA);
134	hbcutb[ibin]->Fill(B);
135	break;
136	}
137	}
138	}
139
140
141
142
143
144	// superimpose cut histograms
145	// new TCanvas();
146	// hnpartcutb[0]->Draw("");
147	for (Int_t icentr=0; icentr<nbins;icentr++) {
148	hnpartcutb[icentr]->SetLineColor(icentr);
149	hnpartcutb[icentr]->Draw("same");
150	}
151
152	for (Int_t icentr=0; icentr<nbins;icentr++) {
153	new TCanvas();
154	hnpartcutb[icentr]->SetLineColor(icentr);
155	hnpartcutb[icentr]->SetStats(1);
156	hnpartcutb[icentr]->Draw("");
157	gNpart->SetPoint(icentr,Float_t(icentr),hnpartcutb[icentr]->GetMean());
158	gNpart->SetPointError(icentr,0,hnpartcutb[icentr]->GetRMS());
159	}
160	cout << endl;
161
162	for (Int_t icentr=0; icentr<nbins;icentr++) {
163	new TCanvas();
164	hncollcutb[icentr]->SetLineColor(icentr);
165	hncollcutb[icentr]->SetStats(1);
166	hncollcutb[icentr]->Draw("");
167	gNcoll->SetPoint(icentr,Float_t(icentr),hncollcutb[icentr]->GetMean());
168	gNcoll->SetPointError(icentr,0,hncollcutb[icentr]->GetRMS());
169	}
170	cout << endl;
171
172	for (Int_t icentr=0; icentr<nbins;icentr++) {
173	new TCanvas();
174	htaacutb[icentr]->SetLineColor(icentr);
175	htaacutb[icentr]->SetStats(1);
176	htaacutb[icentr]->Draw("");
177	gtAA->SetPoint(icentr,Float_t(icentr),htaacutb[icentr]->GetMean());
178	gtAA->SetPointError(icentr,0,htaacutb[icentr]->GetRMS());
179	}
180
181	for (Int_t icentr=0; icentr<nbins;icentr++)
182	cout<<icentr<<" \| "<<setprecision(3)<<centPercent[icentr]<<" \| "
183	<<Multbin->At(icentr-1)<<" \| "<<Multbin->At(icentr)<<" \| "
184	<<hnpartcutb[icentr]->GetMean()<<" \| " <<hnpartcutb[icentr]->GetRMS()<< " \| "
185	<<hncollcutb[icentr]->GetMean()<<" \| " <<hncollcutb[icentr]->GetRMS()<< " \| "
186	<<hbcutb[icentr]->GetMean() <<" \| " <<hbcutb[icentr]->GetRMS() << " \| "
187	<<htaacutb[icentr]->GetMean() <<" \| " <<htaacutb[icentr]->GetRMS() << " \| "<<endl;
188
189
190	//TString suffixhisto=Form("/home/atoia/GlauberNtuple/62/GlauberMC_PbPb_histoB_sigma%d_mind%d_r%d_a%d_%d.root",mysignn,mymind,myr,mya,nbins);
191	//TString suffixhisto=Form("/home/atoia/GlauberNtuple/GlauberMC_PbPb_histoB_sigma%d_mind%d_r%d_a%d_%d.root",mysignn,mymind,myr,mya,nbins);
192	//const Char_t* filehistoname=suffixhisto.Data();
193	//TFile*filefinal=new TFile(filehistoname,"recreate");
194	// --FOR TEST --
195	//TFile*filefinal=new TFile("GlauberInfo.root","recreate");
196
197	// gNpart->Write();
198	// gNcoll->Write();
199	// gtAA->Write();
200
201	// hNpart->Write();
202	// hNcoll->Write();
203	// hB ->Write();
204	// htAA ->Write();
205
206	// for (int i=0; i<nbins; i++) {
207	// hnpartcutb[i]->Write();
208	// hncollcutb[i]->Write();
209	// htaacutb[i] ->Write();
210	// }
211
212	}
213
214
215	void getCentrality(TH1 *histNch, Float_t ff=1.0)
216	// histNch - histo of multiplicity distribution (better with binsize=1)x
217	// ff fraction of accepted events. All losses are assumed to occur in most
218	// peripheral bin
219	{
220
221	//double sum= histNch->GetEntries() - histNch->GetBinContent(1);
222	double sum= histNch->Integral();
223	int nbinsx=histNch->GetNbinsX();
224	double frac=0.;
225	int ic=0;
226	//for (int ib=1;ib<=nbinsx;ib++){
227	for (int ib=nbinsx;ib>0;ib--){
228	frac += histNch->GetBinContent(ib)/sum100.ff;
229	if(frac > centPercent[ic]){
230	binUp->SetAt(ib,ic);
231	Multbin->SetAt(histNch->GetBinCenter(ib),ic);
232	cout<<" centrality="<<centPercent[ic]<<" ncoll <="<< histNch->GetBinCenter(ib) <<endl;
233	//cout<<" centrality="<<centPercent[ic]<<" impact parameter <="<< histNch->GetBinCenter(ib) <<endl;
234	ic++;
235	}
236	if(ic==nbins) break;
237	}
238
239	printf(" \n float binUp[%i] = {",nbins);
240	// cout <<" \n float multCent[nbins] = {";
241
242	for (int ic=nbins-1; ic>-1; ic--){
243	cout<< binUp->At(ic);
244	if (ic!=0) cout<<", ";
245	}
246	cout<<"};\n"<<endl;
247
248
249	printf(" \n float multCent[%i] = {",nbins);
250	// cout <<" \n float multCent[nbins] = {";
251
252	for (int ic=nbins-1; ic>-1; ic--){
253	cout<< histNch->GetBinCenter(binUp->At(ic));
254	if (ic!=0) cout<<", ";
255	}
256	cout<<"};\n"<<endl;
257	}