#include "LxEff.h"
#include "Lx.h"
#include <cmath>
#include <set>

//#define LXDRAW

#ifdef LXDRAW
#include "LxDraw.h"
//#define LXDRAWSAVE
//#define LXDRAWWAIT
#endif//LXDRAW

using namespace std;

LxEff::LxEff(LxFinderTriplet& owner) : finder(owner)
{
}

#ifdef CLUSTER_MODE
void LxEff::AnalyzeNotMatched(ofstream& out, list<LxPoint*> hits[LXSTATIONS][LXLAYERS])
{
  out << "AnalyzeNotMatched:" << endl;

  for (Int_t stNum = LXFIRSTSTATION; stNum < LXSTATIONS - 1; ++stNum)
  {
    out << "Segments #" << stNum << endl;
    list<LxPoint*> lHits[LXLAYERS] = hits[stNum];
    list<LxPoint*> rHits[LXLAYERS] = hits[stNum + 1];

    Double_t maxDx = 0;
    Double_t maxDy = 0;
    Double_t maxDtx = 0;
    Double_t maxDty = 0;
    Double_t minDx = 1000;
    Double_t minDy = 1000;
    Double_t minDtx = 1000;
    Double_t minDty = 1000;

    for (list<LxPoint*>::iterator l0 = lHits[0].begin(); l0 != lHits[0].end(); ++l0)
    {
      for (list<LxPoint*>::iterator l1 = lHits[1].begin(); l1 != lHits[1].end(); ++l1)
      {
        for (list<LxPoint*>::iterator l2 = lHits[2].begin(); l2 != lHits[2].end(); ++l2)
        {
          for (list<LxPoint*>::iterator r0 = rHits[0].begin(); r0 != rHits[0].end(); ++r0)
          {
            for (list<LxPoint*>::iterator r1 = rHits[1].begin(); r1 != rHits[1].end(); ++r1)
            {
              for (list<LxPoint*>::iterator r2 = rHits[2].begin(); r2 != rHits[2].end(); ++r2)
              {
                LxPoint* lPoints[LXLAYERS] = { *l0, *l1, *l2 };
                LxPoint* rPoints[LXLAYERS] = { *r0, *r1, *r2 };
                LxRay* rays[LXLAYERS * LXLAYERS] = {};

                for (Int_t i = 0; i < LXLAYERS; ++i)
                {
                  for (Int_t j = 0; j < LXLAYERS; ++j)
                    rays[j * LXLAYERS + i] = new LxRay(rPoints[j], lPoints[i], j * LXLAYERS + i);
                }

                for (Int_t i = 0; i < LXLAYERS * LXLAYERS; ++i)
                {
                  LxRay* ray1 = rays[i];
                  Double_t tx1 = ray1->tx;
                  Double_t ty1 = ray1->ty;
                  Double_t z1 = ray1->source->z;
                  Double_t x1 = ray1->source->x;
                  Double_t y1 = ray1->source->y;

                  if (minDx > ray1->source->dx)
                    minDx = ray1->source->dx;

                  if (minDy > ray1->source->dy)
                    minDy = ray1->source->dy;

                  if (minDtx > ray1->dtx)
                    minDtx = ray1->dtx;

                  if (minDty > ray1->dty)
                    minDty = ray1->dty;

                  for (Int_t j = i; j < LXLAYERS * LXLAYERS; ++j)
                  {
                    LxRay* ray2 = rays[j];
                    Double_t tx2 = ray2->tx;
                    Double_t ty2 = ray2->ty;
                    Double_t z2 = ray2->source->z;
                    Double_t diffZ = z1 - z2;
                    Double_t x2 = ray2->source->x + tx2 * diffZ;
                    Double_t y2 = ray2->source->y + ty2 * diffZ;

                    Double_t dtx = abs(tx2 - tx1);
                    Double_t dty = abs(ty2 - ty1);
                    Double_t dx = abs(x2 - x1);
                    Double_t dy = abs(y2 - y1);

                    if (maxDx < dx)
                      maxDx = dx;

                    if (maxDy < dy)
                      maxDy = dy;

                    if (maxDtx < dtx)
                      maxDtx = dtx;

                    if (maxDty < dty)
                      maxDty = dty;
                  }
                }

                for (Int_t i = 0; i < LXLAYERS * LXLAYERS; ++i)
                  delete rays[i];
              }
            }
          }
        }
      }
    }

    out << 4.0 * sqrt(minDx * minDx + finder.caSpace.stations[stNum + 1]->clusterXLimit2) - maxDx << " " <<
        4.0 * sqrt(minDy * minDy + finder.caSpace.stations[stNum + 1]->clusterYLimit2) - maxDy << " " <<
        4.0 * sqrt(minDtx * minDtx + finder.caSpace.stations[stNum + 1]->clusterTxLimit2) - maxDtx << " " <<
        4.0 * sqrt(minDty * minDty + finder.caSpace.stations[stNum + 1]->clusterTyLimit2) - maxDty << endl << endl;
  }
}

