FairRoot/PandaRoot
RhoCalculationTools.cxx
Go to the documentation of this file.
1 // //
3 // RhoCalculationTools //
4 // //
5 // Toolset to calculate things centrally. //
6 // //
7 // Authors: //
8 // Ralf Kliemt, HIM/GSI Feb.2013 //
9 // //
11 
12 #include "RhoCalculationTools.h"
13 #include "FairRun.h"
14 #include "FairRunSim.h"
15 #include "FairRunAna.h"
16 #include "FairField.h"
17 
18 Double_t RhoCalculationTools::fBz=0.;
19 Int_t RhoCalculationTools::fVerbose=0;
20 Bool_t RhoCalculationTools::fBzSet=kFALSE;
21 
22 Double_t RhoCalculationTools::GetBz ( const TVector3& pos )
23 {
24  // If field was forced, we return that.
25  if ( kTRUE==fBzSet ) { return fBz; }
26  // Read magnetic filed strength in z direction [kGs]
27  double pnt[3], Bf[3];
28  pnt[0]=pos.X();
29  pnt[1]=pos.Y();
30  pnt[2]=pos.Z();
31  // retrieve the field from the framework
32  if( FairRun::Instance()->IsAna() )
33  {
34  FairRunAna::Instance()->GetField()->GetFieldValue ( pnt, Bf ); //[kGs]
35  } else {
36  FairRunSim::Instance()->GetField()->GetFieldValue ( pnt, Bf ); //[kGs]
37  }
38 
39  return Bf[2];
40 }
41 
42 
43 
44 
45 
46 
47 /*
48 Convert the Cov matrices from x,y,z,px,py,pz ==> px,py,pz,x,y,z
49 or vice versa*/
50 TMatrixD RhoCalculationTools::GetConverted6(TMatrixD tmpMat)
51 {
52  TMatrixD cMat(6,6);
53  for(Int_t i=0; i<6; i++) {
54  for(Int_t j=0; j<6; j++) {
55  if(i>=3) {
56  if(j>=3) { cMat[i-3][j-3] = tmpMat[i][j]; }
57  else { cMat[i-3][j+3] = tmpMat[i][j]; }
58  } else {
59  if(j>=3) { cMat[i+3][j-3] = tmpMat[i][j]; }
60  else { cMat[i+3][j+3] = tmpMat[i][j]; }
61  }
62  }
63  }
64  return cMat;
65 }
66 /*
67 Convert the Cov matrices from x,y,z,px,py,pz ==> px,py,pz,x,y,z
68 or vice versa*/
69 TMatrixDSym RhoCalculationTools::GetConverted6(TMatrixDSym tmpMat)
70 {
71  TMatrixDSym cMat(6);
72  for(Int_t i=0; i<6; i++) {
73  for(Int_t j=0; j<6; j++) {
74  if(i>=3) {
75  if(j>=3) { cMat[i-3][j-3] = tmpMat[i][j]; }
76  else { cMat[i-3][j+3] = tmpMat[i][j]; }
77  } else {
78  if(j>=3) { cMat[i+3][j-3] = tmpMat[i][j]; }
79  else { cMat[i+3][j+3] = tmpMat[i][j]; }
80  }
81  }
82  }
83  return cMat;
84 }
85 //convert cov matrices x,y,z,px,py,pz,E (7x7) ==> px,py,pz,E,x,y,z (7x7)
86 TMatrixDSym RhoCalculationTools::GetCovMat1(TMatrixD tmpMat)
87 {
88  TMatrixDSym cMat(7);
89  for(Int_t i=0; i<7; i++) {
90  for(Int_t j=0; j<7; j++) {
91  if(i>=3) {
92  if(j>=3) { cMat[i-3][j-3] = tmpMat[i][j]; }
93  else { cMat[i-3][j+3] = tmpMat[i][j]; }
94  } else {
95  if(j>=3) { cMat[i+4][j-3] = tmpMat[i][j]; }
96  else { cMat[i+4][j+4] = tmpMat[i][j]; }
97  }
98  }
99  }
100  return cMat;
101 }
102 
103 // Convert the Cov matrices from px,py,pz,E,x,y,z ==> x,y,z,px,py,pz,E
104 TMatrixD RhoCalculationTools::GetConverted7(TMatrixD tmpMat)
105 {
106  TMatrixD cMat(7,7);
107  for(Int_t i=0; i<7; i++) {
108  for(Int_t j=0; j<7; j++) {
109  if(i>=4) {
110  if(j>=4) { cMat[i-4][j-4] = tmpMat[i][j]; }
111  else { cMat[i-4][j+3] = tmpMat[i][j]; }
112  } else {
113  if(j>=4) { cMat[i+3][j-4] = tmpMat[i][j]; }
114  else { cMat[i+3][j+3] = tmpMat[i][j]; }
115  }
116  }
117  }
118  return cMat;
119 }
120 
121 // Convert x,y,z,px,py,pz,E (7x7) ==> px,py,pz,x,y,z (6x6)
122 TMatrixDSym RhoCalculationTools::GetCovMat(TMatrixD tmpMat)
123 {
124  TMatrixDSym cMat(6);
125  for(Int_t i=0; i<6; i++) {
126  for(Int_t j=0; j<6; j++) {
127  if(i>=3) {
128  if(j>=3) { cMat[i-3][j-3] = tmpMat[i][j]; }
129  else { cMat[i-3][j+3] = tmpMat[i][j]; }
130  } else {
131  if(j>=3) { cMat[i+3][j-3] = tmpMat[i][j]; }
132  else { cMat[i+3][j+3] = tmpMat[i][j]; }
133  }
134  }
135  }
136  return cMat;
137 }
138 
139 //converting 7x7 cov matrix to 6x6 matrix
140 TMatrixDSym RhoCalculationTools::GetFitError(TMatrixDSym e)
141 {
142  TMatrixDSym err(6);
143  for(unsigned i=0; i<3; ++i) {
144  for(unsigned j=i; j<3; ++j) {
145  err[i][j] = err[j][i] = e[i][j];
146  err[3+i][3+j] = err[3+j][3+i] = e[4+i][4+j];
147  }
148  }
149  for(unsigned i=0; i<3; ++i) {
150  for(unsigned j=0; j<3; ++j) {
151  err[i][3+j] = err[3+j][i] = e[i][4+j];
152  }
153  }
154  return err;
155 }
156 
157 // converting 6x6 cov matrices to 7x7
158 TMatrixD RhoCalculationTools::GetFitError(TLorentzVector p,TMatrixD e)
159 {
160 // Error(6x6,e) ==> Error(7x7,output(hsm)) using Momentum(p).
161  TMatrixD hsm(7,7);
162  for(unsigned i=0; i<3; ++i) {
163  for(unsigned j=i; j<3; ++j) {
164  hsm[i][j] = hsm[j][i] = e[i][j];
165  hsm[4+i][4+j] = hsm[4+j][4+i] = e[3+i][3+j];
166  }
167  }
168  for(unsigned i=0; i<3; ++i) {
169  for(unsigned j=0; j<3; ++j) {
170  hsm[i][4+j] = hsm[4+j][i] = e[i][3+j];
171  }
172  }
173  double invE = 1./p.E();
174  hsm[0][3] = hsm[3][0] = (p.X()*hsm[0][0]+p.Y()*hsm[0][1]+p.Z()*hsm[0][2])*invE;
175  hsm[1][3] = hsm[3][1] = (p.X()*hsm[0][1]+p.Y()*hsm[1][1]+p.Z()*hsm[1][2])*invE;
176  hsm[2][3] = hsm[3][2] = (p.X()*hsm[0][2]+p.Y()*hsm[1][2]+p.Z()*hsm[2][2])*invE;
177  hsm[3][3] = (p.X()*p.X()*hsm[0][0]+p.Y()*p.Y()*hsm[1][1]+p.Z()*p.Z()*hsm[2][2]
178  +2.0*p.X()*p.Y()*hsm[0][1]
179  +2.0*p.X()*p.Z()*hsm[0][2]
180  +2.0*p.Y()*p.Z()*hsm[1][2])*invE*invE;
181  hsm[3][4] = hsm[4][3] = (p.X()*hsm[0][4]+p.Y()*hsm[1][4]+p.Z()*hsm[2][4])*invE;
182  hsm[3][5] = hsm[5][3] = (p.X()*hsm[0][5]+p.Y()*hsm[1][5]+p.Z()*hsm[2][5])*invE;
183  hsm[3][6] = hsm[6][3] = (p.X()*hsm[0][6]+p.Y()*hsm[1][6]+p.Z()*hsm[2][6])*invE;
184  return hsm;
185 }
186 
187 
188 
189 
190 // Bool_t RhoCalculationTools::FillHelixParams(RhoCandidate* cand, Bool_t skipcov)
191 // {
192 // if (0==cand) {
193 // Error("FillHelixParams():", "RhoCandidate pointer is zero.");
194 // return kFALSE;
195 // }
196 // // write helix parameters & helix cov to RhoCandidate/TFitParams
197 // Float_t helixparams[5];
198 // TMatrixD helixcov(5,5);
199 // Bool_t rc = P7toHelix(cand->Pos(), cand->P4(), cand->GetCharge(), cand->Cov7(), helixparams, helixcov, skipcov);
200 // if (!rc) {Warning("FillHelixParams()","P7toHelix failed"); return kFALSE;}
201 //
202 // cand->SetHelixParms(helixparams);
203 //
204 // if (fVerbose>2) {
205 // std::cout<<"calculated helix Params:"
206 // <<"\n\tD0 ="<<helixparams[0]<<"\t cm"
207 // <<"\n\tPhi0 ="<<helixparams[1]<<"\t rad"
208 // <<"\n\tRho ="<<helixparams[2]<<"\t 1/cm"
209 // <<"\n\tZ0 ="<<helixparams[3]<<"\t cm"
210 // <<"\n\tcotTh ="<<helixparams[4]<<"\t "
211 // <<std::endl;
212 // }
213 //
214 // if(!skipcov){
215 // Float_t rhohelixcov[15];
216 // rhohelixcov[0] = helixcov[0][0];
217 // rhohelixcov[1] = helixcov[1][0];
218 // rhohelixcov[2] = helixcov[2][0];
219 // rhohelixcov[3] = helixcov[3][0];
220 // rhohelixcov[4] = helixcov[4][0];
221 // rhohelixcov[5] = helixcov[1][1];
222 // rhohelixcov[6] = helixcov[2][1];
223 // rhohelixcov[7] = helixcov[3][1];
224 // rhohelixcov[8] = helixcov[4][1];
225 // rhohelixcov[9] = helixcov[2][2];
226 // rhohelixcov[10] = helixcov[3][2];
227 // rhohelixcov[11] = helixcov[4][2];
228 // rhohelixcov[12] = helixcov[3][3];
229 // rhohelixcov[13] = helixcov[4][3];
230 // rhohelixcov[14] = helixcov[4][4];
231 // cand->SetHelixCov(rhohelixcov);
232 // }
233 //
234 //
235 // return kTRUE;
236 // }
237 
238 Bool_t RhoCalculationTools::P7toHelix(const TVector3& pos, const TLorentzVector& p4, const Double_t Q,
239  const TMatrixD& cov77, Float_t* helixparams, TMatrixD& helixCov, Bool_t skipcov)
240 {
241  // Convert from fourmomentum (vx,vy,vz,px,py,pz,e)
242  // to RHO helix parameters (D0,Phi0,rho(omega),Z0,tan(dip))
243  // Assuming vx,vy,vz give the POCA to the z axis.
244  if(p4.Perp()< 1e-9) {Warning("P7toHelix","Too small transverse momentum: %g",p4.Perp()); return kFALSE;}
245 // Double_t pnt[3], Bf[3];
246 // pnt[0]=pos.X();
247 // pnt[1]=pos.Y();
248 // pnt[2]=pos.Z();
249 // FairRunAna::Instance()->GetField()->GetFieldValue(pnt, Bf); //[kGs]
250 // //Double_t B = sqrt(Bf[0]*Bf[0]+Bf[1]*Bf[1]+Bf[2]*Bf[2]);
251 // Double_t B = Bf[2]; // assume field in z only
252 // //Double_t B = 20.;
253 
254  Double_t B = GetBz(pos);
255  if(fVerbose>1) { printf("P7ToHelix: BField is %g kGs\n",B); }
256  Double_t qBc = -0.000299792458*B*Q;//Mind factor from momenta being in GeV
257  //Double_t pti=1/p4.Perp();
258  //Double_t dfi=(pos.Phi()-p4.Phi())*TMath::RadToDeg();
259  //if(dfi>180)dfi-=360;
260  //if(dfi<-180)dfi+=360;
261  //Double_t sign=-Q*((dfi>0)?1:-1); // TODO get a decent D0 sign!!!
262  //helixparams[0]=sign*pos.Perp(); //D0
263  //helixparams[1]=p4.Phi(); //phi0
264  //helixparams[2]=qBc*pti; //omega=rho=1/R[cm]=-2.998*B[kGs]*Q[e]/p_perp[GeV/c]
265  //helixparams[3]=pos.Z(); //z0
266  //helixparams[4]=p4.Pz()/p4.Perp(); //lambda(averey)=cot(theta)=tan(lambda(geane))
267  if(fVerbose>1) { printf("P7ToHelix: Charge is %g e-\n",Q); }
268  if(fVerbose>1) { printf("P7ToHelix: QBc is %g \n",qBc); }
269 
270  const double xp = pos.X();
271  const double yp = pos.Y();
272  const double zp = pos.Z();
273  const double px = p4.Px();
274  const double py = p4.Py();
275  const double pz = p4.Pz();
276  const double phip = TMath::ATan2(py,px);
277  const double pti = 1/TMath::Sqrt(px*px+py*py);
278  //const double phip = p4.Phi();
279  //const double pti = 1/p4.Perp();
280  if(fVerbose>1) { printf("P7ToHelix: P_t is %g GeV/c\n",p4.Perp()); }
281  if(fVerbose>1) { printf("P7ToHelix: P_t^-1 is %g c/GeV\n",pti); }
282 
283  // get rho
284  const double rho = qBc*pti;
285  if(fVerbose>1) { printf("P7ToHelix: rho is %g cm^-1\n",rho); }
286 
287  // get tan(dip)
288  const double tanDip=pz*pti;
289  const double R0 = 1./rho;
290  if(fVerbose>1) { printf("P7ToHelix: tanDip is %g \n",tanDip); }
291  if(fVerbose>1) { printf("P7ToHelix: R0 = rho^-1 is %g cm\n",R0); }
292 
293  //circle center
294  const double xc = xp - py/qBc;
295  const double yc = yp + px/qBc;
296  //const double xc = xp - R0*p4.Py()*pti;
297  //const double yc = yp + R0*p4.Px()*pti;
298  //const double xc = xp - R0*TMath::Sin(phip);
299  //const double yc = yp + R0*TMath::Cos(phip);
300  const double RCsq = xc*xc+yc*yc;
301  const double RC = TMath::Sqrt(xc*xc+yc*yc);
302 
303  //get phi0 at doca
304  const double phi0 = -TMath::ATan2(xc,yc);
305  if(fVerbose>1) { printf("P7ToHelix: phi0 is %g, phiP is %g, DeltaPhi = %g \n",phi0,phip,phip-phi0); }
306 
307  //get D0
308  const double D0 = RC - TMath::Abs(R0);
309  //const double x0 = D0*TMath::Cos(phi0);
310  //const double y0 = D0*TMath::Sin(phi0);
311  if(fVerbose>1) { printf("P7ToHelix: D0 is %g cm\n",D0); }
312 
313  //get z0
314  const double z0 = zp - pz*(phip-phi0)/qBc;
315  //const double z0 = zp - tanDip*R0*(phip-phi0);
316  if(fVerbose>1) { printf("P7ToHelix: z0 is %g cm\n",z0); }
317 
318  helixparams[0]=D0;
319  helixparams[1]=phi0;
320  helixparams[2]=rho;
321  helixparams[3]=z0;
322  helixparams[4]=tanDip; // == lambda
323  if(fVerbose>1) {
324  for(int ai=0; ai<5; ai++) {
325  std::cout<<"helixparams["<<ai<<"]="<<helixparams[ai]<<std::endl;
326  }
327  }
328 
329  if(!skipcov) {
330  TMatrixD jacobian(5,7);
331 
332  jacobian[0][0] = xc/RC; // dD0 / dvx
333  jacobian[0][1] = yc/RC; // dD0 /dvy
334  jacobian[0][2] = 0.; // dD0 /dvz
335  jacobian[0][3] = yc/(qBc*RC)-px*pti/fabs(qBc); // dD0 /dpx
336  jacobian[0][4] = -1.*xc/(qBc*RC)-py*pti/fabs(qBc); // dD0 /dpy
337  jacobian[0][5] = 0.; // dD0 /dpz
338  jacobian[0][6] = 0.; // dD0 /de
339 
340  jacobian[1][0] = yc/(RCsq); // dPhi0 /dvx
341  jacobian[1][1] = -1.*xc/(RCsq); // dPhi0 /dvy
342  jacobian[1][2] = 0.; // dPhi0 /dvz
343  jacobian[1][3] = -1.*xc/(RCsq*qBc); // dPhi0 /dpx
344  jacobian[1][4] = -1.*yc/(RCsq*qBc); // dPhi0 /dpy
345  jacobian[1][5] = 0.; // dPhi0 /dpz
346  jacobian[1][6] = 0.; // dPhi0 /de
347 
348  jacobian[2][0] = 0.; // drho /dvx
349  jacobian[2][1] = 0.; // drho /dvy
350  jacobian[2][2] = 0.; // drho /dvz
351  jacobian[2][3] = -1.*rho*px*pti*pti; // drho /dpx
352  jacobian[2][4] = -1.*rho*py*pti*pti; // drho /dpy
353  jacobian[2][5] = 0.; // drho /dpz
354  jacobian[2][6] = 0.; // drho /de
355 
356  jacobian[3][0] = pz*yc/(RCsq*qBc); // dZ0 /dvx
357  jacobian[3][1] = -1.*pz*xc/(RCsq*qBc); // dZ0 /dvy
358  jacobian[3][2] = 1.; // dZ0 /dvz
359  jacobian[3][3] = pz/qBc * (py*pti*pti + xc/(qBc*RCsq)); // dZ0 /dpx
360  jacobian[3][4] = -1.*pz/qBc * (px*pti*pti - yc/(qBc*RCsq)); // dZ0 /dpy
361  jacobian[3][5] = (phi0-phip)/qBc; // dZ0 /dpz
362  jacobian[3][6] = 0.; // dZ0 /de
363 
364  jacobian[4][0] = 0.; // dtLam /dvx
365  jacobian[4][1] = 0.; // dtLam /dvy
366  jacobian[4][2] = 0.; // dtLam /dvz
367  jacobian[4][3] = -1.*tanDip*px*pti*pti; // dtLam /dpx
368  jacobian[4][4] = -1.*tanDip*py*pti*pti; // dtLam /dpy
369  jacobian[4][5] = pti; // dtLam /dpz
370  jacobian[4][6] = 0.; // dtLam /de
371 
372  TMatrixD tempmat(jacobian,TMatrixD::kMult,cov77);
373  TMatrixD covrho(tempmat,TMatrixD::kMultTranspose,jacobian);
374  helixCov=covrho;
375 
376  if (fVerbose>2) {
377  std::cout<<"cov77: ";
378  cov77.Print();
379  //std::cout<<"sigmas: "; sigmas.Print();
380  std::cout<<"jacobian: ";
381  jacobian.Print();
382  std::cout<<"covrho (D0,Phi0,rho,Z0,tanDip): ";
383  covrho.Print();
384  if(fVerbose>1) {
385  std::cout<<"helixparams[0] = D0 \t= ("<<helixparams[0]<<" \t+- "<<sqrt(helixCov[0][0])<<") cm"<<std::endl;
386  std::cout<<"helixparams[1] = Phi0 \t= ("<<helixparams[1]<<" \t+- "<<sqrt(helixCov[1][1])<<") rad"<<std::endl;
387  std::cout<<"helixparams[2] = rho \t= ("<<helixparams[2]<<" \t+- "<<sqrt(helixCov[2][2])<<") 1/cm"<<std::endl;
388  std::cout<<"helixparams[3] = Z0 \t= ("<<helixparams[3]<<" \t+- "<<sqrt(helixCov[3][3])<<") cm"<<std::endl;
389  std::cout<<"helixparams[4] = tanDip\t= ("<<helixparams[4]<<" \t+- "<<sqrt(helixCov[4][4])<<")"<<std::endl;
390  }
391  }
392  } // skip cov or not
393  return kTRUE;
394 }
395 
396 Bool_t RhoCalculationTools::P7toPRG(const TVector3& pos, const TLorentzVector& p4, const Double_t Q, const TMatrixD& cov77,
397  const TVector3& expPoint, Float_t* helixparams, TMatrixD& helixCov, TMatrixD& jacobian, Bool_t skipcov)
398 {
399  // Convert from fourmomentum (vx,vy,vz,px,py,pz,e)
400  // to PRG helix parameters (epsilon,z0,theta,phi0,rho)
401  // These parameters should NOT be written into the RhoCandidate!
402 
403  // check for enough transverse momentum
404  const double pt = p4.Perp(); // transverse momentum
405  if(pt < 1e-9) {
406  Warning("P7toPRG","Too small transverse momentum: %g",pt);
407  return kFALSE;
408  }
409 
410  // get field
411  const double B = GetBz(pos); // [kGs]
412  const double qBc = -0.000299792458*B*Q; // [GeV/cm] using B/[kGs]
413  //const double qBc = -0.299792458*B*Q; // momenta being in GeV length in cm
414  if(fVerbose>1) { printf("P7toPRG: BField is %g kGs\n",B); }
415  if(fVerbose>1) { printf("P7toPRG: Charge is %g e-\n",Q); }
416  if(fVerbose>1) { printf("P7toPRG: QBc is %g \n",qBc); }
417 
418  // move system to expansion point
419  const double xp = pos.X()-expPoint.X(); // reconstructed point
420  const double yp = pos.Y()-expPoint.Y(); // reconstructed point
421  const double zp = pos.Z()-expPoint.Z(); // reconstructed point
422  const double px = p4.Px(); // reconstructed momentum
423  const double py = p4.Py(); // reconstructed momentum
424  const double pz = p4.Pz(); // reconstructed momentum
425  const double p2 = px*px+py*py+pz*pz;
426  const double p2i = (p2==0.)?0.:1./p2;
427  const double pti = 1/pt;
428 
429  if(fVerbose>1) { printf("P7toPRG: x is [%8g,%8g,%8g] cm \n",xp,yp,zp); }
430  if(fVerbose>1) { printf("P7toPRG: p is [%8g,%8g,%8g] GeV/c \n",px,py,pz); }
431  // phi_p
432  const double phip = p4.Phi();
433  if(fVerbose>1) { printf("P7toPRG: phi is %g \n",phip); }
434 
435 
436 
437  //check if daughter is charged or neutral
438  if(fabs(Q)<1e-6){
439 
440  /*
441  * edited by J. Puetz
442  * reference for the helixparameters for neutral particles
443  * "Vertex reconstruction by means of the method of Kalman filtering" from Rolf Luchsinger and Christoph Grab
444  * doi:10.1016/0010-4655(93)90055-H
445  */
446 
447  double r = TMath::Sqrt(xp*xp+yp*yp);
448  double eta = phip - TMath::ACos(xp/r);
449  double tht = p4.Theta();
450 
451  helixparams[0] = r*sin(eta);
452  helixparams[1] = zp - r*cos(eta)/tan(tht);
453  helixparams[2] = tht;
454  helixparams[3] = phip;
455  helixparams[4] = TMath::Sqrt(p2);
456 
457 
458  if(!skipcov){
459 
460  jacobian[0][0] = xp*sin(eta)/r + TMath::Sqrt(1-xp*xp/(r*r)) * cos(eta); //d(r*sin(eta))/dxp
461  jacobian[0][1] = yp*sin(eta)/r - xp*yp*cos(eta)/(r*r*TMath::Sqrt(1-xp*xp/(r*r))); //d(r*sin(eta))/dyp
462  jacobian[0][2] = 0.; //d(r*sin(eta))/dzp
463  jacobian[0][3] = 0.; //d(r*sin(eta))/dpx
464  jacobian[0][4] = 0.; //d(r*sin(eta))/dpy
465  jacobian[0][5] = 0.; //d(r*sin(eta))/dpz
466  jacobian[0][6] = 0.; //d(r*sin(eta))/dE
467 
468  jacobian[1][0] = TMath::Sqrt(1-xp*xp/(r*r))*sin(eta)/tan(tht) - xp/r * cos(eta)/tan(tht); //d(helixparam[1])/dxp
469  jacobian[1][1] = -(yp/r*cos(eta)/tan(tht)+xp*yp/(r*r)*sin(eta)/tan(tht)/TMath::Sqrt(1-xp*xp/(r*r))); //d(helixparam[1])/dyp
470  jacobian[1][2] = 1.; //d(helixparam[1])/dzp
471  jacobian[1][3] = 0.; //d(helixparam[1])/dpx
472  jacobian[1][4] = 0.; //d(helixparam[1])/dpy
473  jacobian[1][5] = 0.; //d(helixparam[1])/dpz
474  jacobian[1][6] = 0.; //d(helixparam[1])/dE
475 
476  jacobian[2][0] = 0.; // dTheta /dxp
477  jacobian[2][1] = 0.; // dTheta /dyp
478  jacobian[2][2] = 0.; // dTheta /dzp
479  jacobian[2][3] = px*pz*pti*p2i; // dTheta /dpx
480  jacobian[2][4] = py*pz*pti*p2i; // dTheta /dpy
481  jacobian[2][5] = -pt*p2i; // dTheta /dpz
482  jacobian[2][6] = 0.; // dTheta /dE
483 
484  jacobian[3][0] = 0.; // dPhi /dxp
485  jacobian[3][1] = 0.; // dPhi /dyp
486  jacobian[3][2] = 0.; // dPhi /dzp
487  jacobian[3][3] = -py*pti; // dPhi /dpx
488  jacobian[3][4] = px*pti; // dPhi /dpy
489  jacobian[3][5] = 0.; // dPhi /dpz
490  jacobian[3][6] = 0.; // dPhi /dE
491 
492  jacobian[4][0] = 0.; // dp /dxp
493  jacobian[4][1] = 0.; // dp /dyp
494  jacobian[4][4] = 0.; // dp /dzp
495  jacobian[4][3] =px*TMath::Sqrt(p2i); // dp /dpx
496  jacobian[4][4] =py*TMath::Sqrt(p2i);; // dp /dpy
497  jacobian[4][5] =pz*TMath::Sqrt(p2i);; // dp /dpz
498  jacobian[4][6] = 0.; // dp /dE
499  }
500  } //end neutral particle
501 
502  else{
503  // get rho & R0
504 
505  const double rho = qBc*pti;
506  const double R0 = 1./rho; // signed radius
507  if(fVerbose>1) { printf("P7toPRG: rho is %g cm^-1\n",rho); }
508  if(fVerbose>1) { printf("P7toPRG: P_t is %g GeV/c\n",pt); }
509  if(fVerbose>1) { printf("P7toPRG: P_t^-1 is %g c/GeV\n",pti); }
510  if(fVerbose>1) { printf("P7toPRG: R0 = rho^-1 is %g cm\n",R0); }
511 
512  // get theta
513  const double theta=p4.Theta();
514  if(fVerbose>1) { printf("P7toPRG: theta is %g \n",theta); }
515  //const double tanDip = TMath::Tan(0.5*TMath::Pi() - theta); //unused?
516 
517  //circle center in x-y projection
518  const double xc = xp - py/qBc;
519  const double yc = yp + px/qBc;
520  const double RCsq = xc*xc + yc*yc;
521  const double RC = TMath::Sqrt(RCsq);
522 
523  //get epsilon
524  //epsilon is the distance of closest approach with its sign defined as being positive
525  //if the origin lays at the lefthand side of the moving particle
526  const double epsilon = Q*(RC - TMath::Abs(R0));
527  if(fVerbose>1) { printf("P7toPRG: epsilon is %g cm\n",epsilon); }
528 
529  //get phi0 at doca
530  const double direction = TMath::Sign(1.,Q);
531  double phi0 = TMath::ATan2(yc,xc);
532  if (yc < 0) { phi0 += TMath::Pi(); }
533  else { phi0 -= TMath::Pi(); }
534  phi0 -= direction*0.5*TMath::Pi();
535  //get phi0 into [-pi,pi]
536  if (phi0 < -TMath::Pi()) { phi0 += TMath::TwoPi(); }
537  else if (phi0 > TMath::Pi()) { phi0 -= TMath::TwoPi(); }
538 
539  if(fVerbose>1) { printf("P7toPRG: phi0 = %.4g, \tphiP = %.4g, \tDeltaPhi = %.4g \tepsilon=%.4g, \tQ=%g\n",phi0,phip,phip-phi0, epsilon, Q); }
540  if(fVerbose>1) {
541  double xnew=epsilon*sin(phi0);
542  double ynew=-epsilon*cos(phi0);
543  printf("P7toPRG: xnew = %.4g, \tynew = %.4g, \t(x^2+y^2) = %.4g \tepsilon^2=%.4g\n",xnew,ynew,xnew*xnew+ynew*ynew, epsilon*epsilon);
544  }
545 
546  //get z0
547  const double z0 = zp - pz*(phip-phi0)/qBc;
548  //const double z0 = zp - tanDip*R0*(phip-phi0);
549  if(fVerbose>1) { printf("P7toPRG: z0 is %g cm from zp=%.3g cm\n",z0,zp); }
550 
551  helixparams[0]=epsilon;
552  helixparams[1]=z0;
553  helixparams[2]=theta;
554  helixparams[3]=phi0;
555  helixparams[4]=rho;
556  if(fVerbose>0 && skipcov) {
557  for(int ai=0; ai<5; ai++) {
558  std::cout<<"helixparams["<<ai<<"]="<<helixparams[ai]<<std::endl;
559  }
560  }
561 
562  if(!skipcov) {
563  //TMatrixD jacobian(5,7);
564 
565  jacobian[0][0] = Q*xc/RC; // dEpsilon / dvx
566  jacobian[0][1] = Q*yc/RC; // dEpsilon /dvy
567  jacobian[0][2] = 0.; // dEpsilon /dvz
568  jacobian[0][3] = Q*(yc/(qBc*RC)-TMath::Sign(pti,R0)*px/qBc); // dEpsilon /dpx
569  jacobian[0][4] = Q*(-xc/(qBc*RC)-TMath::Sign(pti,R0)*py/qBc); // dEpsilon /dpy
570  jacobian[0][5] = 0.; // dEpsilon /dpz
571  jacobian[0][6] = 0.; // dEpsilon /de
572 
573  jacobian[1][0] = -pz*yc/(RCsq*qBc); // dZ0 /dvx
574  jacobian[1][1] = pz*xc/(RCsq*qBc); // dZ0 /dvy
575  jacobian[1][2] = 1.; // dZ0 /dvz
576  jacobian[1][3] = (py*pti*pti + xc/(qBc*RCsq))*pz/qBc; // dZ0 /dpx
577  jacobian[1][4] = -(px*pti*pti - yc/(qBc*RCsq))*pz/qBc; // dZ0 /dpy
578  jacobian[1][5] = -(phip-phi0)/qBc; // dZ0 /dpz
579  jacobian[1][6] = 0.; // dZ0 /de
580  // ok
581 
582  jacobian[2][0] = 0.; // dTheta /dvx
583  jacobian[2][1] = 0.; // dTheta /dvy
584  jacobian[2][2] = 0.; // dTheta /dvz
585  jacobian[2][3] = px*pz*pti*p2i; // dTheta /dpx
586  jacobian[2][4] = py*pz*pti*p2i; // dTheta /dpy
587  jacobian[2][5] = -pt*p2i; // dTheta /dpz
588  jacobian[2][6] = 0.; // dTheta /de
589  //ok
590 
591  jacobian[3][0] = -yc/(RCsq); // dPhi0 /dvx
592  jacobian[3][1] = xc/(RCsq); // dPhi0 /dvy
593  jacobian[3][2] = 0.; // dPhi0 /dvz
594  jacobian[3][3] = xc/(RCsq*qBc); // dPhi0 /dpx
595  jacobian[3][4] = yc/(RCsq*qBc); // dPhi0 /dpy
596  jacobian[3][5] = 0.; // dPhi0 /dpz
597  jacobian[3][6] = 0.; // dPhi0 /de
598 
599  jacobian[4][0] = 0.; // drho /dvx
600  jacobian[4][1] = 0.; // drho /dvy
601  jacobian[4][2] = 0.; // drho /dvz
602  jacobian[4][3] = -rho*px*pti*pti; // drho /dpx
603  jacobian[4][4] = -rho*py*pti*pti; // drho /dpy
604  jacobian[4][5] = 0.; // drho /dpz
605  jacobian[4][6] = 0.; // drho /de
606  //ok
607  }
608  }//end charged particle
609 
610  TMatrixD tempmat(jacobian,TMatrixD::kMult,cov77);
611  TMatrixD covrho(tempmat,TMatrixD::kMultTranspose,jacobian);
612  helixCov=covrho;
613 
614  if (fVerbose>2) {
615  std::cout<<"cov77: ";
616  cov77.Print();
617  //std::cout<<"sigmas: "; sigmas.Print();
618  std::cout<<"jacobian: ";
619  jacobian.Print();
620  }
621  if (fVerbose>0) {
622  std::cout<<"covrho (epsilon,Z0,Theta,Phi0,rho): ";
623  covrho.Print();
624  std::cout<<"helixparams[0] = epsilon\t= ("<<helixparams[0]<<" \t+- "<<sqrt(helixCov[0][0])<<") cm"<<std::endl;
625  std::cout<<"helixparams[1] = Z0 \t= ("<<helixparams[1]<<" \t+- "<<sqrt(helixCov[1][1])<<") cm"<<std::endl;
626  std::cout<<"helixparams[2] = Theta \t= ("<<helixparams[2]<<" \t+- "<<sqrt(helixCov[2][2])<<") rad"<<std::endl;
627  std::cout<<"helixparams[3] = Phi0 \t= ("<<helixparams[3]<<" \t+- "<<sqrt(helixCov[3][3])<<") rad"<<std::endl;
628  std::cout<<"helixparams[4] = rho\t= ("<<helixparams[4]<<" \t+- "<<sqrt(helixCov[4][4])<<") 1/cm"<<std::endl;
629  }
630  return kTRUE;
631 }
632 
633 
634 // Bool_t RhoCalculationTools::SDtoHelix(FairTrackParH* par, RhoCandidate* cand, Bool_t skipcov)
635 // {
636 // // Adding the helix parameters to the candidate (TFitParams)
637 // // It is nice to know at analysis stage
638 // // rho helix: (D0,Phi0,rho(omega),Z0,tan(dip))
639 // // fair helix: (q/p,lambda, phi, y_perp, z_perp)
640 // Double_t Q=par->GetQ();
641 // if(0==Q) return kFALSE;
642 // // Double_t pnt[3], Bf[3];
643 // // pnt[0]=par->GetX();
644 // // pnt[1]=par->GetY();
645 // // pnt[2]=par->GetZ();
646 // // FairRunAna::Instance()->GetField()->GetFieldValue(pnt, Bf); //[kGs]
647 // // //Double_t B = sqrt(Bf[0]*Bf[0]+Bf[1]*Bf[1]+Bf[2]*Bf[2]);
648 // // Double_t B = Bf[2];
649 //
650 // Double_t B = GetBz(par->GetPosition());
651 //
652 // Double_t qBc = -0.000299792458*B*Q;
653 // Double_t icL = 1. / cos(par->GetLambda()); // inverted for practical reasons (better to multiply than to divide)
654 // Double_t icLs = icL*icL;
655 // Float_t helixparams[5];
656 // helixparams[0]=par->GetY_sc() ; //D0
657 // helixparams[1]=par->GetMomentum().Phi()+TMath::Pi(); //phi0
658 // helixparams[2]=qBc/(par->GetMomentum().Perp()); //omega=rho=1/R[cm]=-2.998e-4*B[kGs]*Q[e]/p_perp[GeV/c]
659 // helixparams[3]=par->GetZ_sc()*icL; //z0
660 // helixparams[4]=tan(par->GetLambda()); //lambda(averey)=cot(theta)=tan(lambda(geane))
661 // cand->SetHelixParms(helixparams);
662 //
663 // if(skipcov) return kTRUE; // stop here when skipping cov matrix calculation
664 //
665 // Double_t fairhelixcov[15];
666 // par->GetCov(fairhelixcov);
667 //
668 // // in the poca to z axis yperp=D0, x_perp^2+z_perp^2 = z_perp/cos(Lambda)= Z0
669 // Float_t rhohelixcov[15];
670 // rhohelixcov[0] = fairhelixcov[12]; // sigma^2 D0
671 // rhohelixcov[1] = fairhelixcov[10]; // cov D0 - Phi0
672 // rhohelixcov[2] = fairhelixcov[3] * qBc * icL; // cov D0 - rho
673 // rhohelixcov[3] = fairhelixcov[13] * icL; // cov D0 - Z0
674 // rhohelixcov[4] = fairhelixcov[7] * icLs; // cov D0 - tan(dip)
675 // rhohelixcov[5] = fairhelixcov[9]; // sigma^2 Phi0
676 // rhohelixcov[6] = fairhelixcov[2] * qBc * icL; // cov Phi0 - rho
677 // rhohelixcov[7] = fairhelixcov[11] * icL; // cov Phi0 - Z0
678 // rhohelixcov[8] = fairhelixcov[6] * icLs; // cov Phi0 - tan(dip)
679 // rhohelixcov[9] = fairhelixcov[0] * qBc * qBc * icLs; // sigma^2 rho
680 // rhohelixcov[10] = fairhelixcov[4] * qBc * icLs; // cov rho - Z0
681 // rhohelixcov[11] = fairhelixcov[1] * qBc * icL * icLs; // cov rho - tan(dip)
682 // rhohelixcov[12] = fairhelixcov[14] * icLs; // sigma^2 Z0
683 // rhohelixcov[13] = fairhelixcov[8] * icL * icLs; //cov Z0 - tan(dip)
684 // rhohelixcov[14] = fairhelixcov[5] * icLs * icLs; // sigma^2 tan(dip) - from
685 // cand->SetHelixCov(rhohelixcov);
686 //
687 // return kTRUE;
688 // }
689 
690 // --------------------------------------------
691 Bool_t RhoCalculationTools::P6FromTrajectory( TVectorD& mom6, TMatrixDSym& cov6, RhoCandidate* cand, Double_t z, Double_t ztolerance=-1. )
692 {
693  // here we transform from Pos&P4 of the RhoCandidate to DecayTreeFitter mom6
694  TVector3 pos = cand->Pos();
695  TLorentzVector mom = cand->P4();
696  TMatrixD cov7 = cand->Cov7();
697  if(ztolerance>0. && fabs(cand->Pos().Z()-z) > ztolerance)
698  {
699  // input/output is (x0,y0,z0,px0,py0,pz0)
700  // transporting on helix to (x1,y1,z,px1,py1,pz0) (z from input, pz stays)
701  const double dz = z-pos.Z();
702  const double px = mom.Px();
703  const double py = mom.Py();
704  const double pzi = 1./mom.Pz();
705  //const double A = dz*pzi;
706  const double B = GetBz(pos); // [kGs]
707  const double a = -0.000299792458*B*cand->Charge(); // [GeV/cm] using B/[kGs]
708  const double ai = 1./a;
709  const double rho = a/mom.P();
710  const double rhoTimesS = a*dz*pzi;
711  const double L=sin(rhoTimesS);
712  //const double K=cos(rhoTimesS);
713  const double K=sqrt(1-L*L);
714 
715  TMatrixD J(7,7);
716  for(int i=0;i<7;i++) J[i][i]=1.;
717  //J[2][0] = pzi*(px*K-py*L); // dx/dz0
718  //J[2][1] = pzi*(py*K+px*L); // dy/dz0
719  //J[2][3] = pzi*a*(-px*L-py*K); // dpx/dz0
720  //J[2][4] = pzi*a*(-py*L+px*K); // dpy/dz0
721  J[0][3] = L*ai; // dx/dpx0
722  J[1][3] = (1.-K)*ai; // dy/dpx0
723  J[0][4] = (K-1.)*ai; // dx/dpy0
724  J[1][4] = L*ai; // dy/dpy0
725  J[3][3] = K; // dpx/dpx0
726  J[4][3] = L; // dpy/dpx0
727  J[3][4] = -L; // dpx/dpy0
728  J[4][4] = K; // dpy/dpy0
729  J[3][5] = rho*pzi*(px*L+py*K); // dpx/dpz0
730  J[4][5] = rho*pzi*(py*L-px*K); // dpy/dpz0
731 
732  double x1 = pos.X() + (L*px-(1.-K)*py)*ai;
733  double y1 = pos.Y() + (L*py+(1.-K)*px)*ai;
734 
735  double px1 = px*K-py*L;
736  double py1 = py*K+px*L;
737  TMatrixD Jcov(J,TMatrixD::kMult,cov7);
738  TMatrixD newcov(Jcov,TMatrixD::kMultTranspose,J);
739 
740  pos.SetXYZ(x1,y1,z);
741  mom.SetX(px1);
742  mom.SetY(py1);
743  cov7=newcov;
744  }
745 
746  mom6[0] = pos.x();
747  mom6[1] = pos.y();
748  mom6[2] = pos.z();
749  mom6[3] = mom.Px();
750  mom6[4] = mom.Py();
751  mom6[5] = mom.Pz();
752 
753  for(int i=0;i<6;i++){
754  for(int j=0;j<6;j++){
755  cov6[i][j]=cov7[i][j];
756  }
757  }
758  return kTRUE;
759 };
760 
761 Double_t RhoCalculationTools::StateFromTrajectory( TVectorD& state, TMatrixDSym& cov, RhoCandidate* cand,
762  double , double , double vz, double ztolerance=-1. ) // vx vy //[R.K.03/2017] unused variable(s)
763 {
764  //std::cout<<"RhoCalculationTools::StateFromTrajectory()"<<std::endl;
765  // here we transform from Pos&P4 of the RhoCandidate to DecayTreeFitter state
766  // (x,y,tx,ty,qoverp)
767 
768  //-----
769 // TVector3 pos0=cand->GetPosition();
770 // TVector3 mom0=cand->GetMomentum();
771 // double dx0 = vx-pos0.x(); // connect cand position with Vertex V
772 // double dy0 = vy-pos0.y();
773 // double dz0 = vz-pos0.z();
774 // double px0 = mom0.x();
775 // double py0 = mom0.y();
776 // double pz0 = mom0.z();
777 // double pt0 = mom0.Pt();
778 // double pt0i = 1./pt0;
779 // double cf0 = px0*pt0i;
780 // double sf0 = py0*pt0i;
781 // double tl0 = pz0*pt0i;
782 // double B0 = GetBz(pos0); // [kGs]
783 // double a0 = -0.000299792458*B0*cand->Charge(); // [GeV/cm] using B/[kGs]
784 // double a0i = 1./a0;
785 // double R0 = pt0*a0i;
786 // double dd = dx0*sf0 - dy0*cf0 + R0*(tl0*tl0 + 1);
787 // double ddi = 1./dd;
788 // double ddiq = ddi*ddi;
789 // double dal = ddi*(dx0*cf0 + dy0*sf0 + dz0*tl0);
790 // double dz1 = tl0*R0*dal; // z shift on helix until poca of V
791 // double myz = cand->Pos().Z() + dz1;
792  double myz = vz; //requested z
793  //-----
794 
795  //if(ztolerance>0. && fabs(dz0) > ztolerance) TransportToZ(cand,myz);
796  if(ztolerance>0. && fabs(cand->Pos().Z()-vz) > ztolerance) TransportToZ(cand,vz);
797  TVector3 pos=cand->Pos();
798  TLorentzVector mom = cand->P4();
799  //TVector3 mom=cand->GetMomentum();
800  double px=mom.Px();
801  double py=mom.Py();
802  double pz=mom.Pz();
803  double pzi=1./pz;
804  double p2=px*px+py*py+pz*pz;
805  double p2i=1./p2;
806  double p=sqrt(p2);
807 
808  if (pz == 0. || p == 0.) {std::cout<<"RhoCalculationTools::StateFromTrajectory(): Zero protection: pz="<<pz<<" Gev/c p="<<p<<" GeV/c"<<std::endl; return kFALSE;}
809  state[0] = pos.x();
810  state[1] = pos.y();
811  state[2] = px*pzi;
812  state[3] = py*pzi;
813  state[4] = cand->Charge()/p;
814 
815  TMatrixD J(5,7);
816  J[0][0] = 1.;
817  J[1][1] = 1.;
818  J[2][3] = pzi;
819  J[2][5] = -pzi*state[2];
820  J[3][4] = pzi;
821  J[3][5] = -pzi*state[3];
822  J[4][3] = -px*state[4]*p2i;
823  J[4][4] = -py*state[4]*p2i;
824  J[4][5] = -pz*state[4]*p2i;
825 
826  TMatrixD cov7 = cand->Cov7();
827  TMatrixD tempmat(J,TMatrixD::kMult,cov7);
828  TMatrixD cov5(tempmat,TMatrixD::kMultTranspose,J);
829  for(int i=0;i<5;i++){
830  for(int j=0;j<5;j++){
831  cov[i][j]=cov5[i][j];
832  }
833  }
834 
835 //std::cout<<"RhoCalculationTools::StateFromTrajectory() particle = "<<cand->PdgCode()<<" Q = "<<cand->Charge()<<std::endl;
836 //std::cout<<"RhoCalculationTools::StateFromTrajectory() p7 (x,y,z,px,py,pz,E) = ("<<pos.x()<<","<<pos.y()<<","<<pos.z()<<","<<px<<","<<py<<","<<pz<<","<<cand->E()<<")"<<std::endl;
837 //std::cout<<"RhoCalculationTools::StateFromTrajectory() state (x,y,(z),tx,ty,q/p) = ("<<state[0]<<","<<state[1]<<","<<myz<<","<<state[2]<<","<<state[3]<<","<<state[4]<<")"<<std::endl;
838 
839  return myz;
840 
841 };
842 
843 
844 void RhoCalculationTools::TransportToZ(RhoCandidate* cand, Double_t z)
845 {
846  //std::cout<<"RhoCalculationTools::TransportToZ(): z="<<z<<std::endl;
847  // input/output is (x0,y0,z0,px0,py0,pz0)
848  // transporting on helix to (x1,y1,z,px1,py1,pz0) (z from input, pz stays)
849  TVector3 pos = cand->Pos();
850  TLorentzVector mom = cand->P4();
851  TMatrixD cov7 = cand->Cov7();
852  const double dz = z-pos.Z();
853  const double px = mom.Px();
854  const double py = mom.Py();
855  const double pzi = 1./mom.Pz();
856  //const double A = dz*pzi;
857  const double B = GetBz(pos); // [kGs]
858  const double a = -0.000299792458*B*cand->Charge(); // [GeV/cm] using B/[kGs]
859  const double ai = 1./a;
860  const double rho = a/mom.P();
861  const double rhoTimesS = a*dz*pzi;
862  const double L=sin(rhoTimesS);
863  //const double K=cos(rhoTimesS);
864  const double K=sqrt(1-L*L);
865 
866  TMatrixD J(7,7);
867  for(int i=0;i<7;i++) J[i][i]=1.;
868  //J[2][0] = pzi*(px*K-py*L); // dx/dz0
869  //J[2][1] = pzi*(py*K+px*L); // dy/dz0
870  //J[2][3] = pzi*a*(-px*L-py*K); // dpx/dz0
871  //J[2][4] = pzi*a*(-py*L+px*K); // dpy/dz0
872  J[0][3] = L*ai; // dx/dpx0
873  J[1][3] = (1.-K)*ai; // dy/dpx0
874  J[0][4] = (K-1.)*ai; // dx/dpy0
875  J[1][4] = L*ai; // dy/dpy0
876  J[3][3] = K; // dpx/dpx0
877  J[4][3] = L; // dpy/dpx0
878  J[3][4] = -L; // dpx/dpy0
879  J[4][4] = K; // dpy/dpy0
880  J[3][5] = rho*pzi*(px*L+py*K); // dpx/dpz0
881  J[4][5] = rho*pzi*(py*L-px*K); // dpy/dpz0
882 
883  double x1 = pos.X() + (L*px-(1.-K)*py)*ai;
884  double y1 = pos.Y() + (L*py+(1.-K)*px)*ai;
885 
886  double px1 = px*K-py*L;
887  double py1 = py*K+px*L;
888  TMatrixD Jcov(J,TMatrixD::kMult,cov7);
889  TMatrixD newcov(Jcov,TMatrixD::kMultTranspose,J);
890 
891  //std::cout<<"RhoCalculationTools::TransportToZ(z="<<z<<") with K="<<std::setprecision(9)<<K<<", L="<<L<<", a="<<a<<", rho="<<rho<<", rhoTimesS="<<rhoTimesS<<std::endl;
892  //std::cout<<"RhoCalculationTools::TransportToZ(): original pos = ( "
893  // <<std::setw(6)<< pos.X() <<", "<<std::setw(6)<< pos.Y() <<", "<<std::setw(6)<< pos.Z()
894  // <<") mom = ("<<std::setw(6)<< px <<", "<<std::setw(6)<< py <<", "<<std::setw(6)<<std::setw(6)<< mom.Pz() <<", "<< mom.E()<<")"<<std::endl;
895 
896  pos.SetXYZ(x1,y1,z);
897  mom.SetX(px1);
898  mom.SetY(py1);
899  //std::cout<<"RhoCalculationTools::TransportToZ(): updated pos = ( "
900  // <<std::setw(6)<< pos.X() <<", "<<std::setw(6)<< pos.Y() <<", "<<std::setw(6)<< pos.Z()
901  // <<") mom = ("<<std::setw(6)<< px1 <<", "<<std::setw(6)<< py1 <<", "<<std::setw(6)<<std::setw(6)<< mom.Pz() <<", "<< mom.E()<<")"<<std::endl;
902 
903  //std::cout<<"RhoCalculationTools::TransportToZ(): original cov:"<<std::endl; cov7.Print();
904  //std::cout<<"RhoCalculationTools::TransportToZ(): jacobian: "<<std::endl; J.Print();
905  //std::cout<<"RhoCalculationTools::TransportToZ(): updated cov: "<<std::endl; newcov.Print();
906 
907  cand->SetPos(pos);
908  cand->SetP4(mom);
909  cand->SetCov7(newcov);
910  //std::cout<<"RhoCalculationTools::TransportToZ(): done."<<std::endl;
911 }
912 
913 
914 
915 
916 
917 
918 
919 
920 
921 
922 
923 
924 
925 //______________________________________________________________________________
926 void RhoCalculationTools::PrintMatrix(TMatrixT<double> m)
927 {
928  // Print the matrix as a table of elements.
929  // By default the format "%11.4g" is used to print one element.
930  // One can specify an alternative format with eg
931  // option ="f= %6.2f "
932  if (!m.IsValid()) {
933  Error("Print","Matrix is invalid");
934  return;
935  }
936  const Int_t ncols = m.GetNcols();
937  const Int_t nrows = m.GetNrows();
938  const Int_t collwb = m.GetColLwb();
939  const Int_t rowlwb = m.GetRowLwb();
940  const Int_t nelems = m.GetNoElements();
941  //build format
942  const char *format = "%11.4g ";
943  //if (option) {
944  // const char *f = strstr(option,"f=");
945  // if (f) format = f+2;
946  //}
947  char topbar[25];
948  const char *topbarstart="-----";
949 
950  snprintf(topbar,25,format,123.456789);
951  Int_t nch = strlen(topbar);//+1;
952  if (nch > 18) nch = 18;
953 
954  const char *ftopbar=Form(" %s%dd |","%",nch-5);
955  std::cout<<std::endl<<nrows<<"x"<<ncols<<" matrix is as follows"<<std::endl;
956  //Int_t nk = 5+nch*ncols;
957  //for (Int_t i = 0; i < nk; i++) topbar[i] = '-';
958  //topbar[nk] = 0;
959  for (Int_t i = 0; i < nch; i++) topbar[i] = '-';
960  topbar[nch] = 0;
961 
962  std::cout<<std::endl<<std::endl<<" |"<<std::flush;
963  for (Int_t j = 1; j <= ncols; j++)
964  {
965  std::cout<<Form(ftopbar,j+collwb-1)<<std::flush;
966  }
967  std::cout<<std::endl<<topbarstart;
968  for (Int_t j = 1; j <= ncols; j++) std::cout<<topbar;
969  std::cout<<std::endl;
970 
971  if (nelems <= 0) return;
972  for (Int_t i = 1; i <= nrows; i++)
973  {
974  std::cout<<Form("%4d |",i+rowlwb-1)<<std::flush;
975  for (Int_t j = 1; j <= ncols; j++)
976  {
977  std::cout<<Form(format,m(i+rowlwb-1,j+collwb-1))<<std::flush;
978  }
979  std::cout<<std::endl;
980  }
981  std::cout<<std::endl;
982 }
983 
984 //______________________________________________________________________________
985 void RhoCalculationTools::PrintMatrix(TMatrixTSym<double> m)
986 {
987  if (!m.IsValid()) {
988  Error("Print","Matrix is invalid");
989  return;
990  }
991  const Int_t ncols = m.GetNcols();
992  const Int_t nrows = m.GetNrows();
993  const Int_t collwb = m.GetColLwb();
994  const Int_t rowlwb = m.GetRowLwb();
995  const Int_t nelems = m.GetNoElements();
996  //build format
997  const char *format = "%11.4g ";
998  //if (option) {
999  // const char *f = strstr(option,"f=");
1000  // if (f) format = f+2;
1001  //}
1002  char topbar[25];
1003  const char *topbarstart="-----";
1004  snprintf(topbar,25,format,123.456789);
1005  Int_t nch = strlen(topbar);//+1;
1006  if (nch > 18) nch = 18;
1007  if (nch < 7) nch = 7;
1008  const char *ftopbar=Form(" %s%dd |","%",nch-5);
1009  std::cout<<std::endl<<nrows<<"x"<<ncols<<" matrix is as follows"<<std::endl;
1010  for (Int_t i = 0; i < nch; i++) topbar[i] = '-';
1011  topbar[nch] = 0;
1012 
1013  std::cout<<std::endl<<std::endl<<" SYM |"<<std::flush;
1014  for (Int_t j = 1; j <= ncols; j++)
1015  {
1016  std::cout<<Form(ftopbar,j+collwb-1)<<std::flush;
1017  }
1018  std::cout<<std::endl<<topbarstart;
1019  for (Int_t j = 1; j <= ncols; j++) std::cout<<topbar;
1020  std::cout<<std::endl;
1021 
1022  if (nelems <= 0) return;
1023  for (Int_t i = 1; i <= nrows; i++)
1024  {
1025  std::cout<<Form("%4d |",i+rowlwb-1)<<std::flush;
1026  for (Int_t j = 1; j <= ncols; j++)
1027  {
1028  //std::cout<<Form(format,m(i+rowlwb-1,j+collwb-1))<<std::flush;
1029  std::cout<< (i==j ? "\e[1m" : "\e[0m") <<Form(format,m(i+rowlwb-1,j+collwb-1))<<"\e[0m"<<std::flush;
1030  }
1031  std::cout<<std::endl;
1032  }
1033  std::cout<<std::endl;
1034 }
1035 //______________________________________________________________________________
1036 void RhoCalculationTools::PrintMatrix(RhoError m)
1037 {
1038  if (!m.IsValid()) {
1039  Error("Print","Matrix is invalid");
1040  return;
1041  }
1042  const Int_t ncols = m.GetNcols();
1043  const Int_t nrows = m.GetNrows();
1044  const Int_t collwb = m.GetColLwb();
1045  const Int_t rowlwb = m.GetRowLwb();
1046  const Int_t nelems = m.GetNoElements();
1047  //build format
1048  const char *format = "%11.4g ";
1049  //if (option) {
1050  // const char *f = strstr(option,"f=");
1051  // if (f) format = f+2;
1052  //}
1053  char topbar[25];
1054  const char *topbarstart ="-----";
1055  snprintf(topbar,25,format,123.456789);
1056  Int_t nch = strlen(topbar);//+1;
1057  if (nch > 18) nch = 18;
1058  if (nch < 7) nch = 7;
1059  const char *ftopbar=Form(" %s%dd |","%",nch-5);
1060  std::cout<<std::endl<<nrows<<"x"<<ncols<<" matrix is as follows"<<std::endl;
1061  for (Int_t i = 0; i < nch; i++) topbar[i] = '-';
1062  topbar[nch] = 0;
1063 
1064  std::cout<<std::endl<<std::endl<<" RHO |"<<std::flush;
1065  for (Int_t j = 1; j <= ncols; j++)
1066  {
1067  std::cout<<Form(ftopbar,j+collwb-1)<<std::flush;
1068  }
1069  std::cout<<std::endl<<topbarstart;
1070  for (Int_t j = 1; j <= ncols; j++) std::cout<<topbar;
1071  std::cout<<std::endl;
1072 
1073  if (nelems <= 0) return;
1074  for (Int_t i = 1; i <= nrows; i++)
1075  {
1076  std::cout<<Form("%4d |",i+rowlwb-1)<<std::flush;
1077  for (Int_t j = 1; j <= ncols; j++)
1078  {
1079  //std::cout<<Form(format,m(i+rowlwb-1,j+collwb-1))<<std::flush;
1080  std::cout<< (i==j ? "\e[1m" : "\e[0m") <<Form(format,m(i+rowlwb-1,j+collwb-1))<<"\e[0m"<<std::flush;
1081  }
1082  std::cout<<std::endl;
1083  }
1084  std::cout<<std::endl;
1085 }
RhoCalculationTools::P6FromTrajectory
static Bool_t P6FromTrajectory(TVectorD &mom6, TMatrixDSym &cov6, RhoCandidate *cand, double z, double ztolerance)
Definition: RhoCalculationTools.cxx:691
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
z0
Double_t z0
Definition: checkhelixhit.C:62
RhoCalculationTools::GetConverted6
static TMatrixDSym GetConverted6(TMatrixDSym)
Definition: RhoCalculationTools.cxx:69
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
r
double r
Definition: RiemannTest.C:14
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
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
RhoCandidate::Pos
TVector3 Pos() const
Definition: RhoCandidate.h:186
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
pnt
TClonesArray * pnt
Definition: checkdedx_helixhit.C:12
mom
Double_t mom
Definition: plot_dirc.C:14
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
p
Double_t p
Definition: anasim.C:58
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
RhoCandidate::SetP4
void SetP4(Double_t mass, const TVector3 &p3)
Definition: RhoCandidate.cxx:373
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
RhoCalculationTools::P7toPRG
static Bool_t P7toPRG(const TVector3 &pos, const TLorentzVector &p4, const Double_t Q, const TMatrixD &cov77, const TVector3 &expPoint, Float_t *helixparams, TMatrixD &helixCov, TMatrixD &jacobian, Bool_t skipcov=kFALSE)
Definition: RhoCalculationTools.cxx:396
a
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
RhoCalculationTools::P7toHelix
static Bool_t P7toHelix(const TVector3 &pos, const TLorentzVector &p4, const Double_t Q, const TMatrixD &cov77, Float_t *helixparams, TMatrixD &helixCov, Bool_t skipcov=kFALSE)
Calculator functions.
Definition: RhoCalculationTools.cxx:238
CAMath::ATan2
static T ATan2(const T &y, const T &x)
tht
TString tht(TString pts, TString exts="px py pz")
Definition: invexp.C:168
RhoCalculationTools::PrintMatrix
static void PrintMatrix(TMatrixT< double > m)
Definition: RhoCalculationTools.cxx:926
RhoCandidate::P4
TLorentzVector P4() const
Definition: RhoCandidate.h:195
z
Double_t z
Definition: createRootFscGeometryFile.C:146
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
RhoCalculationTools::GetConverted7
static TMatrixD GetConverted7(TMatrixD)
Definition: RhoCalculationTools.cxx:104
p2
TPad * p2
Definition: hist-t7.C:117
RhoCandidate::Cov7
TMatrixD Cov7() const
Definition: RhoCandidate.cxx:1643
CAMath::TwoPi
static float TwoPi()
Definition: PndCAMath.h:61
RhoCalculationTools::GetBz
static Double_t GetBz(const TVector3 &position)
Return the magnetic field along the z-axis in kGs
Definition: RhoCalculationTools.cxx:22
RhoCalculationTools::GetFitError
static TMatrixDSym GetFitError(TMatrixDSym)
Definition: RhoCalculationTools.cxx:140
RhoCandidate::Charge
Double_t Charge() const
Definition: RhoCandidate.h:184
RhoCalculationTools::GetCovMat1
static TMatrixDSym GetCovMat1(TMatrixD)
Definition: RhoCalculationTools.cxx:86
RhoCalculationTools::StateFromTrajectory
static Double_t StateFromTrajectory(TVectorD &state, TMatrixDSym &cov, RhoCandidate *track, double vx, double vy, double vz, double ztolerance)
Definition: RhoCalculationTools.cxx:761
eta
TParticlePDG * eta
Definition: Pnd_Hc_etaee7G.C:74
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
RhoCalculationTools::TransportToZ
static void TransportToZ(RhoCandidate *cand, Double_t z=0)
Definition: RhoCalculationTools.cxx:844
RhoCandidate::SetCov7
void SetCov7(const TMatrixD &cov7)
Definition: RhoCandidate.cxx:413
RhoCalculationTools::GetCovMat
static TMatrixDSym GetCovMat(TMatrixD)
Definition: RhoCalculationTools.cxx:122
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
pz
double pz[39]
Definition: pipisigmas.h:14
Generated on Sun Aug 18 2019 03:01:56 for FairRoot/PandaRoot by   doxygen 1.8.8