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

AliGenerator*  CreateGenerator();

void fastGenQuarkonia(Int_t nev = 1, char* filename = "galice.root")
{
//  Update data base
    AliPDG::AddParticlesToPdgDataBase();
    
// Create the fast tracker
    AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
    AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
    AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
    AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
    acc->SetBackground(0);
    eff->SetBackground(0);
    res->SetBackground(0);  
    acc->Init(); 
    eff->Init(); 
    res->Init(); 
    AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
    tracker->AddResponse(trigeff);
    tracker->AddResponse(acc);
    tracker->AddResponse(eff);
    tracker->AddResponse(res);
    tracker->Init();
    tracker->Dump();
//  Histos
    TH1F* massHU = new TH1F("massHU", "Mass Spectrum:                ", 500, 5, 15.);
    TH1F* massHS = new TH1F("massHS", "Mass Spectrum Smeared:        ", 500, 5, 15.);
    TH2F* etaH   = new TH2F("etaH",   "eta vs etas        ", 100, 0., 5., 100, 0., 5.);
    TH1F* etadH  = new TH1F("etaH",   "eta vs etas        ", 100, -1., 1.);
//
//                        Construction
//
    AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
    
    rl->SetCompressionLevel(2);
    rl->SetNumberOfEventsPerFile(nev);
    rl->LoadKinematics("RECREATE");
    rl->MakeTree("E");
    gAlice->SetRunLoader(rl);

//  Create stack
    rl->MakeStack();
    AliStack* stack      = rl->Stack();
 
//  Header
    AliHeader* header = rl->GetHeader();
//
//  Create and Initialize Generator
    AliGenerator *gener = CreateGenerator();
    gener->Init();
    gener->SetStack(stack);
    
//
//                        Event Loop
//
    Int_t iev;
     
    for (iev = 0; iev < nev; iev++) {

	printf("\n \n Event number %d \n \n", iev);
	

//  Initialize event
	header->Reset(0,iev);
	rl->SetEventNumber(iev);
	stack->Reset();
	rl->MakeTree("K");
//  Generate event
	gener->Generate();
		
//  Analysis
	Int_t npart = stack->GetNprimary();
	for (Int_t part=0; part<npart; part++) {
	    TParticle *MPart = stack->Particle(part);
	    Int_t mpart  = MPart->GetPdgCode();
	    if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
	    Int_t ch1 = MPart->GetFirstDaughter();
	    Int_t ch2 = MPart->GetLastDaughter();
	    
	    TParticle *muon1 = stack->Particle(ch1);
	    TParticle *muon2 = stack->Particle(ch2);
	    
	    Float_t theta1 = muon1->Theta();
	    Float_t thetad1 = theta1 * 180./TMath::Pi();
	    Float_t eta1   = muon1->Eta();
	    Float_t pt1    = muon1->Pt();
	    Float_t phi1   = muon1->Phi();
	    Float_t phid1  = phi1 * 180./TMath::Pi() - 180.;
	    Float_t p1     = muon1->P();

	    if (thetad1 > 9. || thetad1 < 2.) continue;

	    Float_t theta2 = muon2->Theta();
	    Float_t thetad2 = theta2 * 180./TMath::Pi();
	    Float_t eta2   = muon2->Eta();
	    Float_t pt2    = muon2->Pt();
	    Float_t phi2   = muon2->Phi();
	    Float_t phid2  = phi2 * 180./TMath::Pi() - 180.;
	    Float_t p2     = muon2->P();


	    if (thetad2 > 9. || thetad2 < 2.) continue;	    
	    
	    Float_t dphi   = phi1 - phi2;
	    Float_t deta   = eta1 - eta2;
	    Float_t m      =  TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));

	    massHU->Fill(m);
// Smear
	    // the mu+
	    Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
	    Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH;
	    
	    res->SetCharge(1);
	    eff->SetCharge(1);
	    acc->SetCharge(1);
	    res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);

	    
	    Float_t effp     = eff->Evaluate(pt1, thetad1, phid1);
	    Float_t accp     = acc->Evaluate(pt1, thetad1, phid1);
	    trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH);
	    thetas1 *= TMath::Pi()/180.;
	    phis1 *= TMath::Pi()/180.;
	    Float_t etas;
	    
	    if (TMath::Tan(thetas1/2.) > 0.) {
		etas  = -TMath::Log(TMath::Tan(thetas1/2.));
		etaH->Fill(etas, eta1);
		etadH->Fill(etas-eta1);
	    }
	    
	    // the mu- 
	    res->SetCharge(-1);
	    acc->SetCharge(-1);
	    eff->SetCharge(-1);
	    res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
	    Float_t effn     = eff->Evaluate(pt2, thetad2, phid2);
	    Float_t accn     = acc->Evaluate(pt2, thetad2, phid2);
	    trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH);
	    thetas2 *= TMath::Pi()/180.;
	    phis2 *= TMath::Pi()/180.;
	    if (TMath::Tan(thetas2/2.) > 0) {
		etas  = -TMath::Log(TMath::Tan(thetas2/2.));
		etaH->Fill(etas, eta2);
		etadH->Fill(etas-eta2);
	    }
	    
	    dphi   = phis1 - phis2;
	    etas1  = - TMath::Log(TMath::Abs(TMath::Tan(thetas1/2.)+1e-4));
	    etas2  = - TMath::Log(TMath::Abs(TMath::Tan(thetas2/2.)+1e-4));	    
	    deta   = etas1 - etas2;
	    pts1   = ps1 * TMath::Sin(thetas1);
	    pts2   = ps2 * TMath::Sin(thetas2);	    
	    
	    m      =  TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
	    Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
	    printf("Mass %f\n", m);
	    if (pts1 > 0. && pts2 > 0.)
	    massHS->Fill(m, wgt);
	}
	
