//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class PndTpcRecoTester
//      see PndTpcRecoTester.hh for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Sebastian Neubert    TUM            (original author)
//
//
//-----------------------------------------------------------

// Panda Headers ----------------------

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

// C/C++ Headers ----------------------
#include <iostream>
#include <map>
#include <algorithm>

// Collaborating Class Headers --------
#include "PndTpcCluster.h"
#include "Kalman.h"
#include "PndTpcClusterRadius.h"
#include "PndTpcConfMapRecoHit.h"
#include "PndTpcConfMapFit.h"
#include "DetPlane.h"
#include "PndTpcConfTrackFinder.h"
#include "TrackCand.h"
#include "Track.h"
#include "TCovEllipse.h"
#include "TGraph.h"
#include "TFile.h"
#include "TTree.h"
#include "TRandom.h"
#include "TH1I.h"
#include "TF2.h"
#include "TH2D.h"
#include "TLorentzVector.h"
#include "TPolyMarker3D.h"
#include "TPolyLine3D.h"
#include "DebugLogger.h"
#include "dbgstream.h"
#include "TMatrixT.h"
#include "TVectorD.h"
#include "StdDiscriminantFcn.h"
#include "PndTpcRiemannHit.h"
#include "PndTpcRiemannTrack.h"
#include "PndTpcRiemannTrackFinder.h"
#include "PndTpcRiemannHTCorrelator.h"
#include "PndTpcProximityHTCorrelator.h"
#include "TMath.h"
// Class Member definitions -----------

ClassImp(PndTpcRecoTester)

void
PndTpcRecoTester::testConfMap(){
  std::vector<PndTpcCluster*> cll; // clusterlist
  maketoytracks(2,cll);
  
  TGraph* confhits=new TGraph(cll.size());
  TGraph* hits=new TGraph(cll.size());
  
  for(int i=0;i<cll.size();++i){
    hits->SetPoint(i,cll[i]->pos().X(),cll[i]->pos().Y());
    PndTpcConfMapRecoHit hit(cll[i]);
    confhits->SetPoint(i,hit.getXcf(),
		       hit.getHitCoord(DetPlane())[0][0]);
  }

  cll.clear();
  
  maketoytracks(2,cll,true);
  TGraph* confhits2=new TGraph(cll.size());
  TGraph* hits2=new TGraph(cll.size());
  
  for(int i=0;i<cll.size();++i){
    hits2->SetPoint(i,cll[i]->pos().X(),cll[i]->pos().Y());
    PndTpcConfMapRecoHit hit(cll[i]);
    confhits2->SetPoint(i,hit.getXcf(),
		       hit.getHitCoord(DetPlane())[0][0]);
  }


  TFile* outfile=new TFile("confhits.root","RECREATE");
  confhits->SetTitle("Conformal Hits");
  confhits->Write("confhits");
  hits->SetTitle("Hits");
  hits->Write("hits");

  confhits2->SetTitle("Conformal Hits from IP");
  confhits2->Write("confhitsIP");
  hits2->SetTitle("Hits from IP");
  hits2->Write("hitsIP");

  outfile->Close();
  return;
}

void
PndTpcRecoTester::testConfMapFit(){
  PndTpcConfMapFit myfit;
  TMatrixT<double> cov(2,2);
  //cov[0][0]=10;
  //cov[1][1]=10;
  //cov[0][1]=-1;
  //cov[1][0]=-1;
  //myfit.setCov(cov);
  myfit.Print();
  PndTpcCluster mycl(TVector3(5,5,3),1,0);
  PndTpcConfMapRecoHit hit(&mycl);
  
  
 
  std::cout<<"Residual:"<<std::endl;
  hit.residualVector(&myfit,myfit.getState()).Print();

  Kalman myfitter;
  double chi2=myfitter.getChi2Hit(&hit,&myfit);
  std::cout<<"Chi2="<<chi2<<std::endl;

}

