From bc29c1f7bada29955c721650c21841831004fc22 Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 12 Jun 2003 07:48:53 +0000 Subject: [PATCH] First commit. --- FASTSIM/fastGenQuarkonia.C | 212 +++++++++++++++++++++++++++++++++++++ 1 file changed, 212 insertions(+) create mode 100644 FASTSIM/fastGenQuarkonia.C diff --git a/FASTSIM/fastGenQuarkonia.C b/FASTSIM/fastGenQuarkonia.C new file mode 100644 index 00000000000..52e7af30bf4 --- /dev/null +++ b/FASTSIM/fastGenQuarkonia.C @@ -0,0 +1,212 @@ +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.); +// +// Construction +// +// Output file + TFile* file = new TFile(filename, "recreate"); +// Create stack + AliStack* stack = new AliStack(10000); + stack->MakeTree(0, filename); + +// Create Header + AliHeader* header = new AliHeader(); +// Create Header Tree + TTree* treeE = new TTree("TE","Headers"); + treeE->Branch("Header", "AliHeader", &header, 4000, 0); + treeE->Write(); +// +// Create and Initialize Generator + AliGenerator *gener = CreateGenerator(); + AliPythia* pyth = AliPythia::Instance(); + 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); + stack->Reset(); + stack->BeginEvent(iev); + +// Generate event + gener->Generate(); + +// Analysis + Int_t npart = stack->GetNprimary(); +// printf("Analyse %d Particles\n", npart); + for (Int_t part=0; partParticle(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))); + + printf("Mass %f %f %f\n", m, phid1, phid2); + 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.; + + // 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.; + + dphi = phis1 - phis2; + etas1 = - TMath::Log(TMath::Tan(thetas1/2.)); + etas2 = - TMath::Log(TMath::Tan(thetas2/2.)); + 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; + Float_t wgt = accn * accp * effp * effn ; + printf("Mass %f\n", m); + if (pts1 > 0. && pts2 > 0.) + massHS->Fill(m, wgt); + } + +// Finish event + header->SetNprimary(stack->GetNprimary()); + header->SetNtrack(stack->GetNtrack()); +// I/O +// +// stack->FinishEvent(); +// header->SetStack(stack); +// treeE->Fill(); +// (stack->TreeK())->Write(0,TObject::kOverwrite); + } // event loop + TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700); + massHU->Draw(); + massHS->Draw("same"); + +// +// Termination +// Generator + gener->FinishRun(); +// Header + treeE->Write(0,TObject::kOverwrite); + delete treeE; treeE = 0; +// Stack + stack->FinishRun(); +// Write file + gener->Write(); + file->Write(); +} + + +AliGenerator* CreateGenerator() +{ + AliGenParam *gener + = new AliGenParam(200000, AliGenMUONlib::kUpsilon, ""); + + gener->SetMomentumRange(0,999); + gener->SetPtRange(0,100.); + gener->SetPhiRange(-180, 180); +// gener->SetYRange(2.5,4); + gener->SetYRange(-6. , 6.); + gener->SetCutOnChild(0); + 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; +} + + + + + + + + + + + + -- 2.39.3