Reweighting resonance decays. More...

#include <GReWeightResonanceDecay.h>

Inheritance diagram for genie::rew::GReWeightResonanceDecay:

Public Member Functions
	GReWeightResonanceDecay ()

	~GReWeightResonanceDecay ()

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...

void	RewNue (bool tf)

void	RewNuebar (bool tf)

void	RewNumu (bool tf)

void	RewNumubar (bool tf)

void	RewCC (bool tf)

void	RewNC (bool tf)

Public Member Functions inherited from genie::rew::GReWeightI
virtual	~GReWeightI ()

Private Member Functions
void	Init (void)

double	RewBR (const EventRecord &event)

double	RewThetaDelta2Npi (const EventRecord &event)

Private Attributes
double	fBR1gammaTwkDial

double	fBR1etaTwkDial

double	fThetaDelta2NpiTwkDial

bool	fRewNue
	reweight nu_e? More...

bool	fRewNuebar
	reweight nu_e_bar? More...

bool	fRewNumu
	reweight nu_mu? More...

bool	fRewNumubar
	reweight nu_mu_bar? More...

bool	fRewCC
	reweight CC? More...

bool	fRewNC
	reweight NC? More...

map< int, TH1D * >	fMpBR1gammaDef

map< int, TH1D * >	fMpBR1etaDef

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

Detailed Description

Reweighting resonance decays.

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

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

Apr 26, 2010

Definition at line 43 of file GReWeightResonanceDecay.h.

Constructor & Destructor Documentation

GReWeightResonanceDecay::GReWeightResonanceDecay ( )

Definition at line 39 of file GReWeightResonanceDecay.cxx.

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

GReWeightResonanceDecay::~GReWeightResonanceDecay ( )

Definition at line 45 of file GReWeightResonanceDecay.cxx.

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

Member Function Documentation

double GReWeightResonanceDecay::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 100 of file GReWeightResonanceDecay.cxx.

 {  
   Interaction * interaction = event.Summary();
 
   bool is_res = interaction->ProcInfo().IsResonant();
   if(!is_res) return 1.;
 
   bool is_cc  = interaction->ProcInfo().IsWeakCC();
   bool is_nc  = interaction->ProcInfo().IsWeakNC();
   if(is_cc && !fRewCC) return 1.;
   if(is_nc && !fRewNC) return 1.;
 
   int nupdg = interaction->InitState().ProbePdg();
   if(nupdg==kPdgNuMu     && !fRewNumu   ) return 1.;
   if(nupdg==kPdgAntiNuMu && !fRewNumubar) return 1.;
   if(nupdg==kPdgNuE      && !fRewNue    ) return 1.;
   if(nupdg==kPdgAntiNuE  && !fRewNuebar ) return 1.;
 
   double wght = 
     this->RewBR(event) * 
     this->RewThetaDelta2Npi(event);
 
   return wght;
 }

void GReWeightResonanceDecay::Init ( void )

private