void
PndTpcRecoTester::testBasicConfFit(){
  PndTpcConfMapRecoHit hit1(0.5,0.51);
  PndTpcConfMapRecoHit hit2(0.7,0.68);
  PndTpcConfMapRecoHit hit3(1.2,1.22);
  PndTpcConfMapRecoHit hit4(1.2,1.24);
  PndTpcConfMapRecoHit hit5(1.2,1.26);
  PndTpcConfMapRecoHit hit6(1.2,1.28);

  Track* trk=new Track(new PndTpcConfMapFit());
  TMatrixT<double> par(2,1);
  par[0][0]=1;
  par[1][0]=0;
  trk->getTrackRep(0)->setState(par);
  trk->getTrackRep(0)->Print();
  
   Kalman myfitter;
  double chi2=myfitter.getChi2Hit(&hit1,trk->getTrackRep(0));
  std::cout<<"Chi2 for hit1="<<chi2<<std::endl;

  trk->addHit(&hit1,2,0);
  myfitter.continueTrack(trk);
  trk->getTrackRep(0)->Print();

  
  chi2=myfitter.getChi2Hit(&hit2,trk->getTrackRep(0));
  std::cout<<"Chi2 for hit2="<<chi2<<std::endl;

  trk->addHit(&hit2,2,1);
  myfitter.continueTrack(trk);
  trk->getTrackRep(0)->Print();
  
  chi2=myfitter.getChi2Hit(&hit3,trk->getTrackRep(0));
  std::cout<<"Chi2 for hit3="<<chi2<<std::endl;
  chi2=myfitter.getChi2Hit(&hit4,trk->getTrackRep(0));
 std::cout<<"Chi2 for hit4="<<chi2<<std::endl;
  chi2=myfitter.getChi2Hit(&hit5,trk->getTrackRep(0));
 std::cout<<"Chi2 for hit5="<<chi2<<std::endl;
  chi2=myfitter.getChi2Hit(&hit6,trk->getTrackRep(0));
 std::cout<<"Chi2 for hit6="<<chi2<<std::endl;
}

int
PndTpcRecoTester::testConfFinder(){

  std::vector<PndTpcCluster*> cll; // clusterlist
  maketoytracks(2,cll);

  //std::cout<<cll.size()<<" clusters produced."<<std::endl;

  PndTpcConfTrackFinder _finder;
  _finder.configure(1,1,2, // proximity cuts
		    3,     // chi2cut
		    4,     // numhitsforfit 
		    false,
		    0.05,0.05, // merge cuts
		    4.);   // proximity merge cut
  
  std::vector<TrackCand*> cands;
  _finder.buildTracks(cll,cands);

  //std::cout<<cands.size()<<" track candidates found."<<std::endl;
  for(unsigned int i=0; i<cands.size(); ++i){
    std::cout<<" candidate"<<i<<" with "<<cands[i]->getNHits()<<" hits and radius "<<1./cands[i]->getCurv()<<std::endl;
  }

  // draw graphs
  std::cout<< _finder.getTrkListRef().size()<<" tracks used in finder"<<std::endl;
  TGraph* confhits=new TGraph();
  TGraph* zshits=new TGraph();
  unsigned int hitcounter=0;
  std::vector<Track*>& trks=_finder.getTrkListRef(); // beware using reference!
  for(unsigned int i=0;i<trks.size();++i){// loop over tracks
    if(trks[i]==NULL)continue;
    if(dynamic_cast<PndTpcConfMapFit*>(trks[i]->getTrackRep(0))->isRotated())
      std::cout<<"Track "<<i<<" is rotated"<<std::endl;
    unsigned int nh=trks[i]->getNumHits();
    for(unsigned int ih=0;ih<nh;++ih){//loop over hits in track
      PndTpcConfMapRecoHit* conf=dynamic_cast<PndTpcConfMapRecoHit*>(trks[i]->getHit(ih));
      confhits->SetPoint(hitcounter,conf->getXcf(),
			 conf->getHitCoord(DetPlane())[0][0]);
      zshits->SetPoint(hitcounter++,conf->s(),
			 conf->z());
    }// end loop over hits in track

  }// end loop over tracks
  TFile* outfile=new TFile("confhits.root","RECREATE");
  confhits->Write("confhits");
  zshits->Write("zshits");
  outfile->Close();
  return cands.size();
}

void
PndTpcRecoTester::testRotated(){
  PndTpcConfMapRecoHit hit1(new PndTpcCluster(TVector3(5,5,3),1,0));
  std::cout<<hit1.getHitCoord(DetPlane())[0][0]<<" "
	   <<hit1.getXcf()<<std::endl;
  hit1.setRotated();
  std::cout<<hit1.getHitCoord(DetPlane())[0][0]<<" "
	   <<hit1.getXcf()<<std::endl;
  hit1.setRotated(false);
  std::cout<<hit1.getHitCoord(DetPlane())[0][0]<<" "
	   <<hit1.getXcf()<<std::endl;
}


