LBNEPrimaryMessenger.cc

Go to the documentation of this file.

#include "LBNEPrimaryMessenger.hh"

#include "LBNEPrimaryGeneratorAction.hh"
#include "G4UIdirectory.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWith3VectorAndUnit.hh"
#include "G4UIcmdWithADouble.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithABool.hh"
#include "G4UIcmdWithoutParameter.hh"
#include "G4ios.hh"
#include "globals.hh"
#include "Randomize.hh"
#include "LBNERunManager.hh"

LBNEPrimaryMessenger::LBNEPrimaryMessenger(LBNEPrimaryGeneratorAction* RA)
  :fPrimaryAction (RA)
{
   LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
   
   if(theRunManager->GetVerboseLevel() > 0)
   {
      G4cout << "LBNEPrimaryMessenger Constructor Called." << G4endl;
   }
  fDirectory = new G4UIdirectory("/LBNE/generator/");

  fBeamOffsetXCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamOffsetX",this);
  fBeamOffsetXCmd->SetGuidance("Set the X offset of the proton beam");
  fBeamOffsetXCmd->SetParameterName("beamOffsetX", false);
  fBeamOffsetXCmd->SetUnitCategory("Length");
  fBeamOffsetXCmd->SetDefaultValue (0.);
  fBeamOffsetXCmd->SetDefaultUnit ("mm");
  fBeamOffsetXCmd->SetUnitCandidates ("mm cm");
  
  fBeamOffsetYCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamOffsetY",this);
  fBeamOffsetYCmd->SetGuidance("Set the Y offset of the proton beam");
  fBeamOffsetYCmd->SetParameterName("beamOffsetY", false);
  fBeamOffsetYCmd->SetUnitCategory("Length");
  fBeamOffsetYCmd->SetDefaultValue (0.);
  fBeamOffsetYCmd->SetDefaultUnit ("mm");
  fBeamOffsetYCmd->SetUnitCandidates ("mm cm");
 
  fBeamSigmaXCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamSigmaX",this);
  fBeamSigmaXCmd->SetGuidance("Set the X Sigma of the proton beam");
  fBeamSigmaXCmd->SetParameterName("beamSigmaX", false);
  fBeamSigmaXCmd->SetUnitCategory("Length");
  fBeamSigmaXCmd->SetDefaultValue (fPrimaryAction->GetBeamSigmaX());
  fBeamSigmaXCmd->SetDefaultUnit ("mm");
  fBeamSigmaXCmd->SetUnitCandidates ("mm cm");
  
  fBeamSigmaYCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamSigmaY",this);
  fBeamSigmaYCmd->SetGuidance("Set the Y Sigma of the proton beam");
  fBeamSigmaYCmd->SetParameterName("beamSigmaY", false);
  fBeamSigmaYCmd->SetUnitCategory("Length");
  fBeamSigmaYCmd->SetDefaultValue (fPrimaryAction->GetBeamSigmaY());
  fBeamSigmaYCmd->SetDefaultUnit ("mm");
  fBeamSigmaYCmd->SetUnitCandidates ("mm cm");
 
  fBeamMaxValXCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamMaxValX",this);
  fBeamMaxValXCmd->SetGuidance("Set the maximum value in Y that a proton can have on target");
  fBeamMaxValXCmd->SetParameterName("beamMaxValX", false);
  fBeamMaxValXCmd->SetUnitCategory("Length");
  fBeamMaxValXCmd->SetDefaultValue (fPrimaryAction->GetBeamMaxValX());
  fBeamMaxValXCmd->SetDefaultUnit ("mm");
  fBeamMaxValXCmd->SetUnitCandidates ("mm cm");
  
  fBeamMaxValYCmd = 
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamMaxValY",this);
  fBeamMaxValYCmd->SetGuidance("Set the maximum value in Y that a proton can have on target");
  fBeamMaxValYCmd->SetParameterName("beamMaxValY", false);
  fBeamMaxValYCmd->SetUnitCategory("Length");
  fBeamMaxValYCmd->SetDefaultValue (fPrimaryAction->GetBeamMaxValY());
  fBeamMaxValYCmd->SetDefaultUnit ("mm");
  fBeamMaxValYCmd->SetUnitCandidates ("mm cm");
 
  fBeamThetaCmd =
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamTheta",this);
  fBeamThetaCmd->SetGuidance("Set the angle (theta) of the proton beam");
  fBeamThetaCmd->SetParameterName("beamTheta", false);
  fBeamThetaCmd->SetUnitCategory("Angle");
  
  fBeamPhiCmd =
    new G4UIcmdWithADoubleAndUnit("/LBNE/generator/beamPhi",this);
  fBeamPhiCmd->SetGuidance("Set the angle (phi) of the proton beam.");
  fBeamPhiCmd->SetParameterName("beamPhi", false);
  fBeamPhiCmd->SetUnitCategory("Angle");

  fCorrectForAngleCmd = 
    new G4UIcmdWithABool("/LBNE/generator/correctForAngle", this);
  fCorrectForAngleCmd->SetGuidance("If true, beam x/y position is corrected");
  fCorrectForAngleCmd->SetGuidance("to hit center of target using angle of");
  fCorrectForAngleCmd->SetGuidance("beam. Offsets in x and y specified via");
  fCorrectForAngleCmd->SetGuidance("messenger are still respected.");
  
  fProtonMomentum  = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/protonMomentum",this);
  fProtonMomentum->SetGuidance("Primary proton momentum delivered by the Fermilab Main Injector ");
  fProtonMomentum->SetParameterName("protonMomentum",true);
  fProtonMomentum->SetDefaultValue (120.0*CLHEP::GeV);
  fProtonMomentum->SetDefaultUnit ("GeV");
  fProtonMomentum->SetUnitCandidates ("GeV");
  fProtonMomentum->AvailableForStates(G4State_Idle);
    
  fBeamOnTargetCmd = 
    new G4UIcmdWithABool("/LBNE/generator/beamOnTarget", this);
  fBeamOnTargetCmd->SetGuidance("If true, forces beam to hit the center");
  fBeamOnTargetCmd->SetGuidance("of target. Any x or y offsets supplied");
  fBeamOnTargetCmd->SetGuidance("via messenger are ignored");
  
  fBeamBetaFunctionX = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/beamBetaFunctionX", this);
  fBeamBetaFunctionX->SetGuidance(
   "LBNE beam line X beta function at the the target (MCZERO) delived by the Fermilab Main Injector ");
  fBeamBetaFunctionX->SetParameterName("beamBetaFunctionX",true);
  fBeamBetaFunctionX->SetDefaultValue (64.842); // ugly, but the primary geenrator class not yet instantiated. 
  fBeamBetaFunctionX->SetDefaultUnit ("m");
  fBeamBetaFunctionX->SetUnitCandidates ("m");
  fBeamBetaFunctionX->AvailableForStates(G4State_Idle);
  
  fBeamBetaFunctionY = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/beamBetaFunctionY", this);
  fBeamBetaFunctionX->SetGuidance(
   "LBNE beam line Y beta function at the the target (MCZERO) delived by the Fermilab Main Injector ");
  fBeamBetaFunctionY->SetParameterName("beamBetaFunctionY",true);
  fBeamBetaFunctionY->SetDefaultValue (64.842); // ugly, but the primary geenrator class not yet instantiated. 
  fBeamBetaFunctionY->SetDefaultUnit ("m");
  fBeamBetaFunctionY->SetUnitCandidates ("m");
  fBeamBetaFunctionY->AvailableForStates(G4State_Idle);

  fBeamEmittanceX = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/beamEmittanceX", this);
  fBeamEmittanceX->SetGuidance(
   "LBNE beam line X Emittance, in pi mm mRad, Fermilab units. ");
  fBeamEmittanceX->SetParameterName("beamEmittanceX",true);
  fBeamEmittanceX->SetDefaultValue (20.); // ugly, but the primary geenrator class not yet instantiated. 
  fBeamEmittanceX->SetDefaultUnit ("mm"); // a lie, but define a unit will only bring more confusion. 
  fBeamEmittanceX->SetUnitCandidates ("mm");
  fBeamEmittanceX->AvailableForStates(G4State_Idle);
  
  fBeamEmittanceY = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/beamEmittanceY", this);
  fBeamEmittanceY->SetGuidance(
   "LBNE beam line Y Emittance, in pi mm mRad, Fermilab units. ");
  fBeamEmittanceY->SetParameterName("beamEmittanceY",true);
  fBeamEmittanceY->SetDefaultValue (20.); // ugly, but the primary geenrator class not yet instantiated. 
  fBeamEmittanceY->SetDefaultUnit ("mm"); // a lie, but define a unit will only bring more confusion. 
  fBeamEmittanceY->SetUnitCandidates ("mm");
  fBeamEmittanceY->AvailableForStates(G4State_Idle);
  
  fRadiusAnnularBeam = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/annularBeamRadius", this);
  fRadiusAnnularBeam->SetGuidance(
   "LBNF Annular beam option. Setting the radius of the annulus. Width set by XY beam size. ");
  fRadiusAnnularBeam->SetParameterName("Annular",true);
  fRadiusAnnularBeam->SetDefaultValue (-1.); // Default is no annulus.  
  fRadiusAnnularBeam->SetDefaultUnit ("mm"); // a lie, but define a unit will only bring more confusion. 
  fRadiusAnnularBeam->SetUnitCandidates ("mm");
  fRadiusAnnularBeam->AvailableForStates(G4State_Idle);
  
    
  fUseGeantino  = new G4UIcmdWithoutParameter("/LBNE/primary/useGeantino",this);
  fUseGeantino->SetGuidance("Using a Geantino at the Primary, to study absorption");
  fUseGeantino->AvailableForStates(G4State_Idle);
    
  fUseMuonGeantino  = new G4UIcmdWithoutParameter("/LBNE/primary/useMuonGeantino",this);
  fUseMuonGeantino->SetGuidance("Using a muon at the Primary, to study absorption, with magnetic field effect ");
  fUseMuonGeantino->AvailableForStates(G4State_Idle);
  
  fUseChargedGeantino  = new G4UIcmdWithoutParameter("/LBNE/primary/useChargedGeantino",this);
  fUseChargedGeantino->SetGuidance("Using a particle with charge, but no other physics (I hope)  magnetic field effect ");
  fUseChargedGeantino->AvailableForStates(G4State_Idle);
  
  fGeantinoOpeningAngle  = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/geantinoOpeningAngle",this);
  fGeantinoOpeningAngle->SetGuidance("Polar angle generating the geantino (or mu geantino)  ");
  fGeantinoOpeningAngle->SetParameterName("GeantinoOpeningAngle",true);
  fGeantinoOpeningAngle->SetDefaultValue (0.005*CLHEP::radian);
  fGeantinoOpeningAngle->SetDefaultUnit("radian");
  fGeantinoOpeningAngle->SetUnitCandidates("radian");
  fGeantinoOpeningAngle->AvailableForStates(G4State_Idle);
   
  fGeantinoOpeningAngleMin  = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/geantinoOpeningAngleMin",this);
  fGeantinoOpeningAngleMin->SetGuidance("Minimum Polar angle generating the geantino (or mu geantino)  ");
  fGeantinoOpeningAngleMin->SetParameterName("GeantinoOpeningAngleMin",true);
  fGeantinoOpeningAngleMin->SetDefaultValue (0.);
  fGeantinoOpeningAngleMin->SetDefaultUnit("radian");
  fGeantinoOpeningAngleMin->SetUnitCandidates("radian");
  fGeantinoOpeningAngleMin->AvailableForStates(G4State_Idle);
   
  fGeantinoZOrigin  = new G4UIcmdWithADoubleAndUnit("/LBNE/primary/geantinoZOrigin",this);
  fGeantinoZOrigin->SetGuidance("Z origin  generating the geantino (or mu geantino) (in mm) ");
  fGeantinoZOrigin->SetParameterName("GeantinoOpeningAngle",true);
  fGeantinoZOrigin->SetDefaultValue (-515.);
  fGeantinoZOrigin->SetDefaultUnit ("mm");
  fGeantinoZOrigin->SetUnitCandidates ("mm cm m");
  fGeantinoZOrigin->AvailableForStates(G4State_Idle);
  
  fGeantinoZOriginSigma  = 
     new G4UIcmdWithADoubleAndUnit("/LBNE/primary/geantinoSigmaZOrigin",this);
  fGeantinoZOriginSigma->SetGuidance("Z origin  longitudinal spread generating the geantino (or mu geantino) (in mm) ");
  fGeantinoZOriginSigma->SetParameterName("GeantinoSigmaZOrigin",true);
  fGeantinoZOriginSigma->SetDefaultValue (100.);
  fGeantinoZOriginSigma->SetDefaultUnit ("mm");
  fGeantinoZOriginSigma->SetUnitCandidates ("mm cm m");
  fGeantinoZOriginSigma->AvailableForStates(G4State_Idle);

  fMuGeantinoZOriginScan  = 
     new G4UIcmdWith3VectorAndUnit("/LBNE/primary/mugeantinoScanZOrigin",this);
  fMuGeantinoZOriginScan->SetGuidance("Z Scan on the origin of the neutrino ");
  fMuGeantinoZOriginScan->SetParameterName("mugeantinoScanZOriginInit", "mugeantinoScanZOriginStep", "mugeantinoScanZOriginFlag", true);
  fMuGeantinoZOriginScan->SetDefaultValue (G4ThreeVector(0., 0.000, 0.)); // Z start, z Step, flag (0 = no, 1 = scan ) 
  fMuGeantinoZOriginScan->SetDefaultUnit ("mm");
  fMuGeantinoZOriginScan->SetUnitCandidates ("mm cm m");
  fMuGeantinoZOriginScan->AvailableForStates(G4State_Idle);

  fMuGeantinoAngleScan  = 
     new G4UIcmdWith3VectorAndUnit("/LBNE/primary/mugeantinoScanAngle",this);
  fMuGeantinoAngleScan->SetGuidance("Angle Scan on the origin of the neutrino "); // in Y, actually.. 
  fMuGeantinoAngleScan->SetParameterName("mugeantinoScanAlngleInit", "mugeantinoScanAlngleStep", "mugeantinoScanAlngleFlag", true);
  fMuGeantinoAngleScan->SetDefaultValue (G4ThreeVector(0.01, 0.0001, 0.)); //  
  fMuGeantinoAngleScan->SetDefaultUnit ("radian");
  fMuGeantinoAngleScan->SetUnitCandidates ("radian");
  fMuGeantinoAngleScan->AvailableForStates(G4State_Idle);

  fMuGeantinoMomentumScan  = 
     new G4UIcmdWith3VectorAndUnit("/LBNE/primary/mugeantinoScanMomentum",this);
  fMuGeantinoMomentumScan->SetGuidance("Momentum Scan on the origin of the neutrino ");
  fMuGeantinoMomentumScan->SetParameterName("mugeantinoScanMomentumInit", "mugeantinoScanMomentumStep","mugeantinoScanMomentumFlag", true);
  fMuGeantinoMomentumScan->SetDefaultValue (G4ThreeVector(10., 0.000, 0.)); // Z start, z Step, flag (0 = no, 1 = scan ) 
  fMuGeantinoMomentumScan->SetDefaultUnit ("GeV");
  fMuGeantinoMomentumScan->SetUnitCandidates ("MeV GeV");
  fMuGeantinoMomentumScan->AvailableForStates(G4State_Idle);

  fMuGeantinoPtScan  = 
     new G4UIcmdWith3VectorAndUnit("/LBNE/primary/mugeantinoScanPt",this);
  fMuGeantinoPtScan->SetGuidance("Pt Scan on the origin of the neutrino ");
  fMuGeantinoPtScan->SetParameterName("mugeantinoScanPtInit", "mugeantinoScanPtStep","mugeantinoScanPtFlag", true);
  fMuGeantinoPtScan->SetDefaultValue (G4ThreeVector(10., 0.000, 0.)); // Z start, z Step, flag (0 = no, 1 = scan ) 
  fMuGeantinoPtScan->SetDefaultUnit ("GeV");
  fMuGeantinoPtScan->SetUnitCandidates ("MeV GeV");
  fMuGeantinoPtScan->AvailableForStates(G4State_Idle);

  fMuGeantinoYOriginScan  = 
     new G4UIcmdWith3VectorAndUnit("/LBNE/primary/mugeantinoScanYOrigin",this);
  fMuGeantinoYOriginScan->SetGuidance("Y Scan on the origin of the neutrino ");
  fMuGeantinoYOriginScan->SetParameterName("mugeantinoScanZOriginInit", "mugeantinoScanZOriginStep","mugeantinoScanZOriginFlag", true);
  fMuGeantinoYOriginScan->SetDefaultValue (G4ThreeVector(0., 0.000, 0.)); // Z start, z Step, flag (0 = no, 1 = scan ) 
  fMuGeantinoYOriginScan->SetDefaultUnit ("mm");
  fMuGeantinoYOriginScan->SetUnitCandidates ("mm cm m");
  fMuGeantinoYOriginScan->AvailableForStates(G4State_Idle);

  fInputFlukaFileName = new G4UIcmdWithAString("/LBNE/run/InputFlukaFileName",this);
  std::string aGuideFluka("Set the input fluka ntuple file name. \n");
  aGuideFluka +=  std::string("Proton beam will be off, we will generate hadrons from Fluka file \n");
  aGuideFluka +=  std::string("Note: No checks are done on the consistency of Fluka vs this Geant geometry \n");
  fInputFlukaFileName->SetGuidance(aGuideFluka.c_str());
  fInputFlukaFileName->SetParameterName("extNtupleFileName",true);
  fInputFlukaFileName->SetDefaultValue (G4String("none"));
  fInputFlukaFileName->AvailableForStates(G4State_PreInit,G4State_Idle);
  
}

