Public Member Functions | Private Attributes | List of all members
LBNEPrimaryGeneratorAction Class Reference

#include <LBNEPrimaryGeneratorAction.hh>

Inheritance diagram for LBNEPrimaryGeneratorAction:

Public Member Functions

 LBNEPrimaryGeneratorAction ()
 
 ~LBNEPrimaryGeneratorAction ()
 
void GeneratePrimaries (G4Event *anEvent)
 
void GenerateG4ProtonBeam (G4Event *anEvent)
 
void GenerateBeamFromInput (G4Event *anEvent)
 
void GenerateFlukaPrimary (G4Event *anEvent)
 
void Geantino (G4Event *anEvent)
 
G4ParticleGun * GetParticleGun ()
 
void SetProtonBeam ()
 
G4bool OpenNtuple (G4String ntupleName)
 
G4bool OpenNtupleFLUKAASCII (G4String ntupleName)
 
void CloseNtuple ()
 
G4int GetNoOfPrimaries ()
 
G4ThreeVector GetProtonOrigin ()
 
G4ThreeVector GetProtonMomentum ()
 
G4ThreeVector GetProtonIntVertex ()
 
G4ThreeVector GetParticlePosition ()
 
G4ThreeVector GetParticleMomentum ()
 
G4double GetWeight ()
 
G4int GetTgen ()
 
G4int GetParticleType ()
 
G4int GetProtonNumber ()
 
void SetCorrectForAngle (G4bool aBool)
 
void SetBeamOnTarget (G4bool aBool)
 
void SetBeamOffsetX (G4double x)
 
void SetBeamOffsetY (G4double y)
 
void SetBeamSigmaX (G4double x)
 
void SetBeamSigmaY (G4double y)
 
void SetBeamMaxValX (G4double x)
 
void SetBeamMaxValY (G4double y)
 
void SetBeamTheta (G4double theta)
 
void SetBeamPhi (G4double phi)
 
void SetProtonMomentum (G4double p)
 
void SetUseCourantSniderParams (G4bool t)
 
void SetUseJustSigmaCoord (G4bool t)
 
G4bool GetUseCourantSniderParams ()
 
G4bool GetUseJustSigmaCoord ()
 
void SetBeamEmittanceX (G4double t)
 
void SetBeamEmittanceY (G4double t)
 
void SetBeamBetaFunctionX (G4double t)
 
void SetBeamBetaFunctionY (G4double t)
 
void SetRadiusAnnularBeam (G4double t)
 
double GetBeamSigmaX () const
 
double GetBeamSigmaY () const
 
double GetBeamOffsetX () const
 
double GetBeamOffsetY () const
 
double GetBeamOffsetZ () const
 
double GetBeamMaxValX () const
 
double GetBeamMaxValY () const
 
double GetBeamEmittanceX () const
 
double GetBeamBetaFunctionX () const
 
double GetBeamEmittanceY () const
 
double GetBeamBetaFunctionY () const
 
double GetBeamAngleTheta () const
 
void SetPolarAngleGeantino (double r)
 
void SetPolarAngleGeantinoMin (double r)
 
void SetUseGeantino (bool t)
 
void SetUseMuonGeantino (bool t)
 
void SetUseChargedGeantino (bool t)
 
void SetZOriginGeantino (double v)
 
void SetSigmaZOriginGeantino (double v)
 
bool GetUseMuonGeantino () const
 
bool GetUseGeantino () const
 
bool GetUseChargedGeantino () const
 
double GetBetaFunctionvsBeamEnergy (double pz)
 
void SetZOriginForMuGeantinoScan (G4ThreeVector v)
 
void SetYOriginForMuGeantinoScan (G4ThreeVector v)
 
void SetMomentumForMuGeantinoScan (G4ThreeVector v)
 
void SetPtForMuGeantinoScan (G4ThreeVector v)
 
void SetAngleForMuGeantinoScan (G4ThreeVector v)
 
G4ThreeVector GetZOriginForMuGeantinoScan () const
 
G4ThreeVector GetYOriginForMuGeantinoScan () const
 
G4ThreeVector GetMomentumForMuGeantinoScan () const
 
G4ThreeVector GetPtForMuGeantinoScan () const
 
G4ThreeVector GetAngleForMuGeantinoScan () const
 

Private Attributes

LBNERunManager * fRunManager
 
LBNEPrimaryMessenger * fPrimaryMessenger
 
G4ParticleGun * fParticleGun
 
TFile * fInputFile
 
TTree * fInputTree
 
std::ifstream fFlukaASCII
 
size_t maxLengthForFlukaLine
 
char * lineForFluka
 
G4int fProtonN
 
G4int fNoOfPrimaries
 
G4int fCurrentPrimaryNo
 
G4double fProtonMomentumMag
 
G4ThreeVector fProtonOrigin
 
G4ThreeVector fProtonMomentum
 
G4ThreeVector fProtonIntVertex
 
G4ThreeVector fParticleMomentum
 
G4ThreeVector fParticlePosition
 
G4int fTgen
 
G4int fType
 
G4int fMuGeantinoCnt
 
G4double fWeight
 
G4bool fCorrectForAngle
 
G4bool fBeamOnTarget
 
G4double fBeamOffsetX
 
G4double fBeamOffsetY
 
G4double fBeamOffsetZ
 
G4double fBeamMaxValX
 
G4double fBeamMaxValY
 
G4double fBeamSigmaX
 
G4double fBeamSigmaY
 
G4double fBeamAngleTheta
 
G4double fBeamAnglePhi
 
G4bool fUseCourantSniderParams
 
G4bool fUseJustSigmaCoord
 
G4double fRadiusAnnularBeam
 
G4double fSigmaForAnnularbeam
 
G4double fBeamEmittanceX
 
G4double fBeamEmittanceY
 
G4double fBeamBetaFunctionX
 
G4double fBeamBetaFunctionY
 
G4double fBeamBetaFunctionAt120
 
G4double fBeamBetaFunctionAt80
 
G4double fBeamBetaFunctionAt60
 
G4double fBeamAlphaFunctionX
 
G4double fBeamAlphaFunctionY
 
G4double fBeamBetaFuncGenX
 
G4double fBeamBetaFuncGenY
 
G4bool fUseGeantino
 
G4bool fUseMuonGeantino
 
G4bool fUseChargedGeantino
 
G4double fZOriginGeantino
 
G4double fSigmaZOriginGeantino
 
G4double fPolarAngleGeantino
 
G4double fPolarAngleGeantinoMin
 
G4ThreeVector fZOriginForMuGeantino
 
G4ThreeVector fYOriginForMuGeantino
 
G4ThreeVector fMomentumForMuGeantino
 
G4ThreeVector fPtForMuGeantino
 
G4ThreeVector fAngleForMuGeantino
 
TF1 * theModBesselFunc1rstKind
 
TH1F * histAnnularRadial
 

Detailed Description

Definition at line 20 of file LBNEPrimaryGeneratorAction.hh.

Constructor & Destructor Documentation

LBNEPrimaryGeneratorAction::LBNEPrimaryGeneratorAction ( )

Definition at line 42 of file LBNEPrimaryGeneratorAction.cc.

