00001 #include "RecoEgamma/EgammaTools/interface/ggPFTracks.h"
00002
00003
00004
00005 ggPFTracks::ggPFTracks(
00006 edm::Handle<BeamSpot>& beamSpotHandle
00007 ):
00008 beamSpotHandle_(beamSpotHandle),
00009 isConv_(false)
00010 {
00011
00012 }
00013
00014 ggPFTracks::~ggPFTracks(){}
00015
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
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
00049
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
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
00092 if(convRMin<SLConvRMin && c_index>-1){
00093
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
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
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
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){
00124
00125 if(convRMin<SLConvRMin && c_index>-1){
00126 if(convRMin<=15)Zint.second=sigmaPix;
00127 else if(convRMin>15 && convRMin<=60)Zint.second=sigmaTib;
00128 else Zint.second=sigmaTob;
00129 }
00130
00131 if(SLConvRMin<convRMin && SLc_index>-1){
00132 if(SLConvRMin<39)Zint.second=EBLR;
00133 else Zint.second=EBHR;
00134 }
00135 }
00136 else{
00137
00138 if(convRMin<SLConvRMin && c_index>-1){
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
00145 if(SLConvRMin<convRMin && SLc_index>-1){
00146 if(SLConvRMin<39)Zint.second=EELR;
00147 else Zint.second=EEHR;
00148 }
00149 }
00150 return Zint;
00151 }
00152
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()){
00162 ZProj=gsfTrackProj(gsf);
00163 return ZProj;
00164 }
00165
00166 float minR=210;
00167 int SLind=-1;
00168 if(SLconversions.size()>0){
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);
00184 float Zerr=((-1*(cos(theta)*cos(theta))/(sin(theta)* sin(theta))-1)*tkR*thetErr);
00185
00186
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
00196
00197 std::pair<float,float> ggPFTracks::gsfTrackProj(
00198 reco::GsfTrackRef gsf
00199 ){
00200
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
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
00251 float sigZProj=TkPjZ.second;
00252 float sigSCPoint=SCZ.second;
00253
00254 if(gsf.isNonnull()){combZ=TkPjZ; return combZ;}
00255
00256 float Z=((SCZ.first/(sigSCPoint*sigSCPoint))+ (TkPjZ.first/( sigZProj* sigZProj)))/(1/(sigSCPoint * sigSCPoint)+ 1/(sigZProj * sigZProj));
00257
00258 float sigZ=sqrt((sigSCPoint * sigSCPoint)+ (sigZProj * sigZProj));
00259 combZ.first=Z; combZ.second=sigZ;
00260 return combZ;
00261 }
00262
00263
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
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
00288 float sigZProj=TkPjZ.second;
00289 float sigSCPoint=SCZ.second;
00290
00291 float Z=((SCZ.first/(sigSCPoint*sigSCPoint))+ (TkPjZ.first/(sigZProj * sigZProj)))/(1/(sigSCPoint*sigSCPoint)+ 1/(sigZProj * sigZProj));
00292
00293 float sigZ=sqrt((sigSCPoint*sigSCPoint)+ (sigZProj * sigZProj));
00294 combZ.first=Z; combZ.second=sigZ;
00295 }
00296 else if(pairConv.size()>0){
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{
00303 combZ.first=beamspot.Z();
00304 combZ.second=0;
00305 hasSL=false;
00306 }
00307 return combZ;
00308 }