Add version number as a global variable
[u/mrichter/AliRoot.git] / EVGEN / AliGenFLUKAsource.cxx
CommitLineData
fe4da5cc 1#include "AliGenFLUKAsource.h"
2#include "AliGenMUONlib.h"
3#include "AliMC.h"
4#include "AliRun.h"
5#include <TDirectory.h>
6#include <TFile.h>
7#include <TTree.h>
8#include <stdlib.h>
9 ClassImp(AliGenFLUKAsource)
10 AliGenFLUKAsource::AliGenFLUKAsource()
11 :AliGenerator(-1)
12{
13 //
14 fName="FLUKA";
15 fTitle="FLUKA Boundary Source";
16 // Read in all particle types by default
17 fIkine=6;
18 // Set maximum admitted age of particles to 1.e-05 by default
19 fAgeMax=1.e-05;
20 // Do not add weight
21 fAddWeight=1.;
22 // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
23 // from specific boundary source to the chamber (fZshift=4.5;),else there
24 // is no need to shift as it reads boundary source for the whole volume of
25 // the Muon Arm; the default value corresponds to boundary source for the
26 // whole volume of the MUON Arm
27 fZshift=0;
28 // Set the default file
29 fFileName="flukasource.root";
30
31 fTreeFluka=0;
32//
33// Read all particles
34 fNpart=-1;
35
36}
37
38AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
39 :AliGenerator(npart)
40{
41 //
42 fName="FLUKA";
43 fTitle="FLUKA Boundary Source";
44 // Read in all particle types by default
45 fIkine=6;
46 // Set maximum admitted age of particles to 1.e-05 by default
47 fAgeMax=1.e-05;
48 // Do not add weight
49 fAddWeight=1.;
50 // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
51 // from specific boundary source to the chamber (fZshift=4.5;),else there
52 // is no need to shift as it reads boundary source for the whole volume of
53 // the Muon Arm; the default value corresponds to boundary source for the
54 // whole volume of the MUON Arm
55 fZshift=0;
56 // Set the default file
57 fFileName="flukasource.root";
58
59 fTreeFluka=0;
60
61}
62
63//____________________________________________________________
64AliGenFLUKAsource::~AliGenFLUKAsource()
65{
66 delete fTreeFluka;
67}
68
69//____________________________________________________________
70void AliGenFLUKAsource::FlukaInit()
71{
72//
73// reset the existing file environment and open a new root file if
74// the pointer to the Fluka tree is null
75
76 TFile *File=0;
77 if (fTreeFluka==0) {
78 if (!File) {
79 File = new TFile(fFileName);
80 File->cd();
81 cout<<"I have opened "<<fFileName<<" file "<<endl;
82 }
83// get the tree address in the Fluka boundary source file
84 fTreeFluka = (TTree*)gDirectory->Get("h1");
85 } else {
86 File = fTreeFluka->GetCurrentFile();
87 File->cd();
88 }
89
90 TTree *h2=fTreeFluka;
91//Set branch addresses
92 h2->SetBranchAddress("Ip",&Ip);
93 h2->SetBranchAddress("Ipp",&Ipp);
94 h2->SetBranchAddress("Xi",&Xi);
95 h2->SetBranchAddress("Yi",&Yi);
96 h2->SetBranchAddress("Zi",&Zi);
97 h2->SetBranchAddress("Px",&Px);
98 h2->SetBranchAddress("Py",&Py);
99 h2->SetBranchAddress("Pz",&Pz);
100 h2->SetBranchAddress("Ekin",&Ekin);
101 h2->SetBranchAddress("Zv",&Zv);
102 h2->SetBranchAddress("Rv",&Rv);
103 h2->SetBranchAddress("Itra",&Itra);
104 h2->SetBranchAddress("Igas",&Igas);
105 h2->SetBranchAddress("Wgt",&Wgt);
106 h2->SetBranchAddress("Etag",&Etag);
107 h2->SetBranchAddress("Ptg",&Ptg);
108 h2->SetBranchAddress("Age",&Age);
109}
110
111//____________________________________________________________
112void AliGenFLUKAsource::Generate()
113{
114
115 AliMC* pMC = AliMC::GetMC();
116
117 const Int_t ifluge[28]={14, 15, 3, 2, 4, 4, 1, 13, 25, 5, 6, 10, 8, 9,
118 11, 12,18, 26, 16, 21, 19, 20, 7, 16,16,0,0,0};
119 Float_t polar[3]= {0,0,0};
120 //
121 Float_t origin[3];
122 Float_t p[3];
123 Float_t prwn;
124 Float_t wgt, fwgt;
125 Float_t phi;
126 char name[100];
127 Float_t amass, charge, tlife;
128 Int_t itrtyp;
129 Int_t iwgt;
130 Int_t i, j, part, nt;
131 static Int_t irwn;
132 //
133 Float_t random[2];
134 //
135 FlukaInit();
136 TTree *h2=fTreeFluka;
137 Int_t nentries = (Int_t) h2->GetEntries();
299e476a 138 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
fe4da5cc 139
140 // loop over number of particles
141 Int_t nb=0;
142 for (i=0; i<fNpart;i++) {
143 Int_t ev=pMC->CurrentEvent();
144 Int_t entry=fNpart*(ev-1)+i;
145 nb = (Int_t)h2->GetEvent(entry);
146 if (irwn > nentries) {
147 printf("No more entries in the FLUKA boundary source file\n");
299e476a 148 TFile *File=0;
149 // Get AliRun object or create it
150 if (!gAlice) {
151 gAlice = (AliRun*)File->Get("gAlice");
152 if (gAlice) printf("AliRun object found on file\n");
153 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
154 }
155 TTree *fAli=gAlice->TreeK();
156 if (fAli) File =fAli->GetCurrentFile();
157 File->cd();
158 printf("Generate - I'm out \n");
fe4da5cc 159 return;
160 }
161 if (Ip > 28 || Ip < 0) {
162 irwn++;
163 continue;
164 }
165
166 if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
167 irwn++;
168 continue;
169 } else if (fIkine == 9) {
170 if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
171 irwn++;
172 continue;
173 }
174 }
175
176
177 irwn++;
178 printf("\n Particle type: %f \n \n ", Ip);
179
180 if (Ip ==7){
181 prwn=Ekin;
182 part=1;
183 } else if (Ip == 8) {
184 prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
185 part=13;
186 } else {
187 part=ifluge[int(Ip)-1];
188 pMC->Gfpart(part, name, itrtyp,
189 amass, charge, tlife);
190 prwn=sqrt(Ekin*Ekin + 2.*amass);
191 }
192 origin[0]=Xi;
193 origin[1]=Yi;
194 origin[2]=Zi;
195
196 p[0]=Px*prwn;
197 p[1]=Py*prwn;
198 p[2]=Pz*prwn;
199
200 //handle particle weight correctly
201 wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
202 iwgt=Int_t(wgt);
203 fwgt=wgt-Float_t(iwgt);
204 pMC->Rndm(random,2);
205 if (random[0] < fwgt) iwgt++;
206 if (part==1 && iwgt>100) iwgt=100;
207 Int_t nstack=0;
208 for (j=0; j<iwgt; j++) {
209 gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
210 pMC->Rndm(random,2);
211 phi=2*random[1]*TMath::Pi();
212 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
213 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
214 p[0]=pn1;
215 p[1]=pn2;
216 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
217 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
218 origin[0]=on1;
219 origin[1]=on2;
220 nstack++;
221 }
222 if (nstack == 0) continue;
223 }
224
225 TFile *File=0;
226// Get AliRun object or create it
227 if (!gAlice) {
228 gAlice = (AliRun*)File->Get("gAlice");
229 if (gAlice) printf("AliRun object found on file\n");
230 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
231 }
232 TTree *fAli=gAlice->TreeK();
233 if (fAli) File =fAli->GetCurrentFile();
234 File->cd();
235}
236
237
238
239
240
241
242
243
244