CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4Core/HelpfulWatchers/src/MonopoleSteppingAction.cc

Go to the documentation of this file.
00001 #include "SimG4Core/HelpfulWatchers/interface/MonopoleSteppingAction.h"
00002 
00003 #include "SimG4Core/Notification/interface/BeginOfJob.h"
00004 #include "SimG4Core/Notification/interface/BeginOfRun.h"
00005 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
00006 #include "SimG4Core/Physics/interface/G4Monopole.hh"
00007 
00008 #include "FWCore/Framework/interface/ESHandle.h"
00009 #include "FWCore/Framework/interface/EventSetup.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 
00012 #include "MagneticField/Engine/interface/MagneticField.h"
00013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00014 
00015 #include "G4Track.hh"
00016 #include "G4Run.hh"
00017 #include "G4Event.hh"
00018 #include "G4ParticleTable.hh"
00019 
00020 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00021 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00022 
00023 MonopoleSteppingAction::MonopoleSteppingAction(edm::ParameterSet const & p) :
00024   actOnTrack (false), bZ(0) {
00025   mode = p.getUntrackedParameter<bool>("ChangeFromFirstStep",true);
00026   edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction set mode for"
00027                                    << " start at first step to " << mode;
00028 }
00029 
00030 MonopoleSteppingAction::~MonopoleSteppingAction() {}
00031 
00032 void MonopoleSteppingAction::update(const BeginOfJob* job) {
00033 
00034   const edm::EventSetup* iSetup = (*job)();
00035   edm::ESHandle<MagneticField> bFieldH;
00036   iSetup->get<IdealMagneticFieldRecord>().get(bFieldH);
00037   const MagneticField *bField = bFieldH.product();
00038   const GlobalPoint p(0,0,0);
00039   bZ = (bField->inTesla(p)).z();
00040   edm::LogInfo("SimG4CoreWatcher") << "Magnetic Field (X): " 
00041                                    << (bField->inTesla(p)).x() << " Y: "
00042                                    << (bField->inTesla(p)).y() << " Z: "
00043                                    << bZ;
00044 }
00045 
00046 void MonopoleSteppingAction::update(const BeginOfRun* ) {
00047 
00048   G4ParticleTable * partTable = G4ParticleTable::GetParticleTable();
00049   magCharge = CLHEP::e_SI/CLHEP::fine_structure_const * 0.5;
00050   for (int ii= 0; ii < partTable->size(); ii++) {
00051     G4ParticleDefinition * particle = partTable->GetParticle(ii);
00052     std::string particleName = (particle->GetParticleName()).substr(0,8);
00053     if (strcmp(particleName.c_str(),"Monopole") == 0) {
00054       magCharge = CLHEP::e_SI*((G4Monopole*)(particle))->MagneticCharge();
00055       pdgCode.push_back(particle->GetPDGEncoding());
00056     }
00057   }
00058   edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction Finds "
00059                                    << pdgCode.size() << " candidates";
00060   for (unsigned int ii=0; ii<pdgCode.size(); ++ii) {
00061     edm::LogInfo("SimG4CoreWatcher") << "PDG Code[" << ii << "] = "
00062                                      << pdgCode[ii];
00063   }
00064   cMevToJ   = CLHEP::e_SI/CLHEP::eV;
00065   cMeVToKgMByS = CLHEP::e_SI*CLHEP::meter/(CLHEP::eV*CLHEP::c_light*CLHEP::second);
00066   cInMByS   = CLHEP::c_light*CLHEP::second/CLHEP::meter;
00067   edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction Constants:" 
00068                                    << " MevToJoules  = " << cMevToJ 
00069                                    << ", MevToKgm/s  = " << cMeVToKgMByS
00070                                    << ", c in m/s    = " << cInMByS 
00071                                    << ", mag. charge = " << magCharge;
00072 }
00073 
00074 void MonopoleSteppingAction::update(const BeginOfTrack * trk) {
00075 
00076   actOnTrack = false;
00077   if (pdgCode.size() > 0) {
00078     const G4Track * aTrack = (*trk)();
00079     int code = aTrack->GetDefinition()->GetPDGEncoding();
00080     if (std::count(pdgCode.begin(),pdgCode.end(),code) > 0) {
00081       actOnTrack = true;
00082       eStart     = aTrack->GetTotalEnergy();
00083       pxStart    = aTrack->GetMomentum().x();
00084       pyStart    = aTrack->GetMomentum().y();
00085       pzStart    = aTrack->GetMomentum().z();
00086       dirxStart  = aTrack->GetMomentumDirection().x();
00087       diryStart  = aTrack->GetMomentumDirection().y();
00088       dirzStart  = aTrack->GetMomentumDirection().z();
00089       LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction Track " 
00090                                    << code << " Flag " << actOnTrack
00091                                    << " (px,py,pz,E) = (" << pxStart/GeV
00092                                    << ", " << pyStart/GeV << ", "
00093                                    <<pzStart/GeV <<", " <<eStart/GeV <<")";
00094     }
00095   }
00096 }
00097 
00098 void MonopoleSteppingAction::update(const G4Step* aStep) {
00099 
00100   if (actOnTrack) {
00101     double eT, pT, pZ, tStep;
00102     G4Track* aTrack = aStep->GetTrack();
00103     G4ThreeVector initialPosition(0,0,0);
00104     if (mode) {
00105       tStep = aTrack->GetGlobalTime();
00106       eT    = eStart*std::sqrt(dirxStart*dirxStart+diryStart*diryStart);
00107       pT    = std::sqrt(pxStart*pxStart+pyStart*pyStart);
00108       pZ    = pzStart;
00109     } else {
00110       G4ThreeVector dirStep = aTrack->GetMomentumDirection();
00111       double        lStep   = aTrack->GetStepLength();
00112       double        xStep   = aTrack->GetPosition().x()-lStep*dirStep.x();
00113       double        yStep   = aTrack->GetPosition().y()-lStep*dirStep.y();
00114       double        zStep   = aTrack->GetPosition().z()-lStep*dirStep.z();
00115       double        vStep   = aTrack->GetVelocity();
00116       initialPosition       = G4ThreeVector(xStep,yStep,zStep);
00117       tStep                 = (vStep > 0. ? (lStep/vStep) : 0.);
00118       double        eStep   = aTrack->GetTotalEnergy();
00119       eT                    = eStep*dirStep.perp();
00120       pT                    = aTrack->GetMomentum().perp();
00121       pZ                    = aTrack->GetMomentum().z();
00122     }
00123     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: tStep " <<tStep
00124                                  << " eT " << eT << " pT " << pT << " pZ " 
00125                                  << pZ;
00126     double eT1 = eT*cMevToJ;
00127     double pT1 = pT*cMeVToKgMByS;
00128     double pZ1 = pZ*cMeVToKgMByS;
00129     double ts1 = tStep/CLHEP::second;
00130     double eT2 = (eT1 > 0. ? (1./eT1) : 0.);
00131     double fac0 = magCharge*bZ*cInMByS;
00132     double fac2 = pZ1*cInMByS*eT2;
00133     double fac1 = fac0*cInMByS*ts1*eT2 + fac2;
00134     double fact1 = (eT1/fac0)*(std::sqrt(1+fac1*fac1)-std::sqrt(1+fac2*fac2));
00135     double fact2 = (pT1/(magCharge*bZ)*(asinh(fac1)-asinh(fac2)));
00136     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Factor eT = "
00137                                  << eT << " " << eT1 << " " << eT2
00138                                  << ", pT = " << pT << " " << pT1
00139                                  << ", pZ = " << pZ << " " << pZ1
00140                                  << ", time = " << tStep << " " << ts1
00141                                  << ", Factors " << fac0 << " " << fac1 
00142                                  << " " << fac2 << ", Final ones " 
00143                                  << fact1 << " " << fact2;
00144     G4ThreeVector ez(0,0,1), et(std::sqrt(0.5),std::sqrt(0.5),0);
00145     G4ThreeVector displacement = (fact1*ez + fact2*et)*CLHEP::m;
00146     
00147     LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Initial " 
00148                                  << initialPosition << " Displacement "
00149                                  << displacement << " Final "
00150                                  << (initialPosition+displacement);
00151     aTrack->SetPosition(initialPosition+displacement);
00152   }
00153 }