X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=FASTSIM%2FuncorrBg.C;h=98f773ac57538f2d2d175ddda0f8cb90f6035513;hb=9175f0df33470ddd543ecb35f0a08ff2624a5436;hp=fa08eb8093c0b7a591398ee8f9ea7038c642f071;hpb=d0768bbb4606d347dba0cb401ffe7f2036f0b33c;p=u%2Fmrichter%2FAliRoot.git diff --git a/FASTSIM/uncorrBg.C b/FASTSIM/uncorrBg.C index fa08eb8093c..98f773ac575 100644 --- a/FASTSIM/uncorrBg.C +++ b/FASTSIM/uncorrBg.C @@ -1,109 +1,133 @@ -void uncorrBg(Int_t nev = 1000000) +#if !defined(__CINT__) || defined(__MAKECINT__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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.) { -// b: 0.00261 (0.0375 per unit of rapidity per event) -// c: 1.029 (1.35 per unit of rapidity per event) // -// rate of central collisions 400 Hz -// central events 4 10^8 -// Get the single muon pT spectra // -// - Float_t scaleC = 200.; - Float_t scaleB = 200.; - Float_t scaleD = 0.34 * 1.362e-3 * 1.e3; - // - Float_t bbx[15] = - { - 0.21, 0.55, 0.65, 0.65, 0.5, - 0.40, 0.26, 0.19, 0.11, 0.1, - 0.075, 0.045, 0.035, 0.02, 0.017 - }; + AliFastGlauber* glauber = new AliFastGlauber(); + glauber->Init(1); + + Float_t lumi = 5.e26; // cm^-2 s^-1 + Float_t time = 1.e6; // s + 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("-------------------------------------------------------------------\n"); + 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("-------------------------------------------------------------------\n"); + +// +// + Float_t ptMinCut = 3.; // GeV + Float_t etamin = 2.543; + Float_t etar = 1.457; + Float_t ptUp = 20.; // GeV + Float_t dpt = 0.01; // GeV +// +// For b = 0 +// (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 + +// // // Fast response // - AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff(); - AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc(); - AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff(); - AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes(); + AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff(); + AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc(); + AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff(); + AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes(); acc->SetBackground(1.); eff->SetBackground(1.); res->SetBackground(1.); - acc->Init(); - eff->Init(); - res->Init(); + + acc ->Init(); + eff ->Init(); + res ->Init(); + trigeff->Init(); + +/* +// To be improved +// AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker"); - tracker->AddResponse(trigeff); tracker->AddResponse(acc); tracker->AddResponse(eff); tracker->AddResponse(res); - tracker->Init(); - tracker->Dump(); + AliFastDetector* trigger = new AliFastDetector("Trigger", "Muon Trigger"); + trigger->AddResponse(trigeff); + + AliFastDetector* spectro = new AliFastDetector("Spectro", "Muon Spectrometer"); + spectro->AddSubdetector(tracker, ""); + spectro->AddSubdetector(trigger, ""); + spectro->Init(); +*/ // // Heavy Flavors // - f = new TFile("HVQinc.root"); - TH1F* ptBB = (TH1F*) f->Get("hPtCorra"); - TH1F* ptCC = (TH1F*) f->Get("hpta"); - TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700); - + 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., 20.); - ptBBHf->SetParameter(0, 2.59961e-03); - ptBBHf->SetParameter(1, 2.41); - ptBBHf->SetParameter(2, 3.075); + TF1* ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., ptUp); + 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, 20.); - 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->Draw("ac"); + TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp); + 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->Draw("ac"); - /* - - TF1* ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.); - ptCCHf->SetParameter(0, 0.249); - ptCCHf->SetParameter(1, 1.15); - ptCCHf->SetParameter(2, 3.33); -// ptCCHf->Draw("ac"); - - TF1* ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.); - ptCCLf->SetParameter(0, 1.125); - ptCCLf->SetParameter(1, 0.525); - ptCCLf->SetParameter(2, 2.42); -// ptCCLf->Draw("ac"); -*/ - + 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., 20.); + TF1* ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp); ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125); ptBf->SetParameter(1, 6.27); ptBf->SetParameter(2, 3.2); -// ptBf->Draw("ac"); // // pi/K -> mu // - f->Close(); - f = new TFile("pikmu.root"); - TH2F* etaptPiK = (TH2F*) f->Get("etaptH"); - TAxis* etaAxis = etaptPiK->GetXaxis(); - TAxis* ptAxis = etaptPiK->GetYaxis(); - TH1F* ptPiK = (TH1F*) f->Get("ptH3"); -// ptAxis = ptPiK->GetXaxis(); - + TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root"); + TH2F* etaptPiK = (TH2F*) f->Get("etaptH"); + TAxis* etaAxis = etaptPiK->GetXaxis(); + TAxis* ptAxis = etaptPiK->GetYaxis(); +// // // Book histograms TH1F* massBBH = new TH1F("massBBH", "Mass Spectrum: b-b ", 150, 0., 15.); @@ -112,27 +136,29 @@ void uncorrBg(Int_t nev = 1000000) TH1F* massDDH = new TH1F("massDDH", "Mass Spectrum: decay-decay", 150, 0., 15.); TH1F* massBDH = new TH1F("massBDH", "Mass Spectrum: decay-b ", 150, 0., 15.); TH1F* massCDH = new TH1F("massCDH", "Mass Spectrum: decay-c ", 150, 0., 15.); - TH1F* ptCH = new TH1F("ptCH", "pt Spectrum mu from c", 20., 0., 10.); - 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* ptBH3 = new TH1F("ptBH3", "pt Spectrum mu from b", 15., 0., 15.); - for (Int_t i = 0; i < 15; i++) - { - ptBH3->SetBinContent(i+1, bbx[i]); - } - ptBH3->Draw(); - - + TH1F* ptCH = new TH1F("ptCH", "pt Spectrum mu from c", 20, 0., 10.); + 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 // Int_t iev; for (iev = 0; iev < nev; iev++) { // +// Collision geometry +// + Float_t b; + b = glauber->GetRandomImpactParameter(bmin, bmax); + Double_t nbinary = glauber->Binaries(b); + Float_t scaleC = scaleC0 * nbinary; + Float_t scaleB = scaleB0 * nbinary; + Float_t scaleD = scaleD0 * nbinary; +// // pT - Float_t pT1 = 20. * gRandom->Rndm(); - Float_t pT2 = 20. * gRandom->Rndm(); + Float_t pT1 = ptUp * gRandom->Rndm(); + Float_t pT2 = ptUp * gRandom->Rndm(); // // phi Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi(); @@ -140,16 +166,19 @@ void uncorrBg(Int_t nev = 1000000) Float_t dphi = phi1 - phi2; // // eta - Float_t eta1 = 1.457 * gRandom->Rndm() + 2.543; - Float_t eta2 = 1.457 * gRandom->Rndm() + 2.543; + Float_t eta1 = etar * gRandom->Rndm() + etamin; + Float_t eta2 = etar * gRandom->Rndm() + etamin; 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 +// Smearing (to be improved) Float_t dm = m * 0.01; m += gRandom->Gaus(0., dm); // @@ -198,22 +227,17 @@ void uncorrBg(Int_t nev = 1000000) etaBin = etaAxis->FindBin(eta2); ptBin = ptAxis ->FindBin(pT2); wgtD2 = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD; - - /* - ptBin = ptAxis ->FindBin(pT1); - wgtD1 = ptPiK->GetBinContent(ptBin) * scaleD; - ptBin = ptAxis ->FindBin(pT2); - wgtD2 = ptPiK->GetBinContent(ptBin) * scaleD; - */ - + // -// efficiencies +// Efficiencies // Float_t theta1 = 2. * TMath::ATan(TMath::Exp(-eta1)) * 180./TMath::Pi(); Float_t theta2 = 2. * TMath::ATan(TMath::Exp(-eta2)) * 180./TMath::Pi(); Float_t phid1 = phi1 * 180./TMath::Pi(); Float_t phid2 = phi2 * 180./TMath::Pi(); - + Float_t p1 = pT1/TMath::Sin(theta1 * TMath::Pi()/180.); + Float_t p2 = pT2/TMath::Sin(theta2 * TMath::Pi()/180.); + res->SetCharge(1); eff->SetCharge(1); acc->SetCharge(1); @@ -225,24 +249,18 @@ void uncorrBg(Int_t nev = 1000000) 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 ( -// (ptMax > 6. && ptMin > 3.) || -// (ptMax < 6. && ptMin > (6. - 0.5 * ptMax)) -// ) -// { - if (ptMin > 3.) { + if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) { massBBH->Fill(m, wgtB1 * wgtB2 / 4. * effA); massCCH->Fill(m, wgtC1 * wgtC2 / 4. * effA); massBCH->Fill(m, wgtC1 * wgtB2 / 4. * effA); @@ -252,6 +270,8 @@ void uncorrBg(Int_t nev = 1000000) 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 @@ -265,12 +285,10 @@ void uncorrBg(Int_t nev = 1000000) ptBH->Fill(pt, wgtB1); ptDH->Fill(pt, wgtD1); } - } - + } // bins } // event loop -// Float_t eff = 0.6 * 1.0; - Float_t evtWgt = 1. / Float_t(nev) * 4.e8; - + + Float_t evtWgt = events / Float_t(nev); massBBH->Scale(evtWgt); massCCH->Scale(evtWgt); @@ -327,5 +345,8 @@ void uncorrBg(Int_t nev = 1000000) 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(); }