42  :
43  fRunManager(0),
44  fPrimaryMessenger(0),
45  fParticleGun(0),
46 
47  fInputFile(0),
48  fInputTree(0),
49  lineForFluka(0),
50 
51  fProtonN(0),
52  fNoOfPrimaries(0),
53  fCurrentPrimaryNo(0),
54 
55  //fTunnelPos(0),
56 
57  fProtonMomentumMag(120.0*CLHEP::GeV),
58  fProtonOrigin(0),
59  fProtonMomentum(0),
60  fProtonIntVertex(0),
61  fParticleMomentum(0),
62  fParticlePosition(0),
63 
64  fTgen(-99),
65  fType(-99),
66  fMuGeantinoCnt(0),
67  fWeight(1.0),
68  fCorrectForAngle(true),
69  fBeamOnTarget(false),
70  fBeamOffsetX(0.0),
71  fBeamOffsetY(0.0),
72  fBeamOffsetZ(-3.6*CLHEP::m), // beam divergence assumed to be very small...
73  // If old way of defining the beam!.. See below.
74  // Should place well in front of the baffle,
75  // hopefully not too much air in the way..
76  fBeamMaxValX(1.0*CLHEP::m), // No truncation of the beam by default.
77  fBeamMaxValY(1.0*CLHEP::m), // No truncation of the beam by default.
78  fBeamSigmaX(1.3*CLHEP::mm),
79  fBeamSigmaY(1.3*CLHEP::mm),
80  fBeamAngleTheta(0.0),
81  fBeamAnglePhi(0.0),
82  fUseCourantSniderParams(true),
83  fUseJustSigmaCoord(false),
84  fRadiusAnnularBeam(-1.0), // default is the usual
85  fSigmaForAnnularbeam(-1.0),
86  fBeamEmittanceX(20.), // In pi mm mRad (Fermi units of emittance)
87  fBeamEmittanceY(20.),
88  fBeamBetaFunctionX(64.842), // in meter, assuming a 120 GeV beam
89  fBeamBetaFunctionY(64.842), // in meter, assuming a 120 GeV beam // John Jonstone, e-mail, Oct 11
90  fBeamBetaFunctionAt120(64.842), // in meter, assuming a 120 GeV beam For the 0.7 MW option!....
91  fBeamBetaFunctionAt80(43.228), // in meter, assuming a 120 GeV beam
92  fBeamBetaFunctionAt60(32.421), // in meter, assuming a 120 GeV beam
93 
94  fUseGeantino(false),
95  fUseMuonGeantino(false),
96  fUseChargedGeantino(false),
97  fZOriginGeantino(-515.), // Upstream of target.
98  fSigmaZOriginGeantino(100.), // Upstream of target.
99  fPolarAngleGeantino(.005),
100  fPolarAngleGeantinoMin(0.),
101  fZOriginForMuGeantino(G4ThreeVector(0., 0., 0.)),
102  fYOriginForMuGeantino(G4ThreeVector(0., 0., 0.)),
103  fMomentumForMuGeantino(G4ThreeVector(7.5, 0., 0.)),
104  fAngleForMuGeantino(G4ThreeVector(0.01, 0., 0.))
105 {
106  fRunManager =(LBNERunManager*)LBNERunManager::GetRunManager();
107  fPrimaryMessenger = new LBNEPrimaryMessenger(this);
108  G4int n_particle = 1;
109  fParticleGun = new G4ParticleGun(n_particle);
110  //
111  // New default situation: for the 1.2 MW option..
112  //
113  fBeamSigmaX = 1.7*CLHEP::mm; fBeamSigmaY = 1.7*CLHEP::mm; fBeamBetaFunctionX = 110.9; fBeamBetaFunctionY = 110.9;
114  //
115  // Test the estimate fo the beta function
116  //
117 // std::cerr << " Pz BetaFunc " << std::endl;
118 // for (int iPz=0; iPz !=100; iPz++)
119 // std::cerr << " " << (50.0 + iPz*1.0) << " " <<
120 // GetBetaFunctionvsBeamEnergy((50.0 + iPz*1.0)) << std:: endl;
121 //
122  // Other CourantSnider functions set down below.
123 }
run_cnn_1class.m
m
Definition: run_cnn_1class.py:75
LBNEPrimaryGeneratorAction::fPrimaryMessenger
LBNEPrimaryMessenger * fPrimaryMessenger
Definition: LBNEPrimaryGeneratorAction.hh:102
LBNEPrimaryGeneratorAction::~LBNEPrimaryGeneratorAction ( )

Definition at line 125 of file LBNEPrimaryGeneratorAction.cc.

126 {
127  delete fPrimaryMessenger;
128  delete fParticleGun;
129  if (lineForFluka != 0) free(lineForFluka);
130 }
LBNEPrimaryGeneratorAction::fPrimaryMessenger
LBNEPrimaryMessenger * fPrimaryMessenger
Definition: LBNEPrimaryGeneratorAction.hh:102

Member Function Documentation

void LBNEPrimaryGeneratorAction::CloseNtuple ( )

Definition at line 324 of file LBNEPrimaryGeneratorAction.cc.

325 {
326  if (fFlukaASCII.is_open()) fFlukaASCII.close();
327  if(!fInputFile) return;
328 
329  if(fInputFile && fInputFile -> IsOpen())
330  {
331  fInputFile->Close();
332  if(fInputFile->IsOpen())
333  {
334  std::cout << " PROBLEM: Failed to close Input Ntuple " << fInputFile -> GetName() << std::endl;
335  }
336  else
337  {
338  std::cout << " Sucessfully closed Input Ntuple " << fInputFile -> GetName() << std::endl;
339  }
340  }
341 
342  std::cout << "LBNEPrimaryGeneratorAction::CloseNtuple() - " << std::endl;
343  std::cout << " Used " << fProtonN << " protons from input file" << endl;
344 
345 
346 
347 
348  fCurrentPrimaryNo=0;
349 }
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void LBNEPrimaryGeneratorAction::Geantino ( G4Event *  anEvent)

Definition at line 504 of file LBNEPrimaryGeneratorAction.cc.

505 {
506 
507 
508 
509 // If nothing else is set, use a proton beam
510 
511  const double dr = fPolarAngleGeantinoMin +
512  (fPolarAngleGeantino - fPolarAngleGeantinoMin)*G4RandFlat::shoot();
513 // Back door to study tile of Horn1 via magentic effect. Take a point source at = 0.
514 // double dPhi = M_PI/2.;
515 // while (((dPhi > M_PI/4.) && (dPhi < 3.0*M_PI/4.)) ||
516 // ((dPhi > (M_PI + M_PI/4.)) && (dPhi < (2.0*M_PI - M_PI/4.))))
517 // dPhi = 2.0*M_PI*G4RandFlat::shoot();
518 
519  // Back door to test limits in case of Helium tube misalignments.
520  const double dPhi = 2.0*M_PI*G4RandFlat::shoot();
521 // const double dPhi = 0.;
522 // const double dPhi = -M_PI/2. + 0.07893; // Shooting in a fixed direction.. Temporary
523 // const double dPhi = -M_PI/2.; // Shooting in a fixed direction.. Temporary
524  double dx = dr*std::cos(dPhi);
525 // if (dx < 0.) dx *=-1.0; // Studying the right side (x > 0.);
526  double dy = dr*std::sin(dPhi);
527  double dz = sqrt(1.0 - (dx*dx + dy*dy));
528  G4ThreeVector direction(dx, dy, dz);
529 // std::cerr << " ... dr " << dr << " direction " << direction << std::endl;
530 //
531 // Store it a proton internally to lbne d/s ... Messy:
532  fProtonN = fCurrentPrimaryNo;
533 
534  fTgen=0;
535 // const double x = G4RandGauss::shoot(0.0, 1.3);
536 // const double y = G4RandGauss::shoot(0.0, 1.3);
537 // Back door to study tile of Horn1 via magentic effect. Take a point source at = 0.
538  double x = G4RandGauss::shoot(0.0, 1.3e-10);
539  double y = G4RandGauss::shoot(0.0, 1.3e-10);
540 // x = 1.5; y = 0.; // off center, to study the height
541 // Back door to study effect of overlap
542 // const double dPhiPos = 2.0*M_PI*G4RandFlat::shoot();
543 // const double x = G4RandGauss::shoot(0.0, 0.00001) + 26.0 * std::sin(dPhiPos);
544 // const double y = G4RandGauss::shoot(0.0, 0.00001) + 26.0 * std::cos(dPhiPos);
545 // const double x = 0.040*CLHEP::mm;
546  double z = fZOriginGeantino + G4RandGauss::shoot(0.0, fSigmaZOriginGeantino);
547  //
548  // November 2014: scan in either Z, Y, angle or momentum.
549  //
550  bool doScanInZPos = (fZOriginForMuGeantino[2] > 0.001);
551  bool doScanInYPos = (fYOriginForMuGeantino[2] > 0.001);
552  bool doScanInMomentum = (fMomentumForMuGeantino[2] > 0.001);
553  bool doScanInAngle = (fAngleForMuGeantino[2] > 0.001);
554  if (doScanInZPos || doScanInYPos || doScanInMomentum || doScanInAngle) {
555  const int evtId = fMuGeantinoCnt;
556  direction[0] = 0.;
557  dy = fAngleForMuGeantino[0];
558  if (doScanInAngle) dy += (evtId)*fAngleForMuGeantino[1];
559  dz = std::sqrt(1.0 - dy*dy);
560  direction = G4ThreeVector(0., dy, dz);
561  x = 0.;
562  y = fYOriginForMuGeantino[0];
563  if (doScanInYPos) y += (evtId)*fYOriginForMuGeantino[1];
564  z = fZOriginForMuGeantino[0];
565  if (doScanInZPos) z += (evtId)*fZOriginForMuGeantino[1];
566  double pMu = fMomentumForMuGeantino[0];
567  if (doScanInMomentum) pMu += (evtId)*fMomentumForMuGeantino[1];
568  fProtonMomentumMag = std::sqrt(pMu*pMu + (105.658*105.658*CLHEP::MeV*CLHEP::MeV)); // misnomer..
569  if (fPtForMuGeantino[2] > 1.) { // Per request from Jim H. and Alberto, scan either in Z position of momentum,
570  // but at a fixed Pt
571  const double Pt = fPtForMuGeantino[0];
572  const double pTotal = std::sqrt(pMu*pMu + Pt*Pt); // pMu is Pz, here..
573  direction = G4ThreeVector(0., Pt/pTotal, pMu/pTotal);
574  fProtonMomentumMag = std::sqrt(pTotal*pTotal + (105.658*105.658*CLHEP::MeV*CLHEP::MeV)); // misnomer..
575  }
576  fMuGeantinoCnt++;
577 // std::cerr << " fAngleForMuGeantino .. " << fAngleForMuGeantino << " dy " << dy << std::endl;
578  }
579 // std::cerr << " Z Particle position.. .. " << z << std::endl;
580  fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z));
581  fParticleGun->SetParticleMomentumDirection(direction);
582 // std::cerr << " Muon Vertex " << G4ThreeVector(x, y, z) << " direction " << direction << std::endl;
583  if (fUseMuonGeantino || fUseChargedGeantino) {
584  if (fUseChargedGeantino)
585  fParticleGun->SetParticleEnergy(fProtonMomentumMag);
586  else
587  fParticleGun->SetParticleEnergy(fProtonMomentumMag - (105.658*CLHEP::MeV));
588 // std::cerr << " I am indeed calling for muon .. and quit " << std::endl; exit(2);
589  }
590  // back door use of the proton momentum data card.,
591  // or see above, for scanning on pMu..This is the kinetic energy of the gun, not the total energy. .
592  fParticleGun->GeneratePrimaryVertex(anEvent);
593 }
y
double y
Definition: GapWidth_module.cc:109
z
double z
Definition: GapWidth_module.cc:108
train.x
list x
Definition: train.py:276
void LBNEPrimaryGeneratorAction::GenerateBeamFromInput ( G4Event *  anEvent)

