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