/* *===================================================== * * CBM Level 1 Reconstruction * * Authors: I.Kisel, S.Gorbunov * * e-mail : ikisel@kip.uni-heidelberg.de * *===================================================== * * Finds tracks using the Cellular Automaton algorithm * */ #include "L1Algo/L1Algo.h" #include "L1Algo/L1TrackPar.h" #include "L1Algo/L1Branch.h" #include "L1Algo/L1Extrapolation.h" #include "L1Algo/L1Filtration.h" #include "L1Algo/L1AddMaterial.h" #include "TStopwatch.h" #include #include #include using std::cout; using std::endl; using std::pair; using std::vector; using std::map; void L1Algo::CATrackFinder() { #ifdef DRAW InitL1Draw(); #endif #ifdef DEBUG_CA InitL1Histo(); #endif //cout.setf(ios::fixed); //cout.setf(ios::showpoint); cout.precision(6); //static double //stat_time = 0.0, stat_time_find = 0.0, stat_time_fit = 0.0, //stat_time_fit_bra = 0.0, stat_time_fit_fin = 0.0, //stat_time_find_10 = 0.0, stat_time_find_20 = 0.0, stat_time_find_30 = 0.0, stat_time_find_40 = 0.0, //stat_time_find_50 = 0.0, //stat_time_find_11 = 0.0, stat_time_find_21 = 0.0, stat_time_find_31 = 0.0; //static unsigned int stat_N = 0; //double //time_find = 0.0, time_fit = 0.0, time_fit_bra = 0.0, time_fit_fin = 0.0, //time_find_10 = 0.0, time_find_20 = 0.0, time_find_30 = 0.0, time_find_40 = 0.0, time_find_50 = 0.0, //time_find_11 = 0.0, time_find_21 = 0.0, time_find_31 = 0.0; const int FIRSTCASTATION = 0; TStopwatch c_time; //cout << endl << "===> CA Track Finder ..." << endl; vTracks.clear(); vRecoHits.clear(); // Arrays of duplets const int MaxArrSize = 2000000; static unsigned short int vecA0[MaxArrSize], vecB0[MaxArrSize], vecA[MaxArrSize], vecB[MaxArrSize], *mrDuplets_hits = vecA, // right hits of middle-right duplets *mrDuplets_start = vecA0,// first right hit in mrDuplets_hits, array paralel to middle hits *lmDuplets_hits = vecB, // same for left-middle duplets *lmDuplets_start = vecB0; // static bool vecFakeA[MaxArrSize], vecFakeB[MaxArrSize], *vFake_m = vecFakeA, *vFake_l = vecFakeB; static pair vecTripA0[MaxArrSize], vecTripB0[MaxArrSize], *vTriplets0_m = vecTripA0, *vTriplets0_l = vecTripB0; pair dup_init(0,0); //(*vA).clear(); (*vB).clear(); for( vector::iterator is= vSFlag.begin(); is!=vSFlag.end(); ++is) SetFUnUsed(*is); for( vector::iterator is= vSFlagB.begin(); is!=vSFlagB.end(); ++is) SetFUnUsed(*is); fvec vINF = 10.0; unsigned int stat_max_n_trip = 0, stat_max_trip_size = 0, stat_max_n_dup = 0, stat_max_BranchSize = 0, stat_max_n_branches = 0; #ifdef DRAW Draw1(); #endif //for (int isec = 0; isec <= 0; isec++){ // D0 !!! //for (int isec = 0; isec <= 1; isec++){ for (int isec = 0; isec <= 2; isec++){ // all finder //for (int isec = 0; isec <= fTrackingLevel; isec++){ // unsigned int stat_n_trip = 0; //cout << "isec " << isec << " : "; //cout<<" total STS hits "<::iterator is= vSFlag.begin(); is!=vSFlag.end(); ++is) SetFUnUsedD(*is); for( vector::iterator is= vSFlagB.begin(); is!=vSFlagB.end(); ++is) SetFUnUsedD(*is); //TStopwatch c_time_find; TRACK_CHI2_CUT = 10.0; // SG double Pick_m = 3.0; double Pick_r = 5.0; double Pick_gather = 5.0; //!!! if (isec == 0){ Pick_m = 2.0; } //double MaxInvMom = 1.0/1.; double MaxInvMom = 1.0/0.5; //if (isec == 1) MaxInvMom = 1.0/0.2; if (isec == 1) MaxInvMom = 1.0/0.1; //if (isec == 2) MaxInvMom = 1.0/0.2; if (isec == 2) MaxInvMom = 1.0/0.1; if (isec == fTrackingLevel ) MaxInvMom = 1./fMomentumCutOff; int nlevel[NStations]; for (int il = 0; il < NStations; ++il) nlevel[il] = 0; double SigmaTargetX, SigmaTargetY; // target constraint static L1FieldValue targB _fvecalignment; fvec targX= 0, targY= 0, targZ= 0; if (isec <= 1){// target targB = vtxFieldValue; if (isec ==-1) SigmaTargetX = SigmaTargetY = 0.01; // target else SigmaTargetX = SigmaTargetY = 0.1; } else{//use outer radius of the 1st station as a constraint L1Station &st = vStations[0]; SigmaTargetX = SigmaTargetY = 10;//st.Rmax[0]; targZ = 0.;//-1.; st.fieldSlice.GetFieldValue( 0, 0, targB ); } static L1XYMeasurementInfo TargetXYInfo _fvecalignment; TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX; TargetXYInfo.C10 = 0; TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY; //cout<<"Triplets found " << endl; //create triplets and find neighbors vTriplets.clear(); int nh = StsHitsStopIndex[NStations-1] - StsHitsStartIndex[NStations-1] + 1; L1StsHit *vStsHits_l = &(vStsHits[0]) + StsHitsStartIndex[NStations-1]; for (int ih =0 ; ih= 0; istal--){// //start downstream chambers for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--){// //start downstream chambers //cout << "istal " << istal << " : " << endl; TripStartIndex[istal] = TripStopIndex[istal] = vTriplets.size(); int nstaltriplets = 0; int istam = istal+1; int istar = istal+2; unsigned int nDuplets_lm = 0; L1Station &stal = vStations[istal]; L1Station &stam = vStations[istam]; L1Station &star = vStations[istar]; fvec qp1 = MaxInvMom; if( isec>=2 ) qp1 = 1; fvec zl = stal.z[0]; fvec dzli = 1./(zl-targZ[0]); fvec Infqp = MaxInvMom/3.*MaxInvMom/3.; const int Portion = 16; const int MaxPortionHit = Portion; const int MaxPortionDup = 10000; const int MaxPortionTrip = 10000; static L1FieldRegion fld_1[2000]; static L1FieldRegion fld_2[20000]; static L1TrackPar T_1[20000]; static L1TrackPar T_2[20000]; static L1TrackPar *T_3 = T_1; static unsigned short int hitsl_2[MaxPortionDup]; static unsigned short int hitsl_3[MaxPortionTrip]; static unsigned short int *hitsl_1 = hitsl_3; static unsigned short int hitsm_2[10000]; static unsigned short int hitsm_3[10000]; static unsigned short int hitsr_3[10000]; static fvec u_front[20000], u_back[20000]; static fvec x_minusV[20000], x_plusV[20000], y_minusV[20000], y_plusV[20000]; static fscal *x_minus = (fscal*) x_minusV; static fscal *x_plus = (fscal*) x_plusV; static fscal *y_minus = (fscal*) y_minusV; static fscal *y_plus = (fscal*) y_plusV; int stat_nhits=0, stat_nduplets=0, stat_n_triplet_c=0, stat_n_triplets=0; L1StsHit *vStsHits_l = &(vStsHits[0]) + StsHitsStartIndex[istal]; L1StsHit *vStsHits_m = &(vStsHits[0]) + StsHitsStartIndex[istam]; L1StsHit *vStsHits_r = 0; int NHits_l = StsHitsStopIndex[istal] - StsHitsStartIndex[istal] + 1; int NHits_m = StsHitsStopIndex[istam] - StsHitsStartIndex[istam] + 1; int NHits_r = 0; if( istar < NStations ){ vStsHits_r = &(vStsHits[0]) + StsHitsStartIndex[istar]; NHits_r = StsHitsStopIndex[istar] - StsHitsStartIndex[istar] + 1; } //hit on the middle plane to start int start_mhit = 0; for( int start_lh=0; start_lhNHits_l ) end_lh = NHits_l; int n1 = 0, n1_V=0, n1_4=0; // prepare the portion of left hits //TStopwatch c_time_find_10; for (int ilh=start_lh; ilh 0 ){ for( int i1_4=n1_4-1; n1_4= MaxPortionDup ) break; //if( imh>=10000 ) cout<<"imh="< y_plus[i1] ) break; if ( (x < x_minus[i1] ) || (x > x_plus[i1]) ) continue; lmDuplets_hits[nDuplets_lm++]=imh; if (istar >= NStations) continue; if( mrDuplets_start[imh] >= mrDuplets_start[imh+1] ) continue; L1TrackPar &T1 = T_1[i1_V]; L1TrackPar &T2 = T_2[n2_V]; L1FieldRegion &f1 = fld_1[i1_V]; L1FieldRegion &f2 = fld_2[n2_V]; u_front[n2_V][n2_4] = x; u_back [n2_V][n2_4] = v; for( fvec *i=&T2.x, *j=&T1.x; i<=&T2.NDF; i++,j++ ) (*i)[n2_4]=(*j)[i1_4]; for( fvec *i=&f2.cx0, *j=&f1.cx0; i<=&f2.z0; i++,j++ ) (*i)[n2_4]=(*j)[i1_4]; hitsl_2[n2] = hitsl_1[i1]; hitsm_2[n2] = imh; n2++; n2_V = n2/fvecLen; n2_4 = n2 % fvecLen; } //cout<0 ){ L1TrackPar &T2 = T_2[n2_V]; L1FieldRegion &f2 = fld_2[n2_V]; for( int i2_4=n2_4-1; n2_4= MaxPortionTrip ) break; int irh = mrDuplets_hits[irh_index]; //if( irh>=10000 ) cout<<"irh="< y_plus [i2] ) break; if ((xr < x_minus[i2]) || (xr > x_plus[i2])) continue; L1TrackPar &T3 = T_3[n3_V]; L1TrackPar &T2 = T_2[i2_V]; hitsl_3[n3] = hitsl_2[i2]; hitsm_3[n3] = hitsm_2[i2]; hitsr_3[n3] = irh; if( n3_4==0 ){ u_front[n3_V] = vStsStrips[hitr.f]; u_back [n3_V] = vStsStripsB[hitr.b]; for( fvec *i=&T3.x, *j=&T2.x; i<=&T3.NDF; i++,j++ ) *i=(*j)[i2_4]; }else{ u_front[n3_V][n3_4] = vStsStrips[hitr.f]; u_back [n3_V][n3_4] = vStsStripsB[hitr.b]; for( fvec *i=&T3.x, *j=&T2.x; i<=&T3.NDF; i++,j++ ) (*i)[n3_4]=(*j)[i2_4]; } n3++; n3_V = n3/fvecLen; n3_4 = n3 % fvecLen; } }// i2 if( n3_4>0 ){ for( int i3_4=n3_4-1; n3_4 TRACK_CHI2_CUT * (T_3[i3_V].NDF[i3_4]-5)) continue; //cout<<510<<" "<1 ) continue; fscal ty = fabs(T_3[i3_V].ty[i3_4]); if( ty>1 ) continue; fscal qp = MaxInvMom + T_3[i3_V].qp[i3_4]; if( qp<0 ) qp = 0; if( qp> MaxInvMom*2 ) qp = MaxInvMom*2; fscal Cqp = 5.*sqrt(fabs(T_3[i3_V].C44[i3_4])); //cout<<520<255 ) Cqp= 255; if( Cqp>20 ) Cqp= 20; qp = (unsigned char) qp; //cout<<"qp "<60000) //cout<<"ihitl="<60000) cout<<"trip.GetLHit()="<Fill(nDuplets_lm); #endif if( stat_max_n_dup *pttmp = vTriplets0_m; vTriplets0_m = vTriplets0_l; vTriplets0_l = pttmp; bool * Ftmp = vFake_m; vFake_m = vFake_l; vFake_l = Ftmp; //cout << "*** stal " << istal << " : " << nstaltriplets << endl; //cout<<" supressed hits = "<= min_level; ilev--){ //for (int ilev = NStations-3; ilev >= 0; ilev--){ //TStopwatch c_time_fit_bra; static vector vtrackcandidate; vtrackcandidate.clear(); //how many levels to check int nlevel = (NStations-2)-ilev+1; int ntracks = 0; for( int istaF = 0; istaF<= NStations-3-ilev; istaF++ ){ //cout<<"istaF="<"<GetLevel(); int jst= GetFStation( vSFlag[vStsHits[first_trip->GetLHit()].f] ); int fl = GetFUsed( vSFlag[vStsHits[first_trip->GetLHit()].f] | vSFlagB[vStsHits[first_trip->GetLHit()].b] ); if ( lv != ilev) continue; if (ilev == 0){ //remove pile-up -> collect only MAPS triplets if ( jst != 0) continue; } if( fl ) continue; // if used /* //check the level if (first_trip->GetLevel() != ilev) continue; if (ilev == 0){ //remove pile-up -> collect only MAPS triplets if ( GetFStation(vSFlag[vStsHits[first_trip->GetLHit()].f]) != 0) continue; } //check if all hits are not used => clear start if( GetFUsed( vSFlag[vStsHits[first_trip->GetLHit()].f] | vSFlagB[vStsHits[first_trip->GetLHit()].b] ) ) continue; // if used */ //cout<GetLHit()); unsigned char curr_L = 1; fscal curr_chi2 = first_trip->GetChi2(); //fscal curr_Qp = first_trip->GetQpOrig(MaxInvMom); //cout << "Triplet Qp " << trip.GetQpOrig(MaxInvMom) << " mom " << 1.0/trip.GetQpOrig(MaxInvMom) << endl; L1Branch best_tr = curr_tr; fscal best_chi2 = curr_chi2; unsigned char best_L = curr_L; int NCalls = 1; //cout<<"CAFindTrack, triplet N "<=0; i--){ if( fabs(T.qp[0])>2. ) break; L1StsHit &hit = vStsHits[best_tr.StsHits[i]]; ista = vSFlag[hit.f]/4; cout<=0; ista-- ){ L1Station &sta = vStations[ista]; L1Extrapolate( T, sta.z, qp0, fld ); fscal dxm_est = ( Pick_gather*sqrt(fabs(T.C00+sta.XYInfo.C00)) )[0]; fscal dym_est = ( Pick_gather*sqrt(fabs(T.C11+sta.XYInfo.C11)) )[0]; fscal x_minus = T.x[0] - dxm_est; fscal x_plus = T.x[0] + dxm_est; fscal y_minus = T.y[0] - dym_est; fscal y_plus = T.y[0] + dym_est; double D2 = 1.e8; int ibest = -1; for( int ih = StsHitsStartIndex[ista]; ih<=StsHitsStopIndex[ista]; ih++ ){ L1StsHit &hit = vStsHits[ih]; if( GetFUsed( vSFlag[hit.f] | vSFlagB[hit.b] ) ) continue; // if used fscal x = vStsStrips[hit.f]; fscal v = vStsStripsB[hit.b]; fscal y = (sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v)[0]; if (y < y_minus) continue; if (y > y_plus ) break; if ((x < x_minus) || (x > x_plus)) continue; double dx = x-T.x[0]; double dy = y-T.y[0]; double d2 = dx*dx + dy*dy; if( d2>D2 ) continue; D2 = d2; ibest = ih; } if( ibest <0 ) break; best_tr.StsHits.push_back(ibest); best_L++; nGathered ++; L1StsHit &hit = vStsHits[ibest]; fvec x = vStsStrips[hit.f]; fvec v = vStsStripsB[hit.b]; fvec y = sta.yInfo.cos_phi*x + sta.yInfo.sin_phi*v; L1Filter( T, sta.frontInfo, x ); L1Filter( T, sta.backInfo, v ); L1AddMaterial( T, sta.materialInfo, qp0 ); fB0 = fB1; fB1 = fB2; fz0 = fz1; fz1 = fz2; sta.fieldSlice.GetFieldValue( x, y, fB2 ); fz2 = sta.z; fld.Set( fB2, fz2, fB1, fz1, fB0, fz0 ); } } //if( nGathered>0 ) cout< TRACK_CHI2_CUT *ndf ) continue; //if ( fabs(T.qp[0]) > MaxInvMom ) continue; if( best_L < min_level+3 ) continue; if( fGhostSuppression ){//suppress ghost if( best_L <3 ) continue; if( best_L == 3 ){ //if( isec == 2 ) continue; // too short secondary track if( isec==2 && istaF != 0 ) continue; // too short non-MAPS track if( (isec<2)&&(best_chi2>5.0) ) continue; } } // track candidate is OK, store best_tr.Set( istaF, best_L, best_chi2, first_trip->GetQpOrig(MaxInvMom)); vtrackcandidate.push_back(best_tr); ntracks++; } //cout<<" level "< vptrackcandidate; vptrackcandidate.clear(); for (vector::iterator trIt = vtrackcandidate.begin(); trIt != vtrackcandidate.end(); ++trIt){ vptrackcandidate.push_back(&(*trIt)); } //if (isec <= 1) sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePChi2); //else //sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePMomentum); ntracks = 0; for (vector::iterator trIt = vptrackcandidate.begin(); trIt != vptrackcandidate.end(); ++trIt){ //count number of used hits L1Branch *tr = *trIt; //cout << "Branch isec " << isec << " Nhits " << tr->StsHits.size() << " Quality " << tr->Quality << " chi2 " << tr-> chi2 << " momentum " << tr->Momentum << endl; int nused = 0; for (vector::iterator phIt = tr->StsHits.begin(); phIt != tr->StsHits.end(); ++phIt){ L1StsHit &h = vStsHits[*phIt]; if ( GetFUsed( vSFlag[h.f] | vSFlagB[h.b] ) ) nused++; } if (nused != 0) continue; //======================================================= // gather MAPS hits using the Kalman filter // was here, skipped for a moment //======================================================= //store track for (vector::iterator phitIt = tr->StsHits.begin(); phitIt != tr->StsHits.end(); ++phitIt){ L1StsHit &h = vStsHits[*phitIt]; SetFUsed( vSFlag[h.f] ); SetFUsed( vSFlagB[h.b] ); vRecoHits.push_back(*phitIt); } //store track L1Track t; t.NHits = tr->StsHits.size(); t.Momentum = tr->Momentum; for( int i=0; i<6; i++) t.TFirst[i] = t.TLast[i]=0; for( int i=0; i<15; i++ ) t.CFirst[i] = t.CLast[i] = 10; t.chi2 = 100; vTracks.push_back(t); ntracks++; } //c_time_fit_fin.Stop(); //time_fit_fin += double(c_time_fit_fin.CpuTime()); unsigned int Bsize=0, nb=0; for (vector::iterator trIt = vptrackcandidate.begin(); trIt != vptrackcandidate.end(); ++trIt){ L1Branch *tr = *trIt; Bsize+= sizeof(L1Branch) + sizeof(unsigned short int) * tr->StsHits.size(); nb++; } if( Bsize>stat_max_BranchSize ) stat_max_BranchSize = Bsize; if( nb>stat_max_n_branches ) stat_max_n_branches = nb; //for (vector::iterator ih = vStsHits.begin(); ih != vStsHits.end(); ++ih){ //if( ih->strip_f->used || ih->strip_b->used) ih->used = 1; //} //cout<< " track candidates stored: "<< ntracks < ... CA Track Finder " << endl << endl; #endif #ifdef DRAW DrawAsk(); #endif } //1234567890123456789012345678901234567890123456789012345678901234567890 void L1Algo::CAFindTrack( int ista, const L1Triplet* curr_trip, L1Branch& best_tr, unsigned char &best_L, fscal &best_chi2, L1Branch &curr_tr, unsigned char &curr_L, fscal &curr_chi2, int &NCalls ){ /**************************************************************** * * * The routine performs recusrsive search for tracks * * * * I. Kisel 06.03.05 * * * ***************************************************************/ // cout << "CA track finder " << curr_trip->level << endl; if (curr_trip->GetLevel() == 0){ // the last triplet -> check and store if ((curr_L > best_L )|| ((curr_L == best_L )&&(curr_chi2level != 0 //store hits & triplets of the current track L1Branch temp_tr = curr_tr; unsigned char temp_L = curr_L; fscal temp_chi2 = curr_chi2; int offset = TripStartIndex[ista+1]; int first_neighbour = offset + curr_trip->GetFirstNeighbour(); int end_neighbour = first_neighbour + curr_trip->GetNNeighbours(); for(int index = first_neighbour; index < end_neighbour; ++index){ L1Triplet *new_trip = &vTriplets[index]; if (curr_trip->GetLevel() != new_trip->GetLevel()+1) continue; fscal qp1 = curr_trip->GetQp(); fscal qp2 = new_trip->GetQp(); fscal dqp = fabs(qp1 - qp2); fscal Cqp = curr_trip->Cqp; Cqp += new_trip->Cqp; if ( dqp > Cqp ) continue; //bad neighbour //no used hits allowed if ( GetFUsed( vSFlag[vStsHits[new_trip->GetLHit()].f] | vSFlagB[vStsHits[new_trip->GetLHit()].b] ) ) { if((curr_L > best_L)|| (curr_L == best_L)&&(curr_chi2GetLHit()); curr_L+= 1; dqp = dqp/Cqp*5.; curr_chi2 += dqp*dqp; if( new_trip->GetLevel()==0 ){ unsigned short int ihitm = new_trip->GetMHit(); unsigned short int ihitr = new_trip->GetRHit(); if( !GetFUsed( vSFlag[vStsHits[ihitm].f] | vSFlagB[vStsHits[ihitm].b] ) ){ curr_tr.StsHits.push_back(ihitm); curr_L+= 1; } if( !GetFUsed( vSFlag[vStsHits[ihitr].f] | vSFlagB[vStsHits[ihitr].b] ) ){ curr_tr.StsHits.push_back(ihitr); curr_L+= 1; } } if( curr_chi2 <= TRACK_CHI2_CUT * (curr_L) ){ // go further CAFindTrack(ista+1,new_trip, best_tr, best_L, best_chi2, curr_tr, curr_L, curr_chi2,NCalls); NCalls++; } }//good triplet added }//first_triplet --- last_triplet }//level != 0 //cout << "Track finder finished " << endl; return; } //========================================================