Definition at line 601 of file LBNEPrimaryGeneratorAction.cc.

602 {
603  /*
604  Fluka and Mars input variables:
605  FLUKA MARS
606  -----------------------------------------------------------------------------------------------
607  TYPE TYPE - type of particle (see LBNEAnalysis::GetParticleName())
608  X, Y, Z X,Y,Z - particle coordinates
609  PX, PY, PZ PX, PY, PZ - particle momentum
610  WEIGHT WEIGHT - particle weight
611  GENER GENER - particle generation
612  PROTVX, PROTVY, PROTVZ PVTXX, PVTXY, PVTXZ - proton interaction vertex
613  PROTX, PROTY, PROTZ PROTX, PROTY, PROTZ - proton initial coordinates
614  PROTPX, PROTPY, PROTPZ PROTPX, PROTPY, PROTPZ - proton initial momentum
615  MOMPX,MOMPY,MOMPZ PPX, PPY, PPZ - ???
616  MOMTYPE PTYPE - ???
617  */
618 
619  G4Exception("LBNEPrimaryGeneratorAction::GenerateBeamFromInput", " ",
620  FatalErrorInArgument, " Input Ntuple not yet supported ");
621  //
622  //Need to create a new Gun each time
623  //so Geant v4.9 doesn't complain
624  //about momentum not match KE
625  //
626  if(fParticleGun){ delete fParticleGun; fParticleGun = 0;}
627  fParticleGun = new G4ParticleGun(1);
628  fParticleGun->SetParticleEnergy(0.0*CLHEP::GeV);
629 
630  G4double x0,y0,z0,px,py,pz;
631  G4String particleName;
632  fInputTree->GetEntry(fCurrentPrimaryNo);
633 
634  fProtonN = fInputTree->GetLeaf("event")->GetValue();
635 
636  x0 = fInputTree->GetLeaf("x")->GetValue()*CLHEP::cm;
637  y0 = fInputTree->GetLeaf("y")->GetValue()*CLHEP::cm;
638  z0 = fBeamOffsetZ; // fixed translation!!! Input file not yet supported anyways...
639 // z0 = fInputTree->GetLeaf("z")->GetValue()*CLHEP::cm+fLBNEData->GetTargetZ0(0)+fLBNEData->GetExtraFlukaNumiTargetZShift();
640  //z0 = fInputTree->GetLeaf("z")->GetValue()*CLHEP::cm+fLBNEData->TargetZ0+fLBNEData->GetExtraFlukaNumiTargetZShift();
641 
642  px = fInputTree->GetLeaf("px")->GetValue()*CLHEP::GeV;
643  py = fInputTree->GetLeaf("py")->GetValue()*CLHEP::GeV;
644  pz = fInputTree->GetLeaf("pz")->GetValue()*CLHEP::GeV;
645 
646  fParticlePosition=G4ThreeVector(x0,y0,z0);
647  fParticleMomentum=G4ThreeVector(px,py,pz);
648 
649  fWeight = fInputTree->GetLeaf("weight")->GetValue();
650  fType = G4int(fInputTree->GetLeaf("type")->GetValue());
651  fTgen = G4int(fInputTree->GetLeaf("gener")->GetValue());
652  particleName=LBNEParticleCode::AsString(LBNEParticleCode::IntToEnum(fType));
653  fProtonOrigin = G4ThreeVector(fInputTree->GetLeaf("protx")->GetValue()*CLHEP::cm,
654  fInputTree->GetLeaf("proty")->GetValue()*CLHEP::cm,
655  fInputTree->GetLeaf("protz")->GetValue()*CLHEP::cm);
656  fProtonMomentum = G4ThreeVector(fInputTree->GetLeaf("protpx")->GetValue()*CLHEP::cm,
657  fInputTree->GetLeaf("protpy")->GetValue()*CLHEP::cm,
658  fInputTree->GetLeaf("protpz")->GetValue()*CLHEP::cm);
659 
660 /*
661  if (fLBNEData->GetUseMarsInput()){
662  fProtonIntVertex = G4ThreeVector(fInputTree->GetLeaf("pvtxx")->GetValue()*CLHEP::cm,
663  fInputTree->GetLeaf("pvtxy")->GetValue()*CLHEP::cm,
664  fInputTree->GetLeaf("pvtxz")->GetValue()*CLHEP::cm);
665  }
666  else if (fLBNEData->GetUseFlukaInput()){
667  fProtonIntVertex = G4ThreeVector(fInputTree->GetLeaf("protvx")->GetValue()*CLHEP::cm,
668  fInputTree->GetLeaf("protvy")->GetValue()*CLHEP::cm,
669  fInputTree->GetLeaf("protvz")->GetValue()*CLHEP::cm);
670  }
671 */
672  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
673  fParticleGun->SetParticleDefinition(particleTable->FindParticle(particleName));
674  //G4double mass=particleTable->FindParticle(particleName)->GetPDGMass();
675 
676  //fParticleGun->SetParticleEnergy(sqrt(mass*CLHEP::mass+px*px+py*py+pz*pz)-mass);
677  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
678  fParticleGun->SetParticleMomentum(G4ThreeVector(px,py,pz));
679  fParticleGun->GeneratePrimaryVertex(anEvent);
680 
681 }
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
LBNEParticleCode::IntToEnum
LBNEParticleCode_t IntToEnum(G4int particleInt)
Definition: LBNEParticleCode.cc:135
LBNEParticleCode::AsString
G4String AsString(LBNEParticleCode_t pCode)
Definition: LBNEParticleCode.cc:10
void LBNEPrimaryGeneratorAction::GenerateFlukaPrimary ( G4Event *  anEvent)

Definition at line 685 of file LBNEPrimaryGeneratorAction.cc.

