HLT/global/AliHLTGlobalTrackResidualsComponent.cxx

   1 // $Id$
   2 //**************************************************************************
   3 //* This file is property of and copyright by the ALICE HLT Project        *
   4 //* ALICE Experiment at CERN, All rights reserved.                         *
   5 //*                                                                        *
   6 //* Primary Authors: Timur Pocheptsov <Timur.Pocheptsov@cern.ch>           *
   7 //*                  for The ALICE HLT Project.                            *
   8 //*                                                                        *
   9 //* Permission to use, copy, modify and distribute this software and its   *
  10 //* documentation strictly for non-commercial purposes is hereby granted   *
  11 //* without fee, provided that the above copyright notice appears in all   *
  12 //* copies and that both the copyright notice and this permission notice   *
  13 //* appear in the supporting documentation. The authors make no claims     *
  14 //* about the suitability of this software for any purpose. It is          *
  15 //* provided "as is" without express or implied warranty.                  *
  16 //**************************************************************************
  17
  18 /// @file  AliHLTGlobalTrackResidualsComponent.cxx
  19 /// @author Timur Pocheptsov
  20 /// @date
  21 /// @brief A histogramming component for plotting the Y and Z track residual
  22 ///
  23
  24 #if __GNUC__>= 3
  25 using namespace std;
  26 #endif
  27
  28 #include <algorithm>
  29
  30 #include <TMath.h>
  31
  32 #include "AliHLTGlobalTrackResidualsComponent.h"
  33 #include "AliHLTTPCClusterDataFormat.h"
  34 #include "AliHLTTPCSpacePointData.h"
  35 #include "AliHLTGlobalBarrelTrack.h"
  36 #include "AliHLTTPCDefinitions.h"
  37 #include "AliHLTDataTypes.h"
  38
  39 ClassImp(AliHLTGlobalTrackResidualsComponent)
  40
  41 //_______________________________________________________________________________________________
  42
  43 AliHLTGlobalTrackResidualsComponent::AliHLTGlobalTrackResidualsComponent()
  44                                 : AliHLTProcessor(),
  45                                   fResY("y_residuals", "y residuals", kNBins, -1., 1.),
  46                                   fResZ("z_residuals", "z residuals", kNBins, -1., 1.),
  47                                   fSortedX()
  48 {
  49   //Ctor.
  50   fResY.SetMarkerStyle(8);
  51   fResY.SetMarkerSize(0.4);
  52   fResY.SetXTitle("Y [cm]");
  53   fResY.SetDirectory(0);
  54
  55   fResZ.SetMarkerStyle(8);
  56   fResZ.SetMarkerSize(0.4);
  57   fResZ.SetXTitle("Z [cm]");
  58   fResZ.SetDirectory(0);
  59
  60   CleanClusters();
  61 }
  62
  63 //_______________________________________________________________________________________________
  64 const char * AliHLTGlobalTrackResidualsComponent::GetComponentID()
  65 {
  66   //Component's name.
  67   return "GlobalTrackResiduals";
  68 }
  69
  70 //_______________________________________________________________________________________________
  71 void AliHLTGlobalTrackResidualsComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  72 {
  73   //Possible input data types
  74   list.clear();
  75   list.push_back(AliHLTTPCDefinitions::fgkClustersDataType|kAliHLTDataOriginTPC);
  76   list.push_back(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
  77 }
  78
  79 //_______________________________________________________________________________________________
  80 AliHLTComponentDataType AliHLTGlobalTrackResidualsComponent::GetOutputDataType()
  81 {
  82   //Output's data type(s).
  83   return kAliHLTMultipleDataType;
  84 }
  85
  86 //_______________________________________________________________________________________________
  87 int AliHLTGlobalTrackResidualsComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
  88 {
  89   //Output's data type(s).
  90   tgtList.clear();
  91   tgtList.push_back(kAliHLTDataTypeTNtuple|kAliHLTDataOriginTPC);
  92   tgtList.push_back(kAliHLTDataTypeHistogram|kAliHLTDataOriginTPC);
  93   return tgtList.size();
  94 }
  95
  96 //_______________________________________________________________________________________________
  97 void AliHLTGlobalTrackResidualsComponent::GetOutputDataSize(unsigned long& constBase, double& inputMultiplier)
  98 {
  99   //Approximate output histograms sizes.
 100   constBase  = sizeof(TH1F) * 2; //Size of histogram objects.
 101   constBase += sizeof(Float_t) * kNBins * 2; //Size of memory, allocated by 1D histograms.
 102   //Forget about strings (name and title), just multiply by 2.
 103   constBase *= 2;
 104   inputMultiplier = 1;
 105 }
 106
 107 //_______________________________________________________________________________________________
 108 AliHLTComponent* AliHLTGlobalTrackResidualsComponent::Spawn()
 109 {
 110   //Create the component.
 111   return new AliHLTGlobalTrackResidualsComponent;
 112 }
 113
 114 //_______________________________________________________________________________________________
 115 int AliHLTGlobalTrackResidualsComponent::DoInit(int /*argc*/, const char** /*argv*/)
 116 {
 117   //(Re)Initialize component.
 118   ResetHistograms();
 119   return 0;
 120 }
 121
 122 //_______________________________________________________________________________________________
 123 int AliHLTGlobalTrackResidualsComponent::DoDeinit()
 124 {
 125   //DoNothing will be better name.
 126   return 0;
 127 }
 128
 129 //_______________________________________________________________________________________________
 130 int AliHLTGlobalTrackResidualsComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/)
 131 {
 132   //Process global barrel tracks and clusters, calculate residuals.
 133   if (GetFirstInputBlock(kAliHLTDataTypeSOR) || GetFirstInputBlock(kAliHLTDataTypeEOR))
 134     return 0;
 135
 136   //Read input data, find residuals, fill histgrams.
 137   ProcessBlocks();
 138
 139   //Do output now.
 140   PushBack(&fResY, kAliHLTDataTypeHistogram | kAliHLTDataOriginTPC, 0);
 141   PushBack(&fResZ, kAliHLTDataTypeHistogram | kAliHLTDataOriginTPC, 0);
 142
 143   return 0;
 144 }
 145
 146 //_______________________________________________________________________________________________
 147 void AliHLTGlobalTrackResidualsComponent::ProcessBlocks()
 148 {
 149   //1. Read cluster blocks.
 150   ReadClusterBlocks();
 151   //2. Loop over merged tracks, calculate residuals, fill histogramms.
 152   std::vector<AliHLTGlobalBarrelTrack> ts;
 153
 154   Int_t totalTracks = 0;
 155   const AliHLTComponentBlockData * i = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
 156
 157   for (; i; i = GetNextInputBlock()) {
 158     if (i->fDataType != (kAliHLTDataTypeTrack | kAliHLTDataOriginTPC))
 159       continue;
 160
 161     ts.clear();
 162     AliHLTGlobalBarrelTrack::ConvertTrackDataArray((AliHLTTracksData*)i->fPtr, i->fSize, ts);
 163
 164     totalTracks += Int_t(ts.size());
 165
 166     std::vector<AliHLTGlobalBarrelTrack>::size_type j = 0, e = ts.size();
 167     for (; j != e; ++j) {
 168       fSortedX.clear();
 169       SortHitsX(ts[j]);
 170       FillResiduals(ts[j]);
 171     }
 172
 173     HLTDebug("TrackResiduals found %d tracks", totalTracks);
 174   }
 175 }
 176
 177 //_______________________________________________________________________________________________
 178 void AliHLTGlobalTrackResidualsComponent::ReadClusterBlocks()
 179 {
 180   //Loop over blocks, find cluster blocks, extract space points.
 181   CleanClusters();
 182
 183   Int_t totalSpacePoints = 0;
 184   const AliHLTComponentBlockData* iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkClustersDataType);
 185
 186   for (; iter; iter = GetNextInputBlock()) {
 187     if (iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType)
 188       continue;
 189
 190     const AliHLTUInt8_t minSlice     = AliHLTTPCDefinitions::GetMinSliceNr(*iter);
 191     const AliHLTUInt8_t minPartition = AliHLTTPCDefinitions::GetMinPatchNr(*iter);
 192
 193     const AliHLTTPCClusterData * clusterData = (AliHLTTPCClusterData *)iter->fPtr;
 194     const Int_t nSpacepoint = Int_t(clusterData->fSpacePointCnt);
 195     totalSpacePoints += nSpacepoint;
 196
 197     //This part is from AliHLTTPCTrackHistoComponent. Logic is not clear -
 198     //is it possible that I can have two blocks with same minSlice and minPartition???
 199     //and one of them with 0 spacepoint?
 200     if (nSpacepoint) {
 201       HLTDebug("TrackResiduals component found %d spacepoints in slice %d partition %d", nSpacepoint, minSlice, minPartition);
 202       fClustersArray[minSlice][minPartition] = (AliHLTTPCSpacePointData*)clusterData->fSpacePoints;
 203       fNSpacePoints[minSlice][minPartition]  = nSpacepoint;
 204     }
 205   }
 206
 207   HLTDebug("TrackResiduals found %d spacepoints", totalSpacePoints);
 208 }
 209
 210 namespace {
 211
 212 void Rotate(Float_t* xy, Float_t alpha);
 213 Bool_t CmpX(const std::pair<Float_t, UInt_t>& rhs, const std::pair<Float_t, UInt_t>& lhs);
 214
 215 }
 216
 217 //_______________________________________________________________________________________________
 218 void AliHLTGlobalTrackResidualsComponent::SortHitsX(const AliHLTGlobalBarrelTrack& gt)
 219 {
 220   //Extract hits' Xs for the track gt, sort them.
 221   fSortedX.clear();
 222
 223   const UInt_t * hitnum = gt.GetPoints();
 224   Int_t prevSlice = -1;
 225   Float_t rotAngle = 0.f;
 226
 227   for (UInt_t i = 0; i < gt.GetNumberOfPoints(); ++i) {
 228     const UInt_t idTrack    = hitnum[i];
 229     const UInt_t pos        = idTrack & 0x3fffff;
 230     const Int_t  sliceTrack = (idTrack >> 25) & 0x7f;
 231     const UInt_t patchTrack = (idTrack >> 22) & 0x7;
 232
 233     if (!fClustersArray[sliceTrack][patchTrack])
 234       continue;
 235
 236     //The following conditional is from the original code.
 237     if (sliceTrack > 36 || patchTrack > 5) {
 238       HLTError("Corrupted TPC cluster Id: slice %d, patch %d, cluster %d", sliceTrack, patchTrack, idTrack);
 239       continue;
 240     }
 241
 242     if (fNSpacePoints[sliceTrack][patchTrack] <= pos) {
 243       HLTError("Space point array out of boundaries!");
 244       continue;
 245     }
 246
 247     if (sliceTrack != prevSlice) {
 248       if (prevSlice != -1)
 249         prevSlice < sliceTrack ? rotAngle += 0.349066 : rotAngle -= 0.349066;
 250       prevSlice = sliceTrack;
 251     }
 252
 253     Float_t clusterXY[] = {fClustersArray[sliceTrack][patchTrack][pos].fX,
 254                            fClustersArray[sliceTrack][patchTrack][pos].fY};
 255
 256     Rotate(clusterXY, rotAngle);
 257
 258     fSortedX.push_back(std::pair<Float_t, UInt_t>(clusterXY[0], i));
 259   }
 260
 261   std::sort(fSortedX.begin(), fSortedX.end(), CmpX);
 262 }
 263
 264 //_______________________________________________________________________________________________
 265 void AliHLTGlobalTrackResidualsComponent::FillResiduals(const AliHLTGlobalBarrelTrack& gt)
 266 {
 267   //Find residuals using clusters and helix approximation.
 268   const UInt_t * hitnum = gt.GetPoints();
 269   AliExternalTrackParam track(gt);
 270   Int_t prevSlice = -1;
 271   Float_t rotAngle = 0.f;
 272
 273   std::vector<std::pair<Float_t, UInt_t> >::size_type i = 0;
 274   for (; i < fSortedX.size(); ++i) {
 275     const UInt_t idTrack = hitnum[fSortedX[i].second];
 276     const UInt_t pos = idTrack & 0x3fffff;
 277     const Int_t sliceTrack = (idTrack >> 25) & 0x7f;
 278     const UInt_t patchTrack = (idTrack >> 22) & 0x7;
 279
 280     if(!fClustersArray[sliceTrack][patchTrack])
 281       continue;
 282
 283     //The following conditionals are from the original code.
 284     if (sliceTrack > 36 || patchTrack > 5) {
 285       HLTError("Corrupted TPC cluster Id: slice %d, patch %d, cluster %d", sliceTrack, patchTrack, idTrack);
 286       continue;
 287     }
 288
 289     if (fNSpacePoints[sliceTrack][patchTrack] <= pos) {
 290       HLTError("Space point array out of boundaries!");
 291       continue;
 292     }
 293
 294     if (sliceTrack != prevSlice) {
 295       if (prevSlice != -1)
 296         prevSlice < sliceTrack ? rotAngle += 0.349066 : rotAngle -= 0.349066;
 297       prevSlice = sliceTrack;
 298     }
 299
 300     if (track.PropagateTo(fSortedX[i].first, GetBz())) {
 301       Float_t clusterXYZ[] = {fClustersArray[sliceTrack][patchTrack][pos].fX,
 302                               fClustersArray[sliceTrack][patchTrack][pos].fY,
 303                               fClustersArray[sliceTrack][patchTrack][pos].fZ};
 304       Rotate(clusterXYZ, rotAngle);
 305
 306       fResY.Fill(clusterXYZ[1] - track.GetY());
 307       fResZ.Fill(clusterXYZ[2] - track.GetZ());
 308     } else
 309       break;
 310   }
 311 }
 312
 313 //_______________________________________________________________________________________________
 314 void AliHLTGlobalTrackResidualsComponent::CleanClusters()
 315 {
 316   //Set pointers and counters to zero.
 317   for (int i = 0; i < 36; ++i) {
 318     for (int j = 0; j < 6; ++j) {
 319       fClustersArray[i][j] = 0;
 320       fNSpacePoints[i][j]  = 0;
 321     }
 322   }
 323 }
 324
 325 //_______________________________________________________________________________________________
 326 void AliHLTGlobalTrackResidualsComponent::ResetHistograms()
 327 {
 328   //Set default values.
 329   fResY.Reset();
 330   fResZ.Reset();
 331
 332   fResY.SetBins(kNBins, -1., 1.);
 333   fResZ.SetBins(kNBins, -1., 1.);
 334 }
 335
 336 namespace
 337 {
 338
 339 //_______________________________________________________________________________________________
 340 void Rotate(Float_t* xy, Float_t alpha)
 341 {
 342   //From hits's local to track's _local_.
 343   const Float_t cosA = TMath::Cos(alpha);
 344   const Float_t sinA = TMath::Sin(alpha);
 345   const Float_t xPrim = xy[0] * cosA - xy[1] * sinA;
 346   const Float_t yPrim = xy[0] * sinA + xy[1] * cosA;
 347   xy[0] = xPrim;
 348   xy[1] = yPrim;
 349 }
 350
 351 //_______________________________________________________________________________________________
 352 Bool_t CmpX(const std::pair<Float_t, UInt_t>& rhs, const std::pair<Float_t, UInt_t>& lhs)
 353 {
 354   //Sort "hits" along x.
 355   return rhs.first < lhs.first;
 356 }
 357
 358
 359 }