Fixed memory leak
[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),
38 fIPSide(0),
39 fPcToTrack(0),
2078b0f7 40 fNnAside(0),
41 fEnAside(0),
42 fNnCside(0),
43 fEnCside(0),
44 fNpAside(0),
45 fEtapAside(0),
46 fNpCside(0),
47 fEtapCside(0)
f1d1affa 48{
49// Default constructor
f1d1affa 50}
51
1c56e311 52AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
53 AliGenReader(reader),
54 fStartEvent(0),
55 fNcurrent(0),
56 fNparticle(0),
57 fTreeNtuple(0),
58 fIPSide(0),
59 fPcToTrack(0),
2078b0f7 60 fNnAside(0),
61 fEnAside(0),
62 fNnCside(0),
63 fEnCside(0),
64 fNpAside(0),
65 fEtapAside(0),
66 fNpCside(0),
67 fEtapCside(0)
1c56e311 68{
69 // Copy Constructor
70 reader.Copy(*this);
71}
f1d1affa 72 // -----------------------------------------------------------------------------------
73AliGenReaderEMD::~AliGenReaderEMD()
74{
75 delete fTreeNtuple;
76}
77
78// -----------------------------------------------------------------------------------
79AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
80{
81// Assignment operator
82 rhs.Copy(*this);
83 return *this;
84}
85
86// -----------------------------------------------------------------------------------
87void AliGenReaderEMD::Copy(TObject&) const
88{
89 //
90 // Copy
91 //
92 Fatal("Copy","Not implemented!\n");
93}
94
95// -----------------------------------------------------------------------------------
96void AliGenReaderEMD::Init()
97{
98//
99// Reset the existing file environment and open a new root file
100
101 TFile *pFile=0;
102 if (!pFile) {
103 pFile = new TFile(fFileName);
104 pFile->cd();
105 printf("\n %s file opened\n\n", fFileName);
106 }
107 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
108 fNcurrent = fStartEvent;
109
110 TTree *Ntu=fTreeNtuple;
111 //
112 // Set branch addresses
113 // **** neutrons
2078b0f7 114 Ntu->SetBranchAddress("Nleft",&fNnAside);
115 Ntu->SetBranchAddress("Eleft",&fEnAside);
116 Ntu->SetBranchAddress("Pxl", fPxnAside);
117 Ntu->SetBranchAddress("Pyl", fPynAside);
118 Ntu->SetBranchAddress("Pzl", fPznAside);
119 Ntu->SetBranchAddress("Nright",&fNnCside);
120 Ntu->SetBranchAddress("Eright",&fEnCside);
121 Ntu->SetBranchAddress("Pxr", fPxnCside);
122 Ntu->SetBranchAddress("Pyr", fPynCside);
123 Ntu->SetBranchAddress("Pzr", fPznCside);
f1d1affa 124 // **** protons
2078b0f7 125 Ntu->SetBranchAddress("Nleft_p",&fNpAside);
126 Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
127 Ntu->SetBranchAddress("Pxl_p", fPxpAside);
128 Ntu->SetBranchAddress("Pyl_p", fPypAside);
129 Ntu->SetBranchAddress("Pzl_p", fPzpAside);
130 Ntu->SetBranchAddress("Nright_p",&fNpCside);
131 Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
132 Ntu->SetBranchAddress("Pxr_p", fPxpCside);
133 Ntu->SetBranchAddress("Pyr_p", fPypCside);
134 Ntu->SetBranchAddress("Pzr_p", fPzpCside);
f1d1affa 135}
136
137// -----------------------------------------------------------------------------------
138Int_t AliGenReaderEMD::NextEvent()
139{
140 // Read the next event
141 Int_t nTracks=0;
2078b0f7 142 fNparticle = 0;
f1d1affa 143
144 TFile* pFile = fTreeNtuple->GetCurrentFile();
145 pFile->cd();
146
147
148 Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
149 if(fNcurrent < nentries) {
150 fTreeNtuple->GetEvent(fNcurrent);
2078b0f7 151 printf("\n *** Event %d read ***\n",fNcurrent);
152 //printf("\t A side-> %d neutrons and %d protons emitted\n", fNnAside, fNpAside);
153 //printf("\t C side-> %d neutrons and %d protons emitted\n\n", fNnCside, fNpCside);
f1d1affa 154 //
2078b0f7 155 // **** IPSIde = 0->both sides, 1->only Cside, 2->only A side
156 // #### fPcToTrack = 0->neutrons, 1->protons
157 if(fPcToTrack==0){ // neutrons
158 if(fIPSide==0){
159 printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
160 nTracks = fNnCside+fNnAside;
161 }
162 else if(fIPSide==1){
163 printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
164 nTracks = fNnCside;
f1d1affa 165 }
2078b0f7 166 else if(fIPSide==2){
167 printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
168 nTracks = fNnAside;
f1d1affa 169 }
170 }
2078b0f7 171 else if(fPcToTrack==1){ //protons
172 if(fIPSide==0){
173 printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
174 nTracks = fNpCside+fNpAside;
f1d1affa 175 }
2078b0f7 176 else if(fIPSide==1){
177 printf("\t Tracking %d protons emitted on C side\n", fNpCside);
178 nTracks = fNpCside;
179 }
180 else if(fIPSide==2){
181 printf("\t Tracking %d protons emitted on A side\n", fNpAside);
182 nTracks = fNpAside;
f1d1affa 183 }
184 }
2078b0f7 185 fNcurrent++;
f1d1affa 186 return nTracks;
187 }
188
189 return 0;
190}
191
192// -----------------------------------------------------------------------------------
193TParticle* AliGenReaderEMD::NextParticle()
194{
195 // Read the next particle
2078b0f7 196 Float_t p[4]={0.,0.,0.,0.};
197
ef064345 198 Int_t ipart=0;
199 if(fPcToTrack==0) ipart = kNeutron;
200 else if(fPcToTrack==1) ipart = kProton;
201 Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
202
2078b0f7 203 if(fPcToTrack==0){
204 if(fIPSide==0){
205 if(fNparticle<fNnCside){
206 p[0] = fPxnCside[fNparticle];
207 p[1] = fPynCside[fNparticle];
208 p[2] = fPznCside[fNparticle];
209 }
210 else{
211 p[0] = fPxnAside[fNparticle-fNnCside];
212 p[1] = fPynAside[fNparticle-fNnCside];
213 p[2] = fPznAside[fNparticle-fNnCside];
214 }
215 }
216 else if(fIPSide==1){
217 p[0] = fPxnCside[fNparticle];
218 p[1] = fPynCside[fNparticle];
219 p[2] = fPznCside[fNparticle];
220 }
221 else if(fIPSide==2){
222 p[0] = fPxnAside[fNparticle-fNpCside];
223 p[1] = fPynAside[fNparticle-fNpCside];
224 p[2] = fPznAside[fNparticle-fNpCside];
225 }
226 }
227 else if(fPcToTrack==1){
228 if(fIPSide==0){
229 if(fNparticle<fNnCside){
230 p[0] = fPxpCside[fNparticle];
231 p[1] = fPypCside[fNparticle];
232 p[2] = fPzpCside[fNparticle];
233 }
234 else{
235 p[0] = fPxpAside[fNparticle];
236 p[1] = fPypAside[fNparticle];
237 p[2] = fPzpAside[fNparticle];
238 }
239 }
240 else if(fIPSide==1){
241 p[0] = fPxpCside[fNparticle];
242 p[1] = fPypCside[fNparticle];
243 p[2] = fPzpCside[fNparticle];
244 }
245 else if(fIPSide==2){
246 p[0] = fPxpAside[fNparticle];
247 p[1] = fPypAside[fNparticle];
248 p[2] = fPzpAside[fNparticle];
249 }
250 }
f1d1affa 251 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
252 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
2078b0f7 253 //printf(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, p[0],p[1],p[2]);
f1d1affa 254
255 if(p[3]<=amass)
256 Warning("Generate","Particle %d E = %f mass = %f \n",ipart,p[3],amass);
257
258 TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
259 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
2078b0f7 260 particle->SetBit(kTransportBit);
f1d1affa 261 fNparticle++;
262 return particle;
263}