// ESgen.cc // Contact person: Nuno Barros // See ESgen.hh for more details //———————————————————————// /// /// Generates an neutrino-elastic scattering event, based on the /// cross-section as function of neutrino energy and the electron's /// recoil energy. Allow for variations in the weak mixing angle and /// the possibility of a neutrino magnetic moment /// /// J. Formaggio (UW) -02/09/2005 /// /// Converted to Geant4+GLG4Sim+RAT by Bill Seligman (07-Feb-2006). /// I'm following the code structure of the IBD classes: /// RATVertexGen_ES handles the G4-related tasks of constructing an /// event, while this class deals with the physics of the /// cross-section. Some of the code (the flux in particular) is copied /// from IBDgen. /// /// N. Barros (LIP) - 15/12/2010 /// Updated the code to use a separate ES cross section class ESCrossSec /// Added a few options to make this generator more general /// /// N. Barros (LIP) - 18/2/2011 /// The class was strongly rewritten to keep the possibility to used for it's original purpose /// The cross section calculations have been completely moved out /// and can be now used independently. /// Only the kinematics have been kept here. /// /// N. Barros (TUD) - 12-Dec-2011 /// /// Updated for easier integration with the SolarGen generator. /// Can still be used as part of the combo generator, /// but more options are necessary to be passed through messenger in that case. /// The decision for the input neutrino direction comes from VertexGen_ES and therefore this class is /// simply a wrapper for the cross section and the neutrino energy. /// /// N. Barros (TUD) - 02-Jul-2012 /// /// Refactored the code to be more dependent on the generators that use it and therefore avoid duplication. /// Removed the usage of Geant4 default streamers (G4cout) to use RAT Log mechanism instead /// Corrected some minor bugs related with the calculation of the recoil energy. /// Corrected a problem leading to underestimated event rate from continuous neutrino spectra. /// Removed attached messenger, as it is no longer needed. /// /// N. Barros (TUD) - 02-Aug-2012 /// /// - Changed the code so that the ES generator can (again) be used as part of the combo generator. /// - The functionality for the solar generator was maintained. /// /// N. Barros (UPenn) - 07-May-2014 /// /// - Renamed G4DEBUG to RATDEBUG to comply with the new naming scheme /// - Solved bug in the calculation of the rate per target for continuous spectra. /// /// J. Caravaca (UCBerkeley) - 14-Jan-2019 /// /// - Correct implementation of trapezoidal rule for calculation of xsec-weighted fluxes /// //////////////////////////////////////////////////////////////////// // - RAT includes #include #include #include #include // - G4 includes #include #include #include #include // - CLHEP includes #include // - ROOT includes #include #include #include using namespace CLHEP; namespace RAT { // This class is a helper to take care of the type of spectrum that is going to be used and // reads the RATDB entries accordingly. ESgen::ESgen() : fNuType(""), fNuFlavor(""),fXS(NULL), fNuSpectrum(NULL), fFluxMax(0.), fGenLoaded(false),fSpectrumRndm(0),fDBName("SOLAR") { // Initialize pointers #ifdef RATDEBUG debug << "[ESgen]::In constructor with gentype [" << fGenType << "]." << newline; #endif fMassElectron = electron_mass_c2; // Load the default generator // Later, depending on the options passed this can be reloaded // Data needed to Load the stuff: type of spectrum (solar or ibd, flux and flavor) if (fNuType.length()>0 && fNuFlavor.length()> 0) { #ifdef RATDEBUG detail << "[ESgen]::ESgen : Loading generator with nutype [" << fNuType << " and flavour " << fNuFlavor << "]." << newline; #endif LoadGenerator(); } #ifdef RATDEBUG debug << "[ESgen]::ESgen :Leaving constructor." << newline; #endif } void ESgen::LoadGenerator() { // Check if the generator is already loaded. // If it is, do nothing if (fGenLoaded) return; if (fNuType.length() ==0 || fNuFlavor.length() == 0) { #ifdef RATDEBUG debug << "[ESgen]::LoadGenerator : Loading failed because not all parameters set [" << fGenType << ","<< fNuType << "," << fNuFlavor << "]." << newline; #endif fGenLoaded = false; return; } // Check what type of ES generator we are dealing with // Depending on type the parameters and cross section pointers // are initialized differently // The parameters are taken from the database, depending of the option passed. // The original test with IBD data should still work DBLinkPtr linkdb; // Solar generator // The nu type is obtained from the job options (it defaults to pep) #ifdef RATDEBUG detail << "[ESgen]::LoadGenerator : Loading generator loaded." << newline; #endif linkdb = DB::Get()->GetLink(fDBName,fNuType); fEnuMin = linkdb->GetD("emin"); fEnuMax = linkdb->GetD("emax"); // fTotalFlux = linkdb->GetD("flux"); fEnuTbl = linkdb->GetDArray("spec_e"); fFluxTbl = linkdb->GetDArray("spec_flux"); // fNuSpectrum = new TGraph(fEnuTbl.size(),&fEnuTbl[0],&fFluxTbl[0]); // initialize the cross-section if (fXS != 0) { delete fXS; } fXS = new ESCrossSec(fNuFlavor); // To sample neutrino energy need to scale flux by total // cross section at that neutrino energy std::vector csScaledFluxTbl(fFluxTbl.size(),0); for (size_t i=0;iSigma(fEnuTbl[i]); } // If random sampler hasn't been initialized yet, lets do it now if (!fSpectrumRndm) { // Be7 is always a particular case due to its discrete nature // The last parameter is set to 1 to disallow interpolations if (fNuType == "be7") { fSpectrumRndm = new CLHEP::RandGeneral(&csScaledFluxTbl[0],fFluxTbl.size(),1); } else { // set interpolation bit to 0 to allow for interpolations in continuous // spectra fSpectrumRndm = new CLHEP::RandGeneral(&csScaledFluxTbl[0],fFluxTbl.size(),0); } } #ifdef RATDEBUG detail << "[ESgen]::LoadGenerator : ES generator loaded for [" << fDBName << "::" << fNuType << "] database input." << newline; #endif // If it reaches this point without failing then everything should be fine fGenLoaded = true; } ESgen::~ESgen() { #ifdef RATDEBUG detail << "[ESgen]:: In destructor." << newline; #endif if ( fXS != 0 ) { delete fXS; fXS = 0; } if ( fNuSpectrum != 0) { delete fNuSpectrum; fNuSpectrum = 0; } if (fSpectrumRndm) { delete fSpectrumRndm; fSpectrumRndm = 0; } } void ESgen::GenerateEvent(const G4ThreeVector& theNeutrino, G4LorentzVector& nu_incoming, G4LorentzVector &electron) { // Check if the generator has been loaded successfully // For now just throw something that can be caught at an upper level. // Need to define a set of specific exceptions if (!fGenLoaded) { G4Exception("[ESgen]::GenerateEvent","ArgError",FatalErrorInArgument,"Vertex generation called but it seems that it is not ready yet."); } ///! ///! Have to be careful with the line neutrino types (pep) ///! and even more careful with the double line (be7) ///! #ifdef RATDEBUG detail << "[ESgen]::GenerateEvent : Starting with direction ( " << theNeutrino.x() << " , " << theNeutrino.y() << " , " << theNeutrino.z() << " ) " << newline; detail << "[ESgen]::GenerateEvent : Current state is:" << newline; Show(); #endif // Throw values against a cross-section. //bool passed=false; double Enu, Te; // Updated sampler (orders of magnitude faster) // Given the neutrino energy, use the differential cross section // shape and sample from it. Enu = SampleNuEnergy()*MeV; Te = SampleRecoilEnergy(Enu)*MeV; // from the incoming neutrino we have already the initial direction. // The final electron direction will follow that. // Now we have : // - the initial direction of the neutrino (unnormalized) // - the energy of the neutrino // - the recoil energy of the electron // build the 4-momentum vector of the neutrino // We will only use the neutrino initial momentum as a baseline to add up to the electron direction #ifdef RATDEBUG detail << "[ESgen]::GenerateEvent : Neutrino incoming energy Enu=" << Enu/MeV << " MeV." << newline; detail << "[ESgen]::GenerateEvent : Electron recoil kinetic energy Te=" << Te/MeV << " MeV." << newline; #endif G4double theta_e = acos(sqrt((Te*(fMassElectron+Enu)*(fMassElectron+Enu))/(2*fMassElectron*Enu*Enu + Enu*Enu*Te))); G4double tot_Ee = Te + fMassElectron; G4double p_e = sqrt(tot_Ee*tot_Ee - fMassElectron*fMassElectron); #ifdef RATDEBUG detail << "[ESgen]::GenerateEvent : Electron recoil angle w.r.t. neutrino direction : " << theta_e << " rad ( " << TMath::RadToDeg()*theta_e << " degrees)." << newline; #endif // We have theta. Now randomly generate phi G4double phi_e = twopi*G4UniformRand(); // This is the incoming neutrino information nu_incoming.setVect(theNeutrino*Enu); nu_incoming.setE(Enu); // compute electron 4-momentum G4ThreeVector e_mom(p_e*theNeutrino); // Rotation from nu direction into electron direction G4ThreeVector rotation_axis = theNeutrino.orthogonal(); rotation_axis.rotate(phi_e,theNeutrino); e_mom.rotate(theta_e,rotation_axis); electron.setVect(e_mom); electron.setE(tot_Ee); #ifdef RATDEBUG detail << "[ESgen]::GenerateEvent : Electron energy-momentum ( " << electron.x() << " , " << electron.y() << " , " << electron.z() << " ) " << newline; #endif // TODO: If we want to keep track of the outgoing neutrino have to add it // here as well and store the information in a new argument G4LorentzVector variable // For the moment we only pass the incoming neutrino information back. } void ESgen::Reset() { // Reset the falg dependent objects. // After this method a call to LoadGenerator should always follow if (fNuSpectrum) { delete fNuSpectrum; fNuSpectrum = 0; } if (fSpectrumRndm) { delete fSpectrumRndm; fSpectrumRndm = 0; } LoadGenerator(); } void ESgen::Show() { detail << "Elastic Scattering Settings:" << newline; detail << "NuType : " << fNuType.c_str() << newline; } // // If we change the neutrino type we should reload the generator // to force it to reload the spectra from the database void ESgen::SetNuType(const G4String &nutype) { if (fGenType != nutype ) { #ifdef RATDEBUG detail << "[ESgen]::SetNuType : Switching to neutrino type [ " << nutype << " ]." << newline; #endif fNuType = nutype; fGenLoaded = false; LoadGenerator(); } } // // If we change the neutrino flavor we should reload the generator // to force it to reload the spectra from the database void ESgen::SetNuFlavor(const G4String &nuflavor) { if (fNuFlavor != nuflavor ) { #ifdef RATDEBUG detail << "[ESgen]::SetNuFlavor : Switching to neutrino flavor [ " << nuflavor << " ]." << newline; #endif fNuFlavor = nuflavor; fGenLoaded = false; LoadGenerator(); } } // This function samples the energy spectrum of the chosen neutrino and // decides from it the proper energy. // Keep in mind that pep is always the same, but be7 is a *very* special case G4double ESgen::SampleRecoilEnergy(G4double Enu) { G4double Te = 0.0; // Get the shape of the differential cross section. TGraph *dsigmadt = fXS->DrawdSigmadT(Enu); // Interpolate between the discrete points CLHEP::RandGeneral *rndm = new CLHEP::RandGeneral(dsigmadt->GetY(),dsigmadt->GetN(),0); Te = rndm->shoot()*(dsigmadt->GetX())[dsigmadt->GetN()-1]; #ifdef RATDEBUG debug << "[ESgen]::SampleRecoilEnergy : Random \'shoot\' for Enu = " << Enu << " returned a recoil of " << Te << newline; #endif delete rndm; delete dsigmadt; return Te; } // This function samples the energy spectrum of the chosen neutrino and // decides from it the proper energy. // Keep in mind that pep is always the same, but be7 is a *very* special case G4double ESgen::SampleNuEnergy() { G4double Enu = 0.0; G4double tmp = 0.0; if (fNuType == G4String("pep")) { // This is a line flux. Only 1 energy possible. Enu = fEnuMin; } else if ( fNuType == G4String("be7") ) { // The neutrino energy of be7 must follow the branching ratio. // Otherwise the flux is not properly sampled. // Use CLHEP RandGeneral generator, with interpolation disabled to build // a discrete distribution of two states. // The random generator will return either 0 or 0.5 tmp = fSpectrumRndm->shoot(); Enu = (tmp<0.5)?fEnuMin:fEnuMax; } else { // Continuous distributions will return a random number between 0 and 1 // Following the spectrum shape. // Scale it to the energy of the spectrum double scale = fEnuMax - fEnuMin; Enu = fEnuMin + fSpectrumRndm->shoot() * scale; } #ifdef RATDEBUG detail << "[ESgen]::SampleNuEnergy : Sampler resulting energy Enu=" << Enu << " MeV." << newline; #endif return Enu; } G4double ESgen::GetRatePerTarget() { if (!fGenLoaded) { G4Exception("ESgen::GetRatePerTarget","404",FatalErrorInArgument,"Trying to get rate before initializing all parameters."); } G4double intRate = 0.0; // Have to deal with the special cases separately if (fNuType == "pep") { intRate = fXS->Sigma(fEnuMin); } else if (fNuType == "be7") { intRate = fXS->Sigma(fEnuMin)*fNuSpectrum->GetY()[0]; intRate += fXS->Sigma(fEnuMax)*fNuSpectrum->GetY()[1]; } else { // For some unknown reason ROOT TGraph::Integral() does not // like to calculate integrals of adjacent points // so I do it by hand using a trapezoidal rule // The input spectra have a fine enough grain that the error // is negligible. for (G4int ip = 0; ip < fNuSpectrum->GetN()-1; ++ip) { G4double de = fNuSpectrum->GetX()[ip+1]-fNuSpectrum->GetX()[ip]; G4double integ = (0.5*de)*(fNuSpectrum->GetY()[ip+1]*fXS->Sigma((fNuSpectrum->GetX())[ip+1]) + fNuSpectrum->GetY()[ip]*fXS->Sigma((fNuSpectrum->GetX())[ip])); intRate += integ; } } #ifdef RATDEBUG detail << "[ESgen]::GetRatePerTarget : Total spectrum weighted cross section " << intRate*1e-42 << " cm^2 " << newline; detail << "[ESgen]::GetRatePerTarget : Total flux " << fTotalFlux << newline; detail << "[ESgen]::GetRatePerTarget : Final rate per target : " << intRate*1e-42*fTotalFlux << " events/target/s." << newline; #endif // don't forget the scale factor from the cross section // nor the total flux return intRate*1e-42*fTotalFlux; } } // namespace RAT