/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
GtMuoneratorTool Class Reference

Generate cosmic muons. More...

#include <GtMuoneratorTool.h>

Inheritance diagram for GtMuoneratorTool:
Inheritance graph
[legend]
Collaboration diagram for GtMuoneratorTool:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 GtMuoneratorTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~GtMuoneratorTool ()
virtual StatusCode initialize ()
virtual StatusCode finalize ()
virtual StatusCode mutate (HepMC::GenEvent &event)

Static Public Member Functions

static const InterfaceID & interfaceID ()

Private Member Functions

StatusCode ReadMuons ()
double GetMuPlusMinusRatio (double momInGeV)
int GetPIDFromMomentum (double momInGeV)

Private Attributes

std::string m_whichSite
std::string m_muonFileName
std::string m_ratioFileName
std::string m_volume
bool m_rotation
TH1F * m_pmratio
std::vector< double > m_muonE
std::vector< double > m_muonPhi
std::vector< double > m_muonTheta
double m_Sx
double m_Sy
double m_Sz
double m_Smax
CLHEP::Hep3Vector m_detectorDim
CLHEP::Hep3Vector m_detectorCenter
double m_dphi
IRndmGenSvc * m_rgs
Rndm::Numbers m_uni

Detailed Description

Generate cosmic muons.

Based on old Generators/Muon

bv@bnl.gov Wed May 5 16:39:04 2010

Definition at line 27 of file GtMuoneratorTool.h.


Constructor & Destructor Documentation

GtMuoneratorTool::GtMuoneratorTool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 29 of file GtMuoneratorTool.cc.

    : 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 ( ) [virtual]

Definition at line 53 of file GtMuoneratorTool.cc.

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

Member Function Documentation

StatusCode GtMuoneratorTool::initialize ( ) [virtual]

Definition at line 62 of file GtMuoneratorTool.cc.

{
    // 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();
}
StatusCode GtMuoneratorTool::finalize ( ) [virtual]

Definition at line 344 of file GtMuoneratorTool.cc.

{
    return StatusCode::SUCCESS;
}
StatusCode GtMuoneratorTool::mutate ( HepMC::GenEvent event) [virtual]

Implements IHepMCEventMutator.

Definition at line 349 of file GtMuoneratorTool.cc.

{
    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;
}
StatusCode GtMuoneratorTool::ReadMuons ( ) [private]

Definition at line 131 of file GtMuoneratorTool.cc.

{
    // 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()
double GtMuoneratorTool::GetMuPlusMinusRatio ( double  momInGeV) [private]

Definition at line 105 of file GtMuoneratorTool.cc.

{
    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) [private]

Definition at line 116 of file GtMuoneratorTool.cc.

{
    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+
}

Member Data Documentation

std::string GtMuoneratorTool::m_whichSite [private]

Definition at line 46 of file GtMuoneratorTool.h.

std::string GtMuoneratorTool::m_muonFileName [private]

Definition at line 47 of file GtMuoneratorTool.h.

std::string GtMuoneratorTool::m_ratioFileName [private]

Definition at line 48 of file GtMuoneratorTool.h.

std::string GtMuoneratorTool::m_volume [private]

Definition at line 49 of file GtMuoneratorTool.h.

Definition at line 50 of file GtMuoneratorTool.h.

TH1F* GtMuoneratorTool::m_pmratio [private]

Definition at line 53 of file GtMuoneratorTool.h.

std::vector<double> GtMuoneratorTool::m_muonE [private]

Definition at line 54 of file GtMuoneratorTool.h.

std::vector<double> GtMuoneratorTool::m_muonPhi [private]

Definition at line 54 of file GtMuoneratorTool.h.

std::vector<double> GtMuoneratorTool::m_muonTheta [private]

Definition at line 54 of file GtMuoneratorTool.h.

double GtMuoneratorTool::m_Sx [private]

Definition at line 55 of file GtMuoneratorTool.h.

double GtMuoneratorTool::m_Sy [private]

Definition at line 55 of file GtMuoneratorTool.h.

double GtMuoneratorTool::m_Sz [private]

Definition at line 55 of file GtMuoneratorTool.h.

double GtMuoneratorTool::m_Smax [private]

Definition at line 55 of file GtMuoneratorTool.h.

CLHEP::Hep3Vector GtMuoneratorTool::m_detectorDim [private]

Definition at line 56 of file GtMuoneratorTool.h.

CLHEP::Hep3Vector GtMuoneratorTool::m_detectorCenter [private]

Definition at line 57 of file GtMuoneratorTool.h.

double GtMuoneratorTool::m_dphi [private]

Definition at line 58 of file GtMuoneratorTool.h.

IRndmGenSvc* GtMuoneratorTool::m_rgs [private]

Definition at line 60 of file GtMuoneratorTool.h.

Rndm::Numbers GtMuoneratorTool::m_uni [private]

Definition at line 61 of file GtMuoneratorTool.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:20:55 for GenMuon by doxygen 1.7.4