Add the number of local boards
[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
28 #include "AliGenReaderEMD.h"
29
30 ClassImp(AliGenReaderEMD)
31
32 AliGenReaderEMD::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)
47 {
48 // Default constructor
49 }
50
51 AliGenReaderEMD::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 }
71   // -----------------------------------------------------------------------------------
72 AliGenReaderEMD::~AliGenReaderEMD()
73 {
74     delete fTreeNtuple;
75 }
76
77 // -----------------------------------------------------------------------------------
78 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
79 {
80 // Assignment operator
81     rhs.Copy(*this);
82     return *this;
83 }
84
85 // -----------------------------------------------------------------------------------
86 void AliGenReaderEMD::Copy(TObject&) const
87 {
88     //
89     // Copy 
90     //
91     Fatal("Copy","Not implemented!\n");
92 }
93
94 // -----------------------------------------------------------------------------------
95 void 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);
115     Ntu->SetBranchAddress("Pxl",  fPxnLeft);
116     Ntu->SetBranchAddress("Pyl",  fPynLeft);
117     Ntu->SetBranchAddress("Pzl",  fPznLeft);
118     Ntu->SetBranchAddress("Nright",&fNnRight);
119     Ntu->SetBranchAddress("Eright",&fEnRight);
120     Ntu->SetBranchAddress("Pxr",   fPxnRight);
121     Ntu->SetBranchAddress("Pyr",   fPynRight);
122     Ntu->SetBranchAddress("Pzr",   fPznRight);
123     // **** protons
124     Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
125     Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
126     Ntu->SetBranchAddress("Pxl_p",  fPxpLeft);
127     Ntu->SetBranchAddress("Pyl_p",  fPypLeft);
128     Ntu->SetBranchAddress("Pzl_p",  fPzpLeft);
129     Ntu->SetBranchAddress("Nright_p",&fNpRight);
130     Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
131     Ntu->SetBranchAddress("Pxr_p",   fPxpRight);
132     Ntu->SetBranchAddress("Pyr_p",   fPypRight);
133     Ntu->SetBranchAddress("Pzr_p",   fPzpRight);
134 }
135
136 // -----------------------------------------------------------------------------------
137 Int_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){
158             printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
159             nTracks    = fNnRight;
160           }
161           else if(fPcToTrack==1){
162             printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
163             nTracks    = fNpRight;
164           }
165         }
166         else if(fIPSide==1){
167           if(fPcToTrack==0){
168             printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
169             nTracks    = fNnLeft;
170           }
171           else if(fPcToTrack==1){
172             printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
173             nTracks    = fNpLeft;
174           }
175         }
176         fNparticle = 0;
177         return nTracks;
178     }
179
180     return 0;
181 }
182
183 // -----------------------------------------------------------------------------------
184 TParticle* AliGenReaderEMD::NextParticle() 
185 {
186     // Read the next particle
187     Float_t p[4];
188     Int_t ipart = kNeutron;
189     Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
190     p[0] = fPxnRight[fNparticle];
191     p[1] = fPynRight[fNparticle];
192     p[2] = fPznRight[fNparticle];
193     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
194     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
195     //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
196     
197     if(p[3]<=amass) 
198       Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
199
200     TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1, 
201         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
202     //fNcurrent++;
203     fNparticle++;
204     return particle;
205 }