Reference Manual

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaTools/interface/ConversionFinder.h"
00002 #include "DataFormats/Math/interface/deltaR.h"
00003 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00004 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "TMath.h"
00007 
00008 typedef math::XYZTLorentzVector LorentzVector;
00009 
00010 bool ConversionFinder::isFromConversion
00011  ( const ConversionInfo & convInfo, double maxAbsDist, double maxAbsDcot )
00012  {
00013   if ( (std::abs(convInfo.dist())<maxAbsDist) && (std::abs(convInfo.dcot())<maxAbsDcot) )
00014     return true ;
00015   return false ;
00016  }
00017 
00018 //-----------------------------------------------------------------------------
00019 ConversionFinder::ConversionFinder() {}
00020 
00021 //-----------------------------------------------------------------------------
00022 ConversionFinder::~ConversionFinder() {}
00023 
00024 //-----------------------------------------------------------------------------
00025 ConversionInfo ConversionFinder::getConversionInfo(const reco::GsfElectron& gsfElectron,
00026                                                                 const edm::Handle<reco::TrackCollection>& ctftracks_h,
00027                                                                 const edm::Handle<reco::GsfTrackCollection>& gsftracks_h,
00028                                                                 const double bFieldAtOrigin,
00029                                                                 const double minFracSharedHits) {
00030 
00031   std::vector<ConversionInfo> temp = getConversionInfos(*gsfElectron.core(),ctftracks_h,gsftracks_h,bFieldAtOrigin,minFracSharedHits) ;
00032   return findBestConversionMatch(temp);
00033 
00034 }
00035 //-----------------------------------------------------------------------------
00036 ConversionInfo ConversionFinder::getConversionInfo(const reco::GsfElectronCore& gsfElectron,
00037                                                     const edm::Handle<reco::TrackCollection>& ctftracks_h,
00038                                                     const edm::Handle<reco::GsfTrackCollection>& gsftracks_h,
00039                                                     const double bFieldAtOrigin,
00040                                                     const double minFracSharedHits) {
00041 
00042   std::vector<ConversionInfo> temp = getConversionInfos(gsfElectron,ctftracks_h,gsftracks_h,bFieldAtOrigin,minFracSharedHits) ;
00043   return findBestConversionMatch(temp);
00044 
00045 }
00046 
00047 
00048 //-----------------------------------------------------------------------------
00049 std::vector<ConversionInfo> ConversionFinder::getConversionInfos(const reco::GsfElectronCore& gsfElectron,
00050                                                                 const edm::Handle<reco::TrackCollection>& ctftracks_h,
00051                                                                 const edm::Handle<reco::GsfTrackCollection>& gsftracks_h,
00052                                                                 const double bFieldAtOrigin,
00053                                                                 const double minFracSharedHits) {
00054 
00055 
00056 
00057   using namespace reco;
00058   using namespace std;
00059   using namespace edm;
00060 
00061 
00062   //get the track collections
00063   const TrackCollection    *ctftracks   = ctftracks_h.product();
00064   const GsfTrackCollection *gsftracks   = gsftracks_h.product();
00065 
00066   //get the references to the gsf and ctf tracks that are made
00067   //by the electron
00068   const reco::TrackRef    el_ctftrack   = gsfElectron.ctfTrack();
00069   const reco::GsfTrackRef el_gsftrack   = gsfElectron.gsfTrack();
00070 
00071   //protect against the wrong collection being passed to the function
00072   if(el_ctftrack.isNonnull() && el_ctftrack.id() != ctftracks_h.id())
00073     throw cms::Exception("ConversionFinderError") << "ProductID of ctf track collection does not match ProductID of electron's CTF track! \n";
00074   if(el_gsftrack.isNonnull() && el_gsftrack.id() != gsftracks_h.id())
00075     throw cms::Exception("ConversionFinderError") << "ProductID of gsf track collection does not match ProductID of electron's GSF track! \n";
00076 
00077   //make p4s for the electron's tracks for use later
00078   LorentzVector el_ctftrack_p4;
00079   if(el_ctftrack.isNonnull() && gsfElectron.ctfGsfOverlap() > minFracSharedHits)
00080     el_ctftrack_p4 = LorentzVector(el_ctftrack->px(), el_ctftrack->py(), el_ctftrack->pz(), el_ctftrack->p());
00081   LorentzVector el_gsftrack_p4(el_gsftrack->px(), el_gsftrack->py(), el_gsftrack->pz(), el_gsftrack->p());
00082 
00083   //the electron's CTF track must share at least 45% of the inner hits
00084   //with the electron's GSF track
00085   int ctfidx = -999.;
00086   int gsfidx = -999.;
00087   if(el_ctftrack.isNonnull() && gsfElectron.ctfGsfOverlap() > minFracSharedHits)
00088     ctfidx = static_cast<int>(el_ctftrack.key());
00089 
00090   gsfidx = static_cast<int>(el_gsftrack.key());
00091 
00092 
00093   //these vectors are for those candidate partner tracks that pass our cuts
00094   vector<ConversionInfo> v_candidatePartners;
00095   //track indices required to make references
00096   int ctftk_i = 0;
00097   int gsftk_i = 0;
00098 
00099 
00100   //loop over the CTF tracks and try to find the partner track
00101   for(TrackCollection::const_iterator ctftk = ctftracks->begin();
00102       ctftk != ctftracks->end(); ctftk++, ctftk_i++) {
00103 
00104     if((ctftk_i == ctfidx))
00105       continue;
00106 
00107     //candidate track's p4
00108     LorentzVector ctftk_p4 = LorentzVector(ctftk->px(), ctftk->py(), ctftk->pz(), ctftk->p());
00109 
00110     //apply quality cuts to remove bad tracks
00111     if(ctftk->ptError()/ctftk->pt() > 0.05)
00112       continue;
00113     if(ctftk->numberOfValidHits() < 5)
00114       continue;
00115 
00116     if(el_ctftrack.isNonnull() &&
00117        gsfElectron.ctfGsfOverlap() > minFracSharedHits &&
00118        fabs(ctftk_p4.Pt() - el_ctftrack->pt())/el_ctftrack->pt() < 0.2)
00119       continue;
00120 
00121     //use the electron's CTF track, if not null, to search for the partner track
00122     //look only in a cone of 0.5 to save time, and require that the track is opp. sign
00123     if(el_ctftrack.isNonnull() && gsfElectron.ctfGsfOverlap() > minFracSharedHits &&
00124        deltaR(el_ctftrack_p4, ctftk_p4) < 0.5 &&
00125        (el_ctftrack->charge() + ctftk->charge() == 0) ) {
00126 
00127       ConversionInfo convInfo = getConversionInfo((const reco::Track*)(el_ctftrack.get()), &(*ctftk), bFieldAtOrigin);
00128 
00129       //need to add the track reference information for completeness
00130       //because the overloaded fnc above does not make a trackRef
00131       int deltaMissingHits = ctftk->trackerExpectedHitsInner().numberOfHits() - el_ctftrack->trackerExpectedHitsInner().numberOfHits();
00132       convInfo = ConversionInfo(convInfo.dist(),
00133                                 convInfo.dcot(),
00134                                 convInfo.radiusOfConversion(),
00135                                 convInfo.pointOfConversion(),
00136                                 TrackRef(ctftracks_h, ctftk_i),
00137                                 GsfTrackRef() ,
00138                                 deltaMissingHits,
00139                                 0);
00140 
00141       v_candidatePartners.push_back(convInfo);
00142 
00143     }//using the electron's CTF track
00144 
00145 
00146     //now we check using the electron's gsf track
00147     if(deltaR(el_gsftrack_p4, ctftk_p4) < 0.5 &&
00148        (el_gsftrack->charge() + ctftk->charge() == 0) &&
00149        el_gsftrack->ptError()/el_gsftrack->pt() < 0.25) {
00150 
00151       int deltaMissingHits = ctftk->trackerExpectedHitsInner().numberOfHits() - el_gsftrack->trackerExpectedHitsInner().numberOfHits();
00152       ConversionInfo convInfo = getConversionInfo((const reco::Track*)(el_gsftrack.get()), &(*ctftk), bFieldAtOrigin);
00153       convInfo = ConversionInfo(convInfo.dist(),
00154                                 convInfo.dcot(),
00155                                 convInfo.radiusOfConversion(),
00156                                 convInfo.pointOfConversion(),
00157                                 TrackRef(ctftracks_h, ctftk_i),
00158                                 GsfTrackRef(),
00159                                 deltaMissingHits,
00160                                 1);
00161 
00162       v_candidatePartners.push_back(convInfo);
00163     }//using the electron's GSF track
00164 
00165   }//loop over the CTF track collection
00166 
00167 
00168   //------------------------------------------------------ Loop over GSF collection ----------------------------------//
00169   for(GsfTrackCollection::const_iterator gsftk = gsftracks->begin();
00170       gsftk != gsftracks->end(); gsftk++, gsftk_i++) {
00171 
00172     //reject the electron's own gsfTrack
00173     if(gsfidx == gsftk_i)
00174       continue;
00175 
00176     LorentzVector gsftk_p4 = LorentzVector(gsftk->px(), gsftk->py(), gsftk->pz(), gsftk->p());
00177 
00178     //apply quality cuts to remove bad tracks
00179     if(gsftk->ptError()/gsftk->pt() > 0.5)
00180       continue;
00181     if(gsftk->numberOfValidHits() < 5)
00182       continue;
00183 
00184     if(fabs(gsftk->pt() - el_gsftrack->pt())/el_gsftrack->pt() < 0.25)
00185       continue;
00186 
00187     //try using the electron's CTF track first if it exists
00188     //look only in a cone of 0.5 around the electron's track
00189         //require opposite sign
00190     if(el_ctftrack.isNonnull() && gsfElectron.ctfGsfOverlap() > minFracSharedHits &&
00191        deltaR(el_ctftrack_p4, gsftk_p4) < 0.5 &&
00192        (el_ctftrack->charge() + gsftk->charge() == 0)) {
00193 
00194       int deltaMissingHits = gsftk->trackerExpectedHitsInner().numberOfHits() - el_ctftrack->trackerExpectedHitsInner().numberOfHits();
00195       ConversionInfo convInfo = getConversionInfo((const reco::Track*)(el_ctftrack.get()), (const reco::Track*)(&(*gsftk)), bFieldAtOrigin);
00196       //fill the Ref info
00197       convInfo = ConversionInfo(convInfo.dist(),
00198                                 convInfo.dcot(),
00199                                 convInfo.radiusOfConversion(),
00200                                 convInfo.pointOfConversion(),
00201                                 TrackRef(),
00202                                 GsfTrackRef(gsftracks_h, gsftk_i),
00203                                 deltaMissingHits,
00204                                 2);
00205       v_candidatePartners.push_back(convInfo);
00206 
00207     }
00208 
00209     //use the electron's gsf track
00210     if(deltaR(el_gsftrack_p4, gsftk_p4) < 0.5 &&
00211        (el_gsftrack->charge() + gsftk->charge() == 0) &&
00212        (el_gsftrack->ptError()/el_gsftrack_p4.pt() < 0.5)) {
00213       ConversionInfo convInfo = getConversionInfo((const reco::Track*)(el_gsftrack.get()), (const reco::Track*)(&(*gsftk)), bFieldAtOrigin);
00214       //fill the Ref info
00215 
00216       int deltaMissingHits = gsftk->trackerExpectedHitsInner().numberOfHits() - el_gsftrack->trackerExpectedHitsInner().numberOfHits();
00217       convInfo = ConversionInfo(convInfo.dist(),
00218                                 convInfo.dcot(),
00219                                 convInfo.radiusOfConversion(),
00220                                 convInfo.pointOfConversion(),
00221                                 TrackRef(),
00222                                 GsfTrackRef(gsftracks_h, gsftk_i),
00223                                 deltaMissingHits,
00224                                 3);
00225 
00226       v_candidatePartners.push_back(convInfo);
00227     }
00228   }//loop over the gsf track collection
00229 
00230 
00231   return v_candidatePartners;
00232 
00233 }
00234 
00235 
00236 //-------------------------------------------------------------------------------------
00237 ConversionInfo ConversionFinder::getConversionInfo(const reco::Track *el_track,
00238                                                    const reco::Track *candPartnerTk,
00239                                                    const double bFieldAtOrigin) {
00240 
00241   using namespace reco;
00242 
00243   //now calculate the conversion related information
00244   LorentzVector el_tk_p4(el_track->px(), el_track->py(), el_track->pz(), el_track->p());
00245   double elCurvature = -0.3*bFieldAtOrigin*(el_track->charge()/el_tk_p4.pt())/100.;
00246   double rEl = fabs(1./elCurvature);
00247   double xEl = -1*(1./elCurvature - el_track->d0())*sin(el_tk_p4.phi());
00248   double yEl = (1./elCurvature - el_track->d0())*cos(el_tk_p4.phi());
00249 
00250 
00251   LorentzVector cand_p4 = LorentzVector(candPartnerTk->px(), candPartnerTk->py(),candPartnerTk->pz(), candPartnerTk->p());
00252   double candCurvature = -0.3*bFieldAtOrigin*(candPartnerTk->charge()/cand_p4.pt())/100.;
00253   double rCand = fabs(1./candCurvature);
00254   double xCand = -1*(1./candCurvature - candPartnerTk->d0())*sin(cand_p4.phi());
00255   double yCand = (1./candCurvature - candPartnerTk->d0())*cos(cand_p4.phi());
00256 
00257   double d = sqrt(pow(xEl-xCand, 2) + pow(yEl-yCand , 2));
00258   double dist = d - (rEl + rCand);
00259   double dcot = 1./tan(el_tk_p4.theta()) - 1./tan(cand_p4.theta());
00260 
00261   //get the point of conversion
00262   double xa1 = xEl   + (xCand-xEl) * rEl/d;
00263   double xa2 = xCand + (xEl-xCand) * rCand/d;
00264   double ya1 = yEl   + (yCand-yEl) * rEl/d;
00265   double ya2 = yCand + (yEl-yCand) * rCand/d;
00266 
00267   double x=.5*(xa1+xa2);
00268   double y=.5*(ya1+ya2);
00269   double rconv = sqrt(pow(x,2) + pow(y,2));
00270   double z = el_track->dz() + rEl*el_track->pz()*TMath::ACos(1-pow(rconv,2)/(2.*pow(rEl,2)))/el_track->pt();
00271 
00272   math::XYZPoint convPoint(x, y, z);
00273 
00274   //now assign a sign to the radius of conversion
00275   float tempsign = el_track->px()*x + el_track->py()*y;
00276   tempsign = tempsign/fabs(tempsign);
00277   rconv = tempsign*rconv;
00278 
00279   //return an instance of ConversionInfo, but with a NULL track refs
00280   return ConversionInfo(dist, dcot, rconv, convPoint, TrackRef(), GsfTrackRef(), -9999, -9999);
00281 
00282 }
00283 
00284 //-------------------------------------------------------------------------------------
00285 const reco::Track* ConversionFinder::getElectronTrack(const reco::GsfElectron& electron, const float minFracSharedHits) {
00286 
00287   if(electron.closestCtfTrackRef().isNonnull() &&
00288      electron.shFracInnerHits() > minFracSharedHits)
00289     return (const reco::Track*)electron.closestCtfTrackRef().get();
00290 
00291   return (const reco::Track*)(electron.gsfTrack().get());
00292 }
00293 
00294 //------------------------------------------------------------------------------------
00295 
00296 //takes in a vector of candidate conversion partners
00297 //and arbitrates between them returning the one with the
00298 //smallest R=sqrt(dist*dist + dcot*dcot)
00299 ConversionInfo ConversionFinder::arbitrateConversionPartnersbyR(const std::vector<ConversionInfo>& v_convCandidates) {
00300 
00301   if(v_convCandidates.size() == 1)
00302     return v_convCandidates.at(0);
00303 
00304   ConversionInfo arbitratedConvInfo = v_convCandidates.at(0);
00305   double R = sqrt(pow(arbitratedConvInfo.dist(),2) + pow(arbitratedConvInfo.dcot(),2));
00306 
00307   for(unsigned int i = 1; i < v_convCandidates.size(); i++) {
00308     ConversionInfo temp = v_convCandidates.at(i);
00309     double temp_R = sqrt(pow(temp.dist(),2) + pow(temp.dcot(),2));
00310     if(temp_R < R) {
00311       R = temp_R;
00312       arbitratedConvInfo = temp;
00313     }
00314 
00315   }
00316 
00317   return arbitratedConvInfo;
00318 
00319  }
00320 
00321 //------------------------------------------------------------------------------------
00322 ConversionInfo ConversionFinder::findBestConversionMatch(const std::vector<ConversionInfo>& v_convCandidates)
00323  {
00324   using namespace std;
00325 
00326   if(v_convCandidates.size() == 0)
00327     return   ConversionInfo(-9999.,-9999.,-9999.,
00328                             math::XYZPoint(-9999.,-9999.,-9999),
00329                             reco::TrackRef(), reco::GsfTrackRef(),
00330                             -9999, -9999);
00331 
00332 
00333   if(v_convCandidates.size() == 1)
00334     return v_convCandidates.at(0);
00335 
00336   vector<ConversionInfo> v_0;
00337   vector<ConversionInfo> v_1;
00338   vector<ConversionInfo> v_2;
00339   vector<ConversionInfo> v_3;
00340   //loop over the candidates
00341   for(unsigned int i = 1; i < v_convCandidates.size(); i++) {
00342     ConversionInfo temp = v_convCandidates.at(i);
00343 
00344     if(temp.flag() == 0) {
00345       bool isConv = false;
00346       if(fabs(temp.dist()) < 0.02 &&
00347          fabs(temp.dcot()) < 0.02 &&
00348          temp.deltaMissingHits() < 3 &&
00349          temp.radiusOfConversion() > -2)
00350         isConv = true;
00351       if(sqrt(pow(temp.dist(),2) + pow(temp.dcot(),2)) < 0.05  &&
00352          temp.deltaMissingHits() < 2 &&
00353          temp.radiusOfConversion() > -2)
00354         isConv = true;
00355 
00356       if(isConv)
00357         v_0.push_back(temp);
00358     }
00359 
00360     if(temp.flag() == 1) {
00361 
00362       if(sqrt(pow(temp.dist(),2) + pow(temp.dcot(),2)) < 0.05  &&
00363          temp.deltaMissingHits() < 2 &&
00364          temp.radiusOfConversion() > -2)
00365         v_1.push_back(temp);
00366     }
00367     if(temp.flag() == 2) {
00368 
00369       if(sqrt(pow(temp.dist(),2) + pow(temp.dcot()*temp.dcot(),2)) < 0.05 &&
00370          temp.deltaMissingHits() < 2 &&
00371          temp.radiusOfConversion() > -2)
00372         v_2.push_back(temp);
00373 
00374     }
00375     if(temp.flag() == 3) {
00376 
00377       if(sqrt(temp.dist()*temp.dist() + temp.dcot()*temp.dcot()) < 0.05
00378          && temp.deltaMissingHits() < 2
00379          && temp.radiusOfConversion() > -2)
00380         v_3.push_back(temp);
00381 
00382     }
00383 
00384   }//candidate conversion loop
00385 
00386   //now do some arbitration
00387 
00388   //give preference to conversion partners found in the CTF collection
00389   //using the electron's CTF track
00390   if(v_0.size() > 0)
00391     return arbitrateConversionPartnersbyR(v_0);
00392 
00393   if(v_1.size() > 0)
00394     return arbitrateConversionPartnersbyR(v_1);
00395 
00396   if(v_2.size() > 0)
00397     return arbitrateConversionPartnersbyR(v_2);
00398 
00399   if(v_3.size() > 0)
00400     return arbitrateConversionPartnersbyR(v_3);
00401 
00402 
00403   //if we get here, we didn't find a candidate conversion partner that
00404   //satisfied even the loose selections
00405   //return the the closest partner by R
00406   return arbitrateConversionPartnersbyR(v_convCandidates);
00407 
00408  }
00409 
00410 
00411 
00412 
00413 //------------------------------------------------------------------------------------
00414 
00415 //------------------------------------------------------------------------------------
00416 // Exists here for backwards compatibility only. Provides only the dist and dcot
00417 std::pair<double, double> ConversionFinder::getConversionInfo(LorentzVector trk1_p4,
00418                                                               int trk1_q, float trk1_d0,
00419                                                               LorentzVector trk2_p4,
00420                                                               int trk2_q, float trk2_d0,
00421                                                               float bFieldAtOrigin) {
00422 
00423 
00424   double tk1Curvature = -0.3*bFieldAtOrigin*(trk1_q/trk1_p4.pt())/100.;
00425   double rTk1 = fabs(1./tk1Curvature);
00426   double xTk1 = -1.*(1./tk1Curvature - trk1_d0)*sin(trk1_p4.phi());
00427   double yTk1 = (1./tk1Curvature - trk1_d0)*cos(trk1_p4.phi());
00428 
00429 
00430   double tk2Curvature = -0.3*bFieldAtOrigin*(trk2_q/trk2_p4.pt())/100.;
00431   double rTk2 = fabs(1./tk2Curvature);
00432   double xTk2 = -1.*(1./tk2Curvature - trk2_d0)*sin(trk2_p4.phi());
00433   double yTk2 = (1./tk2Curvature - trk2_d0)*cos(trk2_p4.phi());
00434 
00435 
00436   double dist = sqrt(pow(xTk1-xTk2, 2) + pow(yTk1-yTk2 , 2));
00437   dist = dist - (rTk1 + rTk2);
00438 
00439   double dcot = 1./tan(trk1_p4.theta()) - 1./tan(trk2_p4.theta());
00440 
00441   return std::make_pair(dist, dcot);
00442 
00443  }
00444 
00445 
00446 //-------------------------------------- Also for backwards compatibility reasons  ---------------------------------------------------------------
00447 ConversionInfo ConversionFinder::getConversionInfo(const reco::GsfElectron& gsfElectron,
00448                                                    const edm::Handle<reco::TrackCollection>& track_h,
00449                                                    const double bFieldAtOrigin,
00450                                                    const double minFracSharedHits) {
00451 
00452 
00453   using namespace reco;
00454   using namespace std;
00455   using namespace edm;
00456 
00457 
00458   const TrackCollection *ctftracks = track_h.product();
00459   const reco::TrackRef el_ctftrack = gsfElectron.closestCtfTrackRef();
00460   int ctfidx = -999.;
00461   int flag = -9999.;
00462   if(el_ctftrack.isNonnull() && gsfElectron.shFracInnerHits() > minFracSharedHits) {
00463     ctfidx = static_cast<int>(el_ctftrack.key());
00464     flag = 0;
00465   } else
00466     flag = 1;
00467 
00468 
00469   /*
00470     determine whether we're going to use the CTF track or the GSF track
00471     using the electron's CTF track to find the dist, dcot has been shown
00472     to reduce the inefficiency
00473   */
00474   const reco::Track* el_track = getElectronTrack(gsfElectron, minFracSharedHits);
00475   LorentzVector el_tk_p4(el_track->px(), el_track->py(), el_track->pz(), el_track->p());
00476 
00477 
00478   int tk_i = 0;
00479   double mindcot = 9999.;
00480   //make a null Track Ref
00481   TrackRef candCtfTrackRef = TrackRef() ;
00482 
00483 
00484   for(TrackCollection::const_iterator tk = ctftracks->begin();
00485       tk != ctftracks->end(); tk++, tk_i++) {
00486     //if the general Track is the same one as made by the electron, skip it
00487     if((tk_i == ctfidx))
00488       continue;
00489 
00490     LorentzVector tk_p4 = LorentzVector(tk->px(), tk->py(),tk->pz(), tk->p());
00491 
00492     //look only in a cone of 0.5
00493     double dR = deltaR(el_tk_p4, tk_p4);
00494     if(dR>0.5)
00495       continue;
00496 
00497 
00498     //require opp. sign -> Should we use the majority logic??
00499     if(tk->charge() + el_track->charge() != 0)
00500       continue;
00501 
00502     double dcot = fabs(1./tan(tk_p4.theta()) - 1./tan(el_tk_p4.theta()));
00503     if(dcot < mindcot) {
00504       mindcot = dcot;
00505       candCtfTrackRef = reco::TrackRef(track_h, tk_i);
00506     }
00507   }//track loop
00508 
00509 
00510   if(!candCtfTrackRef.isNonnull())
00511     return ConversionInfo(-9999.,-9999.,-9999.,
00512                           math::XYZPoint(-9999.,-9999.,-9999),
00513                           reco::TrackRef(), reco::GsfTrackRef(),
00514                           -9999, -9999);
00515 
00516 
00517 
00518   //now calculate the conversion related information
00519   double elCurvature = -0.3*bFieldAtOrigin*(el_track->charge()/el_tk_p4.pt())/100.;
00520   double rEl = fabs(1./elCurvature);
00521   double xEl = -1*(1./elCurvature - el_track->d0())*sin(el_tk_p4.phi());
00522   double yEl = (1./elCurvature - el_track->d0())*cos(el_tk_p4.phi());
00523 
00524 
00525   LorentzVector cand_p4 = LorentzVector(candCtfTrackRef->px(), candCtfTrackRef->py(),candCtfTrackRef->pz(), candCtfTrackRef->p());
00526   double candCurvature = -0.3*bFieldAtOrigin*(candCtfTrackRef->charge()/cand_p4.pt())/100.;
00527   double rCand = fabs(1./candCurvature);
00528   double xCand = -1*(1./candCurvature - candCtfTrackRef->d0())*sin(cand_p4.phi());
00529   double yCand = (1./candCurvature - candCtfTrackRef->d0())*cos(cand_p4.phi());
00530 
00531   double d = sqrt(pow(xEl-xCand, 2) + pow(yEl-yCand , 2));
00532   double dist = d - (rEl + rCand);
00533   double dcot = 1./tan(el_tk_p4.theta()) - 1./tan(cand_p4.theta());
00534 
00535   //get the point of conversion
00536   double xa1 = xEl   + (xCand-xEl) * rEl/d;
00537   double xa2 = xCand + (xEl-xCand) * rCand/d;
00538   double ya1 = yEl   + (yCand-yEl) * rEl/d;
00539   double ya2 = yCand + (yEl-yCand) * rCand/d;
00540 
00541   double x=.5*(xa1+xa2);
00542   double y=.5*(ya1+ya2);
00543   double rconv = sqrt(pow(x,2) + pow(y,2));
00544   double z = el_track->dz() + rEl*el_track->pz()*TMath::ACos(1-pow(rconv,2)/(2.*pow(rEl,2)))/el_track->pt();
00545 
00546   math::XYZPoint convPoint(x, y, z);
00547 
00548   //now assign a sign to the radius of conversion
00549   float tempsign = el_track->px()*x + el_track->py()*y;
00550   tempsign = tempsign/fabs(tempsign);
00551   rconv = tempsign*rconv;
00552 
00553   int deltaMissingHits = -9999;
00554   deltaMissingHits = candCtfTrackRef->trackerExpectedHitsInner().numberOfHits() - el_track->trackerExpectedHitsInner().numberOfHits();
00555   return ConversionInfo(dist, dcot, rconv, convPoint, candCtfTrackRef, GsfTrackRef(), deltaMissingHits, flag);
00556 
00557  }
00558 
00559 //-------------------------------------------------------------------------------------
/data/refman/pasoursint/CMSSW_5_3_3/src/RecoEgamma/EgammaTools/src/ConversionFinder.cc
