/[PAMELA software]/PamVMC/src/PamVMCPrimaryGenerator.cxx
ViewVC logotype

Annotation of /PamVMC/src/PamVMCPrimaryGenerator.cxx

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.5 - (hide annotations) (download)
Fri Jun 12 18:39:32 2009 UTC (15 years, 5 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v1r0
Changes since 1.2: +17 -19 lines 
- Introduced user-defined names of output files and random seeds number.
Users can do it use options of PamVMCApplication constructor:
PamVMCApplication(const char* name,  const char *title, const char*
filename="pamtest", Int_t seed=0).
The Random object that I use is TRandom3 object which has astronomical
large period (in case of default initialization 0). All random generators
in the code use this object by calling of gRandom singleton which keeps
it.

- Corrected TOF digitization routine. No problems with TDC hits due to
hadronic interactions anymore.

- Some small changes was done to compile code under Root 5.23. +
geant4_vmc v. 2.6 without any warnings

- Some classes of PamG4RunConfiguartion was changed for geant4_vmc v.
2.6.Some obsolete classes was deleted as soon as developers implemented
regions.

- Navigation was changed from "geomRootToGeant4" to "geomRoot", because on
VMC web page written that as soon as Geant4 has no option ONLY/MANY
translation of overlapped geometry to Geant4 through VGM could be wrong.
I'd like to stay with Root navigation:
http://root.cern.ch/root/vmc/Geant4VMC.html. This should be default
option.

- New Tracker digitization routine written by Sergio was implemented

- PamVMC again became compatible with geant4_vmc v.2.5 and ROOT 5.20.
 The problem was that ROOT developers introduced in TVirtualMC class a new
method SetMagField and new base class:TVirtualMagField from which
user-defined classes shoukd be derived
1 nikolas 1.1 // $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03
2 nikolas 1.2
3 nikolas 1.1
4     #include <TVirtualMC.h>
5     #include <TVirtualMCStack.h>
6     #include <TPDGCode.h>
7     #include <TDatabasePDG.h>
8     #include <TParticlePDG.h>
9     #include <TVector3.h>
10     #include <TMath.h>
11 pam-rm2 1.5 #include <Riostream.h>
12 nikolas 1.1
13     #include "PamVMCPrimaryGenerator.h"
14    
15 pam-rm2 1.5 using namespace TMath;
16 nikolas 1.2
17     ClassImp(PamVMCPrimary)
18    
19     PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
20     {
21     a.fPDG=b.fPDG;
22     a.fX0=b.fX0;
23     a.fY0=b.fY0;
24     a.fZ0=b.fZ0;
25     a.fTHETA=b.fTHETA;
26     a.fPHI=b.fPHI;
27     a.fP0=b.fP0;
28     a.fGOOD=b.fGOOD;
29    
30     return a;
31     }
32    
33    
34 nikolas 1.1 ClassImp(PamVMCPrimaryGenerator)
35    
36     PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)
37     : TObject(),
38     fStack(stack),
39 nikolas 1.2 fevno(0),
40     fmass(0.),
41     fcharge(0.),
42 pam-rm2 1.5 frandom(0)
43 nikolas 1.1 {
44     // Standard constructor
45 pam-rm2 1.5
46     ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
47     ftheta->SetNpx(1000);
48    
49 nikolas 1.2 fprimColl = new TClonesArray("PamVMCPrimary");
50     fprim.fPDG=kProton;
51     fprim.fX0=1.;
52     fprim.fY0=1.;
53     fprim.fZ0=130.;
54     fprim.fTHETA=0.;
55     fprim.fPHI=0.;
56     fprim.fP0=1.; //1GV
57 nikolas 1.1
58 nikolas 1.2 }
59 nikolas 1.1
60     PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
61     : TObject(),
62     fStack(0),
63 nikolas 1.2 fevno(0),
64     fmass(0.),
65     fcharge(0.),
66     fprimColl(0),
67 pam-rm2 1.5 frandom(0)
68 nikolas 1.1 {
69 pam-rm2 1.5 // Default constructor
70 nikolas 1.2 //Default primary proton
71 pam-rm2 1.5 ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
72     ftheta->SetNpx(1000);
73    
74 nikolas 1.2 fprim.fPDG=kProton;
75     fprim.fX0=1.;
76     fprim.fY0=1.;
77     fprim.fZ0=130.;
78     fprim.fTHETA=0.;
79     fprim.fPHI=0.;
80     fprim.fP0=1.; //1GV
81 nikolas 1.1 }
82    
83     PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
84     {
85 nikolas 1.2 // Destructor
86 pam-rm2 1.5 delete ftheta;
87 nikolas 1.2 delete fprimColl;
88 nikolas 1.1 }
89    
90     // private methods
91    
92    
93     void PamVMCPrimaryGenerator::GeneratePrimary()
94     {
95     // Add one primary particle to the user stack (derived from TVirtualMCStack).
96    
97     // Track ID (filled by stack)
98     Int_t ntr;
99    
100     // Option: to be tracked
101     Int_t toBeDone = 1;
102    
103     // Particle type
104 nikolas 1.2 Int_t pdg = fprim.fPDG;
105    
106     Double_t fvx, fvy, fvz;
107     fvx=fprim.fX0;
108     fvy=fprim.fY0;
109     fvz=fprim.fZ0;
110    
111 nikolas 1.1 // Position
112 nikolas 1.2
113 nikolas 1.1 Double_t tof = 0.;
114    
115     // Energy (in GeV)
116 nikolas 1.2 Double_t kinEnergy = MomentumToKinE(fprim.fP0);
117     Double_t e = fmass + kinEnergy;
118 nikolas 1.1
119     // Particle momentum
120 nikolas 1.2 Double_t px, py, pz;
121    
122     px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI);
123     py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
124     pz = -fprim.fP0*Cos(fprim.fTHETA);
125 nikolas 1.1
126     // Polarization
127     TVector3 polar;
128    
129     // Add particle to stack
130 nikolas 1.2 fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof,
131 nikolas 1.1 polar.X(), polar.Y(), polar.Z(),
132     kPPrimary, ntr, 1., 0);
133 nikolas 1.2
134     PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
135    
136     *pc = fprim;
137 nikolas 1.1 }
138    
139    
140 nikolas 1.2 void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
141     fprim.fPDG=pdg;
142     //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
143     fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
144     fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
145     }
146 nikolas 1.1
147 nikolas 1.2 void PamVMCPrimaryGenerator::SetMomentum(
148     Double_t px, Double_t py, Double_t pz)
149     {
150     fprim.fP0= Sqrt(px*px+py*py+pz*pz);
151     fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz);
152     fprim.fPHI=ATan(py/px);
153 nikolas 1.1 }
154 nikolas 1.2
155     void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
156     {
157     if(isEnergy) {
158 pam-rm2 1.5 fprim.fP0=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
159 nikolas 1.2 } else{
160 pam-rm2 1.5 fprim.fP0=frandom->Uniform(PEmin,PEmax);
161 nikolas 1.2 }
162 nikolas 1.1
163 nikolas 1.2 }
164 nikolas 1.1
165 nikolas 1.2 void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
166     {
167     Double_t alpha = 1.+gamma; //integral spectral index
168     if(alpha==0.){
169 pam-rm2 1.5 fprim.fP0=Exp(Log(PEmin)+frandom->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
170 nikolas 1.2 } else {
171     if(PEmin==0.) PEmin=1.E-10;
172 pam-rm2 1.5 fprim.fP0=Power((frandom->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
173 nikolas 1.2 }
174 nikolas 1.1
175 nikolas 1.2 if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
176 nikolas 1.1
177 nikolas 1.2 }
  ViewVC Help
Powered by ViewVC 1.1.23