From 509201a15cd7cda9ad526fd12fa8563b0965d027 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Fri, 13 Jan 2023 17:19:09 +0000
Subject: [PATCH] fixes in CbmKFTrackQa
---
reco/KF/KFQA/CbmKFTrackQa.cxx | 104 ++++++++++++++++++++++------------
reco/KF/KFQA/CbmKFTrackQa.h | 5 ++
2 files changed, 72 insertions(+), 37 deletions(-)
diff --git a/reco/KF/KFQA/CbmKFTrackQa.cxx b/reco/KF/KFQA/CbmKFTrackQa.cxx
index 79a4ef8e83..e8cd23e981 100644
--- a/reco/KF/KFQA/CbmKFTrackQa.cxx
+++ b/reco/KF/KFQA/CbmKFTrackQa.cxx
@@ -22,6 +22,7 @@
#include "KFParticleMatch.h"
#include "KFParticleTopoReconstructor.h"
#include "KFTopoPerformance.h"
+#include "L1Field.h"
//ROOT headers
#include "TClonesArray.h"
@@ -37,6 +38,7 @@
#include "CbmGlobalTrack.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
+#include "CbmMCEventList.h"
#include "CbmMuchTrack.h"
#include "CbmRichRing.h"
#include "CbmStsTrack.h"
@@ -56,6 +58,7 @@ using std::vector;
CbmKFTrackQa::CbmKFTrackQa(const char* name, Int_t iVerbose, TString outFileName)
: FairTask(name, iVerbose)
+ , fMcEventListBranchName("MCEventList.")
, fStsTrackBranchName("StsTrack")
, fGlobalTrackBranchName("GlobalTrack")
, fRichBranchName("RichRing")
@@ -272,6 +275,13 @@ InitStatus CbmKFTrackQa::Init()
return kERROR;
}
+ // Get mc event list
+ fMcEventList = static_cast<CbmMCEventList*>(ioman->GetObject("fMcEventListBranchName"));
+ if (!fMcEventList) {
+ Warning("CbmKFTrackQa::Init", "mc event list not found!");
+ return kERROR;
+ }
+
// Get sts tracks
fStsTrackArray = (TClonesArray*) ioman->GetObject(fStsTrackBranchName);
if (fStsTrackArray == 0) {
@@ -512,25 +522,30 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/)
Int_t stsTrackIndex = globalTrack->GetStsTrackIndex(); //for STS histos
CbmStsTrack* cbmStsTrack = (CbmStsTrack*) fStsTrackArray->At(stsTrackIndex);
- int stsTrackMCIndex = trackMatch[stsTrackIndex];
- Double_t* stsHistoData = new Double_t[NStsHisto];
- stsHistoData[0] = cbmStsTrack->GetNofHits(); //NHits
- stsHistoData[1] = cbmStsTrack->GetChiSq() / cbmStsTrack->GetNDF(); //Chi2/NDF
- stsHistoData[2] = TMath::Prob(cbmStsTrack->GetChiSq(), cbmStsTrack->GetNDF()); //prob
- if (stsTrackMCIndex > -1) {
+ double stsHistoData[NStsHisto] = {
+ (double) cbmStsTrack->GetNofHits(), //NHits
+ cbmStsTrack->GetChiSq() / cbmStsTrack->GetNDF(), //Chi2/NDF
+ TMath::Prob(cbmStsTrack->GetChiSq(), cbmStsTrack->GetNDF()) //prob
+ };
+
+ for (int iH = 0; iH < NStsHisto; iH++) {
+ hStsHisto[0][iH]->Fill(stsHistoData[iH]);
+ }
+
+ int stsTrackMCIndex = trackMatch[stsTrackIndex];
+ if (stsTrackMCIndex >= 0) {
int iDir = GetHistoIndex(mcTracks[stsTrackMCIndex].PDG());
if (iDir < 3) {
for (int iH = 0; iH < NStsHisto; iH++) {
hStsHisto[iDir][iH]->Fill(stsHistoData[iH]);
- hStsHisto[0][iH]->Fill(stsHistoData[iH]);
}
}
}
- else
- for (int iH = 0; iH < NStsHisto; iH++) //ghost
+ else { //ghost
+ for (int iH = 0; iH < NStsHisto; iH++) {
hStsHisto[7][iH]->Fill(stsHistoData[iH]);
- delete[] stsHistoData;
-
+ }
+ }
Int_t muchIndex = globalTrack->GetMuchTrackIndex(); //for MuCh histos
if (muchIndex > -1) {
@@ -692,37 +707,52 @@ void CbmKFTrackQa::Exec(Option_t* /*opt*/)
if (fTofHitArray && fTofHitMatchArray) {
Int_t tofIndex = globalTrack->GetTofHitIndex(); // for tof histo
Int_t stsIndex = globalTrack->GetStsTrackIndex();
+ // std::cout << "tofIndex: " << tofIndex << " stsIndex: " << stsIndex << std::endl;
if (tofIndex > -1 && stsIndex > -1) {
CbmTofHit* tofHit = (CbmTofHit*) fTofHitArray->At(tofIndex);
CbmMatch* tofHitMatch = (CbmMatch*) fTofHitMatchArray->At(tofIndex);
CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatchArray->At(stsIndex);
-
- Double_t l = globalTrack->GetLength();
- Double_t time = tofHit->GetTime();
- Double_t q = stsPar->GetQp() > 0 ? 1. : -1.;
- Double_t m2 = p * p * (1. / ((l / time / 29.9792458) * (l / time / 29.9792458)) - 1.);
-
- hTofHisto2D[0][0]->Fill(p * q, m2);
- hTofHisto2D[0][1]->Fill(eloss * 1e6, m2);
-
- if (tofHit && tofHitMatch && stsMatch) {
- Int_t tofMCPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
- Int_t tofMCEventNo = tofHitMatch->GetMatchedLink().GetEntry();
- CbmTofPoint* tofPoint = (CbmTofPoint*) fTofPoints->Get(0, tofMCEventNo, tofMCPointIndex);
- Int_t tofMCTrackId = tofPoint->GetTrackID();
-
- int iHitCategory = -1;
- if (stsTrackMCIndex < 0) iHitCategory = 8; // ghost sts track + any tof hit
- else if (tofMCTrackId < 0)
- iHitCategory = 9; // normal sts track + ghost tof hit
- else if (stsTrackMCIndex != tofMCTrackId)
- iHitCategory = 7; // mismatched sts track and tof hit
- else
- iHitCategory = GetHistoIndex(pdg[stsTrackIndex]);
-
- hTofHisto2D[iHitCategory][0]->Fill(p * q, m2);
- hTofHisto2D[iHitCategory][1]->Fill(eloss * 1e6, m2);
+ CbmLink tofLink = tofHitMatch->GetMatchedLink();
+
+ if (tofLink.GetFile() >= 0 && tofLink.GetEntry() >= 0) { // the mc link is needed for the event time
+
+ double eventTime = fMcEventList->GetEventTime(tofLink);
+
+ Double_t l = globalTrack->GetLength();
+ Double_t time = tofHit->GetTime() - eventTime;
+ Double_t q = stsPar->GetQp() > 0 ? 1. : -1.;
+ // Double_t beta = l / time / 29.9792458;
+ // std::cout << " l: " << l << " time: " << time << " beta: " << beta << std::endl;
+ Double_t m2 = p * p * (1. / ((l / time / 29.9792458) * (l / time / 29.9792458)) - 1.);
+ // std::cout << "p: " << p << " m2: " << m2 << std::endl;
+
+ hTofHisto2D[0][0]->Fill(p * q, m2);
+ hTofHisto2D[0][1]->Fill(eloss * 1e6, m2);
+
+ if (tofHit && tofHitMatch && stsMatch) {
+ CbmTofPoint* tofPoint = (CbmTofPoint*) fTofPoints->Get(tofLink);
+ Int_t tofMCTrackId = tofPoint->GetTrackID();
+ // Double_t tofMCTime = tofPoint->GetTime();
+ // Double_t tofMCl = tofPoint->GetLength();
+ // Double_t MCm2 = p * p * (1. / (( tofMCl/ tofMCTime / 29.9792458) * (tofMCl / tofMCTime / 29.9792458)) - 1.);
+ // hTofHisto2D[0][0]->Fill(p * q, MCm2);
+ // hTofHisto2D[0][1]->Fill(eloss * 1e6, MCm2);
+
+ // std::cout << " mcl: " << tofMCl << " mc time: " << tofMCTime << " mc beta: " << tofMCl/tofMCTime/29.9792458 << std::endl;
+
+ int iHitCategory = -1;
+ if (stsTrackMCIndex < 0) iHitCategory = 8; // ghost sts track + any tof hit
+ else if (tofMCTrackId < 0)
+ iHitCategory = 9; // normal sts track + ghost tof hit
+ else if (stsTrackMCIndex != tofMCTrackId)
+ iHitCategory = 7; // mismatched sts track and tof hit
+ else
+ iHitCategory = GetHistoIndex(pdg[stsTrackIndex]);
+
+ hTofHisto2D[iHitCategory][0]->Fill(p * q, m2);
+ hTofHisto2D[iHitCategory][1]->Fill(eloss * 1e6, m2);
+ }
}
}
}
diff --git a/reco/KF/KFQA/CbmKFTrackQa.h b/reco/KF/KFQA/CbmKFTrackQa.h
index 70bd360bff..03353ef419 100644
--- a/reco/KF/KFQA/CbmKFTrackQa.h
+++ b/reco/KF/KFQA/CbmKFTrackQa.h
@@ -24,6 +24,7 @@ class TH1F;
class TH2F;
class TObject;
class CbmMCDataArray;
+class CbmMCEventList;
class CbmKFTrackQa : public FairTask {
public:
@@ -55,6 +56,8 @@ private:
int GetHistoIndex(int pdg);
//names of input branches
+
+ TString fMcEventListBranchName;
TString fStsTrackBranchName;
TString fGlobalTrackBranchName;
TString fRichBranchName;
@@ -70,6 +73,8 @@ private:
TString fMuchTrackMatchBranchName;
//input branches
+
+ CbmMCEventList* fMcEventList;
TClonesArray* fStsTrackArray;
TClonesArray* fGlobalTrackArray;
TClonesArray* fRichRingArray;
--
GitLab