]>
Commit | Line | Data |
---|---|---|
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 | ||
19 | void 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++) { | |
a4103cd4 | 125 | e_dNdpT[ipoint]+=hInvMassSub[ipoint]->GetBinError(ibin)*hInvMassSub[ipoint]->GetBinError(ibin); |
642ca698 | 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++) { | |
a4103cd4 | 233 | e_NumberOfResonances+=hInvMassSub[NumberOfPoints]->GetBinError(ibin)*hInvMassSub[NumberOfPoints]->GetBinError(ibin); |
642ca698 | 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 | } |