gtestAxialFormFactor.cxx File Reference

#include <string>
#include <sstream>
#include <TFile.h>
#include <TTree.h>
#include "Framework/Algorithm/Algorithm.h"
#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Physics/QuasiElastic/XSection/AxialFormFactorModelI.h"
#include "Physics/QuasiElastic/XSection/AxialFormFactor.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/StringUtils.h"

Include dependency graph for gtestAxialFormFactor.cxx:

Go to the source code of this file.

Classes
struct	tdata

Macros
#define	MAX_COEF   3
	Program used for testing / debugging the axial form factor. More...
#define	Q2_LIN   4
#define	Q2_LOG_MIN   0.1
#define	Q2_LOG_MAX   100
#define	N_TREE   100

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	CalculateFormFactor (AxialFormFactor axff, TTree *affnt, tdata params, AlgFactory *algf, Registry *r, const Registry *gc, Interaction *interaction, double t0)
bool	IncrementCoefficients (double *coefmin, double *coefmax, double *coefinc, int kmaxinc, tdata params, AlgFactory *algf, Registry *r, const Registry *gc)
int	main (int argc, char **argv)

Variables
string	kDefOptEvFilePrefix = "test.axialff"
string	gOptEvFilePrefix
bool	gOptDoLogarithmicQ2
double	gOptCoeffMin [MAX_COEF] = {0.}
double	gOptCoeffMax [MAX_COEF] = {0.}
double	gOptCoeffInc [MAX_COEF] = {0.}
int	gOptKmaxInc = 0

Macro Definition Documentation

#define MAX_COEF   3

Program used for testing / debugging the axial form factor.

gtestAxialFormFactor

Author

Aaron Meyer <asmeyer2012 uchicago.edu> based off testElFormFactors by Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Created:

August 20, 2013

gtestAxialFormFactor [-l] [-o output_filename_prefix] [-c min,max,inc[,min,max,inc...]]

[] denotes an optionan argument -l switch which evaluates over logarithmic Q^2 -o output filename prefix for root file output -c allows scanning over z-expansion coefficients – must give a multiple of 3 numbers as arguments: min, max, increment – will scan over at most MAX_COEF coefficients and fill ntuple tree with all entries – can change MAX_COEF and recompile as necessary – to skip scanning over a coefficient, put max < min for that coefficient note that Kmax in UserPhysicsOptions should match the number of scanned coeffs

License:

Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 57 of file gtestAxialFormFactor.cxx.

Referenced by main().

#define N_TREE   100

Definition at line 62 of file gtestAxialFormFactor.cxx.

Referenced by CalculateFormFactor().

#define Q2_LIN   4

Definition at line 59 of file gtestAxialFormFactor.cxx.

Referenced by CalculateFormFactor().

#define Q2_LOG_MAX   100

Definition at line 61 of file gtestAxialFormFactor.cxx.

Referenced by CalculateFormFactor().

#define Q2_LOG_MIN   0.1

Definition at line 60 of file gtestAxialFormFactor.cxx.

Referenced by CalculateFormFactor().

Function Documentation

void CalculateFormFactor	(	AxialFormFactor	axff,
		TTree *	affnt,
		tdata	params,
		AlgFactory *	algf,
		Registry *	r,
		const Registry *	gc,
		Interaction *	interaction,
		double	t0
	)

Definition at line 227 of file gtestAxialFormFactor.cxx.

References genie::AxialFormFactor::Calculate(), genie::AxialFormFactor::FA(), genie::AlgFactory::GetAlgorithm(), gOptDoLogarithmicQ2, genie::Interaction::KinePtr(), genie::constants::kPi0Mass, LOG, N_TREE, pERROR, tdata::pFA, tdata::pmod, tdata::pQ2, tdata::pz, genie::Kinematics::q2(), genie::utils::kinematics::Q2(), Q2_LIN, Q2_LOG_MAX, Q2_LOG_MIN, genie::AxialFormFactor::SetModel(), and genie::Kinematics::SetQ2().

Referenced by main().

