Event reweighting utilities. More...

Functions
double	MeanFreePathWeight (int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double mfp_scale_factor, bool interacted, double nRpi=0.5, double nRnuc=1.0, double NR=3, double R0=1.4)

double	MeanFreePathWeight (double prob_def, double prob_twk, bool interacted)

double	FZoneWeight (int pdgc, const TLorentzVector &vtx, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double fz_scale_factor, bool interacted, double nRpi=0.5, double nRnuc=1.0, double NR=3, double R0=1.4)

double	FateFraction (genie::rew::GSyst_t syst, double kinE, double frac_scale_factor=1.)

double	WhichFateFractionScaleFactor (genie::rew::GSyst_t syst, double kinE, double fate_frac)

bool	HadronizedByAGKY (const EventRecord &event)

bool	HadronizedByAGKYPythia (const EventRecord &event)

TLorentzVector	Hadronic4pLAB (const EventRecord &event)

double	AGKYWeight (int pdgc, double xF, double pT2)

int	Sign (double twkdial)

Detailed Description

Event reweighting utilities.

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

Sep 09, 2009

Function Documentation

double genie::utils::rew::AGKYWeight	(	int	pdgc,
		double	xF,
		double	pT2
	)

Definition at line 299 of file GReWeightUtils.cxx.

 {
   return 1.0;
 }

double genie::utils::rew::FateFraction	(	genie::rew::GSyst_t	syst,
		double	kinE,
		double	frac_scale_factor = `1.`
	)