void LxEff::AnalyzeNotMatched2(ofstream& out, list<LxPoint*> hits[LXSTATIONS][LXLAYERS])
{
  out << "AnalyzeNotMatched2:" << endl;

  for (Int_t stNum = LXFIRSTSTATION; stNum < LXSTATIONS - 1; ++stNum)
  {
    out << "Segments #" << stNum << " has ";
    list<LxPoint*> lHits[LXLAYERS] = hits[stNum];
    list<LxPoint*> rHits[LXLAYERS] = hits[stNum + 1];

    Int_t maxOccupiedLevels = 0;
    Int_t maxNumPointsWithRays = 0;

    for (list<LxPoint*>::iterator l0 = lHits[0].begin(); l0 != lHits[0].end(); ++l0)
    {
      for (list<LxPoint*>::iterator l1 = lHits[1].begin(); l1 != lHits[1].end(); ++l1)
      {
        for (list<LxPoint*>::iterator l2 = lHits[2].begin(); l2 != lHits[2].end(); ++l2)
        {
          for (list<LxPoint*>::iterator r0 = rHits[0].begin(); r0 != rHits[0].end(); ++r0)
          {
            for (list<LxPoint*>::iterator r1 = rHits[1].begin(); r1 != rHits[1].end(); ++r1)
            {
              for (list<LxPoint*>::iterator r2 = rHits[2].begin(); r2 != rHits[2].end(); ++r2)
              {
                LxPoint* lPoints[LXLAYERS] = { *l0, *l1, *l2 };
                LxPoint* rPoints[LXLAYERS] = { *r0, *r1, *r2 };
                Int_t occupiedLevels = 0;
                Int_t numPointsWithRays = 0;

                for (Int_t i = 0; i < LXLAYERS; ++i)
                {
                  LxPoint* rPoint = rPoints[i];

                  if (!rPoint->rays.empty())
                    ++numPointsWithRays;

                  for (Int_t j = 0; j < LXLAYERS; ++j)
                  {
                    LxPoint* lPoint = lPoints[j];

                    for (list<LxRay*>::iterator k = rPoint->rays.begin(); k != rPoint->rays.end(); ++k)
                    {
                      LxRay* ray = *k;

                      if (ray->end == lPoint)
                      {
                        ++occupiedLevels;
                        break;
                      }
                    }
                  }
                }

                if (occupiedLevels > maxOccupiedLevels)
                  maxOccupiedLevels = occupiedLevels;

                if (numPointsWithRays > maxNumPointsWithRays)
                  maxNumPointsWithRays = numPointsWithRays;
              }
            }
          }
        }
      }
    }// for (list<LxPoint*>::iterator l0 = lHits[0].begin(); l0 != lHits[0].end(); ++l0)

    out << maxOccupiedLevels << " of occupied levels and " << maxNumPointsWithRays << " points with rays" << endl << endl;
  }
}

void LxEff::DumpNotMatchedMC(list<LxPoint*> hits[LXSTATIONS][LXLAYERS])
{
  AnalyzeNotMatched(finder.incomplete_events, hits);
  AnalyzeNotMatched2(finder.incomplete_events, hits);
  finder.incomplete_events << "Dumping not matched MC track" << endl;

  for (Int_t i = 0; i < LXSTATIONS; ++i)
  {
    finder.incomplete_events << "S" << i << " ";

    for (Int_t j = 0; j < LXLAYERS; ++j)
    {
      finder.incomplete_events << "L" << j << " ";

      for (list<LxPoint*>::iterator k = hits[i][j].begin(); k != hits[i][j].end(); ++k)
      {
        LxPoint* point = *k;
        finder.incomplete_events << "[";

        if (point->track)
          finder.incomplete_events << "t";
        else
        {
          if (point->leftClusteredRay)
            finder.incomplete_events << "r";

          finder.incomplete_events << "/";

          if (point->rightClusteredRay)
            finder.incomplete_events << "r";

        }

        finder.incomplete_events << "] ";
      }

      finder.incomplete_events << " ";
    }

    finder.incomplete_events << " ";
  }

  finder.incomplete_events << endl << endl;
}
#endif//CLUSTER_MODE

