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

In This Package:

MpNeutron.cc
Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 #include "MuonProphet.h"
00009 #include "../functions/NeutronPDF.h"
00010 #include "TRandom3.h"
00011 
00012 StatusCode MuonProphet::spallationNeutronsInit()
00013 {
00014   m_nEnergyPDF=0;
00015   m_nAnglePDF=0;
00016   m_nLateralPDF=0;
00017 
00022 
00023   // Avarage muon energy is from the estimation of Mengyun Guan. (Ph.D. thesis)
00024   double avaMuonE;
00025   if(m_site == Site::kDayaBay) {
00026     avaMuonE = 55.3*CLHEP::GeV;
00027   } else if (m_site == Site::kLingAo) {
00028     avaMuonE = 61.4*CLHEP::GeV;
00029   } else if (m_site == Site::kFar ) {
00030     avaMuonE = 140.3*CLHEP::GeV;
00031   }
00032   
00037   m_nEnergyPDF = new TF1("m_nEnergyPDF",NeutronEnergyPDF,0.001,4,1);
00038   m_nEnergyPDF->SetNpx(200);  
00039 
00040   m_nEnergyPDF->SetParameter(0,avaMuonE/CLHEP::GeV);
00041   if(!m_nEnergyPDF) {
00042     error()<<"Failed to get NeutronEnergyPDF"<<endreq;
00043     return StatusCode::FAILURE;
00044   }
00045 
00049   m_nAnglePDF = new TF1("m_nAnglePDF",NeutronAngularDis,-1,1,1);
00050   m_nAnglePDF->SetParameter(0,avaMuonE/CLHEP::GeV);
00051   if(!m_nAnglePDF) {
00052     error()<<"Failed to get NeutronAngularDis"<<endreq;
00053     return StatusCode::FAILURE;
00054   }
00055 
00056   
00058   m_nLateralPDF = new TF1("m_nLateralPDF",LateralProfile,0,100,1);
00059   m_nLateralPDF->SetParameter(0,0);
00060   if(!m_nLateralPDF) {
00061     error()<<"Failed to get LateralProfile"<<endreq;
00062     return StatusCode::FAILURE;
00063   }
00064 
00065   return StatusCode::SUCCESS;
00066 }
00067 
00068 
00070 int MuonProphet::spallationNeutrons(HepMC::GenEvent& event) 
00071 {
00072   debug()<<"Spallation Neutron production"<<endreq;
00073 
00074   //--------------------------------------------------------
00075   unsigned int NSeg = m_crucialSegments.size();
00076   
00077   double yield;
00078   if(m_neutronYield>=0) {
00079     yield = m_neutronYield;
00080   } else {
00081     yield = NeutronYield( m_muon->momentum().e()/CLHEP::GeV )*(CLHEP::cm2/CLHEP::gram);
00082   }
00083 
00084   double tracklength;
00085   double multiplicity;
00086 
00087   int nNewTrack = 0;
00088 
00090   HepMC::GenEvent * bkgEvt=new HepMC::GenEvent;
00091 
00092   for (unsigned int idx=0; idx<NSeg; ++idx) 
00093     {
00094       double prob;
00095 
00101       
00102       tracklength = m_crucialSegments[idx].first.second - m_crucialSegments[idx].first.first;
00103       multiplicity = NeutronMulti( m_muon->momentum().e()/CLHEP::GeV,
00104                                    tracklength );
00105       
00106       prob = yield
00107         * tracklength
00108         * m_crucialSegments[idx].second->density()
00109         / multiplicity;
00110       
00111       debug()<<"Yield= "<<NeutronYield( m_muon->momentum().e()/CLHEP::GeV )<<"cm2/g "
00112              <<"TrackLength= "<< tracklength/CLHEP::m <<"m "
00113              <<"Density= "<<m_crucialSegments[idx].second->density()/(CLHEP::g/CLHEP::cm3)<<"g/cm3 "
00114              <<"Multiplicity= "<< multiplicity <<endreq;
00115       debug()<<"prob= "<<prob<<endreq;
00116 
00117       
00118 
00119       if( m_rnd.Rndm() < prob ) 
00120         {
00122           nNewTrack += generateNeutrons(bkgEvt,
00123                                         m_crucialSegments[idx].first.first,   // tick one
00124                                         m_crucialSegments[idx].first.second,  // tick two
00125                                         multiplicity);
00126         }
00127     }  
00128   
00130   // drag vertices out from background event
00131 
00132   HepMC::GenVertex *bkgVtx;
00133   while( !bkgEvt->vertices_empty() )
00134     {
00135     bkgVtx = *(bkgEvt->vertices_begin());
00136     
00137     // add to muon event
00138     // This will automatically transfer the ownership from the old event to new event.
00139     // So the original event can be safely deleted now.
00140     event.add_vertex(bkgVtx);
00141     }
00142  
00143   // delete the bkgEvt
00144   delete bkgEvt;
00145 
00146   return nNewTrack;
00147   
00148 }
00149 
00150 
00152 // tickBegin:  In fact it is a double, indicate the first intersection with the volume
00153 // tickEnd:    Indecate the last intersection with the volume. See m_point, m_vec and m_unit
00154 
00155 int MuonProphet::generateNeutrons(HepMC::GenEvent *bkgEvt,
00156                                   ISolid::Tick tickBegin,
00157                                   ISolid::Tick tickEnd,
00158                                   double multiplicity)
00159 {
00160   int nNeutron = m_rnd.Poisson(multiplicity);
00161   
00162   debug()<<"nNeutron= "<<nNeutron<<endreq;
00163   
00166 
00168   double z;
00169   double r;
00170   double phi;
00171   
00173   double ke;
00174   double eng;
00175   double momentum;
00176   double massN = 939.565346*CLHEP::MeV;
00177   double phiN;
00178   double costh;
00179   double sinth;
00180 
00181   for(int idx =0; idx<nNeutron; ++idx)  
00182     { 
00184       z = m_rnd.Rndm() * (tickEnd - tickBegin) + tickBegin;
00185       r = m_nLateralPDF->GetRandom() * CLHEP::cm ;
00186       phi = m_rnd.Rndm() * 2 * 3.14159265358979323846;
00187 
00188       ke = m_nEnergyPDF->GetRandom() * CLHEP::GeV ;      
00189       eng = ke + massN;
00190       momentum = sqrt (eng*eng - massN*massN);
00191       phiN = m_rnd.Rndm() * 2 * 3.14159265358979323846;
00193       costh = m_nAnglePDF->GetRandom();
00194       sinth = sqrt(1-costh*costh);
00195 
00197 
00198       debug()<<"z "<<z/CLHEP::m<<"m"<<endreq;
00199       debug()<<"r "<<r/CLHEP::cm<<"cm"<<endreq;
00200       debug()<<"ke "<<ke/CLHEP::GeV<<"GeV"<<endreq;
00201       debug()<<"costh "<<costh<<endreq;
00202 
00204       HepGeom::Point3D<double> genPoint;
00205       
00206       genPoint.setX( r*cos(phi) );
00207       genPoint.setY( r*sin(phi) );
00208       genPoint.setZ( z );
00209 
00211       HepGeom::Point3D<double> genVec;
00212       
00213       genVec.setX ( momentum*sinth*cos(phiN) );
00214       genVec.setY ( momentum*sinth*sin(phiN) );
00215       genVec.setZ ( momentum*costh );
00216 
00218       HepGeom::Point3D<double> muonDir(m_unit.X(),
00219                                        m_unit.Y(),
00220                                        m_unit.Z());
00221 
00222       // need a y axis which is perpendicular to muonSlope
00223       HepGeom::Point3D<double> y( m_unit.X(),
00224                                   m_unit.Y(),
00225                                   0);
00226       y.rotateZ(3.14159265358979323846/2);
00227       if(y.r()<=0) y.setX(1);
00228       
00230       HepGeom::Point3D<double> genPointGlobal = muonDir * genPoint.r();
00231       debug()<<m_unit<<endreq;
00232       debug()<<muonDir<<endreq;
00233       debug()<<genPointGlobal<<endreq;
00234 
00235       // rotation
00236       genPointGlobal.rotate(genPoint.theta(), y);
00237       genPointGlobal.rotate(genPoint.phi(),   muonDir);
00238 
00239       // transform
00240       genPointGlobal.setX(genPointGlobal.x() + m_point.X());
00241       genPointGlobal.setY(genPointGlobal.y() + m_point.Y());
00242       genPointGlobal.setZ(genPointGlobal.z() + m_point.Z());
00243 
00245       HepGeom::Point3D<double> genVectorGlobal = muonDir * genVec.r();
00246       
00247       genVectorGlobal.rotate(genVec.theta(),  y);
00248       genVectorGlobal.rotate(genVec.phi(),    muonDir);
00249       
00250 
00252       double t0 = m_muon->production_vertex()->position().t();
00253       debug()<<"t0= "<<t0<<endreq;
00254       t0 += tickEnd / (3e8 * CLHEP::m/CLHEP::s);
00255       debug()<<"dt= "<<tickEnd / (3e8 * CLHEP::m/CLHEP::s)<<endreq;
00256       debug()<<"t0= "<<t0<<endreq;
00257       HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(genPointGlobal.x(),
00258                                                                         genPointGlobal.y(),
00259                                                                         genPointGlobal.z(),
00260                                                                         t0));
00261       
00262       HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(genVectorGlobal.x(),
00263                                                                               genVectorGlobal.y(),
00264                                                                               genVectorGlobal.z(),
00265                                                                               eng),
00266                                                                               2112,
00267                                                                               1/*=status*/);
00269       vertex->add_particle_out(particle);
00270       bkgEvt->add_vertex(vertex);
00271     } 
00272 
00273   return nNeutron;
00274 }
00275 
00276 
| 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