685  {
686  if (!fFlukaASCII.is_open()) {
687  G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ",
688  FatalErrorInArgument, " Fluka Input file is not open, unexpected and fatal " );
689  }
690  fFlukaASCII.getline(lineForFluka, maxLengthForFlukaLine);
691  if (!fFlukaASCII.good()) {
692  G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ",
693  RunMustBeAborted, " Fluka Input file no longer accessible " );
694  return;
695  }
696  if (fFlukaASCII.eof()) {
697  fFlukaASCII.close();
698  G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ",
699  RunMustBeAborted, " Fluka Input file at EOF, no more events to generate " );
700  return;
701  }
702  std::string llStr(lineForFluka);
703  std::istringstream llStrStr(llStr);
704  G4ThreeVector beamPosition(0., 0., 0.);
705  for (size_t k=0; k !=3; k++) { llStrStr >> beamPosition[k]; beamPosition[k] *= 1.0*CLHEP::cm; }
706  fParticleGun->SetParticlePosition(beamPosition);
707  G4ThreeVector partMomentum(0.,0.,0.);
708  for (size_t k=0; k !=3; k++) { llStrStr >> partMomentum[k]; partMomentum[k] *= 1.0*CLHEP::GeV; }
709  G4ThreeVector beamDirection(partMomentum);
710  double pTotTmp = beamDirection.mag();
711  for (size_t k=0; k !=3; k++) beamDirection[k] /= pTotTmp;
712  fParticleGun->SetParticleMomentumDirection(beamDirection);
713  int particleID = 0;
714  llStrStr >> particleID;
715  G4ParticleTable *pTable=G4ParticleTable::GetParticleTable();
716  G4ParticleDefinition* pDef = pTable->FindParticle(particleID);
717  fParticleGun->SetParticleDefinition(pDef);
718  const double pMass = pDef->GetPDGMass();
719  const double eTot = std::sqrt(pMass*pMass + pTotTmp*pTotTmp);
720  fParticleGun->SetParticleEnergy(eTot - pMass);
721  // Note: In the original version of fluka user routine fluscw.f, the weight was determined as:
722  //
723  // if ((totp.gt.30.) .or. (totp .lt. pmin1)) then ! calculate weight
724  // mywei = 1.
725  // else
726  // mywei = 30./totp
727  // endif
728  // Where totp is the total momentum of the particle.
729  // Not sure if this is still relevant..!
730  //
731  fWeight = 1.0;
732  fParticleGun->GeneratePrimaryVertex(anEvent);
733 
734 
735 }
string
std::string string
Definition: nybbler.cc:12
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
void LBNEPrimaryGeneratorAction::GenerateG4ProtonBeam ( G4Event *  anEvent)

Definition at line 417 of file LBNEPrimaryGeneratorAction.cc.

418 {
419 
420 // If nothing else is set, use a proton beam
421  G4double x0 = 0.;
422  G4double y0 = 0.;
423  G4double z0 = 0.;
424 
425  do
426  {
427  x0 = G4RandGauss::shoot(0.0, fBeamSigmaX);
428  }
429  while(std::abs(x0) > fBeamMaxValX);
430  do
431  {
432  y0 = G4RandGauss::shoot(0.0, fBeamSigmaY);
433  }
434  while(std::abs(y0) > fBeamMaxValY);
435 
436  const double xFromEmitt = x0;
437  const double yFromEmitt = y0;
438  z0 = fBeamOffsetZ;
439 
440  if(!fBeamOnTarget){
441  // b.c. if fBeamOnTarget, all offsets are ignored.
442  x0 += fBeamOffsetX;
443  y0 += fBeamOffsetY;
444  }
445 
446  G4double dx=0.;
447  G4double dy=0.;
448  G4double dz=0.;
449 
450  if(fabs(fBeamAngleTheta) > 1e-6){
451  dx += sin(fBeamAngleTheta)*cos(fBeamAnglePhi);
452  dy += sin(fBeamAngleTheta)*sin(fBeamAnglePhi);
453  }
454  dz = std::sqrt(1.0 - (dx*dx + dy*dy));
455  G4ThreeVector direction(dx, dy, dz);
456 
457  if(fCorrectForAngle || fBeamOnTarget){
458  x0 += (dx/dz)*z0;
459  y0 += (dy/dz)*z0;
460  }
461  double phaseX = 0.;
462  double phaseY = 0.;
463  if (fUseCourantSniderParams) { // The tails in the dx, dy distribution not quite Gaussian,
464  // but this is small error.
465  phaseX = 2.0*M_PI*G4RandFlat::shoot();
466  const double amplX = xFromEmitt/std::cos(phaseX);
467  dx += (std::sin(phaseX)*amplX - fBeamAlphaFunctionX*xFromEmitt)/(fBeamBetaFuncGenX*CLHEP::m); // in radian (mm/mm)
468  phaseY = 2.0*M_PI*G4RandFlat::shoot();
469  const double amplY = yFromEmitt/std::cos(phaseY);
470  dy += (std::sin(phaseY)*amplY - fBeamAlphaFunctionY*yFromEmitt)/(fBeamBetaFuncGenY*CLHEP::m);
471  dz = std::sqrt(1.0 - (dx*dx + dy*dy));
472  direction = G4ThreeVector(dx, dy, dz);
473  }
474  if (fRadiusAnnularBeam > 1.0e-6) {
475  const double r = histAnnularRadial->GetRandom();
476  const double phi = 2.0*M_PI*G4RandFlat::shoot();
477  x0 = r*std::cos(phi);
478  y0 = r*std::sin(phi);
479  }
480 
481 // if (std::abs(y0) > 5.0) {
482 // std::cerr << " !!!!!.... Anomalously large y0 " << y0 << " dy " << dy << " Phase y " << phaseY << std::endl;
483 // }
484 
485  fProtonN = fCurrentPrimaryNo;
486  fProtonOrigin = G4ThreeVector(x0, y0, z0);
487  fProtonMomentum = G4ThreeVector(dx*fProtonMomentumMag, dy*fProtonMomentumMag, dz*fProtonMomentumMag);
488  fProtonIntVertex = G4ThreeVector(-9999.,-9999.,-9999.);
489  fWeight=1.; //for primary protons set weight and tgen
490  fTgen=0;
491 
492  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
493  fParticleGun->SetParticleMomentumDirection(direction);
494  //
495  // Backdoor to create lots of beam particle, records data on where we are on the target..
496  //
497  // G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
498  // fParticleGun->SetParticleDefinition(particleTable->FindParticle("geantino"));
499  fParticleGun->GeneratePrimaryVertex(anEvent);
500 }
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
run_cnn_1class.m
m
Definition: run_cnn_1class.py:75
plot_model.r
r
Definition: plot_model.py:57
void LBNEPrimaryGeneratorAction::GeneratePrimaries ( G4Event *  anEvent)

Definition at line 370 of file LBNEPrimaryGeneratorAction.cc.

371 {
372  int verboseLevel = G4EventManager::GetEventManager()->GetVerboseLevel();
373  if(verboseLevel > 1)
374  {
375  std::cout << "Event " << anEvent->GetEventID() << ": LBNEPrimaryGeneratorAction::GeneratePrimaries() Called." << std::endl;
376  }
377 
378 
379  G4int totNoPrim=fRunManager->GetNumberOfEvents();
380  if (totNoPrim>20)
381  {
382  if (fCurrentPrimaryNo%(totNoPrim/20)==0)
383  G4cout<<"Processing particles #: "
384  <<fCurrentPrimaryNo<<" to "<< fCurrentPrimaryNo+(totNoPrim/20) <<G4endl;
385  }
386 
387  {
388 /*
389  if (fLBNEData->GetUseFlukaInput() || fLBNEData->GetUseMarsInput())
390  {
391  LBNEPrimaryGeneratorAction::GenerateBeamFromInput(anEvent);
392  }
393 */
394  if (fUseGeantino || fUseMuonGeantino || fUseChargedGeantino) {
395  LBNEPrimaryGeneratorAction::Geantino(anEvent);
396  }
397  else if (fFlukaASCII.is_open()) {
398  LBNEPrimaryGeneratorAction::GenerateFlukaPrimary(anEvent);
399  } else {
400  LBNEPrimaryGeneratorAction::GenerateG4ProtonBeam(anEvent);
401  }
402 
403  }
404 
405  //fCurrentPrimaryNo++;
406  ++fCurrentPrimaryNo;
407 }
LBNERunManager::GetNumberOfEvents
G4int GetNumberOfEvents()
Definition: LBNERunManager.hh:26
LBNEPrimaryGeneratorAction::GenerateG4ProtonBeam
void GenerateG4ProtonBeam(G4Event *anEvent)
Definition: LBNEPrimaryGeneratorAction.cc:417
LBNEPrimaryGeneratorAction::GenerateFlukaPrimary
void GenerateFlukaPrimary(G4Event *anEvent)
Definition: LBNEPrimaryGeneratorAction.cc:685
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
G4ThreeVector LBNEPrimaryGeneratorAction::GetAngleForMuGeantinoScan ( ) const
inline

