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