#include <MonopoleSteppingAction.h>

Inheritance diagram for MonopoleSteppingAction:

Public Member Functions
	MonopoleSteppingAction (edm::ParameterSet const &p)
void	update (const G4Step *)
	This routine will be called when the appropriate signal arrives.
void	update (const BeginOfTrack *)
	This routine will be called when the appropriate signal arrives.
void	update (const BeginOfRun *)
	This routine will be called when the appropriate signal arrives.
void	update (const BeginOfJob *)
	This routine will be called when the appropriate signal arrives.
	~MonopoleSteppingAction ()
Private Attributes
bool	actOnTrack
double	bZ
double	cInMByS
double	cMevToJ
double	cMeVToKgMByS
double	dirxStart
double	diryStart
double	dirzStart
double	eStart
double	magCharge
bool	mode
std::vector< int >	pdgCode
double	pxStart
double	pyStart
double	pzStart

Detailed Description

Definition at line 15 of file MonopoleSteppingAction.h.

Constructor & Destructor Documentation

MonopoleSteppingAction::MonopoleSteppingAction ( edm::ParameterSet const & p )

Definition at line 23 of file MonopoleSteppingAction.cc.

References edm::ParameterSet::getUntrackedParameter(), and mode.

                                                                        :
  actOnTrack (false), bZ(0) {
  mode = p.getUntrackedParameter<bool>("ChangeFromFirstStep",true);
  edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction set mode for"
                                   << " start at first step to " << mode;
}

MonopoleSteppingAction::~MonopoleSteppingAction ( )

Definition at line 30 of file MonopoleSteppingAction.cc.

{}

Member Function Documentation

