21 std::vector<PhotonMCTruth>
result;
30 const int ELECTRON_FLAV = 1;
31 const int PIZERO_FLAV = 2;
32 const int PHOTON_FLAV = 3;
37 std::vector<SimTrack*> photonTracks;
39 std::vector<SimTrack> trkFromConversion;
40 std::vector<ElectronMCTruth> electronsFromConversions;
43 fill(theSimTracks, theSimVertices);
46 if (!theSimTracks.empty()) {
50 std::vector<SimTrack>::const_iterator iFirstSimTk = theSimTracks.begin();
51 if (!(*iFirstSimTk).noVertex()) {
52 iPV = (*iFirstSimTk).vertIndex();
54 int vtxId = (*iFirstSimTk).vertIndex();
55 primVtx = theSimVertices[vtxId];
57 partType1 = (*iFirstSimTk).type();
68 if (iFirstSimTk != theSimTracks.end()) {
69 if ((*iFirstSimTk).vertIndex() == iPV) {
70 partType2 = (*iFirstSimTk).type();
84 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
85 if ((*iSimTk).noVertex())
92 if ((*iSimTk).vertIndex() == iPV) {
94 if ((*iSimTk).type() == 22) {
97 photonTracks.push_back(&(const_cast<SimTrack&>(*iSimTk)));
106 }
else if (npv == 1) {
109 ievflav = ELECTRON_FLAV;
110 }
else if (partType1 == 111) {
112 ievflav = PIZERO_FLAV;
113 }
else if (partType1 == 22) {
115 ievflav = PHOTON_FLAV;
117 }
else if (npv == 2) {
120 ievflav = ELECTRON_FLAV;
121 }
else if (partType1 == 111 && partType2 == 111) {
123 ievflav = PIZERO_FLAV;
124 }
else if (partType1 == 22 && partType2 == 22) {
126 ievflav = PHOTON_FLAV;
132 int isAconversion = 0;
133 int phoMotherType = -1;
134 int phoMotherVtxIndex = -1;
135 int phoMotherId = -1;
136 if (ievflav == PHOTON_FLAV || ievflav == PIZERO_FLAV || ievtype == PYTHIA) {
143 for (std::vector<SimTrack>::const_iterator iPhoTk = theSimTracks.begin(); iPhoTk != theSimTracks.end();
145 trkFromConversion.clear();
146 electronsFromConversions.clear();
148 if ((*iPhoTk).type() != 22)
150 int photonVertexIndex = (*iPhoTk).vertIndex();
151 int phoTrkId = (*iPhoTk).trackId();
157 if (
vertex.parentIndex() != -1) {
158 unsigned motherGeantId =
vertex.parentIndex();
162 phoMotherType = phoMotherId == -1 ? 0 : theSimTracks[phoMotherId].type();
164 if (phoMotherType == 111 || phoMotherType == 221 || phoMotherType == 331) {
170 for (std::vector<SimTrack>::const_iterator iEleTk = theSimTracks.begin(); iEleTk != theSimTracks.end();
172 if ((*iEleTk).noVertex())
174 if ((*iEleTk).vertIndex() == iPV)
176 if (
std::abs((*iEleTk).type()) != 11)
179 int vertexId = (*iEleTk).vertIndex();
184 if (
vertex.parentIndex() != -1) {
185 unsigned motherGeantId =
vertex.parentIndex();
194 std::vector<CLHEP::Hep3Vector> bremPos;
195 std::vector<CLHEP::HepLorentzVector> pBrem;
196 std::vector<float> xBrem;
198 if (theSimTracks[
motherId].trackId() == (*iPhoTk).trackId()) {
202 trkFromConversion.push_back((*iEleTk));
204 float remainingEnergy = trLast.
momentum().e();
208 (*iEleTk).momentum().y(),
209 (*iEleTk).momentum().z(),
210 (*iEleTk).momentum().e());
212 unsigned int eleId = (*iEleTk).trackId();
213 int eleVtxIndex = (*iEleTk).vertIndex();
219 for (std::vector<SimTrack>::const_iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end();
221 if ((*iSimTk).noVertex())
223 if ((*iSimTk).vertIndex() == iPV)
228 int vertexId1 = (*iSimTk).vertIndex();
229 const SimVertex& vertex1 = theSimVertices[vertexId1];
235 if ((vertexId1 == vertexId2) && ((*iSimTk).type() == (*iEleTk).type()) && trLast.
type() == 22) {
244 float eLoss = remainingEnergy - ((*iSimTk).momentum() + trLast.
momentum()).
e();
255 if (theSimTracks[
motherId].trackId() == eleId) {
257 eleId = (*iSimTk).trackId();
258 remainingEnergy = (*iSimTk).momentum().e();
259 motherMomentum = (*iSimTk).momentum();
261 pBrem.push_back(CLHEP::HepLorentzVector(trLast.
momentum().px(),
265 bremPos.push_back(CLHEP::HepLorentzVector(vertex1.
position().x(),
269 xBrem.push_back(
eLoss);
282 CLHEP::HepLorentzVector tmpEleMom(primEleMom.px(), primEleMom.py(), primEleMom.pz(), primEleMom.e());
283 CLHEP::HepLorentzVector tmpVtxPos(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
285 tmpEleMom, eleVtxIndex, bremPos, pBrem, xBrem, tmpVtxPos, const_cast<SimTrack&>(*iEleTk)));
297 CLHEP::HepLorentzVector phoMotherMom(0., 0., 0., 0.);
298 CLHEP::HepLorentzVector phoMotherVtx(0., 0., 0., 0.);
300 if (phoMotherId >= 0) {
301 phoMotherVtxIndex = theSimTracks[phoMotherId].vertIndex();
302 const SimVertex& motherVtx = theSimVertices[phoMotherVtxIndex];
306 phoMotherMom.setPx(theSimTracks[phoMotherId].momentum().
x());
307 phoMotherMom.setPy(theSimTracks[phoMotherId].momentum().
y());
308 phoMotherMom.setPz(theSimTracks[phoMotherId].momentum().
z());
309 phoMotherMom.setE(theSimTracks[phoMotherId].momentum().
t());
312 phoMotherVtx.setX(motherVtxPosition.x());
313 phoMotherVtx.setY(motherVtxPosition.y());
314 phoMotherVtx.setZ(motherVtxPosition.z());
315 phoMotherVtx.setT(motherVtxPosition.t());
318 if (!electronsFromConversions.empty()) {
324 int convVtxId = electronsFromConversions[0].vertexInd();
325 const SimVertex& convVtx = theSimVertices[convVtxId];
331 CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
332 (*iPhoTk).momentum().py(),
333 (*iPhoTk).momentum().pz(),
334 (*iPhoTk).momentum().e());
335 CLHEP::HepLorentzVector tmpVertex(vtxPosition.x(), vtxPosition.y(), vtxPosition.z(), vtxPosition.t());
336 CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
347 electronsFromConversions));
352 CLHEP::HepLorentzVector vtxPosition(0., 0., 0., 0.);
353 CLHEP::HepLorentzVector tmpPhoMom((*iPhoTk).momentum().px(),
354 (*iPhoTk).momentum().py(),
355 (*iPhoTk).momentum().pz(),
356 (*iPhoTk).momentum().e());
357 CLHEP::HepLorentzVector tmpPrimVtx(primVtxPos.x(), primVtxPos.y(), primVtxPos.z(), primVtxPos.t());
368 electronsFromConversions));
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
const math::XYZTLorentzVectorD & position() const
const math::XYZTLorentzVectorD & momentum() const
int type() const
particle type (HEP PDT convension)
void fill(const std::vector< SimTrack > &theSimTracks, const std::vector< SimVertex > &theSimVertices)
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
std::tuple< layerClusterToCaloParticle, caloParticleToLayerCluster > association
std::map< unsigned, unsigned > geantToIndex_
Abs< T >::type abs(const T &t)
if(threadIdxLocalY==0 &&threadIdxLocalX==0)