Compare trk-cand vs. trk-fit –—————————
(II) missed tracks are defined on hits information GHOST -------------------------------------------------------------------------------------------------------------------------------------------------------————— check repeated assignment to one MC trk (maybe one MC trk rec. trk was assigned to several REC trk) ------------------—————
END check repeated assignment to one MC trk ------------------------------------------------------------------------------------------------------————
END GHOST-----------------------------------------------------------------------------------------------------------------------------------------------------------———
MISSED ------------------------------------------------------------------------------------------------------------------------------------------------------------------——
154 Double_t glXrecLMD, glYrecLMD, glZrecLMD, glThetarecLMD, glPhirecLMD;
155 Double_t glXrec, glYrec, glZrec, glThetarec, glPhirec, glMomrec;
156 Double_t glXmc, glYmc, glZmc, glThetamc, glPhimc, glMommc;
157 Double_t glXmcLMD, glYmcLMD, glZmcLMD, glThetamcLMD, glPhimcLMD, glMommcLMD;
163 int glNumDoubleMChits;
165 int glModule, glHalf;
166 FairRootManager* ioman = FairRootManager::Instance();
168 double glEvTime = ioman->GetEventTime();
176 const int nParticles =
fMCTracks->GetEntriesFast();
177 const int numTrk =
fRecTracks->GetEntriesFast();
178 const int nRecHits =
fRecHits->GetEntriesFast();
179 const int nMCHits =
fMCHits->GetEntriesFast();
185 cout <<
"%%%%%! Event #" <<
fEventNr <<
" has " << nParticles
186 <<
" true particles, " <<
" out of it " << nRecHits <<
" hits, "
187 << nTrkCandidates <<
" trk-cands, " << numTrk <<
" tracks and "
188 << nGeaneTrks <<
" geane Trks!" << endl;
193 for (Int_t iN = 0; iN < nParticles; iN++) {
198 if (motherid < 0 &&
fabs(mcID) < 1e5) {
208 int MCtksSIMhits[nParticles];
209 int MCDoubleHits[nParticles];
210 for (
int imctrk = 0; imctrk < nParticles; imctrk++) {
211 MCtksSIMhits[imctrk] = 0;
212 MCDoubleHits[imctrk] = 0;
214 for (
int imc = 0; imc < nMCHits; imc++) {
216 int MCtrk = MCPoint->GetTrackID();
217 MCtksSIMhits[MCtrk]++;
223 int MCtksREChits[nParticles];
226 for (
int imctrk = 0; imctrk < nParticles; imctrk++) {
227 MCtksREChits[imctrk] = 0;
230 cout <<
" ** ALL REChits are made from MChits: " << endl;
232 for (
int irec = 0; irec < nRecHits; irec++) {
235 int MCtrkidbot, MCtrkidtop;
238 int MCtrkid = MCPointBot->GetTrackID();
240 MCtksREChits[MCtrkid]++;
241 MCtrkidbot = MCtrkid;
243 cout <<
" " << MCtrkid <<
"(MChitID=" << mcrefbot <<
",z="
244 << MCPointBot->GetZ() <<
")";
246 int mcreftop = myHit->GetRefIndex();
249 int MCtrkid = MCPointTop->GetTrackID();
251 MCtksREChits[MCtrkid]++;
252 MCtrkidtop = MCtrkid;
254 cout <<
" " << MCtrkid <<
"(MChitID=" << mcreftop <<
",z="
255 << MCPointTop->GetZ() <<
")";
257 if (mcreftop > 0 && mcrefbot > 0) {
258 if (MCtrkidtop == MCtrkidbot) {
259 MCDoubleHits[MCtrkidbot]++;
262 cout <<
" REChit No." << irec <<
"contain MCid: " << MCtrkidtop
263 <<
", " << MCtrkidbot <<
" !" << endl;
277 for (Int_t iN = 0; iN < nGeaneTrks; iN++) {
282 glThetarecLMD = -9999;
299 glThetamcLMD = -9999;
302 trkRECStatus = -9999;
307 glNumDoubleMChits = -9999;
312 TVector3 PosRecLMD(fFittedTrkP.GetX(), fFittedTrkP.GetY(),
314 TVector3 MomRecLMD(fFittedTrkP.GetPx(), fFittedTrkP.GetPy(),
315 fFittedTrkP.GetPz());
316 MomRecLMD *=
fPbeam / MomRecLMD.Mag();
322 const int numPts = trkcand->
GetNHits();
331 TVector3 MomRecCandLMD = pos2 - posSeed;
332 MomRecCandLMD *=
fPbeam / MomRecCandLMD.Mag();
341 int MCtrk[nParticles];
342 int RECtrkMCid[nGeaneTrks];
343 for (
int nk = 0;
nk < nParticles;
nk++)
345 bool goodTrk[nGeaneTrks];
346 bool ghostTrk[nGeaneTrks];
347 for (Int_t iN = 0; iN < nGeaneTrks; iN++) {
349 ghostTrk[iN] =
false;
354 for (Int_t iN = 0; iN < nGeaneTrks; iN++) {
355 FairTrackParH *fRes = (FairTrackParH*)
fRecBPTracks->At(iN);
356 TVector3 PosRec = fRes->GetPosition();
357 Double_t lyambda = fRes->GetLambda();
359 cout <<
"GEANE didn't propagate trk No." << iN <<
"!" << endl;
367 const int Ntrkcandhits = trkcand->
GetNHits();
369 int MCtrkID[Ntrkcandhits];
370 for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++) {
371 MCtrkID[iHit] = 9999;
378 cout <<
" *** REChits in trk (with " << Ntrkcandhits
379 <<
" hits) are made from MChits: " << endl;
380 for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++) {
382 cout <<
"No." << iHit <<
": ";
391 int MCtrkid = MCPointBot->GetTrackID();
395 cout <<
" " << MCtrkid <<
"!!!";
399 MCtrkID[iHit] = MCtrkid;
401 cout <<
" " << MCtrkid;
405 cout <<
" Ooops, mcrefbot = " << mcrefbot;
407 int mcreftop = myHit->GetRefIndex();
413 int MCtrkid = MCPointTop->GetTrackID();
417 cout <<
" " << MCtrkid <<
"!!!";
421 MCtrkID[iHit] = MCtrkid;
423 cout <<
" " << MCtrkid;
427 cout <<
" Ooops, mcreftop = " << mcreftop;
440 for (Int_t
n = 0;
n < Ntrkcandhits;
n++) {
443 for (Int_t
m =
n + 1;
m < Ntrkcandhits;
m++)
444 if (MCtrkID[
m] < x) {
450 MCtrkID[k] = MCtrkID[
n];
456 Int_t prevID = MCtrkID[0];
457 for (
int nk = 0;
nk < nParticles;
nk++) {
458 if (MCtrkID[0] ==
nk)
459 MCtrk[
nk] = MCtrk[
nk] + 1;
462 for (Int_t
n = 1;
n < Ntrkcandhits;
n++) {
463 if (MCtrkID[
n] > 9998)
465 if (prevID < MCtrkID[
n]) {
468 for (
int nk = 0;
nk < nParticles;
nk++) {
469 if (MCtrkID[n] ==
nk)
470 MCtrk[
nk] = MCtrk[
nk] + 1;
479 RECtrkMCid[iN] = MCtrkID[0];
481 if (MCtrkID[0] < 9999)
488 vector<int> countMC_IDs(diffIDs);
489 for (
int kn = 0; kn < diffIDs; kn++) {
495 for (Int_t n = 0; n < Ntrkcandhits; n++) {
497 if (MCtrkID[n] > 9998) {
503 if (prevID < MCtrkID[n]) {
507 countMC_IDs[diffCount]++;
510 int maxID = countMC_IDs[0];
512 for (
int kn = 0; kn < diffIDs; kn++) {
514 cout <<
"countMC_IDs[" << kn <<
"]=" << countMC_IDs[kn]
515 <<
" 0.65*Ntrkcandhits = " << 0.65 * Ntrkcandhits << endl;
516 if (countMC_IDs[kn] > 0.65 * Ntrkcandhits) {
518 ghostTrk[iN] =
false;
526 if (countMC_IDs[kn] > maxID) {
527 maxID = countMC_IDs[kn];
533 for (Int_t n = 0; n < Ntrkcandhits; n++) {
534 if (diffCount == posIDmax)
535 RECtrkMCid[iN] = prevID;
536 if (prevID < MCtrkID[n]) {
590 int RECassigMC[nGeaneTrks];
591 for (Int_t iN = 0; iN < nGeaneTrks; iN++)
594 for (Int_t iN = 0; iN < nGeaneTrks; iN++) {
595 int mc_i = RECtrkMCid[iN];
596 for (Int_t jN = (iN + 1); jN < nGeaneTrks; jN++) {
597 int mc_j = RECtrkMCid[jN];
601 cout <<
"GHOST?!: rec.trks #" << iN <<
" and #" << jN
602 <<
"were assigned to one MC trk =|" << endl;
608 int goodRecII = 0, ghostRecI = 0, ghostRecII = 0;
609 for (Int_t iN = 0; iN < nGeaneTrks; iN++) {
612 if (RECassigMC[iN] > 1) {
621 cout <<
"... RECtrk#" << iN <<
" was defined as GOOD ..." << endl;
627 cout <<
"... RECtrk#" << iN <<
" was defined as GHOST ..." << endl;
634 cout <<
"Ooops, RECtrk isn't GOOD and isn't GHOST!!!" << endl;
638 FairTrackParH *fRes = (FairTrackParH*)
fRecBPTracks->At(iN);
639 Double_t lyambda = fRes->GetLambda();
645 TVector3 MomRecBP = fRes->GetMomentum();
646 TVector3 PosBP = fRes->GetPosition();
652 glThetarec = thetaBP;
654 glMomrec = MomRecBP.Mag();
662 TVector3 errMomBP(errPx, errPy, errPz);
664 Double_t err_lyambda = fRes->GetDLambda();
671 TVector3 PosRecLMD(fFittedTrkP.GetX(), fFittedTrkP.GetY(),
673 TVector3 MomRecLMD(fFittedTrkP.GetPx(), fFittedTrkP.GetPy(),
674 fFittedTrkP.GetPz());
675 MomRecLMD *=
fPbeam / MomRecLMD.Mag();
676 glXrecLMD = PosRecLMD.X();
677 glYrecLMD = PosRecLMD.Y();
678 glZrecLMD = PosRecLMD.Z();
679 glThetarecLMD = MomRecLMD.Theta();
680 glPhirecLMD = MomRecLMD.Phi();
681 trkRECStatus = trkType;
690 auto const& digi_info(lmd_geo_helper.getHitLocationInfo(sensorID));
691 glModule = digi_info.module;
692 glHalf = digi_info.detector_half;
693 glTrkTime = myHit->GetTimeStamp();
708 int MCidforREC = RECtrkMCid[iN];
711 if (MCidforREC > 9998) {
721 glThetamcLMD = -9999;
725 glNumDoubleMChits = -9999;
733 glThetamc = MomMC.Theta();
734 glPhimc = MomMC.Phi();
735 glMommc = MomMC.Mag();
745 glNumMChits = MCtksSIMhits[MCidforREC];
746 glNumDoubleMChits = MCDoubleHits[MCidforREC];
751 int mcreftop = myHit->GetRefIndex();
753 cout <<
"mcrefbot = " << mcrefbot <<
" mcreftop = " << mcreftop << endl;
764 if (MCPointHit != NULL) {
767 double pxTrue = MCPointHit->GetPx();
768 double pyTrue = MCPointHit->GetPy();
769 double pzTrue = MCPointHit->GetPz();
770 TVector3 MomMClmd(pxTrue, pyTrue, pzTrue);
771 glXmcLMD = PosMClmd.X();
772 glYmcLMD = PosMClmd.Y();
773 glZmcLMD = PosMClmd.Z();
774 glThetamcLMD = MomMClmd.Theta();
775 glPhimcLMD = MomMClmd.Phi();
776 glMommcLMD = MomMClmd.Mag();
789 glThetamcLMD = -9999;
793 glNumDoubleMChits = -9999;
805 trkqlmd->
SetLMDpoint(glXrecLMD, glYrecLMD, glZrecLMD);
806 trkqlmd->
SetLMDdir(glThetarecLMD, glPhirecLMD);
809 trkqlmd->
SetIPmom(glThetarec, glPhirec, glMomrec);
811 trkqlmd->
SetIPerrmom(err_lyambda, err_phi, errMomBP.Mag());
814 trkqlmd->
SetMCmom(glThetamc, glPhimc, glMommc);
819 trkqlmd->
SetMCmomLMD(glThetamcLMD, glPhimcLMD, glMommcLMD);
834 if ((nParticles - goodRecII) > 0) {
836 int nMCmissedTrkSearch = 0;
837 int nMCmissedLossHits = 0;
838 for (
int imc = 0; imc < nParticles; imc++) {
840 for (
int irec = 0; irec < nGeaneTrks; irec++) {
841 int mc_comp = RECtrkMCid[irec];
842 if (mc_comp == imc && goodTrk[irec])
852 glThetamc = MomMC.Theta();
853 glPhimc = MomMC.Phi();
854 glMommc = MomMC.Mag();
858 int minHits = MCDoubleHits[imc] + 2;
859 if (MCDoubleHits[imc] == 4)
863 if (MCtksREChits[imc] > minHits) {
865 cout <<
" --- MCtrk#" << imc
866 <<
" was defined as MISSED during trk-search (#MChits="
867 << MCtksREChits[imc] <<
" with limit>" << minHits <<
")"
869 if (MCDoubleHits[imc] > 0)
870 cout <<
" NB:this MCtrk contains " << MCDoubleHits[imc]
871 <<
" double hits!" << endl;
874 nMCmissedTrkSearch++;
876 if (MCtksREChits[imc] > 0) {
881 cout <<
" --- MCtrk#" << imc
882 <<
" was defined as MISSED due to little amount of hits (#MChits="
883 << MCtksREChits[imc] <<
" with limit>" << minHits <<
")"
887 for (
int imhc = 0; imhc < nMCHits; imhc++) {
889 int MCtrkID = MCPoint->GetTrackID();
891 if (MCtrkID == imc) {
899 double pxTrue = MCPointHit->GetPx();
900 double pyTrue = MCPointHit->GetPy();
901 double pzTrue = MCPointHit->GetPz();
902 TVector3 MomMClmd(pxTrue, pyTrue, pzTrue);
903 glXmcLMD = PosMClmd.X();
904 glYmcLMD = PosMClmd.Y();
905 glZmcLMD = PosMClmd.Z();
906 glThetamcLMD = MomMClmd.Theta();
907 glPhimcLMD = MomMClmd.Phi();
908 glMommcLMD = MomMClmd.Mag();
914 glThetamcLMD = -9999;
923 glThetamcLMD = -9999;
939 glThetarecLMD = -9999;
954 TClonesArray& clref = *
fTrackQ;
955 Int_t size = clref.GetEntriesFast();
960 trkqlmd->
SetLMDpoint(glXrecLMD, glYrecLMD, glZrecLMD);
961 trkqlmd->
SetLMDdir(glThetarecLMD, glPhirecLMD);
964 trkqlmd->
SetIPmom(glThetarec, glPhirec, glMomrec);
966 trkqlmd->
SetMCmom(glThetamc, glPhimc, glMommc);
974 trkqlmd->
SetMCmomLMD(glThetamcLMD, glPhimcLMD, glMommcLMD);
991 cout <<
"PndLmdTrkQTask::Exec END!" << endl;
TClonesArray * fRecBPTracks
void SetEvRECMulti(int tot)
TVector3 GetPosition() const
void SetMCmom(double theta, double phi, double mom)
void SetMCpoint(double x, double y, double z)
PndTrackCandHit GetSortedHit(UInt_t i)
void SetTrkRecStatus(int st)
static PndLmdGeometryHelper & getInstance()
void SetLMDdir(double theta, double phi)
TVector3 GetMomentum() const
void SetIPpoint(double x, double y, double z)
void SetMCpointLMD(double x, double y, double z)
void SetIPerrpoint(double errx, double erry, double errz)
void SetMCmomLMD(double theta, double phi, double mom)
void SetSecondary(int sec)
void SetIPerrmom(double errtheta, double errphi, double errmom)
TClonesArray * fRecCandTracks
friend F32vec4 fabs(const F32vec4 &a)
void SetEvTime(double evtm)
void SetTrkTime(double trktm)
TClonesArray * fRecTracks
void SetSumEvPDG(int sumid)
void SetLMDchi2(double chi2)
TVector3 GetPosition() const
void SetLMDpoint(double x, double y, double z)
Int_t GetSecondMCHit() const
void SetNumMChits(int num)
void SetNumDoubleMChits(int num)
void SetEvMCMulti(int tot)
Int_t GetRefIndex() const
TVector3 GetStartVertex() const
Int_t GetSensorID() const
Int_t GetMotherID() const
void SetIPmom(double theta, double phi, double mom)
FairTrackParP GetParamFirst()