{
  const AxialFormFactorModelI * dipole =
      dynamic_cast<const AxialFormFactorModelI *> (
        algf->GetAlgorithm("genie::DipoleAxialFormFactorModel", "Default"));
  const AxialFormFactorModelI * zexp =
      dynamic_cast<const AxialFormFactorModelI *> (
        algf->GetAlgorithm("genie::ZExpAxialFormFactorModel", "Default"));

  interaction->KinePtr()->SetQ2(0.);
  double t     = interaction->KinePtr()->q2();
  double tcut  = 9.0 * TMath::Power(constants::kPi0Mass, 2);
  double Q2    = 0.;
                                                                                            
  for(int iq=0; iq<N_TREE; iq++) { 
                                                                                            
   if (gOptDoLogarithmicQ2)
   {
     Q2 = TMath::Exp(  TMath::Log(double(Q2_LOG_MIN))
                     + (double(iq+1)/double(N_TREE))
                     * TMath::Log(double(Q2_LOG_MAX)/double(Q2_LOG_MIN)));
     //Q2 = TMath::Exp(((iq+1)*double(Q2_RANGE_LOG/N_TREE))*TMath::Log(10.)); // logarithmic
   } else { Q2 = (iq+1)*double(Q2_LIN)/double(N_TREE) ; } // linear 

   interaction->KinePtr()->SetQ2(Q2);
   *params.pQ2 = Q2;
                                                                                            
   // calculate z parameter used in expansion
   t = interaction->KinePtr()->q2();
   double znum  = TMath::Sqrt(tcut - t) - TMath::Sqrt(tcut - t0);
   double zden  = TMath::Sqrt(tcut - t) + TMath::Sqrt(tcut - t0);
   double zpar = znum/zden;
                                                                                            
   if (zpar != zpar) LOG("testAxialFormFactor", pERROR) << "Undefined expansion parameter";
   *params.pz = zpar;

   axff.SetModel(dipole);
   axff.Calculate(interaction);
   *params.pmod = 0;
   *params.pFA  = axff.FA();
   affnt->Fill();
                                                                                            
   axff.SetModel(zexp);
   axff.Calculate(interaction);
   *params.pmod = 1;
   *params.pFA = axff.FA();
   affnt->Fill();
   
  }
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

bool IncrementCoefficients	(	double *	coefmin,
		double *	coefmax,
		double *	coefinc,
		int	kmaxinc,
		tdata	params,
		AlgFactory *	algf,
		Registry *	r,
		const Registry *	gc
	)

Definition at line 282 of file gtestAxialFormFactor.cxx.

References genie::AlgFactory::ForceReconfiguration(), LOG, tdata::pcAn, pERROR, and genie::Registry::Set().

Referenced by main().

{

  if (kmaxinc < 1)
  {
    LOG("testAxialFormFactor",pERROR) << "No coefficients to increment";
    return false;
  } else {

  ostringstream alg_key;
  int ip = -1;
  bool stopflag;
  do
  {
    if (ip > -1)
    { // a previous iteration went over max
      params.pcAn[ip] = coefmin[ip];
      r->Set(alg_key.str(),params.pcAn[ip]);
    }
    stopflag = true;
    alg_key.str("");                                      // reset sstream
    ip++;                                                 // increment index
    if (ip == kmaxinc) return false;                      // stop if gone too far
    alg_key << "QEL-Z_A" << ip+1;                         // get new name
    params.pcAn[ip] += coefinc[ip];                       // increment parameter
    r->Set(alg_key.str(),params.pcAn[ip]);
    if (params.pcAn[ip] > coefmax[ip]) stopflag=false;    

  } while (! stopflag); // loop


  } // if kmaxinc < 1

  algf->ForceReconfiguration();
  return true;

}

int main	(	int	argc,
		char **	argv
	)

Definition at line 98 of file gtestAxialFormFactor.cxx.

References CalculateFormFactor(), genie::AlgConfigPool::FindRegistry(), genie::AlgFactory::ForceReconfiguration(), genie::AlgFactory::GetAlgorithm(), genie::Registry::GetBool(), genie::Registry::GetBoolDef(), GetCommandLineArgs(), genie::Registry::GetDouble(), genie::Registry::GetDoubleDef(), genie::Registry::GetInt(), genie::Registry::GetIntDef(), genie::AlgConfigPool::GlobalParameterList(), gOptCoeffInc, gOptCoeffMax, gOptCoeffMin, gOptEvFilePrefix, gOptKmaxInc, IncrementCoefficients(), genie::AlgFactory::Instance(), genie::AlgConfigPool::Instance(), genie::kPdgNuMu, genie::kPdgProton, genie::kPdgTgtFe56, LOG, MAX_COEF, tdata::pcAn, tdata::pFA, pINFO, tdata::pKmax, tdata::pmod, tdata::pQ2, pWARN, tdata::pz, genie::utils::kinematics::Q2(), genie::Interaction::QELCC(), genie::Registry::Set(), and genie::Registry::UnLock().

{

  GetCommandLineArgs(argc,argv);

  if (gOptKmaxInc > 0)
  {
    LOG("testAxialFormFactor", pINFO) << "Found coefficient ranges:";
    for (int ik=0;ik<gOptKmaxInc;ik++)
    {
      LOG("testAxialFormFactor", pINFO) << "( min:"  << gOptCoeffMin[ik]  \
                                        << " , max:" << gOptCoeffMax[ik]  \
                                        << " , inc:" << gOptCoeffInc[ik] << " )";
    }
  }

  TTree * affnt = new TTree("affnt","Axial Form Factor Test Tree");
  AxialFormFactor axff;
  AlgFactory * algf = AlgFactory::Instance();
  AlgConfigPool * confp = AlgConfigPool::Instance();
  const Registry * gc = confp->GlobalParameterList();

  AlgId id("genie::ZExpAxialFormFactorModel","Default");
  const Algorithm * alg = algf->GetAlgorithm(id);
  Registry * r = confp->FindRegistry(alg);

  int Kmax;
  int mod;
  double Q2;
  double z;
  double FA;
  double cAn[MAX_COEF];
  affnt = new TTree("axff","Axial Form Factor Test");
  affnt->Branch("Kmax", &Kmax, "Kmax/I"); // number of coefficients
  affnt->Branch("mod",  &mod,  "mod/I" ); // model
  affnt->Branch("Q2",   &Q2,   "Q2/D"  ); // Q^2
  affnt->Branch("z" ,   &z ,   "z/D"   ); // z
  affnt->Branch("FA",   &FA,   "FA/D"  ); // F_A axial form factor
  affnt->Branch("cAn",   cAn,  "cAn[Kmax]/D"); // z-expansion coefficients

  // load all parameters into a struct to make passing them as arguments easier
  tdata params; //defined tdata struct above
  params.pKmax = &Kmax;
  params.pmod  = &mod;
  params.pQ2   = &Q2;
  params.pz    = &z;
  params.pFA   = &FA;
  params.pcAn  = cAn;

  // flag for single output or loop
  bool do_once = (gOptKmaxInc < 1);
  bool do_Q4limit = r->GetBoolDef("QEL-Q4limit", gc->GetBool("QEL-Q4limit"));

  if (gOptKmaxInc < 1) Kmax = r->GetIntDef("QEL-Kmax", gc->GetInt("QEL-Kmax")); 
  else                 Kmax = gOptKmaxInc;
  //if (do_Q4limit) Kmax = Kmax+4;
  if (do_Q4limit) LOG("testAxialFormFactor",pWARN) \
   << "Q4limit specified, note that coefficients will be altered";

  // initialize to zero
  for (int ip=0;ip<MAX_COEF;ip++)
  {
    params.pcAn[ip] = 0.;
  }

  // Get T0 from registry value
  double t0 = r->GetDoubleDef("QEL-T0", gc->GetDouble("QEL-T0")); 

  // make sure coefficients are getting incremented correctly
  if (MAX_COEF < gOptKmaxInc ||
     (gOptKmaxInc == 0 && MAX_COEF < Kmax) )
  {
    LOG("testAxialFormFactor",pWARN) \
      << "Too many coefficient ranges, reducing to " << MAX_COEF;
    Kmax        = MAX_COEF;
    gOptKmaxInc = MAX_COEF;
  }
  LOG("testAxialFormFactor",pWARN) << "pKmax = " << *params.pKmax;

  // set up interaction
  Interaction * interaction = 
         Interaction::QELCC(kPdgTgtFe56,kPdgProton,kPdgNuMu,0.02);

  if (do_once)
  { // calculate once and be done with it

    // pre-load parameters
    ostringstream alg_key;
    for (int ip=1;ip<Kmax+1;ip++)
    {
      alg_key.str("");
      alg_key << "QEL-Z_A" << ip;
      params.pcAn[ip-1] =
        r->GetDoubleDef(alg_key.str(), gc->GetDouble(alg_key.str())); 
      LOG("testAxialFormFactor",pINFO) << "Loading " << alg_key.str() \
       << " : " << params.pcAn[ip-1];
    }

    CalculateFormFactor(axff,affnt,params,algf,r,gc,interaction,t0);//,t0,Kmax);
  } else {
    // override and control coefficient values
    r->UnLock();
    r->Set("QEL-Kmax",*params.pKmax);

    ostringstream alg_key;
    for (int ip=1;ip<gOptKmaxInc+1;ip++)
    {
      alg_key.str("");
      alg_key << "QEL-Z_A" << ip;
      r->Set(alg_key.str(),gOptCoeffMin[ip-1]);
      params.pcAn[ip-1] = gOptCoeffMin[ip-1];
      LOG("testAxialFormFactor",pWARN) << "Set parameter: " << params.pcAn[ip-1];
    }
    algf->ForceReconfiguration();

    do
    { CalculateFormFactor(axff,affnt,params,algf,r,gc,interaction,t0); }//,t0,Kmax); }
    while (IncrementCoefficients(gOptCoeffMin,gOptCoeffMax,gOptCoeffInc,gOptKmaxInc,
                                 params,algf,r,gc));

  }

  TFile f((gOptEvFilePrefix + ".root").c_str(),"recreate");
  affnt->Write();
  f.Close();

  return 0;
}

Variable Documentation

double gOptCoeffInc[MAX_COEF] = {0.}

Definition at line 94 of file gtestAxialFormFactor.cxx.

Referenced by main().

double gOptCoeffMax[MAX_COEF] = {0.}

Definition at line 93 of file gtestAxialFormFactor.cxx.

Referenced by main().

double gOptCoeffMin[MAX_COEF] = {0.}

Definition at line 92 of file gtestAxialFormFactor.cxx.

Referenced by main().

bool gOptDoLogarithmicQ2

Definition at line 91 of file gtestAxialFormFactor.cxx.

Referenced by CalculateFormFactor().

string gOptEvFilePrefix

Definition at line 90 of file gtestAxialFormFactor.cxx.

int gOptKmaxInc = 0

Definition at line 95 of file gtestAxialFormFactor.cxx.

Referenced by main().

string kDefOptEvFilePrefix = "test.axialff"

Definition at line 89 of file gtestAxialFormFactor.cxx.