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