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

Diff of /PamVMC/src/PamVMCPrimaryGenerator.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by nikolas, Thu Jun 28 07:16:57 2007 UTC revision 1.2 by nikolas, Wed Feb 18 17:38:17 2009 UTC
# Line 1  Line 1 
1  // $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03  // $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03
2  //  
3    
4  #include <TVirtualMC.h>  #include <TVirtualMC.h>
5  #include <TVirtualMCStack.h>  #include <TVirtualMCStack.h>
# Line 11  Line 11 
11    
12  #include "PamVMCPrimaryGenerator.h"  #include "PamVMCPrimaryGenerator.h"
13    
14    using TMath::Sqrt;
15    using TMath::Sin;
16    using TMath::Cos;
17    using TMath::ATan;
18    using TMath::Log;
19    using TMath::Power;
20    using TMath::Exp;
21    
22    ClassImp(PamVMCPrimary)
23    
24    PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
25    {
26      a.fPDG=b.fPDG;
27      a.fX0=b.fX0;
28      a.fY0=b.fY0;
29      a.fZ0=b.fZ0;
30      a.fTHETA=b.fTHETA;
31      a.fPHI=b.fPHI;
32      a.fP0=b.fP0;
33      a.fGOOD=b.fGOOD;
34    
35      return a;
36    }
37    
38    
39  ClassImp(PamVMCPrimaryGenerator)  ClassImp(PamVMCPrimaryGenerator)
40    
41  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)
42    : TObject(),    : TObject(),
43      fStack(stack),      fStack(stack),
44      fPdg(kProton),      fevno(0),
45      fKinEnergy(100.0),//100 GeV      fmass(0.),
46      fDirX(0.),      fcharge(0.),
47      fDirY(0.),      frnd(0)
     fDirZ(-1.),  
     fPolAngle(0.),  
     fNofPrimaries(1)  
48  {  {
49  // Standard constructor  // Standard constructor
50      fprimColl = new TClonesArray("PamVMCPrimary");
51      frnd = new TRandom3(0);
52    
53  }    fprim.fPDG=kProton;
54      fprim.fX0=1.;
55      fprim.fY0=1.;
56      fprim.fZ0=130.;
57      fprim.fTHETA=0.;
58      fprim.fPHI=0.;
59      fprim.fP0=1.; //1GV
60    
61    }
62    
63  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()  PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
64    : TObject(),    : TObject(),
65      fStack(0),      fStack(0),
66      fPdg(0),      fevno(0),
67      fKinEnergy(0.),      fmass(0.),
68      fDirX(0.),      fcharge(0.),
69      fDirY(0.),      fprimColl(0),
70      fDirZ(0.),      frnd(0)
     fPolAngle(0.),  
     fNofPrimaries(0)  
71  {      {    
72      frnd = new TRandom3(0);
73  // Default constructor  // Default constructor
74      //Default primary proton
75      fprim.fPDG=kProton;
76      fprim.fX0=1.;
77      fprim.fY0=1.;
78      fprim.fZ0=130.;
79      fprim.fTHETA=0.;
80      fprim.fPHI=0.;
81      fprim.fP0=1.; //1GV
82  }  }
83    
84  PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()  PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
85  {  {
86  // Destructor    // Destructor
87      delete frnd;
88      delete fprimColl;
89  }  }
90    
91  // private methods  // private methods
# Line 62  void PamVMCPrimaryGenerator::GeneratePri Line 103  void PamVMCPrimaryGenerator::GeneratePri
103    Int_t toBeDone = 1;    Int_t toBeDone = 1;
104    
105    // Particle type    // Particle type
106    Int_t pdg  = fPdg;    Int_t pdg  = fprim.fPDG;
107    TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fPdg);    
108      Double_t fvx, fvy, fvz;
109      fvx=fprim.fX0;
110      fvy=fprim.fY0;
111      fvz=fprim.fZ0;
112    
113    // Position    // Position
114    Double_t vx  = 0.;    
   Double_t vy  = 0.;  
   Double_t vz =  110.;  
115    Double_t tof = 0.;    Double_t tof = 0.;
116    
117    // Energy (in GeV)    // Energy (in GeV)
118    Double_t kinEnergy = fKinEnergy;    Double_t kinEnergy = MomentumToKinE(fprim.fP0);    
119    Double_t mass = particlePDG->Mass();          Double_t e  = fmass + kinEnergy;
   Double_t e  = mass + kinEnergy;  
120    
121    // Particle momentum    // Particle momentum
122    Double_t p0, px, py, pz;    Double_t  px, py, pz;
123    p0 = sqrt(e*e - mass*mass);    
124    px = p0 * fDirX;    px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI);
125    py = p0 * fDirY;    py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
126    pz = p0 * fDirZ;    pz = -fprim.fP0*Cos(fprim.fTHETA);
127    
128    // Polarization    // Polarization
129    TVector3 polar;    TVector3 polar;
 //  if ( fPdg == 50000050 ) {  
 //    TVector3 normal (1., 0., 0.);  
 //    TVector3 kphoton = TVector3(fDirX, fDirY, fDirZ);  
 //    TVector3 product = normal.Cross(kphoton);  
 //    Double_t modul2  = product*product;  
   
 //    TVector3 e_perpend (0., 0., 1.);  
 //    if (modul2 > 0.) e_perpend = (1./sqrt(modul2))*product;  
 //    TVector3 e_paralle = e_perpend.Cross(kphoton);  
   
 ///    polar =   TMath::Cos(fPolAngle*TMath::DegToRad())*e_paralle  
   //           + TMath::Sin(fPolAngle*TMath::DegToRad())*e_perpend;  
   //  }  
   //else  
   //  Warning("GeneratePrimary",  
   //          "The primary particle is not an opticalphoton");  
130    
131    // Add particle to stack    // Add particle to stack
132    fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, vx, vy, vz, tof,    fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof,
133                     polar.X(), polar.Y(), polar.Z(),                     polar.X(), polar.Y(), polar.Z(),
134                     kPPrimary, ntr, 1., 0);                     kPPrimary, ntr, 1., 0);
135    
136      PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
137    
138      *pc = fprim;
139  }  }
140    
 // public methods  
141    
142  void PamVMCPrimaryGenerator::GeneratePrimaries()  void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
143  {        fprim.fPDG=pdg;
144  // Fill the user stack (derived from TVirtualMCStack) with primary particle    //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
145      fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
146      fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
147    }
148    
149    for (Int_t i=0; i<fNofPrimaries; i++) GeneratePrimary();    void PamVMCPrimaryGenerator::SetMomentum(
150                                  Double_t px, Double_t py, Double_t pz)
151    {
152      fprim.fP0= Sqrt(px*px+py*py+pz*pz);
153      fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz);
154      fprim.fPHI=ATan(py/px);
155  }  }
156                                        
157    void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
158    {
159      if(isEnergy) {
160        fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
161      } else{
162        fprim.fP0=frnd->Uniform(PEmin,PEmax);
163      }
164    
165  void PamVMCPrimaryGenerator::SetDirection(  }
                               Double_t dirX, Double_t dirY, Double_t dirZ)  
 {  
 // Set normalized direction  
166    
167    Double_t norm = TMath::Sqrt(dirX*dirX + dirY*dirY + dirZ*dirZ);  void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
168      {
169    fDirX = dirX/norm;    Double_t alpha = 1.+gamma; //integral spectral index
170    fDirY = dirY/norm;      if(alpha==0.){
171    fDirZ = dirZ/norm;      fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
172  }                                        } else {
173        if(PEmin==0.) PEmin=1.E-10;
174        fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
175      }
176    
177      if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
178    
179    }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23