]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenReaderEMD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Beam crossing angle and divergence (Chiara)
[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.cxx
CommitLineData
f1d1affa 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>
f1d1affa 27#include "AliGenReaderEMD.h"
2078b0f7 28#include "AliStack.h"
29
f1d1affa 30
706938e6 31ClassImp(AliGenReaderEMD)
f1d1affa 32
1c56e311 33AliGenReaderEMD::AliGenReaderEMD():
34 fStartEvent(0),
35 fNcurrent(0),
36 fNparticle(0),
37 fTreeNtuple(0),
1c56e311 38 fPcToTrack(0),
a16b785f 39 fOffset(0),
2078b0f7 40 fNnAside(0),
41 fEnAside(0),
a16b785f 42 fnPDGCode(0),
2078b0f7 43 fNnCside(0),
44 fEnCside(0),
45 fNpAside(0),
46 fEtapAside(0),
a16b785f 47 fpPDGCode(0),
2078b0f7 48 fNpCside(0),
a16b785f 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)
f1d1affa 75{
a16b785f 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");
f1d1affa 102}
103
a16b785f 104
1c56e311 105AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
106 AliGenReader(reader),
107 fStartEvent(0),
108 fNcurrent(0),
109 fNparticle(0),
110 fTreeNtuple(0),
1c56e311 111 fPcToTrack(0),
a16b785f 112 fOffset(0),
2078b0f7 113 fNnAside(0),
114 fEnAside(0),
a16b785f 115 fnPDGCode(0),
2078b0f7 116 fNnCside(0),
117 fEnCside(0),
118 fNpAside(0),
119 fEtapAside(0),
a16b785f 120 fpPDGCode(0),
2078b0f7 121 fNpCside(0),
a16b785f 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)
1c56e311 148{
149 // Copy Constructor
5eebf646 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 }
1c56e311 172 reader.Copy(*this);
173}
f1d1affa 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();
a16b785f 207 printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
f1d1affa 208 }
209 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
210 fNcurrent = fStartEvent;
211
212 TTree *Ntu=fTreeNtuple;
213 //
214 // Set branch addresses
215 // **** neutrons
2078b0f7 216 Ntu->SetBranchAddress("Nleft",&fNnAside);
217 Ntu->SetBranchAddress("Eleft",&fEnAside);
a16b785f 218 Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
2078b0f7 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);
f1d1affa 227 // **** protons
2078b0f7 228 Ntu->SetBranchAddress("Nleft_p",&fNpAside);
229 Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
a16b785f 230 Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
2078b0f7 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);
a16b785f 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);
f1d1affa 299}
300
301// -----------------------------------------------------------------------------------
302Int_t AliGenReaderEMD::NextEvent()
303{
304 // Read the next event
305 Int_t nTracks=0;
a16b785f 306 fNparticle = 0; fOffset=0;
f1d1affa 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);
a16b785f 315 if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
f1d1affa 316 //
a16b785f 317 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
318 nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
f1d1affa 319 }
a16b785f 320 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
321 nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
322 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
f1d1affa 323 }
2078b0f7 324 fNcurrent++;
a16b785f 325 printf("\t #### Putting %d particles in the stack\n", nTracks);
067ae68e 326 /*if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t %d+%d neutrons, %d+%d protons\n",
a16b785f 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,
067ae68e 330 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);*/
f1d1affa 331 return nTracks;
332 }
333
334 return 0;
335}
336
337// -----------------------------------------------------------------------------------
338TParticle* AliGenReaderEMD::NextParticle()
339{
340 // Read the next particle
2078b0f7 341 Float_t p[4]={0.,0.,0.,0.};
a16b785f 342 int pdgCode=0;
2078b0f7 343
a16b785f 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]);
2078b0f7 351 }
a16b785f 352 else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
2078b0f7 353 p[0] = fPxnCside[fNparticle];
354 p[1] = fPynCside[fNparticle];
355 p[2] = fPznCside[fNparticle];
a16b785f 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]);
2078b0f7 358 }
a16b785f 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]);
2078b0f7 365 }
a16b785f 366 else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
2078b0f7 367 p[0] = fPxpCside[fNparticle];
368 p[1] = fPypCside[fNparticle];
369 p[2] = fPzpCside[fNparticle];
a16b785f 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]);
2078b0f7 372 }
a16b785f 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]);
2078b0f7 419 }
a16b785f 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
2078b0f7 453 }
a16b785f 454
f1d1affa 455 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
a16b785f 456 Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
f1d1affa 457 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
f1d1affa 458
067ae68e 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]);
a16b785f 465
466 TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
f1d1affa 467 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
067ae68e 468 if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
f1d1affa 469 fNparticle++;
470 return particle;
471}