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