]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FASTSIM/uncorrBg.C
Removing warnings in MONITOR
[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C
index 091c1e2001716fc33835d0213b6d8f37a29b451a..98f773ac57538f2d2d175ddda0f8cb90f6035513 100644 (file)
@@ -1,3 +1,20 @@
+#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.)
 {
 //
@@ -8,7 +25,7 @@ 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;
@@ -17,25 +34,27 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
     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
     
 //
 //
@@ -74,31 +93,28 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
 //
 //  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);
@@ -107,8 +123,7 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
 //
 //  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();    
@@ -125,6 +140,7 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
     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
 //
@@ -155,8 +171,11 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
        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)
@@ -230,15 +249,15 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
        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.) {
@@ -251,6 +270,8 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
            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
@@ -324,5 +345,8 @@ void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
     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();
  
 }