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