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 27 28 #include "G4ExcitedStringDecay.hh" 28 #include "G4ExcitedStringDecay.hh" 29 #include "G4SystemOfUnits.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4KineticTrack.hh" 30 #include "G4KineticTrack.hh" 31 #include "G4LundStringFragmentation.hh" << 32 #include "G4HadronicInteractionRegistry.hh" << 33 31 34 #include "G4SampleResonance.hh" 32 #include "G4SampleResonance.hh" 35 33 36 //#define debug_G4ExcitedStringDecay 34 //#define debug_G4ExcitedStringDecay 37 //#define debug_G4ExcitedStringCorr 35 //#define debug_G4ExcitedStringCorr 38 36 39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V << 37 40 : G4VStringFragmentation(), theStringDecay(p << 38 G4ExcitedStringDecay::G4ExcitedStringDecay() : G4VStringFragmentation(),theStringDecay(0) >> 39 {} >> 40 >> 41 >> 42 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay * aStringDecay) >> 43 : G4VStringFragmentation(), >> 44 theStringDecay(aStringDecay) >> 45 {} >> 46 >> 47 >> 48 G4ExcitedStringDecay::G4ExcitedStringDecay(const G4ExcitedStringDecay &) >> 49 : G4VStringFragmentation(), >> 50 theStringDecay(0) 41 { 51 { 42 if(!ptr) { << 52 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible"); 43 G4HadronicInteraction* p = << 53 } 44 G4HadronicInteractionRegistry::Instance( << 54 45 theStringDecay = static_cast<G4VLongitudin << 46 if(!theStringDecay) { theStringDecay = new << 47 } << 48 SetModelName(theStringDecay->GetModelName()) << 49 } << 50 55 51 G4ExcitedStringDecay::~G4ExcitedStringDecay() 56 G4ExcitedStringDecay::~G4ExcitedStringDecay() 52 {} << 57 { >> 58 } >> 59 >> 60 >> 61 const G4ExcitedStringDecay & G4ExcitedStringDecay::operator=(const G4ExcitedStringDecay &) >> 62 { >> 63 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable"); >> 64 return *this; >> 65 } >> 66 >> 67 >> 68 int G4ExcitedStringDecay::operator==(const G4ExcitedStringDecay &) const >> 69 { >> 70 return 0; >> 71 } >> 72 >> 73 >> 74 int G4ExcitedStringDecay::operator!=(const G4ExcitedStringDecay &) const >> 75 { >> 76 return 1; >> 77 } >> 78 53 79 54 G4KineticTrackVector *G4ExcitedStringDecay::Fr 80 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString(const G4ExcitedString &theString) 55 { 81 { 56 return theStringDecay->FragmentString(theStr << 82 if ( theStringDecay == NULL ) theStringDecay=new G4LundStringFragmentation(); >> 83 return theStringDecay->FragmentString(theString); 57 } 84 } 58 85 >> 86 59 G4KineticTrackVector *G4ExcitedStringDecay::Fr 87 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings(const G4ExcitedStringVector * theStrings) 60 { 88 { 61 G4LorentzVector KTsum(0.,0.,0.,0.); 89 G4LorentzVector KTsum(0.,0.,0.,0.); 62 90 63 #ifdef debug_G4ExcitedStringDecay 91 #ifdef debug_G4ExcitedStringDecay 64 G4cout<<G4endl; 92 G4cout<<G4endl; 65 G4cout<<"--------------------------- G4Excit 93 G4cout<<"--------------------------- G4ExcitedStringDecay ----------------------"<<G4endl; 66 G4cout<<"Hadronization of Excited Strings: t 94 G4cout<<"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<G4endl; 67 #endif 95 #endif 68 96 69 for ( unsigned int astring=0; astring < theS 97 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 70 // for ( unsigned int astring=0; astring < 1 98 // for ( unsigned int astring=0; astring < 1; astring++) 71 { 99 { 72 // G4cout<<"theStrings->operator[](astring 100 // G4cout<<"theStrings->operator[](astring)->IsExcited() "<<" "<<astring<<" "<<theStrings->operator[](astring)->IsExcited()<<G4endl; 73 if ( theStrings->operator[](astring)->IsEx 101 if ( theStrings->operator[](astring)->IsExcited() ) 74 {KTsum+= theStrings->operator[](astri 102 {KTsum+= theStrings->operator[](astring)->Get4Momentum();} 75 else {KTsum+=theStrings->operator[](astrin 103 else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();} 76 } 104 } 77 105 78 G4LorentzRotation toCms( -1 * KTsum.boostVec 106 G4LorentzRotation toCms( -1 * KTsum.boostVector() ); 79 G4LorentzRotation toLab(toCms.inverse()); 107 G4LorentzRotation toLab(toCms.inverse()); 80 G4LorentzVector Ptmp; 108 G4LorentzVector Ptmp; 81 KTsum=G4LorentzVector(0.,0.,0.,0.); 109 KTsum=G4LorentzVector(0.,0.,0.,0.); 82 110 83 for ( unsigned int astring=0; astring < theS 111 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 84 // for ( unsigned int astring=0; astring < 1 112 // for ( unsigned int astring=0; astring < 1; astring++) 85 { 113 { 86 if ( theStrings->operator[](astring)->IsEx 114 if ( theStrings->operator[](astring)->IsExcited() ) 87 { 115 { 88 Ptmp=toCms * theStrings->operator[](astr 116 Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum(); 89 theStrings->operator[](astring)->GetLeft 117 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp); 90 118 91 Ptmp=toCms * theStrings->operator[](astr 119 Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum(); 92 theStrings->operator[](astring)->GetRigh 120 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp); 93 121 94 KTsum+= theStrings->operator[](astring)- 122 KTsum+= theStrings->operator[](astring)->Get4Momentum(); 95 } 123 } 96 else 124 else 97 { 125 { 98 Ptmp=toCms * theStrings->operator[](astr 126 Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum(); 99 theStrings->operator[](astring)->GetKine 127 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp); 100 KTsum+= theStrings->operator[](astring)- 128 KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum(); 101 } 129 } 102 } 130 } 103 131 104 G4SampleResonance BrW; 132 G4SampleResonance BrW; 105 const G4ParticleDefinition* TrackDefinition= 133 const G4ParticleDefinition* TrackDefinition=0; 106 134 107 G4KineticTrackVector * theResult = new G4Kin 135 G4KineticTrackVector * theResult = new G4KineticTrackVector; 108 G4int attempts(0); 136 G4int attempts(0); 109 G4bool success=false; 137 G4bool success=false; 110 G4bool NeedEnergyCorrector=false; 138 G4bool NeedEnergyCorrector=false; 111 do { 139 do { 112 #ifdef debug_G4ExcitedStringDecay 140 #ifdef debug_G4ExcitedStringDecay 113 G4cout<<"New try No "<<attempts<<" to 141 G4cout<<"New try No "<<attempts<<" to hadronize strings"<<G4endl; 114 #endif 142 #endif 115 143 116 std::for_each(theResult->begin() , theResult 144 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack()); 117 theResult->clear(); 145 theResult->clear(); 118 146 119 attempts++; 147 attempts++; 120 148 121 G4LorentzVector KTsecondaries(0.,0.,0.,0.); 149 G4LorentzVector KTsecondaries(0.,0.,0.,0.); 122 NeedEnergyCorrector=false; 150 NeedEnergyCorrector=false; 123 151 124 for ( unsigned int astring=0; astring < theS 152 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 125 // for ( unsigned int astring=0; astri 153 // for ( unsigned int astring=0; astring < 1; astring++) // Proj. Last Str. Decay for FTF 126 // for ( unsigned int astring=1; astri 154 // for ( unsigned int astring=1; astring < 2; astring++) // Proj. Last Str. Decay for QGS 127 // for ( unsigned int astring=0; astri 155 // for ( unsigned int astring=0; astring < 1; astring++) // For testing purposes 128 { 156 { 129 #ifdef debug_G4ExcitedStringDecay 157 #ifdef debug_G4ExcitedStringDecay 130 G4cout<<"String No "<<astring+1<<" E 158 G4cout<<"String No "<<astring+1<<" Excited? "<<theStrings->operator[](astring)->IsExcited()<<G4endl; 131 159 132 G4cout<<"String No "<<astring+1<<" 4 160 G4cout<<"String No "<<astring+1<<" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum() 133 <<" "<<theStrings->operator[](astrin 161 <<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl; 134 #endif 162 #endif 135 163 136 G4KineticTrackVector * generatedKine 164 G4KineticTrackVector * generatedKineticTracks = NULL; 137 if ( theStrings->operator[](astring)->IsEx 165 if ( theStrings->operator[](astring)->IsExcited() ) 138 { 166 { 139 #ifdef debug_G4ExcitedStringDecay 167 #ifdef debug_G4ExcitedStringDecay 140 G4cout<<"Fragment String with par 168 G4cout<<"Fragment String with partons: " 141 <<theStrings->operator[](as 169 <<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode() <<" " 142 <<theStrings->operator[](as 170 <<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<" " 143 <<"Direction "<<theStrings- 171 <<"Direction "<<theStrings->operator[](astring)->GetDirection()<<G4endl; 144 #endif 172 #endif 145 generatedKineticTracks=FragmentString 173 generatedKineticTracks=FragmentString(*theStrings->operator[](astring)); 146 #ifdef debug_G4ExcitedStringDecay 174 #ifdef debug_G4ExcitedStringDecay 147 G4cout<<"(G4ExcitedStringDecay) N 175 G4cout<<"(G4ExcitedStringDecay) Number of produced hadrons = " 148 <<generatedKineticTracks->s 176 <<generatedKineticTracks->size()<<G4endl; 149 #endif 177 #endif 150 } else { 178 } else { 151 #ifdef debug_G4ExcitedStringDecay 179 #ifdef debug_G4ExcitedStringDecay 152 G4cout<<" GetTrack from the Str 180 G4cout<<" GetTrack from the String"<<G4endl; 153 #endif 181 #endif 154 G4LorentzVector Mom=theStrings->o 182 G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum(); 155 G4KineticTrack * aTrack= new G4Ki 183 G4KineticTrack * aTrack= new G4KineticTrack( 156 theStrings->operat 184 theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(), 157 theStrings->operat 185 theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(), 158 G4ThreeVector(0), 186 G4ThreeVector(0), Mom); 159 187 160 aTrack->SetPosition(theStrings->o 188 aTrack->SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition()); 161 189 162 #ifdef debug_G4ExcitedStringDecay 190 #ifdef debug_G4ExcitedStringDecay 163 G4cout<<" A particle stored in 191 G4cout<<" A particle stored in the track is "<<aTrack->GetDefinition()->GetParticleName()<<G4endl; 164 #endif 192 #endif 165 193 166 generatedKineticTracks = new G4KineticT 194 generatedKineticTracks = new G4KineticTrackVector; 167 generatedKineticTracks->push_back(aTrac 195 generatedKineticTracks->push_back(aTrack); 168 } 196 } 169 197 170 if (generatedKineticTracks == nullptr || g << 198 if (generatedKineticTracks->size() == 0) 171 { 199 { 172 // G4cerr << "G4VPartonStringMode 200 // G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl; 173 delete generatedKineticTracks; << 201 // continue; 174 generatedKineticTracks = nullptr; << 202 success=false; NeedEnergyCorrector=false; break; 175 continue; << 176 } 203 } 177 204 178 G4LorentzVector KTsum1(0.,0.,0.,0.); 205 G4LorentzVector KTsum1(0.,0.,0.,0.); 179 for ( unsigned int aTrack=0; aTrack< 206 for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++) 180 { 207 { 181 #ifdef debug_G4ExcitedStringDecay 208 #ifdef debug_G4ExcitedStringDecay 182 G4cout<<"Prod part No. "<<aTrack+ 209 G4cout<<"Prod part No. "<<aTrack+1<<" " 183 <<(*generatedKineticTracks) 210 <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" " 184 <<(*generatedKineticTracks) 211 <<(*generatedKineticTracks)[aTrack]->Get4Momentum() 185 <<(*generatedKineticTracks) 212 <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl; 186 #endif 213 #endif 187 // --------------- Sampling mass 214 // --------------- Sampling mass of unstable hadronic resonances ---------------- 188 TrackDefinition = (*generatedKine 215 TrackDefinition = (*generatedKineticTracks)[aTrack]->GetDefinition(); 189 216 190 if (TrackDefinition->IsShortLived()) << 217 if(TrackDefinition->IsShortLived()) 191 { 218 { 192 G4double NewTrackMass = << 219 G4double NewTrackMass = BrW.SampleMass( TrackDefinition, 193 BrW.SampleMass( TrackDefiniti << 220 TrackDefinition->GetPDGMass() + 5.0*TrackDefinition->GetPDGWidth() ); 194 BrW.GetMinimu << 195 TrackDefiniti << 196 G4LorentzVector Tmp=G4LorentzVe 221 G4LorentzVector Tmp=G4LorentzVector((*generatedKineticTracks)[aTrack]->Get4Momentum()); 197 Tmp.setE(std::sqrt(sqr(NewTrack 222 Tmp.setE(std::sqrt(sqr(NewTrackMass) + Tmp.vect().mag2())); 198 223 199 (*generatedKineticTracks)[aTrac 224 (*generatedKineticTracks)[aTrack]->Set4Momentum(Tmp); 200 225 201 #ifdef debug_G4ExcitedStringDec 226 #ifdef debug_G4ExcitedStringDecay 202 G4cout<<"Resonance *** "<<aTrac 227 G4cout<<"Resonance *** "<<aTrack+1<<" " 203 <<(*generatedKineticTrack 228 <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" " 204 <<(*generatedKineticTrack 229 <<(*generatedKineticTracks)[aTrack]->Get4Momentum() 205 <<(*generatedKineticTrack 230 <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl; 206 #endif 231 #endif 207 } 232 } 208 //------------------------------- 233 //------------------------------------------------------------------------------- 209 234 210 theResult->push_back(generatedKin 235 theResult->push_back(generatedKineticTracks->operator[](aTrack)); 211 KTsum1+= (*generatedKineticTracks 236 KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum(); 212 } 237 } 213 KTsecondaries+=KTsum1; 238 KTsecondaries+=KTsum1; 214 239 215 #ifdef debug_G4ExcitedStringDecay 240 #ifdef debug_G4ExcitedStringDecay 216 G4cout << "String secondaries(" <<ge 241 G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ")"<<G4endl 217 <<"Init string momentum: "< 242 <<"Init string momentum: "<< theStrings->operator[](astring)->Get4Momentum()<<G4endl 218 <<"Final hadrons momentum: "< 243 <<"Final hadrons momentum: "<< KTsum1 << G4endl; 219 #endif 244 #endif 220 245 221 if ( KTsum1.e() > 0 && << 246 if ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion ) 222 std::abs((KTsum1.e()-theString << 223 { 247 { 224 NeedEnergyCorrector=true; 248 NeedEnergyCorrector=true; 225 } 249 } 226 250 227 #ifdef debug_G4ExcitedStringDecay 251 #ifdef debug_G4ExcitedStringDecay 228 G4cout<<"NeedEnergyCorrection yes/no 252 G4cout<<"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<G4endl; 229 #endif 253 #endif 230 254 231 // clean up 255 // clean up 232 delete generatedKineticTracks; 256 delete generatedKineticTracks; 233 success=true; 257 success=true; 234 } 258 } 235 259 236 if ( NeedEnergyCorrector ) success=EnergyAnd 260 if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum); 237 } while (!success && (attempts < 100)); /* 261 } while (!success && (attempts < 100)); /* Loop checking, 07.08.2015, A.Ribon */ 238 262 239 for ( unsigned int aTrack=0; aTrack<theResul 263 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++) 240 { 264 { 241 Ptmp=(*theResult)[aTrack]->Get4Momentum(); 265 Ptmp=(*theResult)[aTrack]->Get4Momentum(); 242 Ptmp.transform( toLab); 266 Ptmp.transform( toLab); 243 (*theResult)[aTrack]->Set4Momentum(Ptmp); 267 (*theResult)[aTrack]->Set4Momentum(Ptmp); 244 } 268 } 245 269 246 #ifdef debug_G4ExcitedStringDecay 270 #ifdef debug_G4ExcitedStringDecay 247 G4cout<<"End of the strings fragmentation (G 271 G4cout<<"End of the strings fragmentation (G4ExcitedStringDecay)"<<G4endl; 248 272 249 G4LorentzVector KTsum1(0.,0.,0.,0.); 273 G4LorentzVector KTsum1(0.,0.,0.,0.); 250 274 251 for ( unsigned int aTrack=0; aTrack<theResul 275 for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++) 252 { 276 { 253 G4cout << " corrected tracks .. " << (*t 277 G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName() 254 <<" " << (*theResult)[aTrack]->Get4Mome 278 <<" " << (*theResult)[aTrack]->Get4Momentum() 255 <<" " << (*theResult)[aTrack]->Get4Mome 279 <<" " << (*theResult)[aTrack]->Get4Momentum().mag()<< G4endl; 256 KTsum1+= (*theResult)[aTrack]->Get4Momen 280 KTsum1+= (*theResult)[aTrack]->Get4Momentum(); 257 } 281 } 258 282 259 G4cout << "Needcorrector/success " << NeedEn 283 G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success 260 << ", Corrected total 4 momentum " < 284 << ", Corrected total 4 momentum " << KTsum1 << G4endl; 261 if ( ! success ) G4cout << "failed to correc 285 if ( ! success ) G4cout << "failed to correct E/p" << G4endl; 262 286 263 G4cout<<"End of the Hadronization (G4Excited 287 G4cout<<"End of the Hadronization (G4ExcitedStringDecay)"<<G4endl; 264 #endif 288 #endif 265 289 266 if (!success) << 290 if(!success) 267 { 291 { 268 if (theResult->size() != 0) << 292 if(theResult->size() != 0) 269 { 293 { 270 std::for_each(theResult->begin() , theRe 294 std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack()); 271 theResult->clear(); 295 theResult->clear(); 272 delete theResult; theResult=0; 296 delete theResult; theResult=0; 273 } 297 } 274 for ( unsigned int astring=0; astring < th 298 for ( unsigned int astring=0; astring < theStrings->size(); astring++) 275 // for ( unsigned int astring=0; astring < 299 // for ( unsigned int astring=0; astring < 1; astring++) // Need more correct. For testing purposes. 276 { 300 { 277 if ( theStrings->operator[](astring)->Is 301 if ( theStrings->operator[](astring)->IsExcited() ) 278 { 302 { 279 Ptmp=theStrings->operator[](astring)-> 303 Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum(); 280 Ptmp.transform( toLab); 304 Ptmp.transform( toLab); 281 theStrings->operator[](astring)->GetLe 305 theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp); 282 306 283 Ptmp=theStrings->operator[](astring)-> 307 Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum(); 284 Ptmp.transform( toLab); 308 Ptmp.transform( toLab); 285 theStrings->operator[](astring)->GetRi 309 theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp); 286 } 310 } 287 else 311 else 288 { 312 { 289 Ptmp=theStrings->operator[](astring)-> 313 Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum(); 290 Ptmp.transform( toLab); 314 Ptmp.transform( toLab); 291 theStrings->operator[](astring)->GetKi 315 theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp); 292 } 316 } 293 } 317 } 294 } 318 } 295 return theResult; 319 return theResult; 296 } 320 } 297 321 298 322 299 G4bool G4ExcitedStringDecay:: << 323 G4bool G4ExcitedStringDecay::EnergyAndMomentumCorrector 300 EnergyAndMomentumCorrector(G4KineticTrackVecto << 324 (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom) 301 { 325 { 302 const int nAttemptScale = 500; 326 const int nAttemptScale = 500; 303 const double ErrLimit = 1.E-5; 327 const double ErrLimit = 1.E-5; 304 if (Output->empty()) return TRUE; 328 if (Output->empty()) return TRUE; 305 G4LorentzVector SumMom; 329 G4LorentzVector SumMom; 306 G4double SumMass = 0; 330 G4double SumMass = 0; 307 G4double TotalCollisionMass = Total 331 G4double TotalCollisionMass = TotalCollisionMom.m(); 308 332 309 std::vector<G4double> HadronMass; G4double 333 std::vector<G4double> HadronMass; G4double HadronM(0.); 310 334 311 #ifdef debug_G4ExcitedStringCorr 335 #ifdef debug_G4ExcitedStringCorr 312 G4cout<<G4endl<<"EnergyAndMomentumCorrecto 336 G4cout<<G4endl<<"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<G4endl; 313 #endif 337 #endif 314 // Calculate sum hadron 4-momenta and summ 338 // Calculate sum hadron 4-momenta and summing hadron mass 315 unsigned int cHadron; 339 unsigned int cHadron; 316 for (cHadron = 0; cHadron < Output->size() << 340 for(cHadron = 0; cHadron < Output->size(); cHadron++) 317 { 341 { 318 SumMom += Output->operator[](cHadron) 342 SumMom += Output->operator[](cHadron)->Get4Momentum(); 319 HadronM=Output->operator[](cHadron)->G 343 HadronM=Output->operator[](cHadron)->Get4Momentum().mag(); HadronMass.push_back(HadronM); 320 SumMass += Output->operator[](cHadron) 344 SumMass += Output->operator[](cHadron)->Get4Momentum().mag(); //GetDefinition()->GetPDGMass(); 321 } 345 } 322 346 323 #ifdef debug_G4ExcitedStringCorr 347 #ifdef debug_G4ExcitedStringCorr 324 G4cout<<"Sum part mom "<<SumMom<<" "<<SumM 348 G4cout<<"Sum part mom "<<SumMom<<" "<<SumMom.mag()<<G4endl 325 <<"Sum str mom "<<TotalCollisionMom 349 <<"Sum str mom "<<TotalCollisionMom<<" "<<TotalCollisionMom.mag()<<G4endl; 326 G4cout<<"SumMass TotalCollisionMass "<<Sum 350 G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl; 327 #endif 351 #endif 328 352 329 // Cannot correct a single particle 353 // Cannot correct a single particle 330 if (Output->size() < 2) return FALSE; 354 if (Output->size() < 2) return FALSE; 331 355 332 if (SumMass > TotalCollisionMass) return F 356 if (SumMass > TotalCollisionMass) return FALSE; 333 SumMass = SumMom.m2(); 357 SumMass = SumMom.m2(); 334 if (SumMass < 0) return FALSE; 358 if (SumMass < 0) return FALSE; 335 SumMass = std::sqrt(SumMass); 359 SumMass = std::sqrt(SumMass); 336 360 337 // Compute c.m.s. hadron velocity and boos 361 // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s. 338 // G4ThreeVector Beta = -SumMom.boostVecto 362 // G4ThreeVector Beta = -SumMom.boostVector(); 339 G4ThreeVector Beta = -TotalCollisionMom.bo 363 G4ThreeVector Beta = -TotalCollisionMom.boostVector(); 340 Output->Boost(Beta); 364 Output->Boost(Beta); 341 365 342 // Scale total c.m.s. hadron energy (hadro 366 // Scale total c.m.s. hadron energy (hadron system mass). 343 // It should be equal interaction mass 367 // It should be equal interaction mass 344 G4double Scale = 1; 368 G4double Scale = 1; 345 G4int cAttempt = 0; 369 G4int cAttempt = 0; 346 G4double Sum = 0; 370 G4double Sum = 0; 347 G4bool success = false; 371 G4bool success = false; 348 for (cAttempt = 0; cAttempt < nAttemptScal << 372 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++) 349 { 373 { 350 Sum = 0; 374 Sum = 0; 351 for (cHadron = 0; cHadron < Output->size << 375 for(cHadron = 0; cHadron < Output->size(); cHadron++) 352 { 376 { 353 HadronM = HadronMass.at(cHadron); 377 HadronM = HadronMass.at(cHadron); 354 G4LorentzVector HadronMom = Output->op 378 G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum(); 355 HadronMom.setVect(Scale*HadronMom.vect 379 HadronMom.setVect(Scale*HadronMom.vect()); 356 G4double E = std::sqrt(HadronMom.vect( 380 G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(HadronM)); 357 381 //sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass())); 358 HadronMom.setE(E); 382 HadronMom.setE(E); 359 Output->operator[](cHadron)->Set4Momen 383 Output->operator[](cHadron)->Set4Momentum(HadronMom); 360 Sum += E; 384 Sum += E; 361 } 385 } 362 Scale = TotalCollisionMass/Sum; 386 Scale = TotalCollisionMass/Sum; 363 #ifdef debug_G4ExcitedStringCorr 387 #ifdef debug_G4ExcitedStringCorr 364 G4cout << "Scale-1=" << Scale -1 388 G4cout << "Scale-1=" << Scale -1 365 << ", TotalCollisionMass=" << To 389 << ", TotalCollisionMass=" << TotalCollisionMass 366 << ", Sum=" << Sum 390 << ", Sum=" << Sum 367 << G4endl; 391 << G4endl; 368 #endif 392 #endif 369 if (std::fabs(Scale - 1) <= ErrLimit) 393 if (std::fabs(Scale - 1) <= ErrLimit) 370 { 394 { 371 success = true; 395 success = true; 372 break; 396 break; 373 } 397 } 374 } 398 } 375 399 376 #ifdef debug_G4ExcitedStringCorr 400 #ifdef debug_G4ExcitedStringCorr 377 if (!success) << 401 if(!success) 378 { 402 { 379 G4cout << "G4ExcitedStringDecay::EnergyA 403 G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl; 380 G4cout << " Scale not unity at end of 404 G4cout << " Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl; 381 G4cout << " Number of secondaries: " < 405 G4cout << " Number of secondaries: " << Output->size() << G4endl; 382 G4cout << " Wanted total energy: " << 406 G4cout << " Wanted total energy: " << TotalCollisionMom.e() << G4endl; 383 G4cout << " Increase number of attempt 407 G4cout << " Increase number of attempts or increase ERRLIMIT"<<G4endl; 384 // throw G4HadronicException(__FILE__, _ 408 // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct..."); 385 } 409 } 386 #endif 410 #endif 387 411 388 // Compute c.m.s. interaction velocity and 412 // Compute c.m.s. interaction velocity and KTV back boost 389 Beta = TotalCollisionMom.boostVector(); 413 Beta = TotalCollisionMom.boostVector(); 390 Output->Boost(Beta); 414 Output->Boost(Beta); 391 415 392 return success; 416 return success; 393 } 417 } 394 << 395 418