test/minimalFittingExample/main.cc

0001 #include <ConstField.h>
0002 #include <Exception.h>
0003 #include <FieldManager.h>
0004 #include <KalmanFitterRefTrack.h>
0005 #include <StateOnPlane.h>
0006 #include <Track.h>
0007 #include <TrackPoint.h>
0008
0009 #include <MaterialEffects.h>
0010 #include <RKTrackRep.h>
0011 #include <TGeoMaterialInterface.h>
0012
0013 #include <EventDisplay.h>
0014
0015 #include <HelixTrackModel.h>
0016 #include <MeasurementCreator.h>
0017
0018 #include <TDatabasePDG.h>
0019 #include <TEveManager.h>
0020 #include <TGeoManager.h>
0021 #include <TRandom.h>
0022 #include <TVector3.h>
0023 #include <vector>
0024
0025 #include "TDatabasePDG.h"
0026 #include <TMath.h>
0027
0028
0029
0030
0031 int main() {
0032
0033   gRandom->SetSeed(14);
0034
0035   // init MeasurementCreator
0036   genfit::MeasurementCreator measurementCreator;
0037
0038
0039   // init geometry and mag. field
0040   new TGeoManager("Geometry", "Geane geometry");
0041   TGeoManager::Import("genfitGeom.root");
0042   genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., 15.)); // 15 kGauss
0043   genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
0044
0045
0046   // init event display
0047   genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
0048
0049
0050   // init fitter
0051   genfit::AbsKalmanFitter* fitter = new genfit::KalmanFitterRefTrack();
0052
0053   // main loop
0054   for (unsigned int iEvent=0; iEvent<100; ++iEvent){
0055
0056     // true start values
0057     TVector3 pos(0, 0, 0);
0058     TVector3 mom(1.,0,0);
0059     mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
0060     mom.SetTheta(gRandom->Uniform(0.4*TMath::Pi(),0.6*TMath::Pi()));
0061     mom.SetMag(gRandom->Uniform(0.2, 1.));
0062
0063
0064     // helix track model
0065     const int pdg = 13;               // particle pdg code
0066     const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
0067     genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
0068     measurementCreator.setTrackModel(helix);
0069
0070
0071     unsigned int nMeasurements = gRandom->Uniform(5, 15);
0072
0073
0074     // smeared start values
0075     const bool smearPosMom = true;     // init the Reps with smeared pos and mom
0076     const double posSmear = 0.1;     // cm
0077     const double momSmear = 3. /180.*TMath::Pi();     // rad
0078     const double momMagSmear = 0.1;   // relative
0079
0080     TVector3 posM(pos);
0081     TVector3 momM(mom);
0082     if (smearPosMom) {
0083       posM.SetX(gRandom->Gaus(posM.X(),posSmear));
0084       posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
0085       posM.SetZ(gRandom->Gaus(posM.Z(),posSmear));
0086
0087       momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
0088       momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
0089       momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
0090     }
0091     // approximate covariance
0092     TMatrixDSym covM(6);
0093     double resolution = 0.01;
0094     for (int i = 0; i < 3; ++i)
0095       covM(i,i) = resolution*resolution;
0096     for (int i = 3; i < 6; ++i)
0097       covM(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
0098
0099
0100     // trackrep
0101     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
0102
0103     // smeared start state
0104     genfit::MeasuredStateOnPlane stateSmeared(rep);
0105     stateSmeared.setPosMomCov(posM, momM, covM);
0106
0107
0108     // create track
0109     TVectorD seedState(6);
0110     TMatrixDSym seedCov(6);
0111     stateSmeared.get6DStateCov(seedState, seedCov);
0112     genfit::Track fitTrack(rep, seedState, seedCov);
0113
0114
0115     // create random measurement types
0116     std::vector<genfit::eMeasurementType> measurementTypes;
0117     for (unsigned int i = 0; i < nMeasurements; ++i)
0118       measurementTypes.push_back(genfit::eMeasurementType(gRandom->Uniform(8)));
0119
0120
0121     // create smeared measurements and add to track
0122     try{
0123       for (unsigned int i=0; i<measurementTypes.size(); ++i){
0124         std::vector<genfit::AbsMeasurement*> measurements = measurementCreator.create(measurementTypes[i], i*5.);
0125         fitTrack.insertPoint(new genfit::TrackPoint(measurements, &fitTrack));
0126       }
0127     }
0128     catch(genfit::Exception& e){
0129       std::cerr<<"Exception, next track"<<std::endl;
0130       std::cerr << e.what();
0131       continue;
0132     }
0133
0134     //check
0135     fitTrack.checkConsistency();
0136
0137     // do the fit
0138     fitter->processTrack(&fitTrack);
0139
0140     //check
0141     fitTrack.checkConsistency();
0142
0143
0144     if (iEvent < 1000) {
0145       // add track to event display
0146       display->addEvent(&fitTrack);
0147     }
0148
0149
0150
0151   }// end loop over events
0152
0153   delete fitter;
0154
0155   // open event display
0156   display->open();
0157
0158 }
0159
0160