Coverity fix.
[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
88cb7938 16/* $Id$ */
675e9664 17
18// Read background particles from a FLUKA boundary source file
19// This is a very special generator that works for background studies for the muon-spectrometer.
20// The input files come from FLUKA simulations.
21// Files can be chained.
22// Author: andreas.morsch@cern.ch
23
2524c56f 24#include "AliGenFLUKAsource.h"
ac3faee4 25#include <stdlib.h>
26
27#include <TDatabasePDG.h>
28#include <TPDGCode.h>
9e77b997 29#include <RVersion.h>
ac3faee4 30#include <TChain.h>
31#include <TFile.h>
32#include <TTree.h>
5d12ce38 33#include <TVirtualMC.h>
116cbefd 34
fe4da5cc 35#include "AliRun.h"
5c3fd7ea 36
198bb1c7 37ClassImp(AliGenFLUKAsource)
38
39AliGenFLUKAsource::AliGenFLUKAsource()
1c56e311 40 :AliGenerator(-1),
41 fIkine(6),
42 fAgeMax(1.e-5),
43 fAddWeight(1.),
44 fZshift(0.),
45 fFrac(0.),
46 fSourceId(-1),
47 fFileName(0),
48 fTreeChain(0),
49 fTreeFluka(0),
50 fIp(0.),
51 fIpp(0.),
52 fXi(0.),
53 fYi(0.),
54 fZi(0.),
55 fPx(0.),
56 fPy(0.),
57 fPz(0.),
58 fEkin(0.),
59 fZv(0.),
60 fRv(0.),
61 fItra(0.),
62 fIgas(0.),
63 fWgt(0.),
64 fEtag(0.),
65 fPtg(0.),
66 fAge(0.)
fe4da5cc 67{
f87cfe57 68 // Constructor
fe4da5cc 69 fName="FLUKA";
70 fTitle="FLUKA Boundary Source";
71 // Read in all particle types by default
fe4da5cc 72 // Set maximum admitted age of particles to 1.e-05 by default
fe4da5cc 73 // Do not add weight
fe4da5cc 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
fe4da5cc 79 // Set the default file
886b6f73 80 fTreeChain = new TChain("h1");
fe4da5cc 81//
82// Read all particles
83 fNpart=-1;
fe4da5cc 84}
85
86AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
1c56e311 87 :AliGenerator(npart),
88 fIkine(6),
89 fAgeMax(1.e-5),
90 fAddWeight(1.),
91 fZshift(0.),
92 fFrac(0.),
93 fSourceId(-1),
94 fFileName(""),
95 fTreeChain(new TChain("h1")),
96 fTreeFluka(0),
97 fIp(0.),
98 fIpp(0.),
99 fXi(0.),
100 fYi(0.),
101 fZi(0.),
102 fPx(0.),
103 fPy(0.),
104 fPz(0.),
105 fEkin(0.),
106 fZv(0.),
107 fRv(0.),
108 fItra(0.),
109 fIgas(0.),
110 fWgt(0.),
111 fEtag(0.),
112 fPtg(0.),
113 fAge(0.)
fe4da5cc 114{
f87cfe57 115 // Constructor
675e9664 116 fName = "FLUKA";
117 fTitle = "FLUKA Boundary Source";
fe4da5cc 118}
119
120//____________________________________________________________
121AliGenFLUKAsource::~AliGenFLUKAsource()
122{
f87cfe57 123// Destructor
886b6f73 124 if (fTreeFluka) delete fTreeFluka;
125 if (fTreeChain) delete fTreeChain;
fe4da5cc 126}
127
886b6f73 128void AliGenFLUKAsource::AddFile(const Text_t *filename)
fe4da5cc 129{
f87cfe57 130// Add a file to the chain
886b6f73 131 fTreeChain->Add(filename);
fe4da5cc 132
886b6f73 133}
134
fe4da5cc 135
886b6f73 136//____________________________________________________________
137void AliGenFLUKAsource::FlukaInit()
138{
f87cfe57 139// Set branch addresses of data entries
886b6f73 140 TChain *h2=fTreeChain;
fe4da5cc 141//Set branch addresses
f87cfe57 142 h2->SetBranchAddress("Ip",&fIp);
143 h2->SetBranchAddress("Ipp",&fIpp);
144 h2->SetBranchAddress("Xi",&fXi);
145 h2->SetBranchAddress("Yi",&fYi);
146 h2->SetBranchAddress("Zi",&fZi);
147 h2->SetBranchAddress("Px",&fPx);
148 h2->SetBranchAddress("Py",&fPy);
149 h2->SetBranchAddress("Pz",&fPz);
150 h2->SetBranchAddress("Ekin",&fEkin);
151 h2->SetBranchAddress("Zv",&fZv);
152 h2->SetBranchAddress("Rv",&fRv);
153 h2->SetBranchAddress("Itra",&fItra);
154 h2->SetBranchAddress("Igas",&fIgas);
155 h2->SetBranchAddress("Wgt",&fWgt);
156 h2->SetBranchAddress("Etag",&fEtag);
157 h2->SetBranchAddress("Ptg",&fPtg);
158 h2->SetBranchAddress("Age",&fAge);
fe4da5cc 159}
160
161//____________________________________________________________
162void AliGenFLUKAsource::Generate()
f87cfe57 163{
164// Generate one event
fe4da5cc 165
f87cfe57 166 const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
1578254f 167 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
168 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
169 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
170 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
171 0,kNuMu,kNuMuBar};
886b6f73 172 Float_t polar[3]= {0,0,0};
fe4da5cc 173 //
886b6f73 174 Float_t origin[3];
175 Float_t p[3];
176 Float_t prwn;
177 Float_t wgt, fwgt;
178 Float_t phi;
d3ee5d3d 179 Float_t amass;
886b6f73 180 Int_t iwgt;
181 Int_t i, j, part, nt;
e7bafa6d 182 static Int_t irwn=0;
886b6f73 183 //
184 Float_t random[2];
fe4da5cc 185 //
886b6f73 186 FlukaInit();
187 TChain *h2=fTreeChain;
188 Int_t nentries = (Int_t) h2->GetEntries();
189 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
09d89034 190
191
fe4da5cc 192 // loop over number of particles
886b6f73 193 Int_t nb=0;
cfce8870 194 Int_t ev=gMC->CurrentEvent();
886b6f73 195 for (i=0; i<fNpart;i++) {
196 Int_t entry=fNpart*(ev-1)+i;
197 nb = (Int_t)h2->GetEvent(entry);
198 if (irwn > nentries) {
199 printf("No more entries in the FLUKA boundary source file\n");
f87cfe57 200 TFile *pFile=0;
886b6f73 201 // Get AliRun object or create it
202 if (!gAlice) {
f87cfe57 203 gAlice = (AliRun*)pFile->Get("gAlice");
886b6f73 204 if (gAlice) printf("AliRun object found on file\n");
205 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
206 }
33c3c91a 207 TTree *fAli=AliRunLoader::Instance()->TreeK();
f87cfe57 208 if (fAli) pFile =fAli->GetCurrentFile();
209 pFile->cd();
886b6f73 210 printf("Generate - I'm out \n");
211 return;
212 }
09d89034 213
214 Int_t ifip = Int_t(fIp);
215
f87cfe57 216
217 if (fSourceId != -1 && fIgas !=fSourceId) {
218 irwn++;
219 continue;
220 }
221
09d89034 222 if (ifip > 28 || ifip < 0) {
886b6f73 223 irwn++;
224 continue;
225 }
226
09d89034 227 if ((ifip != fIkine && fIkine != kAll && fIkine != kCharged
228 && fIkine != 10) || fAge > fAgeMax){
886b6f73 229 irwn++;
230 continue;
09d89034 231 } else if (fIkine == kCharged) {
232 if (ifip == 7 || ifip == 8 || fAge > fAgeMax) {
886b6f73 233 irwn++;
234 continue;
235 }
09d89034 236 } else if (fIkine == kNoNeutron) {
237 if (ifip == 8 || fAge > fAgeMax) {
886b6f73 238 irwn++;
239 continue;
240 }
241 }
fe4da5cc 242
fe4da5cc 243
886b6f73 244 irwn++;
02769dba 245//
09d89034 246// PDG code from FLUKA particle type (ifip)
247 part=kIfluge[int(ifip)-1];
02769dba 248//
249// Calculate momentum from kinetic energy and mass of the particle
2ec105c5 250#if ROOT_VERSION_CODE > 197895
b4d4098d 251 amass = gMC->ParticleMass(part);
2ec105c5 252#else
253 amass = (TDatabasePDG::Instance())->GetParticle(part)->Mass();
254#endif
f87cfe57 255 prwn=fEkin*sqrt(1. + 2.*amass/fEkin);
02769dba 256
257
f87cfe57 258 origin[0]=fYi;
259 origin[1]=fXi;
260 origin[2]=fZi;
886b6f73 261
f87cfe57 262 p[0]=fPy*prwn;
263 p[1]=fPx*prwn;
264 p[2]=fPz*prwn;
fe4da5cc 265
886b6f73 266 //handle particle weight correctly
f87cfe57 267 wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
886b6f73 268 iwgt=Int_t(wgt);
269 fwgt=wgt-Float_t(iwgt);
65fb704d 270 Rndm(random,2);
886b6f73 271 if (random[0] < fwgt) iwgt++;
272 if (part==1 && iwgt>100) iwgt=100;
273 Int_t nstack=0;
274 for (j=0; j<iwgt; j++) {
642f15cf 275 PushTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
65fb704d 276 Rndm(random,2);
886b6f73 277 phi=2*random[1]*TMath::Pi();
278 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
279 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
280 p[0]=pn1;
281 p[1]=pn2;
282 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
283 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
284 origin[0]=on1;
285 origin[1]=on2;
286 nstack++;
287 }
288 if (nstack == 0) continue;
fe4da5cc 289 }
886b6f73 290
f87cfe57 291 TFile *pFile=0;
fe4da5cc 292// Get AliRun object or create it
293 if (!gAlice) {
f87cfe57 294 gAlice = (AliRun*)pFile->Get("gAlice");
fe4da5cc 295 if (gAlice) printf("AliRun object found on file\n");
296 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
297 }
33c3c91a 298 TTree *fAli=AliRunLoader::Instance()->TreeK();
f87cfe57 299 if (fAli) pFile =fAli->GetCurrentFile();
300 pFile->cd();
fe4da5cc 301}
302
303
fe4da5cc 304
305
306
307