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