//-----------------------------------------------------------
// 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 "TpcHelMatchTaskPhilipp.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 "TH3.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TMultiGraph.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TMath.h"

#include "GFAbsRecoHit.h"
#include "GFTrack.h"
#include "GFTrackCand.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


//-----------------------------------------------------------
TpcHelMatchTaskPhilipp::TpcHelMatchTaskPhilipp()
  : fFopiEvBranchName("FopiEvent"),
    fTpcTrackFitBranchName("TpcTrackPostFit"),    
    fHelHitBranchName("HeliHit"),
    fHelTrackBranchName("HeliTrack"),
    fPlaTrackBranchName("PlawaTrack"),
    fPersistence(kTRUE)
{	

  fRad = 3.1415926/180.;
}

//-----------------------------------------------------------
TpcHelMatchTaskPhilipp::~TpcHelMatchTaskPhilipp(){
  delete fPar;
}


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

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

  // prepare TPC tree branch for output
  TPCdata = new Double_t[11];
  TPCdataVec = new TVectorD(11);
  fTPCdata = new TClonesArray("TVectorD");
  ioman->Register("TPCdata","TPCdata",fTPCdata,fPersistence);

  // prepare Helitron/PLAWA tree branch for output
  Heldata = new Double_t[12];
  HeldataVec = new TVectorD(12);
  fHeldata = new TClonesArray("TVectorD");
  ioman->Register("Heldata","Heldata",fHeldata,fPersistence);

  // prepare Helitron point tree branch for output
  Helpoint = new Double_t[3];
  HelpointVec = new TVectorD(3);
  fHelpoint = new TClonesArray("TVectorD");
  ioman->Register("Helpoint","Helpoint",fHelpoint,fPersistence);

  // prepare TPC/Helitron/PLAWA Matching tree branch for output
  MatchingPars = new Double_t[23];
  MatchingParsVec = new TVectorD(23);
  fMatchingPars = new TClonesArray("TVectorD");
  ioman->Register("TPCHelMatch","TPCHelMatch",fMatchingPars,fPersistence);


  // prepare combined-tracks tree branch
  CombHit = new Double_t[1000];
  CombHitVec = new TVectorD(1000);
  fCombHit = new TClonesArray("TVectorD");
  ioman->Register("CombHit","CombHit",fCombHit,fPersistence);

  return kSUCCESS;
}

//-----------------------------------------------------------
void TpcHelMatchTaskPhilipp::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");

  // Get Tpc digitisation parameter container
  fPar= (TpcDigiPar*) db->getContainer("TpcDigiPar");
  if (!fPar ) Fatal("SetParContainers", "TpcDigiPar not found");
	
}

