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 #include <ConstField.h>
2 #include <Exception.h>
3 #include <FieldManager.h>
4 #include <KalmanFitterRefTrack.h>
5 #include <StateOnPlane.h>
6 #include <Track.h>
7 #include <TrackPoint.h>
8 
9 #include <MaterialEffects.h>
10 #include <RKTrackRep.h>
11 #include <TGeoMaterialInterface.h>
12 
13 #include <EventDisplay.h>
14 
15 #include <HelixTrackModel.h>
16 #include <MeasurementCreator.h>
17 
18 #include <TDatabasePDG.h>
19 #include <TEveManager.h>
20 #include <TGeoManager.h>
21 #include <TRandom.h>
22 #include <TVector3.h>
23 #include <vector>
24 
25 #include "TDatabasePDG.h"
26 #include <TMath.h>
27 
28 
29 
30 
31 int main() {
32 
33  gRandom->SetSeed(14);
34 
35  // init MeasurementCreator
36  genfit::MeasurementCreator measurementCreator;
37 
38 
39  // init geometry and mag. field
40  new TGeoManager("Geometry", "Geane geometry");
41  TGeoManager::Import("genfitGeom.root");
42  genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., 15.)); // 15 kGauss
44 
45 
46  // init event display
48 
49 
50  // init fitter
52 
53  // main loop
54  for (unsigned int iEvent=0; iEvent<100; ++iEvent){
55 
56  // true start values
57  TVector3 pos(0, 0, 0);
58  TVector3 mom(1.,0,0);
59  mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
60  mom.SetTheta(gRandom->Uniform(0.4*TMath::Pi(),0.6*TMath::Pi()));
61  mom.SetMag(gRandom->Uniform(0.2, 1.));
62 
63 
64  // helix track model
65  const int pdg = 13; // particle pdg code
66  const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
67  genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
68  measurementCreator.setTrackModel(helix);
69 
70 
71  unsigned int nMeasurements = gRandom->Uniform(5, 15);
72 
73 
74  // smeared start values
75  const bool smearPosMom = true; // init the Reps with smeared pos and mom
76  const double posSmear = 0.1; // cm
77  const double momSmear = 3. /180.*TMath::Pi(); // rad
78  const double momMagSmear = 0.1; // relative
79 
80  TVector3 posM(pos);
81  TVector3 momM(mom);
82  if (smearPosMom) {
83  posM.SetX(gRandom->Gaus(posM.X(),posSmear));
84  posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
85  posM.SetZ(gRandom->Gaus(posM.Z(),posSmear));
86 
87  momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
88  momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
89  momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
90  }
91  // approximate covariance
92  TMatrixDSym covM(6);
93  double resolution = 0.01;
94  for (int i = 0; i < 3; ++i)
95  covM(i,i) = resolution*resolution;
96  for (int i = 3; i < 6; ++i)
97  covM(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
98 
99 
100  // trackrep
101  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
102 
103  // smeared start state
104  genfit::MeasuredStateOnPlane stateSmeared(rep);
105  stateSmeared.setPosMomCov(posM, momM, covM);
106 
107 
108  // create track
109  TVectorD seedState(6);
110  TMatrixDSym seedCov(6);
111  stateSmeared.get6DStateCov(seedState, seedCov);
112  genfit::Track fitTrack(rep, seedState, seedCov);
113 
114 
115  // create random measurement types
116  std::vector<genfit::eMeasurementType> measurementTypes;
117  for (unsigned int i = 0; i < nMeasurements; ++i)
118  measurementTypes.push_back(genfit::eMeasurementType(gRandom->Uniform(8)));
119 
120 
121  // create smeared measurements and add to track
122  try{
123  for (unsigned int i=0; i<measurementTypes.size(); ++i){
124  std::vector<genfit::AbsMeasurement*> measurements = measurementCreator.create(measurementTypes[i], i*5.);
125  fitTrack.insertPoint(new genfit::TrackPoint(measurements, &fitTrack));
126  }
127  }
128  catch(genfit::Exception& e){
129  std::cerr<<"Exception, next track"<<std::endl;
130  std::cerr << e.what();
131  continue;
132  }
133 
134  //check
135  fitTrack.checkConsistency();
136 
137  // do the fit
138  fitter->processTrack(&fitTrack);
139 
140  //check
141  fitTrack.checkConsistency();
142 
143 
144  if (iEvent < 1000) {
145  // add track to event display
146  display->addEvent(&fitTrack);
147  }
148 
149 
150 
151  }// end loop over events
152 
153  delete fitter;
154 
155  // open event display
156  display->open();
157 
158 }
159 
160