]> git.uio.no Git - u/mrichter/AliRoot.git/blame - FASTSIM/uncorrBg.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Track references deleted in AliModule
[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C
CommitLineData
d0768bbb 1void uncorrBg(Int_t nev = 1000000)
2{
3// b: 0.00261 (0.0375 per unit of rapidity per event)
4// c: 1.029 (1.35 per unit of rapidity per event)
5//
6// rate of central collisions 400 Hz
7// central events 4 10^8
8// Get the single muon pT spectra
9//
10//
11 Float_t scaleC = 200.;
12 Float_t scaleB = 200.;
13 Float_t scaleD = 0.34 * 1.362e-3 * 1.e3;
14
15//
16 Float_t bbx[15] =
17 {
18 0.21, 0.55, 0.65, 0.65, 0.5,
19 0.40, 0.26, 0.19, 0.11, 0.1,
20 0.075, 0.045, 0.035, 0.02, 0.017
21 };
22//
23// Fast response
24//
25 AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
26 AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
27 AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
28 AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
29 acc->SetBackground(1.);
30 eff->SetBackground(1.);
31 res->SetBackground(1.);
32 acc->Init();
33 eff->Init();
34 res->Init();
35 AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
36 tracker->AddResponse(trigeff);
37 tracker->AddResponse(acc);
38 tracker->AddResponse(eff);
39 tracker->AddResponse(res);
40 tracker->Init();
41 tracker->Dump();
42
43
44//
45// Heavy Flavors
46//
47 f = new TFile("HVQinc.root");
48 TH1F* ptBB = (TH1F*) f->Get("hPtCorra");
49 TH1F* ptCC = (TH1F*) f->Get("hpta");
50 TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700);
51
52 TF1* ptBBLf = new TF1("ptBBLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 3.);
53 ptBBLf->SetParameter(0, 4.46695e-03);
54 ptBBLf->SetParameter(1, 1.60242e+00);
55 ptBBLf->SetParameter(2, 2.24948e+00);
56
57 TF1* ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., 20.);
58 ptBBHf->SetParameter(0, 2.59961e-03);
59 ptBBHf->SetParameter(1, 2.41);
60 ptBBHf->SetParameter(2, 3.075);
61
62 TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, 20.);
63 ptCCHf->SetParameter(0, 6.72360e-01);
64 ptCCHf->SetParameter(1, 7.06043e-01);
65 ptCCHf->SetParameter(2, 2.74240e+00);
66 ptCCHf->SetParameter(3, 8.45018e-03);
67// ptCCHf->Draw("ac");
68
69 TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5);
70 ptCCLf->SetParameter(0, 1.40260e+00);
71 ptCCLf->SetParameter(1, 3.75762e-01);
72 ptCCLf->SetParameter(2, 1.54345e+00);
73 ptCCLf->SetParameter(3, 2.49806e-01);
74// ptCCLf->Draw("ac");
75 /*
76
77 TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
78 ptCCHf->SetParameter(0, 0.249);
79 ptCCHf->SetParameter(1, 1.15);
80 ptCCHf->SetParameter(2, 3.33);
81// ptCCHf->Draw("ac");
82
83 TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
84 ptCCLf->SetParameter(0, 1.125);
85 ptCCLf->SetParameter(1, 0.525);
86 ptCCLf->SetParameter(2, 2.42);
87// ptCCLf->Draw("ac");
88*/
89
90
91 TF1* ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
92 ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125);
93 ptBf->SetParameter(1, 6.27);
94 ptBf->SetParameter(2, 3.2);
95// ptBf->Draw("ac");
96//
97// pi/K -> mu
98//
99 f->Close();
100 f = new TFile("pikmu.root");
101 TH2F* etaptPiK = (TH2F*) f->Get("etaptH");
102 TAxis* etaAxis = etaptPiK->GetXaxis();
103 TAxis* ptAxis = etaptPiK->GetYaxis();
104 TH1F* ptPiK = (TH1F*) f->Get("ptH3");
105// ptAxis = ptPiK->GetXaxis();
106
107//
108// Book histograms
109 TH1F* massBBH = new TH1F("massBBH", "Mass Spectrum: b-b ", 150, 0., 15.);
110 TH1F* massCCH = new TH1F("massCCH", "Mass Spectrum: c-c ", 150, 0., 15.);
111 TH1F* massBCH = new TH1F("massBCH", "Mass Spectrum: b-c ", 150, 0., 15.);
112 TH1F* massDDH = new TH1F("massDDH", "Mass Spectrum: decay-decay", 150, 0., 15.);
113 TH1F* massBDH = new TH1F("massBDH", "Mass Spectrum: decay-b ", 150, 0., 15.);
114 TH1F* massCDH = new TH1F("massCDH", "Mass Spectrum: decay-c ", 150, 0., 15.);
115 TH1F* ptCH = new TH1F("ptCH", "pt Spectrum mu from c", 20., 0., 10.);
116 TH1F* ptBH = new TH1F("ptBH", "pt Spectrum mu from b", 20., 0., 10.);
117 TH1F* ptDH = new TH1F("ptDH", "pt Spectrum mu from pi/K", 20., 0., 10.);
118 TH1F* ptBH2 = new TH1F("ptBH2", "pt Spectrum mu from b", 20., 0., 10.);
119 TH1F* ptBH3 = new TH1F("ptBH3", "pt Spectrum mu from b", 15., 0., 15.);
120 for (Int_t i = 0; i < 15; i++)
121 {
122 ptBH3->SetBinContent(i+1, bbx[i]);
123 }
124 ptBH3->Draw();
125
126
127//
128// Event Loop
129//
130 Int_t iev;
131 for (iev = 0; iev < nev; iev++) {
132//
133// pT
134 Float_t pT1 = 20. * gRandom->Rndm();
135 Float_t pT2 = 20. * gRandom->Rndm();
136//
137// phi
138 Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
139 Float_t phi2 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
140 Float_t dphi = phi1 - phi2;
141//
142// eta
143 Float_t eta1 = 1.457 * gRandom->Rndm() + 2.543;
144 Float_t eta2 = 1.457 * gRandom->Rndm() + 2.543;
145 Float_t deta = eta1 - eta2;
146//
147// invariant mass
148 Float_t m2 = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi));
149 Float_t m = TMath::Sqrt(m2);
150
151//
152// Smearing
153 Float_t dm = m * 0.01;
154 m += gRandom->Gaus(0., dm);
155//
156// Weights
157//
158// Heavy Flavour
159//
160 Int_t ibin;
161 Float_t wgtB1, wgtB2;
162 Float_t wgtC1, wgtC2;
163
164 if (pT1 > 1.5) {
165 wgtC1 = ptCCHf->Eval(pT1) * scaleC;
166 } else {
167 wgtC1 = ptCCLf->Eval(pT1) * scaleC;
168 }
169 if (pT2 > 1.5) {
170 wgtC2 = ptCCHf->Eval(pT2) * scaleC;
171 } else {
172 wgtC2 = ptCCLf->Eval(pT2) * scaleC;
173 }
174
175
176 if (pT1 > 3.) {
177 wgtB1 = ptBBHf->Eval(pT1) * scaleB;
178 } else {
179 wgtB1 = ptBBLf->Eval(pT1) * scaleB;
180 }
181 if (pT2 > 3.) {
182 wgtB2 = ptBBHf->Eval(pT2) * scaleB;
183 } else {
184 wgtB2 = ptBBLf->Eval(pT2) * scaleB;
185 }
186
187
188//
189// Weight for decays
190//
191 Int_t etaBin, ptBin;
192 Float_t wgtD1, wgtD2;
193
194 etaBin = etaAxis->FindBin(eta1);
195 ptBin = ptAxis ->FindBin(pT1);
196 wgtD1 = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
197
198 etaBin = etaAxis->FindBin(eta2);
199 ptBin = ptAxis ->FindBin(pT2);
200 wgtD2 = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
201
202 /*
203 ptBin = ptAxis ->FindBin(pT1);
204 wgtD1 = ptPiK->GetBinContent(ptBin) * scaleD;
205 ptBin = ptAxis ->FindBin(pT2);
206 wgtD2 = ptPiK->GetBinContent(ptBin) * scaleD;
207 */
208
209//
210// efficiencies
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();
216
217 res->SetCharge(1);
218 eff->SetCharge(1);
219 acc->SetCharge(1);
220 Float_t eff1 = eff->Evaluate(pT1, theta1, phid1);
221 Float_t acc1 = acc->Evaluate(pT1, theta1, phid1);
222 Float_t tri1 = trigeff->Evaluate(1, pT1, theta1, phid1);
223 res->SetCharge(-1);
224 eff->SetCharge(-1);
225 acc->SetCharge(-1);
226 Float_t eff2 = eff->Evaluate(pT2, theta2, phid2);
227 Float_t acc2 = acc->Evaluate(pT2, theta2, phid2);
228 Float_t tri2 = trigeff->Evaluate(1, pT2, theta2, phid2);
229
230 Float_t effA = eff1 * eff2 * acc1 * acc2 * tri1 * tri2;
231
232 Float_t ptMax = pT1;
233 Float_t ptMin = pT2;
234 if (pT2 > pT1) {
235 Float_t ptMax = pT2;
236 Float_t ptMin = pT1;
237 }
238
239
240// if (
241// (ptMax > 6. && ptMin > 3.) ||
242// (ptMax < 6. && ptMin > (6. - 0.5 * ptMax))
243// )
244// {
245 if (ptMin > 3.) {
246 massBBH->Fill(m, wgtB1 * wgtB2 / 4. * effA);
247 massCCH->Fill(m, wgtC1 * wgtC2 / 4. * effA);
248 massBCH->Fill(m, wgtC1 * wgtB2 / 4. * effA);
249 massBCH->Fill(m, wgtC2 * wgtB1 / 4. * effA);
250 massDDH->Fill(m, wgtD1 * wgtD2 / 4. * effA);
251 massBDH->Fill(m, wgtB1 * wgtD2 / 4. * effA);
252 massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA);
253 massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA);
254 massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA);
255 }
256 //
257 // pT - Spectra
258 //
259 for (Int_t ipt = 0; ipt < 20; ipt++)
260 {
261 Float_t pt = 0.5 * ipt;
262 ptBH2->Fill(pT1, wgtB1);
263 if (pT1 > pt) {
264 ptCH->Fill(pt, wgtC1);
265 ptBH->Fill(pt, wgtB1);
266 ptDH->Fill(pt, wgtD1);
267 }
268 }
269
270 } // event loop
271// Float_t eff = 0.6 * 1.0;
272 Float_t evtWgt = 1. / Float_t(nev) * 4.e8;
273
274
275 massBBH->Scale(evtWgt);
276 massCCH->Scale(evtWgt);
277 massBCH->Scale(evtWgt);
278 massDDH->Scale(evtWgt);
279 massBDH->Scale(evtWgt);
280 massCDH->Scale(evtWgt);
281
282 TH1F * massALH = new TH1F(*massCDH);
283 massALH->Add(massBDH);
284 massALH->Add(massDDH);
285 massALH->Add(massBCH);
286 massALH->Add(massCCH);
287 massALH->Add(massBBH);
288
289 TCanvas *c0 = new TCanvas("c0","Canvas 1",400,10,600,700);
290 massCCH->SetLineColor(4);
291 massCCH->SetMinimum(1.);
292 massCCH->SetMaximum(1.e4);
293 massCCH->SetXTitle("m_{#mu#mu} [GeV]");
294 massCCH->SetYTitle("Entries/100 MeV /10^{6} s");
295 massCCH->Draw("");
296 massALH->SetLineColor(3);
297 massALH->Draw("same");
298 massBBH->SetLineColor(6);
299 massBBH->Draw("same");
300
301 TCanvas *c2 = new TCanvas("c2","Canvas 3",400,10,600,700);
302 massDDH->SetLineColor(2);
303 massDDH->SetMinimum(1.e2);
304 massDDH->SetMaximum(1.e6);
305 massDDH->SetXTitle("m_{#mu#mu} [GeV]");
306 massDDH->SetYTitle("Entries/100 MeV /10^{6} s");
307 massDDH->Draw("");
308 massALH->SetLineColor(3);
309 massALH->Draw("same");
310 massCCH->SetLineColor(4);
311 massCCH->Draw("same");
312 massBBH->SetLineColor(6);
313 massBBH->Draw("same");
314
315 TCanvas *c3 = new TCanvas("c3","Canvas 4",400,10,600,700);
316 ptCH->Scale(1./float(nev));
317 ptBH->Scale(1./float(nev));
318 ptDH->Scale(1./float(nev));
319 ptCH->SetLineColor(4);
320 ptBH->SetLineColor(6);
321 ptDH->SetLineColor(2);
322 ptCH->SetXTitle("p_{T}^{min} [GeV]");
323 ptCH->SetYTitle("<n>_{#mu}/event");
324
325 ptDH->Draw();
326 ptBH->Draw("same");
327 ptCH->Draw("same");
328 TCanvas *c4 = new TCanvas("c4","Canvas 5",400,10,600,700);
329 ptBH2->Draw();
330
331}