18 #include "TClonesArray.h"
22 #include "FairRootManager.h"
24 #include "FairRuntimeDb.h"
26 #include "FairMultiLinkedData.h"
40 FairTask(
"Creates PndMC test"), fMCInfoBranchName(
"MCTrackInfo"), fRecoInfoBranchName(
"RecoTrackInfo"), fTrackBranchName(trackBranchName), fIdealTrackBranchName(idealBranchName), fPndTrackOrTrackCand(pndTrackData), fEventNr(0) {
53 ioman = FairRootManager::Instance();
55 std::cout <<
"-E- PndTrackingQATask::Init: "
56 <<
"RootManager not instantiated!" << std::endl;
67 std::cout <<
"-E- PndTrackingQATask::Init " <<
"no track branch " <<
fTrackBranchName << std::endl;
71 std::cout <<
"-E- PndTrackingQATask::Init " <<
"no MC track branch " << std::endl;
76 std::cout <<
"-E- PndTrackingQATask::Init " <<
"no ideal track branch " <<
fIdealTrackBranchName << std::endl;
81 fMCTrackInfo =
new TClonesArray(
"PndTrackingQualityMCInfo");
122 fIdealTracksPerEvent =
new TH1I(
"fIdealTracksPerEvent",
"Ideal Tracks per Event", 1000, -0.5,999.5);
124 fIdealPHisto =
new TH1I(
"fIdealPHisto",
"Ideal total Momentum", 1500, -0.5,14.5);
126 fIdealPtHisto =
new TH1I(
"fIdealPtHisto",
"Ideal Tansversal Momentum", 1500, -0.5,14.5);
128 fIdealPlHisto =
new TH1I(
"fIdealPlHisto",
"Ideal Longitudinal Momentum", 1500, -0.5,14.5);
130 fPHisto =
new TH1D(
"fPHisto",
"Momentum Resolution", 1000, -1, 1);
131 fPHisto->GetXaxis()->SetTitle(
"p^{RECO} - p^{MC} / GeV");
132 fPHisto->GetYaxis()->SetTitle(
"counts");
133 fPRelHisto =
new TH1D(
"fPRelHisto",
"Relative Momentum Resolution", 1000, -1, 1);
134 fPRelHisto->GetXaxis()->SetTitle(
"(p^{RECO} - p^{MC}) / p^{MC}");
136 fPtHisto =
new TH1D(
"fPtHisto",
"Transverse Momentum Resolution", 1000, -1, 1);
137 fPtHisto->GetXaxis()->SetTitle(
"p_{t}^{RECO} - p_{t}^{MC} / GeV");
138 fPtHisto->GetYaxis()->SetTitle(
"counts");
139 fPtRelHisto =
new TH1D(
"fPtRelHisto",
"Relative Transverse Momentum Resolution", 1000, -1, 1);
140 fPtRelHisto->GetXaxis()->SetTitle(
"(p_{t}^{RECO} - p_{t}^{MC}) / p_{t}^{MC}");
142 fPlHisto =
new TH1D(
"fPlHisto",
"Longitudinal Momentum Resolution", 1000, -1, 1);
143 fPlHisto->GetXaxis()->SetTitle(
"p_{l}^{RECO} - p_{l}^{MC} / GeV");
144 fPlHisto->GetYaxis()->SetTitle(
"counts");
145 fPlRelHisto =
new TH1D(
"fPlRelHisto",
"Relative Longitudinal Momentum Resolution", 1000, -1, 1);
146 fPlRelHisto->GetXaxis()->SetTitle(
"(p_{l}^{RECO} - p_{l}^{MC}) / p_{l}^{MC}");
149 fQualyHisto =
new TH1I(
"fQualyHisto",
"Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
153 fQualyHisto_all =
new TH1I(
"fQualyHisto_all",
"Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
154 fQualyHisto_pos =
new TH1I(
"fQualyHisto_pos",
"Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
155 fQualyHisto_neg =
new TH1I(
"fQualyHisto_neg",
"Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
156 fQualyHisto_mc =
new TH1I(
"fQualyHisto_mc",
"Quality of Trackfinding;;Counts", 26, -15.5, 10.5);
159 fQualyHisto_rel_all =
new TH1D(
"fQualyHisto_rel_all",
"Quality of Trackfinding;;Relative", 26, -15.5, 10.5);
164 fQualyHisto_rel_possible =
new TH1D(
"fQualyHisto_rel_possible",
"Quality of Trackfinding;;Relative", 26, -15.5, 10.5);
196 FairRuntimeDb*
rtdb = FairRun::Instance()->GetRuntimeDb();
211 std::cout <<
"----- Event " <<
fEventNr <<
" ------" << std::endl;
223 std::map<Int_t, TVector3> recoPMap = qaAna.
GetP();
246 for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++)
248 Int_t mcTrackId = iter->first;
251 if (idealTrackId < 0){
252 std::cout <<
"-W- PndTrackingQATask::Exec no idealTrack for mcTrack " << mcTrackId << std::endl;
255 Int_t trackQuality = iter->second;
259 if (idealtrack == 0){
260 std::cout <<
"-E- No ideal track found for idealTrackId " << idealTrackId << std::endl;
263 if (mcTrackId == -1){
264 std::cout <<
"-W- PndTrackingQATask::Exec mcTrackId == -1" << std::endl;
271 std::cout <<
"-E- PndTrackingQATask::Exec mcMyTrack == 0" << std::endl;
294 for(
int itrk = 0; itrk <
fRecoTrackInfo->GetEntriesFast(); itrk++) {
298 if (idealTrackId < 0){
299 std::cout <<
"-W- PndTrackingQATask::Exec no idealTrack for recoTrack " << itrk << std::endl;
310 for (std::map<Int_t, Int_t>::iterator iter = qualiMap.begin(); iter != qualiMap.end(); iter++) {
311 Int_t mcTrackId = iter->first;
312 Int_t trackQuality = iter->second;
314 if (mcTrackId == -1){
315 std::cout <<
"-W- PndTrackingQATask::Exec mcTrackId == -1" << std::endl;
319 TVector3 recoMomentum = recoPMap[mcTrackId];
324 std::cout <<
"-E- PndTrackingQATask::Exec mcMyTrack == 0" << std::endl;
333 fTuple->
Column(
"McTrackFoundNTimes", (Int_t) mcFoundMap[mcTrackId]);
335 fTuple->
Column(
"Mc_TrackQuality", (Int_t) mcStatusMap[mcTrackId] - 6);
359 for (
size_t branchIndex = 0; branchIndex <
fBranchNames.size(); branchIndex++){
360 result += trackData->GetLinksWithType(
ioman->GetBranchId(
fBranchNames[branchIndex])).GetNLinks();
370 for(std::map<Int_t, Int_t>::iterator iter = trackQualifikation.begin(); iter != trackQualifikation.end(); iter++){
372 if (iter->second > 0){
384 for(std::map<Int_t, Int_t>::iterator iter = trackMCStatus.begin(); iter != trackMCStatus.end(); iter++){
392 for (
std::map<Int_t,
std::map<
TString, std::pair<Double_t, Int_t> > >::iterator iterTracks = efficiencies.begin(); iterTracks != efficiencies.end(); iterTracks++ ){
393 std::map<TString, std::pair<Double_t, Int_t> > branchEfficiency = iterTracks->second;
394 for (
std::map<
TString, std::pair<Double_t, Int_t> >::iterator iterBranch = branchEfficiency.begin(); iterBranch != branchEfficiency.end(); iterBranch++ ){
395 fMapEfficiencies[iterBranch->first]->Fill(iterBranch->second.second, iterBranch->second.first);
400 for (std::map<Int_t, Double_t>::iterator iter = map.begin(); iter != map.end(); iter++ ){
401 histo->Fill(iter->second);
408 Int_t allTracksWithHits = 0;
409 Int_t allPossibleTracksWithHits = 0;
410 Int_t allTracksWithHitsNotFound = 0;
427 Double_t allFound = fullyFound + partiallyFound + spuriousFound;
431 allTracksWithHits += mcAtLeastThreePrim;
432 allTracksWithHits += mcAtLeastThreeSec;
433 allTracksWithHits += mcPossiblePrim;
434 allTracksWithHits += mcPossibleSec;
436 allTracks = allTracksWithHits + lessThanThreePrim;
438 allPossibleTracksWithHits += mcPossiblePrim;
439 allPossibleTracksWithHits += mcPossibleSec;
441 allTracksWithHitsNotFound += atLeastThreePrim;
442 allTracksWithHitsNotFound += atLeastThreeSec;
443 allTracksWithHitsNotFound += possiblePrim;
444 allTracksWithHitsNotFound += possibleSec;
447 if (relative == kTRUE){
448 divisor = allTracks / 100.0;
460 if (mcLessThanThreePrim > 0){
461 if (relative == kTRUE){
462 divisor = mcLessThanThreePrim / 100.0;
469 if (mcAtLeastThreePrim > 0) {
470 if (relative == kTRUE){
471 divisor = mcAtLeastThreePrim / 100.0;
478 if (mcAtLeastThreeSec > 0){
479 if (relative == kTRUE){
480 divisor = mcAtLeastThreeSec / 100.0;
487 if (mcPossiblePrim > 0){
488 if (relative == kTRUE){
489 divisor = mcPossiblePrim / 100.0;
496 if (mcPossibleSec > 0){
497 if (relative == kTRUE){
498 divisor = mcPossibleSec / 100.0;
504 if (relative == kTRUE){
550 Int_t allTracksWithHits = 0;
551 Int_t allPossibleTracksWithHits = 0;
552 Int_t allTracksWithHitsNotFound = 0;
566 allTracksWithHits += mcAtLeastThreePrim;
567 allTracksWithHits += mcAtLeastThreeSec;
568 allTracksWithHits += mcPossiblePrim;
569 allTracksWithHits += mcPossibleSec;
571 allTracks = allTracksWithHits + lessThanThreePrim;
573 allPossibleTracksWithHits += mcPossiblePrim;
574 allPossibleTracksWithHits += mcPossibleSec;
576 allTracksWithHitsNotFound += atLeastThreePrim;
577 allTracksWithHitsNotFound += atLeastThreeSec;
578 allTracksWithHitsNotFound += possiblePrim;
579 allTracksWithHitsNotFound += possibleSec;
589 std::cout <<
"fQualyHisto: All Tracks: " << allTracks << std::endl
590 <<
" Primary Tracks < 3 hits: " << mcLessThanThreePrim <<
": Not Found: " << lessThanThreePrim << std::endl
591 <<
" All Tracks with hits: " << allTracksWithHits <<
". Not Found: " << allTracksWithHitsNotFound << std::endl
592 <<
" Primary Tracks with >= 3 hits, but not a possible track: " << mcAtLeastThreePrim <<
". Not Found: " << atLeastThreePrim << std::endl
593 <<
" Secondary Tracks with >= 3 hits, but not a possible track: " << mcAtLeastThreeSec <<
". Not Found: " << atLeastThreeSec << std::endl
594 <<
" Primary Tracks possible: " << mcPossiblePrim <<
". Not Found: " << possiblePrim << std::endl
595 <<
" Secondary Tracks possible: " << mcPossibleSec <<
". Not Found: " << possibleSec << std::endl
596 <<
" All Possible Tracks with hits: " << allPossibleTracksWithHits << std::endl
598 <<
" FullyFound: " << fullyFound <<
" "
599 << fullyFound / allTracksWithHits * 100.0 <<
"% (of all tracks), "
600 << fullyFound / allPossibleTracksWithHits * 100.0 <<
"% (of all tracks possible)"
602 <<
" PartlyFound: " << partiallyFound <<
" "
603 << partiallyFound / allTracksWithHits * 100.0 <<
"% "
604 << partiallyFound / allPossibleTracksWithHits * 100.0 <<
"% "
606 <<
" Spurious: " << spuriousFound <<
" "
607 << spuriousFound / allTracksWithHits * 100.0 <<
"% "
608 << spuriousFound / allPossibleTracksWithHits * 100.0 <<
"% "
610 <<
" Ghosts: " << ghosts <<
" "
611 << ghosts / allTracksWithHits * 100.0 <<
"% "
612 << ghosts / allPossibleTracksWithHits * 100.0 <<
"% " << std::endl;
629 std::cout <<
"Finish finished!" << std::endl;
665 Int_t nofmvdpixpoint = idealtrkcand->
GetNHitsDet(FairRootManager::Instance()->GetBranchId(
"MVDHitsPixel"));
666 Int_t nofmvdstrpoint = idealtrkcand->
GetNHitsDet(FairRootManager::Instance()->GetBranchId(
"MVDHitsStrip"));
668 Int_t nofgempoint = idealtrkcand->
GetNHitsDet(FairRootManager::Instance()->GetBranchId(
"GEMHit"));
669 Int_t nofftspoint = idealtrkcand->
GetNHitsDet(FairRootManager::Instance()->GetBranchId(
"FTSHit"));
671 int nofsttskewpoint = 0, nofsttparalpoint = 0;
673 for(
size_t ihit = 0; ihit < idealtrkcand->
GetNHits(); ihit++) {
675 Int_t hitID = idealcandhit.
GetHitId();
676 Int_t detID = idealcandhit.
GetDetId();
678 if(detID != FairRootManager::Instance()->GetBranchId(
"STTHit"))
continue;
682 if(tube->
IsSkew()) nofsttskewpoint++;
683 else nofsttparalpoint++;
686 PndTrackingQualityMCInfo info(nofmvdpixpoint, nofmvdstrpoint, nofsttparalpoint, nofsttskewpoint, nofgempoint, nofftspoint);
687 std::vector<FairLink> mcTracks = idealtrack->GetSortedMCTracks();
688 if (mcTracks.size() > 0){
689 PndMCTrack* myMCTrack = (
PndMCTrack*)FairRootManager::Instance()->GetCloneOfLinkData(mcTracks[0]);
690 if (myMCTrack !=
nullptr){
888 for(
int imctrk = 0; imctrk <
fMCTrackInfo->GetEntriesFast(); imctrk++) {
892 double tmpeff = 0., tmppur = 0.;
893 int tmptruerecotrackid = -1;
897 for(
int itrk = 0; itrk <
fRecoTrackInfo->GetEntriesFast(); itrk++) {
900 if(mctrackid != mctrackid0)
continue;
911 tmprecoinfo = recoinfo;
913 if(tmprecoinfo == NULL)
continue;
virtual void FillQualyHisto(std::map< Int_t, Int_t > trackQualifikation, Int_t nGhosts)
TClonesArray * fRecoTrackInfo
UInt_t GetNHitsDet(Int_t detId)
std::map< Int_t, Int_t > GetMCTrackFound()
void SetAssoRecoTrackID(int asso)
void InitializeHistograms()
void LabelQualyHistogram(TH1 *)
static const int kMcAtLeastThreeSec
std::map< Int_t, Double_t > GetPlResolution()
void AnalyseEvent(TClonesArray *recoTrackInfo)
void SetMCTrackID(Int_t mctrackid)
TClonesArray * fMCTrackInfo
virtual void MapToHist(std::map< Int_t, Double_t >, TH1 *)
void SetMomentumLast(TVector3 mom)
static const int kLessThanThreePrim
void SetHitsBranchNames(std::vector< TString > names)
PndTrackCandHit GetSortedHit(UInt_t i)
static const int kPossiblePrim
static const int kAtLeastThreeSec
PndTrackingQATask(TString trackBranchName, TString idealBranchName, Bool_t pndTrackData=kTRUE)
std::map< Int_t, Double_t > GetPResolution()
std::map< Int_t, Double_t > GetPResolutionRel()
TVector3 GetMomentum() const
void SetRecoTrackID(int recotrkid)
TClonesArray * fIdealTrack
static const int kPartiallyFound
static const int kMcPossibleSec
void SetMCTrackInfo(PndTrackingQualityMCInfo *info)
TClonesArray * fSttTubeArray
virtual void FillEfficiencies(std::map< Int_t, std::map< TString, std::pair< Double_t, Int_t > > > efficiencies)
Bool_t fPndTrackOrTrackCand
TH1 * fQualyHisto_rel_possible
virtual void SetQualyHisto(TH1 *histo, Bool_t relative, Int_t base=1)
Int_t GetMCInfoIdFromIdealTrackId(int idealtrackid)
Int_t GetAssoRecoTrackID() const
TClonesArray * fSttHitArray
TH1 * fIdealTracksPerEvent
static const int kPossibleSec
TH1 * fQualyHisto_rel_all
std::map< Int_t, Int_t > GetTrackMCStatus()
PndTrackingQualityMCInfo GetMCInfoFromIdealTrack(PndTrack *idealtrack)
virtual void SetParContainers()
static const int kAtLeastThreePrim
static PndTrackFunctor * make_PndTrackFunctor(std::string functorName)
TString fPossibleTrackFunctorName
std::map< Int_t, TVector3 > GetP()
FairTrackParP GetParamLast()
TClonesArray * FillTubeArray()
std::map< Int_t, std::map< TString, std::pair< Double_t, Int_t > > > GetEfficiencies()
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
void SetIdealTrackId(int idealid)
std::map< Int_t, Double_t > GetPlResolutionRel()
std::map< Int_t, Double_t > GetPtResolution()
void SetMomentumFirst(TVector3 mom)
void SetVerbose(Int_t val)
static const int kMcLessThanThreePrim
std::map< int, int > fMCInfoIdIdealId
static const int kMcAllTracksWithHits
std::map< Int_t, Int_t > GetTrackQualification()
PndTrackFunctor * fPossibleTrackFunctor
static const int kFullyFound
std::vector< TString > fBranchNames
virtual ~PndTrackingQATask()
TString fIdealTrackBranchName
void SetMomentum(TVector3 val)
TString fMCInfoBranchName
void AssociateRecoTracksToMCTracks()
std::map< TString, TH2 * > fMapEfficiencies
void SetIsPrimary(Bool_t val)
Int_t GetIdealTrackIdFromMCTrackId(int mctrackid)
static const int kMcAtLeastThreePrim
void SetPositionLast(TVector3 pos)
void SetPositionFirst(TVector3 pos)
TTree * GetInternalTree()
virtual InitStatus Init()
std::map< Int_t, Double_t > GetPtResolutionRel()
TVector3 GetStartVertex() const
virtual void Exec(Option_t *opt)
PndTrackCand * GetTrackCandPtr()
Int_t GetMotherID() const
void SetMCQuality(int mcquality)
virtual Int_t GetSumOfAllValidMCHits(FairMultiLinkedData *trackData)
TString fRecoInfoBranchName
void AddHitsBranchName(TString name)
Adds branch names of detector data which should be taken into account in the analysis.
void PrintTrackQualityMap(Bool_t detailedInfo=kFALSE)
virtual void FillMCStatus(std::map< Int_t, Int_t > trackMCStatus)
void SetQuality(int quality)
void SetVertex(TVector3 val)
static const int kMcPossiblePrim
PndGeoSttPar * fSttParameters
FairTrackParP GetParamFirst()
static const int kSpuriousFound
static const int kNotFound
Int_t GetIdealTrackIdFromRecoTrackId(int trackid)