Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
main.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file main.cc
1 //
2 // Write fit results to tree, read again, compare.
3 //
4 
5 #include <ConstField.h>
6 #include <Exception.h>
7 #include <FieldManager.h>
8 #include <KalmanFitterRefTrack.h>
9 #include <StateOnPlane.h>
10 #include <Track.h>
11 #include <TrackPoint.h>
12 
13 #include <MaterialEffects.h>
14 #include <RKTrackRep.h>
15 #include <TGeoMaterialInterface.h>
16 
17 #include <EventDisplay.h>
18 
19 #include <HelixTrackModel.h>
20 #include <MeasurementCreator.h>
21 
22 #include <TDatabasePDG.h>
23 #include <TEveManager.h>
24 #include <TGeoManager.h>
25 #include <TRandom.h>
26 #include <TVector3.h>
27 #include <vector>
28 
29 #include "TDatabasePDG.h"
30 #include <TMath.h>
31 #include <TFile.h>
32 #include <TTree.h>
33 
34 
35 #define FILENAME "/tmp/streamerTest.root"
36 
37 constexpr bool verbose = false;
38 
40 {
41  TFile *f = TFile::Open(FILENAME, "RECREATE");
42  f->cd();
43  genfit::Track *t = new genfit::Track();
44  try {
45  t->checkConsistency();
46  } catch (genfit::Exception&) {
47  return false;
48  }
49 
50  t->Write("direct");
51  f->Close();
52  delete t;
53  delete f;
54 
55  f = TFile::Open(FILENAME, "READ");
56  t = (genfit::Track*)f->Get("direct");
57 
58  bool result = false;
59  try {
60  t->checkConsistency();
61  result = true;
62  } catch (genfit::Exception&) {
63  result = false;
64  }
65  delete t;
66  delete f;
67  return result;
68 }
69 
70 
71 int main() {
72  if (!emptyTrackTest()) {
73  std::cout << "emptyTrackTest failed." << std::endl;
74  return 1;
75  }
76 
77 
78  // prepare output tree for Tracks
79  // std::unique_ptr<genfit::Track> fitTrack(new genfit::Track());
80  genfit::Track* fitTrack = new genfit::Track();
81  TVectorD stateFinal;
82  TMatrixDSym covFinal;
83  genfit::DetPlane planeFinal;
84 
85  TFile* fOut = new TFile(FILENAME, "RECREATE");
86  fOut->cd();
87  TTree* tResults = new TTree("tResults", "results from track fit");
88  // it is important to set the splitLevel to -1, because some of the genfit classes use custom streamers
89  tResults->Branch("gfTrack", "genfit::Track", &fitTrack, 32000, -1);
90  tResults->Branch("stateFinal", &stateFinal);
91  tResults->Branch("covFinal", &covFinal, 32000, -1);
92  tResults->Branch("planeFinal", &planeFinal, 32000, -1);
93 
94 
95 
96 
97 
98  gRandom->SetSeed(14);
99 
100  // init MeasurementCreator
101  genfit::MeasurementCreator measurementCreator;
102 
103 
104  // init geometry and mag. field
105  new TGeoManager("Geometry", "Geane geometry");
106  TGeoManager::Import("genfitGeom.root");
107  genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., 15.)); // 15 kGauss
109 
110 
111  // init fitter
113 
114  unsigned int nEvents(100);
115 
116  // main loop
117  for (unsigned int iEvent=0; iEvent<nEvents; ++iEvent){
118 
119  // true start values
120  TVector3 pos(0, 0, 0);
121  TVector3 mom(1.,0,0);
122  mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
123  mom.SetTheta(gRandom->Uniform(0.4*TMath::Pi(),0.6*TMath::Pi()));
124  mom.SetMag(gRandom->Uniform(0.2, 1.));
125 
126 
127  // helix track model
128  const int pdg = 13; // particle pdg code
129  const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
130  genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
131  measurementCreator.setTrackModel(helix);
132 
133 
134  unsigned int nMeasurements = gRandom->Uniform(5, 15);
135 
136 
137  // smeared start values
138  const bool smearPosMom = true; // init the Reps with smeared pos and mom
139  const double posSmear = 0.1; // cm
140  const double momSmear = 3. /180.*TMath::Pi(); // rad
141  const double momMagSmear = 0.1; // relative
142 
143  TVector3 posM(pos);
144  TVector3 momM(mom);
145  if (smearPosMom) {
146  posM.SetX(gRandom->Gaus(posM.X(),posSmear));
147  posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
148  posM.SetZ(gRandom->Gaus(posM.Z(),posSmear));
149 
150  momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
151  momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
152  momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
153  }
154  // approximate covariance
155  TMatrixDSym covM(6);
156  double resolution = 0.01;
157  for (int i = 0; i < 3; ++i)
158  covM(i,i) = resolution*resolution;
159  for (int i = 3; i < 6; ++i)
160  covM(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
161 
162 
163  // trackrep
164  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
165 
166  // smeared start state
167  genfit::MeasuredStateOnPlane stateSmeared(rep);
168  rep->setPosMomCov(stateSmeared, posM, momM, covM);
169 
170 
171  // create track
172  TVectorD seedState(6);
173  TMatrixDSym seedCov(6);
174  rep->get6DStateCov(stateSmeared, seedState, seedCov);
175  fitTrack = new genfit::Track(rep, seedState, seedCov);
176 
177 
178  // create random measurement types
179  std::vector<genfit::eMeasurementType> measurementTypes;
180  for (unsigned int i = 0; i < nMeasurements; ++i)
181  measurementTypes.push_back(genfit::eMeasurementType(gRandom->Uniform(8)));
182 
183 
184  // create smeared measurements and add to track
185  try{
186  for (unsigned int i=0; i<measurementTypes.size(); ++i){
187  std::vector<genfit::AbsMeasurement*> measurements = measurementCreator.create(measurementTypes[i], i*5.);
188  genfit::TrackPoint* tp = new genfit::TrackPoint(measurements, fitTrack);
189  // test scatterers
190  genfit::ThinScatterer* sc = new genfit::ThinScatterer(genfit::SharedPlanePtr(new genfit::DetPlane(TVector3(1,1,1), TVector3(1,1,1))),
191  genfit::Material(1,2,3,4,5));
192  tp->setScatterer(sc);
193 
194  fitTrack->insertPoint(tp);
195  }
196  }
197  catch(genfit::Exception& e){
198  std::cerr<<"Exception, next track"<<std::endl;
199  std::cerr << e.what();
200  delete fitTrack;
201  fitTrack = 0;
202  continue;
203  }
204 
205  fitTrack->checkConsistency();
206 
207  // do the fit
208  fitter->processTrack(fitTrack);
209 
210  fitTrack->checkConsistency();
211 
212 
213  stateFinal.ResizeTo(fitTrack->getFittedState().getState());
214  stateFinal = fitTrack->getFittedState().getState();
215  covFinal.ResizeTo(fitTrack->getFittedState().getCov());
216  covFinal = fitTrack->getFittedState().getCov();
217  planeFinal = *(fitTrack->getFittedState().getPlane());
218 
219  switch (iEvent % 5) {
220  case 0:
221  break;
222  case 1:
223  fitTrack->prune("FL");
224  break;
225  case 2:
226  fitTrack->prune("W");
227  break;
228  case 3:
229  fitTrack->prune("RC");
230  break;
231  case 4:
232  fitTrack->prune("U");
233  break;
234  }
235 
236  tResults->Fill();
237 
238  //fitTrack->Print();
239 
240  delete fitTrack;
241  fitTrack = 0;
242 
243  }// end loop over events
244 
245  delete fitter;
246 
247 
248 
249 
250  fOut->Write();
251  delete fOut;
252 
253  fOut = TFile::Open(FILENAME, "READ");
254  fOut->GetObject("tResults", tResults);
255  TVectorD* pState = 0;
256  tResults->SetBranchAddress("stateFinal", &pState);
257  TMatrixDSym* pMatrix = 0;
258  tResults->SetBranchAddress("covFinal", &pMatrix);
259  genfit::DetPlane* plane = 0;
260  tResults->SetBranchAddress("planeFinal", &plane);
261  tResults->SetBranchAddress("gfTrack", &fitTrack);
262 
263  int fail(0);
264 
265  for (Long_t nEntry = 0; nEntry < tResults->GetEntries(); ++nEntry) {
266  tResults->GetEntry(nEntry);
267 
268  try {
269  fitTrack->checkConsistency();
270  } catch (genfit::Exception& e) {
271  std::cout << e.getExcString() << std::endl;
272  return 1;
273  }
274 
275  try {
276  if (*pState == fitTrack->getFittedState().getState() &&
277  *pMatrix == fitTrack->getFittedState().getCov() &&
278  *plane == *(fitTrack->getFittedState().getPlane())) {
279  // track ok
280  }
281  else {
282  if (verbose) {
283  std::cout << "stored track not equal, small differences can occur if some info has been pruned." << std::endl;
284  pState->Print();
285  fitTrack->getFittedState().getState().Print();
286  pMatrix->Print();
287  fitTrack->getFittedState().getCov().Print();
288  plane->Print();
289  fitTrack->getFittedState().getPlane()->Print();
290  }
291  ++fail;
292  //return 1;
293  }
294  }
295  catch (genfit::Exception& e) {
296  if (verbose) {
297  std::cerr << e.what();
298  }
299  return 1;
300  }
301  }
302  std::cout << nEvents - fail << " stored tracks are identical to fitted tracks, as far as tested." << std::endl;
303  std::cout << fail << " tracks were not identical to fitted tracks, as far as tested." << std::endl;
304  delete fitTrack;
305  std::cout << "deleteing didn't segfault" << std::endl;
306 
307  return 0;
308 }