Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4ChipsHyperonElasticXS.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/cross_sections/src/G4ChipsHyperonElasticXS.cc (Version 11.3.0) and /processes/hadronic/cross_sections/src/G4ChipsHyperonElasticXS.cc (Version 9.5.p1)


  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 //                                                
 29 // G4 Physics class: G4ChipsHyperonElasticXS f    
 30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 5-    
 31 // The last update: M.V. Kossov, CERN/ITEP (Mo    
 32 //                                                
 33 // -------------------------------------------    
 34 // Short description: Interaction cross-sectio    
 35 // Class extracted from CHIPS and integrated i    
 36 // -------------------------------------------    
 37 //                                                
 38                                                   
 39 #include "G4ChipsHyperonElasticXS.hh"             
 40 #include "G4SystemOfUnits.hh"                     
 41 #include "G4DynamicParticle.hh"                   
 42 #include "G4ParticleDefinition.hh"                
 43 #include "G4Lambda.hh"                            
 44 #include "G4SigmaPlus.hh"                         
 45 #include "G4SigmaMinus.hh"                        
 46 #include "G4SigmaZero.hh"                         
 47 #include "G4XiMinus.hh"                           
 48 #include "G4XiZero.hh"                            
 49 #include "G4OmegaMinus.hh"                        
 50 #include "G4Nucleus.hh"                           
 51 #include "G4ParticleTable.hh"                     
 52 #include "G4NucleiProperties.hh"                  
 53 #include "G4IonTable.hh"                          
 54 #include "G4Exp.hh"                               
 55 #include "G4Log.hh"                               
 56 #include "G4Pow.hh"                               
 57                                                   
 58 // factory                                        
 59 #include "G4CrossSectionFactory.hh"               
 60 //                                                
 61 G4_DECLARE_XS_FACTORY(G4ChipsHyperonElasticXS)    
 62                                                   
 63 G4ChipsHyperonElasticXS::G4ChipsHyperonElastic    
 64 {                                                 
 65   lPMin=-8.;  //Min tabulatedLogarithmMomentum    
 66   lPMax= 8.;  //Max tabulatedLogarithmMomentum    
 67   dlnP=(lPMax-lPMin)/nLast;// LogStep inTable     
 68   onlyCS=true;//Flag toCalculOnlyCS(not Si/Bi)    
 69   lastSIG=0.; //Last calculated cross section     
 70   lastLP=-10.;//LastLog(mom_of IncidentHadron)    
 71   lastTM=0.;  //Last t_maximum                    
 72   theSS=0.;   //TheLastSqSlope of 1st difr.Max    
 73   theS1=0.;   //TheLastMantissa of 1st difrMax    
 74   theB1=0.;   //TheLastSlope of 1st difructMax    
 75   theS2=0.;   //TheLastMantissa of 2nd difrMax    
 76   theB2=0.;   //TheLastSlope of 2nd difructMax    
 77   theS3=0.;   //TheLastMantissa of 3d difr.Max    
 78   theB3=0.;   //TheLastSlope of 3d difruct.Max    
 79   theS4=0.;   //TheLastMantissa of 4th difrMax    
 80   theB4=0.;   //TheLastSlope of 4th difructMax    
 81   lastTZ=0;   // Last atomic number of the tar    
 82   lastTN=0;   // Last # of neutrons in the tar    
 83   lastPIN=0.; // Last initialized max momentum    
 84   lastCST=0;  // Elastic cross-section table      
 85   lastPAR=0;  // ParametersForFunctionCalculat    
 86   lastSST=0;  // E-dep ofSqardSlope of 1st dif    
 87   lastS1T=0;  // E-dep of mantissa of 1st dif.    
 88   lastB1T=0;  // E-dep of the slope of 1st dif    
 89   lastS2T=0;  // E-dep of mantissa of 2nd difr    
 90   lastB2T=0;  // E-dep of the slope of 2nd dif    
 91   lastS3T=0;  // E-dep of mantissa of 3d difr.    
 92   lastB3T=0;  // E-dep of the slope of 3d difr    
 93   lastS4T=0;  // E-dep of mantissa of 4th difr    
 94   lastB4T=0;  // E-dep of the slope of 4th dif    
 95   lastN=0;    // The last N of calculated nucl    
 96   lastZ=0;    // The last Z of calculated nucl    
 97   lastP=0.;   // LastUsed inCrossSection Momen    
 98   lastTH=0.;  // Last threshold momentum          
 99   lastCS=0.;  // Last value of the Cross Secti    
100   lastI=0;    // The last position in the DAMD    
101 }                                                 
102                                                   
103 G4ChipsHyperonElasticXS::~G4ChipsHyperonElasti    
104 {                                                 
105   std::vector<G4double*>::iterator pos;           
106   for (pos=CST.begin(); pos<CST.end(); pos++)     
107   { delete [] *pos; }                             
108   CST.clear();                                    
109   for (pos=PAR.begin(); pos<PAR.end(); pos++)     
110   { delete [] *pos; }                             
111   PAR.clear();                                    
112   for (pos=SST.begin(); pos<SST.end(); pos++)     
113   { delete [] *pos; }                             
114   SST.clear();                                    
115   for (pos=S1T.begin(); pos<S1T.end(); pos++)     
116   { delete [] *pos; }                             
117   S1T.clear();                                    
118   for (pos=B1T.begin(); pos<B1T.end(); pos++)     
119   { delete [] *pos; }                             
120   B1T.clear();                                    
121   for (pos=S2T.begin(); pos<S2T.end(); pos++)     
122   { delete [] *pos; }                             
123   S2T.clear();                                    
124   for (pos=B2T.begin(); pos<B2T.end(); pos++)     
125   { delete [] *pos; }                             
126   B2T.clear();                                    
127   for (pos=S3T.begin(); pos<S3T.end(); pos++)     
128   { delete [] *pos; }                             
129   S3T.clear();                                    
130   for (pos=B3T.begin(); pos<B3T.end(); pos++)     
131   { delete [] *pos; }                             
132   B3T.clear();                                    
133   for (pos=S4T.begin(); pos<S4T.end(); pos++)     
134   { delete [] *pos; }                             
135   S4T.clear();                                    
136   for (pos=B4T.begin(); pos<B4T.end(); pos++)     
137   { delete [] *pos; }                             
138   B4T.clear();                                    
139 }                                                 
140                                                   
141 void                                              
142 G4ChipsHyperonElasticXS::CrossSectionDescripti    
143 {                                                 
144     outFile << "G4ChipsHyperonElasticXS provid    
145             << "section for hyperon nucleus sc    
146             << "momentum. The cross section is    
147             << "CHIPS parameterization of cros    
148 }                                                 
149                                                   
150 G4bool G4ChipsHyperonElasticXS::IsIsoApplicabl    
151          const G4Element*,                        
152          const G4Material*)                       
153 {                                                 
154   return true;                                    
155 }                                                 
156                                                   
157 // The main member function giving the collisi    
158 // Make pMom in independent units ! (Now it is    
159 G4double G4ChipsHyperonElasticXS::GetIsoCrossS    
160                const G4Isotope*,                  
161                const G4Element*,                  
162                const G4Material*)                 
163 {                                                 
164   G4double pMom=Pt->GetTotalMomentum();           
165   G4int tgN = A - tgZ;                            
166   G4int pdg = Pt->GetDefinition()->GetPDGEncod    
167                                                   
168   return GetChipsCrossSection(pMom, tgZ, tgN,     
169 }                                                 
170                                                   
171 G4double G4ChipsHyperonElasticXS::GetChipsCros    
172 {                                                 
173                                                   
174   G4bool fCS = false;                             
175   G4double pEn=pMom;                              
176                                                   
177   onlyCS=fCS;                                     
178                                                   
179   G4bool in=false;                   // By def    
180   lastP   = 0.;                      // New mo    
181   lastN   = tgN;                     // The la    
182   lastZ   = tgZ;                     // The la    
183   lastI   = (G4int)colN.size();      // Size o    
184   if(lastI) for(G4int i=0; i<lastI; ++i) // Lo    
185   {                                  // The nu    
186     if(colN[i]==tgN && colZ[i]==tgZ) // Isotop    
187     {                                             
188       lastI=i;                                    
189       lastTH =colTH[i];              // Last T    
190       if(pEn<=lastTH)                             
191       {                                           
192         return 0.;                   // Energy    
193       }                                           
194       lastP  =colP [i];              // Last M    
195       lastCS =colCS[i];              // Last C    
196       //  if(std::fabs(lastP/pMom-1.)<toleranc    
197       if(lastP == pMom)              // Do not    
198       {                                           
199         CalculateCrossSection(fCS,-1,i,pPDG,la    
200         return lastCS*millibarn;     // Use th    
201       }                                           
202       in = true;                       // This    
203       // Momentum pMom is in IU ! @@ Units        
204       lastCS=CalculateCrossSection(fCS,-1,i,pP    
205       if(lastCS<=0. && pEn>lastTH)    // Corre    
206       {                                           
207         lastTH=pEn;                               
208       }                                           
209       break;                           // Go o    
210     }                                             
211   } // End of attampt to find the nucleus in D    
212   if(!in)                            // This n    
213   {                                               
214     //!!The slave functions must provide cross    
215     lastCS=CalculateCrossSection(fCS,0,lastI,p    
216     if(lastCS<=0.)                                
217     {                                             
218       lastTH = 0; //ThresholdEnergy(tgZ, tgN);    
219       if(pEn>lastTH)                              
220       {                                           
221         lastTH=pEn;                               
222       }                                           
223     }                                             
224     colN.push_back(tgN);                          
225     colZ.push_back(tgZ);                          
226     colP.push_back(pMom);                         
227     colTH.push_back(lastTH);                      
228     colCS.push_back(lastCS);                      
229     return lastCS*millibarn;                      
230   } // End of creation of the new set of param    
231   else                                            
232   {                                               
233     colP[lastI]=pMom;                             
234     colCS[lastI]=lastCS;                          
235   }                                               
236   return lastCS*millibarn;                        
237 }                                                 
238                                                   
239 // Calculation of total elastic cross section     
240 // F=0 - create AMDB, F=-1 - read&update AMDB,    
241 G4double G4ChipsHyperonElasticXS::CalculateCro    
242                                              G    
243 {                                                 
244   G4double pMom=pIU/GeV;                // All    
245   onlyCS=CS;                            // Fla    
246   lastLP=G4Log(pMom);                // Make a    
247   if(F)                                 // Thi    
248   {                                               
249     if(F<0)                             // the    
250     {                                             
251       lastPIN = PIN[I];                 // Max    
252       lastPAR = PAR[I];                 // Poi    
253                                                   
254       lastCST = CST[I];                 // Poi    
255       lastSST = SST[I];                 // Poi    
256       lastS1T = S1T[I];                 // Poi    
257       lastB1T = B1T[I];                 // Poi    
258       lastS2T = S2T[I];                 // Poi    
259       lastB2T = B2T[I];                 // Poi    
260       lastS3T = S3T[I];                 // Poi    
261       lastB3T = B3T[I];                 // Poi    
262       lastS4T = S4T[I];                 // Poi    
263       lastB4T = B4T[I];                 // Poi    
264     }                                             
265     if(lastLP>lastPIN && lastLP<lPMax)            
266     {                                             
267       lastPIN=GetPTables(lastLP,lastPIN,PDG,tg    
268       PIN[I]=lastPIN;                   // Rem    
269     }                                             
270   }                                               
271   else                                  // Thi    
272   {                                               
273     lastPAR = new G4double[nPoints];    // All    
274     lastPAR[nLast]=0;                   // Ini    
275     lastCST = new G4double[nPoints];    // All    
276     lastSST = new G4double[nPoints];    // All    
277     lastS1T = new G4double[nPoints];    // All    
278     lastB1T = new G4double[nPoints];    // All    
279     lastS2T = new G4double[nPoints];    // All    
280     lastB2T = new G4double[nPoints];    // All    
281     lastS3T = new G4double[nPoints];    // All    
282     lastB3T = new G4double[nPoints];    // All    
283     lastS4T = new G4double[nPoints];    // All    
284     lastB4T = new G4double[nPoints];    // All    
285     lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,    
286     PIN.push_back(lastPIN);             // Fil    
287     PAR.push_back(lastPAR);             // Fil    
288     CST.push_back(lastCST);             // Fil    
289     SST.push_back(lastSST);             // Fil    
290     S1T.push_back(lastS1T);             // Fil    
291     B1T.push_back(lastB1T);             // Fil    
292     S2T.push_back(lastS2T);             // Fil    
293     B2T.push_back(lastB2T);             // Fil    
294     S3T.push_back(lastS3T);             // Fil    
295     B3T.push_back(lastB3T);             // Fil    
296     S4T.push_back(lastS4T);             // Fil    
297     B4T.push_back(lastB4T);             // Fil    
298   } // End of creation/update of the new set o    
299   // =-----------= NOW Update (if necessary) a    
300   if(lastLP>lastPIN && lastLP<lPMax)              
301   {                                               
302     lastPIN = GetPTables(lastLP,lastPIN,PDG,tg    
303   }                                               
304   if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, p    
305   if(lastLP>lPMin && lastLP<=lastPIN)   // Lin    
306   {                                               
307     if(lastLP==lastPIN)                           
308     {                                             
309       G4double shift=(lastLP-lPMin)/dlnP+.0000    
310       G4int    blast=static_cast<int>(shift);     
311       if(blast<0 || blast>=nLast)G4cout<<"G4QH    
312       lastSIG = lastCST[blast];                   
313       if(!onlyCS)                       // Ski    
314       {                                           
315         theSS  = lastSST[blast];                  
316         theS1  = lastS1T[blast];                  
317         theB1  = lastB1T[blast];                  
318         theS2  = lastS2T[blast];                  
319         theB2  = lastB2T[blast];                  
320         theS3  = lastS3T[blast];                  
321         theB3  = lastB3T[blast];                  
322         theS4  = lastS4T[blast];                  
323         theB4  = lastB4T[blast];                  
324       }                                           
325     }                                             
326     else                                          
327     {                                             
328       G4double shift=(lastLP-lPMin)/dlnP;         
329       G4int    blast=static_cast<int>(shift);     
330       if(blast<0)   blast=0;                      
331       if(blast>=nLast) blast=nLast-1;             
332       shift-=blast;                               
333       G4int lastL=blast+1;                        
334       G4double SIGL=lastCST[blast];               
335       lastSIG= SIGL+shift*(lastCST[lastL]-SIGL    
336       if(!onlyCS)                       // Ski    
337       {                                           
338         G4double SSTL=lastSST[blast];             
339         theSS=SSTL+shift*(lastSST[lastL]-SSTL)    
340         G4double S1TL=lastS1T[blast];             
341         theS1=S1TL+shift*(lastS1T[lastL]-S1TL)    
342         G4double B1TL=lastB1T[blast];             
343         theB1=B1TL+shift*(lastB1T[lastL]-B1TL)    
344         G4double S2TL=lastS2T[blast];             
345         theS2=S2TL+shift*(lastS2T[lastL]-S2TL)    
346         G4double B2TL=lastB2T[blast];             
347         theB2=B2TL+shift*(lastB2T[lastL]-B2TL)    
348         G4double S3TL=lastS3T[blast];             
349         theS3=S3TL+shift*(lastS3T[lastL]-S3TL)    
350         G4double B3TL=lastB3T[blast];             
351         theB3=B3TL+shift*(lastB3T[lastL]-B3TL)    
352         G4double S4TL=lastS4T[blast];             
353         theS4=S4TL+shift*(lastS4T[lastL]-S4TL)    
354         G4double B4TL=lastB4T[blast];             
355         theB4=B4TL+shift*(lastB4T[lastL]-B4TL)    
356       }                                           
357     }                                             
358   }                                               
359   else lastSIG=GetTabValues(lastLP, PDG, tgZ,     
360   if(lastSIG<0.) lastSIG = 0.;                    
361   return lastSIG;                                 
362 }                                                 
363                                                   
364 // It has parameter sets for all tZ/tN/PDG, us    
365 G4double G4ChipsHyperonElasticXS::GetPTables(G    
366                                                   
367 {                                                 
368   // @@ At present all nA==pA ---------> Each     
369   static const G4double pwd=2727;                 
370   const G4int n_hypel=33;                // #o    
371   //                        -0- -1- -2-  -3-      
372   G4double hyp_el[n_hypel]={1.,.002,.12,.0557,    
373                             .001,8.,.055,3.64,    
374                             1.e10,8.5e8,1.e10,    
375   //                        -15--16- -17- -18-    
376   //                         -26-  -27- -28-      
377   //AR-04Jun2014  if(PDG!=3222 && PDG>3000 &&     
378   if(PDG>3000 && PDG<3335)                        
379   {                                               
380     // -- Total pp elastic cross section cs &     
381     //p2=p*p;p3=p2*p;sp=sqrt(p);p2s=p2*sp;lp=l    
382     //CS=2.865/p2s/(1+.0022/p2s)+(18.9+.6461*d    
383     //   par(0)       par(7)     par(1) par(2)    
384     //dl2=lp-5., s1=(74.+3.*dl2*dl2)/(1+3.4/p4    
385     //     par(8) par(9) par(10)        par(11    
386     // b1=8.*p**.055/(1.+3.64/p3); s2=5.e-5+40    
387     // par(15) par(16)  par(17)     par(18) pa    
388     // s3=5.e-5+1.e10/(p4*p4+8.5e8*p2+1.e10);     
389     //  par(24) par(25)     par(26)  par(27) p    
390     //                                            
391     if(lastPAR[nLast]!=pwd) // A unique flag t    
392     {                                             
393       if ( tgZ == 1 && tgN == 0 )                 
394       {                                           
395         for (G4int ip=0; ip<n_hypel; ip++) las    
396       }                                           
397       else                                        
398       {                                           
399         G4double a=tgZ+tgN;                       
400         G4double sa=std::sqrt(a);                 
401         G4double ssa=std::sqrt(sa);               
402         G4double asa=a*sa;                        
403         G4double a2=a*a;                          
404         G4double a3=a2*a;                         
405         G4double a4=a3*a;                         
406         G4double a5=a4*a;                         
407         G4double a6=a4*a2;                        
408         G4double a7=a6*a;                         
409         G4double a8=a7*a;                         
410         G4double a9=a8*a;                         
411         G4double a10=a5*a5;                       
412         G4double a12=a6*a6;                       
413         G4double a14=a7*a7;                       
414         G4double a16=a8*a8;                       
415         G4double a17=a16*a;                       
416         //G4double a20=a16*a4;                    
417         G4double a32=a16*a16;                     
418         // Reaction cross-section parameters (    
419         lastPAR[0]=4./(1.+22/asa);                
420         lastPAR[1]=2.36*asa/(1.+a*.055/ssa);      
421         lastPAR[2]=(1.+.00007*a3/ssa)/(1.+.002    
422         lastPAR[3]=1.76*a/ssa+.00003*a3;          
423         lastPAR[4]=(.03+200./a3)/(1.+1.E5/a3/s    
424         lastPAR[5]=5.;                            
425         lastPAR[6]=0.;                            
426         lastPAR[7]=0.;                            
427         lastPAR[8]=0.;                            
428         // @@ the differential cross-section i    
429         if(a<6.5)                                 
430         {                                         
431           G4double a28=a16*a12;                   
432           // The main pre-exponent      (pel_s    
433           lastPAR[ 9]=4000*a;                     
434           lastPAR[10]=1.2e7*a8+380*a17;           
435           lastPAR[11]=.7/(1.+4.e-12*a16);         
436           lastPAR[12]=2.5/a8/(a4+1.e-16*a32);     
437           lastPAR[13]=.28*a;                      
438           lastPAR[14]=1.2*a2+2.3;                 
439           lastPAR[15]=3.8/a;                      
440           // The main slope             (pel_s    
441           lastPAR[16]=.01/(1.+.0024*a5);          
442           lastPAR[17]=.2*a;                       
443           lastPAR[18]=9.e-7/(1.+.035*a5);         
444           lastPAR[19]=(42.+2.7e-11*a16)/(1.+.1    
445           // The main quadratic         (pel_s    
446           lastPAR[20]=2.25*a3;                    
447           lastPAR[21]=18.;                        
448           lastPAR[22]=2.4e-3*a8/(1.+2.6e-4*a7)    
449           lastPAR[23]=3.5e-36*a32*a8/(1.+5.e-1    
450           // The 1st max pre-exponent   (pel_q    
451           lastPAR[24]=1.e5/(a8+2.5e12/a16);       
452           lastPAR[25]=8.e7/(a12+1.e-27*a28*a28    
453           lastPAR[26]=.0006*a3;                   
454           // The 1st max slope          (pel_q    
455           lastPAR[27]=10.+4.e-8*a12*a;            
456           lastPAR[28]=.114;                       
457           lastPAR[29]=.003;                       
458           lastPAR[30]=2.e-23;                     
459           // The effective pre-exponent (pel_s    
460           lastPAR[31]=1./(1.+.0001*a8);           
461           lastPAR[32]=1.5e-4/(1.+5.e-6*a12);      
462           lastPAR[33]=.03;                        
463           // The effective slope        (pel_s    
464           lastPAR[34]=a/2;                        
465           lastPAR[35]=2.e-7*a4;                   
466           lastPAR[36]=4.;                         
467           lastPAR[37]=64./a3;                     
468           // The gloria pre-exponent    (pel_u    
469           lastPAR[38]=1.e8*G4Exp(.32*asa);        
470           lastPAR[39]=20.*G4Exp(.45*asa);         
471           lastPAR[40]=7.e3+2.4e6/a5;              
472           lastPAR[41]=2.5e5*G4Exp(.085*a3);       
473           lastPAR[42]=2.5*a;                      
474           // The gloria slope           (pel_u    
475           lastPAR[43]=920.+.03*a8*a3;             
476           lastPAR[44]=93.+.0023*a12;              
477         }                                         
478         else                                      
479         {                                         
480           G4double p1a10=2.2e-28*a10;             
481           G4double r4a16=6.e14/a16;               
482           G4double s4a16=r4a16*r4a16;             
483           // a24                                  
484           // a36                                  
485           // The main pre-exponent      (peh_s    
486           lastPAR[ 9]=4.5*G4Pow::GetInstance()    
487           lastPAR[10]=.06*G4Pow::GetInstance()    
488           lastPAR[11]=.6*a/(1.+2.e15/a16);        
489           lastPAR[12]=.17/(a+9.e5/a3+1.5e33/a3    
490           lastPAR[13]=(.001+7.e-11*a5)/(1.+4.4    
491           lastPAR[14]=(p1a10*p1a10+2.e-29)/(1.    
492           // The main slope             (peh_s    
493           lastPAR[15]=400./a12+2.e-22*a9;         
494           lastPAR[16]=1.e-32*a12/(1.+5.e22/a14    
495           lastPAR[17]=1000./a2+9.5*sa*ssa;        
496           lastPAR[18]=4.e-6*a*asa+1.e11/a16;      
497           lastPAR[19]=(120./a+.002*a2)/(1.+2.e    
498           lastPAR[20]=9.+100./a;                  
499           // The main quadratic         (peh_s    
500           lastPAR[21]=.002*a3+3.e7/a6;            
501           lastPAR[22]=7.e-15*a4*asa;              
502           lastPAR[23]=9000./a4;                   
503           // The 1st max pre-exponent   (peh_q    
504           lastPAR[24]=.0011*asa/(1.+3.e34/a32/    
505           lastPAR[25]=1.e-5*a2+2.e14/a16;         
506           lastPAR[26]=1.2e-11*a2/(1.+1.5e19/a1    
507           lastPAR[27]=.016*asa/(1.+5.e16/a16);    
508           // The 1st max slope          (peh_q    
509           lastPAR[28]=.002*a4/(1.+7.e7/G4Pow::    
510           lastPAR[29]=2.e6/a6+7.2/G4Pow::GetIn    
511           lastPAR[30]=11.*a3/(1.+7.e23/a16/a8)    
512           lastPAR[31]=100./asa;                   
513           // The 2nd max pre-exponent   (peh_s    
514           lastPAR[32]=(.1+4.4e-5*a2)/(1.+5.e5/    
515           lastPAR[33]=3.5e-4*a2/(1.+1.e8/a8);     
516           lastPAR[34]=1.3+3.e5/a4;                
517           lastPAR[35]=500./(a2+50.)+3;            
518           lastPAR[36]=1.e-9/a+s4a16*s4a16;        
519           // The 2nd max slope          (peh_s    
520           lastPAR[37]=.4*asa+3.e-9*a6;            
521           lastPAR[38]=.0005*a5;                   
522           lastPAR[39]=.002*a5;                    
523           lastPAR[40]=10.;                        
524           // The effective pre-exponent (peh_u    
525           lastPAR[41]=.05+.005*a;                 
526           lastPAR[42]=7.e-8/sa;                   
527           lastPAR[43]=.8*sa;                      
528           lastPAR[44]=.02*sa;                     
529           lastPAR[45]=1.e8/a3;                    
530           lastPAR[46]=3.e32/(a32+1.e32);          
531           // The effective slope        (peh_u    
532           lastPAR[47]=24.;                        
533           lastPAR[48]=20./sa;                     
534           lastPAR[49]=7.e3*a/(sa+1.);             
535           lastPAR[50]=900.*sa/(1.+500./a3);       
536         }                                         
537         // Parameter for lowEnergyNeutrons        
538         lastPAR[51]=1.e15+2.e27/a4/(1.+2.e-18*    
539       }                                           
540       lastPAR[nLast]=pwd;                         
541       // and initialize the zero element of th    
542       G4double lp=lPMin;                          
543       G4bool memCS=onlyCS;                        
544       onlyCS=false;                               
545       lastCST[0]=GetTabValues(lp, PDG, tgZ, tg    
546       onlyCS=memCS;                               
547       lastSST[0]=theSS;                           
548       lastS1T[0]=theS1;                           
549       lastB1T[0]=theB1;                           
550       lastS2T[0]=theS2;                           
551       lastB2T[0]=theB2;                           
552       lastS3T[0]=theS3;                           
553       lastB3T[0]=theB3;                           
554       lastS4T[0]=theS4;                           
555       lastB4T[0]=theB4;                           
556     }                                             
557     if(LP>ILP)                                    
558     {                                             
559       G4int ini = static_cast<int>((ILP-lPMin+    
560       if(ini<0) ini=0;                            
561       if(ini<nPoints)                             
562       {                                           
563         G4int fin = static_cast<int>((LP-lPMin    
564         if(fin>=nPoints) fin=nLast;               
565         if(fin>=ini)                              
566         {                                         
567           G4double lp=0.;                         
568           for(G4int ip=ini; ip<=fin; ip++)        
569           {                                       
570             lp=lPMin+ip*dlnP;                     
571             G4bool memCS=onlyCS;                  
572             onlyCS=false;                         
573             lastCST[ip]=GetTabValues(lp, PDG,     
574             onlyCS=memCS;                         
575             lastSST[ip]=theSS;                    
576             lastS1T[ip]=theS1;                    
577             lastB1T[ip]=theB1;                    
578             lastS2T[ip]=theS2;                    
579             lastB2T[ip]=theB2;                    
580             lastS3T[ip]=theS3;                    
581             lastB3T[ip]=theB3;                    
582             lastS4T[ip]=theS4;                    
583             lastB4T[ip]=theB4;                    
584           }                                       
585           return lp;                              
586         }                                         
587         else G4cout<<"*Warning*G4ChipsHyperonE    
588                    <<", Z="<<tgZ<<", N="<<tgN<    
589                    <<" > ILP="<<ILP<<" nothing    
590       }                                           
591       else G4cout<<"*Warning*G4ChipsHyperonEla    
592                  <<", Z="<<tgZ<<", N="<<tgN<<"    
593                  <<" > ILP="<<ILP<<", lPMax="<    
594     }                                             
595   } else {                                        
596     // G4cout<<"*Error*G4ChipsHyperonElasticXS    
597     //       <<", N="<<tgN<<", while it is def    
598     // throw G4QException("G4ChipsHyperonElast    
599     G4ExceptionDescription ed;                    
600     ed << "PDG = " << PDG << ", Z = " << tgZ <    
601        << ", while it is defined only for Hype    
602     G4Exception("G4ChipsHyperonElasticXS::GetP    
603                 FatalException, ed);              
604   }                                               
605   return ILP;                                     
606 }                                                 
607                                                   
608 // Returns Q2=-t in independent units (MeV^2)     
609 G4double G4ChipsHyperonElasticXS::GetExchangeT    
610 {                                                 
611   static const G4double GeVSQ=gigaelectronvolt    
612   static const G4double third=1./3.;              
613   static const G4double fifth=1./5.;              
614   static const G4double sevth=1./7.;              
615   //AR-04Jun2014  if(PDG==3222 || PDG<3000 ||     
616   if(PDG<3000 || PDG>3334)G4cout<<"*Warning*G4    
617   if(onlyCS)G4cout<<"*Warning*G4ChipsHyperonEl    
618   if(lastLP<-4.3) return lastTM*GeVSQ*G4Unifor    
619   G4double q2=0.;                                 
620   if(tgZ==1 && tgN==0)                // ===>     
621   {                                               
622     G4double E1=lastTM*theB1;                     
623     G4double R1=(1.-G4Exp(-E1));                  
624     G4double E2=lastTM*theB2;                     
625     G4double R2=(1.-G4Exp(-E2*E2*E2));            
626     G4double E3=lastTM*theB3;                     
627     G4double R3=(1.-G4Exp(-E3));                  
628     G4double I1=R1*theS1/theB1;                   
629     G4double I2=R2*theS2;                         
630     G4double I3=R3*theS3;                         
631     G4double I12=I1+I2;                           
632     G4double rand=(I12+I3)*G4UniformRand();       
633     if     (rand<I1 )                             
634     {                                             
635       G4double ran=R1*G4UniformRand();            
636       if(ran>1.) ran=1.;                          
637       q2=-G4Log(1.-ran)/theB1;                    
638     }                                             
639     else if(rand<I12)                             
640     {                                             
641       G4double ran=R2*G4UniformRand();            
642       if(ran>1.) ran=1.;                          
643       q2=-G4Log(1.-ran);                          
644       if(q2<0.) q2=0.;                            
645       q2=G4Pow::GetInstance()->powA(q2,third)/    
646     }                                             
647     else                                          
648     {                                             
649       G4double ran=R3*G4UniformRand();            
650       if(ran>1.) ran=1.;                          
651       q2=-G4Log(1.-ran)/theB3;                    
652     }                                             
653   }                                               
654   else                                            
655   {                                               
656     G4double a=tgZ+tgN;                           
657     G4double E1=lastTM*(theB1+lastTM*theSS);      
658     G4double R1=(1.-G4Exp(-E1));                  
659     G4double tss=theSS+theSS; // for future so    
660     G4double tm2=lastTM*lastTM;                   
661     G4double E2=lastTM*tm2*theB2;                 
662     if(a>6.5)E2*=tm2;                             
663     G4double R2=(1.-G4Exp(-E2));                  
664     G4double E3=lastTM*theB3;                     
665     if(a>6.5)E3*=tm2*tm2*tm2;                     
666     G4double R3=(1.-G4Exp(-E3));                  
667     G4double E4=lastTM*theB4;                     
668     G4double R4=(1.-G4Exp(-E4));                  
669     G4double I1=R1*theS1;                         
670     G4double I2=R2*theS2;                         
671     G4double I3=R3*theS3;                         
672     G4double I4=R4*theS4;                         
673     G4double I12=I1+I2;                           
674     G4double I13=I12+I3;                          
675     G4double rand=(I13+I4)*G4UniformRand();       
676     if(rand<I1)                                   
677     {                                             
678       G4double ran=R1*G4UniformRand();            
679       if(ran>1.) ran=1.;                          
680       q2=-G4Log(1.-ran)/theB1;                    
681       if(std::fabs(tss)>1.e-7) q2=(std::sqrt(t    
682     }                                             
683     else if(rand<I12)                             
684     {                                             
685       G4double ran=R2*G4UniformRand();            
686       if(ran>1.) ran=1.;                          
687       q2=-G4Log(1.-ran)/theB2;                    
688       if(q2<0.) q2=0.;                            
689       if(a<6.5) q2=G4Pow::GetInstance()->powA(    
690       else      q2=G4Pow::GetInstance()->powA(    
691     }                                             
692     else if(rand<I13)                             
693     {                                             
694       G4double ran=R3*G4UniformRand();            
695       if(ran>1.) ran=1.;                          
696       q2=-G4Log(1.-ran)/theB3;                    
697       if(q2<0.) q2=0.;                            
698       if(a>6.5) q2=G4Pow::GetInstance()->powA(    
699     }                                             
700     else                                          
701     {                                             
702       G4double ran=R4*G4UniformRand();            
703       if(ran>1.) ran=1.;                          
704       q2=-G4Log(1.-ran)/theB4;                    
705       if(a<6.5) q2=lastTM-q2;                     
706     }                                             
707   }                                               
708   if(q2<0.) q2=0.;                                
709   if(!(q2>=-1.||q2<=1.))G4cout<<"*NAN*G4QHyEla    
710   if(q2>lastTM)                                   
711   {                                               
712     q2=lastTM;                                    
713   }                                               
714   return q2*GeVSQ;                                
715 }                                                 
716                                                   
717 // Returns B in independent units (MeV^-2) (al    
718 G4double G4ChipsHyperonElasticXS::GetSlope(G4i    
719 {                                                 
720   static const G4double GeVSQ=gigaelectronvolt    
721   if(onlyCS)G4cout<<"*Warning*G4ChipsHyperonEl    
722   if(lastLP<-4.3) return 0.;          // S-wav    
723   //AR-04Jun2014  if(PDG==3222 || PDG<3000 ||     
724   if(PDG<3000 || PDG>3334)                        
725   {                                               
726     // G4cout<<"*Error*G4ChipsHyperonElasticXS    
727     //       <<", N="<<tgN<<", while it is def    
728     // throw G4QException("G4ChipsHyperonElast    
729     G4ExceptionDescription ed;                    
730     ed << "PDG = " << PDG << ", Z = " << tgZ <    
731        << ", while it is defined only for Hype    
732     G4Exception("G4ChipsHyperonElasticXS::GetS    
733                 FatalException, ed);              
734   }                                               
735   if(theB1<0.) theB1=0.;                          
736   if(!(theB1>=-1.||theB1<=1.)) G4cout<<"*NAN*G    
737   return theB1/GeVSQ;                             
738 }                                                 
739                                                   
740 // Returns half max(Q2=-t) in independent unit    
741 G4double G4ChipsHyperonElasticXS::GetHMaxT()      
742 {                                                 
743   static const G4double HGeVSQ=gigaelectronvol    
744   return lastTM*HGeVSQ;                           
745 }                                                 
746                                                   
747 // lastLP is used, so calculating tables, one     
748 G4double G4ChipsHyperonElasticXS::GetTabValues    
749                                                   
750 {                                                 
751   //AR-04Jun2014  if(PDG==3222 || PDG<3000 ||     
752   if(PDG<3000 || PDG>3334) G4cout<<"*Warning*G    
753                                                   
754                                                   
755   //AR-24Apr2018 Switch to allow transuranic e    
756   const G4bool isHeavyElementAllowed = true;      
757   if(tgZ<0 || ( !isHeavyElementAllowed && tgZ>    
758   {                                               
759     G4cout<<"*Warning*G4QHyperonElastCS::GetTa    
760     return 0.;                                    
761   }                                               
762   G4int iZ=tgZ-1; // Z index                      
763   if(iZ<0)                                        
764   {                                               
765     iZ=0;         // conversion of the neutron    
766     tgZ=1;                                        
767     tgN=0;                                        
768   }                                               
769   G4double p=G4Exp(lp);              // moment    
770   G4double sp=std::sqrt(p);             // sqr    
771   G4double p2=p*p;                                
772   G4double p3=p2*p;                               
773   G4double p4=p3*p;                               
774   if ( tgZ == 1 && tgN == 0 ) // Hyperon+P        
775   {                                               
776     G4double dl2=lp-lastPAR[9];                   
777     theSS=lastPAR[32];                            
778     theS1=(lastPAR[10]+lastPAR[11]*dl2*dl2)/(1    
779           (lastPAR[13]/p2+lastPAR[14]*p)/(p4+l    
780     theB1=lastPAR[16]*G4Pow::GetInstance()->po    
781     theS2=lastPAR[19]+lastPAR[20]/(p4+lastPAR[    
782     theB2=lastPAR[22]+lastPAR[23]/(p4+lastPAR[    
783     theS3=lastPAR[25]+lastPAR[26]/(p4*p4+lastP    
784     theB3=lastPAR[29]+lastPAR[30]/(p4+lastPAR[    
785     theS4=0.;                                     
786     theB4=0.;                                     
787     // Returns the total elastic pim-p cross-s    
788     G4double dp=lp-lastPAR[4];                    
789     return lastPAR[0]/(lastPAR[1]+p2*(lastPAR[    
790            lastPAR[6]/p2)/(1.+lastPAR[7]/sp+la    
791   }                                               
792   else                                            
793   {                                               
794     G4double p5=p4*p;                             
795     G4double p6=p5*p;                             
796     G4double p8=p6*p2;                            
797     G4double p10=p8*p2;                           
798     G4double p12=p10*p2;                          
799     G4double p16=p8*p8;                           
800     //G4double p24=p16*p8;                        
801     G4double dl=lp-5.;                            
802     G4double a=tgZ+tgN;                           
803     G4double pah=G4Pow::GetInstance()->powA(p,    
804     G4double pa=pah*pah;                          
805     G4double pa2=pa*pa;                           
806     if(a<6.5)                                     
807     {                                             
808       theS1=lastPAR[9]/(1.+lastPAR[10]*p4*pa)+    
809             (lastPAR[13]*dl*dl+lastPAR[14])/(1    
810       theB1=(lastPAR[16]+lastPAR[17]*p2)/(p4+l    
811       theSS=lastPAR[20]/(1.+lastPAR[21]/p2)+la    
812       theS2=lastPAR[24]/(pa/p2+lastPAR[25]/p4)    
813       theB2=lastPAR[27]*G4Pow::GetInstance()->    
814       theS3=lastPAR[31]/(pa*p+lastPAR[32]/pa)+    
815       theB3=lastPAR[34]/(p3+lastPAR[35]/p6)+la    
816       theS4=p2*(pah*lastPAR[38]*G4Exp(-pah*las    
817                 lastPAR[40]/(1.+lastPAR[41]*G4    
818       theB4=lastPAR[43]*pa/p2/(1.+pa*lastPAR[4    
819     }                                             
820     else                                          
821     {                                             
822       theS1=lastPAR[9]/(1.+lastPAR[10]/p4)+las    
823             lastPAR[13]/(p5+lastPAR[14]/p16);     
824       theB1=(lastPAR[15]/p8+lastPAR[19])/(p+la    
825             lastPAR[17]/(1.+lastPAR[18]/p4);      
826       theSS=lastPAR[21]/(p4/G4Pow::GetInstance    
827       theS2=lastPAR[24]/p4/(G4Pow::GetInstance    
828       theB2=lastPAR[28]/G4Pow::GetInstance()->    
829       theS3=lastPAR[32]/G4Pow::GetInstance()->    
830             lastPAR[33]/(1.+lastPAR[34]/p6);      
831       theB3=lastPAR[37]/p8+lastPAR[38]/p2+last    
832       theS4=(lastPAR[41]/p4+lastPAR[46]/p)/(1.    
833             (lastPAR[43]+lastPAR[44]*dl*dl)/(1    
834       theB4=lastPAR[47]/(1.+lastPAR[48]/p)+las    
835     }                                             
836     // Returns the total elastic (n/p)A cross-    
837     G4double dlp=lp-lastPAR[5]; // ax             
838     //         p1               p2          p3    
839     return (lastPAR[0]*dlp*dlp+lastPAR[1])/(1.    
840   }                                               
841   return 0.;                                      
842 } // End of GetTableValues                        
843                                                   
844 // Returns max -t=Q2 (GeV^2) for the momentum     
845 G4double G4ChipsHyperonElasticXS::GetQ2max(G4i    
846                                                   
847 {                                                 
848   static const G4double mLamb= G4Lambda::Lambd    
849   static const G4double mLa2= mLamb*mLamb;        
850   G4double pP2=pP*pP;                             
851   if(tgZ || tgN>-1)                               
852   {                                               
853     G4double mt=G4ParticleTable::GetParticleTa    
854                                                   
855     G4double dmt=mt+mt;                           
856     G4double mds=dmt*std::sqrt(pP2+mLa2)+mLa2+    
857     return dmt*dmt*pP2/mds;                       
858   }                                               
859   else                                            
860   {                                               
861     // G4cout<<"*Error*G4ChipsHyperonElasticXS    
862     //       <<tgN<<", while it is defined onl    
863     // throw G4QException("G4ChipsHyperonElast    
864     G4ExceptionDescription ed;                    
865     ed << "PDG = " << PDG << ", Z = " << tgZ <    
866        << ", while it is defined only for p pr    
867     G4Exception("G4ChipsHyperonElasticXS::GetQ    
868                 FatalException, ed);              
869     return 0;                                     
870   }                                               
871 }                                                 
872