HLT/TPCLib/AliHLTTPCTrackGeometry.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE HLT Project        *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  for The ALICE HLT Project.                            *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /// @file   AliHLTTPCTrackGeometry.cxx
  20 /// @author Matthias Richter
  21 /// @date   2011-05-20
  22 /// @brief  Desciption of a track by a sequence of track points
  23 ///
  24
  25 #include "AliHLTTPCTrackGeometry.h"
  26 #include "AliHLTTPCTransform.h"
  27 #include "AliHLTTPCSpacePointData.h"
  28 #include "AliHLTTPCClusterDataFormat.h"
  29 #include "AliHLTTPCSpacePointContainer.h"
  30 #include "AliHLTTPCDefinitions.h"
  31 #include "AliHLTComponent.h"
  32 #include "AliHLTGlobalBarrelTrack.h"
  33 #include "TMath.h"
  34 #include "TH2F.h"
  35 #include <memory>
  36
  37 /** ROOT macro for the implementation of ROOT specific class methods */
  38 ClassImp(AliHLTTPCTrackGeometry)
  39
  40 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
  41   : AliHLTTrackGeometry()
  42   , fRawTrackPoints()
  43 {
  44   /// standard constructor
  45 }
  46
  47 AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
  48   : AliHLTTrackGeometry(src)
  49   , fRawTrackPoints(src.fRawTrackPoints)
  50 {
  51   /// copy constructor
  52 }
  53
  54 AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
  55 {
  56   /// assignment operator
  57   AliHLTTrackGeometry::operator=(src);
  58   fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
  59   return *this;
  60 }
  61
  62 AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
  63 {
  64   /// destructor
  65 }
  66
  67 float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
  68 {
  69   /// alpha of the plane
  70   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
  71   float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
  72   if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
  73   return alpha;
  74 }
  75
  76 float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
  77 {
  78   /// radial distance from global {0,0,0}
  79   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
  80   UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
  81   Int_t row=AliHLTTPCTransform::GetFirstRow(partition)+number;
  82   return AliHLTTPCTransform::Row2X(row);
  83 }
  84
  85 float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
  86 {
  87   /// theta of the plane
  88   return 0.0;
  89 }
  90
  91 bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
  92 {
  93   /// check bounds in u and v coordinate
  94   float r=GetPlaneR(planeId);
  95   if (r<AliHLTTPCTransform::GetFirstRow(0)) return false;
  96
  97   // TODO: check if the pad width needs to be considered here
  98   return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
  99 }
 100
 101 int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
 102 {
 103   /// calculate the track points, expects the global magnetic field to be initialized
 104   AliHLTGlobalBarrelTrack bt(track);
 105   return CalculateTrackPoints(bt);
 106 }
 107
 108 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
 109 {
 110   /// calculate the track points, expects the global magnetic field to be initialized
 111   int iResult=0;
 112   int firstpadrow=0;
 113   for (;
 114        firstpadrow<AliHLTTPCTransform::GetNRows() &&
 115          AliHLTTPCTransform::Row2X(firstpadrow)+AliHLTTPCTransform::GetPadLength(firstpadrow)<track.GetX();
 116        firstpadrow++);
 117   if (firstpadrow>=AliHLTTPCTransform::GetNRows()) return 0;
 118   iResult=CalculateTrackPoints(track, firstpadrow, 1);
 119   if (iResult>=0 && firstpadrow>0)
 120     iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
 121   return iResult;
 122 }
 123
 124 int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
 125 {
 126   /// calculate the track points, expects the global magnetic field to be initialized
 127   float offsetAlpha=0.0;
 128   for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCTransform::GetNRows(); padrow+=step) {
 129     float x=AliHLTTPCTransform::Row2X(padrow);
 130     float y=0.0;
 131     float z=0.0;
 132
 133     int maxshift=9;
 134     int shift=0;
 135     int result=0;
 136     do {
 137       // start calculation of crossing points with padrow planes in the slice of the first point
 138       // plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
 139       // is out of bounds
 140       if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
 141       float pointAlpha=TMath::ATan(y/x);
 142       if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
 143         offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
 144         result=0;
 145       }
 146     } while (result==0 && shift++<maxshift);
 147     if (result<1) continue;
 148     float planealpha=track.GetAlpha()-offsetAlpha;
 149     if (planealpha<0) planealpha+=TMath::TwoPi();
 150     int slice=int(9*planealpha/TMath::Pi());
 151     //if (z<0) slice+=18;
 152     int partition=AliHLTTPCTransform::GetPatch(padrow);
 153     int row=padrow-AliHLTTPCTransform::GetFirstRow(partition);
 154     UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
 155     if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
 156       HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f), expected %f (%.0f)", id, slice, GetPlaneAlpha(id), 180*GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180*planealpha/TMath::Pi());
 157     }
 158     if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
 159       Float_t rpt[3]={0.,y,z}; // row pad time
 160       AliHLTTPCTransform::LocHLT2Raw(rpt, slice, padrow);
 161       // FIXME: there is a mismatch in the definition of the pad coordinate
 162       // should be with respect to middle of pad, that's why the offset of
 163       // 0.5 has been applied when calling the AliHLTTPCClusterTransformation
 164       // and for conversion to AliTPCclusterMI in AliHLTTPCClusterAccessHLTOUT
 165       // AliHLTTPCTransform::LocHLT2Raw seems to define this shift in the
 166       // opposite direction
 167       rpt[1]+=1.;
 168       fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
 169     }
 170   }
 171   return 0;
 172 }
 173
 174 int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
 175 {
 176   /// find the track point which can be associated to a spacepoint with coordinates and id
 177   UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
 178   UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
 179   int row=AliHLTTPCTransform::GetPadRow(spacepoint[0]);
 180   bool bSpecialRow=row==30 || row==90 || row==139;
 181   if (row<AliHLTTPCTransform::GetFirstRow(partition) || row>AliHLTTPCTransform::GetLastRow(partition)) {
 182     HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
 183     return -EINVAL;
 184   }
 185
 186   // find the crossing point of the track with the padrow plane where
 187   // the spacepoint is
 188   // 1) calculate plane id from slice, partition and row (within partition)
 189   row-=AliHLTTPCTransform::GetFirstRow(partition);
 190   UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
 191   const AliHLTTrackPoint* point=GetTrackPoint(id);
 192   // track might be outside the partition and cross the central membrane
 193   // search in the other half of the TPC
 194   if (!point && slice<18) {
 195     // search in the neighboring partition on the C side
 196     id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
 197     point=GetTrackPoint(id);
 198   } else if (!point && slice>=18) {
 199     // search in the neighboring partition on the A side
 200     id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
 201     point=GetTrackPoint(id);
 202   }
 203
 204   // search in the neighboring partition, this takes account for rows
 205   // 30, 90, and 139 which are partly in one and the other partition
 206   if (!point && bSpecialRow) {
 207     row+=AliHLTTPCTransform::GetFirstRow(partition);
 208     row-=AliHLTTPCTransform::GetFirstRow(partition-1);
 209     id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
 210     point=GetTrackPoint(id);
 211     if (!point && slice<18) {
 212       // search in the neighboring partition on the C side
 213       id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
 214       point=GetTrackPoint(id);
 215     } else if (!point && slice>=18) {
 216       // search in the neighboring partition on the A side
 217       id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
 218       point=GetTrackPoint(id);
 219     }
 220   }
 221
 222   if (point) {
 223     planeId=id;
 224     if (point->HaveAssociatedSpacePoint()) {
 225       if (GetVerbosity()>2) HLTInfo("descarding spacepoint 0x%08x z=%f y=%f z=%f: track point 0x%08x already occupied", spacepoint[0], spacepoint[1], spacepoint[2], planeId);
 226       return 0; // already occupied
 227     }
 228     float maxdy=2.;
 229     float maxdz=2.;
 230     if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
 231       if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x y %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetU(), maxdy);
 232       return -ENOENT;
 233     }
 234     if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
 235       if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x z %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetV(), maxdz);
 236       return -ENOENT;
 237     }
 238     return 1;
 239   }
 240   return -ENOENT;
 241 }
 242
 243 AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
 244 {
 245   /// create a collection of all points
 246   std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
 247   if (!spacepoints.get()) return NULL;
 248
 249   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
 250   unsigned i=0;
 251   while (i<trackPoints.size()) {
 252     // allocate buffer for all points, even though the buffer might not be filled
 253     // completely because of a partition change
 254     int nofPoints=trackPoints.size()-i;
 255     int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
 256     AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
 257     if (!pBuffer) return NULL;
 258     AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
 259     pClusterData->fSpacePointCnt=0;
 260     AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
 261     int currentSlice=-1;
 262     int currentPartition=-1;
 263     for (; i<trackPoints.size(); i++) {
 264       if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
 265       AliHLTUInt32_t planeId=trackPoints[i].GetId();
 266       int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
 267       int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
 268       int number=AliHLTTPCSpacePointData::GetNumber(planeId);
 269       if ((currentSlice>=0 && currentSlice!=slice) || (currentPartition>=0 && currentPartition!=partition)) {
 270         // change of partition or slice, need to go to next block
 271         // 2011-07-26 currently all spacepoints go into one block, if separated
 272         // blocks per partition are needed one has to leave the inner loop here
 273         // and set the data block specification below
 274         // Caution: not tested, only the last block seems to make it through
 275         //break;
 276       }
 277       currentSlice=slice;
 278       currentPartition=partition;
 279       pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
 280       pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
 281       pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
 282       pClusters[pClusterData->fSpacePointCnt].fID=planeId;
 283       pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCTransform::GetFirstRow(partition)+number;
 284       pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
 285       pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
 286       pClusters[pClusterData->fSpacePointCnt].fCharge=0;
 287       pClusters[pClusterData->fSpacePointCnt].fQMax=0;
 288       pClusters[pClusterData->fSpacePointCnt].fUsed=0;
 289       pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
 290       pClusterData->fSpacePointCnt++;
 291     }
 292     AliHLTComponentBlockData bd;
 293     AliHLTComponent::FillBlockData(bd);
 294     bd.fPtr=pBuffer;
 295     bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
 296     AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
 297     bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
 298     spacepoints->AddInputBlock(&bd);
 299   }
 300
 301   return spacepoints.release();
 302 }
 303
 304 int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
 305 {
 306   // fill residual histogram
 307   if (!histo || !points) return -EINVAL;
 308   const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
 309   for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
 310        trackpoint!=trackPoints.end(); trackpoint++) {
 311     if (!trackpoint->HaveAssociatedSpacePoint()) continue;
 312     for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
 313          sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
 314     AliHLTUInt32_t spacepointId=sp->fId;
 315     vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
 316     if (rawpoint==fRawTrackPoints.end()) {
 317       HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
 318       continue;
 319     }
 320     if (!points->Check(spacepointId)) {
 321       //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
 322       continue;
 323     }
 324     float value=0.;
 325     if (coordinate==0) {
 326       value=rawpoint->GetU()-points->GetY(spacepointId);
 327       histo->Fill(GetPlaneR(trackpoint->GetId()), value);
 328     } else {
 329       value=rawpoint->GetV()-points->GetZ(spacepointId);
 330       //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
 331       histo->Fill(rawpoint->GetV(), value);
 332     }
 333     }
 334   }
 335   return 0;
 336 }