29 for(
int iv=0; iv<(
int)WrkVrtSet->size(); iv++) {
31 <<
"= "<<(*WrkVrtSet)[iv].vertex[0]
32 <<
", "<<(*WrkVrtSet)[iv].vertex[1]
33 <<
", "<<(*WrkVrtSet)[iv].vertex[2]
34 <<
" NTrk="<<(*WrkVrtSet)[iv].selTrk.size()
35 <<
" is good="<<std::boolalpha<<(*WrkVrtSet)[iv].Good<<std::noboolalpha
36 <<
" Chi2="<<(*WrkVrtSet)[iv].chi2
37 <<
" Mass="<<(*WrkVrtSet)[iv].vertexMom.M()
38 <<
" detached="<<(*WrkVrtSet)[iv].detachedTrack
39 <<
" proj.dist="<<(*WrkVrtSet)[iv].projectedVrt
41 for(
int kk=0;
kk<(
int)(*WrkVrtSet)[iv].selTrk.size();
kk++) {
42 msg(MSG::INFO)<<
", "<<(*WrkVrtSet)[iv].selTrk[
kk];}
43 for(
int kk=0;
kk<(
int)(*WrkVrtSet)[iv].selTrk.size();
kk++) {
44 msg(MSG::INFO)<<
", "<<
momAtVrt((*WrkVrtSet)[iv].trkAtVrt[
kk]).Perp();}
46 if((*WrkVrtSet)[iv].Good)nGoodV++;
55 TVector3 SV_PV( SV.x()-PV.
x(), SV.y()-PV.
y(), SV.z()-PV.
z() );
56 return Direction.Vect().Unit()*SV_PV.Unit();
60 const std::vector<double>& secVrtErr,
double& signif)
63 double distx = primVrt.
x()- secVrt.x();
64 double disty = primVrt.
y()- secVrt.y();
65 double distz = primVrt.
z()- secVrt.z();
69 primCovMtx(0,0) += secVrtErr[0];
70 primCovMtx(0,1) += secVrtErr[1];
71 primCovMtx(1,0) += secVrtErr[1];
72 primCovMtx(1,1) += secVrtErr[2];
73 primCovMtx(0,2) += secVrtErr[3];
74 primCovMtx(2,0) += secVrtErr[3];
75 primCovMtx(1,2) += secVrtErr[4];
76 primCovMtx(2,1) += secVrtErr[4];
77 primCovMtx(2,2) += secVrtErr[5];
81 signif = distx*wgtMtx(0,0)*distx
82 +disty*wgtMtx(1,1)*disty
83 +distz*wgtMtx(2,2)*distz
84 +2.*distx*wgtMtx(0,1)*disty
85 +2.*distx*wgtMtx(0,2)*distz
86 +2.*disty*wgtMtx(1,2)*distz;
87 signif=std::sqrt(std::abs(signif));
88 if( signif!=signif ) signif = 0.;
89 return std::sqrt(distx*distx+disty*disty+distz*distz);
93 const std::vector<double>& secVrtErr,
double& signif)
96 double distx = primVrt.
x()- secVrt.x();
97 double disty = primVrt.
y()- secVrt.y();
100 AmgSymMatrix(3) primCovMtx=primVrt.covariancePosition();
102 covMtx(0,0) = primCovMtx(0,0) + secVrtErr[0];
103 covMtx(0,1) = primCovMtx(0,1) + secVrtErr[1];
104 covMtx(1,0) = primCovMtx(1,0) + secVrtErr[1];
105 covMtx(1,1) = primCovMtx(1,1) + secVrtErr[2];
109 signif = distx*wgtMtx(0,0)*distx
110 +disty*wgtMtx(1,1)*disty
111 +2.*distx*wgtMtx(0,1)*disty;
112 signif=std::sqrt(std::abs(signif));
113 if( signif!=signif ) signif = 0.;
114 return std::sqrt(distx*distx+disty*disty);
119 const Amg::Vector3D & vrt2,
const std::vector<double> & vrtErr2)
122 double distx = vrt1.x()- vrt2.x();
123 double disty = vrt1.y()- vrt2.y();
124 double distz = vrt1.z()- vrt2.z();
127 primCovMtx(0,0) = vrtErr1[0]+vrtErr2[0];
128 primCovMtx(0,1) = primCovMtx(1,0) = vrtErr1[1]+vrtErr2[1];
129 primCovMtx(1,1) = vrtErr1[2]+vrtErr2[2];
130 primCovMtx(0,2) = primCovMtx(2,0) = vrtErr1[3]+vrtErr2[3];
131 primCovMtx(1,2) = primCovMtx(2,1) = vrtErr1[4]+vrtErr2[4];
132 primCovMtx(2,2) = vrtErr1[5]+vrtErr2[5];
137 distx*wgtMtx(0,0)*distx
138 +disty*wgtMtx(1,1)*disty
139 +distz*wgtMtx(2,2)*distz
140 +2.*distx*wgtMtx(0,1)*disty
141 +2.*distx*wgtMtx(0,2)*distz
142 +2.*disty*wgtMtx(1,2)*distz;
143 signif=std::sqrt(std::abs(signif));
144 if(signif != signif) signif = 0.;
150 double dx = Vrt1.x()- Vrt2.x();
151 double dy = Vrt1.y()- Vrt2.y();
152 double dz = Vrt1.z()- Vrt2.z();
153 return std::sqrt(
dx*
dx+
dy*
dy*dz*dz);
161 double DirX=SecVrt.x(), DirY=SecVrt.y();
162 double Covar = DirX*VrtErr[0]*DirX
163 +2.*DirX*VrtErr[1]*DirY
164 +DirY*VrtErr[2]*DirY;
165 Covar /= DirX*DirX + DirY*DirY;
166 Covar=std::sqrt(std::abs(Covar));
167 if(Covar != Covar) Covar = 0.;
177 double massP,
double massPi )
180 double ap1i=std::abs(TrkAtVrt[0][2]);
double ap2i=std::abs(TrkAtVrt[1][2]);
185 double px = phi1.
cs * theta1.
sn * ap1i
186 + phi2.
cs * theta2.
sn * ap2i;
187 double py = phi1.
sn * theta1.
sn * ap1i
188 + phi2.
sn * theta2.
sn * ap2i;
189 double pz = theta1.
cs * ap1i
191 double ee= (ap1i > ap2i) ?
192 (std::sqrt(ap1i*ap1i+massP*massP)+std::sqrt(ap2i*ap2i+massPi*massPi)):
193 (std::sqrt(ap2i*ap2i+massP*massP)+std::sqrt(ap1i*ap1i+massPi*massPi));
195 return test>0 ? std::sqrt(
test) : 0.;
203 double api=1./std::abs(inpTrk[2]);
206 double px = phi.cs * theta.sn * api;
207 double py = phi.sn * theta.sn * api;
208 double pz = theta.cs * api;
219 if( IBLexp==0 )
return -1;
228 if( BLexp==0 )
return -1;
245 truthParticleLinkAcc (
"truthParticleLink");
248 truthParticleLinkAcc(*TP);
249 if( !tplink.
isValid() )
return 0;
251 truthMatchProbabilityAcc (
"truthMatchProbability" );
252 if( truthMatchProbabilityAcc( *TP ) < 0.75 )
return 0;
254 if( (*tplink)->hasProdVtx()){
255 if( (*tplink)->prodVtx()->nIncomingParticles()==1){
256 int PDGID1=0, PDGID2=0, PDGID3=0;
264 if((*tplink)->prodVtx()->perp()>1.)
return 1.;
265 if(noBC1 && noBC2 && noBC3)
return 0;
273 if(PDGID<=0)
return 1;
274 if(PDGID>600 && PDGID<4000)noBC=1;
275 if(PDGID<400 || PDGID>5600)noBC=1;
276 if(PDGID==513 || PDGID==523 || PDGID==533 || PDGID==543)noBC=1;
277 if(PDGID==5114 || PDGID==5214 || PDGID==5224 || PDGID==5314 || PDGID==5324)noBC=1;
296 truthParticleLinkAcc (
"truthParticleLink");
299 truthParticleLinkAcc( *TP );
306 truthParticleLinkAcc (
"truthParticleLink");
309 truthParticleLinkAcc( *TP );
310 if( !tplink.
isValid() )
return 1;
317 const EventContext& ctx = Gaudi::Hive::currentContext();
318 if(Vrt.
fitVertex.perp()<20.)
return 1.e9;
327 someLayer = volume->associatedLayer(Vrt.
fitVertex);
330 nextLayerP=someLayer;
343 if(nextLayerP){ extrapParP =
m_extrapolator->extrapolate(ctx, pseudoVrtPart,
345 if(nextLayerN){ extrapParN =
m_extrapolator->extrapolate(ctx, pseudoVrtPart,
348 float distanceP=1.e9, distanceN=1.e9;
351 if(distanceP==1.e9 && distanceN==1.e9)
return 1.e9;
357 std::vector<double> pntCovar={1.e-2,0.,1.e-2,0.,0.,4.e-2};
367 std::vector<double> estimation(3,0.);
368 int ntsel=selTrk.size();
369 for(
int i=0;
i<ntsel-1;
i++){
370 for(
int j=
i+1; j<ntsel; j++){
371 int k = selTrk[
i]<selTrk[j] ? selTrk[
i]*nTrk+selTrk[j] : selTrk[j]*nTrk+selTrk[
i];
372 estimation[0]+=vrt.at(
k)[0];
373 estimation[1]+=vrt[
k][1];
374 estimation[2]+=vrt[
k][2];
376 estimation[0] /= ntsel*(ntsel-1)/2;
377 estimation[1] /= ntsel*(ntsel-1)/2;
378 estimation[2] /= ntsel*(ntsel-1)/2;