Class moved from AliGenCocktail.h
[u/mrichter/AliRoot.git] / EVGEN / AliGenFLUKAsource.cxx
CommitLineData
4c039060 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$Log$
e7bafa6d 18Revision 1.8 2000/02/14 14:49:38 morsch
19Correct particle type for gamma and neutrons
20More consistent calculation of momentum from kin. energy and mass
21
02769dba 22Revision 1.7 1999/11/03 17:43:20 fca
23New version from G.Martinez & A.Morsch
24
886b6f73 25Revision 1.6 1999/09/29 09:24:12 fca
26Introduction of the Copyright and cvs Log
27
4c039060 28*/
29
fe4da5cc 30#include "AliGenFLUKAsource.h"
31#include "AliGenMUONlib.h"
32#include "AliMC.h"
33#include "AliRun.h"
c5970876 34#include "AliPDG.h"
fe4da5cc 35#include <TDirectory.h>
36#include <TFile.h>
37#include <TTree.h>
38#include <stdlib.h>
39 ClassImp(AliGenFLUKAsource)
40 AliGenFLUKAsource::AliGenFLUKAsource()
41 :AliGenerator(-1)
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;
886b6f73 62 fTreeChain = new TChain("h1");
fe4da5cc 63//
64// Read all particles
65 fNpart=-1;
66
67}
68
69AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
70 :AliGenerator(npart)
71{
72 //
73 fName="FLUKA";
74 fTitle="FLUKA Boundary Source";
75 // Read in all particle types by default
76 fIkine=6;
77 // Set maximum admitted age of particles to 1.e-05 by default
78 fAgeMax=1.e-05;
79 // Do not add weight
80 fAddWeight=1.;
81 // Shift the z-coordinate of the impact point by 4.5 cm only if it reads
82 // from specific boundary source to the chamber (fZshift=4.5;),else there
83 // is no need to shift as it reads boundary source for the whole volume of
84 // the Muon Arm; the default value corresponds to boundary source for the
85 // whole volume of the MUON Arm
86 fZshift=0;
87 // Set the default file
88 fFileName="flukasource.root";
89
90 fTreeFluka=0;
886b6f73 91 fTreeChain = new TChain("h1");
fe4da5cc 92}
93
94//____________________________________________________________
95AliGenFLUKAsource::~AliGenFLUKAsource()
96{
886b6f73 97 if (fTreeFluka) delete fTreeFluka;
98 if (fTreeChain) delete fTreeChain;
99 // if (fFileName) delete fFileName;
fe4da5cc 100}
101
886b6f73 102void AliGenFLUKAsource::AddFile(const Text_t *filename)
fe4da5cc 103{
886b6f73 104 fTreeChain->Add(filename);
fe4da5cc 105
886b6f73 106}
107
fe4da5cc 108
886b6f73 109//____________________________________________________________
110void AliGenFLUKAsource::FlukaInit()
111{
112 TChain *h2=fTreeChain;
fe4da5cc 113//Set branch addresses
114 h2->SetBranchAddress("Ip",&Ip);
115 h2->SetBranchAddress("Ipp",&Ipp);
116 h2->SetBranchAddress("Xi",&Xi);
117 h2->SetBranchAddress("Yi",&Yi);
118 h2->SetBranchAddress("Zi",&Zi);
119 h2->SetBranchAddress("Px",&Px);
120 h2->SetBranchAddress("Py",&Py);
121 h2->SetBranchAddress("Pz",&Pz);
122 h2->SetBranchAddress("Ekin",&Ekin);
123 h2->SetBranchAddress("Zv",&Zv);
124 h2->SetBranchAddress("Rv",&Rv);
125 h2->SetBranchAddress("Itra",&Itra);
126 h2->SetBranchAddress("Igas",&Igas);
127 h2->SetBranchAddress("Wgt",&Wgt);
128 h2->SetBranchAddress("Etag",&Etag);
129 h2->SetBranchAddress("Ptg",&Ptg);
130 h2->SetBranchAddress("Age",&Age);
131}
132
133//____________________________________________________________
134void AliGenFLUKAsource::Generate()
886b6f73 135{
136 AliMC* gMC = AliMC::GetMC();
fe4da5cc 137
886b6f73 138 const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
1578254f 139 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
140 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
141 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
142 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
143 0,kNuMu,kNuMuBar};
886b6f73 144 Float_t polar[3]= {0,0,0};
fe4da5cc 145 //
886b6f73 146 Float_t origin[3];
147 Float_t p[3];
148 Float_t prwn;
149 Float_t wgt, fwgt;
150 Float_t phi;
151 char name[100];
152 Float_t amass, charge, tlife;
153 Int_t itrtyp;
154 Int_t iwgt;
155 Int_t i, j, part, nt;
e7bafa6d 156 static Int_t irwn=0;
886b6f73 157 //
158 Float_t random[2];
fe4da5cc 159 //
886b6f73 160 FlukaInit();
161 TChain *h2=fTreeChain;
162 Int_t nentries = (Int_t) h2->GetEntries();
163 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
fe4da5cc 164
165 // loop over number of particles
886b6f73 166 Int_t nb=0;
cfce8870 167 Int_t ev=gMC->CurrentEvent();
886b6f73 168 for (i=0; i<fNpart;i++) {
169 Int_t entry=fNpart*(ev-1)+i;
170 nb = (Int_t)h2->GetEvent(entry);
171 if (irwn > nentries) {
172 printf("No more entries in the FLUKA boundary source file\n");
173 TFile *File=0;
174 // Get AliRun object or create it
175 if (!gAlice) {
176 gAlice = (AliRun*)File->Get("gAlice");
177 if (gAlice) printf("AliRun object found on file\n");
178 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
179 }
180 TTree *fAli=gAlice->TreeK();
181 if (fAli) File =fAli->GetCurrentFile();
182 File->cd();
183 printf("Generate - I'm out \n");
184 return;
185 }
186 if (Ip > 28 || Ip < 0) {
187 irwn++;
188 continue;
189 }
190
191 if ((Ip != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) || Age > fAgeMax){
192 irwn++;
193 continue;
194 } else if (fIkine == 9) {
195 if (Ip == 7 || Ip == 8 || Age > fAgeMax) {
196 irwn++;
197 continue;
198 }
199 } else if (fIkine == 10) {
200 if (Ip == 8 || Age > fAgeMax) {
201 irwn++;
202 continue;
203 }
204 }
fe4da5cc 205
fe4da5cc 206
886b6f73 207 irwn++;
02769dba 208//
209// PDG code from FLUKA particle type (Ip)
210 part=ifluge[int(Ip)-1];
211//
212// Calculate momentum from kinetic energy and mass of the particle
213 gMC->Gfpart(part, name, itrtyp,
214 amass, charge, tlife);
215 prwn=Ekin*sqrt(1. + 2.*amass/Ekin);
216
217
886b6f73 218 origin[0]=Yi;
219 origin[1]=Xi;
220 origin[2]=Zi;
221
222 p[0]=Py*prwn;
223 p[1]=Px*prwn;
224 p[2]=Pz*prwn;
fe4da5cc 225
886b6f73 226 //handle particle weight correctly
227 wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
228 iwgt=Int_t(wgt);
229 fwgt=wgt-Float_t(iwgt);
cfce8870 230 gMC->Rndm(random,2);
886b6f73 231 if (random[0] < fwgt) iwgt++;
232 if (part==1 && iwgt>100) iwgt=100;
233 Int_t nstack=0;
234 for (j=0; j<iwgt; j++) {
235 gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,Age,"Primary",nt);
236 gMC->Rndm(random,2);
237 phi=2*random[1]*TMath::Pi();
238 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
239 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
240 p[0]=pn1;
241 p[1]=pn2;
242 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
243 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
244 origin[0]=on1;
245 origin[1]=on2;
246 nstack++;
247 }
248 if (nstack == 0) continue;
fe4da5cc 249 }
886b6f73 250
fe4da5cc 251 TFile *File=0;
252// Get AliRun object or create it
253 if (!gAlice) {
254 gAlice = (AliRun*)File->Get("gAlice");
255 if (gAlice) printf("AliRun object found on file\n");
256 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
257 }
258 TTree *fAli=gAlice->TreeK();
259 if (fAli) File =fAli->GetCurrentFile();
260 File->cd();
261}
262
263
264
265
266
267
268
269
270