MpSpallation.cc

Go to the documentation of this file.

#include "MuonProphet.h"


int MuonProphet::spallationProducts(HepMC::GenEvent& event, int idx)
{
  if(m_muon == 0) {
    error()<<"No muon track found"<<endreq;
    // Since no muon track, it should not reach this point.
    return -1;
  }
  
  HepMC::ThreeVector bkgPoint;
  bool bkgGenerated = false;
  bkgGenerated = generateOneBkg(idx, 
                                bkgPoint);

  if( !bkgGenerated ) return 0;

  HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;
  if(m_genTools[idx]->mutate(*bkgEvt).isFailure()) {
    error()<<"Failed to run background gentool "<<m_genToolNames[idx]<<endreq;
    return -1;
  }
  
  // debug
#if 0
  {
    HepMC::GenEvent::particle_iterator p;
    HepMC::GenEvent::particle_iterator particleEnd = bkgEvt->particles_end();
    for ( p = bkgEvt->particles_begin(); p!=particleEnd; ++p) {
      debug()<<*(*p)<<endreq;
    }
    HepMC::GenEvent::vertex_iterator v;
    HepMC::GenEvent::vertex_iterator vertexEnd = bkgEvt->vertices_end();
    for ( v = bkgEvt->vertices_begin(); v!=vertexEnd; ++v) {
      debug()<<*(*v)<<endreq;
    }
  }
#endif

  // change the vertex position first
  setPosition(*bkgEvt,bkgPoint);

  // change the time of the vertex
  // Get muon vertex time 
  double t0 = m_muon->production_vertex()->position().t();
  // to do: A delta should be added to t0 to account for flight time of muon
  //        and flight time of muon's daughter particle which reached the 
  //        production vertex of the background.
  
  setTime(*bkgEvt,t0,m_genLifetimes[idx]);

  // 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 0;
}

bool MuonProphet::generateOneBkg(int idx,  /* used to find the property of one background */
                                 HepMC::ThreeVector& bkgPoint)
{
  if( !m_crossWater ) {
    return false;
  }

  double eMuon = m_muon->momentum().e();

  // Mineral oil and LS densities are assumed to be the same.
  double density = m_mineralOil->density();
  double yield = m_genYields[idx] ;
  double alpha = 0.77;

  double yieldAtThisEnergy = yield * pow( (eMuon / m_genYieldMeasuredAt[idx]), alpha);

  double ProbThres = yieldAtThisEnergy * density * m_trkLengthInWater;

  debug()<<m_genToolNames[idx]<<" eMuon= "<<eMuon/CLHEP::GeV<<"GeV "
         <<" trkLengthInWater= "<<m_trkLengthInWater/CLHEP::meter<<"m "
         <<"ProbThres= "<<ProbThres<<endreq;

  if( ProbThres > 1 || ProbThres < 0 ) {
    debug()<<"Yield too high or too low. Invalid probability "<<ProbThres<<endreq;
  }

  if( m_rnd.Rndm() > ProbThres) {
    return false;
  }

  
  // The position of background in the cylinder coordinate system of muon track
  // Muon direction is z.
  double r,phi,z;

  z = m_rnd.Rndm() * (m_tickWaterMax - m_tickWaterMin) + m_tickWaterMin;
  r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
  phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;


  // try point in coordinate system of intersect 1st
  HepGeom::Point3D<double> tryPoint;

  tryPoint.setX( r*cos(phi) );
  tryPoint.setY( r*sin(phi) );
  tryPoint.setZ( z );

  /* This is a muon track shown below.
   *
   *            intersect 1st (tick1)  intersect 2nd (tick2)
   * Primary vertex -----1--------C----------2------>
   *                        z     | theta
   *                              |  
   *                              | r
   *                              |
   *                              |
   *                           tryPoint
   *
   */

  // direction of muon
  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);

  // tryPoint in ows coordinate system
  HepGeom::Point3D<double> tryPointGlobal = muonDir * tryPoint.r();

  // rotation
  tryPointGlobal.rotate(tryPoint.theta(), y);
  tryPointGlobal.rotate(tryPoint.phi(),   muonDir);

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

  debug()<<tryPoint<<endreq;
  debug()<<tryPointGlobal<<endreq;

  Gaudi::XYZPoint p(tryPointGlobal.x(),
                    tryPointGlobal.y(),
                    tryPointGlobal.z());
  
  bool inside = false;
  for (unsigned int idx = 0; idx<m_ADs.size(); ++idx) 
    {
      if(m_ADs[idx]->isInside(p)) inside = true;
    }
  
  if(!inside) return false;
  
  // done
  bkgPoint.setX(tryPointGlobal.x());
  bkgPoint.setY(tryPointGlobal.y());
  bkgPoint.setZ(tryPointGlobal.z());

  return true;
}