evgen::MUSUN Class Reference

module to produce single or multiple specified particles in the detector More...

Inheritance diagram for evgen::MUSUN:

Public Member Functions
	MUSUN (fhicl::ParameterSet const &pset)
void	produce (art::Event &evt)
void	beginJob ()
void	beginRun (art::Run &run)
void	endRun (art::Run &run)
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	SampleOne (unsigned int i, simb::MCTruth &mct, CLHEP::HepRandomEngine &engine)
void	initialization (double theta1, double theta2, double phi1, double phi2, int figflag, double s_hor, double s_ver1, double s_ver2, double &FI)
void	sampling (double &E, double &theta, double &phi, double &dep, CLHEP::HepRandomEngine &engine)

Private Attributes
CLHEP::HepRandomEngine &	fEngine
	art-managed random-number engine More...
int	fPDG
	PDG code of particles to generate. More...
double	fChargeRatio
	Charge ratio of particle / anti-particle. More...
std::string	fInputDir
	Input Directory. More...
std::string	fInputFile1
	Input File 1. More...
std::string	fInputFile2
	Input File 2. More...
std::string	fInputFile3
	Input File 3. More...
double	fCavernAngle
	Angle of the detector from the North to the East. More...
double	fRockDensity
double	fEmin
	Minimum Kinetic Energy (GeV) More...
double	fEmax
	Maximum Kinetic Energy (GeV) More...
double	fThetamin
	Minimum theta. More...
double	fThetamax
	Maximum theta. More...
double	fPhimin
	Minimum phi. More...
double	fPhimax
	Maximum phi. More...
int	figflag
	If want sampled from sphere or parallelepiped. More...
double	fXmin
	Minimum X position. More...
double	fYmin
	Minimum Y position. More...
double	fZmin
	Minimum Z position. More...
double	fXmax
	Maximum X position. More...
double	fYmax
	Maximum Y position. More...
double	fZmax
	Maximum Z position. More...
double	fT0
	Central t position (ns) in world coordinates. More...
double	fSigmaT
	Variation in t position (ns) More...
int	fTDist
	How to distribute t (gaus, or uniform) More...
int	PdgCode
double	Energy
double	phi
double	theta
double	dep
double	Time
double	Momentum
double	px0
double	py0
double	pz0
double	x0
double	y0
double	z0
double	cx
double	cy
double	cz
double	FI = 0.
double	s_hor = 0.
double	s_ver1 = 0.
double	s_ver2 = 0.
double	spmu [121][62][51]
double	fnmu [32401]
double	depth [360][91]
double	fmu [360][91]
double	the1
double	the2
double	ph1
double	ph2
double	se = 0.
double	st = 0.
double	sp = 0.
double	sd = 0.
unsigned int	NEvents = 0
TTree *	fTree

Static Private Attributes
static const int	kGAUS = 1

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

module to produce single or multiple specified particles in the detector

Definition at line 190 of file MUSUN_module.cc.

Constructor & Destructor Documentation

evgen::MUSUN::MUSUN ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 322 of file MUSUN_module.cc.

    : art::EDProducer{pset}
    // create a default random engine; obtain the random seed from NuRandomService,
    // unless overridden in configuration with key "Seed"
    , fEngine(art::ServiceHandle<rndm::NuRandomService> {}->createEngine(*this, pset, "Seed"))
    , fPDG{pset.get<int>("PDG")}
    , fChargeRatio{pset.get<double>("ChargeRatio")}
    , fInputDir{pset.get<std::string>("InputDir")}
    , fInputFile1{pset.get<std::string>("InputFile1")}
    , fInputFile2{pset.get<std::string>("InputFile2")}
    , fInputFile3{pset.get<std::string>("InputFile3")}
    , fCavernAngle{pset.get<double>("CavernAngle")}
    , fRockDensity{pset.get<double>("RockDensity")}
    , fEmin{pset.get<double>("Emin")}
    , fEmax{pset.get<double>("Emax")}
    , fThetamin{pset.get<double>("Thetamin")}
    , fThetamax{pset.get<double>("Thetamax")}
    , fPhimin{pset.get<double>("Phimin")}
    , fPhimax{pset.get<double>("Phimax")}
    , figflag{pset.get<int>("igflag")}
    , fXmin{pset.get<double>("Xmin")}
    , fYmin{pset.get<double>("Ymin")}
    , fZmin{pset.get<double>("Zmin")}
    , fXmax{pset.get<double>("Xmax")}
    , fYmax{pset.get<double>("Ymax")}
    , fZmax{pset.get<double>("Zmax")}
    , fT0{pset.get<double>("T0")}
    , fSigmaT{pset.get<double>("SigmaT")}
    , fTDist{pset.get<int>("TDist")}
  {
    produces<std::vector<simb::MCTruth>>();
    produces<sumdata::RunData, art::InRun>();
  }

Member Function Documentation

