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
ClassImp(PndTrackingQATask)
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)