test/measurementFactoryExample/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 <TrackCand.h>
0008 #include <TrackPoint.h>
0009
0010 #include <MeasurementProducer.h>
0011 #include <MeasurementFactory.h>
0012
0013 #include "mySpacepointDetectorHit.h"
0014 #include "mySpacepointMeasurement.h"
0015
0016 #include <MaterialEffects.h>
0017 #include <RKTrackRep.h>
0018 #include <TGeoMaterialInterface.h>
0019
0020 #include <EventDisplay.h>
0021
0022 #include <HelixTrackModel.h>
0023 #include <MeasurementCreator.h>
0024
0025 #include <TDatabasePDG.h>
0026 #include <TEveManager.h>
0027 #include <TGeoManager.h>
0028 #include <TRandom.h>
0029 #include <TVector3.h>
0030 #include <vector>
0031
0032 #include "TDatabasePDG.h"
0033 #include <TMath.h>
0034
0035
0036
0037
0038 int main() {
0039
0040   gRandom->SetSeed(14);
0041
0042   // init MeasurementCreator
0043   genfit::MeasurementCreator measurementCreator;
0044
0045
0046   // init geometry and mag. field
0047   new TGeoManager("Geometry", "Geane geometry");
0048   TGeoManager::Import("genfitGeom.root");
0049   genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., 15.)); // 15 kGauss
0050   genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
0051
0052
0053   // init event display
0054   genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
0055
0056
0057   // init fitter
0058   genfit::AbsKalmanFitter* fitter = new genfit::KalmanFitterRefTrack();
0059
0060
0061   TClonesArray myDetectorHitArray("genfit::mySpacepointDetectorHit");
0062
0063   // init the factory
0064   int myDetId(1);
0065   genfit::MeasurementFactory<genfit::AbsMeasurement> factory;
0066   genfit::MeasurementProducer<genfit::mySpacepointDetectorHit, genfit::mySpacepointMeasurement> myProducer(&myDetectorHitArray);
0067   factory.addProducer(myDetId, &myProducer);
0068
0069
0070   // main loop
0071   for (unsigned int iEvent=0; iEvent<100; ++iEvent){
0072
0073     myDetectorHitArray.Clear();
0074
0075     //TrackCand
0076     genfit::TrackCand myCand;
0077
0078     // true start values
0079     TVector3 pos(0, 0, 0);
0080     TVector3 mom(1.,0,0);
0081     mom.SetPhi(gRandom->Uniform(0.,2*TMath::Pi()));
0082     mom.SetTheta(gRandom->Uniform(0.4*TMath::Pi(),0.6*TMath::Pi()));
0083     mom.SetMag(gRandom->Uniform(0.2, 1.));
0084
0085
0086     // helix track model
0087     const int pdg = 13;               // particle pdg code
0088     const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
0089     genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
0090     measurementCreator.setTrackModel(helix);
0091
0092
0093     unsigned int nMeasurements = gRandom->Uniform(5, 15);
0094
0095     // covariance
0096     double resolution = 0.01;
0097     TMatrixDSym cov(3);
0098     for (int i = 0; i < 3; ++i)
0099       cov(i,i) = resolution*resolution;
0100
0101     for (unsigned int i=0; i<nMeasurements; ++i) {
0102       // "simulate" the detector
0103       TVector3 currentPos = helix->getPos(i*2.);
0104       currentPos.SetX(gRandom->Gaus(currentPos.X(), resolution));
0105       currentPos.SetY(gRandom->Gaus(currentPos.Y(), resolution));
0106       currentPos.SetZ(gRandom->Gaus(currentPos.Z(), resolution));
0107
0108       // Fill the TClonesArray and the TrackCand
0109       // In a real experiment, you detector code would deliver mySpacepointDetectorHits and fill the TClonesArray.
0110       // The patternRecognition would create the TrackCand.
0111       new(myDetectorHitArray[i]) genfit::mySpacepointDetectorHit(currentPos, cov);
0112       myCand.addHit(myDetId, i);
0113     }
0114
0115
0116     // smeared start values (would come from the pattern recognition)
0117     const bool smearPosMom = true;   // init the Reps with smeared pos and mom
0118     const double posSmear = 0.1;     // cm
0119     const double momSmear = 3. /180.*TMath::Pi();     // rad
0120     const double momMagSmear = 0.1;   // relative
0121
0122     TVector3 posM(pos);
0123     TVector3 momM(mom);
0124     if (smearPosMom) {
0125       posM.SetX(gRandom->Gaus(posM.X(),posSmear));
0126       posM.SetY(gRandom->Gaus(posM.Y(),posSmear));
0127       posM.SetZ(gRandom->Gaus(posM.Z(),posSmear));
0128
0129       momM.SetPhi(gRandom->Gaus(mom.Phi(),momSmear));
0130       momM.SetTheta(gRandom->Gaus(mom.Theta(),momSmear));
0131       momM.SetMag(gRandom->Gaus(mom.Mag(), momMagSmear*mom.Mag()));
0132     }
0133
0134     // initial guess for cov
0135     TMatrixDSym covSeed(6);
0136     for (int i = 0; i < 3; ++i)
0137       covSeed(i,i) = resolution*resolution;
0138     for (int i = 3; i < 6; ++i)
0139       covSeed(i,i) = pow(resolution / nMeasurements / sqrt(3), 2);
0140
0141
0142     // set start values and pdg to cand
0143     myCand.setPosMomSeedAndPdgCode(posM, momM, pdg);
0144     myCand.setCovSeed(covSeed);
0145
0146
0147     // create track
0148     genfit::Track fitTrack(myCand, factory, new genfit::RKTrackRep(pdg));
0149
0150
0151     // do the fit
0152     try{
0153       fitter->processTrack(&fitTrack);
0154     }
0155     catch(genfit::Exception& e){
0156       std::cerr << e.what();
0157       std::cerr << "Exception, next track" << std::endl;
0158       continue;
0159     }
0160
0161     fitTrack.checkConsistency();
0162
0163
0164     if (iEvent < 1000) {
0165       // add track to event display
0166       display->addEvent(&fitTrack);
0167     }
0168
0169
0170   }// end loop over events
0171
0172   delete fitter;
0173
0174   // open event display
0175   display->open();
0176
0177 }
0178
0179