In This Package:

GtMuoneratorTool.cc

Go to the documentation of this file.

/* Generate cosmic-muons in a GenTool
 *   
 *   Brett Viren
 *   Wed May  5 16:47:41 2010
 *
 *   Disclaimer: this refactors the app from Generators/Muon into a
 *   tool. The essentials themselves is unchanged.
 */

#include "GtMuoneratorTool.h"

#include "GaudiKernel/IRndmGenSvc.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"

#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"

#include "TH1F.h"
#include "TFile.h"

#include <fstream>
#include <cassert>

using namespace std;
using namespace CLHEP;

GtMuoneratorTool::GtMuoneratorTool(const std::string& type,
                                   const std::string& name, 
                                   const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_pmratio(0)
    , m_Sx(0), m_Sy(0), m_Sz(0), m_Smax(0)
    , m_dphi(0)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("WhichSite",m_whichSite="",
                    "Which site (DYB, LA, Mid, Far,SAB)");
    declareProperty("Rotation",m_rotation=false,
                    "Set to true to assume rotation for the given site");
    declareProperty("MuonFile",m_muonFileName="",
                    "Name of file giving muons to sample");
    declareProperty("RatioFile",m_ratioFileName="",
                    "Name of file giving mu+/mu- ratio histogram");
    declareProperty("Volume",m_volume="rock",
                    "Volume name (rock, RPC, ADE)");
    info() << "Creating GtMuoneratorTool(\"" << type << "\", \"" << name << "\")" << endreq;
}


GtMuoneratorTool::~GtMuoneratorTool()
{
    info() << "Destroying GtMuoneratorTool site=" << m_whichSite
           << " MuonFile=" << m_muonFileName
           << " RatioFile=" << m_ratioFileName
           << " Volume=" << m_volume
           << endreq;
}

