genie::rew::GReWeightFGM Class Reference

Reweighting the Fermi Gas nuclear model. More...

#include <GReWeightFGM.h>

Inheritance diagram for genie::rew::GReWeightFGM:

Public Member Functions
	GReWeightFGM ()
	~GReWeightFGM ()
bool	IsHandled (GSyst_t syst)
	does the current weight calculator handle the input nuisance param? More...
void	SetSystematic (GSyst_t syst, double val)
	update the value for the specified nuisance param More...
void	Reset (void)
	set all nuisance parameters to default values More...
void	Reconfigure (void)
	propagate updated nuisance parameter values to actual MC, etc More...
double	CalcWeight (const EventRecord &event)
	calculate a weight for the input event using the current nuisance param values More...
Public Member Functions inherited from genie::rew::GReWeightI
virtual	~GReWeightI ()

Private Member Functions
void	Init (void)
double	RewCCQEPauliSupViaKF (const EventRecord &event)
double	RewCCQEMomDistroFGtoSF (const EventRecord &event)

Private Attributes
double	fKFTwkDial
double	fMomDistroTwkDial
const NuclearModelI *	fFG
const NuclearModelI *	fSF
map< int, TH1D * >	fMapFGn
map< int, TH1D * >	fMapFGp
map< int, TH1D * >	fMapSFn
map< int, TH1D * >	fMapSFp

Additional Inherited Members
Protected Member Functions inherited from genie::rew::GReWeightI
	GReWeightI ()

Detailed Description

Reweighting the Fermi Gas nuclear model.

Author

Jim Dobson <J.Dobson07 imperial.ac.uk> Imperial College London

Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

Apr 26, 2010

Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 46 of file GReWeightFGM.h.

Constructor & Destructor Documentation

GReWeightFGM::GReWeightFGM

(

)

Definition at line 41 of file GReWeightFGM.cxx.

                           :
GReWeightI()
{
  this->Init();
}

GReWeightFGM::~GReWeightFGM

(

)

Definition at line 47 of file GReWeightFGM.cxx.

{
#ifdef _G_REWEIGHT_FGM_DEBUG_
  fTestFile->cd();
  fTestNtp ->Write();
  fTestFile->Close();
  delete fTestFile;
#endif
}

Member Function Documentation

double GReWeightFGM::CalcWeight ( const EventRecord &  event )

virtual

calculate a weight for the input event using the current nuisance param values

Implements genie::rew::GReWeightI.

Definition at line 97 of file GReWeightFGM.cxx.

{ 
  double wght =
    this->RewCCQEPauliSupViaKF   (event) *
    this->RewCCQEMomDistroFGtoSF (event);
 
  return wght;
}

void GReWeightFGM::Init ( void  )

private

Definition at line 251 of file GReWeightFGM.cxx.

{
  fKFTwkDial        = 0.;
  fMomDistroTwkDial = 0.;

  AlgFactory * algf = AlgFactory::Instance();

  fFG = dynamic_cast<const NuclearModelI*> (
    algf->GetAlgorithm("genie::FGMBodekRitchie","Default"));
  fSF = dynamic_cast<const NuclearModelI*> (
    algf->GetAlgorithm("genie::SpectralFunc","Default"));

#ifdef _G_REWEIGHT_FGM_DEBUG_
  fTestFile = new TFile("./fgm_reweight_test.root","recreate");
  fTestNtp  = new TNtupleD("testntp","","Q2:wght");
#endif
}

bool GReWeightFGM::IsHandled ( GSyst_t  syst )

virtual

does the current weight calculator handle the input nuisance param?

Implements genie::rew::GReWeightI.

Definition at line 57 of file GReWeightFGM.cxx.

{
  switch(syst) {
    case ( kSystNucl_CCQEPauliSupViaKF   ) : 
    case ( kSystNucl_CCQEMomDistroFGtoSF ) : 
        return true;
        break;
    default:
        return false;
        break;
  }
  return false;
}

void GReWeightFGM::Reconfigure ( void  )

virtual

propagate updated nuisance parameter values to actual MC, etc

