/* *==================================================================== * * CBM Level 1 Reconstruction * * Authors: I.Kisel, S.Gorbunov * * e-mail : ikisel@kip.uni-heidelberg.de * *==================================================================== * * Read Event information to L1 * *==================================================================== */ #include "CbmL1.h" #include "L1Algo/L1Algo.h" #include "CbmKF.h" #include "CbmMvdHitMatch.h" #include #include using std::cout; using std::endl; using std::vector; struct TmpHitSort{ int ID; double z; int StsHit; CbmL1MCPoint *MCPoint; bool eff; static bool compareIDz( const TmpHitSort &a, const TmpHitSort &b ) { return ( a.ID < b.ID ) || ( ( a.ID == b.ID ) && (a.z < b.z) ); } }; struct TmpHit{ int iStripF, iStripB; int indStripF, indStripB; int iSector, iStation; int ExtIndex; bool isStrip; double u_front, u_back; double x, y; int iMC; static bool Compare( const TmpHit &a, const TmpHit &b ) { return ( a.iStation < b.iStation ) || ( ( a.iStation == b.iStation ) && ( a.y < b.y ) ); } }; struct TmpStrip{ fscal u; int iStation; int iSector, iStrip; bool isStrip; int effIndex; }; void CbmL1::ReadEvent() { vMCPoints.clear(); vMCTracks.clear(); vRTracks.clear(); vHitMCRef.clear(); vHitStore.clear(); algo->vStsHits.clear(); algo->vStsStrips.clear(); algo->vStsStripsB.clear(); algo->vSFlag.clear(); algo->vSFlagB.clear(); vector tmpHits; vector tmpStrips; vector tmpStripsB; // produce Sts hits from space points for(int i = 0; i < NStation; i++){ algo->StsHitsStartIndex[i] = -1; algo->StsHitsStopIndex[i] = -2; } // get MVD hits int nMvdHits=0; if( listMvdHits ){ Int_t nEnt = listMvdHits->GetEntries(); for (int j=0; j At(j); TmpHit th; { CbmMvdHit *mh = (CbmMvdHit*) listMvdHits->At(j); th.ExtIndex = -(1+j); th.iStation = mh->GetStationNr() - 1; th.iSector = 0; th.isStrip = 0; th.iStripF = j; th.iStripB = -1; if( th.iStripF<0 ) continue; if( th.iStripF>=0 && th.iStripB>=0 ) th.isStrip = 1; if( th.iStripB <0 ) th.iStripB = th.iStripF; TVector3 pos, err; mh->Position(pos); mh->PositionError(err); th.x = pos.X(); th.y = pos.Y(); L1Station &st = algo->vStations[th.iStation]; th.u_front = th.x*st.frontInfo.cos_phi[0] + th.y*st.frontInfo.sin_phi[0]; th.u_back = th.x* st.backInfo.cos_phi[0] + th.y* st.backInfo.sin_phi[0]; } th.iMC=-1; int iMC = -1; if( listMvdHitMatches ){ CbmMvdHitMatch *match = (CbmMvdHitMatch*) listMvdHitMatches->At(j); if( match){ iMC = match->GetPointId(); } } if( listStsPts && iMC>=0 ){ CbmL1MCPoint MC; if( ! ReadMCPoint( &MC, iMC, 1 ) ){ MC.iStation = th.iStation; vMCPoints.push_back( MC ); th.iMC = vMCPoints.size()-1; } } //if( h.MC_Point >=0 ) // DEBUG !!!! { tmpHits.push_back(th); nMvdHits++; } } } int nStsHits=0; if( listStsHits ){ Int_t nEnt = listStsHits->GetEntries(); //cout<<"nEnt = "<GetEntries(); int negF=0; for (int j=0; j < nEnt; j++ ) { CbmStsHit *sh = (CbmStsHit*) listStsHits->At(j); TmpHit th; { CbmStsHit *mh = (CbmStsHit*) listStsHits->At(j); th.ExtIndex = j; th.iStation = NMvdStations + mh->GetStationNr() - 1; th.iSector = mh->GetSectorNr(); th.isStrip = 0; th.iStripF = mh->GetDigi(0); th.iStripB = mh->GetDigi(1); if( th.iStripF<0 ){ negF++; continue;} if( th.iStripF>=0 && th.iStripB>=0 ) th.isStrip = 1; if( th.iStripB <0 ) th.iStripB = th.iStripF; th.iStripF += nMvdHits; th.iStripB += nMvdHits; TVector3 pos, err; mh->Position(pos); mh->PositionError(err); th.x = pos.X(); th.y = pos.Y(); L1Station &st = algo->vStations[th.iStation]; th.u_front = th.x*st.frontInfo.cos_phi[0] + th.y*st.frontInfo.sin_phi[0]; th.u_back = th.x* st.backInfo.cos_phi[0] + th.y* st.backInfo.sin_phi[0]; } th.iMC=-1; int iMC = sh->GetRefIndex(); if( listStsPts && iMC>=0 && iMC=0 ) // DEBUG !!!! { tmpHits.push_back(th); nStsHits++; } } //cout<<"negF="<1.e-4 ) continue; th.indStripB = is; } if( th.indStripF<0 ){ TmpStrip tmp; tmp.iStation = th.iStation; tmp.iSector = th.iSector; tmp.iStrip = th.iStripF; tmp.u = th.u_front; tmp.isStrip = th.isStrip; tmpStrips.push_back(tmp); th.indStripF = NStrips++; } if( th.indStripB<0 ){ TmpStrip tmp; tmp.iStation = th.iStation; tmp.iSector = th.iSector; tmp.iStrip = th.iStripB; tmp.isStrip = th.isStrip; tmp.u = th.u_back; tmpStripsB.push_back(tmp); th.indStripB = NStripsB++; } } //cout <<"N strips = "<= NStrips ) continue; if( th.indStripB <0 || th.indStripB >= NStripsB ) continue; if( th.iMC<0 ) continue; if( gRandom->Uniform(1)>fDetectorEfficiency ){ tmpStrips [th.indStripF].effIndex = -2; tmpStripsB[th.indStripB].effIndex = -2; } } } else { // strip unefficiency for( int i=0; iUniform(1)>fDetectorEfficiency ) tmpStrips[i].effIndex = -2; } for( int i=0; iUniform(1)>fDetectorEfficiency ) ts.effIndex = -2; } } Int_t NEffStrips = 0, NEffStripsB = 0; for( int i=0; ivStsStrips.push_back(ts.u); algo->vSFlag.push_back(flag); } } for( int i=0; ivStsStripsB.push_back(ts.u); algo->vSFlagB.push_back(flag); } } vector vUnEffHitMCRef; vector vUnEffHitStore; int nEffHits = 0; for( int i=0; i= NStrips ) continue; if( th.indStripB <0 || th.indStripB >= NStripsB ) continue; TmpStrip &tsF = tmpStrips [th.indStripF]; TmpStrip &tsB = tmpStripsB[th.indStripB]; if( tsF.effIndex <0 || tsB.effIndex <0 ){ vUnEffHitStore.push_back(s); vUnEffHitMCRef.push_back(th.iMC); continue; } L1StsHit h; h.f = tsF.effIndex; h.b = tsB.effIndex; algo->vStsHits.push_back(h); int sta = th.iStation; if (algo->StsHitsStartIndex[sta] < 0) algo->StsHitsStartIndex[sta] = nEffHits; algo->StsHitsStopIndex[sta] = nEffHits; nEffHits ++; vHitStore.push_back(s); vHitMCRef.push_back(th.iMC); } //cout<<" sorting MC points.."< vTmp; for ( int i=0; i::iterator i= vTmp.begin(); i!=vTmp.end(); ++i) { if( vMCTracks.empty() || i->ID!= ID ) { //cout<ID; CbmL1MCTrack T; T.ID = ID; if( listMCTracks ) { T.q = i->MCPoint->q; T.mass = i->MCPoint->mass; T.pdg = i->MCPoint->pdg; CbmMCTrack *MCTrack = (CbmMCTrack*) listMCTracks->At( T.ID ); T.mother_ID = MCTrack->GetMotherId(); TVector3 vr; MCTrack->GetStartVertex(vr); TLorentzVector vp; MCTrack->Get4Momentum(vp); T.px = vp.Px(); T.py = vp.Py(); T.pz = vp.Pz(); T.x = vr.X(); T.y = vr.Y(); T.z = vr.Z(); T.p = sqrt( fabs(T.px*T.px + T.py*T.py + T.pz*T.pz )); } vMCTracks.push_back( T ); if ( T.mother_ID ==-1 ){ // vertex track if ( PrimVtx.MC_ID == 999 || fabs(T.z-Vtxcurr.MC_z)>1.e-7 ){// new vertex if( nvtrackscurr > nvtracks ){ PrimVtx = Vtxcurr; nvtracks = nvtrackscurr; } Vtxcurr.MC_x = T.x; Vtxcurr.MC_y = T.y; Vtxcurr.MC_z = T.z; Vtxcurr.MC_ID = T.mother_ID; nvtrackscurr = 1; }else nvtrackscurr++; } } vMCTracks.back().Points.push_back(i->MCPoint); if( i->eff ) vMCTracks.back().StsHits.push_back(i->StsHit); } if( nvtrackscurr > nvtracks ) PrimVtx = Vtxcurr; //cout<::iterator it = vMCTracks.begin(); it != vMCTracks.end(); ++it){ CbmL1MCTrack &T = *it; int nPoints = T.Points.size(); T.nMCContStations = 0; int istaold = -1, ncont=0; for( int ih=0; ihvStsHits[jh]; int ista = algo->vSFlag[h.f]/4; if (ista - istaold == 1) ncont++; else { if( T.nHitContStationsAt(iPoint); if ( !pt ) return 1; pt->Position(xyzI); pt->Momentum(PI); pt->PositionOut(xyzO); pt->MomentumOut(PO); mcID = pt->GetTrackID(); }else if( listStsPts ){ CbmStsPoint *pt = (CbmStsPoint*) listStsPts->At(iPoint); if ( !pt ) return 1; pt->Position(xyzI); pt->Momentum(PI); pt->PositionOut(xyzO); pt->MomentumOut(PO); mcID = pt->GetTrackID(); } TVector3 xyz = .5*(xyzI + xyzO ); TVector3 P = .5*(PI + PO ); MC->x = xyz.X(); MC->y = xyz.Y(); MC->z = xyz.Z(); MC->px = P.X(); MC->py = P.Y(); MC->pz = P.Z(); MC->p = sqrt( fabs( MC->px*MC->px + MC->py*MC->py + MC->pz*MC->pz ) ); MC->ID = mcID; CbmMCTrack *MCTrack = (CbmMCTrack*) listMCTracks->At( MC->ID ); if ( !MCTrack ) return 1; MC->pdg = MCTrack->GetPdgCode(); MC->mother_ID = MCTrack->GetMotherId(); if ( abs(MC->pdg) >= 10000 ) return 1; MC->q = TDatabasePDG::Instance()->GetParticle(MC->pdg)->Charge()/3.0; MC->mass = TDatabasePDG::Instance()->GetParticle(MC->pdg)->Mass(); return 0; }