LBNEPrimaryGeneratorAction.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
// $Id
//----------------------------------------------------------------------

#include <iostream>
#include <fstream>
#include <iomanip>

#include "LBNEPrimaryGeneratorAction.hh"

#include <math.h>
#include <map>
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ThreeVector.hh"
#include "globals.hh"
#include "Randomize.hh"
#include "G4UImanager.hh"

#include "LBNEPrimaryMessenger.hh"
#include "LBNERunManager.hh"
#include "LBNEParticleCode.hh"
#include "LBNEAnalysis.hh"
#include "LBNEVolumePlacements.hh"


#include "TROOT.h"
#include <TFile.h>
#include <TTree.h>
#include <TLeaf.h>
#include <TSystem.h>
#include "TH1.h"
#include "TRandom3.h"
#include "TF1.h"
#include "TMath.h"
#include "Math/IFunction.h"

using namespace std;

LBNEPrimaryGeneratorAction::LBNEPrimaryGeneratorAction() :
    fRunManager(0),
    fPrimaryMessenger(0),
    fParticleGun(0),

    fInputFile(0),
    fInputTree(0),
    lineForFluka(0), 

    fProtonN(0),
    fNoOfPrimaries(0),
    fCurrentPrimaryNo(0),

    //fTunnelPos(0),

    fProtonMomentumMag(120.0*CLHEP::GeV),
    fProtonOrigin(0),
    fProtonMomentum(0),
    fProtonIntVertex(0),
    fParticleMomentum(0),
    fParticlePosition(0),

    fTgen(-99),
    fType(-99),
    fMuGeantinoCnt(0),
    fWeight(1.0),
    fCorrectForAngle(true),
    fBeamOnTarget(false),
    fBeamOffsetX(0.0),
    fBeamOffsetY(0.0),
    fBeamOffsetZ(-3.6*CLHEP::m), // beam divergence assumed to be very small... 
                          // If old way of defining the beam!.. See below. 
                          // Should place well in front of the baffle, 
                          // hopefully not too much air in the way.. 
    fBeamMaxValX(1.0*CLHEP::m),  // No truncation of the beam by default. 
    fBeamMaxValY(1.0*CLHEP::m), // No truncation of the beam by default. 
    fBeamSigmaX(1.3*CLHEP::mm),
    fBeamSigmaY(1.3*CLHEP::mm),
    fBeamAngleTheta(0.0),
    fBeamAnglePhi(0.0),
    fUseCourantSniderParams(true),
    fUseJustSigmaCoord(false),
    fRadiusAnnularBeam(-1.0), // default is the usual 
    fSigmaForAnnularbeam(-1.0),
    fBeamEmittanceX(20.), // In pi mm mRad (Fermi units of emittance) 
    fBeamEmittanceY(20.),
    fBeamBetaFunctionX(64.842), // in meter, assuming a 120 GeV beam 
    fBeamBetaFunctionY(64.842), // in meter, assuming a 120 GeV beam // John Jonstone, e-mail, Oct 11
    fBeamBetaFunctionAt120(64.842), // in meter, assuming a 120 GeV beam For the 0.7 MW option!.... 
    fBeamBetaFunctionAt80(43.228), // in meter, assuming a 120 GeV beam 
    fBeamBetaFunctionAt60(32.421), // in meter, assuming a 120 GeV beam 
    
    fUseGeantino(false),
    fUseMuonGeantino(false),
    fUseChargedGeantino(false),
    fZOriginGeantino(-515.), // Upstream of target.
    fSigmaZOriginGeantino(100.), // Upstream of target.
    fPolarAngleGeantino(.005),
    fPolarAngleGeantinoMin(0.),
    fZOriginForMuGeantino(G4ThreeVector(0., 0., 0.)),
    fYOriginForMuGeantino(G4ThreeVector(0., 0., 0.)),
    fMomentumForMuGeantino(G4ThreeVector(7.5, 0., 0.)),
    fAngleForMuGeantino(G4ThreeVector(0.01, 0., 0.))
{
   fRunManager       =(LBNERunManager*)LBNERunManager::GetRunManager();
   fPrimaryMessenger = new LBNEPrimaryMessenger(this);
   G4int n_particle = 1;
   fParticleGun = new G4ParticleGun(n_particle);
   //
   // New default situation: for the 1.2 MW option.. 
   //
    fBeamSigmaX = 1.7*CLHEP::mm; fBeamSigmaY = 1.7*CLHEP::mm; fBeamBetaFunctionX = 110.9; fBeamBetaFunctionY = 110.9;
   //
   // Test the estimate fo the beta function 
   //
//   std::cerr << " Pz BetaFunc " << std::endl;
//   for (int iPz=0; iPz !=100; iPz++) 
//     std::cerr << " " << (50.0 + iPz*1.0) << " " <<
//               GetBetaFunctionvsBeamEnergy((50.0 + iPz*1.0)) << std:: endl;
//             
   // Other CourantSnider functions set down below. 
}
//---------------------------------------------------------------------------------------
LBNEPrimaryGeneratorAction::~LBNEPrimaryGeneratorAction()
{
  delete fPrimaryMessenger;
  delete fParticleGun;
  if (lineForFluka != 0) free(lineForFluka);
}
//---------------------------------------------------------------------------------------
void LBNEPrimaryGeneratorAction::SetProtonBeam()
{

   // Make sure we have a consistent way to define the beam 

   if (fUseCourantSniderParams && fUseJustSigmaCoord) {
     G4Exception("LBNEPrimaryGeneratorAction::SetProtonBeam", " ", 
                FatalErrorInArgument,
    "Attempting to refdefined both CourantSnider parameters and beam spot size.  Please choose one!..."  );
   } 
   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
 
  if (fUseGeantino) {
    fParticleGun->SetParticleDefinition(particleTable->FindParticle("geantino"));
  }
  else if (fUseMuonGeantino) {
    fParticleGun->SetParticleDefinition(particleTable->FindParticle("mu+"));
  }
  else if (fUseChargedGeantino) {
    fParticleGun->SetParticleDefinition(particleTable->FindParticle("chargedgeantino"));
  }
  else fParticleGun->SetParticleDefinition(particleTable->FindParticle("proton"));
  
  const double eBeam = std::sqrt(fProtonMomentumMag*fProtonMomentumMag + 0.938272*CLHEP::GeV*0.938272*CLHEP::GeV); 
  const double pEkin = eBeam - 0.938272*CLHEP::GeV;
  fParticleGun->SetParticleEnergy(pEkin);
  // This needs to be implemented via data cards!!!
  G4ThreeVector beamPosition(0., 0., fBeamOffsetZ);
  
  fParticleGun->SetParticlePosition(beamPosition);
  G4ThreeVector beamDirection(1,0,0);
  beamDirection.setTheta(fBeamAngleTheta);
  beamDirection.setPhi(fBeamAnglePhi);
  G4cout << "Beam direction of " << beamDirection[0] << " " <<
  beamDirection[1] << " " << beamDirection[2] << endl;
  fParticleGun->SetParticleMomentumDirection(beamDirection);
  
  fCurrentPrimaryNo=0;
  
//  LBNEVolumePlacements* volP=LBNEVolumePlacements::Instance();
//  bool use1p2MW = volP->GetUse1p2MW();
  // Upgrade the beam sigma April 17: also the beam emittance!
//  if (use1p2MW) { fBeamSigmaX = 1.7*CLHEP::mm; fBeamSigmaY = 1.7*CLHEP::mm; fBeamBetaFunctionX = 110.9; fBeamBetaFunctionY = 110.9;}


  G4String spaces = "   ";
  std::cout << spaces << "Configuring the Proton beam..." << std::endl
            << spaces << "   Momentum          = " << fParticleGun->GetParticleMomentum()/CLHEP::GeV << " GeV/c" << std::endl
            << spaces << "   Kinetic Energy    = " << fParticleGun->GetParticleEnergy()/CLHEP::GeV << " GeV" << std::endl
            << spaces << "   Beam Offset, X    = " << fBeamOffsetX/CLHEP::mm << " mm" << std::endl
            << spaces << "   Beam Offset, Y    = " << fBeamOffsetY/CLHEP::mm << " mm" << std::endl
            << spaces << "   Beam Offset, Z    = " << fBeamOffsetZ/CLHEP::m << " m" << std::endl
            << spaces << "   Beam Sigma, X     = " << fBeamSigmaX/CLHEP::mm << " mm" << std::endl
            << spaces << "   Beam Sigma, Y     = " << fBeamSigmaY/CLHEP::mm << " mm" << std::endl
            << spaces << "   Beam Max. d|X|    = " << fBeamMaxValX/CLHEP::mm << " mm" << std::endl
            << spaces << "   Beam Max, d|Y|    = " << fBeamMaxValY/CLHEP::mm << " mm" << std::endl
            << spaces << "   Direction         = " << fParticleGun->GetParticleMomentumDirection() << std::endl;

  if (fUseCourantSniderParams && (!fUseMuonGeantino) && (!fUseChargedGeantino)) {
     
     const double gammaLorentz = eBeam/(0.938272*CLHEP::GeV);
     const double betaLorentz = std::sqrt(1.0 - 1.0/(gammaLorentz*gammaLorentz));
     const double betaGammaLorentz = gammaLorentz*betaLorentz;
     
     const double betaFuncStarX = (std::abs(fProtonMomentumMag - 120.*CLHEP::GeV) < 0.1) ? 
                              fBeamBetaFunctionX : 
                              GetBetaFunctionvsBeamEnergy(fProtonMomentumMag/CLHEP::GeV); // in meters, not mm 
     const double betaFuncStarY = (std::abs(fProtonMomentumMag - 120.*CLHEP::GeV) < 0.1) ? 
                              fBeamBetaFunctionY : 
                              GetBetaFunctionvsBeamEnergy(fProtonMomentumMag/CLHEP::GeV);
                                 
     fBeamAlphaFunctionX = -(fBeamOffsetZ/CLHEP::m)/betaFuncStarX;
     fBeamAlphaFunctionY = -(fBeamOffsetZ/CLHEP::m)/betaFuncStarY;
     
     fBeamBetaFuncGenX = betaFuncStarX + ((fBeamOffsetZ/CLHEP::m)*(fBeamOffsetZ/CLHEP::m))/betaFuncStarX;// in meters
     fBeamBetaFuncGenY = betaFuncStarY + ((fBeamOffsetZ/CLHEP::m)*(fBeamOffsetZ/CLHEP::m))/betaFuncStarY;
     
     fBeamSigmaX = std::sqrt( 1.0e3 * fBeamBetaFuncGenX/CLHEP::m * fBeamEmittanceX/ (6.0*betaGammaLorentz)); // in mm 
     fBeamSigmaY = std::sqrt(  1.0e3 * fBeamBetaFuncGenY/CLHEP::m * fBeamEmittanceY/ (6.0*betaGammaLorentz));
     std::cout << spaces << "Re-Configuring the Proton beam using Courant-Snider parameters " << std::endl;
     std::cout << spaces << "  Beta function X at Point of creation " << fBeamBetaFuncGenX << " m " << std::endl;
     std::cout << spaces << "  Beta function Y at Point of creation " << fBeamBetaFuncGenY << " m " << std::endl;
     std::cout << spaces << "  Alpha function X at Point of creation " << fBeamAlphaFunctionX << std::endl;
     std::cout << spaces << "  Alpha function Y at Point of creation " << fBeamAlphaFunctionY << std::endl;
     std::cout << spaces << "  Sigma X at Point of creation " << fBeamSigmaX << " mm " << std::endl;
     std::cout << spaces << "  Sigma Y at Point of creation " << fBeamSigmaY << " mm " << std::endl;
     
  }
  if (fRadiusAnnularBeam > 1.0e-6) {
    // We will be using some of the ROOT math utilities, namely generating 
    // random number from a histogram. But we need to control the seed....
    // The following statement may cause havoc in other part of the code where 
    // The root Random generator is called.. 
    //
    int startSeed = CLHEP::HepRandom::getTheSeed();
    startSeed += 1267; // a bit arbitrary. 
    gRandom = new TRandom3(0);
    gRandom->SetSeed(startSeed); 
    fSigmaForAnnularbeam = std::sqrt(fBeamSigmaX*fBeamSigmaX + fBeamSigmaY*fBeamSigmaY);
    std::cout << spaces << " But.. but we will be using an annular beam, defined as follow " << std::endl;
    std::cout << spaces << " f(r; R, sigma) = d^2N/drdTheta =  " << std::endl; 
    std::cout << spaces << " = 1/(2*pi*sigma^2) * exp[(-r^2-R^2)/(2*sigma^2)] * I0(rR/sigma^2) " << std::endl;
    std::cout << spaces << " with R = " << fRadiusAnnularBeam << " and sigma = " 
                        << fSigmaForAnnularbeam << std::endl;
    std::cout << spaces << " This sigma is derived from Sigma X and Sigma Y, stated above  " << std::endl;
    // Prepare the distirbution.    
    theModBesselFunc1rstKind = new TF1("I_0", "TMath::BesselI0(x)", 0., 10.);
    histAnnularRadial = new TH1F("AN1", "Annular beam, radial", 1000, 0., 6.*fRadiusAnnularBeam);
    const double fSigSq = fSigmaForAnnularbeam*fSigmaForAnnularbeam;
    const double RSq = fRadiusAnnularBeam*fRadiusAnnularBeam;
    const double stepR = 6.*fRadiusAnnularBeam/1000.;
    for (int kR=0.; kR != 1000; kR++) {
      const double r = 0.49*stepR + kR*stepR;
      const double f = (1.0/(2.0*M_PI*fSigSq)) * std::exp((-r*r -RSq)/(2.*fSigSq)) * 
                            theModBesselFunc1rstKind->Eval(r*fRadiusAnnularBeam/fSigSq);
      histAnnularRadial->Fill(r, r*f);
    }
  }
  
}
G4bool LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII(G4String ntupleName)
{
   if (fFlukaASCII.is_open())  fFlukaASCII.close();
   fFlukaASCII.open(ntupleName.c_str());
   if (!fFlukaASCII.is_open()) {
      G4Exception("LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII", " ", 
                FatalErrorInArgument, " Fluka Input file is not open, unexpected and fatal " );
    }
  maxLengthForFlukaLine = 1000;
  lineForFluka = static_cast<char*>(malloc(maxLengthForFlukaLine));  
  fFlukaASCII.getline(lineForFluka, maxLengthForFlukaLine);
  if (!fFlukaASCII.good()) {
   G4Exception("LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII", " ", 
                FatalErrorInArgument, " Fluka Input file no longer accessible " );
  }
  if (fFlukaASCII.eof()) {
   fFlukaASCII.close();
   G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ", 
                FatalErrorInArgument, " Fluka Input file is empty....  " );
  }
  //
  // The first line if the Header for R.. 
  // 
  return true;

}
//---------------------------------------------------------------------------------------
G4bool LBNEPrimaryGeneratorAction::OpenNtuple(G4String ntupleName)
{
   G4String message("Input Ntuple"); 
   message += ntupleName + G4String("Not yet supported"); 
   G4Exception("LBNEPrimaryGeneratorAction::OpenNtuple", " ", 
                FatalErrorInArgument, message.c_str());
    /*          
   fCurrentPrimaryNo=0;
   
   G4bool fIsOpen=false;
   fInputFile=new TFile(ntupleName,"READ");
   if (!fInputFile->IsZombie())
   {
      
      //fInputTree=(TTree*)(fInputFile->Get("h3"));
      fInputTree=(TTree*)(fInputFile->Get((fLBNEData->GetInputNtpTreeName()).c_str()));
   
      if(!fInputTree)
      {
         //G4cout<<"LBNEPrimaryGeneratorAction: Can't find tree "<< G4endl;
         G4cout<<"LBNEPrimaryGeneratorAction: Can't find tree with name "
               << fLBNEData->GetInputNtpTreeName() << G4endl;
      }

      fCurrentPrimaryNo=0;
      G4int entries = fInputTree->GetEntries();
      if(fLBNEData->GetNEvents() > 0 && fLBNEData->GetNEvents() < entries)
         fNoOfPrimaries = fLBNEData->GetNEvents();
      else   
         fNoOfPrimaries = fInputTree->GetEntries();
      
      
     fIsOpen=true;
   }
   else
   {
      G4cout<<"LBNEPrimaryGeneratorAction: Input (fluka/mars) root file doesn't exist"
            <<"   Aborting run"<<G4endl;
      fIsOpen=false;
   }
   return fIsOpen;
   */
   return false;
}
//---------------------------------------------------------------------------------------
void LBNEPrimaryGeneratorAction::CloseNtuple()
{
   if (fFlukaASCII.is_open()) fFlukaASCII.close();
   if(!fInputFile) return;

   if(fInputFile && fInputFile -> IsOpen())
   {
      fInputFile->Close();
      if(fInputFile->IsOpen())
      {
         std::cout << "   PROBLEM: Failed to close Input Ntuple " << fInputFile -> GetName() << std::endl;
      }
      else 
      {
         std::cout << "   Sucessfully closed Input Ntuple " << fInputFile -> GetName() << std::endl;
      }
   }

   std::cout << "LBNEPrimaryGeneratorAction::CloseNtuple() - " << std::endl;
   std::cout << "   Used " << fProtonN << " protons from input file" << endl;




   fCurrentPrimaryNo=0;
}

double LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy(double pz){
 
                         // quadratic interpolation, based on the e-mail from 
                                     // John Jostone. 
        G4Exception("LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy", "", FatalException, 
          " With the current 1.2 MW option, the use of Courand Snider params and a beam energy != 120 GeV is prohibited  "); 

 // Define and slove a set of quadratic equation. 
 const double betaSecond80 = (fBeamBetaFunctionAt120 + 2.0*fBeamBetaFunctionAt60 - 
                                   3.0*fBeamBetaFunctionAt80)/2400.;
 const double betaPrime80 = (1.0/40.)*(fBeamBetaFunctionAt120 - fBeamBetaFunctionAt80 - betaSecond80*1600.); 
 
 const double dp = pz - 80.;
 std::cerr << " Using LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy... fBeamBetaFunctionAt120 " <<
           fBeamBetaFunctionAt120  << " At 80 " << fBeamBetaFunctionAt80 
                                  << " fBeamBetaFunctionAt60 " << fBeamBetaFunctionAt60 << std::endl;
 return (fBeamBetaFunctionAt80 + betaPrime80*dp + betaSecond80*dp*dp); 
}
//---------------------------------------------------------------------------------------
void LBNEPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
   int verboseLevel = G4EventManager::GetEventManager()->GetVerboseLevel();
   if(verboseLevel > 1)
   {
      std::cout << "Event " << anEvent->GetEventID() << ": LBNEPrimaryGeneratorAction::GeneratePrimaries() Called." << std::endl;
   }
   
   
   G4int totNoPrim=fRunManager->GetNumberOfEvents();
   if (totNoPrim>20)
   {
      if (fCurrentPrimaryNo%(totNoPrim/20)==0)
         G4cout<<"Processing particles #: "
               <<fCurrentPrimaryNo<<" to "<< fCurrentPrimaryNo+(totNoPrim/20) <<G4endl;
   }
   
   {
/* 
      if (fLBNEData->GetUseFlukaInput() || fLBNEData->GetUseMarsInput()) 
      {
         LBNEPrimaryGeneratorAction::GenerateBeamFromInput(anEvent);
      }
*/      
      if (fUseGeantino || fUseMuonGeantino || fUseChargedGeantino) {
        LBNEPrimaryGeneratorAction::Geantino(anEvent);
      }
      else if (fFlukaASCII.is_open()) {
         LBNEPrimaryGeneratorAction::GenerateFlukaPrimary(anEvent); 
      } else {
         LBNEPrimaryGeneratorAction::GenerateG4ProtonBeam(anEvent);
      }
      
   }   
   
   //fCurrentPrimaryNo++;
   ++fCurrentPrimaryNo;
}



