]> 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
removing unnecessary include file
[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C
CommitLineData
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 18void 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}