HLT/TPCLib/AliHLTTPCTrackHistoComponent.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: Gaute Ovrebekk <ovrebekk@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   AliHLTTPCTrackHistoComponent.cxx
  20     @author Gaute Ovrebekk, Matthias Richter
  21     @date
  22     @brief  The TPC conformal mapping tracker component.
  23 */
  24
  25
  26 #if __GNUC__>= 3
  27 using namespace std;
  28 #endif
  29
  30 #include "AliHLTTPCTrackHistoComponent.h"
  31 #include "AliHLTTPCTransform.h"
  32 #include "AliHLTTPCClusterDataFormat.h"
  33 #include "AliHLTTPCTrackletDataFormat.h"
  34 #include "AliHLTTPCMemHandler.h"
  35 #include "AliHLTTPCDefinitions.h"
  36 #include <TFile.h>
  37 #include <TString.h>
  38 #include "TNtuple.h"
  39 #include "TObjString.h"
  40 #include "TObjArray.h"
  41 #include "AliHLTTPCTrackArray.h"
  42 #include "AliHLTTPCTrack.h"
  43
  44 //#include "AliHLTTPC.h"
  45 //#include <stdlib.h>
  46 //#include <cerrno>
  47
  48 AliHLTTPCTrackHistoComponent gAliHLTTPCTrackHistoComponent;
  49
  50 /** ROOT macro for the implementation of ROOT specific class methods */
  51 ClassImp(AliHLTTPCTrackHistoComponent)
  52
  53 AliHLTTPCTrackHistoComponent::AliHLTTPCTrackHistoComponent()
  54 :
  55 fClusters(NULL),
  56   fTracks(NULL),
  57   fTracksArray(NULL)
  58 {
  59
  60   // see header file for class documentation
  61   // or
  62   // refer to README to build package
  63   // or
  64   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  65
  66 }
  67
  68 AliHLTTPCTrackHistoComponent::~AliHLTTPCTrackHistoComponent()
  69 {
  70   // see header file for class documentation
  71 }
  72
  73 // Public functions to implement AliHLTComponent's interface.
  74 // These functions are required for the registration process
  75
  76 const char* AliHLTTPCTrackHistoComponent::GetComponentID()
  77 {
  78   // see header file for class documentation
  79
  80   return "TPCTrackHisto";
  81 }
  82
  83 void AliHLTTPCTrackHistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  84 {
  85   // see header file for class documentation
  86   list.clear();
  87   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
  88   list.push_back( AliHLTTPCDefinitions::fgkTrackSegmentsDataType );
  89   list.push_back( AliHLTTPCDefinitions::fgkTracksDataType );
  90 }
  91
  92 AliHLTComponentDataType AliHLTTPCTrackHistoComponent::GetOutputDataType()
  93 {
  94   // see header file for class documentation
  95   return kAliHLTDataTypeTNtuple;
  96
  97 }
  98
  99 void AliHLTTPCTrackHistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 100 {
 101   // see header file for class documentation
 102   // XXX TODO: Find more realistic values.
 103   constBase = 0;
 104   inputMultiplier = 1;
 105 }
 106
 107 AliHLTComponent* AliHLTTPCTrackHistoComponent::Spawn()
 108 {
 109   // see header file for class documentation
 110   return new AliHLTTPCTrackHistoComponent;
 111 }
 112
 113 int AliHLTTPCTrackHistoComponent::DoInit( int argc, const char** argv )
 114 {
 115
 116   fClusters = new TNtuple("fCluster", "fCluster", "charge:qmax:residualY:residualZ:used:event");
 117   fTracks = new TNtuple("fTracks", "fTracks", "pt:eta:psi:nclusters:event");
 118   fTracksArray=new AliHLTTPCTrackArray();
 119
 120   int iResult=0;
 121   TString configuration="";
 122   TString argument="";
 123   for (int i=0; i<argc && iResult>=0; i++) {
 124     argument=argv[i];
 125     if (!configuration.IsNull()) configuration+=" ";
 126     configuration+=argument;
 127   }
 128
 129   if (!configuration.IsNull()) {
 130     iResult=Configure(configuration.Data());
 131   }
 132   return iResult;
 133 }
 134
 135 int AliHLTTPCTrackHistoComponent::DoDeinit()
 136 {
 137   // see header file for class documentation
 138
 139   delete fClusters;
 140   delete fTracks;
 141   delete fTracksArray;
 142
 143   return 0;
 144 }
 145
 146 int AliHLTTPCTrackHistoComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/)
 147 {
 148   if(!fTracksArray){fTracksArray=new AliHLTTPCTrackArray();}
 149
 150   const AliHLTComponentBlockData* iter = NULL;
 151
 152   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
 153     return 0;
 154
 155   Int_t TotalTrack = 0;
 156
 157   //Reading Merged Tracks
 158   for ( iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkTracksDataType); iter != NULL; iter = GetNextInputBlock() ) {
 159     if(iter->fDataType!=AliHLTTPCDefinitions::fgkTracksDataType){continue;}
 160     ReadTracks(iter,TotalTrack);
 161   }
 162
 163   //Reading Tracks form slice
 164   for ( iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkTrackSegmentsDataType); iter != NULL; iter = GetNextInputBlock() ) {
 165     if(iter->fDataType!=AliHLTTPCDefinitions::fgkTrackSegmentsDataType){continue;}
 166     ReadTracks(iter,TotalTrack);
 167   }
 168
 169   int TotalSpacePoint = 0;
 170   int nClustersUsed=0;
 171
 172   for ( iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkClustersDataType); iter != NULL; iter = GetNextInputBlock() ) {
 173
 174     if(iter->fDataType!=AliHLTTPCDefinitions::fgkClustersDataType){continue;}
 175
 176     AliHLTUInt8_t slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
 177     AliHLTUInt8_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
 178
 179     //HLTDebug ( "Input Data - TPC cluster - Slice/Patch: %d/%d.", slice, patch );
 180     const AliHLTTPCClusterData* clusterData = (const AliHLTTPCClusterData*) iter->fPtr;
 181     Int_t nSpacepoint = (Int_t) clusterData->fSpacePointCnt;
 182     TotalSpacePoint += nSpacepoint;
 183     //HLTInfo("TrackHisto found %d Spacepoints in slice %d patch %d", nSpacepoint, slice, patch);
 184     AliHLTTPCSpacePointData *clusters = (AliHLTTPCSpacePointData*) clusterData->fSpacePoints;
 185
 186     if (fClustersArray[slice][patch]!=NULL) {
 187       delete(fClustersArray[slice][patch]);
 188       fClustersArray[slice][patch]=NULL;
 189     }
 190     Int_t arraysize=nSpacepoint*sizeof(AliHLTTPCSpacePointData);
 191     fClustersArray[slice][patch] = (AliHLTTPCSpacePointData*)new Byte_t[arraysize];
 192     if (fClustersArray[slice][patch]) {
 193       memcpy(fClustersArray[slice][patch], clusters, arraysize);
 194       fNcl[slice][patch]=nSpacepoint;
 195     } else {
 196       fNcl[slice][patch]=nSpacepoint;
 197       HLTError ( "Memory allocation failed!" );
 198     }
 199
 200     for(int i=0;i<nSpacepoint;i++){
 201       UInt_t idCluster = clusters[i].fID;
 202       Int_t sliceCl = (idCluster>>25) & 0x7f;
 203       Int_t patchCl = (idCluster>>22) & 0x7;
 204       UInt_t pos = idCluster&0x3fffff;
 205       Int_t used = 0;
 206       Float_t resy = 0, resz = 0;
 207       for(UInt_t id=0;id<fTrackClusterID[sliceCl][patchCl].size();id++){
 208         if(fTrackClusterID[sliceCl][patchCl][id]==pos){
 209           clusters[i].fUsed=kTRUE;
 210           nClustersUsed++;
 211           used=1;
 212           FillResidual(pos,sliceCl,patchCl,resy,resz);
 213         }
 214       }
 215       if(used==1){
 216         fClusters->Fill(clusters[i].fCharge,clusters[i].fQMax,resy,resz,used,GetEventId());
 217       }
 218       else{
 219         fClusters->Fill(clusters[i].fCharge,clusters[i].fQMax,-100,-100,used,GetEventId());
 220       }
 221     }
 222   }
 223
 224   HLTInfo("TrackHisto found %d Spacepoints",TotalSpacePoint);
 225   HLTInfo("TrackHisto found %d Tracks",TotalTrack);
 226
 227   PushHisto();
 228
 229   delete fTracksArray;
 230   fTracksArray=NULL;
 231
 232   for(UInt_t i=0;i<36;i++){
 233     for(UInt_t j=0;j<6;j++){
 234       fTrackClusterID[i][j].clear();
 235     }
 236   }
 237
 238   return 0;
 239 }
 240
 241  int AliHLTTPCTrackHistoComponent::Configure(const char* arguments)
 242  {
 243
 244    int iResult=0;
 245    if (!arguments) return iResult;
 246
 247    TString allArgs=arguments;
 248    TString argument;
 249
 250    TObjArray* pTokens=allArgs.Tokenize(" ");
 251    if (pTokens) {
 252      for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 253        argument=((TObjString*)pTokens->At(i))->GetString();
 254        if (argument.IsNull()) continue;
 255
 256        HLTError("unknown argument %s", argument.Data());
 257        iResult=-EINVAL;
 258        break;
 259      }
 260      delete pTokens;
 261    }
 262
 263    return iResult;
 264  }
 265
 266 void AliHLTTPCTrackHistoComponent::ReadTracks(const AliHLTComponentBlockData* iter,Int_t &tt){
 267
 268   AliHLTUInt8_t slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
 269   //AliHLTUInt8_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
 270
 271   //HLTDebug ( "Input Data - TPC cluster - Slice/Patch: %d/%d.", slice, patch );
 272   AliHLTTPCTrackletData* trackData = (AliHLTTPCTrackletData*) iter->fPtr;
 273   AliHLTUInt32_t nTracks = trackData->fTrackletCnt;
 274   fTracksArray->FillTracksChecked(trackData->fTracklets,trackData->fTrackletCnt,iter->fSize,slice,true);
 275   tt += nTracks;
 276   //HLTInfo("TrackHisto found %d Tracks in slice %d patch %d", nTracks, slice, patch);
 277   AliHLTTPCTrackSegmentData *tracks = (AliHLTTPCTrackSegmentData*) trackData->fTracklets;
 278
 279   for(AliHLTUInt32_t i=0;i<nTracks;i++){
 280     UInt_t nHits = tracks->fNPoints;
 281
 282     fTracks->Fill(tracks->fPt,tracks->fPsi,tracks->fTgl,nHits,GetEventId());
 283
 284     const UInt_t *hitnum = tracks->fPointIDs;
 285     for(UInt_t h=0; h<nHits; h++){
 286       UInt_t idTrack = hitnum[h];
 287       Int_t sliceTrack = (idTrack>>25) & 0x7f;
 288       Int_t patchTrack = (idTrack>>22) & 0x7;
 289       UInt_t pos = idTrack&0x3fffff;
 290       fTrackClusterID[sliceTrack][patchTrack].push_back(pos);
 291     }
 292     UChar_t *tmpP = (UChar_t*)tracks;
 293     tmpP += sizeof(AliHLTTPCTrackSegmentData)+tracks->fNPoints*sizeof(UInt_t);
 294     tracks = (AliHLTTPCTrackSegmentData*)tmpP;
 295   }
 296 }
 297
 298 void AliHLTTPCTrackHistoComponent::PushHisto(){
 299
 300     AliHLTUInt32_t fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification(0,35,0,5);
 301     PushBack( (TObject*) fTracks,kAliHLTDataTypeTNtuple, fSpecification);
 302     PushBack( (TObject*) fClusters,kAliHLTDataTypeTNtuple, fSpecification);
 303
 304 }
 305 void AliHLTTPCTrackHistoComponent::FillResidual( UInt_t pos,AliHLTUInt8_t slice,AliHLTUInt8_t patch,Float_t& resy,Float_t& resz){
 306
 307   AliHLTTPCSpacePointData *cl =  &fClustersArray[slice][patch][pos];
 308   if(!cl){return;}
 309
 310   AliHLTTPCTrack *gtrack = NULL;
 311
 312   for(int i;i<fTracksArray->GetNTracks();i++){
 313     AliHLTTPCTrack *tt = fTracksArray->GetCheckedTrack(i);
 314     UInt_t *hitnum =tt->GetHitNumbers();
 315     Int_t nHits = tt->GetNHits();
 316     for(Int_t h=0; h<nHits; h++){
 317       UInt_t id=hitnum[h];
 318       Int_t Tslice = (id>>25) & 0x7f;
 319       Int_t Tpatch = (id>>22) & 0x7;
 320       UInt_t Tpos = id&0x3fffff;
 321       if(Tslice==slice && Tpatch==patch && Tpos==pos) {
 322         gtrack = tt;
 323         break;
 324       }
 325     }
 326   }
 327
 328   if(!gtrack){return;}
 329
 330   Int_t tslice = gtrack->GetSector();
 331   Double_t radius = gtrack->GetRadius();      // radius
 332   Double_t kappa = gtrack->GetKappa();        // curvature = 1/R , signed
 333   Double_t lambda = atan( gtrack->GetTgl() ); // dipAngle lambda
 334
 335   // ------------------------------------
 336   // ++ Get first/last point of the track
 337
 338   Double_t xyzL[3];      // lastpoint of track
 339   Double_t xyzF[3];      // firstpoint of track
 340
 341   xyzF[0] = gtrack->GetFirstPointX();
 342   xyzF[1] = gtrack->GetFirstPointY();
 343   xyzF[2] = gtrack->GetFirstPointZ();
 344
 345   xyzL[0] = gtrack->GetLastPointX();
 346   xyzL[1] = gtrack->GetLastPointY();
 347   xyzL[2] = gtrack->GetLastPointZ();
 348
 349   // --------------------------
 350   // ++ Calculate length of the track
 351
 352   Double_t s = 0.;       // length of the track
 353   if (  AliHLTTPCTransform::GetBFieldValue() == 0. || kappa == 0 )
 354     s = sqrt ( (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]) );
 355   else {
 356     // Calculate the length of the track. If it is to flat in in s,z plane use sxy, otherwise use sz
 357     if (fabs(lambda) > 0.05){
 358       // length of track calculated out of z
 359       s = fabs( (xyzL[2] - xyzF[2]) / sin(lambda) ); // length of track calculated out of z
 360     }
 361     else {
 362       Double_t d = (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]);
 363       // length of track calculated out of xy
 364       s = fabs ( acos( 0.5 * (2 - (d / (radius*radius)))) / ( kappa * cos(lambda) ) );
 365     }
 366   }
 367
 368   gtrack->Rotate(tslice,kTRUE);
 369
 370   Double_t padrows = 0;
 371
 372   Float_t xyzC[3];       // cluster tmp
 373   Float_t xyzTtmp[3];    // track tmp
 374
 375   xyzC[0] = cl->fX;
 376   xyzC[1] = cl->fY;
 377   xyzC[2] = cl->fZ;
 378
 379   Int_t padrow = AliHLTTPCTransform::GetPadRow(cl->fX);
 380
 381   xyzTtmp[0] = gtrack->GetFirstPointX();
 382
 383   if(gtrack->GetCrossingPoint(padrow,xyzTtmp)) {
 384     // ----------------------
 385        // ++ Calculate Residuals
 386
 387        Float_t deltaY = ( xyzC[1] - xyzTtmp[1] );
 388        Float_t deltaZ = ( xyzC[2] - xyzTtmp[2] );
 389
 390        resy = deltaY;
 391        resz = deltaZ;
 392   }
 393   else{
 394     resy = -1000;
 395     resz = -1000;
 396   }
 397 }