void
PndTpcRecoTester::testConfFinderStat(){
  TH1I* histo=new TH1I("nfound","nfound",30,0,30);
  for(int ev=0;ev<1000;++ev){
    histo->Fill(testConfFinder());
  }
  TFile* outfile=new TFile("toystat.root","RECREATE");
  histo->Write();
  outfile->Close();
}


void
PndTpcRecoTester::testConfTrackMerge(int nevents, bool split, bool merge){
  
  
  PndTpcConfTrackFinder _finder;
  _finder.configure(1,1,2, // proximity cuts
		    3,     // chi2cut
		    3,     // numhitsforfit 
		    merge, // do merge
		    0.05,0.05, // merge cuts
		    4.);   // proximity merge cut
  if(split)_finder.setMaxHitsInTrack(20);
  std::vector<TMatrixT<double> > confresults;
  std::vector<TMatrixT<double> > confcov;
  _finder.setConfFitResults(&confresults,&confcov);

  
  

  // make event
  std::vector<double> rv;
  std::vector<TVector2> dirv;
  std::vector<double> dipv;
  std::vector<TVector2> startv;
  //maketoyevent(2,rv,dirv,dipv);
  rv.push_back(50);
  rv.push_back(51);
  TVector2 dir(1,0.3);dir=dir.Unit();
  dirv.push_back(dir);
  TVector2 dir2(1.1,0.6);dir2=dir2.Unit();
  dirv.push_back(dir2);
  dipv.push_back(1);
  dipv.push_back(1);
  TVector2 start1(0,0);startv.push_back(start1);
  TVector2 start2(0,0);startv.push_back(start2);


  for(int evts=0; evts<nevents; ++evts){
    std::vector<PndTpcCluster*> cll; // clusterlist
    maketoyhits(2,rv,dirv,dipv,startv,cll);
    std::vector<TrackCand*> cands;
    std::vector<TrackCand*> cands2;
    _finder.buildTracks(cll,cands);
  }

  TFile* outfile=new TFile("toystat.root","RECREATE");
  TTree* tr=new TTree("stattree","stattree");
  double p1;
  double p2;
  tr->Branch("p1",&p1,"p1/D"); 
  tr->Branch("p2",&p2,"p2/D"); 
  std::cout<<confcov.size()<<std::endl;
  for(int i=0;i<confresults.size();++i){
    p1=confresults[i][0][0];
    p2=confresults[i][1][0];
    tr->Fill();
    TCovEllipse* el=new TCovEllipse(confcov[i],p1,p2);
    confcov[i].Print();
    std::cout<<"p1="<<p1<<"  p2="<<p2<<std::endl;
    char name[20];sprintf(name,"el%i",i);
    el->Write(name);
  }

  tr->Write();
  outfile->Close();
}

//************************************************************************
//****************** discrimant analysis *********************************

// param[0]=m[0], param[1]=m[1], param[2]=cov[0][0], param[3]=cov[1][1],
// param[4]=cov[1][0]=cov[0][1]
double ftf2(Double_t *x, Double_t *params)
{
  TVectorT<double> m(2);
  m[0]=params[0];
  m[1]=params[1];
  TMatrixT<double> cov(2,2);
  cov[0][0]=params[2];
  cov[1][1]=params[3];
  cov[0][1]=params[4];
  cov[1][0]=cov[0][1];
  StdDiscriminantFcn dis(m,cov);
  TVectorT<double> xv(2,x);
  return dis(xv);
}

double dis2(Double_t *x, Double_t *params)
{
  // class1 
  TVectorT<double> m1(2);
  m1[0]=params[0];
  m1[1]=params[1];
  TMatrixT<double> cov1(2,2);
  cov1[0][0]=params[2];
  cov1[1][1]=params[3];
  cov1[0][1]=params[4];
  cov1[1][0]=cov1[0][1];
  StdDiscriminantFcn dis1(m1,cov1);

  // class2
  TVectorT<double> m2(2);
  m2[0]=params[5];
  m2[1]=params[6];
  TMatrixT<double> cov2(2,2);
  cov2[0][0]=params[7];
  cov2[1][1]=params[8];
  cov2[0][1]=params[9];
  cov2[1][0]=cov2[0][1];
  StdDiscriminantFcn dis2(m2,cov2);

  TVectorT<double> xv(2,x);
  return dis1(xv)-dis2(xv);
}