Definition at line 355 of file GReWeightResonanceDecay.cxx.

 {
   this->RewNue    (true);
   this->RewNuebar (true);
   this->RewNumu   (true);
   this->RewNumubar(true);
   this->RewCC     (true);
   this->RewNC     (true);
 
   fBR1gammaTwkDial       = 0.;
   fBR1etaTwkDial         = 0.;
   fThetaDelta2NpiTwkDial = 0.;
 
   const int respdgarray[] = {
     kPdgP33m1232_DeltaM, kPdgP33m1232_Delta0, kPdgP33m1232_DeltaP, kPdgP33m1232_DeltaPP,
     kPdgS31m1620_DeltaM, kPdgS31m1620_Delta0, kPdgS31m1620_DeltaP, kPdgS31m1620_DeltaPP,
     kPdgD33m1700_DeltaM, kPdgD33m1700_Delta0, kPdgD33m1700_DeltaP, kPdgD33m1700_DeltaPP,
     kPdgF35m1905_DeltaM, kPdgF35m1905_Delta0, kPdgF35m1905_DeltaP, kPdgF35m1905_DeltaPP,
     kPdgP31m1910_DeltaM, kPdgP31m1910_Delta0, kPdgP31m1910_DeltaP, kPdgP31m1910_DeltaPP,
     kPdgP33m1920_DeltaM, kPdgP33m1920_Delta0, kPdgP33m1920_DeltaP, kPdgP33m1920_DeltaPP,
     kPdgF37m1950_DeltaM, kPdgF37m1950_Delta0, kPdgF37m1950_DeltaP, kPdgF37m1950_DeltaPP,
     kPdgP11m1440_N0, kPdgP11m1440_NP,    
     kPdgD13m1520_N0, kPdgD13m1520_NP,   
     kPdgS11m1535_N0, kPdgS11m1535_NP,  
     kPdgS11m1650_N0, kPdgS11m1650_NP,  
     kPdgD15m1675_N0, kPdgD15m1675_NP,    
     kPdgF15m1680_N0, kPdgF15m1680_NP,     
     kPdgD13m1700_N0, kPdgD13m1700_NP,  
     kPdgP11m1710_N0, kPdgP11m1710_NP,      
     kPdgP13m1720_N0, kPdgP13m1720_NP,
     0
   };
 
   // init
   const int    kNW   = 30;
   const double kWmin = 0.9;
   const double kWmax = 5.0;
   unsigned int ires=0;
   int respdg = 0;
   while((respdg = respdgarray[ires++])) {
     fMpBR1gammaDef [respdg] = new TH1D("","",kNW,kWmin,kWmax);
     fMpBR1etaDef   [respdg] = new TH1D("","",kNW,kWmin,kWmax);
     fMpBR1gammaDef [respdg] -> SetDirectory(0);
     fMpBR1etaDef   [respdg] -> SetDirectory(0);
   }
 
   // find corresponding decay channels and store default BR
   PDGLibrary * pdglib = PDGLibrary::Instance();
   ires=0;
   while((respdg = respdgarray[ires++])) {
     TParticlePDG * res = pdglib->Find(respdg);
     if(!res) continue;
     for(int j=0; j<res->NDecayChannels(); j++) {
         TDecayChannel * dch = res->DecayChannel(j);
         int ngamma = 0;
         int neta   = 0;
         double br = dch->BranchingRatio();
         double mt = 0.;
         for(int k=0; k<dch->NDaughters(); k++) {
           int dpdg = dch->DaughterPdgCode(k);
           if(dpdg==kPdgGamma) ngamma++; 
           if(dpdg==kPdgEta  ) neta++; 
           mt += 0;
         }//decay channel f/s particles
         bool is_1gamma = (ngamma==1);
         bool is_1eta   = (neta  ==1);
         for(int ibin = 1; ibin <= fMpBR1gammaDef[respdg]->GetNbinsX(); ibin++) {
           double W = fMpBR1gammaDef[respdg]->GetBinLowEdge(ibin) + 
                      fMpBR1gammaDef[respdg]->GetBinWidth(ibin);
           bool is_allowed = (W>mt);
           if(is_allowed && is_1gamma) { fMpBR1gammaDef[respdg]->Fill(W, br); }
           if(is_allowed && is_1eta  ) { fMpBR1etaDef  [respdg]->Fill(W, br); }
         }//W bins
     }//decay channels
   }//resonances
 
 #ifdef _G_REWEIGHT_RESDEC_DEBUG_
   fTestFile = new TFile("./resdec_reweight_test.root","recreate");
   fTestNtp  = new TNtupleD("testntp","","costheta:wght");
 #endif
 }

bool GReWeightResonanceDecay::IsHandled ( GSyst_t syst )

virtual

does the current weight calculator handle the input nuisance param?

Implements genie::rew::GReWeightI.

Definition at line 55 of file GReWeightResonanceDecay.cxx.

 {
   switch(syst) {
      case ( kRDcyTwkDial_BR1gamma         ) :
      case ( kRDcyTwkDial_BR1eta           ) :
      case ( kRDcyTwkDial_Theta_Delta2Npi  ) :
         return true;
         break;
      default:
         return false;
         break;
   }
   return false;
 }

void GReWeightResonanceDecay::Reconfigure ( void )

virtual

propagate updated nuisance parameter values to actual MC, etc

Implements genie::rew::GReWeightI.

Definition at line 95 of file GReWeightResonanceDecay.cxx.

 {
 
 }

void GReWeightResonanceDecay::Reset ( void )

virtual

set all nuisance parameters to default values

Implements genie::rew::GReWeightI.

Definition at line 88 of file GReWeightResonanceDecay.cxx.

 {
   fBR1gammaTwkDial       = 0.;
   fBR1etaTwkDial         = 0.;
   fThetaDelta2NpiTwkDial = 0.;
 }

double GReWeightResonanceDecay::RewBR ( const EventRecord & event )

private

