G4ElasticHadrNucleusHE.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //
00030 //  The generator of high energy hadron-nucleus elastic scattering
00031 //  The hadron kinetic energy T > 1 GeV
00032 //  N.  Starkov 2003.
00033 //
00034 //  19.05.04 Variant for G4 6.1: The 'ApplyYourself' was changed
00035 //  19.11.05 The HE elastic scattering on proton is added (N.Starkov)
00036 //  16.11.06 The low energy boundary is shifted to T = 400 MeV (N.Starkov)
00037 //  23.11.06 General cleanup, ONQ0=3, use pointer instead of particle name (VI)
00038 //  02.05.07 Scale sampled t as p^2 (VI)
00039 //  15.05.07 Redesign and cleanup (V.Ivanchenko)
00040 //  17.05.07 cleanup (V.Grichine)
00041 //  19.04.12 Fixed reproducibility violation (A.Ribon)
00042 //  12.06.12 Fixed warnings of shadowed variables (A.Ribon)
00043 //
00044 
00045 #include  "G4ElasticHadrNucleusHE.hh"
00046 #include  "G4PhysicalConstants.hh"
00047 #include  "G4SystemOfUnits.hh"
00048 #include  "Randomize.hh"
00049 #include  "G4ios.hh"
00050 #include  "G4ParticleTable.hh"
00051 #include  "G4IonTable.hh"
00052 #include  "G4Proton.hh"
00053 #include  "G4NistManager.hh"
00054 
00055 #include <cmath>
00056 
00057 using namespace std;
00058 
00060 //
00061 //
00062 
00063 G4ElasticData::G4ElasticData(const G4ParticleDefinition* p, 
00064                              G4int Z, G4double AWeight, G4double* eGeV)
00065 { 
00066   hadr     = p;
00067   massGeV  = p->GetPDGMass()/GeV;
00068   mass2GeV2= massGeV*massGeV;
00069   AtomicWeight = G4int(AWeight + 0.5);
00070 
00071   DefineNucleusParameters(AWeight);
00072 
00073   limitQ2 = 35./(R1*R1);     //  (GeV/c)^2
00074 
00075   G4double dQ2 = limitQ2/(ONQ2 - 1.);
00076 
00077   TableQ2[0] = 0.0;
00078 
00079   for(G4int ii = 1; ii < ONQ2; ii++) 
00080   {
00081     TableQ2[ii] = TableQ2[ii-1]+dQ2;
00082   }
00083 
00084   massA  = AWeight*amu_c2/GeV;
00085   massA2 = massA*massA; 
00086 
00087   for(G4int kk = 0; kk < NENERGY; kk++) 
00088   {
00089     dnkE[kk] = 0;
00090     G4double elab = eGeV[kk] + massGeV;
00091     G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
00092     G4double Q2m  = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
00093 
00094     if(Z == 1 && p == G4Proton::Proton()) Q2m *= 0.5;
00095 
00096     maxQ2[kk] = std::min(limitQ2, Q2m);
00097     TableCrossSec[ONQ2*kk] = 0.0;
00098 
00099 //    G4cout<<" kk  eGeV[kk] "<<kk<<"  "<<eGeV[kk]<<G4endl;
00100   }
00101 }
00102 
00104 //
00105 //
00106 
00107 void G4ElasticData::DefineNucleusParameters(G4double A)
00108 {
00109   switch (AtomicWeight)
00110     {
00111     case 207:
00112     case 208:
00113       R1       = 20.5;
00114 //      R1       = 17.5;
00115 //      R1       = 21.3;    
00116       R2       = 15.74;
00117 //      R2       = 10.74;
00118 
00119       Pnucl    = 0.4;
00120       Aeff     = 0.7;
00121       break;
00122     case 237:
00123     case 238:
00124       R1       = 21.7;    
00125       R2       = 16.5;
00126       Pnucl    = 0.4;
00127       Aeff     = 0.7;
00128       break;
00129     case 90:
00130     case 91:
00131       R1    = 16.5*1.0;
00132       R2    = 11.62;
00133       Pnucl = 0.4;
00134       Aeff  = 0.7;
00135       break;
00136     case 58:
00137     case 59:
00138       R1    = 15.0*1.05;
00139       R2    = 9.9;
00140       Pnucl = 0.45;
00141       Aeff  = 0.85;
00142       break;
00143     case 48:
00144     case 47:
00145       R1    = 14.0;
00146       R2    = 9.26;
00147       Pnucl = 0.31;
00148       Aeff  = 0.75;
00149       break;
00150     case 40:
00151     case 41:
00152       R1    = 13.3;
00153       R2    = 9.26;
00154       Pnucl = 0.31;
00155       Aeff  = 0.75;
00156       break;
00157     case 28:
00158     case 29:
00159       R1    = 12.0;
00160       R2    = 7.64;
00161       Pnucl = 0.253;
00162       Aeff  = 0.8;
00163       break;
00164     case 16:
00165       R1    = 10.50;
00166       R2    = 5.5;
00167       Pnucl = 0.7;
00168       Aeff  = 0.98;
00169       break;
00170     case 12:
00171       R1    = 9.3936;
00172       R2    = 4.63;
00173       Pnucl = 0.7;
00174 //      Pnucl = 0.5397;
00175       Aeff  = 1.0;
00176       break;
00177     case 11:
00178       R1    = 9.0;
00179       R2    = 5.42;
00180       Pnucl = 0.19;
00181 //      Pnucl = 0.5397;
00182       Aeff  = 0.9;
00183       break;
00184     case 9:
00185       R1    = 9.9;
00186       R2    = 6.5;
00187       Pnucl = 0.690;
00188       Aeff  = 0.95;
00189       break;
00190     case 4:
00191       R1    = 5.3;   
00192       R2    = 3.7;
00193       Pnucl = 0.4;
00194       Aeff  = 0.75;
00195       break;
00196     case 1:
00197       R1    = 4.5;   
00198       R2    = 2.3;
00199       Pnucl = 0.177;
00200       Aeff  = 0.9;
00201       break;
00202     default:
00203       R1    = 4.45*std::pow(A - 1.,0.309)*0.9;
00204       R2    = 2.3 *std::pow(A, 0.36);
00205 
00206       if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*A;
00207       else                 Pnucl = 0.4;
00208 
00209 //G4cout<<" Deault: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
00210 //      <<Aeff<<"  "<<Pnucl<<G4endl;
00211 
00212       if(A >= 100)               Aeff = 0.7;
00213       else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*A;
00214       else                       Aeff = 0.9;
00215       break;
00216     }
00217 //G4cout<<" Result: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
00218 //      <<Aeff<<"  "<<Pnucl<<G4endl;
00219 }
00220 
00222 //
00223 //  The constructor for the generating of events
00224 //
00225 
00226 G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE(const G4String& name)
00227   : G4HadronElastic(name)
00228 {
00229   dQ2 = hMass = hMass2 =  hLabMomentum = hLabMomentum2 = MomentumCM = HadrEnergy 
00230     = R1 = R2 = Pnucl = Aeff = HadrTot = HadrSlope = HadrReIm = TotP = DDSect2
00231     = DDSect3 = ConstU = FmaxT = Slope1 = Slope2 = Coeff1 = Coeff2 = MaxTR 
00232     = Slope0 = Coeff0 = aAIm = aDIm = Dtot11 = 0.0;
00233   NumbN = iHadrCode = iHadron = 0;
00234 
00235   verboseLevel = 0;
00236   plabLowLimit = 20.0*MeV;
00237   lowestEnergyLimit = 0.0;
00238   //Description();
00239 
00240   MbToGeV2  =  2.568;
00241   sqMbToGeV =  1.602;
00242   Fm2ToGeV2 =  25.68;
00243   GeV2      =  GeV*GeV;
00244   protonM   =  proton_mass_c2/GeV;
00245   protonM2  =  protonM*protonM;
00246 
00247   BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
00248   BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
00249   BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
00250   BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
00251   BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
00252   BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
00253   BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
00254 
00255   Binom();
00256   // energy in GeV
00257   Energy[0] = 0.4;
00258   Energy[1] = 0.6;
00259   Energy[2] = 0.8;
00260   LowEdgeEnergy[0] = 0.0;
00261   LowEdgeEnergy[1] = 0.5;
00262   LowEdgeEnergy[2] = 0.7;
00263   G4double e = 1.0;
00264   G4double f = std::pow(10.,0.1);
00265   for(G4int i=3; i<NENERGY; i++) {
00266     Energy[i] = e;
00267     LowEdgeEnergy[i] = e/f;
00268     e *= f*f;
00269   }
00270   nistManager = G4NistManager::Instance();
00271 
00272   // PDG code for hadrons
00273   G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
00274                         3122,3222,3112,3212,3312,3322,3334,
00275                         -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
00276   // internal index 
00277   G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
00278                         0,0,0,0,0,0,0,
00279                         1,7,1,1,1,1,1,1,1};
00280 
00281   G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
00282                         0,0,0,0,0,0,0,
00283                         2,2,2,2,2,2,2,2,2};
00284 
00285   for(G4int j=0; j<NHADRONS; j++) 
00286   {
00287     HadronCode[j]  = ic[j];
00288     HadronType[j]  = id[j];
00289     HadronType1[j] = id1[j];
00290 
00291     for(G4int k = 0; k < 93; k++) { SetOfElasticData[j][k] = 0; }
00292   } 
00293 }
00294 
00295 
00296 void G4ElasticHadrNucleusHE::Description() const
00297 {
00298   char* dirName = getenv("G4PhysListDocDir");
00299   if (dirName) {
00300     std::ofstream outFile;
00301     G4String outFileName = GetModelName() + ".html";
00302     G4String pathName = G4String(dirName) + "/" + outFileName;
00303     outFile.open(pathName);
00304     outFile << "<html>\n";
00305     outFile << "<head>\n";
00306 
00307     outFile << "<title>Description of G4ElasticHadrNucleusHE Model</title>\n";
00308     outFile << "</head>\n";
00309     outFile << "<body>\n";
00310 
00311     outFile << "G4ElasticHadrNucleusHE is a hadron-nucleus elastic scattering\n"
00312             << "model developed by N. Starkov which uses a Glauber model\n"
00313             << "parameterization to calculate the final state.  It is valid\n"
00314             << "for all hadrons with incident energies above 1 GeV.\n";
00315 
00316     outFile << "</body>\n";
00317     outFile << "</html>\n";
00318     outFile.close();
00319   }
00320 }
00321 
00322 
00324 //
00325 //
00326 
00327 G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE()
00328 {
00329   for(G4int j = 0; j < NHADRONS; j++) 
00330   {
00331     for(G4int k = 0; k < 93; k++) 
00332     {
00333       if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k];
00334     }
00335   } 
00336 }
00337 
00339 //
00340 //
00341 
00342 G4double 
00343 G4ElasticHadrNucleusHE::SampleInvariantT(const G4ParticleDefinition* p,
00344                                          G4double inLabMom, 
00345                                          G4int Z, G4int N)
00346 {
00347   G4double plab  = inLabMom/GeV;   // (GeV/c)
00348   G4double Q2 = 0;
00349 
00350   iHadrCode = p->GetPDGEncoding();
00351 
00352   NumbN = N;
00353 
00354   if(verboseLevel > 1)
00355   {
00356     G4cout<< " G4ElasticHadrNucleusHE::SampleT: " 
00357           << " for " << p->GetParticleName() 
00358           << " at Z= " << Z << " A= " << N
00359           << " plab(GeV)= " << plab
00360           << G4endl;
00361   }
00362   G4int idx;
00363 
00364   for( idx = 0 ; idx < NHADRONS; idx++) 
00365   {
00366     if(iHadrCode == HadronCode[idx]) break;
00367   }
00368 
00369   // Hadron is not in the list
00370 
00371   if( idx >= NHADRONS ) return Q2;
00372 
00373   iHadron = HadronType[idx];
00374   iHadrCode = HadronCode[idx];
00375 
00376   if(Z==1)
00377     {
00378       hMass  = p->GetPDGMass()/GeV;
00379       hMass2 = hMass*hMass;
00380 
00381       G4double T = sqrt(plab*plab+hMass2)-hMass;
00382 
00383       if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
00384 
00385       if (verboseLevel>1)
00386         G4cout<<"  Proton : Q2  "<<Q2<<G4endl;
00387     }
00388   else
00389     {
00390       G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
00391 
00392       // Construct elastic data
00393       if(!ElD1) 
00394         {
00395           G4double AWeight = nistManager->GetAtomicMassAmu(Z);
00396           ElD1 = new  G4ElasticData(p, Z, AWeight, Energy);
00397           SetOfElasticData[idx][Z] = ElD1;
00398     
00399           if(verboseLevel > 1)
00400             {
00401               G4cout<< " G4ElasticHadrNucleusHE::SampleT:  new record " << idx
00402                     << " for " << p->GetParticleName() << " Z= " << Z
00403                     << G4endl;
00404             }
00405         }  
00406       hMass          = ElD1->massGeV;
00407       hMass2         = ElD1->mass2GeV2;
00408       G4double M     = ElD1->massA;
00409       G4double M2    = ElD1->massA2;
00410       G4double plab2 = plab*plab;
00411       G4double Q2max = 4.*plab2*M2/
00412         (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
00413 
00414       // sample scattering
00415       G4double T = sqrt(plab2+hMass2)-hMass;
00416 
00417       if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
00418 
00419       if(verboseLevel > 1)
00420         G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< "  t/tmax= " << Q2/Q2max <<G4endl;
00421     }
00422   return  Q2*GeV2;
00423 }
00424 
00426 //
00427 //
00428 
00429 G4double 
00430 G4ElasticHadrNucleusHE::SampleT(const G4ParticleDefinition* p,
00431                                 G4double inLabMom, 
00432                                 G4int Z, G4int N)
00433 {
00434   return SampleInvariantT(p, inLabMom, Z, N);
00435 }
00436 
00438 //
00439 //
00440 
00441 G4double G4ElasticHadrNucleusHE::
00442                           HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z, 
00443                                             G4double plab, G4double tmax)
00444 {
00445   G4double LineFq2[ONQ2];
00446 
00447   G4double Rand = G4UniformRand();
00448 
00449   G4int      iNumbQ2 = 0;
00450   G4double   Q2 = 0.0;
00451 
00452   G4double ptot2 = plab*plab;
00453   G4double ekin  = std::sqrt(hMass2 + ptot2) - hMass;
00454 
00455   if(verboseLevel > 1)
00456     G4cout<<"Q2_2: ekin  plab  "<<ekin<<"    "<<plab<<"  tmax "<<tmax<<G4endl;
00457 
00458   // Find closest energy bin
00459   G4int NumbOnE; 
00460   for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ ) 
00461   {
00462     if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
00463   }
00464   G4double* dNumbQ2 = pElD->TableQ2;
00465 
00466   G4int index = NumbOnE*ONQ2;
00467 
00468   // Select kinematics for node energy
00469   G4double T     = Energy[NumbOnE];
00470   hLabMomentum2  = T*(T + 2.*hMass);
00471   G4double Q2max = pElD->maxQ2[NumbOnE];
00472   G4int length   = pElD->dnkE[NumbOnE];
00473 
00474   // Build vector
00475   if(length == 0) 
00476     {
00477       R1    = pElD->R1;
00478       R2    = pElD->R2;
00479       Aeff  = pElD->Aeff;
00480       Pnucl = pElD->Pnucl;
00481       hLabMomentum = std::sqrt(hLabMomentum2);
00482  
00483       DefineHadronValues(Z);
00484 
00485       if(verboseLevel >0)
00486         {
00487           G4cout<<"1  plab  T  "<<plab<<"  "<<T<<"  sigTot  B  ReIm  "
00488                 <<HadrTot<<"  "<<HadrSlope<<"  "<<HadrReIm<<G4endl;
00489           G4cout<<"  R1  R2  Aeff  p  "<<R1<<"  "<<R2<<"  "<<Aeff<<"  "
00490                 <<Pnucl<<G4endl;
00491         }
00492 
00493       //pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
00494 
00495       if(verboseLevel > 1)
00496         G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE 
00497               << " length= " << length 
00498               << " Q2max= " << Q2max 
00499               << " ekin= " << ekin <<G4endl;
00500     
00501       pElD->TableCrossSec[index] = 0;
00502 
00503 
00504       dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
00505 
00506       GetHeavyFq2(Z, NumbN, LineFq2);  //  %%%%%%%%%%%%%%%%%%%%%%%%%
00507 
00508       for(G4int ii=0; ii<ONQ2; ii++)
00509         {
00510           //if(verboseLevel > 2)
00511           //  G4cout<<"  ii LineFq2  "<<ii<<"  "<<LineFq2[ii]/LineFq2[ONQ2-1]
00512           //    <<"  dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
00513 
00514           pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
00515         }
00516     
00517       pElD->dnkE[NumbOnE] = ONQ2;
00518       length = ONQ2;
00519     } 
00520 
00521   G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
00522 
00523   for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ ) 
00524     {
00525       if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
00526     }
00527   Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
00528 
00529   if(tmax < Q2max) Q2 *= tmax/Q2max;
00530 
00531   if(verboseLevel > 1)
00532     G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2 
00533           << " rand= " << Rand << G4endl;
00534   
00535   return Q2;
00536 }       
00537 
00539 //
00540 //  The randomization of one dimensional array 
00541 //
00542 G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q,
00543                                          G4double * F, G4double ranUni)
00544 {
00545   G4double ranQ2;
00546 
00547   G4double F1  = F[kk-2];
00548   G4double F2  = F[kk-1];
00549   G4double F3  = F[kk];
00550   G4double X1  = Q[kk-2];
00551   G4double X2  = Q[kk-1];
00552   G4double X3  = Q[kk];
00553 
00554   if(verboseLevel > 2) 
00555     G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3 
00556            << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
00557 
00558   if(kk == 1 || kk == 0)
00559   {
00560      F1  = F[0]; 
00561      F2  = F[1];
00562      F3  = F[2];
00563      X1  = Q[0];
00564      X2  = Q[1];
00565      X3  = Q[2];
00566   }
00567 
00568   G4double F12 = F1*F1;
00569   G4double F22 = F2*F2;
00570   G4double F32 = F3*F3;
00571 
00572   G4double D0  = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
00573 
00574   if(verboseLevel > 2) 
00575     G4cout << "       X1= " << X1 << " F1= " << F1 << "  D0= " 
00576            << D0 << G4endl; 
00577 
00578   if(std::abs(D0) < 0.00000001)
00579     { 
00580       ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
00581     }
00582   else    
00583     {
00584       G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
00585       G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
00586       G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
00587                    -X1*F2*F32-X2*F3*F12-X3*F1*F22;
00588       ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
00589     }
00590   return ranQ2;         //  MeV^2
00591 }
00592 
00594 //
00595 //
00596 G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Z, G4int Nucleus, G4double* LineF) 
00597 {
00598   G4int ii, jj, aSimp;
00599   G4double curQ2, curSec;
00600   G4double curSum = 0.0;
00601   G4double totSum = 0.0;
00602 
00603   G4double ddQ2 = dQ2/20;
00604   G4double Q2l  = 0;
00605 
00606   LineF[0] = 0;
00607   for(ii = 1; ii<ONQ2; ii++)
00608     {
00609       curSum = 0;
00610       aSimp  = 4;   
00611 
00612       for(jj = 0; jj<20; jj++)
00613         {
00614           curQ2 = Q2l+jj*ddQ2;
00615 
00616           curSec  = HadrNucDifferCrSec(Z, Nucleus, curQ2);
00617           curSum += curSec*aSimp;
00618 
00619           if(aSimp > 3) aSimp = 2;
00620           else          aSimp = 4;
00621 
00622           if(jj == 0 && verboseLevel>2)
00623             G4cout<<"  Q2  "<<curQ2<<"  AIm  "<<aAIm<<"  DIm  "<<aDIm
00624                   <<"  Diff  "<<curSec<<"  totSum  "<<totSum<<G4endl;
00625         }
00626 
00627       Q2l    += dQ2;
00628       curSum *= ddQ2/2.3;   //  $$$$$$$$$$$$$$$$$$$$$$$
00629       totSum += curSum;
00630 
00631       LineF[ii] = totSum;
00632         
00633       if (verboseLevel>2)
00634         G4cout<<"  GetHeavy: Q2  dQ2  totSum  "<<Q2l<<"  "<<dQ2<<"  "<<totSum
00635               <<"  curSec  "
00636               <<curSec<<"  totSum  "<< totSum<<"  DTot "
00637               <<curSum<<G4endl;
00638     }      
00639   return totSum;
00640 }
00641 
00643 //
00644 //
00645 
00646 G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus, 
00647                                              G4double Q2)
00648 {
00649   // Scattering of proton
00650   if(Z == 1) 
00651   {
00652     G4double SqrQ2  = std::sqrt(Q2);
00653     G4double valueConstU = 2.*(hMass2 + protonM2) - Q2;
00654 
00655     G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2))
00656       + Coeff0*(1.-std::exp(-Slope0*Q2))
00657       + Coeff2/Slope2*std::exp(Slope2*valueConstU)*(std::exp(Slope2*Q2)-1.)
00658       + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
00659 
00660     return y;
00661   }
00662 
00663   // The preparing of probability function  
00664 
00665   G4double prec = Nucleus > 208  ?  1.0e-7 : 1.0e-6;
00666 
00667   G4double    Stot     = HadrTot*MbToGeV2;     //  Gev^-2
00668   G4double    Bhad     = HadrSlope;         //  GeV^-2
00669   G4double    Asq      = 1+HadrReIm*HadrReIm;
00670   G4double    Rho2     = std::sqrt(Asq);
00671 
00672 //  if(verboseLevel >1)
00673     G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<"  "<<HadrSlope<<"  "
00674       <<HadrReIm<<G4endl;
00675 
00676   if(verboseLevel > 1) {
00677     G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad 
00678            <<"  Im "<<HadrReIm 
00679            << " Asq= " << Asq << G4endl;
00680     G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
00681   }
00682   G4double    R12      = R1*R1;
00683   G4double    R22      = R2*R2;
00684   G4double    R12B     = R12+2*Bhad;
00685   G4double    R22B     = R22+2*Bhad;
00686 
00687   G4double    Norm     = (R12*R1-Pnucl*R22*R2); // HP->Aeff;
00688 
00689   G4double    R13      = R12*R1/R12B;
00690   G4double    R23      = Pnucl*R22*R2/R22B;
00691   G4double    Unucl    = Stot/twopi/Norm*R13;
00692   G4double    UnucRho2 = -Unucl*Rho2;
00693 
00694   G4double    FiH      = std::asin(HadrReIm/Rho2);
00695   G4double    NN2      = R23/R13;
00696 
00697   if(verboseLevel > 2) 
00698     G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2 
00699            << " Norm= " << Norm << G4endl;
00700 
00701   G4double    dddd;
00702  
00703   G4double    Prod0    = 0;
00704   G4double    N1       = -1.0;
00705   //G4double    Tot0     = 0;
00706   G4double    exp1;
00707 
00708   G4double    Prod3 ;
00709   G4double    exp2  ;
00710   G4double    N4, N5, N2, Prod1, Prod2;
00711   G4int    i1, i2, j1, j2;
00712 
00713   for(i1 = 1; i1<= Nucleus; i1++) 
00714     {
00715       N1    = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
00716       Prod1 = 0;
00717       //Tot0  = 0;
00718       N2    = -1;
00719 
00720       for(i2 = 1; i2<=Nucleus; i2++) 
00721         {
00722           N2    = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
00723           Prod2 = 0; 
00724           N5    = -1/NN2;
00725           for(j2=0; j2<= i2; j2++) 
00726             {
00727               Prod3 = 0;
00728               exp2  = 1/(j2/R22B+(i2-j2)/R12B);
00729               N5    = -N5*NN2;
00730               N4    = -1/NN2;
00731               for(j1=0; j1<=i1; j1++) 
00732                 {
00733                   exp1  = 1/(j1/R22B+(i1-j1)/R12B);
00734                   dddd  = exp1+exp2;
00735                   N4    = -N4*NN2;
00736                   Prod3 = Prod3+N4*exp1*exp2*
00737                     (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][j1];
00738                }                                   // j1
00739               Prod2 = Prod2 +Prod3*N5*SetBinom[i2][j2];
00740             }                                      // j2
00741           Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
00742 
00743           if (std::fabs(Prod2*N2/Prod1)<prec) break;
00744         }                                         // i2
00745       Prod0   = Prod0 + Prod1*N1;
00746       if(std::fabs(N1*Prod1/Prod0) < prec) break;
00747 
00748     }                                           // i1
00749 
00750   Prod0 *= 0.25*pi/MbToGeV2;  //  This is in mb
00751 
00752   if(verboseLevel>1) 
00753     G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus 
00754            <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
00755   return Prod0;
00756 }
00757 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00758 G4double G4ElasticHadrNucleusHE::
00759    HadrNucDifferCrSec(G4int Z, G4int , G4double aQ2)
00760 {
00761 //   ------ All external kinematical variables are in MeV -------
00762 //            ------ but internal in GeV !!!!  ------
00763 
00764   G4int NWeight = int( nistManager->GetAtomicMassAmu(Z) + 0.5 ); 
00765 
00766   G4double    theQ2 = aQ2;   
00767 
00768   // Scattering of proton
00769   if(NWeight == 1) 
00770   {
00771     G4double SqrQ2  = std::sqrt(aQ2);
00772     G4double valueConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
00773 
00774     G4double MaxT = 4*MomentumCM*MomentumCM;
00775 
00776      BoundaryTL[0] = MaxT;
00777      BoundaryTL[1] = MaxT;
00778      BoundaryTL[3] = MaxT;
00779      BoundaryTL[4] = MaxT;
00780      BoundaryTL[5] = MaxT;
00781 
00782     G4double dSigPodT;
00783 
00784     dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
00785                  (
00786                   Coeff1*std::exp(-Slope1*SqrQ2)+
00787                   Coeff2*std::exp( Slope2*(valueConstU)+aQ2)+
00788                   (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+
00789                  +Coeff0*std::exp(-Slope0*aQ2)
00790 //                +0.1*(1-std::fabs(CosTh))
00791                   )/16/3.1416*2.568;
00792 
00793     return dSigPodT;
00794   }
00795 
00796     G4double    Stot     = HadrTot*MbToGeV2; 
00797     G4double    Bhad     = HadrSlope; 
00798     G4double    Asq      = 1+HadrReIm*HadrReIm;
00799     G4double    Rho2     = std::sqrt(Asq);
00800     G4double    Pnuclp   = 0.001;
00801                 Pnuclp   = Pnucl;
00802     G4double    R12      = R1*R1;
00803     G4double    R22      = R2*R2;
00804     G4double    R12B     = R12+2*Bhad;
00805     G4double    R22B     = R22+2*Bhad;
00806     G4double    R12Ap    = R12+20;
00807     G4double    R22Ap    = R22+20;
00808     G4double    R13Ap    = R12*R1/R12Ap;
00809     G4double    R23Ap    = R22*R2/R22Ap*Pnuclp;
00810     G4double    R23dR13  = R23Ap/R13Ap;
00811     G4double    R12Apd   = 2/R12Ap;
00812     G4double    R22Apd   = 2/R22Ap;
00813     G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
00814 
00815     G4double DDSec1p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4));
00816     G4double DDSec2p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/
00817                              std::sqrt((R12+R22)/2)/4));
00818     G4double DDSec3p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4));
00819 
00820     G4double    Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;
00821     G4double    Normp    = (R12*R1-Pnuclp*R22*R2)*Aeff;
00822     G4double    R13      = R12*R1/R12B;
00823     G4double    R23      = Pnucl*R22*R2/R22B;
00824     G4double    Unucl    = Stot/twopi/Norm*R13;
00825     G4double    UnuclScr = Stot/twopi/Normp*R13Ap;
00826     G4double    SinFi    = HadrReIm/Rho2;
00827     G4double    FiH      = std::asin(SinFi);
00828     G4double    N        = -1;
00829     G4double    N2       = R23/R13;
00830 
00831     G4double    ImElasticAmpl0  = 0;
00832     G4double    ReElasticAmpl0  = 0;
00833 
00834     G4double    exp1;
00835     G4double    N4;
00836     G4double    Prod1, Tot1=0, medTot, DTot1, DmedTot;
00837     G4int       i;
00838 
00839     for( i=1; i<=NWeight; i++)
00840     {
00841       N       = -N*Unucl*(NWeight-i+1)/i*Rho2;
00842       N4      = 1;
00843       Prod1   = std::exp(-theQ2/i*R12B/4)/i*R12B;
00844       medTot  = R12B/i;
00845 
00846        for(G4int l=1; l<=i; l++)
00847        {
00848          exp1    = l/R22B+(i-l)/R12B;
00849          N4      = -N4*(i-l+1)/l*N2;
00850          Prod1   = Prod1+N4/exp1*std::exp(-theQ2/exp1/4);
00851          medTot  = medTot+N4/exp1;
00852        }  // end l
00853 
00854       ReElasticAmpl0  = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
00855       ImElasticAmpl0  = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
00856       Tot1            = Tot1+medTot*N*std::cos(FiH*i);
00857       if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break;
00858     }      // i
00859 
00860     ImElasticAmpl0 = ImElasticAmpl0*pi/2.568;   // The amplitude in mB
00861     ReElasticAmpl0 = ReElasticAmpl0*pi/2.568;   // The amplitude in mB
00862     Tot1           = Tot1*twopi/2.568;
00863 
00864     G4double C1 = R13Ap*R13Ap/2*DDSec1p;
00865     G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
00866     G4double C3 = R23Ap*R23Ap/2*DDSec3p;
00867 
00868     G4double N1p  = 1;
00869 
00870     G4double Din1 = 0.5;     
00871 
00872     Din1  = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap-
00873                  C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+
00874                  C3*R22Ap/2*std::exp(-theQ2/8*R22Ap));
00875 
00876     DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
00877 
00878     G4double exp1p;
00879     G4double exp2p;
00880     G4double exp3p;
00881     G4double N2p;
00882     G4double Din2, BinCoeff;
00883 
00884     BinCoeff = 1;
00885 
00886     for( i = 1; i<= NWeight-2; i++)
00887     {
00888       N1p     = -N1p*UnuclScr*(NWeight-i-1)/i*Rho2;
00889       N2p     = 1;
00890       Din2    = 0;
00891       DmedTot = 0;
00892         for(G4int l = 0; l<=i; l++) 
00893         {
00894           if(l == 0)      BinCoeff = 1;
00895           else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
00896 
00897           exp1  = l/R22B+(i-l)/R12B;
00898           exp1p = exp1+R12Apd;
00899           exp2p = exp1+R12ApdR22Ap;
00900           exp3p = exp1+R22Apd;
00901 
00902           Din2  = Din2 + N2p*BinCoeff*
00903                   (C1/exp1p*std::exp(-theQ2/4/exp1p)-
00904                    C2/exp2p*std::exp(-theQ2/4/exp2p)+
00905                    C3/exp3p*std::exp(-theQ2/4/exp3p));
00906 
00907           DmedTot = DmedTot + N2p*BinCoeff*
00908             (C1/exp1p-C2/exp2p+C3/exp3p);
00909 
00910           N2p   = -N2p*R23dR13;
00911         }     // l
00912 
00913         Din1  = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
00914         DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
00915  
00916         if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
00917     }           //  i
00918 
00919     Din1 = -Din1*NWeight*(NWeight-1)
00920                  /2/pi/Normp/2/pi/Normp*16*pi*pi;
00921 
00922     DTot1 = DTot1*NWeight*(NWeight-1)
00923                  /2/pi/Normp/2/pi/Normp*16*pi*pi;
00924 
00925     DTot1 *= 5;   //  $$$$$$$$$$$$$$$$$$$$$$$$ 
00926 //     Din1  *= 0.2;    //   %%%%%%%%%%%%%%%%%%%%%%%   proton
00927 //     Din1 *= 0.05;    //   %%%%%%%%%%%%%%%%%%%%%%%  pi+
00928 //  ----------------  dSigma/d|-t|,  mb/(GeV/c)^-2  -----------------
00929 
00930     G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
00931                            (ImElasticAmpl0+Din1)*
00932                            (ImElasticAmpl0+Din1))*2/4/pi;
00933 
00934     Tot1   = Tot1-DTot1;
00935      //  Tott1  = Tot1*1.0;
00936     Dtot11 = DTot1;
00937     aAIm   = ImElasticAmpl0;
00938     aDIm   = Din1;
00939 
00940     return DiffCrSec2*1.0;  //  dSig/d|-t|,  mb/(GeV/c)^-2
00941 }   // function
00942 //  ##############################################
00943 
00945 //
00946 //
00947 
00948 void  G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
00949 {
00950   HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
00951 
00952   G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
00953   G4double sqrS  = std::sqrt(sHadr);
00954   G4double Ecm   = 0.5*(sHadr-hMass2+protonM2)/sqrS;
00955   MomentumCM     = std::sqrt(Ecm*Ecm-protonM2);
00956   
00957   if(verboseLevel>2)
00958     G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron 
00959            << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
00960            << " plab= " << hLabMomentum   
00961            <<"  E - m  "<<HadrEnergy - hMass<< G4endl;
00962 
00963   G4double TotN = 0.0;
00964   G4double logE = std::log(HadrEnergy);
00965   G4double logS = std::log(sHadr);
00966            TotP = 0.0;
00967 
00968   switch (iHadron)
00969     {
00970     case 0:                  //  proton, neutron
00971     case 6:
00972 
00973       if(hLabMomentum > 10)
00974         TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); //  mb
00975 
00976       else
00977         {
00978 // ==================  neutron  ================
00979 
00981 
00982 
00983           if( hLabMomentum > 1.4 )
00984             TotN = 33.3+15.2*(hLabMomentum2-1.35)/
00985               (std::pow(hLabMomentum,2.37)+0.95);
00986                 
00987           else if(hLabMomentum > 0.8)
00988             {
00989               G4double A0 = logE + 0.0513;
00990               TotN = 33.0 + 25.5*A0*A0;  
00991             }
00992           else 
00993             {
00994               G4double A0 = logE - 0.2634;  // log(1.3)
00995               TotN = 33.0 + 30.*A0*A0*A0*A0;
00996             }
00997 //  =================  proton  ===============
00998 //       else if(iHadrCode == 2212) 
00999           {
01000             if(hLabMomentum >= 1.05)
01001               {
01002                 TotP = 39.0+75.*(hLabMomentum-1.2)/
01003                   (hLabMomentum2*hLabMomentum+0.15);
01004               }
01005 
01006             else if(hLabMomentum >= 0.7)
01007               {
01008                  G4double A0 = logE + 0.3147;
01009                  TotP = 23.0 + 40.*A0*A0;
01010               }
01011             else 
01012               {
01013                 TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
01014               }
01015           }
01016         }
01017 
01018 //      HadrTot = 0.5*(82*TotP+126*TotN)/104;  //  $$$$$$$$$$$$$$$$$$
01019       HadrTot = 0.5*(TotP+TotN);
01020 //  ...................................................
01021       //  Proton slope
01022       if(hLabMomentum >= 2.)       HadrSlope = 5.44 + 0.88*logS;
01023 
01024       else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
01025 
01026       else                         HadrSlope = 1.5;
01027 
01028 //  ...................................................
01029       if(hLabMomentum >= 1.2)
01030          HadrReIm  = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
01031        
01032       else if(hLabMomentum >= 0.6) 
01033         HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
01034           (std::pow(3*hLabMomentum,2.2)+1);     
01035 
01036       else 
01037         HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
01038 //  ...................................................
01039       DDSect2   = 2.2;                              //mb*GeV-2
01040       DDSect3   = 0.6;                               //mb*GeV-2
01041       //  ================== lambda  ==================
01042       if( iHadrCode == 3122)
01043         {
01044           HadrTot   *= 0.88;
01045           HadrSlope *=0.85;
01046         }
01047       //  ================== sigma +  ==================
01048       else if( iHadrCode == 3222)
01049         {
01050           HadrTot   *=0.81;
01051           HadrSlope *=0.85;
01052         }
01053       //  ================== sigma 0,-  ==================
01054       else if(iHadrCode == 3112 || iHadrCode == 3212 )
01055         {
01056           HadrTot   *=0.88;
01057           HadrSlope *=0.85;
01058         }
01059       //  ===================  xi  =================
01060       else if( iHadrCode == 3312 || iHadrCode == 3322 )
01061         {
01062           HadrTot   *=0.77;
01063           HadrSlope *=0.75;
01064         }
01065       //  =================  omega  =================
01066       else if( iHadrCode == 3334)
01067         {
01068           HadrTot   *=0.78;
01069           HadrSlope *=0.7;
01070         }
01071 
01072       break;
01073 //  ===========================================================
01074     case 1:              //   antiproton
01075     case 7:              //   antineutron
01076 
01077       HadrTot   = 5.2+5.2*logE + 123.2/sqrS;     //  mb
01078       HadrSlope = 8.32+0.57*logS;                //(GeV/c)^-2
01079 
01080       if( HadrEnergy < 1000 )
01081         HadrReIm  = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
01082       else
01083         HadrReIm  = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
01084 
01085       DDSect2   = 11;                            //mb*(GeV/c)^-2
01086       DDSect3   = 3;                             //mb*(GeV/c)^-2
01087       //  ================== lambda  ==================
01088       if( iHadrCode == -3122)
01089         {
01090           HadrTot   *= 0.88;
01091           HadrSlope *=0.85;
01092         }
01093       //  ================== sigma +  ==================
01094       else if( iHadrCode == -3222)
01095         {
01096           HadrTot   *=0.81;
01097           HadrSlope *=0.85;
01098         }
01099       //  ================== sigma 0,-  ==================
01100       else if(iHadrCode == -3112 || iHadrCode == -3212 )
01101         {
01102           HadrTot   *=0.88;
01103           HadrSlope *=0.85;
01104         }
01105       //  ===================  xi  =================
01106       else if( iHadrCode == -3312 || iHadrCode == -3322 )
01107         {
01108           HadrTot   *=0.77;
01109           HadrSlope *=0.75;
01110         }
01111       //  =================  omega  =================
01112       else if( iHadrCode == -3334)
01113         {
01114           HadrTot   *=0.78;
01115           HadrSlope *=0.7;
01116         }
01117 
01118       break;
01119 //  -------------------------------------------
01120     case 2:             //   pi plus, pi minus
01121     case 3:
01122 
01123       if(hLabMomentum >= 3.5)
01124         TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb
01125 //  =========================================
01126       else if(hLabMomentum >= 1.15)
01127         {
01128           G4double x = (hLabMomentum - 2.55)/0.55; 
01129           G4double y = (hLabMomentum - 1.47)/0.225;
01130           TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
01131         }
01132 //  =========================================
01133       else if(hLabMomentum >= 0.4)
01134         {
01135         TotP  = 88*(logE+0.2877)*(logE+0.2877)+14.0;
01136         }
01137 //  =========================================
01138       else 
01139         {
01140           G4double x = (hLabMomentum - 0.29)/0.085;
01141           TotP = 20. + 180.*std::exp(-x*x);
01142         }
01143 //  -------------------------------------------
01144 
01145       if(hLabMomentum >= 3.0 )
01146         TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb
01147 
01148       else if(hLabMomentum >= 1.3) 
01149         {
01150           G4double x = (hLabMomentum - 2.1)/0.4;
01151           G4double y = (hLabMomentum - 1.4)/0.12;
01152           TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) + 
01153                                               1.5*std::exp(-y*y);
01154         }
01155       else if(hLabMomentum >= 0.65)
01156         {
01157           G4double x = (hLabMomentum - 0.72)/0.06;
01158           G4double y = (hLabMomentum - 1.015)/0.075;
01159           TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
01160         }
01161       else if(hLabMomentum >= 0.37)
01162         {
01163           G4double x = std::log(hLabMomentum/0.48);
01164           TotN = 26. + 110.*x*x;
01165         }
01166       else 
01167         {
01168           G4double x = (hLabMomentum - 0.29)/0.07;
01169           TotN = 28.0 + 40.*std::exp(-x*x);
01170         }
01171       HadrTot = (TotP+TotN)/2;
01172 //  ........................................
01173       HadrSlope = 7.28+0.245*logS;        // GeV-2
01174       HadrReIm  = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
01175 
01176       DDSect2   = 0.7;                               //mb*GeV-2
01177       DDSect3   = 0.27;                              //mb*GeV-2
01178 
01179       break;
01180 //  ==========================================================
01181     case 4:            //  K plus
01182 
01183       HadrTot   = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55);  // mb
01184       if(HadrEnergy>100) HadrSlope = 15.0;
01185       else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;   // GeV-2
01186 
01187       HadrReIm  = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
01188       DDSect2   = 0.7;                             //mb*GeV-2
01189       DDSect3   = 0.21;                            //mb*GeV-2
01190       break;
01191 //  =========================================================
01192      case 5:              //   K minus
01193 
01194        HadrTot   = 10+1.8*logE + 25./sqrS; // mb
01195        HadrSlope = 6.98+0.127*logS;         // GeV-2
01196        HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
01197        DDSect2   = 0.7;                             //mb*GeV-2
01198        DDSect3   = 0.27;                            //mb*GeV-2
01199        break;
01200   }   
01201 //  =========================================================
01202   if(verboseLevel>2)
01203     G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
01204            << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
01205            << " DDSect3= " << DDSect3 << G4endl;
01206 
01207   if(Z != 1) return;
01208 
01209   // Scattering of protons
01210 
01211   Coeff0 = Coeff1 = Coeff2 = 0.0;
01212   Slope0 = Slope1 = 1.0;
01213   Slope2 = 5.0;
01214 
01215   // data for iHadron=0
01216   static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
01217   static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
01218   static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
01219   static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
01220   static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
01221       
01222   // data for iHadron=6,7
01223   static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
01224   static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
01225   static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
01226   static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
01227   static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
01228 
01229   // data for iHadron=1
01230   static const G4double EnP[2]={1.5,4.0};
01231   static const G4double C0P[2]={0.001,0.0005};
01232   static const G4double C1P[2]={0.003,0.001};
01233   static const G4double B0P[2]={2.5,4.5};
01234   static const G4double B1P[2]={1.0,4.0};
01235 
01236   // data for iHadron=2
01237   static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
01238   static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
01239   static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
01240   static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
01241   static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
01242 
01243   // data for iHadron=3
01244   static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
01245   static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
01246   static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
01247   static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
01248   static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
01249 
01250   // data for iHadron=4
01251   static const G4double EnK[4]={1.4,2.33,3.0,5.0};
01252   static const G4double C0K[4]={0.0,0.0,0.0,0.0};
01253   static const G4double C1K[4]={0.01,0.007,0.005,0.003};
01254   static const G4double B0K[4]={1.5,2.0,3.8,3.8};
01255   static const G4double B1K[4]={1.6,1.6,1.6,1.6};
01256 
01257   // data for iHadron=5
01258   static const G4double EnKM[2]={1.4,4.0};
01259   static const G4double C0KM[2]={0.006,0.002};
01260   static const G4double C1KM[2]={0.00,0.00};
01261   static const G4double B0KM[2]={2.5,3.5};
01262   static const G4double B1KM[2]={1.6,1.6};
01263 
01264   switch(iHadron)
01265     {
01266     case 0 :
01267 
01268       if(hLabMomentum <BoundaryP[0])
01269         InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
01270 
01271       Coeff2 = 0.8/hLabMomentum/hLabMomentum;
01272       break; 
01273 
01274     case  6 :
01275 
01276       if(hLabMomentum < BoundaryP[1])
01277         InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
01278 
01279       Coeff2 = 0.8/hLabMomentum/hLabMomentum;
01280       break; 
01281 
01282     case 1 :
01283     case 7 :
01284       if(hLabMomentum <  BoundaryP[2])
01285         InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
01286       break; 
01287 
01288     case 2 :
01289 
01290       if(hLabMomentum < BoundaryP[3])
01291         InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
01292 
01293       Coeff2 = 0.02/hLabMomentum;
01294       break; 
01295 
01296     case 3 :
01297 
01298       if(hLabMomentum < BoundaryP[4])
01299         InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
01300 
01301       Coeff2 = 0.02/hLabMomentum;
01302       break;
01303  
01304     case 4 :
01305 
01306       if(hLabMomentum < BoundaryP[5])
01307         InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
01308 
01309       if(hLabMomentum < 1) Coeff2 = 0.34;
01310       else  Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
01311       break; 
01312 
01313     case 5 :
01314       if(hLabMomentum < BoundaryP[6])
01315         InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
01316 
01317       if(hLabMomentum < 1) Coeff2 = 0.01;
01318       else  Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
01319       break; 
01320     }
01321 
01322   if(verboseLevel > 2) 
01323     G4cout<<"  HadrVal : Plasb  "<<hLabMomentum
01324           <<"  iHadron  "<<iHadron<<"  HadrTot  "<<HadrTot<<G4endl;
01325 }
01326 
01327 //  =====================================================
01328 void  G4ElasticHadrNucleusHE::
01329        GetKinematics(const G4ParticleDefinition * aHadron,
01330                            G4double MomentumH)
01331 {
01332   if (verboseLevel>1)
01333     G4cout<<"1  GetKin.: HadronName MomentumH "
01334           <<aHadron->GetParticleName()<<"  "<<MomentumH<<G4endl;
01335 
01336   DefineHadronValues(1);
01337 
01338   G4double Sh     = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
01339 
01340   ConstU = 2*protonM2+2*hMass2-Sh;
01341 
01342   G4double MaxT = 4*MomentumCM*MomentumCM;
01343 
01344   BoundaryTL[0] = MaxT; //2.0;
01345   BoundaryTL[1] = MaxT;
01346   BoundaryTL[3] = MaxT;
01347   BoundaryTL[4] = MaxT;
01348   BoundaryTL[5] = MaxT;
01349 
01350   G4int NumberH=0;
01351 
01352   while(iHadrCode!=HadronCode[NumberH]) NumberH++;
01353 
01354   NumberH = HadronType1[NumberH];   
01355 
01356   if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
01357   else MaxTR = BoundaryTG[NumberH];
01358 
01359   if (verboseLevel>1)
01360     G4cout<<"3  GetKin. : NumberH  "<<NumberH
01361           <<"  Bound.P[NumberH] "<<BoundaryP[NumberH]
01362           <<"  Bound.TL[NumberH] "<<BoundaryTL[NumberH]
01363           <<"  Bound.TG[NumberH] "<<BoundaryTG[NumberH]
01364           <<"  MaxT MaxTR "<<MaxT<<"  "<<MaxTR<<G4endl;
01365 
01366 //     GetParametersHP(aHadron, MomentumH);
01367 }
01368 //  ============================================================
01369 G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
01370 {
01371   G4double Fdistr=0;
01372   G4double SqrQ2 = std::sqrt(Q2);
01373  
01374   Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
01375     /HadrSlope*(1-std::exp(-HadrSlope*Q2))
01376     + Coeff0*(1-std::exp(-Slope0*Q2))
01377     + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
01378     + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
01379 
01380   if (verboseLevel>1)
01381     G4cout<<"Old:  Coeff0 Coeff1 Coeff2 "<<Coeff0<<"  "
01382           <<Coeff1<<"  "<<Coeff2<<"  Slope Slope0 Slope1 Slope2 "
01383           <<HadrSlope<<"  "<<Slope0<<"  "<<Slope1<<"  "<<Slope2
01384           <<"  Fdistr "<<Fdistr<<G4endl;
01385   return Fdistr;
01386 }
01387 //  +++++++++++++++++++++++++++++++++++++++
01388 G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran)
01389 {
01390   G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
01391   G4double Q2=0;
01392 
01393   FmaxT = GetFt(MaxTR);
01394   delta = GetDistrFun(DDD0)-Ran;
01395 
01396   while(std::fabs(delta) > 0.0001)
01397     {
01398       if(delta>0) 
01399       {
01400         DDD2 = DDD0;
01401         DDD0 = (DDD0+DDD1)*0.5;
01402       }
01403       else if(delta<0)
01404       {
01405         DDD1 = DDD0; 
01406         DDD0 = (DDD0+DDD2)*0.5;
01407       }
01408       delta = GetDistrFun(DDD0)-Ran;
01409     }
01410 
01411   Q2 = DDD0;
01412 
01413   return Q2;
01414 }
01415 //  ++++++++++++++++++++++++++++++++++++++++++
01416 G4double G4ElasticHadrNucleusHE::
01417        HadronProtonQ2(const G4ParticleDefinition * p,
01418                             G4double inLabMom)
01419 {
01420 
01421   hMass         = p->GetPDGMass()/GeV;
01422   hMass2        = hMass*hMass;
01423   hLabMomentum  = inLabMom;
01424   hLabMomentum2 = hLabMomentum*hLabMomentum;
01425   HadrEnergy    = sqrt(hLabMomentum2+hMass2); 
01426 
01427   G4double Rand = G4UniformRand();
01428 
01429   GetKinematics(p, inLabMom);
01430 
01431   G4double Q2 = GetQ2(Rand);
01432 
01433   return Q2;
01434 }
01435 
01436 //  ===========================================
01437 
01439 //
01440 //  
01441 
01442 void G4ElasticHadrNucleusHE::Binom()
01443 {
01444   G4int N, M;
01445   G4double  Fact1, J;
01446 
01447   for(N = 0; N < 240; N++)
01448   {
01449     J = 1;
01450 
01451       for( M = 0; M <= N; M++ )
01452       {
01453           Fact1 = 1;
01454 
01455           if ( ( N > 0 ) && ( N > M ) && M > 0 )
01456           {
01457               J *= G4double(N-M+1)/G4double(M);
01458               Fact1 = J;
01459           }
01460           SetBinom[N][M] = Fact1;
01461       }
01462   }
01463   return;
01464 }
01465 
01466 
01467 //
01468 //
01470 

Generated on Mon May 27 17:48:06 2013 for Geant4 by  doxygen 1.4.7