void evgen::MUSUN::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 360 of file MUSUN_module.cc.

  {
    // Make the Histograms....
    art::ServiceHandle<art::TFileService const> tfs;
    /*
    hPDGCode      = tfs->make<TH1D>("hPDGCode"     ,"PDG Code"      ,30 , -15           , 15             );
    hPositionX    = tfs->make<TH1D>("hPositionX"   ,"Position (cm)" ,500, ( fXmin - 10 ), ( fXmax + 10 ) );
    hPositionY    = tfs->make<TH1D>("hPositionY"   ,"Position (cm)" ,500, ( fYmin - 10 ), ( fYmax + 10 ) );
    hPositionZ    = tfs->make<TH1D>("hPositionZ"   ,"Position (cm)" ,500, ( fZmin - 10 ), ( fZmax + 10 ) );
    hTime         = tfs->make<TH1D>("hTime"        ,"Time (s)"      ,500, 0    , 1e6       );
    hMomentumHigh = tfs->make<TH1D>("hMomentumHigh","Momentum (GeV)",500, fEmin, fEmax     );
    hMomentum     = tfs->make<TH1D>("hMomentum"    ,"Momentum (GeV)",500, fEmin, fEmin+1e3 );
    hEnergyHigh   = tfs->make<TH1D>("hEnergyHigh"  ,"Energy (GeV)"  ,500, fEmin, fEmax     );
    hEnergy       = tfs->make<TH1D>("hEnergy"      ,"Energy (GeV)"  ,500, fEmin, fEmin+1e3 );
    hDepth        = tfs->make<TH1D>("hDepth"       ,"Depth (m)"     ,800, 0    , 14000     );

    hDirCosineX = tfs->make<TH1D>("hDirCosineX","Normalised Direction cosine",500, -1, 1 );
    hDirCosineY = tfs->make<TH1D>("hDirCosineY","Normalised Direction cosine",500, -1, 1 );
    hDirCosineZ = tfs->make<TH1D>("hDirCosineZ","Normalised Direction cosine",500, -1, 1 );

    hTheta      = tfs->make<TH1D>("hTheta"     ,"Angle (degrees)",500, 0, 90  );
    hPhi        = tfs->make<TH1D>("hPhi"       ,"Angle (degrees)",500, 0, 365 );
    */
    fTree = tfs->make<TTree>("Generator", "analysis tree");
    fTree->Branch("particleID", &PdgCode, "particleID/I");
    fTree->Branch("energy", &Energy, "energy/D");
    fTree->Branch("time", &Time, "Time/D");
    fTree->Branch("posX", &x0, "posX/D");
    fTree->Branch("posY", &y0, "posY/D");
    fTree->Branch("posZ", &z0, "posZ/D");
    fTree->Branch("cosX", &cx, "cosX/D");
    fTree->Branch("cosY", &cy, "cosY/D");
    fTree->Branch("cosZ", &cz, "cosZ/D");
    fTree->Branch("theta", &theta, "theta/D");
    fTree->Branch("phi", &phi, "phi/D");
    fTree->Branch("depth", &dep, "dep/D");
    /*
    fCryos = tfs->make<TTree>("CryoSizes","cryo tree");
    fCryos->Branch("NumTPCs" , &NumTPCs , "NumTPCs/I"       );
    fCryos->Branch("TPCMinX" , &TPCMinX , "TPCMinX[NumTPCs]/D");
    fCryos->Branch("TPCMaxX" , &TPCMaxX , "TPCMaxX[NumTPCs]/D");
    fCryos->Branch("TPCMinY" , &TPCMinY , "TPCMinY[NumTPCs]/D");
    fCryos->Branch("TPCMaxY" , &TPCMaxY , "TPCMaxY[NumTPCs]/D");
    fCryos->Branch("TPCMinZ" , &TPCMinZ , "TPCMinZ[NumTPCs]/D");
    fCryos->Branch("TPCMaxZ" , &TPCMaxZ , "TPCMaxZ[NumTPCs]/D");
    fCryos->Branch("CryoSize", &CryoSize, "CryoSize[6]/D"   );
    fCryos->Branch("DetHall" , &DetHall , "DetHall[6]/D"    );
    */
  }

void evgen::MUSUN::beginRun ( art::Run &  run )

virtual

Reimplemented from art::EDProducer.