Definition at line 125 of file GReWeightResonanceDecay.cxx.

 {
   bool tweaked_br1gamma = (TMath::Abs(fBR1gammaTwkDial) > controls::kASmallNum);
   bool tweaked_br1eta   = (TMath::Abs(fBR1etaTwkDial)   > controls::kASmallNum);
 
   bool tweaked = (tweaked_br1gamma || tweaked_br1eta);
   if(!tweaked) return 1.;
 
   double wght = 1.;
 
   GSystUncertainty * uncertainty = GSystUncertainty::Instance();
 
   LOG("ReW", pDEBUG) << "Checking resonance decay mode.";
 
   GHepParticle * p = 0;
   TIter iter(&event);
   while((p=(GHepParticle*)iter.Next())) {
     if(utils::res::IsBaryonResonance(p->Pdg())) {
       int fd = p->FirstDaughter();
       int ld = p->LastDaughter();
       int ngamma = 0;
       int neta   = 0;
       for(int i=fd; i<=ld; i++) {
         int dpdg = event.Particle(i)->Pdg();
         if(dpdg==kPdgGamma) ngamma++; 
         if(dpdg==kPdgEta  ) neta++; 
       }//i
       bool is_1gamma = (ngamma==1);
       bool is_1eta   = (neta  ==1);
 
       if(is_1gamma) {
          LOG("ReW", pDEBUG) << "A resonance -> X + 1gamma event";
       }
       if(is_1eta) {
          LOG("ReW", pDEBUG) << "A resonance -> X + 1eta event";
       }
 
       //
       // For Resonance -> X + gamma, tweak the branching ratio as:
       // BR_{tweaked} = BR_{default} * (1 + dial * fractional_error)
       // so that, basically, weight = 1 + dial * fractional_error.
       // For all other decay modes adjust the branching ratio so that Sum{BR} = 1.
       //
       if(tweaked_br1gamma) {
         double frerr = uncertainty->OneSigmaErr(kRDcyTwkDial_BR1gamma);
         double dial  = fBR1gammaTwkDial;
         double w = (1. + dial*frerr);
         w = TMath::Max(0.,w);
         TH1D * brfw  = fMpBR1gammaDef[p->Pdg()];
         double mass = p->P4()->M();
         mass = TMath::Min(mass, brfw->GetXaxis()->GetXmax());
         int ibin = brfw->FindBin(mass);
         double brdef = brfw->GetBinContent(ibin);
         double brtwk = brdef*w;
         if(brtwk>1) {
          brtwk = 1.;
          w = brtwk/brdef;
         }
         if(is_1gamma) {
          wght *= w;
         } else {         
          wght *= ((1-brtwk)/(1-brdef));
         }
       }
       // Similarly for Resonance -> X + eta
       //
       if(tweaked_br1eta) {
         double frerr = uncertainty->OneSigmaErr(kRDcyTwkDial_BR1eta);
         double dial  = fBR1etaTwkDial;
         double w = (1. + dial*frerr);
         w = TMath::Max(0.,w);
         TH1D * brfw  = fMpBR1etaDef[p->Pdg()];
         double mass = p->P4()->M();
         mass = TMath::Min(mass, brfw->GetXaxis()->GetXmax());
         int ibin = brfw->FindBin(mass);
         double brdef = brfw->GetBinContent(ibin);
         double brtwk = brdef*w;
         if(brtwk>1) {
          brtwk = 1.;
          w = brtwk/brdef;
         }
         if(is_1eta) {
          wght *= w;
         } else {         
          wght *= ((1-brtwk)/(1-brdef));
         }
       }
       // Similarly for other modes with tweaked BRs
       //
       //...
       //...
 
     }//res?
   }//p
 
   return wght;
 }

void genie::rew::GReWeightResonanceDecay::RewCC ( bool tf )

inline

Definition at line 61 of file GReWeightResonanceDecay.h.