Definition at line 213 of file LBNEPrimaryGeneratorAction.hh.

213 { return fAngleForMuGeantino;}
double LBNEPrimaryGeneratorAction::GetBeamAngleTheta ( ) const
inline

Definition at line 97 of file LBNEPrimaryGeneratorAction.hh.

97 { return fBeamAngleTheta; }
double LBNEPrimaryGeneratorAction::GetBeamBetaFunctionX ( ) const
inline

Definition at line 92 of file LBNEPrimaryGeneratorAction.hh.

92 { return fBeamBetaFunctionX;}
double LBNEPrimaryGeneratorAction::GetBeamBetaFunctionY ( ) const
inline

Definition at line 95 of file LBNEPrimaryGeneratorAction.hh.

95 { return fBeamBetaFunctionY;}
double LBNEPrimaryGeneratorAction::GetBeamEmittanceX ( ) const
inline

Definition at line 91 of file LBNEPrimaryGeneratorAction.hh.

91 { return fBeamEmittanceX;}
double LBNEPrimaryGeneratorAction::GetBeamEmittanceY ( ) const
inline

Definition at line 94 of file LBNEPrimaryGeneratorAction.hh.

94 { return fBeamEmittanceY;}
double LBNEPrimaryGeneratorAction::GetBeamMaxValX ( ) const
inline

Definition at line 88 of file LBNEPrimaryGeneratorAction.hh.

88 { return fBeamMaxValX; }
double LBNEPrimaryGeneratorAction::GetBeamMaxValY ( ) const
inline

Definition at line 89 of file LBNEPrimaryGeneratorAction.hh.

89 { return fBeamMaxValY; }
double LBNEPrimaryGeneratorAction::GetBeamOffsetX ( ) const
inline

Definition at line 85 of file LBNEPrimaryGeneratorAction.hh.

85 { return fBeamOffsetX; }
double LBNEPrimaryGeneratorAction::GetBeamOffsetY ( ) const
inline

Definition at line 86 of file LBNEPrimaryGeneratorAction.hh.

86 { return fBeamOffsetY; }
double LBNEPrimaryGeneratorAction::GetBeamOffsetZ ( ) const
inline

Definition at line 87 of file LBNEPrimaryGeneratorAction.hh.

87 { return fBeamOffsetZ; }
double LBNEPrimaryGeneratorAction::GetBeamSigmaX ( ) const
inline

Definition at line 83 of file LBNEPrimaryGeneratorAction.hh.

83 { return fBeamSigmaX; }
double LBNEPrimaryGeneratorAction::GetBeamSigmaY ( ) const
inline

Definition at line 84 of file LBNEPrimaryGeneratorAction.hh.

84 { return fBeamSigmaY; }
double LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy ( double  pz)

Definition at line 351 of file LBNEPrimaryGeneratorAction.cc.

351  {
352 
353  // quadratic interpolation, based on the e-mail from
354  // John Jostone.
355  G4Exception("LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy", "", FatalException,
356  " With the current 1.2 MW option, the use of Courand Snider params and a beam energy != 120 GeV is prohibited ");
357 
358  // Define and slove a set of quadratic equation.
359  const double betaSecond80 = (fBeamBetaFunctionAt120 + 2.0*fBeamBetaFunctionAt60 -
360  3.0*fBeamBetaFunctionAt80)/2400.;
361  const double betaPrime80 = (1.0/40.)*(fBeamBetaFunctionAt120 - fBeamBetaFunctionAt80 - betaSecond80*1600.);
362 
363  const double dp = pz - 80.;
364  std::cerr << " Using LBNEPrimaryGeneratorAction::GetBetaFunctionvsBeamEnergy... fBeamBetaFunctionAt120 " <<
365  fBeamBetaFunctionAt120 << " At 80 " << fBeamBetaFunctionAt80
366  << " fBeamBetaFunctionAt60 " << fBeamBetaFunctionAt60 << std::endl;
367  return (fBeamBetaFunctionAt80 + betaPrime80*dp + betaSecond80*dp*dp);
368 }
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
G4ThreeVector LBNEPrimaryGeneratorAction::GetMomentumForMuGeantinoScan ( ) const
inline

Definition at line 211 of file LBNEPrimaryGeneratorAction.hh.

211 { return fMomentumForMuGeantino;}
G4int LBNEPrimaryGeneratorAction::GetNoOfPrimaries ( )
inline

Definition at line 39 of file LBNEPrimaryGeneratorAction.hh.

39 { return fNoOfPrimaries; }
G4ParticleGun* LBNEPrimaryGeneratorAction::GetParticleGun ( )
inline

Definition at line 32 of file LBNEPrimaryGeneratorAction.hh.

32 {return fParticleGun;};
G4ThreeVector LBNEPrimaryGeneratorAction::GetParticleMomentum ( )
inline

Definition at line 52 of file LBNEPrimaryGeneratorAction.hh.

