MuonProphet.h

Go to the documentation of this file.

#ifndef _MUON_PROPHET_H_
#define _MUON_PROPHET_H_ 1

#include "MuonProphet/MpMuonFate.h"
#include "MuonProphet/MpTriggerStat.h"
#include "RectangularBox.h"
#include "RpcPlane.h"
#include "DetDesc/Material.h"
#include "DetDesc/ILVolume.h"
#include "DetDesc/IDetectorElement.h"
#include "DetDesc/IGeometryInfo.h"

#include "GaudiAlg/GaudiTool.h"
#include "GenTools/IHepMCEventMutator.h"
#include "HepMC/GenEvent.h"
#include "Conventions/Site.h"
#include "CLHEP/Units/SystemOfUnits.h"

// random
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/IRndmGenSvc.h"

// root
#include "TF1.h"
#include "TRandom3.h"

using namespace std;

class MuonProphet : public GaudiTool,
                    virtual public IHepMCEventMutator
{
 public:

  MuonProphet(const std::string& type,
              const std::string& name,
              const IInterface* parent);
  virtual ~MuonProphet();
  
  virtual StatusCode initialize();
  virtual StatusCode finalize();
  
  // HepMCEventMutator interface
  virtual StatusCode mutate(HepMC::GenEvent& event);
  
 private:
  // Predict the fate of a muon, alter the status of the muon track.
  // At this moment, fast simulation is enforced to all muon track.
  MpMuonFate predictMuonFate(HepMC::GenEvent& event);
  // get all geometry related parameters
  StatusCode geometryCalculationInit();
  StatusCode geometryCalculation(HepMC::GenParticle * pMuon);

  // Triggered by RPC?
  MpTriggerStat triggeredByRpc(HepMC::GenParticle * pMuon);
  // Triggered by water pool?
  MpTriggerStat triggeredByWaterPool(HepMC::GenParticle * pMuon);
  // to do: I need a more complicate function to show the 
  // relation between a muon track and water pool.
  // path length in water pool
  // double pathLengthInWaterPool();

  // For a given muon track, add spallation products, return number of tracks added,
  // and negtive for failure.
  int spallationProducts(HepMC::GenEvent& event, int idx);
  bool generateOneBkg(int idx,  /* used to find the property of one background */
                      HepMC::ThreeVector& bkgPoint);

  StatusCode spallationNeutronsInit();
  int spallationNeutrons(HepMC::GenEvent& event);
  int generateNeutrons(HepMC::GenEvent *bkgEvt,
                       ISolid::Tick tickBegin, // In fact it is a double, indicate the first intersection with the volume
                       ISolid::Tick tickEnd,   // Indecate the last intersection with the volume. See m_point, m_vec and m_unit
                       double multiplicity);

  StatusCode printout(HepMC::GenEvent& event); // printout all vertex and particle info
  // Put a spallation background to position
  void setPosition(HepMC::GenEvent& event, HepMC::ThreeVector& position);
  // Set a spallation background to right time
  void setTime(HepMC::GenEvent& event, double t0, double lifetime);


  // A switch to turn on or off the whole module
  bool m_active;

  // My muon track
  HepMC::GenParticle * m_muon;

  // name of them
  vector<std::string>         m_genToolNames;
  // id for them. An ID for each of them for easier and faster operation
  vector<int>                 m_genId;
  // the genTool (generator)
  vector<IHepMCEventMutator*> m_genTools;
  /*
   *  About spallation background yield, I take the convention from KamLAND.
   *  arXiv:0907.0066
   *  The first set of parameters are also from it.
   *
   *  Yield = N_bkg / (Muon_Rate * Detector_Live_Time * LS_density * Track_length)  
   *  [Yield] = 1/[g/cm^2]
   *
   *  The yield is also assumed to obey a power-law function of muon energy.
   *  (E_mu)^alpha
   *  alpha = 0.77 in the first version
   *  
   *  15 Oct. 2009
   */
  // Spallation background yield
  vector<double>              m_genYields;
  // The energy at which the yield is measured
  vector<double>              m_genYieldMeasuredAt;
  // Lifetime for each background
  vector<double>              m_genLifetimes;


  // turn on or off neutron production
  bool m_actNeutron;
  // Neutron must have its own generator. 
  // Neutron is special. They should be generated everywhere and with initial momentum.
  
  // Neutron yield. Normally it will pick a NeutronYield function.
  // However sometimes it is very convinient to set this manually.
  // <0:  use Neutron Yield function (default)
  // >=0: use the value of this parameter.
  double m_neutronYield;


  // Site
  // See Conventions/Site.h for conventions
  string m_siteName;
  Site::Site_t m_site;

  TRandom3 m_rnd;

  RectangularBox * m_waterPool;

  string m_topDetectorElementName;
  IGeometryInfo * m_topGeoInfo;
  ILVolume * m_topLVolume;
  vector<IGeometryInfo *> m_ADs;

  ILVolume::Intersections m_intersections;
  ILVolume::Intersections m_crucialSegments;

  Material * m_mixGas;      
  Material * m_owsWater;    
  Material * m_iwsWater;    
  Material * m_rock;        
  Material * m_mineralOil;  

  // m_point primary vertex
  // m_vec   momentum vector
  // m_unit  unit vector of momentum direction
  Gaudi::XYZPoint   m_point;
  Gaudi::XYZVector  m_vec;
  Gaudi::XYZVector  m_unit;


  bool m_crossWater;
  /*  Not in use
  HepGeom::Point3D<double> m_crossWaterA;
  HepGeom::Point3D<double> m_crossWaterB;
  */  

  // Track length in water
  double m_trkLengthInWater;
  // The value of t when muon track intersects with water pool. 
  // See m_point for conventions.
  ISolid::Tick  m_tickWaterMin;
  ISolid::Tick  m_tickWaterMax;

  // Longer than that will be considered to have high trigger rate
  // For example 20cm.
  double m_trkLengthInWaterThres;

  /*
  // Trigger efficiency in the high trigger rate
  // Even a track length is longer than 20 cm, there is still a chance
  // it will miss trigger. Like travel through some dead material,
  // or it is scattered before it enters water pool
  // to do: However the inefficiency is correlated with other system, not in use.
  */
  double m_waterPoolTriggerEff;

  RpcPlane * m_rpcPlane;                       
  // does a track intersect with rpc and if yes, where it is?
  bool m_crossRpc;
  HepGeom::Point3D<double> m_rpcIntersect;
  
  //
  MpTriggerStat m_poolTriggerStatus;
  MpTriggerStat m_rpcTriggerStatus;

  
  //

  TF1* m_nEnergyPDF;
  TF1* m_nAnglePDF;
  TF1* m_nLateralPDF;  
  
};


#endif // _MUON_PROPHET_H_