In This Package:

MpNeutron.cc

Go to the documentation of this file.

#include "MuonProphet.h"
#include "../functions/NeutronPDF.h"
#include "TRandom3.h"

StatusCode MuonProphet::spallationNeutronsInit()
{
  m_nEnergyPDF=0;
  m_nAnglePDF=0;
  m_nLateralPDF=0;


  // Avarage muon energy is from the estimation of Mengyun Guan. (Ph.D. thesis)
  double avaMuonE;
  if(m_site == Site::kDayaBay) {
    avaMuonE = 55.3*CLHEP::GeV;
  } else if (m_site == Site::kLingAo) {
    avaMuonE = 61.4*CLHEP::GeV;
  } else if (m_site == Site::kFar ) {
    avaMuonE = 140.3*CLHEP::GeV;
  }
  
  m_nEnergyPDF = new TF1("m_nEnergyPDF",NeutronEnergyPDF,0.001,4,1);
  m_nEnergyPDF->SetNpx(200);  

  m_nEnergyPDF->SetParameter(0,avaMuonE/CLHEP::GeV);
  if(!m_nEnergyPDF) {
    error()<<"Failed to get NeutronEnergyPDF"<<endreq;
    return StatusCode::FAILURE;
  }

  m_nAnglePDF = new TF1("m_nAnglePDF",NeutronAngularDis,-1,1,1);
  m_nAnglePDF->SetParameter(0,avaMuonE/CLHEP::GeV);
  if(!m_nAnglePDF) {
    error()<<"Failed to get NeutronAngularDis"<<endreq;
    return StatusCode::FAILURE;
  }

  
  m_nLateralPDF = new TF1("m_nLateralPDF",LateralProfile,0,100,1);
  m_nLateralPDF->SetParameter(0,0);
  if(!m_nLateralPDF) {
    error()<<"Failed to get LateralProfile"<<endreq;
    return StatusCode::FAILURE;
  }

  return StatusCode::SUCCESS;
}


int MuonProphet::spallationNeutrons(HepMC::GenEvent& event) 
{
  debug()<<"Spallation Neutron production"<<endreq;

  //--------------------------------------------------------
  unsigned int NSeg = m_crucialSegments.size();
  
  double yield;
  if(m_neutronYield>=0) {
    yield = m_neutronYield;
  } else {
    yield = NeutronYield( m_muon->momentum().e()/CLHEP::GeV )*(CLHEP::cm2/CLHEP::gram);
  }

  double tracklength;
  double multiplicity;

  int nNewTrack = 0;

  HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;

  for (unsigned int idx=0; idx<NSeg; ++idx) 
    {
      double prob;

      
      tracklength = m_crucialSegments[idx].first.second - m_crucialSegments[idx].first.first;
      multiplicity = NeutronMulti( m_muon->momentum().e()/CLHEP::GeV,
                                   tracklength );
      
      prob = yield
        * tracklength
        * m_crucialSegments[idx].second->density()
        / multiplicity;
      
      debug()<<"Yield= "<<NeutronYield( m_muon->momentum().e()/CLHEP::GeV )<<"cm2/g "
             <<"TrackLength= "<< tracklength/CLHEP::m <<"m "
             <<"Density= "<<m_crucialSegments[idx].second->density()/(CLHEP::g/CLHEP::cm3)<<"g/cm3 "
             <<"Multiplicity= "<< multiplicity <<endreq;
      debug()<<"prob= "<<prob<<endreq;

      

      if( m_rnd.Rndm() < prob ) 
        {
          nNewTrack += generateNeutrons(bkgEvt,
                                        m_crucialSegments[idx].first.first,   // tick one
                                        m_crucialSegments[idx].first.second,  // tick two
                                        multiplicity);
        }
    }  
  
  // drag vertices out from background event

  HepMC::GenVertex *bkgVtx;
  while( !bkgEvt->vertices_empty() )
    {
    bkgVtx = *(bkgEvt->vertices_begin());
    
    // add to muon event
    // This will automatically transfer the ownership from the old event to new event.
    // So the original event can be safely deleted now.
    event.add_vertex(bkgVtx);
    }
 
  // delete the bkgEvt
  delete bkgEvt;

  return nNewTrack;
  
}


