Simple/LxEff.cxx

Go to the documentation of this file.

#include "Lx.h"
#ifdef MAKE_EFF_CALC
//#include <cmath>
#include <iostream>
 
using namespace std;
 
void LxFinder::MatchMCToReco() {
  static Int_t referenceMCTracks      = 0;
  static Int_t reconstructedRefTracks = 0;
  bool hasPositiveSignalMuon          = false;
  bool hasNegativeSignalMuon          = false;
  hasSignalInEvent                    = false;
 
  for (vector<LxMCTrack>::iterator i = MCTracks.begin(); i != MCTracks.end();
       ++i) {
    LxMCTrack& mcTrack = *i;
    Int_t pdgCode      = mcTrack.pdg;
    map<LxTrack*, Int_t>
      recoTracks;  // Mapped value is the number of common points in MC and reconstructed tracks.
    memset(mcTrack.hitsOnStations, 0, sizeof(mcTrack.hitsOnStations));
    bool isSignal = true;
    //list<LxPoint*> mcTrackHits[LXSTATIONS][LXLAYERS];
    bool mcTrackMCPoints[LXSTATIONS][LXLAYERS];
    memset(mcTrackMCPoints, 0, sizeof(mcTrackMCPoints));
 
    for (vector<LxMCPoint*>::iterator j = mcTrack.Points.begin();
         j != mcTrack.Points.end();
         ++j) {
      LxMCPoint* pMCPoint                                             = *j;
      mcTrackMCPoints[pMCPoint->stationNumber][pMCPoint->layerNumber] = true;
 
      for (list<LxPoint*>::iterator k = pMCPoint->lxPoints.begin();
           k != pMCPoint->lxPoints.end();
           ++k) {
        LxPoint* point = *k;
        mcTrack.hitsOnStations[point->layer->station->stationNumber]
                              [point->layer->layerNumber] = true;
        //mcTrackHits[point->layer->station->stationNumber][point->layer->layerNumber].push_back(point);
        LxTrack* track = point->track;
 
        if (0 == track) continue;
 
        if (track->matched) continue;
 
        if (track->clone) continue;
 
        map<LxTrack*, int>::iterator l = recoTracks.find(track);
 
        if (l != recoTracks.end())
          ++(l->second);
        else
          recoTracks[track] = 1;
      }  // for (list<LxPoint*>::iterator k = pMCPoint->lxPoints.begin(); k != pMCPoint->lxPoints.end(); ++k)
    }  // for (vector<LxMCPoint*>::iterator j = mcTrack.Points.begin(); j != mcTrack.Points.end(); ++j)
 
    bool enoughHits          = true;
    bool enoughMCPoints      = true;
    mcTrack.stationsWithHits = 0;
    mcTrack.layersWithHits   = 0;
    Int_t mcpCount           = 0;
 
    for (Int_t j = 0; j < LXSTATIONS; ++j) {
      Int_t hitsOnSt = 0;
      Int_t mcpOnSt  = 0;
 
      for (Int_t k = 0; k < LXLAYERS; ++k) {
        if (mcTrack.hitsOnStations[j][k]) {
          ++hitsOnSt;
          ++mcTrack.layersWithHits;
        }
 
        if (mcTrackMCPoints[j][k]) {
          ++mcpOnSt;
          ++mcpCount;
        }
      }  // for (Int_t k = 0; k < LXLAYERS; ++k)
 
      if (hitsOnSt < LXLAYERS) {
        if (j < caSpace.stationsInAlgo) enoughHits = false;
      } else
        ++mcTrack.stationsWithHits;
 
      if (j < caSpace.stationsInAlgo && mcpOnSt < LXLAYERS)
        enoughMCPoints = false;
    }  // for (Int_t j = 0; j < LXSTATIONS; ++j)
 
    scaltype pt2 = mcTrack.px * mcTrack.px + mcTrack.py * mcTrack.py;
 
    if (!enoughMCPoints || mcTrack.mother_ID >= 0
        || (pdgCode != 13 && pdgCode != -13)
        || (pPtCut && (mcTrack.p < 3.0 || pt2 < 1.0)))
      isSignal = false;
    else
      ++referenceMCTracks;
 
    if (!isSignal) continue;
 
    if (-13 == pdgCode)
      hasPositiveSignalMuon = true;
    else
      hasNegativeSignalMuon = true;
 
    LxTrack* matchTrack = 0;
    Int_t matchedPoints = 0;
 
    for (map<LxTrack*, Int_t>::iterator j = recoTracks.begin();
         j != recoTracks.end();
         ++j) {
      if (0 == matchTrack
          || matchedPoints < j->second + j->first->restoredPoints) {
        matchTrack    = j->first;
        matchedPoints = j->second + j->first->restoredPoints;
      }
    }
 
    if (0 == matchTrack) continue;
 
    if (matchedPoints < mcpCount * 0.7) continue;
 
    ++reconstructedRefTracks;
  }  // for (vector<LxMCTrack>::iterator i = MCTracks.begin(); i != MCTracks.end(); ++i)
 
  if (hasPositiveSignalMuon && hasNegativeSignalMuon) {
    hasSignalInEvent = true;
    ++signalCounter;
  }
 
  Double_t result = 100 * reconstructedRefTracks;
  result /= referenceMCTracks;
  cout << "LxFinder::MatchMCToReco(): efficiency = " << result << "% ("
       << reconstructedRefTracks << "/" << referenceMCTracks << ")" << endl;
}  // void LxFinder::MatchMCToReco()
 
