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

In This Package:

MpSpallation.cc
Go to the documentation of this file.
00001 
00010 #include "MuonProphet.h"
00011 
00012 
00015 int MuonProphet::spallationProducts(HepMC::GenEvent& event, int idx)
00016 {
00017   if(m_muon == 0) {
00018     error()<<"No muon track found"<<endreq;
00019     // Since no muon track, it should not reach this point.
00020     return -1;
00021   }
00022   
00025   HepMC::ThreeVector bkgPoint;
00026   bool bkgGenerated = false;
00027   bkgGenerated = generateOneBkg(idx, 
00028                                 bkgPoint);
00029 
00031   if( !bkgGenerated ) return 0;
00032 
00034   HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;
00035   if(m_genTools[idx]->mutate(*bkgEvt).isFailure()) {
00036     error()<<"Failed to run background gentool "<<m_genToolNames[idx]<<endreq;
00037     return -1;
00038   }
00039   
00040   // debug
00041 #if 0
00042   {
00043     HepMC::GenEvent::particle_iterator p;
00044     HepMC::GenEvent::particle_iterator particleEnd = bkgEvt->particles_end();
00045     for ( p = bkgEvt->particles_begin(); p!=particleEnd; ++p) {
00046       debug()<<*(*p)<<endreq;
00047     }
00048     HepMC::GenEvent::vertex_iterator v;
00049     HepMC::GenEvent::vertex_iterator vertexEnd = bkgEvt->vertices_end();
00050     for ( v = bkgEvt->vertices_begin(); v!=vertexEnd; ++v) {
00051       debug()<<*(*v)<<endreq;
00052     }
00053   }
00054 #endif
00055 
00057   // change the vertex position first
00058   setPosition(*bkgEvt,bkgPoint);
00059 
00060   // change the time of the vertex
00061   // Get muon vertex time 
00062   double t0 = m_muon->production_vertex()->position().t();
00063   // to do: A delta should be added to t0 to account for flight time of muon
00064   //        and flight time of muon's daughter particle which reached the 
00065   //        production vertex of the background.
00066   
00067   setTime(*bkgEvt,t0,m_genLifetimes[idx]);
00068 
00070   // drag vertices out from background event 
00071   HepMC::GenVertex *bkgVtx;
00072   while( !bkgEvt->vertices_empty() ){
00073 
00074     bkgVtx = *(bkgEvt->vertices_begin());
00075 
00076     // add to muon event
00077     // This will automatically transfer the ownership from the old event to new event.
00078     // So the original event can be safely deleted now.
00079     event.add_vertex(bkgVtx);
00080   }
00081 
00082   // delete the bkgEvt
00083   delete bkgEvt;
00084 
00085   return 0;
00086 }
00087 
00090 bool MuonProphet::generateOneBkg(int idx,  /* used to find the property of one background */
00091                                  HepMC::ThreeVector& bkgPoint)
00092 {
00095   if( !m_crossWater ) {
00096     return false;
00097   }
00098 
00100   double eMuon = m_muon->momentum().e();
00101 
00102   // Mineral oil and LS densities are assumed to be the same.
00103   double density = m_mineralOil->density();
00104   double yield = m_genYields[idx] ;
00105   double alpha = 0.77;
00106 
00108   double yieldAtThisEnergy = yield * pow( (eMuon / m_genYieldMeasuredAt[idx]), alpha);
00109 
00111   double ProbThres = yieldAtThisEnergy * density * m_trkLengthInWater;
00112 
00113   debug()<<m_genToolNames[idx]<<" eMuon= "<<eMuon/CLHEP::GeV<<"GeV "
00114          <<" trkLengthInWater= "<<m_trkLengthInWater/CLHEP::meter<<"m "
00115          <<"ProbThres= "<<ProbThres<<endreq;
00116 
00117   if( ProbThres > 1 || ProbThres < 0 ) {
00118     debug()<<"Yield too high or too low. Invalid probability "<<ProbThres<<endreq;
00119   }
00120 
00121   if( m_rnd.Rndm() > ProbThres) {
00123     return false;
00124   }
00125 
00129   
00130   // The position of background in the cylinder coordinate system of muon track
00131   // Muon direction is z.
00132   double r,phi,z;
00133 
00135   z = m_rnd.Rndm() * (m_tickWaterMax - m_tickWaterMin) + m_tickWaterMin;
00136   r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
00137   phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;
00138 
00139 
00140   // try point in coordinate system of intersect 1st
00141   HepGeom::Point3D<double> tryPoint;
00142 
00143   tryPoint.setX( r*cos(phi) );
00144   tryPoint.setY( r*sin(phi) );
00145   tryPoint.setZ( z );
00146 
00148   /* This is a muon track shown below.
00149    *
00150    *            intersect 1st (tick1)  intersect 2nd (tick2)
00151    * Primary vertex -----1--------C----------2------>
00152    *                        z     | theta
00153    *                              |  
00154    *                              | r
00155    *                              |
00156    *                              |
00157    *                           tryPoint
00158    *
00159    */
00160 
00161   // direction of muon
00162   HepGeom::Point3D<double> muonDir(m_unit.X(),
00163                                    m_unit.Y(),
00164                                    m_unit.Z());
00165  
00166   // need a y axis which is perpendicular to muonSlope
00167   HepGeom::Point3D<double> y( m_unit.X(),
00168                               m_unit.Y(),
00169                               0);
00170   y.rotateZ(3.14159265358979323846/2);
00171   if(y.r()<=0) y.setX(1);
00172 
00173   // tryPoint in ows coordinate system
00174   HepGeom::Point3D<double> tryPointGlobal = muonDir * tryPoint.r();
00175 
00176   // rotation
00177   tryPointGlobal.rotate(tryPoint.theta(), y);
00178   tryPointGlobal.rotate(tryPoint.phi(),   muonDir);
00179 
00180   // transform                                                                                                              
00181   tryPointGlobal.setX(tryPointGlobal.x() + m_point.X());
00182   tryPointGlobal.setY(tryPointGlobal.y() + m_point.Y());
00183   tryPointGlobal.setZ(tryPointGlobal.z() + m_point.Z());
00184 
00185   debug()<<tryPoint<<endreq;
00186   debug()<<tryPointGlobal<<endreq;
00187 
00190   Gaudi::XYZPoint p(tryPointGlobal.x(),
00191                     tryPointGlobal.y(),
00192                     tryPointGlobal.z());
00193   
00194   bool inside = false;
00195   for (unsigned int idx = 0; idx<m_ADs.size(); ++idx) 
00196     {
00197       if(m_ADs[idx]->isInside(p)) inside = true;
00198     }
00199   
00200   if(!inside) return false;
00201   
00202   // done
00203   bkgPoint.setX(tryPointGlobal.x());
00204   bkgPoint.setY(tryPointGlobal.y());
00205   bkgPoint.setZ(tryPointGlobal.z());
00206 
00207   return true;
00208 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:21:29 for MuonProphet by doxygen 1.7.4