Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Historic fragment from M.Komogorov; clean-u 26 // Historic fragment from M.Komogorov; clean-up still necessary @@@ 27 << 28 #include "G4ExcitedStringDecay.hh" 27 #include "G4ExcitedStringDecay.hh" 29 #include "G4SystemOfUnits.hh" << 30 #include "G4KineticTrack.hh" 28 #include "G4KineticTrack.hh" 31 #include "G4LundStringFragmentation.hh" << 32 #include "G4HadronicInteractionRegistry.hh" << 33 << 34 #include "G4SampleResonance.hh" << 35 << 36 //#define debug_G4ExcitedStringDecay << 37 //#define debug_G4ExcitedStringCorr << 38 29 39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V << 30 G4ExcitedStringDecay::G4ExcitedStringDecay() 40 : G4VStringFragmentation(), theStringDecay(p << 41 { 31 { 42 if(!ptr) { << 32 theStringDecay=NULL; 43 G4HadronicInteraction* p = << 44 G4HadronicInteractionRegistry::Instance( << 45 theStringDecay = static_cast<G4VLongitudin << 46 if(!theStringDecay) { theStringDecay = new << 47 } << 48 SetModelName(theStringDecay->GetModelName()) << 49 } 33 } 50 34 51 G4ExcitedStringDecay::~G4ExcitedStringDecay() << 35 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay * aStringDecay) >> 36 : theStringDecay(aStringDecay) 52 {} 37 {} 53 38 54 G4KineticTrackVector *G4ExcitedStringDecay::Fr << 39 G4ExcitedStringDecay::G4ExcitedStringDecay(const G4ExcitedStringDecay &) : G4VStringFragmentation() 55 { 40 { 56 return theStringDecay->FragmentString(theStr << 41 } 57 } << 58 42 59 G4KineticTrackVector *G4ExcitedStringDecay::Fr << 43 G4ExcitedStringDecay::~G4ExcitedStringDecay() 60 { 44 { 61 G4LorentzVector KTsum(0.,0.,0.,0.); << 45 } 62 << 63 #ifdef debug_G4ExcitedStringDecay << 64 G4cout<<G4endl; << 65 G4cout<<"--------------------------- G4Excit << 66 G4cout<<"Hadronization of Excited Strings: t << 67 #endif << 68 << 69 for ( unsigned int astring=0; astring < theS << 70 // for ( unsigned int astring=0; astring < 1 << 71 { << 72 // G4cout<<"theStrings->operator[](astring << 73 if ( theStrings->operator[](astring)->IsEx << 74 {KTsum+= theStrings->operator[](astri << 75 else {KTsum+=theStrings->operator[](astrin << 76 } << 77 << 78 G4LorentzRotation toCms( -1 * KTsum.boostVec << 79 G4LorentzRotation toLab(toCms.inverse()); << 80 G4LorentzVector Ptmp; << 81 KTsum=G4LorentzVector(0.,0.,0.,0.); << 82 << 83 for ( unsigned int astring=0; astring < theS << 84 // for ( unsigned int astring=0; astring < 1 << 85 { << 86 if ( theStrings->operator[](astring)->IsEx << 87 { << 88 Ptmp=toCms * theStrings->operator[](astr << 89 theStrings->operator[](astring)->GetLeft << 90 << 91 Ptmp=toCms * theStrings->operator[](astr << 92 theStrings->operator[](astring)->GetRigh << 93 << 94 KTsum+= theStrings->operator[](astring)- << 95 } << 96 else << 97 { << 98 Ptmp=toCms * theStrings->operator[](astr << 99 theStrings->operator[](astring)->GetKine << 100 KTsum+= theStrings->operator[](astring)- << 101 } << 102 } << 103 << 104 G4SampleResonance BrW; << 105 const G4ParticleDefinition* TrackDefinition= << 106 << 107 G4KineticTrackVector * theResult = new G4Kin << 108 G4int attempts(0); << 109 G4bool success=false; << 110 G4bool NeedEnergyCorrector=false; << 111 do { << 112 #ifdef debug_G4ExcitedStringDecay << 113 G4cout<<"New try No "<<attempts<<" to << 114 #endif << 115 << 116 std::for_each(theResult->begin() , theResult << 117 theResult->clear(); << 118 << 119 attempts++; << 120 << 121 G4LorentzVector KTsecondaries(0.,0.,0.,0.); << 122 NeedEnergyCorrector=false; << 123 << 124 for ( unsigned int astring=0; astring < theS << 125 // for ( unsigned int astring=0; astri << 126 // for ( unsigned int astring=1; astri << 127 // for ( unsigned int astring=0; astri << 128 { << 129 #ifdef debug_G4ExcitedStringDecay << 130 G4cout<<"String No "<<astring+1<<" E << 131 << 132 G4cout<<"String No "<<astring+1<<" 4 << 133 <<" "<<theStrings->operator[](astrin << 134 #endif << 135 << 136 G4KineticTrackVector * generatedKine << 137 if ( theStrings->operator[](astring)->IsEx << 138 { << 139 #ifdef debug_G4ExcitedStringDecay << 140 G4cout<<"Fragment String with par << 141 <<theStrings->operator[](as << 142 <<theStrings->operator[](as << 143 <<"Direction "<<theStrings- << 144 #endif << 145 generatedKineticTracks=FragmentString << 146 #ifdef debug_G4ExcitedStringDecay << 147 G4cout<<"(G4ExcitedStringDecay) N << 148 <<generatedKineticTracks->s << 149 #endif << 150 } else { << 151 #ifdef debug_G4ExcitedStringDecay << 152 G4cout<<" GetTrack from the Str << 153 #endif << 154 G4LorentzVector Mom=theStrings->o << 155 G4KineticTrack * aTrack= new G4Ki << 156 theStrings->operat << 157 theStrings->operat << 158 G4ThreeVector(0), << 159 << 160 aTrack->SetPosition(theStrings->o << 161 << 162 #ifdef debug_G4ExcitedStringDecay << 163 G4cout<<" A particle stored in << 164 #endif << 165 << 166 generatedKineticTracks = new G4KineticT << 167 generatedKineticTracks->push_back(aTrac << 168 } << 169 << 170 if (generatedKineticTracks == nullptr || g << 171 { << 172 // G4cerr << "G4VPartonStringMode << 173 delete generatedKineticTracks; << 174 generatedKineticTracks = nullptr; << 175 continue; << 176 } << 177 << 178 G4LorentzVector KTsum1(0.,0.,0.,0.); << 179 for ( unsigned int aTrack=0; aTrack< << 180 { << 181 #ifdef debug_G4ExcitedStringDecay << 182 G4cout<<"Prod part No. "<<aTrack+ << 183 <<(*generatedKineticTracks) << 184 <<(*generatedKineticTracks) << 185 <<(*generatedKineticTracks) << 186 #endif << 187 // --------------- Sampling mass << 188 TrackDefinition = (*generatedKine << 189 << 190 if (TrackDefinition->IsShortLived()) << 191 { << 192 G4double NewTrackMass = << 193 BrW.SampleMass( TrackDefiniti << 194 BrW.GetMinimu << 195 TrackDefiniti << 196 G4LorentzVector Tmp=G4LorentzVe << 197 Tmp.setE(std::sqrt(sqr(NewTrack << 198 << 199 (*generatedKineticTracks)[aTrac << 200 << 201 #ifdef debug_G4ExcitedStringDec << 202 G4cout<<"Resonance *** "<<aTrac << 203 <<(*generatedKineticTrack << 204 <<(*generatedKineticTrack << 205 <<(*generatedKineticTrack << 206 #endif << 207 } << 208 //------------------------------- << 209 << 210 theResult->push_back(generatedKin << 211 KTsum1+= (*generatedKineticTracks << 212 } << 213 KTsecondaries+=KTsum1; << 214 << 215 #ifdef debug_G4ExcitedStringDecay << 216 G4cout << "String secondaries(" <<ge << 217 <<"Init string momentum: "< << 218 <<"Final hadrons momentum: "< << 219 #endif << 220 << 221 if ( KTsum1.e() > 0 && << 222 std::abs((KTsum1.e()-theString << 223 { << 224 NeedEnergyCorrector=true; << 225 } << 226 << 227 #ifdef debug_G4ExcitedStringDecay << 228 G4cout<<"NeedEnergyCorrection yes/no << 229 #endif << 230 << 231 // clean up << 232 delete generatedKineticTracks; << 233 success=true; << 234 } << 235 << 236 if ( NeedEnergyCorrector ) success=EnergyAnd << 237 } while (!success && (attempts < 100)); /* << 238 << 239 for ( unsigned int aTrack=0; aTrack<theResul << 240 { << 241 Ptmp=(*theResult)[aTrack]->Get4Momentum(); << 242 Ptmp.transform( toLab); << 243 (*theResult)[aTrack]->Set4Momentum(Ptmp); << 244 } << 245 << 246 #ifdef debug_G4ExcitedStringDecay << 247 G4cout<<"End of the strings fragmentation (G << 248 << 249 G4LorentzVector KTsum1(0.,0.,0.,0.); << 250 << 251 for ( unsigned int aTrack=0; aTrack<theResul << 252 { << 253 G4cout << " corrected tracks .. " << (*t << 254 <<" " << (*theResult)[aTrack]->Get4Mome << 255 <<" " << (*theResult)[aTrack]->Get4Mome << 256 KTsum1+= (*theResult)[aTrack]->Get4Momen << 257 } << 258 << 259 G4cout << "Needcorrector/success " << NeedEn << 260 << ", Corrected total 4 momentum " < << 261 if ( ! success ) G4cout << "failed to correc << 262 << 263 G4cout<<"End of the Hadronization (G4Excited << 264 #endif << 265 46 266 if (!success) << 47 const G4ExcitedStringDecay & G4ExcitedStringDecay::operator=(const G4ExcitedStringDecay &) 267 { << 48 { 268 if (theResult->size() != 0) << 49 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable"); 269 { << 50 return *this; 270 std::for_each(theResult->begin() , theRe << 51 } 271 theResult->clear(); << 272 delete theResult; theResult=0; << 273 } << 274 for ( unsigned int astring=0; astring < th << 275 // for ( unsigned int astring=0; astring < << 276 { << 277 if ( theStrings->operator[](astring)->Is << 278 { << 279 Ptmp=theStrings->operator[](astring)-> << 280 Ptmp.transform( toLab); << 281 theStrings->operator[](astring)->GetLe << 282 52 283 Ptmp=theStrings->operator[](astring)-> << 53 int G4ExcitedStringDecay::operator==(const G4ExcitedStringDecay &) const 284 Ptmp.transform( toLab); << 54 { 285 theStrings->operator[](astring)->GetRi << 55 return 0; 286 } << 287 else << 288 { << 289 Ptmp=theStrings->operator[](astring)-> << 290 Ptmp.transform( toLab); << 291 theStrings->operator[](astring)->GetKi << 292 } << 293 } << 294 } << 295 return theResult; << 296 } 56 } 297 57 >> 58 int G4ExcitedStringDecay::operator!=(const G4ExcitedStringDecay &) const >> 59 { >> 60 return 1; >> 61 } 298 62 299 G4bool G4ExcitedStringDecay:: 63 G4bool G4ExcitedStringDecay:: 300 EnergyAndMomentumCorrector(G4KineticTrackVecto 64 EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom) 301 { << 65 { 302 const int nAttemptScale = 500; 66 const int nAttemptScale = 500; 303 const double ErrLimit = 1.E-5; 67 const double ErrLimit = 1.E-5; 304 if (Output->empty()) return TRUE; << 68 if (Output->empty()) >> 69 return TRUE; 305 G4LorentzVector SumMom; 70 G4LorentzVector SumMom; 306 G4double SumMass = 0; << 71 G4double SumMass = 0; 307 G4double TotalCollisionMass = Total 72 G4double TotalCollisionMass = TotalCollisionMom.m(); 308 << 73 if( !(TotalCollisionMass<1) && !(TotalCollisionMass>-1) ) 309 std::vector<G4double> HadronMass; G4double << 74 { 310 << 75 std::cout << "TotalCollisionMomentum = "<<TotalCollisionMom<<G4endl; 311 #ifdef debug_G4ExcitedStringCorr << 76 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay received nan mass..."); 312 G4cout<<G4endl<<"EnergyAndMomentumCorrecto << 77 } 313 #endif << 314 // Calculate sum hadron 4-momenta and summ 78 // Calculate sum hadron 4-momenta and summing hadron mass 315 unsigned int cHadron; 79 unsigned int cHadron; 316 for (cHadron = 0; cHadron < Output->size() << 80 for(cHadron = 0; cHadron < Output->size(); cHadron++) 317 { 81 { 318 SumMom += Output->operator[](cHadron) 82 SumMom += Output->operator[](cHadron)->Get4Momentum(); 319 HadronM=Output->operator[](cHadron)->G << 83 if( !(SumMom<1) && !(SumMom>-1) ) 320 SumMass += Output->operator[](cHadron) << 84 { >> 85 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::EnergyAndMomentumCorrector() received nan momentum..."); >> 86 } >> 87 SumMass += Output->operator[](cHadron)->GetDefinition()->GetPDGMass(); >> 88 if( !(SumMass<1) && !(SumMass>-1) ) >> 89 { >> 90 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::EnergyAndMomentumCorrector() received nan mass..."); >> 91 } 321 } 92 } 322 << 323 #ifdef debug_G4ExcitedStringCorr << 324 G4cout<<"Sum part mom "<<SumMom<<" "<<SumM << 325 <<"Sum str mom "<<TotalCollisionMom << 326 G4cout<<"SumMass TotalCollisionMass "<<Sum << 327 #endif << 328 << 329 // Cannot correct a single particle 93 // Cannot correct a single particle 330 if (Output->size() < 2) return FALSE; 94 if (Output->size() < 2) return FALSE; 331 << 95 332 if (SumMass > TotalCollisionMass) return F 96 if (SumMass > TotalCollisionMass) return FALSE; 333 SumMass = SumMom.m2(); 97 SumMass = SumMom.m2(); 334 if (SumMass < 0) return FALSE; 98 if (SumMass < 0) return FALSE; 335 SumMass = std::sqrt(SumMass); 99 SumMass = std::sqrt(SumMass); 336 100 337 // Compute c.m.s. hadron velocity and boos << 101 // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s. 338 // G4ThreeVector Beta = -SumMom.boostVecto << 102 G4ThreeVector Beta = -SumMom.boostVector(); 339 G4ThreeVector Beta = -TotalCollisionMom.bo << 340 Output->Boost(Beta); 103 Output->Boost(Beta); 341 104 342 // Scale total c.m.s. hadron energy (hadro 105 // Scale total c.m.s. hadron energy (hadron system mass). 343 // It should be equal interaction mass 106 // It should be equal interaction mass 344 G4double Scale = 1; 107 G4double Scale = 1; 345 G4int cAttempt = 0; 108 G4int cAttempt = 0; 346 G4double Sum = 0; 109 G4double Sum = 0; 347 G4bool success = false; 110 G4bool success = false; 348 for (cAttempt = 0; cAttempt < nAttemptScal << 111 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++) 349 { 112 { 350 Sum = 0; 113 Sum = 0; 351 for (cHadron = 0; cHadron < Output->size << 114 for(cHadron = 0; cHadron < Output->size(); cHadron++) 352 { 115 { 353 HadronM = HadronMass.at(cHadron); << 354 G4LorentzVector HadronMom = Output->op 116 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum(); 355 HadronMom.setVect(Scale*HadronMom.vect 117 HadronMom.setVect(Scale*HadronMom.vect()); 356 G4double E = std::sqrt(HadronMom.vect( << 118 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass())); 357 << 358 HadronMom.setE(E); 119 HadronMom.setE(E); 359 Output->operator[](cHadron)->Set4Momen 120 Output->operator[](cHadron)->Set4Momentum(HadronMom); 360 Sum += E; 121 Sum += E; 361 } 122 } 362 Scale = TotalCollisionMass/Sum; << 123 Scale = TotalCollisionMass/Sum; 363 #ifdef debug_G4ExcitedStringCorr << 124 #ifdef debug_G4ExcitedStringDecay 364 G4cout << "Scale-1=" << Scale -1 125 G4cout << "Scale-1=" << Scale -1 365 << ", TotalCollisionMass=" << To << 126 << ", TotalCollisionMass=" << TotalCollisionMass 366 << ", Sum=" << Sum << 127 << ", Sum=" << Sum 367 << G4endl; << 128 << G4endl; 368 #endif << 129 #endif 369 if (std::fabs(Scale - 1) <= ErrLimit) 130 if (std::fabs(Scale - 1) <= ErrLimit) 370 { 131 { 371 success = true; 132 success = true; 372 break; 133 break; 373 } 134 } 374 } 135 } 375 << 136 #ifdef debug_G4ExcitedStringDecay 376 #ifdef debug_G4ExcitedStringCorr << 137 if(!success) 377 if (!success) << 378 { 138 { 379 G4cout << "G4ExcitedStringDecay::EnergyA 139 G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl; 380 G4cout << " Scale not unity at end of 140 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl; 381 G4cout << " Number of secondaries: " < 141 G4cout << " Number of secondaries: " << Output->size() << G4endl; 382 G4cout << " Wanted total energy: " << 142 G4cout << " Wanted total energy: " << TotalCollisionMom.e() << G4endl; 383 G4cout << " Increase number of attempt 143 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl; 384 // throw G4HadronicException(__FILE__, _ << 144 // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct..."); 385 } 145 } 386 #endif << 146 #endif 387 << 388 // Compute c.m.s. interaction velocity and 147 // Compute c.m.s. interaction velocity and KTV back boost 389 Beta = TotalCollisionMom.boostVector(); 148 Beta = TotalCollisionMom.boostVector(); 390 Output->Boost(Beta); 149 Output->Boost(Beta); 391 << 150 /* // Uzhi >> 151 G4cout<<"Number of produced hadrons // correct E and P "<<Output->size()<<G4endl; // Uzhi >> 152 for(cHadron = 0; cHadron < Output->size(); cHadron++) >> 153 { >> 154 G4cout<<cHadron<<" "<<Output->operator[](cHadron)->Get4Momentum()<<" "<<Output->operator[](cHadron)->Get4Momentum().mag()<<Output->operator[](cHadron)->GetDefinition()->GetParticleName()<<G4endl; >> 155 } >> 156 */ // Uzhi 392 return success; 157 return success; 393 } << 158 } 394 << 395 159