double dis3(Double_t *x, Double_t *params)
{
  double d1=ftf2(x, params);
  double d2=ftf2(x, &(params[5]));
  double d3=ftf2(x, &(params[10]));
  double m1= d1>d2 ? d1 : d2;
  return m1>d3 ? m1 : d3;
}

double dis3_cut(Double_t *x, Double_t *params)
{
  double d=dis3(x,params);
  return params[15]>d ? params[15] : d;
}

void
PndTpcRecoTester::testStdDis(double* par)
{
  TF2* f=new TF2("dis",&ftf2,0,4,0,4,5);
  double pars[16]={2,2,0.4,0.8,0.5,
		   5,1,0.4,0.8,0.5,
		   3,4,0.1,0.2,0.1,-2};
  if(par!=NULL)for(int i=0;i<16;++i)pars[i]=par[i];
  f->SetParameters(pars);
  TMatrixT<double> cov(2,2);
  cov[0][0]=pars[2];
  cov[1][1]=pars[3];
  cov[0][1]=pars[4];
  cov[1][0]=cov[0][1];
  TCovEllipse* el=new TCovEllipse(cov,pars[0],pars[1]);


  TF2* f2=new TF2("dis2",&dis2,-5,10,-5,10,10);
  f2->SetParameters(pars);
  TMatrixT<double> cov2(2,2);
  cov2[0][0]=pars[7];
  cov2[1][1]=pars[8];
  cov2[0][1]=pars[9];
  cov2[1][0]=cov2[0][1];
  TCovEllipse* el2=new TCovEllipse(cov2,pars[5],pars[6]);

  
  TF2* f3=new TF2("dis3",&dis3,-5,10,-5,10,15);
  f3->SetParameters(pars);
  TMatrixT<double> cov3(2,2);
  cov3[0][0]=pars[12];
  cov3[1][1]=pars[13];
  cov3[0][1]=pars[14];
  cov3[1][0]=cov2[0][1];
  TCovEllipse* el3=new TCovEllipse(cov3,pars[10],pars[11]);

  TF2* f4=new TF2("dis3_cut",&dis3_cut,-5,10,-5,10,16);
  f4->SetParameters(pars);

  TFile* file=new TFile("StdDis.root","RECREATE");
  el->Write("covel");
  el2->Write("covel2");
  el3->Write("covel3");
  f->Write();
  f2->Write();
  f3->Write();
  f4->Write();
  file->Close();
}

// ***********************************************************************
// *************** Riemann tracker ***************************************
void
PndTpcRecoTester::testRiemannMap(){

  std::vector<PndTpcCluster*> cll; // clusterlist
  maketoytracks(1,cll);

  std::vector<PndTpcRiemannHit*> rhits;
  int nhits=cll.size();
  for(int i=0;i<nhits;++i){
    rhits.push_back(new PndTpcRiemannHit(cll[i]));
  }
  
  TPolyMarker3D* line=new TPolyMarker3D(nhits);

  PndTpcRiemannTrack track;

  for(int j=0;j<nhits;++j){
    TVector3 x=rhits[j]->x();
    //x.Print();
    
    track.addHit(rhits[j]);

    TVector3 r=x-TVector3(0,0,0.5);
    //std::cout<<"R="<<r.Mag()<<std::endl;
    line->SetPoint(j,x.X(),x.Y(),x.Z());
  }

  track.av().Print();
  track.refit();
  track.szFit();
  track.n().Print();
  std::cout<<"abs(n)="<<TMath::Sqrt(track.n().Norm2Sqr())<<std::endl;
  std::cout<<"r="<<track.r()<<std::endl;
  std::cout<<"m="<<track.m()<<std::endl;
  std::cout<<"dip="<<track.dip()<<std::endl;
  std::cout<<"sign="<<track.sign()<<std::endl;
  std::cout<<"Center:";
  track.orig().Print();
  
  double c=track.c();
  TVectorD o(3);
  o[0]=rhits[0]->x().X();
  o[1]=rhits[0]->x().Y();
  o[2]=rhits[0]->x().Z();
  TVectorD e=o+track.n();

  TGraph* sz=new TGraph(nhits);
  for(int j=0;j<nhits;++j){
    rhits[j]->calcPosOnTrk(&track);
    sz->SetPoint(j,rhits[j]->s(),rhits[j]->z());
  }

  TPolyLine3D* arr=new TPolyLine3D(2);
  arr->SetPoint(0,o[0],o[1],o[2]);
  arr->SetPoint(1,e[0],e[1],e[2]);

  TFile* file=new TFile("Riemann.root","RECREATE");
  line->Write("rieTrack");
  arr->Write("arrow");
  sz->Write("sz");
  file->Close();
    

}

