//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Select good combinations of TPC and Helitron Tracks
//
// Environment:
//      Processing of data from FOPI experiment S339
//
// Author List:
//      Paul Buehler                        SMI/OEAW
//
//-----------------------------------------------------------

#include "THPTrackFilter.h"

#include "TVectorD.h"
#include "TVector3.h"

#include "GFTrack.h"
#include <algorithm>

//-----------------------------------------------------------
// function for sorting of 3-dim points by radial distance
bool mapsorter(TVector3 h1, TVector3 h2) {
  return(h1.Z() < h2.Z());
}

//-----------------------------------------------------------
// uses matching pars and TPC/Helitron track fit to select best fitting
// TPC and Helitron track combinations
int THPTrackFilter (
  TClonesArray* matchpars, TClonesArray* tpcheltracks,
  Double_t fMaxChi2, std::vector<TVector3> *map)
{
   // initialisations
  bool fVerbose = kTRUE;
  TVector3 entry;
  TVectorD *MatchingParsVec = new TVectorD(14);
  
  // reset map
  map->clear(); 

  // fill map
  // [0]: index matchpars
  // [1]: index tpcheltracks
  // [2]: chi2
  Int_t NumMatchPars = matchpars->GetEntries();
  if (NumMatchPars == 0) return 0;
  
  for (Int_t ii=0; ii<NumMatchPars; ii++) {
    MatchingParsVec = (TVectorD*) matchpars->At(ii);
    Double_t *MatchingPars = MatchingParsVec->GetMatrixArray();
    Int_t trkindex = (Int_t)MatchingPars[0];
    if (trkindex >= 0) {
      GFTrack       *trk = (GFTrack*) tpcheltracks->At(trkindex);
      
      // decide which tracks are to be considered
      GFAbsTrackRep *rep = trk->getCardinalRep();
      if (trk->getRedChiSqu()<fMaxChi2 && rep->getStatusFlag()==0) {
      // if (trk->getRedChiSqu()<fMaxChi2) {
        entry = TVector3(ii,trkindex,rep->getRedChiSqu());
        map->push_back(entry);
      }
    }
  }
  // sort map according to chi2
  std::sort(map->begin(),map->end(),mapsorter);
  
  // clean up map and keep only best fitting tracks
  // . specific TPC and Helitron track can only be combined once
  // . select combination with minimum chi2  
  Int_t cc = 0;
  while (map->size() > cc) {
    // TPC and Helitron track indices corresponding to combination with
    // minimal chi2
    Double_t *MP = ((TVectorD*)matchpars->At((*map)[cc].X()))->GetMatrixArray();
    Int_t iiref = MP[1];
    Int_t jjref = MP[2];
  
    // eliminate all combinations containing either TPC track iiref or
    // Helitron track jjref
    Int_t pp = cc+1;
    while (map->size() > pp) {
      Double_t *MP2 = ((TVectorD*)matchpars->At((*map)[pp].X()))->GetMatrixArray();
      Int_t iiloc = MP2[1];
      Int_t jjloc = MP2[2];
      if (iiloc==iiref || jjloc==jjref) {
        map->erase(map->begin()+pp);
      } else {
        pp++;
      }
    }
    cc++;
  }
  // delete MatchingParsVec;
  
  return map->size();
  
}

//-----------------------------------------------------------
int THPTrackFilter (
  TClonesArray* matchpars,
  std::vector<TVector3> *map)
{
   // initialisations
  bool fVerbose = kTRUE;
  TVector3 entry;
  TVectorD *MatchingParsVec = new TVectorD(14);
  
  // reset map
  map->clear(); 

  // fill map
  // [0]: index matchpars
  // [1]: index tpcheltracks
  // [2]: DTheta
  Int_t NumMatchPars = matchpars->GetEntries();
  if (NumMatchPars == 0) return 0;
  
  for (Int_t ii=0; ii<NumMatchPars; ii++) {
    MatchingParsVec = (TVectorD*) matchpars->At(ii);
    Double_t *MatchingPars = MatchingParsVec->GetMatrixArray();
    Int_t trkindex = MatchingPars[0];
    if (trkindex >= 0) {
      Double_t DTheta = fabs(MatchingPars[5]-MatchingPars[9]);
      entry = TVector3(ii,trkindex,DTheta);
      map->push_back(entry);
    }
  }
  // sort map according to DTheta
  sort (map->begin(),map->end(),mapsorter);
  
  // clean up map and keep only best fitting tracks
  // . specific TPC and Helitron track can only be combined once
  // . select combination with minimum DTheta  
  Int_t cc = 0;
  while (map->size() > cc) {
    // TPC and Helitron track indices corresponding to combination with
    // minimal DTheta
    Double_t *MP = ((TVectorD*)matchpars->At((*map)[cc].X()))->GetMatrixArray();
    Int_t iiref = MP[1];
    Int_t jjref = MP[2];
  
    // eliminate all combinations containing either TPC track iiref or
    // Helitron track jjref
    Int_t pp = cc+1;
    while (map->size() > pp) {
      Double_t *MP2 = ((TVectorD*)matchpars->At((*map)[pp].X()))->GetMatrixArray();
      Int_t iiloc = MP2[1];
      Int_t jjloc = MP2[2];
      if (iiloc==iiref || jjloc==jjref) {
        map->erase(map->begin()+pp);
      } else {
        pp++;
      }
    }
    cc++;
  }
  // delete MatchingParsVec;
  
  return map->size();
  
}

//-----------------------------------------------------------