b and c contribution scaled by 2.
[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C
1 void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
2 {
3 //
4 //
5 //
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 //    
36     Float_t scaleC0 = 2. * 1.35 * ptUp / dpt;
37     Float_t scaleB0 = 2. * 1.35 * ptUp / dpt;
38     Float_t scaleD0 = 1.35 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6
39     
40 //
41 //
42 //  Fast response
43 //    
44     AliFastMuonTriggerEff  *trigeff = new AliFastMuonTriggerEff();
45     AliFastMuonTrackingAcc *acc     = new AliFastMuonTrackingAcc();
46     AliFastMuonTrackingEff *eff     = new AliFastMuonTrackingEff();
47     AliFastMuonTrackingRes *res     = new AliFastMuonTrackingRes();
48     acc->SetBackground(1.);
49     eff->SetBackground(1.);
50     res->SetBackground(1.);  
51
52     acc    ->Init(); 
53     eff    ->Init(); 
54     res    ->Init(); 
55     trigeff->Init();
56     
57 /*
58 //  To be improved
59 //
60     AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
61     tracker->AddResponse(acc);
62     tracker->AddResponse(eff);
63     tracker->AddResponse(res);
64
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 */
73             
74 //
75 //  Heavy Flavors
76 //
77     
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
83     TF1*  ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., ptUp);
84     ptBBHf->SetParameter(0, 2.59961e-03);
85     ptBBHf->SetParameter(1, 2.41);
86     ptBBHf->SetParameter(2, 3.075);
87
88     TF1*  ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp);
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);
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);
99     
100     TF1*  ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp);
101     ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125);
102     ptBf->SetParameter(1, 6.27);
103     ptBf->SetParameter(2, 3.2);
104 //
105 //  pi/K -> mu
106 //
107     f->Close();
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 //
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.);    
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.);    
125 //
126 // Event Loop
127 //
128     Int_t iev;
129     for (iev = 0; iev < nev; iev++) {
130 //
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 //
140 // pT
141         Float_t pT1 = ptUp * gRandom->Rndm();
142         Float_t pT2 = ptUp * gRandom->Rndm();
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
150         Float_t eta1 = etar * gRandom->Rndm() + etamin;
151         Float_t eta2 = etar * gRandom->Rndm() + etamin; 
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 //
159 // Smearing (to be improved)
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;
208         
209 //
210 //      Efficiencies
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();
216         Float_t p1     = pT1/TMath::Sin(theta1 * TMath::Pi()/180.);
217         Float_t p2     = pT2/TMath::Sin(theta2 * TMath::Pi()/180.);
218         
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;
233
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         
241         if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) {
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             }
264         } // bins
265     } // event loop
266
267     Float_t evtWgt = events / Float_t(nev);
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 }