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 // Author: Alexei Sytov 27 // Co-author: Gianfranco Paterno (testing) 28 // Using the key points of G4BaierKatkov and d 29 // partially described in L. Bandiera et al. E 30 31 #include "G4CoherentPairProduction.hh" 32 33 #include "Randomize.hh" 34 #include "G4TouchableHistory.hh" 35 #include "G4TouchableHandle.hh" 36 #include "G4SystemOfUnits.hh" 37 38 #include "G4Track.hh" 39 #include "G4Gamma.hh" 40 #include "G4Electron.hh" 41 #include "G4Positron.hh" 42 43 #include "G4ParticleDefinition.hh" 44 #include "G4ProcessManager.hh" 45 #include "G4EmProcessSubType.hh" 46 #include "G4TransportationManager.hh" 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 50 G4CoherentPairProduction::G4CoherentPairProduc 51 52 G4VDiscreteProcess(aName) 53 { 54 SetProcessSubType(fCoherentPairProduction) 55 } 56 57 //....oooOO0OOooo........oooOO0OOooo........oo 58 59 G4double G4CoherentPairProduction::GetMeanFree 60 G4double, 61 G4ForceCon 62 { 63 //current logical volume 64 G4LogicalVolume* crystallogic; 65 66 //momentum direction and coordinates (see 67 G4ThreeVector momentumDirectionGamma,xyzGa 68 //angle of the photon in the local referen 69 G4double txGamma0 = 0, tyGamma0 = 0; 70 71 *condition = NotForced; 72 73 //model activation 74 G4bool modelTrigger = false; 75 76 //photon energy 77 G4double eGamma = aTrack.GetTotalEnergy(); 78 79 //energy cut, at the beginning, to not che 80 if(eGamma > ModelMinPrimaryEnergy()) 81 { 82 //current logical volume 83 crystallogic = aTrack.GetVolume()->Get 84 85 //the model works only in the G4Region 86 if(crystallogic->GetRegion()->GetName( 87 { 88 fCrystalData->SetGeometryParameter 89 90 //the momentum direction of the ph 91 momentumDirectionGamma = 92 (aTrack.GetTouchableHandle()-> 93 GetTopTransform().NetRotat 94 95 //the coordinates of the photon in 96 xyzGamma0 = 97 aTrack.GetTouchableHandle()->G 98 GetTopTransform().Transfor 99 100 // the coordinates of the photon i 101 //a channel (elementary periodic c 102 xyzGamma = fCrystalData->Coordinat 103 104 //angle of the photon in the local 105 //(!!! ONLY FORWARD DIRECTION, mom 106 txGamma0 = std::atan(momentumDirec 107 tyGamma0 = std::atan(momentumDirec 108 109 //recalculate angle into the latti 110 G4double angle = fCrystalData->Ang 111 if (fCrystalData->GetModel()==2) 112 { 113 angle = std::sqrt(angle*angle+ 114 } 115 116 //Applies the parameterisation not 117 //above low energy limit and below 118 modelTrigger = (momentumDirectionG 119 std::abs(angle) < 120 } 121 } 122 123 if(modelTrigger) 124 { 125 //execute the model 126 127 G4double x=0.,y=0.,z=0.;// the coordin 128 //in the co-rotating reference sys 129 //a channel (elementary periodic c 130 G4double tx0=0.,ty0=0.; // the angles 131 // in the local reference system o 132 G4double txPreStep0=0.,tyPreStep0=0.; 133 // in the co-rotating reference sy 134 //a channel (elementary periodic c 135 136 G4ThreeVector scatteringAnglesAndEnerg 137 138 //coordinates in Runge-Kutta calculati 139 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0. 140 //angles in Runge-Kutta calculations 141 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,t 142 //variables in Runge-Kutta calculation 143 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4= 144 //simulation step along z (internal st 145 G4double dz=0.,dzd3=0.,dzd8=0.;//dzd3 146 //simulation step along the momentum d 147 G4double momentumDirectionStep; 148 //effective simulation step (taking in 149 G4double effectiveStep=0.; 150 151 // Baier-Katkov variables 152 G4double dzMeV=0.; //step in MeV^-1 153 G4double axt=0.,ayt=0.; //charged part 154 G4double vxin=0.,vyin=0.;//the angles 155 G4double vxno=0.,vyno=0.;//the angles 156 157 G4double dzmod=0.; 158 G4double fa1=0.,faseBefore=0.,faseBefo 159 G4double faseAfter=0.,fa2dfaseBefore2= 160 161 G4double skJ=0, skIx=0., skIy=0.; 162 G4double sinfa1=0.,cosfa1=0.; 163 164 //2-vector is needed for an initial pa 165 //vector of 2-vectors is an initial pa 166 167 //collection of etotal for a single pa 168 CLHEP::Hep2Vector twoVectorEtotal(0.,0 169 170 //collection of x for a single pair 171 CLHEP::Hep2Vector twoVectorX(0.,0.); 172 173 //collection of y for a single pair 174 CLHEP::Hep2Vector twoVectorY(0.,0.); 175 176 //collection of tx for a single pair 177 CLHEP::Hep2Vector twoVectorTX(0.,0.); 178 179 //collection of tx for a single pair 180 CLHEP::Hep2Vector twoVectorTY(0.,0.); 181 182 fullVectorEtotal.clear(); 183 fullVectorX.clear(); 184 fullVectorY.clear(); 185 fullVectorTX.clear(); 186 fullVectorTY.clear(); 187 fPairProductionCDFdz.clear(); 188 fPairProductionCDFdz.push_back(0.);//0 189 190 const G4double charge[2] = {-1.,1.}; / 191 const G4String particleName[2] = {"e-" 192 193 // the coordinates of a charged partic 194 //a channel (elementary periodic cell) 195 G4ThreeVector xyzparticle = xyzGamma;/ 196 197 //the idea of pair production simulati 198 //since the matrix element of these pr 199 //to radiation: sample the pairs, calc 200 //probabilities using Baier-Katkov ana 201 202 //cycle by sampling e+- pairs 203 for(G4int i=0; i<fNMCPairs;i++) 204 { 205 //pair energy uniform sampling 206 G4double etotal = fMass + fPPKinet 207 G4UniformRand()* 208 209 210 G4double phi = CLHEP::twopi*G4Unif 211 212 //the probability of the productio 213 //per distance 214 G4double probabilityPPdz = 0.; 215 216 //cycle e- and e+ within single pa 217 for(G4int j=0; j<2;j++) 218 { 219 if(j==1){etotal=eGamma-etotal; 220 twoVectorEtotal[j]=etotal; 221 222 //Baier-Katkov input 223 //intermediate variables to re 224 G4double e2 = etotal*etotal; 225 G4double gammaInverse2 = fMass 226 //normalization coefficient 227 G4double coefNorm = CLHEP::fin 228 //G4double phi = CLHEP::twopi* 229 G4double om = eGamma; 230 G4double eprime=om-etotal; //E 231 G4double eprime2 = eprime*epri 232 G4double e2pluseprime2 =e2+epr 233 G4double omprime=etotal*om/epr 234 G4double omprimed2=omprime/2; 235 236 //difference vs G4BaierKatkov: 237 G4double coefNorme2deprime2 = 238 239 G4double gammaInverse2om = gam 240 241 //initialize intermediate inte 242 G4double fa=0.,ss=0.,sc=0.,ssx 243 244 //End of Baier-Katkov input 245 246 G4bool fbreak = false;//flag o 247 248 //set fCrystalData parameters 249 fCrystalData->SetParticlePrope 250 251 252 //needed just to setup the cor 253 //since later it may be change 254 fCrystalData->CoordinatesFromB 255 256 //coordinate sampling: random 257 //in the interaction point 258 if(j==0) 259 { 260 x = fCrystalData->GetChann 261 y = fCrystalData->GetChann 262 } 263 else 264 { 265 x=twoVectorX[0]; 266 y=twoVectorY[0]; 267 } 268 twoVectorX[j] = x; 269 twoVectorY[j] = y; 270 //definite z as a coordinate o 271 //interaction point is taking 272 //of the position of pair prod 273 z = xyzGamma.z(); 274 275 //angles of the photon in the 276 //angular distribution center 277 G4double tx = fCrystalData->An 278 G4double ty = tyGamma0; 279 G4double momentumDirectionZGam 280 281 282 283 //angle sampling: depends on a 284 //defined by the Lindhard angl 285 //to 1/gamma 286 287 //range of MC integration on a 288 G4double paramParticleAngle = 289 290 G4double axangle=0.; 291 if (fCrystalData->GetModel()== 292 { 293 axangle = std::abs(tx); 294 } 295 else if (fCrystalData->GetMod 296 { 297 axangle = std::sqrt(tx*tx+ 298 } 299 300 if(axangle>fCrystalData->GetLi 301 { 302 paramParticleAngle+=axangl 303 -std 304 305 306 } 307 else 308 { 309 paramParticleAngle+=fCryst 310 } 311 312 313 //ONLY forward direction 314 if (paramParticleAngle>CLHEP:: 315 316 G4double rho=1.; 317 G4double rhocut=CLHEP::halfpi/ 318 319 G4double norm=std::atan(rhocut 320 CLHEP::pi*para 321 322 323 //distribution with long tails 324 //after a strong single scatte 325 //at ellipsescale < 1 => half 326 do 327 { 328 rho = std::sqrt(std::tan(C 329 } 330 while (rho>rhocut); 331 332 //normalization coefficient fo 333 G4double angleNormCoef = (1.+r 334 335 tx+=charge[j]*paramParticleAng 336 twoVectorTX[j] = tx; 337 ty+=charge[j]*paramParticleAng 338 twoVectorTY[j] = ty; 339 340 G4double zalongGamma = 0;//nec 341 //depending on the traject 342 //starting the trajectory 343 //here we don't care about the 344 //the trajectory is very short 345 //in Baier-Katkov will be extr 346 for(G4int k=0; k<fNTrajectoryS 347 { 348 //back to the local refere 349 txPreStep0 = fCrystalData- 350 tyPreStep0 = ty; 351 352 dz = fCrystalData->GetSimu 353 dzd3=dz/3; 354 dzd8=dz/8; 355 356 //trajectory calculation: 357 //Runge-Cutt "3/8" 358 //fCrystalData->GetCurv(z) 359 //of the radius on x; GetC 360 361 //first step 362 kvx1=fCrystalData->Ex(x,y) 363 x1=x+tx*dzd3; 364 tx1=tx+(kvx1-fCrystalData- 365 if (fCrystalData->GetModel 366 { 367 kvy1=fCrystalData->Ey( 368 y1=y+ty*dzd3; 369 ty1=ty+kvy1*dzd3; 370 } 371 372 //second step 373 kvx2=fCrystalData->Ex(x1,y 374 x2=x-tx*dzd3+tx1*dz; 375 tx2=tx-(kvx1-fCrystalData- 376 (kvx2-fCrystalData-> 377 if (fCrystalData->GetModel 378 { 379 kvy2=fCrystalData->Ey( 380 y2=y-ty*dzd3+ty1*dz; 381 ty2=ty-kvy1*dzd3+kvy2* 382 } 383 384 //third step 385 kvx3=fCrystalData->Ex(x2,y 386 x3=x+(tx-tx1+tx2)*dz; 387 tx3=tx+(kvx1-kvx2+kvx3- 388 fCrystalData-> 389 if (fCrystalData->GetModel 390 { 391 kvy3=fCrystalData->Ey( 392 y3=y+(ty-ty1+ty2)*dz; 393 ty3=ty+(kvy1-kvy2+kvy3 394 } 395 396 //fourth step 397 kvx4=fCrystalData->Ex(x3,y 398 x4=x+(tx+3.*tx1+3.*tx2+tx3 399 tx4=tx+(kvx1+3.*kvx2+3.*kv 400 fCrystalData->GetCur 401 if (fCrystalData->GetModel 402 { 403 kvy4=fCrystalData->Ey( 404 y4=y+(ty+3.*ty1+3.*ty2 405 ty4=ty+(kvy1+3.*kvy2+3 406 } 407 else 408 { 409 y4 =y+ty*dz; 410 ty4=ty; 411 } 412 413 x=x4; 414 tx=tx4; 415 y=y4; 416 ty=ty4; 417 418 z+=dz*fCrystalData->GetCor 419 //("central plane/axis 420 421 xyzparticle = fCrystalData 422 x=xyzparticle.x(); 423 y=xyzparticle.y(); 424 z=xyzparticle.z(); 425 426 momentumDirectionStep = 427 dz*std::sqrt(1+std::po 428 zalongGamma += dz/momentum 429 430 //default scattering and e 431 scatteringAnglesAndEnergyL 432 433 if(fIncoherentScattering) 434 { 435 //calculate separately 436 for (G4int ii = 0; ii 437 { 438 //effective step ta 439 effectiveStep = mo 440 fC 441 //Coulomb scatteri 442 //(both multiple a 443 scatteringAnglesAn 444 fCrystalData-> 445 446 447 } 448 //electron scattering 449 scatteringAnglesAndEne 450 fCrystalData->MinI 451 fCrystalData->Elec 452 momentumDirectionS 453 tx += scatteringAngles 454 ty += scatteringAngles 455 } 456 457 //To avoid backward direct 458 if(std::abs(tx)>CLHEP::hal 459 std::abs(ty)>CLHEP::ha 460 { 461 G4cout << "Warning: pa 462 "skipping th 463 fbreak = true; 464 break; 465 } 466 467 //**********Baier-Katkov s 468 469 //back to the local refere 470 tx0 = fCrystalData->AngleX 471 ty0 = ty; 472 473 dzMeV=momentumDirectionSte 474 475 // accelerations 476 axt=(tx0-scatteringAnglesA 477 ayt=(ty0-scatteringAnglesA 478 479 //the angles vs the photon 480 vxin = tx0-txGamma0; 481 vyin = ty0-tyGamma0; 482 //the angles vs the photon 483 vxno = vxin-scatteringAngl 484 vyno = vyin-scatteringAngl 485 486 //phase difference before 487 faseBefore=omprimed2*(gamm 488 489 faseBeforedz = faseBefore* 490 faseBeforedzd2 = faseBefor 491 fa+=faseBeforedz; // 492 fa1=fa-faseBeforedzd2;// 493 dzmod=2*std::sin(faseBefor 494 495 //phi''/faseBefore^2 496 fa2dfaseBefore2 = omprime* 497 498 //phase difference after s 499 faseAfter=omprimed2*(gamma 500 501 skJ=1/faseAfter-1/faseBefo 502 skIx=vxin/faseAfter-vxno/f 503 504 skIy=vyin/faseAfter-vyno/f 505 506 507 sinfa1 = std::sin(fa1); 508 cosfa1 = std::cos(fa1); 509 510 ss+=sinfa1*skJ;//sum sin i 511 sc+=cosfa1*skJ;//sum cos i 512 ssx+=sinfa1*skIx;// sum si 513 ssy+=sinfa1*skIy;// sum si 514 scx+=cosfa1*skIx;// sum co 515 scy+=cosfa1*skIy;// sum co 516 } 517 518 //only of the trajectory cycle 519 if(!fbreak) 520 { 521 G4double i2=ssx*ssx+scx*sc 522 G4double j2=ss*ss+sc*sc;// 523 524 probabilityPPdz += coefNor 525 (i2*e2p 526 } 527 } 528 529 //filling the CDF of probabilities 530 fPairProductionCDFdz.push_back(fPa 531 //**********Baier-Katkov end 532 533 //accumulation of initial paramete 534 fullVectorEtotal.push_back(twoVect 535 fullVectorX.push_back(twoVectorX); 536 fullVectorY.push_back(twoVectorY); 537 fullVectorTX.push_back(twoVectorTX 538 fullVectorTY.push_back(twoVectorTY 539 } 540 541 //photon mean free path 542 //fPairProductionCDFdz.back() = full p 543 //simulated for the current photon alo 544 G4double lMeanFreePath = 1/fPairProduc 545 546 fEffectiveLrad = 7.*lMeanFreePath/9.;/ 547 548 return lMeanFreePath; 549 } 550 else 551 { 552 //dummy process, does not occur 553 return DBL_MAX; 554 } 555 556 } 557 558 //....oooOO0OOooo........oooOO0OOooo........oo 559 560 G4VParticleChange* G4CoherentPairProduction::P 561 con 562 { 563 //example with no physical sense 564 aParticleChange.Initialize(aTrack); 565 //G4LogicalVolume* aLV = aTrack.GetVolume( 566 567 const G4ParticleDefinition* chargedParticl 568 {G4Electron::Electron(),G4Positron::Po 569 570 // the coordinates of the photon in the lo 571 G4ThreeVector xyzGamma0 = 572 aTrack.GetTouchableHandle()->GetHistor 573 GetTopTransform().TransformPoint(a 574 575 // the coordinates of the photon in the co 576 //a channel (elementary periodic cell) 577 G4ThreeVector xyzGamma = fCrystalData->Coo 578 579 //global time 580 G4double tGlobalGamma = aTrack.GetGlobalTi 581 582 G4double ksi1 = G4UniformRand()*fPairProdu 583 584 //randomly choosing the pair to be produce 585 //according to the probabilities calculate 586 G4int ipair = FindVectorIndex(fPairProduct 587 //a pair produced 588 589 // the coordinates of a charged particle i 590 //a channel (elementary periodic cell) 591 G4ThreeVector xyzparticle; 592 //cycle e- and e+ within single pair 593 for(G4int j=0; j<2;j++) 594 { 595 xyzparticle.set(fullVectorX[ipair][j], 596 597 //in the local reference system of the 598 G4ThreeVector newParticleCoordinateXYZ 599 fCrystalData->CoordinatesFromL 600 //the same in the global reference sys 601 newParticleCoordinateXYZ = 602 aTrack.GetTouchableHandle()->GetHi 603 GetTopTransform().Inverse().Tr 604 605 //back to the local reference system o 606 G4double tx0 = fCrystalData->AngleXFro 607 G4double ty0 = fullVectorTY[ipair][j]; 608 609 G4double momentumDirectionZ = 1./ 610 std::sqr 611 612 613 //momentum direction vector of the cha 614 //in the local reference system of the 615 G4ThreeVector momentumDirectionParticl 616 617 618 //the same in the global reference sys 619 momentumDirectionParticle = 620 (aTrack.GetTouchableHandle()->GetH 621 momentumDirectionParticle; 622 623 G4DynamicParticle* chargedParticle = 624 new G4DynamicParticle(chargedParti 625 momentumDirectio 626 fullVectorEtotal 627 628 // Create the track for the secondary 629 G4Track* secondaryTrack = new G4Track( 630 631 632 secondaryTrack->SetTouchableHandle(aSt 633 secondaryTrack->SetParentID(aTrack.Get 634 635 //generation of a secondary charged pa 636 aParticleChange.AddSecondary(secondary 637 } 638 639 //killing the photon 640 aParticleChange.ProposeTrackStatus(fStopAn 641 642 return &aParticleChange; 643 } 644 645 //....oooOO0OOooo........oooOO0OOooo........oo 646 647 G4int G4CoherentPairProduction::FindVectorInde 648 { 649 auto iteratorbegin = myvector.begin(); 650 auto iteratorend = myvector.end(); 651 652 //vector index (for non precise values low 653 auto loweriterator = std::lower_bound(iter 654 //return the index of the vector element 655 return (G4int)std::distance(iteratorbegin, 656 } 657 658 //....oooOO0OOooo........oooOO0OOooo........oo 659 660 void G4CoherentPairProduction::Input(const G4M 661 const G4S 662 const G4S 663 { 664 //initializing the class with containing a 665 //the crystal material and crystal lattice 666 //Channeling scattering and ionization pro 667 fCrystalData = new G4ChannelingFastSimCrys 668 //setting all the crystal material and lat 669 fCrystalData->SetMaterialProperties(crysta 670 } 671 672 //....oooOO0OOooo........oooOO0OOooo........oo 673 674 void G4CoherentPairProduction::Input(const G4C 675 { 676 //setting the class with containing all 677 //the crystal material and crystal lattice 678 //Channeling scattering and ionization pro 679 //fCrystalData = new G4ChannelingFastSimCr 680 681 fCrystalData = const_cast<G4ChannelingFast 682 } 683 684 //....oooOO0OOooo........oooOO0OOooo........oo 685 686 void G4CoherentPairProduction::ProcessDescript 687 { 688 out << " Coherent pair production"; 689 G4VDiscreteProcess::ProcessDescription(out 690 } 691 692 //....oooOO0OOooo........oooOO0OOooo........oo 693