Use gMC and not pMC everywhere
[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
1578254f 117 const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
118 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
119 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
120 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
121 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
122 0,kNuMu,kNuMuBar};
fe4da5cc 123 Float_t polar[3]= {0,0,0};
124 //
125 Float_t origin[3];
126 Float_t p[3];
127 Float_t prwn;
128 Float_t wgt, fwgt;
129 Float_t phi;
130 char name[100];
131 Float_t amass, charge, tlife;
132 Int_t itrtyp;
133 Int_t iwgt;
134 Int_t i, j, part, nt;
135 static Int_t irwn;
136 //
137 Float_t random[2];
138 //
139 FlukaInit();
140 TTree *h2=fTreeFluka;
141 Int_t nentries = (Int_t) h2->GetEntries();
299e476a 142 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
fe4da5cc 143
144 // loop over number of particles
145 Int_t nb=0;
146 for (i=0; i<fNpart;i++) {
cfce8870 147 Int_t ev=gMC->CurrentEvent();
fe4da5cc 148 Int_t entry=fNpart*(ev-1)+i;
149 nb = (Int_t)h2->GetEvent(entry);
150 if (irwn > nentries) {
151 printf("No more entries in the FLUKA boundary source file\n");
299e476a 152 TFile *File=0;
153 // Get AliRun object or create it
154 if (!gAlice) {
155 gAlice = (AliRun*)File->Get("gAlice");
156 if (gAlice) printf("AliRun object found on file\n");
157 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
158 }
159 TTree *fAli=gAlice->TreeK();
160 if (fAli) File =fAli->GetCurrentFile();
161 File->cd();
162 printf("Generate - I'm out \n");
fe4da5cc 163 return;
164 }
165 if (Ip > 28 || Ip < 0) {
166 irwn++;
167 continue;
168 }
169
170 if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
171 irwn++;
172 continue;
173 } else if (fIkine == 9) {
174 if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
175 irwn++;
176 continue;
177 }
178 }
179
180
181 irwn++;
182 printf("\n Particle type: %f \n \n ", Ip);
183
184 if (Ip ==7){
185 prwn=Ekin;
186 part=1;
187 } else if (Ip == 8) {
188 prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
189 part=13;
190 } else {
191 part=ifluge[int(Ip)-1];
cfce8870 192 gMC->Gfpart(part, name, itrtyp,
fe4da5cc 193 amass, charge, tlife);
194 prwn=sqrt(Ekin*Ekin + 2.*amass);
195 }
196 origin[0]=Xi;
197 origin[1]=Yi;
198 origin[2]=Zi;
199
200 p[0]=Px*prwn;
201 p[1]=Py*prwn;
202 p[2]=Pz*prwn;
203
204 //handle particle weight correctly
205 wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
206 iwgt=Int_t(wgt);
207 fwgt=wgt-Float_t(iwgt);
cfce8870 208 gMC->Rndm(random,2);
fe4da5cc 209 if (random[0] < fwgt) iwgt++;
210 if (part==1 && iwgt>100) iwgt=100;
211 Int_t nstack=0;
212 for (j=0; j<iwgt; j++) {
213 gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
cfce8870 214 gMC->Rndm(random,2);
fe4da5cc 215 phi=2*random[1]*TMath::Pi();
216 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
217 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
218 p[0]=pn1;
219 p[1]=pn2;
220 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
221 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
222 origin[0]=on1;
223 origin[1]=on2;
224 nstack++;
225 }
226 if (nstack == 0) continue;
227 }
228
229 TFile *File=0;
230// Get AliRun object or create it
231 if (!gAlice) {
232 gAlice = (AliRun*)File->Get("gAlice");
233 if (gAlice) printf("AliRun object found on file\n");
234 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
235 }
236 TTree *fAli=gAlice->TreeK();
237 if (fAli) File =fAli->GetCurrentFile();
238 File->cd();
239}
240
241
242
243
244
245
246
247
248