//---------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------

//---------------------------------------------------------------------------------------
void LBNEPrimaryGeneratorAction::GenerateG4ProtonBeam(G4Event* anEvent)
{

// If nothing else is set, use a proton beam
    G4double x0 = 0.;
    G4double y0 = 0.; 
    G4double z0 = 0.;
    
    do 
    {
       x0 = G4RandGauss::shoot(0.0, fBeamSigmaX);
    }
    while(std::abs(x0) > fBeamMaxValX);
    do 
    {
       y0 = G4RandGauss::shoot(0.0, fBeamSigmaY);
    }
    while(std::abs(y0) > fBeamMaxValY);
    
    const double xFromEmitt = x0;
    const double yFromEmitt = y0;
    z0 = fBeamOffsetZ;

    if(!fBeamOnTarget){
      // b.c. if fBeamOnTarget, all offsets are ignored.
      x0 += fBeamOffsetX;
      y0 += fBeamOffsetY;
    }

    G4double dx=0.;
    G4double dy=0.; 
    G4double dz=0.;
    
    if(fabs(fBeamAngleTheta) > 1e-6){
      dx += sin(fBeamAngleTheta)*cos(fBeamAnglePhi);
      dy += sin(fBeamAngleTheta)*sin(fBeamAnglePhi);
    }
    dz = std::sqrt(1.0 - (dx*dx + dy*dy));
    G4ThreeVector direction(dx, dy, dz);
        
    if(fCorrectForAngle || fBeamOnTarget){
       x0 += (dx/dz)*z0;
       y0 += (dy/dz)*z0;
    }
    double phaseX = 0.;
    double phaseY = 0.;
    if (fUseCourantSniderParams) { // The tails in the dx, dy distribution not quite Gaussian, 
                                   // but this is small error. 
        phaseX = 2.0*M_PI*G4RandFlat::shoot();
        const double amplX = xFromEmitt/std::cos(phaseX);
        dx += (std::sin(phaseX)*amplX - fBeamAlphaFunctionX*xFromEmitt)/(fBeamBetaFuncGenX*CLHEP::m); // in radian (mm/mm) 
        phaseY = 2.0*M_PI*G4RandFlat::shoot();
        const double amplY = yFromEmitt/std::cos(phaseY);
        dy += (std::sin(phaseY)*amplY - fBeamAlphaFunctionY*yFromEmitt)/(fBeamBetaFuncGenY*CLHEP::m);
        dz = std::sqrt(1.0 - (dx*dx + dy*dy));
        direction = G4ThreeVector(dx, dy, dz);
    }
    if (fRadiusAnnularBeam > 1.0e-6) {
      const double r = histAnnularRadial->GetRandom();
      const double phi = 2.0*M_PI*G4RandFlat::shoot();
      x0 = r*std::cos(phi);
      y0 = r*std::sin(phi);
    }

//    if (std::abs(y0) > 5.0) {
//      std::cerr << " !!!!!.... Anomalously large y0 " << y0 << " dy " << dy << " Phase y " << phaseY << std::endl;
//    }

    fProtonN = fCurrentPrimaryNo;
    fProtonOrigin   = G4ThreeVector(x0, y0, z0);
    fProtonMomentum = G4ThreeVector(dx*fProtonMomentumMag, dy*fProtonMomentumMag, dz*fProtonMomentumMag);
    fProtonIntVertex = G4ThreeVector(-9999.,-9999.,-9999.);
    fWeight=1.; //for primary protons set weight and tgen
    fTgen=0;
    
    fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
    fParticleGun->SetParticleMomentumDirection(direction);
    //
    // Backdoor to create lots of beam particle, records data on where we are on the target.. 
    //
    // G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
    // fParticleGun->SetParticleDefinition(particleTable->FindParticle("geantino"));
    fParticleGun->GeneratePrimaryVertex(anEvent);
}
//
// Used to begug geometry, and/or, Radiography of the detector. 
//
void LBNEPrimaryGeneratorAction::Geantino(G4Event* anEvent)
{

   
    
// If nothing else is set, use a proton beam
    
    const double dr = fPolarAngleGeantinoMin + 
                      (fPolarAngleGeantino - fPolarAngleGeantinoMin)*G4RandFlat::shoot();
// Back door to study tile of Horn1 via magentic effect. Take a point source at = 0.
//     double dPhi  = M_PI/2.;
//     while (((dPhi >  M_PI/4.) && (dPhi <  3.0*M_PI/4.)) ||
//            ((dPhi > (M_PI + M_PI/4.)) && (dPhi < (2.0*M_PI - M_PI/4.))))
//        dPhi = 2.0*M_PI*G4RandFlat::shoot();
        
    // Back door to test limits in case of Helium tube misalignments. 
    const double dPhi = 2.0*M_PI*G4RandFlat::shoot();
//      const double dPhi = 0.;
//    const double dPhi = -M_PI/2. + 0.07893; // Shooting in a fixed direction.. Temporary 
//    const double dPhi = -M_PI/2.; // Shooting in a fixed direction.. Temporary 
    double dx = dr*std::cos(dPhi);
//    if (dx < 0.) dx *=-1.0; // Studying the right side (x > 0.);
    double dy = dr*std::sin(dPhi);
    double dz = sqrt(1.0 - (dx*dx + dy*dy));
    G4ThreeVector direction(dx, dy, dz);
//    std::cerr << " ... dr  " << dr << " direction " << direction << std::endl;
//
// Store it a proton internally to lbne d/s  ... Messy: 
    fProtonN = fCurrentPrimaryNo;
    
    fTgen=0;
//    const double x =  G4RandGauss::shoot(0.0, 1.3);  
//    const double y =  G4RandGauss::shoot(0.0, 1.3);  
// Back door to study tile of Horn1 via magentic effect. Take a point source at = 0.
    double x =  G4RandGauss::shoot(0.0, 1.3e-10);  
    double y =  G4RandGauss::shoot(0.0, 1.3e-10);
//    x  = 1.5; y = 0.; // off center, to study the height     
// Back door to study effect of overlap
//    const double dPhiPos = 2.0*M_PI*G4RandFlat::shoot();
//    const double x =  G4RandGauss::shoot(0.0, 0.00001) + 26.0 * std::sin(dPhiPos);  
//    const double y =  G4RandGauss::shoot(0.0, 0.00001) + 26.0 * std::cos(dPhiPos);
//    const double x =  0.040*CLHEP::mm;
    double z =  fZOriginGeantino +  G4RandGauss::shoot(0.0, fSigmaZOriginGeantino);
    //
    // November 2014: scan in either Z, Y, angle or momentum. 
    //
    bool doScanInZPos = (fZOriginForMuGeantino[2] > 0.001);  
    bool doScanInYPos = (fYOriginForMuGeantino[2] > 0.001);  
    bool doScanInMomentum = (fMomentumForMuGeantino[2] > 0.001);  
    bool doScanInAngle = (fAngleForMuGeantino[2] > 0.001);
    if (doScanInZPos || doScanInYPos ||  doScanInMomentum ||  doScanInAngle) {
      const int evtId = fMuGeantinoCnt;
      direction[0] = 0.;
      dy = fAngleForMuGeantino[0];      
      if (doScanInAngle) dy += (evtId)*fAngleForMuGeantino[1];
      dz = std::sqrt(1.0 - dy*dy);
      direction = G4ThreeVector(0., dy, dz);
      x = 0.;
      y = fYOriginForMuGeantino[0];
      if (doScanInYPos) y += (evtId)*fYOriginForMuGeantino[1];
      z = fZOriginForMuGeantino[0];
      if (doScanInZPos) z += (evtId)*fZOriginForMuGeantino[1];
      double pMu = fMomentumForMuGeantino[0];
      if (doScanInMomentum) pMu += (evtId)*fMomentumForMuGeantino[1];
      fProtonMomentumMag = std::sqrt(pMu*pMu + (105.658*105.658*CLHEP::MeV*CLHEP::MeV)); // misnomer.. 
      if (fPtForMuGeantino[2] > 1.) { // Per request from Jim H. and Alberto, scan either in Z position of momentum, 
                                      // but at a fixed Pt 
        const double Pt = fPtForMuGeantino[0];
        const double pTotal = std::sqrt(pMu*pMu + Pt*Pt); // pMu is Pz, here..
        direction = G4ThreeVector(0., Pt/pTotal, pMu/pTotal);
        fProtonMomentumMag = std::sqrt(pTotal*pTotal + (105.658*105.658*CLHEP::MeV*CLHEP::MeV)); // misnomer.. 
      }
       fMuGeantinoCnt++;
//      std::cerr << " fAngleForMuGeantino .. " << fAngleForMuGeantino << " dy " << dy << std::endl;
    }
//    std::cerr << " Z Particle position.. .. " <<  z << std::endl;
    fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z));
    fParticleGun->SetParticleMomentumDirection(direction);
