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