[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCMerger.cxx

//$Id$
// Original: AliHLTMerger.cxx,v 1.16 2005/06/14 10:55:21 cvetan 

/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        * 
 * ALICE Experiment at CERN, All rights reserved.                         *
 *                                                                        *
 * Primary Authors: Uli Frankenfeld, maintained by                          *
 *                  Matthias Richter <Matthias.Richter@ift.uib.no>        *
 *                  for The ALICE HLT Project.                            *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/** @file   AliHLTTPCMerger.cxx
    @author Uli Frankenfeld, maintained by Matthias Richter
    @date   
    @brief  The HLT TPC merger base class
*/

#include "AliHLTTPCLogging.h"
#include "AliHLTTPCMerger.h"
#include "AliHLTTPCTrack.h"
#include "AliHLTTPCTrackSegmentData.h"
#include "AliHLTTPCTransform.h"
#include "AliHLTTPCTrackArray.h"

#ifdef use_root //use root ntuple for slow merge
#include <TNtuple.h>
#include <TTree.h>
#include <TFile.h>
#endif

#if __GNUC__ >= 3
using namespace std;
#endif

ClassImp(AliHLTTPCMerger)

AliHLTTPCMerger::AliHLTTPCMerger()
  :
  fCurrentTracks(0),
  fSlice(0), // TODO its not quite clear were this member is set
  fVertex(NULL),
  f2Global(false), // TODO: initialization was missing at all, assuming false for the moment
  fMaxY(0.0),
  fMaxZ(0.0),
  fMaxKappa(0.0),
  fMaxPsi(0.0),
  fMaxTgl(0.0),
  fTrackType(0),
  fInTrack(NULL),
  fNIn(0),
  fOutTrack(NULL)
{
  //Default constructor
  fInTrack=0;
  fOutTrack=0;
  fCurrentTracks=0;
  fNIn=0;
}

AliHLTTPCMerger::~AliHLTTPCMerger()
{
  //Destructor
  DeleteArray();
}

void AliHLTTPCMerger::InitMerger(Int_t ntrackarrays,const Char_t *tracktype)
{
  //Used to setup all arrays
  
  if(strcmp(tracktype,"AliHLTTPCTrack")==0) fTrackType='t';
  else if(strcmp(tracktype,"AliHLTTPCConfMapTrack")==0) fTrackType='c';
#ifdef INCLUDE_TPC_HOUGH
  else if(strcmp(tracktype,"AliHLTTPCHoughTrack")==0) fTrackType='h';
#endif
  else
    LOG(AliHLTTPCLog::kError,"AliHLTTPCMerger::AliHLTTPCMerger","Track types")
      <<"Unknown tracktype"<<ENDLOG;
  SetArray(ntrackarrays);
  fCurrentTracks=0;

}

void AliHLTTPCMerger::DeleteArray()
{
  //delete arrays
  for(Int_t i=0; i<fNIn;i++)
    {
      if(!fInTrack[i]) continue;
      delete fInTrack[i];
      fInTrack[i]=0;
    }
  if(fInTrack)
    delete[] fInTrack;
  if(fOutTrack)
    delete fOutTrack;
  fInTrack=0;
  fOutTrack=0;
}

void AliHLTTPCMerger::SetArray(Int_t nin)
{ 
  //set arrays
  DeleteArray();//Make sure arrays are cleaned 
  
  fNIn = nin;
  fInTrack = new AliHLTTPCTrackArray*[fNIn];
  for(Int_t i=0; i<fNIn;i++)
    {
#ifdef INCLUDE_TPC_HOUGH
      if(fTrackType=='h')
	fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
      else
#endif
	fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCTrack");
      
    }
#ifdef INCLUDE_TPC_HOUGH
  if(fTrackType=='h')
    fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
  else
#endif
    fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCTrack");
}

void AliHLTTPCMerger::Reset()
{ 
  //reset
  for(Int_t i=0; i<fNIn;i++)
    {
      fInTrack[i]->Reset();
    }
  fOutTrack->Reset();
}

void AliHLTTPCMerger::FillTracks(Int_t ntracks, AliHLTTPCTrackSegmentData* tr)
{
  //Read tracks from shared memory (or memory)
  AliHLTTPCTrackArray *destination = GetInTracks(fCurrentTracks);
  if(Is2Global())
    destination->FillTracksChecked(tr, ntracks, 0, fSlice);
  else
    destination->FillTracksChecked(tr, ntracks, 0);
}

void AliHLTTPCMerger::AddAllTracks()
{ 
  //add all tracks
  for(Int_t i=0; i<GetNIn();i++)
    {
      AliHLTTPCTrackArray *in = GetInTracks(i);
      AliHLTTPCTrackArray *out = GetOutTracks();
      out->AddTracks(in);
    }
}

void AliHLTTPCMerger::SortGlobalTracks(AliHLTTPCTrack **tracks, Int_t ntrack)
{ 
  //sort global tracks
  AliHLTTPCTrack **tmp = new AliHLTTPCTrack*[ntrack]; 
  for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
  Int_t *t = new Int_t[ntrack];
  for(Int_t i=0;i<ntrack;i++) t[i]=-1;
  
  for(Int_t j=0;j<ntrack;j++)
    {
      Double_t minr=300;
      Int_t    mini=0;
      for(Int_t i=0;i<ntrack;i++)
	{
	  if(!tracks[i]) continue;
	  Double_t rr=pow(tracks[i]->GetFirstPointX(),2)+pow(tracks[i]->GetFirstPointY(),2);
	  Double_t r=sqrt(rr);
	  if(r<minr){
	    minr=r;
	    mini=i;
	  }
	}
      t[j]=mini;
      tracks[mini]=0;
    }
  for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
  delete[] t;
  delete[] tmp;
}

void AliHLTTPCMerger::SortTracks(AliHLTTPCTrack **tracks, Int_t ntrack) const
{ 
  //sort tracks
  AliHLTTPCTrack **tmp = new  AliHLTTPCTrack*[ntrack];
  for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
  Int_t *t = new Int_t[ntrack];
  for(Int_t i=0;i<ntrack;i++) t[i]=-1;
  
  for(Int_t j=0;j<ntrack;j++)
    {
      Double_t minx=300; 
      Int_t    mini=0;
      for(Int_t i=0;i<ntrack;i++)
	{
	  if(!tracks[i]) continue;
	  if(tracks[i]->GetFirstPointX()<minx)
	    {
	      minx=tracks[i]->GetFirstPointX();
	      mini=i;
	    }     
	}
      t[j]=mini;  
      tracks[mini]=0;
    }
  for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
  delete[] t;
  delete[] tmp;
}

void AliHLTTPCMerger::AddTrack(AliHLTTPCTrackArray *mergedtrack,AliHLTTPCTrack *track)
{ 
  // add tracks
  AliHLTTPCTrack *t[1];
  t[0] = track;
  MultiMerge(mergedtrack,t,1);
}

AliHLTTPCTrack * AliHLTTPCMerger::MergeTracks(AliHLTTPCTrackArray *mergedtrack,AliHLTTPCTrack *t0,AliHLTTPCTrack *t1)
{ 
  //merge tracks
  AliHLTTPCTrack *t[2];
  t[0] = t0; 
  t[1] = t1;
  SortTracks(t,2);
  return MultiMerge(mergedtrack,t,2);
}

AliHLTTPCTrack * AliHLTTPCMerger::MultiMerge(AliHLTTPCTrackArray *mergedtracks,AliHLTTPCTrack **tracks, Int_t ntrack)
{
  //multi merge the tracks
  //check npoints
  Int_t nps = 0;
  for(Int_t i=0;i<ntrack;i++)
    {
      nps+=tracks[i]->GetNHits();
    }
  if(nps>AliHLTTPCTransform::GetNRows())
    {
      LOG(AliHLTTPCLog::kWarning,"AliHLTTPCMerger::MultiMerge","Adding Points")
	<<AliHLTTPCLog::kDec<<"Too many Points: "<<nps<<ENDLOG;
      return 0;
    }
  
  //create new track
  AliHLTTPCTrack *newtrack = mergedtracks->NextTrack();
  //copy points
  UInt_t * nn = new UInt_t[AliHLTTPCTransform::GetNRows()];
  nps = 0;
  
  //  for(Int_t i=0;i<ntrack;i++){
  for(Int_t i=ntrack-1;i>=0;i--)
    {
      memcpy(&nn[nps],tracks[i]->GetHitNumbers(),tracks[i]->GetNHits()*sizeof(UInt_t));
      nps+=tracks[i]->GetNHits();
    }
  AliHLTTPCTrack *tpf=tracks[0];
  AliHLTTPCTrack *tpl=tracks[ntrack-1];
  AliHLTTPCTrack *best = tpf;
  if(tpf->GetNHits()<tpl->GetNHits() && Is2Global())
    best = tpl;//Best means = most points and therefore best fit (in global case)
  
  newtrack->SetNHits(nps);
  newtrack->SetHits(nps,nn);
  newtrack->SetFirstPoint(tpf->GetFirstPointX(),tpf->GetFirstPointY(),tpf->GetFirstPointZ());
  newtrack->SetLastPoint(tpl->GetLastPointX(),tpl->GetLastPointY(),tpl->GetLastPointZ());
  newtrack->SetPt(best->GetPt());
  newtrack->SetPsi(best->GetPsi());
  newtrack->SetTgl(best->GetTgl());
  newtrack->SetCharge(tpf->GetCharge());
  delete [] nn;
  return newtrack;
}

void* AliHLTTPCMerger::GetNtuple(char *varlist) const
{ 
  //get ntuple
#ifdef use_root
  TNtuple* nt = new TNtuple("ntuple","ntuple",varlist);
  return (void*) nt;
#else
  return 0;
#endif
}

void* AliHLTTPCMerger::GetNtuple() const
{ 
  //get ntuple
#ifdef use_root
  TNtuple* nt = new TNtuple("ntuple","ntuple",
			    "dx:dy:dz:dk:dpsi:dtgl:dq:disx:disy:disz:dis:n0:n1:diff:drx:dry:drz");
  return (void*) nt;
#else
  return 0;
#endif
}

Bool_t AliHLTTPCMerger::WriteNtuple(char *filename, void* nt) const
{ 
  //write ntuple
#ifdef use_root
  TNtuple *ntuple=(TNtuple *) nt;
  TFile *f = new TFile(filename,"RECREATE");
  ntuple->Write();
  f->Close();
  delete ntuple; 
  return kTRUE; 
#else
  return kFALSE;
#endif
}

void AliHLTTPCMerger::FillNtuple(void *nt,AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
{ 
  //fill ntuple
  Float_t data[17];
  if(outertrack->IsPoint()&&innertrack->IsPoint())
    {
      data[0] =Float_t(innertrack->GetPointX()-outertrack->GetPointX());
      data[1] =Float_t(innertrack->GetPointY()-outertrack->GetPointY());
      data[2] =Float_t(innertrack->GetPointZ()-outertrack->GetPointZ());
      data[3] =Float_t(innertrack->GetKappa()-outertrack->GetKappa());
      Double_t psi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
      if(psi>AliHLTTPCTransform::Pi()) psi-=AliHLTTPCTransform::TwoPi();
      else if(psi<-AliHLTTPCTransform::Pi()) psi+=AliHLTTPCTransform::TwoPi();
      data[4] =Float_t(psi);
      data[5] =Float_t(innertrack->GetTgl()-outertrack->GetTgl());
      data[6] =Float_t(innertrack->GetCharge()-outertrack->GetCharge());
      data[7] =Float_t(innertrack->GetLastPointX()-outertrack->GetFirstPointX());
      data[8] =Float_t(innertrack->GetLastPointY()-outertrack->GetFirstPointY());
      data[9] =Float_t(innertrack->GetLastPointZ()-outertrack->GetFirstPointZ());
      data[10] =sqrt(pow(data[7],2)+pow(data[8],2)+pow(data[9],2));
      data[11]= outertrack->GetNHits();
      data[12]= innertrack->GetNHits();
      data[13] = Float_t(TrackDiff(innertrack,outertrack));
      data[14]=0;
      data[15]=0;
      data[16]=0;
#ifdef use_root
      TNtuple *ntuple = (TNtuple *) nt;
      ntuple->Fill(data);
#endif
    }
}

void AliHLTTPCMerger::FillNtuple(void *nt,Float_t *data) const
{ 
  //fill ntuple
#ifdef use_root
  TNtuple *ntuple = (TNtuple *) nt;
  ntuple->Fill(data);
#endif
}

Double_t AliHLTTPCMerger::GetAngle(Double_t a1,Double_t a2)
{ 
  //get angle
  Double_t da = a1 - a2 + 4*AliHLTTPCTransform::Pi();
  da = fmod(da,AliHLTTPCTransform::TwoPi());
  if(da>AliHLTTPCTransform::Pi()) da = AliHLTTPCTransform::TwoPi()-da;
  return da;
}

void AliHLTTPCMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxkappa, Double_t maxpsi, Double_t maxtgl)
{ 
  //set parameters for merger
  fMaxY = maxy;
  fMaxZ = maxz;
  fMaxKappa = maxkappa;
  fMaxPsi = maxpsi;
  fMaxTgl = maxtgl;
}

Bool_t AliHLTTPCMerger::IsTrack(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
{
  //is track to be merged
  if(innertrack->GetCharge()!=outertrack->GetCharge()) return kFALSE;
  if( (!innertrack->IsPoint()) || (!outertrack->IsPoint()) )  return kFALSE; 
  if(innertrack->GetNHits()+outertrack->GetNHits()>AliHLTTPCTransform::GetNRows()) return kFALSE;
  
  if(fabs(innertrack->GetPointY()-outertrack->GetPointY()) >fMaxY) return kFALSE;
  if(fabs(innertrack->GetPointZ()-outertrack->GetPointZ()) >fMaxZ) return kFALSE;
  if(fabs(innertrack->GetKappa()-outertrack->GetKappa())   >fMaxKappa) return kFALSE;
  if(GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi()) >fMaxPsi) return kFALSE;
  if(fabs(innertrack->GetTgl()-outertrack->GetTgl()) >fMaxTgl) return kFALSE;
  //if no rejection up to this point: merge!!
  return kTRUE;
}

Bool_t AliHLTTPCMerger::IsRTrack(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
{ 
  //same as IsTrack
  return IsTrack(innertrack,outertrack);
}

Double_t AliHLTTPCMerger::TrackDiff(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
{ 
  //return track difference
  Double_t diff =-1;
  Double_t x[4],y[4],z[4],dy[4],dz[4];
  AliHLTTPCTrack *tracks[2]; 
  
  tracks[0] = innertrack;
  tracks[1] = outertrack;
  SortGlobalTracks(tracks,2);
  innertrack = tracks[0]; 
  outertrack = tracks[1];
  
  x[0] = innertrack->GetFirstPointX();
  x[1] = innertrack->GetLastPointX();
  x[2] = outertrack->GetFirstPointX();
  x[3] = outertrack->GetLastPointX();
  
  y[0] = innertrack->GetFirstPointY();
  y[1] = innertrack->GetLastPointY();
  y[2] = outertrack->GetFirstPointY();
  y[3] = outertrack->GetLastPointY();

  z[0] = innertrack->GetFirstPointZ();
  z[1] = innertrack->GetLastPointZ();
  z[2] = outertrack->GetFirstPointZ();
  z[3] = outertrack->GetLastPointZ();

  
  outertrack->CalculatePoint(x[0]);
  if(!outertrack->IsPoint()) return diff;
  dy[0] = fabs(y[0] - outertrack->GetPointY());
  dz[0] = fabs(z[0] - outertrack->GetPointZ());
  
  outertrack->CalculatePoint(x[1]);
  if(!outertrack->IsPoint()) return diff;
  dy[1] = fabs(y[1] - outertrack->GetPointY());
  dz[1] = fabs(z[1] - outertrack->GetPointZ());
  
  innertrack->CalculatePoint(x[2]);
  if(!innertrack->IsPoint()) return diff;
  dy[2] = fabs(y[2] - innertrack->GetPointY());
  dz[2] = fabs(z[2] - innertrack->GetPointZ());
  
  innertrack->CalculatePoint(x[3]);
  if(!innertrack->IsPoint()) return diff;
  dy[3] = fabs(y[3] - innertrack->GetPointY());
  dz[3] = fabs(z[3] - innertrack->GetPointZ());

  diff=0;
  for(Int_t i=0;i<4;i++)
    diff+=sqrt(dy[i]*dy[i]+dz[i]*dz[i]);
  return diff; 
}

void AliHLTTPCMerger::PrintDiff(AliHLTTPCTrack *innertrack,AliHLTTPCTrack *outertrack)
{ 
  // print difference
  if(!innertrack->IsPoint()||!outertrack->IsPoint())
    {
      LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","No Points")<<ENDLOG;
      //cerr<<"AliHLTTPCMerger::PrintDiff: No Points"<<endl;
      //cerr<<"---------------------------"<<endl;
      return;
    } 
  
  Double_t dx = innertrack->GetPointX()-outertrack->GetPointX();
  Double_t dy = innertrack->GetPointY()-outertrack->GetPointY();
  Double_t dz = innertrack->GetPointZ()-outertrack->GetPointZ();
  Double_t dk = innertrack->GetKappa()-outertrack->GetKappa();
  Double_t dpsi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
  if(dpsi>AliHLTTPCTransform::Pi()) dpsi-=AliHLTTPCTransform::TwoPi();
  else if(dpsi<-AliHLTTPCTransform::Pi())dpsi+=AliHLTTPCTransform::TwoPi();
  //Double_t dpsi = GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi());
  Double_t dtgl= innertrack->GetTgl()-outertrack->GetTgl();
  Double_t dq =innertrack->GetCharge()-outertrack->GetCharge();
  
  LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","Points") <<"dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<" dk: "<<dk<<" dpsi: "<<dpsi<<" dtgl: "<<dtgl<<" dq: "<<dq<<ENDLOG;
  //fprintf(stderr,"dx: %4f dy: %4f dz: %4f dk: %4f dpsi: %4f dtgl: %4f dq: %4f\n",dx,dy,dz,dk,dpsi,dtgl,dq);
  //cerr<<"---------------------------"<<endl;
  
}

void AliHLTTPCMerger::Print()
{
  // print some infos
  for(Int_t i=0; i<fNIn; i++)
    {
      AliHLTTPCTrackArray *ttt= GetInTracks(i);
      for(Int_t j =0;j<ttt->GetNTracks();j++)
	{
	  AliHLTTPCTrack *track=ttt->GetCheckedTrack(j);
	  if(!track) continue;
	  track->CalculateHelix();
	  //      Double_t angle = atan2(track->GetLastPointY(),track->GetLastPointX());
	  //      if(angle<0) angle+=AliHLTTPCTransform::Pi();
	  if(track->CalculatePoint(135))
	    //      if(!track->CalculateEdgePoint(angle)) cerr<<"**************"<<endl;     
	    //      if(track->CalculatePoint(track->GetLastPointX()))
	    //      if(track->CalculatePoint(0))
	    {
	      //      PrintTrack(track);
	      //      track->CalculateReferencePoint(AliHLTTPCTransform::Pi()/180.);
	      track->CalculateReferencePoint(0.001);
	      Float_t dx=(float)track->GetPointX()-track->GetPointX();
	      Float_t dy=(float)track->GetPointY()-track->GetPointY();
	      Float_t dz=(float)track->GetPointZ()-track->GetPointZ();
	      LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::Print","RefPoint") <<"npt: "<<track->GetNHits()<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<ENDLOG;
	      
	      //fprintf(stderr,"npt: %3d dx: %8.5f dy: %8.5f dz: %8.5f\n",track->GetNHits(),dx,dy,dz);
	      //cerr<<"---------------------------"<<endl;
	    }
	}  
    }
}

void AliHLTTPCMerger::PrintTrack(AliHLTTPCTrack *track)
{ 
  //print track info
  fprintf(stderr,"npt: %3d pt: %.2f psi: %.2f tgl: %5.2f q: %2d\n",
	  track->GetNHits(),track->GetPt(),track->GetPsi(),
	  track->GetTgl(),track->GetCharge());
  fprintf(stderr,
	  "x1: %6.2f y1: %6.2f z1: %6.2f xl: %6.2f yl: %6.2f zl: %6.2f\n",
	  track->GetFirstPointX(),track->GetFirstPointY(),track->GetFirstPointZ(),
	  track->GetLastPointX(),track->GetLastPointY(),track->GetLastPointZ());
  if(track->IsPoint())
    {
      fprintf(stderr,
	      "R: %.2f Xc: %.2f Yc: %.2f Xp: %.2f Yp: %.2f Zp: %.2f Psip: %.2f\n",
	      track->GetRadius(),track->GetCenterX(),track->GetCenterY(),
	      track->GetPointX(),track->GetPointY(),track->GetPointZ(),
	      track->GetPointPsi());
    }
}
Commit	Line	Data
a6c02c85	1	//$Id$
4aa41877	2	// Original: AliHLTMerger.cxx,v 1.16 2005/06/14 10:55:21 cvetan
a6c02c85	3
2a083ac4	4	/**************************************************************************
9be2600f	5	* This file is property of and copyright by the ALICE HLT Project *
9be2600f	6	* ALICE Experiment at CERN, All rights reserved. *
2a083ac4	7	* *
9be2600f	8	* Primary Authors: Uli Frankenfeld, maintained by *
	9	* Matthias Richter <Matthias.Richter@ift.uib.no> *
	10	* for The ALICE HLT Project. *
2a083ac4	11	* *
	12	* Permission to use, copy, modify and distribute this software and its *
	13	* documentation strictly for non-commercial purposes is hereby granted *
	14	* without fee, provided that the above copyright notice appears in all *
	15	* copies and that both the copyright notice and this permission notice *
	16	* appear in the supporting documentation. The authors make no claims *
	17	* about the suitability of this software for any purpose. It is *
	18	* provided "as is" without express or implied warranty. *
	19	**************************************************************************/
	20
	21	/** @file AliHLTTPCMerger.cxx
	22	@author Uli Frankenfeld, maintained by Matthias Richter
	23	@date
	24	@brief The HLT TPC merger base class
	25	*/
a6c02c85	26
a6c02c85	27	#include "AliHLTTPCLogging.h"
	28	#include "AliHLTTPCMerger.h"
	29	#include "AliHLTTPCTrack.h"
	30	#include "AliHLTTPCTrackSegmentData.h"
	31	#include "AliHLTTPCTransform.h"
	32	#include "AliHLTTPCTrackArray.h"
	33
	34	#ifdef use_root //use root ntuple for slow merge
	35	#include <TNtuple.h>
	36	#include <TTree.h>
	37	#include <TFile.h>
	38	#endif
	39
a6c02c85	40	#if __GNUC__ >= 3
	41	using namespace std;
	42	#endif
	43
	44	ClassImp(AliHLTTPCMerger)
	45
	46	AliHLTTPCMerger::AliHLTTPCMerger()
1204d026	47	:
	48	fCurrentTracks(0),
	49	fSlice(0), // TODO its not quite clear were this member is set
	50	fVertex(NULL),
	51	f2Global(false), // TODO: initialization was missing at all, assuming false for the moment
	52	fMaxY(0.0),
	53	fMaxZ(0.0),
	54	fMaxKappa(0.0),
	55	fMaxPsi(0.0),
	56	fMaxTgl(0.0),
	57	fTrackType(0),
	58	fInTrack(NULL),
	59	fNIn(0),
	60	fOutTrack(NULL)
a6c02c85	61	{
	62	//Default constructor
	63	fInTrack=0;
	64	fOutTrack=0;
	65	fCurrentTracks=0;
	66	fNIn=0;
	67	}
	68
	69	AliHLTTPCMerger::~AliHLTTPCMerger()
	70	{
	71	//Destructor
	72	DeleteArray();
	73	}
	74
a6e0ebfe	75	void AliHLTTPCMerger::InitMerger(Int_t ntrackarrays,const Char_t *tracktype)
a6c02c85	76	{
	77	//Used to setup all arrays
	78
	79	if(strcmp(tracktype,"AliHLTTPCTrack")==0) fTrackType='t';
	80	else if(strcmp(tracktype,"AliHLTTPCConfMapTrack")==0) fTrackType='c';
	81	#ifdef INCLUDE_TPC_HOUGH
	82	else if(strcmp(tracktype,"AliHLTTPCHoughTrack")==0) fTrackType='h';
	83	#endif
	84	else
	85	LOG(AliHLTTPCLog::kError,"AliHLTTPCMerger::AliHLTTPCMerger","Track types")
	86	<<"Unknown tracktype"<<ENDLOG;
	87	SetArray(ntrackarrays);
	88	fCurrentTracks=0;
	89
	90	}
	91
	92	void AliHLTTPCMerger::DeleteArray()
	93	{
	94	//delete arrays
	95	for(Int_t i=0; i<fNIn;i++)
	96	{
	97	if(!fInTrack[i]) continue;
	98	delete fInTrack[i];
	99	fInTrack[i]=0;
	100	}
	101	if(fInTrack)
	102	delete[] fInTrack;
	103	if(fOutTrack)
	104	delete fOutTrack;
	105	fInTrack=0;
	106	fOutTrack=0;
	107	}
	108
	109	void AliHLTTPCMerger::SetArray(Int_t nin)
	110	{
	111	//set arrays
	112	DeleteArray();//Make sure arrays are cleaned
	113
	114	fNIn = nin;
	115	fInTrack = new AliHLTTPCTrackArray*[fNIn];
	116	for(Int_t i=0; i<fNIn;i++)
	117	{
	118	#ifdef INCLUDE_TPC_HOUGH
	119	if(fTrackType=='h')
	120	fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
	121	else
	122	#endif
	123	fInTrack[i] = new AliHLTTPCTrackArray("AliHLTTPCTrack");
	124
	125	}
	126	#ifdef INCLUDE_TPC_HOUGH
	127	if(fTrackType=='h')
	128	fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCHoughTrack");
	129	else
	130	#endif
	131	fOutTrack= new AliHLTTPCTrackArray("AliHLTTPCTrack");
	132	}
	133
	134	void AliHLTTPCMerger::Reset()
	135	{
	136	//reset
	137	for(Int_t i=0; i<fNIn;i++)
	138	{
	139	fInTrack[i]->Reset();
140	}
141	fOutTrack->Reset();
142	}
143
144	void AliHLTTPCMerger::FillTracks(Int_t ntracks, AliHLTTPCTrackSegmentData* tr)
145	{
146	//Read tracks from shared memory (or memory)
147	AliHLTTPCTrackArray *destination = GetInTracks(fCurrentTracks);
148	if(Is2Global())
febeb705	149	destination->FillTracksChecked(tr, ntracks, 0, fSlice);
a6c02c85	150	else
febeb705	151	destination->FillTracksChecked(tr, ntracks, 0);
a6c02c85	152	}
	153
	154	void AliHLTTPCMerger::AddAllTracks()
	155	{
	156	//add all tracks
	157	for(Int_t i=0; i<GetNIn();i++)
	158	{
	159	AliHLTTPCTrackArray *in = GetInTracks(i);
	160	AliHLTTPCTrackArray *out = GetOutTracks();
	161	out->AddTracks(in);
	162	}
	163	}
	164
	165	void AliHLTTPCMerger::SortGlobalTracks(AliHLTTPCTrack **tracks, Int_t ntrack)
	166	{
	167	//sort global tracks
	168	AliHLTTPCTrack *tmp = new AliHLTTPCTrack[ntrack];
	169	for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
	170	Int_t *t = new Int_t[ntrack];
	171	for(Int_t i=0;i<ntrack;i++) t[i]=-1;
	172
	173	for(Int_t j=0;j<ntrack;j++)
	174	{
	175	Double_t minr=300;
	176	Int_t mini=0;
	177	for(Int_t i=0;i<ntrack;i++)
	178	{
	179	if(!tracks[i]) continue;
	180	Double_t rr=pow(tracks[i]->GetFirstPointX(),2)+pow(tracks[i]->GetFirstPointY(),2);
	181	Double_t r=sqrt(rr);
	182	if(r<minr){
	183	minr=r;
	184	mini=i;
	185	}
	186	}
	187	t[j]=mini;
	188	tracks[mini]=0;
	189	}
	190	for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
	191	delete[] t;
	192	delete[] tmp;
	193	}
	194
	195	void AliHLTTPCMerger::SortTracks(AliHLTTPCTrack **tracks, Int_t ntrack) const
	196	{
	197	//sort tracks
	198	AliHLTTPCTrack *tmp = new AliHLTTPCTrack[ntrack];
	199	for(Int_t i=0;i<ntrack;i++) tmp[i] = tracks[i];
	200	Int_t *t = new Int_t[ntrack];
	201	for(Int_t i=0;i<ntrack;i++) t[i]=-1;
	202
	203	for(Int_t j=0;j<ntrack;j++)
	204	{
	205	Double_t minx=300;
	206	Int_t mini=0;
	207	for(Int_t i=0;i<ntrack;i++)
	208	{
	209	if(!tracks[i]) continue;
	210	if(tracks[i]->GetFirstPointX()<minx)
	211	{
	212	minx=tracks[i]->GetFirstPointX();
	213	mini=i;
	214	}
	215	}
216	t[j]=mini;
217	tracks[mini]=0;
218	}
219	for(Int_t i=0;i<ntrack;i++) tracks[i] = tmp[t[i]];
220	delete[] t;
221	delete[] tmp;
222	}
223
224	void AliHLTTPCMerger::AddTrack(AliHLTTPCTrackArray mergedtrack,AliHLTTPCTrack track)
225	{
226	// add tracks
227	AliHLTTPCTrack *t[1];
228	t[0] = track;
229	MultiMerge(mergedtrack,t,1);
230	}
231
232	AliHLTTPCTrack * AliHLTTPCMerger::MergeTracks(AliHLTTPCTrackArray mergedtrack,AliHLTTPCTrack t0,AliHLTTPCTrack *t1)
233	{
234	//merge tracks
235	AliHLTTPCTrack *t[2];
236	t[0] = t0;
237	t[1] = t1;
238	SortTracks(t,2);
239	return MultiMerge(mergedtrack,t,2);
240	}
241
242	AliHLTTPCTrack * AliHLTTPCMerger::MultiMerge(AliHLTTPCTrackArray mergedtracks,AliHLTTPCTrack *tracks, Int_t ntrack)
243	{
244	//multi merge the tracks
245	//check npoints
246	Int_t nps = 0;
247	for(Int_t i=0;i<ntrack;i++)
248	{
249	nps+=tracks[i]->GetNHits();
250	}
251	if(nps>AliHLTTPCTransform::GetNRows())
252	{
253	LOG(AliHLTTPCLog::kWarning,"AliHLTTPCMerger::MultiMerge","Adding Points")
254	<<AliHLTTPCLog::kDec<<"Too many Points: "<<nps<<ENDLOG;
255	return 0;
256	}
257
258	//create new track
259	AliHLTTPCTrack *newtrack = mergedtracks->NextTrack();
260	//copy points
261	UInt_t * nn = new UInt_t[AliHLTTPCTransform::GetNRows()];
262	nps = 0;
263
264	// for(Int_t i=0;i<ntrack;i++){
265	for(Int_t i=ntrack-1;i>=0;i--)
266	{
267	memcpy(&nn[nps],tracks[i]->GetHitNumbers(),tracks[i]->GetNHits()*sizeof(UInt_t));
268	nps+=tracks[i]->GetNHits();
269	}
270	AliHLTTPCTrack *tpf=tracks[0];
271	AliHLTTPCTrack *tpl=tracks[ntrack-1];
272	AliHLTTPCTrack *best = tpf;
273	if(tpf->GetNHits()<tpl->GetNHits() && Is2Global())
274	best = tpl;//Best means = most points and therefore best fit (in global case)
275
276	newtrack->SetNHits(nps);
277	newtrack->SetHits(nps,nn);
278	newtrack->SetFirstPoint(tpf->GetFirstPointX(),tpf->GetFirstPointY(),tpf->GetFirstPointZ());
279	newtrack->SetLastPoint(tpl->GetLastPointX(),tpl->GetLastPointY(),tpl->GetLastPointZ());
280	newtrack->SetPt(best->GetPt());
281	newtrack->SetPsi(best->GetPsi());
282	newtrack->SetTgl(best->GetTgl());
283	newtrack->SetCharge(tpf->GetCharge());
284	delete [] nn;
285	return newtrack;
286	}
287
288	void* AliHLTTPCMerger::GetNtuple(char *varlist) const
289	{
290	//get ntuple
291	#ifdef use_root
292	TNtuple* nt = new TNtuple("ntuple","ntuple",varlist);
293	return (void*) nt;
294	#else
295	return 0;
296	#endif
297	}
298
299	void* AliHLTTPCMerger::GetNtuple() const
300	{
301	//get ntuple
302	#ifdef use_root
303	TNtuple* nt = new TNtuple("ntuple","ntuple",
304	"dx:dy:dz:dk:dpsi:dtgl:dq:disx:disy:disz:dis:n0:n1:diff:drx:dry:drz");
305	return (void*) nt;
306	#else
307	return 0;
308	#endif
309	}
310
311	Bool_t AliHLTTPCMerger::WriteNtuple(char filename, void nt) const
312	{
313	//write ntuple
314	#ifdef use_root
315	TNtuple ntuple=(TNtuple ) nt;
316	TFile *f = new TFile(filename,"RECREATE");
317	ntuple->Write();
318	f->Close();
319	delete ntuple;
320	return kTRUE;
321	#else
322	return kFALSE;
323	#endif
324	}
325
326	void AliHLTTPCMerger::FillNtuple(void nt,AliHLTTPCTrack innertrack,AliHLTTPCTrack *outertrack)
327	{
328	//fill ntuple
329	Float_t data[17];
330	if(outertrack->IsPoint()&&innertrack->IsPoint())
331	{
332	data[0] =Float_t(innertrack->GetPointX()-outertrack->GetPointX());
333	data[1] =Float_t(innertrack->GetPointY()-outertrack->GetPointY());
334	data[2] =Float_t(innertrack->GetPointZ()-outertrack->GetPointZ());
335	data[3] =Float_t(innertrack->GetKappa()-outertrack->GetKappa());
336	Double_t psi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
337	if(psi>AliHLTTPCTransform::Pi()) psi-=AliHLTTPCTransform::TwoPi();
338	else if(psi<-AliHLTTPCTransform::Pi()) psi+=AliHLTTPCTransform::TwoPi();
339	data[4] =Float_t(psi);
340	data[5] =Float_t(innertrack->GetTgl()-outertrack->GetTgl());
341	data[6] =Float_t(innertrack->GetCharge()-outertrack->GetCharge());
342	data[7] =Float_t(innertrack->GetLastPointX()-outertrack->GetFirstPointX());
343	data[8] =Float_t(innertrack->GetLastPointY()-outertrack->GetFirstPointY());
344	data[9] =Float_t(innertrack->GetLastPointZ()-outertrack->GetFirstPointZ());
345	data[10] =sqrt(pow(data[7],2)+pow(data[8],2)+pow(data[9],2));
346	data[11]= outertrack->GetNHits();
347	data[12]= innertrack->GetNHits();
348	data[13] = Float_t(TrackDiff(innertrack,outertrack));
349	data[14]=0;
350	data[15]=0;
351	data[16]=0;
352	#ifdef use_root
353	TNtuple ntuple = (TNtuple ) nt;
354	ntuple->Fill(data);
355	#endif
356	}
357	}
358
359	void AliHLTTPCMerger::FillNtuple(void nt,Float_t data) const
360	{
361	//fill ntuple
362	#ifdef use_root
363	TNtuple ntuple = (TNtuple ) nt;
364	ntuple->Fill(data);
365	#endif
366	}
367
368	Double_t AliHLTTPCMerger::GetAngle(Double_t a1,Double_t a2)
369	{
370	//get angle
371	Double_t da = a1 - a2 + 4*AliHLTTPCTransform::Pi();
372	da = fmod(da,AliHLTTPCTransform::TwoPi());
373	if(da>AliHLTTPCTransform::Pi()) da = AliHLTTPCTransform::TwoPi()-da;
374	return da;
375	}
376
377	void AliHLTTPCMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxkappa, Double_t maxpsi, Double_t maxtgl)
378	{
379	//set parameters for merger
380	fMaxY = maxy;
381	fMaxZ = maxz;
382	fMaxKappa = maxkappa;
383	fMaxPsi = maxpsi;
384	fMaxTgl = maxtgl;
385	}
386
387	Bool_t AliHLTTPCMerger::IsTrack(AliHLTTPCTrack innertrack,AliHLTTPCTrack outertrack)
388	{
389	//is track to be merged
390	if(innertrack->GetCharge()!=outertrack->GetCharge()) return kFALSE;
391	if( (!innertrack->IsPoint()) \|\| (!outertrack->IsPoint()) ) return kFALSE;
392	if(innertrack->GetNHits()+outertrack->GetNHits()>AliHLTTPCTransform::GetNRows()) return kFALSE;
393
394	if(fabs(innertrack->GetPointY()-outertrack->GetPointY()) >fMaxY) return kFALSE;
395	if(fabs(innertrack->GetPointZ()-outertrack->GetPointZ()) >fMaxZ) return kFALSE;
396	if(fabs(innertrack->GetKappa()-outertrack->GetKappa()) >fMaxKappa) return kFALSE;
397	if(GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi()) >fMaxPsi) return kFALSE;
398	if(fabs(innertrack->GetTgl()-outertrack->GetTgl()) >fMaxTgl) return kFALSE;
399	//if no rejection up to this point: merge!!
400	return kTRUE;
401	}
402
403	Bool_t AliHLTTPCMerger::IsRTrack(AliHLTTPCTrack innertrack,AliHLTTPCTrack outertrack)
404	{
405	//same as IsTrack
406	return IsTrack(innertrack,outertrack);
407	}
408
409	Double_t AliHLTTPCMerger::TrackDiff(AliHLTTPCTrack innertrack,AliHLTTPCTrack outertrack)
410	{
411	//return track difference
412	Double_t diff =-1;
413	Double_t x[4],y[4],z[4],dy[4],dz[4];
414	AliHLTTPCTrack *tracks[2];
415
416	tracks[0] = innertrack;
417	tracks[1] = outertrack;
418	SortGlobalTracks(tracks,2);
419	innertrack = tracks[0];
420	outertrack = tracks[1];
421
422	x[0] = innertrack->GetFirstPointX();
423	x[1] = innertrack->GetLastPointX();
424	x[2] = outertrack->GetFirstPointX();
425	x[3] = outertrack->GetLastPointX();
426
427	y[0] = innertrack->GetFirstPointY();
428	y[1] = innertrack->GetLastPointY();
429	y[2] = outertrack->GetFirstPointY();
430	y[3] = outertrack->GetLastPointY();
431
432	z[0] = innertrack->GetFirstPointZ();
433	z[1] = innertrack->GetLastPointZ();
434	z[2] = outertrack->GetFirstPointZ();
435	z[3] = outertrack->GetLastPointZ();
436
437
438	outertrack->CalculatePoint(x[0]);
439	if(!outertrack->IsPoint()) return diff;
440	dy[0] = fabs(y[0] - outertrack->GetPointY());
441	dz[0] = fabs(z[0] - outertrack->GetPointZ());
442
443	outertrack->CalculatePoint(x[1]);
444	if(!outertrack->IsPoint()) return diff;
445	dy[1] = fabs(y[1] - outertrack->GetPointY());
446	dz[1] = fabs(z[1] - outertrack->GetPointZ());
447
448	innertrack->CalculatePoint(x[2]);
449	if(!innertrack->IsPoint()) return diff;
450	dy[2] = fabs(y[2] - innertrack->GetPointY());
451	dz[2] = fabs(z[2] - innertrack->GetPointZ());
452
453	innertrack->CalculatePoint(x[3]);
454	if(!innertrack->IsPoint()) return diff;
455	dy[3] = fabs(y[3] - innertrack->GetPointY());
456	dz[3] = fabs(z[3] - innertrack->GetPointZ());
457
458	diff=0;
459	for(Int_t i=0;i<4;i++)
460	diff+=sqrt(dy[i]dy[i]+dz[i]dz[i]);
461	return diff;
462	}
463
464	void AliHLTTPCMerger::PrintDiff(AliHLTTPCTrack innertrack,AliHLTTPCTrack outertrack)
465	{
466	// print difference
467	if(!innertrack->IsPoint()\|\|!outertrack->IsPoint())
468	{
469	LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","No Points")<<ENDLOG;
470	//cerr<<"AliHLTTPCMerger::PrintDiff: No Points"<<endl;
471	//cerr<<"---------------------------"<<endl;
472	return;
473	}
474
475	Double_t dx = innertrack->GetPointX()-outertrack->GetPointX();
476	Double_t dy = innertrack->GetPointY()-outertrack->GetPointY();
477	Double_t dz = innertrack->GetPointZ()-outertrack->GetPointZ();
478	Double_t dk = innertrack->GetKappa()-outertrack->GetKappa();
479	Double_t dpsi= innertrack->GetPointPsi() - outertrack->GetPointPsi();
480	if(dpsi>AliHLTTPCTransform::Pi()) dpsi-=AliHLTTPCTransform::TwoPi();
481	else if(dpsi<-AliHLTTPCTransform::Pi())dpsi+=AliHLTTPCTransform::TwoPi();
482	//Double_t dpsi = GetAngle(innertrack->GetPointPsi(),outertrack->GetPointPsi());
483	Double_t dtgl= innertrack->GetTgl()-outertrack->GetTgl();
484	Double_t dq =innertrack->GetCharge()-outertrack->GetCharge();
485
486	LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::PrintDiff","Points") <<"dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<" dk: "<<dk<<" dpsi: "<<dpsi<<" dtgl: "<<dtgl<<" dq: "<<dq<<ENDLOG;
487	//fprintf(stderr,"dx: %4f dy: %4f dz: %4f dk: %4f dpsi: %4f dtgl: %4f dq: %4f\n",dx,dy,dz,dk,dpsi,dtgl,dq);
488	//cerr<<"---------------------------"<<endl;
489
490	}
491
492	void AliHLTTPCMerger::Print()
493	{
494	// print some infos
495	for(Int_t i=0; i<fNIn; i++)
496	{
497	AliHLTTPCTrackArray *ttt= GetInTracks(i);
498	for(Int_t j =0;j<ttt->GetNTracks();j++)
499	{
500	AliHLTTPCTrack *track=ttt->GetCheckedTrack(j);
501	if(!track) continue;
502	track->CalculateHelix();
503	// Double_t angle = atan2(track->GetLastPointY(),track->GetLastPointX());
504	// if(angle<0) angle+=AliHLTTPCTransform::Pi();
505	if(track->CalculatePoint(135))
506	// if(!track->CalculateEdgePoint(angle)) cerr<<"**************"<<endl;
507	// if(track->CalculatePoint(track->GetLastPointX()))
508	// if(track->CalculatePoint(0))
509	{
510	// PrintTrack(track);
511	// track->CalculateReferencePoint(AliHLTTPCTransform::Pi()/180.);
512	track->CalculateReferencePoint(0.001);
513	Float_t dx=(float)track->GetPointX()-track->GetPointX();
514	Float_t dy=(float)track->GetPointY()-track->GetPointY();
515	Float_t dz=(float)track->GetPointZ()-track->GetPointZ();
516	LOG(AliHLTTPCLog::kInformational,"AliHLTTPCMerger::Print","RefPoint") <<"npt: "<<track->GetNHits()<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<ENDLOG;
517
518	//fprintf(stderr,"npt: %3d dx: %8.5f dy: %8.5f dz: %8.5f\n",track->GetNHits(),dx,dy,dz);
519	//cerr<<"---------------------------"<<endl;
520	}
521	}
522	}
523	}
524
525	void AliHLTTPCMerger::PrintTrack(AliHLTTPCTrack *track)
526	{
527	//print track info
528	fprintf(stderr,"npt: %3d pt: %.2f psi: %.2f tgl: %5.2f q: %2d\n",
529	track->GetNHits(),track->GetPt(),track->GetPsi(),
530	track->GetTgl(),track->GetCharge());
531	fprintf(stderr,
532	"x1: %6.2f y1: %6.2f z1: %6.2f xl: %6.2f yl: %6.2f zl: %6.2f\n",
533	track->GetFirstPointX(),track->GetFirstPointY(),track->GetFirstPointZ(),
534	track->GetLastPointX(),track->GetLastPointY(),track->GetLastPointZ());
535	if(track->IsPoint())
536	{
537	fprintf(stderr,
538	"R: %.2f Xc: %.2f Yc: %.2f Xp: %.2f Yp: %.2f Zp: %.2f Psip: %.2f\n",
539	track->GetRadius(),track->GetCenterX(),track->GetCenterY(),
540	track->GetPointX(),track->GetPointY(),track->GetPointZ(),
541	track->GetPointPsi());
542	}
543	}