FASTSIM/uncorrBg.C

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