// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // $Id: G4QTauNuclearCrossSection.cc,v 1.2 2010/06/02 09:08:25 mkossov Exp $ // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ // // // G4 Physics class: G4QTauNuclearCrossSection for gamma+A cross sections // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-Oct-03 // // **************************************************************************************** // ********** This CLASS is temporary moved from the photolepton_hadron directory ********* // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ****** // **************************************************************************************** //========================================================================================= // Short description: reaction cross-sections for tau-nuclear reactions, which // are integrals over virtual equivalent photons photons. The tau-nuclear // reactions do not exist in GHAD, so by the present physics lists it is not // simulated at all. // -------------------------------------------------------------------------------- ///#define debug #define edebug //#define pdebug //#define ppdebug //#define tdebug //#define sdebug #include "G4QTauNuclearCrossSection.hh" // Initialization of the G4bool G4QTauNuclearCrossSection::onlyCS=true;// Flag to calculate only CS G4double G4QTauNuclearCrossSection::lastSig=0.;// Last calculated cross section G4int G4QTauNuclearCrossSection::lastL=0; // Last used in cross section TheLastBin G4double G4QTauNuclearCrossSection::lastE=0.; // Last used in cross section TheEnergy G4int G4QTauNuclearCrossSection::lastF=0; // Last used in cross section TheFirstBin G4double G4QTauNuclearCrossSection::lastG=0.; // Last used in cross section TheGamma G4double G4QTauNuclearCrossSection::lastH=0.; // LastValue of theHighEnergy A-dependence G4double* G4QTauNuclearCrossSection::lastJ1=0; // Pointer to the LastArray of J1 function G4double* G4QTauNuclearCrossSection::lastJ2=0; // Pointer to the LastArray of J2 function G4double* G4QTauNuclearCrossSection::lastJ3=0; // Pointer to the LastArray of J3 function G4int G4QTauNuclearCrossSection::lastPDG=0; // The last PDG code of the projectile G4int G4QTauNuclearCrossSection::lastN=0; // The last N of calculated nucleus G4int G4QTauNuclearCrossSection::lastZ=0; // The last Z of calculated nucleus G4double G4QTauNuclearCrossSection::lastP=0.; // Last used in cross section Momentum G4double G4QTauNuclearCrossSection::lastTH=0.; // Last threshold momentum G4double G4QTauNuclearCrossSection::lastCS=0.; // Last value of the Cross Section G4int G4QTauNuclearCrossSection::lastI=0; // The last position in the DAMDB std::vector* G4QTauNuclearCrossSection::J1 = new std::vector; std::vector* G4QTauNuclearCrossSection::J2 = new std::vector; std::vector* G4QTauNuclearCrossSection::J3 = new std::vector; // Returns Pointer to the G4VQCrossSection class G4VQCrossSection* G4QTauNuclearCrossSection::GetPointer() { static G4QTauNuclearCrossSection theCrossSection; //**Static body of the Cross Section** return &theCrossSection; } G4QTauNuclearCrossSection::~G4QTauNuclearCrossSection() { G4int lens=J1->size(); for(G4int i=0; isize(); for(G4int i=0; isize(); for(G4int i=0; i colPDG;// Vector of the projectile PDG code static std::vector colN; // Vector of N for calculated nuclei (isotops) static std::vector colZ; // Vector of Z for calculated nuclei (isotops) static std::vector colP; // Vector of last momenta for the reaction static std::vector colTH; // Vector of energy thresholds for the reaction static std::vector colCS; // Vector of last cross sections for the reaction // ***---*** End of the mandatory Static Definitions of the Associative Memory ***---*** G4double pEn=std::sqrt(pMom*pMom+mtu2)-mtu; // ==> tau-/tau+ kinEnergy #ifdef pdebug G4cout<<"G4QTNCS::GetCS:>>> f="< CS=0"< New (inDB) Calculated CS="<lastTH) // Correct the threshold { #ifdef pdebug G4cout<<"G4QTNCS::GetCS: New T="< Threshold="< Threshold="< CS=0 } else if(std::fabs(lastP/pMom-1.)m+(m^2+2lm)/2M // CHIPS - Direct GEANT //G4double mT= G4QPDGCode(111).GetNuclMass(Z,N,0); G4double mT= 0.; if(G4NucleiProperties::IsInStableTable(A,Z)) mT = G4NucleiProperties::GetNuclearMass(A,Z)/MeV; else return 0.; // If it is not in the Table of Stable Nuclei, then the Threshold=0 // --------- Splitting thresholds G4double mP= infEn; if(A>1&&Z&&G4NucleiProperties::IsInStableTable(A-1,Z-1)) mP = G4NucleiProperties::GetNuclearMass(A-1,Z-1)/MeV; // ResNucMass for a proton G4double mN= infEn; if(A>1&&G4NucleiProperties::IsInStableTable(A-1,Z)) mN = G4NucleiProperties::GetNuclearMass(A-1,Z)/MeV; // ResNucMass for a neutron G4double mA= infEn; if(A>4&&Z>1&&G4NucleiProperties::IsInStableTable(A-4,Z-2)) mA = G4NucleiProperties::GetNuclearMass(A-4,Z-2)/MeV; // ResNucMass for an alpha G4double dP= mP +mProt - mT; G4double dN= mN +mNeut - mT; G4double dA= mA +mAlph - mT; #ifdef pdebug G4cout<<"G4TauNucCS::ThreshEn: mP="<size(); if(J3->size()) G4cout<<", p="<<(*J3)[0]; G4cout<GetKineticEnergy()/MeV; // Energy of the Tau-lepton onlyCS=CS; // Flag to calculate only CS (not Si/Bi) G4double TotEnergy2=Momentum*Momentum+mtu2; G4double TotEnergy=std::sqrt(TotEnergy2); // Total energy of the muon lastE=TotEnergy-mtu; // Kinetic energy of the muon #ifdef pdebug G4cout<<"G4QElectronNucCS::CalcCS: P="< G4QTauNucCS::CalcCS: Old CS=0 as lastE="< G4QTauNucCS::CalcCS: CS=0 as lastE="< CREATE { lastJ1 = new G4double[nE]; // Allocate memory for the new J1 function lastJ2 = new G4double[nE]; // Allocate memory for the new J2 function lastJ3 = new G4double[nE]; // Allocate memory for the new J3 function lastF = GetFunctions(A,lastJ1,lastJ2,lastJ3);//newZeroPos and J-functions filling lastH = alop*A*(1.-.072*std::log(A)); // like lastSP of G4PhotonuclearCrossSection #ifdef pdebug G4cout<<"==>G4QTaNCS::CalcCS:lJ1="<size(); if(sync!=I) G4cerr<<"***G4QTauNuclearCS::CalcCS: PDG=15, S="<push_back(lastJ1); J2->push_back(lastJ2); J3->push_back(lastJ3); colF.push_back(lastF); colH.push_back(lastH); } // End of creation of the new set of parameters } // End of parameters udate // ============================== NOW Calculate the Cross Section ===================== if (lastE<=lastTH || lastE<=EMi) // Check that muKiE is higher than ThreshE { lastE=0.; lastG=0.; lastSig=0.; #ifdef pdebug G4cout<<"---> G4QTauNucCS::CalcCS:CS=0 as T="<(shift); #ifdef pdebug G4cout<<"-->G4QTauNuclearCS::CalcCrossSect:LOGfit b="<G4QTauNucCS::CCS:LOGex="<(a+.499); // Make the round integer of the atomic number G4double ai=iA; if(a!=ai) a=ai; for(G4int i=0; i get from Tab { for(G4int k=0; k must be calculated { G4int k=0; // !! To be good for different compilers !! for(k=1; k=nN) k=nN-1; // Extrapolation from the last bin (U) G4int k1=k-1; G4double xi=A[k1]; G4double b=(a-xi)/(A[k]-xi); for(G4int m=0; m Make them general. static const G4int nE=336; // !! If change this, change it in GetFunctions() (*.hh) !! static const G4int mL=nE-1; static const G4double EMi=2.0612; // Minimum Energy static const G4double EMa=50000.; // Maximum Energy static const G4double lEMi=std::log(EMi); // Minimum logarithmic Energy static const G4double lEMa=std::log(EMa); // Maximum logarithmic Energy static const G4double dlnE=(lEMa-lEMi)/mL; // Logarithmic step in Energy static const G4double mtu=1777.; // Mass of muon in MeV static const G4double lmtu=std::log(mtu); // Log of muon mass G4double phLE=0.; // Prototype of the log(nu=E_gamma) G4double Y[nE]; // Prepare the array for randomization #ifdef debug G4cout<<"G4QTauNuclCrossSect::GetExchanEn: B="<leE" <lE="<eE"); } if(std::fabs(d)=imax) G4cerr<<"G4QTauNucCS::SolveTheE:"<"<Use bigMax ln(eE)=" < Q2max="<