Coverity fix.
[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.cxx
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
31 ClassImp(AliGenReaderEMD)
32
33 AliGenReaderEMD::AliGenReaderEMD():
34     fStartEvent(0),
35     fNcurrent(0),  
36     fNparticle(0), 
37     fTreeNtuple(0),
38     fIPSide(0),
39     fPcToTrack(0),
40     fNnAside(0),
41     fEnAside(0),
42     fNnCside(0),
43     fEnCside(0),
44     fNpAside(0),
45     fEtapAside(0),
46     fNpCside(0),
47     fEtapCside(0)
48 {
49 // Default constructor
50 }
51
52 AliGenReaderEMD::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),
60     fNnAside(0),
61     fEnAside(0),
62     fNnCside(0),
63     fEnCside(0),
64     fNpAside(0),
65     fEtapAside(0),
66     fNpCside(0),
67     fEtapCside(0)
68 {
69     // Copy Constructor
70     reader.Copy(*this);
71 }
72   // -----------------------------------------------------------------------------------
73 AliGenReaderEMD::~AliGenReaderEMD()
74 {
75     delete fTreeNtuple;
76 }
77
78 // -----------------------------------------------------------------------------------
79 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
80 {
81 // Assignment operator
82     rhs.Copy(*this);
83     return *this;
84 }
85
86 // -----------------------------------------------------------------------------------
87 void AliGenReaderEMD::Copy(TObject&) const
88 {
89     //
90     // Copy 
91     //
92     Fatal("Copy","Not implemented!\n");
93 }
94
95 // -----------------------------------------------------------------------------------
96 void 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
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);
124     // **** protons
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);
135 }
136
137 // -----------------------------------------------------------------------------------
138 Int_t AliGenReaderEMD::NextEvent() 
139 {
140     // Read the next event  
141     Int_t nTracks=0;
142     fNparticle = 0;
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);
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);
154         //
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;
165           }
166           else if(fIPSide==2){
167             printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
168             nTracks    = fNnAside;
169           }
170         }
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;
175           }
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;
183           }
184         }
185         fNcurrent++;
186         return nTracks;
187     }
188
189     return 0;
190 }
191
192 // -----------------------------------------------------------------------------------
193 TParticle* AliGenReaderEMD::NextParticle() 
194 {
195     // Read the next particle
196     Float_t p[4]={0.,0.,0.,0.};
197     
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
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     } 
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);
253     //printf(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, p[0],p[1],p[2]);
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.);
260     particle->SetBit(kTransportBit);
261     fNparticle++;
262     return particle;
263 }