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