1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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)
22 #include <TParticle.h>
24 #include <TVirtualMC.h>
25 #include <TDatabasePDG.h>
27 #include "AliGenReaderEMD.h"
31 ClassImp(AliGenReaderEMD)
33 AliGenReaderEMD::AliGenReaderEMD():
49 // Default constructor
52 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
72 // -----------------------------------------------------------------------------------
73 AliGenReaderEMD::~AliGenReaderEMD()
78 // -----------------------------------------------------------------------------------
79 AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
81 // Assignment operator
86 // -----------------------------------------------------------------------------------
87 void AliGenReaderEMD::Copy(TObject&) const
92 Fatal("Copy","Not implemented!\n");
95 // -----------------------------------------------------------------------------------
96 void AliGenReaderEMD::Init()
99 // Reset the existing file environment and open a new root file
103 pFile = new TFile(fFileName);
105 printf("\n %s file opened\n\n", fFileName);
107 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
108 fNcurrent = fStartEvent;
110 TTree *Ntu=fTreeNtuple;
112 // Set branch addresses
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);
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);
137 // -----------------------------------------------------------------------------------
138 Int_t AliGenReaderEMD::NextEvent()
140 // Read the next event
144 TFile* pFile = fTreeNtuple->GetCurrentFile();
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);
155 // **** IPSIde = 0->both sides, 1->only Cside, 2->only A side
156 // #### fPcToTrack = 0->neutrons, 1->protons
157 if(fPcToTrack==0){ // neutrons
159 printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
160 nTracks = fNnCside+fNnAside;
163 printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
167 printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
171 else if(fPcToTrack==1){ //protons
173 printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
174 nTracks = fNpCside+fNpAside;
177 printf("\t Tracking %d protons emitted on C side\n", fNpCside);
181 printf("\t Tracking %d protons emitted on A side\n", fNpAside);
192 // -----------------------------------------------------------------------------------
193 TParticle* AliGenReaderEMD::NextParticle()
195 // Read the next particle
196 Float_t p[4]={0.,0.,0.,0.};
199 if(fPcToTrack==0) ipart = kNeutron;
200 else if(fPcToTrack==1) ipart = kProton;
201 Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
205 if(fNparticle<fNnCside){
206 p[0] = fPxnCside[fNparticle];
207 p[1] = fPynCside[fNparticle];
208 p[2] = fPznCside[fNparticle];
211 p[0] = fPxnAside[fNparticle-fNnCside];
212 p[1] = fPynAside[fNparticle-fNnCside];
213 p[2] = fPznAside[fNparticle-fNnCside];
217 p[0] = fPxnCside[fNparticle];
218 p[1] = fPynCside[fNparticle];
219 p[2] = fPznCside[fNparticle];
222 p[0] = fPxnAside[fNparticle-fNpCside];
223 p[1] = fPynAside[fNparticle-fNpCside];
224 p[2] = fPznAside[fNparticle-fNpCside];
227 else if(fPcToTrack==1){
229 if(fNparticle<fNnCside){
230 p[0] = fPxpCside[fNparticle];
231 p[1] = fPypCside[fNparticle];
232 p[2] = fPzpCside[fNparticle];
235 p[0] = fPxpAside[fNparticle];
236 p[1] = fPypAside[fNparticle];
237 p[2] = fPzpAside[fNparticle];
241 p[0] = fPxpCside[fNparticle];
242 p[1] = fPypCside[fNparticle];
243 p[2] = fPzpCside[fNparticle];
246 p[0] = fPxpAside[fNparticle];
247 p[1] = fPypAside[fNparticle];
248 p[2] = fPzpAside[fNparticle];
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]);
256 Warning("Generate","Particle %d E = %f mass = %f \n",ipart,p[3],amass);
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);