void LxFinder::MatchRecoToMC() {
  static Int_t recoTracksCount    = 0;
  static Int_t matchedTracksCount = 0;
  LxTrack* signalMuPlus           = 0;
  LxTrack* signalMuMinus          = 0;
  list<LxTrack*> bgrMuPlusTracks;
  list<LxTrack*> bgrMuMinusTracks;
 
  for (list<LxTrack*>::iterator i = caSpace.tracks.begin();
       i != caSpace.tracks.end();
       ++i) {
    LxTrack* recoTrack = *i;
 
    if (recoTrack->clone) continue;
 
    ++recoTracksCount;
 
    map<LxMCTrack*, Int_t> mcTracks;
 
    for (Int_t j = 0; j < LXSTATIONS * LXLAYERS; ++j) {
      LxPoint* recoPoint = recoTrack->points[j];
 
      if (0 == recoPoint) continue;
 
      for (list<LxMCPoint*>::iterator k = recoPoint->mcPoints.begin();
           k != recoPoint->mcPoints.end();
           ++k) {
        LxMCPoint* mcPoint = *k;
        LxMCTrack* mcTrack = mcPoint->track;
 
        if (0 == mcTrack) continue;
 
        map<LxMCTrack*, Int_t>::iterator mcIter = mcTracks.find(mcTrack);
 
        if (mcIter == mcTracks.end())
          mcTracks[mcTrack] = 1;
        else
          ++(mcIter->second);
      }
    }  // for (Int_t j = 0; j < LXSTATIONS * LXLAYERS; ++j)
 
    LxMCTrack* bestMatch = 0;
    Int_t matchedPoints  = 0;
 
    for (map<LxMCTrack*, Int_t>::iterator j = mcTracks.begin();
         j != mcTracks.end();
         ++j) {
      if (0 == bestMatch || j->second > matchedPoints) {
        bestMatch     = j->first;
        matchedPoints = j->second;
      }
    }
 
    if (0 == bestMatch) continue;
 
    Int_t recoPoints = 0;
 
    for (Int_t j = 0; j < LXSTATIONS * LXLAYERS; ++j) {
      LxPoint* recoPoint = recoTrack->points[j];
 
      if (0 != recoPoint) ++recoPoints;
    }
 
    if (matchedPoints < 0.7 * recoPoints) continue;
 
    ++matchedTracksCount;
 
    LxRay* ray0     = recoTrack->rays[0];
    LxPoint* point0 = ray0->end;
    scaltype x0     = point0->x - point0->z * ray0->tx;
    scaltype y0     = point0->y - point0->z * ray0->ty;
 
    if (0 != bestMatch && bestMatch->mother_ID < 0
        && (-13 == bestMatch->pdg || 13 == bestMatch->pdg)) {
      signalXAtZ0->Fill(x0);
      signalYAtZ0->Fill(y0);
 
      if (-13 == bestMatch->pdg)
        signalMuPlus = recoTrack;
      else
        signalMuMinus = recoTrack;
    } else {
      bgrXAtZ0->Fill(x0);
      bgrYAtZ0->Fill(y0);
      scaltype charge = ray0->tx - point0->x / point0->z;
 
      if (charge > 0)
        bgrMuPlusTracks.push_back(recoTrack);
      else if (charge < 0)
        bgrMuMinusTracks.push_back(recoTrack);
    }
  }  // for (list<LxTrack*>::iterator i = caSpace.tracks.begin(); i != caSpace.tracks.end(); ++i)
 
  Double_t result = 100 * matchedTracksCount;
  result /= recoTracksCount;
  cout << "LxFinder::MatchRecoToMC(): efficiency = " << result << "% ("
       << matchedTracksCount << "/" << recoTracksCount << ")" << endl;
 
  if (0 != signalMuPlus && 0 != signalMuMinus) {
    scaltype x1     = signalMuPlus->rays[0]->end->x;
    scaltype y1     = signalMuPlus->rays[0]->end->y;
    scaltype z1     = signalMuPlus->rays[0]->end->z;
    scaltype x2     = signalMuMinus->rays[0]->end->x;
    scaltype y2     = signalMuMinus->rays[0]->end->y;
    scaltype z2     = signalMuMinus->rays[0]->end->z;
    scaltype deltaX = x1 - x2;
    scaltype deltaY = y1 - y2;
    scaltype deltaZ = z1 - z2;
    scaltype d      = sqrt(deltaX * deltaX + deltaY * deltaY);
    signalInterTracksDistanceOn1st->Fill(d);
    scaltype aSq  = x1 * x1 + y1 * y1 + z1 * z1;
    scaltype a    = sqrt(aSq);
    scaltype bSq  = x2 * x2 + y2 * y2 + z2 * z2;
    scaltype b    = sqrt(bSq);
    scaltype cSq  = deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ;
    scaltype cosC = (aSq + bSq - cSq) / (2 * a * b);
    scaltype C    = acos(cosC);
    C *= 180;
    C /= 3.14159265;
    signalInterTracksAngle->Fill(C);
    signalInterTrackCorrDA->Fill(d, C);
  }  // if (0 != signalMuPlus && 0 != signalMuMinus)
 
  scaltype bgrMuDistSq = 0;
  scaltype bgrMuAngle  = 0;
 
  for (list<LxTrack*>::iterator i = bgrMuPlusTracks.begin();
       i != bgrMuPlusTracks.end();
       ++i) {
    LxTrack* mupt = *i;
    scaltype x1   = mupt->rays[0]->end->x;
    scaltype y1   = mupt->rays[0]->end->y;
    scaltype z1   = mupt->rays[0]->end->z;
 
    for (list<LxTrack*>::iterator j = bgrMuMinusTracks.begin();
         j != bgrMuMinusTracks.end();
         ++j) {
      LxTrack* mumt   = *j;
      scaltype x2     = mumt->rays[0]->end->x;
      scaltype y2     = mumt->rays[0]->end->y;
      scaltype z2     = mumt->rays[0]->end->z;
      scaltype deltaX = x1 - x2;
      scaltype deltaY = y1 - y2;
      scaltype deltaZ = z1 - z2;
      scaltype dSq    = deltaX * deltaX + deltaY * deltaY;
 
      if (dSq > bgrMuDistSq) bgrMuDistSq = dSq;
 
      scaltype aSq  = x1 * x1 + y1 * y1 + z1 * z1;
      scaltype a    = sqrt(aSq);
      scaltype bSq  = x2 * x2 + y2 * y2 + z2 * z2;
      scaltype b    = sqrt(bSq);
      scaltype cSq  = deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ;
      scaltype cosC = (aSq + bSq - cSq) / (2 * a * b);
      scaltype C    = acos(cosC);
      C *= 180;
      C /= 3.14159265;
 
      if (C > bgrMuAngle) bgrMuAngle = C;
    }  // for (list<LxTrack*>::iterator j = bgrMuMinusTracks.begin(); j != bgrMuMinusTracks.end(); ++j)
  }  // for (list<LxTrack*>::iterator i = bgrMuPlusTracks.begin(); i != bgrMuPlusTracks.end(); ++i)
 
  if (bgrMuDistSq > 0) bgrInterTracksDistanceOn1st->Fill(sqrt(bgrMuDistSq));
 
  if (bgrMuAngle > 0) bgrInterTracksAngle->Fill(bgrMuAngle);
 
  if (bgrMuDistSq > 0 && bgrMuAngle > 0)
    bgrInterTrackCorrDA->Fill(sqrt(bgrMuDistSq), bgrMuAngle);
}  // void LxFinder::MatchRecoToMC()
 
#endif  //MAKE_EFF_CALC