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 #include "G4ExcitedStringDecay.hh"
00029 #include "G4SystemOfUnits.hh"
00030 #include "G4KineticTrack.hh"
00031
00032 G4ExcitedStringDecay::G4ExcitedStringDecay() : G4VStringFragmentation(),
00033 theStringDecay(0)
00034 {}
00035
00036 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay * aStringDecay)
00037 : G4VStringFragmentation(),
00038 theStringDecay(aStringDecay)
00039 {}
00040
00041 G4ExcitedStringDecay::G4ExcitedStringDecay(const G4ExcitedStringDecay &)
00042 : G4VStringFragmentation(),
00043 theStringDecay(0)
00044 {
00045 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible");
00046 }
00047
00048 G4ExcitedStringDecay::~G4ExcitedStringDecay()
00049 {
00050 }
00051
00052 const G4ExcitedStringDecay & G4ExcitedStringDecay::operator=(const G4ExcitedStringDecay &)
00053 {
00054 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable");
00055 return *this;
00056 }
00057
00058 int G4ExcitedStringDecay::operator==(const G4ExcitedStringDecay &) const
00059 {
00060 return 0;
00061 }
00062
00063 int G4ExcitedStringDecay::operator!=(const G4ExcitedStringDecay &) const
00064 {
00065 return 1;
00066 }
00067
00068 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString
00069 (const G4ExcitedString &theString)
00070 {
00071 if ( theStringDecay == NULL )
00072
00073 theStringDecay=new G4LundStringFragmentation();
00074
00075 return theStringDecay->FragmentString(theString);
00076 }
00077
00078 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings
00079 (const G4ExcitedStringVector * theStrings)
00080 {
00081 G4LorentzVector KTsum(0.,0.,0.,0.);
00082
00083
00084 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
00085 {
00086 KTsum+= theStrings->operator[](astring)->Get4Momentum();
00087 }
00088
00089 G4KineticTrackVector * theResult = new G4KineticTrackVector;
00090 G4int attempts(0);
00091 G4bool success=false;
00092 G4bool NeedEnergyCorrector=false;
00093 do {
00094
00095 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
00096 theResult->clear();
00097
00098 attempts++;
00099
00100 G4LorentzVector KTsecondaries(0.,0.,0.,0.);
00101 NeedEnergyCorrector=false;
00102
00103 for ( unsigned int astring=0; astring < theStrings->size(); astring++)
00104 {
00105
00106
00107 G4KineticTrackVector * generatedKineticTracks = NULL;
00108 if ( theStrings->operator[](astring)->IsExcited() )
00109 {
00110
00111 generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
00112 } else {
00113 generatedKineticTracks = new G4KineticTrackVector;
00114 generatedKineticTracks->push_back(theStrings->operator[](astring)->GetKineticTrack());
00115 }
00116
00117 if (generatedKineticTracks == NULL)
00118 {
00119 G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl;
00120 continue;
00121 }
00122
00123 G4LorentzVector KTsum1(0.,0.,0.,0.);
00124 for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
00125 {
00126
00127 theResult->push_back(generatedKineticTracks->operator[](aTrack));
00128 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
00129 }
00130 KTsecondaries+=KTsum1;
00131
00132
00133
00134 if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
00135 {
00136
00137
00138 NeedEnergyCorrector=true;
00139 }
00140
00141
00142 delete generatedKineticTracks;
00143 }
00144
00145
00146 success=true;
00147
00148 if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
00149
00150 } while(!success && (attempts < 10));
00151
00152
00153 #ifdef debug_ExcitedStringDecay
00154 G4LorentzVector KTsum1=0;
00155 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
00156 {
00157 G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
00158 <<" " << (*theResult)[aTrack]->Get4Momentum() << G4endl;
00159 KTsum1+= (*theResult)[aTrack]->Get4Momentum();
00160 }
00161
00162 G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success << ", Corrected total 4 momentum " << KTsum1 << G4endl;
00163 if ( ! success ) G4cout << "failed to correct E/p" << G4endl;
00164 #endif
00165
00166 return theResult;
00167 }
00168
00169 G4bool G4ExcitedStringDecay::EnergyAndMomentumCorrector
00170 (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)
00171 {
00172 const int nAttemptScale = 500;
00173 const double ErrLimit = 1.E-5;
00174 if (Output->empty())
00175 return TRUE;
00176 G4LorentzVector SumMom;
00177 G4double SumMass = 0;
00178 G4double TotalCollisionMass = TotalCollisionMom.m();
00179
00180
00181
00182 unsigned int cHadron;
00183 for(cHadron = 0; cHadron < Output->size(); cHadron++)
00184 {
00185 SumMom += Output->operator[](cHadron)->Get4Momentum();
00186 SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass();
00187 }
00188
00189
00190
00191
00192 if (Output->size() < 2) return FALSE;
00193
00194 if (SumMass > TotalCollisionMass) return FALSE;
00195 SumMass = SumMom.m2();
00196 if (SumMass < 0) return FALSE;
00197 SumMass = std::sqrt(SumMass);
00198
00199
00200 G4ThreeVector Beta = -SumMom.boostVector();
00201 Output->Boost(Beta);
00202
00203
00204
00205 G4double Scale = 1;
00206 G4int cAttempt = 0;
00207 G4double Sum = 0;
00208 G4bool success = false;
00209 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
00210 {
00211 Sum = 0;
00212 for(cHadron = 0; cHadron < Output->size(); cHadron++)
00213 {
00214 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
00215 HadronMom.setVect(Scale*HadronMom.vect());
00216 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
00217 HadronMom.setE(E);
00218 Output->operator[](cHadron)->Set4Momentum(HadronMom);
00219 Sum += E;
00220 }
00221 Scale = TotalCollisionMass/Sum;
00222 #ifdef debug_G4ExcitedStringDecay
00223 G4cout << "Scale-1=" << Scale -1
00224 << ", TotalCollisionMass=" << TotalCollisionMass
00225 << ", Sum=" << Sum
00226 << G4endl;
00227 #endif
00228 if (std::fabs(Scale - 1) <= ErrLimit)
00229 {
00230 success = true;
00231 break;
00232 }
00233 }
00234 #ifdef debug_G4ExcitedStringDecay
00235 if(!success)
00236 {
00237 G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl;
00238 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
00239 G4cout << " Number of secondaries: " << Output->size() << G4endl;
00240 G4cout << " Wanted total energy: " << TotalCollisionMom.e() << G4endl;
00241 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl;
00242
00243 }
00244 #endif
00245
00246 Beta = TotalCollisionMom.boostVector();
00247 Output->Boost(Beta);
00248
00249 return success;
00250 }