[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderHLT.cxx

// $Id$

//___________________________________________________________________________
/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// File reader for HLT tracks ESD                                          //
//                                                                         //
// loizides@ikf.uni-frankfurt.de                                           //
/////////////////////////////////////////////////////////////////////////////

#include <Riostream.h>
#include <TMath.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
#include <TFile.h>
#include <AliESD.h>
#include <AliESDtrack.h>
#include <AliESDHLTtrack.h>
#include <AliL3Track.h>
#include <AliL3Vertex.h>
#include <AliKalmanTrack.h>
#include <AliJetEventParticles.h>
#include "AliJetParticlesReaderHLT.h"

ClassImp(AliJetParticlesReaderHLT)

AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(Bool_t bMapper, const Char_t* esdfilename) :
  AliJetParticlesReaderESD(0,esdfilename),
  fTrackerType(bMapper),
  fMinHits(0),
  fMinWeight(0)
{
  //constructor
}

/********************************************************************/
  
AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(
				      Bool_t bMapper,
                                      TObjArray* dirs,
                                      const Char_t* esdfilename) :
  AliJetParticlesReaderESD(0,dirs,esdfilename),
  fTrackerType(bMapper),
  fMinHits(0),
  fMinWeight(0)

{
  //constructor
}

/********************************************************************/

AliJetParticlesReaderHLT::~AliJetParticlesReaderHLT()
{
  //desctructor
}

Int_t AliJetParticlesReaderHLT::ReadESD(AliESD* esd)
{
  //Reads one ESD

  if (esd == 0)
   {
     Error("ReadESD","ESD is NULL");
     return kFALSE;
   }

  Float_t mf = esd->GetMagneticField(); 
  if (mf <= 0.0)  
  {
     Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
     return kFALSE;
  }
  AliKalmanTrack::SetMagneticField(mf/10.);

  Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
  if((!fOwner) || (fEventParticles==0)) 
    fEventParticles = new AliJetEventParticles();

  Int_t ntr=0;
  if(fTrackerType){
    ntr =esd->GetNumberOfHLTHoughTracks();
    Info("ReadESD","Found %d conformal tracks.",ntr);
  } else {
    ntr=esd->GetNumberOfHLTHoughTracks();
    Info("ReadESD","Found %d hough tracks.",ntr);
  }
  fEventParticles->Reset(ntr);

  TString headdesc="";
  headdesc+="Run ";
  headdesc+=esd->GetRunNumber();
  headdesc+=": Ev ";
  headdesc+=esd->GetEventNumber();
  fEventParticles->SetHeader(headdesc);

  Double_t vertexpos[3];//vertex position
  const AliESDVertex* kvertex = esd->GetVertex();
  if (kvertex == 0)
   {
     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
     vertexpos[0] = 0.0;
     vertexpos[1] = 0.0;
     vertexpos[2] = 0.0;
   }
  else
   {
     kvertex->GetXYZ(vertexpos);
   }

  fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
  //cout << vertexpos[0] << " " << vertexpos[1] << " " << vertexpos[2] << endl;

  AliL3Track l3;
  AliL3Vertex v;
  v.SetX(vertexpos[0]);
  v.SetY(vertexpos[1]);
  v.SetZ(vertexpos[2]);
  
  Double_t xc=0.,yc=0.,zc=0.;
  for (Int_t i = 0;i<ntr; i++) {
    AliESDHLTtrack *kesdtrack;
    if(fTrackerType){
      kesdtrack=esd->GetHLTConfMapTrack(i);
    } else {
      kesdtrack=esd->GetHLTHoughTrack(i);
    }

    if (kesdtrack == 0)
      {
        Error("ReadESD","Can not get track %d", i);
        continue;
      }

    //const Float_t kpid=kesdtrack->GetPID();
    const Int_t knhits=kesdtrack->GetNHits();
    const Int_t kweight=kesdtrack->GetWeight();
    //cout << i << " " << kweight << " " << knhits << endl;
    if((fMinHits>0) && (knhits<fMinHits)) continue;    
    if(kweight>1000) continue; //avoid ghosts 
    if((fMinWeight>0) && (kweight<fMinWeight)) continue;

    Float_t px=kesdtrack->GetPx();
    Float_t py=kesdtrack->GetPy();
    Float_t pz=kesdtrack->GetPz();

    if(0&&fTrackerType){
      //if(!kesdtrack->ComesFromMainVertex()) continue;
      cout << kesdtrack->GetPx() << " " << kesdtrack->GetPy() << " " << kesdtrack->GetPt() << endl;
      l3.SetFirstPoint(kesdtrack->GetFirstPointX(),kesdtrack->GetFirstPointY(),kesdtrack->GetFirstPointZ());
      l3.SetLastPoint(kesdtrack->GetLastPointX(),kesdtrack->GetLastPointY(),kesdtrack->GetLastPointZ());
      l3.SetCharge(kesdtrack->GetCharge());
      l3.SetPt(kesdtrack->GetPt());
      l3.SetTgl(kesdtrack->GetTgl());
      l3.SetPsi(kesdtrack->GetPsi());
      l3.CalculateHelix();
      l3.GetClosestPoint(&v,xc,yc,zc);
      if(TMath::Abs(zc)>10.) continue;
      l3.SetFirstPoint(vertexpos[0],vertexpos[1],vertexpos[2]);
      //l3.CalculateHelix();
      l3.UpdateToFirstPoint();
      px=l3.GetPx();
      py=l3.GetPy();
      pz=l3.GetPz();
    }

    const Float_t kpt=TMath::Sqrt(px*px+py*py);
    //cout << px << " " << py << " " << " " << kpt << endl;

    const Float_t kp=TMath::Sqrt(pz*pz+kpt*kpt);
    const Float_t keta=0.5*TMath::Log((kp+pz+1e-30)/(kp-pz+1e-30)); 
    const Float_t kphi=TMath::Pi()+TMath::ATan2(-py,-px);

    if(IsAcceptedParticle(kpt,kphi,keta))
      fEventParticles->AddParticle(px,py,pz,kp,i,kesdtrack->GetMCid(),knhits,kpt,kphi,keta);
  } //loop over tracks

  return kTRUE;
}

#if 0
  SetChi2(0.);
  if(t.GetNHits()==1)
    SetNumberOfClusters(0);
  else
    SetNumberOfClusters(t.GetNHits());
  SetLabel(t.GetMCid());
  SetMass(0.13957);

  fdEdx=0;
  fAlpha = fmod((t.GetSector()+0.5)*(2*TMath::Pi()/18),2*TMath::Pi());
  if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
  else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();

  //First the emiision angle
  Double_t psi = t.GetPsi()-(t.GetSector()+0.5)*(2*TMath::Pi()/18);

  //Then local x,y coordinates
  Double_t radius = t.GetPt()*GetConvConst();
  Double_t xhit = 82.97; //Position at first TPC padrow
  Double_t trackphi0 = psi + (-t.GetCharge())*TMath::Pi()/2;
  Double_t x0 = t.GetFirstPointX()*TMath::Cos(fAlpha) + t.GetFirstPointY()*TMath::Sin(fAlpha);
  Double_t y0 = t.GetFirstPointY()*TMath::Cos(fAlpha) - t.GetFirstPointX()*TMath::Sin(fAlpha);
  Double_t centerx = radius *  cos(trackphi0) + x0;
  Double_t centery = radius *  sin(trackphi0) + y0;
  Double_t aa = (xhit - centerx)*(xhit - centerx);
  Double_t r2 = radius*radius;
  if(aa > r2) throw "AliITStrackV2: conversion failed !\n";
  Double_t aa2 = sqrt(r2 - aa);
  Double_t y1 = centery + aa2;
  Double_t y2 = centery - aa2;
  Double_t yhit = y1;
  if(fabs(y2) < fabs(y1)) yhit = y2;

  //Local z coordinate
  Double_t angle1 = atan2((yhit - centery),(xhit - centerx));
  if(angle1 < 0) angle1 += 2.*TMath::Pi();
  Double_t angle2 = atan2((x0-centery),(y0-centerx));
  if(angle2 < 0) angle2 += 2.*TMath::Pi();
  Double_t diffangle = angle1 - angle2;
  diffangle = fmod(diffangle,2.*TMath::Pi());
  if(((-t.GetCharge())*diffangle) < 0) diffangle = diffangle - (-t.GetCharge())*2.*TMath::Pi();
  Double_t stot = fabs(diffangle)*radius;
  Double_t zhit;
  if(t.GetNHits()==1)
    zhit = zvertex + stot*t.GetTgl();
  else
    zhit = t.GetFirstPointZ() + stot*t.GetTgl();

  //Local sine of track azimuthal angle
  if((-t.GetCharge())<0)
    radius = -radius;
  Double_t sinbeta = -1.*(centerx - xhit)/radius;

  //Filling of the track paramaters
  fX=xhit;
  fP0=yhit;
  fP1=zhit;
  fP2=sinbeta;
  fP3=t.GetTgl();
  fP4=1./radius;
#endif
Commit	Line	Data
5d60c8f9	1	// $Id$
	2
	3	//___________________________________________________________________________
	4	/////////////////////////////////////////////////////////////////////////////
	5	// //
	6	// File reader for HLT tracks ESD //
	7	// //
	8	// loizides@ikf.uni-frankfurt.de //
	9	/////////////////////////////////////////////////////////////////////////////
	10
	11	#include <Riostream.h>
	12	#include <TMath.h>
	13	#include <TString.h>
	14	#include <TObjString.h>
	15	#include <TTree.h>
	16	#include <TFile.h>
	17	#include <AliESD.h>
	18	#include <AliESDtrack.h>
	19	#include <AliESDHLTtrack.h>
75a0c43e	20	#include <AliL3Track.h>
46a140bc	21	#include <AliL3Vertex.h>
04a02430	22	#include <AliKalmanTrack.h>
5d60c8f9	23	#include <AliJetEventParticles.h>
	24	#include "AliJetParticlesReaderHLT.h"
	25
	26	ClassImp(AliJetParticlesReaderHLT)
	27
	28	AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(Bool_t bMapper, const Char_t* esdfilename) :
d9a296b1	29	AliJetParticlesReaderESD(0,esdfilename),
04a02430	30	fTrackerType(bMapper),
	31	fMinHits(0),
	32	fMinWeight(0)
5d60c8f9	33	{
	34	//constructor
	35	}
	36
	37	/********************************************************************/
	38
	39	AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(
	40	Bool_t bMapper,
	41	TObjArray* dirs,
	42	const Char_t* esdfilename) :
d9a296b1	43	AliJetParticlesReaderESD(0,dirs,esdfilename),
04a02430	44	fTrackerType(bMapper),
	45	fMinHits(0),
	46	fMinWeight(0)
	47
5d60c8f9	48	{
	49	//constructor
	50	}
	51
	52	/********************************************************************/
	53
	54	AliJetParticlesReaderHLT::~AliJetParticlesReaderHLT()
	55	{
	56	//desctructor
	57	}
	58
	59	Int_t AliJetParticlesReaderHLT::ReadESD(AliESD* esd)
	60	{
	61	//Reads one ESD
	62
	63	if (esd == 0)
	64	{
	65	Error("ReadESD","ESD is NULL");
	66	return kFALSE;
	67	}
	68
04a02430	69	Float_t mf = esd->GetMagneticField();
	70	if (mf <= 0.0)
	71	{
	72	Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
	73	return kFALSE;
	74	}
	75	AliKalmanTrack::SetMagneticField(mf/10.);
	76
a86edc4e	77	Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
5d60c8f9	78	if((!fOwner) \|\| (fEventParticles==0))
	79	fEventParticles = new AliJetEventParticles();
	80
04a02430	81	Int_t ntr=0;
	82	if(fTrackerType){
	83	ntr =esd->GetNumberOfHLTHoughTracks();
	84	Info("ReadESD","Found %d conformal tracks.",ntr);
	85	} else {
	86	ntr=esd->GetNumberOfHLTHoughTracks();
	87	Info("ReadESD","Found %d hough tracks.",ntr);
	88	}
	89	fEventParticles->Reset(ntr);
	90
	91	TString headdesc="";
	92	headdesc+="Run ";
	93	headdesc+=esd->GetRunNumber();
	94	headdesc+=": Ev ";
	95	headdesc+=esd->GetEventNumber();
	96	fEventParticles->SetHeader(headdesc);
	97
5d60c8f9	98	Double_t vertexpos[3];//vertex position
	99	const AliESDVertex* kvertex = esd->GetVertex();
	100	if (kvertex == 0)
	101	{
	102	Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
	103	vertexpos[0] = 0.0;
	104	vertexpos[1] = 0.0;
	105	vertexpos[2] = 0.0;
	106	}
	107	else
	108	{
	109	kvertex->GetXYZ(vertexpos);
	110	}
	111
	112	fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
74e5a625	113	//cout << vertexpos[0] << " " << vertexpos[1] << " " << vertexpos[2] << endl;
5d60c8f9	114
75a0c43e	115	AliL3Track l3;
46a140bc	116	AliL3Vertex v;
46a140bc	117	v.SetX(vertexpos[0]);
c5368c62	118	v.SetY(vertexpos[1]);
c5368c62	119	v.SetZ(vertexpos[2]);
78a77a8f	120
46a140bc	121	Double_t xc=0.,yc=0.,zc=0.;
04a02430	122	for (Int_t i = 0;i<ntr; i++) {
	123	AliESDHLTtrack *kesdtrack;
	124	if(fTrackerType){
	125	kesdtrack=esd->GetHLTConfMapTrack(i);
	126	} else {
	127	kesdtrack=esd->GetHLTHoughTrack(i);
	128	}
5d60c8f9	129
04a02430	130	if (kesdtrack == 0)
5d60c8f9	131	{
	132	Error("ReadESD","Can not get track %d", i);
	133	continue;
	134	}
	135
46a140bc	136	//const Float_t kpid=kesdtrack->GetPID();
	137	const Int_t knhits=kesdtrack->GetNHits();
	138	const Int_t kweight=kesdtrack->GetWeight();
	139	//cout << i << " " << kweight << " " << knhits << endl;
	140	if((fMinHits>0) && (knhits<fMinHits)) continue;
	141	if(kweight>1000) continue; //avoid ghosts
	142	if((fMinWeight>0) && (kweight<fMinWeight)) continue;
	143
	144	Float_t px=kesdtrack->GetPx();
	145	Float_t py=kesdtrack->GetPy();
	146	Float_t pz=kesdtrack->GetPz();
	147
78a77a8f	148	if(0&&fTrackerType){
75a0c43e	149	//if(!kesdtrack->ComesFromMainVertex()) continue;
78a77a8f	150	cout << kesdtrack->GetPx() << " " << kesdtrack->GetPy() << " " << kesdtrack->GetPt() << endl;
75a0c43e	151	l3.SetFirstPoint(kesdtrack->GetFirstPointX(),kesdtrack->GetFirstPointY(),kesdtrack->GetFirstPointZ());
	152	l3.SetLastPoint(kesdtrack->GetLastPointX(),kesdtrack->GetLastPointY(),kesdtrack->GetLastPointZ());
	153	l3.SetCharge(kesdtrack->GetCharge());
	154	l3.SetPt(kesdtrack->GetPt());
	155	l3.SetTgl(kesdtrack->GetTgl());
	156	l3.SetPsi(kesdtrack->GetPsi());
	157	l3.CalculateHelix();
46a140bc	158	l3.GetClosestPoint(&v,xc,yc,zc);
	159	if(TMath::Abs(zc)>10.) continue;
	160	l3.SetFirstPoint(vertexpos[0],vertexpos[1],vertexpos[2]);
74e5a625	161	//l3.CalculateHelix();
75a0c43e	162	l3.UpdateToFirstPoint();
46a140bc	163	px=l3.GetPx();
	164	py=l3.GetPy();
	165	pz=l3.GetPz();
75a0c43e	166	}
75a0c43e	167
46a140bc	168	const Float_t kpt=TMath::Sqrt(pxpx+pypy);
78a77a8f	169	//cout << px << " " << py << " " << " " << kpt << endl;
78a77a8f	170
46a140bc	171	const Float_t kp=TMath::Sqrt(pzpz+kptkpt);
	172	const Float_t keta=0.5*TMath::Log((kp+pz+1e-30)/(kp-pz+1e-30));
	173	const Float_t kphi=TMath::Pi()+TMath::ATan2(-py,-px);
5d60c8f9	174
04a02430	175	if(IsAcceptedParticle(kpt,kphi,keta))
46a140bc	176	fEventParticles->AddParticle(px,py,pz,kp,i,kesdtrack->GetMCid(),knhits,kpt,kphi,keta);
04a02430	177	} //loop over tracks
5d60c8f9	178
5d60c8f9	179	return kTRUE;
5d60c8f9	180	}
78a77a8f	181
	182	#if 0
	183	SetChi2(0.);
	184	if(t.GetNHits()==1)
	185	SetNumberOfClusters(0);
	186	else
	187	SetNumberOfClusters(t.GetNHits());
	188	SetLabel(t.GetMCid());
	189	SetMass(0.13957);
	190
	191	fdEdx=0;
	192	fAlpha = fmod((t.GetSector()+0.5)(2TMath::Pi()/18),2*TMath::Pi());
	193	if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
	194	else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
	195
	196	//First the emiision angle
	197	Double_t psi = t.GetPsi()-(t.GetSector()+0.5)(2TMath::Pi()/18);
	198
	199	//Then local x,y coordinates
	200	Double_t radius = t.GetPt()*GetConvConst();
	201	Double_t xhit = 82.97; //Position at first TPC padrow
	202	Double_t trackphi0 = psi + (-t.GetCharge())*TMath::Pi()/2;
	203	Double_t x0 = t.GetFirstPointX()TMath::Cos(fAlpha) + t.GetFirstPointY()TMath::Sin(fAlpha);
	204	Double_t y0 = t.GetFirstPointY()TMath::Cos(fAlpha) - t.GetFirstPointX()TMath::Sin(fAlpha);
	205	Double_t centerx = radius * cos(trackphi0) + x0;
	206	Double_t centery = radius * sin(trackphi0) + y0;
	207	Double_t aa = (xhit - centerx)*(xhit - centerx);
	208	Double_t r2 = radius*radius;
	209	if(aa > r2) throw "AliITStrackV2: conversion failed !\n";
	210	Double_t aa2 = sqrt(r2 - aa);
	211	Double_t y1 = centery + aa2;
	212	Double_t y2 = centery - aa2;
	213	Double_t yhit = y1;
	214	if(fabs(y2) < fabs(y1)) yhit = y2;
	215
	216	//Local z coordinate
	217	Double_t angle1 = atan2((yhit - centery),(xhit - centerx));
	218	if(angle1 < 0) angle1 += 2.*TMath::Pi();
	219	Double_t angle2 = atan2((x0-centery),(y0-centerx));
	220	if(angle2 < 0) angle2 += 2.*TMath::Pi();
	221	Double_t diffangle = angle1 - angle2;
	222	diffangle = fmod(diffangle,2.*TMath::Pi());
	223	if(((-t.GetCharge())diffangle) < 0) diffangle = diffangle - (-t.GetCharge())2.*TMath::Pi();
	224	Double_t stot = fabs(diffangle)*radius;
	225	Double_t zhit;
	226	if(t.GetNHits()==1)
	227	zhit = zvertex + stot*t.GetTgl();
	228	else
	229	zhit = t.GetFirstPointZ() + stot*t.GetTgl();
	230
	231	//Local sine of track azimuthal angle
	232	if((-t.GetCharge())<0)
	233	radius = -radius;
	234	Double_t sinbeta = -1.*(centerx - xhit)/radius;
	235
	236	//Filling of the track paramaters
	237	fX=xhit;
	238	fP0=yhit;
	239	fP1=zhit;
	240	fP2=sinbeta;
	241	fP3=t.GetTgl();
	242	fP4=1./radius;
	243	#endif
	244