gtestMuELoss.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <string>
#include <vector>
#include <TFile.h>
#include <TNtuple.h>
#include "Framework/AlgorithmAlgFactory.h"
#include "Framework/Conventions/Units.h"
#include "Physics/MuonEnergyLoss/MuELossI.h"
#include "Physics/MuonEnergyLoss/MuELMaterial.h"
#include "Physics/MuonEnergyLoss/MuELProcess.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"

Go to the source code of this file.

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

Variables
string	gOptMaterials

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 153 of file gtestMuELoss.cxx.

{
  LOG("test", pNOTICE) << "Parsing command line arguments";

  CmdLnArgParser parser(argc,argv);

  if ( parser.OptionExists('m') ) {
    LOG("test", pINFO) << "Reading material ids";
    gOptMaterials = parser.ArgAsString('m');
  } else {
    LOG("test", pINFO) << "Unspecified material ids - Exiting";
    exit(1);
  }
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 56 of file gtestMuELoss.cxx.

{
  GetCommandLineArgs(argc, argv);

  const int N = 14;
  double E[N] = {1,5,10,15,20,30,50,100,200,500, 1000,2000, 5000, 9000}; //GeV

  // split the comma separated list of materials
  vector<string> mtv = utils::str::Split(gOptMaterials,  ",");
  
  // get muon energy loss algorithms
  AlgFactory * algf = AlgFactory::Instance();

  const MuELossI * betheBloch = 
         dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
                         "genie::mueloss::BetheBlochModel","Default"));

  const MuELossI * petrukhinShestakov = 
         dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
                 "genie::mueloss::PetrukhinShestakovModel","Default"));

  const MuELossI * kokoulinPetroukhin = 
         dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
                  "genie::mueloss::KokoulinPetrukhinModel","Default"));

  const MuELossI * bezroukovBugaev = 
         dynamic_cast<const MuELossI *> (algf->GetAlgorithm(
                     "genie::mueloss::BezrukovBugaevModel","Default"));

  assert ( betheBloch         );
  assert ( petrukhinShestakov );
  assert ( kokoulinPetroukhin );
  assert ( bezroukovBugaev    );

  double myunits_conversion = units::GeV/(units::g/units::cm2); 
  string myunits_name       = " GeV/(gr/cm^2)";

  // open a ROOT file and define the output ntuple.
  TFile froot("./genie-mueloss.root", "RECREATE");
  TNtuple muntp("muntp","muon dE/dx", "material:E:ion:brem:pair:pnucl");
  
  //loop over materials
  vector<string>::iterator iter;
  for(iter = mtv.begin(); iter != mtv.end(); ++iter) {

    MuELMaterial_t mt = (MuELMaterial_t) atoi(iter->c_str());

     LOG("test", pINFO) 
        << "---------- Computing/Printing muon energy losses in "
                       << MuELMaterial::AsString(mt) << " ----------";

     // loop over energies
     for(int i=0; i<N; i++)  {
      
       // ionization 
       double ion = betheBloch->dE_dx(E[i],mt) / myunits_conversion;

       LOG("test", pINFO) 
         << "Process: " << MuELProcess::AsString(betheBloch->Process())
         << ", Model: " << betheBloch->Id().Key() 
         << " : \n -dE/dx(E=" << E[i] << ") = " << ion << myunits_name;

       // bremsstrahlung
       double brem = petrukhinShestakov->dE_dx(E[i],mt) / myunits_conversion;

       LOG("test", pINFO) 
         << "Process: " << MuELProcess::AsString(petrukhinShestakov->Process())
         << ", Model: " << petrukhinShestakov->Id().Key() 
         << " : \n -dE/dx(E=" << E[i] << ") = " << brem << myunits_name;

       // e-e+ pair production
       double pair = kokoulinPetroukhin->dE_dx(E[i],mt) / myunits_conversion;

       LOG("test", pINFO) 
         << "Process: " << MuELProcess::AsString(kokoulinPetroukhin->Process())
         << ", Model: " << kokoulinPetroukhin->Id().Key() 
         << " : \n -dE/dx(E=" << E[i] << ") = " << pair << myunits_name;

       // photonuclear interactions
       double pnucl = bezroukovBugaev->dE_dx(E[i],mt) / myunits_conversion;

       LOG("test", pINFO) 
         << "Process: " << MuELProcess::AsString(bezroukovBugaev->Process())
         << ", Model: " << bezroukovBugaev->Id().Key() 
         << " : \n -dE/dx(E=" << E[i] << ") = " << pnucl << myunits_name
         << "\n\n";

       muntp.Fill( (int)mt,E[i],ion,brem,pair,pnucl);
    }//e
  }//m

  muntp.Write();
  froot.Close();

  return 0;
}

Variable Documentation

string gOptMaterials

Definition at line 53 of file gtestMuELoss.cxx.