G4QProtonElasticCrossSection Class Reference

#include <G4QProtonElasticCrossSection.hh>

Inheritance diagram for G4QProtonElasticCrossSection:

G4VQCrossSection

Public Member Functions

 ~G4QProtonElasticCrossSection ()
virtual G4double GetCrossSection (G4bool fCS, G4double pMom, G4int tgZ, G4int tgN, G4int pPDG=2212)
G4double CalculateCrossSection (G4bool CS, G4int F, G4int I, G4int pPDG, G4int Z, G4int N, G4double pP)
G4double GetSlope (G4int tZ, G4int tN, G4int pPDG)
G4double GetExchangeT (G4int tZ, G4int tN, G4int pPDG)
G4double GetHMaxT ()

Static Public Member Functions

static G4VQCrossSectionGetPointer ()

Protected Member Functions

 G4QProtonElasticCrossSection ()

Detailed Description

Definition at line 47 of file G4QProtonElasticCrossSection.hh.


Constructor & Destructor Documentation

G4QProtonElasticCrossSection::G4QProtonElasticCrossSection (  )  [protected]

Definition at line 100 of file G4QProtonElasticCrossSection.cc.

00101 {
00102 }

G4QProtonElasticCrossSection::~G4QProtonElasticCrossSection (  ) 

Definition at line 104 of file G4QProtonElasticCrossSection.cc.

00105 {
00106   std::vector<G4double*>::iterator pos;
00107   for (pos=CST.begin(); pos<CST.end(); pos++)
00108   { delete [] *pos; }
00109   CST.clear();
00110   for (pos=PAR.begin(); pos<PAR.end(); pos++)
00111   { delete [] *pos; }
00112   PAR.clear();
00113   for (pos=SST.begin(); pos<SST.end(); pos++)
00114   { delete [] *pos; }
00115   SST.clear();
00116   for (pos=S1T.begin(); pos<S1T.end(); pos++)
00117   { delete [] *pos; }
00118   S1T.clear();
00119   for (pos=B1T.begin(); pos<B1T.end(); pos++)
00120   { delete [] *pos; }
00121   B1T.clear();
00122   for (pos=S2T.begin(); pos<S2T.end(); pos++)
00123   { delete [] *pos; }
00124   S2T.clear();
00125   for (pos=B2T.begin(); pos<B2T.end(); pos++)
00126   { delete [] *pos; }
00127   B2T.clear();
00128   for (pos=S3T.begin(); pos<S3T.end(); pos++)
00129   { delete [] *pos; }
00130   S3T.clear();
00131   for (pos=B3T.begin(); pos<B3T.end(); pos++)
00132   { delete [] *pos; }
00133   B3T.clear();
00134   for (pos=S4T.begin(); pos<S4T.end(); pos++)
00135   { delete [] *pos; }
00136   S4T.clear();
00137   for (pos=B4T.begin(); pos<B4T.end(); pos++)
00138   { delete [] *pos; }
00139   B4T.clear();
00140 }


Member Function Documentation

G4double G4QProtonElasticCrossSection::CalculateCrossSection ( G4bool  CS,
G4int  F,
G4int  I,
G4int  pPDG,
G4int  Z,
G4int  N,
G4double  pP 
) [virtual]

Implements G4VQCrossSection.

Definition at line 286 of file G4QProtonElasticCrossSection.cc.

References G4cout, and G4endl.

Referenced by GetCrossSection().

