Class Reference for E1039 Core & Analysis Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GenFitExtrapolator.cxx
Go to the documentation of this file.
1 #include "GenFitExtrapolator.h"
2 #include "GFField.h"
3 
4 #include <phool/recoConsts.h>
5 #include <phfield/PHFieldUtility.h>
6 #include <phfield/PHFieldConfig_v3.h>
7 #include <phfield/PHField.h>
8 
9 #include <GenFit/FieldManager.h>
10 #include <GenFit/MaterialEffects.h>
11 #include <GenFit/TGeoMaterialInterface.h>
12 #include <GenFit/MeasuredStateOnPlane.h>
13 #include <GenFit/SharedPlanePtr.h>
14 #include <GenFit/RKTrackRep.h>
15 
16 //ROOT
17 #include <TGeoManager.h>
18 #include <TMath.h>
19 
20 //
21 //#include <G4SystemOfUnits.hh>
22 
23 #include <memory>
24 #include <cassert>
25 
26 //#define _DEBUG_ON
27 
28 namespace
29 {
30  //static flag to indicate the initialized has been done
31  static bool inited = false;
32 
33  //Simple swimming settings
34  static int NSTEPS_TARGET = 100;
35  static int NSTEPS_SHIELDING = 50;
36  static int NSTEPS_FMAG = 100;
37 
38  static double FMAG_LENGTH;
39  static double Z_UPSTREAM;
40 
41  //initialize global variables
42  void initGlobalVariables()
43  {
44  if(!inited)
45  {
46  inited = true;
47  recoConsts* rc = recoConsts::instance();
48 
49  NSTEPS_TARGET = rc->get_IntFlag("NSTEPS_TARGET");
50  NSTEPS_SHIELDING = rc->get_IntFlag("NSTEPS_SHIELDING");
51  NSTEPS_FMAG = rc->get_IntFlag("NSTEPS_FMAG");
52 
53  FMAG_LENGTH = rc->get_DoubleFlag("FMAG_LENGTH");
54  Z_UPSTREAM = rc->get_DoubleFlag("Z_UPSTREAM");
55  }
56  }
57 }
58 
59 static const double c_light = 2.99792458e+8 * m/s;
60 
61 using namespace std;
62 
63 GenFitExtrapolator::GenFitExtrapolator():
64  pos_i(TVector3()), mom_i(TVector3()), cov_i(TMatrixDSym(5)),
65  pos_f(TVector3()), mom_f(TVector3()), cov_f(TMatrixDSym(5)),
66  jac_sd2sc(TMatrixD(5,5)), jac_sc2sd(TMatrixD(5,5))
67 {
68  initGlobalVariables();
69 }
70 
71 GenFitExtrapolator::~GenFitExtrapolator()
72  {}
73 
74 bool GenFitExtrapolator::init(const PHField* field, const TGeoManager *geom)
75 {
76  iParType = 1;
77 
78  assert(field);
79  SQGenFit::GFField *fieldMap = new SQGenFit::GFField(field);
80  genfit::FieldManager::getInstance()->init(fieldMap);
81 
82  _tgeo_manager = const_cast<TGeoManager*>(geom);
83 
84 #ifdef _DEBUG_ON
85  double z_test = 1000;
86  LogInfo("");
87  {
88  double p[4] = {0, 0, z_test*cm, 0};
89  double B[3] = {0, 0, 0};
90  field->GetFieldValue(p, B);
91  cout << "PHField (CLHEP) at Z = " << z_test << endl;
92  cout << B[0] << ", " << B[1] << ", " << B[2] << endl;
93  }
94  {
95  genfit::AbsBField *f = genfit::FieldManager::getInstance()->getField();
96  TVector3 H = f->get(TVector3(0,0,z_test));
97  H *= kilogauss/tesla;
98  cout << "genfit::AbsBField (tesla) at Z = " << z_test << endl;
99  H.Print();
100 
101  z_test = 250;
102  H = f->get(TVector3(0,0,z_test));
103  H *= kilogauss/tesla;
104  cout << "genfit::AbsBField (tesla) at Z = " << z_test << endl;
105  H.Print();
106  }
107 #endif
108 
109  _tgeo_manager->Export("GenFitExtrapolatorGeom.root");
110 
111  genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
112 
113  return true;
114 }
115 
116 void GenFitExtrapolator::setInitialStateWithCov(
117  double z_in,
118  TMatrixD& state_in,
119  TMatrixD& cov_in)
120 {
121 #ifdef _DEBUG_ON
122  LogInfo("z_in: ") << z_in << endl;
123  state_in.Print();
124  cov_in.Print();
125 #endif
126 
127  convertSVtoMP(z_in, state_in, mom_i, pos_i);
128 
129  if(state_in[0][0] > 0)
130  {
131  setParticleType(1);
132  }
133  else
134  {
135  setParticleType(-1);
136  }
137 
138  TMatrixD cov_sd(5, 5);
139  for(int i = 0; i < 5; i++)
140  {
141  for(int j = 0; j < 5; j++)
142  {
143  cov_sd[i][j] = cov_in[i][j];
144 
145  if(i == 0) cov_sd[i][j] = double(iParType)*cov_sd[i][j];
146  if(j == 0) cov_sd[i][j] = double(iParType)*cov_sd[i][j];
147  }
148  }
149 
151  TRSDSC(iParType, mom_i, pos_i);
152  TMatrixD jac_sd2sc_T = jac_sd2sc;
153  jac_sd2sc_T.T();
154 
155  TMatrixD cov_sc = jac_sd2sc*cov_sd*jac_sd2sc_T;
156  for(int i = 0; i < 5; i++)
157  {
158  for(int j = 0; j < 5; j++)
159  {
160  cov_i[i][j] = cov_sc[i][j];
161  }
162  }
163 }
164 
165 void GenFitExtrapolator::setParticleType(int type) {
166  iParType = type;
167 }
168 
169 void GenFitExtrapolator::getFinalStateWithCov(TMatrixD& state_out, TMatrixD& cov_out) {
170  convertMPtoSV(mom_f, pos_f, state_out);
171 
172  TMatrixD cov_sc(5, 5);
173  for(int i = 0; i < 5; i++)
174  {
175  for(int j = 0; j < 5; j++)
176  {
177  cov_sc[i][j] = cov_f[i][j];
178 
179  if(direction == GenFitExtrapolator::Backward)
180  {
181  if(i == 1) cov_sc[i][j] = -cov_sc[i][j];
182  if(j == 1) cov_sc[i][j] = -cov_sc[i][j];
183  if(i == 3) cov_sc[i][j] = -cov_sc[i][j];
184  if(j == 3) cov_sc[i][j] = -cov_sc[i][j];
185  }
186  }
187  }
188 
189  //convert from SC to SD error matrix
190  TRSCSD(iParType, mom_f, pos_f);
191  TMatrixD jac_sc2sd_T = jac_sc2sd;
192  jac_sc2sd_T.T();
193 
194  cov_out = jac_sc2sd*cov_sc*jac_sc2sd_T;
195 
196  for(int i = 0; i < 5; i++)
197  {
198  for(int j = 0; j < 5; j++)
199  {
200  if(i == 0) cov_out[i][j] = double(iParType)*cov_out[i][j];
201  if(j == 0) cov_out[i][j] = double(iParType)*cov_out[i][j];
202  }
203  }
204 }
205 
206 void GenFitExtrapolator::getPropagator(TMatrixD& prop) {
207 // if(fabs(pos_i[2] - pos_f[2]) < 1E-3)
208 // {
209 // prop.UnitMatrix();
210 // return;
211 // }
212 //
213 // for(int i = 0; i < 5; i++)
214 // {
215 // for(int j = 0; j < 5; j++)
216 // {
217 // prop[i][j] = g4eProp[i][j];
218 //
219 // if(i == 0) prop[i][j] = double(iParType)*prop[i][j];
220 // if(j == 0) prop[i][j] = double(iParType)*prop[i][j];
221 // if(direction == GenFitExtrapolator::Backward)
222 // {
223 // if(i == 1) prop[i][j] = -prop[i][j];
224 // if(j == 1) prop[i][j] = -prop[i][j];
225 // if(i == 3) prop[i][j] = -prop[i][j];
226 // if(j == 3) prop[i][j] = -prop[i][j];
227 // }
228 // }
229 // }
230 //
231 // prop = jac_sc2sd*prop*jac_sd2sc;
232 }
233 
234 bool GenFitExtrapolator::extrapolateTo(double z_out) {
235 
237  if(pos_i[2] > 2400 || pos_i[2] < -600 || z_out > 2400 || z_out < -600)
238  {
239  return false;
240  }
241 
243  if(fabs(pos_i[2] - z_out) < 1E-3)
244  {
245  mom_f = mom_i;
246  pos_f = pos_i;
247  cov_f = cov_i;
248 
249  return true;
250  }
251 
252  int pid = iParType > 0 ? -13 : 13;
253 #ifdef _DEBUG_ON
254  cout
255  << "extrapolateToPlane: "
256  << "From: " << pos_i.Z()
257  << " -> " << z_out
258  << endl;
259 #endif
260 
261  auto up_rep = std::unique_ptr<genfit::AbsTrackRep> (new genfit::RKTrackRep(pid));
262  auto rep = up_rep.get();
263 
264  genfit::SharedPlanePtr destPlane(new genfit::DetPlane(TVector3(0, 0, z_out), TVector3(0, 0, 1)));
265 
266  std::unique_ptr<genfit::MeasuredStateOnPlane> currentState = std::unique_ptr < genfit::MeasuredStateOnPlane > (new genfit::MeasuredStateOnPlane(rep));
267  //rep->setPosMomCov(*currentState, pos_i, mom_i, cov_i);
268  currentState->setPosMom(pos_i, mom_i);
269  currentState->setCov(cov_i);
270 
271  try {
272  travelLength = rep->extrapolateToPlane(*currentState, destPlane);
273  } catch (...) {
274 #ifdef _DEBUG_ON
275  cout << "extrapolateToPlane failed!" << endl;
276 #endif
277  return false;
278  }
279 
280  currentState->getPosMom(pos_f, mom_f);
281  cov_f = currentState->getCov();
282 
283 
284  return true;
285 }
286 
287 //int GenFitExtrapolator::propagate() {
288 //
289 // return -1;
290 //}
291 
292 double GenFitExtrapolator::extrapolateToIP() {
293 
294  //Store the steps on each point
295  TVector3 mom[NSTEPS_FMAG + NSTEPS_TARGET + 1];
296  TVector3 pos[NSTEPS_FMAG + NSTEPS_TARGET + 1];
297 
298  //Step size in FMAG/target area, unit is cm.
299  //FIXME Units
300  double step_fmag = FMAG_LENGTH/NSTEPS_FMAG;//*cm;
301  double step_target = fabs(Z_UPSTREAM)/NSTEPS_TARGET;//*cm;
302 
303  //Start from FMAG face downstream
304  extrapolateTo(FMAG_LENGTH);
305  pos[0] = pos_f;
306  mom[0] = mom_f;
307 
308  //Now make the real swimming
309  int iStep = 1;
310  for(; iStep <= NSTEPS_FMAG; ++iStep)
311  {
312  pos_i = pos[iStep-1];
313  mom_i = mom[iStep-1];
314 
315  extrapolateTo((pos_i[2] - step_fmag));
316 
317  pos[iStep] = pos_f;
318  mom[iStep] = mom_f;
319  }
320 
321  for(; iStep < NSTEPS_FMAG+NSTEPS_TARGET+1; ++iStep)
322  {
323  pos_i = pos[iStep-1];
324  mom_i = mom[iStep-1];
325 
326  extrapolateTo((pos_i[2] - step_target));
327 
328  pos[iStep] = pos_f;
329  mom[iStep] = mom_f;
330  }
331 
332  //Find the one step with minimum DCA
333  double dca_min = 1E6;
334  for(int i = 0; i < NSTEPS_FMAG+NSTEPS_TARGET+1; ++i)
335  {
336  double dca = sqrt(pos[i][0]*pos[i][0] + pos[i][1]*pos[i][1]);
337  if(dca < dca_min)
338  {
339  dca_min = dca;
340  iStep = i;
341  }
342  }
343 
344  return pos[iStep][2];
345 }
346 void GenFitExtrapolator::convertSVtoMP(double z, TMatrixD& state, TVector3& mom, TVector3& pos)
347 {
348  double p = fabs(1. / state[0][0]);
349  double pz = p / sqrt(1. + state[1][0] * state[1][0] + state[2][0] * state[2][0]);
350  double px = pz * state[1][0];
351  double py = pz * state[2][0];
352 
353  double x = state[3][0];
354  double y = state[4][0];
355 
356  pos.SetXYZ(x, y, z);
357  mom.SetXYZ(px, py, pz);
358 }
359 
360 void GenFitExtrapolator::convertMPtoSV(TVector3& mom, TVector3& pos, TMatrixD& state) {
361  TVector3 mom_gev = mom;
362  TVector3 pos_cm = pos;
363 
364  state[0][0] = double(iParType) / mom_gev.Mag();
365  state[1][0] = mom_gev.x() / mom_gev.z();
366  state[2][0] = mom_gev.y() / mom_gev.z();
367 
368  state[3][0] = pos_cm.x();
369  state[4][0] = pos_cm.y();
370 }
371 
372 void GenFitExtrapolator::TRSDSC(int charge, TVector3 mom_input, TVector3 pos_input)
373 {
375  TVector3 mom(mom_input[0], mom_input[1], mom_input[2]);
376  TVector3 pos(pos_input[0], pos_input[1], pos_input[2]);
377 
379  TVector3 DJ(1., 0., 0.);
380  TVector3 DK(0., 1., 0.);
381  TVector3 DI = DJ.Cross(DK);
382 
384  double p_inv = 1./mom.Mag();
385  double vp = mom.Dot(DJ)/mom.Dot(DI);
386  double wp = mom.Dot(DK)/mom.Dot(DI);
387  double lambda = TMath::Pi()/2. - mom.Theta();
388  //double phi = mom.Phi();
389 
390  TVector3 TVW;
391  TVW.SetX(1./sqrt(1. + vp*vp + wp*wp));
392  TVW.SetY(vp*TVW.X());
393  TVW.SetZ(wp*TVW.X());
394 
395  TVector3 TN;
396  for(int i = 0; i < 3; i++)
397  {
398  TN[i] = TVW[0]*DI[i] + TVW[1]*DJ[i] + TVW[2]*DK[i];
399  }
400 
401  double cosl = cos(lambda);
402  double cosl1 = 1./cosl;
403 
404  TVector3 UN(-TN.Y()*cosl1, TN.X()*cosl1, 0.);
405  TVector3 VN(-TN.Z()*UN.Y(), TN.Z()*UN.X(), cosl);
406 
407  double UJ = UN.Dot(DJ);
408  double UK = UN.Dot(DK);
409  double VJ = VN.Dot(DJ);
410  double VK = VN.Dot(DK);
411 
412  jac_sd2sc.Zero();
413  jac_sd2sc[0][0] = 1.;
414  jac_sd2sc[1][1] = TVW[0]*VJ;
415  jac_sd2sc[1][2] = TVW[0]*VK;
416  jac_sd2sc[2][1] = TVW[0]*UJ*cosl1;
417  jac_sd2sc[2][2] = TVW[0]*UK*cosl1;
418  jac_sd2sc[3][3] = UJ;
419  jac_sd2sc[3][4] = UK;
420  jac_sd2sc[4][3] = VJ;
421  jac_sd2sc[4][4] = VK;
422 
423  //Takes cm output kGauss (0.1*tesla)
424  genfit::AbsBField *field = genfit::FieldManager::getInstance()->getField();
425 
426  if(charge != 0 && field)
427  {
428  TVector3 H = field->get(pos);
429  H *= kilogauss/tesla;
430 #ifdef _DEBUG_ON
431  LogInfo("");
432  pos.Print();
433  H.Print();
434 #endif
435  double HA = H.Mag();
436  double HAM = HA*p_inv*tesla*GeV;
437  double HM;
438  if(HA < 1E-6)
439  {
440  HM = 0.;
441  }
442  else
443  {
444  HM = charge/HA;
445  }
446 
447  double Q = -HAM*c_light/(km/ns);
448  double sinz = -H.Dot(UN)*HM;
449  double cosz = H.Dot(VN)*HM;
450 
451  jac_sd2sc[1][3] = -Q*TVW[1]*sinz;
452  jac_sd2sc[1][4] = -Q*TVW[2]*sinz;
453  jac_sd2sc[2][3] = -Q*TVW[1]*cosz*cosl1;
454  jac_sd2sc[2][4] = -Q*TVW[2]*cosz*cosl1;
455  }
456 }
457 
458 void GenFitExtrapolator::TRSCSD(int charge, TVector3 mom_input, TVector3 pos_input)
459 {
461  TVector3 mom(mom_input[0], mom_input[1], mom_input[2]);
462  TVector3 pos(pos_input[0], pos_input[1], pos_input[2]);
463 
465  TVector3 DJ(1., 0., 0.);
466  TVector3 DK(0., 1., 0.);
467  TVector3 DI = DJ.Cross(DK);
468 
470  double p_inv = 1./mom.Mag();
471  //double vp = mom.Dot(DJ)/mom.Dot(DI);
472  //double wp = mom.Dot(DK)/mom.Dot(DI);
473  double lambda = TMath::Pi()/2. - mom.Theta();
474  double phi = mom.Phi();
475 
476  double cosl = cos(lambda);
477  double sinp = sin(phi);
478  double cosp = cos(phi);
479 
480  TVector3 TN(cosl*cosp, cosl*sinp, sin(lambda));
481  TVector3 TVW;
482  TVW.SetX(TN.Dot(DI));
483  TVW.SetY(TN.Dot(DJ));
484  TVW.SetZ(TN.Dot(DK));
485 
486  double T1R = 1./TVW[0];
487  double T2R = T1R*T1R;
488  TVector3 UN(-sinp, cosp, 0.);
489  TVector3 VN(-TN.Z()*UN.Y(), TN.Z()*UN.X(), cosl);
490 
491  double UJ = UN.Dot(DJ);
492  double UK = UN.Dot(DK);
493  double VJ = VN.Dot(DJ);
494  double VK = VN.Dot(DK);
495  double UI = UN.Dot(DI);
496  double VI = VN.Dot(DI);
497 
498  jac_sc2sd.Zero();
499  jac_sc2sd[0][0] = 1.;
500  jac_sc2sd[1][1] = -UK*T2R;
501  jac_sc2sd[1][2] = VK*cosl*T2R;
502  jac_sc2sd[2][1] = UJ*T2R;
503  jac_sc2sd[2][2] = -VJ*cosl*T2R;
504  jac_sc2sd[3][3] = VK*T1R;
505  jac_sc2sd[3][4] = -UK*T1R;
506  jac_sc2sd[4][3] = -VJ*T1R;
507  jac_sc2sd[4][4] = UJ*T1R;
508 
509  genfit::AbsBField *field = genfit::FieldManager::getInstance()->getField();
510  if(charge != 0 && field)
511  {
512  TVector3 H = field->get(pos);
513  H *= kilogauss/tesla;
514 
515  double HA = H.Mag();
516  double HAM = HA*p_inv;
517  double HM;
518  if(HA < 1E-6)
519  {
520  HM = 0.;
521  }
522  else
523  {
524  HM = charge/HA;
525  }
526 
527  double Q = -HAM*c_light/(km/ns);
528  double sinz = -H.Dot(UN)*HM;
529  double cosz = H.Dot(VN)*HM;
530  double T3R = Q*T1R*T1R*T1R;
531 
532  jac_sc2sd[1][3] = -UI*(VK*cosz - UK*sinz)*T3R;
533  jac_sc2sd[1][4] = -VI*(VK*cosz - UK*sinz)*T3R;
534  jac_sc2sd[2][3] = UI*(VJ*cosz - UJ*sinz)*T3R;
535  jac_sc2sd[2][4] = VI*(VJ*cosz - UJ*sinz)*T3R;
536  }
537 }
538 
539 void GenFitExtrapolator::print() {
540  cout << "Propagating " << iParType << ":" << endl;
541  cout << "From "
542  << "(" << pos_i.X() << ", " << pos_i.Y() << ", "<< pos_i.Z() << ") "
543  << " To "
544  << "(" << pos_f.X() << ", " << pos_f.Y() << ", "<< pos_f.Z() << ") "
545  << endl;
546 
547  cout << "Momentum change: From "
548  << "(" << mom_i.X() << ", " << mom_i.Y() << ", "<< mom_i.Z() << ") "
549  << " To "
550  << "(" << mom_f.X() << ", " << mom_f.Y() << ", "<< mom_f.Z() << ") "
551  << endl;
552 
553  cout << "Initial error matrix: " << endl;
554  for(int i = 0; i < 5; i++)
555  {
556  for(int j = 0; j < 5; j++)
557  {
558  cout << cov_i[i][j] << " ";
559  }
560  cout << endl;
561  }
562 
563  cout << "Final error matrix: " << endl;
564  for(int i = 0; i < 5; i++)
565  {
566  for(int j = 0; j < 5; j++)
567  {
568  cout << cov_f[i][j] << " ";
569  }
570  cout << endl;
571  }
572 }
GenFitExtrapolator::extrapolateTo
bool extrapolateTo(double z_out)
Definition: GenFitExtrapolator.cxx:234
GenFitExtrapolator::getPropagator
void getPropagator(TMatrixD &prop)
Get the propagator.
Definition: GenFitExtrapolator.cxx:206
DPGEN::s
const double s
Definition: SQPrimaryParticleGen.C:58
c_light
static const double c_light
Definition: GenFitExtrapolator.cxx:59
GenFitExtrapolator::TRSCSD
void TRSCSD(int charge, TVector3 mom_input, TVector3 pos_input)
Definition: GenFitExtrapolator.cxx:458
PHFlag::get_DoubleFlag
virtual double get_DoubleFlag(const std::string &name) const
Definition: PHFlag.cc:49
PHField
transient DST object for field storage and access
Definition: PHField.h:13
recoConsts::instance
static recoConsts * instance()
Definition: recoConsts.cc:7
GenFitExtrapolator::setParticleType
void setParticleType(int type)
Set particle type.
Definition: GenFitExtrapolator.cxx:165
GenFitExtrapolator::init
bool init(const PHField *field, const TGeoManager *geom)
Initialize geometry and physics.
Definition: GenFitExtrapolator.cxx:74
GenFitExtrapolator::getFinalStateWithCov
void getFinalStateWithCov(TMatrixD &state_out, TMatrixD &cov_out)
Get the final state parameters and covariance.
Definition: GenFitExtrapolator.cxx:169
GenFitExtrapolator::setInitialStateWithCov
void setInitialStateWithCov(double z_in, TMatrixD &state_in, TMatrixD &cov_in)
Set input initial state parameters.
Definition: GenFitExtrapolator.cxx:116
GenFitExtrapolator::convertSVtoMP
void convertSVtoMP(double z, TMatrixD &state, TVector3 &mom, TVector3 &pos)
Transformation between the state vector and the mom/pos.
Definition: GenFitExtrapolator.cxx:346
GenFitExtrapolator::convertMPtoSV
void convertMPtoSV(TVector3 &mom, TVector3 &pos, TMatrixD &state)
Definition: GenFitExtrapolator.cxx:360
GenFitExtrapolator::print
void print()
Debug print.
Definition: GenFitExtrapolator.cxx:539
PHField::GetFieldValue
virtual void GetFieldValue(const double Point[4], double *Bfield) const =0
GenFitExtrapolator::TRSDSC
void TRSDSC(int charge, TVector3 mom_input, TVector3 pos_input)
Definition: GenFitExtrapolator.cxx:372
LogInfo
#define LogInfo(message)
Definition: DbSvc.cc:14
PHFlag::get_IntFlag
virtual int get_IntFlag(const std::string &name) const
Definition: PHFlag.cc:117
GenFitExtrapolator::extrapolateToIP
double extrapolateToIP()
Extrapolate to the primary vertex.
Definition: GenFitExtrapolator.cxx:292