1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
24 #include "AliGenFLUKAsource.h"
27 #include <TDatabasePDG.h>
33 #include <TVirtualMC.h>
37 ClassImp(AliGenFLUKAsource)
39 AliGenFLUKAsource::AliGenFLUKAsource()
70 fTitle="FLUKA Boundary Source";
71 // Read in all particle types by default
72 // Set maximum admitted age of particles to 1.e-05 by default
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
79 // Set the default file
80 fTreeChain = new TChain("h1");
86 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
95 fTreeChain(new TChain("h1")),
117 fTitle = "FLUKA Boundary Source";
120 AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource):
121 AliGenerator(FLUKAsource),
150 FLUKAsource.Copy(*this);
154 //____________________________________________________________
155 AliGenFLUKAsource::~AliGenFLUKAsource()
158 if (fTreeFluka) delete fTreeFluka;
159 if (fTreeChain) delete fTreeChain;
162 void AliGenFLUKAsource::AddFile(const Text_t *filename)
164 // Add a file to the chain
165 fTreeChain->Add(filename);
170 //____________________________________________________________
171 void AliGenFLUKAsource::FlukaInit()
173 // Set branch addresses of data entries
174 TChain *h2=fTreeChain;
175 //Set branch addresses
176 h2->SetBranchAddress("Ip",&fIp);
177 h2->SetBranchAddress("Ipp",&fIpp);
178 h2->SetBranchAddress("Xi",&fXi);
179 h2->SetBranchAddress("Yi",&fYi);
180 h2->SetBranchAddress("Zi",&fZi);
181 h2->SetBranchAddress("Px",&fPx);
182 h2->SetBranchAddress("Py",&fPy);
183 h2->SetBranchAddress("Pz",&fPz);
184 h2->SetBranchAddress("Ekin",&fEkin);
185 h2->SetBranchAddress("Zv",&fZv);
186 h2->SetBranchAddress("Rv",&fRv);
187 h2->SetBranchAddress("Itra",&fItra);
188 h2->SetBranchAddress("Igas",&fIgas);
189 h2->SetBranchAddress("Wgt",&fWgt);
190 h2->SetBranchAddress("Etag",&fEtag);
191 h2->SetBranchAddress("Ptg",&fPtg);
192 h2->SetBranchAddress("Age",&fAge);
195 //____________________________________________________________
196 void AliGenFLUKAsource::Generate()
198 // Generate one event
200 const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
201 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
202 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
203 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
204 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
206 Float_t polar[3]= {0,0,0};
215 Int_t i, j, part, nt;
221 TChain *h2=fTreeChain;
222 Int_t nentries = (Int_t) h2->GetEntries();
223 if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
226 // loop over number of particles
228 Int_t ev=gMC->CurrentEvent();
229 for (i=0; i<fNpart;i++) {
230 Int_t entry=fNpart*(ev-1)+i;
231 nb = (Int_t)h2->GetEvent(entry);
232 if (irwn > nentries) {
233 printf("No more entries in the FLUKA boundary source file\n");
235 // Get AliRun object or create it
237 gAlice = (AliRun*)pFile->Get("gAlice");
238 if (gAlice) printf("AliRun object found on file\n");
239 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
241 TTree *fAli=gAlice->TreeK();
242 if (fAli) pFile =fAli->GetCurrentFile();
244 printf("Generate - I'm out \n");
248 Int_t ifip = Int_t(fIp);
251 if (fSourceId != -1 && fIgas !=fSourceId) {
256 if (ifip > 28 || ifip < 0) {
261 if ((ifip != fIkine && fIkine != kAll && fIkine != kCharged
262 && fIkine != 10) || fAge > fAgeMax){
265 } else if (fIkine == kCharged) {
266 if (ifip == 7 || ifip == 8 || fAge > fAgeMax) {
270 } else if (fIkine == kNoNeutron) {
271 if (ifip == 8 || fAge > fAgeMax) {
280 // PDG code from FLUKA particle type (ifip)
281 part=kIfluge[int(ifip)-1];
283 // Calculate momentum from kinetic energy and mass of the particle
284 #if ROOT_VERSION_CODE > 197895
285 amass = gMC->ParticleMass(part);
287 amass = (TDatabasePDG::Instance())->GetParticle(part)->Mass();
289 prwn=fEkin*sqrt(1. + 2.*amass/fEkin);
300 //handle particle weight correctly
301 wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
303 fwgt=wgt-Float_t(iwgt);
305 if (random[0] < fwgt) iwgt++;
306 if (part==1 && iwgt>100) iwgt=100;
308 for (j=0; j<iwgt; j++) {
309 PushTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
311 phi=2*random[1]*TMath::Pi();
312 Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
313 Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
316 Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
317 Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
322 if (nstack == 0) continue;
326 // Get AliRun object or create it
328 gAlice = (AliRun*)pFile->Get("gAlice");
329 if (gAlice) printf("AliRun object found on file\n");
330 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
332 TTree *fAli=gAlice->TreeK();
333 if (fAli) pFile =fAli->GetCurrentFile();
338 AliGenFLUKAsource& AliGenFLUKAsource::operator=(const AliGenFLUKAsource& rhs)
340 // Assignment operator
346 void AliGenFLUKAsource::Copy(TObject &) const
348 Fatal("Copy","Not implemented!\n");