void MonopoleSteppingAction::update ( const BeginOfJob * ) [virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfJob * >.

Definition at line 32 of file MonopoleSteppingAction.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::bField, bZ, edm::EventSetup::get(), MagneticField::inTesla(), AlCaHLTBitMon_ParallelJobs::p, x, detailsBasic3DVector::y, and z.

                                                         {

  const edm::EventSetup* iSetup = (*job)();
  edm::ESHandle<MagneticField> bFieldH;
  iSetup->get<IdealMagneticFieldRecord>().get(bFieldH);
  const MagneticField *bField = bFieldH.product();
  const GlobalPoint p(0,0,0);
  bZ = (bField->inTesla(p)).z();
  edm::LogInfo("SimG4CoreWatcher") << "Magnetic Field (X): " 
                                   << (bField->inTesla(p)).x() << " Y: "
                                   << (bField->inTesla(p)).y() << " Z: "
                                   << bZ;
}

void MonopoleSteppingAction::update ( const G4Step * ) [virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 98 of file MonopoleSteppingAction.cc.

References actOnTrack, bZ, cInMByS, cMevToJ, cMeVToKgMByS, dirxStart, diryStart, eStart, LogDebug, m, magCharge, mode, pxStart, pyStart, pzStart, edm::second(), and mathSSE::sqrt().

                                                       {

  if (actOnTrack) {
    double eT, pT, pZ, tStep;
    G4Track* aTrack = aStep->GetTrack();
    G4ThreeVector initialPosition(0,0,0);
    if (mode) {
      tStep = aTrack->GetGlobalTime();
      eT    = eStart*std::sqrt(dirxStart*dirxStart+diryStart*diryStart);
      pT    = std::sqrt(pxStart*pxStart+pyStart*pyStart);
      pZ    = pzStart;
    } else {
      G4ThreeVector dirStep = aTrack->GetMomentumDirection();
      double        lStep   = aTrack->GetStepLength();
      double        xStep   = aTrack->GetPosition().x()-lStep*dirStep.x();
      double        yStep   = aTrack->GetPosition().y()-lStep*dirStep.y();
      double        zStep   = aTrack->GetPosition().z()-lStep*dirStep.z();
      double        vStep   = aTrack->GetVelocity();
      initialPosition       = G4ThreeVector(xStep,yStep,zStep);
      tStep                 = (vStep > 0. ? (lStep/vStep) : 0.);
      double        eStep   = aTrack->GetTotalEnergy();
      eT                    = eStep*dirStep.perp();
      pT                    = aTrack->GetMomentum().perp();
      pZ                    = aTrack->GetMomentum().z();
    }
    LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: tStep " <<tStep
                                 << " eT " << eT << " pT " << pT << " pZ " 
                                 << pZ;
    double eT1 = eT*cMevToJ;
    double pT1 = pT*cMeVToKgMByS;
    double pZ1 = pZ*cMeVToKgMByS;
    double ts1 = tStep/CLHEP::second;
    double eT2 = (eT1 > 0. ? (1./eT1) : 0.);
    double fac0 = magCharge*bZ*cInMByS;
    double fac2 = pZ1*cInMByS*eT2;
    double fac1 = fac0*cInMByS*ts1*eT2 + fac2;
    double fact1 = (eT1/fac0)*(std::sqrt(1+fac1*fac1)-std::sqrt(1+fac2*fac2));
    double fact2 = (pT1/(magCharge*bZ)*(asinh(fac1)-asinh(fac2)));
    LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Factor eT = "
                                 << eT << " " << eT1 << " " << eT2
                                 << ", pT = " << pT << " " << pT1
                                 << ", pZ = " << pZ << " " << pZ1
                                 << ", time = " << tStep << " " << ts1
                                 << ", Factors " << fac0 << " " << fac1 
                                 << " " << fac2 << ", Final ones " 
                                 << fact1 << " " << fact2;
    G4ThreeVector ez(0,0,1), et(std::sqrt(0.5),std::sqrt(0.5),0);
    G4ThreeVector displacement = (fact1*ez + fact2*et)*CLHEP::m;
    
    LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction: Initial " 
                                 << initialPosition << " Displacement "
                                 << displacement << " Final "
                                 << (initialPosition+displacement);
    aTrack->SetPosition(initialPosition+displacement);
  }
}

void MonopoleSteppingAction::update ( const BeginOfTrack * ) [virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfTrack * >.

Definition at line 74 of file MonopoleSteppingAction.cc.

References actOnTrack, prof2calltree::count, dirxStart, diryStart, dirzStart, eStart, LogDebug, pdgCode, pxStart, pyStart, and pzStart.

                                                            {

  actOnTrack = false;
  if (pdgCode.size() > 0) {
    const G4Track * aTrack = (*trk)();
    int code = aTrack->GetDefinition()->GetPDGEncoding();
    if (std::count(pdgCode.begin(),pdgCode.end(),code) > 0) {
      actOnTrack = true;
      eStart     = aTrack->GetTotalEnergy();
      pxStart    = aTrack->GetMomentum().x();
      pyStart    = aTrack->GetMomentum().y();
      pzStart    = aTrack->GetMomentum().z();
      dirxStart  = aTrack->GetMomentumDirection().x();
      diryStart  = aTrack->GetMomentumDirection().y();
      dirzStart  = aTrack->GetMomentumDirection().z();
      LogDebug("SimG4CoreWatcher") << "MonopoleSeppingAction Track " 
                                   << code << " Flag " << actOnTrack
                                   << " (px,py,pz,E) = (" << pxStart/GeV
                                   << ", " << pyStart/GeV << ", "
                                   <<pzStart/GeV <<", " <<eStart/GeV <<")";
    }
  }
}

void MonopoleSteppingAction::update ( const BeginOfRun * ) [virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 46 of file MonopoleSteppingAction.cc.

References cInMByS, cMevToJ, cMeVToKgMByS, e_SI, magCharge, pdgCode, and edm::second().

                                                      {

  G4ParticleTable * partTable = G4ParticleTable::GetParticleTable();
  magCharge = CLHEP::e_SI/CLHEP::fine_structure_const * 0.5;
  for (int ii= 0; ii < partTable->size(); ii++) {
    G4ParticleDefinition * particle = partTable->GetParticle(ii);
    std::string particleName = (particle->GetParticleName()).substr(0,8);
    if (strcmp(particleName.c_str(),"Monopole") == 0) {
      magCharge = CLHEP::e_SI*((G4Monopole*)(particle))->MagneticCharge();
      pdgCode.push_back(particle->GetPDGEncoding());
    }
  }
  edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction Finds "
                                   << pdgCode.size() << " candidates";
  for (unsigned int ii=0; ii<pdgCode.size(); ++ii) {
    edm::LogInfo("SimG4CoreWatcher") << "PDG Code[" << ii << "] = "
                                     << pdgCode[ii];
  }
  cMevToJ   = CLHEP::e_SI/CLHEP::eV;
  cMeVToKgMByS = CLHEP::e_SI*CLHEP::meter/(CLHEP::eV*CLHEP::c_light*CLHEP::second);
  cInMByS   = CLHEP::c_light*CLHEP::second/CLHEP::meter;
  edm::LogInfo("SimG4CoreWatcher") << "MonopoleSeppingAction Constants:" 
                                   << " MevToJoules  = " << cMevToJ 
                                   << ", MevToKgm/s  = " << cMeVToKgMByS
                                   << ", c in m/s    = " << cInMByS 
                                   << ", mag. charge = " << magCharge;
}