/[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.2 by nikolas, Wed Feb 18 17:38:17 2009 UTC revision 1.5 by pam-rm2, Fri Jun 12 18:39:32 2009 UTC
# Line 8  Line 8 
8  #include <TParticlePDG.h>  #include <TParticlePDG.h>
9  #include <TVector3.h>  #include <TVector3.h>
10  #include <TMath.h>  #include <TMath.h>
11    #include <Riostream.h>
12    
13  #include "PamVMCPrimaryGenerator.h"  #include "PamVMCPrimaryGenerator.h"
14    
15  using TMath::Sqrt;  using namespace TMath;
 using TMath::Sin;  
 using TMath::Cos;  
 using TMath::ATan;  
 using TMath::Log;  
 using TMath::Power;  
 using TMath::Exp;  
16    
17  ClassImp(PamVMCPrimary)  ClassImp(PamVMCPrimary)
18    
# Line 44  PamVMCPrimaryGenerator::PamVMCPrimaryGen Line 39  PamVMCPrimaryGenerator::PamVMCPrimaryGen
39      fevno(0),      fevno(0),
40      fmass(0.),      fmass(0.),
41      fcharge(0.),      fcharge(0.),
42      frnd(0)      frandom(0)
43  {  {
44  // Standard constructor  // Standard constructor
   fprimColl = new TClonesArray("PamVMCPrimary");  
   frnd = new TRandom3(0);  
45    
46      ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
47      ftheta->SetNpx(1000);
48    
49      fprimColl = new TClonesArray("PamVMCPrimary");
50    fprim.fPDG=kProton;    fprim.fPDG=kProton;
51    fprim.fX0=1.;    fprim.fX0=1.;
52    fprim.fY0=1.;    fprim.fY0=1.;
# Line 67  PamVMCPrimaryGenerator::PamVMCPrimaryGen Line 64  PamVMCPrimaryGenerator::PamVMCPrimaryGen
64      fmass(0.),      fmass(0.),
65      fcharge(0.),      fcharge(0.),
66      fprimColl(0),      fprimColl(0),
67      frnd(0)      frandom(0)
68  {      {    
69    frnd = new TRandom3(0);   // Default constructor
 // Default constructor  
70    //Default primary proton    //Default primary proton
71      ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.);
72      ftheta->SetNpx(1000);
73    
74    fprim.fPDG=kProton;    fprim.fPDG=kProton;
75    fprim.fX0=1.;    fprim.fX0=1.;
76    fprim.fY0=1.;    fprim.fY0=1.;
# Line 84  PamVMCPrimaryGenerator::PamVMCPrimaryGen Line 83  PamVMCPrimaryGenerator::PamVMCPrimaryGen
83  PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()  PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
84  {  {
85  // Destructor  // Destructor
86    delete frnd;    delete ftheta;
87    delete fprimColl;    delete fprimColl;
88  }  }
89    
90  // private methods  // private methods
91    
 #include <Riostream.h>  
92    
93  void PamVMCPrimaryGenerator::GeneratePrimary()  void PamVMCPrimaryGenerator::GeneratePrimary()
94  {      {    
# Line 157  void PamVMCPrimaryGenerator::SetMomentum Line 155  void PamVMCPrimaryGenerator::SetMomentum
155  void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)  void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
156  {  {
157    if(isEnergy) {    if(isEnergy) {
158      fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));      fprim.fP0=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
159    } else{    } else{
160      fprim.fP0=frnd->Uniform(PEmin,PEmax);      fprim.fP0=frandom->Uniform(PEmin,PEmax);
161    }    }
162    
163  }  }
# Line 168  void PamVMCPrimaryGenerator::GenSpe(Doub Line 166  void PamVMCPrimaryGenerator::GenSpe(Doub
166  {  {
167    Double_t alpha = 1.+gamma; //integral spectral index    Double_t alpha = 1.+gamma; //integral spectral index
168    if(alpha==0.){    if(alpha==0.){
169      fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));      fprim.fP0=Exp(Log(PEmin)+frandom->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
170    } else {    } else {
171      if(PEmin==0.) PEmin=1.E-10;      if(PEmin==0.) PEmin=1.E-10;
172      fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);      fprim.fP0=Power((frandom->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
173    }    }
174    
175    if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);    if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23