//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      TPC Helitron/Plawa Matching Routine
//
//
// Environment:
//      Software developed for the Prototype Detector at FOPI
//
// Author List:
//      Paul Buehler          SMI
//
//-----------------------------------------------------------



// This Class' Header ------------------
#include "TpcHelMatchingTask.h"

// C++ headers

// Collaborating Class Headers --------
#include "FairRun.h"
#include "FairRunAna.h"
#include "TpcDigiPar.h" 
#include "TpcCluster.h"
#include "FairRuntimeDb.h"

#include "TVectorD.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TLatex.h"
#include "TLegend.h"

#include "RKTrackRep.h"
#include "GFTrack.h"
#include "GFTrackCand.h"
#include "RecoHits/GFAbsRecoHit.h"
#include "PseudoSpacePoint.h"
#include "GFKalman.h"

#include "FOPIField.h"
#include "FopiEvent.h"
#include "CdcHit.h" 
#include "HeliHit.h" 
#include "HeliTrack.h" 
#include "PlawaTrack.h"
#include "FopiForwardHit.h"

#include "GFException.h"

#define DEBUG 0


//-----------------------------------------------------------
TpcHelMatchingTask::TpcHelMatchingTask()
  : fFopiEvBranchName("FopiEvent"),
    fTpcTrackBranchName("TpcPostFit"),    
    fHelHitBranchName("HeliHit"),
    fHelTrackBranchName("HeliTrack"),
    fPlaTrackBranchName("PlawaTrack"),
    fTpcHelMatchParsBranchName("TpcHelMatchPars"),
    fOutBranchName("TpcHelPreFit"),
    fdrMax(100.),
    fWithPlawaHit(kTRUE),
    fminHPmatch(0.),
    fPersistence(kTRUE),
    fVerbose(kFALSE),
		fSmoothing(kFALSE)
{	

  fRad = 3.1415926/180.;
  fNHeliHitstoAdd = 0;
  fPosErr = TVector3(1.,1.,0.1);
  fMomErr = TVector3(0.1,0.1,0.1);
  fZshift = 0.;
  fResScaleFac = 1;
  fwivScale = kFALSE;
  fHRes = TVector3(-1.,-1.,0.1);
  fPRes = TVector3(1.0,1.0,1.0);

}

//-----------------------------------------------------------
TpcHelMatchingTask::~TpcHelMatchingTask(){
}


//-----------------------------------------------------------
InitStatus
TpcHelMatchingTask::Init() {

  //Get ROOT Manager ------------------
  FairRootManager* ioman= FairRootManager::Instance();
  if(ioman==0)
    {
      Error("TpcHelMatchingTask::Init","RootManager not instantiated!");
      return kERROR;
    }
  
  // Get input ------------------------
  fFopiEventArray=(TClonesArray*) ioman->GetObject(fFopiEvBranchName);
  if(fFopiEventArray==0) {
    Error("TpcHelMatchingTask::Init","FOPI event not found!");
    return kERROR;
  } 
  fprintf(stderr,"fTpcTrackBranchName 1: %s\n",fTpcTrackBranchName.Data());
  
  fTpcTrackArray=(TClonesArray*) ioman->GetObject(fTpcTrackBranchName);
  if(fTpcTrackArray==0) {
    Error("TpcHelMatchingTask::Init","GFTrack array not found!");
    return kERROR;
  }
  fHelHitArray=(TClonesArray*) ioman->GetObject(fHelHitBranchName);
  if(fHelHitArray==0) {
    Error("TpcHelMatchingTask::Init","Helitron Hits array not found!");
    return kERROR;
  }  
  fHelTrackArray=(TClonesArray*) ioman->GetObject(fHelTrackBranchName);
  if(fHelTrackArray==0) {
    Error("TpcHelMatchingTask::Init","Helitron track array not found!");
    return kERROR;
  }  
  fPlaTrackArray=(TClonesArray*) ioman->GetObject(fPlaTrackBranchName);
  if(fPlaTrackArray==0) {
    Error("TpcHelMatchingTask::Init","Plawa track array not found!");
    return kERROR;
  }  

  // prepare trees for output
  MatchingPars = new Double_t[14];
  MatchingParsVec = new TVectorD(14);
  fMatchingPars = new TClonesArray("TVectorD");
  ioman->Register(fTpcHelMatchParsBranchName,"Tpc",fMatchingPars,fPersistence);
  
  fTpcHelTracks = new TClonesArray("GFTrack");
  ioman->Register(fOutBranchName,"Tpc",fTpcHelTracks,fPersistence);

  return kSUCCESS;
}

