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