Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // Historic fragment from M.Komogorov; clean-u 27 28 #include "G4ExcitedStringDecay.hh" 29 #include "G4SystemOfUnits.hh" 30 #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 39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V 40 : G4VStringFragmentation(), theStringDecay(p 41 { 42 if(!ptr) { 43 G4HadronicInteraction* p = 44 G4HadronicInteractionRegistry::Instance( 45 theStringDecay = static_cast<G4VLongitudin 46 if(!theStringDecay) { theStringDecay = new 47 } 48 SetModelName(theStringDecay->GetModelName()) 49 } 50 51 G4ExcitedStringDecay::~G4ExcitedStringDecay() 52 {} 53 54 G4KineticTrackVector *G4ExcitedStringDecay::Fr 55 { 56 return theStringDecay->FragmentString(theStr 57 } 58 59 G4KineticTrackVector *G4ExcitedStringDecay::Fr 60 { 61 G4LorentzVector KTsum(0.,0.,0.,0.); 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 266 if (!success) 267 { 268 if (theResult->size() != 0) 269 { 270 std::for_each(theResult->begin() , theRe 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 283 Ptmp=theStrings->operator[](astring)-> 284 Ptmp.transform( toLab); 285 theStrings->operator[](astring)->GetRi 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 } 297 298 299 G4bool G4ExcitedStringDecay:: 300 EnergyAndMomentumCorrector(G4KineticTrackVecto 301 { 302 const int nAttemptScale = 500; 303 const double ErrLimit = 1.E-5; 304 if (Output->empty()) return TRUE; 305 G4LorentzVector SumMom; 306 G4double SumMass = 0; 307 G4double TotalCollisionMass = Total 308 309 std::vector<G4double> HadronMass; G4double 310 311 #ifdef debug_G4ExcitedStringCorr 312 G4cout<<G4endl<<"EnergyAndMomentumCorrecto 313 #endif 314 // Calculate sum hadron 4-momenta and summ 315 unsigned int cHadron; 316 for (cHadron = 0; cHadron < Output->size() 317 { 318 SumMom += Output->operator[](cHadron) 319 HadronM=Output->operator[](cHadron)->G 320 SumMass += Output->operator[](cHadron) 321 } 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 330 if (Output->size() < 2) return FALSE; 331 332 if (SumMass > TotalCollisionMass) return F 333 SumMass = SumMom.m2(); 334 if (SumMass < 0) return FALSE; 335 SumMass = std::sqrt(SumMass); 336 337 // Compute c.m.s. hadron velocity and boos 338 // G4ThreeVector Beta = -SumMom.boostVecto 339 G4ThreeVector Beta = -TotalCollisionMom.bo 340 Output->Boost(Beta); 341 342 // Scale total c.m.s. hadron energy (hadro 343 // It should be equal interaction mass 344 G4double Scale = 1; 345 G4int cAttempt = 0; 346 G4double Sum = 0; 347 G4bool success = false; 348 for (cAttempt = 0; cAttempt < nAttemptScal 349 { 350 Sum = 0; 351 for (cHadron = 0; cHadron < Output->size 352 { 353 HadronM = HadronMass.at(cHadron); 354 G4LorentzVector HadronMom = Output->op 355 HadronMom.setVect(Scale*HadronMom.vect 356 G4double E = std::sqrt(HadronMom.vect( 357 358 HadronMom.setE(E); 359 Output->operator[](cHadron)->Set4Momen 360 Sum += E; 361 } 362 Scale = TotalCollisionMass/Sum; 363 #ifdef debug_G4ExcitedStringCorr 364 G4cout << "Scale-1=" << Scale -1 365 << ", TotalCollisionMass=" << To 366 << ", Sum=" << Sum 367 << G4endl; 368 #endif 369 if (std::fabs(Scale - 1) <= ErrLimit) 370 { 371 success = true; 372 break; 373 } 374 } 375 376 #ifdef debug_G4ExcitedStringCorr 377 if (!success) 378 { 379 G4cout << "G4ExcitedStringDecay::EnergyA 380 G4cout << " Scale not unity at end of 381 G4cout << " Number of secondaries: " < 382 G4cout << " Wanted total energy: " << 383 G4cout << " Increase number of attempt 384 // throw G4HadronicException(__FILE__, _ 385 } 386 #endif 387 388 // Compute c.m.s. interaction velocity and 389 Beta = TotalCollisionMom.boostVector(); 390 Output->Boost(Beta); 391 392 return success; 393 } 394 395