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