52 { return fParticleMomentum; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetParticlePosition ( )
inline

Definition at line 51 of file LBNEPrimaryGeneratorAction.hh.

51 { return fParticlePosition; }
G4int LBNEPrimaryGeneratorAction::GetParticleType ( )
inline

Definition at line 55 of file LBNEPrimaryGeneratorAction.hh.

55 { return fType; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetProtonIntVertex ( )
inline

Definition at line 46 of file LBNEPrimaryGeneratorAction.hh.

46 { return fProtonIntVertex; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetProtonMomentum ( )
inline

Definition at line 45 of file LBNEPrimaryGeneratorAction.hh.

45 { return fProtonMomentum; }
G4int LBNEPrimaryGeneratorAction::GetProtonNumber ( )
inline

Definition at line 57 of file LBNEPrimaryGeneratorAction.hh.

57 { return fProtonN; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetProtonOrigin ( )
inline

Definition at line 44 of file LBNEPrimaryGeneratorAction.hh.

44 { return fProtonOrigin; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetPtForMuGeantinoScan ( ) const
inline

Definition at line 212 of file LBNEPrimaryGeneratorAction.hh.

212 { return fPtForMuGeantino;}
G4int LBNEPrimaryGeneratorAction::GetTgen ( )
inline

Definition at line 54 of file LBNEPrimaryGeneratorAction.hh.

54 { return fTgen; }
bool LBNEPrimaryGeneratorAction::GetUseChargedGeantino ( ) const
inline

Definition at line 198 of file LBNEPrimaryGeneratorAction.hh.

198 { return fUseChargedGeantino; }
G4bool LBNEPrimaryGeneratorAction::GetUseCourantSniderParams ( )
inline

Definition at line 74 of file LBNEPrimaryGeneratorAction.hh.

74 { return fUseCourantSniderParams;}
bool LBNEPrimaryGeneratorAction::GetUseGeantino ( ) const
inline

Definition at line 197 of file LBNEPrimaryGeneratorAction.hh.

197 { return fUseGeantino; }
G4bool LBNEPrimaryGeneratorAction::GetUseJustSigmaCoord ( )
inline

Definition at line 75 of file LBNEPrimaryGeneratorAction.hh.

75 {return fUseJustSigmaCoord;}
bool LBNEPrimaryGeneratorAction::GetUseMuonGeantino ( ) const
inline

Definition at line 196 of file LBNEPrimaryGeneratorAction.hh.

196 { return fUseMuonGeantino; }
G4double LBNEPrimaryGeneratorAction::GetWeight ( )
inline

Definition at line 53 of file LBNEPrimaryGeneratorAction.hh.

53 { return fWeight; }
G4ThreeVector LBNEPrimaryGeneratorAction::GetYOriginForMuGeantinoScan ( ) const
inline

Definition at line 210 of file LBNEPrimaryGeneratorAction.hh.

210 { return fYOriginForMuGeantino;}
G4ThreeVector LBNEPrimaryGeneratorAction::GetZOriginForMuGeantinoScan ( ) const
inline

Definition at line 209 of file LBNEPrimaryGeneratorAction.hh.

209 { return fZOriginForMuGeantino;}
G4bool LBNEPrimaryGeneratorAction::OpenNtuple ( G4String  ntupleName)

Definition at line 279 of file LBNEPrimaryGeneratorAction.cc.

280 {
281  G4String message("Input Ntuple");
282  message += ntupleName + G4String("Not yet supported");
283  G4Exception("LBNEPrimaryGeneratorAction::OpenNtuple", " ",
284  FatalErrorInArgument, message.c_str());
285  /*
286  fCurrentPrimaryNo=0;
287 
288  G4bool fIsOpen=false;
289  fInputFile=new TFile(ntupleName,"READ");
290  if (!fInputFile->IsZombie())
291  {
292 
293  //fInputTree=(TTree*)(fInputFile->Get("h3"));
294  fInputTree=(TTree*)(fInputFile->Get((fLBNEData->GetInputNtpTreeName()).c_str()));
295 
296  if(!fInputTree)
297  {
298  //G4cout<<"LBNEPrimaryGeneratorAction: Can't find tree "<< G4endl;
299  G4cout<<"LBNEPrimaryGeneratorAction: Can't find tree with name "
300  << fLBNEData->GetInputNtpTreeName() << G4endl;
301  }
302 
303  fCurrentPrimaryNo=0;
304  G4int entries = fInputTree->GetEntries();
305  if(fLBNEData->GetNEvents() > 0 && fLBNEData->GetNEvents() < entries)
306  fNoOfPrimaries = fLBNEData->GetNEvents();
307  else
308  fNoOfPrimaries = fInputTree->GetEntries();
309 
310 
311  fIsOpen=true;
312  }
313  else
314  {
315  G4cout<<"LBNEPrimaryGeneratorAction: Input (fluka/mars) root file doesn't exist"
316  <<" Aborting run"<<G4endl;
317  fIsOpen=false;
318  }
319  return fIsOpen;
320  */
321  return false;
322 }
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
python.root_metadata.message
message
Definition: root_metadata.py:25
G4bool LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII ( G4String  ntupleName)

Definition at line 252 of file LBNEPrimaryGeneratorAction.cc.

253 {
254  if (fFlukaASCII.is_open()) fFlukaASCII.close();
255  fFlukaASCII.open(ntupleName.c_str());
256  if (!fFlukaASCII.is_open()) {
257  G4Exception("LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII", " ",
258  FatalErrorInArgument, " Fluka Input file is not open, unexpected and fatal " );
259  }
260  maxLengthForFlukaLine = 1000;
261  lineForFluka = static_cast<char*>(malloc(maxLengthForFlukaLine));
262  fFlukaASCII.getline(lineForFluka, maxLengthForFlukaLine);
263  if (!fFlukaASCII.good()) {
264  G4Exception("LBNEPrimaryGeneratorAction::OpenNtupleFLUKAASCII", " ",
265  FatalErrorInArgument, " Fluka Input file no longer accessible " );
266  }
267  if (fFlukaASCII.eof()) {
268  fFlukaASCII.close();
269  G4Exception("LBNEPrimaryGeneratorAction::GenerateFlukaPrimary", " ",
270  FatalErrorInArgument, " Fluka Input file is empty.... " );
271  }
272  //
273  // The first line if the Header for R..
274  //
275  return true;
276 
277 }
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
void LBNEPrimaryGeneratorAction::SetAngleForMuGeantinoScan ( G4ThreeVector  v)
inline

Definition at line 207 of file LBNEPrimaryGeneratorAction.hh.

207 {fAngleForMuGeantino = v;}
void LBNEPrimaryGeneratorAction::SetBeamBetaFunctionX ( G4double  t)
inline

Definition at line 79 of file LBNEPrimaryGeneratorAction.hh.

79 {fBeamBetaFunctionX = t;}
void LBNEPrimaryGeneratorAction::SetBeamBetaFunctionY ( G4double  t)
inline

Definition at line 80 of file LBNEPrimaryGeneratorAction.hh.

80 {fBeamBetaFunctionY = t;}
void LBNEPrimaryGeneratorAction::SetBeamEmittanceX ( G4double  t)
inline

Definition at line 77 of file LBNEPrimaryGeneratorAction.hh.

77 {fBeamEmittanceX = t;}
void LBNEPrimaryGeneratorAction::SetBeamEmittanceY ( G4double  t)
inline

Definition at line 78 of file LBNEPrimaryGeneratorAction.hh.

78 {fBeamEmittanceY = t;}
void LBNEPrimaryGeneratorAction::SetBeamMaxValX ( G4double  x)
inline

Definition at line 65 of file LBNEPrimaryGeneratorAction.hh.

65 { fBeamMaxValX = x; }
train.x
list x
Definition: train.py:276
void LBNEPrimaryGeneratorAction::SetBeamMaxValY ( G4double  y)
inline

Definition at line 66 of file LBNEPrimaryGeneratorAction.hh.

66 { fBeamMaxValY = y; }
y
double y
Definition: GapWidth_module.cc:109
void LBNEPrimaryGeneratorAction::SetBeamOffsetX ( G4double  x)
inline

Definition at line 61 of file LBNEPrimaryGeneratorAction.hh.

61 { fBeamOffsetX = x; }
train.x
list x
Definition: train.py:276
void LBNEPrimaryGeneratorAction::SetBeamOffsetY ( G4double  y)
inline

Definition at line 62 of file LBNEPrimaryGeneratorAction.hh.

62 { fBeamOffsetY = y; }
y
double y
Definition: GapWidth_module.cc:109
void LBNEPrimaryGeneratorAction::SetBeamOnTarget ( G4bool  aBool)
inline

Definition at line 60 of file LBNEPrimaryGeneratorAction.hh.

60 { fBeamOnTarget = aBool; }
void LBNEPrimaryGeneratorAction::SetBeamPhi ( G4double  phi)
inline

Definition at line 68 of file LBNEPrimaryGeneratorAction.hh.

68 { fBeamAnglePhi = phi; }
void LBNEPrimaryGeneratorAction::SetBeamSigmaX ( G4double  x)
inline

Definition at line 63 of file LBNEPrimaryGeneratorAction.hh.

63 { fBeamSigmaX = x; }
train.x
list x
Definition: train.py:276
void LBNEPrimaryGeneratorAction::SetBeamSigmaY ( G4double  y)
inline

Definition at line 64 of file LBNEPrimaryGeneratorAction.hh.

64 { fBeamSigmaY = y; }
y
double y
Definition: GapWidth_module.cc:109
void LBNEPrimaryGeneratorAction::SetBeamTheta ( G4double  theta)
inline

Definition at line 67 of file LBNEPrimaryGeneratorAction.hh.

67 { fBeamAngleTheta = theta; }
void LBNEPrimaryGeneratorAction::SetCorrectForAngle ( G4bool  aBool)
inline

Definition at line 59 of file LBNEPrimaryGeneratorAction.hh.

59 { fCorrectForAngle = aBool; }
void LBNEPrimaryGeneratorAction::SetMomentumForMuGeantinoScan ( G4ThreeVector  v)
inline

Definition at line 205 of file LBNEPrimaryGeneratorAction.hh.

205 {fMomentumForMuGeantino = v;}
void LBNEPrimaryGeneratorAction::SetPolarAngleGeantino ( double  r)
inline

Definition at line 188 of file LBNEPrimaryGeneratorAction.hh.

188 {fPolarAngleGeantino=r;}
plot_model.r
r
Definition: plot_model.py:57
void LBNEPrimaryGeneratorAction::SetPolarAngleGeantinoMin ( double  r)
inline

Definition at line 189 of file LBNEPrimaryGeneratorAction.hh.

189 {fPolarAngleGeantinoMin=r;}
plot_model.r
r
Definition: plot_model.py:57
void LBNEPrimaryGeneratorAction::SetProtonBeam ( )

Definition at line 132 of file LBNEPrimaryGeneratorAction.cc.

133 {
134 
135  // Make sure we have a consistent way to define the beam
136 
137  if (fUseCourantSniderParams && fUseJustSigmaCoord) {
138  G4Exception("LBNEPrimaryGeneratorAction::SetProtonBeam", " ",
139  FatalErrorInArgument,
140  "Attempting to refdefined both CourantSnider parameters and beam spot size. Please choose one!..." );
141  }
142  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
143 
144  if (fUseGeantino) {
145  fParticleGun->SetParticleDefinition(particleTable->FindParticle("geantino"));
146  }
147  else if (fUseMuonGeantino) {
148  fParticleGun->SetParticleDefinition(particleTable->FindParticle("mu+"));
149  }
150  else if (fUseChargedGeantino) {
151  fParticleGun->SetParticleDefinition(particleTable->FindParticle("chargedgeantino"));
152  }
153  else fParticleGun->SetParticleDefinition(particleTable->FindParticle("proton"));
154 
155  const double eBeam = std::sqrt(fProtonMomentumMag*fProtonMomentumMag + 0.938272*CLHEP::GeV*0.938272*CLHEP::GeV);
156  const double pEkin = eBeam - 0.938272*CLHEP::GeV;
157  fParticleGun->SetParticleEnergy(pEkin);
158  // This needs to be implemented via data cards!!!
159  G4ThreeVector beamPosition(0., 0., fBeamOffsetZ);
160 
161  fParticleGun->SetParticlePosition(beamPosition);
162  G4ThreeVector beamDirection(1,0,0);
163  beamDirection.setTheta(fBeamAngleTheta);
164  beamDirection.setPhi(fBeamAnglePhi);
165  G4cout << "Beam direction of " << beamDirection[0] << " " <<
166  beamDirection[1] << " " << beamDirection[2] << endl;
167  fParticleGun->SetParticleMomentumDirection(beamDirection);
168 
169  fCurrentPrimaryNo=0;
170 
171 // LBNEVolumePlacements* volP=LBNEVolumePlacements::Instance();
172 // bool use1p2MW = volP->GetUse1p2MW();
173  // Upgrade the beam sigma April 17: also the beam emittance!
174 // if (use1p2MW) { fBeamSigmaX = 1.7*CLHEP::mm; fBeamSigmaY = 1.7*CLHEP::mm; fBeamBetaFunctionX = 110.9; fBeamBetaFunctionY = 110.9;}
175 
176 
177  G4String spaces = " ";
178  std::cout << spaces << "Configuring the Proton beam..." << std::endl
179  << spaces << " Momentum = " << fParticleGun->GetParticleMomentum()/CLHEP::GeV << " GeV/c" << std::endl
180  << spaces << " Kinetic Energy = " << fParticleGun->GetParticleEnergy()/CLHEP::GeV << " GeV" << std::endl
181  << spaces << " Beam Offset, X = " << fBeamOffsetX/CLHEP::mm << " mm" << std::endl
182  << spaces << " Beam Offset, Y = " << fBeamOffsetY/CLHEP::mm << " mm" << std::endl
183  << spaces << " Beam Offset, Z = " << fBeamOffsetZ/CLHEP::m << " m" << std::endl
184  << spaces << " Beam Sigma, X = " << fBeamSigmaX/CLHEP::mm << " mm" << std::endl
185  << spaces << " Beam Sigma, Y = " << fBeamSigmaY/CLHEP::mm << " mm" << std::endl
186  << spaces << " Beam Max. d|X| = " << fBeamMaxValX/CLHEP::mm << " mm" << std::endl
187  << spaces << " Beam Max, d|Y| = " << fBeamMaxValY/CLHEP::mm << " mm" << std::endl
188  << spaces << " Direction = " << fParticleGun->GetParticleMomentumDirection() << std::endl;
189 
190  if (fUseCourantSniderParams && (!fUseMuonGeantino) && (!fUseChargedGeantino)) {
191 
192  const double gammaLorentz = eBeam/(0.938272*CLHEP::GeV);
193  const double betaLorentz = std::sqrt(1.0 - 1.0/(gammaLorentz*gammaLorentz));
194  const double betaGammaLorentz = gammaLorentz*betaLorentz;
195 
196  const double betaFuncStarX = (std::abs(fProtonMomentumMag - 120.*CLHEP::GeV) < 0.1) ?
197  fBeamBetaFunctionX :
198  GetBetaFunctionvsBeamEnergy(fProtonMomentumMag/CLHEP::GeV); // in meters, not mm
199  const double betaFuncStarY = (std::abs(fProtonMomentumMag - 120.*CLHEP::GeV) < 0.1) ?
200  fBeamBetaFunctionY :
201  GetBetaFunctionvsBeamEnergy(fProtonMomentumMag/CLHEP::GeV);
202 
203  fBeamAlphaFunctionX = -(fBeamOffsetZ/CLHEP::m)/betaFuncStarX;
204  fBeamAlphaFunctionY = -(fBeamOffsetZ/CLHEP::m)/betaFuncStarY;
205 
206  fBeamBetaFuncGenX = betaFuncStarX + ((fBeamOffsetZ/CLHEP::m)*(fBeamOffsetZ/CLHEP::m))/betaFuncStarX;// in meters
207  fBeamBetaFuncGenY = betaFuncStarY + ((fBeamOffsetZ/CLHEP::m)*(fBeamOffsetZ/CLHEP::m))/betaFuncStarY;
208 
209  fBeamSigmaX = std::sqrt( 1.0e3 * fBeamBetaFuncGenX/CLHEP::m * fBeamEmittanceX/ (6.0*betaGammaLorentz)); // in mm
210  fBeamSigmaY = std::sqrt( 1.0e3 * fBeamBetaFuncGenY/CLHEP::m * fBeamEmittanceY/ (6.0*betaGammaLorentz));
211  std::cout << spaces << "Re-Configuring the Proton beam using Courant-Snider parameters " << std::endl;
212  std::cout << spaces << " Beta function X at Point of creation " << fBeamBetaFuncGenX << " m " << std::endl;
213  std::cout << spaces << " Beta function Y at Point of creation " << fBeamBetaFuncGenY << " m " << std::endl;
214  std::cout << spaces << " Alpha function X at Point of creation " << fBeamAlphaFunctionX << std::endl;
215  std::cout << spaces << " Alpha function Y at Point of creation " << fBeamAlphaFunctionY << std::endl;
216  std::cout << spaces << " Sigma X at Point of creation " << fBeamSigmaX << " mm " << std::endl;
217  std::cout << spaces << " Sigma Y at Point of creation " << fBeamSigmaY << " mm " << std::endl;
218 
219  }
220  if (fRadiusAnnularBeam > 1.0e-6) {
221  // We will be using some of the ROOT math utilities, namely generating
222  // random number from a histogram. But we need to control the seed....
223  // The following statement may cause havoc in other part of the code where
224  // The root Random generator is called..
225  //
226  int startSeed = CLHEP::HepRandom::getTheSeed();
227  startSeed += 1267; // a bit arbitrary.
228  gRandom = new TRandom3(0);
229  gRandom->SetSeed(startSeed);
230  fSigmaForAnnularbeam = std::sqrt(fBeamSigmaX*fBeamSigmaX + fBeamSigmaY*fBeamSigmaY);
231  std::cout << spaces << " But.. but we will be using an annular beam, defined as follow " << std::endl;
232  std::cout << spaces << " f(r; R, sigma) = d^2N/drdTheta = " << std::endl;
233  std::cout << spaces << " = 1/(2*pi*sigma^2) * exp[(-r^2-R^2)/(2*sigma^2)] * I0(rR/sigma^2) " << std::endl;
234  std::cout << spaces << " with R = " << fRadiusAnnularBeam << " and sigma = "
235  << fSigmaForAnnularbeam << std::endl;
236  std::cout << spaces << " This sigma is derived from Sigma X and Sigma Y, stated above " << std::endl;
237  // Prepare the distirbution.
238  theModBesselFunc1rstKind = new TF1("I_0", "TMath::BesselI0(x)", 0., 10.);
239  histAnnularRadial = new TH1F("AN1", "Annular beam, radial", 1000, 0., 6.*fRadiusAnnularBeam);
240  const double fSigSq = fSigmaForAnnularbeam*fSigmaForAnnularbeam;
241  const double RSq = fRadiusAnnularBeam*fRadiusAnnularBeam;
242  const double stepR = 6.*fRadiusAnnularBeam/1000.;
243  for (int kR=0.; kR != 1000; kR++) {
244  const double r = 0.49*stepR + kR*stepR;
245  const double f = (1.0/(2.0*M_PI*fSigSq)) * std::exp((-r*r -RSq)/(2.*fSigSq)) *
246  theModBesselFunc1rstKind->Eval(r*fRadiusAnnularBeam/fSigSq);
247  histAnnularRadial->Fill(r, r*f);
248  }
249  }
250 
251 }
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
run_cnn_1class.m
m
Definition: run_cnn_1class.py:75
plot_model.r
r
Definition: plot_model.py:57
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void LBNEPrimaryGeneratorAction::SetProtonMomentum ( G4double  p)
inline

Definition at line 69 of file LBNEPrimaryGeneratorAction.hh.

69 {fProtonMomentumMag = p;}
test.p
p
Definition: test.py:223
void LBNEPrimaryGeneratorAction::SetPtForMuGeantinoScan ( G4ThreeVector  v)
inline

Definition at line 206 of file LBNEPrimaryGeneratorAction.hh.

206 {fPtForMuGeantino = v;}
void LBNEPrimaryGeneratorAction::SetRadiusAnnularBeam ( G4double  t)
inline

Definition at line 81 of file LBNEPrimaryGeneratorAction.hh.

81 {fRadiusAnnularBeam = t;}
void LBNEPrimaryGeneratorAction::SetSigmaZOriginGeantino ( double  v)
inline

Definition at line 194 of file LBNEPrimaryGeneratorAction.hh.

194 {fSigmaZOriginGeantino=v;}
void LBNEPrimaryGeneratorAction::SetUseChargedGeantino ( bool  t)
inline

Definition at line 192 of file LBNEPrimaryGeneratorAction.hh.

192 {fUseChargedGeantino=t;}
void LBNEPrimaryGeneratorAction::SetUseCourantSniderParams ( G4bool  t)
inline

Definition at line 71 of file LBNEPrimaryGeneratorAction.hh.

71 {fUseCourantSniderParams = t;}
void LBNEPrimaryGeneratorAction::SetUseGeantino ( bool  t)
inline

Definition at line 190 of file LBNEPrimaryGeneratorAction.hh.

190 {fUseGeantino=t;}
void LBNEPrimaryGeneratorAction::SetUseJustSigmaCoord ( G4bool  t)
inline

Definition at line 72 of file LBNEPrimaryGeneratorAction.hh.

72 {fUseJustSigmaCoord = t;}
void LBNEPrimaryGeneratorAction::SetUseMuonGeantino ( bool  t)
inline

Definition at line 191 of file LBNEPrimaryGeneratorAction.hh.

191 {fUseMuonGeantino=t;}
void LBNEPrimaryGeneratorAction::SetYOriginForMuGeantinoScan ( G4ThreeVector  v)
inline

Definition at line 204 of file LBNEPrimaryGeneratorAction.hh.

204 {fYOriginForMuGeantino = v;}
void LBNEPrimaryGeneratorAction::SetZOriginForMuGeantinoScan ( G4ThreeVector  v)
inline

Definition at line 203 of file LBNEPrimaryGeneratorAction.hh.

203 {fZOriginForMuGeantino = v;}
void LBNEPrimaryGeneratorAction::SetZOriginGeantino ( double  v)
inline

Definition at line 193 of file LBNEPrimaryGeneratorAction.hh.

193 {fZOriginGeantino=v;}

Member Data Documentation

G4ThreeVector LBNEPrimaryGeneratorAction::fAngleForMuGeantino
private

Definition at line 176 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamAlphaFunctionX
private

Definition at line 157 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamAlphaFunctionY
private

Definition at line 158 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamAnglePhi
private

Definition at line 140 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamAngleTheta
private

Definition at line 139 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFuncGenX
private

Definition at line 159 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFuncGenY
private

Definition at line 160 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFunctionAt120
private

Definition at line 150 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFunctionAt60
private

Definition at line 152 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFunctionAt80
private

Definition at line 151 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFunctionX
private

Definition at line 148 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamBetaFunctionY
private

Definition at line 149 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamEmittanceX
private

Definition at line 146 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamEmittanceY
private

Definition at line 147 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamMaxValX
private

Definition at line 135 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamMaxValY
private

Definition at line 136 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamOffsetX
private

Definition at line 132 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamOffsetY
private

Definition at line 133 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamOffsetZ
private

Definition at line 134 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fBeamOnTarget
private

Definition at line 131 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamSigmaX
private

Definition at line 137 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fBeamSigmaY
private

Definition at line 138 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fCorrectForAngle
private

Definition at line 130 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fCurrentPrimaryNo
private

Definition at line 114 of file LBNEPrimaryGeneratorAction.hh.

std::ifstream LBNEPrimaryGeneratorAction::fFlukaASCII
private

Definition at line 108 of file LBNEPrimaryGeneratorAction.hh.

TFile* LBNEPrimaryGeneratorAction::fInputFile
private

Definition at line 106 of file LBNEPrimaryGeneratorAction.hh.

TTree* LBNEPrimaryGeneratorAction::fInputTree
private

Definition at line 107 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fMomentumForMuGeantino
private

Definition at line 174 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fMuGeantinoCnt
private

Definition at line 127 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fNoOfPrimaries
private

Definition at line 113 of file LBNEPrimaryGeneratorAction.hh.

G4ParticleGun* LBNEPrimaryGeneratorAction::fParticleGun
private

Definition at line 104 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fParticleMomentum
private

Definition at line 122 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fParticlePosition
private

Definition at line 123 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fPolarAngleGeantino
private

Definition at line 167 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fPolarAngleGeantinoMin
private

Definition at line 168 of file LBNEPrimaryGeneratorAction.hh.

LBNEPrimaryMessenger* LBNEPrimaryGeneratorAction::fPrimaryMessenger
private

Definition at line 102 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fProtonIntVertex
private

Definition at line 121 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fProtonMomentum
private

Definition at line 120 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fProtonMomentumMag
private

Definition at line 118 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fProtonN
private

Definition at line 112 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fProtonOrigin
private

Definition at line 119 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fPtForMuGeantino
private

Definition at line 175 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fRadiusAnnularBeam
private

Definition at line 144 of file LBNEPrimaryGeneratorAction.hh.

LBNERunManager* LBNEPrimaryGeneratorAction::fRunManager
private

Definition at line 101 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fSigmaForAnnularbeam
private

Definition at line 145 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fSigmaZOriginGeantino
private

Definition at line 166 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fTgen
private

Definition at line 125 of file LBNEPrimaryGeneratorAction.hh.

G4int LBNEPrimaryGeneratorAction::fType
private

Definition at line 126 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fUseChargedGeantino
private

Definition at line 164 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fUseCourantSniderParams
private

Definition at line 142 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fUseGeantino
private

Definition at line 162 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fUseJustSigmaCoord
private

Definition at line 143 of file LBNEPrimaryGeneratorAction.hh.

G4bool LBNEPrimaryGeneratorAction::fUseMuonGeantino
private

Definition at line 163 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fWeight
private

Definition at line 128 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fYOriginForMuGeantino
private

Definition at line 173 of file LBNEPrimaryGeneratorAction.hh.

G4ThreeVector LBNEPrimaryGeneratorAction::fZOriginForMuGeantino
private

Definition at line 172 of file LBNEPrimaryGeneratorAction.hh.

G4double LBNEPrimaryGeneratorAction::fZOriginGeantino
private

Definition at line 165 of file LBNEPrimaryGeneratorAction.hh.

TH1F* LBNEPrimaryGeneratorAction::histAnnularRadial
private

Definition at line 181 of file LBNEPrimaryGeneratorAction.hh.

char* LBNEPrimaryGeneratorAction::lineForFluka
private

Definition at line 110 of file LBNEPrimaryGeneratorAction.hh.

size_t LBNEPrimaryGeneratorAction::maxLengthForFlukaLine
private

Definition at line 109 of file LBNEPrimaryGeneratorAction.hh.

TF1* LBNEPrimaryGeneratorAction::theModBesselFunc1rstKind
private

Definition at line 180 of file LBNEPrimaryGeneratorAction.hh.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11