StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFwdTrackMaker.cxx
1 #include "StFwdTrackMaker/StFwdTrackMaker.h"
2 #include "StFwdTrackMaker/include/Tracker/FwdHit.h"
3 #include "StFwdTrackMaker/include/Tracker/FwdTracker.h"
4 #include "StFwdTrackMaker/include/Tracker/TrackFitter.h"
5 #include "StFwdTrackMaker/include/Tracker/FwdGeomUtils.h"
6 #include "StFwdTrackMaker/include/Tracker/ObjExporter.h"
7 
8 #include "KiTrack/IHit.h"
9 #include "GenFit/Track.h"
10 #include "GenFit/GFRaveVertexFactory.h"
11 
12 #include "TMath.h"
13 
14 #include <limits>
15 #include <map>
16 #include <string>
17 #include <string>
18 #include <vector>
19 
20 #include "StBFChain/StBFChain.h"
21 
22 #include "StEvent/StEvent.h"
23 #include "StEvent/StGlobalTrack.h"
24 #include "StEvent/StHelixModel.h"
25 #include "StEvent/StPrimaryTrack.h"
26 #include "StEvent/StRnDHit.h"
27 #include "StEvent/StRnDHitCollection.h"
28 #include "StEvent/StTrack.h"
29 #include "StEvent/StTrackGeometry.h"
30 #include "StEvent/StTrackNode.h"
31 #include "StEvent/StPrimaryVertex.h"
32 #include "StEvent/StEnumerations.h"
33 #include "StEvent/StTrackDetectorInfo.h"
34 #include "StEvent/StFttPoint.h"
35 #include "StEvent/StFcsHit.h"
36 #include "StEvent/StFcsCluster.h"
37 #include "StEvent/StFttCollection.h"
38 #include "StEvent/StFcsCollection.h"
39 #include "StEvent/StTriggerData.h"
40 #include "StEvent/StFstHitCollection.h"
41 #include "StEvent/StFstHit.h"
42 #include "StEvent/StFwdTrackCollection.h"
43 #include "StChain/StChainOpt.h"
44 
45 #include "StEventUtilities/StEventHelper.h"
46 
47 #include "tables/St_g2t_fts_hit_Table.h"
48 #include "tables/St_g2t_track_Table.h"
49 #include "tables/St_g2t_vertex_Table.h"
50 #include "tables/St_g2t_event_Table.h"
51 
52 #include "StarMagField/StarMagField.h"
53 
54 #include "St_base/StMessMgr.h"
55 #include "StarClassLibrary/StPhysicalHelix.hh"
56 #include "StarClassLibrary/SystemOfUnits.h"
57 
58 #include <SystemOfUnits.h>
59 #include <exception>
60 
61 #include "TROOT.h"
62 #include "TLorentzVector.h"
63 
64 #include "StRoot/StEpdUtil/StEpdGeom.h"
65 #include "StFcsDbMaker/StFcsDb.h"
66 #include "StFstUtil/StFstCollection.h"
67 
68 #include "StEvent/StFwdTrack.h"
69 #include "GenFit/AbsMeasurement.h"
70 
71 
72 FwdSystem* FwdSystem::sInstance = nullptr;
73 TMVA::Reader * BDTCrit2::reader = nullptr;
74 float BDTCrit2::Crit2_RZRatio = -999;
75 float BDTCrit2::Crit2_DeltaRho = -999;
76 float BDTCrit2::Crit2_DeltaPhi = -999;
77 float BDTCrit2::Crit2_StraightTrackRatio = -999;
78 
79 
80 
81 //_______________________________________________________________________________________
82 class GenfitUtils{
83  public:
84 
85  // For now, accept anything we are passed, no matter what it is or how bad it is
86  template<typename T> static bool accept( T ) { return true; }
87 
88 
89 }; // GenfitUtils
90 
91 // Basic sanity cuts on genfit tracks
92 template<> bool GenfitUtils::accept( genfit::Track *track )
93 {
94  // This also gets rid of failed fits (but may need to explicitly
95  // for fit failure...)
96  if (track->getNumPoints() <= 0 ) return false; // fit may have failed
97 
98  auto cardinal = track->getCardinalRep();
99 
100  // Check that the track fit converged
101  auto status = track->getFitStatus( cardinal );
102  if ( !status->isFitConverged() ) {
103  return false;
104  }
105 
106 
107  // Next, check that all points on the track have fitter info
108  // (may be another indication of a failed fit?)
109  for ( auto point : track->getPoints() ) {
110  if ( !point->hasFitterInfo(cardinal) ) {
111  return false;
112  }
113  }
114 
115  // Following line fails with an exception, because some tracks lack
116  // forward update, or prediction in fitter info at the first point
117  //
118  // genfit::KalmanFitterInfo::getFittedState(bool) const of
119  // GenFit/fitters/src/KalmanFitterInfo.cc:250
120 
121  // Fitted state at the first point
122  // const auto &atFirstPoint = track->getFittedState();
123 
124  // Getting the fitted state from a track occasionally fails, because
125  // the first point on the fit doesn't have forward/backward fit
126  // information. So we want the first point with fit info...
127 
128  genfit::TrackPoint* first = nullptr;
129  unsigned int ipoint = 0;
130  for ( ipoint = 0; ipoint < track->getNumPoints(); ipoint++ ) {
131  first = track->getPointWithFitterInfo( ipoint );
132  if ( first ) break;
133  }
134 
135  // No points on the track have fit information
136  if ( !first ) {
137  LOG_WARN << "No fit information on fwd genfit track" << endm;
138  return false;
139  }
140 
141  auto& fittedState= track->getFittedState(ipoint);
142 
143  TVector3 momentum = fittedState.getMom();
144  double pt = momentum.Perp();
145 
146  if (pt < 0.10 ) return false; // below this
147 
148  return true;
149 
150 };
151 
152 //______________________________________________________________________________________
153 class SiRasterizer {
154  public:
155  SiRasterizer() {}
156  SiRasterizer(FwdTrackerConfig &_cfg) { setup(_cfg); }
157  ~SiRasterizer() {}
158  void setup(FwdTrackerConfig &_cfg) {
159  cfg = _cfg;
160  mRasterR = cfg.get<double>("SiRasterizer:r", 3.0);
161  mRasterPhi = cfg.get<double>("SiRasterizer:phi", 0.1);
162  }
163 
164  bool active() {
165  return cfg.get<bool>("SiRasterizer:active", false);
166  }
167 
168  TVector3 raster(TVector3 p0) {
169  TVector3 p = p0;
170  double r = p.Perp();
171  double phi = p.Phi();
172  const double minR = 5.0;
173  // 5.0 is the r minimum of the Si
174  p.SetPerp(minR + (std::floor((r - minR) / mRasterR) * mRasterR + mRasterR / 2.0));
175  p.SetPhi(-TMath::Pi() + (std::floor((phi + TMath::Pi()) / mRasterPhi) * mRasterPhi + mRasterPhi / 2.0));
176  return p;
177  }
178 
179  FwdTrackerConfig cfg;
180  double mRasterR, mRasterPhi;
181 };
182 
183 // Wrapper class around the forward tracker
184 class ForwardTracker : public ForwardTrackMaker {
185  public:
186  // Replaces original initialization. Config file and hitloader
187  // will be provided by the maker.
188  void initialize( TString geoCache, bool genHistograms ) {
189  nEvents = 1; // only process single event
190 
191  // Create the forward system...
192  FwdSystem::sInstance = new FwdSystem();
193 
194  // make our quality plotter
195  mQualityPlotter = new QualityPlotter(mConfig);
196  mQualityPlotter->makeHistograms(mConfig.get<size_t>("TrackFinder:nIterations", 1));
197 
198  // initialize the track fitter
199  mTrackFitter = new TrackFitter(mConfig, geoCache);
200  mTrackFitter->setGenerateHistograms(genHistograms);
201  mTrackFitter->setup();
202 
203  ForwardTrackMaker::initialize( geoCache, genHistograms );
204  }
205 
206  void finish() {
207 
208  if ( mGenHistograms ){
209  mQualityPlotter->finish();
210  writeEventHistograms();
211  }
212 
213  if (FwdSystem::sInstance){
214  delete FwdSystem::sInstance;
215  FwdSystem::sInstance = 0;
216  }
217  if (mQualityPlotter){
218  delete mQualityPlotter;
219  mQualityPlotter = 0;
220  }
221  if (mTrackFitter){
222  delete mTrackFitter;
223  mTrackFitter= 0;
224  }
225  }
226 };
227 
228 //________________________________________________________________________
229 StFwdTrackMaker::StFwdTrackMaker() : StMaker("fwdTrack"), mGenHistograms(false), mGenTree(false), mForwardTracker(nullptr), mForwardData(nullptr){
230  SetAttr("useFtt",1); // Default Ftt on
231  SetAttr("useFst",1); // Default Fst on
232  SetAttr("useFcs",1); // Default Fcs on
233  SetAttr("config", "config.xml"); // Default configuration file (user may override before Init())
234  SetAttr("fillEvent",1); // fill StEvent
235 };
236 
237 int StFwdTrackMaker::Finish() {
238 
239  auto prevDir = gDirectory;
240  if ( mGenHistograms ) {
241 
242  // output file name
243  string name = mFwdConfig.get<string>("Output:url", "fwdTrackerOutput.root");
244  LOG_INFO << "Saving StFwdTrackMaker Histograms to ROOT file: " << name << endm;
245  TFile *fOutput = new TFile(name.c_str(), "RECREATE");
246  fOutput->cd();
247 
248  fOutput->mkdir("StFwdTrackMaker");
249  fOutput->cd("StFwdTrackMaker");
250  for (auto nh : mHistograms) {
251  nh.second->SetDirectory(gDirectory);
252  nh.second->Write();
253  }
254  fOutput->cd("");
255  }
256 
257  mForwardTracker->finish();
258 
259  prevDir->cd();
260 
261  if (mGenTree) {
262  mTreeFile->cd();
263  mTree->Write();
264  mTreeFile->Write();
265  }
266  return kStOk;
267 }
268 
269 void StFwdTrackMaker::LoadConfiguration() {
270  if (mConfigFile.length() < 5){
271  LOG_INFO << "Forward Tracker is using default config for ";
272  if ( defaultConfig == defaultConfigData ){
273  LOG_INFO << " DATA" << endm;
274  } else {
275  LOG_INFO << " Simulation" << endm;
276  }
277  mFwdConfig.load( defaultConfig, true );
278  } else {
279  LOG_INFO << "Forward Tracker is using config from file : " << mConfigFile << endm;
280  mFwdConfig.load( mConfigFile );
281  }
282  configLoaded = true;
283 }
284 
285 //________________________________________________________________________
286 int StFwdTrackMaker::Init() {
287 
288  if ( !configLoaded ){
289  LoadConfiguration();
290  }
291 
292  if (mGenTree) {
293  mTreeFile = new TFile("fwdtree.root", "RECREATE");
294  mTree = new TTree("fwd", "fwd tracking tree");
295  mTree->Branch("fttN", &mTreeData. fttN, "fttN/I");
296  mTree->Branch("fttX", &mTreeData. fttX );
297  mTree->Branch("fttY", &mTreeData. fttY );
298  mTree->Branch("fttZ", &mTreeData. fttZ );
299 
300  mTree->Branch("fttTrackId", &mTreeData. fttTrackId );
301  mTree->Branch("fttVolumeId", &mTreeData. fttVolumeId );
302  mTree->Branch("fttPt", &mTreeData. fttPt );
303  mTree->Branch("fttVertexId", &mTreeData. fttVertexId );
304 
305  mTree->Branch("fstN", &mTreeData. fstN, "fstN/I");
306  mTree->Branch("fstX", &mTreeData. fstX );
307  mTree->Branch("fstY", &mTreeData. fstY );
308  mTree->Branch("fstZ", &mTreeData. fstZ );
309  mTree->Branch("fstTrackId", &mTreeData. fstTrackId );
310 
311  mTree->Branch("fcsN", &mTreeData. fcsN, "fcsN/I");
312  mTree->Branch("fcsX", &mTreeData. fcsX );
313  mTree->Branch("fcsY", &mTreeData. fcsY );
314  mTree->Branch("fcsZ", &mTreeData. fcsZ );
315  mTree->Branch("fcsDet", &mTreeData. fcsDet );
316 
317  // mc tracks
318  mTree->Branch("mcN", &mTreeData. mcN, "mcN/I");
319  mTree->Branch("mcPt", &mTreeData. mcPt );
320  mTree->Branch("mcEta", &mTreeData. mcEta );
321  mTree->Branch("mcPhi", &mTreeData. mcPhi );
322  mTree->Branch("mcCharge", &mTreeData. mcCharge );
323  mTree->Branch("mcVertexId", &mTreeData. mcVertexId );
324 
325  // mcverts
326  mTree->Branch("vmcN", &mTreeData. vmcN, "vmcN/I");
327  mTree->Branch("vmcX", &mTreeData. vmcX );
328  mTree->Branch("vmcY", &mTreeData. vmcY );
329  mTree->Branch("vmcZ", &mTreeData. vmcZ );
330 
331  // rcverts
332  mTree->Branch("vrcN", &mTreeData. vrcN, "vrcN/I");
333  mTree->Branch("vrcX", &mTreeData. vrcX );
334  mTree->Branch("vrcY", &mTreeData. vrcY );
335  mTree->Branch("vrcZ", &mTreeData. vrcZ );
336 
337  // rc tracks
338  mTree->Branch("rcN", &mTreeData. rcN, "rcN/I");
339  mTree->Branch("rcPt", &mTreeData. rcPt );
340  mTree->Branch("rcEta", &mTreeData. rcEta );
341  mTree->Branch("rcPhi", &mTreeData. rcPhi );
342  mTree->Branch("rcCharge", &mTreeData. rcCharge );
343  mTree->Branch("rcTrackId", &mTreeData. rcTrackId );
344  mTree->Branch("rcNumFST", &mTreeData. rcNumFST );
345  mTree->Branch("rcNumFTT", &mTreeData. rcNumFTT );
346  mTree->Branch("rcNumPV", &mTreeData. rcNumPV );
347  mTree->Branch("rcQuality", &mTreeData. rcQuality );
348 
349  mTree->Branch("thdN", &mTreeData. thdN, "thdN/I");
350  mTree->Branch("thdX", &mTreeData. thdX );
351  mTree->Branch("thdY", &mTreeData. thdY );
352  mTree->Branch("thaX", &mTreeData. thaX );
353  mTree->Branch("thaY", &mTreeData. thaY );
354  mTree->Branch("thaZ", &mTreeData. thaZ );
355 
356  // track projections
357  mTree->Branch("tprojN", &mTreeData. tprojN, "tprojN/I");
358  mTree->Branch("tprojIdD", &mTreeData. tprojIdD);
359  mTree->Branch("tprojIdT", &mTreeData. tprojIdT);
360  mTree->Branch("tprojX", &mTreeData. tprojX);
361  mTree->Branch("tprojY", &mTreeData. tprojY);
362  mTree->Branch("tprojZ", &mTreeData. tprojZ);
363  mTree->Branch("tprojPx", &mTreeData. tprojPx);
364  mTree->Branch("tprojPy", &mTreeData. tprojPy);
365  mTree->Branch("tprojPz", &mTreeData. tprojPz);
366 
367  std::string path = "TrackFinder.Iteration[0].SegmentBuilder";
368  std::vector<string> paths = mFwdConfig.childrenOf(path);
369 
370  if (mTreeData.saveCrit){
371  for (string p : paths) {
372  string name = mFwdConfig.get<string>(p + ":name", "");
373  mTreeData.Crits[name]; // create the entry
374  mTree->Branch(name.c_str(), &mTreeData.Crits[name]);
375  mTree->Branch((name + "_trackIds").c_str(), &mTreeData.CritTrackIds[name]);
376 
377  if ( name == "Crit2_RZRatio" ){
378  string n = name + "_x1";
379  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
380 
381  n = name + "_y1";
382  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
383 
384  n = name + "_z1";
385  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
386 
387  n = name + "_x2";
388  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
389 
390  n = name + "_y2";
391  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
392 
393  n = name + "_z2";
394  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
395 
396  n = name + "_h1";
397  mTreeData.CritTrackIds[(n)]; mTree->Branch(n.c_str(), &mTreeData.CritTrackIds[n]);
398  n = name + "_h2";
399  mTreeData.CritTrackIds[(n)]; mTree->Branch(n.c_str(), &mTreeData.CritTrackIds[n]);
400  n = name + "_h3";
401  mTreeData.CritTrackIds[(n)]; mTree->Branch(n.c_str(), &mTreeData.CritTrackIds[n]);
402  }
403 
404  if ( name == "Crit2_BDT" ){
405  string n = name + "_DeltaPhi";
406  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
407  n = name + "_DeltaRho";
408  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
409  n = name + "_RZRatio";
410  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
411  n = name + "_StraightTrackRatio";
412  mTreeData.Crits[(n)]; mTree->Branch(n.c_str(), &mTreeData.Crits[n]);
413  }
414  }
415 
416  // Three hit criteria
417  path = "TrackFinder.Iteration[0].ThreeHitSegments";
418  paths = mFwdConfig.childrenOf(path);
419 
420  for (string p : paths) {
421  string name = mFwdConfig.get<string>(p + ":name", "");
422  mTreeData.Crits[name]; // create the entry
423  mTree->Branch(name.c_str(), &mTreeData.Crits[name]);
424  mTree->Branch((name + "_trackIds").c_str(), &mTreeData.CritTrackIds[name]);
425  }
426  }
427 
428  mTree->SetAutoFlush(0);
429  } // gen tree
430 
431 
433  GetDataBase("VmcGeometry");
434 
435 
436 
437  TString geoCache = GetChainOpt()->GetFileOut();
438  if ( geoCache=="" )
439  geoCache = GetChainOpt()->GetFileIn();
440 
441  // Strip out @ symbol
442  geoCache = geoCache.ReplaceAll("@","");
443  // Strip off the last extention in the geoCache
444  geoCache = geoCache( 0, geoCache.Last('.') );
445  // Append geom.root to the extentionless geoCache
446  geoCache+=".geom.root";
447 
448  // create an SiRasterizer in case we need it
449  mSiRasterizer = std::shared_ptr<SiRasterizer>( new SiRasterizer(mFwdConfig));
450  mForwardTracker = std::shared_ptr<ForwardTracker>(new ForwardTracker( ));
451  mForwardTracker->setConfig(mFwdConfig);
452 
453  // only save criteria values if we are generating a tree.
454  mForwardTracker->setSaveCriteriaValues(mGenTree);
455 
456  mForwardData = std::shared_ptr<FwdDataSource>(new FwdDataSource());
457  mForwardTracker->setData(mForwardData);
458  mForwardTracker->initialize( geoCache, mGenHistograms );
459 
460  if ( mGenHistograms ){
461  mHistograms["fwdVertexZ"] = new TH1D("fwdVertexZ", "FWD Vertex (RAVE);z", 1000, -50, 50);
462  mHistograms["fwdVertexXY"] = new TH2D("fwdVertexXY", "FWD Vertex (RAVE);x;y", 100, -1, 1, 100, -1, 1);
463  mHistograms["fwdVertexDeltaZ"] = new TH2D("fwdVertexDeltaZ", "FWD Vertex - MC Vertex;#Delta z", 100, -1, 1, 100, -1, 1);
464 
465  mHistograms["McEventEta"] = new TH1D("McEventEta", ";MC Track Eta", 1000, -5, 5);
466  mHistograms["McEventPt"] = new TH1D("McEventPt", ";MC Track Pt (GeV/c)", 1000, 0, 10);
467  mHistograms["McEventPhi"] = new TH1D("McEventPhi", ";MC Track Phi", 1000, 0, 6.2831852);
468 
469  // these are tracks within 2.5 < eta < 4.0
470  mHistograms["McEventFwdEta"] = new TH1D("McEventFwdEta", ";MC Track Eta", 1000, -5, 5);
471  mHistograms["McEventFwdPt"] = new TH1D("McEventFwdPt", ";MC Track Pt (GeV/c)", 1000, 0, 10);
472  mHistograms["McEventFwdPhi"] = new TH1D("McEventFwdPhi", ";MC Track Phi", 1000, 0, 6.2831852);
473 
474  // create mHistograms
475  mHistograms["nMcTracks"] = new TH1I("nMcTracks", ";# MC Tracks/Event", 1000, 0, 1000);
476  mHistograms["nMcTracksFwd"] = new TH1I("nMcTracksFwd", ";# MC Tracks/Event", 1000, 0, 1000);
477  mHistograms["nMcTracksFwdNoThreshold"] = new TH1I("nMcTracksFwdNoThreshold", ";# MC Tracks/Event", 1000, 0, 1000);
478 
479  mHistograms["nHitsSTGC"] = new TH1I("nHitsSTGC", ";# STGC Hits/Event", 1000, 0, 1000);
480  mHistograms["nHitsFSI"] = new TH1I("nHitsFSI", ";# FSIT Hits/Event", 1000, 0, 1000);
481 
482  mHistograms["stgc_volume_id"] = new TH1I("stgc_volume_id", ";stgc_volume_id", 50, 0, 50);
483  mHistograms["fsi_volume_id"] = new TH1I("fsi_volume_id", ";fsi_volume_id", 50, 0, 50);
484 
485  mHistograms["fsiHitDeltaR"] = new TH1F("fsiHitDeltaR", "FSI; delta r (cm); ", 500, -5, 5);
486  mHistograms["fsiHitDeltaPhi"] = new TH1F("fsiHitDeltaPhi", "FSI; delta phi; ", 500, -5, 5);
487 
488  // there are 4 stgc stations
489  for (int i = 0; i < 4; i++) {
490  mHistograms[TString::Format("stgc%dHitMap", i).Data()] = new TH2F(TString::Format("stgc%dHitMap", i), TString::Format("STGC Layer %d; x (cm); y(cm)", i), 200, -100, 100, 200, -100, 100);
491 
492  mHistograms[TString::Format("stgc%dHitMapPrim", i).Data()] = new TH2F(TString::Format("stgc%dHitMapPrim", i), TString::Format("STGC Layer %d; x (cm); y(cm)", i), 200, -100, 100, 200, -100, 100);
493  mHistograms[TString::Format("stgc%dHitMapSec", i).Data()] = new TH2F(TString::Format("stgc%dHitMapSec", i), TString::Format("STGC Layer %d; x (cm); y(cm)", i), 200, -100, 100, 200, -100, 100);
494  }
495 
496  // There are 3 silicon stations
497  for (int i = 0; i < 3; i++) {
498  mHistograms[TString::Format("fsi%dHitMap", i).Data()] = new TH2F(TString::Format("fsi%dHitMap", i), TString::Format("FSI Layer %d; x (cm); y(cm)", i), 200, -100, 100, 200, -100, 100);
499  mHistograms[TString::Format("fsi%dHitMapZ", i).Data()] = new TH2F(TString::Format("fsi%dHitMapZ", i), TString::Format("FSI Layer %d; x (cm); y(cm)", i), 200, -100, 100, 200, -100, 100);
500 
501  mHistograms[TString::Format("fsi%dHitMapR", i).Data()] = new TH1F(TString::Format("fsi%dHitMapR", i), TString::Format("FSI Layer %d; r (cm); ", i), 500, 0, 50);
502  mHistograms[TString::Format("fsi%dHitMapPhi", i).Data()] = new TH1F(TString::Format("fsi%dHitMapPhi", i), TString::Format("FSI Layer %d; phi; ", i), 320, 0, TMath::Pi() * 2 + 0.1);
503  }
504 
505  } // mGenHistograms
506  LOG_DEBUG << "StFwdTrackMaker::Init" << endm;
507  return kStOK;
508 };
509 
510 TMatrixDSym makeSiCovMat(TVector3 hit, FwdTrackerConfig &xfg) {
511  // we can calculate the CovMat since we know the det info, but in future we should probably keep this info in the hit itself
512 
513  float rSize = xfg.get<float>("SiRasterizer:r", 3.0);
514  float phiSize = xfg.get<float>("SiRasterizer:phi", 0.004);
515 
516  // measurements on a plane only need 2x2
517  // for Si geom we need to convert from cylindrical to cartesian coords
518  TMatrixDSym cm(2);
519  TMatrixD T(2, 2);
520  TMatrixD J(2, 2);
521  const float x = hit.X();
522  const float y = hit.Y();
523  const float R = sqrt(x * x + y * y);
524  const float cosphi = x / R;
525  const float sinphi = y / R;
526  const float sqrt12 = sqrt(12.);
527 
528  const float dr = rSize / sqrt12;
529  const float dphi = (phiSize) / sqrt12;
530 
531  // Setup the Transposed and normal Jacobian transform matrix;
532  // note, the si fast sim did this wrong
533  // row col
534  T(0, 0) = cosphi;
535  T(0, 1) = -R * sinphi;
536  T(1, 0) = sinphi;
537  T(1, 1) = R * cosphi;
538 
539  J(0, 0) = cosphi;
540  J(0, 1) = sinphi;
541  J(1, 0) = -R * sinphi;
542  J(1, 1) = R * cosphi;
543 
544  TMatrixD cmcyl(2, 2);
545  cmcyl(0, 0) = dr * dr;
546  cmcyl(1, 1) = dphi * dphi;
547 
548  TMatrixD r = T * cmcyl * J;
549 
550  // note: float sigmaX = sqrt(r(0, 0));
551  // note: float sigmaY = sqrt(r(1, 1));
552 
553  cm(0, 0) = r(0, 0);
554  cm(1, 1) = r(1, 1);
555  cm(0, 1) = r(0, 1);
556  cm(1, 0) = r(1, 0);
557 
558  TMatrixDSym tamvoc(3);
559  tamvoc( 0, 0 ) = cm(0, 0); tamvoc( 0, 1 ) = cm(0, 1); tamvoc( 0, 2 ) = 0.0;
560  tamvoc( 1, 0 ) = cm(1, 0); tamvoc( 1, 1 ) = cm(1, 1); tamvoc( 1, 2 ) = 0.0;
561  tamvoc( 2, 0 ) = 0.0; tamvoc( 2, 1 ) = 0.0; tamvoc( 2, 2 ) = 0.01*0.01;
562 
563 
564  return tamvoc;
565 }
566 
567 void StFwdTrackMaker::loadFttHits( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
568  LOG_DEBUG << "Loading FTT Hits" << endm;
569  // Get the StEvent handle to see if the rndCollection is available
570  StEvent *event = (StEvent *)GetDataSet("StEvent");
571  string fttFromSource = mFwdConfig.get<string>( "Source:ftt", "" );
572 
573  if (!event){
574  LOG_ERROR << "No StEvent, cannot load Ftt Data" << endm;
575  return;
576  }
577 
578  // Load GEANT hits directly if requested
579  if ( "GEANT" == fttFromSource ) {
580  LOG_DEBUG << "Loading sTGC hits directly from GEANT hits" << endm;
581  loadFttHitsFromGEANT( mcTrackMap, hitMap, count );
582  return;
583  }
584 
585  StFttCollection *col = event->fttCollection();
586  // From Data
587  if ( col || "DATA" == fttFromSource ) {
588  loadFttHitsFromStEvent( mcTrackMap, hitMap, count );
589  return;
590  }
591 } // loadFttHits
592 
593 void StFwdTrackMaker::loadFttHitsFromStEvent( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
594  LOG_DEBUG << "Loading FTT Hits from Data" << endm;
595  StEvent *event = (StEvent *)GetDataSet("StEvent");
596  StFttCollection *col = event->fttCollection();
597 
598  mTreeData.fttN = 0;
599 
600  if ( col && col->numberOfPoints() > 0 ){
601  LOG_DEBUG << "The Ftt Collection has " << col->numberOfPoints() << " points" << endm;
602  TMatrixDSym hitCov3(3);
603  const double sigXY = 0.2; //
604  hitCov3(0, 0) = sigXY * sigXY;
605  hitCov3(1, 1) = sigXY * sigXY;
606  hitCov3(2, 2) = 4; // unused since they are loaded as points on plane
607  for ( auto point : col->points() ){
608 
609  FwdHit *hit = new FwdHit(count++, point->xyz().x()/10.0, point->xyz().y()/10.0, point->xyz().z(), -point->plane(), 0, hitCov3, nullptr);
610  mFttHits.push_back( TVector3( point->xyz().x()/10.0, point->xyz().y()/10.0, point->xyz().z() ) );
611  if ( mGenHistograms ) {
612  mHistograms[TString::Format("stgc%dHitMapSec", point->plane()).Data()]->Fill(point->xyz().x()/10.0, point->xyz().y()/10.0);
613  }
614  // Add the hit to the hit map
615  hitMap[hit->getSector()].push_back(hit);
616 
617  if (mGenTree && (unsigned)mTreeData.fttN < MAX_TREE_ELEMENTS) {
618  LOG_DEBUG << "FttPoint( " << TString::Format( "[plane=%d, quad=%d, nClu=%d]", point->plane(), point->quadrant(), point->nClusters() ) << point->xyz().x()/10.0 << ", " << point->xyz().y()/10.0 << ", " << point->xyz().z() << " )" << endm;
619  mTreeData.fttX.push_back( point->xyz().x()/10.0 );
620  mTreeData.fttY.push_back( point->xyz().y()/10.0 );
621  mTreeData.fttZ.push_back( point->xyz().z() );
622  mTreeData.fttTrackId.push_back( 0 );
623  mTreeData.fttVolumeId.push_back( point->plane() );
624  mTreeData.fttPt.push_back( 0 );
625  mTreeData.fttVertexId.push_back( 0 );
626  mTreeData.fttN++;
627  }
628  }
629 
630  return;
631  } else {
632  LOG_DEBUG << "The Ftt Collection is EMPTY points" << endm;
633  }
634  LOG_DEBUG << "Number of FTT in TTree: " << mTreeData.fttN << endm;
635 }
636 
637 void StFwdTrackMaker::loadFttHitsFromGEANT( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
638  /************************************************************/
639  // STGC Hits
640  St_g2t_fts_hit *g2t_stg_hits = (St_g2t_fts_hit *)GetDataSet("geant/g2t_stg_hit");
641 
642 
643  mTreeData.fttN = 0;
644  if (!g2t_stg_hits){
645  LOG_WARN << "geant/g2t_stg_hit is empty" << endm;
646  return;
647  }
648 
649  // make the Covariance Matrix once and then reuse
650  TMatrixDSym hitCov3(3);
651  const double sigXY = 0.01;
652  hitCov3(0, 0) = sigXY * sigXY;
653  hitCov3(1, 1) = sigXY * sigXY;
654  hitCov3(2, 2) = 1.0; // unused since they are loaded as points on plane
655 
656  int nstg = g2t_stg_hits->GetNRows();
657 
658  LOG_DEBUG << "This event has " << nstg << " stg hits in geant/g2t_stg_hit " << endm;
659  if ( mGenHistograms ) {
660  mHistograms["nHitsSTGC"]->Fill(nstg);
661  }
662 
663 
664  bool filterGEANT = mFwdConfig.get<bool>( "Source:fttFilter", false );
665  for (int i = 0; i < nstg; i++) {
666 
667  g2t_fts_hit_st *git = (g2t_fts_hit_st *)g2t_stg_hits->At(i);
668  if (0 == git)
669  continue; // geant hit
670  int track_id = git->track_p;
671  int volume_id = git->volume_id;
672  int plane_id = (volume_id - 1) / 100; // from 1 - 16. four chambers per station
673 
674  // only use the hits on the front modules
675  if ( volume_id % 2 ==0 )
676  continue;
677 
678  float x = git->x[0] + gRandom->Gaus(0, sigXY); // 100 micron blur according to approx sTGC reso
679  float y = git->x[1] + gRandom->Gaus(0, sigXY); // 100 micron blur according to approx sTGC reso
680  float z = git->x[2];
681 
682 
683 
684  if (mGenTree && (unsigned)mTreeData.fttN < MAX_TREE_ELEMENTS) {
685  mTreeData.fttX.push_back( x );
686  mTreeData.fttY.push_back( y );
687  mTreeData.fttZ.push_back( z );
688  mTreeData.fttTrackId.push_back( track_id );
689  mTreeData.fttVolumeId.push_back( plane_id );
690  mTreeData.fttPt.push_back( mcTrackMap[track_id]->mPt );
691  mTreeData.fttVertexId.push_back( mcTrackMap[track_id]->mStartVertex );
692  mTreeData.fttN++;
693  } else if ( mGenTree ){
694  LOG_WARN << "Truncating hits in TTree output" << endm;
695  }
696 
697  if ( mGenHistograms ){
698  mHistograms["stgc_volume_id"]->Fill(volume_id);
699  }
700 
701  if (plane_id < 4 && plane_id >= 0) {
702  if ( mGenHistograms ){
703  mHistograms[TString::Format("stgc%dHitMap", plane_id).Data()]->Fill(x, y);
704  }
705  } else {
706  continue;
707  }
708 
709  // this rejects GEANT hits with eta -999 - do we understand this effect?
710  if ( filterGEANT ) {
711  if ( mcTrackMap[track_id] && fabs(mcTrackMap[track_id]->mEta) > 5.0 ){
712 
713  if ( mGenHistograms )
714  mHistograms[TString::Format("stgc%dHitMapSec", plane_id).Data()]->Fill(x, y);
715  continue;
716  } else if ( mcTrackMap[track_id] && fabs(mcTrackMap[track_id]->mEta) < 5.0 ){
717  if ( mGenHistograms ) mHistograms[TString::Format("stgc%dHitMapPrim", plane_id).Data()]->Fill(x, y);
718  }
719  }
720 
721  FwdHit *hit = new FwdHit(count++, x, y, z, -plane_id, track_id, hitCov3, mcTrackMap[track_id]);
722 
723  // Add the hit to the hit map
724  hitMap[hit->getSector()].push_back(hit);
725  mFttHits.push_back( TVector3( x, y, z ) );
726 
727  // Add hit pointer to the track
728  if (mcTrackMap[track_id]){
729  mcTrackMap[track_id]->addHit(hit);
730  } else {
731  LOG_ERROR << "Cannot find MC track for GEANT hit (FTT), track_id = " << track_id << endm;
732  }
733  } // loop on hits
734 
735  if (mGenTree){
736  LOG_INFO << "Saving " << mTreeData.fttN << " hits in Tree" << endm;
737  }
738 } // loadFttHits
739 
740 void StFwdTrackMaker::loadFstHits( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
741  StEvent *event = (StEvent *)GetDataSet("StEvent");
742  StFstHitCollection *fstHitCollection = event->fstHitCollection();
743 
744  if ( fstHitCollection ){
745  // reuse this to store cov mat
746  TMatrixDSym hitCov3(3);
747  LOG_DEBUG << "StFstHitCollection is NOT NULL, loading hits" << endm;
748  for ( unsigned int iw = 0; iw < kFstNumWedges; iw++ ){
749  StFstWedgeHitCollection * wc = fstHitCollection->wedge( iw );
750 
751  for ( unsigned int is = 0; is < kFstNumSensorsPerWedge; is++ ){
752  StFstSensorHitCollection * sc = wc->sensor( is );
753 
754  StSPtrVecFstHit fsthits = sc->hits();
755  mTreeData.fstN = 0;
756  for ( unsigned int ih = 0; ih < fsthits.size(); ih++ ){
757  float vR = fsthits[ih]->localPosition(0);
758  float vPhi = fsthits[ih]->localPosition(1);
759  float vZ = fsthits[ih]->localPosition(2);
760 
761  const float dz0 = fabs( vZ - 151.75 );
762  const float dz1 = fabs( vZ - 165.248 );
763  const float dz2 = fabs( vZ - 178.781 );
764 
765  int d = 0 * ( dz0 < 1.0 ) + 1 * ( dz1 < 1.0 ) + 2 * ( dz2 < 1.0 );
766 
767 
768 
769  float x0 = vR * cos( vPhi );
770  float y0 = vR * sin( vPhi );
771  hitCov3 = makeSiCovMat( TVector3( x0, y0, vZ ), mFwdConfig );
772 
773  LOG_DEBUG << "FST HIT: d = " << d << ", x=" << x0 << ", y=" << y0 << ", z=" << vZ << endm;
774  mFstHits.push_back( TVector3( x0, y0, vZ) );
775 
776  FwdHit *hit = new FwdHit(count++, x0, y0, vZ, d, 0, hitCov3, nullptr);
777  // Add the hit to the hit map
778  hitMap[hit->getSector()].push_back(hit);
779 
780  mTreeData.fstX.push_back( x0 );
781  mTreeData.fstY.push_back( y0 );
782  mTreeData.fstZ.push_back( vZ );
783 
784  mTreeData.fstN++;
785  }
786  } // loop is
787  } // loop iw
788  LOG_DEBUG << " FOUND " << mFstHits.size() << " FST HITS" << endm;
789  return;
790  } // fstHitCollection
791 
792  StRnDHitCollection *rndCollection = nullptr;
793  if (nullptr != event) {
794  rndCollection = event->rndHitCollection();
795  }
796  bool siRasterizer = mFwdConfig.get<bool>( "SiRasterizer:active", false );
797  if ( siRasterizer || rndCollection == nullptr ){
798  LOG_DEBUG << "Loading Fst hits from GEANT with SiRasterizer" << endm;
799  loadFstHitsFromGEANT( mcTrackMap, hitMap, count );
800  } else {
801  LOG_DEBUG << "Loading Fst hits from StEvent" << endm;
802  loadFstHitsFromStEvent( mcTrackMap, hitMap, count );
803  }
804 } // loadFstHits
805 
806 void StFwdTrackMaker::loadFstHitsFromStEvent( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
807 
808  // Get the StEvent handle
809  StEvent *event = (StEvent *)GetDataSet("StEvent");
810  if (!event)
811  return;
812 
813  StRnDHitCollection *rndCollection = event->rndHitCollection();
814 
815  const StSPtrVecRnDHit &hits = rndCollection->hits();
816 
817  // we will reuse this to hold the cov mat
818  TMatrixDSym hitCov3(3);
819 
820  mTreeData.fstN = 0;
821  for (unsigned int fsthit_index = 0; fsthit_index < hits.size(); fsthit_index++) {
822  StRnDHit *hit = hits[fsthit_index];
823 
824  if ( hit->layer() > 6 ){
825  // skip sTGC hits here
826  continue;
827  }
828 
829  const StThreeVectorF pos = hit->position();
830 
831  StMatrixF covmat = hit->covariantMatrix();
832 
833  // copy covariance matrix element by element from StMatrixF
834  hitCov3(0,0) = covmat[0][0]; hitCov3(0,1) = covmat[0][1]; hitCov3(0,2) = covmat[0][2];
835  hitCov3(1,0) = covmat[1][0]; hitCov3(1,1) = covmat[1][1]; hitCov3(1,2) = covmat[1][2];
836  hitCov3(2,0) = covmat[2][0]; hitCov3(2,1) = covmat[2][1]; hitCov3(2,2) = covmat[2][2];
837 
838  FwdHit *fhit = new FwdHit(count++, hit->position().x(), hit->position().y(), hit->position().z(), hit->layer(), hit->idTruth(), hitCov3, mcTrackMap[hit->idTruth()]);
839  size_t index = hit->layer()-4;
840  if (mGenHistograms && index < 3 ){
841  ((TH2*)mHistograms[TString::Format("fsi%luHitMapZ", index).Data()]) -> Fill( hit->position().x(), hit->position().y(), hit->position().z() );
842  }
843 
844  // Add the hit to the hit map
845  hitMap[fhit->getSector()].push_back(fhit);
846  mFstHits.push_back( TVector3( hit->position().x(), hit->position().y(), hit->position().z()) );
847 
848  mTreeData.fstX.push_back( hit->position().x() );
849  mTreeData.fstY.push_back( hit->position().y() );
850  mTreeData.fstZ.push_back( hit->position().z() );
851  mTreeData.fstTrackId.push_back( hit->idTruth() );
852 
853  mTreeData.fstN++;
854 
855  }
856 } //loadFstHitsFromStEvent
857 
858 void StFwdTrackMaker::loadFstHitsFromGEANT( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){
859  /************************************************************/
860  // Load FSI Hits from GEANT
861  St_g2t_fts_hit *g2t_fsi_hits = (St_g2t_fts_hit *)GetDataSet("geant/g2t_fsi_hit");
862 
863  if ( !g2t_fsi_hits )
864  return;
865 
866  int nfsi = g2t_fsi_hits->GetNRows();
867 
868 
869  // reuse this to store cov mat
870  TMatrixDSym hitCov3(3);
871 
872  if ( mGenHistograms ) mHistograms["nHitsFSI"]->Fill(nfsi);
873 
874  mTreeData.fstN = 0;
875  for (int i = 0; i < nfsi; i++) {
876 
877  g2t_fts_hit_st *git = (g2t_fts_hit_st *)g2t_fsi_hits->At(i);
878 
879  if (0 == git)
880  continue; // geant hit
881 
882  int track_id = git->track_p;
883  int volume_id = git->volume_id; // 4, 5, 6
884  int d = volume_id / 1000; // disk id
885 
886  int plane_id = d - 4;
887  float x = git->x[0];
888  float y = git->x[1];
889  float z = git->x[2];
890 
891  if (mSiRasterizer->active()) {
892  TVector3 rastered = mSiRasterizer->raster(TVector3(git->x[0], git->x[1], git->x[2]));
893 
894  if ( mGenHistograms ) {
895  mHistograms["fsiHitDeltaR"]->Fill(std::sqrt(x * x + y * y) - rastered.Perp());
896  mHistograms["fsiHitDeltaPhi"]->Fill(std::atan2(y, x) - rastered.Phi());
897  }
898  x = rastered.X();
899  y = rastered.Y();
900  }
901 
902 
903  if ( mGenHistograms ) mHistograms["fsi_volume_id"]->Fill(d);
904 
905  if (plane_id < 3 && plane_id >= 0) {
906 
907  if ( mGenHistograms ) {
908  mHistograms[TString::Format("fsi%dHitMap", plane_id).Data()]->Fill(x, y);
909  mHistograms[TString::Format("fsi%dHitMapR", plane_id).Data()]->Fill(std::sqrt(x * x + y * y));
910  mHistograms[TString::Format("fsi%dHitMapPhi", plane_id).Data()]->Fill(std::atan2(y, x) + TMath::Pi());
911  }
912  } else {
913  continue;
914  }
915 
916  hitCov3 = makeSiCovMat( TVector3( x, y, z ), mFwdConfig );
917  FwdHit *hit = new FwdHit(count++, x, y, z, d, track_id, hitCov3, mcTrackMap[track_id]);
918  mFstHits.push_back( TVector3( x, y, z ) );
919 
920  mTreeData.fstX.push_back( x );
921  mTreeData.fstY.push_back( y );
922  mTreeData.fstZ.push_back( z );
923  mTreeData.fstTrackId.push_back( track_id );
924 
925  mTreeData.fstN++;
926 
927  // Add the hit to the hit map
928  hitMap[hit->getSector()].push_back(hit);
929  }
930 } // loadFstHitsFromGEANT
931 
932 size_t StFwdTrackMaker::loadMcTracks( FwdDataSource::McTrackMap_t &mcTrackMap ){
933 
934 
935  LOG_DEBUG << "Looking for GEANT sim vertex info" << endm;
936  St_g2t_vertex *g2t_vertex = (St_g2t_vertex *)GetDataSet("geant/g2t_vertex");
937 
938  if ( g2t_vertex != nullptr ) {
939  // Set the MC Vertex for track fitting
940  g2t_vertex_st *vert = (g2t_vertex_st*)g2t_vertex->At(0);
941  mForwardTracker->setEventVertex( TVector3( vert->ge_x[0], vert->ge_x[1], vert->ge_x[2] ) );
942  }
943 
944  // Get geant tracks
945  St_g2t_track *g2t_track = (St_g2t_track *)GetDataSet("geant/g2t_track");
946 
947  if (!g2t_track)
948  return 0;
949 
950  size_t nShowers = 0;
951 
952  mTreeData.mcN = 1;
953  LOG_DEBUG << g2t_track->GetNRows() << " mc tracks in geant/g2t_track " << endm;
954  if ( mGenHistograms ) mHistograms["nMcTracks"]->Fill(g2t_track->GetNRows());
955 
956  for (int irow = 0; irow < g2t_track->GetNRows(); irow++) {
957  g2t_track_st *track = (g2t_track_st *)g2t_track->At(irow);
958 
959  if (0 == track)
960  continue;
961 
962  int track_id = track->id;
963  float pt2 = track->p[0] * track->p[0] + track->p[1] * track->p[1];
964  float pt = std::sqrt(pt2);
965  float eta = track->eta;
966  float phi = std::atan2(track->p[1], track->p[0]); //track->phi;
967  int q = track->charge;
968 
969  if (!mcTrackMap[track_id] )
970  mcTrackMap[track_id] = shared_ptr<McTrack>(new McTrack(pt, eta, phi, q, track->start_vertex_p));
971 
972  if (mGenTree && (unsigned)mTreeData.mcN < MAX_TREE_ELEMENTS) {
973  mTreeData.mcPt.push_back( pt );
974  mTreeData.mcEta.push_back( eta );
975  mTreeData.mcPhi.push_back( phi );
976  mTreeData.mcCharge.push_back( q );
977  mTreeData.mcVertexId.push_back( track->start_vertex_p );
978 
979  if (track->is_shower)
980  nShowers++;
981 
982  mTreeData.mcN++;
983  } else if ( mGenTree ) {
984  LOG_WARN << "Truncating Mc tracks in TTree output" << endm;
985  }
986 
987  } // loop on track (irow)
988 
989 
990  // now check the Mc tracks against the McEvent filter
991  size_t nForwardTracks = 0;
992  size_t nForwardTracksNoThreshold = 0;
993  for (auto mctm : mcTrackMap ){
994  if ( mctm.second == nullptr ) continue;
995 
996  if ( mGenHistograms ){
997  mHistograms[ "McEventPt" ] ->Fill( mctm.second->mPt );
998  mHistograms[ "McEventEta" ] ->Fill( mctm.second->mEta );
999  mHistograms[ "McEventPhi" ] ->Fill( mctm.second->mPhi );
1000  }
1001 
1002  if ( mctm.second->mEta > 2.5 && mctm.second->mEta < 4.0 ){
1003 
1004  if ( mGenHistograms ){
1005  mHistograms[ "McEventFwdPt" ] ->Fill( mctm.second->mPt );
1006  mHistograms[ "McEventFwdEta" ] ->Fill( mctm.second->mEta );
1007  mHistograms[ "McEventFwdPhi" ] ->Fill( mctm.second->mPhi );
1008  }
1009 
1010  nForwardTracksNoThreshold++;
1011  if ( mctm.second->mPt > 0.05 )
1012  nForwardTracks++;
1013  }
1014  } // loop on mcTrackMap
1015 
1016  if ( mGenHistograms ) {
1017  mHistograms[ "nMcTracksFwd" ]->Fill( nForwardTracks );
1018  mHistograms[ "nMcTracksFwdNoThreshold" ]->Fill( nForwardTracksNoThreshold );
1019  }
1020 
1021 
1022  return nForwardTracks;
1023 } // loadMcTracks
1024 
1025 void StFwdTrackMaker::loadFcs( ) {
1026  StEvent *stEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
1027  StFcsDb* fcsDb=static_cast<StFcsDb*>(GetDataSet("fcsDb"));
1028  if ( !stEvent || !fcsDb ){
1029  return;
1030  }
1031  StFcsCollection* fcsCol = stEvent->fcsCollection();
1032  if ( !fcsCol ){
1033  return;
1034  }
1035 
1036  StEpdGeom epdgeo;
1037 
1038  mTreeData.fcsN = 0;
1039 
1040  // LOAD ECAL / HCAL CLUSTERS
1041  for ( int idet = 0; idet < 4; idet++ ){
1042  StSPtrVecFcsCluster& clusters = fcsCol->clusters(idet);
1043  int nc=fcsCol->numberOfClusters(idet);
1044  for ( int i = 0; i < nc; i++ ){
1045  StFcsCluster* clu = clusters[i];
1046  StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y());
1047  mFcsClusters.push_back( TVector3( xyz.x(), xyz.y(), xyz.z() - 2 ) );
1048  mFcsClusterEnergy.push_back( clu->energy() );
1049 
1050  mTreeData.fcsX.push_back( xyz.x() );
1051  mTreeData.fcsY.push_back( xyz.y() );
1052  mTreeData.fcsZ.push_back( xyz.z() - 2 );
1053  mTreeData.fcsDet.push_back( idet );
1054  }
1055  }
1056 
1057  // LOAD PRESHOWER HITS (EPD)
1058  for ( int det = 4; det < 6; det ++ ) {
1059 
1060  StSPtrVecFcsHit& hits = stEvent->fcsCollection()->hits(det);
1061  int nh=fcsCol->numberOfHits(det);
1062  for ( int i = 0; i < nh; i++ ){
1063  StFcsHit* hit=hits[i];
1064 
1065  if(det==kFcsPresNorthDetId || det==kFcsPresSouthDetId){ //EPD
1066  double zepd=375.0;
1067  int pp,tt,n;
1068  double x[5],y[5];
1069 
1070  if ( hit->energy() < 0.2 ) continue;
1071  fcsDb->getEPDfromId(det,hit->id(),pp,tt);
1072  epdgeo.GetCorners(100*pp+tt,&n,x,y);
1073  double x0 = (x[0] + x[1] + x[2] + x[3]) / 4.0;
1074  double y0 = (y[0] + y[1] + y[2] + y[3]) / 4.0;
1075  mFcsPreHits.push_back( TVector3( x0, y0, zepd ) );
1076 
1077  mTreeData.fcsX.push_back( x0 );
1078  mTreeData.fcsY.push_back( y0 );
1079  mTreeData.fcsZ.push_back( zepd );
1080  mTreeData.fcsDet.push_back( det );
1081 
1082  }
1083  }
1084  }
1085 } // loadFcs
1086 
1087 
1088 //________________________________________________________________________
1089 int StFwdTrackMaker::Make() {
1090  // START time for measuring tracking
1091  long long itStart = FwdTrackerUtils::nowNanoSecond();
1092 
1093  // Access forward Tracker maps
1094  FwdDataSource::McTrackMap_t &mcTrackMap = mForwardData->getMcTracks();
1095  FwdDataSource::HitMap_t &hitMap = mForwardData->getFttHits();
1096  FwdDataSource::HitMap_t &fsiHitMap = mForwardData->getFstHits();
1097 
1098  // clear vectors for visualization OBJ hits
1099  mFttHits.clear();
1100  mFstHits.clear();
1101  mFcsPreHits.clear();
1102  mFcsClusters.clear();
1103  mFwdTracks.clear();
1104 
1105  // default event vertex
1106  mForwardTracker->setEventVertex( TVector3( 0, 0, 0 ) );
1107 
1108  /**********************************************************************/
1109  // Load MC tracks
1110  size_t nForwardTracks = loadMcTracks( mcTrackMap );
1111  size_t maxForwardTracks = mFwdConfig.get<size_t>( "McEvent.Mult:max", 10000 );
1112  if ( nForwardTracks > maxForwardTracks ){
1113  LOG_WARN << "Skipping event with more than " << maxForwardTracks << " forward tracks" << endm;
1114  return kStOk;
1115  }
1116  LOG_DEBUG << "We have " << nForwardTracks << " forward MC tracks" << endm;
1117 
1118  /**********************************************************************/
1119  // Load sTGC
1120  LOG_DEBUG << ">>StFwdTrackMaker::loadFttHits" << endm;
1121  if ( IAttr("useFtt") ) {
1122  loadFttHits( mcTrackMap, hitMap );
1123  }
1124 
1125 
1126  /**********************************************************************/
1127  // Load FST
1128  LOG_DEBUG << ">>StFwdTrackMaker::loadFstHits" << endm;
1129  if ( IAttr("useFst") ) {
1130  loadFstHits( mcTrackMap, fsiHitMap );
1131  }
1132 
1133  /**********************************************************************/
1134  // Load FCS
1135  LOG_DEBUG << ">>StFwdTrackMaker::loadFcsHits" << endm;
1136  if ( IAttr("useFcs") ) {
1137  loadFcs();
1138  }
1139 
1140  /**********************************************************************/
1141  // Run Track finding + fitting
1142  LOG_DEBUG << ">>START Event Forward Tracking" << endm;
1143  mForwardTracker->doEvent();
1144  LOG_DEBUG << "<<FINISH Event Forward Tracking" << endm;
1145  LOG_DEBUG << "<<Made " << mForwardTracker -> getRecoTracks().size() << " GenFit Tracks" << endm;
1146 
1147 
1148  FitVertex();
1149 
1150  StEvent *stEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
1151 
1152  /**********************************************************************/
1153  // Run Track finding + fitting
1154 
1155  const auto &genfitTracks = mForwardTracker -> globalTracks();
1156  if ( mVisualize /* && genfitTracks.size() > 0 && genfitTracks.size() < 200*/ ) {
1157  const auto &seed_tracks = mForwardTracker -> getRecoTracks();
1158 
1159  ObjExporter woe;
1160  woe.output(
1161  TString::Format( "ev%lu", eventIndex ).Data(),
1162  stEvent,
1163  seed_tracks, genfitTracks, mRaveVertices,
1164  mFttHits, mFstHits, mFcsPreHits, mFcsClusters, mFcsClusterEnergy );
1165  eventIndex++;
1166  LOG_DEBUG << "Done Writing OBJ " << endm;
1167  } else if (mVisualize && genfitTracks.size() == 0) {
1168  LOG_DEBUG << "Skipping visualization, no FWD tracks" << endm;
1169  } else if (mVisualize && genfitTracks.size() >= 20) {
1170  LOG_DEBUG << "Skipping visualization, too many FWD tracks" << endm;
1171  }
1172 
1173  // Fill Track Deltas in ttree for helpful alignment info
1174  FillTrackDeltas();
1175 
1176  LOG_INFO << "Forward tracking on this event took " << (FwdTrackerUtils::nowNanoSecond() - itStart) * 1e-6 << " ms" << endm;
1177 
1178  if ( true && IAttr("fillEvent") ) {
1179 
1180  if (!stEvent) {
1181  LOG_WARN << "No StEvent found. Forward tracks will not be saved" << endm;
1182  return kStWarn;
1183  }
1184 
1185  FillEvent();
1186  } // IAttr FillEvent
1187 
1188  LOG_DEBUG << "Filling fwd Tree for event: " << GetIventNumber() << endm;
1189  FillTTree();
1190  return kStOK;
1191 } // Make
1192 
1193 
1194 StFwdTrack * StFwdTrackMaker::makeStFwdTrack( GenfitTrackResult &gtr, size_t indexTrack ){
1195  LOG_DEBUG << "StFwdTrackMaker::makeStFwdTrack()" << endm;
1196  StFwdTrack *fwdTrack = new StFwdTrack( );
1197 
1198  auto track = gtr.track;
1199  // if FST refit is available save that
1200  if ( gtr.nFST > 0 && gtr.fstTrack != nullptr){
1201  LOG_DEBUG << "\tSave FST refit track since we have FST points" << endm;
1202  track = gtr.fstTrack;
1203  } else if (gtr.nFST > 0 && gtr.fstTrack == nullptr) {
1204  LOG_WARN << "\tFST refit failed even though we have " << gtr.nFST << " FST points" << endm;
1205  }
1206 
1207  // Fit failed beyond use
1208  if ( track == nullptr ){
1209  LOG_DEBUG << "Track is nullptr, not saving StFwdTrack" << endm;
1210  return nullptr;
1211  }
1212 
1213  auto fitStatus = track->getFitStatus();
1214  if ( !fitStatus )
1215  return nullptr;
1216 
1217  // Fill charge and quality info
1218  fwdTrack->setDidFitConverge( fitStatus->isFitConverged() );
1219  fwdTrack->setDidFitConvergeFully( fitStatus->isFitConvergedFully() );
1220  fwdTrack->setNumberOfFailedPoints( fitStatus->getNFailedPoints() );
1221 
1222  fwdTrack->setNumberOfFitPoints( track->getNumPoints() );
1223  fwdTrack->setChi2( fitStatus->getChi2() );
1224  fwdTrack->setNDF( fitStatus->getNdf() );
1225  fwdTrack->setPval( fitStatus->getPVal() );
1226 
1227  auto cr = track->getCardinalRep();
1228  // charge at first point
1229  fwdTrack->setCharge( gtr.charge );
1230 
1231  TVector3 p = cr->getMom( track->getFittedState( 0, cr ));
1232  fwdTrack->setPrimaryMomentum( StThreeVectorD( gtr.momentum.X(), gtr.momentum.Y(), gtr.momentum.Z() ) );
1233  LOG_DEBUG << "Making StFwdTrack with " << TString::Format( "p=(%f, %f, %f)", fwdTrack->momentum().x(), fwdTrack->momentum().y(), fwdTrack->momentum().z() ) << endm;
1234 
1235  int nSeedPoints = 0;
1236  // store the seed points from FTT
1237  for ( auto s : gtr.trackSeed ){
1238  FwdHit * fh = static_cast<FwdHit*>( s );
1239  if (!fh) continue;
1240  float cov[9];
1241  cov[0] = fh->_covmat(0,0); cov[3] = fh->_covmat(1,0); cov[6] = fh->_covmat(2,0);
1242  cov[1] = fh->_covmat(0,1); cov[4] = fh->_covmat(1,1); cov[7] = fh->_covmat(2,1);
1243  cov[2] = fh->_covmat(0,2); cov[5] = fh->_covmat(1,2); cov[8] = fh->_covmat(2,2);
1244 
1245  StFwdTrackSeedPoint p( StThreeVectorD( fh->getX(), fh->getY(), fh->getZ() ), fh->getSector(), fh->getTrackId(), cov );
1246  fwdTrack->mFTTPoints.push_back( p );
1247  nSeedPoints++;
1248  }
1249 
1250  for ( auto s : gtr.fstSeed ){
1251  FwdHit * fh = static_cast<FwdHit*>( s );
1252  if (!fh) continue;
1253  float cov[9];
1254  cov[0] = fh->_covmat(0,0); cov[3] = fh->_covmat(1,0); cov[6] = fh->_covmat(2,0);
1255  cov[1] = fh->_covmat(0,1); cov[4] = fh->_covmat(1,1); cov[7] = fh->_covmat(2,1);
1256  cov[2] = fh->_covmat(0,2); cov[5] = fh->_covmat(1,2); cov[8] = fh->_covmat(2,2);
1257 
1258  StFwdTrackSeedPoint p( StThreeVectorD( fh->getX(), fh->getY(), fh->getZ() ), fh->getSector(), fh->getTrackId(), cov );
1259  fwdTrack->mFSTPoints.push_back( p );
1260  nSeedPoints++;
1261  }
1262 
1263  // set total number of seed points
1264  fwdTrack->setNumberOfSeedPoints( nSeedPoints );
1265 
1266  // compute projections to z-planes of various detectors
1267  vector<float> zPlanes = {
1268  0, // PV TODO, update with event vertex?
1269  151.750000, 165.248001, 178.781006, // FST
1270  280.904999, 303.704987, 326.605011, 349.404999, // FTT
1271  375.0, // EPD
1272  715.0, //ECAL
1273  807.0 // HCAL
1274  };
1275 
1276  // Note: as discussed, after verification storage of the projections
1277  // @ the FST and FTT may no longer be needed, not saved in e.g. MuDst
1278 
1279  // this should always be the case, but being careful
1280  if (gGeoManager) {
1281  FwdGeomUtils fwdGeoUtils( gGeoManager );
1282 
1283  // get the z-locations from geometry model and fallback to the defaults
1284  auto fstZ = fwdGeoUtils.fstZ( {151.750000, 165.248001, 178.781006} );
1285  auto fttZ = fwdGeoUtils.fttZ( {280.904999, 303.704987, 326.605011, 349.404999} );
1286 
1287  // copy new values into the zPlanes vector
1288  std::copy( fstZ.begin(), fstZ.end(), zPlanes.begin()+1 );
1289  std::copy( fttZ.begin(), fttZ.end(), zPlanes.begin()+4 );
1290  }
1291 
1292  // Match these to the z-planes above
1293  const int FST = kFstId;
1294  const int FTT = kFttId;
1295  vector<int> detMap = {
1296  kTpcId,
1297  FST, FST, FST,
1298  FTT, FTT, FTT, FTT,
1299  kFcsPresId,
1300  kFcsWcalId,
1301  kFcsHcalId
1302  };
1303 
1304  size_t zIndex = 0;
1305  int detIndex = 0;
1306  for ( float z : zPlanes ){
1307  detIndex = detMap[ zIndex];
1308  // LOG_DEBUG << "Calculating Projection for detId=" << detIndex << " @ z=" << z << endm;
1309  TVector3 mom(0, 0, 0);
1310  float cov[9];
1311 
1312  TVector3 tv3(0, 0, 0);
1313  if ( detIndex != kFcsHcalId ){
1314  tv3 = ObjExporter::trackPosition( track, z, cov, mom );
1315  } else {
1316  // use a straight line projection to HCAL since GenFit cannot handle long projections
1317  tv3 = ObjExporter::projectAsStraightLine( track, 575.0, 625.0, z, cov, mom );
1318  }
1319  fwdTrack->mProjections.push_back( StFwdTrackProjection( detIndex, StThreeVectorF( tv3.X(), tv3.Y(), tv3.Z() ), StThreeVectorF( mom.X(), mom.Y(), mom.Z() ), cov) );
1320 
1321  // // Add Proj info to TTree
1322  mTreeData.tprojX.push_back( tv3.X() );
1323  mTreeData.tprojY.push_back( tv3.Y() );
1324  mTreeData.tprojZ.push_back( tv3.Z() );
1325  mTreeData.tprojIdD.push_back( detIndex );
1326  mTreeData.tprojIdT.push_back( indexTrack );
1327 
1328  zIndex++;
1329  }
1330 
1331  return fwdTrack;
1332 }
1333 
1334 void StFwdTrackMaker::FillEvent() {
1335  LOG_DEBUG << "StFwdTrackMaker::FillEvent()" << endm;
1336  // Now fill StEvent
1337  StEvent *stEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
1338 
1339  // FillEvent();
1340  StFwdTrackCollection * ftc = stEvent->fwdTrackCollection();
1341  if ( !ftc ){
1342  LOG_INFO << "Creating the StFwdTrackCollection" << endm;
1343  ftc = new StFwdTrackCollection();
1344  stEvent->setFwdTrackCollection( ftc );
1345  }
1346 
1347  mTreeData.tprojN = 0;
1348  mTreeData.tprojX.clear();
1349  mTreeData.tprojY.clear();
1350  mTreeData.tprojZ.clear();
1351  mTreeData.tprojIdD.clear();
1352  mTreeData.tprojIdT.clear();
1353  size_t indexTrack = 0;
1354  for ( auto gtr : mForwardTracker->getTrackResults() ) {
1355  StFwdTrack* fwdTrack = makeStFwdTrack( gtr, indexTrack );
1356  if (nullptr == fwdTrack)
1357  continue;
1358  ftc->addTrack( fwdTrack );
1359  }
1360 
1361 
1362  mTreeData.tprojN = mTreeData.tprojX.size();
1363  LOG_DEBUG << "StFwdTrackCollection has " << ftc->numberOfTracks() << " tracks now" << endm;
1364  ProcessFwdTracks();
1365 }
1366 
1367 void StFwdTrackMaker::FillTrackDeltas(){
1368  LOG_DEBUG << "Filling Track Deltas for Alignment" << endm;
1369  const auto &fittedTracks = mForwardTracker -> getTrackResults();
1370 
1371  for ( size_t i = 0; i < fittedTracks.size(); i++ ){
1372  auto st = fittedTracks[i].trackSeed;
1373  auto gt = fittedTracks[i].track;
1374 
1375  if (fittedTracks[i].isFitConvergedFully == false){
1376  LOG_DEBUG << "Skipping track, failed fit" << endm;
1377  continue;
1378  }
1379 
1380  for ( KiTrack::IHit * hit : st ){
1381  TVector3 htv3(hit->getX(), hit->getY(), hit->getZ());
1382 
1383  auto ttv3 = ObjExporter::trackPosition( gt, htv3.Z() );
1384 
1385  mTreeData.thdX.push_back( (ttv3.X() - htv3.X()) );
1386  mTreeData.thdY.push_back( (ttv3.Y() - htv3.Y()) );
1387  mTreeData.thaX.push_back( htv3.X() );
1388  mTreeData.thaY.push_back( htv3.Y() );
1389  mTreeData.thaZ.push_back( htv3.Z() );
1390  mTreeData.thdN++;
1391  }
1392  }
1393 } //FillTrackDeltas
1394 
1395 void StFwdTrackMaker::FitVertex(){
1396  const auto &genfitTracks = mForwardTracker -> globalTracks();
1397 
1398  if ( genfitTracks.size() >= 2 ){
1399  genfit::GFRaveVertexFactory gfrvf;
1400 
1401  TMatrixDSym bscm(3);
1402  const double bssXY = 2.0;
1403  bscm(0, 0) = bssXY*bssXY;
1404  bscm(1, 1) = bssXY*bssXY;
1405  bscm(2, 2) = 50.5 * 50.5;
1406  gfrvf.setBeamspot( TVector3( 0, 0, 0 ), bscm );
1407  // std::vector< genfit::GFRaveVertex * > vertices;
1408  const auto &genfitTracks = mForwardTracker -> globalTracks();
1409  mRaveVertices.clear();
1410  gfrvf.findVertices( &mRaveVertices, genfitTracks, false );
1411 
1412  LOG_DEBUG << "mRaveVertices.size() = " << mRaveVertices.size() << endm;
1413 
1414  for ( auto vert : mRaveVertices ){
1415  LOG_DEBUG << TString::Format( "RAVE vertex @(%f, %f, %f)\n\n", vert->getPos().X(), vert->getPos().Y(), vert->getPos().Z() ) << endm;
1416  }
1417  }
1418 }
1419 
1420 void StFwdTrackMaker::FillTTree(){
1421 
1422  St_g2t_vertex *g2t_vertex = (St_g2t_vertex *)GetDataSet("geant/g2t_vertex");
1423  if (mGenTree) {
1424 
1425  // VERTICES
1426  if ( g2t_vertex ){
1427  mTreeData.vmcN = g2t_vertex->GetNRows();
1428  if ( (unsigned)mTreeData.vmcN >= MAX_TREE_ELEMENTS ) mTreeData.vmcN = MAX_TREE_ELEMENTS;
1429  LOG_INFO << "Saving " << mTreeData.vmcN << " vertices in TTree" << endm;
1430  for ( int i = 0; i < mTreeData.vmcN; i++ ){
1431  g2t_vertex_st *vert = (g2t_vertex_st*)g2t_vertex->At(i);
1432  mTreeData.vmcX.push_back( vert->ge_x[0] );
1433  mTreeData.vmcY.push_back( vert->ge_x[1] );
1434  mTreeData.vmcZ.push_back( vert->ge_x[2] );
1435  }
1436  }
1437 
1438  // RAVE RECO VERTICES
1439  mTreeData.vrcN = mRaveVertices.size();
1440  if ( (unsigned)mTreeData.vrcN >= MAX_TREE_ELEMENTS ) mTreeData.vrcN = MAX_TREE_ELEMENTS;
1441  LOG_INFO << "Saving " << mTreeData.vrcN << " RAVE vertices in TTree" << endm;
1442  for ( int i = 0; i < mTreeData.vrcN; i++ ) {
1443  auto vert = mRaveVertices[i];
1444  mTreeData.vrcX.push_back( vert->getPos().X() );
1445  mTreeData.vrcY.push_back( vert->getPos().Y() );
1446  mTreeData.vrcZ.push_back( vert->getPos().Z() );
1447  }
1448 
1449 
1450 
1451 
1452  if (mForwardTracker->getSaveCriteriaValues() && mTreeData.saveCrit ) {
1453  for (auto crit : mForwardTracker->getTwoHitCriteria()) {
1454  string name = crit->getName();
1455 
1456  // special, save all hit info for this one
1457 
1458 
1459  if ( name == "Crit2_BDT" ){
1460  mTreeData.Crits["Crit2_BDT_DeltaPhi"].clear();
1461  mTreeData.Crits["Crit2_BDT_DeltaRho"].clear();
1462  mTreeData.Crits["Crit2_BDT_RZRatio"].clear();
1463  mTreeData.Crits["Crit2_BDT_StraightTrackRatio"].clear();
1464 
1465  for (auto kv : mForwardTracker->getCriteriaAllValues(name)) {
1466  mTreeData.Crits["Crit2_BDT_DeltaPhi"].push_back( kv["Crit2_BDT_DeltaPhi"] );
1467  mTreeData.Crits["Crit2_BDT_DeltaRho"].push_back( kv["Crit2_BDT_DeltaRho"] );
1468  mTreeData.Crits["Crit2_BDT_RZRatio"].push_back( kv["Crit2_BDT_RZRatio"] );
1469  mTreeData.Crits["Crit2_BDT_StraightTrackRatio"].push_back( kv["Crit2_BDT_StraightTrackRatio"] );
1470  }
1471 
1472  }
1473 
1474  if ( name == "Crit2_RZRatio" ){
1475  LOG_INFO << "allValues.size() = " << mForwardTracker->getCriteriaAllValues(name).size() << " == " << mForwardTracker->getCriteriaTrackIds(name).size() << endm;
1476  assert( mForwardTracker->getCriteriaAllValues(name).size() == mForwardTracker->getCriteriaTrackIds(name).size() && " Crit lengths must be equal" );
1477  mTreeData.Crits["Crit2_RZRatio_x1"].clear();
1478  mTreeData.Crits["Crit2_RZRatio_y1"].clear();
1479  mTreeData.Crits["Crit2_RZRatio_z1"].clear();
1480  mTreeData.Crits["Crit2_RZRatio_x2"].clear();
1481  mTreeData.Crits["Crit2_RZRatio_y2"].clear();
1482  mTreeData.Crits["Crit2_RZRatio_z2"].clear();
1483 
1484  mTreeData.CritTrackIds["Crit2_RZRatio_h1"].clear();
1485  mTreeData.CritTrackIds["Crit2_RZRatio_h2"].clear();
1486  mTreeData.CritTrackIds["Crit2_RZRatio_h3"].clear();
1487 
1488 
1489  for (auto kv : mForwardTracker->getCriteriaAllValues(name)) {
1490  mTreeData.Crits["Crit2_RZRatio_x1"].push_back( kv["Crit2_RZRatio_x1"] );
1491  mTreeData.Crits["Crit2_RZRatio_y1"].push_back( kv["Crit2_RZRatio_y1"] );
1492  mTreeData.Crits["Crit2_RZRatio_z1"].push_back( kv["Crit2_RZRatio_z1"] );
1493 
1494  mTreeData.Crits["Crit2_RZRatio_x2"].push_back( kv["Crit2_RZRatio_x2"] );
1495  mTreeData.Crits["Crit2_RZRatio_y2"].push_back( kv["Crit2_RZRatio_y2"] );
1496  mTreeData.Crits["Crit2_RZRatio_z2"].push_back( kv["Crit2_RZRatio_z2"] );
1497 
1498  mTreeData.CritTrackIds["Crit2_RZRatio_h1"].push_back( kv["Crit2_RZRatio_h1"] );
1499  mTreeData.CritTrackIds["Crit2_RZRatio_h2"].push_back( kv["Crit2_RZRatio_h2"] );
1500  mTreeData.CritTrackIds["Crit2_RZRatio_h3"].push_back( -1 );
1501  }
1502  }
1503 
1504 
1505  LOG_DEBUG << "Saving Criteria values from " << name << " in TTree" << endm;
1506  mTreeData.Crits[name].clear();
1507  mTreeData.CritTrackIds[name].clear();
1508  // copy by value so ROOT doesnt get lost (uses pointer to vector)
1509  for (float v : mForwardTracker->getCriteriaValues(name)) {
1510  mTreeData.Crits[name].push_back(v);
1511  }
1512  for (int v : mForwardTracker->getCriteriaTrackIds(name)) {
1513  mTreeData.CritTrackIds[name].push_back(v);
1514  }
1515  }
1516 
1517  // three hit criteria
1518  for (auto crit : mForwardTracker->getThreeHitCriteria()) {
1519  string name = crit->getName();
1520 
1521  // special, save all hit info for this one
1522  if ( name == "Crit2_RZRatio" ){
1523  LOG_INFO << "allValues.size() = " << mForwardTracker->getCriteriaAllValues(name).size() << " == " << mForwardTracker->getCriteriaTrackIds(name).size() << endm;
1524  assert( mForwardTracker->getCriteriaAllValues(name).size() == mForwardTracker->getCriteriaTrackIds(name).size() && " Crit lengths must be equal" );
1525 
1526  mTreeData.CritTrackIds["Crit2_RZRatio_h1"].clear();
1527  mTreeData.CritTrackIds["Crit2_RZRatio_h2"].clear();
1528  mTreeData.CritTrackIds["Crit2_RZRatio_h3"].clear();
1529 
1530  for (auto kv : mForwardTracker->getCriteriaAllValues(name)) {
1531  mTreeData.CritTrackIds["Crit2_RZRatio_h1"].push_back( kv["Crit2_RZRatio_h1"] );
1532  mTreeData.CritTrackIds["Crit2_RZRatio_h2"].push_back( kv["Crit2_RZRatio_h2"] );
1533  mTreeData.CritTrackIds["Crit2_RZRatio_h3"].push_back( kv["Crit2_RZRatio_h3"] );
1534  }
1535  }
1536 
1537 
1538  LOG_DEBUG << "Saving Criteria values from " << name << " in TTree" << endm;
1539  mTreeData.Crits[name].clear();
1540  mTreeData.CritTrackIds[name].clear();
1541  // copy by value so ROOT doesnt get lost (uses pointer to vector)
1542  for (float v : mForwardTracker->getCriteriaValues(name)) {
1543  mTreeData.Crits[name].push_back(v);
1544  }
1545  for (int v : mForwardTracker->getCriteriaTrackIds(name)) {
1546  mTreeData.CritTrackIds[name].push_back(v);
1547  }
1548  }
1549 
1550  // clear them
1551  mForwardTracker->clearSavedCriteriaValues();
1552  }
1553 
1554  // SAVE RECO tracks
1555 
1556  mTreeData.rcN = 0;
1557  const auto &fittedTracks = mForwardTracker -> getTrackResults();
1558 
1559  LOG_INFO << "There are " << fittedTracks.size() << " seed tracks to save" << endm;
1560  size_t maxToSave = fittedTracks.size();
1561  if (maxToSave >= 200) {
1562  maxToSave = 0;
1563  LOG_INFO << "More than 200 tracks , not saving unfit tracks" << endm;
1564  }
1565 
1566  for ( size_t i = 0; i < maxToSave; i++ ){
1567  if ( i >= MAX_TREE_ELEMENTS ){
1568  LOG_WARN << "Truncating Reco tracks in TTree output" << endm;
1569  break;
1570  }
1571 
1572  int idt = 0;
1573  double qual = 0;
1574  idt = MCTruthUtils::dominantContribution(fittedTracks[i].trackSeed, qual);
1575 
1576  if ( fittedTracks[i].track == nullptr || fittedTracks[i].trackRep == nullptr ) {
1577  LOG_INFO << "Skip saving null track" << endm;
1578  continue;
1579  }
1580 
1581  if ( fittedTracks[i].isFitConverged == false ){
1582  LOG_INFO << "Skip saving track where fit did not converge" << endm;
1583  continue;
1584  }
1585 
1586 
1587  mTreeData.rcQuality.push_back( qual );
1588  mTreeData.rcTrackId.push_back( idt );
1589 
1590  mTreeData.rcCharge.push_back( fittedTracks[i].charge );
1591  mTreeData.rcPt.push_back( fittedTracks[i].momentum.Pt() );
1592  mTreeData.rcEta.push_back( fittedTracks[i].momentum.Eta() );
1593  mTreeData.rcPhi.push_back( fittedTracks[i].momentum.Phi() );
1594 
1595  mTreeData.rcNumPV.push_back( fittedTracks[i].nPV );
1596  mTreeData.rcNumFTT.push_back( fittedTracks[i].nFTT );
1597  mTreeData.rcNumFST.push_back( fittedTracks[i].nFST );
1598 
1599  mTreeData.rcN ++;
1600  }
1601  LOG_INFO << "Filling TTree" << endm;
1602  mTree->Fill();
1603  } // if mGenTree
1604 }
1605 
1606 
1607 //________________________________________________________________________
1608 void StFwdTrackMaker::Clear(const Option_t *opts) {
1609  LOG_DEBUG << "StFwdTrackMaker::CLEAR" << endm;
1610  mForwardData->clear();
1611 
1612  if (mGenTree){
1613  mTreeData.thdN = mTreeData.fttN = mTreeData.rcN = mTreeData.mcN = mTreeData.vmcN = mTreeData.vrcN = mTreeData.fcsN = 0;
1614  mTreeData.fttX.clear();
1615  mTreeData.fttY.clear();
1616  mTreeData.fttZ.clear();
1617  mTreeData.fttTrackId.clear();
1618  mTreeData.fttVolumeId.clear();
1619  mTreeData.fttPt.clear();
1620  mTreeData.fttVertexId.clear();
1621 
1622  mTreeData.fstX.clear();
1623  mTreeData.fstY.clear();
1624  mTreeData.fstZ.clear();
1625  mTreeData.fstTrackId.clear();
1626 
1627  mTreeData.fcsX.clear();
1628  mTreeData.fcsY.clear();
1629  mTreeData.fcsZ.clear();
1630  mTreeData.fcsDet.clear();
1631 
1632  mTreeData.rcPt.clear();
1633  mTreeData.rcEta.clear();
1634  mTreeData.rcPhi.clear();
1635  mTreeData.rcQuality.clear();
1636  mTreeData.rcTrackId.clear();
1637  mTreeData.rcNumFST.clear();
1638  mTreeData.rcCharge.clear();
1639  mTreeData.rcNumFTT.clear();
1640  mTreeData.rcNumPV.clear();
1641 
1642 
1643  mTreeData.mcPt.clear();
1644  mTreeData.mcEta.clear();
1645  mTreeData.mcPhi.clear();
1646  mTreeData.mcVertexId.clear();
1647  mTreeData.mcCharge.clear();
1648  mTreeData.vmcX.clear();
1649  mTreeData.vmcY.clear();
1650  mTreeData.vmcZ.clear();
1651 
1652  mTreeData.tprojX.clear();
1653  mTreeData.tprojY.clear();
1654  mTreeData.tprojZ.clear();
1655  mTreeData.tprojPx.clear();
1656  mTreeData.tprojPy.clear();
1657  mTreeData.tprojPz.clear();
1658  mTreeData.vrcX.clear();
1659  mTreeData.vrcY.clear();
1660  mTreeData.vrcZ.clear();
1661  mTreeData.thdX.clear();
1662  mTreeData.thdY.clear();
1663  mTreeData.thaX.clear();
1664  mTreeData.thaY.clear();
1665  mTreeData.thaZ.clear();
1666 
1667  mTreeData.thdX.clear();
1668  mTreeData.thdY.clear();
1669  mTreeData.thaX.clear();
1670  mTreeData.thaY.clear();
1671  mTreeData.thaZ.clear();
1672 
1673  mTreeData.fttN = 0;
1674  mTreeData.fstN = 0;
1675  mTreeData.rcN = 0;
1676  mTreeData.mcN = 0;
1677  mTreeData.vmcN = 0;
1678  mTreeData.tprojN = 0;
1679  mTreeData.vrcN = 0;
1680  mTreeData.thdN = 0;
1681  }
1682 }
1683 //________________________________________________________________________
1684 void StFwdTrackMaker::ProcessFwdTracks( ){
1685  // This is an example of how to process fwd track collection
1686  LOG_DEBUG << "StFwdTrackMaker::ProcessFwdTracks" << endm;
1687  StEvent *stEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
1688  StFwdTrackCollection * ftc = stEvent->fwdTrackCollection();
1689  for ( auto fwdTrack : ftc->tracks() ){
1690  LOG_DEBUG << TString::Format("StFwdTrack[ nProjections=%lu, nFTTSeeds=%lu, nFSTSeeds=%lu, mPt=%f ]", fwdTrack->mProjections.size(), fwdTrack->mFTTPoints.size(), fwdTrack->mFSTPoints.size(), fwdTrack->momentum().perp()) << endm;
1691  for ( auto proj : fwdTrack->mProjections ) {
1692  LOG_DEBUG << TString::Format("Proj[ %d, %f, %f, %f ]", proj.mDetId, proj.mXYZ.x(), proj.mXYZ.y(), proj.mXYZ.z() ) << endm;
1693  }
1694  }
1695 }
1696 
1697 
1698 std::string StFwdTrackMaker::defaultConfigIdealSim = R"( <?xml version="1.0" encoding="UTF-8"?> <config> <Output url="fwdTrackMaker_ideal_sim.root" /> <Source ftt="GEANT" /> <TrackFitter refit="true" mcSeed="true" > <Vertex sigmaXY="0.001" sigmaZ="0.01" includeInFit="true" smearMcVertex="true" /> </TrackFitter> </config> )";
1699 
1700 
1701 
1702 std::string StFwdTrackMaker::defaultConfigData = R"( <?xml version="1.0" encoding="UTF-8"?> <config> <Output url="stfwdtrackmaker_data.root" /> <Source ftt="DATA" /> <SiRasterizer r="3.0" phi="0.004" /> <TrackFinder nIterations="1"> <Iteration nPhiSlices="32" > <!-- Options for first iteration --> <SegmentBuilder> <Criteria name="Crit2_RZRatio" min="0" max="1.20" /> <Criteria name="Crit2_DeltaRho" min="-50" max="50.9"/> <Criteria name="Crit2_DeltaPhi" min="0" max="30.0" /> <Criteria name="Crit2_StraightTrackRatio" min="0.01" max="5.85"/> </SegmentBuilder> <ThreeHitSegments> <Criteria name="Crit3_3DAngle" min="0" max="30" /> <Criteria name="Crit3_PT" min="0" max="100" /> <Criteria name="Crit3_ChangeRZRatio" min="0.8" max="1.21" /> <Criteria name="Crit3_2DAngle" min="0" max="30" /> </ThreeHitSegments> </Iteration> <Connector distance="1"/> <SubsetNN active="true" min-hits-on-track="3" > <!-- <InitialTemp>2.1</InitialTemp> --> <!-- <InfTemp>0.1</InfTemp> --> <Omega>0.99</Omega> <StableThreshold>0.001</StableThreshold> </SubsetNN> <HitRemover active="false" /> </TrackFinder> <TrackFitter refitSi="true" mcSeed="false" zeroB="true"> <Vertex sigmaXY="0.01" sigmaZ="0.01" includeInFit="true" smearMcVertex="false" /> </TrackFitter> </config> )";
FwdHit
Definition: FwdHit.h:68
FwdTrackerConfig::childrenOf
std::vector< std::string > childrenOf(std::string path) const
list the paths of children nodes for a given node
Definition: FwdTrackerConfig.h:257
StFcsDb::getStarXYZfromColumnRow
StThreeVectorD getStarXYZfromColumnRow(int det, float col, float row, float FcsZ=-1.0) const
get the STAR frame cooridnates from other way
Definition: StFcsDb.cxx:863
McTrack
Definition: FwdHit.h:35
FwdTrackerConfig::load
void load(std::string filename, bool asString=false)
Main setup routine Loads the given XML file (or string) and maps it.
Definition: FwdTrackerConfig.h:305
StFcsDb::getEPDfromId
void getEPDfromId(int det, int id, int &pp, int &tt)
Get FCS&#39;s EPD map foom EPD mapping.
Definition: StFcsDb.cxx:1484
FwdTrackerConfig::get
T get(std::string path, T dv) const
template function for getting any type that can be converted from string via stringstream ...
Definition: FwdTrackerConfig.h:190
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
StMaker::GetIventNumber
virtual Int_t GetIventNumber() const
Returns the current event number.
Definition: StMaker.cxx:1033
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5