gtestFGPauliBlockSuppr.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*!

\program gtestFGPauliBlockSuppr

\brief   Plot suppression factor due to Pauli-blocking.

\author  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
 University of Liverpool & STFC Rutherford Appleton Laboratory

\created June 20, 2004

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

#include <TFile.h>
#include <TNtuple.h>

#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Physics/NuclearState/NuclearUtils.h"
#include "Physics/NuclearState/FermiMomentumTablePool.h"
#include "Physics/NuclearState/FermiMomentumTable.h"

using namespace genie;

int main(int /*argc*/, char ** /*argv*/)
{
  TNtuple * nt = new TNtuple("nt","","Z:A:nucl:Q2:Rdef:Rm10p:Rp10p");

  const int ntgt = 4;
  const int nnuc = 4;

  int target  [ntgt] = { kPdgTgtDeuterium, kPdgTgtC12, kPdgTgtO16, kPdgTgtFe56 };
  int nucleon [nnuc] = { kPdgProton, kPdgNeutron };

  const int    N        = 3001;
  const double Q2min    = 0.000001*units::GeV2;
  const double Q2max    = 10*units::GeV2;
  const double logQ2min = TMath::Log(Q2min);
  const double logQ2max = TMath::Log(Q2max);
  const double dlogQ2   = (logQ2max-logQ2min)/(N-1);
  const double kMN      = constants::kNucleonMass;
  const double pmax     = 0.5*units::GeV;

  FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
  const FermiMomentumTable * kft  = kftp->GetTable("Default");

  for(int itgt=0; itgt<ntgt; itgt++) {
   for(int inuc=0; inuc<nnuc; inuc++) {

     double kFi = kft->FindClosestKF(target[itgt], nucleon[inuc]);

     int Z = pdg::IonPdgCodeToZ(target[itgt]);
     int A = pdg::IonPdgCodeToA(target[itgt]);

     LOG("test", pNOTICE) 
       << "nuclear target = " << target[itgt] << ", nucleon = " << nucleon[inuc];

     for(int i=0; i<N; i++) {
       double Q2 = TMath::Exp(logQ2min + i*dlogQ2); 
       double Rdef  = utils::nuclear::RQEFG_generic(-1*Q2,kMN,    kFi,    kFi,pmax);
       double Rm10p = utils::nuclear::RQEFG_generic(-1*Q2,kMN,0.9*kFi,0.9*kFi,pmax); // kF -> +10%
       double Rp10p = utils::nuclear::RQEFG_generic(-1*Q2,kMN,1.1*kFi,1.1*kFi,pmax); // kF -> -10%

       LOG("test", pNOTICE) 
         << "Q2 = " << Q2 << " GeV, Rdef = " << Rdef 
         << ", R(kF->0.9kF) = " << Rm10p << ", R(kF->1.1*kF) = " << Rp10p;

       nt->Fill(Z,A,nucleon[inuc],Q2,Rdef,Rm10p,Rp10p);

    }//i
   }//inuc
  }//itgt

  TFile f("./fg_pauli_suppression_factors.root","recreate");
  nt->Write();
  f.Close();

  return 0;
}