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  = 3.;
  28     Float_t etamin    = 2.543;
  29     Float_t etar      = 1.457;
  30     Float_t ptUp      = 20.;   // GeV
  31     Float_t dpt       = 0.01;   // GeV
  32 //
  33 //  For b = 0
  34 //  (factor 1.28 to scale from 10% most central to b=0)
  35 //
  36     Float_t aGlauber = 208. * 208. * 7.03e-4;    // 1/mb
  37
  38     Float_t scaleC0 = aGlauber * ptUp / dpt;
  39     Float_t scaleB0 = aGlauber * ptUp / dpt;
  40     Float_t scaleD0 = 1.28 * etar * ptUp / 1.35; // scaled by 1.35 to match ALICE-INT-2002-6
  41
  42 //
  43 //
  44 //  Fast response
  45 //
  46     AliFastMuonTriggerEff  *trigeff = new AliFastMuonTriggerEff();
  47     AliFastMuonTrackingAcc *acc     = new AliFastMuonTrackingAcc();
  48     AliFastMuonTrackingEff *eff     = new AliFastMuonTrackingEff();
  49     AliFastMuonTrackingRes *res     = new AliFastMuonTrackingRes();
  50     acc->SetBackground(1.);
  51     eff->SetBackground(1.);
  52     res->SetBackground(1.);
  53
  54     acc    ->Init();
  55     eff    ->Init();
  56     res    ->Init();
  57     trigeff->Init();
  58
  59 /*
  60 //  To be improved
  61 //
  62     AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
  63     tracker->AddResponse(acc);
  64     tracker->AddResponse(eff);
  65     tracker->AddResponse(res);
  66
  67     AliFastDetector* trigger = new AliFastDetector("Trigger", "Muon Trigger");
  68     trigger->AddResponse(trigeff);
  69
  70     AliFastDetector* spectro = new AliFastDetector("Spectro", "Muon Spectrometer");
  71     spectro->AddSubdetector(tracker, "");
  72     spectro->AddSubdetector(trigger, "");
  73     spectro->Init();
  74 */
  75
  76 //
  77 //  Heavy Flavors
  78 //
  79
  80     TF1*  ptBBLf = new TF1("ptBBLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 3.);
  81     ptBBLf->SetParameter(0, 4.390e-05);
  82     ptBBLf->SetParameter(1, 1.8706);
  83     ptBBLf->SetParameter(2, 2.6623);
  84
  85     TF1*  ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., ptUp);
  86     ptBBHf->SetParameter(0, 2.5329e-05);
  87     ptBBHf->SetParameter(1, 2.6067);
  88     ptBBHf->SetParameter(2, 3.3821);
  89
  90     TF1*  ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, ptUp);
  91     ptCCHf->SetParameter(0, 4.8234e-03);
  92     ptCCHf->SetParameter(1, 7.5656e-01);
  93     ptCCHf->SetParameter(2, 2.7707e+00);
  94     ptCCHf->SetParameter(3, 2.3362e-02);
  95
  96     TF1*  ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5);
  97     ptCCLf->SetParameter(0, 1.190e-02);
  98     ptCCLf->SetParameter(1, 3.6343e-01);
  99     ptCCLf->SetParameter(2, 1.4689e+00);
 100     ptCCLf->SetParameter(3, 2.5772e-01);
 101
 102     TF1*  ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., ptUp);
 103     ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125);
 104     ptBf->SetParameter(1, 6.27);
 105     ptBf->SetParameter(2, 3.2);
 106 //
 107 //  pi/K -> mu
 108 //
 109     f = new TFile("$(ALICE_ROOT)/FASTSIM/data/pikmu.root");
 110     TH2F*  etaptPiK = (TH2F*) f->Get("etaptH");
 111     TAxis* etaAxis  = etaptPiK->GetXaxis();
 112     TAxis* ptAxis   = etaptPiK->GetYaxis();
 113 //
 114 //
 115 // Book histograms
 116     TH1F* massBBH = new TH1F("massBBH", "Mass Spectrum: b-b        ", 150, 0., 15.);
 117     TH1F* massCCH = new TH1F("massCCH", "Mass Spectrum: c-c        ", 150, 0., 15.);
 118     TH1F* massBCH = new TH1F("massBCH", "Mass Spectrum: b-c        ", 150, 0., 15.);
 119     TH1F* massDDH = new TH1F("massDDH", "Mass Spectrum: decay-decay", 150, 0., 15.);
 120     TH1F* massBDH = new TH1F("massBDH", "Mass Spectrum: decay-b    ", 150, 0., 15.);
 121     TH1F* massCDH = new TH1F("massCDH", "Mass Spectrum: decay-c    ", 150, 0., 15.);
 122     TH1F* ptCH    = new TH1F("ptCH", "pt Spectrum mu from c",          20, 0., 10.);
 123     TH1F* ptBH    = new TH1F("ptBH", "pt Spectrum mu from b",          20, 0., 10.);
 124     TH1F* ptDH    = new TH1F("ptDH", "pt Spectrum mu from pi/K",       20, 0., 10.);
 125     TH1F* ptBH2   = new TH1F("ptBH2", "pt Spectrum mu from b",         20, 0., 10.);
 126     TH1F* costBBH = new TH1F("costBBH", "cos theta*    ", 20, -1, 1.);
 127 //
 128 // Event Loop
 129 //
 130     Int_t iev;
 131     for (iev = 0; iev < nev; iev++) {
 132 //
 133 //  Collision geometry
 134 //
 135         Float_t b;
 136         b = glauber->GetRandomImpactParameter(bmin, bmax);
 137         Double_t nbinary = glauber->Binaries(b);
 138         Float_t  scaleC  = scaleC0 * nbinary;
 139         Float_t  scaleB  = scaleB0 * nbinary;
 140         Float_t  scaleD  = scaleD0 * nbinary;
 141 //
 142 // pT
 143         Float_t pT1 = ptUp * gRandom->Rndm();
 144         Float_t pT2 = ptUp * gRandom->Rndm();
 145 //
 146 // phi
 147         Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
 148         Float_t phi2 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
 149         Float_t dphi = phi1 - phi2;
 150 //
 151 // eta
 152         Float_t eta1 = etar * gRandom->Rndm() + etamin;
 153         Float_t eta2 = etar * gRandom->Rndm() + etamin;
 154         Float_t deta = eta1 - eta2;
 155 //
 156 // invariant mass
 157         Float_t m2   = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi));
 158         Float_t m    = TMath::Sqrt(m2);
 159         Float_t mt2  = pT1 * pT1 + pT2 * pT2 + 2. * pT1 * pT2 * TMath::CosH(deta);
 160         Float_t mt   = TMath::Sqrt(mt2);
 161         Float_t cost = 2. * pT1 * pT2 * TMath::SinH(deta) / m / mt;
 162
 163 //
 164 // Smearing (to be improved)
 165         Float_t dm = m * 0.01;
 166         m += gRandom->Gaus(0., dm);
 167 //
 168 // Weights
 169 //
 170 //      Heavy Flavour
 171 //
 172         Int_t ibin;
 173         Float_t wgtB1, wgtB2;
 174         Float_t wgtC1, wgtC2;
 175
 176         if (pT1 > 1.5) {
 177             wgtC1 = ptCCHf->Eval(pT1) * scaleC;
 178         } else {
 179             wgtC1 = ptCCLf->Eval(pT1) * scaleC;
 180         }
 181         if (pT2 > 1.5) {
 182             wgtC2 = ptCCHf->Eval(pT2) * scaleC;
 183         } else {
 184             wgtC2 = ptCCLf->Eval(pT2) * scaleC;
 185         }
 186
 187
 188         if (pT1 > 3.) {
 189             wgtB1 = ptBBHf->Eval(pT1) * scaleB;
 190         } else {
 191             wgtB1 = ptBBLf->Eval(pT1) * scaleB;
 192         }
 193         if (pT2 > 3.) {
 194             wgtB2 = ptBBHf->Eval(pT2) * scaleB;
 195         } else {
 196             wgtB2 = ptBBLf->Eval(pT2) * scaleB;
 197         }
 198
 199
 200 //
 201 //      Weight  for decays
 202 //
 203         Int_t etaBin, ptBin;
 204         Float_t wgtD1, wgtD2;
 205
 206         etaBin = etaAxis->FindBin(eta1);
 207         ptBin  = ptAxis ->FindBin(pT1);
 208         wgtD1  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
 209
 210         etaBin = etaAxis->FindBin(eta2);
 211         ptBin  = ptAxis ->FindBin(pT2);
 212         wgtD2  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
 213
 214 //
 215 //      Efficiencies
 216 //
 217         Float_t theta1 = 2. * TMath::ATan(TMath::Exp(-eta1)) * 180./TMath::Pi();
 218         Float_t theta2 = 2. * TMath::ATan(TMath::Exp(-eta2)) * 180./TMath::Pi();
 219         Float_t phid1  = phi1 * 180./TMath::Pi();
 220         Float_t phid2  = phi2 * 180./TMath::Pi();
 221         Float_t p1     = pT1/TMath::Sin(theta1 * TMath::Pi()/180.);
 222         Float_t p2     = pT2/TMath::Sin(theta2 * TMath::Pi()/180.);
 223
 224         res->SetCharge(1);
 225         eff->SetCharge(1);
 226         acc->SetCharge(1);
 227         Float_t eff1  = eff->Evaluate(pT1, theta1, phid1);
 228         Float_t acc1  = acc->Evaluate(pT1, theta1, phid1);
 229         Float_t tri1  = trigeff->Evaluate(1, pT1, theta1, phid1);
 230         res->SetCharge(-1);
 231         eff->SetCharge(-1);
 232         acc->SetCharge(-1);
 233         Float_t eff2  = eff->Evaluate(pT2, theta2, phid2);
 234         Float_t acc2  = acc->Evaluate(pT2, theta2, phid2);
 235         Float_t tri2  = trigeff->Evaluate(-1, pT2, theta2, phid2);
 236
 237         Float_t effA   = eff1 * eff2 * acc1 * acc2 * tri1 * tri2;
 238
 239         Float_t ptMax = pT1;
 240         Float_t ptMin = pT2;
 241         if (pT2 > pT1) {
 242             Float_t ptMax = pT2;
 243             Float_t ptMin = pT1;
 244         }
 245
 246         if (ptMin > ptMinCut && p1 > 4. && p2 > 4.) {
 247             massBBH->Fill(m, wgtB1 * wgtB2 / 4. * effA);
 248             massCCH->Fill(m, wgtC1 * wgtC2 / 4. * effA);
 249             massBCH->Fill(m, wgtC1 * wgtB2 / 4. * effA);
 250             massBCH->Fill(m, wgtC2 * wgtB1 / 4. * effA);
 251             massDDH->Fill(m, wgtD1 * wgtD2 / 4. * effA);
 252             massBDH->Fill(m, wgtB1 * wgtD2 / 4. * effA);
 253             massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA);
 254             massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA);
 255             massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA);
 256
 257             costBBH->Fill(cost, wgtB1 * wgtB2 / 4. * effA);
 258         }
 259         //
 260         // pT - Spectra
 261         //
 262         for (Int_t ipt = 0; ipt < 20; ipt++)
 263         {
 264             Float_t pt = 0.5 * ipt;
 265             ptBH2->Fill(pT1, wgtB1);
 266             if (pT1 > pt) {
 267                 ptCH->Fill(pt, wgtC1);
 268                 ptBH->Fill(pt, wgtB1);
 269                 ptDH->Fill(pt, wgtD1);
 270             }
 271         } // bins
 272     } // event loop
 273
 274     Float_t evtWgt = events / Float_t(nev);
 275
 276     massBBH->Scale(evtWgt);
 277     massCCH->Scale(evtWgt);
 278     massBCH->Scale(evtWgt);
 279     massDDH->Scale(evtWgt);
 280     massBDH->Scale(evtWgt);
 281     massCDH->Scale(evtWgt);
 282
 283     TH1F * massALH = new TH1F(*massCDH);
 284     massALH->Add(massBDH);
 285     massALH->Add(massDDH);
 286     massALH->Add(massBCH);
 287     massALH->Add(massCCH);
 288     massALH->Add(massBBH);
 289
 290     TCanvas *c0 = new TCanvas("c0","Canvas 1",400,10,600,700);
 291     massCCH->SetLineColor(4);
 292     massCCH->SetMinimum(1.);
 293     massCCH->SetMaximum(1.e4);
 294     massCCH->SetXTitle("m_{#mu#mu} [GeV]");
 295     massCCH->SetYTitle("Entries/100 MeV /10^{6} s");
 296     massCCH->Draw("");
 297     massALH->SetLineColor(3);
 298     massALH->Draw("same");
 299     massBBH->SetLineColor(6);
 300     massBBH->Draw("same");
 301
 302     TCanvas *c2 = new TCanvas("c2","Canvas 3",400,10,600,700);
 303     massDDH->SetLineColor(2);
 304     massDDH->SetMinimum(1.e2);
 305     massDDH->SetMaximum(1.e6);
 306     massDDH->SetXTitle("m_{#mu#mu} [GeV]");
 307     massDDH->SetYTitle("Entries/100 MeV /10^{6} s");
 308     massDDH->Draw("");
 309     massALH->SetLineColor(3);
 310     massALH->Draw("same");
 311     massCCH->SetLineColor(4);
 312     massCCH->Draw("same");
 313     massBBH->SetLineColor(6);
 314     massBBH->Draw("same");
 315
 316     TCanvas *c3 = new TCanvas("c3","Canvas 4",400,10,600,700);
 317     ptCH->Scale(1./float(nev));
 318     ptBH->Scale(1./float(nev));
 319     ptDH->Scale(1./float(nev));
 320     ptCH->SetLineColor(4);
 321     ptBH->SetLineColor(6);
 322     ptDH->SetLineColor(2);
 323     ptCH->SetXTitle("p_{T}^{min} [GeV]");
 324     ptCH->SetYTitle("<n>_{#mu}/event");
 325
 326     ptDH->Draw();
 327     ptBH->Draw("same");
 328     ptCH->Draw("same");
 329     TCanvas *c4 = new TCanvas("c4","Canvas 5",400,10,600,700);
 330     ptBH2->Draw();
 331
 332     TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700);
 333     costBBH->Draw();
 334
 335 }