Notification/src/TrackWithHistory.cc

0001 #include "SimG4Core/Notification/interface/TrackWithHistory.h"
0002 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0003 #include "SimG4Core/Notification/interface/TrackInformationExtractor.h"
0004 #include "SimG4Core/Notification/interface/GenParticleInfoExtractor.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "G4VProcess.hh"
0009
0010 #include <iostream>
0011
0012 G4ThreadLocal G4Allocator<TrackWithHistory>* fpTrackWithHistoryAllocator = nullptr;
0013
0014 //#define DEBUG
0015
0016 TrackWithHistory::TrackWithHistory(const G4Track* g4trk)
0017     : trackID_(0),
0018       particleID_(0),
0019       parentID_(0),
0020       momentum_(math::XYZVectorD(0., 0., 0.)),
0021       totalEnergy_(0),
0022       vertexPosition_(math::XYZVectorD(0., 0., 0.)),
0023       globalTime_(0),
0024       localTime_(0),
0025       properTime_(0),
0026       creatorProcess_(nullptr),
0027       weight_(0),
0028       storeTrack_(false),
0029       saved_(false) {
0030   if (g4trk != nullptr) {
0031     TrackInformationExtractor extractor;
0032     trackID_ = g4trk->GetTrackID();
0033     particleID_ = G4TrackToParticleID::particleID(g4trk);
0034     parentID_ = g4trk->GetParentID();
0035     momentum_ = math::XYZVectorD(g4trk->GetMomentum().x(), g4trk->GetMomentum().y(), g4trk->GetMomentum().z());
0036     totalEnergy_ = g4trk->GetTotalEnergy();
0037     vertexPosition_ = math::XYZVectorD(g4trk->GetPosition().x(), g4trk->GetPosition().y(), g4trk->GetPosition().z());
0038     globalTime_ = g4trk->GetGlobalTime();
0039     localTime_ = g4trk->GetLocalTime();
0040     properTime_ = g4trk->GetProperTime();
0041     creatorProcess_ = g4trk->GetCreatorProcess();
0042     storeTrack_ = extractor(g4trk).storeTrack();
0043     saved_ = false;
0044     crossedBoundary_ = false;
0045     genParticleID_ = extractGenID(g4trk);
0046     // V.I. weight is computed in the same way as before
0047     // without usage of G4Track::GetWeight()
0048     weight_ = 10000 * genParticleID_;
0049 #ifdef DEBUG
0050     LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother "
0051                                  << parentID_;
0052 #endif
0053   }
0054 }
0055
0056 void TrackWithHistory::checkAtEnd(const G4Track* gt) {
0057   math::XYZVectorD vposdir(gt->GetVertexPosition().x(), gt->GetVertexPosition().y(), gt->GetVertexPosition().z());
0058   math::XYZVectorD vmomdir(
0059       gt->GetVertexMomentumDirection().x(), gt->GetVertexMomentumDirection().y(), gt->GetVertexMomentumDirection().z());
0060   bool ok = true;
0061   double epsilon = 1.e-6;
0062   double eps2 = epsilon * epsilon;
0063   if ((vertexPosition_ - vposdir).Mag2() > eps2) {
0064     edm::LogWarning("TrackInformation") << "TrackWithHistory vertex position check failed"
0065                                         << "\nAt construction: " << vertexPosition_ << "\nAt end:          " << vposdir;
0066     ok = false;
0067   }
0068   math::XYZVectorD dirDiff = momentum_.Unit() - vmomdir;
0069   if (dirDiff.Mag2() > eps2 && momentum_.Unit().R() > eps2) {
0070     edm::LogWarning("TrackInformation") << "TrackWithHistory momentum direction check failed"
0071                                         << "\nAt construction: " << momentum_.Unit()
0072                                         << "\nAt end:          " << vmomdir;
0073     ok = false;
0074   }
0075
0076   if (!ok)
0077     G4Exception("TrackWithHistory::checkAtEnd()", "mc001", FatalException, "check at track end failed");
0078 }
0079
0080 int TrackWithHistory::extractGenID(const G4Track* gt) const {
0081   void* vgprimary = gt->GetDynamicParticle()->GetPrimaryParticle();
0082   if (vgprimary == nullptr)
0083     return -1;
0084   // replace old-style cast with appropriate new-style cast...
0085   G4PrimaryParticle* gprimary = (G4PrimaryParticle*)vgprimary;
0086   GenParticleInfoExtractor ext;
0087   return ext(gprimary).id();
0088 }