25 std::cout <<
" PizeroMCTruthFinder::find " << std::endl;
30 std::vector<PizeroMCTruth>
result;
31 std::vector<PhotonMCTruth> photonsFromPizero;
32 std::vector<ElectronMCTruth> electronsFromPizero;
38 const int ELECTRON_FLAV=1;
39 const int PIZERO_FLAV=2;
40 const int PHOTON_FLAV=3;
45 std::vector<SimTrack> pizeroTracks;
49 fill(theSimTracks, theSimVertices);
54 std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
55 if ( !(*iFirstSimTk).noVertex() ) {
56 iPV = (*iFirstSimTk).vertIndex();
58 int vtxId = (*iFirstSimTk).vertIndex();
59 primVtx = theSimVertices[vtxId];
61 partType1 = (*iFirstSimTk).type();
62 std::cout <<
" Primary vertex id " << iPV <<
" first track type " << (*iFirstSimTk).type() << std::endl;
64 std::cout <<
" First track has no vertex " << std::endl;
75 if ( iFirstSimTk!= theSimTracks.end() ) {
77 if ( (*iFirstSimTk).vertIndex() == iPV) {
78 partType2 = (*iFirstSimTk).type();
79 std::cout <<
" second track type " << (*iFirstSimTk).type() <<
" vertex " << (*iFirstSimTk).vertIndex() << std::endl;
82 std::cout <<
" Only one kine track from Primary Vertex " << std::endl;
88 for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
89 if ( (*iSimTk).noVertex() )
continue;
91 int vertexId = (*iSimTk).vertIndex();
92 SimVertex vertex = theSimVertices[vertexId];
94 std::cout <<
" Particle type " << (*iSimTk).type() <<
" Sim Track ID " << (*iSimTk).trackId() <<
" momentum " << (*iSimTk).momentum() <<
" vertex position " << vertex.
position() << std::endl;
96 if ( (*iSimTk).vertIndex() == iPV ) {
98 if (
std::abs((*iSimTk).type() ) == 111) {
100 std::cout <<
" Found a primary pizero with ID " << (*iSimTk).trackId() <<
" momentum " << (*iSimTk).momentum() << std::endl;
102 pizeroTracks.push_back( *iSimTk );
106 (*iSimTk).momentum().y(),
107 (*iSimTk).momentum().z(),
108 (*iSimTk).momentum().e());
116 std::cout <<
" There are " << npv <<
" particles originating in the PV " << std::endl;
120 }
else if(npv == 1) {
122 ievtype = SINGLE; ievflav = ELECTRON_FLAV;
123 }
else if(partType1 == 111) {
124 ievtype = SINGLE; ievflav = PIZERO_FLAV;
125 }
else if(partType1 == 22) {
126 ievtype = SINGLE; ievflav = PHOTON_FLAV;
128 }
else if(npv == 2) {
130 ievtype = DOUBLE; ievflav = ELECTRON_FLAV;
131 }
else if(partType1 == 111 && partType2 == 111) {
132 ievtype = DOUBLE; ievflav = PIZERO_FLAV;
133 }
else if(partType1 == 22 && partType2 == 22) {
134 ievtype = DOUBLE; ievflav = PHOTON_FLAV;
140 for (std::vector<SimTrack>::iterator iPizTk = pizeroTracks.begin(); iPizTk != pizeroTracks.end(); ++iPizTk){
141 std::cout <<
" Looping on the primary pizero pt " <<
sqrt((*iPizTk).momentum().perp2()) <<
" pizero track ID " << (*iPizTk).trackId() << std::endl;
143 photonsFromPizero.clear();
144 std::cout <<
" mcPhotons.size " << mcPhotons.size() << std::endl;
145 for ( std::vector<PhotonMCTruth>::iterator iPho=mcPhotons.begin(); iPho !=mcPhotons.end(); ++iPho ){
146 int phoVtxIndex = (*iPho).vertexInd();
147 SimVertex phoVtx = theSimVertices[phoVtxIndex];
149 std::cout <<
" photon parent vertex index " << phoParentInd << std::endl;
151 if (phoParentInd == (*iPizTk).trackId()) {
152 std::cout <<
"Matched Photon ID " << (*iPho).trackId() <<
" vtx " << phoParentInd <<
" with pizero " << (*iPizTk).trackId() << std::endl;
153 photonsFromPizero.push_back( *iPho);
157 std::cout <<
" Photon matching the pizero vertex " << photonsFromPizero.size() <<std::endl;
161 CLHEP::HepLorentzVector tmpMom( (*iPizTk).momentum().px(), (*iPizTk).momentum().py(),
162 (*iPizTk).momentum().pz(), (*iPizTk).momentum().e() ) ;
163 CLHEP::HepLorentzVector tmpPos( primVtx.
position().x(), primVtx.
position().y(),
165 result.push_back(
PizeroMCTruth ( tmpMom, photonsFromPizero, tmpPos ) );
171 std::cout <<
" Pizero size " << result.size() << std::endl;
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
PhotonMCTruthFinder * thePhotonMCTruthFinder_
void fill(std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)
const math::XYZTLorentzVectorD & position() const
std::vector< PhotonMCTruth > find(std::vector< SimTrack > simTracks, std::vector< SimVertex > simVertices)