FairRoot/PandaRoot
Static Public Member Functions | Private Member Functions | Static Private Member Functions | Static Private Attributes | List of all members
RhoCalculationTools Class Reference

#include <RhoCalculationTools.h>

Static Public Member Functions

static void SetVerbose (Int_t level)
 Allows debugging output for these utilities. More...
 
static Double_t GetBz (const TVector3 &position)
 Return the magnetic field along the z-axis in kGs More...
 
static void ForceConstantBz (Double_t bz=0.)
 Force a constant B field value for all positions. More...
 
static void UnForceConstantBz ()
 Release the B field to be fetched from the database. More...
 
static TMatrixDSym GetConverted6 (TMatrixDSym)
 
static TMatrixD GetConverted6 (TMatrixD)
 
static TMatrixD GetConverted7 (TMatrixD)
 
static TMatrixDSym GetCovMat (TMatrixD)
 
static TMatrixDSym GetCovMat1 (TMatrixD)
 
static TMatrixDSym GetFitError (TMatrixDSym)
 
static TMatrixD GetFitError (TLorentzVector, TMatrixD)
 
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. More...
 
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)
 
static Bool_t P6FromTrajectory (TVectorD &mom6, TMatrixDSym &cov6, RhoCandidate *cand, double z, double ztolerance)
 
static Double_t StateFromTrajectory (TVectorD &state, TMatrixDSym &cov, RhoCandidate *track, double vx, double vy, double vz, double ztolerance)
 
static void TransportToZ (RhoCandidate *cand, Double_t z=0)
 
static void PrintMatrix (TMatrixT< double > m)
 
static void PrintMatrix (TMatrixTSym< double > m)
 
static void PrintMatrix (RhoError m)
 

Private Member Functions

 RhoCalculationTools ()
 
 ~RhoCalculationTools ()
 

Static Private Member Functions

static std::ostream & bold_on (std::ostream &os)
 
static std::ostream & bold_off (std::ostream &os)
 

Static Private Attributes

static Double_t fBz =0.
 
static Int_t fVerbose =0
 
static Bool_t fBzSet =kFALSE
 

Detailed Description

Definition at line 23 of file RhoCalculationTools.h.

Constructor & Destructor Documentation

RhoCalculationTools::RhoCalculationTools ( )
private
RhoCalculationTools::~RhoCalculationTools ( )
private

Member Function Documentation

static std::ostream& RhoCalculationTools::bold_off ( std::ostream &  os)
inlinestaticprivate

Definition at line 85 of file RhoCalculationTools.h.

86  {
87  return os << "\e[0m";
88  }
static std::ostream& RhoCalculationTools::bold_on ( std::ostream &  os)
inlinestaticprivate

Definition at line 80 of file RhoCalculationTools.h.

81  {
82  return os << "\e[1m";
83  }
static void RhoCalculationTools::ForceConstantBz ( Double_t  bz = 0.)
inlinestatic

Force a constant B field value for all positions.

Definition at line 37 of file RhoCalculationTools.h.

References fBz, and fBzSet.

Referenced by ana_jpsi(), ana_jpsi_task(), ana_jpsi_task_ST(), ana_task_mult(), ana_triple_task(), PndPmtTask::Init(), poormantracks(), prod_ana(), prod_myana(), qa_softtrig(), quickana(), runpmt(), tut_ana_fast(), and tut_ana_task_d0().

37  {
38  fBz=bz;
39  fBzSet=kTRUE;
40  };
Double_t RhoCalculationTools::GetBz ( const TVector3 &  position)
static

Return the magnetic field along the z-axis in kGs

Definition at line 22 of file RhoCalculationTools.cxx.

References fBz, fBzSet, and pnt.

Referenced by DecayTreeFitter::ParticleBase::bFieldOverC(), RhoKinVtxFitter::GetCovariance(), RhoKinHyperonVtxFitter::GetCovariance(), RhoVtxPoca::GetPocaChargedToNeutral(), RhoVtxPoca::GetPocaTwoCharged(), P6FromTrajectory(), P7toHelix(), P7toPRG(), RhoKinVtxFitter::ReadKinMatrix(), RhoKinHyperonVtxFitter::ReadKinMatrix(), RhoKinVtxFitter::ReadMassKinMatrix(), RhoKinFitter::ReadMassKinMatrix(), RhoKinHyperonVtxFitter::ReadMassKinMatrix(), RhoKinVtxFitter::SetOutput(), RhoKinFitter::SetOutput(), RhoKinHyperonVtxFitter::SetOutput(), RhoKinHyperonVtxFitter::TransportToPoca(), RhoKinVtxFitter::TransportToVertex(), RhoKinHyperonVtxFitter::TransportToVertex(), and TransportToZ().

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 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
pnt
TClonesArray * pnt
Definition: checkdedx_helixhit.C:12
TMatrixDSym RhoCalculationTools::GetConverted6 ( TMatrixDSym  tmpMat)
static

Definition at line 69 of file RhoCalculationTools.cxx.

