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 }