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 // 27 28 #include "globals.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4CollisionComposite.hh" 31 #include "G4VCollision.hh" 32 #include "G4CollisionVector.hh" 33 #include "G4KineticTrack.hh" 34 #include "G4KineticTrackVector.hh" 35 #include "G4VCrossSectionSource.hh" 36 #include "G4HadTmpUtil.hh" 37 #include "G4AutoLock.hh" 38 39 const G4int G4CollisionComposite::nPoints = 32 40 41 const G4double G4CollisionComposite::theT[nPoi 42 {.01, .03, .05, .1, .15, .2, .3, .4, .5, .6, . 43 44 G4CollisionComposite::G4CollisionComposite() 45 { 46 G4MUTEXINIT( bufferMutex ); 47 } 48 49 50 G4CollisionComposite::~G4CollisionComposite() 51 { 52 G4MUTEXDESTROY(bufferMutex); 53 std::for_each(components.begin(), components 54 } 55 56 57 G4double G4CollisionComposite::CrossSection(co 58 const G4KineticTrack& trk2) cons 59 { 60 G4double crossSect = 0.; 61 const G4VCrossSectionSource* xSource = GetCr 62 if (xSource != 0) 63 // There is a total cross section for this C 64 { 65 crossSect = xSource->CrossSection(trk1,trk 66 } 67 else 68 { 69 G4AutoLock l(&bufferMutex); 70 // waiting for mutable to enable buffering 71 const_cast<G4CollisionComposite *>(this)-> 72 // G4cerr << "Buffer filled, reying with sq 73 crossSect = BufferedCrossSection(trk1,trk2 74 } 75 return crossSect; 76 } 77 78 79 G4KineticTrackVector* G4CollisionComposite::Fi 80 const G4KineticTrack& trk2) co 81 { 82 std::vector<G4double> cxCache; 83 G4double partialCxSum = 0.0; 84 85 size_t i; 86 for (i=0; i<components.size(); i++) 87 { 88 G4double partialCx; 89 // cout << "comp" << i << " " << components 90 if (components[i]->IsInCharge(trk1,trk2)) 91 { 92 partialCx = components[i]->CrossSection( 93 } 94 else 95 { 96 partialCx = 0.0; 97 } 98 // cout << " cx=" << partialCx << endl; 99 partialCxSum += partialCx; 100 cxCache.push_back(partialCx); 101 } 102 103 G4double random = G4UniformRand()*partialCxS 104 G4double running = 0; 105 for (i=0; i<cxCache.size(); i++) 106 { 107 running += cxCache[i]; 108 if (running > random) 109 { 110 return components[i]->FinalState(trk1, t 111 } 112 } 113 // G4cerr <<"in charge = "<<IsInCharge(trk1, 114 // G4cerr <<"Cross-section = "<<CrossSection( 115 // G4cerr <<"Names = "<<trk1.GetDefinition()- 116 // throw G4HadronicException(__FILE__, __LINE 117 return NULL; 118 } 119 120 121 G4bool G4CollisionComposite::IsInCharge(const 122 const G4KineticTrack& trk2) const 123 { 124 G4bool isInCharge = false; 125 126 // The composite is in charge if any of its 127 128 const G4CollisionVector* comps = GetComponen 129 if (comps) 130 { 131 G4CollisionVector::const_iterator iter; 132 for (iter = comps->begin(); iter != comp 133 { 134 if ( ((*iter))->IsInCharge(trk1,trk2) ) isI 135 } 136 } 137 138 return isInCharge; 139 } 140 141 void G4CollisionComposite:: 142 BufferCrossSection(const G4ParticleDefinition 143 { 144 // check if already buffered 145 size_t i; 146 for(i=0; i<theBuffer.size(); i++) 147 { 148 if(theBuffer[i].InCharge(aP, bP)) return; 149 } 150 // G4cerr << "Buffering for "<<aP->GetPartic 151 152 // buffer the new one. 153 G4CrossSectionBuffer aNewBuff(aP, bP); 154 size_t maxE=nPoints; 155 for(size_t tt=0; tt<maxE; tt++) 156 { 157 G4double aT = theT[tt]*GeV; 158 G4double crossSect = 0; 159 // The total cross-section is summed over 160 161 //A.R. 28-Sep-2012 Fix reproducibility pr 162 // Assign the kinetic ene 163 // two particles, instead 164 G4double atime = 0; 165 G4double btime = 0; 166 G4ThreeVector aPosition(0,0,0); 167 G4ThreeVector bPosition(0,0,0); 168 G4double aM = aP->GetPDGMass(); 169 G4double bM = bP->GetPDGMass(); 170 G4double aE = aM; 171 G4double bE = bM; 172 G4ThreeVector aMom(0,0,0); 173 G4ThreeVector bMom(0,0,0); 174 if ( aM <= bM ) { 175 aE += aT; 176 aMom = G4ThreeVector(0,0,std::sqrt(aE*aE 177 } else { 178 bE += aT; 179 bMom = G4ThreeVector(0,0,std::sqrt(bE*bE 180 } 181 G4LorentzVector a4Momentum(aE, aMom); 182 G4LorentzVector b4Momentum(bE, bMom); 183 G4KineticTrack a(aP, atime, aPosition, a4 184 G4KineticTrack b(bP, btime, bPosition, b4 185 186 for (i=0; i<components.size(); i++) 187 { 188 if(components[i]->IsInCharge(a,b)) 189 { 190 crossSect += components[i]->CrossSection(a, 191 } 192 } 193 G4double sqrts = (a4Momentum+b4Momentum). 194 aNewBuff.push_back(sqrts, crossSect); 195 } 196 theBuffer.push_back(std::move(aNewBuff)); 197 // theBuffer.back().Print(); 198 } 199 200 201 G4double G4CollisionComposite:: 202 BufferedCrossSection(const G4KineticTrack& trk 203 { 204 for(size_t i=0; i<theBuffer.size(); i++) 205 { 206 if(theBuffer[i].InCharge(trk1.GetDefiniti 207 { 208 return theBuffer[i].CrossSection(trk1, 209 } 210 } 211 throw G4HadronicException(__FILE__, __LINE_ 212 return 0; 213 } 214 215