/* *===================================================== * * CBM Level 1 Reconstruction * * Authors: M.Zyzak * * e-mail : * *===================================================== * * Merge Clones * */ #include "L1Algo.h" #include "L1Extrapolation.h" #include // using namespace std; using std::cout; using std::endl; using std::vector; void L1Algo::InvertCholetsky(fvec a[15]) { fvec d[5], uud, u[5][5]; for(int i=0; i<5; i++) { d[i]=0.f; for(int j=0; j<5; j++) u[i][j]=0.; } for(int i=0; i<5; i++) { uud=0.; for(int j=0; j FirstHit; vector LastHit; vector FirstHitIndex; vector LastHitIndex; vector Neighbour; vector TrackChi2; vector IsNext; vector IsUsed; vector< L1Track > vTracksNew; vTracksNew.reserve(vTracks.size()); vector< THitI > vRecoHitsNew; vRecoHitsNew.reserve(vRecoHits.size()); FirstHit.resize(vTracks.size()); LastHit.resize(vTracks.size()); FirstHitIndex.resize(vTracks.size()); LastHitIndex.resize(vTracks.size()); IsUsed.resize(vTracks.size()); TrackChi2.resize(vTracks.size()); Neighbour.resize(vTracks.size()); IsNext.resize(vTracks.size()); THitI start_hit = 0; unsigned short ista = 0; for(unsigned short iTr = 0; iTr < vTracks.size(); iTr++) { FirstHitIndex[iTr] = start_hit; ista = vSFlag[vStsHits[vRecoHits[start_hit]].f]/4; FirstHit[iTr]=ista; start_hit += vTracks[iTr].NHits-1; LastHitIndex[iTr] = start_hit; ista = vSFlag[vStsHits[vRecoHits[start_hit]].f]/4; LastHit[iTr]=ista; start_hit++; IsUsed[iTr] = 0; Neighbour[iTr] = 50000; TrackChi2[iTr] = 100000; IsNext[iTr] = 0; } L1KFTrackFitter(); //KFTrackFitter_simple(); for(int iTr = 0; iTr < static_cast(vTracks.size()); iTr++) { if(static_cast(vTracks[iTr].NHits) > 6) continue; for(int jTr = 0; jTr < static_cast(vTracks.size()); jTr++) { if(iTr == jTr) continue; if(static_cast(vTracks[iTr].NHits) > 6) continue; unsigned short dist = 0; L1TrackPar Tb; L1TrackPar Tf; unsigned short stab=0, staf=0; bool IsNextTemp=0; //if((vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])*(vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4]) // > 9*(vTracks[iTr].CFirst[14]+vTracks[jTr].CFirst[14]) ) continue; if(FirstHit[iTr] > LastHit[jTr]) { dist = FirstHit[iTr] - LastHit[jTr]; stab = FirstHit[iTr]; staf = LastHit[jTr]; IsNextTemp = 1; Tb.x = vTracks[iTr].TFirst[0]; Tb.y = vTracks[iTr].TFirst[1]; Tb.tx = vTracks[iTr].TFirst[2]; Tb.ty = vTracks[iTr].TFirst[3]; Tb.qp = vTracks[iTr].TFirst[4]; Tb.z = vTracks[iTr].TFirst[5]; Tb.C00 = vTracks[iTr].CFirst[0]; Tb.C10 = vTracks[iTr].CFirst[1]; Tb.C11 = vTracks[iTr].CFirst[2]; Tb.C20 = vTracks[iTr].CFirst[3]; Tb.C21 = vTracks[iTr].CFirst[4]; Tb.C22 = vTracks[iTr].CFirst[5]; Tb.C30 = vTracks[iTr].CFirst[6]; Tb.C31 = vTracks[iTr].CFirst[7]; Tb.C32 = vTracks[iTr].CFirst[8]; Tb.C33 = vTracks[iTr].CFirst[9]; Tb.C40 = vTracks[iTr].CFirst[10]; Tb.C41 = vTracks[iTr].CFirst[11]; Tb.C42 = vTracks[iTr].CFirst[12]; Tb.C43 = vTracks[iTr].CFirst[13]; Tb.C44 = vTracks[iTr].CFirst[14]; Tf.x = vTracks[jTr].TLast[0]; Tf.y = vTracks[jTr].TLast[1]; Tf.tx = vTracks[jTr].TLast[2]; Tf.ty = vTracks[jTr].TLast[3]; Tf.qp = vTracks[jTr].TLast[4]; Tf.z = vTracks[jTr].TLast[5]; Tf.C00 = vTracks[jTr].CLast[0]; Tf.C10 = vTracks[jTr].CLast[1]; Tf.C11 = vTracks[jTr].CLast[2]; Tf.C20 = vTracks[jTr].CLast[3]; Tf.C21 = vTracks[jTr].CLast[4]; Tf.C22 = vTracks[jTr].CLast[5]; Tf.C30 = vTracks[jTr].CLast[6]; Tf.C31 = vTracks[jTr].CLast[7]; Tf.C32 = vTracks[jTr].CLast[8]; Tf.C33 = vTracks[jTr].CLast[9]; Tf.C40 = vTracks[jTr].CLast[10]; Tf.C41 = vTracks[jTr].CLast[11]; Tf.C42 = vTracks[jTr].CLast[12]; Tf.C43 = vTracks[jTr].CLast[13]; Tf.C44 = vTracks[jTr].CLast[14]; //std::cout << "!!!!!!! Chi2 !!!!!! "< LastHit[iTr]) { dist = FirstHit[jTr] - LastHit[iTr]; stab = FirstHit[jTr]; staf = LastHit[iTr]; Tb.x = vTracks[jTr].TFirst[0]; Tb.y = vTracks[jTr].TFirst[1]; Tb.tx = vTracks[jTr].TFirst[2]; Tb.ty = vTracks[jTr].TFirst[3]; Tb.qp = vTracks[jTr].TFirst[4]; Tb.z = vTracks[jTr].TFirst[5]; Tb.C00 = vTracks[jTr].CFirst[0]; Tb.C10 = vTracks[jTr].CFirst[1]; Tb.C11 = vTracks[jTr].CFirst[2]; Tb.C20 = vTracks[jTr].CFirst[3]; Tb.C21 = vTracks[jTr].CFirst[4]; Tb.C22 = vTracks[jTr].CFirst[5]; Tb.C30 = vTracks[jTr].CFirst[6]; Tb.C31 = vTracks[jTr].CFirst[7]; Tb.C32 = vTracks[jTr].CFirst[8]; Tb.C33 = vTracks[jTr].CFirst[9]; Tb.C40 = vTracks[jTr].CFirst[10]; Tb.C41 = vTracks[jTr].CFirst[11]; Tb.C42 = vTracks[jTr].CFirst[12]; Tb.C43 = vTracks[jTr].CFirst[13]; Tb.C44 = vTracks[jTr].CFirst[14]; Tf.x = vTracks[iTr].TLast[0]; Tf.y = vTracks[iTr].TLast[1]; Tf.tx = vTracks[iTr].TLast[2]; Tf.ty = vTracks[iTr].TLast[3]; Tf.qp = vTracks[iTr].TLast[4]; Tf.z = vTracks[iTr].TLast[5]; Tf.C00 = vTracks[iTr].CLast[0]; Tf.C10 = vTracks[iTr].CLast[1]; Tf.C11 = vTracks[iTr].CLast[2]; Tf.C20 = vTracks[iTr].CLast[3]; Tf.C21 = vTracks[iTr].CLast[4]; Tf.C22 = vTracks[iTr].CLast[5]; Tf.C30 = vTracks[iTr].CLast[6]; Tf.C31 = vTracks[iTr].CLast[7]; Tf.C32 = vTracks[iTr].CLast[8]; Tf.C33 = vTracks[iTr].CLast[9]; Tf.C40 = vTracks[iTr].CLast[10]; Tf.C41 = vTracks[iTr].CLast[11]; Tf.C42 = vTracks[iTr].CLast[12]; Tf.C43 = vTracks[iTr].CLast[13]; Tf.C44 = vTracks[iTr].CLast[14]; } if(dist == 0) continue; //if(((Tf.qp - Tb.qp)*(Tf.qp - Tb.qp)/(Tb.C44+Tf.C44))[0] > 25*10*7) continue; L1FieldValue fBm, fBb, fBf _fvecalignment; L1FieldRegion fld _fvecalignment; unsigned short stam; vStations[staf].fieldSlice.GetFieldValue( Tf.x, Tf.y, fBf ); vStations[stab].fieldSlice.GetFieldValue( Tb.x, Tb.y, fBb ); if(dist > 1) stam = staf + 1; else stam = staf - 1; fvec zm = vStations[stam].z; fvec xm = 0.5*(Tf.x + Tf.tx*(zm - Tf.z) + Tb.x + Tb.tx*(zm - Tb.z)); fvec ym = 0.5*(Tb.y + Tb.ty*(zm - Tb.z) + Tb.y + Tb.ty*(zm - Tb.z)); vStations[stam].fieldSlice.GetFieldValue( xm, ym, fBm ); fld.Set( fBb, Tb.z, fBm, zm, fBf, Tf.z ); fvec zMiddle = 0.5*(Tb.z + Tf.z); L1Extrapolate( Tf, zMiddle, Tf.qp, fld ); L1Extrapolate( Tb, zMiddle, Tb.qp, fld ); fvec Chi2Tracks = 0.f; fvec rf[5] = {Tf.x,Tf.y,Tf.tx,Tf.ty,Tf.qp}; fvec rb[5] = {Tb.x,Tb.y,Tb.tx,Tb.ty,Tb.qp}; fvec Cf[15] = {Tf.C00,Tf.C10,Tf.C11,Tf.C20,Tf.C21,Tf.C22,Tf.C30,Tf.C31,Tf.C32,Tf.C33,Tf.C40,Tf.C41,Tf.C42,Tf.C43,Tf.C44}; fvec Cb[15] = {Tb.C00,Tb.C10,Tb.C11,Tb.C20,Tb.C21,Tb.C22,Tb.C30,Tb.C31,Tb.C32,Tb.C33,Tb.C40,Tb.C41,Tb.C42,Tb.C43,Tb.C44}; FilterTracks(rf,Cf,rb,Cb,0,0,&Chi2Tracks); if(Chi2Tracks[0] > 50 ) continue; if(Chi2Tracks[0] < TrackChi2[iTr] || Chi2Tracks[0] < TrackChi2[jTr]) { //std::cout << "!!!!!!! Chi2 !!!!!! "<(50000)) { Neighbour[Neighbour[iTr]] = 50000; TrackChi2[Neighbour[iTr]] = 100000; IsNext[Neighbour[iTr]] = 0; } if(Neighbour[jTr] < static_cast(50000)) { Neighbour[Neighbour[jTr]] = 50000; TrackChi2[Neighbour[jTr]] = 100000; IsNext[Neighbour[jTr]] = 0; } Neighbour[iTr] = jTr; Neighbour[jTr] = iTr; TrackChi2[iTr] = Chi2Tracks[0]; TrackChi2[jTr] = Chi2Tracks[0]; IsNext[iTr] = IsNextTemp; IsNext[jTr] = (!IsNextTemp); } } } for(int iTr = 0; iTr < static_cast(vTracks.size()); iTr++) { if(IsUsed[iTr]) continue; vTracksNew.push_back(vTracks[iTr]); if(!IsNext[iTr]) for(THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++) vRecoHitsNew.push_back(vRecoHits[HI]); if(Neighbour[iTr] < 50000) { IsUsed[Neighbour[iTr]] = 1; vTracksNew.back().NHits += vTracks[Neighbour[iTr]].NHits; for(THitI HI = FirstHitIndex[Neighbour[iTr]]; HI <= LastHitIndex[Neighbour[iTr]]; HI++) vRecoHitsNew.push_back(vRecoHits[HI]); } if(IsNext[iTr]) for(THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++) vRecoHitsNew.push_back(vRecoHits[HI]); } vTracks.resize(vTracksNew.size()); for(unsigned short iTr=0; iTr < vTracksNew.size(); iTr++) vTracks[iTr] = vTracksNew[iTr]; vRecoHits.resize(vRecoHitsNew.size()); for(THitI iHi=0; iHi < vRecoHitsNew.size(); iHi++) vRecoHits[iHi] = vRecoHitsNew[iHi]; //std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!! new "<