]>
Commit | Line | Data |
---|---|---|
1 | void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.) | |
2 | { | |
3 | // | |
4 | // | |
5 | // | |
6 | AliFastGlauber* glauber = new AliFastGlauber(); | |
7 | glauber->Init(1); | |
8 | ||
9 | Float_t lumi = 5.e26; // cm^-2 s^-1 | |
10 | Float_t time = 1.e6; // s | |
11 | Float_t rate = lumi/1.e24 * glauber->CrossSection(bmin, bmax); // Hz | |
12 | Float_t fhard = glauber->FractionOfHardCrossSection(bmin, bmax); | |
13 | Float_t fgeo = glauber->CrossSection(bmin, bmax) / glauber->CrossSection(0, 100.); | |
14 | Float_t events = rate * time; | |
15 | ||
16 | printf("-------------------------------------------------------------------\n"); | |
17 | printf("Impact parameter range: %10.3f - %10.3f fm \n", bmin, bmax); | |
18 | printf("Luminosity: %10.3e cm^-2 s^-1 \n", lumi); | |
19 | printf("Rate: %10.3f Hz\n", rate); | |
20 | printf("Fraction of hard cross-section: %10.3f %\n", fhard * 100.); | |
21 | printf("Fraction of geom. cross-section: %10.3f %\n", fgeo * 100.); | |
22 | printf("Events in 10^6 s: %10.3e %\n", events); | |
23 | printf("-------------------------------------------------------------------\n"); | |
24 | ||
25 | // | |
26 | // | |
27 | Float_t ptMinCut = 1.; | |
28 | Float_t etamin = 2.543; | |
29 | Float_t etar = 1.457; | |
30 | Float_t ptUp = 20.; // GeV | |
31 | Float_t dpt = 0.1; // GeV | |
32 | // | |
33 | // For b = 0 | |
34 | // (factor 1.35 to scale from 10% most central to b=0) | |
35 | // | |
36 | Float_t scaleC0 = 1.35 * ptUp / dpt; | |
37 | Float_t scaleB0 = 1.35 * ptUp / dpt; | |
38 | Float_t scaleD0 = 1.35 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6 | |
39 | ||
40 | // | |
41 | // | |
42 | // Fast response | |
43 | // | |
44 | AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff(); | |
45 | AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc(); | |
46 | AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff(); | |
47 | AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes(); | |
48 | acc->SetBackground(1.); | |
49 | eff->SetBackground(1.); | |
50 | res->SetBackground(1.); | |
51 | ||
52 | acc ->Init(); | |
53 | eff ->Init(); | |
54 | res ->Init(); | |
55 | trigeff->Init(); | |
56 | ||
57 | /* | |
58 | // To be improved | |
59 | // | |
60 | AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker"); | |
61 | tracker->AddResponse(acc); | |
62 | tracker->AddResponse(eff); | |
63 | tracker->AddResponse(res); | |
64 | ||
65 | AliFastDetector* trigger = new AliFastDetector("Trigger", "Muon Trigger"); | |
66 | trigger->AddResponse(trigeff); | |
67 | ||
68 | AliFastDetector* spectro = new AliFastDetector("Spectro", "Muon Spectrometer"); | |
69 | spectro->AddSubdetector(tracker, ""); | |
70 | spectro->AddSubdetector(trigger, ""); | |
71 | spectro->Init(); | |
72 | */ | |
73 | ||
74 | // | |
75 | // Heavy Flavors | |
76 | // | |
77 | ||
78 | TF1* ptBBLf = new TF1("ptBBLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 2.); | |
79 | ptBBLf->SetParameter(0, 1.4651e-02); | |
80 | ptBBLf->SetParameter(1, 9.3409e-01); | |
81 | ptBBLf->SetParameter(2, 1.3583e+00); | |
82 | ||
83 | TF1* ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 2., ptUp); | |
84 | ptBBHf->SetParameter(0, 7.7122e-03); | |
85 | ptBBHf->SetParameter(1, 2.38); | |
86 | ptBBHf->SetParameter(2, 3.32); | |
87 | ||
88 | TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp); | |
89 | ptCCHf->SetParameter(0, 8.6675e-01); | |
90 | ptCCHf->SetParameter(1, 8.1384e-01); | |
91 | ptCCHf->SetParameter(2, 2.8933e+00); | |
92 | ptCCHf->SetParameter(3, 1.4865e-02); | |
93 | ||
94 | TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5); | |
95 | ptCCLf->SetParameter(0, 2.4899e+00); | |
96 | ptCCLf->SetParameter(1, 3.8394e-01); | |
97 | ptCCLf->SetParameter(2, 1.5505e+00); | |
98 | ptCCLf->SetParameter(3, 2.4679e-01); | |
99 | ||
100 | TF1* ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp); | |
101 | ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125); | |
102 | ptBf->SetParameter(1, 6.27); | |
103 | ptBf->SetParameter(2, 3.2); | |
104 | // | |
105 | // pi/K -> mu | |
106 | // | |
107 | f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root"); | |
108 | TH2F* etaptPiK = (TH2F*) f->Get("etaptH"); | |
109 | TAxis* etaAxis = etaptPiK->GetXaxis(); | |
110 | TAxis* ptAxis = etaptPiK->GetYaxis(); | |
111 | // | |
112 | // | |
113 | // Book histograms | |
114 | TH1F* massBBH = new TH1F("massBBH", "Mass Spectrum: b-b ", 150, 0., 15.); | |
115 | TH1F* massCCH = new TH1F("massCCH", "Mass Spectrum: c-c ", 150, 0., 15.); | |
116 | TH1F* massBCH = new TH1F("massBCH", "Mass Spectrum: b-c ", 150, 0., 15.); | |
117 | TH1F* massDDH = new TH1F("massDDH", "Mass Spectrum: decay-decay", 150, 0., 15.); | |
118 | TH1F* massBDH = new TH1F("massBDH", "Mass Spectrum: decay-b ", 150, 0., 15.); | |
119 | TH1F* massCDH = new TH1F("massCDH", "Mass Spectrum: decay-c ", 150, 0., 15.); | |
120 | TH1F* ptCH = new TH1F("ptCH", "pt Spectrum mu from c", 20, 0., 10.); | |
121 | TH1F* ptBH = new TH1F("ptBH", "pt Spectrum mu from b", 20, 0., 10.); | |
122 | TH1F* ptDH = new TH1F("ptDH", "pt Spectrum mu from pi/K", 20, 0., 10.); | |
123 | TH1F* ptBH2 = new TH1F("ptBH2", "pt Spectrum mu from b", 20, 0., 10.); | |
124 | // | |
125 | // Event Loop | |
126 | // | |
127 | Int_t iev; | |
128 | for (iev = 0; iev < nev; iev++) { | |
129 | // | |
130 | // Collision geometry | |
131 | // | |
132 | Float_t b; | |
133 | b = glauber->GetRandomImpactParameter(bmin, bmax); | |
134 | Double_t nbinary = glauber->Binaries(b); | |
135 | Float_t scaleC = scaleC0 * nbinary; | |
136 | Float_t scaleB = scaleB0 * nbinary; | |
137 | Float_t scaleD = scaleD0 * nbinary; | |
138 | // | |
139 | // pT | |
140 | Float_t pT1 = ptUp * gRandom->Rndm(); | |
141 | Float_t pT2 = ptUp * gRandom->Rndm(); | |
142 | // | |
143 | // phi | |
144 | Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi(); | |
145 | Float_t phi2 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi(); | |
146 | Float_t dphi = phi1 - phi2; | |
147 | // | |
148 | // eta | |
149 | Float_t eta1 = etar * gRandom->Rndm() + etamin; | |
150 | Float_t eta2 = etar * gRandom->Rndm() + etamin; | |
151 | Float_t deta = eta1 - eta2; | |
152 | // | |
153 | // invariant mass | |
154 | Float_t m2 = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi)); | |
155 | Float_t m = TMath::Sqrt(m2); | |
156 | ||
157 | // | |
158 | // Smearing (to be improved) | |
159 | Float_t dm = m * 0.01; | |
160 | m += gRandom->Gaus(0., dm); | |
161 | // | |
162 | // Weights | |
163 | // | |
164 | // Heavy Flavour | |
165 | // | |
166 | Int_t ibin; | |
167 | Float_t wgtB1, wgtB2; | |
168 | Float_t wgtC1, wgtC2; | |
169 | ||
170 | if (pT1 > 2.0) { | |
171 | wgtC1 = ptCCHf->Eval(pT1) * scaleC; | |
172 | } else { | |
173 | wgtC1 = ptCCLf->Eval(pT1) * scaleC; | |
174 | } | |
175 | if (pT2 > 2.0) { | |
176 | wgtC2 = ptCCHf->Eval(pT2) * scaleC; | |
177 | } else { | |
178 | wgtC2 = ptCCLf->Eval(pT2) * scaleC; | |
179 | } | |
180 | ||
181 | ||
182 | if (pT1 > 2.) { | |
183 | wgtB1 = ptBBHf->Eval(pT1) * scaleB; | |
184 | } else { | |
185 | wgtB1 = ptBBLf->Eval(pT1) * scaleB; | |
186 | } | |
187 | if (pT2 > 2.) { | |
188 | wgtB2 = ptBBHf->Eval(pT2) * scaleB; | |
189 | } else { | |
190 | wgtB2 = ptBBLf->Eval(pT2) * scaleB; | |
191 | } | |
192 | ||
193 | ||
194 | // | |
195 | // Weight for decays | |
196 | // | |
197 | Int_t etaBin, ptBin; | |
198 | Float_t wgtD1, wgtD2; | |
199 | ||
200 | etaBin = etaAxis->FindBin(eta1); | |
201 | ptBin = ptAxis ->FindBin(pT1); | |
202 | wgtD1 = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD; | |
203 | ||
204 | etaBin = etaAxis->FindBin(eta2); | |
205 | ptBin = ptAxis ->FindBin(pT2); | |
206 | wgtD2 = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD; | |
207 | ||
208 | // | |
209 | // Efficiencies | |
210 | // | |
211 | Float_t theta1 = 2. * TMath::ATan(TMath::Exp(-eta1)) * 180./TMath::Pi(); | |
212 | Float_t theta2 = 2. * TMath::ATan(TMath::Exp(-eta2)) * 180./TMath::Pi(); | |
213 | Float_t phid1 = phi1 * 180./TMath::Pi(); | |
214 | Float_t phid2 = phi2 * 180./TMath::Pi(); | |
215 | Float_t p1 = pT1/TMath::Sin(theta1 * TMath::Pi()/180.); | |
216 | Float_t p2 = pT2/TMath::Sin(theta2 * TMath::Pi()/180.); | |
217 | ||
218 | res->SetCharge(1); | |
219 | eff->SetCharge(1); | |
220 | acc->SetCharge(1); | |
221 | Float_t eff1 = eff->Evaluate(pT1, theta1, phid1); | |
222 | Float_t acc1 = acc->Evaluate(pT1, theta1, phid1); | |
223 | Float_t tri1 = trigeff->Evaluate(1, pT1, theta1, phid1); | |
224 | res->SetCharge(-1); | |
225 | eff->SetCharge(-1); | |
226 | acc->SetCharge(-1); | |
227 | Float_t eff2 = eff->Evaluate(pT2, theta2, phid2); | |
228 | Float_t acc2 = acc->Evaluate(pT2, theta2, phid2); | |
229 | Float_t tri2 = trigeff->Evaluate(1, pT2, theta2, phid2); | |
230 | ||
231 | Float_t effA = eff1 * eff2 * acc1 * acc2 * tri1 * tri2; | |
232 | ||
233 | Float_t ptMax = pT1; | |
234 | Float_t ptMin = pT2; | |
235 | if (pT2 > pT1) { | |
236 | Float_t ptMax = pT2; | |
237 | Float_t ptMin = pT1; | |
238 | } | |
239 | ||
240 | if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) { | |
241 | massBBH->Fill(m, wgtB1 * wgtB2 / 4. * effA); | |
242 | massCCH->Fill(m, wgtC1 * wgtC2 / 4. * effA); | |
243 | massBCH->Fill(m, wgtC1 * wgtB2 / 4. * effA); | |
244 | massBCH->Fill(m, wgtC2 * wgtB1 / 4. * effA); | |
245 | massDDH->Fill(m, wgtD1 * wgtD2 / 4. * effA); | |
246 | massBDH->Fill(m, wgtB1 * wgtD2 / 4. * effA); | |
247 | massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA); | |
248 | massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA); | |
249 | massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA); | |
250 | } | |
251 | // | |
252 | // pT - Spectra | |
253 | // | |
254 | for (Int_t ipt = 0; ipt < 20; ipt++) | |
255 | { | |
256 | Float_t pt = 0.5 * ipt; | |
257 | ptBH2->Fill(pT1, wgtB1); | |
258 | if (pT1 > pt) { | |
259 | ptCH->Fill(pt, wgtC1); | |
260 | ptBH->Fill(pt, wgtB1); | |
261 | ptDH->Fill(pt, wgtD1); | |
262 | } | |
263 | } // bins | |
264 | } // event loop | |
265 | ||
266 | Float_t evtWgt = events / Float_t(nev); | |
267 | ||
268 | massBBH->Scale(evtWgt); | |
269 | massCCH->Scale(evtWgt); | |
270 | massBCH->Scale(evtWgt); | |
271 | massDDH->Scale(evtWgt); | |
272 | massBDH->Scale(evtWgt); | |
273 | massCDH->Scale(evtWgt); | |
274 | ||
275 | TH1F * massALH = new TH1F(*massCDH); | |
276 | massALH->Add(massBDH); | |
277 | massALH->Add(massDDH); | |
278 | massALH->Add(massBCH); | |
279 | massALH->Add(massCCH); | |
280 | massALH->Add(massBBH); | |
281 | ||
282 | TCanvas *c0 = new TCanvas("c0","Canvas 1",400,10,600,700); | |
283 | massCCH->SetLineColor(4); | |
284 | massCCH->SetMinimum(1.); | |
285 | massCCH->SetMaximum(1.e4); | |
286 | massCCH->SetXTitle("m_{#mu#mu} [GeV]"); | |
287 | massCCH->SetYTitle("Entries/100 MeV /10^{6} s"); | |
288 | massCCH->Draw(""); | |
289 | massALH->SetLineColor(3); | |
290 | massALH->Draw("same"); | |
291 | massBBH->SetLineColor(6); | |
292 | massBBH->Draw("same"); | |
293 | ||
294 | TCanvas *c2 = new TCanvas("c2","Canvas 3",400,10,600,700); | |
295 | massDDH->SetLineColor(2); | |
296 | massDDH->SetMinimum(1.e2); | |
297 | massDDH->SetMaximum(1.e6); | |
298 | massDDH->SetXTitle("m_{#mu#mu} [GeV]"); | |
299 | massDDH->SetYTitle("Entries/100 MeV /10^{6} s"); | |
300 | massDDH->Draw(""); | |
301 | massALH->SetLineColor(3); | |
302 | massALH->Draw("same"); | |
303 | massCCH->SetLineColor(4); | |
304 | massCCH->Draw("same"); | |
305 | massBBH->SetLineColor(6); | |
306 | massBBH->Draw("same"); | |
307 | ||
308 | TCanvas *c3 = new TCanvas("c3","Canvas 4",400,10,600,700); | |
309 | ptCH->Scale(1./float(nev)); | |
310 | ptBH->Scale(1./float(nev)); | |
311 | ptDH->Scale(1./float(nev)); | |
312 | ptCH->SetLineColor(4); | |
313 | ptBH->SetLineColor(6); | |
314 | ptDH->SetLineColor(2); | |
315 | ptCH->SetXTitle("p_{T}^{min} [GeV]"); | |
316 | ptCH->SetYTitle("<n>_{#mu}/event"); | |
317 | ||
318 | ptDH->Draw(); | |
319 | ptBH->Draw("same"); | |
320 | ptCH->Draw("same"); | |
321 | TCanvas *c4 = new TCanvas("c4","Canvas 5",400,10,600,700); | |
322 | ptBH2->Draw(); | |
323 | ||
324 | } |