00288 {
00289   // *** Begin of Associative Memory DB for acceleration of the cross section calculations
00290   static std::vector <G4double>  PIN;   // Vector of max initialized log(P) in the table
00291   // *** End of Static Definitions (Associative Memory Data Base) ***
00292   G4double pMom=pIU/GeV;                // All calculations are in GeV
00293   onlyCS=CS;                            // Flag to calculate only CS (not Si/Bi)
00294 #ifdef pdebug
00295   G4cout<<"G4QProtonElasticCrossS::CalcCS:->onlyCS="<<onlyCS<<",F="<<F<<",p="<<pIU<<G4endl;
00296 #endif
00297   lastLP=std::log(pMom);                // Make a logarithm of the momentum for calculation
00298   if(F)                                 // This isotope was found in AMDB =>RETRIEVE/UPDATE
00299   {
00300     if(F<0)                             // the AMDB must be loded
00301     {
00302       lastPIN = PIN[I];                 // Max log(P) initialised for this table set
00303       lastPAR = PAR[I];                 // Pointer to the parameter set
00304       lastCST = CST[I];                 // Pointer to the total sross-section table
00305       lastSST = SST[I];                 // Pointer to the first squared slope
00306       lastS1T = S1T[I];                 // Pointer to the first mantissa
00307       lastB1T = B1T[I];                 // Pointer to the first slope
00308       lastS2T = S2T[I];                 // Pointer to the second mantissa
00309       lastB2T = B2T[I];                 // Pointer to the second slope
00310       lastS3T = S3T[I];                 // Pointer to the third mantissa
00311       lastB3T = B3T[I];                 // Pointer to the rhird slope
00312       lastS4T = S4T[I];                 // Pointer to the 4-th mantissa
00313       lastB4T = B4T[I];                 // Pointer to the 4-th slope
00314 #ifdef pdebug
00315       G4cout<<"G4QElasticCS::CalcCS: DB is updated for I="<<I<<",*,PIN4="<<PIN[4]<<G4endl;
00316 #endif
00317     }
00318 #ifdef pdebug
00319     G4cout<<"G4QProtonElasticCrossS::CalcCS:*read*, LP="<<lastLP<<",PIN="<<lastPIN<<G4endl;
00320 #endif
00321     if(lastLP>lastPIN && lastLP<lPMax)
00322     {
00323       lastPIN=GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);// Can update upper logP-Limit in tabs
00324 #ifdef pdebug
00325       G4cout<<"G4QElCS::CalcCS:*updated(I)*,LP="<<lastLP<<"<IN["<<I<<"]="<<lastPIN<<G4endl;
00326 #endif
00327       PIN[I]=lastPIN;                   // Remember the new P-Limit of the tables
00328     }
00329   }
00330   else                                  // This isotope wasn't initialized => CREATE
00331   {
00332     lastPAR = new G4double[nPoints];    // Allocate memory for parameters of CS function
00333     lastPAR[nLast]=0;                   // Initialization for VALGRIND
00334     lastCST = new G4double[nPoints];    // Allocate memory for Tabulated CS function    
00335     lastSST = new G4double[nPoints];    // Allocate memory for Tabulated first sqaredSlope 
00336     lastS1T = new G4double[nPoints];    // Allocate memory for Tabulated first mantissa 
00337     lastB1T = new G4double[nPoints];    // Allocate memory for Tabulated first slope    
00338     lastS2T = new G4double[nPoints];    // Allocate memory for Tabulated second mantissa
00339     lastB2T = new G4double[nPoints];    // Allocate memory for Tabulated second slope   
00340     lastS3T = new G4double[nPoints];    // Allocate memory for Tabulated third mantissa 
00341     lastB3T = new G4double[nPoints];    // Allocate memory for Tabulated third slope    
00342     lastS4T = new G4double[nPoints];    // Allocate memory for Tabulated 4-th mantissa 
00343     lastB4T = new G4double[nPoints];    // Allocate memory for Tabulated 4-th slope    
00344 #ifdef pdebug
00345     G4cout<<"G4QProtonElasticCrossS::CalcCS:*ini*,lastLP="<<lastLP<<",min="<<lPMin<<G4endl;
00346 #endif
00347     lastPIN = GetPTables(lastLP,lPMin,PDG,tgZ,tgN); // Returns the new P-limit for tables
00348 #ifdef pdebug
00349     G4cout<<"G4QProtElCS::CCS:i,Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<",LP"<<lastPIN<<G4endl;
00350 #endif
00351     PIN.push_back(lastPIN);             // Fill parameters of CS function to AMDB
00352     PAR.push_back(lastPAR);             // Fill parameters of CS function to AMDB
00353     CST.push_back(lastCST);             // Fill Tabulated CS function to AMDB    
00354     SST.push_back(lastSST);             // Fill Tabulated first sq.slope to AMDB 
00355     S1T.push_back(lastS1T);             // Fill Tabulated first mantissa to AMDB 
00356     B1T.push_back(lastB1T);             // Fill Tabulated first slope to AMDB    
00357     S2T.push_back(lastS2T);             // Fill Tabulated second mantissa to AMDB 
00358     B2T.push_back(lastB2T);             // Fill Tabulated second slope to AMDB    
00359     S3T.push_back(lastS3T);             // Fill Tabulated third mantissa to AMDB 
00360     B3T.push_back(lastB3T);             // Fill Tabulated third slope to AMDB    
00361     S4T.push_back(lastS4T);             // Fill Tabulated 4-th mantissa to AMDB 
00362     B4T.push_back(lastB4T);             // Fill Tabulated 4-th slope to AMDB    
00363   } // End of creation/update of the new set of parameters and tables
00364   // =--------= NOW Update (if necessary) and Calculate the Cross Section =------------=
00365 #ifdef pdebug
00366   G4cout<<"G4QElCS::CalcCS:?update?,LP="<<lastLP<<",IN="<<lastPIN<<",ML="<<lPMax<<G4endl;
00367 #endif
00368   if(lastLP>lastPIN && lastLP<lPMax)
00369   {
00370     lastPIN = GetPTables(lastLP,lastPIN,PDG,tgZ,tgN);
00371 #ifdef pdebug
00372     G4cout<<"G4QElCS::CalcCS: *updated(O)*, LP="<<lastLP<<" < IN="<<lastPIN<<G4endl;
00373 #endif
00374   }
00375 #ifdef pdebug
00376   G4cout<<"G4QElastCS::CalcCS: lastLP="<<lastLP<<",lPM="<<lPMin<<",lPIN="<<lastPIN<<G4endl;
00377 #endif
00378   if(!onlyCS) lastTM=GetQ2max(PDG, tgZ, tgN, pMom); // Calculate (-t)_max=Q2_max (GeV2)
00379 #ifdef pdebug
00380   G4cout<<"G4QElasticCrosSec::CalcCS:oCS="<<onlyCS<<",-t="<<lastTM<<", p="<<lastLP<<G4endl;
00381 #endif
00382   if(lastLP>lPMin && lastLP<=lastPIN)   // Linear fit is made using precalculated tables
00383   {
00384     if(lastLP==lastPIN)
00385     {
00386       G4double shift=(lastLP-lPMin)/dlnP+.000001; // Log distance from lPMin
00387       G4int    blast=static_cast<int>(shift); // this is a bin number of the lower edge (0)
00388       if(blast<0 || blast>=nLast) G4cout<<"G4QEleastCS::CCS:b="<<blast<<","<<nLast<<G4endl;
00389       lastSIG = lastCST[blast];
00390       if(!onlyCS)                       // Skip the differential cross-section parameters
00391       {
00392         theSS  = lastSST[blast];
00393         theS1  = lastS1T[blast];
00394         theB1  = lastB1T[blast];
00395         theS2  = lastS2T[blast];
00396         theB2  = lastB2T[blast];
00397         theS3  = lastS3T[blast];
00398         theB3  = lastB3T[blast];
00399         theS4  = lastS4T[blast];
00400         theB4  = lastB4T[blast];
00401       }
00402 #ifdef pdebug
00403       G4cout<<"G4QProtonElasticCrossS::CalculateCS:(E) S1="<<theS1<<", B1="<<theB1<<G4endl;
00404 #endif
00405     }
00406     else
00407     {
00408       G4double shift=(lastLP-lPMin)/dlnP;        // a shift from the beginning of the table
00409       G4int    blast=static_cast<int>(shift);    // the lower bin number
00410       if(blast<0)   blast=0;
00411       if(blast>=nLast) blast=nLast-1;            // low edge of the last bin
00412       shift-=blast;                              // step inside the unit bin
00413       G4int lastL=blast+1;                       // the upper bin number
00414       G4double SIGL=lastCST[blast];              // the basic value of the cross-section
00415       lastSIG= SIGL+shift*(lastCST[lastL]-SIGL); // calculated total elastic cross-section
00416 #ifdef pdebug
00417       G4cout<<"G4QElCS::CalcCrossSection: Sig="<<lastSIG<<", P="<<pMom<<", Z="<<tgZ<<", N="
00418             <<tgN<<", PDG="<<PDG<<", onlyCS="<<onlyCS<<G4endl;
00419 #endif
00420       if(!onlyCS)                       // Skip the differential cross-section parameters
00421       {
00422         G4double SSTL=lastSST[blast];           // the low bin of the first squared slope
00423         theSS=SSTL+shift*(lastSST[lastL]-SSTL); // the basic value of the first sq.slope
00424         G4double S1TL=lastS1T[blast];           // the low bin of the first mantissa
00425         theS1=S1TL+shift*(lastS1T[lastL]-S1TL); // the basic value of the first mantissa
00426         G4double B1TL=lastB1T[blast];           // the low bin of the first slope
00427 #ifdef pdebug
00428         G4cout<<"G4QElCS::CalcCrossSection:bl="<<blast<<",ls="<<lastL<<",SL="<<S1TL<<",SU="
00429               <<lastS1T[lastL]<<",BL="<<B1TL<<",BU="<<lastB1T[lastL]<<G4endl;
00430 #endif
00431         theB1=B1TL+shift*(lastB1T[lastL]-B1TL); // the basic value of the first slope
00432         G4double S2TL=lastS2T[blast];           // the low bin of the second mantissa
00433         theS2=S2TL+shift*(lastS2T[lastL]-S2TL); // the basic value of the second mantissa
00434         G4double B2TL=lastB2T[blast];           // the low bin of the second slope
00435         theB2=B2TL+shift*(lastB2T[lastL]-B2TL); // the basic value of the second slope
00436         G4double S3TL=lastS3T[blast];           // the low bin of the third mantissa
00437         theS3=S3TL+shift*(lastS3T[lastL]-S3TL); // the basic value of the third mantissa
00438 #ifdef pdebug
00439         G4cout<<"G4QElCS::CCS: s3l="<<S3TL<<",sh3="<<shift<<",s3h="<<lastS3T[lastL]<<",b="
00440               <<blast<<",l="<<lastL<<G4endl;
00441 #endif
00442         G4double B3TL=lastB3T[blast];           // the low bin of the third slope
00443         theB3=B3TL+shift*(lastB3T[lastL]-B3TL); // the basic value of the third slope
00444         G4double S4TL=lastS4T[blast];           // the low bin of the 4-th mantissa
00445         theS4=S4TL+shift*(lastS4T[lastL]-S4TL); // the basic value of the 4-th mantissa
00446 #ifdef pdebug
00447         G4cout<<"G4QElCS::CCS: s4l="<<S4TL<<",sh4="<<shift<<",s4h="<<lastS4T[lastL]<<",b="
00448               <<blast<<",l="<<lastL<<G4endl;
00449 #endif
00450         G4double B4TL=lastB4T[blast];           // the low bin of the 4-th slope
00451         theB4=B4TL+shift*(lastB4T[lastL]-B4TL); // the basic value of the 4-th slope
00452       }
00453 #ifdef pdebug
00454       G4cout<<"G4QProtonElasticCrossS::CalculateCS:(I) S1="<<theS1<<", B1="<<theB1<<G4endl;
00455 #endif
00456     }
00457   }
00458   else lastSIG=GetTabValues(lastLP, PDG, tgZ, tgN); // Direct calculation beyond the table
00459   if(lastSIG<0.) lastSIG = 0.;                   // @@ a Warning print can be added
00460 #ifdef pdebug
00461   G4cout<<"G4QProtonElasticCrossSection::CalculateCS: END, onlyCS="<<onlyCS<<G4endl;
00462 #endif
00463   return lastSIG;
00464 }

