CMSSW Reference Manual

00001 #include <iostream>
00002 #include <vector>
00003 #include <memory>
00004 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackPairFinder.h"
00005 // Framework
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Utilities/interface/Exception.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 //
00010 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00011 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00012 
00013 
00014 #include "DataFormats/EgammaTrackReco/interface/TrackCaloClusterAssociation.h"
00015 //
00016 #include "CLHEP/Units/PhysicalConstants.h"
00017 #include <TMath.h>
00018 #include <vector>
00019 #include <map>
00020 
00021 //using namespace std;
00022 
00023 ConversionTrackPairFinder::ConversionTrackPairFinder( ){ 
00024 
00025   LogDebug("ConversionTrackPairFinder") << " CTOR  " <<  "\n";  
00026 
00027 }
00028 
00029 ConversionTrackPairFinder::~ConversionTrackPairFinder() {
00030 
00031   LogDebug("ConversionTrackPairFinder") << " DTOR " <<  "\n";  
00032     
00033 }
00034 
00035 
00036 
00037 
00038 std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr>  ConversionTrackPairFinder::run(std::vector<reco::TransientTrack> outInTrk,  
00039                                                                                                         const edm::Handle<reco::TrackCollection>& outInTrkHandle,
00040                                                                                                         const edm::Handle<reco::TrackCaloClusterPtrAssociation>& outInTrackSCAssH, 
00041                                                                                                         std::vector<reco::TransientTrack> inOutTrk, 
00042                                                                                                         const edm::Handle<reco::TrackCollection>& inOutTrkHandle,
00043                                                                                                         const edm::Handle<reco::TrackCaloClusterPtrAssociation>& inOutTrackSCAssH  ) {
00044 
00045   
00046   LogDebug("ConversionTrackPairFinder")  << "ConversionTrackPairFinder::run " <<  "\n";  
00047   
00048   std::vector<reco::TransientTrack>  selectedOutInTk;
00049   std::vector<reco::TransientTrack>  selectedInOutTk;
00050   std::vector<reco::TransientTrack>  allSelectedTk;
00051   std::map<reco::TransientTrack,  reco::CaloClusterPtr>  scTrkAssocMap; 
00052  
00053   bool oneLeg=false;
00054   bool noTrack=false;  
00055   
00056   
00057   
00058   int iTrk=0;
00059   for( std::vector<reco::TransientTrack>::const_iterator  iTk =  outInTrk.begin(); iTk !=  outInTrk.end(); iTk++) {
00060     //std::cout  << " Out In Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum()  << "\n";  
00061     
00062     if ( iTk->numberOfValidHits() <3 ||   iTk->normalizedChi2() > 5000 ) continue; 
00063     
00064     
00065     
00066     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
00067     reco::TrackRef myTkRef= ttt->persistentTrackRef(); 
00068     //std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner momentum " <<  myTkRef->innerMomentum() << "\n";
00069     
00070     edm::Ref<reco::TrackCollection> trackRef(outInTrkHandle, iTrk );
00071     
00072     //std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner momentum " <<  trackRef->innerMomentum() << "\n";
00073 
00074     const reco::CaloClusterPtr  aClus = (*outInTrackSCAssH)[trackRef];
00075 
00076     //    std::cout << "ConversionTrackPairFinder  Reading the OutIn Map  " << *outInTrackSCAss[trackRef] <<  " " << &outInTrackSCAss[trackRef] <<  std::endl;
00077     //    std::cout << "ConversionTrackPairFinder  Out In track belonging to SC with energy " << aClus->energy() << "\n"; 
00078 
00079     scTrkAssocMap[*iTk]= aClus;
00080     selectedOutInTk.push_back(*iTk);
00081     allSelectedTk.push_back(*iTk);
00082 
00083     iTrk++;
00084     
00085   }
00086 
00087 
00088   iTrk=0;
00089   for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  inOutTrk.begin(); iTk !=  inOutTrk.end(); iTk++) {
00090     //   std::cout << " In Out Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
00091     
00092     if ( iTk->numberOfValidHits() <3 ||   iTk->normalizedChi2() >5000 ) continue; 
00093     
00094     
00095     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
00096     reco::TrackRef myTkRef= ttt->persistentTrackRef(); 
00097     //std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner momentum " <<  myTkRef->innerMomentum() << "\n";
00098     
00099     edm::Ref<reco::TrackCollection> trackRef(inOutTrkHandle, iTrk );
00100     
00101     //std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner momentum " <<  trackRef->innerMomentum() << "\n";
00102     
00103     const reco::CaloClusterPtr  aClus = (*inOutTrackSCAssH)[trackRef];
00104 
00105     //    std::cout << "ConversionTrackPairFinder  Filling the InOut Map  " << &(*inOutTrackSCAss[trackRef]) << " " << &inOutTrackSCAss[trackRef] <<  std::endl;
00106     // std::cout << "ConversionTrackPairFinder  In Out  track belonging to SC with energy " << aClus.energy() << "\n"; 
00107     
00108     scTrkAssocMap[*iTk]= aClus;
00109     selectedInOutTk.push_back(*iTk);
00110     allSelectedTk.push_back(*iTk);    
00111 
00112     iTrk++;
00113     
00114   }
00115   
00116 
00117   //  std::cout << " ConversionTrackPairFinder allSelectedTk size " << allSelectedTk.size() << "  scTrkAssocMap  size " <<  scTrkAssocMap.size() << "\n"; 
00118   
00119   // Sort tracks in decreasing number of hits
00120   if(selectedOutInTk.size() > 0)
00121     std::stable_sort(selectedOutInTk.begin(), selectedOutInTk.end(), ByNumOfHits());
00122   if(selectedInOutTk.size() > 0)
00123     std::stable_sort(selectedInOutTk.begin(), selectedInOutTk.end(), ByNumOfHits());
00124   if(allSelectedTk.size() > 0)
00125     std::stable_sort(allSelectedTk.begin(),   allSelectedTk.end(),   ByNumOfHits());
00126   
00127   
00128   
00129   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedOutInTk.begin(); iTk !=  selectedOutInTk.end(); iTk++) {
00130     // std::cout << " Selected Out In  Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
00131   //}
00132   //for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedInOutTk.begin(); iTk !=  selectedInOutTk.end(); iTk++) {
00133     //std::cout << " Selected In Out Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
00134   //}
00135   
00136   
00137   
00138   std::vector<reco::TransientTrack > thePair(2);
00139   std::vector<std::vector<reco::TransientTrack> > allPairs;
00140   std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr > allPairSCAss;
00141 
00142   std::map<reco::TransientTrack,  reco::CaloClusterPtr>::const_iterator iMap1;
00143   std::map<reco::TransientTrack,  reco::CaloClusterPtr>::const_iterator iMap2;
00144 
00145   
00146   
00147   if ( scTrkAssocMap.size() > 2 ){
00148     
00149 
00150     for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00151       for( iMap2 =  iMap1; iMap2 !=   scTrkAssocMap.end(); ++iMap2) {
00152         // consider only tracks associated to the same SC 
00153 
00154         if (  (iMap1->second) != (iMap2->second) ) continue;  
00155         
00156         if (   ((iMap1->first)).charge() *  ((iMap2->first)).charge()  < 0 ) {
00157           
00158           //      std::cout << " ConversionTrackPairFinde All selected from the map First  Track charge " <<   (iMap1->first).charge() << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << " inner momentum " <<  ((iMap1->first)).track().innerMomentum() << " Ass SC " << (iMap1->second)->energy() <<  "\n";  
00159           
00160           //std::cout << " ConversionTrackPairFinde All selected from the map Second  Track charge " <<   ((iMap2->first)).charge() << " Num of RecHits " <<  ((iMap2->first)).recHitsSize() << " inner momentum " <<  ((iMap2->first)).track().innerMomentum() << " Ass SC " << (iMap2->second)->energy()  <<  "\n";  
00161 
00162 
00163 
00164           thePair.clear();
00165           thePair.push_back( iMap1->first );
00166           thePair.push_back( iMap2->first  );
00167           allPairs.push_back ( thePair );       
00168           allPairSCAss[thePair]= iMap1->second; 
00169 
00170         }
00171       }
00172     }
00173 
00174     //    std::cout << " ConversionTrackPairFinder  INTERMIDIATE allPairSCAss size " << allPairSCAss.size() << "\n";
00175 
00176     if ( allPairSCAss.size() == 0) { 
00177       //      std::cout << " All Tracks had the same charge: Need to send out a single track  " <<   "\n";
00178 
00179       for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00180 
00181         thePair.clear();
00182         thePair.push_back(iMap1->first);
00183         allPairs.push_back ( thePair );
00184         allPairSCAss[thePair]= iMap1->second; 
00185         
00186       }
00187 
00188     }
00189 
00190 
00191 
00192 
00193 
00194   } else if (  (scTrkAssocMap.size() ==2) ) {
00195 
00196     iMap1=scTrkAssocMap.begin();
00197     iMap2=scTrkAssocMap.end();
00198     if (  (iMap1->second) == (iMap2->second)  &&
00199           ((iMap1->first).charge() * (iMap2->first).charge() < 0 )   ) {
00200 
00201 
00202       //      std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap1->first).charge() << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << " inner momentum " <<  ((iMap1->first)).track().innerMomentum() << " Ass SC " << (iMap1->second)->energy() <<  "\n";  
00203 
00204       //      std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap2->first).charge() << " Num of RecHits " <<  ((iMap2->first)).recHitsSize() << " inner momentum " <<  ((iMap2->first)).track().innerMomentum() << " Ass SC " << (iMap2->second)->energy() <<  "\n";  
00205       
00206 
00207       thePair.clear();
00208       thePair.push_back( iMap1->first );
00209       thePair.push_back( iMap2->first );
00210       allPairs.push_back ( thePair );   
00211 
00212       allPairSCAss[thePair]= iMap1->second; 
00213       
00214     } else {
00215       oneLeg=true; 
00216 
00217     }
00218     
00219   } else if  (scTrkAssocMap.size() ==1   ) { 
00220     oneLeg=true;  
00221   } else {
00222     noTrack=true;
00223   } 
00224   
00225   
00226   if ( oneLeg ) {
00227     thePair.clear();               
00228 
00229 
00230     iMap1=scTrkAssocMap.begin();   
00231     thePair.push_back(iMap1->first);
00232  
00233     allPairs.push_back ( thePair );
00234     allPairSCAss[thePair]= iMap1->second; 
00235 
00236     //    std::cout << "  WARNING ConversionTrackPairFinder::tracks The candidate has just one leg. Need to find another way to evaltuate the vertex !!! "   << "\n";
00237   }
00238   
00239   if ( noTrack) {
00240     //  std::cout << "  WARNING ConversionTrackPairFinder::tracks case noTrack "   << "\n";  
00241     thePair.clear();  
00242     allPairSCAss.clear();
00243   }
00244   
00245 
00246 
00247 
00248   for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00249         
00250     int nFound=0;
00251     for (  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
00252 
00253       if (  (iMap1->second) == (iPair->second)  )   nFound++;   
00254       
00255     }      
00256     
00257     if ( nFound == 0) {
00258 
00259       thePair.clear();  
00260       thePair.push_back(iMap1->first);
00261       
00262       allPairs.push_back ( thePair );
00263       allPairSCAss[thePair]= iMap1->second; 
00264     }
00265 
00266 
00267   }
00268 
00269 
00270   //  std::cout << " ConversionTrackPairFinder FINAL allPairSCAss size " << allPairSCAss.size() << "\n";
00271  
00272  
00273   // for (  std::map<std::vector<reco::TransientTrack>, const reco::CaloCluster*>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
00274   //std::cout << " ConversionTrackPairFindder FINAL allPairSCAss " << (iPair->first).size() << " SC Energy " << (iPair->second)->energy() << " eta " << (iPair->second)->eta() << " phi " <<  (iPair->second)->phi() << "\n";  
00275   //}
00276 
00277 
00278 
00279 
00280   return allPairSCAss;
00281  
00282 
00283 }
00284 
00285
ConversionTrackPairFinder.cc
