]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenReaderEMD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Protection in copy constructor
[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 for(int i=0; i<70; i++){
151 fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
152 fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
153 if(i<50){
154 fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
155 fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
156 if(i<30){
157 fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
158 fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
159 fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
160 fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
161 fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
162 fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
163 if(i<15){
164 fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
165 fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
166 fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
167 fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
168 }
169 }
170 }
171 }
172 reader.Copy(*this);
173}
174 // -----------------------------------------------------------------------------------
175AliGenReaderEMD::~AliGenReaderEMD()
176{
177 delete fTreeNtuple;
178}
179
180// -----------------------------------------------------------------------------------
181AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
182{
183// Assignment operator
184 rhs.Copy(*this);
185 return *this;
186}
187
188// -----------------------------------------------------------------------------------
189void AliGenReaderEMD::Copy(TObject&) const
190{
191 //
192 // Copy
193 //
194 Fatal("Copy","Not implemented!\n");
195}
196
197// -----------------------------------------------------------------------------------
198void AliGenReaderEMD::Init()
199{
200//
201// Reset the existing file environment and open a new root file
202
203 TFile *pFile=0;
204 if (!pFile) {
205 pFile = new TFile(fFileName);
206 pFile->cd();
207 printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
208 }
209 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
210 fNcurrent = fStartEvent;
211
212 TTree *Ntu=fTreeNtuple;
213 //
214 // Set branch addresses
215 // **** neutrons
216 Ntu->SetBranchAddress("Nleft",&fNnAside);
217 Ntu->SetBranchAddress("Eleft",&fEnAside);
218 Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
219 Ntu->SetBranchAddress("Pxl", fPxnAside);
220 Ntu->SetBranchAddress("Pyl", fPynAside);
221 Ntu->SetBranchAddress("Pzl", fPznAside);
222 Ntu->SetBranchAddress("Nright",&fNnCside);
223 Ntu->SetBranchAddress("Eright",&fEnCside);
224 Ntu->SetBranchAddress("Pxr", fPxnCside);
225 Ntu->SetBranchAddress("Pyr", fPynCside);
226 Ntu->SetBranchAddress("Pzr", fPznCside);
227 // **** protons
228 Ntu->SetBranchAddress("Nleft_p",&fNpAside);
229 Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
230 Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
231 Ntu->SetBranchAddress("Pxl_p", fPxpAside);
232 Ntu->SetBranchAddress("Pyl_p", fPypAside);
233 Ntu->SetBranchAddress("Pzl_p", fPzpAside);
234 Ntu->SetBranchAddress("Nright_p",&fNpCside);
235 Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
236 Ntu->SetBranchAddress("Pxr_p", fPxpCside);
237 Ntu->SetBranchAddress("Pyr_p", fPypCside);
238 Ntu->SetBranchAddress("Pzr_p", fPzpCside);
239 // **** pi+
240 Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
241 Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
242 Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
243 Ntu->SetBranchAddress("Pxl_pp", fPxppAside);
244 Ntu->SetBranchAddress("Pyl_pp", fPyppAside);
245 Ntu->SetBranchAddress("Pzl_pp", fPzppAside);
246 Ntu->SetBranchAddress("Nright_pp",&fNppCside);
247 Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
248 Ntu->SetBranchAddress("Pxr_pp", fPxppCside);
249 Ntu->SetBranchAddress("Pyr_pp", fPyppCside);
250 Ntu->SetBranchAddress("Pzr_pp", fPzppCside);
251 // **** pi-
252 Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
253 Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
254 Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
255 Ntu->SetBranchAddress("Pxl_pm", fPxpmAside);
256 Ntu->SetBranchAddress("Pyl_pm", fPypmAside);
257 Ntu->SetBranchAddress("Pzl_pm", fPzpmAside);
258 Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
259 Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
260 Ntu->SetBranchAddress("Pxr_pm", fPxpmCside);
261 Ntu->SetBranchAddress("Pyr_pm", fPypmCside);
262 Ntu->SetBranchAddress("Pzr_pm", fPzpmCside);
263 // **** pi0
264 Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
265 Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
266 Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
267 Ntu->SetBranchAddress("Pxl_p0", fPxp0Aside);
268 Ntu->SetBranchAddress("Pyl_p0", fPyp0Aside);
269 Ntu->SetBranchAddress("Pzl_p0", fPzp0Aside);
270 Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
271 Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
272 Ntu->SetBranchAddress("Pxr_p0", fPxp0Cside);
273 Ntu->SetBranchAddress("Pyr_p0", fPyp0Cside);
274 Ntu->SetBranchAddress("Pzr_p0", fPzp0Cside);
275 // **** eta
276 Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
277 Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
278 Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
279 Ntu->SetBranchAddress("Pxl_et", fPxetaAside);
280 Ntu->SetBranchAddress("Pyl_et", fPyetaAside);
281 Ntu->SetBranchAddress("Pzl_et", fPzetaAside);
282 Ntu->SetBranchAddress("Nright_et",&fNetaCside);
283 Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
284 Ntu->SetBranchAddress("Pxr_et", fPxetaCside);
285 Ntu->SetBranchAddress("Pyr_et", fPyetaCside);
286 Ntu->SetBranchAddress("Pzr_et", fPzetaCside);
287 // **** omega
288 Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
289 Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
290 Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
291 Ntu->SetBranchAddress("Pxl_om", fPxomegaAside);
292 Ntu->SetBranchAddress("Pyl_om", fPyomegaAside);
293 Ntu->SetBranchAddress("Pzl_om", fPzomegaAside);
294 Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
295 Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
296 Ntu->SetBranchAddress("Pxr_om", fPxomegaCside);
297 Ntu->SetBranchAddress("Pyr_om", fPyomegaCside);
298 Ntu->SetBranchAddress("Pzr_om", fPzomegaCside);
299}
300
301// -----------------------------------------------------------------------------------
302Int_t AliGenReaderEMD::NextEvent()
303{
304 // Read the next event
305 Int_t nTracks=0;
306 fNparticle = 0; fOffset=0;
307
308 TFile* pFile = fTreeNtuple->GetCurrentFile();
309 pFile->cd();
310
311
312 Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
313 if(fNcurrent < nentries) {
314 fTreeNtuple->GetEvent(fNcurrent);
315 if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
316 //
317 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
318 nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
319 }
320 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
321 nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
322 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
323 }
324 fNcurrent++;
325 printf("\t #### Putting %d particles in the stack\n", nTracks);
326 /*if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t %d+%d neutrons, %d+%d protons\n",
327 fNnAside,fNnCside, fNpAside,fNpCside);
328 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
329 fNppAside,fNppCside,fNpmAside,fNpmCside,
330 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);*/
331 return nTracks;
332 }
333
334 return 0;
335}
336
337// -----------------------------------------------------------------------------------
338TParticle* AliGenReaderEMD::NextParticle()
339{
340 // Read the next particle
341 Float_t p[4]={0.,0.,0.,0.};
342 int pdgCode=0;
343
344 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
345 if(fNparticle<fNnAside){
346 p[0] = fPxnAside[fNparticle];
347 p[1] = fPynAside[fNparticle];
348 p[2] = fPznAside[fNparticle];
349 pdgCode = fnPDGCode;
350// printf(" pc%d n sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
351 }
352 else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
353 p[0] = fPxnCside[fNparticle];
354 p[1] = fPynCside[fNparticle];
355 p[2] = fPznCside[fNparticle];
356 pdgCode = fnPDGCode;
357// printf(" pc%d n sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
358 }
359 else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
360 p[0] = fPxpAside[fNparticle];
361 p[1] = fPypAside[fNparticle];
362 p[2] = fPzpAside[fNparticle];
363 pdgCode = fpPDGCode;
364// printf(" pc%d p sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
365 }
366 else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
367 p[0] = fPxpCside[fNparticle];
368 p[1] = fPypCside[fNparticle];
369 p[2] = fPzpCside[fNparticle];
370 pdgCode = fpPDGCode;
371// printf(" pc%d p sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
372 }
373 fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
374 } //**********************************************************************************************
375 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
376 if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
377 p[0] = fPxppAside[fNparticle];
378 p[1] = fPyppAside[fNparticle];
379 p[2] = fPzppAside[fNparticle];
380 pdgCode = fppPDGCode;
381// printf(" pc%d pi+ sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
382 }
383 if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
384 p[0] = fPxppCside[fNparticle];
385 p[1] = fPyppCside[fNparticle];
386 p[2] = fPzppCside[fNparticle];
387 pdgCode = fppPDGCode;
388// printf(" pc%d pi+ sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
389 }
390 if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
391 p[0] = fPxpmAside[fNparticle];
392 p[1] = fPypmAside[fNparticle];
393 p[2] = fPzpmAside[fNparticle];
394 pdgCode = fpmPDGCode;
395// printf(" pc%d pi- sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
396 }
397 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
398 p[0] = fPxpmCside[fNparticle];
399 p[1] = fPypmCside[fNparticle];
400 p[2] = fPzpmCside[fNparticle];
401 pdgCode = fpmPDGCode;
402// printf(" pc%d pi- sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
403 }
404 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
405 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
406 p[0] = fPxp0Aside[fNparticle];
407 p[1] = fPyp0Aside[fNparticle];
408 p[2] = fPzp0Aside[fNparticle];
409 pdgCode = fp0PDGCode;
410// printf(" pc%d pi0 sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
411 }
412 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside &&
413 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
414 p[0] = fPxp0Cside[fNparticle];
415 p[1] = fPyp0Cside[fNparticle];
416 p[2] = fPzp0Cside[fNparticle];
417 pdgCode = fp0PDGCode;
418// printf(" pc%d pi0 sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
419 }
420 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
421 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
422 p[0] = fPxetaAside[fNparticle];
423 p[1] = fPyetaAside[fNparticle];
424 p[2] = fPzetaAside[fNparticle];
425 pdgCode = fetaPDGCode;
426// printf(" pc%d eta sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
427 }
428 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside &&
429 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
430 p[0] = fPxetaCside[fNparticle];
431 p[1] = fPyetaCside[fNparticle];
432 p[2] = fPzetaCside[fNparticle];
433 pdgCode = fetaPDGCode;
434// printf(" pc%d eta sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
435 }
436 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
437 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
438 p[0] = fPxomegaAside[fNparticle];
439 p[1] = fPyomegaAside[fNparticle];
440 p[2] = fPzomegaAside[fNparticle];
441 pdgCode = fomegaPDGCode;
442// printf(" pc%d omega sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
443 }
444 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
445 && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
446 p[0] = fPxomegaCside[fNparticle];
447 p[1] = fPyomegaCside[fNparticle];
448 p[2] = fPzomegaCside[fNparticle];
449 pdgCode = fomegaPDGCode;
450// printf(" pc%d omega sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
451 }
452
453 }
454
455 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
456 Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
457 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
458
459 if(p[3]<=amass){
460 Warning("Generate","Particle %d E = %f GeV mass = %f GeV ",pdgCode,p[3],amass);
461 }
462
463 //printf(" Pc %d: PDGcode %d p(%1.2f, %1.2f, %1.2f, %1.3f)\n",
464 // fNparticle,pdgCode,p[0], p[1], p[2], p[3]);
465
466 TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
467 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
468 if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
469 fNparticle++;
470 return particle;
471}