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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "G4VPartonStringModel.hh"
00041 #include "G4ios.hh"
00042 #include "G4ShortLivedConstructor.hh"
00043
00044 #include "G4ParticleTable.hh"
00045 #include "G4IonTable.hh"
00046
00047
00048 G4VPartonStringModel::G4VPartonStringModel(const G4String& modelName)
00049 : G4VHighEnergyGenerator(modelName),
00050 stringFragmentationModel(0),
00051 theThis(0)
00052 {
00053
00054 G4ShortLivedConstructor ShortLived;
00055 ShortLived.ConstructParticle();
00056 }
00057
00058 G4VPartonStringModel::~G4VPartonStringModel()
00059 {
00060 }
00061
00062 G4KineticTrackVector * G4VPartonStringModel::Scatter(const G4Nucleus &theNucleus,
00063 const G4DynamicParticle &aPrimary)
00064 {
00065 G4ExcitedStringVector * strings = NULL;
00066
00067 G4DynamicParticle thePrimary=aPrimary;
00068
00069 G4LorentzRotation toZ;
00070 G4LorentzVector Ptmp=thePrimary.Get4Momentum();
00071 toZ.rotateZ(-1*Ptmp.phi());
00072 toZ.rotateY(-1*Ptmp.theta());
00073 thePrimary.Set4Momentum(toZ*Ptmp);
00074 G4LorentzRotation toLab(toZ.inverse());
00075
00076 G4int attempts = 0, maxAttempts=20;
00077 while ( strings == NULL )
00078 {
00079 if (attempts++ > maxAttempts )
00080 {
00081 throw G4HadronicException(__FILE__, __LINE__, "G4VPartonStringModel::Scatter(): fails to generate strings");
00082 }
00083
00084 theThis->Init(theNucleus,thePrimary);
00085
00086 strings = GetStrings();
00087 }
00088
00089 G4double stringEnergy(0);
00090 G4LorentzVector SumStringMom(0.,0.,0.,0.);
00091
00092 for ( unsigned int astring=0; astring < strings->size(); astring++)
00093 {
00094
00095 stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
00096 stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
00097 (*strings)[astring]->LorentzRotate(toLab);
00098 SumStringMom+=(*strings)[astring]->Get4Momentum();
00099 }
00100
00101 G4double InvMass=SumStringMom.mag();
00102
00103
00104 #ifdef debug_PartonStringModel
00105
00106 G4V3DNucleus * fancynucleus=theThis->GetWoundedNucleus();
00107
00108
00109 G4int hits(0), charged_hits(0);
00110 G4ThreeVector hitNucleonMomentum(0.,0.,0.);
00111 G4Nucleon * theCurrentNucleon = fancynucleus->StartLoop() ? fancynucleus->GetNextNucleon() : NULL;
00112 while( theCurrentNucleon )
00113 {
00114 if(theCurrentNucleon->AreYouHit())
00115 {
00116 hits++;
00117 hitNucleonMomentum += theCurrentNucleon->Get4Momentum().vect();
00118 if ( theCurrentNucleon->GetDefinition() == G4Proton::Proton() ) ++charged_hits;
00119 }
00120 theCurrentNucleon = fancynucleus->GetNextNucleon();
00121 }
00122
00123 G4int initialZ=fancynucleus->GetCharge();
00124 G4int initialA=fancynucleus->GetMassNumber();
00125 G4double initial_mass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ,initialA);
00126 G4double final_mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(initialZ-charged_hits, initialA-hits);
00127 G4cout << "G4VPSM: strE, hit nucleons, Primary, SumStringE, nucleus intial, nucleus final, excitation estimate "
00128 << stringEnergy << " " << hits << ", " << Ptmp.e() << ", "<< SumStringMom.e() << ", "
00129 << initial_mass<< ", " << final_mass<< ", " << Ptmp.e() + initial_mass - final_mass - stringEnergy << G4endl;
00130 G4cout << "momentum balance " << thePrimary.GetMomentum() + hitNucleonMomentum - SumStringMom.vect()<< G4endl;
00131 #endif
00132
00133
00134
00135 G4KineticTrackVector * theResult = 0;
00136 G4double SumMass(0.);
00137 attempts = 0;
00138 maxAttempts=100;
00139 do
00140 {
00141 attempts++;
00142 if(theResult != 0)
00143 {
00144 std::for_each(theResult->begin(), theResult->end(), DeleteKineticTrack());
00145 delete theResult;
00146 }
00147 theResult = stringFragmentationModel->FragmentStrings(strings);
00148 if(attempts > maxAttempts ) break;
00149
00150
00151
00152 SumMass=0.;
00153
00154 for ( unsigned int i=0; i < theResult->size(); i++)
00155 {
00156 SumMass+=(*theResult)[i]->GetDefinition()->GetPDGMass();
00157
00158
00159
00160 }
00161
00162
00163 } while(SumMass > InvMass);
00164
00165 std::for_each(strings->begin(), strings->end(), DeleteString() );
00166 delete strings;
00167
00168 return theResult;
00169 }
00170
00171 void G4VPartonStringModel::ModelDescription(std::ostream& outFile) const
00172 {
00173 outFile << GetModelName() << " has no description yet.\n";
00174 }