gtestRewght.cxx File Reference

#include <string>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include "Framework/EventGen/EventRecord.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Ntuple/NtpMCTreeHeader.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/Messenger/Messenger.h"
#include "Tools/ReWeight/GReWeightI.h"
#include "Tools/ReWeight/GSystSet.h"
#include "Tools/ReWeight/GReWeight.h"
#include "Tools/ReWeight/GReWeightNuXSecCCQE.h"
#include "Tools/ReWeight/GReWeightNuXSecCCQEvec.h"
#include "Tools/ReWeight/GReWeightNuXSecCCRES.h"
#include "Tools/ReWeight/GReWeightNuXSecNCRES.h"
#include "Tools/ReWeight/GReWeightNuXSecDIS.h"
#include "Tools/ReWeight/GReWeightNuXSecCOH.h"
#include "Tools/ReWeight/GReWeightNonResonanceBkg.h"
#include "Tools/ReWeight/GReWeightFGM.h"
#include "Tools/ReWeight/GReWeightDISNuclMod.h"
#include "Tools/ReWeight/GReWeightResonanceDecay.h"
#include "Tools/ReWeight/GReWeightFZone.h"
#include "Tools/ReWeight/GReWeightINuke.h"
#include "Tools/ReWeight/GReWeightAGKY.h"
#include "Framework/Utils/CmdLnArgParser.h"

Include dependency graph for gtestRewght.cxx:

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
int	main (int argc, char **argv)

Variables
int	gOptNEvt
string	gOptInpFilename

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

int main	(	int	argc,
		char **	argv
	)

Definition at line 67 of file gtestRewght.cxx.

References genie::NtpMCEventRecord::Clear(), genie::NtpMCEventRecord::event, GetCommandLineArgs(), gOptInpFilename, gOptNEvt, LOG, and pNOTICE.

{
  GetCommandLineArgs (argc, argv);

  // open the ROOT file and get the TTree & its header
  TTree *           tree = 0;
  NtpMCTreeHeader * thdr = 0;

  TFile file(gOptInpFilename.c_str(),"READ");

  tree = dynamic_cast <TTree *>           ( file.Get("gtree")  );
  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );

  if(!tree) return 1;

  LOG("test", pNOTICE) << "Input tree header: " << *thdr;

  int nev = (gOptNEvt > 0) ?
        TMath::Min(gOptNEvt, (int)tree->GetEntries()) :
        (int) tree->GetEntries();

  LOG("test", pNOTICE) << "Will process " << nev << " events";

  //
  // Create a GReWeight object and add to it a set of 
  // weight calculators
  //

  GReWeight rw;

  rw.AdoptWghtCalc( "xsec_ccqe",       new GReWeightNuXSecCCQE      );
  rw.AdoptWghtCalc( "xsec_ccqe_vec",   new GReWeightNuXSecCCQEvec   );
  rw.AdoptWghtCalc( "xsec_ccres",      new GReWeightNuXSecCCRES     );
  rw.AdoptWghtCalc( "xsec_ncres",      new GReWeightNuXSecNCRES     );
  rw.AdoptWghtCalc( "xsec_nonresbkg",  new GReWeightNonResonanceBkg );
  rw.AdoptWghtCalc( "xsec_dis",        new GReWeightNuXSecDIS       );
  rw.AdoptWghtCalc( "xsec_coh",        new GReWeightNuXSecCOH       );
  rw.AdoptWghtCalc( "nuclear_qe",      new GReWeightFGM             );
  rw.AdoptWghtCalc( "nuclear_dis",     new GReWeightDISNuclMod      );
  rw.AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay  );
  rw.AdoptWghtCalc( "hadro_fzone",     new GReWeightFZone           );
  rw.AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke           );
  rw.AdoptWghtCalc( "hadro_agky",      new GReWeightAGKY            );

  //
  // Create a list of systematic params (more to be found at GSyst.h)
  // set non-default values and re-configure.
  // Weight calculators included above must be able to handle the tweaked params.
  // Each tweaking dial t modifies a physics parameter p as:
  // p_{tweaked} = p_{default} ( 1 + t * dp/p )
  // So setting a tweaking dial to +/-1 modifies a physics quantity
  // by +/- 1sigma.
  // Default fractional errors are defined in GSystUncertainty
  // and can be overriden.
  //

  GSystSet & syst = rw.Systematics();

  syst.Set(kXSecTwkDial_NormCCQE,        +1.0);
  syst.Set(kXSecTwkDial_MaCCQEshape,     +1.0);
  syst.Set(kXSecTwkDial_NormCCRES,       -1.0);
  syst.Set(kXSecTwkDial_VecFFCCQEshape,  -1.0);
  syst.Set(kXSecTwkDial_MaCCRESshape,    -1.0);
  syst.Set(kXSecTwkDial_MvCCRESshape,    +0.5);
  syst.Set(kXSecTwkDial_NormNCRES,       +1.0);
  syst.Set(kXSecTwkDial_MaNCRESshape,    -0.7);
  syst.Set(kXSecTwkDial_MvNCRESshape,    +0.3);
  syst.Set(kXSecTwkDial_RvpCC1pi,        +0.5);
  syst.Set(kXSecTwkDial_RvnCC1pi,        +0.5);
  syst.Set(kXSecTwkDial_MaCOHpi,         -0.5);
  syst.Set(kINukeTwkDial_MFP_pi,         +1.0);
  syst.Set(kINukeTwkDial_MFP_N,          -1.0);
  syst.Set(kINukeTwkDial_FrPiProd_pi,    -0.7);
  syst.Set(kHadrAGKYTwkDial_xF1pi,       -1.0);
  syst.Set(kHadrAGKYTwkDial_pT1pi,       +1.0);
  syst.Set(kHadrNuclTwkDial_FormZone,    +1.0);
  syst.Set(kRDcyTwkDial_Theta_Delta2Npi, +1.0);

  rw.Reconfigure();

  //
  // Concrete weight calculators can be retrieved and fine-tuned.
  // For example:

  GReWeightNuXSecCCQE * rwccqe = 
    dynamic_cast<GReWeightNuXSecCCQE *> (rw.WghtCalc("xsec_ccqe"));
  rwccqe -> RewNue    (false); 
  rwccqe -> RewNuebar (false); 
  rwccqe -> RewNumubar(false); 

  //
  // Event loop
  //

  NtpMCEventRecord * mcrec = 0;
  tree->SetBranchAddress("gmcrec", &mcrec);

  for(int i = 0; i < nev; i++) {
    tree->GetEntry(i);

    EventRecord & event = *(mcrec->event);
    LOG("test", pNOTICE) << event;

    double wght = rw.CalcWeight(event);
    LOG("test", pNOTICE) << "Overall weight = " << wght;

    mcrec->Clear();
  }

  file.Close();

  LOG("test", pNOTICE)  << "Done!";
  return 0;
}

Variable Documentation

string gOptInpFilename

Definition at line 64 of file gtestRewght.cxx.

int gOptNEvt

Definition at line 63 of file gtestRewght.cxx.