HLT/TRD/AliHLTTRDTrackerV1Component.cxx

   1 // $Id: AliHLTTRDTrackerV1Component.cxx 23618 2008-01-29 13:07:38Z hristov $
   2
   3 /**************************************************************************
   4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   5  *                                                                        *
   6  * Authors: Matthias Richter <Matthias.Richter@ift.uib.no>                *
   7  *          Timm Steinbeck <timm@kip.uni-heidelberg.de>                   *
   8  *          for The ALICE Off-line 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   AliHLTTRDTrackerV1Component.cxx
  20     @author Timm Steinbeck, Matthias Richter
  21     @date
  22     @brief  A TRDTrackerV1 processing component for the HLT. */
  23
  24 #if __GNUC__ >= 3
  25 using namespace std;
  26 #endif
  27
  28 #include "AliHLTTRDTrackerV1Component.h"
  29 #include "AliHLTTRDDefinitions.h"
  30 #include "AliHLTTRDCluster.h"
  31 #include "AliHLTTRDTrack.h"
  32 #include "AliHLTTRDUtils.h"
  33
  34 #include "TFile.h"
  35 #include "TChain.h"
  36
  37 #include "AliGeomManager.h"
  38 #include "AliCDBManager.h"
  39 #include "AliESDEvent.h"
  40 #include "AliMagF.h"
  41 #include "AliESDfriend.h"
  42
  43 #include "AliTRDcalibDB.h"
  44 #include "AliTRDReconstructor.h"
  45 #include "AliTRDtrackerV1.h"
  46 #include "AliTRDrecoParam.h"
  47
  48 #include <cstdlib>
  49 #include <cerrno>
  50 #include <string>
  51
  52 #ifdef HAVE_VALGRIND_CALLGRIND_H
  53 #include <valgrind/callgrind.h>
  54 #else
  55 #define CALLGRIND_START_INSTRUMENTATION do { } while (0)
  56 #define CALLGRIND_STOP_INSTRUMENTATION do { } while (0)
  57 #endif
  58
  59 ClassImp(AliHLTTRDTrackerV1Component);
  60
  61 AliHLTTRDTrackerV1Component::AliHLTTRDTrackerV1Component():
  62   AliHLTProcessor(),
  63   fOutputPercentage(100), // By default we copy to the output exactly what we got as input
  64   fStrorageDBpath("local://$ALICE_ROOT/OCDB"),
  65   fCDB(NULL),
  66   fGeometryFileName(""),
  67   fTracker(NULL),
  68   fRecoParam(NULL),
  69   fReconstructor(NULL)
  70 {
  71   // Default constructor
  72
  73   fGeometryFileName = getenv("ALICE_ROOT");
  74   fGeometryFileName += "/HLT/TRD/geometry.root";
  75 }
  76
  77 AliHLTTRDTrackerV1Component::~AliHLTTRDTrackerV1Component()
  78 {
  79   // Destructor
  80 }
  81
  82 const char* AliHLTTRDTrackerV1Component::GetComponentID()
  83 {
  84   // Return the component ID const char *
  85   return "TRDTrackerV1"; // The ID of this component
  86 }
  87
  88 void AliHLTTRDTrackerV1Component::GetInputDataTypes( vector<AliHLTComponent_DataType>& list)
  89 {
  90   // Get the list of input data
  91   list.clear(); // We do not have any requirements for our input data type(s).
  92   list.push_back( AliHLTTRDDefinitions::fgkClusterDataType );
  93 }
  94
  95 AliHLTComponent_DataType AliHLTTRDTrackerV1Component::GetOutputDataType()
  96 {
  97   // Get the output data type
  98   return AliHLTTRDDefinitions::fgkClusterDataType;
  99 }
 100
 101 void AliHLTTRDTrackerV1Component::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 102 {
 103   // Get the output data size
 104   constBase = 0;
 105   inputMultiplier = ((double)fOutputPercentage)/100.0;
 106 }
 107
 108 // Spawn function, return new instance of this class
 109 AliHLTComponent* AliHLTTRDTrackerV1Component::Spawn()
 110 {
 111   return new AliHLTTRDTrackerV1Component;
 112 };
 113
 114
 115 int AliHLTTRDTrackerV1Component::DoInit( int argc, const char** argv )
 116 {
 117   // perform initialization. We check whether our relative output size is specified in the arguments.
 118   fOutputPercentage = 100;
 119   int i = 0;
 120   char* cpErr;
 121
 122
 123   Int_t iRecoParamType = -1; // default will be the low flux
 124   Int_t iNtimeBins = -1;     // number of time bins for the tracker to use
 125   Int_t iMagneticField = -1; // magnetic field: 0==OFF and 1==ON
 126   Bool_t bHLTMode = kTRUE, bWriteClusters = kFALSE;
 127
 128   while ( i < argc )
 129     {
 130       HLTDebug("argv[%d] == %s", i, argv[i] );
 131       if ( !strcmp( argv[i], "output_percentage" ) )
 132         {
 133           if ( i+1>=argc )
 134             {
 135               HLTError("Missing output_percentage parameter");
 136               return ENOTSUP;
 137             }
 138           HLTDebug("argv[%d+1] == %s", i, argv[i+1] );
 139           fOutputPercentage = strtoul( argv[i+1], &cpErr, 0 );
 140           if ( *cpErr )
 141             {
 142               HLTError("Cannot convert output_percentage parameter '%s'", argv[i+1] );
 143               return EINVAL;
 144             }
 145           HLTInfo("Output percentage set to %lu %%", fOutputPercentage );
 146           i += 2;
 147
 148         }
 149
 150       else if ( !strcmp( argv[i], "-NTimeBins" ) )
 151         {
 152           if ( i+1>=argc )
 153             {
 154               HLTError("Missing -NTimeBins parameter");
 155               return ENOTSUP;
 156             }
 157           HLTDebug("Arguments", "argv[%d+1] == %s", i, argv[i+1] );
 158           iNtimeBins = strtoul( argv[i+1], &cpErr, 0 );
 159           if ( *cpErr )
 160             {
 161               HLTError("Wrong Argument. Cannot convert -NTimeBins parameter '%s'", argv[i+1] );
 162               return EINVAL;
 163             }
 164           i += 2;
 165
 166         }
 167
 168       else if ( strcmp( argv[i], "-cdb" ) == 0)
 169         {
 170           if ( i+1 >= argc )
 171             {
 172               HLTError( "Missing -cdb argument");
 173               return ENOTSUP;
 174             }
 175           fStrorageDBpath = argv[i+1];
 176           HLTInfo("DB storage is %s", fStrorageDBpath.c_str() );
 177           i += 2;
 178
 179         }
 180
 181       else if ( strcmp( argv[i], "-geometry" ) == 0)
 182         {
 183           if ( i+1 >= argc )
 184             {
 185               HLTError("Missing -geometry argument");
 186               return ENOTSUP;
 187             }
 188           fGeometryFileName = argv[i+1];
 189           HLTInfo("GeomFile storage is %s",
 190                   fGeometryFileName.c_str() );
 191           i += 2;
 192
 193         }
 194
 195       // the flux parametrizations
 196       else if ( strcmp( argv[i], "-lowflux" ) == 0)
 197         {
 198           iRecoParamType = 0;
 199           HLTDebug("Low flux reco selected.");
 200           i++;
 201
 202         }
 203
 204       else if ( strcmp( argv[i], "-highflux" ) == 0)
 205         {
 206           iRecoParamType = 1;
 207           HLTDebug("Low flux reco selected.");
 208           i++;
 209         }
 210
 211       else if ( strcmp( argv[i], "-cosmics" ) == 0)
 212         {
 213           iRecoParamType = 2;
 214           HLTDebug("Cosmic test reco selected.");
 215           i++;
 216         }
 217
 218       else if ( strcmp( argv[i], "-magnetic_field_ON" ) == 0)
 219         {
 220           iMagneticField = 1;
 221           i++;
 222         }
 223       else if ( strcmp( argv[i], "-magnetic_field_OFF" ) == 0)
 224         {
 225           iMagneticField = 0;
 226           i++;
 227         }
 228       else if ( strcmp( argv[i], "-offlineMode" ) == 0)
 229         {
 230           bHLTMode=kFALSE;
 231           HLTDebug("Using standard offline tracking.");
 232           i++;
 233         }
 234       else {
 235         HLTError("Unknown option '%s'", argv[i] );
 236         return EINVAL;
 237       }
 238
 239     }
 240
 241   // THE "REAL" INIT COMES HERE
 242   // offline condition data base
 243   fCDB = AliCDBManager::Instance();
 244   if (!fCDB)
 245     {
 246       HLTError("Could not get CDB instance", "fCDB 0x%x", fCDB);
 247       return -1;
 248     }
 249   else
 250     {
 251       fCDB->SetRun(0); // THIS HAS TO BE RETRIEVED !!!
 252       fCDB->SetDefaultStorage(fStrorageDBpath.c_str());
 253       HLTDebug("CDB instance", "fCDB 0x%x", fCDB);
 254     }
 255
 256   // check if the N of time bins make sense
 257   if (iNtimeBins <= 0)
 258     {
 259       HLTError("Sorry. Tracker needs number of time bins. At the moment you have to provide it with -NTimeBins <value>. The simulation always had 24 and the real data 30. Take your pick. Make sure the information is correct. Ask offline to implement how to propagate this information into clusters/cluster tree.");
 260       return -1;
 261     }
 262
 263   if (iNtimeBins < 24 || iNtimeBins > 30)
 264     {
 265       HLTWarning("The number of time bins seems to be strange = %d. But okay. Let's try it...", iNtimeBins);
 266     }
 267
 268   HLTDebug("The number of time bins = %d.", iNtimeBins);
 269   AliTRDtrackerV1::SetNTimeBins(iNtimeBins);
 270
 271   // !!!! THIS IS IMPORTANT
 272   // init alifield map - temporarly via parameter - should come from a DB or DCS ?
 273   // !!!!
 274   if (iMagneticField < 0)
 275     {
 276       iMagneticField = 0;
 277       HLTWarning("No magnetic field switch stated. Use -magnetic_field_ON or -magnetic_field_OFF flag. Defaulting to OFF = NO MAGNETIC FIELD");
 278     }
 279
 280   if (!TGeoGlobalMagField::Instance()->IsLocked()) {
 281     if (iMagneticField == 0)
 282       {
 283         // magnetic field OFF
 284         AliMagF* field = new AliMagF("Maps","Maps",2,0.,0., 10.,AliMagF::k5kGUniform);
 285         TGeoGlobalMagField::Instance()->SetField(field);
 286         HLTDebug("Magnetic field is OFF.");
 287       }
 288
 289     if (iMagneticField == 1)
 290       {
 291         // magnetic field ON
 292         AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
 293         TGeoGlobalMagField::Instance()->SetField(field);
 294         HLTDebug("Magnetic field is ON.");
 295       }
 296   }
 297   else {
 298     HLTError("Magnetic field is already set and locked, cannot redefine it." );
 299   }
 300
 301   // reconstruction parameters
 302   if (iRecoParamType < 0 || iRecoParamType > 2)
 303     {
 304       HLTWarning("No reco param selected. Use -lowflux -highflux -cosmics flags. Defaulting to low flux.");
 305       iRecoParamType = 0;
 306     }
 307
 308   if (iRecoParamType == 0)
 309     {
 310       fRecoParam = AliTRDrecoParam::GetLowFluxParam();
 311       HLTDebug("Low flux params init.");
 312     }
 313
 314   if (iRecoParamType == 1)
 315     {
 316       fRecoParam = AliTRDrecoParam::GetHighFluxParam();
 317       HLTDebug("High flux params init.");
 318     }
 319
 320   if (iRecoParamType == 2)
 321     {
 322       fRecoParam = AliTRDrecoParam::GetCosmicTestParam();
 323       HLTDebug("Cosmic Test params init.");
 324     }
 325
 326   if (fRecoParam == 0)
 327     {
 328       HLTError("No reco params initialized. Sniffing big trouble!");
 329       return -1;
 330     }
 331
 332   fReconstructor = new AliTRDReconstructor();
 333   //   fRecoParam->SetChi2Y(.1);
 334   //   fRecoParam->SetChi2Z(5.);
 335   fReconstructor->SetRecoParam(fRecoParam);
 336   // write clusters [cw] = true
 337   // track seeding (stand alone tracking) [sa] = true
 338   // PID method in reconstruction (NN) [nn] = true
 339   // write online tracklets [tw] = false
 340   // drift gas [ar] = false
 341   // sl_tr_0 = StreamLevel_task_Level
 342   //  fReconstructor->SetOption("sa,!cw,hlt,sl_tr_0");
 343   TString recoOptions="sa,sl_tr_0";
 344
 345   if (bWriteClusters)
 346     {
 347       recoOptions += ",cw";
 348     }
 349   else
 350     {
 351       recoOptions += ",!cw";
 352     }
 353   if (bHLTMode)
 354     recoOptions += ",hlt";
 355
 356   fReconstructor->SetOption(recoOptions.Data());
 357   HLTDebug("Reconstructor options are: %s",recoOptions.Data());
 358
 359   if((AliGeomManager::GetGeometry()) == NULL){
 360
 361     if ( TFile::Open(fGeometryFileName.c_str())) {
 362       AliGeomManager::LoadGeometry(fGeometryFileName.c_str());
 363     }
 364     else {
 365       HLTError("Cannot load geometry from file %s",fGeometryFileName.c_str());
 366       return EINVAL;
 367     }
 368   }
 369   else
 370     HLTInfo("Geometry Already Loaded");
 371
 372   // create the tracker
 373   fTracker = new AliTRDtrackerV1();
 374   fTracker->SetReconstructor(fReconstructor);
 375   HLTDebug("TRDTracker at 0x%x", fTracker);
 376
 377   if (fTracker == 0)
 378     {
 379       HLTError("Unable to create the tracker!");
 380       return -1;
 381     }
 382
 383   return 0;
 384 }
 385
 386 int AliHLTTRDTrackerV1Component::DoDeinit()
 387 {
 388   // Deinitialization of the component
 389
 390   fTracker->SetClustersOwner(kFALSE);
 391   delete fTracker;
 392   fTracker = 0x0;
 393
 394   // We need to set clusters in Reconstructor to null to prevent from
 395   // double deleting, since we delete TClonesArray by ourself in DoEvent.
 396   fReconstructor->SetClusters(0x0);
 397   delete fReconstructor;
 398   fReconstructor = 0x0;
 399
 400   AliTRDcalibDB::Terminate();
 401
 402   return 0;
 403 }
 404
 405 int AliHLTTRDTrackerV1Component::DoEvent( const AliHLTComponentEventData& evtData,
 406                                           const AliHLTComponentBlockData* blocks,
 407                                           AliHLTComponent_TriggerData& /*trigData*/,
 408                                           AliHLTUInt8_t* outputPtr,
 409                                           AliHLTUInt32_t& size,
 410                                           vector<AliHLTComponent_BlockData>& outputBlocks )
 411 {
 412   // Process an event
 413
 414   HLTDebug("NofBlocks %lu", evtData.fBlockCnt );
 415
 416   AliHLTUInt32_t totalSize = 0, offset = 0;
 417   AliHLTUInt32_t dBlockSpecification = 0;
 418
 419   vector<AliHLTComponent_DataType> expectedDataTypes;
 420   GetInputDataTypes(expectedDataTypes);
 421   if (evtData.fEventID == 1)
 422     CALLGRIND_START_INSTRUMENTATION;
 423   for ( unsigned long iBlock = 0; iBlock < evtData.fBlockCnt; iBlock++ )
 424     {
 425       const AliHLTComponentBlockData &block = blocks[iBlock];
 426       offset = totalSize;
 427       AliHLTComponentDataType inputDataType = block.fDataType;
 428       Bool_t correctDataType = kFALSE;
 429
 430       for(UInt_t i = 0; i < expectedDataTypes.size(); i++)
 431         if( expectedDataTypes.at(i) == inputDataType)
 432           correctDataType = kTRUE;
 433       if (!correctDataType)
 434         {
 435           HLTDebug( "Block # %i/%i; Event 0x%08LX (%Lu) Wrong received datatype: %s - Skipping",
 436                     iBlock, evtData.fBlockCnt-1,
 437                     evtData.fEventID, evtData.fEventID,
 438                     DataType2Text(inputDataType).c_str());
 439           continue;
 440         }
 441       else {
 442           HLTDebug("We get the right data type: Block # %i/%i; Event 0x%08LX (%Lu) Received datatype: %s",
 443                    iBlock, evtData.fBlockCnt-1,
 444                    evtData.fEventID, evtData.fEventID,
 445                    DataType2Text(inputDataType).c_str());
 446       }
 447
 448
 449       TClonesArray* clusterArray = new TClonesArray("AliTRDcluster"); // would be nice to allocate memory for all clusters here.
 450       AliHLTTRDUtils::ReadClusters(clusterArray, block.fPtr, block.fSize);
 451       HLTDebug("TClonesArray of clusters: nbEntries = %i", clusterArray->GetEntriesFast());
 452       fTracker->LoadClusters(clusterArray);
 453
 454
 455
 456       // maybe it is not so smart to create it each event? clear is enough ?
 457       AliESDEvent *esd = new AliESDEvent();
 458       esd->CreateStdContent();
 459       fTracker->Clusters2Tracks(esd);
 460
 461       Int_t nTracks = esd->GetNumberOfTracks();
 462       HLTInfo("Number of tracks  == %d ==", nTracks);
 463
 464       TClonesArray* trdTracks = fTracker->GetListOfTracks();
 465
 466       if (trdTracks)
 467         totalSize += TransportTracks(trdTracks, outputPtr, outputBlocks, offset, dBlockSpecification);
 468       else {
 469           HLTDebug("Bad array trdTracks = 0x%x", trdTracks);
 470       }
 471
 472       HLTDebug("totalSize: %i", totalSize);
 473
 474 //       if ( totalSize > allocSize )
 475 //      {
 476 //        HLTError("Too much data; Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
 477 //        totalSize, size );
 478 //        return EMSGSIZE;
 479 //      }
 480
 481       //here we are deleting clusters (but not the TClonesArray itself)
 482       fTracker->UnloadClusters();
 483
 484       AliTRDReconstructor::SetClusters(0x0);
 485       delete esd;
 486       clusterArray->Delete();
 487       delete clusterArray;
 488
 489       }
 490
 491   size = totalSize;
 492   HLTDebug("Event is done. size written to the output is %i", size);
 493   return 0;
 494 }
 495
 496
 497 /**
 498  * Transport tracks to the next component
 499  * Return Numbers of bytes written to the output
 500  */
 501 //============================================================================
 502 AliHLTUInt32_t AliHLTTRDTrackerV1Component::TransportTracks(TClonesArray *inTracksArray, AliHLTUInt8_t* output,
 503                                                     vector<AliHLTComponent_BlockData>& outputBlocks, AliHLTUInt32_t inOffset, AliHLTUInt32_t inSpec)
 504 {
 505   Int_t nbTracks=inTracksArray->GetEntriesFast();
 506   if (nbTracks>0)
 507     {
 508       HLTDebug("We have an output array: pointer to inTracksArray = 0x%x, nbEntries = %i", inTracksArray, nbTracks);
 509
 510         // Using low-level interface
 511         // with interface classes
 512         AliHLTUInt32_t addedSize = AliHLTTRDUtils::AddTracksToOutput(inTracksArray, output);
 513
 514         // Fill block
 515         AliHLTComponentBlockData bd;
 516         FillBlockData( bd );
 517
 518         bd.fPtr = output;
 519         bd.fOffset = inOffset;
 520         bd.fSize = addedSize;
 521         bd.fSpecification = inSpec;
 522         bd.fDataType = AliHLTTRDDefinitions::fgkTRDSATracksDataType;
 523         outputBlocks.push_back( bd );
 524         HLTDebug("BD fPtr 0x%x, fOffset %i, fSize %i, fSpec 0x%x", bd.fPtr, bd.fOffset, bd.fSize, bd.fSpecification);
 525
 526         return addedSize;
 527
 528       }
 529   return 0;
 530
 531 }