Definition at line 144 of file GReWeightUtils.cxx.

 {
   double fate_frac = 0.0;
       
   INukeHadroData * hd = INukeHadroData::Instance();
 
   // convert to MeV and      
   double ke = kinE / units::MeV;         
   ke = TMath::Max(INukeHadroData::fMinKinEnergy,   ke);
   ke = TMath::Min(INukeHadroData::fMaxKinEnergyHA, ke);
          
   switch (syst) {
 
     //
     // pions
     //
 
     case (genie::rew::kINukeTwkDial_FrCEx_pi) :
     {
       fate_frac = hd->Frac(kPdgPiP, kIHAFtCEx, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrElas_pi) :
     {
       fate_frac = hd->Frac(kPdgPiP, kIHAFtElas, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrInel_pi) :
     {
       fate_frac = hd->Frac(kPdgPiP, kIHAFtInelas, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrAbs_pi) :
     {  
       fate_frac = hd->Frac(kPdgPiP, kIHAFtAbs, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrPiProd_pi) :
     {
       fate_frac = hd->Frac(kPdgPiP, kIHAFtPiProd,  ke);
     }
     break;
  
     //
     // nucleons
     //
 
     case (genie::rew::kINukeTwkDial_FrCEx_N) :
     {
       fate_frac = hd->Frac(kPdgProton, kIHAFtCEx, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrElas_N) :
     {
       fate_frac = hd->Frac(kPdgProton, kIHAFtElas, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrInel_N) :
     {
       fate_frac = hd->Frac(kPdgProton, kIHAFtInelas, ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrAbs_N) :
     {
       fate_frac = hd->Frac(kPdgProton, kIHAFtAbs,    ke);
     }
     break;
 
     case (genie::rew::kINukeTwkDial_FrPiProd_N) :
     {
       fate_frac = hd->Frac(kPdgProton, kIHAFtPiProd,  ke);
     }   
     break;
 
     default:
     {
       LOG("ReW", pDEBUG)
         << "Have reached default case and assigning fraction{fate} = 0";
       fate_frac = 0;
     }
     break;
 
   } // hadron_fate?
 
   fate_frac *= frac_scale_factor;
                  
   return fate_frac;
 }

double genie::utils::rew::FZoneWeight	(	int	pdgc,
		const TLorentzVector &	vtx,
		const TLorentzVector &	x4,
		const TLorentzVector &	p4,
		double	A,
		double	Z,
		double	fz_scale_factor,
		bool	interacted,
		double	nRpi = `0.5`,
		double	nRnuc = `1.0`,
		double	NR = `3`,
		double	R0 = `1.4`
	)

Definition at line 81 of file GReWeightUtils.cxx.

 {
    // Calculate hadron start assuming tweaked formation zone
    TLorentzVector fz    = x4 - vtx;
    TLorentzVector fztwk = fz_scale_factor*fz;
    TLorentzVector x4twk = x4 + fztwk - fz;
 
    LOG("ReW", pDEBUG)  << "Formation zone = "<< fz.Vect().Mag() << " fm";
 
    // Get nominal survival probability.
    double pdef = utils::intranuke::ProbSurvival(
       pdgc,x4,p4,A,Z,1.,nRpi,nRnuc,NR,R0);
    LOG("ReW", pDEBUG)  << "Survival probability (nominal) = "<< pdef;
    if(pdef<=0) return 1.; 
    if(pdef>=1.){
      LOG("ReW", pERROR) 
          <<  "Default formation zone takes hadron outside "
          <<  "nucleus so cannot reweight!" ;
      return 1.;
    }
 
    // Get tweaked survival probability.
    double ptwk = utils::intranuke::ProbSurvival(
       pdgc,x4twk,p4,A,Z,1.,nRpi,nRnuc,NR,R0);
    if(ptwk<=0) return 1.;
    LOG("ReW", pDEBUG)  << "Survival probability (tweaked) = "<< ptwk;
 
    // Calculate weight
    double w_fz = utils::rew::MeanFreePathWeight(pdef, ptwk, interacted);
    LOG("ReW", pDEBUG)  
        << "Particle weight for formation zone tweak = "<< ptwk;
    return w_fz;  
 }

TLorentzVector genie::utils::rew::Hadronic4pLAB ( const EventRecord & event )

Definition at line 278 of file GReWeightUtils.cxx.

 {
   GHepParticle * nu = event.Probe();                    // incoming v
   GHepParticle * N  = event.HitNucleon();               // struck nucleon
   GHepParticle * l  = event.FinalStatePrimaryLepton();  // f/s primary lepton
  
   assert(nu);                                                              
   assert(N);
   assert(l);
 
   // Compute the Final State Hadronic System 4p (PX = Pv + PN - Pl)
 
   const TLorentzVector & p4nu = *(nu->P4());
   const TLorentzVector & p4N  = *(N ->P4());
   const TLorentzVector & p4l  = *(l ->P4());
 
   TLorentzVector pX4 = p4nu + p4N - p4l;
 
   return pX4;
 }

bool genie::utils::rew::HadronizedByAGKY ( const EventRecord & event )

Definition at line 254 of file GReWeightUtils.cxx.

 { 
   Interaction * interaction = event.Summary();
   assert(interaction);                                                     
 
   bool is_dis  = interaction->ProcInfo().IsDeepInelastic();
   bool charm   = interaction->ExclTag().IsCharmEvent();
   bool by_agky = is_dis && !charm;
 
   return by_agky;
 }

bool genie::utils::rew::HadronizedByAGKYPythia ( const EventRecord & event )

Definition at line 266 of file GReWeightUtils.cxx.

 { 
 // Check whether the event was hadronized by AGKY/KNO or AGKY/PYTHIA
    
   GHepStatus_t prefragm = kIStDISPreFragmHadronicState;
   bool found_string  = (event.FindParticle(kPdgString,  prefragm, 0) != 0);
   bool found_cluster = (event.FindParticle(kPdgCluster, prefragm, 0) != 0);
   bool handled_by_pythia = found_string || found_cluster;
 
   return handled_by_pythia;
 }

double genie::utils::rew::MeanFreePathWeight	(	int	pdgc,
		const TLorentzVector &	x4,
		const TLorentzVector &	p4,
		double	A,
		double	Z,
		double	mfp_scale_factor,
		bool	interacted,
		double	nRpi = `0.5`,
		double	nRnuc = `1.0`,
		double	NR = `3`,
		double	R0 = `1.4`
	)

Definition at line 49 of file GReWeightUtils.cxx.

 {
    LOG("ReW", pINFO)  
      << "Calculating mean free path weight: "
      << "A = " << A << ", Z = " << Z << ", mfp_scale = " << mfp_scale_factor
      << ", interacted = " << interacted;
    LOG("ReW", pDEBUG)  
      << "nR_pion = " << nRpi << ", nR_nucleon = " << nRnuc 
      << ", NR = " << NR << ", R0 = " << R0;
 
    // Get the nominal survival probability
    double pdef = utils::intranuke::ProbSurvival(
       pdgc,x4,p4,A,Z,1.,nRpi,nRnuc,NR,R0);
    LOG("ReW", pINFO)  << "Probability(default mfp) = " << pdef;      
    if(pdef<=0) return 1.;
 
    // Get the survival probability for the tweaked mean free path
    double ptwk = utils::intranuke::ProbSurvival(
       pdgc,x4,p4,A,Z,mfp_scale_factor,nRpi,nRnuc,NR,R0);
    LOG("ReW", pINFO)  << "Probability(tweaked mfp) = " << ptwk;      
    if(ptwk<=0) return 1.;
 
    // Calculate weight
    double w_mfp = utils::rew::MeanFreePathWeight(pdef, ptwk, interacted);
    LOG("ReW", pINFO)  << "Mean free path weight = " << w_mfp;     
    return w_mfp;
 }

double genie::utils::rew::MeanFreePathWeight	(	double	prob_def,
		double	prob_twk,
		bool	interacted
	)

Definition at line 119 of file GReWeightUtils.cxx.

 {
 // Returns a weight to account for a change in hadron mean free path inside
 // insidea nuclear medium.
 //
 // Inputs:
 //   pdef : nominal survival probability
 //   ptwk : survival probability for the tweaked mean free path
 //   interacted : flag indicating whether the hadron interacted or escaped
 //
 // See utils::intranuke::ProbSurvival() for the calculation of probabilities.
 //
   double w_mfp = 1.;
 
   if(interacted) {
      w_mfp = (1-pdef>0) ?  (1-ptwk)  /   (1-pdef)  : 1;
   } else {
      w_mfp = (pdef>0) ? ptwk / pdef : 1;
   }
   w_mfp = TMath::Max(0.,w_mfp);
 
   return w_mfp;
 }

int genie::utils::rew::Sign ( double twkdial )

Definition at line 304 of file GReWeightUtils.cxx.

 {
   if(twkdial < 0.) return -1;
   if(twkdial > 0.) return +1;
   return 0;
 }

double genie::utils::rew::WhichFateFractionScaleFactor	(	genie::rew::GSyst_t	syst,
		double	kinE,
		double	fate_frac
	)

Definition at line 241 of file GReWeightUtils.cxx.

 {  
   double fate_frac_nominal = FateFraction(syst,kinE,1.);
 
   if(TMath::Abs(fate_frac-fate_frac_nominal) < kASmallNum) return 0;
 
   if(fate_frac_nominal <= 0) { return -99999; }
 
   double scale = TMath::Max(0.,fate_frac)/fate_frac_nominal;
   return scale;
 }

Functions

Detailed Description

Function Documentation