61 { fRewCC = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewCC

bool fRewCC

reweight CC?

Definition: GReWeightResonanceDecay.h:78

void genie::rew::GReWeightResonanceDecay::RewNC ( bool tf )

inline

Definition at line 62 of file GReWeightResonanceDecay.h.

62 { fRewNC = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewNC

bool fRewNC

reweight NC?

Definition: GReWeightResonanceDecay.h:79

void genie::rew::GReWeightResonanceDecay::RewNue ( bool tf )

inline

Definition at line 57 of file GReWeightResonanceDecay.h.

57 { fRewNue = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewNue

bool fRewNue

reweight nu_e?

Definition: GReWeightResonanceDecay.h:74

void genie::rew::GReWeightResonanceDecay::RewNuebar ( bool tf )

inline

Definition at line 58 of file GReWeightResonanceDecay.h.

58 { fRewNuebar = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewNuebar

bool fRewNuebar

reweight nu_e_bar?

Definition: GReWeightResonanceDecay.h:75

void genie::rew::GReWeightResonanceDecay::RewNumu ( bool tf )

inline

Definition at line 59 of file GReWeightResonanceDecay.h.

59 { fRewNumu = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewNumu

bool fRewNumu

reweight nu_mu?

Definition: GReWeightResonanceDecay.h:76

void genie::rew::GReWeightResonanceDecay::RewNumubar ( bool tf )

inline

Definition at line 60 of file GReWeightResonanceDecay.h.

60 { fRewNumubar = tf; }

tf

Definition: tf_graph.h:23

genie::rew::GReWeightResonanceDecay::fRewNumubar

bool fRewNumubar

reweight nu_mu_bar?

Definition: GReWeightResonanceDecay.h:77

double GReWeightResonanceDecay::RewThetaDelta2Npi ( const EventRecord & event )

private

Definition at line 223 of file GReWeightResonanceDecay.cxx.

 {
   bool tweaked = (TMath::Abs(fThetaDelta2NpiTwkDial) > controls::kASmallNum);
   if(!tweaked) return 1.;
 
   bool is_Delta_1pi = false;
   int ir = -1; // resonance position
   int ip = -1; // pion position
   int i  =  0;
   GHepParticle * p = 0;
   TIter iter(&event);
   while((p=(GHepParticle*)iter.Next())) {
     bool is_Deltapp = (p->Pdg()==kPdgP33m1232_DeltaPP);
     if(is_Deltapp) {
       ir = i;
       int fd = p->FirstDaughter();
       int ld = p->LastDaughter();
       int nd = 1 + ld - fd;
       if(nd==2) {
         int fpdg = event.Particle(fd)->Pdg();
         int lpdg = event.Particle(ld)->Pdg();
         if(fpdg==kPdgProton && lpdg==kPdgPiP) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgProton && fpdg==kPdgPiP) { is_Delta_1pi = true; ip = fd; }
       }
     }
 
     bool is_Deltap = (p->Pdg()==kPdgP33m1232_DeltaP);
     if(is_Deltap) {
       ir = i;
       int fd = p->FirstDaughter();
       int ld = p->LastDaughter();
       int nd = 1 + ld - fd;
       if(nd==2) {
         int fpdg = event.Particle(fd)->Pdg();
         int lpdg = event.Particle(ld)->Pdg();
         if(fpdg==kPdgProton && lpdg==kPdgPi0) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgProton && fpdg==kPdgPi0) { is_Delta_1pi = true; ip = fd; }
         if(fpdg==kPdgNeutron && lpdg==kPdgPiP) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgNeutron && fpdg==kPdgPiP) { is_Delta_1pi = true; ip = fd; }
       }
     }
 
     bool is_Delta0 = (p->Pdg()==kPdgP33m1232_Delta0);
     if(is_Delta0) {
       ir = i;
       int fd = p->FirstDaughter();
       int ld = p->LastDaughter();
       int nd = 1 + ld - fd;
       if(nd==2) {
         int fpdg = event.Particle(fd)->Pdg();
         int lpdg = event.Particle(ld)->Pdg();
         if(fpdg==kPdgProton && lpdg==kPdgPiM) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgProton && fpdg==kPdgPiM) { is_Delta_1pi = true; ip = fd; }
         if(fpdg==kPdgNeutron && lpdg==kPdgPi0) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgNeutron && fpdg==kPdgPi0) { is_Delta_1pi = true; ip = fd; }
       }
     }
 
     bool is_DeltaM = (p->Pdg()==kPdgP33m1232_DeltaM);
     if(is_DeltaM) {
       ir = i;
       int fd = p->FirstDaughter();
       int ld = p->LastDaughter();
       int nd = 1 + ld - fd;
       if(nd==2) {
         int fpdg = event.Particle(fd)->Pdg();
         int lpdg = event.Particle(ld)->Pdg();
         if(fpdg==kPdgNeutron && lpdg==kPdgPiM) { is_Delta_1pi = true; ip = ld; }
         if(lpdg==kPdgNeutron && fpdg==kPdgPiM) { is_Delta_1pi = true; ip = fd; }
       }
     }
 
     if(is_Delta_1pi) break;
     i++;
   }
 
   if(!is_Delta_1pi) return 1.;
 
   LOG("ReW", pDEBUG) << "A Delta++ -> p pi+ event:";
   LOG("ReW", pDEBUG) << "Resonance is at position: " << ir;
   LOG("ReW", pDEBUG) << "Pion is at position: " << ip;
 //LOG("ReW", pDEBUG) << event;
 
   // Get Delta and pi+ 4-momentum vectors
   TLorentzVector p4res(*event.Particle(ir)->P4());
   TLorentzVector p4pip(*event.Particle(ip)->P4());
 
   // Boost pi+ to the Delta CM
   TVector3 bv = -1*p4res.BoostVector();
   p4pip.Boost(bv);
 
   //
   // Calculate weight.
   // Angular distribution for pi+, given 2 possible sub-states mi=1/2,3/2
   // is: Wpi(theta) = 1 - p(3/2)*P2(costheta) + p(1/2)*P2(costheta).
   // mi is the projection of Delta angular momentum along quantization axis.
   // P2(costheta) is the 2nd order Legendre polynomial.
   // Isotropy implies equal populations: p(3/2)=p(1/2)=0.5
   // R/S predict: p(3/2)=0.75, p(1/2)=0.25
   // The above follow the notation used in Sam Zeller's invited talk in a 
   // MINOS Physics Simulation mtg on Sep 28th, 2007.
   // Sam's talk cited a paper by G.Garvey and O.Lalakulich which I haven't
   // been able to locate.
   //
 
   double p32iso   = 0.50;
   double p12iso   = 0.50;
   double p32rs    = 0.75;
   double p12rs    = 0.25;
   double costheta = p4pip.Vect().CosTheta();
   double P2       = 0.5 * (3.*costheta*costheta - 1.);
   double Wiso     = 1 - p32iso * P2 + p12iso * P2; // = 1.0
   double Wrs      = 1 - p32rs  * P2 + p12rs  * P2;
   double dial     = fThetaDelta2NpiTwkDial;
   double Wdef     = Wiso;
   double Wtwk     = dial*Wrs + (1-dial)*Wiso;
 
   double wght = 1.;
   if(Wdef>0. && Wtwk>0.) {
      wght = Wtwk/Wdef;
   }
 
   LOG("ReW", pDEBUG) 
        << "Pion Cos(ThetaCM) = " << costheta << ", weight = " << wght;
 
 #ifdef _G_REWEIGHT_RESDEC_DEBUG_
   fTestNtp->Fill(costheta,wght);
 #endif
 
   return wght;
 }

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

