00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "G4PrimaryTransformer.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4Event.hh"
00033 #include "G4PrimaryVertex.hh"
00034 #include "G4ParticleDefinition.hh"
00035 #include "G4DynamicParticle.hh"
00036 #include "G4Track.hh"
00037 #include "G4ThreeVector.hh"
00038 #include "G4DecayProducts.hh"
00039 #include "G4UnitsTable.hh"
00040 #include "G4ios.hh"
00041 #include "Randomize.hh"
00042
00043 G4PrimaryTransformer::G4PrimaryTransformer()
00044 :verboseLevel(0),trackID(0),unknown(0),unknownParticleDefined(false)
00045 {
00046 particleTable = G4ParticleTable::GetParticleTable();
00047 CheckUnknown();
00048 }
00049
00050 G4PrimaryTransformer::~G4PrimaryTransformer()
00051 {;}
00052
00053 void G4PrimaryTransformer::CheckUnknown()
00054 {
00055 unknown = particleTable->FindParticle("unknown");
00056 if(unknown)
00057 { unknownParticleDefined = true; }
00058 else
00059 { unknownParticleDefined = false; }
00060 }
00061
00062 G4TrackVector* G4PrimaryTransformer::GimmePrimaries(G4Event* anEvent,G4int trackIDCounter)
00063 {
00064 trackID = trackIDCounter;
00065
00066
00067 for( size_t ii=0; ii<TV.size();ii++)
00068 { delete TV[ii]; }
00069 TV.clear();
00070 G4PrimaryVertex* nextVertex = anEvent->GetPrimaryVertex();
00071 while(nextVertex)
00072 {
00073 GenerateTracks(nextVertex);
00074 nextVertex = nextVertex->GetNext();
00075 }
00076 return &TV;
00077 }
00078
00079 void G4PrimaryTransformer::GenerateTracks(G4PrimaryVertex* primaryVertex)
00080 {
00081 G4double X0 = primaryVertex->GetX0();
00082 G4double Y0 = primaryVertex->GetY0();
00083 G4double Z0 = primaryVertex->GetZ0();
00084 G4double T0 = primaryVertex->GetT0();
00085 G4double WV = primaryVertex->GetWeight();
00086
00087 #ifdef G4VERBOSE
00088 if(verboseLevel>2) {
00089 primaryVertex->Print();
00090 } else if (verboseLevel==1) {
00091 G4cout << "G4PrimaryTransformer::PrimaryVertex ("
00092 << X0 / mm << "(mm),"
00093 << Y0 / mm << "(mm),"
00094 << Z0 / mm << "(mm),"
00095 << T0 / nanosecond << "(nsec))" << G4endl;
00096 }
00097 #endif
00098
00099 G4PrimaryParticle* primaryParticle = primaryVertex->GetPrimary();
00100 while( primaryParticle != 0 )
00101 {
00102 GenerateSingleTrack( primaryParticle, X0, Y0, Z0, T0, WV );
00103 primaryParticle = primaryParticle->GetNext();
00104 }
00105 }
00106
00107 void G4PrimaryTransformer::GenerateSingleTrack
00108 (G4PrimaryParticle* primaryParticle,
00109 G4double x0,G4double y0,G4double z0,G4double t0,G4double wv)
00110 {
00111 static G4ParticleDefinition* optPhoton = 0;
00112 static G4int nWarn = 0;
00113 if(!optPhoton) optPhoton = particleTable->FindParticle("opticalphoton");
00114
00115 G4ParticleDefinition* partDef = GetDefinition(primaryParticle);
00116 if(!IsGoodForTrack(partDef))
00117
00118 {
00119 #ifdef G4VERBOSE
00120 if(verboseLevel>2)
00121 {
00122 G4cout << "Primary particle (PDGcode " << primaryParticle->GetPDGcode()
00123 << ") --- Ignored" << G4endl;
00124 }
00125 #endif
00126 G4PrimaryParticle* daughter = primaryParticle->GetDaughter();
00127 while(daughter)
00128 {
00129 GenerateSingleTrack(daughter,x0,y0,z0,t0,wv);
00130 daughter = daughter->GetNext();
00131 }
00132 }
00133
00134
00135 else
00136 {
00137
00138 #ifdef G4VERBOSE
00139 if(verboseLevel>1)
00140 {
00141 G4cout << "Primary particle (" << partDef->GetParticleName()
00142 << ") --- Transfered with momentum " << primaryParticle->GetMomentum()
00143 << G4endl;
00144 }
00145 #endif
00146 G4DynamicParticle* DP =
00147 new G4DynamicParticle(partDef,
00148 primaryParticle->GetMomentumDirection(),
00149 primaryParticle->GetKineticEnergy());
00150 if(partDef==optPhoton && primaryParticle->GetPolarization().mag2()==0.)
00151 {
00152 if(nWarn<10)
00153 {
00154 G4Exception("G4PrimaryTransformer::GenerateSingleTrack","ZeroPolarization",JustWarning,
00155 "Polarization of the optical photon is null. Random polarization is assumed.");
00156 G4cerr << "This warning message is issued up to 10 times." << G4endl;
00157 nWarn++;
00158 }
00159
00160 G4double angle = G4UniformRand() * 360.0*deg;
00161 G4ThreeVector normal (1., 0., 0.);
00162 G4ThreeVector kphoton = DP->GetMomentumDirection();
00163 G4ThreeVector product = normal.cross(kphoton);
00164 G4double modul2 = product*product;
00165
00166 G4ThreeVector e_perpend (0., 0., 1.);
00167 if (modul2 > 0.) e_perpend = (1./std::sqrt(modul2))*product;
00168 G4ThreeVector e_paralle = e_perpend.cross(kphoton);
00169
00170 G4ThreeVector polar = std::cos(angle)*e_paralle + std::sin(angle)*e_perpend;
00171 DP->SetPolarization(polar.x(),polar.y(),polar.z());
00172 }
00173 else
00174 {
00175 DP->SetPolarization(primaryParticle->GetPolX(),
00176 primaryParticle->GetPolY(),
00177 primaryParticle->GetPolZ());
00178 }
00179 if(primaryParticle->GetProperTime()>0.0)
00180 { DP->SetPreAssignedDecayProperTime(primaryParticle->GetProperTime()); }
00181
00182
00183 G4double pmas = primaryParticle->GetMass();
00184 if(pmas>=0.)
00185 { DP->SetMass(pmas); }
00186
00187
00188 if (primaryParticle->GetCharge()<DBL_MAX) {
00189 if (partDef->GetAtomicNumber() <0) {
00190 DP->SetCharge(primaryParticle->GetCharge());
00191 } else {
00192
00193 G4int iz = partDef->GetAtomicNumber();
00194 G4int iq = static_cast<int>(primaryParticle->GetCharge()/eplus);
00195 G4int n_e = iz - iq;
00196 if (n_e>0) DP->AddElectron(0,n_e);
00197 }
00198 }
00199
00200 SetDecayProducts( primaryParticle, DP );
00201
00202 DP->SetPrimaryParticle(primaryParticle);
00203
00204 if(partDef->GetPDGEncoding()==0 && primaryParticle->GetPDGcode()!=0)
00205 {
00206 DP->SetPDGcode(primaryParticle->GetPDGcode());
00207 }
00208
00209 if(!CheckDynamicParticle(DP))
00210 {
00211 delete DP;
00212 return;
00213 }
00214
00215 G4Track* track = new G4Track(DP,t0,G4ThreeVector(x0,y0,z0));
00216
00217 trackID++;
00218 track->SetTrackID(trackID);
00219 primaryParticle->SetTrackID(trackID);
00220
00221 track->SetParentID(0);
00222
00223 track->SetWeight(wv*(primaryParticle->GetWeight()));
00224
00225 TV.push_back( track );
00226
00227 }
00228 }
00229
00230 void G4PrimaryTransformer::SetDecayProducts
00231 (G4PrimaryParticle* mother, G4DynamicParticle* motherDP)
00232 {
00233 G4PrimaryParticle* daughter = mother->GetDaughter();
00234 if(!daughter) return;
00235 G4DecayProducts* decayProducts = (G4DecayProducts*)(motherDP->GetPreAssignedDecayProducts() );
00236 if(!decayProducts)
00237 {
00238 decayProducts = new G4DecayProducts(*motherDP);
00239 motherDP->SetPreAssignedDecayProducts(decayProducts);
00240 }
00241 while(daughter)
00242 {
00243 G4ParticleDefinition* partDef = GetDefinition(daughter);
00244 if(!IsGoodForTrack(partDef))
00245 {
00246 #ifdef G4VERBOSE
00247 if(verboseLevel>2)
00248 {
00249 G4cout << " >> Decay product (PDGcode " << daughter->GetPDGcode()
00250 << ") --- Ignored" << G4endl;
00251 }
00252 #endif
00253 SetDecayProducts(daughter,motherDP);
00254 }
00255 else
00256 {
00257 #ifdef G4VERBOSE
00258 if(verboseLevel>1)
00259 {
00260 G4cout << " >> Decay product (" << partDef->GetParticleName()
00261 << ") --- Attached with momentum " << daughter->GetMomentum()
00262 << G4endl;
00263 }
00264 #endif
00265 G4DynamicParticle*DP
00266 = new G4DynamicParticle(partDef,daughter->GetMomentum());
00267 DP->SetPrimaryParticle(daughter);
00268
00269 if(daughter->GetProperTime()>0.0)
00270 { DP->SetPreAssignedDecayProperTime(daughter->GetProperTime()); }
00271
00272 if (daughter->GetCharge()<DBL_MAX) {
00273 DP->SetCharge(daughter->GetCharge());
00274 }
00275 G4double pmas = daughter->GetMass();
00276 if(pmas>=0.)
00277 { DP->SetMass(pmas); }
00278 decayProducts->PushProducts(DP);
00279 SetDecayProducts(daughter,DP);
00280
00281 if(!CheckDynamicParticle(DP))
00282 {
00283 delete DP;
00284 return;
00285 }
00286 }
00287 daughter = daughter->GetNext();
00288 }
00289 }
00290
00291 void G4PrimaryTransformer::SetUnknnownParticleDefined(G4bool vl)
00292 {
00293 unknownParticleDefined = vl;
00294 if(unknownParticleDefined && !unknown)
00295 { G4cerr << "unknownParticleDefined cannot be set true because G4UnknownParticle is not defined in the physics list."
00296 << G4endl << "Command ignored." << G4endl;
00297 unknownParticleDefined = false;
00298 }
00299 }
00300
00301 G4bool G4PrimaryTransformer::CheckDynamicParticle(G4DynamicParticle*DP)
00302 {
00303 if(IsGoodForTrack(DP->GetDefinition())) return true;
00304 G4DecayProducts* decayProducts = (G4DecayProducts*)(DP->GetPreAssignedDecayProducts());
00305 if(decayProducts && decayProducts->entries()>0) return true;
00306 G4cerr << G4endl
00307 << "G4PrimaryTransformer: a shortlived primary particle is found" << G4endl
00308 << " without any valid decay table nor pre-assigned decay mode." << G4endl;
00309 G4Exception("G4PrimaryTransformer","InvalidPrimary",JustWarning,
00310 "This primary particle will be ignored.");
00311 return false;
00312 }
00313
00314 G4ParticleDefinition* G4PrimaryTransformer::GetDefinition(G4PrimaryParticle*pp)
00315 {
00316 G4ParticleDefinition* partDef = pp->GetG4code();
00317 if(!partDef) partDef = particleTable->FindParticle(pp->GetPDGcode());
00318 if(unknownParticleDefined && ((!partDef)||partDef->IsShortLived())) partDef = unknown;
00319 return partDef;
00320 }
00321
00322 G4bool G4PrimaryTransformer::IsGoodForTrack(G4ParticleDefinition* pd)
00323 {
00324 if(!pd)
00325 { return false; }
00326 else if(!(pd->IsShortLived()))
00327 { return true; }
00328
00329
00330 else if(pd->GetDecayTable())
00331 { return true; }
00332 return false;
00333 }
00334