not needed as no one start the session from RICH
[u/mrichter/AliRoot.git] / FASTSIM / fastGenQuarkonia.C
CommitLineData
bc29c1f7 1AliGenerator* CreateGenerator();
2
3void 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.);
29//
30// Construction
31//
32// Output file
33 TFile* file = new TFile(filename, "recreate");
34// Create stack
35 AliStack* stack = new AliStack(10000);
36 stack->MakeTree(0, filename);
37
38// Create Header
39 AliHeader* header = new AliHeader();
40// Create Header Tree
41 TTree* treeE = new TTree("TE","Headers");
42 treeE->Branch("Header", "AliHeader", &header, 4000, 0);
43 treeE->Write();
44//
45// Create and Initialize Generator
46 AliGenerator *gener = CreateGenerator();
47 AliPythia* pyth = AliPythia::Instance();
48 gener->SetStack(stack);
49
50//
51// Event Loop
52//
53 Int_t iev;
54
55 for (iev = 0; iev < nev; iev++) {
56
57 printf("\n \n Event number %d \n \n", iev);
58
59// Initialize event
60 header->Reset(0,iev);
61 stack->Reset();
62 stack->BeginEvent(iev);
63
64// Generate event
65 gener->Generate();
66
67// Analysis
68 Int_t npart = stack->GetNprimary();
69// printf("Analyse %d Particles\n", npart);
70 for (Int_t part=0; part<npart; part++) {
71 TParticle *MPart = stack->Particle(part);
72 Int_t mpart = MPart->GetPdgCode();
73 if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
74 Int_t ch1 = MPart->GetFirstDaughter();
75 Int_t ch2 = MPart->GetLastDaughter();
76
77 TParticle *muon1 = stack->Particle(ch1);
78 TParticle *muon2 = stack->Particle(ch2);
79
80 Float_t theta1 = muon1->Theta();
81 Float_t thetad1 = theta1 * 180./TMath::Pi();
82 Float_t eta1 = muon1->Eta();
83 Float_t pt1 = muon1->Pt();
84 Float_t phi1 = muon1->Phi();
85 Float_t phid1 = phi1 * 180./TMath::Pi() - 180.;
86 Float_t p1 = muon1->P();
87
88 if (thetad1 > 9. || thetad1 < 2.) continue;
89
90 Float_t theta2 = muon2->Theta();
91 Float_t thetad2 = theta2 * 180./TMath::Pi();
92 Float_t eta2 = muon2->Eta();
93 Float_t pt2 = muon2->Pt();
94 Float_t phi2 = muon2->Phi();
95 Float_t phid2 = phi2 * 180./TMath::Pi() - 180.;
96 Float_t p2 = muon2->P();
97
98
99 if (thetad2 > 9. || thetad2 < 2.) continue;
100
101 Float_t dphi = phi1 - phi2;
102 Float_t deta = eta1 - eta2;
103 Float_t m = TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
104
105 printf("Mass %f %f %f\n", m, phid1, phid2);
106 massHU->Fill(m);
107// Smear
108 // the mu+
109 Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
110 Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH;
111
112 res->SetCharge(1);
113 eff->SetCharge(1);
114 acc->SetCharge(1);
115 res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
116 Float_t effp = eff->Evaluate(pt1, thetad1, phid1);
117 Float_t accp = acc->Evaluate(pt1, thetad1, phid1);
118 trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH);
119 thetas1 *= TMath::Pi()/180.;
120 phis1 *= TMath::Pi()/180.;
121
122 // the mu-
123 res->SetCharge(-1);
124 acc->SetCharge(-1);
125 eff->SetCharge(-1);
126 res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
127 Float_t effn = eff->Evaluate(pt2, thetad2, phid2);
128 Float_t accn = acc->Evaluate(pt2, thetad2, phid2);
129 trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH);
130 thetas2 *= TMath::Pi()/180.;
131 phis2 *= TMath::Pi()/180.;
132
133 dphi = phis1 - phis2;
134 etas1 = - TMath::Log(TMath::Tan(thetas1/2.));
135 etas2 = - TMath::Log(TMath::Tan(thetas2/2.));
136 deta = etas1 - etas2;
137 pts1 = ps1 * TMath::Sin(thetas1);
138 pts2 = ps2 * TMath::Sin(thetas2);
139
140 m = TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
141// Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
142 Float_t wgt = accn * accp * effp * effn ;
143 printf("Mass %f\n", m);
144 if (pts1 > 0. && pts2 > 0.)
145 massHS->Fill(m, wgt);
146 }
147
148// Finish event
149 header->SetNprimary(stack->GetNprimary());
150 header->SetNtrack(stack->GetNtrack());
151// I/O
152//
153// stack->FinishEvent();
154// header->SetStack(stack);
155// treeE->Fill();
156// (stack->TreeK())->Write(0,TObject::kOverwrite);
157 } // event loop
158 TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
159 massHU->Draw();
160 massHS->Draw("same");
161
162//
163// Termination
164// Generator
165 gener->FinishRun();
166// Header
167 treeE->Write(0,TObject::kOverwrite);
168 delete treeE; treeE = 0;
169// Stack
170 stack->FinishRun();
171// Write file
172 gener->Write();
173 file->Write();
174}
175
176
177AliGenerator* CreateGenerator()
178{
179 AliGenParam *gener
180 = new AliGenParam(200000, AliGenMUONlib::kUpsilon, "");
181
182 gener->SetMomentumRange(0,999);
183 gener->SetPtRange(0,100.);
184 gener->SetPhiRange(-180, 180);
185// gener->SetYRange(2.5,4);
186 gener->SetYRange(-6. , 6.);
187 gener->SetCutOnChild(0);
188 gener->SetChildThetaRange(2.0,9.);
189 gener->SetForceDecay(kDiMuon);
190
191 AliDecayer *decayer = new AliDecayerPythia();
192 decayer->SetForceDecay(kAll);
193 decayer->Init();
194 gener->SetDecayer(decayer);
195 gener->Init();
196
197 gener->Draw();
198
199 return gener;
200}
201
202
203
204
205
206
207
208
209
210
211
212