16#include "MdcGeom/Constants.h"
17#include "MdcGeom/BesAngle.h"
18#include "MdcData/MdcHit.h"
19#include "MdcData/MdcHitOnTrack.h"
20#include "MdcGeom/MdcDetector.h"
21#include "MdcTrkRecon/MdcTrackListBase.h"
22#include "MdcTrkRecon/MdcMap.h"
23#include "MdcTrkRecon/MdcTrack.h"
24#include "TrkBase/TrkErrCode.h"
25#include "TrkBase/TrkRecoTrk.h"
26#include "TrkBase/TrkFit.h"
27#include "TrkBase/TrkFitStatus.h"
28#include "TrkBase/TrkHitList.h"
29#include "TrkBase/TrkExchangePar.h"
30#include "MdcRecoUtil/Pdt.h"
31#include "MdcGeom/BesAngle.h"
32#include "MdcData/MdcRecoHitOnTrack.h"
33#include "MdcRawEvent/MdcDigi.h"
34#include "CLHEP/Matrix/SymMatrix.h"
35#include "CLHEP/Vector/ThreeVector.h"
36#include "TrkFitter/TrkContextEv.h"
37#include "TrkBase/TrkRep.h"
38#include "Identifier/MdcID.h"
39#include "Identifier/Identifier.h"
40#include "AIDA/IHistogram1D.h"
41#include "AIDA/IHistogram2D.h"
42#include "CLHEP/Geometry/Point3D.h"
43#include "CLHEP/Geometry/Vector3D.h"
44#include "TrkBase/TrkDifTraj.h"
46#ifndef ENABLE_BACKWARDS_COMPATIBILITY
51#ifndef ENABLE_BACKWARDS_COMPATIBILITY
60#include "GaudiKernel/NTuple.h"
78 for (
int itrack = 0; itrack <
nTrack(); itrack++) {
82 track->
storeTrack(trackId, trackList, hitList, tkStat);
85 HepAListDeleteAll(*
this);
94 std::cout<<
"nTrack "<<
nTrack() << std::endl;
95 for (
int itrack = 0; itrack <
nTrack(); itrack++) {
97 if (atrack == NULL)
continue;
116 std::cout <<
"=======Print before arbitrateHits=======" << std::endl;
120 std::vector<MdcTrack*> trksToKill;
121 trksToKill.reserve(4);
126 int* usedInTrackNum =
new int [
nTrack()];
130 int *refitTrack =
new int [
nTrack()];
131 for (
int i = 0; i <
nTrack(); i++) {
137 for (itrack = 0; itrack <
nTrack(); itrack++) {
139 if (atrack == 0)
continue;
141 trkXRef[itrack] = atrack;
144 for (itrack = 0; itrack <
nTrack(); itrack++) {
146 if (8 ==
tkParam.
lPrint) std::cout<<
"arbitrate track No."<<itrack<< std::endl;
148 if (atrack == 0)
continue;
155 assert (hitList != 0);
157 for (
int ii = 0; ii <
nTrack(); ii++) usedInTrackNum[ii] = 0;
162 int maxGapLength = 0;
166 int nDeleteInLayer[43];
167 for(
int i=0;i<43;i++){
171 if(8 ==
tkParam.
lPrint) std::cout<<
"--arbitrate--"<<std::endl;
173 int nUsed = ihit->hit()->nUsedHits();
175 std::cout<<
"nUsed="<<nUsed<<
":";
176 ihit->hit()->printAll(std::cout);
179 double deltaChi = -999;
180 ihit->getFitStuff(deltaChi);
181 std::cout<<
"deltaChi="<<deltaChi<<std::endl;
183 int layer = ihit->layerNumber();
184 nHitInLayer[layer]++;
186 if (!ihit->isActive()) {
191 nDeleteInLayer[layer]++;
193 std::cout<<
"=remove above inactive "<<std::endl;
197 if(ihit == hitList->
end())
break;
203 bool wasUsed =
false;
204 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q =
205 ihit->hit()->getUsedHits();
207 if ( !i->isActive() )
continue;
209 int id = recoTrk->
id();
210 if (
id == aRecoTrk.
id())
continue;
212 idMap.
get(
id, index);
214 usedInTrackNum[index]++;
216 std::cout<<
" track "<<itrack<<
"&" <<index
217 <<
" shared hits "<<usedInTrackNum[index]<<
":";
218 ihit->printAll(std::cout);
222 if (wasUsed) nPrev++;
228 for (
int i=0;i<43;i++){
231 std::cout<<i<<
" nHitInLayer "<<nHitInLayer[i]
232 <<
" nDeleteInLayer "<<nDeleteInLayer[i]<<std::endl;
235 if(nHitInLayer[i]>0 && (nHitInLayer[i]-nDeleteInLayer[i])==0) {
239 cout <<
"rec hits have been deleted in this layer"<<std::endl;
243 }
else if(nHitInLayer[i]==0){
252 if(testGap>=3){ nGapGE3++; }
253 if(testGap>maxGapLength) maxGapLength=testGap;
260 bool toBeDeleted =
false;
262 if(
tkParam.
lPrint>1) std::cout<<
"arbitrateHits tkNo:"<<itrack<<
" nGapGE2= "<<nGapGE2 <<
" nGapGE3= "<<nGapGE3 <<
" maxGapLength= "<<maxGapLength<<std::endl;
268 <<
" Killing tkNo " << itrack << endl;
277 cout <<
"arbitrateHits: nGapGE2 "<<nGapGE2<<
" >= "<<
tkParam.
nGapGE2 <<
" Killing tkNo " << itrack << endl;
283 cout <<
"arbitrateHits: nGapGE3 "<<nGapGE3<<
" >= "<<
tkParam.
nGapGE3 <<
" Killing tkNo " << itrack << endl;
289 cout <<
"arbitrateHits: maxGapLength "<<maxGapLength<<
" >= "<<
tkParam.
maxGapLength<<
" Killing tkNo " << itrack << endl;
296 delete &(atrack->
track());
298 trksToKill.push_back(atrack);
306 for (
int ii = 0; ii <
nTrack(); ii++) {
308 std::cout<<
"tk:"<<itrack<<
"&"<<ii
309 <<
" shared "<<usedInTrackNum[ii]<<
" hits "<< std::endl;
311 if (usedInTrackNum[ii] > nMost) {
312 nMost = usedInTrackNum[ii];
318 if (trackMost == trackOld) {
319 std::cout <<
"ErrMsg(error) MdcTrackListBase:"
320 <<
"Something ghastly happened in MdcTrackListBase::arbitrateHits"
324 trackOld = trackMost;
329 double groupDiff = 0.0;
338 std::cout<<
"track "<<trackMost<<
" shared "<<nMost<<
" hits > Cut nOverlap "
350 if(8 ==
tkParam.
lPrint) std::cout<<
"Go back through hits, looking up overlap hits"<< std::endl;
354 int nUsed = ihit->hit()->nUsedHits();
357 std::cout<<
"--hit go back, nUsed="<<nUsed<<
":";
358 ihit->hit()->printAll(std::cout);
362 if (nUsed < 2) {
continue; }
365 if (!ihit->isActive()) {
366 if (8 ==
tkParam.
lPrint){ std::cout<<
"act=0 continue"<<std::endl; }
371 std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator>
q = ihit->hit()->getUsedHits();
372 while (
q.first!=
q.second) {
377 if (!otherHot->
isActive())
continue;
380 if ( &aRecoTrk == otherTrack)
continue;
381 int otherId = otherTrack->
id();
382 long otherIndex = -1;
383 idMap.
get(otherId, otherIndex); assert(otherIndex >= 0);
386 if (lGroupHits && otherIndex != trackMost)
continue;
390 std::cout<<
"group hits "<< std::endl;
396 int aDof = tkFit->
nActive() - 5;
399 if (aDof <= 0) {groupDiff = 999;}
400 else if (otherDof <= 0) {groupDiff = -999;}
402 groupDiff += ihit->resid(0) * ihit->resid(0) * ihit->weight() /
407 theseHits[nFound] =
const_cast<TrkHitOnTrk*
>(ihit.get());
408 thoseHits[nFound] = otherHot;
414 std::cout<<
"handle hits individually"<< std::endl;
417 if (fabs(ihit->resid(0)) > fabs(otherHot->
resid(0)) ) {
422 const_cast<TrkHitOnTrk*
>(ihit.get())->setActivity(0);
425 std::cout<<
"dorp hit ";
426 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
431 refitTrack[otherIndex] = 1;
435 std::cout<<
"inactive hit on other track";
436 const_cast<TrkHitOnTrk*
>(ihit.get())->print(std::cout);
442 if (dropThisHit == 1)
break;
447 if (lGroupHits && nFound == nMost || nFound == nPrev) {
449 std::cout<<
"we've found all of the shared hits on this track,Quit"<<std::endl;
459 cout <<
"nGroup: " << nMost <<
" groupDiff: " << groupDiff << endl;
460 cout <<
"Track: " << aRecoTrk.
id() <<
" nHit: "
461 << hitList->
nHit() <<
" nActive: "
462 << tkFit->
nActive() <<
" chisq/dof: " <<
465 cout <<
"Track: "<< othTrack.
id() <<
" nHit: " <<
466 othTrack.
hits()->
nHit() <<
" nActive: " <<
472 if (groupDiff > 0.0) {
475 for (
int ii = 0; ii < nMost; ii++) {
481 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on this track, No."<<aRecoTrk.
id()<< std::endl;
484 refitTrack[trackMost] = 1;
485 for (
int ii = 0; ii < nMost; ii++) {
491 if (8 ==
tkParam.
lPrint) std::cout<<
"inactive hits on other track "<< std::endl;
504 idMap.
get(aRecoTrk.
id(), index); assert (index >= 0);
506 if (lRefit || refitTrack[index] == 1) {
508 std::cout<<
"after group ,refit track"<<aRecoTrk.
id()<< std::endl;
510 fitResult = hitList->
fit();
514 fitResult.
print(std::cerr);
522 int nDOF = tkFit->
nActive() - 5;
524 chisqperDOF = tkFit->
chisq() / nDOF;
526 chisqperDOF = tkFit->
chisq();
531 double tem2 = (float) hitList->
nHit() - tkFit->
nActive();
548 cout <<
"Refitting track " << aRecoTrk.
id() <<
" success = "
549 << fitResult.
success() <<
"\n";
552 if (fitResult.
failure() || badFit ) {
558 cout <<
"fitResult.failure? "<<fitResult.
failure()
559 <<
" badFit? "<<badFit <<
" Killing tkNo " << itrack << endl;
561 delete &(atrack->
track());
563 trksToKill.push_back(atrack);
568 if (lGroupHits)
goto restart;
571 if (8 ==
tkParam.
lPrint) std::cout<<
"end of loop over tracks"<< std::endl;
574 for (
int itk = 0; itk < (int)trksToKill.size(); itk++) {
576 if (8 ==
tkParam.
lPrint) std::cout<<
"remode dead track No."<<itk<< std::endl;
578 if (8 ==
tkParam.
lPrint) std::cout<<
"---end of arbitrateHits"<< std::endl;
580 delete [] usedInTrackNum;
581 delete [] refitTrack;
605MdcTrackListBase::transferTrack() {
ObjectVector< RecMdcHit > RecMdcHitCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
HepGeom::Point3D< double > HepPoint3D
HepGeom::Vector3D< double > HepVector3D
bool get(const K &theKey, V &theAnswer) const
void put(const K &, const V &)
void newParams(const MdcTrackParams &tkPar)
MdcTrackListBase(const MdcTrackParams &tkPar)
virtual ~MdcTrackListBase()
void store(RecMdcTrackCol *, RecMdcHitCol *)
void remove(MdcTrack *atrack)
void setTrack(TrkRecoTrk *trk)
void storeTrack(int trackId, RecMdcTrackCol *trackList, RecMdcHitCol *hitList, int tkStat)
virtual double chisq() const =0
void print(std::ostream &ostr) const
virtual void addHistory(const TrkErrCode &status, const char *modulename)
virtual int nActive() const =0
hot_iterator begin() const
bool removeHit(const TrkFundHit *theHit)
void setActivity(bool turnOn)
double resid(bool exclude=false) const
TrkRecoTrk * parentTrack() const
const TrkFundHit * hit() const
const TrkFit * fitResult() const
virtual void printAll(std::ostream &) const
const TrkFitStatus * status() const