Definition at line 414 of file MUSUN_module.cc.

  {
    // Check fcl parameters were set correctly
    if (fThetamax > 90.5)
      throw cet::exception("MUSUNGen")
        << "\nThetamax has to be less than " << M_PI / 2 << ", but was entered as " << fThetamax
        << ", this causes an error so leaving program now...\n\n";
    if (fThetamin < 0)
      throw cet::exception("MUSUNGen")
        << "\nThetamin has to be more than 0, but was entered as " << fThetamin
        << ", this causes an error so leaving program now...\n\n";
    if (fThetamax < fThetamin)
      throw cet::exception("MUSUNGen") << "\nMinimum angle is bigger than maximum angle....causes "
                                          "an error so leaving program now....\n\n";
    if (fPhimax > 360.5)
      throw cet::exception("MUSUNGen")
        << "\nPhimax has to be less than " << 2 * M_PI << ", but was entered as " << fPhimax
        << ", this cause an error so leaving program now...\n\n";
    if (fPhimin < 0)
      throw cet::exception("MUSUNGen")
        << "\nPhimin has to be more than 0, but was entered as " << fPhimin
        << ", this causes an error so leaving program now...\n\n";
    if (fPhimax < fPhimin)
      throw cet::exception("MUSUNGen") << "\nMinimum angle is bigger than maximum angle....causes "
                                          "an error so leaving program now....\n\n";
    if (fEmax < fEmin)
      throw cet::exception("MUSUNGen")
        << "\nMinimum energy is bigger than maximum energy....causes an error so leaving program "
           "now....\n\n";

    // grab the geometry object to see what geometry we are using
    art::ServiceHandle<geo::Geometry const> geo;
    run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()));

    // area of the horizontal plane of the parallelepiped
    s_hor = (fZmax - fZmin) * (fXmax - fXmin);
    // area of the vertical plane of the parallelepiped, perpendicular to z-axis
    s_ver1 = (fXmax - fXmin) * (fYmax - fYmin);
    // area of the vertical plane of the parallelepiped, perpendicular to x-axis
    s_ver2 = (fZmax - fZmin) * (fYmax - fYmin);

    //std::cout << s_hor << " " << s_ver1 << " " << s_ver2 << std::endl;

    initialization(fThetamin, fThetamax, fPhimin, fPhimax, figflag, s_hor, s_ver1, s_ver2, FI);

    std::cout << "Material - SURF rock" << std::endl;
    std::cout << "Density = " << fRockDensity << " g/cm^3" << std::endl;
    std::cout << "Parameters for muon spectrum are from LVD best fit" << std::endl;
    std::cout << "Muon energy range = " << fEmin << " - " << fEmax << " GeV" << std::endl;
    std::cout << "Zenith angle range = " << fThetamin << " - " << fThetamax << " degrees"
              << std::endl;
    std::cout << "Azimuthal angle range = " << fPhimin << " - " << fPhimax << " degrees"
              << std::endl;
    std::cout << "Global intensity = " << FI
              << " (cm^2 s)^(-1) or s^(-1) (for muons on the surface)" << std::endl;
    /*
    NumTPCs = geo->NTPC(0);
    std::cout << "There are " << NumTPCs << " in cryostat 0" << std::endl;
    for (unsigned int c=0; c<geo->Ncryostats(); c++) {
      const geo::CryostatGeo& cryostat=geo->Cryostat(c);
      geo->CryostatBoundaries( CryoSize, 0 );
      std::cout << "Cryo bounds " << CryoSize[0] << " "<< CryoSize[1] << " "<< CryoSize[2] << " "<< CryoSize[3] << " "<< CryoSize[4] << " "<< CryoSize[5] << std::endl;
      for (unsigned int t=0; t<cryostat.NTPC(); t++) {
        geo::TPCID id;
        id.Cryostat=c;
        id.TPC=t;
        id.isValid=true;
        const geo::TPCGeo& tpc=cryostat.TPC(id);
        TPCMinX[t] = tpc.MinX();
        TPCMaxX[t] = tpc.MaxX();
        TPCMinY[t] = tpc.MinY();
        TPCMaxY[t] = tpc.MaxY();
        TPCMinZ[t] = tpc.MinZ();
        TPCMaxZ[t] = tpc.MaxZ();
        std::cout << t << "\t" << TPCMinX[t] << " " << TPCMaxX[t] << " " << TPCMinY[t] << " " << TPCMaxY[t] << " " << TPCMinZ[t] << " " << TPCMaxZ[t] << std::endl;
      }
    }
    fCryos -> Fill();
    */
    return;
  }

void evgen::MUSUN::endRun ( art::Run &  run )

virtual

Reimplemented from art::EDProducer.

Definition at line 500 of file MUSUN_module.cc.

  {
    std::cout << "\n\nNumber of muons = " << NEvents << std::endl;
    std::cout << "Mean muon energy = " << se / NEvents << " GeV" << std::endl;
    std::cout << "Mean zenith angle (theta) = " << st / NEvents << " degrees" << std::endl;
    std::cout << "Mean azimuthal angle (phi)= " << sp / NEvents << " degrees" << std::endl;
    std::cout << "Mean slant depth = " << sd / NEvents << " m w.e." << std::endl;
  }

void evgen::MUSUN::initialization ( double  theta1, double  theta2, double  phi1, double  phi2, int  figflag, double  s_hor, double  s_ver1, double  s_ver2, double &  FI  )

private