//  Finish event
//      I/O
//	
	stack->FinishEvent();
	header->SetStack(stack);
	rl->TreeE()->Fill();
	header->SetNprimary(stack->GetNprimary());
	header->SetNtrack(stack->GetNtrack());  
    } // event loop
    TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
    massHU->Draw();
    massHS->Draw("same");

    TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
    etaH->Draw();
    
//
//                         Termination
//  Generator
    gener->FinishRun();
//  Stack
    stack->FinishRun();
//  Write file
    rl->WriteHeader("OVERWRITE");
    gener->Write();
    rl->Write();
}


AliGenerator*  CreateGenerator()
{
    AliGenParam *gener 
	= new AliGenParam(2000, AliGenMUONlib::kUpsilon, "");
        
    gener->SetMomentumRange(0,999);
    gener->SetPtRange(0,100.);
    gener->SetPhiRange(-180, 180);
    gener->SetYRange(2.5,4);
    gener->SetCutOnChild(1);
    gener->SetChildThetaRange(2.0,9.);
    gener->SetForceDecay(kDiMuon);

    AliDecayer *decayer = new AliDecayerPythia();
    decayer->SetForceDecay(kAll);
    decayer->Init();
    gener->SetDecayer(decayer);
    gener->Init();
    
    gener->Draw();
    
    return gener;
}
Commit	Line	Data
bc29c1f7	1	AliGenerator* CreateGenerator();
	2
	3	void fastGenQuarkonia(Int_t nev = 1, char* filename = "galice.root")
	4	{
	5	// Update data base
	6	AliPDG::AddParticlesToPdgDataBase();
	7
	8	// Create the fast tracker
	9	AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
	10	AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
	11	AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
	12	AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
	13	acc->SetBackground(0);
	14	eff->SetBackground(0);
	15	res->SetBackground(0);
	16	acc->Init();
	17	eff->Init();
	18	res->Init();
	19	AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
	20	tracker->AddResponse(trigeff);
	21	tracker->AddResponse(acc);
	22	tracker->AddResponse(eff);
	23	tracker->AddResponse(res);
	24	tracker->Init();
	25	tracker->Dump();
	26	// Histos
	27	TH1F* massHU = new TH1F("massHU", "Mass Spectrum: ", 500, 5, 15.);
	28	TH1F* massHS = new TH1F("massHS", "Mass Spectrum Smeared: ", 500, 5, 15.);
0cc695c9	29	TH2F* etaH = new TH2F("etaH", "eta vs etas ", 100, 0., 5., 100, 0., 5.);
0cc695c9	30	TH1F* etadH = new TH1F("etaH", "eta vs etas ", 100, -1., 1.);
bc29c1f7	31	//
	32	// Construction
	33	//
0cc695c9	34	AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
	35
	36	rl->SetCompressionLevel(2);
	37	rl->SetNumberOfEventsPerFile(nev);
	38	rl->LoadKinematics("RECREATE");
	39	rl->MakeTree("E");
	40	gAlice->SetRunLoader(rl);
	41
bc29c1f7	42	// Create stack
0cc695c9	43	rl->MakeStack();
	44	AliStack* stack = rl->Stack();
	45
	46	// Header
	47	AliHeader* header = rl->GetHeader();
bc29c1f7	48	//
	49	// Create and Initialize Generator
	50	AliGenerator *gener = CreateGenerator();
0cc695c9	51	gener->Init();
bc29c1f7	52	gener->SetStack(stack);
	53
	54	//
	55	// Event Loop
	56	//
	57	Int_t iev;
	58
	59	for (iev = 0; iev < nev; iev++) {
	60
	61	printf("\n \n Event number %d \n \n", iev);
	62
0cc695c9	63
bc29c1f7	64	// Initialize event
bc29c1f7	65	header->Reset(0,iev);
0cc695c9	66	rl->SetEventNumber(iev);
bc29c1f7	67	stack->Reset();
0cc695c9	68	rl->MakeTree("K");
bc29c1f7	69	// Generate event
	70	gener->Generate();
	71
	72	// Analysis
	73	Int_t npart = stack->GetNprimary();
bc29c1f7	74	for (Int_t part=0; part<npart; part++) {
	75	TParticle *MPart = stack->Particle(part);
	76	Int_t mpart = MPart->GetPdgCode();
	77	if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
	78	Int_t ch1 = MPart->GetFirstDaughter();
	79	Int_t ch2 = MPart->GetLastDaughter();
	80
	81	TParticle *muon1 = stack->Particle(ch1);
	82	TParticle *muon2 = stack->Particle(ch2);
	83
	84	Float_t theta1 = muon1->Theta();
	85	Float_t thetad1 = theta1 * 180./TMath::Pi();
	86	Float_t eta1 = muon1->Eta();
	87	Float_t pt1 = muon1->Pt();
	88	Float_t phi1 = muon1->Phi();
	89	Float_t phid1 = phi1 * 180./TMath::Pi() - 180.;
	90	Float_t p1 = muon1->P();
	91
	92	if (thetad1 > 9. \|\| thetad1 < 2.) continue;
	93
	94	Float_t theta2 = muon2->Theta();
	95	Float_t thetad2 = theta2 * 180./TMath::Pi();
	96	Float_t eta2 = muon2->Eta();
	97	Float_t pt2 = muon2->Pt();
	98	Float_t phi2 = muon2->Phi();
	99	Float_t phid2 = phi2 * 180./TMath::Pi() - 180.;
	100	Float_t p2 = muon2->P();
	101
	102
	103	if (thetad2 > 9. \|\| thetad2 < 2.) continue;
	104
	105	Float_t dphi = phi1 - phi2;
	106	Float_t deta = eta1 - eta2;
	107	Float_t m = TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
	108
bc29c1f7	109	massHU->Fill(m);
	110	// Smear
	111	// the mu+
	112	Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
	113	Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH;
	114
	115	res->SetCharge(1);
	116	eff->SetCharge(1);
	117	acc->SetCharge(1);
	118	res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
0cc695c9	119
0cc695c9	120
bc29c1f7	121	Float_t effp = eff->Evaluate(pt1, thetad1, phid1);
	122	Float_t accp = acc->Evaluate(pt1, thetad1, phid1);
	123	trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH);
	124	thetas1 *= TMath::Pi()/180.;
	125	phis1 *= TMath::Pi()/180.;
0cc695c9	126	Float_t etas;
	127
	128	if (TMath::Tan(thetas1/2.) > 0.) {
	129	etas = -TMath::Log(TMath::Tan(thetas1/2.));
	130	etaH->Fill(etas, eta1);
	131	etadH->Fill(etas-eta1);
	132	}
	133
bc29c1f7	134	// the mu-
	135	res->SetCharge(-1);
	136	acc->SetCharge(-1);
	137	eff->SetCharge(-1);
	138	res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
	139	Float_t effn = eff->Evaluate(pt2, thetad2, phid2);
	140	Float_t accn = acc->Evaluate(pt2, thetad2, phid2);
	141	trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH);
	142	thetas2 *= TMath::Pi()/180.;
	143	phis2 *= TMath::Pi()/180.;
0cc695c9	144	if (TMath::Tan(thetas2/2.) > 0) {
	145	etas = -TMath::Log(TMath::Tan(thetas2/2.));
	146	etaH->Fill(etas, eta2);
	147	etadH->Fill(etas-eta2);
	148	}
	149
bc29c1f7	150	dphi = phis1 - phis2;
0cc695c9	151	etas1 = - TMath::Log(TMath::Abs(TMath::Tan(thetas1/2.)+1e-4));
0cc695c9	152	etas2 = - TMath::Log(TMath::Abs(TMath::Tan(thetas2/2.)+1e-4));
bc29c1f7	153	deta = etas1 - etas2;
	154	pts1 = ps1 * TMath::Sin(thetas1);
	155	pts2 = ps2 * TMath::Sin(thetas2);
	156
	157	m = TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
0cc695c9	158	Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
bc29c1f7	159	printf("Mass %f\n", m);
	160	if (pts1 > 0. && pts2 > 0.)
	161	massHS->Fill(m, wgt);
	162	}
	163
	164	// Finish event
