nvartrk=7 in ESE task
[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.);
0cc695c9 29 TH2F* etaH = new TH2F("etaH", "eta vs etas ", 100, 0., 5., 100, 0., 5.);
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
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;
bc77d57e 113 Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH, trigeffnA;
bc29c1f7 114
115 res->SetCharge(1);
116 eff->SetCharge(1);
117 acc->SetCharge(1);
118 res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
0cc695c9 119
120
bc29c1f7 121 Float_t effp = eff->Evaluate(pt1, thetad1, phid1);
122 Float_t accp = acc->Evaluate(pt1, thetad1, phid1);
bc77d57e 123 trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH, trigeffnA);
bc29c1f7 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);
bc77d57e 141 trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH, trigeffnA);
bc29c1f7 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));
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
166//
0cc695c9 167 stack->FinishEvent();
168 header->SetStack(stack);
169 rl->TreeE()->Fill();
170 header->SetNprimary(stack->GetNprimary());
171 header->SetNtrack(stack->GetNtrack());
bc77d57e 172 rl->WriteKinematics("OVERWRITE");
bc29c1f7 173 } // event loop
174 TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
175 massHU->Draw();
176 massHS->Draw("same");
0cc695c9 177
178 TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
179 etaH->Draw();
bc29c1f7 180
181//
182// Termination
183// Generator
184 gener->FinishRun();
bc29c1f7 185// Stack
186 stack->FinishRun();
187// Write file
0cc695c9 188 rl->WriteHeader("OVERWRITE");
bc29c1f7 189 gener->Write();
0cc695c9 190 rl->Write();
bc29c1f7 191}
192
193
194AliGenerator* CreateGenerator()
195{
196 AliGenParam *gener
0cc695c9 197 = new AliGenParam(2000, AliGenMUONlib::kUpsilon, "");
bc29c1f7 198
199 gener->SetMomentumRange(0,999);
200 gener->SetPtRange(0,100.);
55bc562c 201 gener->SetPhiRange(0., 360.);
0cc695c9 202 gener->SetYRange(2.5,4);
203 gener->SetCutOnChild(1);
bc29c1f7 204 gener->SetChildThetaRange(2.0,9.);
205 gener->SetForceDecay(kDiMuon);
206
207 AliDecayer *decayer = new AliDecayerPythia();
208 decayer->SetForceDecay(kAll);
209 decayer->Init();
210 gener->SetDecayer(decayer);
211 gener->Init();
212
213 gener->Draw();
214
215 return gener;
216}
217
218
219
220
221
222
223
224
225
226
227
228