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