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