Implements genie::rew::GReWeightI.

Definition at line 92 of file GReWeightFGM.cxx.

{
  
}

void GReWeightFGM::Reset ( void  )

virtual

set all nuisance parameters to default values

Implements genie::rew::GReWeightI.

Definition at line 86 of file GReWeightFGM.cxx.

{
  fKFTwkDial        = 0.;
  fMomDistroTwkDial = 0.;
}

double GReWeightFGM::RewCCQEMomDistroFGtoSF ( const EventRecord &  event )

private

Definition at line 174 of file GReWeightFGM.cxx.

{
  bool momdistro_tweaked = (TMath::Abs(fMomDistroTwkDial) > controls::kASmallNum);
  if(!momdistro_tweaked) return 1.;

  bool is_qe = event.Summary()->ProcInfo().IsQuasiElastic();
  bool is_cc = event.Summary()->ProcInfo().IsWeakCC();
  if(!is_qe || !is_cc) return 1.;

  GHepParticle * tgtnucleus = event.TargetNucleus();
  if(!tgtnucleus) return 1.; // scattering off free-nucleon 

  GHepParticle * hitnucleon = event.HitNucleon();
  if(!hitnucleon) return 1.;

  const double kPmax = 0.5;
  double p = hitnucleon->P4()->Vect().Mag();
  if(p > kPmax) return 1.;

  TH1D * hfg = 0;
  TH1D * hsf = 0;

  int tgtpdg = tgtnucleus -> Pdg();
  int nucpdg = hitnucleon -> Pdg();

  map<int, TH1D*> & mapfg = pdg::IsNeutron(nucpdg) ? fMapFGn : fMapFGp;
  map<int, TH1D*> & mapsf = pdg::IsNeutron(nucpdg) ? fMapSFn : fMapSFp;

  map<int, TH1D*>::iterator it;
  it = mapfg.find(tgtpdg);
  if(it != mapfg.end()) { hfg = it->second; }
  it = mapsf.find(tgtpdg);
  if(it != mapsf.end()) { hsf = it->second; }

  bool have_weight_func = (hfg!=0) && (hsf!=0);
  if(!have_weight_func) {
     const int kNEv  = 20000;
     const int kNP   = 500;
     hfg = new TH1D("","",kNP,0.,kPmax);
     hsf = new TH1D("","",kNP,0.,kPmax);
     hfg -> SetDirectory(0);
     hsf -> SetDirectory(0);
     const Target & tgt = event.Summary()->InitState().Tgt();
     bool ok = true;
     for(int iev=0; iev<kNEv; iev++) {
       ok = fFG->GenerateNucleon(tgt);
       if(!ok) {
         delete hfg;
         delete hsf;
         return 1.;
       }
       hfg->Fill(fFG->Momentum());
     }//fg
     for(int iev=0; iev<kNEv; iev++) {
       ok = fSF->GenerateNucleon(tgt);
       if(!ok) {
         delete hfg;
         delete hsf;
         return 1.;
       }
       hsf->Fill(fSF->Momentum());
     }//sf
     hfg->Scale(1. / hfg->Integral("width"));
     hsf->Scale(1. / hsf->Integral("width"));
     mapfg.insert(map<int,TH1D*>::value_type(tgtpdg,hfg));
     mapsf.insert(map<int,TH1D*>::value_type(tgtpdg,hsf));
  }//create & store momentum distributions


  double f_fg = hfg->GetBinContent( hfg->FindBin(p) );
  double f_sf = hsf->GetBinContent( hsf->FindBin(p) );
  double dial = fMomDistroTwkDial;
  double wght = (f_sf * dial + f_fg * (1-dial)) / f_fg;

  return wght;
}

double GReWeightFGM::RewCCQEPauliSupViaKF ( const EventRecord &  event )

private

Definition at line 106 of file GReWeightFGM.cxx.

