New beam pipe geometry: compensator in ZDC code
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderHLT.cxx
CommitLineData
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
26ClassImp(AliJetParticlesReaderHLT)
27
28AliJetParticlesReaderHLT::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
39AliJetParticlesReaderHLT::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
54AliJetParticlesReaderHLT::~AliJetParticlesReaderHLT()
55{
56 //desctructor
57}
58
59Int_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;
117 v.SetX(vertexpos[0]);
c5368c62 118 v.SetY(vertexpos[1]);
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 }
167
46a140bc 168 const Float_t kpt=TMath::Sqrt(px*px+py*py);
78a77a8f 169 //cout << px << " " << py << " " << " " << kpt << endl;
170
46a140bc 171 const Float_t kp=TMath::Sqrt(pz*pz+kpt*kpt);
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;
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)*(2*TMath::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)*(2*TMath::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