12ea5e503f3e45a64fb5c7d22f6e62ea5fa6be32
[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
33   // -----------------------------------------------------------------------------------
34 AliGenReaderEMD::AliGenReaderEMD() 
35 {
36 // Default constructor
37     fStartEvent  = 0;
38     fTreeNtuple  = 0;
39     fIPSide      = 0;
40     fPcToTrack = 0;
41 }
42
43   // -----------------------------------------------------------------------------------
44 AliGenReaderEMD::~AliGenReaderEMD()
45 {
46     delete fTreeNtuple;
47 }
48
49 // -----------------------------------------------------------------------------------
50 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
51 {
52 // Assignment operator
53     rhs.Copy(*this);
54     return *this;
55 }
56
57 // -----------------------------------------------------------------------------------
58 void AliGenReaderEMD::Copy(TObject&) const
59 {
60     //
61     // Copy 
62     //
63     Fatal("Copy","Not implemented!\n");
64 }
65
66 // -----------------------------------------------------------------------------------
67 void AliGenReaderEMD::Init() 
68 {
69 //
70 // Reset the existing file environment and open a new root file
71     
72     TFile *pFile=0;
73     if (!pFile) {
74         pFile = new TFile(fFileName);
75         pFile->cd();
76         printf("\n %s file opened\n\n", fFileName);
77     }
78     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
79     fNcurrent = fStartEvent;
80
81     TTree *Ntu=fTreeNtuple;
82     //
83     // Set branch addresses
84     // **** neutrons
85     Ntu->SetBranchAddress("Nleft",&fNnLeft);
86     Ntu->SetBranchAddress("Eleft",&fEnLeft);
87     Ntu->SetBranchAddress("Pxl",  &fPxnLeft);
88     Ntu->SetBranchAddress("Pyl",  &fPynLeft);
89     Ntu->SetBranchAddress("Pzl",  &fPznLeft);
90     Ntu->SetBranchAddress("Nright",&fNnRight);
91     Ntu->SetBranchAddress("Eright",&fEnRight);
92     Ntu->SetBranchAddress("Pxr",   &fPxnRight);
93     Ntu->SetBranchAddress("Pyr",   &fPynRight);
94     Ntu->SetBranchAddress("Pzr",   &fPznRight);
95     // **** protons
96     Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
97     Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
98     Ntu->SetBranchAddress("Pxl_p",  &fPxpLeft);
99     Ntu->SetBranchAddress("Pyl_p",  &fPypLeft);
100     Ntu->SetBranchAddress("Pzl_p",  &fPzpLeft);
101     Ntu->SetBranchAddress("Nright_p",&fNpRight);
102     Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
103     Ntu->SetBranchAddress("Pxr_p",   &fPxpRight);
104     Ntu->SetBranchAddress("Pyr_p",   &fPypRight);
105     Ntu->SetBranchAddress("Pzr_p",   &fPzpRight);
106 }
107
108 // -----------------------------------------------------------------------------------
109 Int_t AliGenReaderEMD::NextEvent() 
110 {
111     // Read the next event  
112     Int_t nTracks=0;
113     
114     TFile* pFile = fTreeNtuple->GetCurrentFile();
115     pFile->cd();
116     
117
118     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
119     if(fNcurrent < nentries) {
120         fTreeNtuple->GetEvent(fNcurrent);
121         printf("\n *** Event %d read ***\n\n",fNcurrent);
122         fNcurrent++;
123         //printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
124         //printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
125         //
126         // **** IPSIde          =0->RIGHT side, =1->LEFT side  
127         // #### fPcToTrack      =0->neutrons, =1->protons
128         if(fIPSide==0){
129           if(fPcToTrack==0){
130             printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
131             nTracks    = fNnRight;
132           }
133           else if(fPcToTrack==1){
134             printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
135             nTracks    = fNpRight;
136           }
137         }
138         else if(fIPSide==1){
139           if(fPcToTrack==0){
140             printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
141             nTracks    = fNnLeft;
142           }
143           else if(fPcToTrack==1){
144             printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
145             nTracks    = fNpLeft;
146           }
147         }
148         fNparticle = 0;
149         return nTracks;
150     }
151
152     return 0;
153 }
154
155 // -----------------------------------------------------------------------------------
156 TParticle* AliGenReaderEMD::NextParticle() 
157 {
158     // Read the next particle
159     Float_t p[4];
160     Int_t ipart = kNeutron;
161     Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
162     p[0] = fPxnRight[fNparticle];
163     p[1] = fPynRight[fNparticle];
164     p[2] = fPznRight[fNparticle];
165     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
166     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
167     //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
168     
169     if(p[3]<=amass) 
170       Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
171
172     TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1, 
173         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
174     //fNcurrent++;
175     fNparticle++;
176     return particle;
177 }