void LxEff::AnalyzeNotMatchedMC(ofstream& out, list<LxPoint*> hits[LXSTATIONS][LXLAYERS])
{
  out << "Not matched MC track:" << endl;
  out << "All triplets:   ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
  {
    Int_t tripletsCount = 0;

    for (list<LxPoint*>::iterator j = hits[i][1].begin(); j != hits[i][1].end(); ++j)
    {
      LxPoint* point = *j;
      tripletsCount += point->triplets.size();
    }

    out << "[" << tripletsCount << "]";
  }

  out << endl << "True triplets:  ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
  {
    Int_t trueTripletsCount = 0;

    for (list<LxPoint*>::iterator j = hits[i][1].begin(); j != hits[i][1].end(); ++j)
    {
      LxPoint* point = *j;

      for (list<LxTriplet*>::iterator k = point->triplets.begin(); k != point->triplets.end(); ++k)
      {
        LxTriplet* triplet = *k;
        LxPoint* lPoint = triplet->left;
        LxPoint* rPoint = triplet->right;
        bool trueLPoint = false;
        bool trueRPoint = false;

        for (list<LxPoint*>::iterator l = hits[i][0].begin(); l != hits[i][0].end(); ++l)
        {
          if (*l == lPoint)
            trueLPoint = true;
        }

        for (list<LxPoint*>::iterator l = hits[i][2].begin(); l != hits[i][2].end(); ++l)
        {
          if (*l == rPoint)
            trueRPoint = true;
        }

        if (trueLPoint && trueRPoint)
          ++trueTripletsCount;
      }
    }

    out << "[" << trueTripletsCount << "]";
  }

  out << endl << "Have neighbors: ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
  {
    out << "[";

    if (0 == i)
      out << " ";
    else
    {
      set<LxTriplet*> neighbours;

      for (list<LxPoint*>::iterator j = hits[i][1].begin(); j != hits[i][1].end(); ++j)
      {
        LxPoint* rPoint = *j;

        for (list<LxTriplet*>::iterator k = rPoint->triplets.begin(); k != rPoint->triplets.end(); ++k)
        {
          LxTriplet* rTriplet = *k;

          for (list<LxPoint*>::iterator l = hits[i - 1][1].begin(); l != hits[i - 1][1].end(); ++l)
          {
            LxPoint* lPoint = *l;

            for (list<LxTriplet*>::iterator m = lPoint->triplets.begin(); m != lPoint->triplets.end(); ++m)
            {
              LxTriplet* lTriplet = *m;

              for (list<pair<LxTriplet*, Double_t> >::iterator n = rTriplet->neighbours.begin(); n != rTriplet->neighbours.end(); ++n)
              {
                pair<LxTriplet*, Double_t>& tc2 = *n;

                if (tc2.first == lTriplet)
                  neighbours.insert(lTriplet);
              }
            }
          }
        }
      }

      out << neighbours.size();
    }

    out << "]";
  }

  Double_t maxCrDx[LXSTATIONS] = {};
  Double_t crDxDx2[LXSTATIONS] = {};
  Double_t maxCrDy[LXSTATIONS] = {};
  Double_t crDyDy2[LXSTATIONS] = {};
  Double_t maxClDx[LXSTATIONS] = {};
  Double_t clDxDx2[LXSTATIONS] = {};
  Double_t maxClDy[LXSTATIONS] = {};
  Double_t clDyDy2[LXSTATIONS] = {};

  for (Int_t i = 0; i < LXSTATIONS; ++i)
  {
    for (list<LxPoint*>::iterator j = hits[i][1].begin(); j != hits[i][1].end(); ++j)
    {
      LxPoint* cPoint = *j;
      Double_t cTx = cPoint->x / cPoint->z;
      Double_t cTy = cPoint->y / cPoint->z;

      for (list<LxPoint*>::iterator k = hits[i][2].begin(); k != hits[i][2].end(); ++k)
      {
        LxPoint* rPoint = *k;
        Double_t crDeltaZ = rPoint->z - cPoint->z;
        Double_t crTx = (rPoint->x - cPoint->x) / crDeltaZ;
        Double_t crTy = (rPoint->y - cPoint->y) / crDeltaZ;
        Double_t crDx = abs(rPoint->x - cPoint->x - cTx * crDeltaZ);
        Double_t crDy = abs(rPoint->y - cPoint->y - cTy * crDeltaZ);

        if (maxCrDx[i] < crDx)
        {
          maxCrDx[i] = crDx;
          crDxDx2[i] = cPoint->dx * cPoint->dx + rPoint->dx * rPoint->dx;
        }

        if (maxCrDy[i] < crDy)
        {
          maxCrDy[i] = crDy;
          crDyDy2[i] = cPoint->dy * cPoint->dy + rPoint->dy * rPoint->dy;
        }

        for (list<LxPoint*>::iterator l = hits[i][0].begin(); l != hits[i][0].end(); ++l)
        {
          LxPoint* lPoint = *l;
          Double_t clDeltaZ = lPoint->z - cPoint->z;
          Double_t clDx = abs(lPoint->x - cPoint->x - crTx * clDeltaZ);
          Double_t clDy = abs(lPoint->y - cPoint->y - crTy * clDeltaZ);

          if (maxClDx[i] < clDx)
          {
            maxClDx[i] = clDx;
            clDxDx2[i] = lPoint->dx * lPoint->dx + cPoint->dx * cPoint->dx + rPoint->dx * rPoint->dx;
          }

          if (maxClDy[i] < clDy)
          {
            maxClDy[i] = clDy;
            clDyDy2[i] = lPoint->dy * lPoint->dy + cPoint->dy * cPoint->dy + rPoint->dy * rPoint->dy;
          }
        }
      }
    }
  }

  out << endl << "Triplet C-->R X dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
    out << "[" << maxCrDx[i] / sqrt(finder.caSpace.stations[i]->xDispRight2 + crDxDx2[i]) << "]";

  out << endl << "Triplet C-->R Y dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
    out << "[" << maxCrDy[i] / sqrt(finder.caSpace.stations[i]->yDispRight2 + crDyDy2[i]) << "]";

  out << endl << "Triplet C-->L X dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
    out << "[" << maxClDx[i] / sqrt(finder.caSpace.stations[i]->xDispRL2 + clDxDx2[i]) << "]";

  out << endl << "Triplet C-->L Y dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS; ++i)
    out << "[" << maxClDy[i] / sqrt(finder.caSpace.stations[i]->yDispRL2 + clDyDy2[i]) << "]";

  Double_t maxDx[LXSTATIONS - 1] = {};
  Double_t dxDx2[LXSTATIONS - 1] = {};
  Double_t maxDy[LXSTATIONS - 1] = {};
  Double_t dyDy2[LXSTATIONS - 1] = {};
  Double_t maxDtx[LXSTATIONS - 1] = {};
  Double_t dtxDtx2[LXSTATIONS - 1] = {};
  Double_t maxDty[LXSTATIONS - 1] = {};
  Double_t dtyDty2[LXSTATIONS - 1] = {};

  for (Int_t i = 1; i < LXSTATIONS; ++i)
  {
    for (list<LxPoint*>::iterator j = hits[i][1].begin(); j != hits[i][1].end(); ++j)
    {
      LxPoint* rPoint = *j;

      for (list<LxPoint*>::iterator k = hits[i - 1][1].begin(); k != hits[i - 1][1].end(); ++k)
      {
        LxPoint* lPoint = *k;
        Double_t deltaZ = lPoint->z - rPoint->z;

        for (list<LxPoint*>::iterator l = hits[i][0].begin(); l != hits[i][0].end(); ++l)
        {
          LxPoint* rPoint0 = *l;

          for (list<LxPoint*>::iterator m = hits[i][2].begin(); m != hits[i][2].end(); ++m)
          {
            LxPoint* rPoint2 = *m;
            Double_t rdz = rPoint2->z - rPoint0->z;
            Double_t rdz2 = rdz * rdz;
            Double_t rtx = (rPoint2->x - rPoint0->x) / rdz;
            Double_t rty = (rPoint2->y - rPoint0->y) / rdz;
            Double_t rdtx2 = (rPoint0->dx2 + rPoint2->dx2) / rdz2;
            Double_t rdty2 = (rPoint0->dy2 + rPoint2->dy2) / rdz2;

#ifdef OUT_DISP_BY_TRIPLET_DIR
            Double_t dx = abs(lPoint->x - rPoint->x - rtx * deltaZ);
            Double_t dy = abs(lPoint->y - rPoint->y - rty * deltaZ);
#else//OUT_DISP_BY_TRIPLET_DIR
            Double_t dx = abs(lPoint->x - rPoint->x - deltaZ * rPoint->x / rPoint->z);
            Double_t dy = abs(lPoint->y - rPoint->y - deltaZ * rPoint->y / rPoint->z);
#endif//OUT_DISP_BY_TRIPLET_DIR

            if (dx > maxDx[i - 1])
            {
              maxDx[i - 1] = dx;
              dxDx2[i - 1] = lPoint->dx2 + rPoint->dx2;
            }

            if (dy > maxDy[i - 1])
            {
              maxDy[i - 1] = dy;
              dyDy2[i - 1] = lPoint->dy2 + rPoint->dy2;
            }

            for (list<LxPoint*>::iterator n = hits[i - 1][0].begin(); n != hits[i - 1][0].end(); ++n)
            {
              LxPoint* lPoint0 = *n;

              for (list<LxPoint*>::iterator o = hits[i - 1][2].begin(); o != hits[i - 1][2].end(); ++o)
              {
                LxPoint* lPoint2 = *o;
                Double_t ldz = lPoint2->z - lPoint0->z;
                Double_t ldz2 = ldz * ldz;
                Double_t ltx = (lPoint2->x - lPoint0->x) / ldz;
                Double_t lty = (lPoint2->y - lPoint0->y) / ldz;
                Double_t ldtx2 = (lPoint0->dx2 + lPoint2->dx2) / ldz2;
                Double_t ldty2 = (lPoint0->dy2 + lPoint2->dy2) / ldz2;

                Double_t dtx = abs(ltx - rtx);
                Double_t dty = abs(lty - rty);

                if (dtx > maxDtx[i - 1])
                {
                  maxDtx[i - 1] = dtx;
                  dtxDtx2[i - 1] = rdtx2 + ldtx2;
                }

                if (dty > maxDty[i - 1])
                {
                  maxDty[i - 1] = dty;
                  dtyDty2[i - 1] = rdty2 + ldty2;
                }
              }
            }
          }
        }
      }
    }
  }

  out << endl << "X dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS - 1; ++i)
#ifdef OUT_DISP_BY_TRIPLET_DIR
    out << "[" << maxDx[i] / sqrt(finder.caSpace.stations[i + 1]->xOutDispTriplet2 + dxDx2[i]) << "]";
#else//OUT_DISP_BY_TRIPLET_DIR
    out << "[" << maxDx[i] / sqrt(finder.caSpace.stations[i + 1]->xOutDispVertex2 + dxDx2[i]) << "]";
#endif//OUT_DISP_BY_TRIPLET_DIR

  out << endl << "Y dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS - 1; ++i)
#ifdef OUT_DISP_BY_TRIPLET_DIR
    out << "[" << maxDy[i] / sqrt(finder.caSpace.stations[i + 1]->yOutDispTriplet2 + dyDy2[i]) << "]";
#else//OUT_DISP_BY_TRIPLET_DIR
  out << "[" << maxDy[i] / sqrt(finder.caSpace.stations[i + 1]->yOutDispVertex2 + dyDy2[i]) << "]";
#endif//OUT_DISP_BY_TRIPLET_DIR

  out << endl << "Tx dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS - 1; ++i)
    out << "[" << maxDtx[i] / sqrt(finder.caSpace.stations[i + 1]->txInterTripletBreak2 + dtxDtx2[i]) << "]";

  out << endl << "Ty dispersions: ";

  for (Int_t i = 0; i < LXSTATIONS - 1; ++i)
    out << "[" << maxDty[i] / sqrt(finder.caSpace.stations[i + 1]->tyInterTripletBreak2 + dtyDty2[i]) << "]";

  out << endl;
}

