#ifndef L1Algo_h #define L1Algo_h 1 // #define TBB // TODO: Doesn't work now. Renew /// Debug features // #define PULLS // triplets pulls // #define TRIP_PERFORMANCE // triplets efficiencies // #define DOUB_PERFORMANCE // doublets efficiencies //#define DRAW // event display #define XXX // time debug //#define COUNTERS // diff counters (hits, doublets, ... ) //#define FIND_GAPED_TRACKS // use triplets with gaps #define EXTEND_TRACKS //#define MERGE_CLONES #include "L1Field.h" #include "L1Station.h" #include "L1StsHit.h" #include "L1Triplet.h" #include "L1Branch.h" #include "L1Track.h" #include "L1TrackPar.h" #include "L1Portion.h" #include "L1HitPoint.h" #include "L1Strip.h" #include "L1Grid.h" #include "L1HitsSortHelper.h" #include "../CbmL1Def.h" #include #include #include #include #include "omp.h" using std::vector; using std::map; #ifdef PULLS #define TRIP_PERFORMANCE class L1AlgoPulls; #endif #ifdef TRIP_PERFORMANCE template class L1AlgoEfficiencyPerformance; #endif #ifdef DOUB_PERFORMANCE template class L1AlgoEfficiencyPerformance; #endif typedef int Tindex; // struct hitCheck { // hitCheck(); // ~hitCheck(); // omp_lock_t Occupied; // int trackCandidateIndex; // fscal Chi2Track; // short int Length; // char ista; // bool UsedByTrack; // }; class L1Algo{ public: L1Algo(int nThreads=1): NStations(0), // number of all detector stations NMvdStations(0), // number of mvd stations vStsStrips(0), // strips positions created from hits. Front strips vStsStripsB(0), // back strips vStsZPos(0), // all possible z-positions of hits vStsHits(0), // hits as a combination of front-, backstrips and z-position vSFlag(0), // information of hits station & using hits in tracks(), vSFlagB(0), CATime(0), // time of trackfinding TripForHit(1000000), TripForHitV(1000000), stripB(1000000), stripF(1000000), NHitsIsecAll(0), NTracksIsecAll(0), vTracks(100000), // reconstructed tracks vRecoHits(500000),// packed hits of reconstructed tracks vStsDontUsedHits_A(1000000), vStsDontUsedHits_B(1000000), vStsDontUsedHits_Buf(1000000), vStsDontUsedHitsxy_A(1000000), vStsDontUsedHitsxy_buf(1000000), vStsDontUsedHitsxy_B(1000000), RealIHit_v(1000000), RealIHit_v_buf(1000000), RealIHit_v_buf2(1000000), sh (), TRACK_CHI2_CUT(10.), TRIPLET_CHI2_CUT(5.), DOUBLET_CHI2_CUT(5.), fNThreads(nThreads), Pick_gather(0.f), PickNeighbour(0.f), // (PickNeighbour < dp/dp_error) => triplets are neighbours MaxInvMom(0.f), // max considered q/p for tracks targX(0.f), targY(0.f), targZ(0.f), // target coor targB(), // field in the target point TargetXYInfo(), // target constraint [cm] vtxFieldRegion(),// really doesn't used vtxFieldValue(), // field at teh vertex position. //vTripletsP(), // container for triplets got in finding fTrackingLevel(0.f), fGhostSuppression(0.f), // really doesn't used fMomentumCutOff(0.f)// really doesn't used { for (int i=0; i<20; i++){ vTracks_local[i].resize(30000); vRecoHits_local[i].resize(2000000); } for (int i=0; i > HitCandidatesSortB; // vector > HitCandidatesSortF; // vector HitCandidatesBi; // vector HitCandidatesFi; L1HitsSortHelper sh; vector TripletsLocalReal[8][5000]; vector TripletsLocalRealAll[8][5000]; vector TripletsLocalRealV[8][5000]; vector TripletsLocalRealAllV[8][5000]; vector CandidatesTrack[20]; vector CandidatesTrack_buf[20]; vector CandidatesTrack_buf2[20]; vector* CandidatesTrack_bufP[20]; vector* CandidatesTrack_buf2P[20]; Tindex portionStopIndex[8]; vector n_g1; int numberCandidateThreadBuf [20]; void Init( const fscal geo[] ); #ifdef TBB2 void SetData( const vector< L1StsHit > & StsHits_, const vector< L1Strip > & StsStrips_, const vector< L1Strip > & StsStripsB_, const vector< fscal > & StsZPos_, const vector< unsigned char > & SFlag_, const vector< unsigned char > & SFlagB_, const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_ ); void PrintHits(); #endif /// The main procedure - find tracks. void CATrackFinder(); /// Track fitting procedures void KFTrackFitter_simple(); // version, which use procedured used during the reconstruction void L1KFTrackFitter(); // version from SIMD-KF benchmark /// ----- Input data ----- // filled in CbmL1::ReadEvent(); void SetNThreads(int n=1) {fNThreads = n;} enum{ MaxNStations = 10 }; int NStations, // number of all detector stations NMvdStations; // number of mvd stations L1Station vStations[MaxNStations] _fvecalignment; // station info const vector< L1Strip > *vStsStrips, // strips positions created from hits. Front strips *vStsStripsB; // back strips const vector< fscal > *vStsZPos; // all possible z-positions of hits const vector< L1StsHit > *vStsHits; // hits as a combination of front-, backstrips and z-position L1Grid vGrid[MaxNStations]; // hits as a combination of front-, backstrips and z-position const vector< unsigned char > *vSFlag, // information of hits station & using hits in tracks; *vSFlagB; const THitI *StsHitsStartIndex, *StsHitsStopIndex; // station-bounders in vStsHits array /// --- data used during finding iterations int isec; // iteration vector< L1StsHit > *vStsHitsUnused; L1StsHitV vStsHitsUnusedV[8]; vector< THitI > *RealIHitP; vector< THitI > *RealIHitPBuf; std::vector< L1HitPoint > *vStsHitPointsUnused; L1HitPointV vStsHitPointsUnusedV[8]; THitI *RealIHit; // index in vStsHits indexed by index in vStsHitsUnused THitI StsHitsUnusedStartIndex[MaxNStations+1], StsHitsUnusedStopIndex[MaxNStations+1]; THitI StsHitsUnusedStartIndexEnd[MaxNStations+1], StsHitsUnusedStopIndexEnd[MaxNStations+1]; vector< vector > TripForHit; vector< vector > TripForHitV; vector stripB; vector stripF; vector stripBUsed[20]; vector stripFUsed[20]; Tindex NHits_l[MaxNStations]; Tindex NHits_l_P[MaxNStations]; /// ----- Output data ----- double CATime; // time of trackfinding // L1Branch* pointer; int NHitsIsecAll; int NTracksIsecAll; vector< L1Track > vTracks; vector< THitI > vRecoHits; vector< L1StsHit > vStsDontUsedHits_A; vector< L1StsHit > vStsDontUsedHits_B; vector< L1StsHit > vStsDontUsedHits_Buf; vector< L1HitPoint > vStsDontUsedHitsxy_A; vector< L1HitPoint > vStsDontUsedHitsxy_buf; vector< L1HitPoint > vStsDontUsedHitsxy_B; vector< L1Track > vTracks_local[20]; vector< THitI > vRecoHits_local[20]; vector RealIHit_v; vector RealIHit_v_buf; vector RealIHit_v_buf2; int TripNPortion[8][5000]; int TripNPortionV[8][5000]; friend class CbmL1; const L1FieldValue& GetVtxFieldValue() const {return vtxFieldValue;} const L1FieldRegion& GetVtxFieldRegion() const {return vtxFieldRegion;} /// ----- Hit-point-strips conversion routines ------ void GetHitCoor(const L1StsHit& _h, fscal &_x, fscal &_y, fscal &_z, const L1Station &sta); void GetHitCoor(const L1StsHit& _h, fscal &_x, fscal &_y, char iS); void StripsToCoor(const fscal &u, const fscal &v, fscal &x, fscal &y, const L1Station &sta) const; // convert strip positions to coordinates void StripsToCoor(const fvec &u, const fvec &v, fvec &x, fvec &y, const L1Station &sta) const; void StripsToCoor(const fvec &u, const fvec &v, const fvec &z, fvec &x, fvec &y, char sta) const; L1HitPoint CreateHitPoint(const L1StsHit &hit, char ista); // full the hit point by hit information. L1HitPointV CreateHitPointV(const L1StsHitV &vStsHitsUnusedV, int ih, int ista, L1HitPointV* Point); private: /// ================================= FUNCTIONAL PART ================================= /// ----- Subroutines used by L1Algo::CATrackFinder() ------ void CAFindTrack(int ista, L1Branch& best_tr, unsigned char &best_L, fscal &best_chi2, const L1Triplet* curr_trip, L1Branch &curr_tr, unsigned char &curr_L, fscal &curr_chi2, unsigned char min_best_l, int &NCalls, L1Branch* new_tr ); /// Fit track /// t - track with hits /// T - track params /// dir - 0 - forward, 1 - backward /// qp0 - momentum for extrapolation /// initialize - should be params ititialized. 1 - yes. void BranchFitterFast(const L1Branch &t, L1TrackPar& T, const bool dir, const fvec qp0 = 0.f, const bool initParams = true); /// Fit track. more precise than FitterFast void BranchFitter(const L1Branch &t, L1TrackPar& T, const bool dir, const fvec qp0 = 0.f, const bool initParams = true); /// Find additional hits for existing track /// t - track with hits /// T - track params /// dir - 0 - forward, 1 - backward /// qp0 - momentum for extrapolation void FindMoreHits(L1Branch &t, L1TrackPar& T, const bool dir, const fvec qp0 = 0.0f); /// Find additional hits for existing track /// return chi2 fscal BranchExtender(L1Branch &t); /// ----- Subroutines used by L1Algo::CAMergeClones() ------ void InvertCholetsky(fvec a[15]); void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]); void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]); void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]); void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5], fvec W[15], fvec *chi2); void CAMergeClones(); /// -- Flags routines -- unsigned char GetFStation( unsigned char flag ){ return flag/4; } bool GetFUsed ( unsigned char flag ){ return (flag&0x02)!=0; } // bool GetFUsedD ( unsigned char flag ){ return (flag&0x01)!=0; } void SetFStation ( unsigned char &flag, unsigned int iStation ){ flag = iStation*4 + (flag%4); } void SetFUsed ( unsigned char &flag ){ flag |= 0x02; } // void SetFUsedD ( unsigned char &flag ){ flag |= 0x01; } void SetFUnUsed ( unsigned char &flag ){ flag &= 0xFC; } // void SetFUnUsedD ( unsigned char &flag ){ flag &= 0xFE; } /// Prepare the portion of left hits data void f10( // input Tindex start_lh, Tindex n1, L1HitPoint *vStsHits_l, // output fvec *u_front, fvec *u_back, fvec *zPos, THitI *hitsl_1 ); /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station. void f11( // input int istal, int istam, Tindex n1_V, fvec *u_front, fvec *u_back, fvec *zPos, // output L1TrackPar *T_1, L1FieldRegion *fld_1 ); /// Find the doublets. Reformat data in the portion of doublets. void f20( // input Tindex n1, L1Station &stam, L1HitPoint *vStsHits_m, L1TrackPar *T_1, // THitI *hitsl_1, // output Tindex &n2, vector &i1_2, #ifdef DOUB_PERFORMANCE vector &hitsl_2, #endif // DOUB_PERFORMANCE vector &hitsm_2 // , // vector &lmDuplets ); /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station. void f11SIMD( // input int istal, int istam, Tindex n1_V, fvec *u_front, fvec *u_back, fvec *zPos, // output L1TrackParV *T_1V, L1FieldRegionV *fld_1V ); /// Find the doublets. Reformat data in the portion of doublets. void f20SIMD( // input Tindex n1, L1Station &stam, L1TrackParV *T_1V, L1TrackParV *T_1VBuf, L1TrackParV *T_1VBuf2, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, // THitI *hitsl_1, // output Tindex &n2, int indexHit, nsL1::vector ::TSimd &LHits3 ); /// Add the middle hits to parameters estimation. Propagate to right station. /// Find the triplets (right hit). Reformat data in the portion of triplets. void f30( // input L1HitPoint *vStsHits_r, L1Station &stam, L1Station &star, int istam, int istar, L1HitPoint *vStsHits_m, L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, Tindex n2, vector &hitsm_2, vector &i1_2, // const vector &mrDuplets, // output Tindex &n3, nsL1::vector::TSimd &T_3, vector &hitsl_3, vector &hitsm_3, vector &hitsr_3, nsL1::vector::TSimd &u_front_3, nsL1::vector::TSimd &u_back_3, nsL1::vector::TSimd &z_Pos_3 ); /// Add the right hits to parameters estimation. void f31( // input Tindex n3_V, L1Station &star, nsL1::vector::TSimd &u_front_3, nsL1::vector::TSimd &u_back_3, nsL1::vector::TSimd &z_Pos_3, // output nsL1::vector::TSimd &T_3 ); /// Refit Triplets. void f32( // input Tindex n3, int istal, nsL1::vector::TSimd &T_3, vector &hitsl_3, vector &hitsm_3, vector &hitsr_3, int nIterations = 0 ); /// Select triplets. Save them into vTriplets. void f4( // input Tindex n3, int istal, int istam, int istar, nsL1::vector::TSimd &T_3, vector &hitsl_3, vector &hitsm_3, vector &hitsr_3, // output Tindex &nstaltriplets, Tindex ip_cur // #ifdef XXX // ,unsigned int &stat_n_trip // #endif ); /// Add the middle hits to parameters estimation. Propagate to right station. /// Find the triplets (right hit). Reformat data in the portion of triplets. void f30SIMD( // input L1Station &stam, L1Station &star, int istam, int istar, L1TrackParV *T_1V, L1TrackParV *T_1VBuf2, L1TrackParV *T_1VBuf, L1FieldRegionV *fld_1V, L1TrackParV &TripletsArray, Tindex n2, // const vector &mrDuplets, // output Tindex &n3, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, nsL1::vector ::TSimd &RHits, nsL1::vector::TSimd &u_front_3V, nsL1::vector::TSimd &u_back_3V, nsL1::vector::TSimd &z_Pos_3V, nsL1::vector ::TSimd &TripletsArrayV, nsL1::vector ::TSimd &LHitsV, nsL1::vector ::TSimd &MHitsV, nsL1::vector ::TSimd &LHits3 ); /// Add the right hits to parameters estimation. void f31SIMD( // input Tindex n3_V, L1Station &star, nsL1::vector::TSimd &u_front_3, nsL1::vector::TSimd &u_back_3, nsL1::vector::TSimd &z_Pos_3, // output L1TrackParV &TripletsArray, nsL1::vector ::TSimd &TripletsArrayV ); /// Refit Triplets. void f32SIMD( // input Tindex n3, int istal, L1TrackParV &TripletsArray, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, nsL1::vector ::TSimd &RHits, int nIterations = 0 ); /// Select triplets. Save them into vTriplets. void f4SIMD( // input Tindex n3, int istal, int istam, int istar, L1TrackParV &TripletsArray, // output Tindex &nstaltriplets, Tindex ip_cur, nsL1::vector ::TSimd &TripletsArrayV, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, nsL1::vector ::TSimd &RHits, nsL1::vector ::TSimd &LHitsV, nsL1::vector ::TSimd &MHitsV, nsL1::vector ::TSimd &LHits3, Tindex n3S // #ifdef XXX // ,unsigned int &stat_n_trip // #endif ); /// Find doublets on station void DupletsStaPort( // input int istal, int istam, Tindex ip, vector& n_g1, Tindex *portionStopIndex, // output L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, // vector &lmDuplets, Tindex &n_2, vector &i1_2, vector &hitsm_2 ); void DupletsStaPortSIMD( // input int istal, int istam, Tindex ip, vector& n_g1, Tindex *portionStopIndex, // output L1TrackParV *T_1V, L1TrackParV *T_1VBuf, L1TrackParV *T_1VBuf2, L1FieldRegionV *fld_1V, // vector &lmDuplets, Tindex &n_2, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, nsL1::vector ::TSimd &LHits3 ); /// Find triplets on station void TripletsStaPort( // input int istal, int istam, int istar, Tindex& nstaltriplets, L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, Tindex &n_2, vector &i1_2, vector &hitsm_2, // const vector &mrDuplets, // output Tindex ip_cur ); /// Find triplets on station void TripletsStaPortSIMD( // input int istal, int istam, int istar, Tindex& nstaltriplets, L1TrackParV *T_1V, L1TrackParV *T_1VBuf2, L1TrackParV *T_1VBuf, L1FieldRegionV *fld_1V, L1TrackParV &TripletsArray, Tindex &n_2, // const vector &mrDuplets, // output Tindex ip_cur, nsL1::vector ::TSimd &MHits, nsL1::vector ::TSimd &LHits, nsL1::vector ::TSimd &RHits, nsL1::vector ::TSimd &TripletsArrayV, nsL1::vector ::TSimd &LHits3 ); /// ------ Subroutines used by L1Algo::KFTrackFitter() ------ void GuessVec( L1TrackPar &t, fvec *xV, fvec *yV, fvec *zV, fvec *Sy, fvec *wV, int NHits, fvec *zCur = 0 ); void FilterFirst( L1TrackPar &track,fvec &x, fvec &y, L1Station &st ); #ifdef TBB enum { nthreads = 3, // number of threads nblocks = 1 // number of stations on one thread }; friend class ParalleledDup; friend class ParalleledTrip; #endif // TBB #ifdef TBB2 public: Tindex thrId; #endif // TBB2 private: /// ================================= DATA PART ================================= /// ----- Different parameters of CATrackFinder ----- Tindex FIRSTCASTATION; //first station used in CA // fNFindIterations - set number of interation for trackfinding // itetation of finding: #ifdef FIND_GAPED_TRACKS enum { fNFindIterations = 4 }; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter4 enum { kFastPrimIter, // primary fast tracks kAllPrimIter, // primary all tracks kAllPrimJumpIter, // primary tracks with jumped triplets kAllSecIter, // secondary all tracks kFastPrimJumpIter, // primary fast tracks with jumped triplets kFastPrimIter2, kAllSecJumpIter // secondary tracks with jumped triplets }; #else enum { fNFindIterations = 1}; enum { kFastPrimIter = 0, // primary fast tracks kAllPrimIter, // primary all tracks kAllSecIter, // secondary all tracks kFastPrimJumpIter, // disabled kAllPrimJumpIter, // disabled kFastPrimIter2, kAllSecJumpIter }; #endif // FIND_GAPED_TRACKS map threadNumberToCpuMap; int fNThreads; float TRACK_CHI2_CUT; float TRIPLET_CHI2_CUT; // = 5.0; // cut for selecting triplets before collecting tracks.per one DoF float DOUBLET_CHI2_CUT; fvec MaxDZ; // correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm) /// parameters which are different for different iterations. Set in the begin of CAL1TrackFinder float Pick_gather; // same for attaching additional hits to track float PickNeighbour; // (PickNeighbour < dp/dp_error) => triplets are neighbours fvec MaxInvMom; // max considered q/p for tracks fvec MaxSlope; // max slope (tx\ty) in prim vertex fvec targX, targY, targZ; // target coor L1FieldValue targB _fvecalignment; // field in the target point L1XYMeasurementInfo TargetXYInfo _fvecalignment; // target constraint [cm] /// standard sizes of the arrays enum { multiCoeff = 1, // central - 1, mbias coeff = 64/4, Portion = 1024/coeff, // portion of left hits MaxPortionDoublets = 10000/5 * 64/2 /coeff /*/ multiCoeff*/*1, MaxPortionTriplets = 10000*5 * 64/2 /coeff /*/ multiCoeff*/*1, MaxNPortion = 40 * coeff / multiCoeff, MaxArrSize = MaxNPortion*MaxPortionDoublets/MaxNStations //200000, // standart size of big arrays // mas be 40000 for normal work in cbmroot! }; L1FieldRegion vtxFieldRegion _fvecalignment;// really doesn't used L1FieldValue vtxFieldValue _fvecalignment; // field at teh vertex position. // vector vTriplets; // container for triplets got in finding // vector vTripletsP; int numPortions[8]; vector *TripletsLocal[MaxNStations-2]; // int TripNumThread; int fTrackingLevel, fGhostSuppression; // really doesn't used float fMomentumCutOff;// really doesn't used /// ----- Debug features ----- #ifdef PULLS L1AlgoPulls* fL1Pulls; #endif #ifdef TRIP_PERFORMANCE L1AlgoEfficiencyPerformance<3>* fL1Eff_triplets; L1AlgoEfficiencyPerformance<3>* fL1Eff_triplets2; #endif #ifdef DOUB_PERFORMANCE L1AlgoEfficiencyPerformance<2>* fL1Eff_doublets; #endif #ifdef DRAW friend class L1AlgoDraw; #endif } _fvecalignment; #endif