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