FairRoot/PandaRoot
RhoKinHyperonVtxFitter.cxx
Go to the documentation of this file.
1 
2 #include <iostream>
3 #include "RhoKinHyperonVtxFitter.h"
4 #include "RhoBase/RhoCandListIterator.h"
5 #include "RhoBase/RhoFactory.h"
6 #include "RhoCalculationTools.h"
7 #include "TDecompLU.h"
8 #include "TMatrixD.h"
9 #include "TMatrixDSym.h"
10 #include "RhoVtxPoca.h"
11 
12 using namespace std;
13 
14 
15 ClassImp(RhoKinHyperonVtxFitter)
16 
17 TBuffer& operator>>(TBuffer& buf, RhoKinHyperonVtxFitter *&obj)
18 {
19  obj = (RhoKinHyperonVtxFitter*) buf.ReadObject(RhoKinHyperonVtxFitter::Class());
20  return buf;
21 }
22 
23 //Include only those constraint which need vertex Info....
24 RhoKinHyperonVtxFitter::RhoKinHyperonVtxFitter( RhoCandidate* b) :
25  RhoFitterBase( b )
26 {
27  fMassConstraint =-1;
28  //fPointConstraint=-1;
29 
30  fMinDChisq=0.01;
31  fNMaxIterations=20;
32  fIterateExact=false;
33 
34 
35 }
36 
37 RhoKinHyperonVtxFitter::~RhoKinHyperonVtxFitter()
38 {
39 }
40 
41 void RhoKinHyperonVtxFitter::AddMassConstraint(double mass)
42 {
43  fMassConstraint = 1;
44  fMass=mass;
45 }
46 
47 // void RhoKinHyperonVtxFitter::AddPointingConstraint(TVector3 pVtx)
48 // {
49 // fPointConstraint = 1;
50 // fpVtx=pVtx;
51 // }
52 //
53 
54 Bool_t RhoKinHyperonVtxFitter::FitNode(RhoCandidate *cand)
55 {
56  SetDaugthersFromComposite(cand);
57  Bool_t check=Compute(cand);
58  SetOutput(cand);
59  return check;
60 }
61 
62 
63 void RhoKinHyperonVtxFitter::SetMatrices()
64 {
65  int nd=fDaughters.size();
66 
67  fNvar=7;
68  fNpart=nd;
69  fNpar =nd*fNvar;
70  fNcon=NumCon;
71  // if(fVerbose) cout << fNcon << "Num " << endl;
72  fNc=0;
73  fNiter=0;
74 
75  al0.ResizeTo(7*nd,1);
76  V_al0.ResizeTo(fNpar,fNpar);
77  al1.ResizeTo(7*nd,1);
78  V_al1.ResizeTo(fNpar,fNpar);
79  vtx_ex.ResizeTo(3,1);
80  vtx_st.ResizeTo(3,1);
81  mPull.ResizeTo(7*nd,1);
82  covC.ResizeTo(7,7);
83 }
84 
85 
86 void RhoKinHyperonVtxFitter::ResetMatrices()
87 {
88  al0.Zero();
89  V_al0.Zero();
90  al1.Zero();
91  V_al1.Zero();
92  mPull.Zero();
93  vtx_ex.Zero();
94  vtx_st.Zero();
95  covC.Zero();
96 }
97 
98 
99 
100 Bool_t RhoKinHyperonVtxFitter::Compute(RhoCandidate* c)
101 {
102 
103  // int nd=fDaughters.size();
104  int nd=fDaughters.size();
105 //cout<<"nd="<<nd<<endl;
106  NumCon=2*nd;
107 
108  if(fMassConstraint >0) {
109  NumCon += 1;
110  }
111 // if(fPointConstraint >0) {
112 // NumCon += 2;
113 // }
114 
115  fNDegreesOfFreedom=NumCon-3;
116 
117 
118 
119  SetMatrices();
120  ResetMatrices();
121  ReadMatrix();
122 
123  TMatrixD cov_al_x(7*nd,3);
124 
125  TVector3 startVtx;
126  //Getting point of closed approach as start vertex point
127  RhoVtxPoca poca; //changed from internal method to class RhoVtxPoca by J.Puetz
128  poca.GetPocaVtx(startVtx, c);
129 
130  vtx_st[0][0]=startVtx.X();
131  vtx_st[1][0]=startVtx.Y();
132  vtx_st[2][0]=startVtx.Z();
133 
134  vtx_ex=vtx_st;
135  if(fVerbose) { cout<<"Initial vertex Position is "<<vtx_ex[0][0]<<" "<<vtx_ex[1][0]<<" "<<vtx_ex[2][0]<<endl; }
136 
137  al1=al0;
138  V_al1=V_al0;
139  TransportToPoca(al0,V_al0,al1,V_al1,vtx_ex);
140  al0=al1;
141  V_al0=V_al1;
142 
143 
144 
145 
146  TMatrixD V_vtx(3,3);
147  V_vtx[0][0] = 9000.;
148  V_vtx[1][1] = 9000.;
149  V_vtx[2][2] = 9000.;
150  // vtx_ex=vtx_st;
151 
152  double ierr =0 ; // used to check inversions
153  TMatrixD chi2(1,1);
154  TMatrixD chi2_1(1,1);
155 
156  chi2[0][0]=2000000.;
157  //double tmp_chiSq = 999;
158 
159  for(Int_t j1=0; j1<fNMaxIterations; ++j1) {
160 
161  fNc=0;
162  if(fMassConstraint >0) {
163  ReadMassKinMatrix();
164  } else {
165  ReadKinMatrix();
166  }
167 
168  TMatrixD mD_t=mD;
169  mD_t.T();
170  // mD_t=mD_t.Transpose(mD);
171 
172 
173  TMatrixD Vd_inv = mD*V_al0*mD_t;
174  if(Vd_inv==0) { continue; }
175  TMatrixD Vd = Vd_inv.Invert(&ierr);
176 
177 
178  // if( ierr != 0 ){
179  // if(fVerbose) cout << "Inversion of constraint-matrix failed! " << endl;
180  // return 0;}
181  // Vd.Print();
182 
183 
184  TMatrixD del_al = al0 - al1;
185  TMatrixD del_x0 = vtx_st - vtx_ex;
186 
187 
188  // Lagrange multiplier
189 
190  //TMatrixD lam0=Vd*md;
191  TMatrixD lam0 = Vd* ( mD*del_al + md);
192  // if(fVerbose) cout << " lam0 calculated" << endl;
193 
194 
195  // Position Derivative matrix ...............
196  TMatrixD mE_t=mE;
197  mE_t.T();
198  //TMatrixD Vx_inv = mE_t*Vd*mE;
199  TMatrixD Vx0_inv=V_vtx.Invert(&ierr);
200  TMatrixD Vx_inv = Vx0_inv + mE_t*Vd*mE;
201  TMatrixD Vx=Vx_inv.Invert(&ierr);
202 // TMatrixD del_x=Vx*Vx0_inv*del_x0- Vx*mE_t*lam0;// by x.song
203 // cout<<" *** *** *** *** ***"<<endl;
204 // cout<<"mE "; mE.Print();
205 // cout<<"mE_t "; mE_t.Print();
206 // cout<<"Vd "; Vd.Print();
207 // cout<<"Vx_inv "; Vx_inv.Print();
208 // cout<<"Vx "; Vx.Print();
209 // cout<<" *******************"<<endl;
210 
211  // New vertex and covariance ........
212  TMatrixD V_vtx_new(3,3);
213  TMatrixD vtx_new(vtx_ex);
214 
215 
216  //vtx_new -= Vx*mE_t*lam0;
217  vtx_new += Vx*Vx0_inv*del_x0 - Vx*mE_t*lam0;
218  //cout << " New vtx calculated" << endl;
219  //vtx_new.Print();
220  V_vtx_new = Vx;
221 
222 
223 
224  // Final Lagarange multiplier ........
225  TMatrixD lam = lam0 + (Vd * mE) * (vtx_new - vtx_ex);
226  // if(fVerbose) cout << " New lam calculated" << endl;
227 
228 
229 
230  // New track parameters.............
231  TMatrixD al_new(al0);
232  al_new -= V_al0*mD_t*lam;
233  // if(fVerbose) cout << " New track param calculated" << endl;
234 
235 
236 
237  //chiSquared
238  TMatrixD lam_t=lam;
239  lam_t.T();
240  TMatrixD chi2_new = lam_t*(mD*(al0 - al_new) + mE*(vtx_st-vtx_ex) + md);
241 
242 
243 
244  // New Covariance Matrix................
245  // TMatrixD V_al_new(V_al0);
246  // V_al_new-=V_al0*mD_t*Vd*mD*V_al0_t;
247 
248 
249  // protect against errors. RK: is that safe to do?
250  if(TMath::IsNaN(chi2_new[0][0])) continue;
251 
252  double deltaChi=chi2_new[0][0]-chi2[0][0];
253 
254  // Check chi^2. If better yes update the values ..............................
255  if (deltaChi>0.1*chi2[0][0]) {continue;}
256  if( chi2_new[0][0] < chi2[0][0] ) {
257  vtx_ex = vtx_new;
258  al1 = al_new;
259  // V_al0 = V_al_new;
260  }
261 
262 
263  // if (j1==0) {chi2=chi2_new;continue;}
264 
265  // If Chi^2 change is small then go out of iteration......................
266  if( (j1+1 == fNMaxIterations) ||
267  (!fIterateExact && ( fabs(chi2[0][0]-chi2_new[0][0])<fMinDChisq) ) ) {
268  chi2 = chi2_new;
269  vtx_ex = vtx_new;
270  al0 =al_new;
271 
272 
273  TMatrixD Vd_new(Vd);
274  Vd_new -= Vd*(mE*Vx*mE_t)*Vd.T();
275  TMatrixD V_al_new(V_al0);
276  V_al_new -= V_al0*(mD_t*Vd_new*mD)*V_al0.T();
277  cov_al_x -= V_al0*mD_t*Vd*mE*Vx; // Vertex-track Correlation
278 
279 
280 
281  //double covdif=(V_al0[6][6]-V_al_new[6][6]);
282  // if (covdif > 0 ) {mPull[0][0] =(al0[6][0]-al_new[6][0])/sqrt(covdif);}
283  mPull[0][0] =(al0[6][0]-al_new[6][0]);
284  fPull=mPull[0][0];
285 
286  V_al0 = V_al_new;
287  V_vtx = V_vtx_new;
288  if(fVerbose) { cout <<"iteration Number " << " " << j1 <<" final." <<endl; }
289  if(fVerbose) { cout <<" chi2 in iteration" << " " << chi2[0][0] << " pull="<<fPull<< endl; }
290  break; // that was the final iteration, stop the loop
291  }
292  chi2 = chi2_new;
293  if(fVerbose) { cout << "iteration Number " << " " << j1 << endl; }
294  if(fVerbose) { cout <<" vertex Position is "<<vtx_new[0][0]<<" "<<vtx_new[1][0]<<" "<<vtx_new[2][0]<<endl; }
295  if(fVerbose) { cout << " chi2 in iteration" << " " << chi2[0][0] << endl; }
296  } // end of iteration-loop
297 
298  // tell that the fit failed if we have no updated chi2
299  if( TMath::IsNaN(chi2[0][0]) || chi2[0][0]==2000000. ) return kFALSE;
300 
301  TMatrixD al_new_vtx(7*nd,1);
302  TMatrixD Va_new_vtx(7*nd,7*nd);
303 
304 
305 
306 
307  // Trivial way for covariance matrix of composite
308 // covC=
309 // for(Int_t k=0;k<nd;k++) {
310 // for(int j = 0; j< nd ; ++j) {
311 // {if(k<=j) CovS(k,j) = V_al_0(k,j);}}
312 // CovC -= CovS;
313 // covC+=V_al0.GetSub(k*7,(k+1)*7-1,k*7,(k+1)*7-1);}
314 
315  GetCovariance(V_al0,cov_al_x,V_vtx,covC);
316 
317 
318  // al_new_vtx=al1;
319  // Va_new_vtx=V_al1;
320 
321 //transport the daughters para and cov to those on the fitted vertex
322 
323  TransportToVertex(al0, V_al0, al_new_vtx, Va_new_vtx, vtx_ex);
324  al0=al_new_vtx;
325  V_al0=Va_new_vtx;
326  fChiSquare=chi2[0][0];
327  // fChi2Diff=chi2_1[0][0]-chi2[0][0];
328 
329  return kTRUE;
330 }
331 
332 
333 
334 
335 
336 
337 //Write output
338 void RhoKinHyperonVtxFitter::SetOutput(RhoCandidate* head)
339 {
340 
341  int nd=fDaughters.size();
342  TMatrixD m(nd,1);
343 
344 
345  double sumA=0;
346  double a;
347  for (int k=0; k<nd; k++) {
348 
349 
350  //skip locked daughters
351  if(fDaughters[k]->IsLocked()) continue;
352  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
353  a = -0.00299792458*bField*fDaughters[k]->GetCharge();
354  sumA += a;
355  TVector3 pos(al0[k*7+4][0],al0[k*7+5][0],al0[k*7+6][0]);
356  TLorentzVector mom4(al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0], al0[k*7+3][0]);
357 // fDaughters[k]->SetP7(pos,mom4);
358 //better to put daugthers with mass hypothesis .......?? VJ
359  TLorentzVector momM;
360  double fM=fDaughters[k]->Mass();
361  momM.SetXYZM(al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0],fM);
362 // momM.SetP4(fM,al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0]);
363  fDaughters[k]->SetP7(pos,momM);
364 
365 
366  //Extend matrix for energy for each candidates if daughters from mass hypothesis 6x6 covariance
367  TMatrixD p1Cov(7,7);
368  TLorentzVector p1=fDaughters[k]->P4();
369  TMatrixD p2Cov(7,7);
370 
371  for(int i=0; i<7; i++) {
372  for (int j=0; j<7; j++) {
373  p1Cov[i][j]= V_al0[k*7+i][k*7+j];
374  }
375  }
376 
377 
378  //Change from px,py,pz,E,x,y,z
379  // to x,y,z,px,py,pz,E
380  for(int i=0; i<7; i++) {
381  for(int j=0; j<7; j++) {
382  if(i>=4) {
383  if(j>=4) {
384  p2Cov[i-4][j-4] = p1Cov[i][j];
385  } else { p2Cov[i-4][j+3] = p1Cov[i][j]; }
386  } else {
387  if(j>=4) {
388  p2Cov[i+3][j-4] = p1Cov[i][j];
389  } else { p2Cov[i+3][j+3] = p1Cov[i][j]; }
390  }
391  }
392  }
393 
394  // create cov with E... check it
395  double invE = 1./al0[k*7+3][0];
396  p2Cov[3][6] = p2Cov[6][3] = (p1.X()*p1Cov[0][0]+p1.Y()*p1Cov[0][1]+p1.Z()*p1Cov[0][2])*invE;
397  p2Cov[4][6] = p2Cov[6][4] = (p1.X()*p1Cov[1][0]+p1.Y()*p1Cov[1][1]+p1.Z()*p1Cov[1][2])*invE;
398  p2Cov[5][6] = p2Cov[6][5] = (p1.X()*p1Cov[2][0]+p1.Y()*p1Cov[2][1]+p1.Z()*p1Cov[2][2])*invE;
399 
400  p2Cov[6][6] = (p1.X()*p1.X()*p1Cov[0][0]+p1.Y()*p1.Y()*p1Cov[1][1]+p1.Z()*p1.Z()*p1Cov[2][2]
401  +2.0*p1.X()*p1.Y()*p1Cov[0][1]
402  +2.0*p1.X()*p1.Z()*p1Cov[0][2]
403  +2.0*p1.Y()*p1.Z()*p1Cov[1][2])*invE*invE;
404 
405  p2Cov[6][0] = p2Cov[0][6] = (p1.X()*p1Cov[0][4]+p1.Y()*p1Cov[1][4]+p1.Z()*p1Cov[2][4])*invE;
406  p2Cov[6][1] = p2Cov[1][6] = (p1.X()*p1Cov[0][5]+p1.Y()*p1Cov[1][5]+p1.Z()*p1Cov[2][5])*invE;
407  p2Cov[6][2] = p2Cov[2][6] = (p1.X()*p1Cov[0][6]+p1.Y()*p1Cov[1][6]+p1.Z()*p1Cov[2][6])*invE;
408 
409  fDaughters[k]->SetCov7(p2Cov); //New covariance matrix with correlations
410  //cout<< " ####### KinVtx daughter cov check... " << endl;
411  //cout<<"p1Cov"; p1Cov.Print();
412  //cout<<"p2Cov"; p2Cov.Print();
413 
414  }
415 
416 
418 //
419 // // For the composite particle ..............................
421  double fpx=0,fpy=0,fpz=0,fe=0;
422  for (int k=0; k<nd; k++) {
423  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
424  a = -0.00299792458*bField*fDaughters[k]->GetCharge();
425 
426  fpx+= al0[k*7+0][0]-a*(vtx_ex[1][0] - al0[k*7+5][0]);
427  fpy+= al0[k*7+1][0]+a*(vtx_ex[0][0] - al0[k*7+4][0]);
428  fpz+= al0[k*7+2][0];
429  fe += al0[k*7+3][0];
430  }
431 
432  TLorentzVector sum(fpx,fpy,fpz,fe);
433  TVector3 vtx(vtx_ex[0][0],vtx_ex[1][0],vtx_ex[2][0]);
434  head->SetP7(vtx,sum); //[ralfk:01.12.11 Try to make it a leaf-by-leaf fit]
435  //We set the decay vertex of the mother! [R.K.]
436 
437  TMatrixD CovV = covC.GetSub(0,2,0,2);
438 
439  head->SetPos(vtx);//P4 is defined here
440  SetDecayVertex(head,vtx,CovV);
441 // SetFourMomentumByDaughters(head);//propagates cov7 from daughters
442  head->SetCov7(covC);//which one to use???
443  //cout<<" KinVtx Cov7: ";covC.Print();
444  if(fVerbose) { cout<<"Final vertex Position is "<<vtx_ex[0][0]<<" "<<vtx_ex[1][0]<<" "<<vtx_ex[2][0]<<endl; }
445  if(fVerbose) { cout<<"Final Momenta are "<<al0[0][0]<<" "<<al1[1][0]<<" "<<al1[2][0]<<endl; }
446 }
447 
448 
449 
450 //Read the input vector
451 void RhoKinHyperonVtxFitter::ReadMatrix()
452 {
453  int nd =fDaughters.size();
454  TMatrixD m(nd,1);
455  for (int k=0; k<nd; k++) {
456  int kN=k*7;
457  //px,py,pz,E,x,y,z
458  TLorentzVector p1=fDaughters[k]->P4(); //4-momentum of the daughter
459  TVector3 p2=fDaughters[k]->Pos(); // position of the daughter
460  al0[kN+0][0]=p1.X();
461  al0[kN+1][0]=p1.Y();
462  al0[kN+2][0]=p1.Z();
463  // al0[kN+3][0]=p1.E();
464  al0[kN+4][0]=p2.X();
465  al0[kN+5][0]=p2.Y();
466  al0[kN+6][0]=p2.Z();
467 
468  double fm=fDaughters[k]->Mass();
469  al0[kN+3][0]=sqrt(al0[kN+0][0]*al0[kN+0][0]+ al0[kN+1][0]*al0[kN+1][0]+al0[kN+2][0]*al0[kN+2][0]+fm*fm);
470 
471  // Read Covariance Matrix .... Can read 6x6 matrices..................
472  TMatrixD p1Cov(7,7);
473  TMatrixD p3Cov(6,6); //Why 6x6 here if 7x7 should be cpoied (below)
474  TMatrixD p2Cov(7,7);
475  TMatrixD p4Cov(7,7);
476  p1Cov=fDaughters[k]->Cov7(); //Cov Matrix x,y,z,px,py,pz,E
477 
478  for (int ii=0; ii<6; ii++) {for(int jj=0; jj<6; jj++) {p3Cov[ii][jj]=p1Cov[ii][jj];}} //test
479 
480  //Extend matrix for energy for each candidates .....6x6 to 7x7
481  TMatrixD J(7,6) ;
482  J.Zero();
483  TMatrixD J_t(6,7);
484  for (int ii=0; ii<6; ii++) {for(int jj=0; jj<6; jj++) {J[ii][jj] = 1;}}
485  for(int i=3; i<6; ++i) {J[6][i] = al0[kN+i-3][0]/al0[kN+3][0];}
486  // p2Cov= J*p3Cov*(J_t.Transpose(J));
487  p2Cov=p1Cov;
488  //Change to px,py,pz,E,x,y,z
489  for(int i=0; i<7; i++) {
490  for(int j=0; j<7; j++) {
491  if(i>=3) {
492  if(j>=3) {
493  p4Cov[i-3][j-3] = p2Cov[i][j];
494  } else { p4Cov[i-3][j+4] = p2Cov[i][j]; }
495  } else {
496  if(j>=3) {
497  p4Cov[i+4][j-3] = p2Cov[i][j];
498  } else { p4Cov[i+4][j+4] = p2Cov[i][j]; }
499  }
500  }
501  }
502 
503  for(int i=0; i<7; i++) {
504  for (int j=0; j<7; j++) {
505  V_al0[k*7+i][k*7+j] = p4Cov[i][j];
506  }
507  }
508  }
509 }
510 
511 
512 //Read Constraint Matrices ... D, E and d
513 //unsigned RhoKinHyperonVtxFitter:: ReadKinMatrix( TMatrixD & mD, TMatrixD & mE, TMatrixD & md)
514 void RhoKinHyperonVtxFitter::ReadKinMatrix()
515 {
516 
517  int nd=fDaughters.size();
518  fNc=0;
519  mD.ResizeTo(fNcon,fNpar);
520  mE.ResizeTo(fNcon,3);
521  md.ResizeTo(fNcon,1);
522  for (int k=0; k<nd; k++) {
523 
524 
525  int kN=k*7;
526  int k2=k*2;
527  double delX = vtx_ex[0][0] - al1[kN+4][0];
528  double delY = vtx_ex[1][0] - al1[kN+5][0];
529  double delZ = vtx_ex[2][0] - al1[kN+6][0];
530  double px = al1[kN+0][0];
531  double py = al1[kN+1][0];
532  double pz = al1[kN+2][0];
533  double ch=fDaughters[k]->GetCharge();
534  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
535  double a = -0.0029979246*ch*bField;
536  double pT_2 = px*px + py*py;
537 
538  double J = a*(delX*px + delY*py)/pT_2;
539  double Rx = delX - 2.*px*(delX*px + delY*py)/pT_2;
540  double Ry = delY - 2.*py*(delX*px + delY*py)/pT_2;
541 
542  // if(fabs(J) > 1) {return 0;}
543  if (J>=1.0 || J<= -1.0) { J = (J>=1.0 ? 0.99 : -0.99);}
544 
545  double S = 1./(pT_2*sqrt(1-(J*J)));
546  // if(fVerbose) cout << ch << "ch" << S << "pt2" << J << "bField" << bField <<endl;
547  double asin_J = asin(J);
548  //charged particle
549  if(ch !=0) {
550  mD[fNc+0+k2][kN+0] = delY ;
551  //if(fVerbose) cout << al0[kN+4][0] << " " << delY << endl;
552  mD[fNc+0+k2][kN+1] = -(delX) ;
553  mD[fNc+0+k2][kN+2] = 0. ;
554  mD[fNc+0+k2][kN+3] = 0. ;
555  mD[fNc+0+k2][kN+4] = py + a*delX ;
556  mD[fNc+0+k2][kN+5] = -px + a*delY ;
557  mD[fNc+0+k2][kN+6] = 0. ;
558 
559  mD[fNc+1+k2][kN+0] = -pz*S*Rx ;
560  mD[fNc+1+k2][kN+1] = -pz*S*Ry ;
561  mD[fNc+1+k2][kN+2] = -asin_J/a ;
562  mD[fNc+1+k2][kN+3] = 0.;
563  mD[fNc+1+k2][kN+4] = px*pz*S;
564  mD[fNc+1+k2][kN+5] = py*pz*S;
565  mD[fNc+1+k2][kN+6] = -1.;
566  } else {
567  //neutral particle
568  mD[fNc+0+k2][kN+0] = delY ;
569  mD[fNc+0+k2][kN+1] = -(delX) ;
570  mD[fNc+0+k2][kN+2] = 0. ;
571  mD[fNc+0+k2][kN+3] = 0. ;
572  mD[fNc+0+k2][kN+4] = py;
573  mD[fNc+0+k2][kN+5] = -px;
574  mD[fNc+0+k2][kN+6] = 0. ;
575 
576  mD[fNc+1+k2][kN+0] = 2*(delX*px+delY*py)*px*pz/(pT_2*pT_2) - pz*delX/(pT_2);
577  mD[fNc+1+k2][kN+1] = 2*(delX*px+delY*py)*py*pz/(pT_2*pT_2) - pz*delY/(pT_2);
578  mD[fNc+1+k2][kN+2] =-(delX*px+delY*py)/(pT_2);
579  mD[fNc+1+k2][kN+3] = 0.;
580  mD[fNc+1+k2][kN+4] = px*pz/pT_2;
581  mD[fNc+1+k2][kN+5] = py*pz/pT_2;
582  mD[fNc+1+k2][kN+6] = -1.;
583  }
584  //DMat_trk.push_back(DMat_tmp);
585  //E jacobian matrix
586  if(ch !=0 ) {
587  mE[fNc+0+k2][0] = -(py + a*delX);
588  mE[fNc+0+k2][1] = (px - a*delY);
589  mE[fNc+0+k2][2] = 0.;
590  mE[fNc+1+k2][0] = -px*pz*S;
591  mE[fNc+1+k2][1] = -py*pz*S;
592  mE[fNc+1+k2][2] = 1.;
593  } else {
594  mE[fNc+0+k2][0] = -py ;
595  mE[fNc+0+k2][1] = px;
596  mE[fNc+0+k2][2] = 0.;
597  mE[fNc+1+k2][0] = -px*pz/pT_2;
598  mE[fNc+1+k2][1] = -py*pz/pT_2;
599  mE[fNc+1+k2][2] = 1.;
600  }
601  if(ch !=0 ) {
602  md[fNc+0+k2][0] = delY*px - delX*py - (a/2.)*(delX*delX + delY*delY);
603  md[fNc+1+k2][0] = delZ - (pz/a)*asin_J;
604  } else {
605  md[fNc+0+k2][0] = delY*px - delX*py;
606  md[fNc+1+k2][0] = delZ - pz*(delX * px + delY * py)/(pT_2);
607  }
608 
609  }
610  fNc +=2*nd;
611 
612 
613 
614 
615 }
616 
617 
618 void RhoKinHyperonVtxFitter::ReadMassKinMatrix()
619 // unsigned RhoKinHyperonVtxFitter:: ReadMassKinMatrix( TMatrixD & mD, TMatrixD & mE, TMatrixD & md)
620 {
621  // if(m_fitIncludingVertex == 0)
622  //{
623  int nd=fDaughters.size();
624 
625  mD.ResizeTo(fNcon,fNpar);
626  mE.ResizeTo(fNcon,3);
627  md.ResizeTo(fNcon,1);
628 
629  double Etot = 0.;
630  double Px = 0.;
631  double Py = 0.;
632  double Pz = 0.;
633 
634  TMatrixD al1p(al1);
635  double a=0; //TODO: is that right to initialize with 0?
636  TMatrixD m(nd,1);
637 
638  /*
639  // Mass Constraint without vertex Info.............................
640 
641  for(unsigned k=0;k<nd;++k){
642  int kN=k*7;
643  TLorentzVector p1=fDaughters[k]->P4();
644  m[k][0]=p1.M();
645  double px = al1p[kN+0][0];
646  double py = al1p[kN+1][0];
647  double pz = al1p[kN+2][0];
648  double E = TMath::Sqrt(px*px+py*py+pz*pz+m[k][0]*m[k][0]);
649 
650  // Etot += E; Px +=px; Py +=py; Pz += pz;}
651  // md[fNc+0][0] = Etot*Etot - Px*Px - Py*Py - Pz*Pz - fMass*fMass ;
652 
653  a = -0.00299792458*2.0*fDaughters[k]->GetCharge();
654  Double_t invE = 1./E;
655 
656  mD[fNc+0][kN+0] = 2.*(Etot*px*invE-Px);
657  mD[fNc+0][kN+1] = 2.*(Etot*py*invE-Py);
658  mD[fNc+0][kN+2] = 2.*(Etot*pz*invE-Pz);
659  // mD[fNc+0][kN+3] = 0.0;
660  mD[fNc+0][kN+3] = 2* m[k][0]*Etot*invE;
661  mD[fNc+0][kN+4] = 2.*Py*a;
662  mD[fNc+0][kN+5] = -2.*Px*a;
663  mD[fNc+0][kN+6] = 0.0;
664 
665  //................Simple....................
666 
667  // mD[fNc+0][kN+0] = -2.*Px;
668  // mD[fNc+0][kN+1] = -2.*Py;;
669  // mD[fNc+0][kN+2] = -2.*Pz;
670  // mD[fNc+0][kN+3] = 2.*Etot;
671  // mD[fNc+0][kN+3] = 2* m[k][0]*Etot*invE;
672  // mD[fNc+0][kN+4] = 2.*(Etot*py*invE-Py)*a;
673  // mD[fNc+0][kN+5] = 2.*(Etot*px*invE-Px)*a;
674  // mD[fNc+0][kN+4] = 0.0;
675  // mD[fNc+0][kN+5] = 0.0;
676  // mD[fNc+0][kN+6] = 0.0;
677  }
678  */
679 
680  // With Vertex Info.............................................
681  for(int k=0; k<nd; ++k) {
682  int kN=k*7;
683  TLorentzVector p1=fDaughters[k]->P4();
684  m[k][0]=p1.M();
685  double delX = vtx_ex[0][0] - al1p[kN+4][0];
686  double delY = vtx_ex[1][0] - al1p[kN+5][0];
687  // double delX=0.;
688  // double delY=0.;
689 
690  al1p[kN+0][0] = al1p[kN+0][0]-a*delY;
691  al1p[kN+1][0] = al1p[kN+1][0]+a*delX;
692  // al1p[kN+2][0] = al1p[kN+2][0];
693  double E = TMath::Sqrt(al1p[kN+0][0]* al1p[kN+0][0]+al1p[kN+1][0]*al1p[kN+1][0]+al1p[kN+2][0]*al1p[kN+2][0]+m[k][0]*m[k][0]);
694  Etot += E;
695 
696  Px += al1p[kN+0][0] ;
697  Py += al1p[kN+1][0];
698  Pz += al1p[kN+2][0];
699  }
700  md[fNc+0][0] = Etot*Etot - Px*Px - Py*Py - Pz*Pz - fMass*fMass ;
701 
702  double sumA=0;
703  for(int k=0; k<nd; ++k) {
704  int kN=k*7;
705  double px = al1p[kN+0][0];
706  double py = al1p[kN+1][0];
707  double pz = al1p[kN+2][0];
708  // double E = al1p[kN+3][0];
709  double E = TMath::Sqrt(px*px+py*py+pz*pz+m[k][0]*m[k][0]);
710 
711  //here there should be implemented the algorithm for neutral particles
712 
713  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
714  a = -0.00299792458*bField*fDaughters[k]->GetCharge();//TODO BField
715  sumA += a;
716  Double_t invE = 1./E;
717 
718  mD[fNc+0][kN+0] = 2.*(Etot*al1p[kN+0][0]*invE-Px);
719  mD[fNc+0][kN+1] = 2.*(Etot*al1p[kN+1][0]*invE-Py);
720  mD[fNc+0][kN+2] = 2.*(Etot*al1p[kN+2][0]*invE-Pz);
721  mD[fNc+0][kN+3] = 0.;
722  mD[fNc+0][kN+4] =-2.*(Etot*al1p[kN+1][0]*invE-Py)*a;
723  mD[fNc+0][kN+5] = 2.*(Etot*al1p[kN+0][0]*invE-Px)*a;
724  mD[fNc+0][kN+6] = 0.;
725  }
726  cout << "Sum A: " << sumA << endl;
727 
728 
729  mE[fNc+0][0] = 2*sumA*Px;
730  mE[fNc+0][1] = -2*sumA*Py;
731  mE[fNc+0][2] = 0.;
732  fNc +=1;
733 }
734 
735 // not used yet?
736 // void RhoKinHyperonVtxFitter::ReadPointingKinMatrix(RhoCandidate* head)
737 // {
738 // //Pass the vertex point
739 // // To be applied on the composite particle
740 //
741 // //int nd=fDaughters.size(); //unused
742 //
743 // mD.ResizeTo(fNcon,fNpar);
744 // mE.ResizeTo(fNcon,3);
745 // md.ResizeTo(fNcon,1);
746 //
747 // //double ch= fHeadOfTree->GetCharge(); //unused
748 // //TLorentzVector p1=fHeadOfTree->P4();
749 // //TVector3 p2=fHeadOfTree->Pos();
750 // TLorentzVector p1=head->P4(); //[ralfk:01.12.11 Try to make it a leaf-by-leaf fit]
751 // TVector3 p2=head->Pos(); //[ralfk:01.12.11 Try to make it a leaf-by-leaf fit]
752 // double delX = p2.X() - fpVtx.X();
753 // double delY = p2.Y() - fpVtx.Y();
754 // double delZ = p2.Z() - fpVtx.Z();
755 //
756 // double px = p1.X();
757 // double py = p1.Y();
758 // double pz = p1.Z();
759 //
760 //
761 // //double bField = TRho::Instance()->GetMagnetField(); //unused, why?
762 // //double a = -0.0029979246*ch*2.0; //TODO: bfield put manually here?
763 // double pT = sqrt(px*px + py*py);
764 // double delT= sqrt(delX*delX + delY*delY);
765 // double p = sqrt(px*px + py*py + pz*pz);
766 // double T = sqrt(delX*delX + delY*delY + delZ*delZ);
767 //
768 //
769 // md[fNc+0][0] = (1-delX/delT)/(delY/delT) - (1-px/pT)/(py/pT);
770 // md[fNc+1][0]= (1-delT/T)/(delZ/T) - (1-(pT/p))/(pz/p);
771 //
772 // mD[fNc+0][0] = -(px/(py*pT) - 1/py);
773 // mD[fNc+0][1] = -(1/pT - pT/(py*py)+ px/(py*py));
774 // mD[fNc+0][2] = 0;
775 // mD[fNc+0][3] = 0;
776 // mD[fNc+0][4] = delX/(delY*delT) - 1/delY;
777 // mD[fNc+0][5] = 1/delT - delT/delY*delY+ delX/(delY*delY);
778 // mD[fNc+0][6] = 0;
779 //
780 // //half angle solution
781 // mD[fNc+1][0] = -(px/pz)*(1/p - 1/pT);
782 // mD[fNc+1][1] = -(py/pz)*(1/p - 1/pT);
783 // mD[fNc+1][2] = -((1/(pz*pz))*(pT - p) + 1/p);
784 // mD[fNc+1][3] = 0;
785 // mD[fNc+1][4] = (delX/delZ)*(1/T - 1/delT);
786 // mD[fNc+1][5] = (delY/delZ)*(1/T - 1/delT);
787 // mD[fNc+1][6] = (1/(delZ*delZ))*(delT - T) + 1/T;
788 //
789 // }
790 
791 /*
792  // better alternate way .........
793  md[fNc+0][0] = py/p*delX/delT-px/pT*delY/delT;
794  md[fNc+1][0]= pz/p*pT/T-pT/p*delZ/T;
795 
796  mD[fNc+0][0] = -((py*(dx*px + dy*py))/(delT*pT*pT*PT));
797  mD[fNc+0][1] = -((py*(dx*px + dy*py))/(delT*pT*pT*PT));
798  mD[fNc+0][2] = 0;
799  mD[fNc+0][3] = 0;
800  mD[fNc+0][4] = (delY*(delX*px + delY*py))/(delT*delT*delT*pT) ;
801  mD[fNc+0][5] = -((delX*(delX*px + delY*py))/(delT*delT*delT*pT)) ;
802  mD[fNc+0][6] = 0;
803 
804  //2nd equation
805  mD[fNc+1][0] = -((px*pz*(delT*pT + dz*pz))/(T*pT*(p*p*p)); //correct
806  mD[fNc+1][1] = -((py*pz*(delT*pT + dz*pz))/(T*pT*(p*p*p));
807  mD[fNc+1][2] = (delT*pT*pT + dz*pT*pz)/(delT*p*p*p) ;
808  mD[fNc+1][3] = 0;
809  mD[fNc+1][4] = (dx*dz*(delT*pT + dz*pz))/(delT*T*T*T*p);
810  mD[fNc+1][5] = (dy*dz*(delT*pT + dz*pz))/(delT*T*T*T*p);
811  mD[fNc+1][6] = (-(delT*delT*pT) - delT*dz*pz)/(delT*delT*delT*p);
812 
813  }
814  */
815 
816 void RhoKinHyperonVtxFitter::TransportToPoca(TMatrixD& a_in, TMatrixD& a_cov_in, TMatrixD& a_out,TMatrixD& a_cov_out ,TMatrixD& xref)
817 {
818  // by Xinying Song
819  // Transport the p7 and cov7 to the point nearest to poca or the vertex from fit
820  // it should be a iteration process
821 
822  int fNDau=fDaughters.size();
823  int nd=fNDau;
824  TMatrixD U(7*nd,7*nd);
825  int kN=0;
826  for(int k=0; k<nd; k++) {
827  kN=7*k;
828 
829  //double m = fDaughters[k]->Mass(); //[R.K. 01/2017] unused variable?
830  //check, if daughter particle is either neutral or charged
831  if (fabs(fDaughters[k]->GetCharge())<1e-6){//begin neutral
832  //Get position, energy and momentum for the daughter particle
833  double px=a_in[kN+0][0];
834  double py=a_in[kN+1][0];
835  double pz=a_in[kN+2][0];
836  //double p2=px*px+py*py+pz*pz; //[R.K. 01/2017] unused variable?
837  //double E= sqrt(m*m+p2);
838  double x=a_in[kN+4][0];
839  double y=a_in[kN+5][0];
840  double z=a_in[kN+6][0];
841  double ptot=sqrt(pz*pz+py*py+px*px);
842 
843  double t=(px*(xref[0][0]-x)+ py*(xref[1][0]-y)+pz*(xref[2][0]-z))/(ptot*ptot);
844 
845  //correct particle position to vertex, keep momentum and energy
846  a_out[kN+0][0] = px;
847  a_out[kN+1][0] = py;
848  a_out[kN+2][0] = pz;
849  a_out[kN+3][0] = a_in[kN+3][0];
850  a_out[kN+4][0] = x+px*t;
851  a_out[kN+5][0] = y+py*t;
852  a_out[kN+6][0] = z+pz*t;
853 
854  U[kN+0][kN+0] = 1.;
855  U[kN+1][kN+1] = 1.;
856  U[kN+2][kN+2] = 1.;
857  U[3+kN][3+kN] = 1.;
858 
859  U[4+kN][4+kN] = 1.;
860  U[4+kN][0+kN] = t;
861 
862  U[5+kN][5+kN] = 1.;
863  U[5+kN][1+kN] = t;
864 
865  U[6+kN][6+kN] = 1.;
866  U[6+kN][2+kN] = t;
867 
868  } // end neutral
869  else{
870  // for the charged track, the reference is Transporting Track parameters and covariance Matrix in magnetic field section III.4
871 
872  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
873  double a = -0.00299792458*bField*fDaughters[k]->GetCharge();
874 
875  //Get position, energy and momentum for the daughter particle
876 
877  double px=a_in[kN+0][0];
878  double py=a_in[kN+1][0];
879  double pz=a_in[kN+2][0];
880  double x=a_in[kN+4][0];
881  double y=a_in[kN+5][0];
882  double z=a_in[kN+6][0];
883 
884  double ptot = sqrt(px*px + py*py + pz*pz);
885  //double E=sqrt(m*m+ptot*ptot);
886  double rho = a/ptot; //1/R with R: the radius of the trajectory
887 
888  TVector3 x0(x,y,z);
889  TVector3 xp(xref[0][0], xref[1][0], xref[2][0]);
890  TVector3 delta= x0-xp;
891  TVector3 p0(px,py,pz);
892  TVector3 h(0,0,1);// assume magnetic field is only in z direction
893  double c = delta.Dot(p0);
894  double b = ptot- rho* delta.Dot( p0.Cross(h) );
895  double delta_h = delta.Dot(h);
896  double p_h = p0.Dot(h);
897  double A = -1.0* ( c - delta_h*p_h)*rho*rho/2.0;
898  // A*s^2 + b*s + c =0
899  double s1, s2,s;
900  if ((b*b-4*A*c)<0) {s1=s2=0;}
901  else{
902  s1=(-1.0*b -sqrt(b*b-4*A*c))/(2.0*A);
903  s2=(-1.0*b +sqrt(b*b-4*A*c))/(2.0*A);
904 
905  }
906 
907  double rs1=rho*s1;
908  double rs2=rho*s2;
909 
910  TVector3 x1(x+px/a*sin(rs1)-py/a*(1-cos(rs1)), y + py/a*sin(rs1)+px/a*(1-cos(rs1)), z+pz/ptot*s1 );
911  TVector3 x2(x+px/a*sin(rs2)-py/a*(1-cos(rs2)), y + py/a*sin(rs2)+px/a*(1-cos(rs2)), z+pz/ptot*s2 );
912  double d1= (x1-xp).Mag();
913  double d2= (x2-xp).Mag();
914  if(d1<=d2) {s=s1;}
915  else{s=s2;}
916  double cos_rho_s=cos(rho*s);
917  double sin_rho_s=sin(rho*s);
918 
919  a_out[kN+0][0] = px*cos_rho_s-py*sin_rho_s;
920  a_out[kN+1][0] = py*cos_rho_s+px*sin_rho_s;
921  a_out[kN+2][0] = pz;
922  a_out[kN+3][0] = a_in[kN+3][0];
923  a_out[kN+4][0] = x + (px*sin_rho_s - py*(1-cos_rho_s))/a;
924  a_out[kN+5][0] = y + (py*sin_rho_s + px*(1-cos_rho_s))/a;
925  a_out[kN+6][0] = z + (pz/ptot)*s;
926 
927 
928  //matrix to calculate the corrected covariance matrix
929  U[kN+0][kN+0] = cos_rho_s;
930  U[kN+0][kN+1] = -sin_rho_s;
931 
932  U[kN+1][kN+0] = sin_rho_s;
933  U[kN+1][kN+1] = cos_rho_s;
934 
935  U[kN+2][kN+2] = 1.;
936  U[3+kN][3+kN] = 1.;
937 
938  U[4+kN][0+kN] = sin_rho_s/a;
939  U[4+kN][1+kN] = -1.*(1.-cos_rho_s)/a;
940  U[4+kN][4+kN] = 1.;
941 
942  U[5+kN][0+kN] = (1.-cos_rho_s)/a;
943  U[5+kN][1+kN] = sin_rho_s/a;
944  U[5+kN][5+kN] = 1.;
945 
946  U[6+kN][2+kN] = s/ptot;
947  U[6+kN][6+kN] = 1.;
948 
949 
950  }// end charged
951 
952  }//end k
953 
954  TMatrixD U_t=U;
955  U_t=U_t.T();
956  a_cov_out = U*a_cov_in*U_t;
957 
958  }
959 
960 
961  void RhoKinHyperonVtxFitter::TransportToVertex(TMatrixD& a_in, TMatrixD& a_cov_in, TMatrixD& a_out, TMatrixD& a_cov_out, TMatrixD& xref)
962  {
963  //modified by Xinying Song, original version is from J.Puetz
964  //different from TransportToPoca(), this is for the daughters after vertex fit. their momentum is transported to the vertex fitted.
965 
966  int fNDau=fDaughters.size();
967  int nd=fNDau;
968  TMatrixD U(7*nd,7*nd);
969  int kN=0;
970  for(int k=0; k<nd; k++) {
971  kN=7*k;
972 
973  //double m = fDaughters[k]->Mass(); //[R.K. 01/2017] unused variable?
974  //check, if daughter particle is either neutral or charged
975  if (fabs(fDaughters[k]->GetCharge())<1e-6){//begin neutral
976 
977  //Get position, energy and momentum for the daughter particle
978  double px=a_in[kN+0][0];
979  double py=a_in[kN+1][0];
980  double pz=a_in[kN+2][0];
981  //double p2=px*px+py*py+pz*pz; //[R.K. 01/2017] unused variable?
982  //double E= sqrt(m*m+p2);
983  double x=a_in[kN+4][0];
984  double y=a_in[kN+5][0];
985  double z=a_in[kN+6][0];
986 
987  //correct particle position to vertex, keep momentum and energy
988  a_out[kN+0][0] = px;
989  a_out[kN+1][0] = py;
990  a_out[kN+2][0] = pz;
991  a_out[kN+3][0] = a_in[kN+3][0];
992  a_out[kN+4][0] = xref[0][0];
993  a_out[kN+5][0] = xref[1][0];
994  a_out[kN+6][0] = xref[2][0];
995  double s=sqrt((xref[0][0]-x)*(xref[0][0]-x)+(xref[1][0]-y)*(xref[1][0]-y) + (xref[2][0]-z)*(xref[2][0]-z));
996  double ptot=sqrt(px*px+py*py+pz*pz);
997  double t=s/ptot;
998 
999  //matrix U corrects the covariant matrix a_cov_in to a_cov_out= U * a_cov_in * U^T
1000  //note that x_new= x0- px*t;
1001 
1002  U[kN+0][kN+0] = 1.;
1003  U[kN+1][kN+1] = 1.;
1004  U[kN+2][kN+2] = 1.;
1005  U[3+kN][3+kN] = 1.;
1006 
1007  U[4+kN][4+kN] = 1.;
1008  U[4+kN][0+kN] = -1.*t;
1009 
1010  U[5+kN][5+kN] = 1.;
1011  U[5+kN][1+kN] = -1.*t;
1012 
1013  U[6+kN][6+kN] = 1.;
1014  U[6+kN][2+kN] = -1.*t;
1015 
1016 
1017  }//end neutral
1018 
1019  else{ // if particle is charged
1020 
1021  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
1022  double a = -0.00299792458*bField*fDaughters[k]->GetCharge();
1023 
1024  // if(fVerbose)cout << "a " << a << endl;
1025 
1026  //Get position, energy and momentum for the daughter particle
1027 
1028  double px=a_in[kN+0][0];
1029  double py=a_in[kN+1][0];
1030  double pz=a_in[kN+2][0];
1031  double x=a_in[kN+4][0];
1032  double y=a_in[kN+5][0];
1033  double z=a_in[kN+6][0];
1034 
1035  double ptot = sqrt(px*px + py*py + pz*pz);
1036  //double E=sqrt(m*m+ptot*ptot);
1037  double rho = a/ptot; //1/R with R: the radius of the trajectory
1038 
1039  //double s=(xref[2][0]-z)*ptot/pz;
1040  double cos_rho_s = 1+a*(py*(xref[0][0]-x)-px*(xref[1][0]-y))/(ptot*ptot-pz*pz);
1041  double sin_rho_s = a*(px*(xref[0][0]-x) + py*(xref[1][0]-y))/(ptot*ptot-pz*pz);
1042 
1043 
1044 
1045  //double s = atan2(sin_rho_s,cos_rho_s);
1046  double s = atan2(sin_rho_s,cos_rho_s)/rho ;
1047  a_out[kN+0][0] = px*cos_rho_s-py*sin_rho_s;
1048  a_out[kN+1][0] = py*cos_rho_s+px*sin_rho_s;
1049  a_out[kN+2][0] = pz;
1050  a_out[kN+3][0] = a_in[kN+3][0];
1051  a_out[kN+4][0] = x + (px*sin_rho_s - py*(1-cos_rho_s))/a;
1052  a_out[kN+5][0] = y + (py*sin_rho_s + px*(1-cos_rho_s))/a;
1053  a_out[kN+6][0] = z + (pz/ptot)*s;
1054 
1055 
1056  //matrix to calculate the corrected covariance matrix
1057  U[kN+0][kN+0] = cos_rho_s;
1058  U[kN+0][kN+1] = -sin_rho_s;
1059 
1060  U[kN+1][kN+0] = sin_rho_s;
1061  U[kN+1][kN+1] = cos_rho_s;
1062 
1063  U[kN+2][kN+2] = 1.;
1064  U[3+kN][3+kN] = 1.;
1065 
1066  U[4+kN][0+kN] = sin_rho_s/a;
1067  U[4+kN][1+kN] = -1.*(1.-cos_rho_s)/a;
1068  U[4+kN][4+kN] = 1.;
1069 
1070  U[5+kN][0+kN] = (1.-cos_rho_s)/a;
1071  U[5+kN][1+kN] = sin_rho_s/a;
1072  U[5+kN][5+kN] = 1.;
1073 
1074  U[6+kN][2+kN] = s/ptot;
1075  U[6+kN][6+kN] = 1.;
1076 
1077 
1078  }// end charged
1079 
1080  }//end k
1081 
1082  TMatrixD U_t=U;
1083  U_t=U_t.T();
1084  a_cov_out = U*a_cov_in*U_t;
1085 
1086  }
1087 
1088 
1089 
1090 
1091 
1092 //void RhoKinHyperonVtxFitter::TransportToVertex(TMatrixD& a_in, TMatrixD& a_cov_in, TMatrixD& a_out, TMatrixD& a_cov_out, TMatrixD& xref)
1093 //{
1094 // //edited by J.Puetz
1095 // //added parametrization for neutral daughter particles
1096 //
1097 // ///Correct one helix and/or track(s) to vertex point for charged and/or neutral particles
1098 // int fNDau=fDaughters.size();
1099 // int nd=fNDau;
1100 // TMatrixD U(7*nd,7*nd);
1101 // int kN=0;
1102 // for(int k=0; k<nd; k++) {
1103 // kN=7*k;
1104 //
1105 // //double m = fDaughters[k]->Mass(); //[R.K. 01/2017] unused variable?
1106 // //check, if daughter particle is either neutral or charged
1107 // if (fabs(fDaughters[k]->GetCharge())<1e-6){//begin neutral
1108 //
1109 // //Get position, energy and momentum for the daughter particle
1110 // double px=a_in[kN+0][0];
1111 // double py=a_in[kN+1][0];
1112 // double pz=a_in[kN+2][0];
1113 // //double p2=px*px+py*py+pz*pz; //[R.K. 01/2017] unused variable?
1114 // //double E= sqrt(m*m+p2);
1115 // double x=a_in[kN+4][0];
1116 // double y=a_in[kN+5][0];
1117 // double z=a_in[kN+6][0];
1118 //
1119 // //correct particle position to vertex, keep momentum and energy
1120 // a_out[kN+0][0] = px;
1121 // a_out[kN+1][0] = py;
1122 // a_out[kN+2][0] = pz;
1123 // a_out[kN+3][0] = a_in[kN+3][0];
1124 // a_out[kN+4][0] = x+px;
1125 // a_out[kN+5][0] = y+py;
1126 // a_out[kN+6][0] = z+pz;
1127 //
1128 // //matrix U corrects the covariant matrix a_cov_in to a_cov_out= U * a_cov_in * U^T
1129 // U[kN+0][kN+0] = 1.;
1130 // U[kN+1][kN+1] = 1.;
1131 // U[kN+2][kN+2] = 1.;
1132 // U[3+kN][3+kN] = 1.;
1133 //
1134 // U[4+kN][4+kN] = 1.;
1135 // U[4+kN][0+kN] = 1.;
1136 //
1137 // U[5+kN][5+kN] = 1.;
1138 // U[5+kN][1+kN] = 1.;
1139 //
1140 // U[6+kN][6+kN] = 1.;
1141 // U[6+kN][2+kN] = 1.;
1142 //
1143 //
1144 // }//end neutral
1145 //
1146 // else{ // if particle is charged
1147 //
1148 // Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
1149 // double a = -0.00299792458*bField*fDaughters[k]->GetCharge();
1150 //
1151 // // if(fVerbose)cout << "a " << a << endl;
1152 //
1153 // //Get position, energy and momentum for the daughter particle
1154 //
1155 // double px=a_in[kN+0][0];
1156 // double py=a_in[kN+1][0];
1157 // double pz=a_in[kN+2][0];
1158 // double x=a_in[kN+4][0];
1159 // double y=a_in[kN+5][0];
1160 // double z=a_in[kN+6][0];
1161 //
1162 // double ptot = sqrt(px*px + py*py + pz*pz);
1163 // //double E=sqrt(m*m+ptot*ptot);
1164 // double rho = a/ptot; //1/R with R: the radius of the trajectory
1165 // double A1 = 1 - pow(pz/ptot,2) - ( (x-xref[0][0])*py - (y-xref[1][0])*px )*rho/ptot ;
1166 // double A2 = (x-xref[0][0])*px + (y-xref[1][0])*py;
1167 // A2 = A2/ptot;
1168 //
1169 //
1170 // double det = sqrt(A1*A1+rho*rho*A2*A2);
1171 // double cos_rho_s = A1/det;
1172 // double sin_rho_s = -rho*A2/det;
1173 //
1174 //
1175 //
1176 // //double s = atan2(sin_rho_s,cos_rho_s);
1177 // double s = atan2(sin_rho_s,cos_rho_s)/rho ;
1178 // a_out[kN+0][0] = px*cos_rho_s-py*sin_rho_s;
1179 // a_out[kN+1][0] = py*cos_rho_s+px*sin_rho_s;
1180 // a_out[kN+2][0] = pz;
1181 // a_out[kN+3][0] = a_in[kN+3][0];
1182 // a_out[kN+4][0] = x + (px*sin_rho_s - py*(1-cos_rho_s))/a;
1183 // a_out[kN+5][0] = y + (py*sin_rho_s + px*(1-cos_rho_s))/a;
1184 // a_out[kN+6][0] = z + (pz/ptot)*s;
1185 //
1186 //
1187 // //matrix to calculate the corrected covariance matrix
1188 // U[kN+0][kN+0] = cos_rho_s;
1189 // U[kN+0][kN+1] = -sin_rho_s;
1190 //
1191 // U[kN+1][kN+0] = sin_rho_s;
1192 // U[kN+1][kN+1] = cos_rho_s;
1193 //
1194 // U[kN+2][kN+2] = 1.;
1195 // U[3+kN][3+kN] = 1.;
1196 //
1197 // U[4+kN][0+kN] = sin_rho_s/a;
1198 // U[4+kN][1+kN] = (1.-cos_rho_s)/a;
1199 // U[4+kN][4+kN] = 1.;
1200 //
1201 // U[5+kN][0+kN] = (1.-cos_rho_s)/a;
1202 // U[5+kN][1+kN] = sin_rho_s/a;
1203 // U[5+kN][5+kN] = 1.;
1204 //
1205 // U[6+kN][2+kN] = s/ptot;
1206 // U[6+kN][6+kN] = 1.;
1207 //
1208 //
1209 // }// end charged
1210 //
1211 // }//end k
1212 //
1213 // TMatrixD U_t=U;
1214 // U_t=U_t.T();
1215 // a_cov_out = U*a_cov_in*U_t;
1216 //
1217 //}
1218 
1219 void RhoKinHyperonVtxFitter::GetCovariance(TMatrixD& a_cov0, TMatrixD& cov_al_x, TMatrixD& V_vtx, TMatrixD& covS)
1220 //The covariance for the vitual particle ( a bit complicated)
1221 {
1222  int fNDau=fDaughters.size();
1223  int nd=fNDau;
1224 
1225 
1226  TMatrixD cov_al_x_temp=cov_al_x;
1227 
1228  TMatrixD pA(4*nd,7*nd);
1229  pA.Zero();
1230  double sumA=0;
1231  double a;
1232  int kN, jN;
1233  for (int k=0; k<nd; k++) {
1234  kN=k*7;
1235  jN=k*4;
1236  Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
1237  a = -0.00299792458*bField*fDaughters[k]->GetCharge();
1238  sumA += a;
1239  pA[jN][kN]=pA[jN+1][kN+1]=pA[jN+2][kN+2]=pA[jN+3][kN+3]=1;
1240  pA[jN+1][kN+4]=-a;
1241  pA[jN][kN+5]=a;
1242  }
1243 
1244  TMatrixD pB(4,3);
1245  TMatrixD pB_t(3,4);
1246  pB.Zero();
1247  pB[0][1]=-sumA;
1248  pB[1][0]=sumA;
1249 
1250  TMatrixD covA_al_xi=pA*cov_al_x;
1251 //TMatrixD covB_al_xi=pB*cov_al_x;
1252  TMatrixD cov_al_xT(3,7*nd);
1253  cov_al_xT = cov_al_x_temp.T();
1254 
1255  TMatrixD pA_t(7*nd,4*nd);
1256 //pA_t=pA.T();
1257 
1258 //TMatrixD covP(4*nd,4*nd);
1259  TMatrixD covP = pA*a_cov0*(pA_t.Transpose(pA));
1260  TMatrixD covA=covA_al_xi*(pB_t.Transpose(pB));
1261  TMatrixD covB=pB*(cov_al_xT)*(pA_t.Transpose(pA));
1262  TMatrixD covBV=pB*(V_vtx)*(pB_t.Transpose(pB));
1263 
1264  TMatrixD covA_al_x(4,3);
1265  TMatrixD SumcovP(4,4);
1266  SumcovP=covBV;
1267 
1268  for (int k=0; k<nd; k++) {
1269  kN=k*7;
1270  jN=k*4;
1271  for (int i=0; i<4; i++) {
1272  for (int j=0; j<4; j++) {
1273  SumcovP[i][j] += covP[jN+i][jN+j];
1274  SumcovP[i][j] += covA[jN+i][j];
1275  SumcovP[i][j] += covB[i][jN+j];
1276  }
1277  for (int jj=0; jj<3; jj++) {
1278  covA_al_x[i][jj] += covA_al_xi[jN+i][jj];
1279  }
1280  }
1281  }
1282 
1283 
1284  TMatrixD covPX(4,3);
1285  covPX = covA_al_x ;
1286  covPX += (pB *V_vtx);
1287  TMatrixD covPX_t(3,4);
1288  covPX_t=covPX_t.Transpose(covPX);
1289 
1290 
1291 
1292  covS.SetSub(0,0,V_vtx);
1293  covS.SetSub(0,3,covPX_t);
1294  covS.SetSub(3,0,covPX);
1295  covS.SetSub(3,3,SumcovP);
1296 
1297 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
RhoCandidate::SetP7
void SetP7(const TVector3 &pos, const TLorentzVector &p4)
Definition: RhoCandidate.cxx:393
x0
Double_t x0
Definition: checkhelixhit.C:70
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
RhoCandidate::SetPos
void SetPos(const TVector3 &pos)
Definition: RhoCandidate.h:235
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
RhoKinHyperonVtxFitter::FitNode
Bool_t FitNode(RhoCandidate *b)
Definition: RhoKinHyperonVtxFitter.cxx:54
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
RhoKinHyperonVtxFitter::AddMassConstraint
void AddMassConstraint(double mass)
Definition: RhoKinHyperonVtxFitter.cxx:41
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
s
TLorentzVector s
Definition: Pnd2DStar.C:50
ptot
Double_t ptot
Definition: checkhelixhit.C:68
RhoKinHyperonVtxFitter::GetCovariance
void GetCovariance(TMatrixD &a_cov0, TMatrixD &cov_al_x, TMatrixD &V_vtx, TMatrixD &covS)
Definition: RhoKinHyperonVtxFitter.cxx:1219
jj
Int_t jj
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:434
a
Int_t a
Definition: anaLmdDigi.C:126
RhoVtxPoca::GetPocaVtx
Double_t GetPocaVtx(TVector3 &vertex, RhoCandidate *composite)
Definition: RhoVtxPoca.cxx:27
RhoKinHyperonVtxFitter::SetOutput
void SetOutput(RhoCandidate *head)
Definition: RhoKinHyperonVtxFitter.cxx:338
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
z
Double_t z
Definition: createRootFscGeometryFile.C:146
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
RhoFitterBase::fDaughters
std::vector< RhoCandidate * > fDaughters
Definition: RhoFitterBase.h:69
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
p2
TPad * p2
Definition: hist-t7.C:117
atan2
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
x
Double_t x
Definition: createRootFscGeometryFile.C:159
RhoKinHyperonVtxFitter::TransportToVertex
void TransportToVertex(TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &)
Definition: RhoKinHyperonVtxFitter.cxx:961
fe
int fe
Definition: anaLmdDigi.C:67
RhoCalculationTools::GetBz
static Double_t GetBz(const TVector3 &position)
Return the magnetic field along the z-axis in kGs
Definition: RhoCalculationTools.cxx:22
RhoFitterBase::fChiSquare
double fChiSquare
Definition: RhoFitterBase.h:74
ClassImp
ClassImp(PndAnaContFact)
p1
TPad * p1
Definition: hist-t7.C:116
t
TTree * t
Definition: bump_analys.C:13
RhoKinHyperonVtxFitter::Compute
Bool_t Compute(RhoCandidate *c)
Definition: RhoKinHyperonVtxFitter.cxx:100
y
Double_t y
Definition: createRootFscGeometryFile.C:159
RhoFitterBase::SetDaugthersFromComposite
void SetDaugthersFromComposite(RhoCandidate *cand)
Definition: RhoFitterBase.cxx:136
RhoKinHyperonVtxFitter::RhoKinHyperonVtxFitter
RhoKinHyperonVtxFitter(RhoCandidate *b)
Definition: RhoKinHyperonVtxFitter.cxx:24
RhoFitterBase::SetDecayVertex
void SetDecayVertex(RhoCandidate *composite, const TVector3 &vtx, const TMatrixD &CovVV)
Definition: RhoFitterBase.cxx:178
RhoCandidate::SetCov7
void SetCov7(const TMatrixD &cov7)
Definition: RhoCandidate.cxx:413
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
pz
double pz[39]
Definition: pipisigmas.h:14
RhoKinHyperonVtxFitter::TransportToPoca
void TransportToPoca(TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &)
Definition: RhoKinHyperonVtxFitter.cxx:816
RhoFitterBase::fNDegreesOfFreedom
int fNDegreesOfFreedom
Definition: RhoFitterBase.h:75
Generated on Wed Apr 3 2019 10:03:00 for FairRoot/PandaRoot by   doxygen 1.8.5