StatusCode GtMuoneratorTool::initialize()
{
    // Set up random numbers
    m_rgs = 0;
    if (service("RndmGenSvc",m_rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }
    if (m_uni.initialize(m_rgs, Rndm::Flat(0,1)).isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }

    // leak this file
    TFile* file = new TFile(m_ratioFileName.c_str(),"read");
    if (file->IsZombie()) {
        warning() << "Failed to open ROOT file \"" << m_ratioFileName << "\""
                  << " will fall back to parametrized ratio"
                  << endreq;
    }
    else {
        m_pmratio = (TH1F*)file->Get("mu_plus_minus_ratio");

        if (!m_pmratio) {
            error() << "Failed to get \"mu_plus_minus_ratio\" TH1F from ROOT file "
                    << m_ratioFileName << endreq;
            return StatusCode::FAILURE;
        }
    }

    return ReadMuons();
}

static int GetMuPlusMinusRatioParam(double momInGeV)
{
    Double_t par[2]={0.5505,0.6745};
    Double_t part1=par[0]/(1+(1.1*momInGeV+1)/115);
    Double_t part2=0.054*par[1]/(1+(1.1*momInGeV+1)/850);
    Double_t part3=(1-par[0])/(1+(1.1*momInGeV-1)/115);
    Double_t part4=0.054*(1-par[1])/(1+(1.1*momInGeV-1)/850);
    return (part1+part2)/(part3+part4);
}

double GtMuoneratorTool::GetMuPlusMinusRatio(double momInGeV)
{
    if (!m_pmratio) return GetMuPlusMinusRatioParam(momInGeV);

    int binno = m_pmratio->FindBin(momInGeV);
    //return the measured value if within range. otherwise return 1.4
    if(binno>=1&&binno<=m_pmratio->GetNbinsX()) 
        return m_pmratio->GetBinContent(binno);
    return 1.4;
}

int GtMuoneratorTool::GetPIDFromMomentum(double momInGeV)
{
    double ratio = GetMuPlusMinusRatio(momInGeV);
    double prob_mu_plus = ratio/(1.+ratio);
    //int number = gRandom->Binomial(1, prob_mu_plus);
    Rndm::Numbers bi;
    if (bi.initialize(m_rgs, Rndm::Binomial(1,prob_mu_plus)).isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }
    int number = bi();
    if (number==0) return 13; //mu-
    return -13;//mu+
}

StatusCode GtMuoneratorTool::ReadMuons()
{
    // djaffe 7mar07. (Corrected 12mar07. x-,y-dimensions were swapped.)
    // set dimensions and detector center based on CD1 release (doc-765,762,720)
    // for DYB, LA, Far and default site
    Double_t height_RPC = 0.50*CLHEP::m ; // RPCs are nominally 50cm above  surface of shield volume (includes 20cm of air above water)
    Double_t thick_RPC  = 0.056*CLHEP::m + 0.030*CLHEP::m ; // 1.4cm/layer*4layers + 1cm spacing*3spaces between layers
    Double_t overlap_RPC= 1.0*CLHEP::m  ; // nominal 1m overlap of RPC with rock around pool
    Double_t extra      = 0.05*CLHEP::m ; // 5cm extra space to make sure muons start outside of all active volumes
    Double_t width_near = 10.0*CLHEP::m ; // width of near pool
    Double_t width_far  = 16.0*CLHEP::m ; // width of far pool
    Double_t length     = 16.0*CLHEP::m ; // length of near and far pool
    Double_t depth      = 10.0*CLHEP::m ; // depth of near and far pool
    Double_t air_thick  =  0.2*CLHEP::m ; // 20cm of air above water
    Double_t height_hall= 15.0*CLHEP::m ; // 15m of air above water, from water surface to hall ceiling
    Double_t top_rock_extra = 2*CLHEP::m; // 2m of rock above hall ceiling
    Double_t side_rock_extra= 4*CLHEP::m; // this will make sure at least there is 2m rock of hall wall
    // for pool 4m of rock
    Double_t bottom_rock_extra = 2*CLHEP::m; // 2m of rock under the pool

    Double_t yDim ;  // depends on near,far
    Double_t xDim ;  // same for near,far
    Double_t zDim; // box height
    // the origin (0,0,0) is always at the center of water pool
    Double_t zOffset; // center of the volume

    // stainless steel
    // <!-- Radius for SST -->
    Double_t ADsstRadius=2500*CLHEP::mm;
    // <!-- Height for SST -->
    Double_t ADsstHeight=5000*CLHEP::mm;
    // AD envelop size
    // <!-- ADE extention beyond SST in radius -->
    Double_t ADadeWall=1.0*CLHEP::cm; // 0.25*m;
    // <!-- ADE head gap above tank -->
    Double_t ADadeHead=1.0*CLHEP::cm; // 1.0*m;
    // <!-- ADE foot gap below tank -->
    Double_t ADadeFoot=1.0*CLHEP::cm;
    // <!-- ADE radius -->
    Double_t ADadeRadius=ADsstRadius+ADadeWall;
    // <!-- ADE height -->
    Double_t ADadeHeight=ADadeFoot+ADsstHeight+ADadeHead;

    if (m_volume == "rock") {    
        //Z
        zDim=depth + height_hall + top_rock_extra + bottom_rock_extra;  // same for near, far;
        // the origin (0,0,0) is always at the center of water pool
        zOffset=(height_hall + top_rock_extra - bottom_rock_extra)/2;   // vertical offset of the certer of z
        //X
        xDim = length + 2.*overlap_RPC + 2.*side_rock_extra; // same for near,far
        //Y
        if (m_whichSite == "DYB" || m_whichSite == "LA") {
            yDim=width_near + 2.*overlap_RPC + 2.*side_rock_extra;
        }
        else if (m_whichSite == "Far") {
            yDim=width_far + 2.*overlap_RPC + 2.*side_rock_extra;
        }
        else {
            // fixme: error - unknown site
        }
    } 
    else if (m_volume == "RPC") {
        //Z 
        zDim=height_RPC + thick_RPC + depth + extra + air_thick;  // same for near,far
        // the origin (0,0,0) is always at the center of water pool
        zOffset=(height_RPC + thick_RPC + extra + air_thick)/2. ; // vertical offset of center of 'detector'
        //X
        xDim = length + 2.*overlap_RPC + 2.*extra;
        //Y
        if (m_whichSite == "DYB" || m_whichSite == "LA") {
            yDim=width_near + 2.*overlap_RPC + 2.*extra;
        }
        else if (m_whichSite == "Far") {
            yDim=width_far + 2.*overlap_RPC + 2.*extra;
        }
        else {
            // fixme: error - unknown site
        }
    }
    else if (m_volume == "ADE") {
        // Z
        zDim=ADadeHeight;
        // the origin (0,0,0) is always at the center of AD
        zOffset=0;
        // X
        xDim= 2*ADadeRadius;
        // Y
        yDim=xDim;
    }
    else {
        cout<<"muon_generator: *** ERROR *** Invalid volume " << m_volume
            << " selected. Valid volume names are rock, RPC and ADE"<<endl;
        assert(0); // this will end the job
    }

    Double_t dphi=0;           // detector angle.


    // get the right detector setting and the right muon flux file 
    ifstream infile;
    infile.open(m_muonFileName.c_str(), ios_base::in);
    if (!infile.is_open()) {
        error() << "Failed to open muon flux file \"" << m_muonFileName << "\"" 
                << endreq;
        return StatusCode::FAILURE;
    }
    debug() << "Opened file \"" << m_muonFileName << "\"" << endreq;

    if(m_whichSite == "DYB") {
        m_detectorDim = Hep3Vector(xDim, yDim, zDim);
        m_detectorCenter = Hep3Vector(0, 0, zOffset );
        if(m_rotation) {
            dphi = 56.7*degree;
        }
    }
    else if(m_whichSite == "CDR") {
        m_detectorDim = Hep3Vector(16.06*CLHEP::m, 10.06*CLHEP::m, 10.0*CLHEP::m );
        m_detectorCenter = Hep3Vector(0, 0, -5*CLHEP::m );
        dphi = 0; 
    }
    else if(m_whichSite == "LA") {
        m_detectorDim = Hep3Vector(xDim, yDim, zDim);
        m_detectorCenter = Hep3Vector(0, 0, zOffset );
        if(m_rotation) {
            dphi = 79.6*degree;
        }
    }
    else if(m_whichSite == "Mid") {
        error()<<" No implemented yet. "<<endl; 
        return StatusCode::FAILURE;
    }
    else if(m_whichSite == "Far") {
        m_detectorDim = Hep3Vector(xDim, yDim, zDim);
        m_detectorCenter = Hep3Vector(0, 0, zOffset );
        if(m_rotation) {
            dphi = 60.5*degree;
        }
    }
    else if(m_whichSite == "SAB") {
        m_detectorDim = Hep3Vector(5.2*CLHEP::m, 5.2*CLHEP::m, 5.2*CLHEP::m );
        m_detectorCenter = Hep3Vector(0, 0, 0 );
        dphi = 0; 
    }
    else {
        // 23mar07 djaffe If no valid site name, then STOP
        error() <<"muon_generator: *** ERROR *** Invalid site " << m_whichSite
                << " selected. Valid site names are Far, Mid, DYB, LA, CDR,SAB"<<endreq;
        return StatusCode::FAILURE;
    }

    // Calculate size of sides to generate on for later.

    // the detector geometry 
    Double_t dimx=0, dimy=0, dimz=0; 
    Double_t xdet=0, ydet=0, zdet=0; 
  
    dimx = m_detectorDim.x()/CLHEP::m; 
    dimy = m_detectorDim.y()/CLHEP::m; 
    dimz = m_detectorDim.z()/CLHEP::m;
    xdet = m_detectorCenter.x()/CLHEP::m; 
    ydet = m_detectorCenter.y()/CLHEP::m; 
    zdet = m_detectorCenter.z()/CLHEP::m; 
    cout<<"The detector dimension is: "
        <<dimx<<"*m, "<<dimy<<"*m, "<<dimz<<"*m"<<endl;
    cout<<"The detector center is: "
        <<xdet<<"*m, "<<ydet<<"*m, "<<zdet<<"*m"<<endl; 
    cout<<"The detector angle is: "<<dphi<<endl;
    cout<<"MAKE SURE the setting is consistent with Geant4 setting!!!"<<endl;

    m_Sz = dimx * dimy;
    m_Sx = dimy * dimz;
    m_Sy = dimx * dimz;
    m_Smax = sqrt(m_Sx*m_Sx + m_Sy*m_Sy + m_Sz*m_Sz);
    m_dphi = dphi;
    
    // read in muons

    //temporary variable to hold data in flux file.
    Double_t energy, phi, theta, energy0, phi0, theta0;  
  
    Double_t flux_of_site;
    // line 1 to line 6 are comments
    char tmp[255];
    for (int lineNo=0; lineNo<4; lineNo++) {      
        infile.getline(tmp, 255);      
    }
    // line 5 contains flux info of that site.
    infile.getline(tmp,255,':');
    infile>>flux_of_site;
    infile.getline(tmp,255);
    // line 6
    infile.getline(tmp,255);

    while (!infile.eof()) {
        infile>>energy>>theta>>phi>>energy0>>theta0>>phi0;   
        m_muonE.push_back(energy*CLHEP::GeV);
        m_muonPhi.push_back(phi);
        m_muonTheta.push_back(theta);
        verbose() << "Loading muon E=" << energy*CLHEP::GeV 
                  << " phi=" << phi
                  << " theta=" << theta 
                  << endreq;
    }

    m_muonE.pop_back();
    m_muonPhi.pop_back();
    m_muonTheta.pop_back();

    return StatusCode::SUCCESS;

} // ReadMuons()


StatusCode GtMuoneratorTool::finalize()
{
    return StatusCode::SUCCESS;
}

StatusCode GtMuoneratorTool::mutate(HepMC::GenEvent& event)
{
    double pS_z=0, pS_x=0, pS_y=0; // projected area in x, y, z direction.
    const double PI = 3.1415926;
    Hep3Vector direction;          // the muon direction.
    double muon_mass = 0.105658*CLHEP::GeV;
    double weight = 1.;
    double arandom=0;
    int ievent = 0;
  
    do {             // pick a muon       
        ievent = (int)(m_uni() * (m_muonE.size()-1));
        debug() << "Pick muon " << ievent << " out of " << m_muonE.size() << endreq;

        // convert coordinates from MUSIC system to digimap system.
        double theta = PI - PI * m_muonTheta[ievent]/180;
        double phi = PI/2 -PI * m_muonPhi[ievent]/180;
        if (phi < 0) phi = phi + 2 * PI;

        // The WC detector may be rotated an angle to better shield
        // neutrons, It is convenient to rotate muon direction instead
        // of rotating detector in Geant4
        phi = phi + m_dphi;
        if(phi > 2 * PI) phi = phi - 2 * PI;
        direction.setX(sin(theta) * cos(phi));
        direction.setY(sin(theta) * sin(phi));
        direction.setZ(cos(theta));
    
        // projected area of the WE surface to this muon.
        pS_z = m_Sz * abs(direction.z()); 
        pS_x = m_Sx * abs(direction.x()); 
        pS_y = m_Sy * abs(direction.y());

        weight = ( pS_z + pS_x + pS_y )/m_Smax; 
        arandom = m_uni();
    } while ( arandom > weight ); // drop this muon.
    
    // calculate the muon momentum and position.
    
    double eMuon = m_muonE[ievent];
    double pMuon = sqrt(eMuon*eMuon - muon_mass*muon_mass);
    Hep3Vector momentum;
    momentum.setX(pMuon * direction.x());
    momentum.setY(pMuon * direction.y());
    momentum.setZ(pMuon * direction.z());

    Hep3Vector position;
    // muon hits the top panel.
    if( arandom < pS_z/m_Smax) {
        position.setX(m_detectorCenter.x() +
                      (m_uni()-0.5) * m_detectorDim.x()); 
        position.setY(m_detectorCenter.y() +
                      (m_uni()-0.5) * m_detectorDim.y()); 
        position.setZ(m_detectorCenter.z() + 0.5 * m_detectorDim.z());

    } 
    // muon hits the X panel.
    else if( arandom < (pS_z + pS_x)/m_Smax) {
        position.setY( m_detectorCenter.y() +
                       (m_uni()-0.5) * m_detectorDim.y() ); 
        position.setZ( m_detectorCenter.z() +
                       (m_uni()-0.5) * m_detectorDim.z() ); 
        if (direction.x() > 0)
            position.setX( m_detectorCenter.x() - 0.5 * m_detectorDim.x() );
        if (direction.x() < 0)
            position.setX( m_detectorCenter.x() + 0.5 * m_detectorDim.x() );
    } 
    // muon hits the Y panel.
    else if( arandom < weight) {
        position.setX( m_detectorCenter.x() +
                       (m_uni()-0.5) * m_detectorDim.x() ); 
        position.setZ( m_detectorCenter.z() +
                       (m_uni()-0.5) * m_detectorDim.z() ); 
        if (direction.y() > 0)
            position.setY( m_detectorCenter.y() - 0.5 * m_detectorDim.y() );
        if (direction.y() < 0)
            position.setY( m_detectorCenter.y() + 0.5 * m_detectorDim.y() );
    }


    int pid = GetPIDFromMomentum(momentum.mag()/CLHEP::GeV);

    HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(position.x(), 
                                                                      position.y(), 
                                                                      position.z(), 0));

    HepMC::FourVector fourmom(momentum.x(), momentum.y(), momentum.z(), eMuon);

    debug() << "momentum.mag()=" << momentum.mag()/CLHEP::GeV << " GeV/c "
            << "pMuon=" << pMuon/CLHEP::GeV << " GeV/c "
            << "eMuon=" << eMuon/CLHEP::GeV << " GeV "
            << "muon_mass=" << muon_mass/CLHEP::MeV << " MeV/c^2 "
            << "sqrt(e^2-p^2)=" << sqrt(eMuon*eMuon-pMuon*pMuon)
            << endreq;

    HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,pid,1/*=status*/);
    particle->set_generated_mass(muon_mass);

    vertex->add_particle_out(particle);
    event.set_signal_process_vertex(vertex);

    return StatusCode::SUCCESS;
}

---