References i.

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixD RhoCalculationTools::GetConverted6 ( TMatrixD  tmpMat)
static

Definition at line 50 of file RhoCalculationTools.cxx.

References i.

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
TMatrixD RhoCalculationTools::GetConverted7 ( TMatrixD  tmpMat)
static

Definition at line 104 of file RhoCalculationTools.cxx.

References i.

Referenced by PndPidCorrelator::GetTrackInfo(), PndAnalysis::Propagator(), and PndAnalysis::ResetCandidate().

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
TMatrixDSym RhoCalculationTools::GetCovMat ( TMatrixD  tmpMat)
static

Definition at line 122 of file RhoCalculationTools.cxx.

References i.

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixDSym RhoCalculationTools::GetCovMat1 ( TMatrixD  tmpMat)
static

Definition at line 86 of file RhoCalculationTools.cxx.

References i.

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixDSym RhoCalculationTools::GetFitError ( TMatrixDSym  e)
static

Definition at line 140 of file RhoCalculationTools.cxx.

References i.

Referenced by PndPidCorrelator::GetTrackInfo(), PndAnalysis::Propagator(), and PndAnalysis::ResetCandidate().

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 }
i
Int_t i
Definition: run_full.C:25
TMatrixD RhoCalculationTools::GetFitError ( TLorentzVector  p,
TMatrixD  e 
)
static

Definition at line 158 of file RhoCalculationTools.cxx.

References i.

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 }
p
Double_t p
Definition: anasim.C:58
i
Int_t i
Definition: run_full.C:25
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
Bool_t RhoCalculationTools::P6FromTrajectory ( TVectorD &  mom6,
TMatrixDSym &  cov6,
RhoCandidate *  cand,
double  z,
double  ztolerance 
)
static

Definition at line 691 of file RhoCalculationTools.cxx.

References a, RhoCandidate::Charge(), RhoCandidate::Cov7(), dz, fabs(), GetBz(), i, L, mom, RhoCandidate::P4(), pos, RhoCandidate::Pos(), sin(), sqrt(), and z.

Referenced by DecayTreeFitter::RecoTrackStateProvider::P6FromTrack().

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 };
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
RhoCandidate::Pos
TVector3 Pos() const
Definition: RhoCandidate.h:186
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
mom
Double_t mom
Definition: plot_dirc.C:14
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
a
Int_t a
Definition: anaLmdDigi.C:126
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
RhoCandidate::Cov7
TMatrixD Cov7() const
Definition: RhoCandidate.cxx:1643
RhoCalculationTools::GetBz
static Double_t GetBz(const TVector3 &position)
Return the magnetic field along the z-axis in kGs
Definition: RhoCalculationTools.cxx:22
RhoCandidate::Charge
Double_t Charge() const
Definition: RhoCandidate.h:184
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
Bool_t RhoCalculationTools::P7toHelix ( const TVector3 &  pos,
const TLorentzVector &  p4,
const Double_t  Q,
const TMatrixD &  cov77,
Float_t *  helixparams,
TMatrixD &  helixCov,
Bool_t  skipcov = kFALSE 
)
static

Calculator functions.

Definition at line 238 of file RhoCalculationTools.cxx.

References CAMath::Abs(), CAMath::ATan2(), Double_t, fabs(), fVerbose, GetBz(), phi0, printf(), pz, sqrt(), CAMath::Sqrt(), and z0.

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 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
z0
Double_t z0
Definition: checkhelixhit.C:62
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
phi0
Double_t phi0
Definition: checkhelixhit.C:60
CAMath::ATan2
static T ATan2(const T &y, const T &x)
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
RhoCalculationTools::GetBz
static Double_t GetBz(const TVector3 &position)
Return the magnetic field along the z-axis in kGs
Definition: RhoCalculationTools.cxx:22
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
pz
double pz[39]
Definition: pipisigmas.h:14
Bool_t RhoCalculationTools::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 
)
static

Definition at line 396 of file RhoCalculationTools.cxx.

References CAMath::Abs(), CAMath::ATan2(), cos(), eta, fabs(), fVerbose, GetBz(), p2, phi0, Pi, printf(), pt(), pz, r, sin(), sqrt(), CAMath::Sqrt(), theta, tht(), CAMath::TwoPi(), and z0.

Referenced by RhoKalmanVtxFitter::CalcPrgParams().

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 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
z0
Double_t z0
Definition: checkhelixhit.C:62
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
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
pt
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
phi0
Double_t phi0
Definition: checkhelixhit.C:60
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
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
p2
TPad * p2
Definition: hist-t7.C:117
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
eta
TParticlePDG * eta
Definition: Pnd_Hc_etaee7G.C:74
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
pz
double pz[39]
Definition: pipisigmas.h:14
void RhoCalculationTools::PrintMatrix ( TMatrixT< double >  m)
static

Definition at line 926 of file RhoCalculationTools.cxx.

References i, and m.

