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 <TrackCand.h>
8 #include <TrackPoint.h>
9 
10 #include <MeasurementProducer.h>
11 #include <MeasurementFactory.h>
12 
15 
16 #include <MaterialEffects.h>
17 #include <RKTrackRep.h>
18 #include <TGeoMaterialInterface.h>
19 
20 #include <EventDisplay.h>
21 
22 #include <HelixTrackModel.h>
23 #include <MeasurementCreator.h>
24 
25 #include <TDatabasePDG.h>
26 #include <TEveManager.h>
27 #include <TGeoManager.h>
28 #include <TRandom.h>
29 #include <TVector3.h>
30 #include <vector>
31 
32 #include "TDatabasePDG.h"
33 #include <TMath.h>
34 
35 
36 
37 
38 int main() {
39 
40  gRandom->SetSeed(14);
41 
42  // init MeasurementCreator
43  genfit::MeasurementCreator measurementCreator;
44 
45 
46  // init geometry and mag. field
47  new TGeoManager("Geometry", "Geane geometry");
48  TGeoManager::Import("genfitGeom.root");
49  genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., 15.)); // 15 kGauss
51 
52 
53  // init event display
55 
56 
57  // init fitter
59 
60 
61  TClonesArray myDetectorHitArray("genfit::mySpacepointDetectorHit");
62 
63  // init the factory
64  int myDetId(1);
67  factory.addProducer(myDetId, &myProducer);
68 
69 
70  // main loop
71  for (unsigned int iEvent=0; iEvent<100; ++iEvent){
72 
73  myDetectorHitArray.Clear();
74 
75  //TrackCand
76  genfit::TrackCand myCand;
77 
78  // true start values
79  TVector3 pos(0, 0, 0);
80  TVector3 mom(1.,0,0);
81  mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
82  mom.SetTheta(gRandom->Uniform(0.4*TMath::Pi(),0.6*TMath::Pi()));
83  mom.SetMag(gRandom->Uniform(0.2, 1.));
84 
85 
86  // helix track model
87  const int pdg = 13; // particle pdg code
88  const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
89  genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
90  measurementCreator.setTrackModel(helix);
91 
92 
93  unsigned int nMeasurements = gRandom->Uniform(5, 15);
94 
95  // covariance
96  double resolution = 0.01;
97  TMatrixDSym cov(3);
98  for (int i = 0; i < 3; ++i)
99  cov(i,i) = resolution*resolution;
100 
101  for (unsigned int i=0; i<nMeasurements; ++i) {
102  // "simulate" the detector
103  TVector3 currentPos = helix->getPos(i*2.);
104  currentPos.SetX(gRandom->Gaus(currentPos.X(), resolution));
105  currentPos.SetY(gRandom->Gaus(currentPos.Y(), resolution));
106  currentPos.SetZ(gRandom->Gaus(currentPos.Z(), resolution));
107 
108  // Fill the TClonesArray and the TrackCand
109  // In a real experiment, you detector code would deliver mySpacepointDetectorHits and fill the TClonesArray.
110  // The patternRecognition would create the TrackCand.
111  new(myDetectorHitArray[i]) genfit::mySpacepointDetectorHit(currentPos, cov);
112  myCand.addHit(myDetId, i);
113  }
114 
115 
116  // smeared start values (would come from the pattern recognition)
117  const bool smearPosMom = true; // init the Reps with smeared pos and mom
118  const double posSmear = 0.1; // cm
119  const double momSmear = 3. /180.*TMath::Pi(); // rad
120  const double momMagSmear = 0.1; // relative
121 
122  TVector3 posM(pos);
123  TVector3 momM(mom);
124  if (smearPosMom) {
125  posM.SetX(gRandom->Gaus(posM.X(),posSmear));
126  posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
127  posM.SetZ(gRandom->Gaus(posM.Z(),posSmear));
128 
129  momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
130  momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
131  momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
132  }
133 
134  // initial guess for cov
135  TMatrixDSym covSeed(6);
136  for (int i = 0; i < 3; ++i)
137  covSeed(i,i) = resolution*resolution;
138  for (int i = 3; i < 6; ++i)
139  covSeed(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
140 
141 
142  // set start values and pdg to cand
143  myCand.setPosMomSeedAndPdgCode(posM, momM, pdg);
144  myCand.setCovSeed(covSeed);
145 
146 
147  // create track
148  genfit::Track fitTrack(myCand, factory, new genfit::RKTrackRep(pdg));
149 
150 
151  // do the fit
152  try{
153  fitter->processTrack(&fitTrack);
154  }
155  catch(genfit::Exception& e){
156  std::cerr << e.what();
157  std::cerr << "Exception, next track" << std::endl;
158  continue;
159  }
160 
161  fitTrack.checkConsistency();
162 
163 
164  if (iEvent < 1000) {
165  // add track to event display
166  display->addEvent(&fitTrack);
167  }
168 
169 
170  }// end loop over events
171 
172  delete fitter;
173 
174  // open event display
175  display->open();
176 
177 }
178 
179