File indexing completed on 2023-03-17 11:17:44
0001 #include "RecoEgamma/EgammaMCTools/interface/PizeroMCTruthFinder.h"
0002 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
0003 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"
0004
0005 #include <algorithm>
0006
0007 PizeroMCTruthFinder::PizeroMCTruthFinder() {
0008 thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
0009 theElectronMCTruthFinder_ = new ElectronMCTruthFinder();
0010 }
0011
0012 PizeroMCTruthFinder::~PizeroMCTruthFinder() {
0013 delete thePhotonMCTruthFinder_;
0014 delete theElectronMCTruthFinder_;
0015 std::cout << "~PizeroMCTruthFinder" << std::endl;
0016 }
0017
0018 std::vector<PizeroMCTruth> PizeroMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
0019 const std::vector<SimVertex>& theSimVertices) {
0020 std::cout << " PizeroMCTruthFinder::find " << std::endl;
0021
0022 std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
0023
0024 std::vector<PizeroMCTruth> result;
0025 std::vector<PhotonMCTruth> photonsFromPizero;
0026 std::vector<ElectronMCTruth> electronsFromPizero;
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039 std::vector<SimTrack> pizeroTracks;
0040 SimVertex primVtx;
0041
0042 fill(theSimTracks, theSimVertices);
0043
0044 int iPV = -1;
0045
0046
0047 std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
0048 if (!(*iFirstSimTk).noVertex()) {
0049 iPV = (*iFirstSimTk).vertIndex();
0050
0051 int vtxId = (*iFirstSimTk).vertIndex();
0052 primVtx = theSimVertices[vtxId];
0053
0054
0055 std::cout << " Primary vertex id " << iPV << " first track type " << (*iFirstSimTk).type() << std::endl;
0056 } else {
0057 std::cout << " First track has no vertex " << std::endl;
0058 }
0059
0060
0061 math::XYZTLorentzVectorD primVtxPos(
0062 primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
0063
0064
0065 iFirstSimTk++;
0066 if (iFirstSimTk != theSimTracks.end()) {
0067 if ((*iFirstSimTk).vertIndex() == iPV) {
0068
0069 std::cout << " second track type " << (*iFirstSimTk).type() << " vertex " << (*iFirstSimTk).vertIndex()
0070 << std::endl;
0071
0072 } else {
0073 std::cout << " Only one kine track from Primary Vertex " << std::endl;
0074 }
0075 }
0076
0077 int npv = 0;
0078
0079 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0080 if ((*iSimTk).noVertex())
0081 continue;
0082
0083 int vertexId = (*iSimTk).vertIndex();
0084 const SimVertex& vertex = theSimVertices[vertexId];
0085
0086 std::cout << " Particle type " << (*iSimTk).type() << " Sim Track ID " << (*iSimTk).trackId() << " momentum "
0087 << (*iSimTk).momentum() << " vertex position " << vertex.position() << std::endl;
0088
0089 if ((*iSimTk).vertIndex() == iPV) {
0090 npv++;
0091 if (std::abs((*iSimTk).type()) == 111) {
0092 std::cout << " Found a primary pizero with ID " << (*iSimTk).trackId() << " momentum " << (*iSimTk).momentum()
0093 << std::endl;
0094
0095 pizeroTracks.push_back(*iSimTk);
0096
0097
0098 math::XYZTLorentzVectorD momentum(
0099 (*iSimTk).momentum().x(), (*iSimTk).momentum().y(), (*iSimTk).momentum().z(), (*iSimTk).momentum().e());
0100 }
0101 }
0102 }
0103
0104 std::cout << " There are " << npv << " particles originating in the PV " << std::endl;
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126 for (std::vector<SimTrack>::iterator iPizTk = pizeroTracks.begin(); iPizTk != pizeroTracks.end(); ++iPizTk) {
0127 std::cout << " Looping on the primary pizero pt " << sqrt((*iPizTk).momentum().perp2()) << " pizero track ID "
0128 << (*iPizTk).trackId() << std::endl;
0129
0130 photonsFromPizero.clear();
0131 std::cout << " mcPhotons.size " << mcPhotons.size() << std::endl;
0132 for (std::vector<PhotonMCTruth>::iterator iPho = mcPhotons.begin(); iPho != mcPhotons.end(); ++iPho) {
0133 int phoVtxIndex = (*iPho).vertexInd();
0134 const SimVertex& phoVtx = theSimVertices[phoVtxIndex];
0135 unsigned int phoParentInd = phoVtx.parentIndex();
0136 std::cout << " photon parent vertex index " << phoParentInd << std::endl;
0137
0138 if (phoParentInd == (*iPizTk).trackId()) {
0139 std::cout << "Matched Photon ID " << (*iPho).trackId() << " vtx " << phoParentInd << " with pizero "
0140 << (*iPizTk).trackId() << std::endl;
0141 photonsFromPizero.push_back(*iPho);
0142 }
0143 }
0144 std::cout << " Photon matching the pizero vertex " << photonsFromPizero.size() << std::endl;
0145
0146
0147 CLHEP::HepLorentzVector tmpMom(
0148 (*iPizTk).momentum().px(), (*iPizTk).momentum().py(), (*iPizTk).momentum().pz(), (*iPizTk).momentum().e());
0149 CLHEP::HepLorentzVector tmpPos(
0150 primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().t());
0151 result.push_back(PizeroMCTruth(tmpMom, photonsFromPizero, tmpPos));
0152
0153 }
0154
0155 std::cout << " Pizero size " << result.size() << std::endl;
0156
0157 return result;
0158 }
0159
0160 void PizeroMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0161 unsigned nVtx = simVertices.size();
0162 unsigned nTks = simTracks.size();
0163
0164
0165 if (nVtx == 0)
0166 return;
0167
0168
0169
0170 for (unsigned it = 0; it < nTks; ++it) {
0171 geantToIndex_[simTracks[it].trackId()] = it;
0172 std::cout << " PizeroMCTruthFinder::fill it " << it << " simTracks[it].trackId() " << simTracks[it].trackId()
0173 << std::endl;
0174 }
0175 }