bc29c1f7	165	// I/O
bc29c1f7	166	//
0cc695c9	167	stack->FinishEvent();
	168	header->SetStack(stack);
	169	rl->TreeE()->Fill();
	170	header->SetNprimary(stack->GetNprimary());
	171	header->SetNtrack(stack->GetNtrack());
bc29c1f7	172	} // event loop
	173	TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
	174	massHU->Draw();
	175	massHS->Draw("same");
0cc695c9	176
	177	TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
	178	etaH->Draw();
bc29c1f7	179
	180	//
	181	// Termination
	182	// Generator
	183	gener->FinishRun();
bc29c1f7	184	// Stack
	185	stack->FinishRun();
	186	// Write file
0cc695c9	187	rl->WriteHeader("OVERWRITE");
bc29c1f7	188	gener->Write();
0cc695c9	189	rl->Write();
bc29c1f7	190	}
	191
	192
	193	AliGenerator* CreateGenerator()
	194	{
	195	AliGenParam *gener
0cc695c9	196	= new AliGenParam(2000, AliGenMUONlib::kUpsilon, "");
bc29c1f7	197
	198	gener->SetMomentumRange(0,999);
	199	gener->SetPtRange(0,100.);
	200	gener->SetPhiRange(-180, 180);
0cc695c9	201	gener->SetYRange(2.5,4);
0cc695c9	202	gener->SetCutOnChild(1);
bc29c1f7	203	gener->SetChildThetaRange(2.0,9.);
	204	gener->SetForceDecay(kDiMuon);
	205
	206	AliDecayer *decayer = new AliDecayerPythia();
	207	decayer->SetForceDecay(kAll);
	208	decayer->Init();
	209	gener->SetDecayer(decayer);
	210	gener->Init();
	211
	212	gener->Draw();
	213
	214	return gener;
	215	}
	216
	217
	218
	219
	220
	221
	222
	223
	224
	225
	226
	227