Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnaMvtxTelescopeHits.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnaMvtxTelescopeHits.C
1 
9 #include "AnaMvtxTelescopeHits.h"
10 
11 #include <phool/getClass.h>
12 #include <phool/phool.h>
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHIODataNode.h>
15 #include <phool/PHNodeIterator.h>
16 
17 #include <phgeom/PHGeomUtility.h>
18 
19 #include <fun4all/Fun4AllServer.h>
20 #include <fun4all/Fun4AllReturnCodes.h>
21 #include <fun4all/PHTFileServer.h>
22 
23 #include <g4main/PHG4TruthInfoContainer.h>
24 #include <g4main/PHG4Particle.h>
25 #include <g4main/PHG4Hit.h>
26 #include <g4main/PHG4HitContainer.h>
27 
28 #include <g4hough/SvtxCluster.h>
29 #include <g4hough/SvtxClusterMap.h>
30 #include <g4hough/SvtxHitMap.h>
31 #include <g4hough/SvtxTrack.h>
32 #include <g4hough/SvtxVertex.h>
33 #include <g4hough/SvtxTrackMap.h>
34 #include <g4hough/SvtxVertexMap.h>
35 
36 #include <g4jets/JetMap.h>
37 #include <g4jets/Jet.h>
38 
39 #include <g4detectors/PHG4CellContainer.h>
40 #include <g4detectors/PHG4CylinderGeomContainer.h>
41 #include <g4detectors/PHG4Cell.h>
42 #include <g4detectors/PHG4CylinderGeom_MAPS.h>
43 #include <g4detectors/PHG4CylinderGeom_Siladders.h>
44 
45 #include <phgenfit/Fitter.h>
46 #include <phgenfit/PlanarMeasurement.h>
47 #include <phgenfit/Track.h>
48 #include <phgenfit/SpacepointMeasurement.h>
49 
50 //GenFit
51 #include <GenFit/FieldManager.h>
52 #include <GenFit/GFRaveConverters.h>
53 #include <GenFit/GFRaveVertex.h>
54 #include <GenFit/GFRaveVertexFactory.h>
55 #include <GenFit/MeasuredStateOnPlane.h>
56 #include <GenFit/RKTrackRep.h>
57 #include <GenFit/StateOnPlane.h>
58 #include <GenFit/Track.h>
59 
60 //Rave
61 #include <rave/Version.h>
62 #include <rave/Track.h>
63 #include <rave/VertexFactory.h>
64 #include <rave/ConstantMagneticField.h>
65 
66 #include <phhepmc/PHHepMCGenEvent.h>
67 #include <HepMC/GenEvent.h>
68 #include <HepMC/GenVertex.h>
69 
70 #include <HFJetTruthGeneration/HFJetDefs.h>
71 
72 #include "TTree.h"
73 #include "TFile.h"
74 #include "TLorentzVector.h"
75 #include "TH1.h"
76 #include "TH2.h"
77 #include "TVectorD.h"
78 #include "TMatrixD.h"
79 #include "TDecompSVD.h"
80 
81 
82 #include <iostream>
83 #include <map>
84 #include <utility>
85 #include <vector>
86 #include <memory>
87 
88 
89 
90 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
91 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
92 #define LogWarning(exp) std::cout<<"WARNING: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
93 
94 
95 AnaMvtxTelescopeHits::AnaMvtxTelescopeHits(const std::string &name,
96  const std::string &ofName) :
97  SubsysReco(name),
98  _clustermap(NULL),
99  _foutname(ofName),
100  _f(NULL)
101 {
102 
103  hlayer = NULL;
104  for (int il = 0; il < 4; il++)
105  {
106  hsize_phi[il] = NULL;
107  hsize_z[il] = NULL;
108  hphiz[il] = NULL;
109  hdphi[il] = NULL;
110  hdz[il] = NULL;
111  }
112 }
113 
114 int AnaMvtxTelescopeHits::Init(PHCompositeNode *topNode)
115 {
116 
117  //-- Open file
118  _f = new TFile(_foutname.c_str(), "RECREATE");
119 
120  //-- Initialize histograms
121  hlayer = new TH1D("hlayer", ";layer", 6, -0.5, 5.5);
122 
123  for (int il = 0; il < 4; il++)
124  {
125  hsize_phi[il] = new TH1D(Form("hsize_phi_l%i", il),
126  ";cluster size #phi [cm]",
127  100, 0, 0.1);
128 
129  hsize_z[il] = new TH1D(Form("hsize_z_l%i", il),
130  ";cluster size z [cm]",
131  100, 0, 0.1);
132 
133  hphiz[il] = new TH2D(Form("hphiz_l%i", il),
134  ";z [cm]; #phi",
135  1000, -10, 10,
136  100, -0.1, 0.1);
137 
138  hdphi[il] = new TH1D(Form("hdphi_l%i", il),
139  "; d#phi [#mum]",
140  2000, -100, 100);
141 
142  hdz[il] = new TH1D(Form("hdz_l%i", il),
143  "; dz [#mum]",
144  2000, -100, 100);
145  } // il
146 
147 
148 
149  return 0;
150 
151 }
152 
153 int AnaMvtxTelescopeHits::InitRun(PHCompositeNode *topNode)
154 {
155 
156  // TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
157 
158  // _fitter = PHGenFit::Fitter::getInstance(tgeo_manager, _mag_field_file_name.data(),
159  // (_reverse_mag_field) ? -1. * _mag_field_re_scaling_factor : _mag_field_re_scaling_factor,
160  // _track_fitting_alg_name, "RKTrackRep", _do_evt_display);
161 
162  // if (!_fitter) {
163  // std::cerr << PHWHERE << std::endl;
164  // return Fun4AllReturnCodes::ABORTRUN;
165  // }
166 
167  //-- set counters
168  _ievent = -1; // since we count at the beginning of the event, start at -1
169 
170 
171  return 0;
172 
173 }
174 
175 int AnaMvtxTelescopeHits::process_event(PHCompositeNode *topNode)
176 {
177 
178  if ( verbosity > VERBOSITY_MORE )
179  std::cout << "AnaMvtxTelescopeHits::process_event()" << std::endl;
180 
181  // count event
182  _ievent++;
183 
184  //------
185  //-- Get the nodes
186  //------
187  int retval = GetNodes(topNode);
188  if ( retval != Fun4AllReturnCodes::EVENT_OK )
189  return retval;
190 
191  if ( verbosity > VERBOSITY_MORE )
192  {
193  std::cout << "-- evnt:" << _ievent << std::endl;
194  std::cout << "-- Nclus:" << _clustermap->size() << std::endl;
195  }
196 
197 
198  //------
199  //-- Loop over clusters
200  //------
201 
202  if ( verbosity > VERBOSITY_MORE )
203  std::cout << "-- Looping over clusters" << std::endl;
204 
205  LyrClusMap _mmap_lyr_clus;
206 
207  for (SvtxClusterMap::Iter iter = _clustermap->begin();
208  iter != _clustermap->end();
209  ++iter)
210  {
211  SvtxCluster *clus = iter->second;
212 
213  hlayer->Fill(clus->get_layer());
214 
215  unsigned int lyr = clus->get_layer();
216  if (lyr < 0 || lyr > 3)
217  continue;
218  hsize_phi[lyr]->Fill(clus->get_phi_size());
219  hsize_z[lyr]->Fill(clus->get_z_size());
220 
221  double phi = TMath::ATan2(clus->get_y(), clus->get_x());
222  hphiz[lyr]->Fill(clus->get_z(), phi);
223 
224  // insert into map
225  _mmap_lyr_clus.insert(std::make_pair(lyr, clus));
226  }
227 
228  if ( verbosity > VERBOSITY_MORE )
229  std::cout << "-- _mmap_lyr_clus.size():" << _mmap_lyr_clus.size() << std::endl;
230 
231  //------
232  //-- Simple tracking
233  //------
234 
235  if ( verbosity > VERBOSITY_MORE )
236  std::cout << "-- Performing simple tracking" << std::endl;
237 
238  // loop over clusters in layer 0
239  LyrClusMap::iterator lyr0_itr;
240  LyrClusMap::iterator lyr_itr;
241  LyrClusMap::iterator lyr3_itr;
242  for ( lyr0_itr = _mmap_lyr_clus.lower_bound(0);
243  lyr0_itr != _mmap_lyr_clus.upper_bound(0);
244  ++lyr0_itr)
245  {
246  SvtxCluster* clus0 = lyr0_itr->second;
247 
248  if ( !clus0 )
249  {
250  std::cout << "WARNING: clus0 not found" << std::endl;
251  continue;
252  }
253 
254  // get (r, phi, z)
255  double r0 = TMath::Sqrt(TMath::Power(clus0->get_x(), 2) + TMath::Power(clus0->get_y(), 2));
256  double p0 = TMath::ATan2(clus0->get_y(), clus0->get_x());
257  double z0 = clus0->get_z();
258 
259  if ( verbosity > VERBOSITY_MORE )
260  {
261  std::cout << "-- clus0:(" << r0 << ", " << p0 << ", " << z0 << ")" << std::endl;
262  }
263 
264  // loop over all clusters in the outer layer
265  for ( lyr3_itr = _mmap_lyr_clus.lower_bound(3);
266  lyr3_itr != _mmap_lyr_clus.upper_bound(3);
267  ++lyr3_itr)
268  {
269  SvtxCluster* clus3 = lyr3_itr->second;
270 
271  if ( !clus3 )
272  {
273  std::cout << "WARNING: clus3 not found" << std::endl;
274  continue;
275  }
276 
277  // get (r, phi, z)
278  double r3 = TMath::Sqrt(TMath::Power(clus3->get_x(), 2) + TMath::Power(clus3->get_y(), 2));
279  double p3 = TMath::ATan2(clus3->get_y(), clus3->get_x());
280  double z3 = clus3->get_z();
281 
282  // calc m, b in (phi, r) space
283  double mpr = CalcSlope(r0, p0, r3, p3);
284  double bpr = CalcIntecept(r0, p0, mpr);
285 
286  // calc m, b in (z, r) space
287  double mzr = CalcSlope(r0, z0, r3, z3);
288  double bzr = CalcIntecept(r0, z0, mzr);
289 
290 
291  if ( verbosity > VERBOSITY_MORE )
292  {
293  std::cout << "-- clus3:(" << r3 << ", " << p3 << ", " << z3 << ")" << std::endl;
294  std::cout << "-- mpr:" << mpr << " bpr:" << bpr << std::endl;
295  std::cout << "-- mzr:" << mzr << " bzr:" << bzr << std::endl;
296  }
297 
298 
299 
300  // loop over all clusters in layer 1 & 2 and calc displacement
301  for (int ilyr = 1; ilyr <= 2; ilyr++)
302  {
303  for ( lyr_itr = _mmap_lyr_clus.lower_bound(ilyr);
304  lyr_itr != _mmap_lyr_clus.upper_bound(ilyr);
305  ++lyr_itr)
306  {
307  SvtxCluster* clus = lyr_itr->second;
308 
309  if ( !clus )
310  continue;
311 
312  // get (r, phi, z)
313  double r = TMath::Sqrt(TMath::Power(clus->get_x(), 2) + TMath::Power(clus->get_y(), 2));
314  double phi = TMath::ATan2(clus->get_y(), clus->get_x());
315  double z = clus->get_z();
316 
317  // calc track projection in phi, z
318  double phi_proj = CalcProjection(r, mpr, bpr);
319  double z_proj = CalcProjection(r, mzr, bzr);
320 
321  // calc proj - true
322  double dphi = (phi_proj - phi) * 1e4;
323  double dz = (z_proj - z) * 1e4;
324 
325  if ( !hdphi[ilyr] || !hdz[ilyr] )
326  {
327  std::cout << "WARNING!! can't find hdphi or hdz for lyr:" << ilyr << std::endl;
328  }
329  else
330  {
331  hdphi[ilyr]->Fill(dphi);
332  hdz[ilyr]->Fill(dz);
333  }
334 
335  if ( verbosity > VERBOSITY_MORE )
336  {
337  std::cout << "------------------" << std::endl;
338  std::cout << "-- clus" << clus->get_layer() << ":(" << r3 << ", " << p3 << ", " << z3 << ")" << std::endl;
339  std::cout << "-- proj:(" << phi_proj << ", " << z_proj << ")" << std::endl;
340  std::cout << "-- dphi [micron]: " << dphi << std::endl;
341  std::cout << "-- dz [micron]: " << dz << std::endl;
342  std::cout << "------------------" << std::endl;
343  }
344 
345  } // lyr1_itr
346  } // ilyr;
347 
348  } // lyr3_itr
349 
350  } //lyr0_itr
351 
352 
353  if ( verbosity > VERBOSITY_MORE )
354  std::cout << "-- Done simple tracking" << std::endl;
355 
356 
357 
358  // //------
359  // //-- Track candidates using all layers
360  // //------
361 
362  // LyrClusMap::iterator lyr_itr;
363  // TrkVec _vec_trk;
364  // ClusVec _trk_candidate;
365 
366  // for ( lyr_itr = _mmap_lyr_clus.lower_bound(0);
367  // lyr_itr != _mmap_lyr_clus.upper_bound(0);
368  // ++lyr_itr)
369  // {
370  // SvtxCluster* clus0 = lyr_itr->second;
371 
372  // }
373 
374 
375  //-- Cleanup
376 
377 
378  //-- End
379  return 0;
380 
381 }
382 
383 int AnaMvtxTelescopeHits::End(PHCompositeNode * topNode)
384 {
385  _f->Write();
386  _f->Close();
387 
388  return 0;
389 }
390 
391 /*
392  * Get all the necessary nodes from the node tree
393  */
394 int AnaMvtxTelescopeHits::GetNodes(PHCompositeNode * topNode)
395 {
396 
397  // Input Svtx Clusters
398  _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
399  if (!_clustermap && _ievent < 2) {
400  std::cout << PHWHERE << " SvtxClusterMap node not found on node tree"
401  << std::endl;
402  return Fun4AllReturnCodes::ABORTEVENT;
403  }
404 
405  return Fun4AllReturnCodes::EVENT_OK;
406 
407 }
408 
409 /*
410  * Simple helper function for calculating the slope between two points
411  */
412 double
413 AnaMvtxTelescopeHits::CalcSlope(double x0, double y0, double x1, double y1)
414 {
415  return (y1 - y0) / (x1 - x0);
416 }
417 
418 /*
419  * Simple helper function for calculating intercept
420  */
421 double
422 AnaMvtxTelescopeHits::CalcIntecept(double x0, double y0, double m)
423 {
424  return y0 - x0 * m;
425 }
426 
427 /*
428  * Simple helper function for calculating projection
429  */
430 double
431 AnaMvtxTelescopeHits::CalcProjection(double x, double m, double b)
432 {
433  return m * x + b;
434 }
435 
436 
437 // /*
438 // *
439 // */
440 // int AnaMvtxTelescopeHits::GetTrackCandidates(LyrClusMap* clusmap,
441 // TrkVec* trkvec)
442 // {
443 
444 // //-- Make a vector to hold all track candidates for a given first cluster
445 // TrkVec trkcand;
446 
447 // //-- Loop over all clusters in the first layer of the cluster map
448 // unsigned int lyr = (clusmap->begin())->first;
449 // for (LyrClusMap::iterator itr = clusmap->lower_bound(lyr);
450 // itr != clusmap->upper_bound(lyr);
451 // ++itr)
452 // {
453 // trkcand.clear();
454 
455 // //-- Find all track candidates for this cluster
456 // ClusVec trk;
457 // trk.insert(itr->second);
458 // LinkClusters(&trk, lyr+1, clusmap, &trkcand);
459 
460 // //-- Choose the best track candidate using this cluster
461 // ClusVec best_trk = ChooseBestTrk(&trkcand);
462 // if ( best_trk.size() > 2 )
463 // trkvec->insert(best_trk);
464 // }
465 
466 // return trkvec->size();
467 // }
468 
469 // /*
470 // *
471 // */
472 // void AnaMvtxTelescopeHits::LinkClusters(ClusVec* trk,
473 // unsigned int next_lyr,
474 // LyrClusMap* clusmap,
475 // TrkVec* trkvec)
476 // {
477 
478 // //-- Check layer
479 // if ( next_lyr > 3 )
480 // {
481 // // add this trk to the trk vector, we're done with this trk
482 // // we only consider tracks with at least 3 clusters
483 // if ( trk->size() > 2 )
484 // trkvec->push_back(*trk);
485 // return;
486 // }
487 
488 // bool found = false;
489 
490 // //-- Search for clusters in the next layer
491 // for ( LyrClusMap::iterator itr = LyrClusMap.lower_bound(next_lyr);
492 // itr != LyrClusMap.upper_bound(next_lyr);
493 // ++itr)
494 // {
495 // // make a copy of this trk candidate and run the recursion
496 // ClusVec new_trk = *trk;
497 // new_trk.insert(itr->second);
498 // LinkClusters(&new_trk, next_lyr + 1, clusmap, trkvec);
499 // found = true;
500 // }
501 
502 // //-- if no cluster found in next layer, try next layer
503 // if ( !found )
504 // {
505 // LinkClusters(*trk, next_lyr + 2, clusmap, trkvec);
506 // }
507 
508 // }
509 
510 // /*
511 // * Choose the best track candidates which all start with the same cluster
512 // */
513 // ClusVec AnaMvtxTelescopeHits::ChooseBestTrk(TrkVec* trkcnd)
514 // {
515 
516 // int best_idx = -1;
517 // double best_chi2 = 100;
518 // for (TrkVec::iterator itr = trkcnd->begin();
519 // itr != trkcnd->end();
520 // ++itr)
521 // {
522 // double chi2 = FitTrk(&itr->second);
523 // if ( chi2 < best_chi2 )
524 // {
525 // best_chi2 = chi2;
526 // best_idx = itr->first;
527 // }
528 // }
529 
530 // if ( best_idx >= 0 && best_chi2 < 5 )
531 // {
532 // return trkcnd->at(best_idx);
533 // }
534 // else
535 // {
536 // return ClusVec;
537 // }
538 // }
539 
540 // /*
541 // *
542 // */
543 // double AnaMvtxTelescopeHits::FitTrk(ClusVec* trk)
544 // {
545 // return 0;
546 // }
547 
548 /*
549  * Simple generalized linear least-squares fitter.
550  * Solve y(X) = beta'*X + error(L) for beta (vector of regression coefs.)
551  * L is the inverse covariance matrix for y.
552  * Least-squares solution involves solving X' L X beta = X' L y
553  *
554  */
555 /*
556 TVectorD
557 AnaMvtxTelescopeHits::SolveGLS(TMatrixD &X, TVectorD &y, TMatrixD &L, TMatrixD &cov)
558 {
559 
560  double tol = 1.e-15;
561  TMatrixD XT(X); XT.T();
562  TMatrixD A = XT * L * X;
563  TVectorD b = XT * L * y;
564 
565  // Now solve A*beta = b using SVD. Decompose A = U S V'.
566  TDecompSVD svd(A);
567  TMatrixD UT = svd.GetU(); UT.T();
568  TMatrixD V = svd.GetV();
569 
570  // Construct Moore-Penrose pseudoinverse of S
571  TVectorD s = svd.GetSig();
572  int n = s.GetNrows();
573  TMatrixD Sd(n, n);
574  for (int i = 0; i < n; i++)
575  Sd(i, i) = s(i) > tol ? 1. / s(i) : 0.;
576 
577  // Result
578  TVectorD beta = V * Sd * UT * b;
579 
580  // and covariance matrix
581  V *= Sd;
582  TMatrixD C(V, TMatrixD::kMultTranspose, V);
583  cov.ResizeTo(C);
584  cov = C;
585 
586  return beta;
587 }
588 */
589