AMSMCEventClassifier.h

Go to the documentation of this file.

 
#include <root.h>
#include "amschain.h"
 
#include <TFile.h>
#include <TH1F.h>
#include <TH2.h>
 
#include <map>
#include <vector>
#include <iostream>
#include <unordered_map>
 
#include "Point.h"
#include "TRDConstants.h"
#include "Utilities.h"
#include "AMSTrdMCTrack.h"
#include "Intersection3D.h"
#include "ClusterFitter.h"
#include "LinearAlgebra.h"
#include "AMSTracker.h"
 
 
 
std::pair<int, double> GetTRDVertexClassifierOutput(AMSEventR * pev, TH2D * h_TrackerL1_XY, bool debugMode = false)
{
    int NLENGTH = 80;
    
    int TRDVertexStatus = -99; 
    
    
    double PrimaryParticleL1HitDistance =  -9;   
    
    int nFittedPrimaryMCTRDtrkInEvent = 0;
    int nPrimaryParticleL1Hit = 0;
 
 
 
    if (pev == NULL)
        return {-999, -90};
    
    
    
    MCEventgR * mc0 = pev->GetPrimaryMC();
    const int mc0_PID = mc0->Particle;  // primary particle type
    int mc0_trkID  = mc0->trkID;
    double mc0_mom = mc0->Momentum;
    int mc0_Nskip  = mc0->Nskip;
    double mc0_x   = mc0->Coo[0];
    double mc0_y   = mc0->Coo[1];
    double mc0_z   = mc0->Coo[2];
    // delete mc0;
    
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\n\nFor primary MC particle in this AMS event:\n");
        printf("g3 PID = %10d,  trkID = %10d, Momentum = %7.2f,  Nskip = %12d, (MCCooX, MCCooY, MCCooZ) = %8.2f  %8.2f  %8.2f\n",
                      mc0_PID,     mc0_trkID,          mc0_mom,     mc0_Nskip,                            mc0_x, mc0_y, mc0_z);
    }
    
    
    
    int nmceventg = pev->nMCEventg(); 
    std::vector<int> vMCtrkIDs = {};
    std::vector<Point::Point2DwID> vMCEventgClustersXY = {};  
    
 
    if (debugMode)
        printf("\nFound %d MCEventgR's in this AMS event.\n", nmceventg);
    
    for(int imc = 0; imc < nmceventg; imc++)
    {
        MCEventgR mc = pev->MCEventg(imc);
        vMCtrkIDs.push_back(mc.trkID);
        
        if ((mc.Coo[2] >= TrkL1_ZPos - TrkL1_ZPosCutoff)  &&  (mc.Coo[2] <= TrkL1_ZPos + TrkL1_ZPosCutoff)  &&    
            TrackerL1::HasL1XYHit(h_TrackerL1_XY, mc.Coo[0], mc.Coo[1]))                               
        {
            vMCEventgClustersXY.emplace_back(mc.Coo[0], mc.Coo[1], mc.trkID, mc.Particle, -1);  
        }
        
        if (debugMode)
        {
            printf("g3 PID = %10d,  trkID = %10d, Momentum = %7.2f,  Nskip = %12d, (MCCooX, MCCooY, MCCooZ) = %8.2f  %8.2f  %8.2f\n",
                      mc.Particle,      mc.trkID,      mc.Momentum,      mc.Nskip,      mc.Coo[0],      mc.Coo[1],      mc.Coo[2]);
        }
    } 
    
    auto MCEventgResultMap = VectorUtilities::GetUniqueElementsAndFrequency(vMCtrkIDs);
    
    if (debugMode)
    {
        printf("Unique trkIDs found:  %zu.\nDetails of trkIDs with nClusters > 1:\n", MCEventgResultMap.size());
        for(const auto& pair : MCEventgResultMap)
        {
            if (pair.second > 1)
                printf("trkID = %10d, nClusters = %10d\n", pair.first, pair.second);
        }
    }
    
    
    
    
 
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\nNow analyzing MC TRD track clusters in this event...");
    }
    
    int ntrdmccluster = pev->nTrdMCCluster();
    std::vector<int> vTRDMCtrkIDs = {};
    
    if (debugMode)
        printf("\nFound %d TrdMCClusterR's in this AMS event.\n", ntrdmccluster);
    
    if (ntrdmccluster == 0)
    {
        TRDVertexStatus = -29;
        return std::make_pair(TRDVertexStatus, PrimaryParticleL1HitDistance);
    }
    
    for(int iTRDMCCluster = 0; iTRDMCCluster < ntrdmccluster; iTRDMCCluster++)
    {
        TrdMCClusterR trdCltr = pev->TrdMCCluster(iTRDMCCluster);        
        vTRDMCtrkIDs.push_back(trdCltr.GtrkID);
 
        if (debugMode)
        {
            printf("g3 PID = %10d, GtrkID = %10d,     Ekin = %7.2f, ProcID = %12d, MCCooX = %8.2f, MCCooY = %8.2f, MCCooZ = %8.2f\n",
               trdCltr.ParticleNo,trdCltr.GtrkID,     trdCltr.Ekin,trdCltr.ProcID, trdCltr.Xgl[0], trdCltr.Xgl[1], trdCltr.Xgl[2]);
        }
    }
    
    std::map<int, int> TRDMCresultMap = VectorUtilities::GetUniqueElementsAndFrequency(vTRDMCtrkIDs);  
    if (debugMode)
        printf("\n\nUnique TRD trkIDs found:  %zu.", TRDMCresultMap.size());
    
    VectorUtilities::filterAndShrinkMap(TRDMCresultMap, nMIN_CLUSTERS_FOR_FIT);    
    unsigned int size_TRDMCresultMap = TRDMCresultMap.size(); 
    
    if (debugMode)
    {
        printf("\nGood TRD trkIDs found  :  %zu.\nDetails of good TRD trkIDs with nClusters > 1:\n", TRDMCresultMap.size());
        for(const auto& pair : TRDMCresultMap)
        {
            if (pair.second < nMIN_CLUSTERS_FOR_FIT) 
                continue;
 
            printf("TRD trkID = %10d, nClusters = %10d\n", pair.first, pair.second);
        }
    }
    
    
    
    
 
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\nNow grouping clusters under separate AMSTrdMCTrack::Tracks...");
    }
    
    std::vector<AMSTrdMCTrack::Track> vMCTracksInEvent = {}; 
 
    for(const auto& pair : TRDMCresultMap)
    {
        if (pair.second < nMIN_CLUSTERS_FOR_FIT) 
            continue;
        
        AMSTrdMCTrack::Track track;
        track.trkID = pair.first;
        
        std::vector<int> vClusterG3PID;
        
        for(int iTRDMCCluster = 0; iTRDMCCluster < ntrdmccluster; iTRDMCCluster++)
        {
            TrdMCClusterR trdCltr = pev->TrdMCCluster(iTRDMCCluster);
            
            if (pair.first != trdCltr.GtrkID)
                continue;
            
            vClusterG3PID.push_back(trdCltr.ParticleNo);
            
            AMSTrdMCTrack::Cluster cluster;
            cluster.x = trdCltr.Xgl[0];
            cluster.y = trdCltr.Xgl[1];
            cluster.z = trdCltr.Xgl[2];
            track.clusters.push_back(cluster);
        }
        
        
        if (VectorUtilities::areAllElementsEqual(vClusterG3PID) == true)
        {
            track.g3PID = vClusterG3PID[0];
            vMCTracksInEvent.push_back(track);
        }
        else if (VectorUtilities::countDifferentFromMajority(vClusterG3PID) == 1)
        {
            if (debugMode)
                printf("\n\nWARNING: In TRD trkID = %10d, nClusters = %10d, one vClusterG3PID element is different from the rest. Ignoring this element and moving on...\n", pair.first, pair.second);
            track.g3PID = VectorUtilities::getMajorityElement(vClusterG3PID);
            vMCTracksInEvent.push_back(track);
        }
        else
        {
            --size_TRDMCresultMap;
            if (debugMode)
            {
                printf("\n\n\n\nERROR! Track clusters do not have the same g3 PID!\n");
                printf("vClusterG3PID.size() = %zu\n", vClusterG3PID.size());
                printf("vClusterG3PID elements:\n");
                VectorUtilities::printVector(vClusterG3PID);
                printf("These clusters won't be used!\n\n\n\n"); 
            }
        }
        
        vClusterG3PID.clear(); vClusterG3PID.resize(0);
    }
    
    if (vMCTracksInEvent.size() != size_TRDMCresultMap)
        printf("\n\n\n\nERROR! Number of tracks saved does not equal unique tracks found!\n\nPlease confirm logic!\n\n\n\n"); 
    
    if (debugMode)
        printf("\n\nSaved %zu tracks in this AMS event for fitting:", vMCTracksInEvent.size());
    
    if (vMCTracksInEvent.size() == 0)
    {
        TRDVertexStatus = -29;
        return std::make_pair(TRDVertexStatus, PrimaryParticleL1HitDistance);
    }
    
    
 
    
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\nNow performing track fits...");
    }
    
    std::vector<std::vector<Line::FitLine2D>> vvFittedTracks = {}; 
    for (const AMSTrdMCTrack::Track& track : vMCTracksInEvent)
    {
        if (debugMode)
            printf("\n\ntrkID = %d, G3PID = %d, nClusters = %zu\n", track.trkID, track.g3PID, track.clusters.size());
        
        for (const AMSTrdMCTrack::Cluster& cluster : track.clusters)
        {
            if (debugMode)
                printf("---> Cluster: %9.3f %9.3f %9.3f\n", cluster.x, cluster.y, cluster.z);
        }
        
        std::vector<Line::FitLine2D> vFitted_Lines2D = ClusterFitter::GetClusterProjectionLinearFits(
        AMSTrdMCTrack::GetTrkPosValues(track, AMSTrdMCTrack::X),
        AMSTrdMCTrack::GetTrkPosValues(track, AMSTrdMCTrack::Y),
        AMSTrdMCTrack::GetTrkPosValues(track, AMSTrdMCTrack::Z),
        track.trkID,
        track.g3PID,
        debugMode,
        fitChi2perNDF_threshold);
        
        if (vFitted_Lines2D.size() == 3)
            vvFittedTracks.push_back(vFitted_Lines2D);
        
        // vFitted_Lines2D.clear(); vFitted_Lines2D.resize(0); /// Not required as it is not a vector of pointers.
    }
    
    unsigned int nFittedTracks = vvFittedTracks.size();
    if (nFittedTracks == 0) 
    {
        TRDVertexStatus = -29;
        return std::make_pair(TRDVertexStatus, PrimaryParticleL1HitDistance);
    }
    
    
    if (debugMode)
        printf("\n\n%d fitted AMSTrdMCTrack::Track found in event.\nVerifying track fit prameters:\n", nFittedTracks);
    
    for (unsigned int iFittedTrack = 0; iFittedTrack < nFittedTracks; iFittedTrack++)
    {
        if (debugMode)
        {
            for (unsigned int iPlane = 0; iPlane < 3; iPlane++)
            {
                printf("Track#: %3d (%5d, %5d) = Fit plane: %s, M = %14.7f, C = %14.7f\n",
                iFittedTrack, vvFittedTracks[iFittedTrack][iPlane].trkID, vvFittedTracks[iFittedTrack][iPlane].g3PID,
                vFitPlanes[iPlane].c_str(), vvFittedTracks[iFittedTrack][iPlane].m, vvFittedTracks[iFittedTrack][iPlane].c);
            }
        }
        
        if (vvFittedTracks[iFittedTrack][0].trkID == 1 && vvFittedTracks[iFittedTrack][0].g3PID == mc0_PID)
            ++nFittedPrimaryMCTRDtrkInEvent;
    }
    if (nFittedPrimaryMCTRDtrkInEvent < 1) 
    {
        TRDVertexStatus = -9;
        return std::make_pair(TRDVertexStatus, PrimaryParticleL1HitDistance);
    }
    
    
    
    
 
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\nNow finding 3D line from fitted plane projections...");
    }
    
    
    for (unsigned int iFittedTrack = 0; iFittedTrack < nFittedTracks; iFittedTrack++)
    {
        const double m1 = vvFittedTracks[iFittedTrack][0].m; const double b1 = vvFittedTracks[iFittedTrack][0].c;  
        const double m2 = vvFittedTracks[iFittedTrack][1].m; const double b2 = vvFittedTracks[iFittedTrack][1].c;  
        const double m3 = vvFittedTracks[iFittedTrack][2].m; const double b3 = vvFittedTracks[iFittedTrack][2].c;  
        const int fit_trkID = vvFittedTracks[iFittedTrack][0].trkID;
        const int fit_g3PID = vvFittedTracks[iFittedTrack][0].g3PID;
        
        if (debugMode)
        {
            printf("\n\n----> For fitted track# %d, (trkID = %d, g3PID = %d), the 2D projections provided are:\n", iFittedTrack, fit_trkID, fit_g3PID);
            printf("XZ:  m1, b1 = (%9.4f, %9.4f)\n", m1, b1);
            printf("YZ:  m2, b2 = (%9.4f, %9.4f)\n", m2, b2);
            printf("XY:  m3, b3 = (%9.4f, %9.4f)\n", m3, b3);
        }
        
        Line::Line3DParametricForm line3D = Line::Getline3DParametricForm(m1, b1, m2, b2, m3, b3, fit_trkID, fit_g3PID, debugMode);
        
        if (!line3D.point.empty() && !line3D.direction.empty())
        {
            if (debugMode)
            {
                // std::cout<<"The point on the line is: (" << line3D.point[0] << ", " << line3D.point[1] << ", " << line3D.point[2] <<")\n";
                // std::cout<<"The direction ratios of the line are: (" << line3D.direction[0] << ", " << line3D.direction[1] << ", " << line3D.direction[2] <<")\n";
            }
        }
        
        else
        {
            if (debugMode)
                std::cout<<"\nInvalid input or no solution exists for the given projections.\n";
            
            continue;
        }
        
        
        auto closestPoint = TrackerL1::FindClosestTrkL1ClusterFrom3DLine(h_TrackerL1_XY, line3D, vMCEventgClustersXY);
        
        if (closestPoint.second == false)
        {
            if (debugMode)
                std::cout<<"\nNo L1 hit found for this fitted track.\n\n";
        }
        
        else
        {
            if (debugMode)
                std::cout<<"\nFitted track potentially intersects with L1. Matching with MCEventgCluster hits...";
            
            if (closestPoint.first.trkID == fit_trkID && closestPoint.first.PID == fit_g3PID)
            {
                if (debugMode)
                    std::cout<<"\nThe closest L1 hit found: ("<<closestPoint.first.x<< ", "<<closestPoint.first.y<<") at "<<closestPoint.first.closedis<<" cm.\n";   
                
                if (fit_trkID ==  mc0_trkID && fit_g3PID == mc0_PID)
                {
                    ++nPrimaryParticleL1Hit;
                    PrimaryParticleL1HitDistance = closestPoint.first.closedis;
                }
            }
            
            else
            {
                if (debugMode)
                    std::cout<<"\nBut no matched L1 MCEventgCluster hit found for this fitted track.\n\n";
            }
        }
    }
    
    if (nPrimaryParticleL1Hit != 1)
        PrimaryParticleL1HitDistance = -5;
 
 
 
 
 
 
 
 
 
 
    int nMAX_POSSIBLE_PAIRS       = Utilities::CalculatePairs(nFittedTracks);        
    int nNEXTTOMAX_POSSIBLE_PAIRS = Utilities::CalculatePairs(nFittedTracks - 1);
    
    if (debugMode)
    {
        std::cerr<<std::string(NLENGTH, '_')<<std::endl;
        printf("\n\nFinding vertices in fitted intersecting tracks...");
        printf("\n\nFor %d fitted tracks in this event.\n", nFittedTracks);
        printf("nMAX_POSSIBLE_PAIRS       = %d\n", nMAX_POSSIBLE_PAIRS);
        printf("nNEXTTOMAX_POSSIBLE_PAIRS = %d\n", nNEXTTOMAX_POSSIBLE_PAIRS);
    }
    
    
    
    int nIntersecting_3Dlines = Intersection3D::GetNIntersectingLinePairs(vvFittedTracks, debugMode);
    if (debugMode)
        printf("\nnIntersecting 3D lines found = %d\n", nIntersecting_3Dlines);
    
    
    
    if (nIntersecting_3Dlines == 0)                  
    {
        if (nFittedTracks == 1)                      
            TRDVertexStatus = 1;
        else                                         
            TRDVertexStatus = 2;
    }
    
    else if (nIntersecting_3Dlines == nMAX_POSSIBLE_PAIRS && nFittedTracks > 1)
    {
        if ((nFittedTracks == 2) &&                                                         
            (vvFittedTracks[0][0].g3PID == -3 || vvFittedTracks[1][0].g3PID == -3) &&       
            (vvFittedTracks[0][0].g3PID * vvFittedTracks[1][0].g3PID < 0))                  
        {
            TRDVertexStatus = 3;
        }
        
        else
            TRDVertexStatus = -19;
    }
    
    else if (nIntersecting_3Dlines == nNEXTTOMAX_POSSIBLE_PAIRS && nFittedTracks > 2)
    {
        
        TRDVertexStatus = 0;
    }
    
    else if (nIntersecting_3Dlines > nNEXTTOMAX_POSSIBLE_PAIRS && nIntersecting_3Dlines < nMAX_POSSIBLE_PAIRS)
    {
        TRDVertexStatus = -2;
    }
    
    else if (nIntersecting_3Dlines > 0 && nIntersecting_3Dlines < nNEXTTOMAX_POSSIBLE_PAIRS && nFittedTracks > 3)
    {
        TRDVertexStatus = -1;
    }
    
    else
    {
        TRDVertexStatus = -39;
        if (debugMode)
                printf("\n\n\n\nWARNING: Unknown logic: unable to assign the TRDVertexStatus!\n\n\n\n"); 
    }
 
    
    if (debugMode)
        printf("\nEvent has been classified as TRDVertexStatus = %d\n", TRDVertexStatus);
    
    
    
    
    vMCtrkIDs.clear(); vMCtrkIDs.resize(0);
    vMCEventgClustersXY.clear(); vMCEventgClustersXY.resize(0);
    vTRDMCtrkIDs.clear(); vTRDMCtrkIDs.resize(0);
    vMCTracksInEvent.clear(); vMCTracksInEvent.resize(0); //std::vector<AMSTrdMCTrack::Track>().swap(vMCTracksInEvent);
    vvFittedTracks.clear(); vvFittedTracks.resize(0);
    
    return std::make_pair(TRDVertexStatus, PrimaryParticleL1HitDistance);
}