Line data Source code
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 "AliGenReadersEMD.h"
28 : #include "AliStack.h"
29 :
30 :
31 6 : ClassImp(AliGenReadersEMD)
32 :
33 0 : AliGenReadersEMD::AliGenReadersEMD():
34 0 : fStartEvent(0),
35 0 : fNcurrent(0),
36 0 : fNparticle(0),
37 0 : fTreeNtuple(0),
38 0 : fPcToTrack(2),
39 0 : fNneu(0),
40 0 : fEneu(0),
41 0 : fPxfrag(0),
42 0 : fPyfrag(0),
43 0 : fPzfrag(0),
44 0 : fAfrag(0),
45 0 : fZfrag(0),
46 0 : fNpro(0),
47 0 : fEpro(0)
48 0 : {
49 : // Std constructor
50 0 : for(int i=0; i<70; i++){
51 0 : fPxneu[i] = fPyneu[i] = fPzneu[i] = 0.;
52 0 : if(i<50){
53 0 : fPxpro[i] = fPypro[i] = fPzpro[i] = 0.;
54 0 : }
55 : }
56 : //
57 0 : if(fPcToTrack==kAll) printf("\n\t *** AliGenReadersEMD will track n, p and fragments \n\n");
58 0 : else if(fPcToTrack==kNucleons) printf("\n\t *** AliGenReadersEMD will track all nucleons\n\n");
59 0 : else if(fPcToTrack==kOnlyNeutrons) printf("\n\t *** AliGenReadersEMD will track only neutrons\n\n");
60 0 : }
61 :
62 :
63 : AliGenReadersEMD::AliGenReadersEMD(const AliGenReadersEMD &reader):
64 0 : AliGenReader(reader),
65 0 : fStartEvent(0),
66 0 : fNcurrent(0),
67 0 : fNparticle(0),
68 0 : fTreeNtuple(0),
69 0 : fPcToTrack(2),
70 0 : fNneu(0),
71 0 : fEneu(0),
72 0 : fPxfrag(0),
73 0 : fPyfrag(0),
74 0 : fPzfrag(0),
75 0 : fAfrag(0),
76 0 : fZfrag(0),
77 0 : fNpro(0),
78 0 : fEpro(0)
79 0 : {
80 : // Copy Constructor
81 0 : for(int i=0; i<70; i++){
82 0 : fPxneu[i] = fPyneu[i] = fPzneu[i] = 0.;
83 0 : if(i<50){
84 0 : fPxpro[i] = fPypro[i] = fPzpro[i] = 0.;
85 0 : }
86 : }
87 0 : reader.Copy(*this);
88 0 : }
89 : // -----------------------------------------------------------------------------------
90 : AliGenReadersEMD::~AliGenReadersEMD()
91 0 : {
92 0 : delete fTreeNtuple;
93 0 : }
94 :
95 : // -----------------------------------------------------------------------------------
96 : AliGenReadersEMD& AliGenReadersEMD::operator=(const AliGenReadersEMD& rhs)
97 : {
98 : // Assignment operator
99 0 : rhs.Copy(*this);
100 0 : return *this;
101 : }
102 :
103 : // -----------------------------------------------------------------------------------
104 : void AliGenReadersEMD::Copy(TObject&) const
105 : {
106 : //
107 : // Copy
108 : //
109 0 : Fatal("Copy","Not implemented!\n");
110 0 : }
111 :
112 : // -----------------------------------------------------------------------------------
113 : void AliGenReadersEMD::Init()
114 : {
115 : //
116 : // Reset the existing file environment and open a new root file
117 :
118 : TFile *pFile=0;
119 0 : if (!pFile) {
120 0 : pFile = new TFile(fFileName);
121 0 : pFile->cd();
122 0 : printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
123 0 : }
124 0 : fTreeNtuple = (TTree*)gDirectory->Get("h2034");
125 0 : fNcurrent = fStartEvent;
126 :
127 0 : TTree *Ntu=fTreeNtuple;
128 : //
129 : // Set branch addresses
130 : // **** neutrons
131 0 : Ntu->SetBranchAddress("N_na50",&fNneu);
132 0 : Ntu->SetBranchAddress("E_na50",&fEneu);
133 0 : Ntu->SetBranchAddress("Pxna50", fPxneu);
134 0 : Ntu->SetBranchAddress("Pyna50", fPyneu);
135 0 : Ntu->SetBranchAddress("Pzna50", fPzneu);
136 0 : Ntu->SetBranchAddress("Px_mhfrag", &fPxfrag);
137 0 : Ntu->SetBranchAddress("Py_mhfrag", &fPyfrag);
138 0 : Ntu->SetBranchAddress("Pz_mhfrag", &fPzfrag);
139 0 : Ntu->SetBranchAddress("A_mhfrag", &fAfrag);
140 0 : Ntu->SetBranchAddress("Z_mhfrag", &fZfrag);
141 0 : Ntu->SetBranchAddress("N_na50_p",&fNpro);
142 0 : Ntu->SetBranchAddress("E_na50_p",&fEpro);
143 0 : Ntu->SetBranchAddress("Pxna50_p", fPxpro);
144 0 : Ntu->SetBranchAddress("Pyna50_p", fPypro);
145 0 : Ntu->SetBranchAddress("Pzna50_p", fPzpro);
146 0 : }
147 :
148 : // -----------------------------------------------------------------------------------
149 : Int_t AliGenReadersEMD::NextEvent()
150 : {
151 : // Read the next event
152 : Int_t nTracks=0;
153 0 : fNparticle = 0;
154 :
155 0 : TFile* pFile = fTreeNtuple->GetCurrentFile();
156 0 : pFile->cd();
157 :
158 0 : Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
159 0 : if(fNcurrent < nentries) {
160 0 : fTreeNtuple->GetEvent(fNcurrent);
161 0 : if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
162 : //
163 0 : if(fPcToTrack==kAll){ // all
164 0 : nTracks = fNneu+fNpro+1;
165 0 : }
166 0 : else if(fPcToTrack==kNucleons){ // nucleons
167 0 : nTracks = fNneu+fNpro;
168 0 : }
169 0 : else if(fPcToTrack==kOnlyNeutrons){ // neutrons
170 0 : nTracks = fNneu;
171 0 : }
172 0 : fNcurrent++;
173 0 : printf("\t #### Putting %d particles in the stack\n", nTracks);
174 0 : return nTracks;
175 : }
176 :
177 0 : return 0;
178 0 : }
179 :
180 : // -----------------------------------------------------------------------------------
181 : TParticle* AliGenReadersEMD::NextParticle()
182 : {
183 : // Read the next particle
184 : Float_t p[4]={0.,0.,0.,0.};
185 : int pdgCode=0;
186 :
187 0 : if(fNparticle<fNneu){
188 0 : p[0] = fPxneu[fNparticle];
189 0 : p[1] = fPyneu[fNparticle];
190 0 : p[2] = fPzneu[fNparticle];
191 : pdgCode = 2112;
192 : // printf(" pc%d n: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
193 0 : }
194 :
195 0 : if(fPcToTrack==kAll || fPcToTrack==kNucleons){
196 0 : if(fNparticle>=fNneu && fNparticle<(fNneu+fNpro)){
197 0 : p[0] = fPxpro[fNparticle];
198 0 : p[1] = fPypro[fNparticle];
199 0 : p[2] = fPzpro[fNparticle];
200 : pdgCode = 2212;
201 : // printf(" pc%d p: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
202 0 : }
203 : }
204 0 : if(fPcToTrack==kAll){
205 0 : if(fNparticle>=(fNneu+fNpro)){
206 0 : p[0] = fPxfrag;
207 0 : p[1] = fPyfrag;
208 0 : p[2] = fPzfrag;
209 0 : pdgCode = 1e10+1e6*fZfrag+10.*fAfrag;
210 : // printf(" pc%d fragment: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
211 0 : }
212 : }
213 :
214 0 : Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
215 0 : Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
216 0 : p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
217 :
218 0 : if(p[3]<=amass){
219 0 : Warning("Generate","Particle %d E = %f GeV mass = %f GeV ",pdgCode,p[3],amass);
220 0 : }
221 :
222 : //printf(" Pc %d: PDGcode %d p(%1.2f, %1.2f, %1.2f, %1.3f)\n",
223 : //fNparticle,pdgCode,p[0], p[1], p[2], p[3]);
224 :
225 0 : TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
226 0 : p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
227 0 : if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
228 0 : fNcurrent++;
229 0 : fNparticle++;
230 0 : return particle;
231 0 : }
232 :
233 : //___________________________________________________________
234 : void AliGenReadersEMD::RewindEvent()
235 : {
236 : // Go back to the first particle of the event
237 0 : fNparticle = 0;
238 0 : }
|