Definition at line 681 of file MUSUN_module.cc.

  {
    //
    //  Read in the data files
    //
    int lineNumber = 0, index = 0;
    char inputLine[10000];
    std::string fROOTfile;

    for (int a = 0; a < 121; ++a)
      for (int b = 0; b < 62; ++b)
        for (int c = 0; c < 50; ++c)
          spmu[a][b][c] = 0;
    for (int a = 0; a < 23401; ++a)
      fnmu[a] = 0;
    for (int a = 0; a < 360; ++a)
      for (int b = 0; b < 91; ++b)
        depth[a][b] = 0;
    for (int a = 0; a < 360; ++a)
      for (int b = 0; b < 91; ++b)
        fmu[a][b] = 0;

    std::ostringstream File1LocStream;
    File1LocStream << fInputDir << fInputFile1;
    std::string File1Loc = File1LocStream.str();
    cet::search_path sp1("FW_SEARCH_PATH");
    if (sp1.find_file(fInputFile1, fROOTfile)) File1Loc = fROOTfile;
    std::ifstream file1(File1Loc.c_str(), std::ios::in);
    if (!file1.good())
      throw cet::exception("MUSUNGen")
        << "\nFile1 " << fInputFile1 << " not found in FW_SEARCH_PATH or at " << fInputDir
        << "\n\n";

    while (file1.good()) {
      //std::cout << "Looking at file 1...." << std::endl;
      file1.getline(inputLine, 9999);
      char* token;
      token = strtok(inputLine, " ");
      while (token != NULL) {
        //std::cout << "While loop file 1..." << std::endl;
        fmu[index][lineNumber] = atof(token);
        token = strtok(NULL, " ");
        index++;
        if (index == 360) {
          //std::cout << "If statement file 1..." << std::endl;
          index = 0;
          lineNumber++;
        }
      }
    }
    file1.close();

    std::ostringstream File2LocStream;
    File2LocStream << fInputDir << fInputFile2;
    std::string File2Loc = File2LocStream.str();
    cet::search_path sp2("FW_SEARCH_PATH");
    if (sp2.find_file(fInputFile2, fROOTfile)) File2Loc = fROOTfile;
    std::ifstream file2(File2Loc.c_str(), std::ios::binary | std::ios::in);
    if (!file2.good())
      throw cet::exception("MUSUNGen")
        << "\nFile2 " << fInputFile2 << " not found in FW_SEARCH_PATH or at " << fInputDir
        << "\n\n";

    int i1 = 0, i2 = 0, i3 = 0;
    float readVal;
    while (file2.good()) {
      //std::cout << "Looking at file 2...." << std::endl;
      file2.read((char*)(&readVal), sizeof(float));
      spmu[i1][i2][i3] = readVal;
      i1++;
      if (i1 == 121) {
        //std::cout << "First if statement file 2..." << std::endl;
        i2++;
        i1 = 0;
      }
      if (i2 == 62) {
        //std::cout << "Second if statement file 2..." << std::endl;
        i3++;
        i2 = 0;
      }
    }
    file2.close();
    for (int i = 0; i < 120; i++)
      for (int j = 0; j < 62; j++)
        for (int k = 0; k < 51; k++)
          spmu[i][j][k] = spmu[i + 1][j][k];
    spmu[1][1][0] = 0.000853544;
    //std::cout << "Set spmu to some value..." << std::endl;

    std::ostringstream File3LocStream;
    File3LocStream << fInputDir << fInputFile3;
    std::string File3Loc = File3LocStream.str();
    cet::search_path sp3("FW_SEARCH_PATH");
    if (sp3.find_file(fInputFile3, fROOTfile)) File3Loc = fROOTfile;
    std::ifstream file3(File3Loc.c_str(), std::ios::in);
    if (!file3.good())
      throw cet::exception("MUSUNGen")
        << "\nFile3 " << fInputFile3 << " not found in FW_SEARCH_PATH or at " << fInputDir
        << "\n\n";

    lineNumber = index = 0;
    while (file3.good()) {
      //std::cout << "Looking at file 3...." << std::endl;
      file3.getline(inputLine, 9999);
      char* token;
      token = strtok(inputLine, " ");
      while (token != NULL) {
        //std::cout << "While loop file 3..." << std::endl;
        depth[index][lineNumber] = atof(token);
        token = strtok(NULL, " ");
        index++;
        if (index == 360) {
          //std::cout << "If statement file 3..." << std::endl;
          index = 0;
          lineNumber++;
        }
      }
    }
    file3.close();

    //
    //  Set up variables
    //

    the1 = theta1;
    the2 = theta2;
    // for c2: c1 and c2 are unused
    //double c1 = cos(M_PI/180.*theta1);
    //double c2 = cos(M_PI/180.*theta2);
    ph1 = M_PI / 180. * phi1;
    ph2 = M_PI / 180. * phi2;
    // for c2: dph is unused
    //double dph = ph2-ph1;

    int ipc = 1;
    double theta = theta1;
    double dc = 1.;
    double sc = 0.;
    int iteration = 0;
    while (theta < theta2 - dc / 2.) {
      theta += dc / 2.;
      double theta0 = M_PI / 180. * theta;
      double cc = cos(theta0);
      double ash = s_hor * cc;
      double asv01 = s_ver1 * sqrt(1. - cc * cc);
      double asv02 = s_ver2 * sqrt(1. - cc * cc);
      int ic1 = (theta + 0.999);
      int ic2 = ic1 + 1;
      if (ic2 > 91) ic2 = 91;
      if (ic1 < 1) ic1 = 1;
      double phi = phi1;
      double dp = 1.;

      while (phi < phi2 - dp / 2.) {
        phi += dp / 2.;
        //  the long side of the cavern is pointing at 14 deg to the north:
        //  double phi0 = M_PI / 180. * (phi + 90. - 14);

        //  Want our co-ord system going from East to South.
        double phi0 = M_PI / 180. * (phi + fCavernAngle);

        double asv1 = asv01 * fabs(cos(phi0));
        double asv2 = asv02 * fabs(sin(phi0));
        double asv0 = ash + asv1 + asv2;
        double fl = 1.;
        if (figflag == 1) fl = asv0;
        int ip1 = (phi + 0.999);
        int ip2 = ip1 + 1;
        if (ip2 > 360) ip2 = 1;
        if (ip1 < 1) ip1 = 360;
        double sp1 = 0.;

        for (int ii = 0; ii < 4; ii++) {
          int iic = ii / 2;
          int iip = ii % 2;
          if (ip1 == 360 && (ii == 1 || ii == 3)) iip = -359;
          if (fmu[ip1 + iip - 1][ic1 + iic - 1] < 0) {
            if (pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 > 1e-6) {
              std::cout << "Looking at fmu [ " << ip1 << " + " << iip << " - 1 (" << ip1 + iip - 1
                        << ") ] [ " << ic1 << " + " << iic << " - 1 (" << ic1 + iic - 1 << ") ] ."
                        << "\nChanging sp1 from " << sp1 << " to "
                        << sp1 + pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4 << "..........."
                        << sp1 << " + 10 ^ (" << fmu[ip1 + iip - 1][ic1 + iic - 1] << ") / 4 "
                        << std::endl;
            }
            sp1 = sp1 + pow(10., fmu[ip1 + iip - 1][ic1 + iic - 1]) / 4;
          }
        }
        /*
              std::cout << iteration<< " time of new sc value! Theta " << theta << ", phi " << phi + dp / 2. << ", sc = " << sc + sp1 * fl * dp * M_PI / 180. * sin(theta0) * dc * M_PI / 180. << " = "
              << sc << " + " << sp1 << " * " << fl << " * " << dp << " * " << M_PI/180 << " * sin(" << theta0 << ") * " << dc << " * " << M_PI/180 << ".....sin(theta)=" << sin(theta) << "\n"
              << std::endl; */
        sc = sc + sp1 * fl * dp * M_PI / 180. * sin(theta0) * dc * M_PI / 180.;
        ++iteration;
        ipc = ipc + 1;
        fnmu[ipc - 1] = sc;
        phi = phi + dp / 2.;
      }

      theta = theta + dc / 2.;
    }
    //std::cout << *FI << " = " << sc << std::endl;
    FI = sc;
    for (int ipc1 = 0; ipc1 < ipc; ipc1++)
      fnmu[ipc1] = fnmu[ipc1] / fnmu[ipc - 1];
  }

