Adding renormalisation of the combinatotial bgk for peripheral events
[u/mrichter/AliRoot.git] / MUON / MUONQuarkoniaPtDistri.C
CommitLineData
4692d494 1//
2// Macro for analysis of 2d InvMass v Pt distribution for the generation of:
3// - Pt distribution of quarkonia,
4// - Mass and sigma of the quarkonia resonance as a function of pT
5// - and invariant mass distributions ( all, background and subtracted)
6// as a function of pT
7// Usual utilisation:
8// QuarkoniaPtDistri("ProdMuonCocktailPer1/MUONmassPlot.root",4,2.9,3.3)
9//
10// MUONmassPlotFile root file form MUONmassPlot_ESD.C
11// NRes=1 bining of the pT distribution
12// ==1, equal to MUONmassPlotFile 2D histo,
13// ==2, bin twice larger
14// InvMass1, InvMass2 Range of the quarkonia resonance
15//
16// Gines MARTINEZ, Subatech, May 04
17//
18
19void MUONQuarkoniaPtDistri(char * MUONmassPlotFile="MUONmassPlot.root", Int_t NRes=4, Float_t InvMass1=2.9, Float_t InvMass2=3.3)
20{
21 TFile MUONmassPlot(MUONmassPlotFile);
22
23 TH2F * hInvMassAll_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassAll_vs_Pt");
24 TH2F * hInvMassBgk_vs_Pt = (TH2F*) MUONmassPlot.Get("hInvMassBgk_vs_Pt");
25 hInvMassAll_vs_Pt->Sumw2();
26 hInvMassBgk_vs_Pt->Sumw2();
27
28 gROOT->cd();
29
30 Int_t Nbins = hInvMassAll_vs_Pt->GetYaxis()->GetNbins();
31 Float_t ptMin = hInvMassAll_vs_Pt->GetYaxis()->GetXmin();
32 Float_t ptMax = hInvMassAll_vs_Pt->GetYaxis()->GetXmax();
33
34 const Int_t NumberOfPoints = Nbins/NRes;
35 Float_t *pT = new Float_t[NumberOfPoints];
36 Float_t *e_pT = new Float_t[NumberOfPoints];
37 Float_t *dNdpT = new Float_t[NumberOfPoints];
38 Float_t *e_dNdpT = new Float_t[NumberOfPoints];
39
40 Float_t *pT1 = new Float_t[NumberOfPoints];
41 Float_t *e_pT1 = new Float_t[NumberOfPoints];
42 Float_t *Cent = new Float_t[NumberOfPoints];
43 Float_t *e_Cent = new Float_t[NumberOfPoints];
44
45 Float_t *pT2 = new Float_t[NumberOfPoints];
46 Float_t *e_pT2 = new Float_t[NumberOfPoints];
47 Float_t *Sigma = new Float_t[NumberOfPoints];
48 Float_t *e_Sigma = new Float_t[NumberOfPoints];
49
642ca698 50 TH1D * hInvMassAll[NumberOfPoints+1];
51 TH1D * hInvMassBgk[NumberOfPoints+1];
52 TH1D * hInvMassSub[NumberOfPoints+1];
4692d494 53
54 char histoname[60];
55
56 Int_t ibin, ibin1,ibin2, ipoint;
57 Int_t integral1, integral2;
642ca698 58 for(ipoint=0; ipoint<=NumberOfPoints; ipoint++) {
59
60 if (ipoint!=NumberOfPoints) {
61 ibin1 = ipoint * NRes;
62 ibin2 = ibin1+ NRes;
63 pT[ipoint]=0.;
64 e_pT[ipoint]=0.;
65 dNdpT[ipoint]=0.;
66 e_dNdpT[ipoint]=0.;
67 for(ibin=ibin1; ibin<ibin2; ibin++) {
68 pT[ipoint]+=hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin);
69 }
70 pT[ipoint]/=NRes;
71 pT1[ipoint]=pT[ipoint];
72 pT2[ipoint]=pT[ipoint];
73
74 e_pT[ipoint] = (hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin2)-hInvMassAll_vs_Pt->GetYaxis()->GetBinCenter(ibin1))/2.;
75 e_pT1[ipoint]=e_pT[ipoint];
76 e_pT2[ipoint]=e_pT[ipoint];
77 }
78 else {
79 ibin1 = 0;
80 ibin2 = Nbins;
4692d494 81 }
4692d494 82
4692d494 83
642ca698 84 if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassAll_%2.0f",100*pT[ipoint]);
85 else sprintf(histoname,"hInvMassAll_full");
4692d494 86 TH1D * nume = hInvMassAll_vs_Pt->ProjectionX("nume",ibin1,ibin2-1);
87 hInvMassAll[ipoint] = new TH1D(*nume);
88 hInvMassAll[ipoint]->SetName(histoname); hInvMassAll[ipoint]->SetTitle(histoname);
89 hInvMassAll[ipoint]->Sumw2();
90
642ca698 91 if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassSub_%2.0f",100*pT[ipoint]);
92 else sprintf(histoname,"hInvMassSub_full");
4692d494 93 hInvMassSub[ipoint] = new TH1D(*nume);
94 hInvMassSub[ipoint]->SetName(histoname); hInvMassSub[ipoint]->SetTitle(histoname);
95 hInvMassSub[ipoint]->Sumw2();
96
642ca698 97 if (ipoint!=NumberOfPoints) sprintf(histoname,"hInvMassBgk_%2.0f",100*pT[ipoint]);
98 else sprintf(histoname,"hInvMassBgk_full");
4692d494 99 TH1D * deno = hInvMassBgk_vs_Pt->ProjectionX("deno",ibin1,ibin2-1);
100 hInvMassBgk[ipoint]= new TH1D(*deno);
101 hInvMassBgk[ipoint]->SetName(histoname);hInvMassBgk[ipoint]->SetTitle(histoname);
102 hInvMassBgk[ipoint]->Sumw2();
103
104 hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
105 hInvMassSub[ipoint]->Fit("gaus","","",InvMass1,InvMass2);
642ca698 106
4692d494 107 integral1 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
108 integral2 = hInvMassSub[ipoint]->FindBin( (hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2) ) );
642ca698 109
110 // Renormalisation to take into account the resonance statistics
111 if ( hInvMassBgk[ipoint]->Integral()>0.) {
112 Float_t renormalisation = (hInvMassBgk[ipoint]->Integral() - hInvMassSub[ipoint]->Integral(integral1,integral2))/hInvMassBgk[ipoint]->Integral();
113 hInvMassBgk[ipoint]->Scale(renormalisation);
114 hInvMassSub[ipoint]->Reset("ICE");
115 hInvMassSub[ipoint]->Add(hInvMassAll[ipoint],+1.);
116 hInvMassSub[ipoint]->Add(hInvMassBgk[ipoint],-1.);
117 }
118 if (ipoint!=NumberOfPoints) {
119 Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(1);
120 e_Cent[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(1);
121 Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParameter(2);
122 e_Sigma[ipoint] = hInvMassSub[ipoint]->GetFunction("gaus")->GetParError(2);
123 dNdpT[ipoint] = hInvMassSub[ipoint]->Integral(integral1,integral2);
124 for(ibin=integral1; ibin<=integral2; ibin++) {
125 e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinContent(ibin)*hInvMassSub[ipoint]->GetBinContent(ibin);
126 }
127 e_dNdpT[ipoint] = TMath::Sqrt( e_dNdpT[ipoint]);
4692d494 128 }
4692d494 129 }
130
131 MUONmassPlot.Close();
132
133 TFile out("QuarkoniaPtDistri.root","RECREATE");
134 TGraphErrors * gePtDistri = new TGraphErrors(NumberOfPoints,pT,dNdpT,e_pT, e_dNdpT);
135 TGraphErrors * geCenter = new TGraphErrors(NumberOfPoints,pT1,Cent,e_pT1, e_Cent);
136 TGraphErrors * geSigma = new TGraphErrors(NumberOfPoints,pT2,Sigma,e_pT2, e_Sigma);
137 gePtDistri->Write();
138 geCenter->Write();
139 geSigma->Write();
140 for(ipoint=0;ipoint<NumberOfPoints; ipoint++) {
141 hInvMassAll[ipoint]->Write();
142 hInvMassBgk[ipoint]->Write();
143 hInvMassSub[ipoint]->Write();
144 }
145 out.Write();
146 out.Close();
147
148 TCanvas *c1 =new TCanvas("c1","c1");
149 TH2D * hframe1 = new TH2D("hframe1","hframe1",10,0,12,10,InvMass1,InvMass2);
150 hframe1->SetXTitle("Transverse Momentum (GeV/c)");
151 hframe1->SetYTitle("Resonance Centroide (GeV/C^{2})");
152 hframe1->Draw();
153 geCenter->SetMarkerStyle(27);
154 geCenter->Draw("p");
155
156 TCanvas *c2 =new TCanvas("c2","c2");
157 TH2D * hframe2 = new TH2D("hframe2","hframe2",10,0,12,10,0.03,0.200);
158 hframe2->SetXTitle("Transverse Momentum (GeV/c)");
159 hframe2->SetYTitle("Resonance Sigma (GeV/C^2)");
160 hframe2->Draw();
161 geSigma->SetMarkerStyle(28);
162 geSigma->Draw("p");
163
164 TCanvas *c3 =new TCanvas("c3","c3");
165 TH2D * hframe3 = new TH2D("hframe3","hframe3",10,0,12,10,0.1,10000);
166 hframe3->SetXTitle("Transverse Momentum (GeV/c)");
167 hframe3->SetYTitle("dN/dpT (GeV^{-1})");
168 hframe3->Draw();
169 c3->SetLogy();
170 hframe3->Draw();
171 gePtDistri->SetMarkerStyle(29);
172 gePtDistri->Draw("p");
173
174 TCanvas *c4 =new TCanvas("c4","c4");
175 c4->SetLogy();
176 c4->Divide(3,3);
177 Int_t i;
178 Float_t izoom1;
179 Float_t izoom2;
180 for (i=1;i<10;i++) {
181 izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
182 izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
183 c4->cd(i);
184 hInvMassAll[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
185 hInvMassAll[i]->SetYTitle("Counts (1/25 MeV)");
186 hInvMassAll[i]->GetXaxis()->SetRange(izoom1,izoom2);
187 hInvMassAll[i]->Draw();
188 hInvMassBgk[i]->SetLineColor(4);
189 hInvMassBgk[i]->Draw("histo,same");
190 }
191
192 TCanvas *c5 =new TCanvas("c5","c5");
193 c5->Divide(3,3);
194 Int_t i;
195 for (i=1;i<10;i++) {
196 c5->cd(i);
197 izoom1 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
198 izoom2 = hInvMassSub[i]->FindBin( (hInvMassSub[i]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[i]->GetFunction("gaus")->GetParameter(2) ) );
199 hInvMassSub[i]->GetXaxis()->SetRange(izoom1,izoom2);
200 hInvMassSub[i]->SetXTitle("Invariant Mass (GeV/C^{2})");
201 hInvMassSub[i]->SetYTitle("Counts (1/25 MeV)");
202 hInvMassSub[i]->Draw();
203 }
204
642ca698 205 TCanvas *c6 =new TCanvas("c6","c6");
206 c6->Divide(2,1);
207 c6->cd(1);
208 izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
209 izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
210 c4->cd(i);
211 hInvMassAll[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
212 hInvMassAll[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
213 hInvMassAll[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
214 hInvMassAll[NumberOfPoints]->Draw();
215 hInvMassBgk[NumberOfPoints]->SetLineColor(4);
216 hInvMassBgk[NumberOfPoints]->Draw("histo,same");
217 c6->cd(2);
218 izoom1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
219 izoom2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+15.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
220 hInvMassSub[NumberOfPoints]->GetXaxis()->SetRange(izoom1,izoom2);
221 hInvMassSub[NumberOfPoints]->SetXTitle("Invariant Mass (GeV/C^{2})");
222 hInvMassSub[NumberOfPoints]->SetYTitle("Counts (1/25 MeV)");
223 hInvMassSub[NumberOfPoints]->Draw();
224
225 //Number of resonances
226
227 integral1 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
228 integral2 = hInvMassSub[NumberOfPoints]->FindBin( (hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ) );
229
230 Int_t NumberOfResonances = hInvMassSub[NumberOfPoints]->Integral(integral1,integral2);
231 Int_t e_NumberOfResonances = 0;
232 for(ibin=integral1; ibin<=integral2; ibin++) {
233 e_NumberOfResonances+=hInvMassSub[NumberOfPoints]->GetBinContent(ibin)*hInvMassSub[NumberOfPoints]->GetBinContent(ibin);
234 }
235 e_NumberOfResonances= TMath::Sqrt( e_NumberOfResonances);
236 printf(">>> Number of resonances is %d+-%d in the inv mass range (2sigma) %f, %f \n",NumberOfResonances,e_NumberOfResonances,
237 hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)-2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) ,
238 hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1)+2.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2) );;
239
240 printf(">>> Centroide is %f +-%f GeV\n", hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(1), hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(1));
241 printf(">>> Sigma is %f +-%f MeV \n", 1000.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParameter(2), 1000.*hInvMassSub[NumberOfPoints]->GetFunction("gaus")->GetParError(2));
242
4692d494 243}