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 // The lust update: M.V. Kossov, CERN/ITEP(Mos 28 // 29 // 30 // G4 Physics class: G4ChipsPionMinusInelastic 31 // Created: M.V. Kossov, CERN/ITEP(Moscow), 20 32 // The last update: M.V. Kossov, CERN/ITEP (Mo 33 // 34 // ------------------------------------------- 35 // Short description: Cross-sections extracted 36 // pion interactions. Original author: M. Koss 37 // ------------------------------------------- 38 // 39 40 #include "G4ChipsPionMinusInelasticXS.hh" 41 #include "G4ChipsPionPlusInelasticXS.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4DynamicParticle.hh" 44 #include "G4ParticleDefinition.hh" 45 #include "G4PionMinus.hh" 46 47 #include "G4Log.hh" 48 #include "G4Exp.hh" 49 50 // factory 51 #include "G4CrossSectionFactory.hh" 52 // 53 G4_DECLARE_XS_FACTORY(G4ChipsPionMinusInelasti 54 55 G4ChipsPionMinusInelasticXS::G4ChipsPionMinusI 56 { 57 // Initialization of the 58 lastLEN=0; // Pointer to lastArray of LowEn 59 lastHEN=0; // Pointer to lastArray of HighEn 60 lastN=0; // The last N of calculated nucle 61 lastZ=0; // The last Z of calculated nucle 62 lastP=0.; // Last used cross section Moment 63 lastTH=0.; // Last threshold momentum 64 lastCS=0.; // Last value of the Cross Sectio 65 lastI=0; // The last position in the DAMDB 66 LEN = new std::vector<G4double*>; 67 HEN = new std::vector<G4double*>; 68 } 69 70 G4ChipsPionMinusInelasticXS::~G4ChipsPionMinus 71 { 72 std::size_t lens=LEN->size(); 73 for(std::size_t i=0; i<lens; ++i) delete[] ( 74 delete LEN; 75 std::size_t hens=HEN->size(); 76 for(std::size_t i=0; i<hens; ++i) delete[] ( 77 delete HEN; 78 } 79 80 void 81 G4ChipsPionMinusInelasticXS::CrossSectionDescr 82 { 83 outFile << "G4ChipsPionMinusInelasticXS pr 84 << "section for pion- nucleus scat 85 << "momentum. The cross section is 86 << "CHIPS parameterization of cros 87 } 88 89 G4bool G4ChipsPionMinusInelasticXS::IsIsoAppli 90 const G4Element*, 91 const G4Material*) 92 { 93 return true; 94 } 95 96 // The main member function giving the collisi 97 // Make pMom in independent units ! (Now it is 98 G4double G4ChipsPionMinusInelasticXS::GetIsoCr 99 const G4Isotope*, 100 const G4Element*, 101 const G4Material*) 102 { 103 G4double pMom=Pt->GetTotalMomentum(); 104 G4int tgN = A - tgZ; 105 106 return GetChipsCrossSection(pMom, tgZ, tgN, 107 } 108 109 G4double G4ChipsPionMinusInelasticXS::GetChips 110 { 111 112 G4bool in=false; // By d 113 if(tgN!=lastN || tgZ!=lastZ) // The 114 { 115 in = false; // By d 116 lastP = 0.; // New 117 lastN = tgN; // The 118 lastZ = tgZ; // The 119 lastI = (G4int)colN.size(); // Size 120 j = 0; // A#0f 121 if(lastI) for(G4int i=0; i<lastI; ++i) // 122 { 123 if(colN[i]==tgN && colZ[i]==tgZ) // Try 124 { 125 lastI=i; // Reme 126 lastTH =colTH[i]; // The 127 if(pMom<=lastTH) 128 { 129 return 0.; // Ener 130 } 131 lastP =colP [i]; // Last 132 lastCS =colCS[i]; // Last 133 in = true; // This 134 // Momentum pMom is in IU ! @@ Units 135 lastCS=CalculateCrossSection(-1,j,-211 136 if(lastCS<=0. && pMom>lastTH) // Corr 137 { 138 lastCS=0.; 139 lastTH=pMom; 140 } 141 break; // Go o 142 } 143 j++; // Incr 144 } 145 if(!in) // This 146 { 147 //!!The slave functions must provide cro 148 lastCS=CalculateCrossSection(0,j,-211,la 149 //if(lastCS>0.) // It 150 //{ 151 152 lastTH = 0; //ThresholdEnergy(tgZ, tgN); 153 colN.push_back(tgN); 154 colZ.push_back(tgZ); 155 colP.push_back(pMom); 156 colTH.push_back(lastTH); 157 colCS.push_back(lastCS); 158 //} // M.K. Presence of H1 with high thr 159 return lastCS*millibarn; 160 } // End of creation of the new set of par 161 else 162 { 163 colP[lastI]=pMom; 164 colCS[lastI]=lastCS; 165 } 166 } // End of parameters udate 167 else if(pMom<=lastTH) 168 { 169 return 0.; // Mome 170 } 171 else // It i 172 { 173 lastCS=CalculateCrossSection(1,j,-211,last 174 lastP=pMom; 175 } 176 return lastCS*millibarn; 177 } 178 179 // The main member function giving the gamma-A 180 G4double G4ChipsPionMinusInelasticXS::Calculat 181 G4int, 182 { 183 static const G4double THmin=27.; // defa 184 static const G4double THmiG=THmin*.001; // m 185 static const G4double dP=10.; // step 186 static const G4double dPG=dP*.001; // step 187 static const G4int nL=105; // A#of 188 static const G4double Pmin=THmin+(nL-1)*dP; 189 static const G4double Pmax=227000.; // maxP 190 static const G4int nH=224; // A#of 191 static const G4double milP=G4Log(Pmin);// Lo 192 static const G4double malP=G4Log(Pmax);// Hi 193 static const G4double dlP=(malP-milP)/(nH-1) 194 static const G4double milPG=G4Log(.001*Pmin) 195 if(F<=0) // This 196 { 197 if(F<0) // This 198 { 199 G4int sync=(G4int)LEN->size(); 200 if(sync<=I) G4cerr<<"*!*G4ChipsPiMinusNu 201 lastLEN=(*LEN)[I]; // Poin 202 lastHEN=(*HEN)[I]; // Poin 203 } 204 else // This 205 { 206 lastLEN = new G4double[nL]; // Allo 207 lastHEN = new G4double[nH]; // Allo 208 // --- Instead of making a separate func 209 G4double P=THmiG; // Tabl 210 for(G4int k=0; k<nL; k++) 211 { 212 lastLEN[k] = CrossSectionLin(targZ, ta 213 P+=dPG; 214 } 215 G4double lP=milPG; 216 for(G4int n=0; n<nH; n++) 217 { 218 lastHEN[n] = CrossSectionLog(targZ, ta 219 lP+=dlP; 220 } 221 // --- End of possible separate function 222 // *** The synchronization check *** 223 G4int sync=(G4int)LEN->size(); 224 if(sync!=I) 225 { 226 G4cerr<<"***G4ChipsPiMinusNuclCS::Calc 227 <<", N="<<targN<<", F="<<F<<G4en 228 //G4Exception("G4PiMinusNuclearCS::Cal 229 } 230 LEN->push_back(lastLEN); // reme 231 HEN->push_back(lastHEN); // reme 232 } // End of creation of the new set of par 233 } // End of parameters udate 234 // =---------------------= NOW the Magic For 235 G4double sigma; 236 if (Momentum<lastTH) return 0.; // It m 237 else if (Momentum<Pmin) // High 238 { 239 sigma=EquLinearFit(Momentum,nL,THmin,dP,la 240 } 241 else if (Momentum<Pmax) // High 242 { 243 G4double lP=G4Log(Momentum); 244 sigma=EquLinearFit(lP,nH,milP,dlP,lastHEN) 245 } 246 else // UHE 247 { 248 G4double P=0.001*Momentum; // Appr 249 sigma=CrossSectionFormula(targZ, targN, P, 250 } 251 if(sigma<0.) return 0.; 252 return sigma; 253 } 254 255 // Calculation formula for piMinus-nuclear ine 256 G4double G4ChipsPionMinusInelasticXS::CrossSec 257 { 258 G4double lP=G4Log(P); 259 return CrossSectionFormula(tZ, tN, P, lP); 260 } 261 262 // Calculation formula for piMinus-nuclear ine 263 G4double G4ChipsPionMinusInelasticXS::CrossSec 264 { 265 G4double P=G4Exp(lP); 266 return CrossSectionFormula(tZ, tN, P, lP); 267 } 268 // Calculation formula for piMinus-nuclear ine 269 G4double G4ChipsPionMinusInelasticXS::CrossSec 270 271 { 272 G4double sigma=0.; 273 if(tZ==1 && !tN) // P 274 { 275 G4double lr=lP+1.27; // 276 G4double LE=1.53/(lr*lr+.0676); // 277 G4double ld=lP-3.5; 278 G4double ld2=ld*ld; 279 G4double p2=P*P; 280 G4double p4=p2*p2; 281 G4double sp=std::sqrt(P); 282 G4double lm=lP+.36; 283 G4double md=lm*lm+.04; 284 G4double lh=lP-.017; 285 G4double hd=lh*lh+.0025; 286 G4double El=(.0557*ld2+2.4+7./sp)/(1.+.7/p 287 G4double To=(.3*ld2+22.3+12./sp)/(1.+.4/p4 288 sigma=(To-El)+.4/md+.01/hd; 289 sigma+=LE*2; // 290 } 291 else if(tZ==1 && tN==1) // 292 { 293 G4double p2=P*P; 294 G4double d=lP-2.7; 295 G4double f=lP+1.25; 296 G4double gg=lP-.017; 297 sigma=(.55*d*d+38.+23./std::sqrt(P))/(1.+. 298 } 299 else if(tZ<97 && tN<152) // G 300 { 301 G4double d=lP-4.2; 302 G4double p2=P*P; 303 G4double p4=p2*p2; 304 G4double a=tN+tZ; // A 305 G4double al=G4Log(a); 306 G4double sa=std::sqrt(a); 307 G4double ssa=std::sqrt(sa); 308 G4double a2=a*a; 309 G4double c=41.*G4Exp(al*.68)*(1.+44./a2)/( 310 G4double f=120*sa/(1.+24./a/ssa); 311 G4double gg=-1.32-al*.043; 312 G4double u=lP-gg; 313 G4double h=al*(.388-.046*al); 314 sigma=(c+d*d)/(1.+.17/p4)+f/(u*u+h*h); 315 } 316 else 317 { 318 G4cerr<<"-Warning-G4ChipsPiMinusNuclearCro 319 sigma=0.; 320 } 321 if(sigma<0.) return 0.; 322 return sigma; 323 } 324 325 G4double G4ChipsPionMinusInelasticXS::EquLinea 326 { 327 if(DX<=0. || N<2) 328 { 329 G4cerr<<"***G4ChipsPionMinusInelasticXS: 330 return Y[0]; 331 } 332 333 G4int N2=N-2; 334 G4double d=(X-X0)/DX; 335 G4int jj=static_cast<int>(d); 336 if (jj<0) jj=0; 337 else if(jj>N2) jj=N2; 338 d-=jj; // excess 339 G4double yi=Y[jj]; 340 G4double sigma=yi+(Y[jj+1]-yi)*d; 341 342 return sigma; 343 } 344