Returned to old version where parameters are given at first point.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderHLT.cxx
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>
20 #include <AliL3Track.h>
21 #include <AliL3Vertex.h>
22 #include <AliKalmanTrack.h>
23 #include <AliJetEventParticles.h>
24 #include "AliJetParticlesReaderHLT.h"
25
26 ClassImp(AliJetParticlesReaderHLT)
27
28 AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(Bool_t bMapper, const Char_t* esdfilename) :
29   AliJetParticlesReaderESD(0,esdfilename),
30   fTrackerType(bMapper),
31   fMinHits(0),
32   fMinWeight(0)
33 {
34   //constructor
35 }
36
37 /********************************************************************/
38   
39 AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(
40                                       Bool_t bMapper,
41                                       TObjArray* dirs,
42                                       const Char_t* esdfilename) :
43   AliJetParticlesReaderESD(0,dirs,esdfilename),
44   fTrackerType(bMapper),
45   fMinHits(0),
46   fMinWeight(0)
47
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
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
77   Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
78   if((!fOwner) || (fEventParticles==0)) 
79     fEventParticles = new AliJetEventParticles();
80
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
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]);
113   //cout << vertexpos[0] << " " << vertexpos[1] << " " << vertexpos[2] << endl;
114
115   AliL3Track l3;
116   AliL3Vertex v;
117   v.SetX(vertexpos[0]);
118   v.SetY(vertexpos[1]);
119   v.SetZ(vertexpos[2]);
120   
121   Double_t xc=0.,yc=0.,zc=0.;
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     }
129
130     if (kesdtrack == 0)
131       {
132         Error("ReadESD","Can not get track %d", i);
133         continue;
134       }
135
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
148     if(0&&fTrackerType){
149       //if(!kesdtrack->ComesFromMainVertex()) continue;
150       cout << kesdtrack->GetPx() << " " << kesdtrack->GetPy() << " " << kesdtrack->GetPt() << endl;
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();
158       l3.GetClosestPoint(&v,xc,yc,zc);
159       if(TMath::Abs(zc)>10.) continue;
160       l3.SetFirstPoint(vertexpos[0],vertexpos[1],vertexpos[2]);
161       //l3.CalculateHelix();
162       l3.UpdateToFirstPoint();
163       px=l3.GetPx();
164       py=l3.GetPy();
165       pz=l3.GetPz();
166     }
167
168     const Float_t kpt=TMath::Sqrt(px*px+py*py);
169     //cout << px << " " << py << " " << " " << kpt << endl;
170
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);
174
175     if(IsAcceptedParticle(kpt,kphi,keta))
176       fEventParticles->AddParticle(px,py,pz,kp,i,kesdtrack->GetMCid(),knhits,kpt,kphi,keta);
177   } //loop over tracks
178
179   return kTRUE;
180 }
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