//-----------------------------------------------------------
void TpcHelMatchingTask::SetParContainers() {

  //  std::cout<<"TpcResidualTask::SetParContainers"<<std::endl;
  std::cout.flush();

  // Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

}

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

//-----------------------------------------------------------
Bool_t TpcHelMatchingTask::matching (Double_t *mpars) {

  // Definition of limits
  Double_t DThetaMax = 10.;
  Double_t DPhiMax = 90.;
  
  Double_t the1 = mpars[9];
  Double_t the2 = mpars[5];
  Double_t phi1 = mpars[10];
  if (phi1 < 0.) phi1 += 360.;
  Double_t phi2 = mpars[7];
  if (phi2 < 0.) phi2 += 360.;
  
  Double_t DThe = the1-the2;
  Double_t DPhi = phi1-phi2;
  if (fabs(DPhi) > 180.) DPhi = 360.-fabs(DPhi);

  if (fVerbose) {
    fprintf(stderr,"phi1/phi2: %f %f\n",phi1,phi2);
    fprintf(stderr,"the1/the2: %f %f\n",the1,the2);
    fprintf(stderr,"DThe/DPhi: %f %f\n",DThe,DPhi);
  }
  return  (fabs(DThe)<DThetaMax && fabs(DPhi)<DPhiMax);

}