void
PndTpcRecoTester::testRiemannFit(unsigned int ntrk, double r){
  
  TH1D* rh=new TH1D("rh","rh",100,r-20,r+20);
  TH1D* diph=new TH1D("diph","diph",40,-0.1,0.1);

  double mean=0;

  PndTpcRiemannTrackFinder finder;
  finder.setMinHitsForFit(4);
  finder.addCorrelator(new PndTpcRiemannHTCorrelator(1E-3,0.01,2.0));
  finder.addCorrelator(new PndTpcProximityHTCorrelator(2.));
  

  for(int i=0;i<ntrk;++i){
    std::vector<PndTpcCluster*> cll; // clusterlist
    std::vector<double> rv;
    std::vector<TVector2> dirv;
    std::vector<double> dipv;
    std::vector<TVector2> startv;
    
    rv.push_back(r);
    dipv.push_back(10);
    
    maketoyevent(1,rv,dirv,dipv,startv);
    maketoyhits(1,rv,dirv,dipv,startv,cll);
    
    std::vector<PndTpcRiemannTrack*> cands;
        
    finder.buildTracks(cll,cands);
    
    int ncands=cands.size();
    for(int i=0;i<ncands;++i){
      PndTpcRiemannTrack* trk=cands.at(i);
      if(trk->r()>10000)continue;
      mean+=trk->r()*100;
      rh->Fill(trk->r()*100);
      diph->Fill(fabs(trk->m())-0.2);
    }
  }
  mean/=(double)ntrk;  
  std::cout<<"mean="<<mean<<std::endl;

  TFile* outf=new TFile("riemannfit.root","RECREATE");
  if(!outf->IsOpen()){std::cout<<"No file"<<std::endl;return;}
  rh->Write();
  diph->Write();
  outf->Close();
}


void
PndTpcRecoTester::testRiemannFinder(unsigned int ntrk){

  std::vector<PndTpcCluster*> cll; // clusterlist
  maketoytracks(ntrk,cll);
  std::vector<PndTpcRiemannTrack*> cands;

  PndTpcRiemannTrackFinder finder;
  finder.setMinHitsForFit(4);
  finder.addCorrelator(new PndTpcRiemannHTCorrelator(1E-3,0.01,2.0));
  finder.addCorrelator(new PndTpcProximityHTCorrelator(2.));
  finder.buildTracks(cll,cands);
  
  
  unsigned int ntrks=cands.size();
  //std::vector<TPolyMarker3D*> marker;
  TFile* file=new TFile("RiemannFinder.root","RECREATE");
  for(int i=0;i<ntrks;++i){
    PndTpcRiemannTrack* trk=cands.at(i);
    
    //abuse TLorentzVector:
    TLorentzVector* plane=new TLorentzVector(trk->n()[0],
					     trk->n()[1],
					     trk->n()[2],
					     trk->c());

    

    trk->szFit();
    unsigned int nhits=trk->getNumHits();
    TPolyMarker3D* maker=new TPolyMarker3D(nhits);
    TPolyMarker3D* clust=new TPolyMarker3D(nhits);
    TGraph* sz=new TGraph(nhits);
    for(int h=0;h<nhits;++h){
      const PndTpcRiemannHit* hit=trk->getHit(h);
      maker->SetPoint(h,hit->x().X(),hit->x().Y(),hit->x().Z());
      sz->SetPoint(h,hit->s(),hit->z());
      clust->SetPoint(h,hit->cluster()->pos().X()/30,hit->cluster()->pos().Y()/30,0);
    }
    char name[200];
    sprintf(name,"track%i",i);
    maker->Write(name);
    sprintf(name,"sztrack%i",i);
    sz->Write(name);
    sprintf(name,"clusters%i",i);
    clust->Write(name);
    sprintf(name,"plane%i",i);
    plane->Write(name);

  }
  file->Close(); 
  
}



// ***********************************************************************
// *************** Private methods ***************************************