//    std::cerr << " Muon Vertex " << G4ThreeVector(x, y, z) << " direction " << direction << std::endl;
    if (fUseMuonGeantino || fUseChargedGeantino) {
         if (fUseChargedGeantino) 
           fParticleGun->SetParticleEnergy(fProtonMomentumMag);
         else 
           fParticleGun->SetParticleEnergy(fProtonMomentumMag -  (105.658*CLHEP::MeV));
//       std::cerr << " I am indeed calling for muon .. and quit " << std::endl; exit(2);
    }
     // back door use of the proton momentum data card.,
                            // or see above, for scanning on pMu..This is the kinetic energy of the gun, not the total energy. . 
    fParticleGun->GeneratePrimaryVertex(anEvent);
}



//---------------------------------------------------------------------------------------
// 
// Old interface... 
// 
void LBNEPrimaryGeneratorAction::GenerateBeamFromInput(G4Event* anEvent)
{
    /*
     Fluka and Mars input variables:
     FLUKA                           MARS                   
     -----------------------------------------------------------------------------------------------
     TYPE                            TYPE                                      - type of particle (see LBNEAnalysis::GetParticleName())
     X, Y, Z                         X,Y,Z                                     - particle coordinates
     PX, PY, PZ                      PX, PY, PZ                                - particle momentum
     WEIGHT                          WEIGHT                                    - particle weight
     GENER                           GENER                                     - particle generation
     PROTVX, PROTVY, PROTVZ          PVTXX, PVTXY, PVTXZ                       - proton interaction vertex 
     PROTX, PROTY, PROTZ             PROTX, PROTY, PROTZ                       - proton initial coordinates
     PROTPX, PROTPY, PROTPZ          PROTPX, PROTPY, PROTPZ                    - proton initial momentum
     MOMPX,MOMPY,MOMPZ               PPX, PPY, PPZ                             - ???
     MOMTYPE                         PTYPE                                     - ???
     */

   G4Exception("LBNEPrimaryGeneratorAction::GenerateBeamFromInput", " ", 
                FatalErrorInArgument, " Input Ntuple not yet supported ");
      //
    //Need to create a new Gun each time
    //so Geant v4.9 doesn't complain
    //about momentum not match KE
    //
     if(fParticleGun){ delete fParticleGun; fParticleGun = 0;}
     fParticleGun = new G4ParticleGun(1);
     fParticleGun->SetParticleEnergy(0.0*CLHEP::GeV);
    
    G4double x0,y0,z0,px,py,pz;
    G4String particleName;
    fInputTree->GetEntry(fCurrentPrimaryNo);

    fProtonN = fInputTree->GetLeaf("event")->GetValue();
    
    x0 = fInputTree->GetLeaf("x")->GetValue()*CLHEP::cm;
    y0 = fInputTree->GetLeaf("y")->GetValue()*CLHEP::cm;
    z0 = fBeamOffsetZ; // fixed translation!!! Input file not yet supported anyways... 
//    z0 = fInputTree->GetLeaf("z")->GetValue()*CLHEP::cm+fLBNEData->GetTargetZ0(0)+fLBNEData->GetExtraFlukaNumiTargetZShift();
    //z0 = fInputTree->GetLeaf("z")->GetValue()*CLHEP::cm+fLBNEData->TargetZ0+fLBNEData->GetExtraFlukaNumiTargetZShift();
    
    px = fInputTree->GetLeaf("px")->GetValue()*CLHEP::GeV;
    py = fInputTree->GetLeaf("py")->GetValue()*CLHEP::GeV;
    pz = fInputTree->GetLeaf("pz")->GetValue()*CLHEP::GeV;
    
    fParticlePosition=G4ThreeVector(x0,y0,z0);
    fParticleMomentum=G4ThreeVector(px,py,pz);
    
    fWeight = fInputTree->GetLeaf("weight")->GetValue();
    fType = G4int(fInputTree->GetLeaf("type")->GetValue());
    fTgen = G4int(fInputTree->GetLeaf("gener")->GetValue());
    particleName=LBNEParticleCode::AsString(LBNEParticleCode::IntToEnum(fType));
    fProtonOrigin   = G4ThreeVector(fInputTree->GetLeaf("protx")->GetValue()*CLHEP::cm,
                                    fInputTree->GetLeaf("proty")->GetValue()*CLHEP::cm,
                                    fInputTree->GetLeaf("protz")->GetValue()*CLHEP::cm);
    fProtonMomentum = G4ThreeVector(fInputTree->GetLeaf("protpx")->GetValue()*CLHEP::cm,
                                    fInputTree->GetLeaf("protpy")->GetValue()*CLHEP::cm,
                                    fInputTree->GetLeaf("protpz")->GetValue()*CLHEP::cm);
    
/*
    if (fLBNEData->GetUseMarsInput()){
      fProtonIntVertex = G4ThreeVector(fInputTree->GetLeaf("pvtxx")->GetValue()*CLHEP::cm,
                                       fInputTree->GetLeaf("pvtxy")->GetValue()*CLHEP::cm,
                                       fInputTree->GetLeaf("pvtxz")->GetValue()*CLHEP::cm);
    }
    else if (fLBNEData->GetUseFlukaInput()){
      fProtonIntVertex = G4ThreeVector(fInputTree->GetLeaf("protvx")->GetValue()*CLHEP::cm,
                                       fInputTree->GetLeaf("protvy")->GetValue()*CLHEP::cm,
                                       fInputTree->GetLeaf("protvz")->GetValue()*CLHEP::cm);
    }
*/    
    G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
    fParticleGun->SetParticleDefinition(particleTable->FindParticle(particleName));
    //G4double mass=particleTable->FindParticle(particleName)->GetPDGMass();
    
    //fParticleGun->SetParticleEnergy(sqrt(mass*CLHEP::mass+px*px+py*py+pz*pz)-mass);
    fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
    fParticleGun->SetParticleMomentum(G4ThreeVector(px,py,pz));
    fParticleGun->GeneratePrimaryVertex(anEvent);

}
//
// New Fluka interface 
//
void LBNEPrimaryGeneratorAction::GenerateFlukaPrimary(G4Event* anEvent) {
  if (!fFlukaASCII.is_open()) {
   G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ", 
                FatalErrorInArgument, " Fluka Input file is not open, unexpected and fatal " );
  }
  fFlukaASCII.getline(lineForFluka, maxLengthForFlukaLine);
  if (!fFlukaASCII.good()) {
   G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ", 
                RunMustBeAborted, " Fluka Input file no longer accessible " );
                return;
  }
  if (fFlukaASCII.eof()) {
   fFlukaASCII.close();
   G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ", 
                RunMustBeAborted, " Fluka Input file at EOF, no more events to generate " );
                return;
  }
  std::string llStr(lineForFluka);
  std::istringstream llStrStr(llStr);  
  G4ThreeVector beamPosition(0., 0., 0.);
  for (size_t k=0; k !=3; k++) { llStrStr >> beamPosition[k]; beamPosition[k] *= 1.0*CLHEP::cm; } 
  fParticleGun->SetParticlePosition(beamPosition);  
  G4ThreeVector partMomentum(0.,0.,0.);
  for (size_t k=0; k !=3; k++) { llStrStr >> partMomentum[k]; partMomentum[k] *= 1.0*CLHEP::GeV; }
  G4ThreeVector beamDirection(partMomentum);
  double pTotTmp = beamDirection.mag();
  for (size_t k=0; k !=3; k++) beamDirection[k] /= pTotTmp;
  fParticleGun->SetParticleMomentumDirection(beamDirection);   
  int particleID = 0;
  llStrStr >> particleID;
  G4ParticleTable *pTable=G4ParticleTable::GetParticleTable();
  G4ParticleDefinition* pDef = pTable->FindParticle(particleID);
  fParticleGun->SetParticleDefinition(pDef);
  const double pMass = pDef->GetPDGMass();
  const double eTot = std::sqrt(pMass*pMass + pTotTmp*pTotTmp);  
  fParticleGun->SetParticleEnergy(eTot - pMass);
  // Note: In the original version of fluka user routine fluscw.f, the weight was determined as: 
  //
  //       if ((totp.gt.30.) .or. (totp .lt. pmin1)) then  ! calculate weight
  //          mywei = 1.
  //       else
  //          mywei = 30./totp
  //       endif
  // Where totp is the total momentum of the particle. 
  // Not sure if this is still relevant..!
  //
  fWeight = 1.0; 
  fParticleGun->GeneratePrimaryVertex(anEvent);

  
}