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 // G4QMDMeanField implementation << 26 // 081120 Add Update by T. Koi 27 // 27 // 28 // Author: Tatsumi Koi (SLAC), 29 March 2007 << 29 // ------------------------------------------- << 30 28 31 #include <map> 29 #include <map> 32 #include <algorithm> 30 #include <algorithm> 33 #include <numeric> 31 #include <numeric> 34 32 35 #include <cmath> << 36 #include <CLHEP/Random/Stat.h> 33 #include <CLHEP/Random/Stat.h> 37 34 38 #include "G4QMDMeanField.hh" 35 #include "G4QMDMeanField.hh" 39 #include "G4QMDParameters.hh" 36 #include "G4QMDParameters.hh" 40 #include "G4Exp.hh" << 37 41 #include "G4Pow.hh" << 42 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh" 43 #include "Randomize.hh" 39 #include "Randomize.hh" 44 40 45 G4QMDMeanField::G4QMDMeanField() 41 G4QMDMeanField::G4QMDMeanField() >> 42 : rclds ( 4.0 ) // distance for cluster judgement >> 43 , epsx ( -20.0 ) // gauss term >> 44 , epscl ( 0.0001 ) // coulomb term >> 45 , irelcr ( 1 ) 46 { 46 { >> 47 47 G4QMDParameters* parameters = G4QMDParamete 48 G4QMDParameters* parameters = G4QMDParameters::GetInstance(); 48 wl = parameters->Get_wl(); 49 wl = parameters->Get_wl(); 49 cl = parameters->Get_cl(); 50 cl = parameters->Get_cl(); 50 rho0 = parameters->Get_rho0(); 51 rho0 = parameters->Get_rho0(); 51 hbc = parameters->Get_hbc(); 52 hbc = parameters->Get_hbc(); 52 gamm = parameters->Get_gamm(); 53 gamm = parameters->Get_gamm(); 53 54 54 cpw = parameters->Get_cpw(); 55 cpw = parameters->Get_cpw(); 55 cph = parameters->Get_cph(); 56 cph = parameters->Get_cph(); 56 cpc = parameters->Get_cpc(); 57 cpc = parameters->Get_cpc(); 57 58 58 c0 = parameters->Get_c0(); 59 c0 = parameters->Get_c0(); 59 c3 = parameters->Get_c3(); 60 c3 = parameters->Get_c3(); 60 cs = parameters->Get_cs(); 61 cs = parameters->Get_cs(); 61 62 62 // distance << 63 // distance 63 c0w = 1.0/4.0/wl; 64 c0w = 1.0/4.0/wl; >> 65 //c3w = 1.0/4.0/wl; //no need 64 c0sw = std::sqrt( c0w ); 66 c0sw = std::sqrt( c0w ); 65 clw = 2.0 / std::sqrt ( 4.0 * pi * wl ); 67 clw = 2.0 / std::sqrt ( 4.0 * pi * wl ); 66 68 67 // graduate << 69 // graduate 68 c0g = - c0 / ( 2.0 * wl ); 70 c0g = - c0 / ( 2.0 * wl ); 69 c3g = - c3 / ( 4.0 * wl ) * gamm; 71 c3g = - c3 / ( 4.0 * wl ) * gamm; 70 csg = - cs / ( 2.0 * wl ); 72 csg = - cs / ( 2.0 * wl ); 71 pag = gamm - 1; 73 pag = gamm - 1; 72 74 73 system = nullptr; // will be set through Se << 74 } 75 } 75 76 >> 77 >> 78 >> 79 G4QMDMeanField::~G4QMDMeanField() >> 80 { >> 81 ; >> 82 } >> 83 >> 84 >> 85 76 void G4QMDMeanField::SetSystem ( G4QMDSystem* 86 void G4QMDMeanField::SetSystem ( G4QMDSystem* aSystem ) 77 { 87 { >> 88 >> 89 //std::cout << "QMDMeanField SetSystem" << std::endl; >> 90 78 system = aSystem; 91 system = aSystem; 79 92 80 G4int n = system->GetTotalNumberOfParticipa 93 G4int n = system->GetTotalNumberOfParticipant(); 81 94 82 pp2.clear(); 95 pp2.clear(); 83 rr2.clear(); 96 rr2.clear(); 84 rbij.clear(); 97 rbij.clear(); 85 rha.clear(); 98 rha.clear(); 86 rhe.clear(); 99 rhe.clear(); 87 rhc.clear(); 100 rhc.clear(); 88 101 89 rr2.resize( n ); 102 rr2.resize( n ); 90 pp2.resize( n ); 103 pp2.resize( n ); 91 rbij.resize( n ); 104 rbij.resize( n ); 92 rha.resize( n ); 105 rha.resize( n ); 93 rhe.resize( n ); 106 rhe.resize( n ); 94 rhc.resize( n ); 107 rhc.resize( n ); 95 108 96 for ( G4int i = 0 ; i < n ; ++i ) << 109 for ( int i = 0 ; i < n ; i++ ) 97 { 110 { 98 rr2[i].resize( n ); 111 rr2[i].resize( n ); 99 pp2[i].resize( n ); 112 pp2[i].resize( n ); 100 rbij[i].resize( n ); 113 rbij[i].resize( n ); 101 rha[i].resize( n ); 114 rha[i].resize( n ); 102 rhe[i].resize( n ); 115 rhe[i].resize( n ); 103 rhc[i].resize( n ); 116 rhc[i].resize( n ); 104 } 117 } 105 118 >> 119 106 ffr.clear(); 120 ffr.clear(); 107 ffp.clear(); 121 ffp.clear(); 108 rh3d.clear(); 122 rh3d.clear(); 109 123 110 ffr.resize( n ); 124 ffr.resize( n ); 111 ffp.resize( n ); 125 ffp.resize( n ); 112 rh3d.resize( n ); 126 rh3d.resize( n ); 113 127 114 Cal2BodyQuantities(); 128 Cal2BodyQuantities(); >> 129 115 } 130 } 116 131 117 void G4QMDMeanField::SetNucleus ( G4QMDNucleus 132 void G4QMDMeanField::SetNucleus ( G4QMDNucleus* aNucleus ) 118 { 133 { >> 134 >> 135 //std::cout << "QMDMeanField SetNucleus" << std::endl; >> 136 119 SetSystem( aNucleus ); 137 SetSystem( aNucleus ); 120 138 121 G4double totalPotential = GetTotalPotential 139 G4double totalPotential = GetTotalPotential(); 122 aNucleus->SetTotalPotential( totalPotential 140 aNucleus->SetTotalPotential( totalPotential ); >> 141 123 aNucleus->CalEnergyAndAngularMomentumInCM() 142 aNucleus->CalEnergyAndAngularMomentumInCM(); >> 143 124 } 144 } 125 145 >> 146 >> 147 126 void G4QMDMeanField::Cal2BodyQuantities() 148 void G4QMDMeanField::Cal2BodyQuantities() 127 { 149 { 128 if ( system->GetTotalNumberOfParticipant() << 129 150 130 for ( G4int j = 1 ; j < system->GetTotalNum << 151 if ( system->GetTotalNumberOfParticipant() < 2 ) return; >> 152 >> 153 for ( G4int j = 1 ; j < system->GetTotalNumberOfParticipant() ; j++ ) 131 { 154 { >> 155 132 G4ThreeVector rj = system->GetParticipan 156 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 133 G4LorentzVector p4j = system->GetPartici 157 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 134 158 135 for ( G4int i = 0 ; i < j ; ++i ) << 159 for ( G4int i = 0 ; i < j ; i++ ) 136 { 160 { >> 161 137 G4ThreeVector ri = system->GetPartici 162 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 138 G4LorentzVector p4i = system->GetPart 163 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 139 164 140 G4ThreeVector rij = ri - rj; 165 G4ThreeVector rij = ri - rj; 141 G4ThreeVector pij = (p4i - p4j).v(); 166 G4ThreeVector pij = (p4i - p4j).v(); 142 G4LorentzVector p4ij = p4i - p4j; 167 G4LorentzVector p4ij = p4i - p4j; 143 G4ThreeVector bij = ( p4i + p4j ).boo 168 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 144 G4double gammaij = ( p4i + p4j ).gamm 169 G4double gammaij = ( p4i + p4j ).gamma(); 145 170 146 G4double eij = ( p4i + p4j ).e(); 171 G4double eij = ( p4i + p4j ).e(); 147 172 148 G4double rbrb = rij*bij; 173 G4double rbrb = rij*bij; >> 174 // G4double bij2 = bij*bij; 149 G4double rij2 = rij*rij; 175 G4double rij2 = rij*rij; 150 G4double pij2 = pij*pij; 176 G4double pij2 = pij*pij; 151 177 152 rbrb = irelcr * rbrb; 178 rbrb = irelcr * rbrb; 153 G4double gamma2_ij = gammaij*gammaij 179 G4double gamma2_ij = gammaij*gammaij; 154 180 >> 181 155 rr2[i][j] = rij2 + gamma2_ij * rbrb*r 182 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 156 rr2[j][i] = rr2[i][j]; 183 rr2[j][i] = rr2[i][j]; 157 184 158 rbij[i][j] = gamma2_ij * rbrb; 185 rbij[i][j] = gamma2_ij * rbrb; 159 rbij[j][i] = - rbij[i][j]; 186 rbij[j][i] = - rbij[i][j]; 160 187 161 pp2[i][j] = pij2 188 pp2[i][j] = pij2 162 + irelcr * ( - G4Pow::GetIn << 189 + irelcr * ( - std::pow ( p4i.e() - p4j.e() , 2 ) 163 + gamma2_ij * G4Pow::GetIns << 190 + gamma2_ij * std::pow ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 164 << 191 >> 192 165 pp2[j][i] = pp2[i][j]; 193 pp2[j][i] = pp2[i][j]; 166 194 167 // Gauss term << 195 // Gauss term 168 196 169 G4double expa1 = - rr2[i][j] * c0w; 197 G4double expa1 = - rr2[i][j] * c0w; 170 198 171 G4double rh1; 199 G4double rh1; 172 if ( expa1 > epsx ) 200 if ( expa1 > epsx ) 173 { 201 { 174 rh1 = G4Exp( expa1 ); << 202 rh1 = std::exp( expa1 ); 175 } 203 } 176 else 204 else 177 { 205 { 178 rh1 = 0.0; 206 rh1 = 0.0; 179 } 207 } 180 208 181 G4int ibry = system->GetParticipant(i 209 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 182 G4int jbry = system->GetParticipant(j 210 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 183 211 >> 212 184 rha[i][j] = ibry*jbry*rh1; 213 rha[i][j] = ibry*jbry*rh1; 185 rha[j][i] = rha[i][j]; 214 rha[j][i] = rha[i][j]; 186 215 187 // Coulomb terms << 216 // Coulomb terms 188 217 189 G4double rrs2 = rr2[i][j] + epscl; 218 G4double rrs2 = rr2[i][j] + epscl; 190 G4double rrs = std::sqrt ( rrs2 ); 219 G4double rrs = std::sqrt ( rrs2 ); 191 220 192 G4int icharge = system->GetParticipan 221 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 193 G4int jcharge = system->GetParticipan 222 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 194 223 195 G4double xerf = 0.0; << 224 G4double erf = 0.0; 196 // T. K. add this protection. 5.8 is 225 // T. K. add this protection. 5.8 is good enough for double 197 if ( rrs*c0sw < 5.8 ) 226 if ( rrs*c0sw < 5.8 ) 198 { << 227 erf = CLHEP::HepStat::erf ( rrs*c0sw ); 199 #if defined WIN32-VC << 200 xerf = CLHEP::HepStat::erf ( rrs*c << 201 #else << 202 xerf = std::erf ( rrs*c0sw ); << 203 #endif << 204 } << 205 else 228 else 206 { << 229 erf = 1.0; 207 xerf = 1.0; << 230 208 } << 231 G4double erfij = erf/rrs; 209 232 210 G4double erfij = xerf/rrs; << 211 233 212 rhe[i][j] = icharge*jcharge * erfij; 234 rhe[i][j] = icharge*jcharge * erfij; >> 235 213 rhe[j][i] = rhe[i][j]; 236 rhe[j][i] = rhe[i][j]; >> 237 214 rhc[i][j] = icharge*jcharge * ( - erf 238 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; >> 239 215 rhc[j][i] = rhc[i][j]; 240 rhc[j][i] = rhc[i][j]; >> 241 216 } // i 242 } // i 217 } // j 243 } // j 218 } 244 } 219 245 >> 246 >> 247 220 void G4QMDMeanField::Cal2BodyQuantities( G4int 248 void G4QMDMeanField::Cal2BodyQuantities( G4int i ) 221 { 249 { >> 250 >> 251 //std::cout << "Cal2BodyQuantities " << i << std::endl; >> 252 222 G4ThreeVector ri = system->GetParticipant( 253 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 223 G4LorentzVector p4i = system->GetParticipan 254 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 224 255 225 for ( G4int j = 0 ; j < system->GetTotalNum << 256 >> 257 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) 226 { 258 { 227 if ( j == i ) { continue; } << 259 if ( j == i ) continue; 228 260 229 G4ThreeVector rj = system->GetParticipan 261 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 230 G4LorentzVector p4j = system->GetPartici 262 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 231 263 232 G4ThreeVector rij = ri - rj; << 264 G4ThreeVector rij = ri - rj; 233 G4ThreeVector pij = (p4i - p4j).v(); << 265 G4ThreeVector pij = (p4i - p4j).v(); 234 G4LorentzVector p4ij = p4i - p4j; << 266 G4LorentzVector p4ij = p4i - p4j; 235 G4ThreeVector bij = ( p4i + p4j ).boostV << 267 G4ThreeVector bij = ( p4i + p4j ).boostVector(); 236 G4double gammaij = ( p4i + p4j ).gamma() << 268 G4double gammaij = ( p4i + p4j ).gamma(); 237 << 269 238 G4double eij = ( p4i + p4j ).e(); << 270 G4double eij = ( p4i + p4j ).e(); 239 << 271 240 G4double rbrb = rij*bij; << 272 G4double rbrb = rij*bij; 241 G4double rij2 = rij*rij; << 273 // G4double bij2 = bij*bij; 242 G4double pij2 = pij*pij; << 274 G4double rij2 = rij*rij; >> 275 G4double pij2 = pij*pij; >> 276 >> 277 rbrb = irelcr * rbrb; >> 278 G4double gamma2_ij = gammaij*gammaij; >> 279 >> 280 /* >> 281 G4double rbrb = 0.0; >> 282 G4double beta2_ij = 0.0; >> 283 G4double rij2 = 0.0; >> 284 G4double pij2 = 0.0; >> 285 >> 286 // >> 287 G4LorentzVector p4ip4j = p4i + p4j; >> 288 G4double eij = p4ip4j.e(); >> 289 >> 290 G4ThreeVector r = ri - rj; >> 291 G4LorentzVector p4 = p4i - p4j; >> 292 >> 293 rbrb = r.x()*p4ip4j.x()/eij >> 294 + r.y()*p4ip4j.y()/eij >> 295 + r.z()*p4ip4j.z()/eij; >> 296 >> 297 beta2_ij = ( p4ip4j.x()*p4ip4j.x() + p4ip4j.y()*p4ip4j.y() + p4ip4j.z()*p4ip4j.z() ) / ( eij*eij ); >> 298 rij2 = r*r; >> 299 pij2 = p4.v()*p4.v(); 243 300 244 rbrb = irelcr * rbrb; 301 rbrb = irelcr * rbrb; 245 G4double gamma2_ij = gammaij*gammaij; << 302 >> 303 G4double gamma2_ij = 1 / ( 1 - beta2_ij ); >> 304 */ 246 305 247 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb 306 rr2[i][j] = rij2 + gamma2_ij * rbrb*rbrb; 248 rr2[j][i] = rr2[i][j]; 307 rr2[j][i] = rr2[i][j]; 249 308 250 rbij[i][j] = gamma2_ij * rbrb; 309 rbij[i][j] = gamma2_ij * rbrb; 251 rbij[j][i] = - rbij[i][j]; 310 rbij[j][i] = - rbij[i][j]; 252 311 253 pp2[i][j] = pij2 312 pp2[i][j] = pij2 254 + irelcr * ( - G4Pow::GetInsta << 313 + irelcr * ( - std::pow ( p4i.e() - p4j.e() , 2 ) 255 + gamma2_ij * G4Pow::GetInstan << 314 + gamma2_ij * std::pow ( ( ( p4i.m2() - p4j.m2() ) / eij ) , 2 ) ); 256 << 315 257 pp2[j][i] = pp2[i][j]; 316 pp2[j][i] = pp2[i][j]; 258 317 259 // Gauss term << 318 // Gauss term 260 319 261 G4double expa1 = - rr2[i][j] * c0w; 320 G4double expa1 = - rr2[i][j] * c0w; 262 321 263 G4double rh1; 322 G4double rh1; 264 if ( expa1 > epsx ) 323 if ( expa1 > epsx ) 265 { 324 { 266 rh1 = G4Exp( expa1 ); << 325 rh1 = std::exp( expa1 ); 267 } 326 } 268 else 327 else 269 { 328 { 270 rh1 = 0.0; 329 rh1 = 0.0; 271 } 330 } 272 331 273 G4int ibry = system->GetParticipant(i)-> 332 G4int ibry = system->GetParticipant(i)->GetBaryonNumber(); 274 G4int jbry = system->GetParticipant(j)-> 333 G4int jbry = system->GetParticipant(j)->GetBaryonNumber(); 275 334 >> 335 276 rha[i][j] = ibry*jbry*rh1; 336 rha[i][j] = ibry*jbry*rh1; 277 rha[j][i] = rha[i][j]; 337 rha[j][i] = rha[i][j]; 278 338 279 // Coulomb terms << 339 // Coulomb terms 280 340 281 G4double rrs2 = rr2[i][j] + epscl; 341 G4double rrs2 = rr2[i][j] + epscl; 282 G4double rrs = std::sqrt ( rrs2 ); 342 G4double rrs = std::sqrt ( rrs2 ); 283 343 284 G4int icharge = system->GetParticipant(i 344 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 285 G4int jcharge = system->GetParticipant(j 345 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 286 346 287 G4double xerf = 0.0; << 347 G4double erf = 0.0; 288 // T. K. add this protection. 5.8 is goo 348 // T. K. add this protection. 5.8 is good enough for double 289 if ( rrs*c0sw < 5.8 ) 349 if ( rrs*c0sw < 5.8 ) 290 { << 350 erf = CLHEP::HepStat::erf ( rrs*c0sw ); 291 #if defined WIN32-VC << 292 xerf = CLHEP::HepStat::erf ( rrs*c0sw << 293 #else << 294 xerf = std::erf ( rrs*c0sw ); << 295 #endif << 296 } << 297 else 351 else 298 { << 352 erf = 1.0; 299 xerf = 1.0; << 300 } << 301 353 302 G4double erfij = xerf/rrs; << 354 G4double erfij = erf/rrs; 303 355 >> 356 304 rhe[i][j] = icharge*jcharge * erfij; 357 rhe[i][j] = icharge*jcharge * erfij; >> 358 305 rhe[j][i] = rhe[i][j]; 359 rhe[j][i] = rhe[i][j]; >> 360 >> 361 // G4double clw; >> 362 306 rhc[i][j] = icharge*jcharge * ( - erfij 363 rhc[i][j] = icharge*jcharge * ( - erfij + clw * rh1 ) / rrs2; 307 rhc[j][i] = rhc[i][j]; << 364 >> 365 rhc[j][i] = rhc[i][j]; >> 366 308 } 367 } >> 368 309 } 369 } 310 370 >> 371 >> 372 311 void G4QMDMeanField::CalGraduate() 373 void G4QMDMeanField::CalGraduate() 312 { 374 { >> 375 313 ffr.resize( system->GetTotalNumberOfPartici 376 ffr.resize( system->GetTotalNumberOfParticipant() ); 314 ffp.resize( system->GetTotalNumberOfPartici 377 ffp.resize( system->GetTotalNumberOfParticipant() ); 315 rh3d.resize( system->GetTotalNumberOfPartic 378 rh3d.resize( system->GetTotalNumberOfParticipant() ); 316 379 317 for ( G4int i = 0 ; i < system->GetTotalNum << 380 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ ) 318 { 381 { 319 G4double rho3 = 0.0; 382 G4double rho3 = 0.0; 320 for ( G4int j = 0 ; j < system->GetTotal << 383 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) 321 { 384 { 322 rho3 += rha[j][i]; 385 rho3 += rha[j][i]; 323 } 386 } 324 rh3d[i] = G4Pow::GetInstance()->powA ( r << 387 rh3d[i] = std::pow ( rho3 , pag ); 325 } 388 } 326 389 327 for ( G4int i = 0 ; i < system->GetTotalNum << 390 >> 391 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i ++ ) 328 { 392 { >> 393 329 G4ThreeVector ri = system->GetParticipan 394 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 330 G4LorentzVector p4i = system->GetPartici 395 G4LorentzVector p4i = system->GetParticipant( i )->Get4Momentum(); 331 396 332 G4ThreeVector betai = p4i.v()/p4i.e(); 397 G4ThreeVector betai = p4i.v()/p4i.e(); 333 398 334 // R-JQMD << 399 // R-JQMD 335 G4double Vi = GetPotential( i ); 400 G4double Vi = GetPotential( i ); 336 G4double p_zero = std::sqrt( p4i.e()*p4i 401 G4double p_zero = std::sqrt( p4i.e()*p4i.e() + 2*p4i.m()*Vi); 337 G4ThreeVector betai_R = p4i.v()/p_zero; 402 G4ThreeVector betai_R = p4i.v()/p_zero; 338 G4double mi_R = p4i.m()/p_zero; 403 G4double mi_R = p4i.m()/p_zero; 339 << 404 // 340 ffr[i] = betai_R; 405 ffr[i] = betai_R; 341 ffp[i] = G4ThreeVector( 0.0 ); 406 ffp[i] = G4ThreeVector( 0.0 ); 342 407 343 for ( G4int j = 0 ; j < system->GetTotal << 408 if ( false ) 344 { 409 { >> 410 ffr[i] = betai; >> 411 mi_R = 1.0; >> 412 } >> 413 >> 414 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j ++ ) >> 415 { >> 416 345 G4ThreeVector rj = system->GetPartici 417 G4ThreeVector rj = system->GetParticipant( j )->GetPosition(); 346 G4LorentzVector p4j = system->GetPart 418 G4LorentzVector p4j = system->GetParticipant( j )->Get4Momentum(); 347 419 348 G4double eij = p4i.e() + p4j.e(); 420 G4double eij = p4i.e() + p4j.e(); 349 421 350 G4int icharge = system->GetParticipan 422 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 351 G4int jcharge = system->GetParticipan 423 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 352 424 353 G4int inuc = system->GetParticipant(i 425 G4int inuc = system->GetParticipant(i)->GetNuc(); 354 G4int jnuc = system->GetParticipant(j 426 G4int jnuc = system->GetParticipant(j)->GetNuc(); 355 427 356 G4double ccpp = c0g * rha[j][i] 428 G4double ccpp = c0g * rha[j][i] 357 + c3g * rha[j][i] * ( r 429 + c3g * rha[j][i] * ( rh3d[j] + rh3d[i] ) 358 + csg * rha[j][i] * jnu 430 + csg * rha[j][i] * jnuc * inuc 359 * ( 1. - 2. * std << 431 * ( 1. - 2. * std::abs( jcharge - icharge ) ) 360 + cl * rhc[j][i]; 432 + cl * rhc[j][i]; 361 ccpp *= mi_R; 433 ccpp *= mi_R; 362 434 >> 435 /* >> 436 std::cout << c0g << " " << c3g << " " << csg << " " << cl << std::endl; >> 437 std::cout << "ccpp " << i << " " << j << " " << ccpp << std::endl; >> 438 std::cout << "rha[j][i] " << rha[j][i] << std::endl; >> 439 std::cout << "rh3d " << rh3d[j] << " " << rh3d[i] << std::endl; >> 440 std::cout << "rhc[j][i] " << rhc[j][i] << std::endl; >> 441 */ >> 442 363 G4double grbb = - rbij[j][i]; 443 G4double grbb = - rbij[j][i]; 364 G4double ccrr = grbb * ccpp / eij; 444 G4double ccrr = grbb * ccpp / eij; 365 445 >> 446 /* >> 447 std::cout << "ccrr " << ccrr << std::endl; >> 448 std::cout << "grbb " << grbb << std::endl; >> 449 */ >> 450 >> 451 366 G4ThreeVector rij = ri - rj; 452 G4ThreeVector rij = ri - rj; 367 G4ThreeVector betaij = ( p4i + p4j ) << 453 G4ThreeVector betaij = ( p4i + p4j ).v()/eij; >> 454 368 G4ThreeVector cij = betaij - betai; 455 G4ThreeVector cij = betaij - betai; 369 456 370 ffr[i] = ffr[i] + 2*ccrr* ( rij + grb 457 ffr[i] = ffr[i] + 2*ccrr* ( rij + grbb*cij ); >> 458 371 ffp[i] = ffp[i] - 2*ccpp* ( rij + grb 459 ffp[i] = ffp[i] - 2*ccpp* ( rij + grbb*betaij ); >> 460 372 } 461 } 373 } 462 } >> 463 >> 464 //std::cout << "gradu 0 " << ffr[0] << " " << ffp[0] << std::endl; >> 465 //std::cout << "gradu 1 " << ffr[1] << " " << ffp[1] << std::endl; >> 466 374 } 467 } 375 468 >> 469 >> 470 376 G4double G4QMDMeanField::GetPotential( G4int i 471 G4double G4QMDMeanField::GetPotential( G4int i ) 377 { 472 { 378 G4int n = system->GetTotalNumberOfParticipa 473 G4int n = system->GetTotalNumberOfParticipant(); 379 474 380 G4double rhoa = 0.0; 475 G4double rhoa = 0.0; 381 G4double rho3 = 0.0; 476 G4double rho3 = 0.0; 382 G4double rhos = 0.0; 477 G4double rhos = 0.0; 383 G4double rhoc = 0.0; 478 G4double rhoc = 0.0; 384 479 >> 480 385 G4int icharge = system->GetParticipant(i)-> 481 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 386 G4int inuc = system->GetParticipant(i)->Get 482 G4int inuc = system->GetParticipant(i)->GetNuc(); 387 483 388 for ( G4int j = 0 ; j < n ; ++j ) << 484 for ( G4int j = 0 ; j < n ; j ++ ) 389 { 485 { 390 G4int jcharge = system->GetParticipant(j 486 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 391 G4int jnuc = system->GetParticipant(j)-> 487 G4int jnuc = system->GetParticipant(j)->GetNuc(); 392 488 393 rhoa += rha[j][i]; 489 rhoa += rha[j][i]; 394 rhoc += rhe[j][i]; 490 rhoc += rhe[j][i]; 395 rhos += rha[j][i] * jnuc * inuc 491 rhos += rha[j][i] * jnuc * inuc 396 * ( 1 - 2 * std::abs ( << 492 * ( 1 - 2 * std::abs ( jcharge - icharge ) ); 397 } 493 } 398 494 399 rho3 = G4Pow::GetInstance()->powA ( rhoa , << 495 rho3 = std::pow ( rhoa , gamm ); >> 496 >> 497 G4double potential = c0 * rhoa >> 498 + c3 * rho3 >> 499 + cs * rhos >> 500 + cl * rhoc; 400 501 401 // return potential << 502 return potential; 402 return c0 * rhoa + c3 * rho3 + cs * rhos + << 403 } 503 } 404 504 >> 505 >> 506 405 G4double G4QMDMeanField::GetTotalPotential() 507 G4double G4QMDMeanField::GetTotalPotential() 406 { 508 { >> 509 407 G4int n = system->GetTotalNumberOfParticipa 510 G4int n = system->GetTotalNumberOfParticipant(); 408 511 409 std::vector < G4double > rhoa ( n , 0.0 ); 512 std::vector < G4double > rhoa ( n , 0.0 ); 410 std::vector < G4double > rho3 ( n , 0.0 ); 513 std::vector < G4double > rho3 ( n , 0.0 ); 411 std::vector < G4double > rhos ( n , 0.0 ); 514 std::vector < G4double > rhos ( n , 0.0 ); 412 std::vector < G4double > rhoc ( n , 0.0 ); 515 std::vector < G4double > rhoc ( n , 0.0 ); >> 516 413 517 414 for ( G4int i = 0 ; i < n ; ++i ) << 518 for ( G4int i = 0 ; i < n ; i ++ ) 415 { 519 { 416 G4int icharge = system->GetParticipant(i 520 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 417 G4int inuc = system->GetParticipant(i)-> 521 G4int inuc = system->GetParticipant(i)->GetNuc(); 418 522 419 for ( G4int j = 0 ; j < n ; ++j ) << 523 for ( G4int j = 0 ; j < n ; j ++ ) 420 { 524 { 421 G4int jcharge = system->GetParticipan 525 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 422 G4int jnuc = system->GetParticipant(j 526 G4int jnuc = system->GetParticipant(j)->GetNuc(); 423 527 424 rhoa[i] += rha[j][i]; 528 rhoa[i] += rha[j][i]; 425 rhoc[i] += rhe[j][i]; 529 rhoc[i] += rhe[j][i]; 426 rhos[i] += rha[j][i] * jnuc * inuc 530 rhos[i] += rha[j][i] * jnuc * inuc 427 * ( 1 - 2 * std::abs ( jcharg << 531 * ( 1 - 2 * std::abs ( jcharge - icharge ) ); 428 } 532 } 429 533 430 rho3[i] = G4Pow::GetInstance()->powA ( r << 534 rho3[i] = std::pow ( rhoa[i] , gamm ); 431 } 535 } 432 536 433 // return potential << 537 G4double potential = c0 * std::accumulate( rhoa.begin() , rhoa.end() , 0.0 ) 434 return c0 * std::accumulate( rhoa.cbegin() << 538 + c3 * std::accumulate( rho3.begin() , rho3.end() , 0.0 ) 435 + c3 * std::accumulate( rho3.cbeg << 539 + cs * std::accumulate( rhos.begin() , rhos.end() , 0.0 ) 436 + cs * std::accumulate( rhos.cbeg << 540 + cl * std::accumulate( rhoc.begin() , rhoc.end() , 0.0 ); 437 + cl * std::accumulate( rhoc.cbeg << 541 >> 542 return potential; >> 543 438 } 544 } 439 545 >> 546 >> 547 440 G4double G4QMDMeanField::calPauliBlockingFacto 548 G4double G4QMDMeanField::calPauliBlockingFactor( G4int i ) 441 { 549 { 442 // i is supposed beyond total number of Par << 443 550 444 G4double pf = 0.0; 551 G4double pf = 0.0; >> 552 // i is supposed beyond total number of Participant() 445 G4int icharge = system->GetParticipant(i)-> 553 G4int icharge = system->GetParticipant(i)->GetChargeInUnitOfEplus(); 446 554 447 for ( G4int j = 0 ; j < system->GetTotalNum << 555 for ( G4int j = 0 ; j < system->GetTotalNumberOfParticipant() ; j++ ) 448 { 556 { >> 557 449 G4int jcharge = system->GetParticipant(j 558 G4int jcharge = system->GetParticipant(j)->GetChargeInUnitOfEplus(); 450 G4int jnuc = system->GetParticipant(j)-> 559 G4int jnuc = system->GetParticipant(j)->GetNuc(); 451 560 452 if ( jcharge == icharge && jnuc == 1 ) 561 if ( jcharge == icharge && jnuc == 1 ) 453 { 562 { 454 G4double expa = -rr2[i][j]*cpw; << 563 >> 564 /* >> 565 std::cout << "Pauli i j " << i << " " << j << std::endl; >> 566 std::cout << "Pauli icharge " << icharge << std::endl; >> 567 std::cout << "Pauli jcharge " << jcharge << std::endl; >> 568 */ >> 569 G4double expa = -rr2[i][j]*cpw; >> 570 >> 571 455 if ( expa > epsx ) 572 if ( expa > epsx ) 456 { 573 { 457 expa = expa - pp2[i][j]*cph; 574 expa = expa - pp2[i][j]*cph; >> 575 /* >> 576 std::cout << "Pauli cph " << cph << std::endl; >> 577 std::cout << "Pauli pp2 " << pp2[i][j] << std::endl; >> 578 std::cout << "Pauli expa " << expa << std::endl; >> 579 std::cout << "Pauli epsx " << epsx << std::endl; >> 580 */ 458 if ( expa > epsx ) 581 if ( expa > epsx ) 459 { 582 { 460 pf = pf + G4Exp ( expa ); << 583 // std::cout << "Pauli phase " << pf << std::endl; >> 584 pf = pf + std::exp ( expa ); 461 } 585 } 462 } 586 } 463 } 587 } >> 588 464 } 589 } 465 590 466 return ( pf - 1.0 ) * cpc; << 591 >> 592 pf = ( pf - 1.0 ) * cpc; >> 593 >> 594 //std::cout << "Pauli pf " << pf << std::endl; >> 595 >> 596 return pf; >> 597 467 } 598 } 468 599 >> 600 >> 601 469 G4bool G4QMDMeanField::IsPauliBlocked( G4int i 602 G4bool G4QMDMeanField::IsPauliBlocked( G4int i ) 470 { 603 { 471 G4bool result = false; 604 G4bool result = false; 472 605 473 if ( system->GetParticipant( i )->GetNuc() 606 if ( system->GetParticipant( i )->GetNuc() == 1 ) 474 { 607 { 475 G4double pf = calPauliBlockingFactor( i 608 G4double pf = calPauliBlockingFactor( i ); 476 G4double rand = G4UniformRand(); 609 G4double rand = G4UniformRand(); 477 if ( pf > rand ) { result = true; } << 610 if ( pf > rand ) result = true; 478 } 611 } 479 612 480 return result; 613 return result; 481 } 614 } 482 615 >> 616 >> 617 483 void G4QMDMeanField::DoPropagation( G4double d 618 void G4QMDMeanField::DoPropagation( G4double dt ) 484 { 619 { >> 620 485 G4double cc2 = 1.0; 621 G4double cc2 = 1.0; 486 G4double cc1 = 1.0 - cc2; 622 G4double cc1 = 1.0 - cc2; 487 G4double cc3 = 1.0 / 2.0 / cc2; 623 G4double cc3 = 1.0 / 2.0 / cc2; 488 624 489 G4double dt3 = dt * cc3; 625 G4double dt3 = dt * cc3; 490 G4double dt1 = dt * ( cc1 - cc3 ); 626 G4double dt1 = dt * ( cc1 - cc3 ); 491 G4double dt2 = dt * cc2; 627 G4double dt2 = dt * cc2; 492 628 493 CalGraduate(); 629 CalGraduate(); 494 630 495 G4int n = system->GetTotalNumberOfParticipa 631 G4int n = system->GetTotalNumberOfParticipant(); 496 632 497 // 1st Step << 633 // 1st Step 498 634 499 std::vector< G4ThreeVector > f0r, f0p; 635 std::vector< G4ThreeVector > f0r, f0p; 500 f0r.resize( n ); 636 f0r.resize( n ); 501 f0p.resize( n ); 637 f0p.resize( n ); 502 638 503 for ( G4int i = 0 ; i < n ; ++i ) << 639 for ( G4int i = 0 ; i < n ; i++ ) 504 { 640 { 505 G4ThreeVector ri = system->GetParticipan 641 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 506 G4ThreeVector p3i = system->GetParticipa 642 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 507 643 508 ri += dt3* ffr[i]; 644 ri += dt3* ffr[i]; 509 p3i += dt3* ffp[i]; 645 p3i += dt3* ffp[i]; 510 646 511 f0r[i] = ffr[i]; 647 f0r[i] = ffr[i]; 512 f0p[i] = ffp[i]; 648 f0p[i] = ffp[i]; 513 649 514 system->GetParticipant( i )->SetPosition 650 system->GetParticipant( i )->SetPosition( ri ); 515 system->GetParticipant( i )->SetMomentum 651 system->GetParticipant( i )->SetMomentum( p3i ); 516 652 517 // we do not need set total momentum by << 653 // we do not need set total momentum by ourselvs 518 } 654 } 519 655 520 // 2nd Step << 656 // 2nd Step 521 << 522 Cal2BodyQuantities(); 657 Cal2BodyQuantities(); 523 CalGraduate(); 658 CalGraduate(); 524 659 525 for ( G4int i = 0 ; i < n ; ++i ) << 660 for ( G4int i = 0 ; i < n ; i++ ) 526 { 661 { 527 G4ThreeVector ri = system->GetParticipan 662 G4ThreeVector ri = system->GetParticipant( i )->GetPosition(); 528 G4ThreeVector p3i = system->GetParticipa 663 G4ThreeVector p3i = system->GetParticipant( i )->GetMomentum(); 529 664 530 ri += dt1* f0r[i] + dt2* ffr[i]; 665 ri += dt1* f0r[i] + dt2* ffr[i]; 531 p3i += dt1* f0p[i] + dt2* ffp[i]; 666 p3i += dt1* f0p[i] + dt2* ffp[i]; 532 667 533 system->GetParticipant( i )->SetPosition 668 system->GetParticipant( i )->SetPosition( ri ); 534 system->GetParticipant( i )->SetMomentum 669 system->GetParticipant( i )->SetMomentum( p3i ); 535 670 536 // we do not need set total momentum by << 671 // we do not need set total momentum by ourselvs 537 } 672 } 538 673 539 Cal2BodyQuantities(); << 674 Cal2BodyQuantities(); >> 675 540 } 676 } 541 677 >> 678 >> 679 542 std::vector< G4QMDNucleus* > G4QMDMeanField::D 680 std::vector< G4QMDNucleus* > G4QMDMeanField::DoClusterJudgment() 543 { 681 { >> 682 >> 683 //std::cout << "MeanField DoClusterJudgemnt" << std::endl; >> 684 544 Cal2BodyQuantities(); 685 Cal2BodyQuantities(); 545 686 546 G4double cpf2 = G4Pow::GetInstance()->A23 ( << 687 G4double cpf2 = std::pow ( 1.5 * pi*pi * std::pow ( 4.0 * pi * wl , -1.5 ) >> 688 , >> 689 2./3. ) 547 * hbc * hbc; 690 * hbc * hbc; 548 G4double rcc2 = rclds*rclds; 691 G4double rcc2 = rclds*rclds; 549 692 550 G4int n = system->GetTotalNumberOfParticipa 693 G4int n = system->GetTotalNumberOfParticipant(); 551 std::vector < G4double > rhoa; 694 std::vector < G4double > rhoa; 552 rhoa.resize ( n ); 695 rhoa.resize ( n ); 553 696 554 for ( G4int i = 0 ; i < n ; ++i ) << 697 for ( G4int i = 0 ; i < n ; i++ ) 555 { 698 { 556 rhoa[i] = 0.0; << 699 rhoa[i] = 0.0; 557 700 558 if ( system->GetParticipant( i )->GetBary << 701 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) 559 { << 702 { 560 for ( G4int j = 0 ; j < n ; ++j ) << 703 for ( G4int j = 0 ; j < n ; j++ ) 561 { << 704 { 562 if ( system->GetParticipant( j )->Get 705 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 563 rhoa[i] += rha[i][j]; 706 rhoa[i] += rha[i][j]; 564 } << 707 } 565 } << 708 } >> 709 >> 710 rhoa[i] = std::pow ( rhoa[i] + 1 , 1.0/3.0 ); 566 711 567 rhoa[i] = G4Pow::GetInstance()->A13 ( rho << 568 } 712 } 569 713 570 // identification of the cluster << 714 // identification of the cluster >> 715 >> 716 std::map < G4int , std::vector < G4int > > cluster_map; 571 std::vector < G4bool > is_already_belong_so 717 std::vector < G4bool > is_already_belong_some_cluster; 572 718 573 // cluster_id participant_id 719 // cluster_id participant_id 574 std::multimap < G4int , G4int > comb_map; 720 std::multimap < G4int , G4int > comb_map; 575 std::multimap < G4int , G4int > assign_map; 721 std::multimap < G4int , G4int > assign_map; 576 assign_map.clear(); 722 assign_map.clear(); 577 723 578 std::vector < G4int > mascl; 724 std::vector < G4int > mascl; 579 std::vector < G4int > num; 725 std::vector < G4int > num; 580 mascl.resize ( n ); 726 mascl.resize ( n ); 581 num.resize ( n ); 727 num.resize ( n ); 582 is_already_belong_some_cluster.resize ( n ) 728 is_already_belong_some_cluster.resize ( n ); 583 729 584 std::vector < G4int > is_assigned_to ( n , 730 std::vector < G4int > is_assigned_to ( n , -1 ); 585 std::multimap < G4int , G4int > clusters; 731 std::multimap < G4int , G4int > clusters; 586 732 587 for ( G4int i = 0 ; i < n ; ++i ) << 733 for ( G4int i = 0 ; i < n ; i++ ) 588 { 734 { 589 mascl[i] = 1; << 735 mascl[i] = 1; 590 num[i] = 1; << 736 num[i] = 1; 591 is_already_belong_some_cluster[i] = false << 737 >> 738 is_already_belong_some_cluster[i] = false; 592 } 739 } 593 740 >> 741 >> 742 G4int nclst = 1; 594 G4int ichek = 1; 743 G4int ichek = 1; >> 744 >> 745 595 G4int id = 0; 746 G4int id = 0; 596 G4int cluster_id = -1; 747 G4int cluster_id = -1; 597 for ( G4int i = 0 ; i < n-1 ; ++i ) << 748 for ( G4int i = 0 ; i < n-1 ; i++ ) 598 { 749 { 599 G4bool hasThisCompany = false; << 600 750 601 if ( system->GetParticipant( i )->GetBary << 751 G4bool hasThisCompany = false; 602 { << 752 // Check only for bryons? 603 G4int j1 = i + 1; << 753 // std::cout << "Check Baryon " << i << std::endl; 604 for ( G4int j = j1 ; j < n ; ++j ) << 754 605 { << 755 if ( system->GetParticipant( i )->GetBaryonNumber() == 1 ) >> 756 { >> 757 >> 758 // if ( is_already_belong_some_cluster[i] != true ) >> 759 // { >> 760 //G4int j1 = ichek + 1; >> 761 G4int j1 = i + 1; >> 762 for ( G4int j = j1 ; j < n ; j++ ) >> 763 { >> 764 606 std::vector < G4int > cluster_partici 765 std::vector < G4int > cluster_participants; 607 if ( system->GetParticipant( j )->Get 766 if ( system->GetParticipant( j )->GetBaryonNumber() == 1 ) 608 { 767 { 609 G4double rdist2 = rr2[ i ][ j ]; << 768 G4double rdist2 = rr2[ i ][ j ]; 610 G4double pdist2 = pp2[ i ][ j ]; << 769 G4double pdist2 = pp2[ i ][ j ]; 611 G4double pcc2 = cpf2 << 770 //G4double rdist2 = rr2[ num[i] ][ num[j] ]; 612 * ( rhoa[ i ] + rhoa[ << 771 //G4double pdist2 = pp2[ num[i] ][ num[j] ]; 613 * ( rhoa[ i ] + rhoa[ << 772 G4double pcc2 = cpf2 614 << 773 * ( rhoa[ i ] + rhoa[ j ] ) 615 // Check phase space: close enough? << 774 * ( rhoa[ i ] + rhoa[ j ] ); 616 if ( rdist2 < rcc2 && pdist2 < pcc2 << 775 617 { << 776 // Check phase space: close enough? 618 if ( is_assigned_to [ j ] == -1 ) << 777 if ( rdist2 < rcc2 && pdist2 < pcc2 ) 619 { << 778 { >> 779 >> 780 /* >> 781 std::cout << "G4QMDRESULT " >> 782 << i << " " << j << " " << id << " " >> 783 << is_assigned_to [ i ] << " " << is_assigned_to [ j ] >> 784 << std::endl; >> 785 */ >> 786 >> 787 if ( is_assigned_to [ j ] == -1 ) >> 788 { 620 if ( is_assigned_to [ i ] == -1 789 if ( is_assigned_to [ i ] == -1 ) 621 { 790 { 622 if ( clusters.size() != 0 ) << 791 if ( clusters.size() != 0 ) 623 { << 792 { 624 id = clusters.rbegin()->fir << 793 id = clusters.rbegin()->first + 1; 625 } << 794 //std::cout << "id is increare " << id << std::endl; 626 else << 795 } 627 { << 796 else 628 id = 0; << 797 { 629 } << 798 id = 0; 630 clusters.insert ( std::multim << 799 } 631 is_assigned_to [ i ] = id; << 800 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , i ) ); 632 clusters.insert ( std::multim << 801 is_assigned_to [ i ] = id; 633 is_assigned_to [ j ] = id; << 802 clusters.insert ( std::multimap<G4int,G4int>::value_type ( id , j ) ); >> 803 is_assigned_to [ j ] = id; 634 } 804 } 635 else 805 else 636 { 806 { 637 clusters.insert ( std::multim << 807 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ i ] , j ) ); 638 is_assigned_to [ j ] = is_ass << 808 is_assigned_to [ j ] = is_assigned_to [ i ]; 639 } 809 } 640 } << 810 } 641 else << 811 else 642 { << 812 { 643 // j is already belong to some << 813 // j is already belong to some cluester 644 if ( is_assigned_to [ i ] == -1 814 if ( is_assigned_to [ i ] == -1 ) 645 { 815 { 646 clusters.insert ( std::multim << 816 clusters.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , i ) ); 647 is_assigned_to [ i ] = is_ass << 817 is_assigned_to [ i ] = is_assigned_to [ j ]; 648 } 818 } 649 else 819 else 650 { 820 { 651 // i has companion << 821 // i has companion 652 if ( is_assigned_to [ i ] != << 822 if ( is_assigned_to [ i ] != is_assigned_to [ j ] ) 653 { << 823 { 654 // move companions to the c << 824 // move companions to the cluster 655 std::multimap< G4int , G4in << 825 // 656 G4int target_cluster_id; << 826 //std::cout << "combine " << is_assigned_to [ i ] << " to " << is_assigned_to [ j ] << std::endl; 657 if ( is_assigned_to [ i ] > << 827 std::multimap< G4int , G4int > clusters_tmp; 658 { << 828 G4int target_cluster_id; 659 target_cluster_id = is_as << 829 if ( is_assigned_to [ i ] > is_assigned_to [ j ] ) 660 } << 830 target_cluster_id = is_assigned_to [ i ]; 661 else << 662 { << 663 target_cluster_id = is_as << 664 } << 665 for ( auto it = clusters.cb << 666 { << 667 if ( it->first == target_ << 668 { << 669 is_assigned_to [ it->se << 670 clusters_tmp.insert(std << 671 } << 672 else 831 else >> 832 target_cluster_id = is_assigned_to [ j ]; >> 833 >> 834 for ( std::multimap< G4int , G4int >::iterator it >> 835 = clusters.begin() ; it != clusters.end() ; it++ ) 673 { 836 { 674 clusters_tmp.insert(std << 837 >> 838 //std::cout << it->first << " " << it->second << " " << target_cluster_id << std::endl; >> 839 if ( it->first == target_cluster_id ) >> 840 { >> 841 //std::cout << "move " << it->first << " " << it->second << std::endl; >> 842 is_assigned_to [ it->second ] = is_assigned_to [ j ]; >> 843 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( is_assigned_to [ j ] , it->second ) ); >> 844 } >> 845 else >> 846 { >> 847 clusters_tmp.insert ( std::multimap<G4int,G4int>::value_type ( it->first , it->second ) ); >> 848 } 675 } 849 } 676 } << 850 677 clusters = std::move(cluste << 851 clusters = clusters_tmp; 678 } << 852 //id = clusters.rbegin()->first; >> 853 //id = target_cluster_id; >> 854 //std::cout << "id " << id << std::endl; >> 855 } 679 } 856 } 680 } << 857 } >> 858 >> 859 //std::cout << "combination " << i << " " << j << std::endl; >> 860 comb_map.insert( std::multimap<G4int,G4int>::value_type ( i , j ) ); >> 861 cluster_participants.push_back ( j ); >> 862 681 863 682 comb_map.insert( std::multimap<G4 << 683 cluster_participants.push_back ( << 684 864 685 if ( assign_map.find( cluster_id << 865 if ( assign_map.find( cluster_id ) == assign_map.end() ) 686 { << 866 { 687 is_already_belong_some_cluster[ 867 is_already_belong_some_cluster[i] = true; 688 assign_map.insert ( std::multim 868 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , i ) ); 689 hasThisCompany = true; 869 hasThisCompany = true; 690 } << 870 } 691 assign_map.insert ( std::multimap << 871 assign_map.insert ( std::multimap<G4int,G4int>::value_type ( cluster_id , j ) ); 692 is_already_belong_some_cluster[j] << 872 is_already_belong_some_cluster[j] = true; 693 } << 873 694 << 874 } 695 if ( ichek == i ) << 875 696 { << 876 if ( ichek == i ) 697 ++ichek; << 877 { 698 } << 878 nclst++; >> 879 ichek++; >> 880 } >> 881 } >> 882 >> 883 if ( cluster_participants.size() > 0 ) >> 884 { >> 885 // cluster , participant >> 886 cluster_map.insert ( std::pair < G4int , std::vector < G4int > > ( i , cluster_participants ) ); 699 } 887 } 700 } << 888 } 701 } << 889 // } 702 if ( hasThisCompany == true ) { ++cluster << 890 } >> 891 if ( hasThisCompany == true ) cluster_id++; 703 } 892 } 704 893 705 // sort << 894 //std::cout << " id " << id << std::endl; 706 // Heavy cluster comes first << 895 707 // size cluster_id << 896 // sort >> 897 // Heavy cluster comes first >> 898 // size cluster_id 708 std::multimap< G4int , G4int > sorted_clust 899 std::multimap< G4int , G4int > sorted_cluster_map; 709 for ( G4int i = 0 ; i <= id ; ++i ) // << << 900 for ( G4int i = 0 ; i <= id ; i++ ) // << "<=" because id is highest cluster nubmer. 710 { 901 { 711 sorted_cluster_map.insert ( std::multimap << 902 >> 903 //std::cout << i << " cluster has " << clusters.count( i ) << " nucleons." << std::endl; >> 904 sorted_cluster_map.insert ( std::multimap<G4int,G4int>::value_type ( (G4int) clusters.count( i ) , i ) ); >> 905 712 } 906 } 713 << 907 714 // create nucleus from divided clusters << 908 >> 909 // create nucleus from devided clusters 715 std::vector < G4QMDNucleus* > result; 910 std::vector < G4QMDNucleus* > result; 716 for ( auto it = sorted_cluster_map.crbegin( << 911 for ( std::multimap < G4int , G4int >::reverse_iterator it >> 912 = sorted_cluster_map.rbegin() ; it != sorted_cluster_map.rend() ; it ++) 717 { 913 { 718 if ( it->first != 0 ) << 914 719 { << 915 //G4cout << "Add Participants to cluseter " << it->second << G4endl; 720 G4QMDNucleus* nucleus = new G4QMDNucleu << 916 721 for ( auto itt = clusters.cbegin(); itt << 917 if ( it->first != 0 ) 722 { << 918 { 723 if ( it->second == itt->first ) << 919 G4QMDNucleus* nucleus = new G4QMDNucleus(); >> 920 for ( std::multimap < G4int , G4int >::iterator itt >> 921 = clusters.begin() ; itt != clusters.end() ; itt ++) 724 { 922 { 725 nucleus->SetParticipant( system->Ge << 923 >> 924 if ( it->second == itt->first ) >> 925 { >> 926 nucleus->SetParticipant( system->GetParticipant ( itt->second ) ); >> 927 //G4cout << "Add Participants " << itt->second << " " << system->GetParticipant ( itt->second )->GetPosition() << G4endl; >> 928 } >> 929 726 } 930 } 727 } << 931 result.push_back( nucleus ); 728 result.push_back( nucleus ); << 932 } 729 } << 933 730 } 934 } 731 935 732 // delete participants from current system << 936 // delete participants from current system 733 for ( auto it = result.cbegin(); it != resu << 937 >> 938 for ( std::vector < G4QMDNucleus* > ::iterator it >> 939 = result.begin() ; it != result.end() ; it++ ) 734 { 940 { 735 system->SubtractSystem ( *it ); << 941 system->SubtractSystem ( *it ); 736 } 942 } 737 943 738 return result; 944 return result; >> 945 739 } 946 } >> 947 >> 948 740 949 741 void G4QMDMeanField::Update() 950 void G4QMDMeanField::Update() 742 { 951 { 743 SetSystem( system ); 952 SetSystem( system ); 744 } 953 } 745 954