Referenced by DecayTreeFitter::InteractionPoint::InteractionPoint(), DecayTreeFitter::InteractionPoint::projectBeamConstraint(), DecayTreeFitter::InternalParticle::projectConstraint(), DecayTreeFitter::ParticleBase::projectGeoConstraint(), DecayTreeFitter::InteractionPoint::projectIPConstraint(), DecayTreeFitter::RecoComposite::projectRecoComposite(), DecayTreeFitter::RecoPhoton::projectRecoConstraint(), DecayTreeFitter::RecoTrack::projectRecoConstraint(), and DecayTreeFitter::KalmanCalculator::updateCov().

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 }
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
void RhoCalculationTools::PrintMatrix ( TMatrixTSym< double >  m)
static

Definition at line 985 of file RhoCalculationTools.cxx.

References i, and 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 }
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
void RhoCalculationTools::PrintMatrix ( RhoError  m)
static

Definition at line 1036 of file RhoCalculationTools.cxx.

References i, and 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 }
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
static void RhoCalculationTools::SetVerbose ( Int_t  level)
inlinestatic

Allows debugging output for these utilities.

Definition at line 28 of file RhoCalculationTools.h.

References fVerbose.

28  {
29  fVerbose = level;
30  }
Double_t RhoCalculationTools::StateFromTrajectory ( TVectorD &  state,
TMatrixDSym &  cov,
RhoCandidate *  track,
double  vx,
double  vy,
double  vz,
double  ztolerance = -1. 
)
static

Definition at line 761 of file RhoCalculationTools.cxx.

References RhoCandidate::Charge(), RhoCandidate::Cov7(), fabs(), i, mom, p, p2, RhoCandidate::P4(), pos, RhoCandidate::Pos(), pz, sqrt(), and TransportToZ().

Referenced by DecayTreeFitter::RecoTrackStateProvider::state().

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 };
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
p
Double_t p
Definition: anasim.C:58
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
mom
Double_t mom
Definition: plot_dirc.C:14
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
p2
TPad * p2
Definition: hist-t7.C:117
RhoCalculationTools::TransportToZ
static void TransportToZ(RhoCandidate *cand, Double_t z=0)
Definition: RhoCalculationTools.cxx:844
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
pz
double pz[39]
Definition: pipisigmas.h:14
void RhoCalculationTools::TransportToZ ( RhoCandidate *  cand,
Double_t  z = 0 
)
static

Definition at line 844 of file RhoCalculationTools.cxx.

References a, RhoCandidate::Charge(), RhoCandidate::Cov7(), dz, GetBz(), i, L, mom, RhoCandidate::P4(), pos, RhoCandidate::Pos(), RhoCandidate::SetCov7(), RhoCandidate::SetP4(), RhoCandidate::SetPos(), sin(), and sqrt().

Referenced by StateFromTrajectory().

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 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
RhoCandidate::SetPos
void SetPos(const TVector3 &pos)
Definition: RhoCandidate.h:235
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
RhoCandidate::Pos
TVector3 Pos() const
Definition: RhoCandidate.h:186
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
mom
Double_t mom
Definition: plot_dirc.C:14
dz
Double_t dz
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:24
RhoCandidate::SetP4
void SetP4(Double_t mass, const TVector3 &p3)
Definition: RhoCandidate.cxx:373
a
Int_t a
Definition: anaLmdDigi.C:126
RhoCandidate::P4
TLorentzVector P4() const
Definition: RhoCandidate.h:195
z
Double_t z
Definition: createRootFscGeometryFile.C:146
RhoCandidate::Cov7
TMatrixD Cov7() const
Definition: RhoCandidate.cxx:1643
RhoCalculationTools::GetBz
static Double_t GetBz(const TVector3 &position)
Return the magnetic field along the z-axis in kGs
Definition: RhoCalculationTools.cxx:22
RhoCandidate::Charge
Double_t Charge() const
Definition: RhoCandidate.h:184
RhoCandidate::SetCov7
void SetCov7(const TMatrixD &cov7)
Definition: RhoCandidate.cxx:413
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
static void RhoCalculationTools::UnForceConstantBz ( )
inlinestatic

Release the B field to be fetched from the database.

Definition at line 42 of file RhoCalculationTools.h.

References fBzSet.

42  {
43  fBzSet=kFALSE;
44  };

Member Data Documentation

Double_t RhoCalculationTools::fBz =0.
staticprivate

Definition at line 93 of file RhoCalculationTools.h.

Referenced by ForceConstantBz(), and GetBz().

Bool_t RhoCalculationTools::fBzSet =kFALSE
staticprivate

Definition at line 95 of file RhoCalculationTools.h.

Referenced by ForceConstantBz(), GetBz(), and UnForceConstantBz().

Int_t RhoCalculationTools::fVerbose =0
staticprivate

Definition at line 94 of file RhoCalculationTools.h.

Referenced by P7toHelix(), P7toPRG(), and SetVerbose().

The documentation for this class was generated from the following files:
Generated on Wed Apr 3 2019 10:03:49 for FairRoot/PandaRoot by   doxygen 1.8.5