// tickBegin:  In fact it is a double, indicate the first intersection with the volume
// tickEnd:    Indecate the last intersection with the volume. See m_point, m_vec and m_unit

int MuonProphet::generateNeutrons(HepMC::GenEvent *bkgEvt,
                                  ISolid::Tick tickBegin,
                                  ISolid::Tick tickEnd,
                                  double multiplicity)
{
  int nNeutron = m_rnd.Poisson(multiplicity);
  
  debug()<<"nNeutron= "<<nNeutron<<endreq;
  

  double z;
  double r;
  double phi;
  
  double ke;
  double eng;
  double momentum;
  double massN = 939.565346*CLHEP::MeV;
  double phiN;
  double costh;
  double sinth;

  for(int idx =0; idx<nNeutron; ++idx)  
    { 
      z = m_rnd.Rndm() * (tickEnd - tickBegin) + tickBegin;
      r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
      phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;

      ke = m_nEnergyPDF->GetRandom() * CLHEP::GeV ;      
      eng = ke + massN;
      momentum = sqrt (eng*eng - massN*massN);
      phiN = m_rnd.Rndm() * 2 * 3.14159265358979323846;
      costh = m_nAnglePDF->GetRandom();
      sinth = sqrt(1-costh*costh);


      debug()<<"z "<<z/CLHEP::m<<"m"<<endreq;
      debug()<<"r "<<r/CLHEP::cm<<"cm"<<endreq;
      debug()<<"ke "<<ke/CLHEP::GeV<<"GeV"<<endreq;
      debug()<<"costh "<<costh<<endreq;

      HepGeom::Point3D<double> genPoint;
      
      genPoint.setX( r*cos(phi) );
      genPoint.setY( r*sin(phi) );
      genPoint.setZ( z );

      HepGeom::Point3D<double> genVec;
      
      genVec.setX ( momentum*sinth*cos(phiN) );
      genVec.setY ( momentum*sinth*sin(phiN) );
      genVec.setZ ( momentum*costh );

      HepGeom::Point3D<double> muonDir(m_unit.X(),
                                       m_unit.Y(),
                                       m_unit.Z());

      // need a y axis which is perpendicular to muonSlope
      HepGeom::Point3D<double> y( m_unit.X(),
                                  m_unit.Y(),
                                  0);
      y.rotateZ(3.14159265358979323846/2);
      if(y.r()<=0) y.setX(1);
      
      HepGeom::Point3D<double> genPointGlobal = muonDir * genPoint.r();
      debug()<<m_unit<<endreq;
      debug()<<muonDir<<endreq;
      debug()<<genPointGlobal<<endreq;

      // rotation
      genPointGlobal.rotate(genPoint.theta(), y);
      genPointGlobal.rotate(genPoint.phi(),   muonDir);

      // transform
      genPointGlobal.setX(genPointGlobal.x() + m_point.X());
      genPointGlobal.setY(genPointGlobal.y() + m_point.Y());
      genPointGlobal.setZ(genPointGlobal.z() + m_point.Z());

      HepGeom::Point3D<double> genVectorGlobal = muonDir * genVec.r();
      
      genVectorGlobal.rotate(genVec.theta(),  y);
      genVectorGlobal.rotate(genVec.phi(),    muonDir);
      

      double t0 = m_muon->production_vertex()->position().t();
      debug()<<"t0= "<<t0<<endreq;
      t0 += tickEnd / (3e8 * CLHEP::m/CLHEP::s);
      debug()<<"dt= "<<tickEnd / (3e8 * CLHEP::m/CLHEP::s)<<endreq;
      debug()<<"t0= "<<t0<<endreq;
      HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(genPointGlobal.x(),
                                                                        genPointGlobal.y(),
                                                                        genPointGlobal.z(),
                                                                        t0));
      
      HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(genVectorGlobal.x(),
                                                                              genVectorGlobal.y(),
                                                                              genVectorGlobal.z(),
                                                                              eng),
                                                                              2112,
                                                                              1/*=status*/);
      vertex->add_particle_out(particle);
      bkgEvt->add_vertex(vertex);
    } 

  return nNeutron;
}

---