Reference Manual

Go to the documentation of this file.
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 #include "TrackingTools/TransientTrack/interface/GsfTransientTrack.h"
00013 
00014 #include "DataFormats/EgammaTrackReco/interface/TrackCaloClusterAssociation.h"
00015 //
00016 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00017 #include <TMath.h>
00018 #include <vector>
00019 #include <map>
00020 
00021 
00022 //using namespace std;
00023 
00024 ConversionTrackPairFinder::ConversionTrackPairFinder( ){ 
00025 
00026   LogDebug("ConversionTrackPairFinder") << " CTOR  " <<  "\n";  
00027 
00028 }
00029 
00030 ConversionTrackPairFinder::~ConversionTrackPairFinder() {
00031 
00032   LogDebug("ConversionTrackPairFinder") << " DTOR " <<  "\n";  
00033     
00034 }
00035 
00036 
00037 
00038 
00039 std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr>  ConversionTrackPairFinder::run(std::vector<reco::TransientTrack> outInTrk,  
00040                                                                                                         const edm::Handle<reco::TrackCollection>& outInTrkHandle,
00041                                                                                                         const edm::Handle<reco::TrackCaloClusterPtrAssociation>& outInTrackSCAssH, 
00042                                                                                                         std::vector<reco::TransientTrack> inOutTrk, 
00043                                                                                                         const edm::Handle<reco::TrackCollection>& inOutTrkHandle,
00044                                                                                                         const edm::Handle<reco::TrackCaloClusterPtrAssociation>& inOutTrackSCAssH  ) 
00045 {
00046 
00047   
00048   LogDebug("ConversionTrackPairFinder")  << "ConversionTrackPairFinder::run " <<  "\n";  
00049   
00050   std::vector<reco::TransientTrack>  selectedOutInTk;
00051   std::vector<reco::TransientTrack>  selectedInOutTk;
00052   std::vector<reco::TransientTrack>  allSelectedTk;
00053   std::map<reco::TransientTrack,  reco::CaloClusterPtr>  scTrkAssocMap; 
00054   std::multimap<int,reco::TransientTrack,std::greater<int> >  auxMap; 
00055  
00056   bool oneLeg=false;
00057   bool noTrack=false;  
00058   
00059   
00060   
00061   int iTrk=0;
00062   for( std::vector<reco::TransientTrack>::iterator  iTk =  outInTrk.begin(); iTk !=  outInTrk.end(); iTk++) {
00063     edm::Ref<reco::TrackCollection> trackRef(outInTrkHandle, iTrk );    
00064     iTrk++;
00065     
00066     if ( iTk->numberOfValidHits() <3 ||   iTk->normalizedChi2() > 5000 ) continue; 
00067     if ( fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
00068          fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
00069          fabs(iTk->impactPointState().globalPosition().z()) > 280 ) continue;
00070     
00071     //    std::cout  << " Out In Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00072     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
00073     reco::TrackRef myTkRef= ttt->persistentTrackRef(); 
00074     //std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00075     //    std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00076 
00077     const reco::CaloClusterPtr  aClus = (*outInTrackSCAssH)[trackRef];
00078 
00079     //    std::cout << "ConversionTrackPairFinder  Reading the OutIn Map  " << *outInTrackSCAss[trackRef] <<  " " << &outInTrackSCAss[trackRef] <<  std::endl;
00080     //    std::cout << "ConversionTrackPairFinder  Out In track belonging to SC with energy " << aClus->energy() << "\n"; 
00081 
00082     int nHits=iTk->recHitsSize();
00083     scTrkAssocMap[*iTk]= aClus;
00084     auxMap.insert(std::pair<int,reco::TransientTrack >(nHits,(*iTk)) );
00085     selectedOutInTk.push_back(*iTk);
00086     allSelectedTk.push_back(*iTk);
00087 
00088 
00089     
00090   }
00091 
00092 
00093   iTrk=0;
00094   for(  std::vector<reco::TransientTrack>::iterator  iTk =  inOutTrk.begin(); iTk !=  inOutTrk.end(); iTk++) {
00095     edm::Ref<reco::TrackCollection> trackRef(inOutTrkHandle, iTrk );
00096     iTrk++;
00097     
00098     if ( iTk->numberOfValidHits() <3 ||   iTk->normalizedChi2() >5000 ) continue; 
00099     if ( fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
00100          fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
00101          fabs(iTk->impactPointState().globalPosition().z()) > 280 ) continue;
00102     
00103     //    std::cout << " In Out Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";   
00104     const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
00105     reco::TrackRef myTkRef= ttt->persistentTrackRef(); 
00106     // std::cout <<  " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00107     //    std::cout <<  " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt  " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00108     
00109     const reco::CaloClusterPtr  aClus = (*inOutTrackSCAssH)[trackRef];
00110 
00111     //    std::cout << "ConversionTrackPairFinder  Filling the InOut Map  " << &(*inOutTrackSCAss[trackRef]) << " " << &inOutTrackSCAss[trackRef] <<  std::endl;
00112     // std::cout << "ConversionTrackPairFinder  In Out  track belonging to SC with energy " << aClus.energy() << "\n"; 
00113     
00114     scTrkAssocMap[*iTk]= aClus;
00115     int nHits=iTk->recHitsSize();
00116     auxMap.insert(std::pair<int,reco::TransientTrack >(nHits,(*iTk)) );
00117     selectedInOutTk.push_back(*iTk);
00118     allSelectedTk.push_back(*iTk);    
00119 
00120 
00121     
00122   }
00123   
00124 
00125   //  std::cout << " ConversionTrackPairFinder allSelectedTk size " << allSelectedTk.size() << "  scTrkAssocMap  size " <<  scTrkAssocMap.size() << "\n"; 
00126   
00127   // Sort tracks in decreasing number of hits
00128   if(selectedOutInTk.size() > 0)
00129     std::stable_sort(selectedOutInTk.begin(), selectedOutInTk.end(), ByNumOfHits());
00130   if(selectedInOutTk.size() > 0)
00131     std::stable_sort(selectedInOutTk.begin(), selectedInOutTk.end(), ByNumOfHits());
00132   if(allSelectedTk.size() > 0)
00133     std::stable_sort(allSelectedTk.begin(),   allSelectedTk.end(),   ByNumOfHits());
00134   
00135   
00136   
00137   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedOutInTk.begin(); iTk !=  selectedOutInTk.end(); iTk++) {
00138     // std::cout << " Selected Out In  Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
00139   //}
00140   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  selectedInOutTk.begin(); iTk !=  selectedInOutTk.end(); iTk++) {
00141   // std::cout << " Selected In Out Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";  
00142   // }
00143   //  for(  std::vector<reco::TransientTrack>::const_iterator  iTk =  allSelectedTk.begin(); iTk !=  allSelectedTk.end(); iTk++) {
00144   // std::cout << " All Selected  Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " chi2 " << iTk->normalizedChi2() <<  " pt " <<  sqrt(iTk->track().innerMomentum().perp2()) << "\n";  
00145   //}
00146   
00147 
00148   
00149   
00150   std::vector<reco::TransientTrack > thePair(2);
00151   std::vector<std::vector<reco::TransientTrack> > allPairs;
00152   std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr > allPairSCAss;
00153   std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr > allPairOrdInPtSCAss;
00154 
00155   std::map<reco::TransientTrack,  reco::CaloClusterPtr>::const_iterator iMap1;
00156   std::map<reco::TransientTrack,  reco::CaloClusterPtr>::const_iterator iMap2;
00157 
00158  for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00159  
00160    //   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;
00161   }
00162 
00163 
00164  std::multimap<int, reco::TransientTrack>::const_iterator iAux;
00165 
00166 
00167  // for( iAux = auxMap.begin(); iAux!= auxMap.end(); ++iAux) {
00168  //  //   std::cout << " Aux Map  " << (iAux->first)  <<" pt " << sqrt(((iAux->second)).track().innerMomentum().Perp2()) << std::endl;
00169  // for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00170  //   if ( (iMap1->first) == (iAux->second) ) std::cout << " ass SC " <<  (iMap1->second)->energy() << std::endl;
00171  // }
00172  // }
00173 
00174   if ( scTrkAssocMap.size() > 2 ){
00175     
00176 
00177     for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00178       for( iMap2 =  iMap1; iMap2 !=   scTrkAssocMap.end(); ++iMap2) {
00179         // consider only tracks associated to the same SC 
00180 
00181         if (  (iMap1->second) != (iMap2->second) ) continue;  
00182         
00183         if (   ((iMap1->first)).charge() *  ((iMap2->first)).charge()  < 0 ) {
00184           
00185           //      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";  
00186           
00187           //  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";  
00188 
00189 
00190 
00191           thePair.clear();
00192           thePair.push_back( iMap1->first );
00193           thePair.push_back( iMap2->first  );
00194           allPairs.push_back ( thePair );       
00195           allPairSCAss[thePair]= iMap1->second; 
00196 
00197         }
00198       }
00199     }
00200 
00201     //    std::cout << " ConversionTrackPairFinder  INTERMIDIATE allPairSCAss size " << allPairSCAss.size() << "\n";
00202 
00203     if ( allPairSCAss.size() == 0) { 
00204       //      std::cout << " All Tracks had the same charge: Need to send out a single track  " <<   "\n";
00205 
00206       for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00207 
00208         thePair.clear();
00209         thePair.push_back(iMap1->first);
00210         allPairs.push_back ( thePair );
00211         allPairSCAss[thePair]= iMap1->second; 
00212         
00213       }
00214 
00215     }
00216 
00217 
00218 
00219 
00220 
00221   } else if (  (scTrkAssocMap.size() ==2) ) {
00222     
00223     iMap1=scTrkAssocMap.begin();//get the first
00224     iMap2=iMap1;
00225     iMap2++;//get the second
00226     if (  (iMap1->second) == (iMap2->second)  ) {
00227       if  ( (iMap1->first).charge() * (iMap2->first).charge() < 0 )  {
00228         
00229         //      std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap1->first).charge() << std::endl;
00230         //std::cout << " Num of RecHits " <<  ((iMap1->first)).recHitsSize() << std::endl;
00231         //      std::cout << " inner pt " <<  sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << std::endl; 
00232         //std::cout << " Ass SC " << (iMap1->second)->energy() <<  "\n";  
00233 
00234         //      std::cout << " ConversionTrackPairFinder Case when  (scTrkAssocMap.size() ==2)  " <<   (iMap2->first).charge() << std::endl;
00235         //      std::cout << " Num of RecHits " <<  ((iMap2->first)).recHitsSize() << std::endl;
00236         //std::cout << " inner pt " <<  sqrt(((iMap2->first)).track().innerMomentum().Perp2()) << std::endl; 
00237         //std::cout << " Ass SC " << (iMap2->second)->energy() <<  "\n";  
00238         
00239         thePair.clear();
00240         thePair.push_back( iMap1->first );
00241         thePair.push_back( iMap2->first );
00242         allPairs.push_back ( thePair ); 
00243         
00244         allPairSCAss[thePair]= iMap1->second; 
00245 
00246 
00247       } else {
00248         //std::cout << " ConversionTrackPairFinder oneLeg case when 2 tracks with same sign Pick up the longest one" << std::endl;
00249         if (  ((iMap1->first)).recHitsSize() > ((iMap2->first)).recHitsSize() ) {
00250           thePair.clear();
00251           thePair.push_back(iMap1->first);
00252           allPairs.push_back ( thePair );
00253           allPairSCAss[thePair]= iMap1->second; 
00254         } else {
00255           thePair.clear();
00256           thePair.push_back(iMap2->first);
00257           allPairs.push_back ( thePair );
00258           allPairSCAss[thePair]= iMap2->second;
00259         }
00260       }
00261       
00262     }
00263 
00264   } else if  (scTrkAssocMap.size() ==1   ) { 
00265     //    std::cout << " ConversionTrackPairFinder oneLeg case when 1 track only " << std::endl;
00266     oneLeg=true;  
00267   } else {
00268     noTrack=true;
00269   } 
00270   
00271   
00272   if ( oneLeg ) {
00273     thePair.clear();               
00274     // std::cout << " ConversionTrackPairFinder oneLeg case charge  " << std::endl;
00275                                                       
00276                                                       
00277     iMap1=scTrkAssocMap.begin();   
00278                                                       
00279     //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";  
00280                                                       
00281     thePair.push_back(iMap1->first);
00282     allPairs.push_back ( thePair );
00283     allPairSCAss[thePair]= iMap1->second; 
00284 
00285     // std::cout << "  WARNING ConversionTrackPairFinder::tracks The candidate has just one leg. Need to find another way to evaltuate the vertex !!! "   << "\n";
00286   }
00287   
00288   if ( noTrack) {
00289     //    std::cout << "  WARNING ConversionTrackPairFinder::tracks case noTrack "   << "\n";  
00290     thePair.clear();  
00291     allPairSCAss.clear();
00292   }
00293   
00294 
00295 
00297   for( iMap1 =   scTrkAssocMap.begin(); iMap1 !=   scTrkAssocMap.end(); ++iMap1) {
00298     
00299     int nFound=0;
00300     for (  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
00301       if (  (iMap1->second) == (iPair->second)  ) nFound++;
00302     }
00303     
00304     if ( nFound == 0) {
00305       //      std::cout << " nFound zero case " << std::endl;
00306       int iList=0;
00307       for( iAux = auxMap.begin(); iAux!= auxMap.end(); ++iAux) {
00308         if ( (iMap1->first) == (iAux->second) && iList==0 )  {
00309           thePair.clear();   
00310           thePair.push_back(iAux->second);
00311           allPairSCAss[thePair]= iMap1->second;         
00312           
00313         }
00314         
00315         iList++;
00316       }  
00317     }
00318     
00319   }
00320 
00321 
00322 
00323 
00324   // order the tracks in the pair in order of decreasing pt 
00325     for (  std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
00326       thePair.clear();
00327       if ( (iPair->first).size() ==2 ) {
00328         if (  sqrt((iPair->first)[0].track().innerMomentum().perp2()) > sqrt((iPair->first)[1].track().innerMomentum().perp2())  ) {
00329           thePair.push_back((iPair->first)[0]);
00330           thePair.push_back((iPair->first)[1]);
00331         } else {
00332           thePair.push_back((iPair->first)[1]);
00333           thePair.push_back((iPair->first)[0]);
00334         }
00335       } else {
00336         thePair.push_back((iPair->first)[0]);
00337       }
00338 
00339       allPairOrdInPtSCAss[thePair]=iPair->second;
00340     }
00341 
00342 
00343     //    std::cout << " ConversionTrackPairFinder FINAL allPairOrdInPtSCAss size " << allPairOrdInPtSCAss.size() << "\n";
00344     // for (  std::map<std::vector<reco::TransientTrack>,  reco::CaloClusterPtr>::const_iterator iPair= allPairOrdInPtSCAss.begin(); iPair!= allPairOrdInPtSCAss.end(); ++iPair ) {
00345     // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss " << (iPair->first).size() << " SC Energy " << (iPair->second)->energy() << " eta " << (iPair->second)->eta() << " phi " <<  (iPair->second)->phi() << "\n";  
00346     // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss (iPair->first).size() " << (iPair->first).size() << std::endl;
00347       
00348     // for ( std::vector<reco::TransientTrack>::const_iterator iTk=(iPair->first).begin(); iTk!= (iPair->first).end(); ++iTk) {
00349     //  std::cout << " ConversionTrackPair ordered track pt " << sqrt(iTk->track().innerMomentum().perp2()) << std::endl;
00350     // }
00351     //}
00352 
00353   
00354 
00355   return allPairOrdInPtSCAss;
00356  
00357 
00358 }
00359 
00360 
00361 
00362 
00363 
00364 
00365 
00366
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoEgamma/EgammaPhotonAlgos/src/ConversionTrackPairFinder.cc