void evgen::MUSUN::produce ( art::Event &  evt )

virtual

unique_ptr allows ownership to be transferred to the art::Event after the put statement

Implements art::EDProducer.

Definition at line 512 of file MUSUN_module.cc.

  {
    ///unique_ptr allows ownership to be transferred to the art::Event after the put statement
    std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);

    simb::MCTruth truth;
    truth.SetOrigin(simb::kSingleParticle);

    ++NEvents;

    SampleOne(NEvents, truth, fEngine);

    MF_LOG_DEBUG("MUSUN") << truth;

    truthcol->push_back(truth);

    evt.put(std::move(truthcol));
  }

void evgen::MUSUN::SampleOne ( unsigned int  i, simb::MCTruth &  mct, CLHEP::HepRandomEngine &  engine  )

private

Definition at line 536 of file MUSUN_module.cc.

  {
    CLHEP::RandFlat flat(engine);
    CLHEP::RandGaussQ gauss(engine);

    Energy = 0;
    theta = 0;
    phi = 0;
    dep = 0;
    Time = 0;

    sampling(Energy, theta, phi, dep, engine);

    theta = theta * M_PI / 180;

    //  changing the angle phi so z-axis is positioned along the long side
    //  of the cavern pointing at 14 deg from the North to the East.
    //  phi += (90. - 14.0);
    //  Want our co-ord rotation going from East to South.
    phi += fCavernAngle;
    if (phi >= 360.) phi -= 360.;
    if (phi < 0) phi += 360.;
    phi *= M_PI / 180.;

    // set track id to -i as these are all primary particles and have id <= 0
    int trackid = -1 * (i + 1);
    std::string primary("primary");

    // Work out whether particle/antiparticle, and mass...
    double m = 0.0;
    PdgCode = fPDG;
    double ChargeCheck = 1. / (1 + fChargeRatio);
    double pdgfire = flat.fire();
    if (pdgfire < ChargeCheck) PdgCode = -PdgCode;

    static TDatabasePDG pdgt;
    TParticlePDG* pdgp = pdgt.GetParticle(PdgCode);
    if (pdgp) m = pdgp->Mass();

    //std::cout << pdgfire << " " << ChargeCheck << " " << PdgCode << " " << m << std::endl;

    // Work out T0...
    if (fTDist == kGAUS) { Time = gauss.fire(fT0, fSigmaT); }
    else {
      Time = fT0 + fSigmaT * (2.0 * flat.fire() - 1.0);
    }

    //  The minus sign above is for y-axis pointing up, so the y-momentum
    //  is always pointing down
    cx = +sin(theta) * sin(phi);
    cy = -cos(theta);
    cz = +sin(theta) * cos(phi);
    Momentum = std::sqrt(Energy * Energy - m * m); // Get momentum
    px0 = Momentum * cx;
    py0 = Momentum * cy;
    pz0 = Momentum * cz;
    TLorentzVector pvec(px0, py0, pz0, Energy);

    //  Muon coordinates
    double sh1 = s_hor * cos(theta);
    double sv1 = s_ver1 * sin(theta) * fabs(cos(phi));
    double sv2 = s_ver2 * sin(theta) * fabs(sin(phi));
    double ss = sh1 + sv1 + sv2;
    double xfl1 = flat.fire();
    if (xfl1 <= sh1 / ss) {
      x0 = (fXmax - fXmin) * flat.fire() + fXmin;
      y0 = fYmax;
      z0 = (fZmax - fZmin) * flat.fire() + fZmin;
    }
    else if (xfl1 <= (sh1 + sv1) / ss) {
      x0 = (fXmax - fXmin) * flat.fire() + fXmin;
      y0 = (fYmax - fYmin) * flat.fire() + fYmin;
      if (cz >= 0)
        z0 = fZmin;
      else
        z0 = fZmax;
    }
    else {
      if (cx >= 0)
        x0 = fXmin;
      else
        x0 = fXmax;
      y0 = (fYmax - fYmin) * flat.fire() + fYmin;
      z0 = (fZmax - fZmin) * flat.fire() + fZmin;
    }
    // Make Lorentz vector for x and t....
    TLorentzVector pos(x0, y0, z0, Time);

    //  Parameters written to the file muons_surf_v2_test*.dat
    //      nmu - muon sequential number
    //      id_part - muon charge (10 - positive, 11 - negative )
    //      Energy  - total muon energy in GeV assuming ultrarelativistic muons
    //      x0, y0, z0 - muon coordinates on the surface of parallelepiped
    //          specified above; x-axis and y-axis are pointing in the
    //          directions such that the angle phi (from the slant depth
    //          distribution files) is measured from x to y. z-axis is
    //          pointing upwards.
    //      cx, cy, cz - direction cosines.

    simb::MCParticle part(trackid, PdgCode, primary);
    part.AddTrajectoryPoint(pos, pvec);

    mct.Add(part);

    theta = theta * 180 / M_PI;
    phi = phi * 180 / M_PI;

    // Sum energies, angles, depth for average outputs.
    se += Energy;
    st += theta;
    sp += phi;
    sd += dep;

    // Fill Histograms.....
    /*
    hPDGCode      ->Fill (PdgCode);
    hPositionX    ->Fill (x0);
    hPositionY    ->Fill (y0);
    hPositionZ    ->Fill (z0);
    hTime         ->Fill (Time);
    hMomentumHigh ->Fill (Momentum);
    hMomentum     ->Fill (Momentum);
    hEnergyHigh   ->Fill (Energy);
    hEnergy       ->Fill (Energy);
    hDepth        ->Fill (dep);
    hDirCosineX   ->Fill (cx);
    hDirCosineY   ->Fill (cy);
    hDirCosineZ   ->Fill (cz);
    hTheta        ->Fill (theta);
    hPhi          ->Fill (phi);
    */
    /*
    // Write event by event outsputs.....
    std::cout << "Theta: " << theta << " Phi: " << phi << " Energy: " << Energy << " Momentum: " << Momentum << std::endl;
    std::cout << "x: " <<  pos.X() << " y: " << pos.Y() << " z: " << pos.Z() << " time: " << pos.T() << std::endl;
    std::cout << "Px: " <<  pvec.Px() << " Py: " << pvec.Py() << " Pz: " << pvec.Pz() << std::endl;
    std::cout << "Normalised..." << cx << " " << cy << " " << cz << std::endl;
    */
    fTree->Fill();
  }

