+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <stdio.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TDatabasePDG.h>
+#include <TRandom.h>
+#include <TStopwatch.h>
+#include "AliFastGlauber.h"
+#include "AliFastMuonTrackingRes.h"
+#include "AliFastMuonTrackingAcc.h"
+#include "AliFastMuonTrackingEff.h"
+#include "AliFastMuonTriggerEff.h"
+#endif
+
void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
{
//
Float_t lumi = 5.e26; // cm^-2 s^-1
Float_t time = 1.e6; // s
- Float_t rate = lumi/1.e24 * glauber->CrossSection(bmin, bmax); // Hz
+ Float_t rate = lumi/1.e27 * glauber->CrossSection(bmin, bmax); // Hz
Float_t fhard = glauber->FractionOfHardCrossSection(bmin, bmax);
Float_t fgeo = glauber->CrossSection(bmin, bmax) / glauber->CrossSection(0, 100.);
Float_t events = rate * time;
printf("Impact parameter range: %10.3f - %10.3f fm \n", bmin, bmax);
printf("Luminosity: %10.3e cm^-2 s^-1 \n", lumi);
printf("Rate: %10.3f Hz\n", rate);
- printf("Fraction of hard cross-section: %10.3f %\n", fhard * 100.);
- printf("Fraction of geom. cross-section: %10.3f %\n", fgeo * 100.);
- printf("Events in 10^6 s: %10.3e %\n", events);
+ printf("Fraction of hard cross-section: %10.3f \n", fhard * 100.);
+ printf("Fraction of geom. cross-section: %10.3f \n", fgeo * 100.);
+ printf("Events in 10^6 s: %10.3e \n", events);
printf("-------------------------------------------------------------------\n");
//
//
- Float_t ptMinCut = 1.;
+ Float_t ptMinCut = 3.; // GeV
Float_t etamin = 2.543;
Float_t etar = 1.457;
- Float_t ptUp = 20.; // GeV
- Float_t dpt = 0.1; // GeV
+ Float_t ptUp = 20.; // GeV
+ Float_t dpt = 0.01; // GeV
//
// For b = 0
-// (factor 1.35 to scale from 10% most central to b=0)
-//
- Float_t scaleC0 = 1.35 * ptUp / dpt;
- Float_t scaleB0 = 1.35 * ptUp / dpt;
- Float_t scaleD0 = 1.35 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6
+// (factor 1.28 to scale from 10% most central to b=0)
+//
+ Float_t aGlauber = 208. * 208. * 7.03e-4; // 1/mb
+
+ Float_t scaleC0 = aGlauber * ptUp / dpt;
+ Float_t scaleB0 = aGlauber * ptUp / dpt;
+ Float_t scaleD0 = 1.28 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6
//
//
//
// Heavy Flavors
//
- f = new TFile("HVQinc.root");
- TH1F* ptBB = (TH1F*) f->Get("hPtCorra");
- TH1F* ptCC = (TH1F*) f->Get("hpta");
TF1* ptBBLf = new TF1("ptBBLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 3.);
- ptBBLf->SetParameter(0, 4.46695e-03);
- ptBBLf->SetParameter(1, 1.60242e+00);
- ptBBLf->SetParameter(2, 2.24948e+00);
+ ptBBLf->SetParameter(0, 4.390e-05);
+ ptBBLf->SetParameter(1, 1.8706);
+ ptBBLf->SetParameter(2, 2.6623);
TF1* ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., ptUp);
- ptBBHf->SetParameter(0, 2.59961e-03);
- ptBBHf->SetParameter(1, 2.41);
- ptBBHf->SetParameter(2, 3.075);
+ ptBBHf->SetParameter(0, 2.5329e-05);
+ ptBBHf->SetParameter(1, 2.6067);
+ ptBBHf->SetParameter(2, 3.3821);
TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp);
- ptCCHf->SetParameter(0, 6.72360e-01);
- ptCCHf->SetParameter(1, 7.06043e-01);
- ptCCHf->SetParameter(2, 2.74240e+00);
- ptCCHf->SetParameter(3, 8.45018e-03);
+ ptCCHf->SetParameter(0, 4.8234e-03);
+ ptCCHf->SetParameter(1, 7.5656e-01);
+ ptCCHf->SetParameter(2, 2.7707e+00);
+ ptCCHf->SetParameter(3, 2.3362e-02);
TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5);
- ptCCLf->SetParameter(0, 1.40260e+00);
- ptCCLf->SetParameter(1, 3.75762e-01);
- ptCCLf->SetParameter(2, 1.54345e+00);
- ptCCLf->SetParameter(3, 2.49806e-01);
+ ptCCLf->SetParameter(0, 1.190e-02);
+ ptCCLf->SetParameter(1, 3.6343e-01);
+ ptCCLf->SetParameter(2, 1.4689e+00);
+ ptCCLf->SetParameter(3, 2.5772e-01);
TF1* ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp);
ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125);
//
// pi/K -> mu
//
- f->Close();
- f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root");
+ TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root");
TH2F* etaptPiK = (TH2F*) f->Get("etaptH");
TAxis* etaAxis = etaptPiK->GetXaxis();
TAxis* ptAxis = etaptPiK->GetYaxis();
TH1F* ptBH = new TH1F("ptBH", "pt Spectrum mu from b", 20, 0., 10.);
TH1F* ptDH = new TH1F("ptDH", "pt Spectrum mu from pi/K", 20, 0., 10.);
TH1F* ptBH2 = new TH1F("ptBH2", "pt Spectrum mu from b", 20, 0., 10.);
+ TH1F* costBBH = new TH1F("costBBH", "cos theta* ", 20, -1, 1.);
//
// Event Loop
//
Float_t deta = eta1 - eta2;
//
// invariant mass
- Float_t m2 = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi));
- Float_t m = TMath::Sqrt(m2);
+ Float_t m2 = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi));
+ Float_t m = TMath::Sqrt(m2);
+ Float_t mt2 = pT1 * pT1 + pT2 * pT2 + 2. * pT1 * pT2 * TMath::CosH(deta);
+ Float_t mt = TMath::Sqrt(mt2);
+ Float_t cost = 2. * pT1 * pT2 * TMath::SinH(deta) / m / mt;
//
// Smearing (to be improved)
acc->SetCharge(-1);
Float_t eff2 = eff->Evaluate(pT2, theta2, phid2);
Float_t acc2 = acc->Evaluate(pT2, theta2, phid2);
- Float_t tri2 = trigeff->Evaluate(1, pT2, theta2, phid2);
+ Float_t tri2 = trigeff->Evaluate(-1, pT2, theta2, phid2);
Float_t effA = eff1 * eff2 * acc1 * acc2 * tri1 * tri2;
Float_t ptMax = pT1;
Float_t ptMin = pT2;
if (pT2 > pT1) {
- Float_t ptMax = pT2;
- Float_t ptMin = pT1;
+ ptMax = pT2;
+ ptMin = pT1;
}
if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) {
massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA);
massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA);
massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA);
+
+ costBBH->Fill(cost, wgtB1 * wgtB2 / 4. * effA);
}
//
// pT - Spectra
ptCH->Draw("same");
TCanvas *c4 = new TCanvas("c4","Canvas 5",400,10,600,700);
ptBH2->Draw();
+
+ TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700);
+ costBBH->Draw();
}