1 #ifndef FPCCDSILICONTRACKING_MarlinTrk_H
2 #define FPCCDSILICONTRACKING_MarlinTrk_H 1
4 #include "marlin/Processor.h"
5 #include <marlin/Global.h>
12 #include <IMPL/TrackImpl.h>
13 #include "ClusterExtended.h"
14 #include "TrackExtended.h"
15 #include "TrackerHitExtended.h"
16 #include "HelixClass_double.h"
17 #include <EVENT/SimTrackerHit.h>
18 #include <EVENT/MCParticle.h>
19 #include <IMPL/TrackImpl.h>
21 #include "MarlinTrk/IMarlinTrack.h"
23 #include <UTIL/BitField64.h>
24 #include <UTIL/ILDConf.h>
30 #include "TStopwatch.h"
32 #include "../src/GetPurity.h"
34 namespace DiagnosticsHistograms {
38 namespace DiagnosticsHistograms2D {
42 using namespace lcio ;
43 using namespace marlin ;
51 class IMarlinTrkSystem ;
55 class LCRelationNavigator ;
229 virtual void init() ;
233 virtual void processRunHeader( LCRunHeader* run ) ;
237 virtual void processEvent( LCEvent * evt ) ;
240 virtual void check( LCEvent * evt ) ;
252 EVENT::LCEvent* _current_event;
254 int _nDivisionsInPhi;
255 int _nDivisionsInTheta;
258 MarlinTrk::HelixFit* _fastfitter;
263 bool _runMarlinTrkDiagnostics;
264 std::string _MarlinTrkDiagnosticsName;
266 bool _MSOn, _ElossOn, _SmoothOn ;
268 float _initialTrackError_d0;
269 float _initialTrackError_phi0;
270 float _initialTrackError_omega;
271 float _initialTrackError_z0;
272 float _initialTrackError_tanL;
274 double _maxChi2PerHit;
275 double _maxChi2PerHit2nd;
277 bool _UseEventDisplay;
278 std::vector<int> _colours;
285 bool _createDiagnosticsHistograms;
289 int _ntriplets, _ntriplets_good, _ntriplets_2MCP, _ntriplets_3MCP, _ntriplets_1MCP_Bad, _ntriplets_bad;
293 LCCollection* GetCollection( LCEvent * evt, std::string colName ) ;
296 LCRelationNavigator* GetRelations( LCEvent * evt, std::string RelName ) ;
305 bool operator()(TrackExtended *a, TrackExtended *b)
const {
306 if ( a == b )
return false;
307 return (a->getChi2()/a->getNDF() < b->getChi2()/b->getNDF() );
312 std::string _VTXHitCollection;
313 std::string _FTDPixelHitCollection;
314 std::string _FTDSpacePointCollection;
315 std::string _SITHitCollection;
316 std::string _siTrkCollection;
318 std::vector< LCCollection* > _colTrackerHits;
319 std::map< LCCollection*, std::string > _colNamesTrackerHits;
321 std::vector<TrackerHitExtendedVec> _sectors;
322 std::vector<TrackerHitExtendedVec> _sectorsFTD;
337 _tracksNHits.resize(maxHits-2);
338 _maxIndex=(maxHits-3);
351 return _tracksNHits.at(std::min(nHits-3, _maxIndex));
355 std::vector< TrackExtendedVec > _tracksNHits;
359 TracksWithNHitsContainer _tracksWithNHitsContainer;
361 int InitialiseVTX(LCEvent * evt);
362 int InitialiseFTD(LCEvent * evt);
363 void ProcessOneSector(
int iSectorPhi,
int iSectorTheta);
364 void ProcessOneSectorVer2(
int iSectorPhi,
int iSectorTheta);
366 TrackExtended * TestTriplet(TrackerHitExtended * outerHit,
367 TrackerHitExtended * middleHit,
368 TrackerHitExtended * innerHit,
369 HelixClass_double & helix,
372 int BuildTrack_KalFit(TrackerHitExtended * outerHit,
373 TrackerHitExtended * middleHit,
374 TrackerHitExtended * innerHit,
375 HelixClass_double & helix,
377 TrackExtended * trackAR);
379 int _useBuildTrackForHighPt;
380 double _cosThetaRangeForBuildTrackForHighPt;
381 double _phiRangeForBuildTrackForHighPt;
382 void getPhiThetaRegionForHighPt(
int* boundaries,TrackExtended* trackAR);
384 void Sorting( TrackExtendedVec & trackVec);
385 void CreateTrack(TrackExtended * trackAR );
386 void AttachRemainingVTXHitsSlow();
387 void AttachRemainingFTDHitsSlow();
388 void AttachRemainingVTXHitsFast();
389 void AttachRemainingFTDHitsFast();
390 void AttachRemainingVTXHitsVeryFast();
391 void TrackingInFTD();
392 int BuildTrackFTD(TrackExtended* trackAR,
int* nLR,
int iS);
393 int AttachHitToTrack(TrackExtended * trackAR, TrackerHitExtended * hit,
int iopt);
394 int AttachHitToTrack_KalFit(TrackExtended * trackAR, TrackerHitExtended * hit);
396 void FinalRefit(LCCollectionVec* trk_col, LCCollectionVec* rel_col);
399 float _chi2WRPhiTriplet;
400 float _chi2WRPhiQuartet;
401 float _chi2WRPhiSeptet;
402 float _chi2WZTriplet;
403 float _chi2WZQuartet;
405 float _minDistCutAttachForFTD;
406 float _minDistCutAttachForVXD;
407 int _minimalLayerToAttach;
408 int _minMissAddition;
411 static const double TWOPI;
419 std::vector<int> _Combinations;
420 std::vector<int> _CombinationsFTD;
422 float _resolutionRPhiVTX;
423 float _resolutionZVTX;
425 float _resolutionRPhiFTD;
426 float _resolutionZFTD;
428 float _resolutionRPhiSIT;
429 float _resolutionZSIT;
431 float _phiCutForMerging;
432 float _tanlambdaCutForMerging;
433 float _angleCutForMerging;
439 float _minDistToDelta;
444 float _chi2FitCut_lowPt;
447 TrackExtendedVec _trackImplVec;
450 float _cutOnD0, _cutOnZ0, _cutOnOmegaVXD, _cutOnOmegaFTD;
451 double _cutOnPtVXD,_cutOnPtFTD;
458 int _max_hits_per_sector;
460 int _nTotalVTXHits,_nTotalFTDHits,_nTotalSITHits;
467 UTIL::BitField64* _encoder;
468 int getDetectorID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::subdet]; }
469 int getSideID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::side]; };
470 int getLayerID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::layer]; };
471 int getModuleID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::module]; };
472 int getSensorID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::sensor]; };
474 int getDetectorID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::subdet]; }
475 int getSideID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::side]; };
476 int getLayerID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::layer]; };
477 int getModuleID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::module]; };
478 int getSensorID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0());
return (*_encoder)[lcio::ILDCellID0::sensor]; };
480 void setupGearGeom(
const gear::GearMgr* gearMgr ) ;
483 unsigned int _nLayersVTX;
485 unsigned int _nLayersSIT;
488 std::vector<float> _zLayerFTD;
490 unsigned int _nlayersFTD;
491 bool _petalBasedFTDWithOverlaps;
494 int _output_track_col_quality;
495 static const int _output_track_col_quality_GOOD;
496 static const int _output_track_col_quality_FAIR;
497 static const int _output_track_col_quality_POOR;
500 float _chi2FitCut_kalman;
501 bool _useClusterRejection;
502 float _minDotOf2Clusters;
504 int getIntersectionEasy(HelixClass_double& helix, TrackerHit* curInmos ,
int layer,
double* isec,
double* ref);
505 int getIntersectionEasyTest(HelixClass_double& helix, TrackerHit* basisTrkhit ,
int layer, std::vector<double> &vec );
506 int CheckTiltOf2Clusters(TrackerHit* A, TrackerHit* B,
int level);
507 float DotOf2Clusters(TrackerHit* A, TrackerHit* B);
508 int KalFit(
int& ndf,
float& Chi2, TrackerHitVec trkHits,TrackerHitVec& hits_in_fit, TrackerHitVec& outliers,
float* par ,
float* epar, HelixClass_double& helix);
509 int getPhiThetaRegion(TrackExtended* trackAR,
int layer,
int* Boundaries);
516 std::vector<double> cosphi;
517 std::vector<double> sinphi;
518 std::vector<double> phi;
519 std::vector<double> phiAtXiMin;
520 std::vector<double> phiAtXiMax;
521 std::vector<double> ladder_incline;
534 std::vector<GeoData_t> geodata;
537 void InitVXDGeometry();
538 void InitSITGeometry();
539 FloatVec _pixelSizeVec;
541 TVector3 LocalToGlobal(TVector3 local,
int layer,
int ladder);
543 void calcTrackParameterOfMCP(MCParticle* pmcp,
double* par);
551 unsigned int xiwidth ;
552 unsigned int zetawidth ;
555 unsigned int quality ;
558 cellid0 = hit->getCellID0();
559 cellid1 = hit->getCellID1();
562 xiwidth = cellid1 & 0x000001ff ;
563 zetawidth = ( cellid1 >> 9 ) & 0x000001ff;
564 nPix = ( cellid1 >> 18) & 0x000003ff;
565 tilt = ( cellid1 >> 28) & 0x00000003;
566 quality = ( cellid1 >> 30) & 0x00000003;
567 const int cellId_trk = hit->getCellID0();
568 UTIL::BitField64 encoder_trk( lcio::ILDCellID0::encoder_string ) ;
569 encoder_trk.setValue(cellId_trk) ;
570 layer = encoder_trk[lcio::ILDCellID0::layer] ;
571 ladder = encoder_trk[lcio::ILDCellID0::module] ;
577 typedef std::map< std::pair< int, int >,
int > RangeMap;
578 RangeMap _phiRangeForTriplet;
580 double getNeededPhiSectors(
double Pt,
int outly ,
int inly);
584 typedef std::map< std::vector<int> , std::vector<int> > RangeMapVer2;
585 RangeMapVer2 _phiRangeForTripletVer2;
586 void getNeededPhiSectorsVer2(
double Pt, std::vector<int> layers, std::vector<double>& phiDiff);
589 float _safetyPhiRange_ratio;
595 float _fudgeFactorForSITsr_rphi;
596 float _fudgeFactorForSITsr_z;
598 int _fudgeThetaRange;
600 std::string _colNameVXDTrackerHitRelations;
601 std::string _colNameSITSpacePointRelations;
602 std::string _colNameFTDSpacePointRelations;
603 std::string _colNameFTDPixelTrackerHitRelations;
604 LCRelationNavigator* _navVXD;
605 LCRelationNavigator* _navSIT;
606 LCRelationNavigator* _navFTDsp;
607 LCRelationNavigator* _navFTDpix;
608 std::string _colNameVXDSimHit;
609 std::string _colNameSITSimHit;
610 std::string _colNameFTDspSimHit;
611 std::string _colNameFTDpixSimHit;
612 LCCollection* _simVXD;
613 LCCollection* _simSIT;
614 LCCollection* _simFTDsp;
615 LCCollection* _simFTDpix;
619 double _nSigmaBuild_phi;
620 double _nSigmaBuild_theta;
632 bool _mydebugIntersection;
634 bool _mydebugRootVXD;
635 bool _mydebugRootFTD;
636 bool _mydebugRootVXDFTD;
637 bool _mydebugVXDHits;
638 bool _mydebugSITHits;
639 bool _mydebugTriplet;
640 int _mydebugTripletMode;
641 bool _mydebugTripletVXD;
642 bool _mydebugTripletFTD;
643 bool _mydebugTripletVXDFTD;
644 bool _mydebugBuildTrack;
645 int _mydebugBuildTrackMode;
646 bool _mydebugAttachVXD;
647 bool _mydebugPrintMCP;
652 TStopwatch InitialiseVTX;
653 TStopwatch InitialiseFTD;
654 TStopwatch ProcessOneSector;
655 TStopwatch TrackingInFTD;
657 TStopwatch CreateTrack;
658 TStopwatch AttachRemainingVXD;
659 TStopwatch AttachRemainingFTD;
660 TStopwatch FinalRefit;
663 InitialiseVTX.Reset();
664 InitialiseFTD.Reset();
665 ProcessOneSector.Reset();
666 TrackingInFTD.Reset();
669 AttachRemainingVXD.Reset();
670 AttachRemainingFTD.Reset();
674 printf(
"inAnEvt : RT=%.3f s, CPU=%.3f s \n",inAnEvt.RealTime(),inAnEvt.CpuTime());
675 printf(
"InitialiseVTX : RT=%.3f s, CPU=%.3f s \n",InitialiseVTX.RealTime(),InitialiseVTX.CpuTime());
676 printf(
"InitialiseFTD : RT=%.3f s, CPU=%.3f s \n",InitialiseFTD.RealTime(),InitialiseFTD.CpuTime());
677 printf(
"ProcessOneSector : RT=%.3f s, CPU=%.3f s \n",ProcessOneSector.RealTime(),ProcessOneSector.CpuTime());
678 printf(
"TrackingInFTD : RT=%.3f s, CPU=%.3f s \n",TrackingInFTD.RealTime(),TrackingInFTD.CpuTime());
679 printf(
"Sorting : RT=%.3f s, CPU=%.3f s \n",Sorting.RealTime(),Sorting.CpuTime());
680 printf(
"CreateTrack : RT=%.3f s, CPU=%.3f s \n",CreateTrack.RealTime(),CreateTrack.CpuTime());
681 printf(
"AttachRemainingVXD: RT=%.3f s, CPU=%.3f s \n",AttachRemainingVXD.RealTime(),AttachRemainingVXD.CpuTime());
682 printf(
"AttachRemainingFTD: RT=%.3f s, CPU=%.3f s \n",AttachRemainingFTD.RealTime(),AttachRemainingFTD.CpuTime());
683 printf(
"FinalRefit : RT=%.3f s, CPU=%.3f s \n",FinalRefit.RealTime(),FinalRefit.CpuTime());
687 TStopwatch _timer2Triplet;
688 TStopwatch _timer2Build;
689 bool _mydebugstopwatch2;
691 float _currentPurity;
693 std::map< MCParticle*, SimTrackerHitVec > _mcpVXD;
694 std::map< MCParticle*, SimTrackerHitVec > _mcpVXDFTD;
695 std::map< MCParticle*, SimTrackerHitVec > _mcpVXDSIT;
696 std::map< MCParticle*, SimTrackerHitVec > _mcpVXDFTDSIT;
697 std::map< MCParticle*, SimTrackerHitVec > _mcpFTD;
698 std::map< MCParticle*, SimTrackerHitVec > _mcpFTDSIT;
699 std::map< MCParticle*, SimTrackerHitVec > _mcpSIT;
702 enum MCPContributions{
714 enum MCPContributions getMCPContribution(IntVec nsub);
718 static int numOfProcesses;
725 SimTrackerHitVec simVecFromTrk;
726 std::pair< MCParticle* , SimTrackerHitVec > mcTruth ;
729 TrackerHitVec truthTrkHits ;
730 enum MCPContributions mcpcont;
742 Triplet(
bool incrementID =
true){
748 detID[0] = detID[1] = detID[2] = -1;
749 layer[0] = layer[1] = layer[2] = -1;
752 simVecFromTrk.clear();
753 mcTruth.first = 0; mcTruth.second.clear();
754 truthTrkHits.clear();
757 mcp_d0 = mcp_z0 = mcp_omega = mcp_phi0 = mcp_tanL = 1e20;
766 enum BuildTrackResult{
775 static int numOfProcesses;
778 BuildTrackResult result;
781 const MCParticle* truthmcp;
782 SimTrackerHitVec simvec;
785 SimTrackerHitVec simVecFromTrk;
789 result = BuildTrackResultSize;
795 simVecFromTrk.clear();
801 bool _availableInBuildedTrack;
803 typedef std::vector< BuildedTrack > BuildedTrackVec;
804 BuildedTrackVec _buildedTrackContainer;
805 void SetBuildedTrack(TrackExtended* trackAR, BuildedTrackVec& btrackvec);
806 void BuildedTrackDebuger1(BuildedTrackVec::iterator begin, BuildedTrackVec::iterator end);
807 void BuildedTrackDebuger2(std::vector<const BuildedTrack*>::iterator begin, std::vector<const BuildedTrack*>::iterator end);
812 std::pair<MCParticle*,SimTrackerHitVec> pair;
814 std::vector<const BuildedTrack*> BuildedTrackVec;
815 std::vector<const BuildedTrack*> BuildedTrack2ndVec;
816 std::vector<const BuildedTrack*> BuildedTrack3rdVec;
818 typedef std::map<MCParticle*,MCP_BuildedTrack> BTMap;
819 BTMap _mcpBuildedTrackContainer;
821 std::vector< std::map<MCParticle*, MCP_BuildedTrack> > _mcpRemainingBTContainerA;
823 std::vector< std::map<MCParticle*, MCP_BuildedTrack> > _mcpRemainingBTContainerB;
829 std::vector< Triplet > _tripletContainer;
833 std::pair<MCParticle*,SimTrackerHitVec> pair;
835 std::vector<const Triplet*> tripvec;
836 std::vector<const Triplet*> trip2ndvec;
837 std::vector<const Triplet*> trip3rdvec;
839 typedef std::map<MCParticle*,MCP_Triplet> TripMap;
840 TripMap _mcpTripletContainer;
843 std::vector< std::map<MCParticle*, MCP_Triplet> > _mcpRemainingTRContainerA;
845 std::vector< std::map<MCParticle*, MCP_Triplet> > _mcpRemainingTRContainerB;
848 void TripletDebugerWithMCPRemain(
int nth, std::vector< std::map<MCParticle*, MCP_Triplet> > A);
851 std::vector<LCRelationNavigator*> _naviVec;
852 IntVec getNHitsInSubDet(SimTrackerHitVec simvec);
853 void TripletDebuger1(std::vector<Triplet>::iterator begin, std::vector<Triplet>::iterator end);
854 void TripletDebuger2(std::vector<const Triplet*>::iterator begin, std::vector<const Triplet*>::iterator end);
859 int FPCCDSiliconTracking_MarlinTrk::Triplet::numOfProcesses = 0;
860 int FPCCDSiliconTracking_MarlinTrk::BuildedTrack::numOfProcesses = 0;
Compare tracks according to their chi2/ndf.
Definition: FPCCDSiliconTracking_MarlinTrk.h:303
void resize(size_t maxHits)
Set the size to allow a maximum of maxHit hits.
Definition: FPCCDSiliconTracking_MarlinTrk.h:336
TrackExtendedVec & getTracksWithNHitsVec(size_t nHits)
Returns the TrackExtendedVec for track with n hits.
Definition: FPCCDSiliconTracking_MarlinTrk.h:348
Definition: FPCCDSiliconTracking_MarlinTrk.h:649
Definition: FPCCDSiliconTracking_MarlinTrk.h:531
Definition: FPCCDSiliconTracking_MarlinTrk.h:773
Definition: FPCCDSiliconTracking_MarlinTrk.h:810
Definition: FPCCDSiliconTracking_MarlinTrk.h:545
MarlinTrk::IMarlinTrkSystem * _trksystem
pointer to the IMarlinTrkSystem instance
Definition: FPCCDSiliconTracking_MarlinTrk.h:262
A helper class to allow good code readability by accessing tracks with N hits.
Definition: FPCCDSiliconTracking_MarlinTrk.h:330
Definition: GetPurity.h:37
=== FPCCDSiliconTracking_MarlinTrk Processor === This processor is based on SiliconTracking_MarlinT...
Definition: FPCCDSiliconTracking_MarlinTrk.h:217
std::string _colNameMCParticles
input MCParticle collection and threshold used for Drawing
Definition: FPCCDSiliconTracking_MarlinTrk.h:299
Definition: FPCCDSiliconTracking_MarlinTrk.h:511
Definition: SiliconTracking_MarlinTrk.cc:90
Definition: FPCCDSiliconTracking_MarlinTrk.h:831
Definition: FPCCDSiliconTracking_MarlinTrk.h:716
int _safetyPhiRange_fix
Extra range in addition to main range used for triplet construction process and needed to find triple...
Definition: FPCCDSiliconTracking_MarlinTrk.h:594