PamVMC/src/PamVMCPrimaryGenerator.cxx

// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 


#include <TVirtualMC.h>
#include <TVirtualMCStack.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TVector3.h>
#include <TMath.h>

#include "PamVMCPrimaryGenerator.h"

using TMath::Sqrt;
using TMath::Sin;
using TMath::Cos;
using TMath::ATan;
using TMath::Log;
using TMath::Power;
using TMath::Exp;

ClassImp(PamVMCPrimary)

PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b)
{
  a.fPDG=b.fPDG;
  a.fX0=b.fX0;
  a.fY0=b.fY0;
  a.fZ0=b.fZ0;
  a.fTHETA=b.fTHETA;
  a.fPHI=b.fPHI;
  a.fP0=b.fP0;
  a.fGOOD=b.fGOOD;

  return a;
}


ClassImp(PamVMCPrimaryGenerator)

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) 
  : TObject(),
    fStack(stack),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    frnd(0)
{
// Standard constructor
  fprimColl = new TClonesArray("PamVMCPrimary");
  frnd = new TRandom3(0);

  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV

} 

PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
  : TObject(),
    fStack(0),
    fevno(0),
    fmass(0.),
    fcharge(0.),
    fprimColl(0),
    frnd(0)
{    
  frnd = new TRandom3(0);
// Default constructor
  //Default primary proton
  fprim.fPDG=kProton;
  fprim.fX0=1.;
  fprim.fY0=1.;
  fprim.fZ0=130.;
  fprim.fTHETA=0.;
  fprim.fPHI=0.;
  fprim.fP0=1.; //1GV
}

PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() 
{
// Destructor
  delete frnd;
  delete fprimColl;
}

// private methods

#include <Riostream.h>

void PamVMCPrimaryGenerator::GeneratePrimary()
{    
// Add one primary particle to the user stack (derived from TVirtualMCStack).
  
  // Track ID (filled by stack)
  Int_t ntr;
 
  // Option: to be tracked
  Int_t toBeDone = 1; 
 
  // Particle type
  Int_t pdg  = fprim.fPDG;
  
  Double_t fvx, fvy, fvz;
  fvx=fprim.fX0;
  fvy=fprim.fY0;
  fvz=fprim.fZ0;

  // Position
  
  Double_t tof = 0.;

  // Energy (in GeV)
  Double_t kinEnergy = MomentumToKinE(fprim.fP0);     
  Double_t e  = fmass + kinEnergy;
 
  // Particle momentum
  Double_t  px, py, pz;
   
  px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI);
  py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI);
  pz = -fprim.fP0*Cos(fprim.fTHETA);
 
  // Polarization
  TVector3 polar;

  // Add particle to stack 
  fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof, 
                   polar.X(), polar.Y(), polar.Z(), 
                   kPPrimary, ntr, 1., 0);

  PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
 
  *pc = fprim;
}


void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){
  fprim.fPDG=pdg;
  //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG);
  fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass();
  fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.;
}

void PamVMCPrimaryGenerator::SetMomentum(
                              Double_t px, Double_t py, Double_t pz) 
{
  fprim.fP0= Sqrt(px*px+py*py+pz*pz);
  fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz);
  fprim.fPHI=ATan(py/px);
}
                                     
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy)
{
  if(isEnergy) {
    fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax));
  } else{
    fprim.fP0=frnd->Uniform(PEmin,PEmax);
  }

}

void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy)
{
  Double_t alpha = 1.+gamma; //integral spectral index
  if(alpha==0.){
    fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin)));
  } else {
    if(PEmin==0.) PEmin=1.E-10;
    fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha);
  }

  if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);

}
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
12			#include "PamVMCPrimaryGenerator.h"
13
14	nikolas	1.2	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	nikolas	1.1	ClassImp(PamVMCPrimaryGenerator)
40
41			PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack)
42			: TObject(),
43			fStack(stack),
44	nikolas	1.2	fevno(0),
45			fmass(0.),
46			fcharge(0.),
47			frnd(0)
48	nikolas	1.1	{
49			// Standard constructor
50	nikolas	1.2	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	nikolas	1.1
61	nikolas	1.2	}
62	nikolas	1.1
63			PamVMCPrimaryGenerator::PamVMCPrimaryGenerator()
64			: TObject(),
65			fStack(0),
66	nikolas	1.2	fevno(0),
67			fmass(0.),
68			fcharge(0.),
69			fprimColl(0),
70			frnd(0)
71	nikolas	1.1	{
72	nikolas	1.2	frnd = new TRandom3(0);
73	nikolas	1.1	// Default constructor
74	nikolas	1.2	//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	nikolas	1.1	}
83
84			PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator()
85			{
86	nikolas	1.2	// Destructor
87			delete frnd;
88			delete fprimColl;
89	nikolas	1.1	}
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	nikolas	1.2	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	nikolas	1.1	// Position
114	nikolas	1.2
115	nikolas	1.1	Double_t tof = 0.;
116
117			// Energy (in GeV)
118	nikolas	1.2	Double_t kinEnergy = MomentumToKinE(fprim.fP0);
119			Double_t e = fmass + kinEnergy;
120	nikolas	1.1
121			// Particle momentum
122	nikolas	1.2	Double_t px, py, pz;
123
124			px = fprim.fP0Sin(fprim.fTHETA)Cos(fprim.fPHI);
125			py = fprim.fP0Sin(fprim.fTHETA)Sin(fprim.fPHI);
126			pz = -fprim.fP0*Cos(fprim.fTHETA);
127	nikolas	1.1
128			// Polarization
129			TVector3 polar;
130
131			// Add particle to stack
132	nikolas	1.2	fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof,
133	nikolas	1.1	polar.X(), polar.Y(), polar.Z(),
134			kPPrimary, ntr, 1., 0);
135	nikolas	1.2
136			PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++);
137
138			*pc = fprim;
139	nikolas	1.1	}
140
141
142	nikolas	1.2	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	nikolas	1.1
149	nikolas	1.2	void PamVMCPrimaryGenerator::SetMomentum(
150			Double_t px, Double_t py, Double_t pz)
151			{
152			fprim.fP0= Sqrt(pxpx+pypy+pz*pz);
153			fprim.fTHETA=ATan(Sqrt(pxpx+pypy)/pz);
154			fprim.fPHI=ATan(py/px);
155	nikolas	1.1	}
156	nikolas	1.2
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	nikolas	1.1
165	nikolas	1.2	}
166	nikolas	1.1
167	nikolas	1.2	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	nikolas	1.1
177	nikolas	1.2	if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0);
178	nikolas	1.1
179	nikolas	1.2	}