void 
PndTpcRecoTester::maketoyevent(int n,std::vector<double>& rv,std::vector<TVector2>& dirv, std::vector<double>& dipv, std::vector<TVector2>& startv){
  bool rflag=false;
  bool dipflag=false;
  if(rv.size()!=0)rflag=true;
  if(dipv.size()!=0)dipflag=true;
  for(int i=0;i<n;++i){
    double r=gRandom->Uniform(30,70);
    double dirstartx;
    double dirstarty;
    gRandom->Circle(dirstartx,dirstarty,15);
    TVector2 dir(dirstartx,dirstarty);
    double dz_dphi=gRandom->Uniform(1.,100.);

    double startr=gRandom->Uniform(0,40);
    double startx; double starty;
    gRandom->Circle(startx,starty,startr);
    TVector2 start(startx,starty);
    std::cout<<"Track: r="<<r<<"  dz/dphi="<<dz_dphi
	     <<"   dz/ds="<<dz_dphi/r<<" dir=("
	     <<dir.X()<<","<<dir.Y()<<")  start=("
	     <<start.X()<<","<<start.Y()<<")"<<std::endl;
    if(!rflag)rv.push_back(r);
    dirv.push_back(dir);
    if(!dipflag)dipv.push_back(dz_dphi);
    startv.push_back(start);
  }
}


void
PndTpcRecoTester::maketoyhits(int ntrk,std::vector<double>& rv,
			   std::vector<TVector2>& dirv, 
			   std::vector<double>& dipv,
			   std::vector<TVector2>& startv,
			   std::vector<PndTpcCluster*>& cll){
  //produce hits
  TGraph* hits=new TGraph();
  TGraph* confhits=new TGraph();
  TGraph* hitsxz=new TGraph();
  TGraph* confhits_shift=new TGraph();

  unsigned int hitcounter=0;
  double sig=0.01; // sigma in space

  for(int t=0;t<ntrk;++t){// track loop
    
    unsigned int hitsintrk=0;

    double r=rv[t];
    std::cout<<"Radius r="<<r<<std::endl;
    TVector2 start=startv[t];
    TVector2 dir=dirv[t];
    TVector2 radial=dir.Unit().Rotate(TMath::Pi()/2)*r;
    TVector2 origin=start+radial;
    std::cout<<"Origin("<<origin.X()<<","<<origin.Y()<<")"<<std::endl;
    
    double dz_dphi=dipv[t];
    double ds2=0.25;
    double dphi2=ds2/(r*r+dz_dphi*dz_dphi); // a hit every 5 mm 
    double dphi=TMath::Sqrt(dphi2);
    double rot=gRandom->Uniform(-1,1)>0 ? 1. : -1.;
    std::cout<<"rot="<<rot<<std::endl;


    for(int i=1; i<60; ++i){// hit loop 1 degrees
      double phi=i*dphi;
      TVector2 point=origin-radial.Rotate(phi*rot);
      double xhit=point.X()+gRandom->Gaus(0,sig);
      double yhit=point.Y()+gRandom->Gaus(0,sig);
      double zhit=phi*dz_dphi+gRandom->Gaus(0,sig);
      double r=TMath::Sqrt(xhit*xhit+yhit*yhit);
      if(r<5) continue;
      hits->SetPoint(hitcounter,xhit,yhit);
      hitsxz->SetPoint(hitcounter,xhit,zhit);
      PndTpcCluster* cl=new PndTpcCluster(TVector3(xhit,yhit,zhit),1,hitcounter);
      PndTpcConfMapRecoHit conf(cl);
      conf.setRotated();
      cll.push_back(cl);
      confhits->SetPoint(hitcounter,conf.getXcf(),conf.getHitCoord(DetPlane())[0][0]);
      conf.setReferencePoint(start.X(),start.Y());
      confhits_shift->SetPoint(hitcounter,conf.getXcf(),conf.getHitCoord(DetPlane())[0][0]);
      ++hitcounter;
      ++hitsintrk;
    }
    //std::cout<<hitsintrk<<" hits in track"<<t<<std::endl;
  }
  TFile* outfile=new TFile("toyhits.root","RECREATE");
  hits->Write("hits");
  hitsxz->Write("hitsxz");
  confhits->Write("confhits");
  confhits_shift->Write("confhits2");
  outfile->Close();

}

void
PndTpcRecoTester::maketoytracks(int n,std::vector<PndTpcCluster*>& cll, bool vertex)
{
  // make events
  std::vector<double> rv;
  std::vector<TVector2> dirv;
  std::vector<double> dipv;
  std::vector<TVector2> startv;
  maketoyevent(n,rv,dirv,dipv,startv);
  if(vertex){
    for(int i=0;i<startv.size();++i)(startv[i]).Set(0.0,0.0);
  }
  maketoyhits(n,rv,dirv,dipv,startv,cll);
}

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