Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:57 

0001 #include <iostream>
0002 #include <vector>
0003 #include <memory>
0004 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackPairFinder.h"
0005 // Framework
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 //
0008 
0009 //
0010 #include <vector>
0011 #include <map>
0012 
0013 //using namespace std;
0014 
0015 ConversionTrackPairFinder::ConversionTrackPairFinder() {
0016   LogDebug("ConversionTrackPairFinder") << " CTOR  "
0017                                         << "\n";
0018 }
0019 
0020 ConversionTrackPairFinder::~ConversionTrackPairFinder() {
0021   LogDebug("ConversionTrackPairFinder") << " DTOR "
0022                                         << "\n";
0023 }
0024 
0025 std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>
0026 ConversionTrackPairFinder::run(const std::vector<reco::TransientTrack>& outInTrk,
0027                                const edm::Handle<reco::TrackCollection>& outInTrkHandle,
0028                                const edm::Handle<reco::TrackCaloClusterPtrAssociation>& outInTrackSCAssH,
0029                                const std::vector<reco::TransientTrack>& _inOutTrk,
0030                                const edm::Handle<reco::TrackCollection>& inOutTrkHandle,
0031                                const edm::Handle<reco::TrackCaloClusterPtrAssociation>& inOutTrackSCAssH) {
0032   std::vector<reco::TransientTrack> inOutTrk = _inOutTrk;
0033 
0034   LogDebug("ConversionTrackPairFinder") << "ConversionTrackPairFinder::run "
0035                                         << "\n";
0036 
0037   std::vector<reco::TransientTrack> selectedOutInTk;
0038   std::vector<reco::TransientTrack> selectedInOutTk;
0039   std::vector<reco::TransientTrack> allSelectedTk;
0040   std::map<reco::TransientTrack, reco::CaloClusterPtr, CompareTwoTracks> scTrkAssocMap;
0041   std::multimap<int, reco::TransientTrack, std::greater<int> > auxMap;
0042 
0043   bool oneLeg = false;
0044   bool noTrack = false;
0045 
0046   int iTrk = 0;
0047   for (std::vector<reco::TransientTrack>::const_iterator iTk = outInTrk.begin(); iTk != outInTrk.end(); iTk++) {
0048     edm::Ref<reco::TrackCollection> trackRef(outInTrkHandle, iTrk);
0049     iTrk++;
0050 
0051     if (iTk->numberOfValidHits() < 3 || iTk->normalizedChi2() > 5000)
0052       continue;
0053     if (fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
0054         fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
0055         fabs(iTk->impactPointState().globalPosition().z()) > 280)
0056       continue;
0057 
0058     //    std::cout  << " Out In Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0059     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0060     reco::TrackRef myTkRef = ttt->persistentTrackRef();
0061     //std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0062     //    std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0063 
0064     const reco::CaloClusterPtr aClus = (*outInTrackSCAssH)[trackRef];
0065 
0066     //    std::cout << "ConversionTrackPairFinder  Reading the OutIn Map  " << *outInTrackSCAss[trackRef] <<  " " << &outInTrackSCAss[trackRef] <<  std::endl;
0067     //    std::cout << "ConversionTrackPairFinder  Out In track belonging to SC with energy " << aClus->energy() << "\n";
0068 
0069     int nHits = iTk->recHitsSize();
0070     scTrkAssocMap[*iTk] = aClus;
0071     auxMap.insert(std::pair<int, reco::TransientTrack>(nHits, (*iTk)));
0072     selectedOutInTk.push_back(*iTk);
0073     allSelectedTk.push_back(*iTk);
0074   }
0075 
0076   iTrk = 0;
0077   for (std::vector<reco::TransientTrack>::const_iterator iTk = inOutTrk.begin(); iTk != inOutTrk.end(); iTk++) {
0078     edm::Ref<reco::TrackCollection> trackRef(inOutTrkHandle, iTrk);
0079     iTrk++;
0080 
0081     if (iTk->numberOfValidHits() < 3 || iTk->normalizedChi2() > 5000)
0082       continue;
0083     if (fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
0084         fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
0085         fabs(iTk->impactPointState().globalPosition().z()) > 280)
0086       continue;
0087 
0088     //    std::cout << " In Out Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0089     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0090     reco::TrackRef myTkRef = ttt->persistentTrackRef();
0091     // std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0092     //    std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0093 
0094     const reco::CaloClusterPtr aClus = (*inOutTrackSCAssH)[trackRef];
0095 
0096     //    std::cout << "ConversionTrackPairFinder  Filling the InOut Map  " << &(*inOutTrackSCAss[trackRef]) << " " << &inOutTrackSCAss[trackRef] <<  std::endl;
0097     // std::cout << "ConversionTrackPairFinder  In Out  track belonging to SC with energy " << aClus.energy() << "\n";
0098 
0099     scTrkAssocMap[*iTk] = aClus;
0100     int nHits = iTk->recHitsSize();
0101     auxMap.insert(std::pair<int, reco::TransientTrack>(nHits, (*iTk)));
0102     selectedInOutTk.push_back(*iTk);
0103     allSelectedTk.push_back(*iTk);
0104   }
0105 
0106   //  std::cout << " ConversionTrackPairFinder allSelectedTk size " << allSelectedTk.size() << "  scTrkAssocMap  size " <<  scTrkAssocMap.size() << "\n";
0107 
0108   // Sort tracks in decreasing number of hits
0109   if (!selectedOutInTk.empty())
0110     std::stable_sort(selectedOutInTk.begin(), selectedOutInTk.end(), ByNumOfHits());
0111   if (!selectedInOutTk.empty())
0112     std::stable_sort(selectedInOutTk.begin(), selectedInOutTk.end(), ByNumOfHits());
0113   if (!allSelectedTk.empty())
0114     std::stable_sort(allSelectedTk.begin(), allSelectedTk.end(), ByNumOfHits());
0115 
0116   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedOutInTk.begin(); iTk !=  selectedOutInTk.end(); iTk++) {
0117   // std::cout << " Selected Out In  Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0118   //}
0119   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedInOutTk.begin(); iTk !=  selectedInOutTk.end(); iTk++) {
0120   // std::cout << " Selected In Out Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0121   // }
0122   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  allSelectedTk.begin(); iTk !=  allSelectedTk.end(); iTk++) {
0123   // std::cout << " All Selected  Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " chi2 " << iTk->normalizedChi2() <<  " pt " <<  sqrt(iTk->track().innerMomentum().perp2()) << "\n";
0124   //}
0125 
0126   std::vector<reco::TransientTrack> thePair(2);
0127   std::vector<std::vector<reco::TransientTrack> > allPairs;
0128 
0129   std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairSCAss;
0130   std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairOrdInPtSCAss;
0131   std::map<reco::TransientTrack, reco::CaloClusterPtr>::const_iterator iMap1;
0132   std::map<reco::TransientTrack, reco::CaloClusterPtr>::const_iterator iMap2;
0133 
0134   for (iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
0135     //   std::cout << " Ass map track  charge " << (iMap1->first).charge()  <<" pt " << sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " SC E  " << (iMap1->second)->energy() << " SC eta " << (iMap1->second)->eta() <<  " SC phi " << (iMap1->second)->phi() << std::endl;
0136   }
0137 
0138   std::multimap<int, reco::TransientTrack>::const_iterator iAux;
0139 
0140   // for( iAux = auxMap.begin(); iAux!= auxMap.end(); ++iAux) {
0141   //  //   std::cout << " Aux Map  " << (iAux->first)  <<" pt " << sqrt(((iAux->second)).track().innerMomentum().Perp2()) << std::endl;
0142   // for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
0143   //   if ( (iMap1->first) == (iAux->second) ) std::cout << " ass SC " <<  (iMap1->second)->energy() << std::endl;
0144   // }
0145   // }
0146 
0147   if (scTrkAssocMap.size() > 2) {
0148     for (iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
0149       for (iMap2 = iMap1; iMap2 != scTrkAssocMap.end(); ++iMap2) {
0150         // consider only tracks associated to the same SC
0151 
0152         if ((iMap1->second) != (iMap2->second))
0153           continue;
0154 
0155         if (((iMap1->first)).charge() * ((iMap2->first)).charge() < 0) {
0156           //      std::cout << " ConversionTrackPairFinde All selected from the map First  Track charge " <<   (iMap1->first).charge() << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << " inner pt " <<  sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap1->second)->energy() <<  "\n";
0157 
0158           //  std::cout << " ConversionTrackPairFinde All selected from the map Second  Track charge " <<   ((iMap2->first)).charge() << " Num of RecHits " <<  ((iMap2->first)).recHitsSize() << " inner pt " <<  sqrt(((iMap2->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap2->second)->energy()  <<  "\n";
0159 
0160           thePair.clear();
0161           thePair.push_back(iMap1->first);
0162           thePair.push_back(iMap2->first);
0163           allPairs.push_back(thePair);
0164           allPairSCAss[thePair] = iMap1->second;
0165         }
0166       }
0167     }
0168 
0169     //    std::cout << " ConversionTrackPairFinder  INTERMIDIATE allPairSCAss size " << allPairSCAss.size() << "\n";
0170 
0171     if (allPairSCAss.empty()) {
0172       //      std::cout << " All Tracks had the same charge: Need to send out a single track  " <<   "\n";
0173 
0174       for (iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
0175         thePair.clear();
0176         thePair.push_back(iMap1->first);
0177         allPairs.push_back(thePair);
0178         allPairSCAss[thePair] = iMap1->second;
0179       }
0180     }
0181 
0182   } else if ((scTrkAssocMap.size() == 2)) {
0183     iMap1 = scTrkAssocMap.begin();  //get the first
0184     iMap2 = iMap1;
0185     iMap2++;  //get the second
0186     if ((iMap1->second) == (iMap2->second)) {
0187       if ((iMap1->first).charge() * (iMap2->first).charge() < 0) {
0188         //  std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap1->first).charge() << std::endl;
0189         //std::cout << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << std::endl;
0190         //  std::cout << " inner pt " <<  sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << std::endl;
0191         //std::cout << " Ass SC " << (iMap1->second)->energy() <<  "\n";
0192 
0193         //  std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap2->first).charge() << std::endl;
0194         //  std::cout << " Num of RecHits " <<  ((iMap2->first)).recHitsSize() << std::endl;
0195         //std::cout << " inner pt " <<  sqrt(((iMap2->first)).track().innerMomentum().Perp2()) << std::endl;
0196         //std::cout << " Ass SC " << (iMap2->second)->energy() <<  "\n";
0197 
0198         thePair.clear();
0199         thePair.push_back(iMap1->first);
0200         thePair.push_back(iMap2->first);
0201         allPairs.push_back(thePair);
0202 
0203         allPairSCAss[thePair] = iMap1->second;
0204 
0205       } else {
0206         //std::cout << " ConversionTrackPairFinder oneLeg case when 2 tracks with same sign Pick up the longest one" << std::endl;
0207         if (((iMap1->first)).recHitsSize() > ((iMap2->first)).recHitsSize()) {
0208           thePair.clear();
0209           thePair.push_back(iMap1->first);
0210           allPairs.push_back(thePair);
0211           allPairSCAss[thePair] = iMap1->second;
0212         } else {
0213           thePair.clear();
0214           thePair.push_back(iMap2->first);
0215           allPairs.push_back(thePair);
0216           allPairSCAss[thePair] = iMap2->second;
0217         }
0218       }
0219     }
0220 
0221   } else if (scTrkAssocMap.size() == 1) {  /// ONly one track in input to the finder
0222     //    std::cout << " ConversionTrackPairFinder oneLeg case when 1 track only " << std::endl;
0223     oneLeg = true;
0224   } else {
0225     noTrack = true;
0226   }
0227 
0228   if (oneLeg) {
0229     thePair.clear();
0230     // std::cout << " ConversionTrackPairFinder oneLeg case charge  " << std::endl;
0231 
0232     iMap1 = scTrkAssocMap.begin();
0233 
0234     //std::cout << " ConversionTrackPairFinder oneLeg case charge  " <<   (iMap1->first).charge() << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << " inner pt " <<  sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap1->second)->energy() <<  "\n";
0235 
0236     thePair.push_back(iMap1->first);
0237     allPairs.push_back(thePair);
0238     allPairSCAss[thePair] = iMap1->second;
0239 
0240     // std::cout << "  WARNING ConversionTrackPairFinder::tracks The candidate has just one leg. Need to find another way to evaltuate the vertex !!! "   << "\n";
0241   }
0242 
0243   if (noTrack) {
0244     //    std::cout << "  WARNING ConversionTrackPairFinder::tracks case noTrack "   << "\n";
0245     thePair.clear();
0246     allPairSCAss.clear();
0247   }
0248 
0249   /// all cases above failed and some track-SC association is still missing
0250   for (iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
0251     int nFound = 0;
0252     for (std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair = allPairSCAss.begin();
0253          iPair != allPairSCAss.end();
0254          ++iPair) {
0255       if ((iMap1->second) == (iPair->second))
0256         nFound++;
0257     }
0258 
0259     if (nFound == 0) {
0260       //      std::cout << " nFound zero case " << std::endl;
0261       int iList = 0;
0262       for (iAux = auxMap.begin(); iAux != auxMap.end(); ++iAux) {
0263         if ((iMap1->first) == (iAux->second) && iList == 0) {
0264           thePair.clear();
0265           thePair.push_back(iAux->second);
0266           allPairSCAss[thePair] = iMap1->second;
0267         }
0268 
0269         iList++;
0270       }
0271     }
0272   }
0273 
0274   // order the tracks in the pair in order of decreasing pt
0275   for (std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair = allPairSCAss.begin();
0276        iPair != allPairSCAss.end();
0277        ++iPair) {
0278     thePair.clear();
0279     if ((iPair->first).size() == 2) {
0280       if (sqrt((iPair->first)[0].track().innerMomentum().perp2()) >
0281           sqrt((iPair->first)[1].track().innerMomentum().perp2())) {
0282         thePair.push_back((iPair->first)[0]);
0283         thePair.push_back((iPair->first)[1]);
0284       } else {
0285         thePair.push_back((iPair->first)[1]);
0286         thePair.push_back((iPair->first)[0]);
0287       }
0288     } else {
0289       thePair.push_back((iPair->first)[0]);
0290     }
0291 
0292     allPairOrdInPtSCAss[thePair] = iPair->second;
0293   }
0294 
0295   //    std::cout << " ConversionTrackPairFinder FINAL allPairOrdInPtSCAss size " << allPairOrdInPtSCAss.size() << "\n";
0296   // for (  std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr>::const_iterator iPair= allPairOrdInPtSCAss.begin(); iPair!= allPairOrdInPtSCAss.end(); ++iPair ) {
0297   // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss " << (iPair->first).size() << " SC Energy " << (iPair->second)->energy() << " eta " << (iPair->second)->eta() << " phi " <<  (iPair->second)->phi() << "\n";
0298   // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss (iPair->first).size() " << (iPair->first).size() << std::endl;
0299 
0300   // for ( std::vector<reco::TransientTrack>::const_iterator iTk=(iPair->first).begin(); iTk!= (iPair->first).end(); ++iTk) {
0301   //    std::cout << " ConversionTrackPair ordered track pt " << sqrt(iTk->track().innerMomentum().perp2()) << std::endl;
0302   // }
0303   //}
0304 
0305   return allPairOrdInPtSCAss;
0306 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!