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