Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHGenFitTrackProjection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHGenFitTrackProjection.cc
1 
7 #include "PHGenFitTrackProjection.h"
8 
9 #include <trackbase_historic/SvtxTrackMap.h>
10 #include <trackbase_historic/SvtxTrack.h>
11 #include <trackbase_historic/SvtxTrackState.h> // for SvtxTrackState
12 
13 #include <calobase/RawTowerGeomContainer.h>
14 #include <calobase/TowerInfoContainerv1.h>
15 #include <calobase/TowerInfo.h>
16 #include <calobase/RawClusterContainer.h>
17 #include <calobase/RawCluster.h>
18 #include <calobase/RawClusterUtility.h>
19 
20 #include <phgenfit/Fitter.h>
21 
22 #include <phgeom/PHGeomUtility.h>
23 
24 #include <phfield/PHFieldUtility.h>
25 
26 #include <fun4all/Fun4AllReturnCodes.h>
27 #include <fun4all/SubsysReco.h> // for SubsysReco
28 
29 #include <phool/getClass.h>
30 #include <phool/phool.h> // for PHWHERE
31 
32 #include <GenFit/AbsTrackRep.h> // for AbsTrackRep
33 #include <GenFit/MeasuredStateOnPlane.h> // for MeasuredStateOnPlane
34 #include <GenFit/RKTrackRep.h>
35 
36 
37 //ROOT
38 #include <TDatabasePDG.h> // for TDatabasePDG
39 #include <TMatrixDSymfwd.h> // for TMatrixDSym
40 #include <TMatrixTSym.h> // for TMatrixTSym
41 #include <TMatrixTUtils.h> // for TMatrixTRow
42 #include <TParticlePDG.h> // for TParticlePDG
43 #include <TVector3.h> // for TVector3
44 
45 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
46 
47 // standard includes
48 #include <cfloat> // for DBL_MAX
49 #include <cmath> // for isnan, atan2, sqrt, NAN
50 #include <cstdlib> // for abs
51 #include <iostream>
52 #include <map> // for _Rb_tree_iterator
53 #include <memory>
54 #include <utility> // for pair
55 #include <vector>
56 
57 class PHField;
58 class TGeoManager;
59 
60 //#define DEBUG
61 
62 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
63 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
64 #define LogWarning(exp) std::cout<<"WARNING: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<std::endl
65 
66 using namespace std;
67 
68 PHGenFitTrackProjection::PHGenFitTrackProjection(const string &name, const int pid_guess) :
69  SubsysReco(name),
70 
71  _fitter(nullptr),
72 
73  _pid_guess(pid_guess),
74 
75  _num_cal_layers(4)
76 {
77  _cal_radii.assign(_num_cal_layers, NAN);
78  _cal_names.push_back("PRES"); // PRES not yet in G4
79  _cal_names.push_back("CEMC");
80  _cal_names.push_back("HCALIN");
81  _cal_names.push_back("HCALOUT");
82  _cal_types.push_back(SvtxTrack::PRES); // PRES not yet in G4
83  _cal_types.push_back(SvtxTrack::CEMC);
84  _cal_types.push_back(SvtxTrack::HCALIN);
85  _cal_types.push_back(SvtxTrack::HCALOUT);
86 }
87 
88 int PHGenFitTrackProjection::Init(PHCompositeNode */*topNode*/) {
89  return Fun4AllReturnCodes::EVENT_OK;
90 }
91 
92 int PHGenFitTrackProjection::InitRun(PHCompositeNode *topNode) {
93  for (int i = 0; i < _num_cal_layers; ++i) {
94  string nodename = "TOWERGEOM_" + _cal_names[i];
95  RawTowerGeomContainer *geo = findNode::getClass<RawTowerGeomContainer>(
96  topNode, nodename.c_str());
97  if (geo)
98  _cal_radii[i] = geo->get_radius();//+0.5*geo->get_thickness();
99  }
100 
101  TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
102 
103 #ifdef DEBUG
104  tgeo_manager->Export("Geo_extract.root");
105 #endif
106 
107  PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
108  _fitter = PHGenFit::Fitter::getInstance(tgeo_manager,field,
109  "DafRef",
110  "RKTrackRep", false);
111 
112 
113  if (!_fitter) {
114  cerr << PHWHERE << endl;
115  return Fun4AllReturnCodes::ABORTRUN;
116  }
117 
118  _fitter->set_verbosity(Verbosity());
119 
120  if (Verbosity() > 0) {
121  cout
122  << "================== PHGenFitTrackProjection::InitRun() ====================="
123  << endl;
124  for (int i = 0; i < _num_cal_layers; ++i) {
125  if (!std::isnan(_cal_radii[i])) {
126  cout << " " << _cal_names[i] << " projection radius: "
127  << _cal_radii[i] << " cm" << endl;
128  }
129  }
130  cout << " projections still curl after the mag field" << endl;
131  cout
132  << " projections start from the vertex momentum vector (M.S. effects could be large)"
133  << endl;
134  cout << " projections don't correct for the slat HCAL geometry" << endl;
135  cout
136  << "==========================================================================="
137  << endl;
138  }
139 
140  return Fun4AllReturnCodes::EVENT_OK;
141 }
142 
143 int PHGenFitTrackProjection::process_event(PHCompositeNode *topNode) {
144  if (Verbosity() > 1)
145  cout << "PHGenFitTrackProjection::process_event -- entered" << endl;
146 
147  //---------------------------------
148  // Get Objects off of the Node Tree
149  //---------------------------------
150 
151  // Pull the reconstructed track information off the node tree...
152  SvtxTrackMap* _g4tracks = findNode::getClass<SvtxTrackMap>(topNode,
153  "SvtxTrackMap");
154  if (!_g4tracks) {
155  cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap." << endl;
156  return Fun4AllReturnCodes::ABORTRUN;
157  }
158 
159  for (int i = 0; i < _num_cal_layers; ++i) {
160 
161  if (std::isnan(_cal_radii[i]))
162  continue;
163 
164  if (Verbosity() > 1)
165  cout << "Projecting tracks into: " << _cal_names[i] << endl;
166 
167  // pull the tower geometry
168  string towergeonodename = "TOWERGEOM_" + _cal_names[i];
169  RawTowerGeomContainer *towergeo = findNode::getClass<
170  RawTowerGeomContainer>(topNode, towergeonodename.c_str());
171  if (!towergeo) {
172  cerr << PHWHERE << " ERROR: Can't find node " << towergeonodename
173  << endl;
174  return Fun4AllReturnCodes::ABORTRUN;
175  }
176 
177  // pull the towers
178  string towernodename = "TOWER_CALIB_" + _cal_names[i];
179  TowerInfoContainer *towerList = findNode::getClass<TowerInfoContainerv1>(
180  topNode, towernodename.c_str());
181 
182  if (!towerList) {
183  cerr << PHWHERE << " ERROR: Can't find node " << towernodename
184  << endl;
185  return Fun4AllReturnCodes::ABORTRUN;
186  }
187 
188  // pull the clusters
189  string clusternodename = "CLUSTER_" + _cal_names[i];
190  RawClusterContainer *clusterList = findNode::getClass<
191  RawClusterContainer>(topNode, clusternodename.c_str());
192 
193  if(_use_poscalib_cemc and
194  _cal_names[i].compare("CEMC") == 0) {
195  std::string nodeName = "CLUSTER_POS_COR_" + _cal_names[i];
196  clusterList = findNode::getClass<RawClusterContainer>( topNode, nodeName.c_str());
197 
198  if(Verbosity() > 1)
199  std::cout << "Grabbing CEMC position recalib clusters"
200  << std::endl;
201  }
202 
203  if (!clusterList) {
204  cerr << PHWHERE << " ERROR: Can't find node " << clusternodename
205  << endl;
206  return Fun4AllReturnCodes::ABORTRUN;
207  }
208 
209  // loop over all tracks
210  for (SvtxTrackMap::Iter iter = _g4tracks->begin();
211  iter != _g4tracks->end(); ++iter) {
212  SvtxTrack *track = iter->second;
213  if(!track) {
214  if(Verbosity() >= 2) LogWarning("!track");
215  continue;
216  }
217 
218 #ifdef DEBUG
219  cout
220  <<__LINE__
221  <<": track->get_charge(): "<<track->get_charge()
222  <<endl;
223 #endif
224  if (Verbosity() > 1)
225  cout << "projecting track id " << track->get_id() << endl;
226 
227  if (Verbosity() > 1) {
228  cout << " track pt = " << track->get_pt() << endl;
229  }
230 
231  std::vector<double> point;
232  point.assign(3, -9999.);
233 
234  auto last_state_iter = --track->end_states();
235 
236  SvtxTrackState * trackstate = last_state_iter->second;
237 
238  if(!trackstate) {
239  if(Verbosity() >= 2) LogWarning("!trackstate");
240  continue;
241  }
242 
243  auto pdg = unique_ptr<TDatabasePDG> (TDatabasePDG::Instance());
244  int reco_charge = track->get_charge();
245  int gues_charge = pdg->GetParticle(_pid_guess)->Charge();
246  if(reco_charge*gues_charge<0) _pid_guess *= -1;
247 #ifdef DEBUG
248  cout
249  <<__LINE__
250  <<": guess charge: " << gues_charge
251  <<": reco charge: " << reco_charge
252  <<": pid: " << _pid_guess
253  <<": pT: " << sqrt(trackstate->get_px()*trackstate->get_px() + trackstate->get_py()*trackstate->get_py())
254  <<endl;
255 #endif
256 
257  auto rep = unique_ptr<genfit::AbsTrackRep> (new genfit::RKTrackRep(_pid_guess));
258 
259  unique_ptr<genfit::MeasuredStateOnPlane> msop80 = nullptr;
260 
261  {
262  TVector3 pos(trackstate->get_x(), trackstate->get_y(), trackstate->get_z());
263  //pos.SetXYZ(0.01,0,0);
264 
265  TVector3 mom(trackstate->get_px(), trackstate->get_py(), trackstate->get_pz());
266  //mom.SetXYZ(1,0,0);
267 
268  TMatrixDSym cov(6);
269  for (int k = 0; k < 6; ++k) {
270  for (int j = 0; j < 6; ++j) {
271  cov[k][j] = trackstate->get_error(k, j);
272  }
273  }
274 
275  msop80 = unique_ptr<genfit::MeasuredStateOnPlane> (new genfit::MeasuredStateOnPlane(rep.get()));
276 
277  msop80->setPosMomCov(pos, mom, cov);
278  }
279 
280 #ifdef DEBUG
281  {
282  double x = msop80->getPos().X();
283  double y = msop80->getPos().Y();
284  double z = msop80->getPos().Z();
285 // double px = msop80->getMom().X();
286 // double py = msop80->getMom().Y();
287  double pz = msop80->getMom().Z();
288  genfit::FieldManager *field_mgr = genfit::FieldManager::getInstance();
289  double Bx=0, By=0, Bz=0;
290  field_mgr->getFieldVal(x,y,z,Bx,By,Bz);
291  cout
292  << __LINE__
293  << ": { " << msop80->getPos().Perp() << ", " << msop80->getPos().Phi() << ", " << msop80->getPos().Eta() << "} @ "
294  //<< "{ " << Bx << ", " << By << ", " << Bz << "}"
295  << "{ " << msop80->getMom().Perp() << ", " << msop80->getMom().Phi() << ", " << pz << "} "
296  <<endl;
297  //msop80->Print();
298  }
299 #endif
300  try {
301  rep->extrapolateToCylinder(*msop80, _cal_radii[i], TVector3(0,0,0), TVector3(0,0,1));
302  //rep->extrapolateToCylinder(*msop80, 5., TVector3(0,0,0), TVector3(0,0,1));
303  } catch (...) {
304  if(Verbosity() >= 2) LogWarning("extrapolateToCylinder failed");
305  continue;
306  }
307 
308 #ifdef DEBUG
309  {
310  cout<<__LINE__<<endl;
311  //msop80->Print();
312  double x = msop80->getPos().X();
313  double y = msop80->getPos().Y();
314  double z = msop80->getPos().Z();
315 // double px = msop80->getMom().X();
316 // double py = msop80->getMom().Y();
317  double pz = msop80->getMom().Z();
318  genfit::FieldManager *field_mgr = genfit::FieldManager::getInstance();
319  double Bx=0, By=0, Bz=0;
320  field_mgr->getFieldVal(x,y,z,Bx,By,Bz);
321  cout
322  << __LINE__
323  << ": { " << msop80->getPos().Perp() << ", " << msop80->getPos().Phi() << ", " << msop80->getPos().Eta() << "} @ "
324  //<< "{ " << Bx << ", " << By << ", " << Bz << "}"
325  << "{ " << msop80->getMom().Perp() << ", " << msop80->getMom().Phi() << ", " << pz << "} "
326  <<endl;
327  }
328 #endif
329 
330  point[0] = msop80->getPos().X();
331  point[1] = msop80->getPos().Y();
332  point[2] = msop80->getPos().Z();
333 
334 #ifdef DEBUG
335  cout
336  <<__LINE__
337  <<": GenFit: {"
338  << point[0] <<", "
339  << point[1] <<", "
340  << point[2] <<" }"
341  <<endl;
342 #endif
343 
344  if (std::isnan(point[0]))
345  continue;
346  if (std::isnan(point[1]))
347  continue;
348  if (std::isnan(point[2]))
349  continue;
350 
351  double x = point[0];
352  double y = point[1];
353  double z = point[2];
354 
355  double phi = atan2(y, x);
356  double eta = asinh(z / sqrt(x * x + y * y));
357 
358  if (Verbosity() > 1) {
359  cout << " initial track phi = " << track->get_phi();
360  cout << ", eta = " << track->get_eta() << endl;
361  cout << " calorimeter phi = " << phi << ", eta = " << eta
362  << endl;
363  }
364 
365  // projection is outside the detector extent
366  // TODO towergeo doesn't make this easy to extract, but this should be
367  // fetched from the node tree instead of hardcoded
368  if (fabs(eta) >= 1.0)
369  continue;
370 
371  // calculate 3x3 tower energy
372  int binphi = towergeo->get_phibin(phi);
373  int bineta = towergeo->get_etabin(eta);
374 
375  double energy_3x3 = 0.0;
376  double energy_5x5 = 0.0;
377  for (int iphi = binphi - 2; iphi <= binphi + 2; ++iphi) {
378  for (int ieta = bineta - 2; ieta <= bineta + 2; ++ieta) {
379 
380  // wrap around
381  int wrapphi = iphi;
382  if (wrapphi < 0) {
383  wrapphi = towergeo->get_phibins() + wrapphi;
384  }
385  if (wrapphi >= towergeo->get_phibins()) {
386  wrapphi = wrapphi - towergeo->get_phibins();
387  }
388 
389  // edges
390  if (ieta < 0)
391  continue;
392  if (ieta >= towergeo->get_etabins())
393  continue;
394  unsigned int towerkey = (ieta << 16U) + wrapphi;
395  TowerInfo* tower = towerList->get_tower_at_key(towerkey);
396  if (tower) {
397 
398  energy_5x5 += tower->get_energy();
399  if (abs(iphi - binphi) <= 1 and abs(ieta - bineta) <= 1)
400  energy_3x3 += tower->get_energy();
401 
402  if (Verbosity() > 1)
403  cout << " tower " << ieta << " " << wrapphi
404  << " energy = " << tower->get_energy()
405  << endl;
406  }
407  }
408  }
409 
410  track->set_cal_energy_3x3(_cal_types[i], energy_3x3);
411  track->set_cal_energy_5x5(_cal_types[i], energy_5x5);
412 
413  // loop over all clusters and find nearest
414  double min_r = DBL_MAX;
415  double min_index = -9999;
416  double min_dphi = NAN;
417  double min_deta = NAN;
418  double min_e = NAN;
419 #ifdef DEBUG
420  double min_cluster_phi = NAN;
421 #endif
422  for (const auto & iterator : clusterList->getClustersMap()) {
423 
424  const RawCluster *cluster = iterator.second;
425 
427  const float cluster_eta = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0,0,0));
428 
429  double dphi = atan2(sin(phi - cluster->get_phi()),
430  cos(phi - cluster->get_phi()));
431  double deta = eta - cluster_eta;
432  double r = sqrt(pow(dphi, 2) + pow(deta, 2));
433 
434  if (r < min_r) {
435  min_index = iterator.first;
436  min_r = r;
437  min_dphi = dphi;
438  min_deta = deta;
439  min_e = cluster->get_energy();
440 #ifdef DEBUG
441  min_cluster_phi = cluster->get_phi();
442 #endif
443  }
444  }
445 
446  if (min_index != -9999) {
447  track->set_cal_dphi(_cal_types[i], min_dphi);
448  track->set_cal_deta(_cal_types[i], min_deta);
449  track->set_cal_cluster_id(_cal_types[i], min_index);
450  track->set_cal_cluster_e(_cal_types[i], min_e);
451 
452 #ifdef DEBUG
453  cout
454  <<__LINE__
455  <<": min_cluster_phi: "<<min_cluster_phi
456  <<endl;
457 #endif
458 
459  if (Verbosity() > 1) {
460  cout << " nearest cluster dphi = " << min_dphi << " deta = "
461  << min_deta << " e = " << min_e << endl;
462  }
463  }
464 
465  } // end track loop
466  } // end calorimeter layer loop
467 
468  if (Verbosity() > 1)
469  cout << "PHGenFitTrackProjection::process_event -- exited" << endl;
470 
471  return Fun4AllReturnCodes::EVENT_OK;
472 }
473 
474 int PHGenFitTrackProjection::End(PHCompositeNode */*topNode*/) {
475  return Fun4AllReturnCodes::EVENT_OK;
476 }