{
  bool kF_tweaked = (TMath::Abs(fKFTwkDial) > controls::kASmallNum);
  if(!kF_tweaked) return 1.;

  bool is_qe = event.Summary()->ProcInfo().IsQuasiElastic();
  bool is_cc = event.Summary()->ProcInfo().IsWeakCC();
  if(!is_qe || !is_cc) return 1.;

  const Target & target = event.Summary()->InitState().Tgt();
  if (!target.IsNucleus()) {
     return 1.; 
  }

  double A = target.A();
  if(A <= 4) {
     return 1.; 
  }

  int target_pdgc         = target.Pdg();
  int struck_nucleon_pdgc = target.HitNucPdg();
  int final_nucleon_pdgc  = pdg::SwitchProtonNeutron(struck_nucleon_pdgc);
             
  FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
  const FermiMomentumTable * kft  = kftp->GetTable("Default");

  // Fermi momentum for initial, final nucleons
  double kFi = kft->FindClosestKF(target_pdgc, struck_nucleon_pdgc);
  double kFf = kft->FindClosestKF(target_pdgc, final_nucleon_pdgc );
 
  // hit nucleon mass
  double Mn = target.HitNucP4Ptr()->M(); // can be off m/shell

  // momentum transfer
  const Kinematics & kine = event.Summary()->Kine();
  bool selected = true;
  double q2 = kine.q2(selected);
  
  const double pmax = 0.5; // GeV

  // default nuclear suppression 
  double R = utils::nuclear::RQEFG_generic(q2, Mn, kFi, kFf, pmax);

  // tweak kF
  GSystUncertainty * uncertainty = GSystUncertainty::Instance();
  double kF_fracerr = uncertainty->OneSigmaErr(kSystNucl_CCQEPauliSupViaKF);

  double kFi_twk = kFi * (1 + fKFTwkDial * kF_fracerr);
  double kFf_twk = kFf * (1 + fKFTwkDial * kF_fracerr);

  // calculate tweaked nuclear suppression factor
  double Rtwk = nuclear::RQEFG_generic(q2, Mn, kFi_twk, kFf_twk, pmax);
  
  // calculate weight (ratio of suppression factors)

  double wght = 1.0;

  if(R>0 && Rtwk>0) {
    wght = Rtwk/R;
  }

#ifdef _G_REWEIGHT_FGM_DEBUG_
  fTestNtp->Fill(-q2,wght);
#endif

  return wght;
}

void GReWeightFGM::SetSystematic ( GSyst_t  syst, double  val  )

virtual

update the value for the specified nuisance param

Implements genie::rew::GReWeightI.

Definition at line 71 of file GReWeightFGM.cxx.

{
  switch(syst) {
    case ( kSystNucl_CCQEPauliSupViaKF   ) : 
        fKFTwkDial = val;
        break;
    case ( kSystNucl_CCQEMomDistroFGtoSF ) : 
        fMomDistroTwkDial = val;
        break;
    default:
        return;
        break;
  }
}

Member Data Documentation

const NuclearModelI* genie::rew::GReWeightFGM::fFG

private

Definition at line 69 of file GReWeightFGM.h.

double genie::rew::GReWeightFGM::fKFTwkDial

private

Definition at line 66 of file GReWeightFGM.h.

map<int, TH1D *> genie::rew::GReWeightFGM::fMapFGn

private

Definition at line 72 of file GReWeightFGM.h.

map<int, TH1D *> genie::rew::GReWeightFGM::fMapFGp

private

Definition at line 73 of file GReWeightFGM.h.

map<int, TH1D *> genie::rew::GReWeightFGM::fMapSFn

private

Definition at line 74 of file GReWeightFGM.h.

map<int, TH1D *> genie::rew::GReWeightFGM::fMapSFp

private

Definition at line 75 of file GReWeightFGM.h.

double genie::rew::GReWeightFGM::fMomDistroTwkDial

private

Definition at line 67 of file GReWeightFGM.h.

const NuclearModelI* genie::rew::GReWeightFGM::fSF

private

Definition at line 70 of file GReWeightFGM.h.

---

The documentation for this class was generated from the following files:

• GENIE/src/ReWeight/GReWeightFGM.h
• GENIE/src/ReWeight/GReWeightFGM.cxx