00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "BesPrimaryGeneratorAction.hh"
00012 #include "BesPrimaryGeneratorMessenger.hh"
00013 #include "globals.hh"
00014 #include "G4Event.hh"
00015 #include "G4ParticleGun.hh"
00016 #include "G4ParticleTable.hh"
00017 #include "G4ParticleDefinition.hh"
00018 #include "G4ParticleMomentum.hh"
00019 #include "Randomize.hh"
00020 #include "G4HEPEvtInterface.hh"
00021 #include "TFile.h"
00022 #include "TH1F.h"
00023
00024 BesPrimaryGeneratorAction::BesPrimaryGeneratorAction()
00025 {
00026 nParticle = 1;
00027 particleName = "pi-";
00028 minCos = -0.8;
00029 maxCos = 0.8;
00030 phiStart = 0.;
00031 phiEnd = 360.;
00032 pMomentum = 1.;
00033 deltaP = 0.;
00034 posX=0;
00035 posY=0;
00036 posZ=0;
00037
00038 particleGun = new G4ParticleGun(1);
00039 isGenbes = true;
00040 HEPEvt = 0;
00041 generatorName = "tester";
00042
00043 messenger = new BesPrimaryGeneratorMessenger(this);
00044
00045 if(generatorName == "cosmic")
00046 {
00047 std::string path = getenv("GENSIMROOT");
00048 G4cout<<"path: "<<path<<G4endl;
00049
00050 path += "/root/";
00051 std::string pFile = path + "ppdc.root";
00052 std::string thetaFile = path + "theta.root";
00053 std::string phiFile = path +"phi.root";
00054
00055 TFile* f1 = new TFile(pFile.c_str());
00056 h1 = (TH1F*)f1->Get("htemp");
00057
00058 TFile* f2 = new TFile(thetaFile.c_str());
00059 h2 = (TH1F*)f2->Get("htemp");
00060
00061 TFile* f3 = new TFile(phiFile.c_str());
00062 h3 = (TH1F*)f3->Get("htemp");
00063
00064
00065
00066
00067 }
00068 }
00069
00070 BesPrimaryGeneratorAction::~BesPrimaryGeneratorAction()
00071 {
00072 delete particleGun;
00073 if(messenger) delete messenger;
00074 if(generatorName=="genbes")delete HEPEvt;
00075 }
00076
00077 void BesPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
00078 {
00079 if(generatorName=="tester")
00080 {
00081 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
00082 G4ParticleDefinition* particle
00083 = particleTable->FindParticle(particleName);
00084 particleGun->SetParticleDefinition(particle);
00085 particleGun->SetParticlePosition(G4ThreeVector(posX,posY,posZ));
00086
00087
00088 for(G4int i=0; i<nParticle; i++)
00089 {
00090 G4double pMag = pMomentum;
00091
00092 if(deltaP>0.)pMag = pMomentum - deltaP*(1.0 - 2.0*G4UniformRand());
00093 pMag = pMag*GeV;
00094
00095 G4double costheta = minCos +(maxCos - minCos)*G4UniformRand();
00096
00097 G4double phi = phiStart+(phiEnd-phiStart)*G4UniformRand();
00098 phi = phi*degree;
00099 G4double sintheta = sqrt(1.-costheta*costheta);
00100
00101 G4ParticleMomentum aMomentum;
00102 aMomentum[0] = pMag*sintheta*cos(phi);
00103 aMomentum[1] = pMag*sintheta*sin(phi);
00104 aMomentum[2] = pMag*costheta;
00105
00106 particleGun->SetParticleMomentum(aMomentum);
00107 particleGun->GeneratePrimaryVertex(anEvent);
00108 }
00109 }
00110 else if(generatorName=="cosmic")
00111 {
00112 G4cout<<"generatorName: "<<generatorName<<G4endl;
00113 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
00114 G4ParticleDefinition* particle
00115 = particleTable->FindParticle(particleName);
00116 G4cout<<"particleName: "<<particleName<<G4endl;
00117
00118 particleGun->SetParticleDefinition(particle);
00119 particleGun->SetParticlePosition(G4ThreeVector(posX,posY,posZ));
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 G4double pMag = h1->GetRandom()*GeV;
00132 G4cout<<"pMag: "<<pMag<<G4endl;
00133
00134
00135
00136
00137
00138 G4double theta =(Double_t)h2->GetRandom();
00139 G4cout<<"theta: "<<theta<<G4endl;
00140 G4double costheta = cos(theta);
00141
00142
00143
00144
00145
00146 G4double phi = (Double_t)h3->GetRandom();
00147 G4cout<<"phi: "<<phi<<G4endl;
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 G4double sintheta = sqrt(1.-costheta*costheta);
00160
00161 G4ParticleMomentum aMomentum;
00162 aMomentum[0] = pMag*sintheta*cos(phi);
00163 aMomentum[1] = pMag*sintheta*sin(phi);
00164 aMomentum[2] = pMag*costheta;
00165
00166 particleGun->SetParticleMomentum(aMomentum);
00167 particleGun->GeneratePrimaryVertex(anEvent);
00168 }
00169
00170 else if(generatorName=="genbes")
00171 {
00172 G4cout<<"genbes called"<<G4endl;
00173 if(isGenbes)
00174 {
00175 isGenbes=false;
00176 HEPEvt=new G4HEPEvtInterface(genbesName);
00177 }
00178 HEPEvt->GeneratePrimaryVertex(anEvent);
00179 }
00180 }
00181