Bugfixes again for HLT.
[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 <AliL3ITStrack.h>
22 #include <AliL3Vertex.h>
23 #include <AliKalmanTrack.h>
24 #include <AliJetEventParticles.h>
25 #include "AliJetParticlesReaderHLT.h"
26
27 ClassImp(AliJetParticlesReaderHLT)
28
29 AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(Bool_t bMapper, const Char_t* esdfilename) :
30   AliJetParticlesReaderESD(0,esdfilename),
31   fTrackerType(bMapper),
32   fMinHits(0),
33   fMinWeight(0)
34 {
35   //constructor
36 }
37
38 /********************************************************************/
39   
40 AliJetParticlesReaderHLT::AliJetParticlesReaderHLT(
41                                       Bool_t bMapper,
42                                       TObjArray* dirs,
43                                       const Char_t* esdfilename) :
44   AliJetParticlesReaderESD(0,dirs,esdfilename),
45   fTrackerType(bMapper),
46   fMinHits(0),
47   fMinWeight(0)
48
49 {
50   //constructor
51 }
52
53 /********************************************************************/
54
55 AliJetParticlesReaderHLT::~AliJetParticlesReaderHLT()
56 {
57   //desctructor
58 }
59
60 Int_t AliJetParticlesReaderHLT::ReadESD(AliESD* esd)
61 {
62   //Reads one ESD
63
64   if (esd == 0)
65    {
66      Error("ReadESD","ESD is NULL");
67      return kFALSE;
68    }
69
70   Float_t mf = esd->GetMagneticField(); 
71   if (mf <= 0.0)  
72   {
73      Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
74      return kFALSE;
75   }
76   AliKalmanTrack::SetMagneticField(mf/10.);
77
78   Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
79   if((!fOwner) || (fEventParticles==0)) 
80     fEventParticles = new AliJetEventParticles();
81
82   Int_t ntr=0;
83   if(fTrackerType){
84     ntr =esd->GetNumberOfHLTHoughTracks();
85     Info("ReadESD","Found %d conformal tracks.",ntr);
86   } else {
87     ntr=esd->GetNumberOfHLTHoughTracks();
88     Info("ReadESD","Found %d hough tracks.",ntr);
89   }
90   fEventParticles->Reset(ntr);
91
92   TString headdesc="";
93   headdesc+="Run ";
94   headdesc+=esd->GetRunNumber();
95   headdesc+=": Ev ";
96   headdesc+=esd->GetEventNumber();
97   fEventParticles->SetHeader(headdesc);
98
99   Double_t vertexpos[3];//vertex position
100   const AliESDVertex* kvertex = esd->GetVertex();
101   if (kvertex == 0)
102    {
103      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
104      vertexpos[0] = 0.0;
105      vertexpos[1] = 0.0;
106      vertexpos[2] = 0.0;
107    }
108   else
109    {
110      kvertex->GetXYZ(vertexpos);
111    }
112
113   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
114   //cout << vertexpos[0] << " " << vertexpos[1] << " " << vertexpos[2] << endl;
115   AliL3Vertex v;
116   //v.SetX(vertexpos[0]);
117   //v.SetY(vertexpos[1]);
118   //v.SetZ(vertexpos[2]);
119
120   for (Int_t i = 0;i<ntr; i++) {
121     AliESDHLTtrack *kesdtrack;
122     if(fTrackerType){
123       kesdtrack=esd->GetHLTConfMapTrack(i);
124     } else {
125       kesdtrack=esd->GetHLTHoughTrack(i);
126     }
127
128     if (kesdtrack == 0)
129       {
130         Error("ReadESD","Can not get track %d", i);
131         continue;
132       }
133
134     //const Float_t kpid=kesdtrack->GetPID();
135     const Int_t knhits=kesdtrack->GetNHits();
136     const Int_t kweight=kesdtrack->GetWeight();
137     //cout << i << " " << kweight << " " << knhits << endl;
138     if((fMinHits>0) && (knhits<fMinHits)) continue;    
139     if(kweight>1000) continue; //avoid ghosts 
140     if((fMinWeight>0) && (kweight<fMinWeight)) continue;
141
142     Float_t px=kesdtrack->GetPx();
143     Float_t py=kesdtrack->GetPy();
144     Float_t pz=kesdtrack->GetPz();
145     //cout << px << " " << py << " " << pz <<  " " << TMath::Sqrt(px*px+py*py) << endl;
146
147     if(0&&fTrackerType){
148 #if 0
149       try {
150         Double_t mom[3];
151         AliL3ITStrack l3(*kesdtrack,vertexpos[2]);
152         if(!l3.GetPxPyPzAt(0,mom)) continue;
153         px=mom[0];
154         py=mom[1];
155         pz=mom[2];
156       } catch (const Char_t *msg) {
157         //Warning("Clusters2Tracks",msg);
158         continue;
159       }
160 #else
161       AliL3Track l3;
162       //if(!kesdtrack->ComesFromMainVertex()) continue;
163       //cout << "Pos: " << kesdtrack->GetFirstPointX() << " " << kesdtrack->GetFirstPointY() << " " << kesdtrack->GetFirstPointZ() << endl;      
164       l3.SetFirstPoint(kesdtrack->GetFirstPointX(),kesdtrack->GetFirstPointY(),kesdtrack->GetFirstPointZ());
165       l3.SetLastPoint(kesdtrack->GetLastPointX(),kesdtrack->GetLastPointY(),kesdtrack->GetLastPointZ());
166       l3.SetCharge(kesdtrack->GetCharge());
167       l3.SetPt(kesdtrack->GetPt());
168       l3.SetTgl(kesdtrack->GetTgl());
169       l3.SetPsi(kesdtrack->GetPsi());
170       l3.CalculateHelix();
171       Double_t xc=0.,yc=0.,zc=0.;
172       l3.GetClosestPoint(&v,xc,yc,zc);
173       if(TMath::Abs(zc)>10.) continue;
174       l3.SetFirstPoint(xc,yc,zc);
175       //cout << "Pos: " << xc << " " << yc << " " << zc << endl;
176       l3.UpdateToFirstPoint();
177       px=l3.GetPx();
178       py=l3.GetPy();
179       pz=l3.GetPz();
180 #endif
181     }
182     const Float_t kpt=TMath::Sqrt(px*px+py*py);
183     //cout << px << " " << py << " " << pz <<  " " << kpt << endl;
184
185     const Float_t kp=TMath::Sqrt(pz*pz+kpt*kpt);
186     const Float_t keta=0.5*TMath::Log((kp+pz+1e-30)/(kp-pz+1e-30)); 
187     const Float_t kphi=TMath::Pi()+TMath::ATan2(-py,-px);
188
189     if(IsAcceptedParticle(kpt,kphi,keta))
190       fEventParticles->AddParticle(px,py,pz,kp,i,kesdtrack->GetMCid(),knhits,kpt,kphi,keta);
191   } //loop over tracks
192
193   return kTRUE;
194 }