//-----------------------------------------------------------
void TpcHelMatchTaskPhilipp::Exec(Option_t* opt) {


	Double_t ps[6] = {0.1,0.25,0.15,0.15,0.15,0.25};  // symbol size
	Double_t const phmatchmin = 0.0;
	Int_t const TpcId = 8;
	char ltxt[200];
	TVector3 fitPos;
	TVector3 tpcMasterPos,tpcMasterMom;
	
  // clean up result buffer
  if(fHeldata==0) 
    Fatal("PndTpcHelMatchingTask::Exec()","No fHeldata");
  fHeldata->Delete();

  if(fHelpoint==0) 
    Fatal("PndTpcHelMatchingTask::Exec()","No fHelpoint");
  fHelpoint->Delete();

  if(fTPCdata==0) 
    Fatal("PndTpcHelMatchingTask::Exec()","No fTPCdata");
  fTPCdata->Delete();

  if(fMatchingPars==0) 
    Fatal("PndTpcHelMatchingTask::Exec()","No fMatchingPars");
  fMatchingPars->Delete();

  if(fCombHit==0) 
    Fatal("TpcHelMatchTaskPhilipp::Exec()","No fCombHit");
  fCombHit->Delete();

  FairRunAna* run = FairRunAna::Instance();
  FOPIField *fMagField = (FOPIField*) run->GetField();
  Double_t Bz = fMagField->GetBz(0.,0.,0.);
  // fprintf(stderr,"Magfield: %f\n",Bz);
  
	// 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();

	printf("run: %i\n", nrun);
	printf("spill: %i\n", nspill);
	printf("event: %i\n", nevent);
  
  // transformation Helitron-to-beam 
  Double_t taroff = fopiEv->GetTaroff();
  Double_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 = fTpcTrackFitArray->GetEntriesFast();
  Int_t NumHelTracks = fopiEv->GetHmul();
  
	//////////////////////////////////////////
	//		TPC			//
	//////////////////////////////////////////

// prepare 
  Int_t *NumTpcTrackHits = new Int_t[NumTpcTracks];
  int tcounter=0;
  TVector3 pos0, mom0;
  TMatrixT<double> cov0;  
  TVector3 pos;

  Double_t tx[NumTpcTracks][200];	//TPC x of track hit
  Double_t ty[NumTpcTracks][200];	//TPC y of track hit
  Double_t tz[NumTpcTracks][200];	//TPC z of track hit
  Double_t tr[NumTpcTracks][200];	//TPC r of track hit

  Double_t *theTPC = new Double_t[3*NumTpcTracks];	//TPC track Theta
  Double_t *phiTPC = new Double_t[3*NumTpcTracks];	//TPC track Phi
  Double_t *momTPC = new Double_t[NumTpcTracks];	//TPC track Momentum
  Double_t *chaTPC = new Double_t[NumTpcTracks];	//TPC track Charge 

  GFDetPlane *plane0 = new GFDetPlane(TVector3(0.,0.,15.),
				      TVector3(1.,0.,0.),
				      TVector3(0.,1.,0.));
  
//analyse TPC tracks
  for (Int_t ttrack=0; ttrack<NumTpcTracks; ttrack++) {
	GFTrack* tpcsingletrack = (GFTrack*) fTpcTrackFitArray->At(ttrack);
	NumTpcTrackHits[ttrack] = tpcsingletrack->getNumHits();
  //get position of TPC-clusters for each track
	for(Int_t ii=0; ii<NumTpcTrackHits[ttrack]; ii++){
		GFAbsRecoHit *tsinglehit= tpcsingletrack->getHit(ii);
		TVectorD raw(3);
		raw = tsinglehit->getRawHitCoord();
		pos.SetXYZ(raw[0],raw[1],raw[2]);
		tx[ttrack][ii]=pos.X();
		ty[ttrack][ii]=pos.Y();
		tz[ttrack][ii]=pos.Z();
	}
  //angles from the track fit are evaluated at plane0
	momTPC[ttrack] = tpcsingletrack->getMom().Mag();
	chaTPC[ttrack] = tpcsingletrack->getCharge();
	// extrapolate track to given z-position (plane0)
	try {
	tpcsingletrack->getPosMomCov(*plane0,pos0,mom0,cov0);
	}
	catch (GFException& e){
	e.what();
	}
	theTPC[3*ttrack]=atan2(sqrt(pow(mom0.X(),2)+pow(mom0.Y(),2)),mom0.Z())/fRad;
	phiTPC[3*ttrack]=atan2(mom0.X(),mom0.Y())/fRad;
	if(phiTPC[3*ttrack]<0){phiTPC[3*ttrack]+=360;}
  // estimate theta and phi using hit points
    	if (NumTpcTrackHits[ttrack] > 1) {
	theTPC[3*ttrack+1]=0;
	phiTPC[3*ttrack+1]=0;
      	for (Int_t ii=0; ii<NumTpcTrackHits[ttrack]; ii++) {
        	GFAbsRecoHit* hit=tpcsingletrack->getHit(ii);
		TVectorD raw(3);
		raw = hit->getRawHitCoord();
        	pos.SetXYZ(raw[0],raw[1],raw[2]);
		theTPC[3*ttrack+1] += acos(tz[ttrack][ii]/sqrt(pow(tx[ttrack][ii],2)+pow(ty[ttrack][ii],2)+pow(tz[ttrack][ii],2)));
		phiTPC[3*ttrack+1] += atan2(ty[ttrack][ii],tx[ttrack][ii])/fRad;
	}
	theTPC[3*ttrack+1] = theTPC[3*ttrack+1]/NumTpcTrackHits[ttrack]/fRad;
	phiTPC[3*ttrack+1] = phiTPC[3*ttrack+1]/NumTpcTrackHits[ttrack];
	if(phiTPC[3*ttrack+1]<0){phiTPC[3*ttrack+1]+=360;}
	}
  //estimate theta and phi using hit pairs
	theTPC[3*ttrack+2]=0;
	phiTPC[3*ttrack+2]=0;
	tcounter=0;
	if(NumTpcTrackHits[ttrack]>1){
	   for(Int_t ii=0; ii<NumTpcTrackHits[ttrack]; ii++){
		for(int jj=ii+1; jj<NumTpcTrackHits[ttrack] ;jj++){
		//theta calculation
			theTPC[3*ttrack+2]+=acos((tz[ttrack][jj]-tz[ttrack][ii])/sqrt(pow(tx[ttrack][jj]-tx[ttrack][ii],2)+pow(ty[ttrack][jj]-ty[ttrack][ii],2)+pow(tz[ttrack][jj]-tz[ttrack][ii],2)));
		//phi calculation
			phiTPC[3*ttrack+2]+=atan2((ty[ttrack][jj]-ty[ttrack][ii]),
						   (tx[ttrack][jj]-tx[ttrack][ii]));
			tcounter++;
		}
	   }
   	   theTPC[3*ttrack+2]=theTPC[ttrack]/tcounter/fRad;
	   phiTPC[3*ttrack+2]=phiTPC[ttrack]/tcounter/fRad;
	   if(phiTPC[3*ttrack+2]<0){phiTPC[ttrack]+=360;}
	}
	else{theTPC[3*ttrack+2]=360;}
  }

	//////////////////////////////////////////
	//		Helitron		//
	//////////////////////////////////////////

//prepare
  Double_t hrsec;
  Int_t hhrc = -1;
  Int_t *NumHelTrackHits = new Int_t[NumHelTracks];

  Double_t hx[NumHelTracks][100];
  Double_t hy[NumHelTracks][100];
  Double_t hz[NumHelTracks][100];

  Double_t *px = new Double_t[NumHelTracks];
  Double_t *py = new Double_t[NumHelTracks];
  Double_t *pz = new Double_t[NumHelTracks];

  Double_t *theHel = new Double_t[2*NumHelTracks];	//Hel track Theta
  Double_t *phiHel = new Double_t[2*NumHelTracks];	//Hel track Phi
  Double_t *momHel = new Double_t[NumHelTracks];	//Hel track momentum
  Double_t *charHel = new Double_t[NumHelTracks];	//Hel track charge
  Double_t *massHel = new Double_t[NumHelTracks];	//Hel track mass
  Double_t *dedxHel = new Double_t[NumHelTracks];	//Hel track dEdx
  Double_t *phmatch = new Double_t[NumHelTracks];	//Quality of matched PLAWA hit


// analyse Helitron tracks	
  Int_t counter = -1;	
  for (Int_t htrack=0; htrack<NumHelTracks; htrack++) {
	HeliTrack *hsingletrack = (HeliTrack*) fHelTrackArray->At(htrack);
	PlawaTrack *psingletrack = (PlawaTrack*) fPlaTrackArray->At(htrack);
	NumHelTrackHits[htrack] = hsingletrack->GetHhmul();
  // momentum and Helitron/PLAWA matching
	theHel[2*htrack] = hsingletrack->GetHthe();
	phiHel[2*htrack] = hsingletrack->GetHphi();
	momHel[htrack] = hsingletrack->GetHmom();
	massHel[htrack] = hsingletrack->GetHmass();
	charHel[htrack] = hsingletrack->GetHchar();
	dedxHel[htrack] = hsingletrack->GetHdedx();	
	phmatch[htrack] = hsingletrack->GetPHmatch();
  // reference sector
	hrsec = hsingletrack->GetHrsec();
	lochit.SetHrsec(hrsec);
  // compute theta and phi from track hits
      	theHel[2*htrack+1]=0;
	phiHel[2*htrack+1]=0;
	for (Int_t ii=0; ii<NumHelTrackHits[htrack]; ii++) {
		hhrc++;
		HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(hhrc);
	// transform into global coordinate system
		lochit.SetPos(0,hsinglehit->GetHhx(),
				hsinglehit->GetHhy(),
				hsinglehit->GetHhz());
	//get position of Helitron tracks
		hx[htrack][ii]=lochit.GetGlobalPos().X();
		hy[htrack][ii]=lochit.GetGlobalPos().Y();
		hz[htrack][ii]=lochit.GetGlobalPos().Z();
//		hr[htrack][ii] = sqrt(pow(hx[htrack][ii],2)+pow(hy[htrack][ii],2));
	//angles
      		theHel[2*htrack+1] += acos(hz[htrack][ii]/sqrt(pow(hx[htrack][ii],2)+pow(hy[htrack][ii],2)+pow(hz[htrack][ii],2)));
		phiHel[2*htrack+1] += atan2(hy[htrack][ii],hx[htrack][ii]);
		Helpoint[0] = acos(hz[htrack][ii]/sqrt(pow(hx[htrack][ii],2)+pow(hy[htrack][ii],2)+pow(hz[htrack][ii],2)))/fRad;
		Helpoint[1] = atan2(hy[htrack][ii],hx[htrack][ii])/fRad;
		Helpoint[2] = hsingletrack->GetPHmatch();

		counter++;
		HelpointVec->SetElements(Helpoint);
		new((*fHelpoint)[counter]) TVectorD(*HelpointVec);
	}
    	theHel[2*htrack+1] = theHel[2*htrack+1]/NumHelTrackHits[htrack]/fRad;
    	phiHel[2*htrack+1] = phiHel[2*htrack+1]/NumHelTrackHits[htrack]/fRad;
	if(phiHel[2*htrack+1]<0){phiHel[2*htrack+1]+=360;}


  // get matching PLAWA hit
	if (phmatch[htrack] > phmatchmin) {
		Double_t phthe = fRad*psingletrack->GetPhthe();
		Double_t phphi = fRad*psingletrack->GetPhphi();
		pz[htrack] = psingletrack->GetPhz();
		px[htrack] = pz[htrack]*tan(phthe)*cos(phphi);
		py[htrack] = pz[htrack]*tan(phthe)*sin(phphi);
		Float_t pr = pz[htrack]*tan(phthe);
		globhit.SetHrsec(hrsec);
		globhit.SetPos(1,pr*cos(phphi),pr*sin(phphi),pz[htrack]);
	}
    
  }

	//////////////////////////////////////////////////////////
	//	Matching of TPC and Helitron Angles		//
	//////////////////////////////////////////////////////////
	Float_t theTHMatch[NumTpcTracks][NumHelTracks];
	for(Int_t ttrack=0; ttrack<NumTpcTracks; ttrack++){
		for(Int_t htrack=0; htrack<NumHelTracks;htrack++){
		   theTHMatch[ttrack][htrack]=0;
		   for(Int_t kk=0; kk<NumHelTrackHits[htrack]; kk++){
			theTHMatch[ttrack][htrack]+=acos((hz[htrack][kk]-tz[ttrack][0])/sqrt(pow(hx[htrack][kk]-tx[ttrack][0],2)+pow(hy[htrack][kk]-ty[ttrack][0],2)+pow(hz[htrack][kk]-tz[ttrack][0],2)));
		   }
		   theTHMatch[ttrack][htrack]=theTHMatch[ttrack][htrack]/NumHelTrackHits[htrack]/fRad;
		}
	}


	//////////////////////////////////////////////////////////
	//	Combining of TPC and Helitron Hits		//
	//////////////////////////////////////////////////////////

  Int_t hitcounter = -1;
  counter= -1;
 
  float thehistmean= -0.5835;
  float themeanerr = 0.0280;
  float thehistsigma = 1.954;
  float phihistmean = 10.42;
  float phimeanerr = 0.07;
  float phihistsigma = 2.822;


  for (Int_t ttrack=0; ttrack<NumTpcTracks; ttrack++){
	GFTrack* tpcsingletrack = (GFTrack*) fTpcTrackFitArray->At(ttrack);
  	hhrc= -1;
	if(theTPC[3*ttrack+2]<0){
	for (Int_t htrack=0; htrack<NumHelTracks; htrack++){
		HeliTrack *hsingletrack = (HeliTrack*) fHelTrackArray->At(htrack);
		PlawaTrack *psingletrack = (PlawaTrack*) fPlaTrackArray->At(htrack); 
		counter++;
  		hitcounter = -1;
		CombHit[0]=nrun; 		//run num
		CombHit[1]=nspill;		//spill num
		CombHit[2]=nevent;		//event num
		CombHit[3]= ttrack;		//TPC track number
		CombHit[4]= htrack;		//Hel track number
		CombHit[5]= NumTpcTracks;	//total number of TPC-tracks
		CombHit[6]= NumHelTracks;	//total number of Hel-tracks
		if (phmatch[htrack] > phmatchmin) {
		  CombHit[7]= 3*NumTpcTrackHits[ttrack] + 3*NumHelTrackHits[htrack] + 3;//num trackhits with matched PLAWA hit
		}
		else{CombHit[7]= 3*NumTpcTrackHits[ttrack] + 3*NumHelTrackHits[htrack];}//num trackhits w/o matched PLAWA hit
		CombHit[8]= phmatch[htrack];	//matched PLAWA hit
		CombHit[9] = theTPC[3*ttrack+2];	//Theta of TPC track cluster pairs
		CombHit[10] = phiTPC[3*ttrack+2];	//phi of TPC track cluster pairs
		CombHit[11]= theHel[2*htrack];		//Theta of Hel track fit
		CombHit[12]= phiHel[2*htrack];		//phi of Hel track fit
		CombHit[13]= momHel[htrack];		//momentum of Hel track
		//TPC clusters
		for(Int_t ii=0; ii<NumTpcTrackHits[ttrack]; ii++){
			hitcounter++;
			GFAbsRecoHit *tsinglehit= tpcsingletrack->getHit(ii);
			TVectorD raw(3);
			raw = tsinglehit->getRawHitCoord();
			pos.SetXYZ(raw[0],raw[1],raw[2]);
			CombHit[3*hitcounter+14]=pos.X();
			CombHit[3*hitcounter+15]=pos.Y();
			CombHit[3*hitcounter+16]=pos.Z();
		}
		//HeliHits
		hrsec = hsingletrack->GetHrsec();
		lochit.SetHrsec(hrsec);
		for(Int_t jj=0; jj<NumHelTrackHits[htrack]; jj++){
			hitcounter++;
			hhrc++;
			HeliHit *hsinglehit = (HeliHit*) fHelHitArray->At(hhrc);
		// transform into global coordinate system
			lochit.SetPos(0,hsinglehit->GetHhx(),
					hsinglehit->GetHhy(),
					hsinglehit->GetHhz());
		//get position of Helitron tracks
			CombHit[3*hitcounter+14]=lochit.GetGlobalPos().X();
			CombHit[3*hitcounter+15]=lochit.GetGlobalPos().Y();
			CombHit[3*hitcounter+16]=lochit.GetGlobalPos().Z();
		}
		if (phmatch[htrack] > phmatchmin) {
			hitcounter++;
			CombHit[3*hitcounter+14]=px[htrack];
			CombHit[3*hitcounter+15]=py[htrack];
			CombHit[3*hitcounter+16]=pz[htrack];
		}
	   CombHitVec->SetElements(CombHit);
	   new((*fCombHit)[counter]) TVectorD(*CombHitVec);
	}
	}
  }

//TPcHelMatchPars
  counter = -1;
  for (Int_t ttrack=0; ttrack<NumTpcTracks; ttrack++) {
    MatchingPars[0] = nrun;			//run number
    MatchingPars[1] = nspill;			//spill number
    MatchingPars[2] = nevent;			//event number
    MatchingPars[3] = ttrack;			//TPC track number
    MatchingPars[4] = theTPC[3*ttrack];		//Theta of TPC track fit
    MatchingPars[5] = theTPC[3*ttrack+1];	//Theta of TPC track clusters (with vertex)
    MatchingPars[6] = theTPC[3*ttrack+2];	//Theta of TPC track cluster pairs (w/o vertex)
    MatchingPars[7] = phiTPC[3*ttrack];		//Phi of TPC track fit
    MatchingPars[8] = phiTPC[3*ttrack+1];	//Phi of TPC track clusters (with vertex)
    MatchingPars[9] = phiTPC[3*ttrack+2];	//Phi of TPC track cluster pairs (w/o vertex)
    MatchingPars[10] = momTPC[ttrack];		//Momentum of TPC tracks
    MatchingPars[11] = NumTpcTrackHits[ttrack];	//Number of clusters per TPC track
    for (Int_t htrack=0; htrack<NumHelTracks; htrack++) {
      MatchingPars[12] = htrack;		//Hel track number
      MatchingPars[13] = theHel[2*htrack];	//Hel track Theta from fit
      MatchingPars[14] = theHel[2*htrack+1];	//Hel track Theta from track points
      MatchingPars[15] = phiHel[2*htrack];	//Hel track Phi from fit
      MatchingPars[16] = phiHel[2*htrack+1];	//Hel track Phi from track points
      MatchingPars[17] = theTHMatch[ttrack][htrack];	//Theta of Helitron tracks 
      MatchingPars[18] = momHel[htrack];	//Momentum of Helitron tracks
      MatchingPars[19] = massHel[htrack];	//Mass of Helitron tracks
      MatchingPars[20] = charHel[htrack];	//Char of Helitron tracks
      MatchingPars[21] = dedxHel[htrack];	//dEdx of Helitron tracks
      MatchingPars[22] = phmatch[htrack];	//Matched Plawa hit
      
      counter++;
      MatchingParsVec->SetElements(MatchingPars);
      new((*fMatchingPars)[counter]) TVectorD(*MatchingParsVec);
    }
  }
//TPCdata
  counter = -1;
  for (Int_t ttrack=0; ttrack<NumTpcTracks; ttrack++) {
    TPCdata[0] = nrun;			//run number
    TPCdata[1] = nspill;		//spill number
    TPCdata[2] = nevent;		//event number
    TPCdata[3] = theTPC[3*ttrack];	//Theta of TPC track fit
    TPCdata[4] = theTPC[3*ttrack+1];	//Theta of TPC track clusters (with vertex)
    TPCdata[5] = theTPC[3*ttrack+2];	//Theta of TPC track cluster pairs (w/o vertex)
    TPCdata[6] = phiTPC[3*ttrack];	//Phi of TPC track fit
    TPCdata[7] = phiTPC[3*ttrack+1];	//Phi of TPC track clusters (with vertex)
    TPCdata[8] = phiTPC[3*ttrack+2];	//Phi of TPC track cluster pairs (w/o vertex)
    TPCdata[9] = momTPC[ttrack];	//Momentum of TPC tracks
    TPCdata[10]=NumTpcTrackHits[ttrack];//Number of clusters per TPC track
    
    counter++;
    TPCdataVec->SetElements(TPCdata);
    new((*fTPCdata)[counter]) TVectorD(*TPCdataVec);
  }
//Heldata
  counter = -1;
  for (Int_t htrack=0; htrack<NumHelTracks; htrack++) {
    Heldata[0] = nrun;			//run number
    Heldata[1] = nspill;		//spill number
    Heldata[2] = nevent;		//event number
    Heldata[3] = theHel[2*htrack];	//Hel track Theta from fit
    Heldata[4] = theHel[2*htrack+1];	//Hel track Theta from track points
    Heldata[5] = phiHel[2*htrack];	//Hel track Phi from fit
    Heldata[6] = phiHel[2*htrack+1];	//Hel track Phi from track points
    Heldata[7] = momHel[htrack];	//Momentum of Helitron tracks
    Heldata[8] = massHel[htrack];	//Mass of Helitron tracks
    Heldata[9] = charHel[htrack];	//Char of Helitron tracks
    Heldata[10] = dedxHel[htrack];	//dEdx of Helitron tracks
    Heldata[11] = phmatch[htrack];	//Matched Plawa hit
    
    counter++;
    HeldataVec->SetElements(Heldata);
    new((*fHeldata)[counter]) TVectorD(*HeldataVec);
  }

 
}

ClassImp(TpcHelMatchTaskPhilipp)

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