Geant4: G4FTFParameters.cc Source 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 // GEANT4 tag $Name:  $
00029 //
00030 
00031 #include <utility>                                        
00032 
00033 #include "G4FTFParameters.hh"
00034 
00035 #include "G4ios.hh"
00036 #include "G4PhysicalConstants.hh"
00037 #include "G4SystemOfUnits.hh"
00038 
00039 #include "G4ParticleDefinition.hh"             // 31 May 2011
00040 
00041 #include "G4Proton.hh"                         // 31 May 2011
00042 #include "G4Neutron.hh"                        // 31 May 2011
00043 
00044 #include "G4PionPlus.hh"                       // 31 May 2011
00045 #include "G4PionMinus.hh"                      // 31 May 2011
00046 #include "G4KaonPlus.hh"                       // 31 May 2011
00047 #include "G4KaonMinus.hh"                      // 31 May 2011
00048 
00049 G4FTFParameters::G4FTFParameters() :
00050   //A.R. 14-Aug-2012 : Coverity fix.
00051   FTFhNcmsEnergy(0.0), 
00052   FTFxsManager(0),
00053   FTFXtotal(0.0), FTFXelastic(0.0), FTFXinelastic(0.0), FTFXannihilation(0.0),
00054   ProbabilityOfAnnihilation(0.0), ProbabilityOfElasticScatt(0.0),
00055   RadiusOfHNinteractions2(0.0), FTFSlope(0.0), 
00056   AvaragePt2ofElasticScattering(0.0), FTFGamma0(0.0),
00057   MagQuarkExchange(0.0), SlopeQuarkExchange(0.0), DeltaProbAtQuarkExchange(0.0),
00058   ProbOfSameQuarkExchange(0.0), ProjMinDiffMass(0.0), ProjMinNonDiffMass(0.0),
00059   ProbabilityOfProjDiff(0.0), TarMinDiffMass(0.0), TarMinNonDiffMass(0.0),
00060   ProbabilityOfTarDiff(0.0), AveragePt2(0.0), ProbLogDistr(0.0),
00061   Pt2kink(0.0),
00062   MaxNumberOfCollisions(0.0), ProbOfInelInteraction(0.0), CofNuclearDestruction(0.0),
00063   R2ofNuclearDestruction(0.0), ExcitationEnergyPerWoundedNucleon(0.0),
00064   DofNuclearDestruction(0.0), Pt2ofNuclearDestruction(0.0), MaxPt2ofNuclearDestruction(0.0) 
00065 {}
00066 
00067 
00068 G4FTFParameters::~G4FTFParameters()
00069 {}
00070 //**********************************************************************************************
00071 G4FTFParameters::G4FTFParameters(const G4ParticleDefinition * particle, 
00072                                                    G4int      theA,
00073                                                    G4int      theZ,
00074                                                    G4double   PlabPerParticle) 
00075 {
00076 
00077     //A.R. 25-Jul-2012 Coverity fix.
00078     FTFXannihilation = 0.0;
00079     FTFhNcmsEnergy = 0.0;
00080     ProbOfSameQuarkExchange = 0.0;
00081 
00082     G4int    ProjectilePDGcode    = particle->GetPDGEncoding();
00083     G4int    ProjectileabsPDGcode = std::abs(ProjectilePDGcode);
00084     G4double ProjectileMass       = particle->GetPDGMass();
00085     G4double ProjectileMass2      =ProjectileMass*ProjectileMass;
00086 
00087     G4int    ProjectileBaryonNumber(0), AbsProjectileBaryonNumber(0);
00088     G4int                               AbsProjectileCharge(0);
00089     G4bool   ProjectileIsNucleus=false;
00090 
00091     if(std::abs(particle->GetBaryonNumber()) > 1)
00092     { // The projectile is a nucleus
00093      ProjectileIsNucleus      =true;
00094      ProjectileBaryonNumber   =particle->GetBaryonNumber();
00095      AbsProjectileBaryonNumber=std::abs(ProjectileBaryonNumber);
00096      AbsProjectileCharge      =(G4int) particle->GetPDGCharge();
00097 
00098      if(ProjectileBaryonNumber > 1)
00099      {      ProjectilePDGcode= 2212; ProjectileabsPDGcode=2212;} // Proton
00100      else { ProjectilePDGcode=-2212; ProjectileabsPDGcode=2212;} // Anti-Proton
00101 
00102      ProjectileMass =G4Proton::Proton()->GetPDGMass();
00103      ProjectileMass2=sqr(ProjectileMass);
00104     } 
00105 
00106     G4double TargetMass     = G4Proton::Proton()->GetPDGMass();
00107     G4double TargetMass2    = TargetMass*TargetMass;
00108 
00109     G4double Plab = PlabPerParticle;
00110     G4double Elab = std::sqrt(Plab*Plab+ProjectileMass2);
00111     G4double KineticEnergy = Elab-ProjectileMass;                 // 31 May 2011
00112 
00113     G4double S=ProjectileMass2 + TargetMass2 + 2.*TargetMass*Elab;
00114 
00115 //G4cout<<"Proj Plab "<<ProjectilePDGcode<<" "<<Plab<<G4endl;
00116 //G4cout<<"Mass KinE "<<ProjectileMass<<" "<<KineticEnergy<<G4endl;
00117 //G4cout<<" A Z "<<theA<<" "<<theZ<<G4endl;
00118 
00119     G4double Ylab,Xtotal,Xelastic,Xannihilation;
00120     G4int NumberOfTargetNucleons;
00121 
00122     Ylab=0.5*std::log((Elab+Plab)/(Elab-Plab));
00123 
00124     G4double ECMSsqr=S/GeV/GeV;
00125     G4double SqrtS  =std::sqrt(S)/GeV;
00126 //G4cout<<"Sqrt(s) "<<SqrtS<<G4endl;
00127 
00128     TargetMass     /=GeV; TargetMass2     /=(GeV*GeV);
00129     ProjectileMass /=GeV; ProjectileMass2 /=(GeV*GeV);
00130 
00131     static G4ChipsComponentXS* _instance = new G4ChipsComponentXS();  // Witek Pokorski
00132     FTFxsManager = _instance;
00133 
00134     Plab/=GeV;
00135 //  G4double LogPlab    = std::log( Plab );
00136 //  G4double sqrLogPlab = LogPlab * LogPlab;
00137 
00138     G4int NumberOfTargetProtons  = theZ; 
00139     G4int NumberOfTargetNeutrons = theA-theZ;
00140 
00141     NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
00142 
00143     if( (ProjectilePDGcode == 2212) || 
00144         (ProjectilePDGcode == 2112)   )    //------Projectile is nucleon --------
00145       {        
00146        G4double XtotPP = FTFxsManager->
00147                   GetTotalElementCrossSection(  particle,KineticEnergy,1,0);
00148        G4ParticleDefinition* Neutron=G4Neutron::Neutron();
00149        G4double XtotPN = FTFxsManager->
00150                   GetTotalElementCrossSection(   Neutron,KineticEnergy,1,0);
00151 
00152 
00153        G4double XelPP  = FTFxsManager->
00154                   GetElasticElementCrossSection(particle,KineticEnergy,1,0);
00155        G4double XelPN  = FTFxsManager->
00156                   GetElasticElementCrossSection(   Neutron,KineticEnergy,1,0);
00157 //G4cout<<"Xs "<<XtotPP/millibarn<<" "<<XelPP/millibarn<<G4endl;
00158 //G4cout<<"Xs "<<XtotPN/millibarn<<" "<<XelPN/millibarn<<G4endl;
00159        if(!ProjectileIsNucleus)
00160        { // Projectile is hadron
00161         Xtotal          = ( NumberOfTargetProtons  * XtotPP + 
00162                             NumberOfTargetNeutrons * XtotPN  ) / NumberOfTargetNucleons;
00163         Xelastic        = ( NumberOfTargetProtons  * XelPP  + 
00164                             NumberOfTargetNeutrons * XelPN   ) / NumberOfTargetNucleons;
00165        } else
00166        { // Projectile is a nucleus
00167         Xtotal  = (
00168                  AbsProjectileCharge                           *NumberOfTargetProtons *XtotPP + 
00169                 (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons*XtotPP + 
00170                ( AbsProjectileCharge                           *NumberOfTargetNeutrons +
00171                 (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons)*XtotPN
00172                    )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
00173 
00174         Xelastic= (
00175                  AbsProjectileCharge                           *NumberOfTargetProtons *XelPP + 
00176                 (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons*XelPP + 
00177                ( AbsProjectileCharge                           *NumberOfTargetNeutrons +
00178                 (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons)*XelPN
00179                    )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
00180       }
00181 
00182        Xannihilation   = 0.;
00183        Xtotal/=millibarn;
00184        Xelastic/=millibarn;
00185       }
00186     else if( ProjectilePDGcode < -1000 )         //------Projectile is anti_baryon --------
00187       {        
00188 
00189        G4double X_a(0.), X_b(0.), X_c(0.), X_d(0.);
00190        G4double MesonProdThreshold=ProjectileMass+TargetMass+(2.*0.14+0.016); // 2 Mpi +DeltaE;
00191 
00192        if(PlabPerParticle < 40.*MeV)
00193        { // Low energy limits. Projectile at rest.
00194         Xtotal=   1512.9;    // mb
00195         Xelastic=  473.2;    // mb
00196         X_a=       625.1;    // mb
00197         X_b=         9.780;  // mb
00198         X_c=        49.989;  // mb
00199         X_d=         6.614;  // mb
00200        }
00201        else
00202        { // Total and elastic cross section of PbarP interactions a'la Arkhipov
00203         G4double LogS=std::log(ECMSsqr/33.0625);
00204         G4double Xasmpt=36.04+0.304*LogS*LogS;    // mb
00205 
00206                  LogS=std::log(SqrtS/20.74);
00207         G4double Basmpt=11.92+0.3036*LogS*LogS;   // GeV^(-2)
00208         G4double R0=std::sqrt(0.40874044*Xasmpt-Basmpt); // GeV^(-1)
00209 
00210         G4double FlowF=SqrtS/
00211         std::sqrt(ECMSsqr*ECMSsqr+ProjectileMass2*ProjectileMass2+TargetMass2*TargetMass2-
00212         2.*ECMSsqr*ProjectileMass2 -2.*ECMSsqr*TargetMass2 -2.*ProjectileMass2*TargetMass2);
00213 
00214         Xtotal=Xasmpt*(1.+13.55*FlowF/R0/R0/R0*
00215                          (1.-4.47/SqrtS+12.38/ECMSsqr-12.43/SqrtS/ECMSsqr)); // mb
00216 
00217         Xasmpt=4.4+0.101*LogS*LogS;    // mb
00218         Xelastic=Xasmpt*(1.+59.27*FlowF/R0/R0/R0*
00219                          (1.-6.95/SqrtS+23.54/ECMSsqr-25.34/SqrtS/ECMSsqr)); // mb
00220 //G4cout<<"Param Xtotal Xelastic "<<Xtotal<<" "<<Xelastic<<G4endl;
00221 //G4cout<<"FlowF "<<FlowF<<" SqrtS "<<SqrtS<<G4endl;
00222 //G4cout<<"Param Xelastic-NaN "<<Xelastic<<" "<<1.5*16.654/pow(ECMSsqr/2.176/2.176,2.2)<<" "<<ECMSsqr<<G4endl;
00223         X_a=25.*FlowF;               // mb, 3-shirts diagram
00224 
00225         if(SqrtS < MesonProdThreshold)
00226         {
00227          X_b=3.13+140.*std::pow(MesonProdThreshold-SqrtS,2.5);// mb anti-quark-quark annihilation
00228          Xelastic-=3.*X_b;  // Xel-X(PbarP->NNbar)
00229         } else
00230         {
00231          X_b=6.8/SqrtS;                                 // mb anti-quark-quark annihilation
00232          Xelastic-=3.*X_b;  // Xel-X(PbarP->NNbar)
00233         }
00234 
00235         X_c=2.*FlowF*sqr(ProjectileMass+TargetMass)/ECMSsqr; // mb rearrangement
00236 
00237 //G4cout<<"Old new Xa "<<35.*FlowF<<" "<<25.*FlowF<<G4endl;
00238 
00239         X_d=23.3/ECMSsqr;                       // mb anti-quark-quark string creation
00240        }
00241 //---------------------------------------------------------------
00242 //G4cout<<"Param Xtotal Xelastic "<<Xtotal<<" "<<Xelastic<<G4endl;
00243 //G4cout<<"Para a b c d "<<X_a<<" "<<X_b<<" "<<X_c<<" "<<X_d<<G4endl;
00244 //G4cout<<"Para a b c d "<<X_a<<" "<<5.*X_b<<" "<<5.*X_c<<" "<<6.*X_d<<G4endl;
00245        G4double Xann_on_P(0.), Xann_on_N(0.);
00246 
00247        if(ProjectilePDGcode == -2212)       // Pbar+P/N
00248        {Xann_on_P=X_a + X_b*5. + X_c*5. + X_d*6.; Xann_on_N=X_a + X_b*4. + X_c*4. + X_d*4.;} 
00249        else if(ProjectilePDGcode == -2112) // NeutrBar+P/N
00250        {Xann_on_P=X_a + X_b*4. + X_c*4. + X_d*4.; Xann_on_N=X_a + X_b*5. + X_c*5. + X_d*6.;} 
00251        else if(ProjectilePDGcode == -3122) // LambdaBar+P/N
00252        {Xann_on_P=X_a + X_b*3. + X_c*3. + X_d*2.; Xann_on_N=X_a + X_b*3. + X_c*3. + X_d*2.;} 
00253        else if(ProjectilePDGcode == -3112) // Sigma-Bar+P/N
00254        {Xann_on_P=X_a + X_b*2. + X_c*2. + X_d*0.; Xann_on_N=X_a + X_b*4. + X_c*4. + X_d*2.;} 
00255        else if(ProjectilePDGcode == -3212) // Sigma0Bar+P/N
00256        {Xann_on_P=X_a + X_b*3. + X_c*3. + X_d*2.; Xann_on_N=X_a + X_b*3. + X_c*3. + X_d*2.;} 
00257        else if(ProjectilePDGcode == -3222) // Sigma+Bar+P/N
00258        {Xann_on_P=X_a + X_b*4. + X_c*4. + X_d*2.; Xann_on_N=X_a + X_b*2. + X_c*2. + X_d*0.;} 
00259        else if(ProjectilePDGcode == -3312) // Xi-Bar+P/N
00260        {Xann_on_P=X_a + X_b*1. + X_c*1. + X_d*0.; Xann_on_N=X_a + X_b*2. + X_c*2. + X_d*0.;} 
00261        else if(ProjectilePDGcode == -3322) // Xi0Bar+P/N
00262        {Xann_on_P=X_a + X_b*2. + X_c*2. + X_d*0.; Xann_on_N=X_a + X_b*1. + X_c*1. + X_d*0.;} 
00263        else if(ProjectilePDGcode == -3334) // Omega-Bar+P/N
00264        {Xann_on_P=X_a + X_b*0. + X_c*0. + X_d*0.; Xann_on_N=X_a + X_b*0. + X_c*0. + X_d*0.;} 
00265        else {G4cout<<"Unknown anti-baryon for FTF annihilation"<<G4endl;}
00266 //---------------------------------------------------------------
00267 
00268 //G4cout<<"Sum          "<<Xann_on_P<<G4endl;
00269 
00270        if(!ProjectileIsNucleus)
00271        { // Projectile is anti-baryon
00272         Xannihilation   = ( NumberOfTargetProtons  * Xann_on_P  + 
00273                             NumberOfTargetNeutrons * Xann_on_N   ) / NumberOfTargetNucleons;
00274        } else
00275        { // Projectile is a nucleus
00276         Xannihilation=(
00277           ( AbsProjectileCharge                           *NumberOfTargetProtons+ 
00278            (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons )*Xann_on_P + 
00279           ( AbsProjectileCharge                           *NumberOfTargetNeutrons+
00280            (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons  )*Xann_on_N
00281                       )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
00282        }
00283 
00284        G4double Xftf=0.;  
00285        MesonProdThreshold=ProjectileMass+TargetMass+(0.14+0.08); // Mpi +DeltaE
00286        if(SqrtS > MesonProdThreshold) {Xftf=36.*(1.-MesonProdThreshold/SqrtS);}
00287 
00288        Xtotal = Xelastic + Xannihilation + Xftf;
00289 /*
00290 G4cout<<"Plab Xtotal, Xelastic  Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
00291 G4cout<<"Plab Xelastic/Xtotal,  Xann/Xin "<<Plab<<" "<<Xelastic/Xtotal<<" "<<Xannihilation/(Xtotal-Xelastic)<<G4endl;
00292 //G4int Uzhi; G4cin>>Uzhi;
00293 */
00294 //---------------------------------------------------------------
00295       }
00296     else if( ProjectilePDGcode ==  211 )    //------Projectile is PionPlus -------
00297       {
00298        G4double XtotPiP = FTFxsManager->
00299                   GetTotalElementCrossSection(  particle,KineticEnergy,1,0); 
00300        G4ParticleDefinition* PionMinus=G4PionMinus::PionMinus();
00301        G4double XtotPiN =  FTFxsManager->
00302                   GetTotalElementCrossSection( PionMinus,KineticEnergy,1,0);
00303            
00304        G4double XelPiP  = FTFxsManager->
00305                   GetElasticElementCrossSection(particle,KineticEnergy,1,0); 
00306        G4double XelPiN  = FTFxsManager->
00307                   GetElasticElementCrossSection( PionMinus,KineticEnergy,1,0);
00308 
00309        Xtotal           = ( NumberOfTargetProtons  * XtotPiP + 
00310                             NumberOfTargetNeutrons * XtotPiN  ) / NumberOfTargetNucleons;
00311        Xelastic         = ( NumberOfTargetProtons  * XelPiP  + 
00312                             NumberOfTargetNeutrons * XelPiN   ) / NumberOfTargetNucleons; 
00313        Xannihilation   = 0.;
00314        Xtotal/=millibarn;
00315        Xelastic/=millibarn;
00316       }
00317     else if( ProjectilePDGcode == -211 )            //------Projectile is PionMinus -------
00318       {
00319        G4double XtotPiP = FTFxsManager->
00320                   GetTotalElementCrossSection(  particle,KineticEnergy,1,0);
00321        G4ParticleDefinition* PionPlus=G4PionPlus::PionPlus();
00322        G4double XtotPiN = FTFxsManager->
00323                   GetTotalElementCrossSection(  PionPlus,KineticEnergy,1,0);
00324            
00325        G4double XelPiP  = FTFxsManager->
00326                   GetElasticElementCrossSection(particle,KineticEnergy,1,0);
00327        G4double XelPiN  = FTFxsManager->
00328                   GetElasticElementCrossSection(  PionPlus,KineticEnergy,1,0);
00329 
00330        Xtotal           = ( NumberOfTargetProtons  * XtotPiP + 
00331                             NumberOfTargetNeutrons * XtotPiN  ) / NumberOfTargetNucleons;
00332        Xelastic         = ( NumberOfTargetProtons  * XelPiP  + 
00333                             NumberOfTargetNeutrons * XelPiN   ) / NumberOfTargetNucleons;
00334        Xannihilation   = 0.;
00335        Xtotal/=millibarn;
00336        Xelastic/=millibarn;
00337       }
00338 
00339     else if( ProjectilePDGcode ==  111 )          //------Projectile is PionZero  -------
00340       {
00341        G4ParticleDefinition* PionPlus=G4PionPlus::PionPlus();
00342        G4double XtotPipP= FTFxsManager->
00343                   GetTotalElementCrossSection(   PionPlus,KineticEnergy,1,0);
00344 
00345        G4ParticleDefinition* PionMinus=G4PionMinus::PionMinus();
00346        G4double XtotPimP= FTFxsManager->
00347                   GetTotalElementCrossSection(   PionMinus,KineticEnergy,1,0);
00348            
00349        G4double XelPipP = FTFxsManager->
00350                   GetElasticElementCrossSection(  PionPlus,KineticEnergy,1,0);
00351        G4double XelPimP = FTFxsManager->
00352                   GetElasticElementCrossSection( PionMinus,KineticEnergy,1,0);
00353 
00354        G4double XtotPiP= (XtotPipP + XtotPimP)/2.;
00355        G4double XtotPiN=XtotPiP;
00356        G4double XelPiP = (XelPipP  + XelPimP )/2.;
00357        G4double XelPiN = XelPiP;
00358 
00359        Xtotal           = ( NumberOfTargetProtons  * XtotPiP + 
00360                             NumberOfTargetNeutrons * XtotPiN  ) / NumberOfTargetNucleons;
00361        Xelastic         = ( NumberOfTargetProtons  * XelPiP  + 
00362                             NumberOfTargetNeutrons * XelPiN   ) / NumberOfTargetNucleons; 
00363        Xannihilation   = 0.;
00364 
00365        Xtotal/=millibarn;
00366        Xelastic/=millibarn;
00367       }
00368     else if( ProjectilePDGcode == 321 )             //------Projectile is KaonPlus -------
00369       {
00370        G4double XtotKP = FTFxsManager->
00371                   GetTotalElementCrossSection(  particle,KineticEnergy,1,0);
00372 
00373        G4ParticleDefinition* KaonMinus=G4KaonMinus::KaonMinus();
00374        G4double XtotKN = FTFxsManager->
00375                   GetTotalElementCrossSection( KaonMinus,KineticEnergy,1,0);
00376 
00377        G4double XelKP  = FTFxsManager->
00378                   GetElasticElementCrossSection(particle,KineticEnergy,1,0);
00379        G4double XelKN  = FTFxsManager->
00380                   GetElasticElementCrossSection( KaonMinus,KineticEnergy,1,0);
00381 
00382        Xtotal          = ( NumberOfTargetProtons  * XtotKP + 
00383                            NumberOfTargetNeutrons * XtotKN  ) / NumberOfTargetNucleons;
00384        Xelastic        = ( NumberOfTargetProtons  * XelKP  + 
00385                            NumberOfTargetNeutrons * XelKN   ) / NumberOfTargetNucleons;
00386        Xannihilation   = 0.;
00387 
00388        Xtotal/=millibarn;
00389        Xelastic/=millibarn;
00390       }
00391     else if( ProjectilePDGcode ==-321 )             //------Projectile is KaonMinus ------
00392       {
00393        G4double XtotKP = FTFxsManager->
00394                   GetTotalElementCrossSection(  particle,KineticEnergy,1,0);
00395 
00396        G4ParticleDefinition* KaonPlus=G4KaonPlus::KaonPlus();
00397        G4double XtotKN = FTFxsManager->
00398                   GetTotalElementCrossSection(  KaonPlus,KineticEnergy,1,0);
00399 
00400        G4double XelKP  = FTFxsManager->
00401                   GetElasticElementCrossSection(particle,KineticEnergy,1,0);
00402 
00403        G4double XelKN  = FTFxsManager->
00404                   GetElasticElementCrossSection(KaonPlus,KineticEnergy,1,0); 
00405 
00406        Xtotal          = ( NumberOfTargetProtons  * XtotKP + 
00407                            NumberOfTargetNeutrons * XtotKN  ) / NumberOfTargetNucleons;
00408        Xelastic        = ( NumberOfTargetProtons  * XelKP  + 
00409                            NumberOfTargetNeutrons * XelKN   ) / NumberOfTargetNucleons;
00410        Xannihilation   = 0.;
00411        
00412        Xtotal/=millibarn;
00413        Xelastic/=millibarn;
00414       }
00415     else if((ProjectilePDGcode == 311) || 
00416             (ProjectilePDGcode == 130) || 
00417             (ProjectilePDGcode == 310))               //Projectile is KaonZero
00418       {
00419        G4ParticleDefinition* KaonPlus=G4KaonPlus::KaonPlus();
00420        G4double XtotKpP= FTFxsManager->
00421                   GetTotalElementCrossSection(    KaonPlus,KineticEnergy,1,0);
00422 
00423        G4ParticleDefinition* KaonMinus=G4KaonMinus::KaonMinus();
00424        G4double XtotKmP= FTFxsManager->
00425                   GetTotalElementCrossSection(   KaonMinus,KineticEnergy,1,0);
00426 
00427        G4double XelKpP = FTFxsManager->
00428                   GetElasticElementCrossSection(  KaonPlus,KineticEnergy,1,0);
00429        G4double XelKmP = FTFxsManager->
00430                   GetElasticElementCrossSection(   KaonMinus,KineticEnergy,1,0);
00431 
00432        G4double XtotKP=(XtotKpP+XtotKmP)/2.;
00433        G4double XtotKN=XtotKP;
00434        G4double XelKP =(XelKpP +XelKmP )/2.; 
00435        G4double XelKN =XelKP;
00436 
00437        Xtotal          = ( NumberOfTargetProtons  * XtotKP + 
00438                            NumberOfTargetNeutrons * XtotKN  ) / NumberOfTargetNucleons;
00439        Xelastic        = ( NumberOfTargetProtons  * XelKP  + 
00440                            NumberOfTargetNeutrons * XelKN   ) / NumberOfTargetNucleons;
00441        Xannihilation   = 0.;
00442 
00443        Xtotal/=millibarn;
00444        Xelastic/=millibarn;
00445       }
00446     else                 //------Projectile is undefined, Nucleon assumed
00447       {
00448        G4ParticleDefinition* Proton=G4Proton::Proton();
00449        G4double XtotPP = FTFxsManager->
00450                   GetTotalElementCrossSection(  Proton,KineticEnergy,1,0);
00451 
00452        G4ParticleDefinition* Neutron=G4Neutron::Neutron();
00453        G4double XtotPN = FTFxsManager->
00454                   GetTotalElementCrossSection( Neutron,KineticEnergy,1,0);
00455 
00456        G4double XelPP  = FTFxsManager->
00457                   GetElasticElementCrossSection(Proton,KineticEnergy,1,0);
00458        G4double XelPN  = FTFxsManager->
00459                   GetElasticElementCrossSection( Neutron,KineticEnergy,1,0);
00460 
00461        Xtotal          = ( NumberOfTargetProtons  * XtotPP + 
00462                            NumberOfTargetNeutrons * XtotPN  ) / NumberOfTargetNucleons;
00463        Xelastic        = ( NumberOfTargetProtons  * XelPP  + 
00464                            NumberOfTargetNeutrons * XelPN   ) / NumberOfTargetNucleons;
00465        Xannihilation   = 0.;
00466 
00467        Xtotal/=millibarn;
00468        Xelastic/=millibarn;
00469       };
00470 
00471 //----------- Geometrical parameters ------------------------------------------------
00472       SetTotalCrossSection(Xtotal);
00473       SetElastisCrossSection(Xelastic);
00474       SetInelasticCrossSection(Xtotal-Xelastic);
00475 
00476 /*
00477 G4cout<<"Plab Xtotal, Xelastic Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
00478 if(Xtotal-Xelastic != 0.)
00479 {
00480   G4cout<<"Plab Xelastic/Xtotal,  Xann/Xin "<<Plab<<" "<<Xelastic/Xtotal<<" "<<Xannihilation/
00481   (Xtotal-Xelastic)<<G4endl;
00482 } else 
00483 {
00484   G4cout<<"Plab Xelastic/Xtotal,  Xann     "<<Plab<<" "<<Xelastic/Xtotal<<" "<<
00485   Xannihilation<<G4endl;
00486 }
00487 //G4int Uzhi; G4cin>>Uzhi;
00488 */
00489 //  // Interactions with elastic and inelastic collisions
00490       SetProbabilityOfElasticScatt(Xtotal, Xelastic);
00491       SetRadiusOfHNinteractions2(Xtotal/pi/10.);
00492 
00493       if(Xtotal-Xelastic == 0.)
00494       {
00495        SetProbabilityOfAnnihilation(0.);
00496       } else
00497       {SetProbabilityOfAnnihilation(Xannihilation/(Xtotal-Xelastic));}
00498 //
00499 //SetProbabilityOfElasticScatt(Xtotal, 0.);
00500 // //==== No elastic scattering ============================
00501 //      SetProbabilityOfElasticScatt(Xtotal, 0.);
00502 //      SetRadiusOfHNinteractions2((Xtotal-Xelastic)/pi/10.);
00503 //      SetProbabilityOfAnnihilation(1.);
00504 //        SetProbabilityOfAnnihilation(0.);
00505 // //=======================================================
00506 
00507 //-----------------------------------------------------------------------------------  
00508 
00509       SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
00510                                                         //      (GeV/c)^(-2))
00511 //G4cout<<"Slope "<<GetSlope()<<G4endl;
00512 //-----------------------------------------------------------------------------------
00513       SetGamma0( GetSlope()*Xtotal/10./2./pi );
00514 
00515 //----------- Parameters of elastic scattering --------------------------------------
00516                                                         // Gaussian parametrization of
00517                                                         // elastic scattering amplitude assumed
00518       SetAvaragePt2ofElasticScattering(1./(Xtotal*Xtotal/16./pi/Xelastic/0.3894)*GeV*GeV);
00519 //G4cout<<"AvaragePt2ofElasticScattering "<<GetAvaragePt2ofElasticScattering()<<G4endl;
00520 //----------- Parameters of excitations ---------------------------------------------
00521 
00522             G4double Xinel=Xtotal-Xelastic;                        // Uzhi 25.04.2012
00523 //G4cout<<"Param ProjectilePDGcode "<<ProjectilePDGcode<<G4endl;
00524            if( ProjectilePDGcode > 1000 )             //------Projectile is baryon --------
00525              {
00526               SetMagQuarkExchange(1.84);//(3.63);
00527               SetSlopeQuarkExchange(0.7);//(1.2);
00528               SetDeltaProbAtQuarkExchange(0.);
00529               if(NumberOfTargetNucleons > 26) {SetProbOfSameQuarkExchange(1.);}
00530               else                            {SetProbOfSameQuarkExchange(0.);}
00531 
00532               SetProjMinDiffMass(1.16);                              // GeV 
00533               SetProjMinNonDiffMass(1.16);                           // GeV 
00534 //            SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab));  // Uzhi 21.05.2012
00535               SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr);        // Uzhi 25.04.2012
00536 
00537               SetTarMinDiffMass(1.16);                               // GeV
00538               SetTarMinNonDiffMass(1.16);                            // GeV 
00539 //            SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab));   // Uzhi 21.05.2012
00540               SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr);         // Uzhi 25.04.2012
00541 //            SetAveragePt2(0.15);                                   // 0.15 GeV^2
00542               SetAveragePt2(0.3);                         // 0.30 GeV^2 Uzhi 21.05.2012
00543 
00544               SetProbLogDistr(0.5);                                  // Uzhi 21.05.2012
00545              }
00546            else if( ProjectilePDGcode < -1000 )  //------Projectile is anti_baryon --------
00547              {
00548               SetMagQuarkExchange(0.);
00549               SetSlopeQuarkExchange(0.);
00550               SetDeltaProbAtQuarkExchange(0.);
00551               SetProbOfSameQuarkExchange(0.);
00552 
00553               SetProjMinDiffMass(ProjectileMass+0.22);               // GeV 
00554               SetProjMinNonDiffMass(ProjectileMass+0.22);            // GeV
00555 //            SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab));  // Uzhi 21.05.2012
00556               SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr);        // Uzhi 25.04.2012
00557 
00558               SetTarMinDiffMass(TargetMass+0.22);                  // GeV
00559               SetTarMinNonDiffMass(TargetMass+0.22);               // GeV
00560 //            SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab));   // Uzhi 21.05.2012
00561               SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr);         // Uzhi 25.04.2012
00562 
00563               SetAveragePt2(0.3);                   // 0.15 GeV^2    // Uzhi 21.05.2012
00564 
00565               SetProbLogDistr(0.5);                                  // Uzhi 21.05.2012
00566              }
00567            else if( ProjectileabsPDGcode == 211 || 
00568                     ProjectilePDGcode ==  111)     //------Projectile is Pion -----------
00569              {
00570               SetMagQuarkExchange(240.); 
00571               SetSlopeQuarkExchange(2.);         
00572               SetDeltaProbAtQuarkExchange(0.56); //(0.35);
00573 
00574               SetProjMinDiffMass(0.5);                               // GeV
00575               SetProjMinNonDiffMass(0.5);                            // GeV 0.3
00576 //              SetProbabilityOfProjDiff(0.);                        // Uzhi 3.06.2012 
00577               SetProbabilityOfProjDiff((6.2-3.7*std::exp(-sqr(SqrtS-7.)/16.))/Xinel*0.);
00578 
00579               SetTarMinDiffMass(1.16);                               // GeV
00580               SetTarMinNonDiffMass(1.16);                            // GeV
00581 //            SetProbabilityOfTarDiff(0.8*std::exp(-0.6*(Ylab-3.))); // Uzhi 3.06.2012
00582               SetProbabilityOfTarDiff((2.+22./ECMSsqr)/Xinel);
00583 
00584               SetAveragePt2(0.3);                                   // GeV^2 7 June 2011
00585               SetProbLogDistr(0.);                                   // Uzhi 21.05.2012
00586 
00587               SetProbLogDistr(1.);
00588              }
00589            else if( (ProjectileabsPDGcode == 321) || 
00590                     (ProjectileabsPDGcode == 311) || 
00591                     (ProjectilePDGcode == 130)    || 
00592                     (ProjectilePDGcode == 310))        //Projectile is Kaon
00593              {
00594               SetMagQuarkExchange(40.);
00595               SetSlopeQuarkExchange(2.25);
00596               SetDeltaProbAtQuarkExchange(0.6);
00597 
00598               SetProjMinDiffMass(0.6);                               // GeV 0.7 0.6
00599               SetProjMinNonDiffMass(0.6);                            // GeV 0.7 0.6
00600 //            SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
00601               SetProbabilityOfProjDiff(0.*4.7/Xinel);                   // Uzhi 5.06.2012
00602 
00603               SetTarMinDiffMass(1.1);                                // GeV
00604               SetTarMinNonDiffMass(1.1);                             // GeV
00605 //            SetProbabilityOfTarDiff(0.45*std::pow(s/GeV/GeV,-0.5));// 40/32 X-dif/X-inel
00606               SetProbabilityOfTarDiff(1.5/Xinel);                    // Uzhi 5.06.2012
00607               SetAveragePt2(0.3);                                    // GeV^2 7 June 2011
00608               SetProbLogDistr(1.);                                   // Uzhi 5.06.2012
00609              }
00610            else                                           //------Projectile is undefined,
00611                                                           //------Nucleon assumed
00612              {
00613 /*                 // Uzhi 6.06.2012
00614               SetMagQuarkExchange(1.85);       // 7 June 2011
00615               SetSlopeQuarkExchange(0.7);      // 7 June 2011
00616               SetDeltaProbAtQuarkExchange(0.); // 7 June 2011
00617 
00618               SetProjMinDiffMass((940.+160.*MeV)/GeV);     // particle->GetPDGMass()
00619               SetProjMinNonDiffMass((940.+160.*MeV)/GeV);  // particle->GetPDGMass()
00620               SetProbabilityOfProjDiff(0.805*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
00621 
00622               SetTarMinDiffMass(1.16);                     // GeV
00623               SetTarMinNonDiffMass(1.16);                  // GeV
00624               SetProbabilityOfTarDiff(0.805*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
00625 */
00626 
00627               SetMagQuarkExchange(0.);
00628               SetSlopeQuarkExchange(0.);
00629               SetDeltaProbAtQuarkExchange(0.);
00630               SetProbOfSameQuarkExchange(0.);
00631 
00632               SetProjMinDiffMass(ProjectileMass+0.22);               // GeV 
00633               SetProjMinNonDiffMass(ProjectileMass+0.22);            // GeV
00634               SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr);        // Uzhi 25.04.2012
00635 
00636               SetTarMinDiffMass(TargetMass+0.22);                    // GeV
00637               SetTarMinNonDiffMass(TargetMass+0.22);                 // GeV
00638               SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr);         // Uzhi 25.04.2012
00639 
00640               SetAveragePt2(0.3);                         // 0.15 GeV^2 Uzhi 21.05.2012
00641               SetProbLogDistr(0.5);                                  // Uzhi 21.05.2012
00642              }
00643 
00644 //           if(theA > 4) SetProbabilityOfProjDiff(0.);     // Uzhi 6.07.2012 Closed
00645 
00646 //G4cout<<"Param Get Min Dif "<<GetProjMinNonDiffMass()<<G4endl;
00647 
00648 // ---------- Set parameters of a string kink -------------------------------
00649              SetPt2Kink(6.*GeV*GeV);
00650              G4double Puubar(1./3.), Pddbar(1./3.), Pssbar(1./3.); // SU(3) symmetry
00651 //           G4double Puubar(0.41 ), Pddbar(0.41 ), Pssbar(0.18 ); // Broken SU(3) symmetry
00652              SetQuarkProbabilitiesAtGluonSplitUp(Puubar, Pddbar, Pssbar);
00653 
00654 // --------- Set parameters of nuclear destruction--------------------
00655 
00656     if( ProjectileabsPDGcode < 1000 )               // Meson projectile
00657     {
00658       SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
00659       SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
00660                               (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
00661 
00662       SetR2ofNuclearDestruction(1.5*fermi*fermi);
00663 
00664       SetDofNuclearDestruction(0.3);
00665       SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
00666                                          (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
00667 //G4cout<<"Parm Pt2 Y "<<(0.035+0.04*std::exp(4.*(Ylab-2.5))/(1.+std::exp(4.*(Ylab-2.5))))<<" "<<Ylab<<G4endl;
00668       SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
00669 
00670       SetExcitationEnergyPerWoundedNucleon(100.*MeV);
00671     } else if( ProjectilePDGcode < -1000 )             // for anti-baryon projectile
00672     {
00673 //G4cout<<"Nucl destruct Anti Bar"<<G4endl;
00674 
00675       SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
00676       SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
00677                               (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
00678 
00679       SetR2ofNuclearDestruction(1.5*fermi*fermi);
00680 
00681       SetDofNuclearDestruction(0.3);
00682       SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
00683                                          (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
00684       SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
00685 
00686       SetExcitationEnergyPerWoundedNucleon(100.*MeV);
00687 
00688       if(Plab < 2.)   // 2 GeV/c
00689       { // For slow anti-baryon we have to garanty putting on mass-shell
00690        SetCofNuclearDestruction(0.);
00691        SetR2ofNuclearDestruction(1.5*fermi*fermi);
00692        SetDofNuclearDestruction(0.01);
00693        SetPt2ofNuclearDestruction(0.035*GeV*GeV);
00694        SetMaxPt2ofNuclearDestruction(0.04*GeV*GeV);
00695 //       SetExcitationEnergyPerWoundedNucleon(0.);   // ?????
00696       }
00697     } else                                        // Projectile baryon assumed
00698     {
00699       SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
00700       SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
00701                               (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
00702 
00703       SetR2ofNuclearDestruction(1.5*fermi*fermi);
00704 
00705       SetDofNuclearDestruction(0.3);
00706       SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
00707                                          (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
00708       SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
00709 
00710       SetExcitationEnergyPerWoundedNucleon(100.*MeV);
00711     }
00712 
00713 //SetCofNuclearDestruction(0.47*std::exp(2.*(Ylab-2.5))/(1.+std::exp(2.*(Ylab-2.5)))); 
00714 //SetPt2ofNuclearDestruction((0.035+0.1*std::exp(4.*(Ylab-3.))/(1.+std::exp(4.*(Ylab-3.))))*GeV*GeV);
00715 
00716 //SetMagQuarkExchange(120.); // 210. PipP
00717 //SetSlopeQuarkExchange(2.0);
00718 //SetDeltaProbAtQuarkExchange(0.6);
00719 //SetProjMinDiffMass(0.7);                    // GeV 1.1
00720 //SetProjMinNonDiffMass(0.7);                 // GeV
00721 //SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
00722 //SetTarMinDiffMass(1.1);                     // GeV
00723 //SetTarMinNonDiffMass(1.1);                  // GeV
00724 //SetProbabilityOfTarDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
00725 //
00726 //SetAveragePt2(0.3);                         // GeV^2
00727 //------------------------------------
00728 //SetProbabilityOfElasticScatt(1.,1.); //(Xtotal, Xelastic);
00729 //SetProbabilityOfProjDiff(1.*0.62*std::pow(s/GeV/GeV,-0.51)); // 0->1
00730 //SetProbabilityOfTarDiff(4.*0.62*std::pow(s/GeV/GeV,-0.51)); // 2->4
00731 //SetAveragePt2(0.3);                              //(0.15);
00732 //SetAvaragePt2ofElasticScattering(0.);
00733 
00734 //SetMaxNumberOfCollisions(Plab,6.); //(4.*(Plab+0.01),Plab); //6.); // ##########
00735 //SetAveragePt2(0.15);
00736 //G4cout<<"Cnd "<<GetCofNuclearDestruction()<<G4endl;
00737 //SetCofNuclearDestruction(0.4);// (0.2); //(0.4);                  
00738 //SetExcitationEnergyPerWoundedNucleon(0.*MeV); //(75.*MeV); 
00739 //SetDofNuclearDestruction(0.);                  
00740 //SetPt2ofNuclearDestruction(0.*GeV*GeV); //(0.168*GeV*GeV); 
00741 //G4cout<<"Pt2 "<<GetPt2ofNuclearDestruction()/GeV/GeV<<G4endl;
00742 //G4int Uzhi; G4cin>>Uzhi;
00743 } 
00744 //**********************************************************************************************