File indexing completed on 2024-04-06 12:24:57
0001 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
0002 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
0003 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
0004
0005
0006 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0007 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0008 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0009 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0010
0011 #include <algorithm>
0012
0013 PhotonMCTruthFinder::PhotonMCTruthFinder() {
0014
0015 }
0016
0017 std::vector<PhotonMCTruth> PhotonMCTruthFinder::find(const std::vector<SimTrack>& theSimTracks,
0018 const std::vector<SimVertex>& theSimVertices) {
0019
0020
0021 std::vector<PhotonMCTruth> result;
0022
0023
0024
0025
0026
0027 const int SINGLE = 1;
0028 const int DOUBLE = 2;
0029 const int PYTHIA = 3;
0030 const int ELECTRON_FLAV = 1;
0031 const int PIZERO_FLAV = 2;
0032 const int PHOTON_FLAV = 3;
0033
0034 int ievtype = 0;
0035 int ievflav = 0;
0036
0037 std::vector<SimTrack*> photonTracks;
0038
0039 std::vector<SimTrack> trkFromConversion;
0040 std::vector<ElectronMCTruth> electronsFromConversions;
0041 SimVertex primVtx;
0042
0043 fill(theSimTracks, theSimVertices);
0044
0045
0046 if (!theSimTracks.empty()) {
0047 int iPV = -1;
0048 int partType1 = 0;
0049 int partType2 = 0;
0050 std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
0051 if (!(*iFirstSimTk).noVertex()) {
0052 iPV = (*iFirstSimTk).vertIndex();
0053
0054 int vtxId = (*iFirstSimTk).vertIndex();
0055 primVtx = theSimVertices[vtxId];
0056
0057 partType1 = (*iFirstSimTk).type();
0058
0059 } else {
0060
0061 }
0062
0063 math::XYZTLorentzVectorD primVtxPos(
0064 primVtx.position().x(), primVtx.position().y(), primVtx.position().z(), primVtx.position().e());
0065
0066
0067 iFirstSimTk++;
0068 if (iFirstSimTk != theSimTracks.end()) {
0069 if ((*iFirstSimTk).vertIndex() == iPV) {
0070 partType2 = (*iFirstSimTk).type();
0071
0072
0073 } else {
0074
0075 }
0076 }
0077
0078
0079
0080 int npv = 0;
0081
0082
0083
0084 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0085 if ((*iSimTk).noVertex())
0086 continue;
0087
0088
0089
0090
0091
0092 if ((*iSimTk).vertIndex() == iPV) {
0093 npv++;
0094 if ((*iSimTk).type() == 22) {
0095
0096
0097 photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));
0098 }
0099 }
0100 }
0101
0102
0103
0104 if (npv >= 3) {
0105 ievtype = PYTHIA;
0106 } else if (npv == 1) {
0107 if (std::abs(partType1) == 11) {
0108 ievtype = SINGLE;
0109 ievflav = ELECTRON_FLAV;
0110 } else if (partType1 == 111) {
0111 ievtype = SINGLE;
0112 ievflav = PIZERO_FLAV;
0113 } else if (partType1 == 22) {
0114 ievtype = SINGLE;
0115 ievflav = PHOTON_FLAV;
0116 }
0117 } else if (npv == 2) {
0118 if (std::abs(partType1) == 11 && std::abs(partType2) == 11) {
0119 ievtype = DOUBLE;
0120 ievflav = ELECTRON_FLAV;
0121 } else if (partType1 == 111 && partType2 == 111) {
0122 ievtype = DOUBLE;
0123 ievflav = PIZERO_FLAV;
0124 } else if (partType1 == 22 && partType2 == 22) {
0125 ievtype = DOUBLE;
0126 ievflav = PHOTON_FLAV;
0127 }
0128 }
0129
0130
0131
0132 int isAconversion = 0;
0133 int phoMotherType = -1;
0134 int phoMotherVtxIndex = -1;
0135 int phoMotherId = -1;
0136 if (ievflav == PHOTON_FLAV || ievflav == PIZERO_FLAV || ievtype == PYTHIA) {
0137
0138
0139
0140
0141
0142
0143 for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end();
0144 ++iPhoTk) {
0145 trkFromConversion.clear();
0146 electronsFromConversions.clear();
0147
0148 if ((*iPhoTk).type() != 22)
0149 continue;
0150 int photonVertexIndex = (*iPhoTk).vertIndex();
0151 int phoTrkId = (*iPhoTk).trackId();
0152
0153
0154
0155 const SimVertex& vertex = theSimVertices[photonVertexIndex];
0156 phoMotherId = -1;
0157 if (vertex.parentIndex() != -1) {
0158 unsigned motherGeantId = vertex.parentIndex();
0159 std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0160 if (association != geantToIndex_.end())
0161 phoMotherId = association->second;
0162 phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
0163
0164 if (phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331) {
0165
0166
0167 }
0168 }
0169
0170 for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end();
0171 ++iEleTk) {
0172 if ((*iEleTk).noVertex())
0173 continue;
0174 if ((*iEleTk).vertIndex() == iPV)
0175 continue;
0176 if (std::abs((*iEleTk).type()) != 11)
0177 continue;
0178
0179 int vertexId = (*iEleTk).vertIndex();
0180 const SimVertex& vertex = theSimVertices[vertexId];
0181 int motherId = -1;
0182
0183
0184 if (vertex.parentIndex() != -1) {
0185 unsigned motherGeantId = vertex.parentIndex();
0186 std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0187 if (association != geantToIndex_.end())
0188 motherId = association->second;
0189
0190
0191
0192
0193
0194 std::vector<CLHEP::Hep3Vector> bremPos;
0195 std::vector<CLHEP::HepLorentzVector> pBrem;
0196 std::vector<float> xBrem;
0197
0198 if (theSimTracks[motherId].trackId() == (*iPhoTk).trackId()) {
0199
0200
0201
0202 trkFromConversion.push_back((*iEleTk));
0203 SimTrack trLast = (*iEleTk);
0204 float remainingEnergy = trLast.momentum().e();
0205
0206
0207 math::XYZTLorentzVectorD primEleMom((*iEleTk).momentum().x(),
0208 (*iEleTk).momentum().y(),
0209 (*iEleTk).momentum().z(),
0210 (*iEleTk).momentum().e());
0211 math::XYZTLorentzVectorD motherMomentum(primEleMom);
0212 unsigned int eleId = (*iEleTk).trackId();
0213 int eleVtxIndex = (*iEleTk).vertIndex();
0214
0215 bremPos.clear();
0216 pBrem.clear();
0217 xBrem.clear();
0218
0219 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end();
0220 ++iSimTk) {
0221 if ((*iSimTk).noVertex())
0222 continue;
0223 if ((*iSimTk).vertIndex() == iPV)
0224 continue;
0225
0226
0227
0228 int vertexId1 = (*iSimTk).vertIndex();
0229 const SimVertex& vertex1 = theSimVertices[vertexId1];
0230 int vertexId2 = trLast.vertIndex();
0231
0232
0233 int motherId = -1;
0234
0235 if ((vertexId1 == vertexId2) && ((*iSimTk).type() == (*iEleTk).type()) && trLast.type() == 22) {
0236
0237
0238
0239
0240
0241
0242
0243
0244 float eLoss = remainingEnergy - ((*iSimTk).momentum() + trLast.momentum()).e();
0245
0246
0247 if (vertex1.parentIndex() != -1) {
0248 unsigned motherGeantId = vertex1.parentIndex();
0249 std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0250 if (association != geantToIndex_.end())
0251 motherId = association->second;
0252
0253
0254
0255 if (theSimTracks[motherId].trackId() == eleId) {
0256
0257 eleId = (*iSimTk).trackId();
0258 remainingEnergy = (*iSimTk).momentum().e();
0259 motherMomentum = (*iSimTk).momentum();
0260
0261 pBrem.push_back(CLHEP::HepLorentzVector(trLast.momentum().px(),
0262 trLast.momentum().py(),
0263 trLast.momentum().pz(),
0264 trLast.momentum().e()));
0265 bremPos.push_back(CLHEP::HepLorentzVector(vertex1.position().x(),
0266 vertex1.position().y(),
0267 vertex1.position().z(),
0268 vertex1.position().t()));
0269 xBrem.push_back(eLoss);
0270 }
0271
0272 } else {
0273
0274 }
0275 }
0276 trLast = (*iSimTk);
0277
0278 }
0279
0280
0281
0282 CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(), primEleMom.py(), primEleMom.pz(), primEleMom.e());
0283 CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0284 electronsFromConversions.push_back(ElectronMCTruth(
0285 tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos, const_cast<SimTrack&>(*iEleTk)));
0286 }
0287
0288 } else {
0289
0290 }
0291
0292 }
0293
0294
0295
0296 math::XYZTLorentzVectorD motherVtxPosition(0., 0., 0., 0.);
0297 CLHEP::HepLorentzVector phoMotherMom(0., 0., 0., 0.);
0298 CLHEP::HepLorentzVector phoMotherVtx(0., 0., 0., 0.);
0299
0300 if (phoMotherId >= 0) {
0301 phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
0302 const SimVertex& motherVtx = theSimVertices[phoMotherVtxIndex];
0303 motherVtxPosition = math::XYZTLorentzVectorD(
0304 motherVtx.position().x(), motherVtx.position().y(), motherVtx.position().z(), motherVtx.position().e());
0305
0306 phoMotherMom.setPx(theSimTracks[phoMotherId].momentum().x());
0307 phoMotherMom.setPy(theSimTracks[phoMotherId].momentum().y());
0308 phoMotherMom.setPz(theSimTracks[phoMotherId].momentum().z());
0309 phoMotherMom.setE(theSimTracks[phoMotherId].momentum().t());
0310
0311
0312 phoMotherVtx.setX(motherVtxPosition.x());
0313 phoMotherVtx.setY(motherVtxPosition.y());
0314 phoMotherVtx.setZ(motherVtxPosition.z());
0315 phoMotherVtx.setT(motherVtxPosition.t());
0316 }
0317
0318 if (!electronsFromConversions.empty()) {
0319
0320 isAconversion = 1;
0321
0322
0323
0324 int convVtxId = electronsFromConversions[0].vertexInd();
0325 const SimVertex& convVtx = theSimVertices[convVtxId];
0326
0327 math::XYZTLorentzVectorD vtxPosition(
0328 convVtx.position().x(), convVtx.position().y(), convVtx.position().z(), convVtx.position().e());
0329
0330
0331 CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
0332 (*iPhoTk).momentum().py(),
0333 (*iPhoTk).momentum().pz(),
0334 (*iPhoTk).momentum().e());
0335 CLHEP::HepLorentzVector tmpVertex(vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t());
0336 CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0337
0338 result.push_back(PhotonMCTruth(isAconversion,
0339 tmpPhoMom,
0340 photonVertexIndex,
0341 phoTrkId,
0342 phoMotherType,
0343 phoMotherMom,
0344 phoMotherVtx,
0345 tmpVertex,
0346 tmpPrimVtx,
0347 electronsFromConversions));
0348
0349 } else {
0350 isAconversion = 0;
0351
0352 CLHEP::HepLorentzVector vtxPosition(0., 0., 0., 0.);
0353 CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
0354 (*iPhoTk).momentum().py(),
0355 (*iPhoTk).momentum().pz(),
0356 (*iPhoTk).momentum().e());
0357 CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
0358
0359 result.push_back(PhotonMCTruth(isAconversion,
0360 tmpPhoMom,
0361 photonVertexIndex,
0362 phoTrkId,
0363 phoMotherType,
0364 phoMotherMom,
0365 phoMotherVtx,
0366 vtxPosition,
0367 tmpPrimVtx,
0368 electronsFromConversions));
0369 }
0370
0371 }
0372
0373 }
0374
0375
0376 }
0377
0378 return result;
0379 }
0380
0381 void PhotonMCTruthFinder::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0382
0383
0384
0385
0386 unsigned nVtx = simVertices.size();
0387 unsigned nTks = simTracks.size();
0388
0389
0390
0391
0392 if (nVtx == 0)
0393 return;
0394
0395
0396
0397 for (unsigned it = 0; it < nTks; ++it) {
0398 geantToIndex_[simTracks[it].trackId()] = it;
0399
0400 }
0401 }