void evgen::MUSUN::sampling ( double &  E, double &  theta, double &  phi, double &  dep, CLHEP::HepRandomEngine &  engine  )

private

Definition at line 900 of file MUSUN_module.cc.

  {
    CLHEP::RandFlat flat(engine);
    CLHEP::RandGaussQ gauss(engine);

#if 0 // this code is disabled for good
    double xfl = flat.fire();
    int loIndex = 0, hiIndex = 32400;
    int i = (loIndex+hiIndex)/2;
    bool foundIndex = false;
    if( xfl < fnmu[loIndex] ) {
      i = loIndex;
      foundIndex = true;
    } else if ( xfl > fnmu[hiIndex] ) {
      i = hiIndex;
      foundIndex = true;
    } else if ( xfl > fnmu[i-1] && xfl <= fnmu[i] )
      foundIndex = true;
    while( !foundIndex ) {
      if( xfl < fnmu[i] )
        hiIndex = i;
      else
        loIndex = i;
      i = (loIndex + hiIndex)/2;

      if( xfl > fnmu[i-1] && xfl <= fnmu[i] )
        foundIndex = true;
    }
#else
    double xfl = flat.fire();
    int i = 0;
    while (xfl > fnmu[i])
      ++i;
#endif
    int ic = (i - 2) / 360;
    int ip = i - 2 - ic * 360;

    xfl = flat.fire();
    theta = the1 + ((double)ic + xfl);
    xfl = flat.fire();
    phi = ph1 + ((double)ip + xfl);
    if (phi > 360) phi = phi - 360;
    dep = depth[ip][ic] * fRockDensity;

    int ic1 = cos(M_PI / 180. * theta) * 50.;
    if (ic1 < 0) ic1 = 0;
    if (ic1 > 50) ic1 = 50;
    int ip1 = dep / 200. - 16;
    if (ip1 < 0) ip1 = 0;
    if (ip1 > 61) ip1 = 61;

    xfl = flat.fire();
#if 0
    loIndex = 0, hiIndex = 120;
    int j = (loIndex+hiIndex)/2;
    foundIndex = false;
    if( xfl < spmu[loIndex][ip1][ic1] ) {
      j = loIndex;
      foundIndex = true;
    } else if ( xfl > spmu[hiIndex][ip1][ic1] ) {
      j = hiIndex;
      foundIndex = true;
    } else if ( xfl > spmu[j-1][ip1][ic1] && xfl <= spmu[j][ip1][ic1] )
      foundIndex = true;
    while( !foundIndex ) {
      if( xfl < spmu[j][ip1][ic1] )
        hiIndex = j;
      else
        loIndex = j;
      j = (loIndex + hiIndex)/2;

      if( xfl > spmu[j-1][ip1][ic1] && xfl <= spmu[j][ip1][ic1] )
        foundIndex = true;
    }
#else
    int j = 0;
    while (xfl > spmu[j][ip1][ic1])
      ++j;
#endif

    double En1 = 0.05 * (j - 1);
    double En2 = 0.05 * (j);
    E = pow(10., En1 + (En2 - En1) * flat.fire());

    return;
  }

