]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - FASTSIM/uncorrBg.C
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Simplifying calling test, explicite functions for each type of test
[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C
... / ...
CommitLineData
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
18void 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}