Fix for bug #73998: AliESDZDC.cxx changes to be committed in trunk and ported to...
[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>
27
28#include "AliGenReaderEMD.h"
29
706938e6 30ClassImp(AliGenReaderEMD)
f1d1affa 31
1c56e311 32AliGenReaderEMD::AliGenReaderEMD():
33 fStartEvent(0),
34 fNcurrent(0),
35 fNparticle(0),
36 fTreeNtuple(0),
37 fIPSide(0),
38 fPcToTrack(0),
39 fNnLeft(0),
40 fEnLeft(0),
41 fNnRight(0),
42 fEnRight(0),
43 fNpLeft(0),
44 fEtapLeft(0),
45 fNpRight(0),
46 fEtapRight(0)
f1d1affa 47{
48// Default constructor
f1d1affa 49}
50
1c56e311 51AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
52 AliGenReader(reader),
53 fStartEvent(0),
54 fNcurrent(0),
55 fNparticle(0),
56 fTreeNtuple(0),
57 fIPSide(0),
58 fPcToTrack(0),
59 fNnLeft(0),
60 fEnLeft(0),
61 fNnRight(0),
62 fEnRight(0),
63 fNpLeft(0),
64 fEtapLeft(0),
65 fNpRight(0),
66 fEtapRight(0)
67{
68 // Copy Constructor
69 reader.Copy(*this);
70}
f1d1affa 71 // -----------------------------------------------------------------------------------
72AliGenReaderEMD::~AliGenReaderEMD()
73{
74 delete fTreeNtuple;
75}
76
77// -----------------------------------------------------------------------------------
78AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
79{
80// Assignment operator
81 rhs.Copy(*this);
82 return *this;
83}
84
85// -----------------------------------------------------------------------------------
86void AliGenReaderEMD::Copy(TObject&) const
87{
88 //
89 // Copy
90 //
91 Fatal("Copy","Not implemented!\n");
92}
93
94// -----------------------------------------------------------------------------------
95void AliGenReaderEMD::Init()
96{
97//
98// Reset the existing file environment and open a new root file
99
100 TFile *pFile=0;
101 if (!pFile) {
102 pFile = new TFile(fFileName);
103 pFile->cd();
104 printf("\n %s file opened\n\n", fFileName);
105 }
106 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
107 fNcurrent = fStartEvent;
108
109 TTree *Ntu=fTreeNtuple;
110 //
111 // Set branch addresses
112 // **** neutrons
113 Ntu->SetBranchAddress("Nleft",&fNnLeft);
114 Ntu->SetBranchAddress("Eleft",&fEnLeft);
2524c56f 115 Ntu->SetBranchAddress("Pxl", fPxnLeft);
116 Ntu->SetBranchAddress("Pyl", fPynLeft);
117 Ntu->SetBranchAddress("Pzl", fPznLeft);
f1d1affa 118 Ntu->SetBranchAddress("Nright",&fNnRight);
119 Ntu->SetBranchAddress("Eright",&fEnRight);
2524c56f 120 Ntu->SetBranchAddress("Pxr", fPxnRight);
121 Ntu->SetBranchAddress("Pyr", fPynRight);
122 Ntu->SetBranchAddress("Pzr", fPznRight);
f1d1affa 123 // **** protons
124 Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
125 Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
2524c56f 126 Ntu->SetBranchAddress("Pxl_p", fPxpLeft);
127 Ntu->SetBranchAddress("Pyl_p", fPypLeft);
128 Ntu->SetBranchAddress("Pzl_p", fPzpLeft);
f1d1affa 129 Ntu->SetBranchAddress("Nright_p",&fNpRight);
130 Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
2524c56f 131 Ntu->SetBranchAddress("Pxr_p", fPxpRight);
132 Ntu->SetBranchAddress("Pyr_p", fPypRight);
133 Ntu->SetBranchAddress("Pzr_p", fPzpRight);
f1d1affa 134}
135
136// -----------------------------------------------------------------------------------
137Int_t AliGenReaderEMD::NextEvent()
138{
139 // Read the next event
140 Int_t nTracks=0;
141
142 TFile* pFile = fTreeNtuple->GetCurrentFile();
143 pFile->cd();
144
145
146 Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
147 if(fNcurrent < nentries) {
148 fTreeNtuple->GetEvent(fNcurrent);
149 printf("\n *** Event %d read ***\n\n",fNcurrent);
150 fNcurrent++;
151 //printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
152 //printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
153 //
154 // **** IPSIde =0->RIGHT side, =1->LEFT side
155 // #### fPcToTrack =0->neutrons, =1->protons
156 if(fIPSide==0){
157 if(fPcToTrack==0){
ef064345 158 printf("\n \t \t Tracking %d neutrons emitted on C side\n\n", fNnRight);
f1d1affa 159 nTracks = fNnRight;
160 }
161 else if(fPcToTrack==1){
ef064345 162 printf("\n \t \t Tracking %d protons emitted on C side\n\n", fNpRight);
f1d1affa 163 nTracks = fNpRight;
164 }
165 }
166 else if(fIPSide==1){
167 if(fPcToTrack==0){
ef064345 168 printf("\n \t \t Tracking %d neutrons emitted on A side\n", fNnLeft);
f1d1affa 169 nTracks = fNnLeft;
170 }
171 else if(fPcToTrack==1){
ef064345 172 printf("\n \t \t Tracking %d protons emitted on A side\n", fNpLeft);
f1d1affa 173 nTracks = fNpLeft;
174 }
175 }
176 fNparticle = 0;
177 return nTracks;
178 }
179
180 return 0;
181}
182
183// -----------------------------------------------------------------------------------
184TParticle* AliGenReaderEMD::NextParticle()
185{
186 // Read the next particle
187 Float_t p[4];
ef064345 188
189 Int_t ipart=0;
190 if(fPcToTrack==0) ipart = kNeutron;
191 else if(fPcToTrack==1) ipart = kProton;
192 Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
193
f1d1affa 194 p[0] = fPxnRight[fNparticle];
195 p[1] = fPynRight[fNparticle];
196 p[2] = fPznRight[fNparticle];
197 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
198 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
199 //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
200
201 if(p[3]<=amass)
202 Warning("Generate","Particle %d E = %f mass = %f \n",ipart,p[3],amass);
203
204 TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
205 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
206 //fNcurrent++;
207 fNparticle++;
208 return particle;
209}