#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 #include #include #include 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 long Tindex; class L1Algo{ public: L1Algo(): NStations(0), // number of all detector stations NMvdStations(0), // number of mvd stations vStsStrips(), // strips positions created from hits. Front strips vStsStripsB(), // back strips vStsZPos(), // all possible z-positions of hits vStsHits(), // hits as a combination of front-, backstrips and z-position vSFlag(), // information of hits station & using hits in tracks(), vSFlagB(), vTracks(), // reconstructed tracks vRecoHits(), // packed hits of reconstructed tracks CATime(0), // time of trackfinding TRACK_CHI2_CUT(10.), TRIPLET_CHI2_CUT(5.), DOUBLET_CHI2_CUT(5.), Pick_gather(0), PickNeighbour(0), // (PickNeighbour < dp/dp_error) => triplets are neighbours MaxInvMom(0), // max considered q/p for tracks targX(0), targY(0), targZ(0), // target coor targB(), // field in the target point TargetXYInfo(), // target constraint [cm] vtxFieldRegion(),// really doesn't used vtxFieldValue(), // field at teh vertex position. vTriplets(), // container for triplets got in finding fTrackingLevel(0), fGhostSuppression(0), // really doesn't used fMomentumCutOff(0)// really doesn't used { vTriplets.reserve(MaxArrSize); for ( int i = 0; i < MaxNStations; i++ ) { fLmDuplets[i].reserve(MaxNHitsPerStation); fLmDupletsG[i].reserve(MaxNHitsPerStation); } fHitsm_2.reserve(MaxPortionDoublets); fI1_2.reserve(MaxPortionDoublets); fHitsmG_2.reserve(MaxPortionDoublets); fI1G_2.reserve(MaxPortionDoublets); fT_3.reserve(MaxPortionTriplets/fvecLen); fHitsl_3.reserve(MaxPortionTriplets); fHitsm_3.reserve(MaxPortionTriplets); fHitsr_3.reserve(MaxPortionTriplets); fU_front3.reserve(MaxPortionTriplets/fvecLen); fU_back3.reserve(MaxPortionTriplets/fvecLen); fZ_pos3.reserve(MaxPortionTriplets/fvecLen); fNeighCands.reserve(8); // ~average is 1-2 for central, up to 5 fTrackcandidate.reserve(MaxNHitsPerStation); } void Init( const fscal geo[] ); 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(); /// 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(); enum{ MaxNStations = 10 }; int NStations, // number of all detector stations NMvdStations; // number of mvd stations L1Station vStations[MaxNStations] _fvecalignment; // station info vector< L1Strip > vStsStrips, // strips positions created from hits. Front strips vStsStripsB; // back strips vector< fscal > vStsZPos; // all possible z-positions of hits 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 vector< unsigned char > vSFlag, // information of hits station & using hits in tracks; vSFlagB; THitI StsHitsStartIndex[MaxNStations+1], StsHitsStopIndex[MaxNStations+1]; // station-bounders in vStsHits array /// --- data used during finding iterations int isec; // iteration vector< L1StsHit > *vStsHitsUnused; std::vector< L1HitPoint > *vStsHitPointsUnused; THitI *RealIHit; // index in vStsHits indexed by index in vStsHitsUnused THitI StsHitsUnusedStartIndex[MaxNStations+1], StsHitsUnusedStopIndex[MaxNStations+1]; /// ----- Output data ----- vector< L1Track > vTracks; // reconstructed tracks vector< THitI > vRecoHits; // packed hits of reconstructed tracks double CATime; // time of trackfinding 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 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; L1HitPoint CreateHitPoint(const L1StsHit &hit, char ista); // full the hit point by hit information. 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 ); void CleanTriplets(bool isClean = 0 #ifdef COUNTERS ,Tindex* nlevels = 0 #endif ); // isClean = true - find neighbours only /// 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., const bool initParams = true); /// Fit track. more precise than FitterFast void BranchFitter(const L1Branch &t, L1TrackPar& T, const bool dir, const fvec qp0 = 0., 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.0); /// 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 ); /// 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, vector &vTriplets_part, Tindex *TripStartIndexH, Tindex *TripStopIndexH // #ifdef XXX // ,unsigned int &stat_n_trip // #endif ); /// Find neighbours of triplets. Calculate level of triplets. void FindNeighbours( const vector& TripStartIndexH, const vector& TripStopIndexH #ifdef COUNTERS ,Tindex* nlevels = 0 #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 ); /// 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 vector *vTriplets_part, Tindex *TripStartIndexH, Tindex *TripStopIndexH ); /// ------ 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 & kAllSecJumpIter 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 = 3 }; enum { kFastPrimIter = 0, // primary fast tracks kAllPrimIter, // primary all tracks kAllSecIter, // secondary all tracks kFastPrimJumpIter, // disabled kAllPrimJumpIter, // disabled kFastPrimIter2, kAllSecJumpIter }; #endif // FIND_GAPED_TRACKS 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! MaxNHitsPerStation = 3000 }; L1FieldRegion vtxFieldRegion _fvecalignment;// really doesn't used L1FieldValue vtxFieldValue _fvecalignment; // field at teh vertex position. vector vTriplets; // container for triplets got in finding Tindex TripStartIndex[MaxNStations-2], TripStopIndex[MaxNStations-2]; // containers for stations bounders in vTriplets. indices are from start to stop-1 int fTrackingLevel, fGhostSuppression; // really doesn't used float fMomentumCutOff;// really doesn't used /// ---- containers for temporary data --- vector fLmDuplets[MaxNStations]; // is exist a doublet started from indexed by left hit vector fLmDupletsG[MaxNStations]; // is exist a doublet started from indexed by left hit /// has to be cleaned after each use vector fHitsm_2; /// middle hits indexed by number of doublets in portion(i2) vector fI1_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2) vector fHitsmG_2; vector fI1G_2; nsL1::vector::TSimd fT_3; vector fHitsl_3, fHitsm_3, fHitsr_3; nsL1::vector::TSimd fU_front3, fU_back3, fZ_pos3; vector fNeighCands; // to save neighbour candidates vector fTrackcandidate; /// ----- 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