//-----------------------------------------------------------
void TpcHelMatchingTask::Exec(Option_t* opt) {
  if (fVerbose) fprintf(stderr,"TpcHelMatchingTask:Exec\n");

  Int_t NumTPCTrackHits;
  Int_t NumHelTrackHits;
  TVector3 pos0, mom0;
  TMatrixD cov0;  
  Int_t hrsec;
  Double_t HelTrackHitsMin = 10.;

  // Definition of Helitron and Plawa hit covariance matrices
  TMatrixD HHcov = TMatrixD(3,3); 
  Double_t HHcovvals[9] = {0.,0.,0.,0.,0.,0.,0.,0.,pow(fHRes.Z(),2)};
  TMatrixD PHcov = TMatrixD(3,3); 
  Double_t PHcovvals[9] = {
    pow(fPRes.X(),2),0.,0.,0.,pow(fPRes.Y(),2),0.,0.,0.,pow(fPRes.Z(),2)};
  PHcov.SetMatrixArray(PHcovvals,"F");

  std::vector<unsigned int> HELcandIDs;
  std::vector<FopiForwardHit> helihits;

  // clean up result buffer
  if(fMatchingPars==0) 
    Fatal("TpcHelMatchingTask::Exec()","No fMatchingPars");
  fMatchingPars->Delete();
  if(fTpcHelTracks==0) 
    Fatal("TpcHelMatchingTask::Exec()","No fTpcHelTracks");
  fTpcHelTracks->Delete();

  //FairRunAna* run = FairRunAna::Instance();
  //FOPIField *fMagField = (FOPIField*) run->GetField();
  //Double_t Bz = fMagField->GetBz(0.,0.,0.);
  
	// get FOPI event
  FopiEvent *fopiEv = (FopiEvent*) fFopiEventArray->At(0);
  if (fopiEv == 0) return;
	
	Int_t nrun = fopiEv->GetRunNum();
	Int_t nspill = fopiEv->GetSpillNum();
	Int_t nevent = fopiEv->GetEvtNum();
  if (fVerbose)
    fprintf(stderr,"Event nrun, nspill, nevent: %i %i %i\n",nrun,nspill,nevent);
  
  // transformation Helitron-to-beam 
  Float_t taroff = fopiEv->GetTaroff();
  Float_t phi0 = fopiEv->GetPhi0();
  TVector3 hoff = fopiEv->GetHoff();
  TVector3 hrot = fopiEv->GetHrot();
  
	// define globhit and lochit
	FopiForwardHit globhit = FopiForwardHit(1,0,taroff,hoff,hrot);
	FopiForwardHit lochit  = FopiForwardHit(0,0,taroff,hoff,hrot);

  // how many tracks are there
  Int_t NumTpcTracks = fTpcTrackArray->GetEntriesFast();
  Int_t NumHelTracks = fopiEv->GetHmul();
	if (fVerbose)
		fprintf(stderr,"Number of Tpc and Helitron tracks: %i %i\n",
			NumTpcTracks,NumHelTracks);
  
  // prepare buffers
  Double_t the2use = 0.;
  Double_t *theTPCs = new Double_t[2*NumTpcTracks];
  Double_t *phiTPCs = new Double_t[2*NumTpcTracks];
  Double_t *momTPCs = new Double_t[NumTpcTracks];
  Double_t *theHels = new Double_t[NumHelTracks];
  Double_t *phiHels = new Double_t[NumHelTracks];
  Double_t *momHels = new Double_t[NumHelTracks];
  Double_t *dEdxHels = new Double_t[NumHelTracks];
  Double_t *HelPlaMats = new Double_t[NumHelTracks];
  
  // define reference plane for TPC track angle estimations
  // angles from the track fit are evaluated at plane0
  GFDetPlane *plane0 = new GFDetPlane(
    TVector3(0.,0.,0.),
    TVector3(1.,0.,0.),
    TVector3(0.,1.,0.));
  
  // start with Helitron tracks
  Int_t ntracks = -1;
  for (Int_t ii=0; ii<NumHelTracks; ii++) {
		HeliTrack *hsingletrack = (HeliTrack*) fHelTrackArray->At(ii);
    PlawaTrack *psingletrack = (PlawaTrack*) fPlaTrackArray->At(ii);

		// momentum and Helitron/PLAWA matching
    momHels[ii]    = hsingletrack->GetHmom();
    dEdxHels[ii]   = hsingletrack->GetHdedx();
    HelPlaMats[ii] = hsingletrack->GetPHmatch();
    
    // get track hit errors
    TVector3 HRes = fHRes;
    if (fHRes.X()<0.) HRes = TVector3(hsingletrack->GetHdx(),HRes.Y(),HRes.Z());
    if (fHRes.Y()<0.) HRes = TVector3(HRes.X(),hsingletrack->GetHdy(),HRes.Z());
    if (fVerbose) fprintf(stderr,"Hit errors: %f %f\n",HRes.X(),HRes.Y());
    
    // reference sector
    hrsec = (Int_t)hsingletrack->GetHrsec();
    lochit.SetHrsec(hrsec);
    globhit.SetHrsec(hrsec);

    if (fVerbose) {
      fprintf(stderr,"....................................................\n");
      fprintf(stderr,"Matching Helitron track (ref.sec. %i)> mass: %f mom: %f charge: %f\n",
        (Int_t)hrsec,hsingletrack->GetHmass(),momHels[ii],
        hsingletrack->GetHchar());
    }
    
    // compute theta and phi from track hits
    HELcandIDs.clear(); 
    HELcandIDs = hsingletrack->GetHelIdHits();
    NumHelTrackHits = HELcandIDs.size();
    if (fVerbose)
      fprintf(stderr,"Number of Helitron track hits: %i\n",NumHelTrackHits);
      
    TVector3 possum = TVector3(0.,0.,0.);
    theHels[ii] = 0.;
    phiHels[ii] = 0.;
    for (Int_t jj=0; jj<NumHelTrackHits; jj++) {
			HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(HELcandIDs[jj]);
      
      // convert to global coordinate system
      lochit.SetPos(0,
        hsinglehit->GetHhx(),hsinglehit->GetHhy(),hsinglehit->GetHhz()+fZshift);

      possum += lochit.GetGlobalPos();

    }
    theHels[ii] = atan2(possum.Perp(),possum.Z())/fRad;
    phiHels[ii] = atan2(possum.Y(),possum.X())/fRad;    
    
    // now loop over all TPC tracks
    if (fVerbose)
      fprintf(stderr,"Helitron> Phi: %8.3f           Theta    : %8.3f\n",
        phiHels[ii],theHels[ii]);
      
    for (Int_t jj=0; jj<NumTpcTracks; jj++) {
      GFTrack* tpcsingletrack = (GFTrack*) fTpcTrackArray->At(jj);
      NumTPCTrackHits = tpcsingletrack->getNumHits();
    
      // estimate theta and phi from track fit
      // extrapolate track to given z-position (plane0)
      //try {
      //  tpcsingletrack->getPosMomCov(*plane0,pos0,mom0,cov0);
      //} catch (GFException& e) {
      //  e.what();
      //}
      pos0 = tpcsingletrack->getPos();
      mom0 = tpcsingletrack->getMom();
      
      if (fVerbose) fprintf(stderr,"TPC mom: %f %f %f\n",mom0.X(),mom0.Y(),mom0.Z());
      momTPCs[jj]   = mom0.Mag();
      theTPCs[2*jj] = atan2(mom0.Perp(),mom0.Z())/fRad;
      phiTPCs[2*jj] = atan2(mom0.Y(),mom0.X())/fRad;
    
      // estimate theta and phi using hit points
      TVector3 pos1, pos2;
      possum = TVector3(0.,0.,0.);
      theTPCs[2*jj+1] = 0.;
      phiTPCs[2*jj+1] = 0.;
      if (NumTPCTrackHits > 1) {
        for (Int_t kk=0; kk<NumTPCTrackHits; kk++) {
          GFAbsRecoHit* hit1=tpcsingletrack->getHit(kk);
	  TVectorD raw(3);
	  raw = hit1->getRawHitCoord();
	  pos1.SetXYZ(raw[0],raw[1],raw[2]);
                  
          for (Int_t ll=kk+1; ll<NumTPCTrackHits; ll++) {
            GFAbsRecoHit* hit2=tpcsingletrack->getHit(ll);
	    raw = hit2->getRawHitCoord();
	    pos2.SetXYZ(raw[0],raw[1],raw[2]);
                    
            possum += (pos2-pos1);
           
          }
        }
        theTPCs[2*jj+1] = atan2(possum.Perp(),possum.Z())/fRad;
        phiTPCs[2*jj+1] = atan2(possum.Y(),possum.X())/fRad;
        
      }

      // save MatchingPars
      MatchingPars[1]  = (Float_t) jj;
      MatchingPars[2]  = (Float_t) ii;
      MatchingPars[3]  = (Float_t) hrsec;
      MatchingPars[4]  = theTPCs[2*jj];
      MatchingPars[5]  = theTPCs[2*jj+1];
      MatchingPars[6]  = phiTPCs[2*jj];
      MatchingPars[7]  = phiTPCs[2*jj+1];
      MatchingPars[8]  = momTPCs[jj];
      MatchingPars[9]  = theHels[ii];
      MatchingPars[10]  = phiHels[ii];
      MatchingPars[11] = momHels[ii];
      MatchingPars[12] = dEdxHels[ii];
      MatchingPars[13] = HelPlaMats[ii];

      // printout
      if (fVerbose) {
        if (fabs(MatchingPars[11])>0.) {
          fprintf(stderr,"TPC     >      %8.3f (%8.3f)           %8.3f (%8.3f)\n",
            MatchingPars[6],MatchingPars[7]+10./MatchingPars[11],
            MatchingPars[4],MatchingPars[5]);
        } else {
          fprintf(stderr,"TPC     >      %8.3f (%8.3f)           %8.3f (%8.3f)\n",
            MatchingPars[6],MatchingPars[7],MatchingPars[4],MatchingPars[5]);
        }
      
      }
        
      RKTrackRep* rep1 = NULL;
      RKTrackRep* rep2 = NULL;
      GFTrack* trk = NULL;

      // use MatchingPars to select matching tracks
      Int_t ismatched = -1;
      if (matching(MatchingPars) && NumHelTrackHits>HelTrackHitsMin) {
        if (fVerbose) fprintf(stderr,"ismatched!!\n");
      
        // compute best theta angle
        possum = TVector3(0.,0.,0.);
        the2use = 0.;
        for (Int_t kk=0; kk<NumTPCTrackHits; kk++) {
          GFAbsRecoHit* hit1=tpcsingletrack->getHit(kk);
	  TVectorD raw(3);
	  raw = hit1->getRawHitCoord();
	  pos1.SetXYZ(raw[0],raw[1],raw[2]);
	  
	  for (Int_t ll=0; ll<NumHelTrackHits; ll++) {
            HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(HELcandIDs[ll]);
      
            // convert to global coordinate system
            lochit.SetPos(0,
              hsinglehit->GetHhx(),hsinglehit->GetHhy(),
              hsinglehit->GetHhz()+fZshift);
            pos2 = TVector3(lochit.GetGlobalPos().X(),
              lochit.GetGlobalPos().Y(),lochit.GetGlobalPos().Z());
        
            possum += (pos2-pos1);
            
          }
        }
        the2use = atan2(possum.Perp(),possum.Z())/fRad;
      
        // create common track
        // estimate momentum from Helitron, TPC, and Dphi
        Double_t Dphi = MatchingPars[10]-MatchingPars[7];
        
        // momentum from Dphi
        Double_t momest = MatchingPars[11];
        if (momest <= 0.) momest = -10./Dphi;
        if (fabs(momest)<0.2) momest *= 0.2/fabs(momest);
        if (fabs(momest)>1.5) momest *= 1.5/fabs(momest);
        
        // what kind of particle is it?
        Int_t pdg1, pdg2;
				
        // use FOPI estimate *********************
				// Int_t pid = hsingletrack->GetHpid();
        // switch (pid) {
        //   case 1: pdg1 = 211;      // pion
        //           if (momest<0) pdg1 = -211;
        //           break;
        //   case 2: pdg1 = 2212;    // proton
        //           break;
        //   default:pdg1 = 211;
        //           if (momest<0) pdg1 = -211;
        // }
				// pdg2 = -pdg1;
        // use FOPI estimate ********************
				
				// use pion and proton hypothesis
				pdg1 = 211;
				if (momest<0) pdg1 = -211;
				pdg2 = 2212;
				
        // create track representation
        // track parameters have to be estimated
        GFAbsRecoHit* hit = tpcsingletrack->getHit(0);

	TVectorD raw(3);
	raw = hit->getRawHitCoord();
	// compute pos0, mom0
	pos0.SetXYZ(raw[0],raw[1],raw[2]);
        
        mom0.SetMagThetaPhi(fabs(momest),
        the2use*fRad,MatchingPars[7]*fRad);
        
        // get it from TPC track fit
        //pos0 = tpcsingletrack->getPos();
        //mom0 = tpcsingletrack->getMom();

        TVector3 momerr(mom0.X()*fMomErr.X(),mom0.Y()*fMomErr.Y(),
          mom0.Z()*fMomErr.Z()); 

        if (fVerbose) {
          fprintf(stderr,"Start values: \n");
          fprintf(stderr,"PDG: %i\n",pdg1);
          fprintf(stderr,"pos0: %f %f %f\n",
            pos0.X(),pos0.Y(),pos0.Z());
          fprintf(stderr,"mom0: %f %f %f\n",
            mom0.X(),mom0.Y(),mom0.Z());
          fprintf(stderr,"the2use: %f\n",the2use);
          fprintf(stderr,"Hel: %f, TPC: %f, used: %f\n\n",
            MatchingPars[11],momTPCs[jj],momest);
        }
        
        // define track representations
        rep1 = new RKTrackRep(pos0,mom0,fPosErr,momerr,pdg1);
        rep2 = new RKTrackRep(pos0,mom0,fPosErr,momerr,pdg2);
      
        // create new track
        trk = new GFTrack(rep1);
        trk->addTrackRep(rep2);
        trk->setCardinalRep(0);
				if (fSmoothing) {
					trk->setSmoothing(kTRUE);
				} else {
					trk->setSmoothing(kFALSE);
				}
      
        // add TPC hit points
        for (Int_t kk=0; kk<NumTPCTrackHits; kk++) {
          GFAbsRecoHit *TPCHit = tpcsingletrack->getHit(kk);
          trk->addHit(TPCHit);
	  TVectorD raw2(3);
	  raw2 = TPCHit->getRawHitCoord();
	            
          if (fVerbose) {
            TVector3 pos;
	    pos.SetXYZ(raw2[0],raw2[1],raw2[2]);
            fprintf(stderr,"TPC pos: %8.3f %8.3f %8.3f\n",
              pos.X(),pos.Y(),pos.Z());
          }
        
        }
      
        // get list of Helitron hits and sort according to 
        helihits.clear();
        for (Int_t kk=0; kk<NumHelTrackHits; kk++) {
          HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(HELcandIDs[kk]);
          lochit.SetPos(0,TVector3(
            hsinglehit->GetHhx(),hsinglehit->GetHhy(),
            hsinglehit->GetHhz()+fZshift));
          helihits.push_back(lochit);
        }
        sort (helihits.begin(),helihits.end(),sortcriterium);
      
        // add Helitron hit points
        PseudoSpacePoint *fwdhits = NULL;
        Int_t HeliHitstoAdd = fNHeliHitstoAdd;
        if (HeliHitstoAdd > NumHelTrackHits) HeliHitstoAdd = NumHelTrackHits;
        
        Int_t numfwdhits = HeliHitstoAdd;
        if (fWithPlawaHit && HelPlaMats[ii]>fminHPmatch) numfwdhits++;
        if (numfwdhits > 0) {
          fwdhits = new PseudoSpacePoint [numfwdhits];
        
          for (Int_t kk=0; kk<HeliHitstoAdd; kk++) {
            
            // determine hit covariance matrix
            Double_t rnom =
                helihits[kk].GetLocalPos().Z()*tan(MatchingPars[9]*fRad);
            Double_t rmes = sqrt(pow(helihits[kk].GetLocalPos().X(),2)+
              pow(helihits[kk].GetLocalPos().Y(),2));
            Double_t dr = fabs(rnom-rmes);
            
            // is this point added?
            if (dr < fdrMax) {
            
              Double_t efac = 1.0;
              if (fwivScale) efac = dr/HRes.Perp();
              HHcovvals[0] = pow(fResScaleFac*efac*HRes.X(),2);
              HHcovvals[4] = pow(fResScaleFac*efac*HRes.Y(),2);
              HHcovvals[8] = pow(fResScaleFac*HRes.Z(),2);
              HHcov.SetMatrixArray(HHcovvals,"F");
              helihits[kk].SetCov(0,HHcov);
              fwdhits[kk] = PseudoSpacePoint(&(helihits[kk]));
            
              trk->addHit(&fwdhits[kk]);
            
              if (fVerbose) {
                TVector3 pos = helihits[kk].GetGlobalPos();
                fprintf(stderr,"HEL pos: %8.3f/%8.3f %8.3f/%8.3f %8.3f/%8.3f : %8.3f\n",
                  pos.X(),fResScaleFac*efac*HRes.X(),
                  pos.Y(),fResScaleFac*efac*HRes.Y(),
                  pos.Z(),fHRes.Z(),dr);
              }

            }
            
          }
          // fprintf(stderr,"Number of hits on track: %i\n",trk->getNumHits());
        
          // add matched Plawa point to track
          fprintf(stderr,"Plawa hit: %i %f\n",fWithPlawaHit,HelPlaMats[ii]);
          if (fWithPlawaHit && HelPlaMats[ii]>fminHPmatch) {
            
            Float_t phthe = fRad*psingletrack->GetPhthe();
            Float_t phphi = fRad*psingletrack->GetPhphi();
            Float_t phz   = psingletrack->GetPhz()+fZshift;
            Float_t pr    = phz*tan(phthe);
            
            globhit.SetPos(1,pr*cos(phphi),pr*sin(phphi),phz);
            globhit.SetCov(1,PHcov);
            fwdhits[HeliHitstoAdd] = PseudoSpacePoint(&globhit);

            trk->addHit(&fwdhits[HeliHitstoAdd]);
        
            if (fVerbose) {
              TVector3 pos = globhit.GetGlobalPos();
              fprintf(stderr,"PLA pos: %8.3f %8.3f %8.3f\n",
                pos.X(),pos.Y(),pos.Z());
            }
          }
        }
        
        // create track on output array
        ntracks++;
        ismatched = ntracks;
        new((*fTpcHelTracks)[ntracks]) GFTrack(*trk);
        
        // delete trk; trk = NULL;
        delete rep2; rep2 = NULL;
        delete rep1; rep1 = NULL;
        delete [] fwdhits; fwdhits = NULL;
        
      }
      
      // save MatchingPars
      MatchingPars[0]  = (Float_t) ismatched;
      
      MatchingParsVec->SetElements(MatchingPars);
      Int_t ncombs = fMatchingPars->GetEntries();
			new((*fMatchingPars)[ncombs]) TVectorD(*MatchingParsVec);
      
    }
  }
  // fprintf(stderr,"number of track combinations: %i / %i\n",
	//   fMatchingPars->GetEntries(),ntracks);
    
  delete [] theTPCs;
  delete [] phiTPCs;
  delete [] momTPCs;
  delete [] theHels;
  delete [] phiHels;
  delete [] momHels;
  delete [] dEdxHels;
  delete [] HelPlaMats;
  delete plane0;
	  
}

ClassImp(TpcHelMatchingTask)

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