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