G4double G4QProtonElasticCrossSection::GetCrossSection ( G4bool  fCS,
G4double  pMom,
G4int  tgZ,
G4int  tgN,
G4int  pPDG = 2212 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 151 of file G4QProtonElasticCrossSection.cc.

References CalculateCrossSection(), G4cout, G4endl, and G4VQCrossSection::ThresholdEnergy().

00153 {
00154   static std::vector <G4int>    colN;  // Vector of N for calculated nuclei (isotops)
00155   static std::vector <G4int>    colZ;  // Vector of Z for calculated nuclei (isotops)
00156   static std::vector <G4double> colP;  // Vector of last momenta for the reaction
00157   static std::vector <G4double> colTH; // Vector of energy thresholds for the reaction
00158   static std::vector <G4double> colCS; // Vector of last cross sections for the reaction
00159   // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---***
00160   G4double pEn=pMom;
00161   onlyCS=fCS;
00162 #ifdef pdebug
00163   G4cout<<"G4QPElCS::GetCS:>>> f="<<fCS<<", p="<<pMom<<", Z="<<tgZ<<"("<<lastZ<<") ,N="
00164         <<tgN<<"("<<lastN<<"), T="<<pEn<<"("<<lastTH<<")"<<",Sz="<<colN.size()<<G4endl;
00165   //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00166 #endif
00167   if(pPDG!=2212)
00168   {
00169     G4cout<<"*Warning*G4QProtonElCS::GetCS: *** Found pPDG="<<pPDG<<" =--=> CS=0"<<G4endl;
00170     //CalculateCrossSection(fCS,-27,j,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00171     return 0.;                         // projectile PDG=0 is a mistake (?!) @@
00172   }
00173   G4bool in=false;                   // By default the isotope must be found in the AMDB
00174   lastP   = 0.;                      // New momentum history (nothing to compare with)
00175   lastN   = tgN;                     // The last N of the calculated nucleus
00176   lastZ   = tgZ;                     // The last Z of the calculated nucleus
00177   lastI   = colN.size();             // Size of the Associative Memory DB in the heap
00178   if(lastI) for(G4int i=0; i<lastI; i++) // Loop over proj/tgZ/tgN lines of DB
00179   {                                  // The nucleus with projPDG is found in AMDB
00180     if(colN[i]==tgN && colZ[i]==tgZ) // Isotope is foind in AMDB
00181     {
00182       lastI=i;
00183       lastTH =colTH[i];              // Last THreshold (A-dependent)
00184 #ifdef pdebug
00185       G4cout<<"G4QElCS::GetCS:*Found* P="<<pMom<<",Threshold="<<lastTH<<",i="<<i<<G4endl;
00186       //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00187 #endif
00188       if(pEn<=lastTH)
00189       {
00190 #ifdef pdebug
00191         G4cout<<"G4QElCS::GetCS:Found T="<<pEn<<" < Threshold="<<lastTH<<",CS=0"<<G4endl;
00192         //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00193 #endif
00194         return 0.;                   // Energy is below the Threshold value
00195       }
00196       lastP  =colP [i];              // Last Momentum  (A-dependent)
00197       lastCS =colCS[i];              // Last CrossSect (A-dependent)
00198       //  if(std::fabs(lastP/pMom-1.)<tolerance) //VI (do not use tolerance)
00199       if(lastP == pMom)              // Do not recalculate
00200       {
00201 #ifdef pdebug
00202         G4cout<<"G4QElCS::GetCS:P="<<pMom<<",CS="<<lastCS*millibarn<<G4endl;
00203 #endif
00204         CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // Update param's only
00205         return lastCS*millibarn;     // Use theLastCS
00206       }
00207       in = true;                       // This is the case when the isotop is found in DB
00208       // Momentum pMom is in IU ! @@ Units
00209 #ifdef pdebug
00210       G4cout<<"G4QElCS::G:UpdateDB P="<<pMom<<",f="<<fCS<<",I="<<lastI<<",i="<<i<<G4endl;
00211 #endif
00212       lastCS=CalculateCrossSection(fCS,-1,i,pPDG,lastZ,lastN,pMom); // read & update
00213 #ifdef pdebug
00214       G4cout<<"G4QElCS::GetCrosSec: *****> New (inDB) Calculated CS="<<lastCS<<G4endl;
00215       //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00216 #endif
00217       if(lastCS<=0. && pEn>lastTH)    // Correct the threshold
00218       {
00219 #ifdef pdebug
00220         G4cout<<"G4QElCS::GetCS: New T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00221 #endif
00222         lastTH=pEn;
00223       }
00224       break;                           // Go out of the LOOP with found lastI
00225     }
00226 #ifdef pdebug
00227     G4cout<<"---G4QElCrossSec::GetCrosSec:pPDG="<<pPDG<<",i="<<i<<",N="<<colN[i]
00228           <<",Z["<<i<<"]="<<colZ[i]<<G4endl;
00229     //CalculateCrossSection(fCS,-27,i,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00230 #endif
00231   } // End of attampt to find the nucleus in DB
00232   if(!in)                            // This nucleus has not been calculated previously
00233   {
00234 #ifdef pdebug
00235     G4cout<<"G4QElCS::GetCrosSec:CalcNew P="<<pMom<<",f="<<fCS<<",lastI="<<lastI<<G4endl;
00236 #endif
00238     lastCS=CalculateCrossSection(fCS,0,lastI,pPDG,lastZ,lastN,pMom);//calculate&create
00239     if(lastCS<=0.)
00240     {
00241       lastTH = ThresholdEnergy(tgZ, tgN); // The Threshold Energy which is now the last
00242 #ifdef pdebug
00243       G4cout<<"G4QElCrossSection::GetCrossSect: NewThresh="<<lastTH<<",T="<<pEn<<G4endl;
00244 #endif
00245       if(pEn>lastTH)
00246       {
00247 #ifdef pdebug
00248         G4cout<<"G4QElCS::GetCS: First T="<<pEn<<"(CS=0) > Threshold="<<lastTH<<G4endl;
00249 #endif
00250         lastTH=pEn;
00251       }
00252     }
00253 #ifdef pdebug
00254     G4cout<<"G4QElCS::GetCrosSec: New CS="<<lastCS<<",lZ="<<lastN<<",lN="<<lastZ<<G4endl;
00255     //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00256 #endif
00257     colN.push_back(tgN);
00258     colZ.push_back(tgZ);
00259     colP.push_back(pMom);
00260     colTH.push_back(lastTH);
00261     colCS.push_back(lastCS);
00262 #ifdef pdebug
00263     G4cout<<"G4QElCS::GetCS:1st,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00264     //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00265 #endif
00266     return lastCS*millibarn;
00267   } // End of creation of the new set of parameters
00268   else
00269   {
00270 #ifdef pdebug
00271     G4cout<<"G4QElCS::GetCS: Update lastI="<<lastI<<G4endl;
00272 #endif
00273     colP[lastI]=pMom;
00274     colCS[lastI]=lastCS;
00275   }
00276 #ifdef pdebug
00277   G4cout<<"G4QElCS::GetCrSec:End,P="<<pMom<<"(MeV),CS="<<lastCS*millibarn<<"(mb)"<<G4endl;
00278   //CalculateCrossSection(fCS,-27,lastI,pPDG,lastZ,lastN,pMom); // DUMMY TEST
00279   G4cout<<"G4QElCS::GetCrSec:***End***, onlyCS="<<onlyCS<<G4endl;
00280 #endif
00281   return lastCS*millibarn;
00282 }

G4double G4QProtonElasticCrossSection::GetExchangeT ( G4int  tZ,
G4int  tN,
G4int  pPDG 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 744 of file G4QProtonElasticCrossSection.cc.

References G4cout, G4endl, and G4UniformRand.

00745 {
00746   static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00747   static const G4double third=1./3.;
00748   static const G4double fifth=1./5.;
00749   static const G4double sevth=1./7.;
00750 #ifdef tdebug
00751   G4cout<<"G4QProtElCS::GetExcT: F="<<onlyCS<<",Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00752 #endif
00753   if(PDG!=2212) G4cout<<"**Warning*G4QProtonElasticCrossSection::GetExT:PDG="<<PDG<<G4endl;
00754   if(onlyCS) G4cout<<"**Warning*G4QProtonElasticCrossSection::GetExchanT:onlyCS=1"<<G4endl;
00755   if(lastLP<-4.3) return lastTM*GeVSQ*G4UniformRand();// S-wave for p<14 MeV/c (kinE<.1MeV)
00756   G4double q2=0.;
00757   if(tgZ==1 && tgN==0)                // ===> p+p=p+p
00758   {
00759 #ifdef tdebug
00760     G4cout<<"G4QElasticCS::GetExchangeT: TM="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",S2="
00761           <<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",GeV2="<<GeVSQ<<G4endl;
00762 #endif
00763     G4double E1=lastTM*theB1;
00764     G4double R1=(1.-std::exp(-E1));
00765 #ifdef tdebug
00766     G4double ts1=-std::log(1.-R1)/theB1;
00767     G4double ds1=std::fabs(ts1-lastTM)/lastTM;
00768     if(ds1>.0001)
00769       G4cout<<"*Warning*G4QElCS::GetExT:1p "<<ts1<<"#"<<lastTM<<",d="<<ds1
00770             <<",R1="<<R1<<",E1="<<E1<<G4endl;
00771 #endif
00772     G4double E2=lastTM*theB2;
00773     G4double R2=(1.-std::exp(-E2*E2*E2));
00774 #ifdef tdebug
00775     G4double ts2=std::pow(-std::log(1.-R2),.333333333)/theB2;
00776     G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00777     if(ds2>.0001)
00778       G4cout<<"*Warning*G4QElCS::GetExT:2p "<<ts2<<"#"<<lastTM<<",d="<<ds2
00779             <<",R2="<<R2<<",E2="<<E2<<G4endl;
00780 #endif
00781     G4double E3=lastTM*theB3;
00782     G4double R3=(1.-std::exp(-E3));
00783 #ifdef tdebug
00784     G4double ts3=-std::log(1.-R3)/theB3;
00785     G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00786     if(ds3>.0001)
00787       G4cout<<"*Warning*G4QElCS::GetExT:3p "<<ts3<<"#"<<lastTM<<",d="<<ds3
00788             <<",R3="<<R1<<",E3="<<E3<<G4endl;
00789 #endif
00790     G4double I1=R1*theS1/theB1;
00791     G4double I2=R2*theS2;
00792     G4double I3=R3*theS3;
00793     G4double I12=I1+I2;
00794     G4double rand=(I12+I3)*G4UniformRand();
00795     if     (rand<I1 )
00796     {
00797       G4double ran=R1*G4UniformRand();
00798       if(ran>1.) ran=1.;
00799       q2=-std::log(1.-ran)/theB1;
00800     }
00801     else if(rand<I12)
00802     {
00803       G4double ran=R2*G4UniformRand();
00804       if(ran>1.) ran=1.;
00805       q2=-std::log(1.-ran);
00806       if(q2<0.) q2=0.;
00807       q2=std::pow(q2,third)/theB2;
00808     }
00809     else
00810     {
00811       G4double ran=R3*G4UniformRand();
00812       if(ran>1.) ran=1.;
00813       q2=-std::log(1.-ran)/theB3;
00814     }
00815   }
00816   else
00817   {
00818     G4double a=tgZ+tgN;
00819 #ifdef tdebug
00820     G4cout<<"G4QElCS::GetExT: a="<<a<<",t="<<lastTM<<",S1="<<theS1<<",B1="<<theB1<<",SS="
00821           <<theSS<<",S2="<<theS2<<",B2="<<theB2<<",S3="<<theS3<<",B3="<<theB3<<",S4="
00822           <<theS4<<",B4="<<theB4<<G4endl;
00823 #endif
00824     G4double E1=lastTM*(theB1+lastTM*theSS);
00825     G4double R1=(1.-std::exp(-E1));
00826     G4double tss=theSS+theSS; // for future solution of quadratic equation (imediate check)
00827 #ifdef tdebug
00828     G4double ts1=-std::log(1.-R1)/theB1;
00829     if(std::fabs(tss)>1.e-7) ts1=(std::sqrt(theB1*(theB1+(tss+tss)*ts1))-theB1)/tss;
00830     G4double ds1=(ts1-lastTM)/lastTM;
00831     if(ds1>.0001)
00832       G4cout<<"*Warning*G4QElCS::GetExT:1a "<<ts1<<"#"<<lastTM<<",d="<<ds1
00833             <<",R1="<<R1<<",E1="<<E1<<G4endl;
00834 #endif
00835     G4double tm2=lastTM*lastTM;
00836     G4double E2=lastTM*tm2*theB2;                   // power 3 for lowA, 5 for HighA (1st)
00837     if(a>6.5)E2*=tm2;                               // for heavy nuclei
00838     G4double R2=(1.-std::exp(-E2));
00839 #ifdef tdebug
00840     G4double ts2=-std::log(1.-R2)/theB2;
00841     if(a<6.5)ts2=std::pow(ts2,third);
00842     else     ts2=std::pow(ts2,fifth);
00843     G4double ds2=std::fabs(ts2-lastTM)/lastTM;
00844     if(ds2>.0001)
00845       G4cout<<"*Warning*G4QElCS::GetExT:2a "<<ts2<<"#"<<lastTM<<",d="<<ds2
00846             <<",R2="<<R2<<",E2="<<E2<<G4endl;
00847 #endif
00848     G4double E3=lastTM*theB3;
00849     if(a>6.5)E3*=tm2*tm2*tm2;                       // power 1 for lowA, 7 (2nd) for HighA
00850     G4double R3=(1.-std::exp(-E3));
00851 #ifdef tdebug
00852     G4double ts3=-std::log(1.-R3)/theB3;
00853     if(a>6.5)ts3=std::pow(ts3,sevth);
00854     G4double ds3=std::fabs(ts3-lastTM)/lastTM;
00855     if(ds3>.0001)
00856       G4cout<<"*Warning*G4QElCS::GetExT:3a "<<ts3<<"#"<<lastTM<<",d="<<ds3
00857             <<",R3="<<R3<<",E3="<<E3<<G4endl;
00858 #endif
00859     G4double E4=lastTM*theB4;
00860     G4double R4=(1.-std::exp(-E4));
00861 #ifdef tdebug
00862     G4double ts4=-std::log(1.-R4)/theB4;
00863     G4double ds4=std::fabs(ts4-lastTM)/lastTM;
00864     if(ds4>.0001)
00865       G4cout<<"*Warning*G4QElCS::GetExT:4a "<<ts4<<"#"<<lastTM<<",d="<<ds4
00866             <<",R4="<<R4<<",E4="<<E4<<G4endl;
00867 #endif
00868     G4double I1=R1*theS1;
00869     G4double I2=R2*theS2;
00870     G4double I3=R3*theS3;
00871     G4double I4=R4*theS4;
00872     G4double I12=I1+I2;
00873     G4double I13=I12+I3;
00874     G4double rand=(I13+I4)*G4UniformRand();
00875 #ifdef tdebug
00876     G4cout<<"G4QElCS::GtExT:1="<<I1<<",2="<<I2<<",3="<<I3<<",4="<<I4<<",r="<<rand<<G4endl;
00877 #endif
00878     if(rand<I1)
00879     {
00880       G4double ran=R1*G4UniformRand();
00881       if(ran>1.) ran=1.;
00882       q2=-std::log(1.-ran)/theB1;
00883       if(std::fabs(tss)>1.e-7) q2=(std::sqrt(theB1*(theB1+(tss+tss)*q2))-theB1)/tss;
00884 #ifdef tdebug
00885       G4cout<<"G4QElCS::GetExT:Q2="<<q2<<",ss="<<tss/2<<",b1="<<theB1<<",t1="<<ts1<<G4endl;
00886 #endif
00887     }
00888     else if(rand<I12)
00889     {
00890       G4double ran=R2*G4UniformRand();
00891       if(ran>1.) ran=1.;
00892       q2=-std::log(1.-ran)/theB2;
00893       if(q2<0.) q2=0.;
00894       if(a<6.5) q2=std::pow(q2,third);
00895       else      q2=std::pow(q2,fifth);
00896 #ifdef tdebug
00897       G4cout<<"G4QElCS::GetExT: Q2="<<q2<<", r2="<<R2<<", b2="<<theB2<<",t2="<<ts2<<G4endl;
00898 #endif
00899     }
00900     else if(rand<I13)
00901     {
00902       G4double ran=R3*G4UniformRand();
00903       if(ran>1.) ran=1.;
00904       q2=-std::log(1.-ran)/theB3;
00905       if(q2<0.) q2=0.;
00906       if(a>6.5) q2=std::pow(q2,sevth);
00907 #ifdef tdebug
00908       G4cout<<"G4QElCS::GetExT:Q2="<<q2<<", r3="<<R2<<", b3="<<theB2<<",t3="<<ts2<<G4endl;
00909 #endif
00910     }
00911     else
00912     {
00913       G4double ran=R4*G4UniformRand();
00914       if(ran>1.) ran=1.;
00915       q2=-std::log(1.-ran)/theB4;
00916       if(a<6.5) q2=lastTM-q2;                    // u reduced for lightA (starts from 0)
00917 #ifdef tdebug
00918       G4cout<<"G4QElCS::GetExT:Q2="<<q2<<",m="<<lastTM<<",b4="<<theB3<<",t4="<<ts3<<G4endl;
00919 #endif
00920     }
00921   }
00922   if(q2<0.) q2=0.;
00923   if(!(q2>=-1.||q2<=1.)) G4cout<<"*NAN*G4QElasticCrossSect::GetExchangeT: -t="<<q2<<G4endl;
00924   if(q2>lastTM)
00925   {
00926 #ifdef tdebug
00927     G4cout<<"*Warning*G4QElasticCrossSect::GetExT:-t="<<q2<<">"<<lastTM<<G4endl;
00928 #endif
00929     q2=lastTM;
00930   }
00931   return q2*GeVSQ;
00932 }

G4double G4QProtonElasticCrossSection::GetHMaxT (  )  [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 960 of file G4QProtonElasticCrossSection.cc.

00961 {
00962   static const G4double HGeVSQ=gigaelectronvolt*gigaelectronvolt/2.;
00963   return lastTM*HGeVSQ;
00964 }

G4VQCrossSection * G4QProtonElasticCrossSection::GetPointer (  )  [static]

Definition at line 143 of file G4QProtonElasticCrossSection.cc.

Referenced by G4QuasiFreeRatios::ChExer(), G4CHIPSElastic::G4CHIPSElastic(), G4CHIPSElasticXS::G4CHIPSElasticXS(), G4QHadronElasticDataSet::GetIsoCrossSection(), G4QIonIonElastic::GetMeanFreePath(), G4QElastic::GetMeanFreePath(), G4QLowEnergy::PostStepDoIt(), G4QIonIonElastic::PostStepDoIt(), G4QElastic::PostStepDoIt(), and G4QuasiFreeRatios::Scatter().

00144 {
00145   static G4QProtonElasticCrossSection theCrossSection;//*StaticBody of theQEl CrossSection*
00146   return &theCrossSection;
00147 }

G4double G4QProtonElasticCrossSection::GetSlope ( G4int  tZ,
G4int  tN,
G4int  pPDG 
) [virtual]

Reimplemented from G4VQCrossSection.

Definition at line 935 of file G4QProtonElasticCrossSection.cc.

References FatalException, G4cout, G4endl, and G4Exception().

00936 {
00937   static const G4double GeVSQ=gigaelectronvolt*gigaelectronvolt;
00938 #ifdef tdebug
00939   G4cout<<"G4QElasticCS::GetSlope:"<<onlyCS<<", Z="<<tgZ<<",N="<<tgN<<",PDG="<<PDG<<G4endl;
00940 #endif
00941   if(onlyCS) G4cout<<"*Warning*G4QProtonElasticCrossSection::GetSlope:onlyCS=true"<<G4endl;
00942   if(lastLP<-4.3) return 0.;          // S-wave for p<14 MeV/c (kinE<.1MeV)
00943   if(PDG!=2212)
00944   {
00945     // G4cout<<"*Error*G4QProtonElasticCrossSection::GetSlope: PDG="<<PDG<<", Z="<<tgZ<<", N="
00946     //       <<tgN<<", while it is defined only for PDG=2212"<<G4endl;
00947     // throw G4QException("G4QProtonElasticCrossSection::GetSlope: pA are implemented");
00948     G4ExceptionDescription ed;
00949     ed << "PDG = " << PDG << ", Z = " << tgZ << ", N = " << tgN
00950        << ", while it is defined only for PDG=2212 (p)" << G4endl;
00951     G4Exception("G4QProtonElasticCrossSection::GetSlope()", "HAD_CHPS_0000",
00952                 FatalException, ed);
00953   }
00954   if(theB1<0.) theB1=0.;
00955   if(!(theB1>=-1.||theB1<=1.))G4cout<<"*NAN*G4QElasticCrossSect::Getslope:"<<theB1<<G4endl;
00956   return theB1/GeVSQ;
00957 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:14 2013 for Geant4 by  doxygen 1.4.7