FASTSIM/uncorrBg.C

   1 void uncorrBg(Int_t nev = 1000000)
   2 {
   3 // b: 0.00261     (0.0375 per unit of rapidity per event)
   4 // c: 1.029       (1.35 per unit of rapidity per event)
   5 //
   6 // rate of central collisions  400 Hz
   7 // central events 4 10^8
   8 // Get the single muon pT spectra
   9 //
  10 //
  11     Float_t scaleC = 200.;
  12     Float_t scaleB = 200.;
  13     Float_t scaleD = 0.34 * 1.362e-3 * 1.e3;
  14
  15 //
  16     Float_t bbx[15] =
  17         {
  18             0.21, 0.55, 0.65, 0.65,  0.5,
  19             0.40, 0.26, 0.19, 0.11,  0.1,
  20             0.075, 0.045, 0.035, 0.02, 0.017
  21         };
  22 //
  23 //  Fast response
  24 //
  25     AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
  26     AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
  27     AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
  28     AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
  29     acc->SetBackground(1.);
  30     eff->SetBackground(1.);
  31     res->SetBackground(1.);
  32     acc->Init();
  33     eff->Init();
  34     res->Init();
  35     AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
  36     tracker->AddResponse(trigeff);
  37     tracker->AddResponse(acc);
  38     tracker->AddResponse(eff);
  39     tracker->AddResponse(res);
  40     tracker->Init();
  41     tracker->Dump();
  42
  43
  44 //
  45 //  Heavy Flavors
  46 //
  47     f = new TFile("HVQinc.root");
  48     TH1F* ptBB = (TH1F*) f->Get("hPtCorra");
  49     TH1F* ptCC = (TH1F*) f->Get("hpta");
  50     TCanvas *c5 = new TCanvas("c5","Canvas 6",400,10,600,700);
  51
  52     TF1*  ptBBLf = new TF1("ptBBLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 3.);
  53     ptBBLf->SetParameter(0, 4.46695e-03);
  54     ptBBLf->SetParameter(1, 1.60242e+00);
  55     ptBBLf->SetParameter(2, 2.24948e+00);
  56
  57     TF1*  ptBBHf = new TF1("ptBBHf", "[0] * x / (1. + (x/[1])**2)**[2]", 3., 20.);
  58     ptBBHf->SetParameter(0, 2.59961e-03);
  59     ptBBHf->SetParameter(1, 2.41);
  60     ptBBHf->SetParameter(2, 3.075);
  61
  62     TF1*  ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 1.5, 20.);
  63     ptCCHf->SetParameter(0, 6.72360e-01);
  64     ptCCHf->SetParameter(1, 7.06043e-01);
  65     ptCCHf->SetParameter(2, 2.74240e+00);
  66     ptCCHf->SetParameter(3, 8.45018e-03);
  67 //    ptCCHf->Draw("ac");
  68
  69     TF1*  ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**([2] + [3] * x)", 0., 1.5);
  70     ptCCLf->SetParameter(0, 1.40260e+00);
  71     ptCCLf->SetParameter(1, 3.75762e-01);
  72     ptCCLf->SetParameter(2, 1.54345e+00);
  73     ptCCLf->SetParameter(3, 2.49806e-01);
  74 //    ptCCLf->Draw("ac");
  75     /*
  76
  77     TF1*  ptCCHf = new TF1("ptCCHf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
  78     ptCCHf->SetParameter(0, 0.249);
  79     ptCCHf->SetParameter(1, 1.15);
  80     ptCCHf->SetParameter(2, 3.33);
  81 //    ptCCHf->Draw("ac");
  82
  83     TF1*  ptCCLf = new TF1("ptCCLf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
  84     ptCCLf->SetParameter(0, 1.125);
  85     ptCCLf->SetParameter(1, 0.525);
  86     ptCCLf->SetParameter(2, 2.42);
  87 //    ptCCLf->Draw("ac");
  88 */
  89
  90
  91     TF1*  ptBf = new TF1("ptBf", "[0] * x / (1. + (x/[1])**2)**[2]", 0., 20.);
  92     ptBf->SetParameter(0, 1.e5 * 0.7 * 1.125);
  93     ptBf->SetParameter(1, 6.27);
  94     ptBf->SetParameter(2, 3.2);
  95 //    ptBf->Draw("ac");
  96 //
  97 //  pi/K -> mu
  98 //
  99     f->Close();
 100     f = new TFile("pikmu.root");
 101     TH2F* etaptPiK = (TH2F*) f->Get("etaptH");
 102     TAxis* etaAxis = etaptPiK->GetXaxis();
 103     TAxis* ptAxis  = etaptPiK->GetYaxis();
 104     TH1F* ptPiK    = (TH1F*) f->Get("ptH3");
 105 //    ptAxis = ptPiK->GetXaxis();
 106
 107 //
 108 // Book histograms
 109     TH1F* massBBH = new TH1F("massBBH", "Mass Spectrum: b-b        ", 150, 0., 15.);
 110     TH1F* massCCH = new TH1F("massCCH", "Mass Spectrum: c-c        ", 150, 0., 15.);
 111     TH1F* massBCH = new TH1F("massBCH", "Mass Spectrum: b-c        ", 150, 0., 15.);
 112     TH1F* massDDH = new TH1F("massDDH", "Mass Spectrum: decay-decay", 150, 0., 15.);
 113     TH1F* massBDH = new TH1F("massBDH", "Mass Spectrum: decay-b    ", 150, 0., 15.);
 114     TH1F* massCDH = new TH1F("massCDH", "Mass Spectrum: decay-c    ", 150, 0., 15.);
 115     TH1F* ptCH    = new TH1F("ptCH", "pt Spectrum mu from c", 20., 0., 10.);
 116     TH1F* ptBH    = new TH1F("ptBH", "pt Spectrum mu from b", 20., 0., 10.);
 117     TH1F* ptDH    = new TH1F("ptDH", "pt Spectrum mu from pi/K", 20., 0., 10.);
 118     TH1F* ptBH2    = new TH1F("ptBH2", "pt Spectrum mu from b", 20., 0., 10.);
 119     TH1F* ptBH3    = new TH1F("ptBH3", "pt Spectrum mu from b", 15., 0., 15.);
 120     for (Int_t i = 0; i < 15; i++)
 121     {
 122         ptBH3->SetBinContent(i+1, bbx[i]);
 123     }
 124     ptBH3->Draw();
 125
 126
 127 //
 128 // Event Loop
 129 //
 130     Int_t iev;
 131     for (iev = 0; iev < nev; iev++) {
 132 //
 133 // pT
 134         Float_t pT1 = 20. * gRandom->Rndm();
 135         Float_t pT2 = 20. * gRandom->Rndm();
 136 //
 137 // phi
 138         Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
 139         Float_t phi2 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
 140         Float_t dphi = phi1 - phi2;
 141 //
 142 // eta
 143         Float_t eta1 = 1.457 * gRandom->Rndm() + 2.543;
 144         Float_t eta2 = 1.457 * gRandom->Rndm() + 2.543;
 145         Float_t deta = eta1 - eta2;
 146 //
 147 // invariant mass
 148         Float_t m2 = 2. * pT1 * pT2 * (TMath::CosH(deta) - TMath::Cos(dphi));
 149         Float_t m  = TMath::Sqrt(m2);
 150
 151 //
 152 // Smearing
 153         Float_t dm = m * 0.01;
 154         m += gRandom->Gaus(0., dm);
 155 //
 156 // Weights
 157 //
 158 //      Heavy Flavour
 159 //
 160         Int_t ibin;
 161         Float_t wgtB1, wgtB2;
 162         Float_t wgtC1, wgtC2;
 163
 164         if (pT1 > 1.5) {
 165             wgtC1 = ptCCHf->Eval(pT1) * scaleC;
 166         } else {
 167             wgtC1 = ptCCLf->Eval(pT1) * scaleC;
 168         }
 169         if (pT2 > 1.5) {
 170             wgtC2 = ptCCHf->Eval(pT2) * scaleC;
 171         } else {
 172             wgtC2 = ptCCLf->Eval(pT2) * scaleC;
 173         }
 174
 175
 176         if (pT1 > 3.) {
 177             wgtB1 = ptBBHf->Eval(pT1) * scaleB;
 178         } else {
 179             wgtB1 = ptBBLf->Eval(pT1) * scaleB;
 180         }
 181         if (pT2 > 3.) {
 182             wgtB2 = ptBBHf->Eval(pT2) * scaleB;
 183         } else {
 184             wgtB2 = ptBBLf->Eval(pT2) * scaleB;
 185         }
 186
 187
 188 //
 189 //      Weight  for decays
 190 //
 191         Int_t etaBin, ptBin;
 192         Float_t wgtD1, wgtD2;
 193
 194         etaBin = etaAxis->FindBin(eta1);
 195         ptBin  = ptAxis ->FindBin(pT1);
 196         wgtD1  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
 197
 198         etaBin = etaAxis->FindBin(eta2);
 199         ptBin  = ptAxis ->FindBin(pT2);
 200         wgtD2  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
 201
 202         /*
 203         ptBin  = ptAxis ->FindBin(pT1);
 204         wgtD1  = ptPiK->GetBinContent(ptBin) * scaleD;
 205         ptBin  = ptAxis ->FindBin(pT2);
 206         wgtD2  = ptPiK->GetBinContent(ptBin) * scaleD;
 207         */
 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
 217         res->SetCharge(1);
 218         eff->SetCharge(1);
 219         acc->SetCharge(1);
 220         Float_t eff1  = eff->Evaluate(pT1, theta1, phid1);
 221         Float_t acc1  = acc->Evaluate(pT1, theta1, phid1);
 222         Float_t tri1  = trigeff->Evaluate(1, pT1, theta1, phid1);
 223         res->SetCharge(-1);
 224         eff->SetCharge(-1);
 225         acc->SetCharge(-1);
 226         Float_t eff2  = eff->Evaluate(pT2, theta2, phid2);
 227         Float_t acc2  = acc->Evaluate(pT2, theta2, phid2);
 228         Float_t tri2  = trigeff->Evaluate(1, pT2, theta2, phid2);
 229
 230         Float_t effA   = eff1 * eff2 * acc1 * acc2 * tri1 * tri2;
 231
 232         Float_t ptMax = pT1;
 233         Float_t ptMin = pT2;
 234         if (pT2 > pT1) {
 235             Float_t ptMax = pT2;
 236             Float_t ptMin = pT1;
 237         }
 238
 239
 240 //      if (
 241 //          (ptMax > 6. && ptMin > 3.) ||
 242 //          (ptMax < 6. && ptMin > (6. - 0.5 * ptMax))
 243 //          )
 244 //      {
 245         if (ptMin > 3.) {
 246             massBBH->Fill(m, wgtB1 * wgtB2 / 4. * effA);
 247             massCCH->Fill(m, wgtC1 * wgtC2 / 4. * effA);
 248             massBCH->Fill(m, wgtC1 * wgtB2 / 4. * effA);
 249             massBCH->Fill(m, wgtC2 * wgtB1 / 4. * effA);
 250             massDDH->Fill(m, wgtD1 * wgtD2 / 4. * effA);
 251             massBDH->Fill(m, wgtB1 * wgtD2 / 4. * effA);
 252             massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA);
 253             massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA);
 254             massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA);
 255         }
 256         //
 257         // pT - Spectra
 258         //
 259         for (Int_t ipt = 0; ipt < 20; ipt++)
 260         {
 261             Float_t pt = 0.5 * ipt;
 262             ptBH2->Fill(pT1, wgtB1);
 263             if (pT1 > pt) {
 264                 ptCH->Fill(pt, wgtC1);
 265                 ptBH->Fill(pt, wgtB1);
 266                 ptDH->Fill(pt, wgtD1);
 267             }
 268         }
 269
 270     } // event loop
 271 //    Float_t eff    = 0.6 * 1.0;
 272     Float_t evtWgt = 1. / Float_t(nev) * 4.e8;
 273
 274
 275     massBBH->Scale(evtWgt);
 276     massCCH->Scale(evtWgt);
 277     massBCH->Scale(evtWgt);
 278     massDDH->Scale(evtWgt);
 279     massBDH->Scale(evtWgt);
 280     massCDH->Scale(evtWgt);
 281
 282     TH1F * massALH = new TH1F(*massCDH);
 283     massALH->Add(massBDH);
 284     massALH->Add(massDDH);
 285     massALH->Add(massBCH);
 286     massALH->Add(massCCH);
 287     massALH->Add(massBBH);
 288
 289     TCanvas *c0 = new TCanvas("c0","Canvas 1",400,10,600,700);
 290     massCCH->SetLineColor(4);
 291     massCCH->SetMinimum(1.);
 292     massCCH->SetMaximum(1.e4);
 293     massCCH->SetXTitle("m_{#mu#mu} [GeV]");
 294     massCCH->SetYTitle("Entries/100 MeV /10^{6} s");
 295     massCCH->Draw("");
 296     massALH->SetLineColor(3);
 297     massALH->Draw("same");
 298     massBBH->SetLineColor(6);
 299     massBBH->Draw("same");
 300
 301     TCanvas *c2 = new TCanvas("c2","Canvas 3",400,10,600,700);
 302     massDDH->SetLineColor(2);
 303     massDDH->SetMinimum(1.e2);
 304     massDDH->SetMaximum(1.e6);
 305     massDDH->SetXTitle("m_{#mu#mu} [GeV]");
 306     massDDH->SetYTitle("Entries/100 MeV /10^{6} s");
 307     massDDH->Draw("");
 308     massALH->SetLineColor(3);
 309     massALH->Draw("same");
 310     massCCH->SetLineColor(4);
 311     massCCH->Draw("same");
 312     massBBH->SetLineColor(6);
 313     massBBH->Draw("same");
 314
 315     TCanvas *c3 = new TCanvas("c3","Canvas 4",400,10,600,700);
 316     ptCH->Scale(1./float(nev));
 317     ptBH->Scale(1./float(nev));
 318     ptDH->Scale(1./float(nev));
 319     ptCH->SetLineColor(4);
 320     ptBH->SetLineColor(6);
 321     ptDH->SetLineColor(2);
 322     ptCH->SetXTitle("p_{T}^{min} [GeV]");
 323     ptCH->SetYTitle("<n>_{#mu}/event");
 324
 325     ptDH->Draw();
 326     ptBH->Draw("same");
 327     ptCH->Draw("same");
 328     TCanvas *c4 = new TCanvas("c4","Canvas 5",400,10,600,700);
 329     ptBH2->Draw();
 330
 331 }