static Double_t InterLinesDist(LxTrack* track1, LxTrack* track2)
{
  LxTriplet* triplet1 = track1->branches[0];
  LxPoint* point1 = triplet1->center;
  Double_t x1 = point1->x;
  Double_t y1 = point1->y;
  Double_t z1 = point1->z;

  //Double_t tx1 = point1->x / point1->z;
  //Double_t delta1 = triplet1->tx - tx1;

  Double_t ax = triplet1->tx;// - 2 * delta1;
  Double_t ay = triplet1->ty;
  Double_t az = 1.0;

  LxTriplet* triplet2 = track2->branches[0];
  LxPoint* point2 = triplet2->center;
  Double_t x2 = point2->x;
  Double_t y2 = point2->y;
  Double_t z2 = point2->z;

  //Double_t tx2 = point2->x / point2->z;
  //Double_t delta2 = triplet2->tx - tx2;

  Double_t bx = triplet2->tx;// - 2 * delta2;
  Double_t by = triplet2->ty;
  Double_t bz = 1.0;

  Double_t A = ay * bz - az * by;
  Double_t B = az * bx - ax * bz;
  Double_t C = ax * by - ay * bx;
  Double_t P = A * x2 + B * y2 + C * z2;

  Double_t ABCLen = sqrt(A * A + B * B + C * C);

  Double_t cosA = A / ABCLen;
  Double_t cosB = B / ABCLen;
  Double_t cosC = C / ABCLen;
  Double_t p = P / ABCLen;
  Double_t interTrackDistance = abs(cosA * x1 + cosB * y1 + cosC * z1 - p);

  return interTrackDistance;
}