virtual

update the value for the specified nuisance param

Implements genie::rew::GReWeightI.

Definition at line 70 of file GReWeightResonanceDecay.cxx.

 {
   switch(syst) {
      case ( kRDcyTwkDial_BR1gamma ) :
         fBR1gammaTwkDial = twk_dial;
         break;
      case ( kRDcyTwkDial_BR1eta ) :
         fBR1etaTwkDial = twk_dial;
         break;
      case ( kRDcyTwkDial_Theta_Delta2Npi  ) :
         fThetaDelta2NpiTwkDial = twk_dial;
         break;
      default:
         return;
         break;
   }
 }

Member Data Documentation

double genie::rew::GReWeightResonanceDecay::fBR1etaTwkDial

private

Definition at line 71 of file GReWeightResonanceDecay.h.

double genie::rew::GReWeightResonanceDecay::fBR1gammaTwkDial

private

Definition at line 70 of file GReWeightResonanceDecay.h.

map<int, TH1D*> genie::rew::GReWeightResonanceDecay::fMpBR1etaDef

private

Definition at line 82 of file GReWeightResonanceDecay.h.

map<int, TH1D*> genie::rew::GReWeightResonanceDecay::fMpBR1gammaDef

private

Definition at line 81 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewCC

private

reweight CC?

Definition at line 78 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewNC

private

reweight NC?

Definition at line 79 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewNue

private

reweight nu_e?

Definition at line 74 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewNuebar

private

reweight nu_e_bar?

Definition at line 75 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewNumu

private

reweight nu_mu?

Definition at line 76 of file GReWeightResonanceDecay.h.

bool genie::rew::GReWeightResonanceDecay::fRewNumubar

private

reweight nu_mu_bar?

Definition at line 77 of file GReWeightResonanceDecay.h.

double genie::rew::GReWeightResonanceDecay::fThetaDelta2NpiTwkDial

private

Definition at line 72 of file GReWeightResonanceDecay.h.

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

GENIE/src/ReWeight/GReWeightResonanceDecay.h
GENIE/src/ReWeight/GReWeightResonanceDecay.cxx

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation