| 1 |
// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 |
| 2 |
|
| 3 |
|
| 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 |
|
| 12 |
#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) |
| 40 |
|
| 41 |
PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) |
| 42 |
: TObject(), |
| 43 |
fStack(stack), |
| 44 |
fevno(0), |
| 45 |
fmass(0.), |
| 46 |
fcharge(0.), |
| 47 |
frnd(0) |
| 48 |
{ |
| 49 |
// 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() |
| 64 |
: TObject(), |
| 65 |
fStack(0), |
| 66 |
fevno(0), |
| 67 |
fmass(0.), |
| 68 |
fcharge(0.), |
| 69 |
fprimColl(0), |
| 70 |
frnd(0) |
| 71 |
{ |
| 72 |
frnd = new TRandom3(0); |
| 73 |
// 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() |
| 85 |
{ |
| 86 |
// Destructor |
| 87 |
delete frnd; |
| 88 |
delete fprimColl; |
| 89 |
} |
| 90 |
|
| 91 |
// private methods |
| 92 |
|
| 93 |
#include <Riostream.h> |
| 94 |
|
| 95 |
void PamVMCPrimaryGenerator::GeneratePrimary() |
| 96 |
{ |
| 97 |
// Add one primary particle to the user stack (derived from TVirtualMCStack). |
| 98 |
|
| 99 |
// Track ID (filled by stack) |
| 100 |
Int_t ntr; |
| 101 |
|
| 102 |
// Option: to be tracked |
| 103 |
Int_t toBeDone = 1; |
| 104 |
|
| 105 |
// Particle type |
| 106 |
Int_t pdg = fprim.fPDG; |
| 107 |
|
| 108 |
Double_t fvx, fvy, fvz; |
| 109 |
fvx=fprim.fX0; |
| 110 |
fvy=fprim.fY0; |
| 111 |
fvz=fprim.fZ0; |
| 112 |
|
| 113 |
// Position |
| 114 |
|
| 115 |
Double_t tof = 0.; |
| 116 |
|
| 117 |
// Energy (in GeV) |
| 118 |
Double_t kinEnergy = MomentumToKinE(fprim.fP0); |
| 119 |
Double_t e = fmass + kinEnergy; |
| 120 |
|
| 121 |
// Particle momentum |
| 122 |
Double_t px, py, pz; |
| 123 |
|
| 124 |
px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI); |
| 125 |
py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI); |
| 126 |
pz = -fprim.fP0*Cos(fprim.fTHETA); |
| 127 |
|
| 128 |
// Polarization |
| 129 |
TVector3 polar; |
| 130 |
|
| 131 |
// Add particle to stack |
| 132 |
fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof, |
| 133 |
polar.X(), polar.Y(), polar.Z(), |
| 134 |
kPPrimary, ntr, 1., 0); |
| 135 |
|
| 136 |
PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++); |
| 137 |
|
| 138 |
*pc = fprim; |
| 139 |
} |
| 140 |
|
| 141 |
|
| 142 |
void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){ |
| 143 |
fprim.fPDG=pdg; |
| 144 |
//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 |
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 |
} |
| 166 |
|
| 167 |
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy) |
| 168 |
{ |
| 169 |
Double_t alpha = 1.+gamma; //integral spectral index |
| 170 |
if(alpha==0.){ |
| 171 |
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 |
} |
Parent Directory
| 
Revision Log