static Double_t RecoTrackChi2(LxTrack* recoTrack)
{
  Double_t breakChi2 = 0;

  for (Int_t j = recoTrack->length - 1; j > 0; --j)
  {
    LxTriplet* rTriplet = recoTrack->branches[j];
    LxTriplet* lTriplet = recoTrack->branches[j - 1];

    LxStation* station = rTriplet->center->layer->station;
    breakChi2 += (lTriplet->tx - rTriplet->tx) * (lTriplet->tx - rTriplet->tx) / (station->txInterTripletBreak2 + lTriplet->dtx2 + rTriplet->dtx2) +
        (lTriplet->ty - rTriplet->ty) * (lTriplet->ty - rTriplet->ty) / (station->tyInterTripletBreak2 + lTriplet->dty2 + rTriplet->dty2);
  }

  breakChi2 /= 2 * (recoTrack->length - 1);

  return breakChi2;
}

void LxEff::CalcRecoEff(bool joinExt)
{
  Double_t result = 0;
  static Int_t signalRecoTracks = 0;
  static Int_t signalMCTracks = 0;
//  static Int_t bgrRecoTracks = 0;
//  static Int_t bgrMCTracks = 0;
  static Int_t mc2recoZeroMatched = 0;
  static Int_t mc2recoIncomplete = 0;

  bool hasPositiveMCTrack = false;
  bool hasNegativeMCTrack = false;

  for (vector<LxMCTrack>::iterator i = finder.MCTracks.begin(); i != finder.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];

    for (vector<LxMCPoint*>::iterator j = mcTrack.Points.begin(); j != mcTrack.Points.end(); ++j)
    {
      LxMCPoint* pMCPoint = *j;

      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;
      }
    }
    
    bool enoughHits = true;
    mcTrack.stationsWithHits = 0;
    mcTrack.layersWithHits = 0;

    for (Int_t j = 0; j < LXSTATIONS; ++j)
    {
      Int_t hitsOnSt = 0;

      for (Int_t k = 0; k < LXLAYERS; ++k)
      {
        if (mcTrack.hitsOnStations[j][k])
        {
          ++hitsOnSt;
          ++mcTrack.layersWithHits;
        }
      }

      if (hitsOnSt < 3)
      {
        if (j < finder.caSpace.stationsInAlgo)
          enoughHits = false;
      }
      else
        ++mcTrack.stationsWithHits;
    }

    Double_t pt2 = mcTrack.px * mcTrack.px + mcTrack.py * mcTrack.py;

    if (!enoughHits || mcTrack.mother_ID >= 0 || (pdgCode != 13 && pdgCode != -13) || (finder.pPtCut && (mcTrack.p < 3.0 || pt2 < 1.0)))
      isSignal = false;
    else
      ++signalMCTracks;

    if (isSignal)
    {
      if (-13 == pdgCode)
        hasPositiveMCTrack = true;
      else
        hasNegativeMCTrack = 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)
      {
        matchTrack = j->first;
        matchedPoints = j->second;
      }
    }

    if (0 == matchTrack)
    {
#ifdef MAKE_EFF_CALC
      //finder.incomplete_events << finder.eventNumber << " matched points: " << matchedPoints << " ; MC points: " << mcTrack.Points.size() << endl;
      //finder.caSpace.CheckArray(finder.incomplete_events, mcTrack);
#endif//MAKE_EFF_CALC
      //cout << "MC track does not have common point with a reconstructed" << endl;
      if (isSignal)
      {
        ++mc2recoZeroMatched;
        AnalyzeNotMatchedMC(finder.incomplete_events, mcTrackHits);
      }

#ifdef CLUSTER_MODE
      if (isSignal)
        DumpNotMatchedMC(mcTrackHits);
#endif//CLUSTER_MODE

      continue;
    }

    Int_t numberOfMuchMCPoints = LXLAYERS * finder.caSpace.stationsInAlgo;

    if (matchedPoints < 0.7 * numberOfMuchMCPoints)
    {
#ifdef MAKE_EFF_CALC
      //finder.incomplete_events << finder.eventNumber << " matched points: " << matchedPoints << " ; MC points: " << mcTrack.Points.size() << endl;
      //finder.caSpace.CheckArray(finder.incomplete_events, mcTrack);
#endif//MAKE_EFF_CALC
      //cout << "MC track have only " << matchedPoints << " common points with reconstructed" << endl;
      if (isSignal)
      {
        ++mc2recoIncomplete;
        AnalyzeNotMatchedMC(finder.incomplete_events, mcTrackHits);
      }

#ifdef CLUSTER_MODE
      if (isSignal)
        DumpNotMatchedMC(mcTrackHits);
#endif//CLUSTER_MODE

      continue;
    }

    matchTrack->mcTrack = &mcTrack;
    matchTrack->matched = true;

    if (!isSignal)
      continue;

