]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenReaderEMD.cxx
move from AliAODTrack to AliVTrack
[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17// Class to read events from external (TNtupla) file
18// Events -> neutron removal by EM dissociation of Pb nuclei
19// Data from RELDIS code (by I. Pshenichov)
20
21#include <TFile.h>
22#include <TParticle.h>
23#include <TTree.h>
24#include <TVirtualMC.h>
25#include <TDatabasePDG.h>
26#include <TPDGCode.h>
27#include "AliGenReaderEMD.h"
28#include "AliStack.h"
29
30
31ClassImp(AliGenReaderEMD)
32
33AliGenReaderEMD::AliGenReaderEMD():
34 fStartEvent(0),
35 fNcurrent(0),
36 fNparticle(0),
37 fTreeNtuple(0),
38 fPcToTrack(0),
39 fOffset(0),
40 fNnAside(0),
41 fEnAside(0),
42 fnPDGCode(0),
43 fNnCside(0),
44 fEnCside(0),
45 fNpAside(0),
46 fEtapAside(0),
47 fpPDGCode(0),
48 fNpCside(0),
49 fEtapCside(0),
50 fNppAside(0),
51 fEtappAside(0),
52 fppPDGCode(0),
53 fNppCside(0),
54 fEtappCside(0),
55 fNpmAside(0),
56 fEtapmAside(0),
57 fpmPDGCode(0),
58 fNpmCside(0),
59 fEtapmCside(0),
60 fNp0Aside(0),
61 fEtap0Aside(0),
62 fp0PDGCode(0),
63 fNp0Cside(0),
64 fEtap0Cside(0),
65 fNetaAside(0),
66 fEtaetaAside(0),
67 fetaPDGCode(0),
68 fNetaCside(0),
69 fEtaetaCside(0),
70 fNomegaAside(0),
71 fEtaomegaAside(0),
72 fomegaPDGCode(0),
73 fNomegaCside(0),
74 fEtaomegaCside(0)
75{
76// Std constructor
77 for(int i=0; i<70; i++){
78 fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
79 fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
80 if(i<50){
81 fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
82 fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
83 if(i<30){
84 fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
85 fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
86 fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
87 fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
88 fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
89 fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
90 if(i<15){
91 fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
92 fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
93 fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
94 fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
95 }
96 }
97 }
98 }
99 if(fPcToTrack==kAll) printf("\n\t *** AliGenReaderEMD will track all produced particles \n\n");
100 else if(fPcToTrack==kNotNucleons) printf("\n\t *** AliGenReaderEMD will track all produced particles except nucleons\n\n");
101 else if(fPcToTrack==kOnlyNucleons) printf("\n\t *** AliGenReaderEMD will track only nucleons\n\n");
102}
103
104
105AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
106 AliGenReader(reader),
107 fStartEvent(0),
108 fNcurrent(0),
109 fNparticle(0),
110 fTreeNtuple(0),
111 fPcToTrack(0),
112 fOffset(0),
113 fNnAside(0),
114 fEnAside(0),
115 fnPDGCode(0),
116 fNnCside(0),
117 fEnCside(0),
118 fNpAside(0),
119 fEtapAside(0),
120 fpPDGCode(0),
121 fNpCside(0),
122 fEtapCside(0),
123 fNppAside(0),
124 fEtappAside(0),
125 fppPDGCode(0),
126 fNppCside(0),
127 fEtappCside(0),
128 fNpmAside(0),
129 fEtapmAside(0),
130 fpmPDGCode(0),
131 fNpmCside(0),
132 fEtapmCside(0),
133 fNp0Aside(0),
134 fEtap0Aside(0),
135 fp0PDGCode(0),
136 fNp0Cside(0),
137 fEtap0Cside(0),
138 fNetaAside(0),
139 fEtaetaAside(0),
140 fetaPDGCode(0),
141 fNetaCside(0),
142 fEtaetaCside(0),
143 fNomegaAside(0),
144 fEtaomegaAside(0),
145 fomegaPDGCode(0),
146 fNomegaCside(0),
147 fEtaomegaCside(0)
148{
149 // Copy Constructor
150 reader.Copy(*this);
151}
152 // -----------------------------------------------------------------------------------
153AliGenReaderEMD::~AliGenReaderEMD()
154{
155 delete fTreeNtuple;
156}
157
158// -----------------------------------------------------------------------------------
159AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
160{
161// Assignment operator
162 rhs.Copy(*this);
163 return *this;
164}
165
166// -----------------------------------------------------------------------------------
167void AliGenReaderEMD::Copy(TObject&) const
168{
169 //
170 // Copy
171 //
172 Fatal("Copy","Not implemented!\n");
173}
174
175// -----------------------------------------------------------------------------------
176void AliGenReaderEMD::Init()
177{
178//
179// Reset the existing file environment and open a new root file
180
181 TFile *pFile=0;
182 if (!pFile) {
183 pFile = new TFile(fFileName);
184 pFile->cd();
185 printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
186 }
187 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
188 fNcurrent = fStartEvent;
189
190 TTree *Ntu=fTreeNtuple;
191 //
192 // Set branch addresses
193 // **** neutrons
194 Ntu->SetBranchAddress("Nleft",&fNnAside);
195 Ntu->SetBranchAddress("Eleft",&fEnAside);
196 Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
197 Ntu->SetBranchAddress("Pxl", fPxnAside);
198 Ntu->SetBranchAddress("Pyl", fPynAside);
199 Ntu->SetBranchAddress("Pzl", fPznAside);
200 Ntu->SetBranchAddress("Nright",&fNnCside);
201 Ntu->SetBranchAddress("Eright",&fEnCside);
202 Ntu->SetBranchAddress("Pxr", fPxnCside);
203 Ntu->SetBranchAddress("Pyr", fPynCside);
204 Ntu->SetBranchAddress("Pzr", fPznCside);
205 // **** protons
206 Ntu->SetBranchAddress("Nleft_p",&fNpAside);
207 Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
208 Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
209 Ntu->SetBranchAddress("Pxl_p", fPxpAside);
210 Ntu->SetBranchAddress("Pyl_p", fPypAside);
211 Ntu->SetBranchAddress("Pzl_p", fPzpAside);
212 Ntu->SetBranchAddress("Nright_p",&fNpCside);
213 Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
214 Ntu->SetBranchAddress("Pxr_p", fPxpCside);
215 Ntu->SetBranchAddress("Pyr_p", fPypCside);
216 Ntu->SetBranchAddress("Pzr_p", fPzpCside);
217 // **** pi+
218 Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
219 Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
220 Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
221 Ntu->SetBranchAddress("Pxl_pp", fPxppAside);
222 Ntu->SetBranchAddress("Pyl_pp", fPyppAside);
223 Ntu->SetBranchAddress("Pzl_pp", fPzppAside);
224 Ntu->SetBranchAddress("Nright_pp",&fNppCside);
225 Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
226 Ntu->SetBranchAddress("Pxr_pp", fPxppCside);
227 Ntu->SetBranchAddress("Pyr_pp", fPyppCside);
228 Ntu->SetBranchAddress("Pzr_pp", fPzppCside);
229 // **** pi-
230 Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
231 Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
232 Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
233 Ntu->SetBranchAddress("Pxl_pm", fPxpmAside);
234 Ntu->SetBranchAddress("Pyl_pm", fPypmAside);
235 Ntu->SetBranchAddress("Pzl_pm", fPzpmAside);
236 Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
237 Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
238 Ntu->SetBranchAddress("Pxr_pm", fPxpmCside);
239 Ntu->SetBranchAddress("Pyr_pm", fPypmCside);
240 Ntu->SetBranchAddress("Pzr_pm", fPzpmCside);
241 // **** pi0
242 Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
243 Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
244 Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
245 Ntu->SetBranchAddress("Pxl_p0", fPxp0Aside);
246 Ntu->SetBranchAddress("Pyl_p0", fPyp0Aside);
247 Ntu->SetBranchAddress("Pzl_p0", fPzp0Aside);
248 Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
249 Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
250 Ntu->SetBranchAddress("Pxr_p0", fPxp0Cside);
251 Ntu->SetBranchAddress("Pyr_p0", fPyp0Cside);
252 Ntu->SetBranchAddress("Pzr_p0", fPzp0Cside);
253 // **** eta
254 Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
255 Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
256 Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
257 Ntu->SetBranchAddress("Pxl_et", fPxetaAside);
258 Ntu->SetBranchAddress("Pyl_et", fPyetaAside);
259 Ntu->SetBranchAddress("Pzl_et", fPzetaAside);
260 Ntu->SetBranchAddress("Nright_et",&fNetaCside);
261 Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
262 Ntu->SetBranchAddress("Pxr_et", fPxetaCside);
263 Ntu->SetBranchAddress("Pyr_et", fPyetaCside);
264 Ntu->SetBranchAddress("Pzr_et", fPzetaCside);
265 // **** omega
266 Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
267 Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
268 Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
269 Ntu->SetBranchAddress("Pxl_om", fPxomegaAside);
270 Ntu->SetBranchAddress("Pyl_om", fPyomegaAside);
271 Ntu->SetBranchAddress("Pzl_om", fPzomegaAside);
272 Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
273 Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
274 Ntu->SetBranchAddress("Pxr_om", fPxomegaCside);
275 Ntu->SetBranchAddress("Pyr_om", fPyomegaCside);
276 Ntu->SetBranchAddress("Pzr_om", fPzomegaCside);
277}
278
279// -----------------------------------------------------------------------------------
280Int_t AliGenReaderEMD::NextEvent()
281{
282 // Read the next event
283 Int_t nTracks=0;
284 fNparticle = 0; fOffset=0;
285
286 TFile* pFile = fTreeNtuple->GetCurrentFile();
287 pFile->cd();
288
289
290 Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
291 if(fNcurrent < nentries) {
292 fTreeNtuple->GetEvent(fNcurrent);
293 if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
294 //
295 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
296 nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
297 }
298 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
299 nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
300 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
301 }
302 fNcurrent++;
303 printf("\t #### Putting %d particles in the stack\n", nTracks);
304 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t %d+%d neutrons, %d+%d protons\n",
305 fNnAside,fNnCside, fNpAside,fNpCside);
306 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
307 fNppAside,fNppCside,fNpmAside,fNpmCside,
308 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);
309 return nTracks;
310 }
311
312 return 0;
313}
314
315// -----------------------------------------------------------------------------------
316TParticle* AliGenReaderEMD::NextParticle()
317{
318 // Read the next particle
319 Float_t p[4]={0.,0.,0.,0.};
320 int pdgCode=0;
321
322 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
323 if(fNparticle<fNnAside){
324 p[0] = fPxnAside[fNparticle];
325 p[1] = fPynAside[fNparticle];
326 p[2] = fPznAside[fNparticle];
327 pdgCode = fnPDGCode;
328// printf(" pc%d n sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
329 }
330 else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
331 p[0] = fPxnCside[fNparticle];
332 p[1] = fPynCside[fNparticle];
333 p[2] = fPznCside[fNparticle];
334 pdgCode = fnPDGCode;
335// printf(" pc%d n sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
336 }
337 else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
338 p[0] = fPxpAside[fNparticle];
339 p[1] = fPypAside[fNparticle];
340 p[2] = fPzpAside[fNparticle];
341 pdgCode = fpPDGCode;
342// printf(" pc%d p sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
343 }
344 else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
345 p[0] = fPxpCside[fNparticle];
346 p[1] = fPypCside[fNparticle];
347 p[2] = fPzpCside[fNparticle];
348 pdgCode = fpPDGCode;
349// printf(" pc%d p sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
350 }
351 fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
352 } //**********************************************************************************************
353 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
354 if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
355 p[0] = fPxppAside[fNparticle];
356 p[1] = fPyppAside[fNparticle];
357 p[2] = fPzppAside[fNparticle];
358 pdgCode = fppPDGCode;
359// printf(" pc%d pi+ sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
360 }
361 if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
362 p[0] = fPxppCside[fNparticle];
363 p[1] = fPyppCside[fNparticle];
364 p[2] = fPzppCside[fNparticle];
365 pdgCode = fppPDGCode;
366// printf(" pc%d pi+ sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
367 }
368 if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
369 p[0] = fPxpmAside[fNparticle];
370 p[1] = fPypmAside[fNparticle];
371 p[2] = fPzpmAside[fNparticle];
372 pdgCode = fpmPDGCode;
373// printf(" pc%d pi- sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
374 }
375 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
376 p[0] = fPxpmCside[fNparticle];
377 p[1] = fPypmCside[fNparticle];
378 p[2] = fPzpmCside[fNparticle];
379 pdgCode = fpmPDGCode;
380// printf(" pc%d pi- sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
381 }
382 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
383 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
384 p[0] = fPxp0Aside[fNparticle];
385 p[1] = fPyp0Aside[fNparticle];
386 p[2] = fPzp0Aside[fNparticle];
387 pdgCode = fp0PDGCode;
388// printf(" pc%d pi0 sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
389 }
390 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside &&
391 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
392 p[0] = fPxp0Cside[fNparticle];
393 p[1] = fPyp0Cside[fNparticle];
394 p[2] = fPzp0Cside[fNparticle];
395 pdgCode = fp0PDGCode;
396// printf(" pc%d pi0 sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
397 }
398 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
399 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
400 p[0] = fPxetaAside[fNparticle];
401 p[1] = fPyetaAside[fNparticle];
402 p[2] = fPzetaAside[fNparticle];
403 pdgCode = fetaPDGCode;
404// printf(" pc%d eta sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
405 }
406 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside &&
407 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
408 p[0] = fPxetaCside[fNparticle];
409 p[1] = fPyetaCside[fNparticle];
410 p[2] = fPzetaCside[fNparticle];
411 pdgCode = fetaPDGCode;
412// printf(" pc%d eta sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
413 }
414 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
415 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
416 p[0] = fPxomegaAside[fNparticle];
417 p[1] = fPyomegaAside[fNparticle];
418 p[2] = fPzomegaAside[fNparticle];
419 pdgCode = fomegaPDGCode;
420// printf(" pc%d omega sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
421 }
422 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
423 && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
424 p[0] = fPxomegaCside[fNparticle];
425 p[1] = fPyomegaCside[fNparticle];
426 p[2] = fPzomegaCside[fNparticle];
427 pdgCode = fomegaPDGCode;
428// printf(" pc%d omega sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
429 }
430
431 }
432
433 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
434 Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
435 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
436
437 if(p[3]<=amass)
438 Warning("Generate","Particle %d E = %f GeV mass = %f GeV/c^2",pdgCode,p[3],amass);
439
440 TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
441 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
442 particle->SetBit(kTransportBit);
443 fNparticle++;
444 return particle;
445}