]>
Commit | Line | Data |
---|---|---|
ec7f2c5f | 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 | ||
a51eecef | 18 | void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.) |
d0768bbb | 19 | { |
d0768bbb | 20 | // |
d0768bbb | 21 | // |
d0768bbb | 22 | // |
a51eecef | 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 | |
8f94e97c | 28 | Float_t rate = lumi/1.e27 * glauber->CrossSection(bmin, bmax); // Hz |
a51eecef | 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); | |
ec7f2c5f | 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); | |
a51eecef | 40 | printf("-------------------------------------------------------------------\n"); |
41 | ||
42 | // | |
43 | // | |
ec7f2c5f | 44 | Float_t ptMinCut = 3.; // GeV |
a51eecef | 45 | Float_t etamin = 2.543; |
46 | Float_t etar = 1.457; | |
ec7f2c5f | 47 | Float_t ptUp = 20.; // GeV |
1b882d7d | 48 | Float_t dpt = 0.01; // GeV |
a51eecef | 49 | // |
50 | // For b = 0 | |
ce839c0a | 51 | // (factor 1.28 to scale from 10% most central to b=0) |
1b882d7d | 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; | |
ce839c0a | 57 | Float_t scaleD0 = 1.28 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6 |
a51eecef | 58 | |
59 | // | |
d0768bbb | 60 | // |
61 | // Fast response | |
62 | // | |
a51eecef | 63 | AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff(); |
64 | AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc(); | |
65 | AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff(); | |
66 | AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes(); | |
d0768bbb | 67 | acc->SetBackground(1.); |
68 | eff->SetBackground(1.); | |
69 | res->SetBackground(1.); | |
a51eecef | 70 | |
71 | acc ->Init(); | |
72 | eff ->Init(); | |
73 | res ->Init(); | |
74 | trigeff->Init(); | |
75 | ||
76 | /* | |
77 | // To be improved | |
78 | // | |
d0768bbb | 79 | AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker"); |
d0768bbb | 80 | tracker->AddResponse(acc); |
81 | tracker->AddResponse(eff); | |
82 | tracker->AddResponse(res); | |
d0768bbb | 83 | |
a51eecef | 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 | */ | |
d0768bbb | 92 | |
93 | // | |
94 | // Heavy Flavors | |
95 | // | |
a51eecef | 96 | |
1b882d7d | 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); | |
d0768bbb | 101 | |
1b882d7d | 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); | |
d0768bbb | 106 | |
a51eecef | 107 | TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp); |
1b882d7d | 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); | |
d0768bbb | 112 | |
113 | TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5); | |
1b882d7d | 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); | |
d0768bbb | 118 | |
a51eecef | 119 | TF1* ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp); |
d0768bbb | 120 | ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125); |
121 | ptBf->SetParameter(1, 6.27); | |
122 | ptBf->SetParameter(2, 3.2); | |
d0768bbb | 123 | // |
124 | // pi/K -> mu | |
125 | // | |
ec7f2c5f | 126 | TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root"); |
a51eecef | 127 | TH2F* etaptPiK = (TH2F*) f->Get("etaptH"); |
128 | TAxis* etaAxis = etaptPiK->GetXaxis(); | |
129 | TAxis* ptAxis = etaptPiK->GetYaxis(); | |
130 | // | |
d0768bbb | 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.); | |
a51eecef | 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.); | |
1b882d7d | 143 | TH1F* costBBH = new TH1F("costBBH", "cos theta* ", 20, -1, 1.); |
d0768bbb | 144 | // |
145 | // Event Loop | |
146 | // | |
147 | Int_t iev; | |
148 | for (iev = 0; iev < nev; iev++) { | |
149 | // | |
a51eecef | 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 | // | |
d0768bbb | 159 | // pT |
a51eecef | 160 | Float_t pT1 = ptUp * gRandom->Rndm(); |
161 | Float_t pT2 = ptUp * gRandom->Rndm(); | |
d0768bbb | 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 | |
a51eecef | 169 | Float_t eta1 = etar * gRandom->Rndm() + etamin; |
170 | Float_t eta2 = etar * gRandom->Rndm() + etamin; | |
d0768bbb | 171 | Float_t deta = eta1 - eta2; |
172 | // | |
173 | // invariant mass | |
1b882d7d | 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; | |
d0768bbb | 179 | |
180 | // | |
a51eecef | 181 | // Smearing (to be improved) |
d0768bbb | 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 | ||
ce839c0a | 193 | if (pT1 > 1.5) { |
d0768bbb | 194 | wgtC1 = ptCCHf->Eval(pT1) * scaleC; |
195 | } else { | |
196 | wgtC1 = ptCCLf->Eval(pT1) * scaleC; | |
197 | } | |
ce839c0a | 198 | if (pT2 > 1.5) { |
d0768bbb | 199 | wgtC2 = ptCCHf->Eval(pT2) * scaleC; |
200 | } else { | |
201 | wgtC2 = ptCCLf->Eval(pT2) * scaleC; | |
202 | } | |
203 | ||
204 | ||
1b882d7d | 205 | if (pT1 > 3.) { |
d0768bbb | 206 | wgtB1 = ptBBHf->Eval(pT1) * scaleB; |
207 | } else { | |
208 | wgtB1 = ptBBLf->Eval(pT1) * scaleB; | |
209 | } | |
1b882d7d | 210 | if (pT2 > 3.) { |
d0768bbb | 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; | |
a51eecef | 230 | |
d0768bbb | 231 | // |
a51eecef | 232 | // Efficiencies |
d0768bbb | 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(); | |
a51eecef | 238 | Float_t p1 = pT1/TMath::Sin(theta1 * TMath::Pi()/180.); |
239 | Float_t p2 = pT2/TMath::Sin(theta2 * TMath::Pi()/180.); | |
240 | ||
d0768bbb | 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); | |
3e1ba7de | 252 | Float_t tri2 = trigeff->Evaluate(-1, pT2, theta2, phid2); |
d0768bbb | 253 | |
254 | Float_t effA = eff1 * eff2 * acc1 * acc2 * tri1 * tri2; | |
a51eecef | 255 | |
d0768bbb | 256 | Float_t ptMax = pT1; |
257 | Float_t ptMin = pT2; | |
258 | if (pT2 > pT1) { | |
ec7f2c5f | 259 | ptMax = pT2; |
260 | ptMin = pT1; | |
d0768bbb | 261 | } |
262 | ||
a51eecef | 263 | if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) { |
d0768bbb | 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); | |
1b882d7d | 273 | |
274 | costBBH->Fill(cost, wgtB1 * wgtB2 / 4. * effA); | |
d0768bbb | 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 | } | |
a51eecef | 288 | } // bins |
d0768bbb | 289 | } // event loop |
a51eecef | 290 | |
291 | Float_t evtWgt = events / Float_t(nev); | |
d0768bbb | 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(); | |
1b882d7d | 348 | |
349 | TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700); | |
350 | costBBH->Draw(); | |
d0768bbb | 351 | |
352 | } |