LBNEPrimaryMessenger::~LBNEPrimaryMessenger()
{
  delete fBeamOffsetXCmd;
  delete fBeamOffsetYCmd;
  delete fBeamSigmaXCmd;
  delete fBeamSigmaYCmd;
  delete fBeamMaxValXCmd;
  delete fBeamMaxValYCmd;
  delete fBeamPhiCmd;
  delete fBeamThetaCmd;
  delete fCorrectForAngleCmd;
  delete fBeamOnTargetCmd;
  delete fDirectory;
  delete fUseGeantino;
  delete fProtonMomentum;
  delete fUseMuonGeantino;
  delete fUseChargedGeantino;
  delete fGeantinoOpeningAngle;
  delete fGeantinoZOrigin;
  delete fGeantinoZOriginSigma;
  delete fMuGeantinoZOriginScan;
  delete fMuGeantinoAngleScan; 
  delete fMuGeantinoMomentumScan; 
  delete fMuGeantinoPtScan; 
  delete fMuGeantinoYOriginScan; 
  delete fInputFlukaFileName;
  
}

void LBNEPrimaryMessenger::SetNewValue(G4UIcommand* cmd, G4String val)
{
   LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

   if(theRunManager->GetVerboseLevel() > 1)
   {
      G4cout << "LBNEPrimaryMessenger::SetNewValue - Setting Parameter value from input macro." << G4endl;
   }


   if(theRunManager->GetVerboseLevel() > 1)
   {
      G4cout << "LBNEPrimaryMessenger::SetNewValue - Done Setting parameter value." << G4endl;
   }

  if(cmd == fBeamOffsetXCmd){
    fPrimaryAction->SetBeamOffsetX(fBeamOffsetXCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamOffsetYCmd){
    fPrimaryAction->SetBeamOffsetY(fBeamOffsetYCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamSigmaXCmd){
    fPrimaryAction->SetBeamSigmaX(fBeamSigmaXCmd->GetNewDoubleValue(val)); 
    fPrimaryAction->SetUseJustSigmaCoord(true); 
    fPrimaryAction->SetUseCourantSniderParams(false);
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamSigmaYCmd){
    fPrimaryAction->SetBeamSigmaY(fBeamSigmaYCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetUseJustSigmaCoord(true); 
    fPrimaryAction->SetUseCourantSniderParams(false);
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamMaxValXCmd){
    fPrimaryAction->SetBeamMaxValX(fBeamMaxValXCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamMaxValYCmd){
    fPrimaryAction->SetBeamMaxValY(fBeamMaxValYCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamThetaCmd){
    fPrimaryAction->SetBeamTheta(fBeamThetaCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamPhiCmd){
    fPrimaryAction->SetBeamPhi(fBeamPhiCmd->GetNewDoubleValue(val));   
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fCorrectForAngleCmd){
    fPrimaryAction->SetCorrectForAngle(fCorrectForAngleCmd->GetNewBoolValue(val));
    fPrimaryAction->SetProtonBeam();
  }
  if(cmd == fBeamOnTargetCmd){
    fPrimaryAction->SetBeamOnTarget(fBeamOnTargetCmd->GetNewBoolValue(val));
    fPrimaryAction->SetProtonBeam();
  }
  if (cmd == fGeantinoOpeningAngle) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetPolarAngleGeantino(cmdWD->GetNewDoubleValue(val));
  } else if (cmd ==  fGeantinoZOrigin ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetZOriginGeantino( cmdWD->GetNewDoubleValue(val));   
  } else if (cmd ==  fGeantinoOpeningAngleMin ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetPolarAngleGeantinoMin(cmdWD->GetNewDoubleValue(val));   
  } else if (cmd ==  fGeantinoZOriginSigma ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetSigmaZOriginGeantino( cmdWD->GetNewDoubleValue(val));   
  } else if (cmd ==  fMuGeantinoZOriginScan ) {
      G4UIcmdWith3VectorAndUnit* cmdWD = dynamic_cast<G4UIcmdWith3VectorAndUnit*> (cmd);
      fPrimaryAction->SetZOriginForMuGeantinoScan( cmdWD->GetNew3VectorValue(val));   
  } else if (cmd ==  fMuGeantinoMomentumScan ) {
      G4UIcmdWith3VectorAndUnit* cmdWD = dynamic_cast<G4UIcmdWith3VectorAndUnit*> (cmd);
      fPrimaryAction->SetMomentumForMuGeantinoScan( cmdWD->GetNew3VectorValue(val));   
  } else if (cmd ==  fMuGeantinoPtScan ) {
      G4UIcmdWith3VectorAndUnit* cmdWD = dynamic_cast<G4UIcmdWith3VectorAndUnit*> (cmd);
      fPrimaryAction->SetPtForMuGeantinoScan( cmdWD->GetNew3VectorValue(val));   
  } else if (cmd ==  fMuGeantinoAngleScan ) {
      G4UIcmdWith3VectorAndUnit* cmdWD = dynamic_cast<G4UIcmdWith3VectorAndUnit*> (cmd);
      fPrimaryAction->SetAngleForMuGeantinoScan( cmdWD->GetNew3VectorValue(val));   
  } else if (cmd ==  fMuGeantinoYOriginScan ) {
      G4UIcmdWith3VectorAndUnit* cmdWD = dynamic_cast<G4UIcmdWith3VectorAndUnit*> (cmd);
      fPrimaryAction->SetYOriginForMuGeantinoScan( cmdWD->GetNew3VectorValue(val));   
   } else if (cmd ==  fUseGeantino ) {
      if ((fPrimaryAction->GetUseMuonGeantino()) || (fPrimaryAction->GetUseChargedGeantino()) ) {
        G4Exception("LBNEPrimaryMessenger", "Inconsistency in particle choice ", FatalException,
                      "Can't use both a muon geantino, and a geantino ");
      }
      fPrimaryAction->SetUseGeantino(true);
   } else if (cmd ==  fUseMuonGeantino ) {
      if ((fPrimaryAction->GetUseGeantino()) || (fPrimaryAction->GetUseChargedGeantino())) {
        G4Exception("LBNEPrimaryMessenger", "Inconsistency in particle choice ", FatalException,
                      "Can't use both a muon (or charged) geantino, and a geantino ");
      }
      fPrimaryAction->SetUseMuonGeantino(true);
   } else if (cmd ==  fUseChargedGeantino ) {
      if ((fPrimaryAction->GetUseGeantino()) || (fPrimaryAction->GetUseMuonGeantino())) {
        G4Exception("LBNEPrimaryMessenger", "Inconsistency in particle choice ", FatalException,
                      "Can't use both a geantino or a muon geantino, and a Charged geantino ");
      }
      fPrimaryAction->SetUseChargedGeantino(true);
   } else if (cmd ==  fProtonMomentum ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      const double pz =  cmdWD->GetNewDoubleValue(val);
      if (std::abs(pz - 120.0) > 0.1) {
        std::string msg("With the current 1.2 MW option, and Proton Momentum != 120 GeV, \n");
        msg += std::string("The calculation of the beam spot size is unreliable.  \n");
        msg += std::string("The user must set beam signa X and sigma Y explicitely via cmd /LBNE/generator/beamSigmaX (Y) ");
        G4Exception("LBNEPrimaryMessenger", "", JustWarning, msg.c_str());
      }
      // 
      //Valid only for 1.2 MW option !!! 
      //
      fPrimaryAction->SetProtonMomentum( pz);
      fPrimaryAction->SetUseCourantSniderParams(false);
      fPrimaryAction->SetProtonBeam();
      
//      fPrimaryAction->SetProtonBeam();
//      double bf = (std::abs(pz - 120.0) < 0.1) ?
//                   110.9 : 
//                   fPrimaryAction->GetBetaFunctionvsBeamEnergy(pz);
//      fPrimaryAction->SetBeamBetaFunctionX(bf);
//      fPrimaryAction->SetBeamBetaFunctionY(bf);
//      fPrimaryAction->SetProtonBeam();
   }
   
   else if (cmd ==  fBeamBetaFunctionX ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetBeamBetaFunctionX(cmdWD->GetNewDoubleValue(val));
      fPrimaryAction->SetUseCourantSniderParams(true);
      fPrimaryAction->SetProtonBeam();
   }
   else if (cmd ==  fBeamBetaFunctionY ) {
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetBeamBetaFunctionY(cmdWD->GetNewDoubleValue(val));
      fPrimaryAction->SetUseCourantSniderParams(true);
      fPrimaryAction->SetProtonBeam();
   }
   else if (cmd ==  fBeamEmittanceX ) {      
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetBeamEmittanceX(cmdWD->GetNewDoubleValue(val));
      fPrimaryAction->SetUseCourantSniderParams(true);
      fPrimaryAction->SetProtonBeam();
   }
   else if (cmd ==  fBeamEmittanceY ) {      
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetBeamEmittanceY(cmdWD->GetNewDoubleValue(val));
      fPrimaryAction->SetUseCourantSniderParams(true);
      fPrimaryAction->SetProtonBeam();
   }
   else if (cmd == fInputFlukaFileName)
   {
      fPrimaryAction->OpenNtupleFLUKAASCII(val.c_str());
   }
   else if (cmd ==  fRadiusAnnularBeam ) {      
      G4UIcmdWithADoubleAndUnit* cmdWD = dynamic_cast<G4UIcmdWithADoubleAndUnit*> (cmd);
      fPrimaryAction->SetRadiusAnnularBeam(cmdWD->GetNewDoubleValue(val));
   }
   
}