Member Data Documentation

double evgen::MUSUN::cx

private

Definition at line 279 of file MUSUN_module.cc.

double evgen::MUSUN::cy

private

Definition at line 279 of file MUSUN_module.cc.

double evgen::MUSUN::cz

private

Definition at line 279 of file MUSUN_module.cc.

double evgen::MUSUN::dep

private

Definition at line 277 of file MUSUN_module.cc.

double evgen::MUSUN::depth[360][91]

private

Definition at line 289 of file MUSUN_module.cc.

double evgen::MUSUN::Energy

private

Definition at line 277 of file MUSUN_module.cc.

double evgen::MUSUN::fCavernAngle

private

Angle of the detector from the North to the East.

Definition at line 232 of file MUSUN_module.cc.

double evgen::MUSUN::fChargeRatio

private

Charge ratio of particle / anti-particle.

Definition at line 225 of file MUSUN_module.cc.

double evgen::MUSUN::fEmax

private

Maximum Kinetic Energy (GeV)

Definition at line 239 of file MUSUN_module.cc.

double evgen::MUSUN::fEmin

private

Minimum Kinetic Energy (GeV)

Definition at line 238 of file MUSUN_module.cc.

CLHEP::HepRandomEngine& evgen::MUSUN::fEngine

private

art-managed random-number engine

Definition at line 222 of file MUSUN_module.cc.

double evgen::MUSUN::FI = 0.

private

Definition at line 282 of file MUSUN_module.cc.

int evgen::MUSUN::figflag

private

If want sampled from sphere or parallelepiped.

Definition at line 246 of file MUSUN_module.cc.

std::string evgen::MUSUN::fInputDir

private

Input Directory.

Definition at line 227 of file MUSUN_module.cc.

std::string evgen::MUSUN::fInputFile1

private

Input File 1.

Definition at line 228 of file MUSUN_module.cc.

std::string evgen::MUSUN::fInputFile2

private

Input File 2.

Definition at line 229 of file MUSUN_module.cc.

std::string evgen::MUSUN::fInputFile3

private

Input File 3.

Definition at line 230 of file MUSUN_module.cc.

double evgen::MUSUN::fmu[360][91]

private

Definition at line 290 of file MUSUN_module.cc.

double evgen::MUSUN::fnmu[32401]

private

