CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoEgamma/EgammaTools/src/ggPFTracks.cc

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaTools/interface/ggPFTracks.h"
00002 //Class by Rishi Patel rpatel@cern.ch for Single Leg (high pt leg of Conversion
00003 //For Vertex Selection, also has functions to do pointing using conversion pairs//that are stored with PFPhotons/PFElectrons
00004 
00005 ggPFTracks::ggPFTracks(
00006                        edm::Handle<BeamSpot>& beamSpotHandle
00007                        ):
00008   beamSpotHandle_(beamSpotHandle),//pass beamspot
00009   isConv_(false)
00010 {
00011   
00012 }
00013 
00014 ggPFTracks::~ggPFTracks(){}
00015 //Fills a Vector of Conversion Tracks (single legs, or track pairs from a PFPhoton) and also Conversion Objects
00016 void ggPFTracks::getPFConvTracks(
00017                                  reco::Photon phot,
00018                                  vector<edm::RefToBase<reco::Track> > &Tracks, 
00019                                  reco::ConversionRefVector &conversions,
00020                                  vector<edm::RefToBase<reco::Track> > &SLTracks, 
00021                                  reco::ConversionRefVector &SLconversions
00022                                  ){
00023   Tracks.clear();
00024   conversions.clear();
00025   SLTracks.clear();
00026   SLconversions.clear();
00027   conversions=phot.conversions();
00028   //loop over conversion paris
00029   for(unsigned int c=0; c<conversions.size(); ++c){
00030     const std::vector<edm::RefToBase<reco::Track> > tracks = conversions[c]->tracks();
00031     for (unsigned int t=0; t<tracks.size(); t++){
00032       Tracks.push_back(tracks[t]);
00033     }
00034   }
00035   
00036   reco::ConversionRefVector  SLConversions=phot.conversionsOneLeg();
00037   
00038   for(unsigned int SLc=0; SLc<SLConversions.size(); ++SLc){
00039     const std::vector<edm::RefToBase<reco::Track> > tracks = SLConversions[SLc]->tracks();
00040     SLconversions.push_back(SLConversions[SLc]);
00041     
00042     for (unsigned int t=0; t<tracks.size(); t++){
00043       SLTracks.push_back(tracks[t]);
00044     }
00045   }
00046   if(SLConversions.size()>0 || conversions.size()>0)isConv_=true;
00047 }
00048 //This function does SuperCluster Pointing using Single Legs and Conversion pairs
00049 //There is a bool that is flagged when it uses SLconv exclusively
00050 std::pair<float,float> ggPFTracks::BeamLineInt( 
00051                                                reco::SuperClusterRef sc, 
00052                                                vector<edm::RefToBase<reco::Track> > &Tracks, 
00053                                                reco::ConversionRefVector &conversions,
00054                                                vector<edm::RefToBase<reco::Track> > &SLTracks, 
00055                                                reco::ConversionRefVector &SLconversions
00056                                                ){
00057   std::pair<float,float> Zint(0,0);
00058   TVector3 beamspot(beamSpotHandle_->position().x(),beamSpotHandle_->position().y(),
00059               beamSpotHandle_->position().z());
00060   TVector3 SCPos(sc->position().x()-beamspot[0], sc->position().y()-beamspot[1], sc->position().z()-beamspot[2]);
00061   //find min Conversion radius for track pairs if they exist. 
00062   float convRMin=130;
00063   float SLConvRMin=130;
00064   int c_index=-1; int SLc_index=-1;
00065   if(conversions.size()>0)
00066     {
00067       for(unsigned int c=0; c<conversions.size(); ++c){
00068         float convR=sqrt(conversions[c]->conversionVertex().x()* conversions[c]->conversionVertex().x() + conversions[c]->conversionVertex().y()* conversions[c]->conversionVertex().y());
00069         if(convRMin>convR){
00070           convRMin=convR;
00071           c_index=c;
00072         }
00073       }
00074     }
00075   
00076   if(SLconversions.size()>0){
00077     for(unsigned int SLc=0; SLc<SLconversions.size(); ++SLc){      
00078       std::vector<math::XYZPointF> innerPos=SLconversions[SLc]->tracksInnerPosition();
00079       for (unsigned int t=0; t<innerPos.size(); t++){
00080         float convR=sqrt( innerPos[t].X()* innerPos[t].X() + innerPos[t].Y()* innerPos[t].Y());
00081         
00082         if(SLConvRMin>convR){
00083           SLc_index=SLc;  
00084           SLConvRMin=convR;
00085           
00086         }
00087       }
00088     }
00089   }
00090   TVector3 TkPos(beamspot[0],beamspot[1],beamspot[2]);
00091   //take the smaller Radius:
00092   if(convRMin<SLConvRMin && c_index>-1){
00093     //point using Conversion Vertex
00094     TkPos.SetXYZ(conversions[c_index]->conversionVertex().x()-beamspot[0],
00095                  conversions[c_index]->conversionVertex().y()-beamspot[1],
00096                  conversions[c_index]->conversionVertex().z()-beamspot[2] );
00097   }
00098   if(SLConvRMin<convRMin && SLc_index>-1){
00099     reco::Vertex conv=SLconversions[SLc_index]->conversionVertex();
00100     TkPos.SetXYZ(conv.x()-beamspot[0],
00101                  conv.y()-beamspot[1],
00102                  conv.z()-beamspot[2] );
00103   }
00104   //Intersection fromt the two points:
00105   float R1=sqrt(SCPos.X()* SCPos.X() + SCPos.Y()*SCPos.Y()); 
00106   float R2=sqrt(TkPos.X()* TkPos.X() + TkPos.Y()*TkPos.Y());
00107   float Z1=SCPos.Z();
00108   float Z2=TkPos.Z();
00109   float slope=(Z1-Z2)/(R1-R2);
00110   Zint.first=Z2 - R2*slope;
00111   //determine error based on Tracking Region for conversion pairs:
00112   float sigmaPix=0.06;
00113   float sigmaTib=0.67;
00114   float sigmaTob=2.04;
00115   float sigmaFwd1=0.18;
00116   float sigmaFwd2=0.61;
00117   float sigmaFwd3=0.99;
00118   //error of SL tracks of Conversion based on EB/EE and R>39 R<39 (4 cat)
00119   float EBLR=0.24;
00120   float EBHR=0.478;
00121   float EELR=0.416;
00122   float EEHR=0.888;
00123   if(sc->eta()<1.4442){//Barrel
00124     //if conversion
00125     if(convRMin<SLConvRMin && c_index>-1){//3 tracking regions
00126       if(convRMin<=15)Zint.second=sigmaPix;
00127       else if(convRMin>15 && convRMin<=60)Zint.second=sigmaTib;
00128       else Zint.second=sigmaTob;
00129     }  
00130     //if SL 
00131     if(SLConvRMin<convRMin && SLc_index>-1){//2 tracking regions
00132       if(SLConvRMin<39)Zint.second=EBLR;
00133       else Zint.second=EBHR;
00134     }
00135   }
00136   else{
00137     //if conversion
00138     if(convRMin<SLConvRMin && c_index>-1){//3 foreward tracking regions
00139       float convz=conversions[c_index]->conversionVertex().z();
00140       if(convz<=50)Zint.second=sigmaFwd1;
00141       else if(convz>50 && convz<=60)Zint.second=sigmaFwd2;
00142       else Zint.second=sigmaFwd3;
00143     }
00144     //if SL
00145     if(SLConvRMin<convRMin && SLc_index>-1){//2 tracking regions based on R
00146       if(SLConvRMin<39)Zint.second=EELR;
00147       else Zint.second=EEHR;
00148     }
00149   }
00150   return Zint;//return intersection at beamline and error in the pointing based on tracking region
00151 }
00152 //This Function does the track Projection (using innermost hit and inner momentum of the Single Track, NOTE: Can also use for GSF tracks
00153 std::pair<float,float> ggPFTracks::TrackProj(
00154                                              bool isEb, 
00155                                              reco::GsfTrackRef gsf,
00156                                              vector<edm::RefToBase<reco::Track> > &SLTracks, 
00157                                              reco::ConversionRefVector &SLconversions
00158                             ){
00159   std::pair<float,float> ZProj(0,0);
00160   
00161   if(gsf.isNonnull()){//if there is a gsf track then use this for track projection Plenty of inner hits
00162     ZProj=gsfTrackProj(gsf);
00163     return ZProj;
00164   }
00165   
00166   float minR=210;
00167   int SLind=-1;
00168   if(SLconversions.size()>0){//find track with min Radius (starts earliest)
00169     for(unsigned int SLc=0; SLc<SLconversions.size(); ++SLc){     
00170       reco::Vertex conv=SLconversions[SLc]->conversionVertex();
00171       float convR=sqrt( conv.x() * conv.x() + conv.y() * conv.y());
00172       if(convR<minR){
00173         minR=convR;
00174         SLind=SLc;
00175       }
00176     }
00177     reco::Vertex conv=SLconversions[SLind]->conversionVertex();
00178     const std::vector<edm::RefToBase<reco::Track> > tracks = SLconversions[SLind]->tracks();
00179     float theta =tracks[0]->theta();
00180     float tkz=conv.z();
00181     float tkR=sqrt( conv.x() * conv.x() + conv.y() * conv.y());
00182     float thetErr=tracks[0]->thetaError();
00183     float Z=tkz-tkR/tan(theta);//track projection
00184     float Zerr=((-1*(cos(theta)*cos(theta))/(sin(theta)* sin(theta))-1)*tkR*thetErr); //projection error based on theta err of track theta derivative of the Zproj
00185     //for early tracks theta error is very small so just hard code an error 
00186     //that is measured by looking at the track Proj resolutin in MC
00187     if(tkR<39 && isEb)Zerr=0.234; 
00188     if(tkR<39 && !isEb)Zerr=0.341;
00189     ZProj.first=Z; ZProj.second=Zerr;
00190     return ZProj;
00191   }
00192   return ZProj;
00193 }
00194 
00195 //this function combines the results from the Track Projection and SC pointing 
00196 //can also use conversion pairs and gsf tracks, even when there is no Single leg
00197 std::pair<float,float> ggPFTracks::gsfTrackProj(
00198                                                 reco::GsfTrackRef gsf
00199                                                 ){
00200   //if there is a gsf track then use this for track projection Plenty of inner hits
00201   
00202   float theta =gsf->innerMomentum().theta();
00203   
00204   float tkz=gsf->innerPosition().Z();
00205   float tkR=sqrt(gsf->innerPosition().X()* gsf->innerPosition().X()+ gsf->innerPosition().Y()* gsf->innerPosition().Y());
00206   float thetErr=gsf->thetaError();
00207   float Z=tkz-tkR/tan(theta);
00208   float Zerr=((-1*(cos(theta)*cos(theta))/(sin(theta)* sin(theta))-1)*tkR*thetErr);
00209   std::pair<float,float> ZProj(0,0);
00210   ZProj.first=Z; ZProj.second=Zerr;
00211   return ZProj;
00212   
00213 }
00214 
00215 std::pair<float,float> ggPFTracks::gsfElectronProj(
00216                                                    reco::GsfElectron gsf
00217                                                    ){
00218   //if there is a gsf electron then use this for track projection Plenty of inner hits
00219   
00220   float theta =gsf.trackMomentumAtVtx().theta();
00221   
00222   float tkz=gsf.trackPositionAtVtx().Z();
00223   float tkR=sqrt(gsf.trackPositionAtVtx().X()* gsf.trackPositionAtVtx().X()+ gsf.trackPositionAtVtx().Y()* gsf.trackPositionAtVtx().Y());
00224   float thetErr=gsf.gsfTrack()->thetaError();
00225   float Z=tkz-tkR/tan(theta);
00226   float Zerr=((-1*(cos(theta)*cos(theta))/(sin(theta)* sin(theta))-1)*tkR*thetErr);
00227   std::pair<float,float> ZProj(0,0);
00228   ZProj.first=Z; ZProj.second=Zerr;
00229   return ZProj;
00230   
00231 }
00232 
00233 std::pair<float,float> ggPFTracks::CombZVtx(
00234                            reco::SuperClusterRef sc, 
00235                            reco::GsfTrackRef gsf,
00236                            vector<edm::RefToBase<reco::Track> > &Tracks, 
00237                            reco::ConversionRefVector &conversions,
00238                            vector<edm::RefToBase<reco::Track> > &SLTracks, 
00239                            reco::ConversionRefVector &SLconversions
00240                            
00241                            ){
00242   std::pair<float, float> combZ(0,0);
00243   bool isEb;
00244   TVector3 beamspot(beamSpotHandle_->position().x(),beamSpotHandle_->position().y(),
00245                     beamSpotHandle_->position().z());
00246   if(fabs(sc->eta())<1.4442)isEb=true;
00247   else isEb=false;
00248   std::pair< float,float> SCZ=BeamLineInt(sc, Tracks, conversions,SLTracks, SLconversions);
00249   std::pair<float, float> TkPjZ=TrackProj(isEb,gsf,SLTracks, SLconversions);
00250   //errors in the two methods
00251   float sigZProj=TkPjZ.second;
00252   float sigSCPoint=SCZ.second;
00253   
00254   if(gsf.isNonnull()){combZ=TkPjZ; return combZ;}//for gsf Tracks just return track Proj
00255   //weighted avg of the two methods, where weights are based on the error
00256   float Z=((SCZ.first/(sigSCPoint*sigSCPoint))+ (TkPjZ.first/( sigZProj* sigZProj)))/(1/(sigSCPoint * sigSCPoint)+ 1/(sigZProj * sigZProj));
00257   //total error sum of the two errors in quadrature
00258   float sigZ=sqrt((sigSCPoint * sigSCPoint)+ (sigZProj * sigZProj));
00259   combZ.first=Z; combZ.second=sigZ;
00260   return combZ;  
00261 }
00262 //this function combines the results from the Track Projection and SC pointing 
00263 //simpler version of the above function, but returns bspot Z when there is no SL
00264 std::pair<float, float> ggPFTracks::SLCombZVtx(
00265                                                reco::Photon phot,
00266                                                bool &hasSL
00267                                                ){
00268   std::pair<float, float> combZ(0,0);
00269   TVector3 beamspot(beamSpotHandle_->position().x(),beamSpotHandle_->position().y(),
00270                     beamSpotHandle_->position().z());
00271   bool isEb=phot.isEB();
00272   //only want to use SL so the conversion pair variables will be dummy for the get function
00273   vector<edm::RefToBase<reco::Track> >convTracks; 
00274   reco::ConversionRefVector pairConv;
00275   
00276   reco::ConversionRefVector SLConversions;
00277   vector<edm::RefToBase<reco::Track> >SLTracks; 
00278   reco::GsfTrackRef gsfdummy;
00279   getPFConvTracks(phot,convTracks,pairConv ,SLTracks,SLConversions);
00280   
00281   if(SLConversions.size()>0){
00282     hasSL=true;
00283     vector<edm::RefToBase<reco::Track> >dummy; 
00284     reco::ConversionRefVector pairdummy;
00285     std::pair<float, float> TkPjZ=TrackProj(isEb,gsfdummy,SLTracks, SLConversions);
00286     std::pair< float,float> SCZ=BeamLineInt(phot.superCluster(), dummy,pairdummy,SLTracks, SLConversions);
00287     //errors in the two methods
00288     float sigZProj=TkPjZ.second;
00289     float sigSCPoint=SCZ.second;
00290     //weighted avg of the two methods, where weights are based on the error
00291     float Z=((SCZ.first/(sigSCPoint*sigSCPoint))+ (TkPjZ.first/(sigZProj * sigZProj)))/(1/(sigSCPoint*sigSCPoint)+ 1/(sigZProj * sigZProj));
00292     //total error sum of the two errors in quadrature
00293     float sigZ=sqrt((sigSCPoint*sigSCPoint)+ (sigZProj * sigZProj));
00294     combZ.first=Z; combZ.second=sigZ;
00295   }
00296   else if(pairConv.size()>0){//else use Conversion pairs if available 
00297     std::pair< float,float> SCZ=BeamLineInt(phot.superCluster(),convTracks, pairConv, SLTracks, SLConversions);   
00298     combZ.first=SCZ.first;
00299     combZ.second=SCZ.second;
00300     hasSL=false;
00301   }
00302   else{//returns beamspot
00303     combZ.first=beamspot.Z();
00304     combZ.second=0;
00305     hasSL=false;
00306   }
00307   return combZ;
00308 }