#ifdef CALC_LINK_WITH_STS_EFF
    if (0 != matchTrack && matchedPoints >= 0.7 * numberOfMuchMCPoints)
      matchTrack->mcTracks.push_back(&mcTrack);
#endif//CALC_LINK_WITH_STS_EFF

    if (!joinExt)
      ++signalRecoTracks;
    else
    {
      if (0 == matchTrack->externalTrack)
      {
        if (0 == mcTrack.externalTrack)
          ++signalRecoTracks;
      }
      else if (matchTrack->externalTrack->track == mcTrack.externalTrack)
        ++signalRecoTracks;
    }
  }
  
  result = 100 * signalRecoTracks;
  result /= signalMCTracks;
  
#ifdef MAKE_EFF_CALC
  //if (foundNow != mcNow)
    //finder.incomplete_events << finder.eventNumber << endl;
#endif//MAKE_EFF_CALC

  cout << "LxEff::CalcRecoEff signal efficiency: " << result << "( " << signalRecoTracks << " / " << signalMCTracks << " )" << endl;
  cout << "LxEff::CalcRecoEff: zero-matched: " << mc2recoZeroMatched << ", incomplete: " << mc2recoIncomplete << endl;

  bool hasPositiveTrack = false;
  bool hasNegativeTrack = false;

  for (list<LxTrack*>::iterator i = finder.caSpace.tracks.begin(); i != finder.caSpace.tracks.end(); ++i)
  {
    LxTrack* recoTrack = *i;

    if (recoTrack->clone)
      continue;

    //if (recoTrack->length > finder.caSpace.stationsInAlgo)
      //continue;

    LxTriplet* firstTriplet = recoTrack->branches[0];
    LxPoint* firstPoint = firstTriplet->center;
    Double_t particleSign = firstTriplet->tx - firstPoint->x / firstPoint->z;
    LxTriplet* secondTriplet = recoTrack->branches[1];
    LxPoint* secondPoint = secondTriplet->center;
    Double_t tx = (secondPoint->x - firstPoint->x) / (secondPoint->z - firstPoint->z);
    Double_t particleSign2 = tx - firstPoint->x / firstPoint->z;
//    Double_t yAtZ0 = firstPoint->y - firstTriplet->ty * firstPoint->z;
//    Double_t breakChi2 = RecoTrackChi2(recoTrack);

    //if (-9 <= yAtZ0 && 9 >= yAtZ0)
    //if (breakChi2 <= 3)
    {
      if (particleSign > 0 && particleSign2 > 0)
        hasPositiveTrack = true;
      else if (particleSign < 0 && particleSign2 < 0)
        hasNegativeTrack = true;
    }
  }

  static Int_t numberOfJpsi = 0;

  if (hasPositiveTrack && hasNegativeTrack)
    ++numberOfJpsi;

  static Int_t numberOfMCJpsi = 0;
  static Int_t numberOfTrueJpsi = 0;

  if (hasPositiveMCTrack && hasNegativeMCTrack)
  {
    ++numberOfMCJpsi;

    if (hasPositiveTrack && hasNegativeTrack)
      ++numberOfTrueJpsi;
  }

  cout << "LxEff::CalcRecoEff: number of J/Psi hypoteses " << numberOfJpsi << endl;
  result = 100 * numberOfTrueJpsi;
  result /= numberOfMCJpsi;
  cout << "LxEff::CalcRecoEff: J/Psi efficiency " << result << " % ( " << numberOfTrueJpsi << " / " << numberOfMCJpsi << " )" << endl;

  return;

  // The code below is mostly dedicated to counting ghosts.
  static Int_t recoTracks = 0;
  static Int_t matchedMCTracks = 0;
  static Int_t completeleyUnmatched = 0;
  static Int_t notEnoughPoints = 0;
  static Int_t notEnoughStations = 0;
  static Int_t incompleteStations[LXSTATIONS] = {};
  static Int_t doesntHaveStsPart = 0;
  static Int_t doesntHaveStsPart1 = 0;
  static Int_t doesntHaveStsPart2 = 0;
  static Int_t unmatchedButHaveStsPart = 0;

  for (list<LxTrack*>::iterator i = finder.caSpace.tracks.begin(); i != finder.caSpace.tracks.end(); ++i)
  {
    //++recoTracks;
    LxTrack* recoTrack = *i;

    if (recoTrack->clone)
      continue;

    map<LxMCTrack*, Int_t> mcTracks;// Mapped value is the number of common points in MC and reconstructed tracks.

#ifdef CLUSTER_MODE
    for (Int_t ii = 0; ii < recoTrack->length; ++ii)
    {
      LxRay* ray = recoTrack->rays[ii];

      for (list<LxPoint*>::iterator j = ray->clusterPoints.begin(); j != ray->clusterPoints.end(); ++j)
      {
        LxPoint* recoPoint = *j;
#else//CLUSTER_MODE
    for (Int_t j = 0; j < finder.caSpace.stationsInAlgo * LXLAYERS; ++j)
    {
      LxPoint* recoPoint = recoTrack->points[j];
#endif//CLUSTER_MODE

      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);
      }
#ifdef CLUSTER_MODE
    }
#endif//CLUSTER_MODE
    }// for (Int_t j = 0; j < (LXSTATIONS - LXFIRSTSTATION) * 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)
    {
      ++completeleyUnmatched;
      cout << "This reconstructed track is not intersected with any MC track in any point" << endl;

      if (!joinExt || 0 != recoTrack->externalTrack)
        ++recoTracks;

      continue;
    }

    ++recoTracks;

    bool unmatched = false;

    if (matchedPoints < 0.7 * finder.caSpace.stationsInAlgo * LXLAYERS)
    {
      ++notEnoughPoints;
      cout << "This reconstructed track match with an MC track in only " << matchedPoints << "points" << endl;
      unmatched = true;
    }

    if (bestMatch->stationsWithHits < finder.caSpace.stationsInAlgo)
    {
      ++notEnoughStations;

#ifdef LXDRAW
      LxDraw drawer;
      drawer.ClearView();
      drawer.DrawInputHits();
      drawer.DrawMCTracks();
      drawer.DrawRays();

#ifdef LXDRAWSAVE
      char drawSaveName[64];
      sprintf(drawSaveName, "MC_Rays_%d_", finder.eventNumber);
      drawer.SaveCanvas(drawSaveName);
#endif//LXDRAWSAVE

#ifdef LXDRAWWAIT
      gPad->WaitPrimitive();
#endif//LXDRAWWAIT
#endif//LXDRAW

      for (Int_t j = 0; j < finder.caSpace.stationsInAlgo; ++j)
      {
        Int_t hitsOnStation = 0;

        for (Int_t k = 0; k < LXLAYERS; ++k)
        {
          if (bestMatch->hitsOnStations[j][k])
            ++hitsOnStation;
        }

        if (hitsOnStation < 3)
          ++incompleteStations[j];
      }

      cout << "The matched MC track has hits not on all stations. Only on: " << bestMatch->stationsWithHits << endl;
      unmatched = true;
    }

    if (0 == bestMatch->externalTrack)
    {
      ++doesntHaveStsPart;

      if (unmatched)
        ++doesntHaveStsPart1;
      else
        ++doesntHaveStsPart2;
    }

    if (!unmatched)
    {
      if (!joinExt)
        ++matchedMCTracks;
      else
      {
        if (0 == recoTrack->externalTrack)
        {
          if (0 == bestMatch->externalTrack)
            ++matchedMCTracks;
          else
            cout << "LxEff::CalcRecoEff: STS track not found" << endl;
        }
        else
        {
          if (bestMatch->externalTrack == recoTrack->externalTrack->track)
            ++matchedMCTracks;
          else
            cout << "LxEff::CalcRecoEff: STS track determined wrongly" << endl;
        }
      }
    }
    else if (0 != bestMatch->externalTrack)
      ++unmatchedButHaveStsPart;
  }// for (list<LxTrack*>::iterator i = finder.caSpace.tracks.begin(); i != finder.caSpace.tracks.end(); ++i)

  result = 100 * matchedMCTracks;
  result /= recoTracks;
  cout << "LxEff::CalcRecoEff matching reconstucted efficiency: " << result << "( " << matchedMCTracks << " / " << recoTracks << " )" << endl;
  cout << "LxEff::CalcRecoEff 0-matched: " << completeleyUnmatched << ", not enough points: " << notEnoughPoints <<
      ", not enough stations: " << notEnoughStations << " : ";

  for (Int_t i = 0; i < finder.caSpace.stationsInAlgo; ++i)
    cout << "[" << incompleteStations[i] << "]";

  cout << endl;
  cout << "LxEff::CalcRecoEff doesn't have an STS part: " << doesntHaveStsPart <<
      ", unmatched: " << doesntHaveStsPart1 << ", matched: " << doesntHaveStsPart2 << endl;
  cout << "LxEff::CalcRecoEff unmatched but have an STS part: " << unmatchedButHaveStsPart << endl;
}// void LxEff::CalcRecoEff(bool joinExt)

