[u/mrichter/AliRoot.git] / FASTSIM / uncorrBg.C

void uncorrBg(Int_t nev = 1000000, Double_t bmin = 0., Double_t bmax = 5.)
{
//
//
//
    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.e24 * 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  = 1.;
    Float_t etamin    = 2.543;
    Float_t etar      = 1.457;
    Float_t ptUp      = 20.;   // GeV
    Float_t dpt       = 0.1;   // 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
    
//
//
//  Fast response
//    
    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(); 
    trigeff->Init();
    
/*
//  To be improved
//
    AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
    tracker->AddResponse(acc);
    tracker->AddResponse(eff);
    tracker->AddResponse(res);

    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");   
    
    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);

    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);

    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);

    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);
    
    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);
//
//  pi/K -> mu
//
    f->Close();
    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.);
    TH1F* massCCH = new TH1F("massCCH", "Mass Spectrum: c-c        ", 150, 0., 15.);
    TH1F* massBCH = new TH1F("massBCH", "Mass Spectrum: b-c        ", 150, 0., 15.);
    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.);    
//
// 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 = ptUp * gRandom->Rndm();
	Float_t pT2 = ptUp * gRandom->Rndm();
//
// phi
	Float_t phi1 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
	Float_t phi2 = 2. * TMath::Pi() * gRandom->Rndm() - TMath::Pi();
	Float_t dphi = phi1 - phi2;
//
// eta
	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);

//
// Smearing (to be improved)
	Float_t dm = m * 0.01;
	m += gRandom->Gaus(0., dm);	
//
// Weights
//
//      Heavy Flavour
//
	Int_t ibin;
	Float_t wgtB1, wgtB2;
	Float_t wgtC1, wgtC2;

	if (pT1 > 1.5) {
	    wgtC1 = ptCCHf->Eval(pT1) * scaleC;
	} else {
	    wgtC1 = ptCCLf->Eval(pT1) * scaleC;
	}
	if (pT2 > 1.5) {
	    wgtC2 = ptCCHf->Eval(pT2) * scaleC;
	} else {
	    wgtC2 = ptCCLf->Eval(pT2) * scaleC;
	}


	if (pT1 > 3.) {
	    wgtB1 = ptBBHf->Eval(pT1) * scaleB;
	} else {
	    wgtB1 = ptBBLf->Eval(pT1) * scaleB;
	}
	if (pT2 > 3.) {
	    wgtB2 = ptBBHf->Eval(pT2) * scaleB;
	} else {
	    wgtB2 = ptBBLf->Eval(pT2) * scaleB;
	}


//
//      Weight  for decays
//
	Int_t etaBin, ptBin;
	Float_t wgtD1, wgtD2;
	
	etaBin = etaAxis->FindBin(eta1);
	ptBin  = ptAxis ->FindBin(pT1);	
	wgtD1  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
	
	etaBin = etaAxis->FindBin(eta2);
	ptBin  = ptAxis ->FindBin(pT2);	
	wgtD2  = etaptPiK->GetBinContent(etaBin, ptBin) * scaleD;
	
//
//      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);
	Float_t eff1  = eff->Evaluate(pT1, theta1, phid1);
	Float_t acc1  = acc->Evaluate(pT1, theta1, phid1);
	Float_t tri1  = trigeff->Evaluate(1, pT1, theta1, phid1);
	res->SetCharge(-1);
	eff->SetCharge(-1);
	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 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;
	}
	
	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);
	    massBCH->Fill(m, wgtC2 * wgtB1 / 4. * effA);
	    massDDH->Fill(m, wgtD1 * wgtD2 / 4. * effA);
	    massBDH->Fill(m, wgtB1 * wgtD2 / 4. * effA);
	    massBDH->Fill(m, wgtB2 * wgtD1 / 4. * effA);
	    massCDH->Fill(m, wgtC1 * wgtD2 / 4. * effA);
	    massCDH->Fill(m, wgtC2 * wgtD1 / 4. * effA);
	}
	//
	// pT - Spectra
	//
	for (Int_t ipt = 0; ipt < 20; ipt++)
	{
	    Float_t pt = 0.5 * ipt;
	    ptBH2->Fill(pT1, wgtB1);		
	    if (pT1 > pt) {
		ptCH->Fill(pt, wgtC1);
		ptBH->Fill(pt, wgtB1);
		ptDH->Fill(pt, wgtD1);
	    }
	} // bins
    } // event loop

    Float_t evtWgt = events / Float_t(nev);
    
    massBBH->Scale(evtWgt);
    massCCH->Scale(evtWgt);
    massBCH->Scale(evtWgt);
    massDDH->Scale(evtWgt);
    massBDH->Scale(evtWgt);
    massCDH->Scale(evtWgt);
    
    TH1F * massALH = new TH1F(*massCDH);
    massALH->Add(massBDH);
    massALH->Add(massDDH);
    massALH->Add(massBCH);    
    massALH->Add(massCCH);
    massALH->Add(massBBH);     

    TCanvas *c0 = new TCanvas("c0","Canvas 1",400,10,600,700);
    massCCH->SetLineColor(4);
    massCCH->SetMinimum(1.);
    massCCH->SetMaximum(1.e4);
    massCCH->SetXTitle("m_{#mu#mu} [GeV]");
    massCCH->SetYTitle("Entries/100 MeV /10^{6} s");
    massCCH->Draw("");
    massALH->SetLineColor(3);
    massALH->Draw("same");
    massBBH->SetLineColor(6);
    massBBH->Draw("same");

    TCanvas *c2 = new TCanvas("c2","Canvas 3",400,10,600,700);
    massDDH->SetLineColor(2);
    massDDH->SetMinimum(1.e2);
    massDDH->SetMaximum(1.e6);
    massDDH->SetXTitle("m_{#mu#mu} [GeV]");
    massDDH->SetYTitle("Entries/100 MeV /10^{6} s");
    massDDH->Draw("");
    massALH->SetLineColor(3);
    massALH->Draw("same");
    massCCH->SetLineColor(4);
    massCCH->Draw("same");
    massBBH->SetLineColor(6);
    massBBH->Draw("same");

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