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