void LxFinderTriplet::MatchRecoToMC()
{
  static Int_t recoTracks = 0;
  static Int_t matchedTracks = 0;
  LxTrack* positiveSignal = 0;
  Double_t positiveDelta = 0;
  Double_t negativeDelta = 0;
  LxTrack* negativeSignal = 0;
  list<LxTrack*> positiveBgrs;
  list<LxTrack*> negativeBgrs;

  for (list<LxTrack*>::iterator i = caSpace.tracks.begin(); i != caSpace.tracks.end(); ++i)
  {
    LxTrack* recoTrack = *i;

    if (recoTrack->clone)
      continue;

    ++recoTracks;

    map<LxMCTrack*, Int_t> mcTracks;// Mapped value is the number of common points in MC and reconstructed tracks.

    for (Int_t j = 0; j < caSpace.stationsInAlgo * 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 (list<LxMCPoint*>::iterator k = recoPoint->mcPoints.begin(); k != recoPoint->mcPoints.end(); ++k)
    }// for (Int_t j = 0; j < caSpace.stationsInAlgo * 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;
      }
    }

    LxTriplet* firstTriplet = recoTrack->branches[0];
    LxPoint* firstPoint = firstTriplet->center;
    Double_t particleSign = firstTriplet->tx - firstPoint->x / firstPoint->z;
    Double_t yAtZ0 = firstPoint->y - firstTriplet->ty * firstPoint->z;
    LxTriplet* secondTriplet = recoTrack->branches[1];
    LxPoint* secondPoint = secondTriplet->center;
    Double_t yAtZ0_2 = firstPoint->y - firstPoint->z * (secondPoint->y - firstPoint->y) / (secondPoint->z - firstPoint->z);
    Double_t breakChi2 = RecoTrackChi2(recoTrack);

    if (0 == bestMatch || matchedPoints < 0.7 * caSpace.stationsInAlgo * LXLAYERS)
    {
      if (particleSign > 0)
        positiveBgrs.push_back(recoTrack);
      else if (particleSign < 0)
        negativeBgrs.push_back(recoTrack);

      bgrYAtZ0->Fill(yAtZ0);
      bgrYAtZ0_2->Fill(yAtZ0_2);
      backgroundChi2Histo->Fill(breakChi2);

      continue;
    }

    ++matchedTracks;

    if (0 > bestMatch->mother_ID && (-13 == bestMatch->pdg || 13 == bestMatch->pdg))
    {
      if (particleSign > 0 && -13 == bestMatch->pdg)
      {
        positiveSignal = recoTrack;
        positiveDelta = particleSign;
      }
      else if (particleSign < 0 && 13 == bestMatch->pdg)
      {
        negativeSignal = recoTrack;
        negativeDelta = particleSign;
      }

      signalYAtZ0->Fill(yAtZ0);
      signalYAtZ0_2->Fill(yAtZ0_2);
      signalChi2Histo->Fill(breakChi2);
    }
    else
    {
      if (particleSign > 0)
        positiveBgrs.push_back(recoTrack);
      else if (particleSign < 0)
        negativeBgrs.push_back(recoTrack);

      bgrYAtZ0->Fill(yAtZ0);
      bgrYAtZ0_2->Fill(yAtZ0_2);
      backgroundChi2Histo->Fill(breakChi2);
    }
  }// for (list<LxTrack*>::iterator i = caSpace.tracks.begin(); i != caSpace.tracks.end(); ++i)

  Double_t result = 100 * matchedTracks;
  result /= recoTracks;
  cout << "LxFinderTriplet::MatchRecoToMC(): efficiency: " << result << " % ( " << matchedTracks << " / " << recoTracks << " )" << endl;

  if (0 != positiveSignal && 0 != negativeSignal)
  {
    signalInterTracksDistance->Fill(InterLinesDist(positiveSignal, negativeSignal));
    signalSignDefect->Fill(positiveDelta + negativeDelta);
  }

  LxTrack* posBgr = 0;
  LxTrack* negBgr = 0;
  Double_t bgrDist = 10000;
  Double_t signDefect = 10000;

  for (list<LxTrack*>::iterator i = positiveBgrs.begin(); i != positiveBgrs.end(); ++i)
  {
    LxTrack* track1 = *i;

    for (list<LxTrack*>::iterator j = negativeBgrs.begin(); j != negativeBgrs.end(); ++j)
    {
      LxTrack* track2 = *j;
      Double_t dist = InterLinesDist(track1, track2);

      if (dist < bgrDist)
      {
        posBgr = track1;
        negBgr = track2;
        bgrDist = dist;
      }

      LxTriplet* firstTriplet1 = track1->branches[0];
      LxPoint* firstPoint1 = firstTriplet1->center;
      Double_t particleSign1 = firstTriplet1->tx - firstPoint1->x / firstPoint1->z;

      LxTriplet* firstTriplet2 = track2->branches[0];
      LxPoint* firstPoint2 = firstTriplet2->center;
      Double_t particleSign2 = firstTriplet2->tx - firstPoint2->x / firstPoint2->z;

      Double_t defect = particleSign1 + particleSign2;

      if (abs(defect) < abs(signDefect))
        signDefect = defect;
    }
  }// for (list<LxTrack*>::iterator i = positiveBgrs.begin(); i != positiveBgrs.end(); ++i)

  if (0 != posBgr && 0 != negBgr)
  {
    bgrInterTracksDistance->Fill(bgrDist);
    bgrSignDefect->Fill(signDefect);
  }
}// void LxFinderTriplet::MatchRecoToMC()