Definition at line 288 of file MUSUN_module.cc.

int evgen::MUSUN::fPDG

private

PDG code of particles to generate.

Definition at line 224 of file MUSUN_module.cc.

double evgen::MUSUN::fPhimax

private

Maximum phi.

Definition at line 244 of file MUSUN_module.cc.

double evgen::MUSUN::fPhimin

private

Minimum phi.

Definition at line 243 of file MUSUN_module.cc.

double evgen::MUSUN::fRockDensity

private

Default rock density is 2.70 g cm-3. If this is changed then the three input files need to be remade. If there is a desire for this contact Vitaly Kudryavtsev at V.Kud.nosp@m.ryav.nosp@m.tsev@.nosp@m.shef.nosp@m..ac.u.nosp@m.k

Definition at line 233 of file MUSUN_module.cc.

double evgen::MUSUN::fSigmaT

private

Variation in t position (ns)

Definition at line 255 of file MUSUN_module.cc.

double evgen::MUSUN::fT0

private

Central t position (ns) in world coordinates.

Definition at line 254 of file MUSUN_module.cc.

int evgen::MUSUN::fTDist

private

How to distribute t (gaus, or uniform)

Definition at line 256 of file MUSUN_module.cc.

double evgen::MUSUN::fThetamax

private

Maximum theta.

Definition at line 242 of file MUSUN_module.cc.

double evgen::MUSUN::fThetamin

private

Minimum theta.

Definition at line 241 of file MUSUN_module.cc.

TTree* evgen::MUSUN::fTree

private

Definition at line 302 of file MUSUN_module.cc.

double evgen::MUSUN::fXmax

private

Maximum X position.

Definition at line 250 of file MUSUN_module.cc.

double evgen::MUSUN::fXmin

private

Minimum X position.

Definition at line 247 of file MUSUN_module.cc.

double evgen::MUSUN::fYmax

private

Maximum Y position.

Definition at line 251 of file MUSUN_module.cc.

double evgen::MUSUN::fYmin

private

Minimum Y position.

Definition at line 248 of file MUSUN_module.cc.

double evgen::MUSUN::fZmax

private

Maximum Z position.

Definition at line 252 of file MUSUN_module.cc.

double evgen::MUSUN::fZmin

private

Minimum Z position.

Definition at line 249 of file MUSUN_module.cc.

const int evgen::MUSUN::kGAUS = 1

staticprivate

Definition at line 220 of file MUSUN_module.cc.

double evgen::MUSUN::Momentum

private

Definition at line 278 of file MUSUN_module.cc.

unsigned int evgen::MUSUN::NEvents = 0

private

Definition at line 299 of file MUSUN_module.cc.

int evgen::MUSUN::PdgCode

private

Definition at line 276 of file MUSUN_module.cc.

double evgen::MUSUN::ph1

private

Definition at line 293 of file MUSUN_module.cc.

double evgen::MUSUN::ph2

private

Definition at line 293 of file MUSUN_module.cc.

double evgen::MUSUN::phi

private

Definition at line 277 of file MUSUN_module.cc.

double evgen::MUSUN::px0

private

Definition at line 278 of file MUSUN_module.cc.

double evgen::MUSUN::py0

private

Definition at line 278 of file MUSUN_module.cc.

double evgen::MUSUN::pz0

private

Definition at line 278 of file MUSUN_module.cc.

double evgen::MUSUN::s_hor = 0.

private

Definition at line 283 of file MUSUN_module.cc.

double evgen::MUSUN::s_ver1 = 0.

private

Definition at line 284 of file MUSUN_module.cc.

double evgen::MUSUN::s_ver2 = 0.

private

Definition at line 285 of file MUSUN_module.cc.

double evgen::MUSUN::sd = 0.

private

Definition at line 297 of file MUSUN_module.cc.

double evgen::MUSUN::se = 0.

private

Definition at line 294 of file MUSUN_module.cc.

double evgen::MUSUN::sp = 0.

private

Definition at line 296 of file MUSUN_module.cc.

double evgen::MUSUN::spmu[121][62][51]

private

Definition at line 287 of file MUSUN_module.cc.

double evgen::MUSUN::st = 0.

private

Definition at line 295 of file MUSUN_module.cc.

double evgen::MUSUN::the1

private

Definition at line 293 of file MUSUN_module.cc.

double evgen::MUSUN::the2

private

Definition at line 293 of file MUSUN_module.cc.

double evgen::MUSUN::theta

private

Definition at line 277 of file MUSUN_module.cc.

double evgen::MUSUN::Time

private

Definition at line 277 of file MUSUN_module.cc.

double evgen::MUSUN::x0

private

Definition at line 279 of file MUSUN_module.cc.

double evgen::MUSUN::y0

private

Definition at line 279 of file MUSUN_module.cc.

double evgen::MUSUN::z0

private

Definition at line 279 of file MUSUN_module.cc.

---

The documentation for this class was generated from the following file:

• larsim/larsim/EventGenerator/MuonPropagation/MUSUN_module.cc