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 #if __GNUC__>= 3
  26 using namespace std;
  27 #endif
  28
  29 #include "AliHLTTPCTrackHistoComponent.h"
  30 #include "AliHLTTPCTransform.h"
  31 #include "AliHLTTPCClusterDataFormat.h"
  32 #include "AliHLTTPCTrackletDataFormat.h"
  33 #include "AliHLTTPCMemHandler.h"
  34 #include "AliHLTTPCDefinitions.h"
  35 #include "AliHLTTPCTrackArray.h"
  36 #include "AliHLTTPCTrack.h"
  37 //#include "AliHLTGlobalBarrelTrack.h"
  38 #include "AliHLTExternalTrackParam.h"
  39 #include "AliHLTDataTypes.h"
  40
  41 #include <TFile.h>
  42 #include <TString.h>
  43 #include "TNtuple.h"
  44 #include "TH1F.h"
  45 #include "TProfile.h"
  46 #include "TObjString.h"
  47 #include "TObjArray.h"
  48
  49
  50 AliHLTTPCTrackHistoComponent gAliHLTTPCTrackHistoComponent;
  51
  52 /** ROOT macro for the implementation of ROOT specific class methods */
  53 ClassImp(AliHLTTPCTrackHistoComponent)
  54
  55 AliHLTTPCTrackHistoComponent::AliHLTTPCTrackHistoComponent()
  56   :
  57   fMinSlice(35),
  58   fMaxSlice(0),
  59   fMinPartition(5),
  60   fMaxPartition(0),
  61   //fReset(0),
  62   fNEvents(0),
  63   fNtotTracks(0),
  64   fNEvtMod(20),
  65   fMeanMultiplicity(NULL),
  66   fMultiplicity(NULL),
  67   fDeDxVsP(NULL),
  68   fClusters(NULL),
  69   fTracks(NULL),
  70   //fNClusterVsXY(NULL),
  71   //fChargeVsXY(NULL),
  72   fTracksArray(NULL)
  73   //fClustersArray(NULL),
  74   //fNSpacePoints(NULL)
  75 {
  76   // see header file for class documentation
  77   // or
  78   // refer to README to build package
  79   // or
  80   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  81 }
  82
  83 AliHLTTPCTrackHistoComponent::~AliHLTTPCTrackHistoComponent(){
  84 // see header file for class documentation
  85 }
  86
  87 // Public functions to implement AliHLTComponent's interface.
  88 // These functions are required for the registration process
  89
  90 const char* AliHLTTPCTrackHistoComponent::GetComponentID(){
  91 // see header file for class documentation
  92
  93   return "TPCTrackHisto";
  94 }
  95
  96 void AliHLTTPCTrackHistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list){
  97 // see header file for class documentation
  98   list.clear();
  99   list.push_back(AliHLTTPCDefinitions::fgkClustersDataType|kAliHLTDataOriginTPC);
 100   list.push_back(AliHLTTPCDefinitions::fgkTrackSegmentsDataType|kAliHLTDataOriginTPC);
 101   //list.push_back(AliHLTTPCDefinitions::fgkTracksDataType);
 102   list.push_back(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
 103   list.push_back(kAliHLTDataTypeTObjArray|kAliHLTDataOriginTPC);
 104 }
 105
 106 AliHLTComponentDataType AliHLTTPCTrackHistoComponent::GetOutputDataType(){
 107 // see header file for class documentation
 108   return kAliHLTMultipleDataType;
 109 }
 110
 111 int AliHLTTPCTrackHistoComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList){
 112 // see header file for class documentation
 113   tgtList.clear();
 114   tgtList.push_back(kAliHLTDataTypeTNtuple|kAliHLTDataOriginTPC);
 115   tgtList.push_back(kAliHLTDataTypeHistogram|kAliHLTDataOriginTPC);
 116   return tgtList.size();
 117 }
 118
 119 void AliHLTTPCTrackHistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ){
 120 // see header file for class documentation
 121
 122   constBase = 0;
 123   inputMultiplier = 1;// XXX TODO: Find more realistic value
 124 }
 125
 126 AliHLTComponent* AliHLTTPCTrackHistoComponent::Spawn(){
 127 // see header file for class documentation
 128   return new AliHLTTPCTrackHistoComponent;
 129 }
 130
 131 int AliHLTTPCTrackHistoComponent::DoInit( int argc, const char** argv ){
 132 // see header file for class documentation
 133
 134   fClusters = new TNtuple("fClusters", "fClusters", "charge:qmax:residualY:residualZ:event");
 135   fTracks   = new TNtuple("fTracks",  "fTracks",  "pt:eta:psi:nclusters:event");
 136   fTracksArray = new AliHLTTPCTrackArray();
 137
 138   fMultiplicity     = new TH1F("fMultiplicity",     "Track multiplicity per event",     1000,           0, 1000);
 139   fMeanMultiplicity = new TH1F("fMeanMultiplicity", "Mean track multiplicity vs. #evt", 10000/fNEvtMod, 0, 10000);
 140   fDeDxVsP          = new TProfile("fDeDxVsP",      "E-deposition per unit length vs. p",100, 0, 100);
 141
 142   int iResult = 0;
 143   TString configuration = "";
 144   TString argument = "";
 145   for(int i=0; i<argc && iResult>=0; i++){
 146       argument = argv[i];
 147       if(!configuration.IsNull()) configuration += " ";
 148       configuration += argument;
 149   }
 150
 151   if(!configuration.IsNull()){
 152      iResult = Configure(configuration.Data());
 153   }
 154   return iResult;
 155 }
 156
 157 int AliHLTTPCTrackHistoComponent::DoDeinit(){
 158 // see header file for class documentation
 159
 160   delete fClusters;
 161   delete fTracks;
 162   delete fTracksArray;
 163
 164   delete fMultiplicity;
 165   delete fMeanMultiplicity;
 166   delete fDeDxVsP;
 167
 168   return 0;
 169 }
 170
 171 int AliHLTTPCTrackHistoComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/){
 172 // see header file for class documentation
 173
 174   if(GetFirstInputBlock(kAliHLTDataTypeSOR) || GetFirstInputBlock(kAliHLTDataTypeEOR)) return 0;
 175
 176   fNEvents++;
 177
 178   if(!fTracksArray) fTracksArray = new AliHLTTPCTrackArray();
 179
 180   const AliHLTComponentBlockData *iter = NULL;
 181
 182 //   //----------------- loop over slice tracks ----------------------//
 183 //
 184 //   for(iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkTrackSegmentsDataType); iter != NULL; iter = GetNextInputBlock()){
 185 //       if(iter->fDataType!=AliHLTTPCDefinitions::fgkTrackSegmentsDataType)continue;
 186 //       ReadTracks(iter,totalTracks);
 187 //   }
 188
 189
 190
 191   //----------------- loop over cluster blocks ---------------------//
 192
 193   Int_t totalSpacePoints = 0;
 194
 195   for(iter = GetFirstInputBlock(AliHLTTPCDefinitions::fgkClustersDataType); iter != NULL; iter = GetNextInputBlock()){
 196
 197       if(iter->fDataType!=AliHLTTPCDefinitions::fgkClustersDataType) continue;
 198
 199       AliHLTUInt8_t minSlice     = AliHLTTPCDefinitions::GetMinSliceNr(*iter);
 200       AliHLTUInt8_t minPartition = AliHLTTPCDefinitions::GetMinPatchNr(*iter);
 201       //HLTDebug("Input Data - TPC cluster - slice/partition: %d/%d.", minSlice, minPartition);
 202
 203       const AliHLTTPCClusterData* clusterData = (const AliHLTTPCClusterData*)iter->fPtr;
 204       Int_t nSpacepoint = (Int_t)clusterData->fSpacePointCnt;
 205       totalSpacePoints += nSpacepoint;
 206       HLTDebug("TrackHisto component found %d spacepoints in slice %d partition %d", nSpacepoint, minSlice, minPartition);
 207
 208       AliHLTTPCSpacePointData *clusters = (AliHLTTPCSpacePointData*)clusterData->fSpacePoints;
 209
 210       if(fClustersArray[minSlice][minPartition] != NULL){
 211          //delete(fClustersArray[minSlice][minPartition]);
 212          fClustersArray[minSlice][minPartition] = NULL;
 213       }
 214
 215       // fill the array with AliHLTTPCSpacePointData pointers
 216       // it will be used in the track loop to access information
 217       // for the used clusters only
 218       fClustersArray[minSlice][minPartition] = clusters;
 219       fNSpacePoints[minSlice][minPartition]  = nSpacepoint;
 220
 221       if(nSpacepoint==0) fClustersArray[minSlice][minPartition] = NULL;
 222
 223   } // end of loop over cluster data blocks
 224
 225   HLTInfo("TrackHisto found %d spacepoints",totalSpacePoints);
 226
 227
 228
 229
 230   //----------------- loop over merged tracks -------------------//
 231
 232   Int_t totalTracks = 0;
 233
 234   for(iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC); iter != NULL; iter = GetNextInputBlock()){
 235       if(iter->fDataType != (kAliHLTDataTypeTrack|kAliHLTDataOriginTPC)) continue;
 236       ReadTracks(iter,totalTracks);
 237   }
 238
 239   HLTInfo("TrackHisto found %d tracks", totalTracks);
 240
 241   fMultiplicity->Fill(totalTracks);
 242
 243   fNtotTracks += totalTracks;
 244   if(fNEvents%fNEvtMod==0){
 245      fMeanMultiplicity->Fill(fNEvents, Float_t(fNtotTracks)/(fNEvtMod));
 246      //HLTInfo("Event number: %d, total tracks accummulated %d", fNEvents, fNtotTracks);
 247      fNtotTracks = 0;
 248   }
 249
 250   PushHisto();
 251
 252   delete fTracksArray; fTracksArray = NULL;
 253   return 0;
 254 }
 255
 256 int AliHLTTPCTrackHistoComponent::Configure(const char* arguments){
 257 // see header file for class documentation
 258
 259    int iResult=0;
 260    if (!arguments) return iResult;
 261
 262    TString allArgs=arguments;
 263    TString argument;
 264
 265    TObjArray* pTokens=allArgs.Tokenize(" ");
 266    if (pTokens) {
 267      for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 268        argument=((TObjString*)pTokens->At(i))->GetString();
 269        if (argument.IsNull()) continue;
 270
 271        HLTError("unknown argument %s", argument.Data());
 272        iResult=-EINVAL;
 273        break;
 274      }
 275      delete pTokens;
 276    }
 277    return iResult;
 278  }
 279
 280 void AliHLTTPCTrackHistoComponent::ReadTracks(const AliHLTComponentBlockData* iter,Int_t &tt){
 281 // see header file for class documentation
 282
 283   AliHLTUInt8_t slice     = AliHLTTPCDefinitions::GetMinSliceNr(*iter);
 284   AliHLTUInt8_t partition = AliHLTTPCDefinitions::GetMinPatchNr(*iter);
 285
 286   if( slice     < fMinSlice )     fMinSlice     = slice;
 287   if( slice     > fMaxSlice )     fMaxSlice     = slice;
 288   if( partition < fMinPartition ) fMinPartition = partition;
 289   if( partition > fMaxPartition ) fMaxPartition = partition;
 290
 291   AliHLTTracksData *trackData = (AliHLTTracksData*)(iter->fPtr);
 292   AliHLTUInt32_t nTracks = trackData->fCount;
 293
 294   AliHLTExternalTrackParam *track = (AliHLTExternalTrackParam*)trackData->fTracklets;
 295   tt+= nTracks;
 296
 297   fTracksArray->FillTracksChecked(trackData->fTracklets,trackData->fCount,iter->fSize,slice,true);
 298
 299   Int_t usedSpacePoints = 0;
 300
 301   for(AliHLTUInt32_t i=0;i<nTracks;i++){
 302
 303     AliHLTTPCTrack * tpcTrack = 0;
 304     tpcTrack = (AliHLTTPCTrack*) fTracksArray->GetTrack(i);
 305     if(!tpcTrack) continue;
 306     Double_t trackLength = GetTrackLength(tpcTrack);
 307
 308     UInt_t nHits = track->fNPoints;
 309     fTracks->Fill( track->fq1Pt, track->fSinPsi, track->fTgl, nHits, GetEventId() );
 310
 311     const UInt_t *hitnum = track->fPointIDs;
 312
 313     Double_t totCharge = 0;
 314     for(UInt_t h=0; h<nHits; h++){
 315
 316       UInt_t idTrack = hitnum[h];
 317       Int_t sliceTrack = (idTrack>>25) & 0x7f;
 318       Int_t patchTrack = (idTrack>>22) & 0x7;
 319       UInt_t pos = idTrack&0x3fffff;
 320
 321       // use the fClustersArray that was filled in the cluster loop
 322       if( !fClustersArray[sliceTrack][patchTrack]  ) continue;
 323       if( fNSpacePoints[sliceTrack][patchTrack]<pos ) HLTError("Space point array out of boundaries!");
 324
 325       Float_t resy = 0., resz = 0.;
 326
 327       FillResidual(pos,sliceTrack,patchTrack,resy,resz);
 328
 329       totCharge += (fClustersArray[sliceTrack][patchTrack])[pos].fCharge;
 330
 331       fClusters->Fill( (fClustersArray[sliceTrack][patchTrack])[pos].fCharge, (fClustersArray[sliceTrack][patchTrack])[pos].fQMax, resy, resz, GetEventId() );
 332
 333       usedSpacePoints++;
 334     }
 335
 336     if(trackLength > 0) fDeDxVsP->Fill(track->fq1Pt*TMath::Sqrt(1.+track->fTgl*track->fTgl), totCharge/trackLength);
 337
 338     UChar_t *tmpP = (UChar_t*)track;
 339     tmpP += sizeof(AliHLTExternalTrackParam)+track->fNPoints*sizeof(UInt_t);
 340     track = (AliHLTExternalTrackParam*)tmpP;
 341   }
 342
 343
 344 /*  //HLTDebug ( "Input Data - TPC cluster - Slice/Patch: %d/%d.", slice, partition );
 345   AliHLTTPCTrackletData* trackData = (AliHLTTPCTrackletData*) iter->fPtr;
 346   AliHLTUInt32_t nTracks = trackData->fTrackletCnt;
 347   fTracksArray->FillTracksChecked(trackData->fTracklets,trackData->fTrackletCnt,iter->fSize,slice,true);
 348
 349   //HLTInfo("TrackHisto found %d Tracks in slice %d partition %d", nTracks, slice, partition);
 350   AliHLTTPCTrackSegmentData *tracks = (AliHLTTPCTrackSegmentData*) trackData->fTracklets;
 351
 352   for(AliHLTUInt32_t i=0;i<nTracks;i++){
 353     UInt_t nHits = tracks->fNPoints;
 354
 355     fTracks->Fill(tracks->fPt,tracks->fPsi,tracks->fTgl,nHits,GetEventId());
 356
 357     const UInt_t *hitnum = tracks->fPointIDs;
 358     for(UInt_t h=0; h<nHits; h++){
 359       UInt_t idTrack = hitnum[h];
 360       Int_t sliceTrack = (idTrack>>25) & 0x7f;
 361       Int_t patchTrack = (idTrack>>22) & 0x7;
 362       UInt_t pos = idTrack&0x3fffff;
 363       fTrackClusterID[sliceTrack][patchTrack].push_back(pos);
 364     }
 365     UChar_t *tmpP = (UChar_t*)tracks;
 366     tmpP += sizeof(AliHLTTPCTrackSegmentData)+tracks->fNPoints*sizeof(UInt_t);
 367     tracks = (AliHLTTPCTrackSegmentData*)tmpP;
 368   } */
 369
 370 }
 371
 372 void AliHLTTPCTrackHistoComponent::PushHisto(){
 373 // see header file for class documentation
 374
 375     AliHLTUInt32_t fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification(fMinSlice,fMaxSlice,fMinPartition,fMaxPartition);
 376
 377     PushBack( (TObject*)fTracks,           kAliHLTDataTypeTNtuple  |kAliHLTDataOriginTPC, fSpecification);
 378     PushBack( (TObject*)fClusters,         kAliHLTDataTypeTNtuple  |kAliHLTDataOriginTPC, fSpecification);
 379     PushBack( (TObject*)fMultiplicity,     kAliHLTDataTypeHistogram|kAliHLTDataOriginTPC, fSpecification);
 380     PushBack( (TObject*)fMeanMultiplicity, kAliHLTDataTypeHistogram|kAliHLTDataOriginTPC, fSpecification);
 381     PushBack( (TObject*)fDeDxVsP,          kAliHLTDataTypeHistogram|kAliHLTDataOriginTPC, fSpecification);
 382 }
 383
 384 void AliHLTTPCTrackHistoComponent::FillResidual(UInt_t pos,AliHLTUInt8_t slice,AliHLTUInt8_t partition,Float_t& resy,Float_t& resz){
 385 // see header file for class documentation
 386
 387   AliHLTTPCSpacePointData *cl =  &fClustersArray[slice][partition][pos];
 388   if(!cl) return;
 389
 390   AliHLTTPCTrack *gtrack = NULL;
 391
 392    for(int i=0;i<fTracksArray->GetNTracks();i++){
 393
 394        AliHLTTPCTrack *tt = fTracksArray->GetCheckedTrack(i);
 395        UInt_t *hitnum =tt->GetHitNumbers();
 396        Int_t nHits = tt->GetNHits();
 397
 398         for(Int_t h=0; h<nHits; h++){
 399             UInt_t id=hitnum[h];
 400              Int_t Tslice = (id>>25) & 0x7f;
 401              Int_t Tpatch = (id>>22) & 0x7;
 402              UInt_t Tpos = id&0x3fffff;
 403              if(Tslice==slice && Tpatch==partition && Tpos==pos){
 404                gtrack = tt;
 405                break;
 406              }
 407          }
 408    }
 409
 410   if(!gtrack) return;
 411
 412   Int_t tslice = gtrack->GetSector();
 413
 414   /*
 415   Double_t radius = gtrack->GetRadius();      // radius
 416   Double_t kappa = gtrack->GetKappa();        // curvature = 1/R , signed
 417   Double_t lambda = atan( gtrack->GetTgl() ); // dipAngle lambda
 418
 419   // ------------------------------------
 420   // ++ Get first/last point of the track
 421
 422   Double_t xyzL[3];      // lastpoint of track
 423   Double_t xyzF[3];      // firstpoint of track
 424
 425   xyzF[0] = gtrack->GetFirstPointX();
 426   xyzF[1] = gtrack->GetFirstPointY();
 427   xyzF[2] = gtrack->GetFirstPointZ();
 428
 429   xyzL[0] = gtrack->GetLastPointX();
 430   xyzL[1] = gtrack->GetLastPointY();
 431   xyzL[2] = gtrack->GetLastPointZ();
 432
 433   // --------------------------
 434   // ++ Calculate length of the track
 435
 436   Double_t s = 0.;       // length of the track
 437   if (  AliHLTTPCTransform::GetBFieldValue() == 0. || kappa == 0 )
 438     s = sqrt ( (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]) );
 439   else {
 440     // Calculate the length of the track. If it is to flat in in s,z plane use sxy, otherwise use sz
 441     if (fabs(lambda) > 0.05){
 442       // length of track calculated out of z
 443       s = fabs( (xyzL[2] - xyzF[2]) / sin(lambda) ); // length of track calculated out of z
 444     }
 445     else {
 446       Double_t d = (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0]) + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]);
 447       // length of track calculated out of xy
 448       s = fabs ( acos( 0.5 * (2 - (d / (radius*radius)))) / ( kappa * cos(lambda) ) );
 449     }
 450   }
 451   */  //???
 452
 453   gtrack->Rotate(tslice,kTRUE);
 454
 455   //Double_t padrows = 0;
 456
 457   Float_t xyzC[3];       // cluster tmp
 458   Float_t xyzTtmp[3];    // track tmp
 459
 460   xyzC[0] = cl->fX;
 461   xyzC[1] = cl->fY;
 462   xyzC[2] = cl->fZ;
 463
 464   Int_t padrow = AliHLTTPCTransform::GetPadRow(cl->fX);
 465
 466   xyzTtmp[0] = gtrack->GetFirstPointX();
 467
 468   if(gtrack->GetCrossingPoint(padrow,xyzTtmp)) {
 469     // ----------------------
 470        // ++ Calculate Residuals
 471
 472        Float_t deltaY = ( xyzC[1] - xyzTtmp[1] );
 473        Float_t deltaZ = ( xyzC[2] - xyzTtmp[2] );
 474
 475        resy = deltaY;
 476        resz = deltaZ;
 477   } else {
 478     resy = -1000;
 479     resz = -1000;
 480   }
 481 }
 482
 483 Double_t AliHLTTPCTrackHistoComponent::GetTrackLength(AliHLTTPCTrack *hltTrack)
 484 {
 485
 486   AliHLTTPCTrack * gtrack = hltTrack;
 487
 488   //Caculate the HLT Track Length
 489
 490   Double_t radius = gtrack->GetRadius();      // radius
 491   Double_t kappa = gtrack->GetKappa();        // curvature = 1/R , signed
 492   Double_t lambda = atan( gtrack->GetTgl() ); // dipAngle lambda
 493
 494   // ------------------------------------
 495   // ++ Get first/last point of the track
 496
 497   Double_t xyzL[3];      // lastpoint of track
 498   Double_t xyzF[3];      // firstpoint of track
 499
 500   xyzF[0] = gtrack->GetFirstPointX();
 501   xyzF[1] = gtrack->GetFirstPointY();
 502   xyzF[2] = gtrack->GetFirstPointZ();
 503
 504   xyzL[0] = gtrack->GetLastPointX();
 505   xyzL[1] = gtrack->GetLastPointY();
 506   xyzL[2] = gtrack->GetLastPointZ();
 507
 508   // --------------------------
 509   // ++ Calculate length of the track
 510
 511   Double_t s = 0.;       // length of the track
 512   if (  AliHLTTPCTransform::GetBFieldValue() == 0. || kappa == 0 )
 513     s = sqrt ( (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0])
 514                + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]) );
 515   else {
 516     // Calculate the length of the track. If it is to flat in in s,z plane use sxy, otherwise use sz
 517     if (fabs(lambda) > 0.05){
 518       // length of track calculated out of z
 519       s = fabs( (xyzL[2] - xyzF[2]) / sin(lambda) ); // length of track calculated out of z
 520   } else {
 521       Double_t d = (xyzL[0] - xyzF[0])*(xyzL[0] - xyzF[0])
 522         + (xyzL[1] - xyzF[1])*(xyzL[1] - xyzF[1]);
 523       // length of track calculated out of xy
 524       s = fabs ( acos( 0.5 * (2 - (d / (radius*radius))))
 525                  / ( kappa * cos(lambda) ) );
 526     }
 527   }
 528   return s;
 529 }