77 void fillVertex(
VKVertex *vk,
int NTRK,
const long int *ich,
78 double xyz0[3],
double (*par0)[3],
double (*inp_Trk5)[5],
79 double (*inp_CovTrk5)[15]) {
83 double xyz[3] = {0., 0., 0.};
90 for (
int tk = 0; tk < NTRK; tk++) {
92 vk->
TrackList[tk] = std::make_unique<VKTrack>(
93 TrkID, &inp_Trk5[tk][0], &inp_CovTrk5[tk][0], vk, vrtForCFT.
wm[tk]);
96 vk->
tmpArr[tk]= std::make_unique< TWRK > ();
107 for(
int tk=0; tk<NTRK; tk++){
117 bool insideGoodVolume =
true;
124 return insideGoodVolume;
127 void protectCurvatureSign(
128 double ini,
double cur,
130 double x =
cur / ini;
135 Wgt[14] *= (0.7 * 0.7) / (
x *
x);
140 int i, jerr, tk,
it=0;
143 double chi2min, chi21s=11., chi22s=10., vShift;
144 double aermd[6],tmpd[6]={0.};
145 double tmpPer[5],tmpCov[15], tmpWgt[15];
146 double VrtMomCov[21],PartMom[4];
147 double cnstRemnants=0., iniCnstRem=1.e-12;
150 const double ConstraintAccuracy=1.e-6;
163 double zeroV[3]={0.,0.,0.};
169 std::vector<int>
index;
173 for(tk=0; tk<NTRK; tk++){
if( vrtForCFT.
indtrkmc[
ic][tk] )
index.push_back(tk); }
179 vk->
ConstraintList.emplace_back(std::make_unique<VKPointConstraint>( NTRK, vrtForCFT.
vrt, vk,
false));
182 vk->
ConstraintList.emplace_back(std::make_unique<VKPointConstraint>( NTRK, vrtForCFT.
vrt, vk,
true));
186 IERR=
cfdinv(tmpd, aermd, -3);
if (IERR) { IERR = -4;
return IERR; }
194 vk->
ConstraintList.emplace_back(std::make_unique<VKPlaneConstraint>( NTRK, vrtForCFT.
Ap, vrtForCFT.
Bp, vrtForCFT.
Cp, vrtForCFT.
Dp, vk));
220 bool extrapolationDone=
false;
221 bool forcedExtrapolation=
false;
223 double cnstRemnantsPrev=1.e20, Chi2Prev=0.;
int countBadStep=0;
227 savedExtrapVertices[
it-1].Set(newVrtXYZ);
232 extrapolationDone=
false;
233 bool insideGoodVolume=checkPosition(vk,newVrtXYZ);
236 if( forcedExtrapolation || vk->
truncatedStep || !insideGoodVolume){
237 insideGoodVolume=
false;
238 while( !insideGoodVolume ){
239 newVrtXYZ[0]=(2.*savedExtrapVertices[
it-1].X + savedExtrapVertices[
it-2].X)/3.;
240 newVrtXYZ[1]=(2.*savedExtrapVertices[
it-1].Y + savedExtrapVertices[
it-2].Y)/3.;
241 newVrtXYZ[2]=(2.*savedExtrapVertices[
it-1].Z + savedExtrapVertices[
it-2].Z)/3.;
242 savedExtrapVertices[
it-1].Set(newVrtXYZ);
243 insideGoodVolume=checkPosition(vk,newVrtXYZ);
244 if(!insideGoodVolume && savedExtrapVertices[
it-1].Dist3D(savedExtrapVertices[
it-2])<5.) {
return -11; }
247 extrapolationDone=
true;
248 forcedExtrapolation=
false;
249 double targV[3]={newVrtXYZ[0],newVrtXYZ[1],newVrtXYZ[2]};
250 for (tk = 0; tk < NTRK; ++tk) {
253 if(std::isnan(tmpCov[14])) {
return -7;}
254 if (std::abs(tmpCov[14])<1.
e-20
260 if(eig5>0 && eig5<1.
e-15 ){
261 tmpCov[0]+=1.e-15; tmpCov[2]+=1.e-15; tmpCov[5]+=1.e-15; tmpCov[9]+=1.e-15; tmpCov[14]+=1.e-15;
262 }
else if(tmpCov[0]>1.e9 || eig5<0.) {
267 tmpCov[1]=0.;tmpCov[3]=0.;tmpCov[6]=0.;tmpCov[10]=0.;
268 tmpCov[4]=0.;tmpCov[7]=0.;tmpCov[11]=0.;
269 tmpCov[8]=0.;tmpCov[12]=0.;
271 tmpCov[0]=std::abs(tmpCov[0]); tmpCov[2]=std::abs(tmpCov[2]);tmpCov[5]=std::abs(tmpCov[5]);
272 tmpCov[9]=std::abs(tmpCov[9]); tmpCov[14]=std::abs(tmpCov[14]);
275 if(tmpCov[0]>1.e9){ IERR=-7;
return IERR; }
276 jerr=
cfInv5(tmpCov, tmpWgt);
277 if(jerr)jerr=
cfdinv(tmpCov, tmpWgt,5);
278 if(jerr){ IERR=-7;
return IERR; }
279 vk->
TrackList[tk]->setCurrent(tmpPer,tmpWgt);
288 for( tk=0; tk<NTRK; tk++) vk->
TrackList[tk]->restoreCurrentWgt();
295 for (tk = 0; tk < NTRK; ++tk){
301 for( tk=0; tk<NTRK; tk++){
311 if (IERR)
return IERR;
314 if(extrapolationDone && chi21s*10.<vk->
Chi2 &&
it>2 && vk->
Chi2>NTRK*5.){
315 double ddstep=savedExtrapVertices[
it-1].Dist3D(savedExtrapVertices[
it-2]);
318 forcedExtrapolation=
true;
323 chi22s = chi21s * 1.01 + 10.;
332 if(jerr!=0)
return -17;
333 cfdcopy( PartMom, &dparst[3], 3);
341 for( tk=0; tk<NTRK; tk++){
351 if ( IERR )
return IERR;
357 if (chi22s > 1e8) {
return -1; }
358 if (chi22s != chi22s) {
return -14; }
365 double chi2df = std::abs(chi21s - chi22s);
371 if((chi2df < PrecLimit) && (vShift < 0.001) && (cnstRemnants<ConstraintAccuracy)){
374 forcedExtrapolation=
true;
381 if ((chi2min*100. < chi22s) && (chi22s>
std::max( (2*NTRK-3)*10., 100.)) && (
it>5)){
383 IERR = -5;
return IERR;
391 cfdcopy( PartMom, &dparst[3], 3);
406 if(
it==1){
if(cnstRemnants>0.){iniCnstRem=cnstRemnants;}
else{iniCnstRem=1.e-12;} }
408 if(cnstRemnants > 0.1 && cnstRemnants/iniCnstRem > 0.1 ){
409 IERR = -6;
return IERR;
415 if(
it>5 && ((cnstRemnants>cnstRemnantsPrev)||cnstRemnants==0.) && chi22s>(Chi2Prev*1.0001)){
416 if(++countBadStep>5)
break;
417 }
else{countBadStep=0;}
418 Chi2Prev=chi22s; cnstRemnantsPrev=cnstRemnants;
424 if( cnstRemnants > ConstraintAccuracy ){
425 IERR = -6;
return IERR;
436 long int *ich,
double xyz0[3],
double (*par0)[3],
437 double (*inp_Trk5)[5],
double (*inp_CovTrk5)[15],
438 double xyzfit[3],
double (*parfs)[3],
double ptot[4],
439 double covf[21],
double &
chi2,
double *chi2tr)
442 std::fill_n(ptot,4,0.);
443 std::fill_n(covf,21,0.);
450 std::unique_ptr<VKVertex> MainVRT(
new VKVertex(*FitCONTROL));
451 fillVertex(MainVRT.get(), NTRK, ich, xyz0, par0, inp_Trk5, inp_CovTrk5);
455 int IERR = fitVertex( MainVRT.get());
456 if(IERR)
return IERR;
463 for (
int tk = 0; tk < NTRK; tk++) {
466 chi2tr[tk] = trk->
Chi2;
469 ptot[0]+=pp[0]; ptot[1]+=pp[1]; ptot[2]+=pp[2]; ptot[3]+=pp[3];
472 double vM2=(ptot[3]-ptot[2])*(ptot[3]+ptot[2]) - ptot[1]*ptot[1] - ptot[0]*ptot[0];
479 double dptot[5], VrtMomCov[21];
485 for(
int ti=0; ti<3+3*NTRK; ti++){
486 for(
int tj=0; tj<=ti; tj++){
490 int activeTrk[
vkalNTrkM]={0}; std::fill_n(activeTrk,NTRK,1);
491 double vrtMass=0., vrtMassError=0.;