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   return  kAliHLTDataTypeTrack | kAliHLTDataOriginTRD;
 100 }
 101
 102 void AliHLTTRDTrackerV1Component::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 103 {
 104   // Get the output data size
 105   constBase = 0;
 106   inputMultiplier = ((double)fOutputPercentage)/100.0;
 107 }
 108
 109 // Spawn function, return new instance of this class
 110 AliHLTComponent* AliHLTTRDTrackerV1Component::Spawn()
 111 {
 112   return new AliHLTTRDTrackerV1Component;
 113 };
 114
 115
 116 int AliHLTTRDTrackerV1Component::DoInit( int argc, const char** argv )
 117 {
 118   // perform initialization. We check whether our relative output size is specified in the arguments.
 119   fOutputPercentage = 100;
 120   int i = 0;
 121   char* cpErr;
 122
 123
 124   Int_t iRecoParamType = -1; // default will be the low flux
 125   Int_t iNtimeBins = -1;     // number of time bins for the tracker to use
 126   Int_t iMagneticField = -1; // magnetic field: 0==OFF and 1==ON
 127   Int_t iPIDmethod = 1;      // 0=LikelyHood(LH) 1=NeuronalNetwork(NN) 2=TruncatedMean(TM)
 128   Bool_t bHLTMode = kTRUE, bWriteClusters = kFALSE;
 129
 130   while ( i < argc )
 131     {
 132       HLTDebug("argv[%d] == %s", i, argv[i] );
 133       if ( !strcmp( argv[i], "output_percentage" ) )
 134         {
 135           if ( i+1>=argc )
 136             {
 137               HLTError("Missing output_percentage parameter");
 138               return ENOTSUP;
 139             }
 140           HLTDebug("argv[%d+1] == %s", i, argv[i+1] );
 141           fOutputPercentage = strtoul( argv[i+1], &cpErr, 0 );
 142           if ( *cpErr )
 143             {
 144               HLTError("Cannot convert output_percentage parameter '%s'", argv[i+1] );
 145               return EINVAL;
 146             }
 147           HLTInfo("Output percentage set to %lu %%", fOutputPercentage );
 148           i += 2;
 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       else if ( strcmp( argv[i], "-cdb" ) == 0)
 167         {
 168           if ( i+1 >= argc )
 169             {
 170               HLTError( "Missing -cdb argument");
 171               return ENOTSUP;
 172             }
 173           fStrorageDBpath = argv[i+1];
 174           HLTInfo("DB storage is %s", fStrorageDBpath.c_str() );
 175           i += 2;
 176         }
 177       else if ( strcmp( argv[i], "-geometry" ) == 0)
 178         {
 179           if ( i+1 >= argc )
 180             {
 181               HLTError("Missing -geometry argument");
 182               return ENOTSUP;
 183             }
 184           fGeometryFileName = argv[i+1];
 185           HLTInfo("GeomFile storage is %s",
 186                   fGeometryFileName.c_str() );
 187           i += 2;
 188         }
 189       // the flux parametrizations
 190       else if ( strcmp( argv[i], "-lowflux" ) == 0)
 191         {
 192           iRecoParamType = 0;
 193           HLTDebug("Low flux reco selected.");
 194           i++;
 195         }
 196       else if ( strcmp( argv[i], "-highflux" ) == 0)
 197         {
 198           iRecoParamType = 1;
 199           HLTDebug("Low flux reco selected.");
 200           i++;
 201         }
 202       else if ( strcmp( argv[i], "-cosmics" ) == 0)
 203         {
 204           iRecoParamType = 2;
 205           HLTDebug("Cosmic test reco selected.");
 206           i++;
 207         }
 208       else if ( strcmp( argv[i], "-magnetic_field_ON" ) == 0)
 209         {
 210           iMagneticField = 1;
 211           i++;
 212         }
 213       else if ( strcmp( argv[i], "-magnetic_field_OFF" ) == 0)
 214         {
 215           iMagneticField = 0;
 216           i++;
 217         }
 218       else if ( strcmp( argv[i], "-offlineMode" ) == 0)
 219         {
 220           bHLTMode=kFALSE;
 221           HLTDebug("Using standard offline tracking.");
 222           i++;
 223         }
 224       else if ( strcmp( argv[i], "-PIDmethod" ) == 0)
 225         {
 226           if ( i+1 >= argc )
 227             {
 228               HLTError("Missing -PIDmethod argument");
 229               return ENOTSUP;
 230             }
 231           if( strcmp(argv[i], "LH") )
 232             iPIDmethod=0;
 233           else if( strcmp(argv[i], "NN") )
 234             iPIDmethod=1;
 235           else if( strcmp(argv[i], "TM") )
 236             iPIDmethod=2;
 237           else {
 238             HLTError("Unknown -PIDmethod argument");
 239             return ENOTSUP;
 240           }
 241           i += 2;
 242         }
 243
 244       else {
 245         HLTError("Unknown option '%s'", argv[i] );
 246         return EINVAL;
 247       }
 248
 249     }
 250
 251   // THE "REAL" INIT COMES HERE
 252   // offline condition data base
 253   fCDB = AliCDBManager::Instance();
 254   if (!fCDB)
 255     {
 256       HLTError("Could not get CDB instance", "fCDB 0x%x", fCDB);
 257       return -1;
 258     }
 259   else
 260     {
 261       fCDB->SetRun(0); // THIS HAS TO BE RETRIEVED !!!
 262       fCDB->SetDefaultStorage(fStrorageDBpath.c_str());
 263       HLTDebug("CDB instance", "fCDB 0x%x", fCDB);
 264     }
 265
 266   // check if the N of time bins make sense
 267   if (iNtimeBins <= 0)
 268     {
 269       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.");
 270       return -1;
 271     }
 272
 273   if (iNtimeBins < 24 || iNtimeBins > 30)
 274     {
 275       HLTWarning("The number of time bins seems to be strange = %d. But okay. Let's try it...", iNtimeBins);
 276     }
 277
 278   HLTDebug("The number of time bins = %d.", iNtimeBins);
 279   AliTRDtrackerV1::SetNTimeBins(iNtimeBins);
 280
 281   // !!!! THIS IS IMPORTANT
 282   // init alifield map - temporarly via parameter - should come from a DB or DCS ?
 283   // !!!!
 284   if (iMagneticField < 0)
 285     {
 286       iMagneticField = 0;
 287       HLTWarning("No magnetic field switch stated. Use -magnetic_field_ON or -magnetic_field_OFF flag. Defaulting to OFF = NO MAGNETIC FIELD");
 288     }
 289
 290   if (!TGeoGlobalMagField::Instance()->IsLocked()) {
 291     if (iMagneticField == 0)
 292       {
 293         // magnetic field OFF
 294         AliMagF* field = new AliMagF("Maps","Maps",2,0.,0., 10.,AliMagF::k5kGUniform);
 295         TGeoGlobalMagField::Instance()->SetField(field);
 296         HLTDebug("Magnetic field is OFF.");
 297       }
 298
 299     if (iMagneticField == 1)
 300       {
 301         // magnetic field ON
 302         AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
 303         TGeoGlobalMagField::Instance()->SetField(field);
 304         HLTDebug("Magnetic field is ON.");
 305       }
 306   }
 307   else {
 308     HLTError("Magnetic field is already set and locked, cannot redefine it." );
 309   }
 310
 311   // reconstruction parameters
 312   if (iRecoParamType < 0 || iRecoParamType > 2)
 313     {
 314       HLTWarning("No reco param selected. Use -lowflux -highflux -cosmics flags. Defaulting to low flux.");
 315       iRecoParamType = 0;
 316     }
 317
 318   if (iRecoParamType == 0)
 319     {
 320       fRecoParam = AliTRDrecoParam::GetLowFluxParam();
 321       HLTDebug("Low flux params init.");
 322     }
 323
 324   if (iRecoParamType == 1)
 325     {
 326       fRecoParam = AliTRDrecoParam::GetHighFluxParam();
 327       HLTDebug("High flux params init.");
 328     }
 329
 330   if (iRecoParamType == 2)
 331     {
 332       fRecoParam = AliTRDrecoParam::GetCosmicTestParam();
 333       HLTDebug("Cosmic Test params init.");
 334     }
 335
 336   if (fRecoParam == 0)
 337     {
 338       HLTError("No reco params initialized. Sniffing big trouble!");
 339       return -1;
 340     }
 341
 342   fReconstructor = new AliTRDReconstructor();
 343   //   fRecoParam->SetChi2Y(.1);
 344   //   fRecoParam->SetChi2Z(5.);
 345   fReconstructor->SetRecoParam(fRecoParam);
 346   // write clusters [cw] = true
 347   // track seeding (stand alone tracking) [sa] = true
 348   // PID method in reconstruction (NN) [nn] = true
 349   // write online tracklets [tw] = false
 350   // drift gas [ar] = false
 351   // sl_tr_0 = StreamLevel_task_Level
 352   //  fReconstructor->SetOption("sa,!cw,hlt,sl_tr_0");
 353   TString recoOptions="sa,sl_tr_0";
 354
 355   if (bWriteClusters)
 356     {
 357       recoOptions += ",cw";
 358     }
 359   else
 360     {
 361       recoOptions += ",!cw";
 362     }
 363   if (bHLTMode)
 364     recoOptions += ",hlt";
 365
 366   switch(iPIDmethod){
 367   case 0: recoOptions += ",!nn"; break;
 368   case 1: recoOptions += ",nn"; break;
 369   case 2: recoOptions += ",!nn"; break;
 370   }
 371
 372   fReconstructor->SetOption(recoOptions.Data());
 373   HLTDebug("Reconstructor options are: %s",recoOptions.Data());
 374
 375   if((AliGeomManager::GetGeometry()) == NULL){
 376
 377     if ( TFile::Open(fGeometryFileName.c_str())) {
 378       AliGeomManager::LoadGeometry(fGeometryFileName.c_str());
 379     }
 380     else {
 381       HLTError("Cannot load geometry from file %s",fGeometryFileName.c_str());
 382       return EINVAL;
 383     }
 384   }
 385   else
 386     HLTInfo("Geometry Already Loaded");
 387
 388   // create the tracker
 389   fTracker = new AliTRDtrackerV1();
 390   fTracker->SetReconstructor(fReconstructor);
 391   HLTDebug("TRDTracker at 0x%x", fTracker);
 392
 393   if (fTracker == 0)
 394     {
 395       HLTError("Unable to create the tracker!");
 396       return -1;
 397     }
 398
 399   return 0;
 400 }
 401
 402 int AliHLTTRDTrackerV1Component::DoDeinit()
 403 {
 404   // Deinitialization of the component
 405
 406   fTracker->SetClustersOwner(kFALSE);
 407   delete fTracker;
 408   fTracker = 0x0;
 409
 410   // We need to set clusters in Reconstructor to null to prevent from
 411   // double deleting, since we delete TClonesArray by ourself in DoEvent.
 412   fReconstructor->SetClusters(0x0);
 413   delete fReconstructor;
 414   fReconstructor = 0x0;
 415
 416   AliTRDcalibDB::Terminate();
 417
 418   return 0;
 419 }
 420
 421 int AliHLTTRDTrackerV1Component::DoEvent( const AliHLTComponentEventData& evtData,
 422                                           const AliHLTComponentBlockData* blocks,
 423                                           AliHLTComponent_TriggerData& /*trigData*/,
 424                                           AliHLTUInt8_t* outputPtr,
 425                                           AliHLTUInt32_t& size,
 426                                           vector<AliHLTComponent_BlockData>& outputBlocks )
 427 {
 428   // Process an event
 429
 430   HLTDebug("NofBlocks %lu", evtData.fBlockCnt );
 431
 432   AliHLTUInt32_t totalSize = 0, offset = 0;
 433   AliHLTUInt32_t dBlockSpecification = 0;
 434
 435   vector<AliHLTComponent_DataType> expectedDataTypes;
 436   GetInputDataTypes(expectedDataTypes);
 437   if (evtData.fEventID > 1)
 438     CALLGRIND_START_INSTRUMENTATION;
 439   for ( unsigned long iBlock = 0; iBlock < evtData.fBlockCnt; iBlock++ )
 440     {
 441       const AliHLTComponentBlockData &block = blocks[iBlock];
 442       AliHLTComponentDataType inputDataType = block.fDataType;
 443       Bool_t correctDataType = kFALSE;
 444
 445       for(UInt_t i = 0; i < expectedDataTypes.size(); i++){
 446         if( expectedDataTypes.at(i) == inputDataType)
 447           correctDataType = kTRUE;
 448       }
 449       if (!correctDataType)
 450         {
 451           HLTDebug( "Block # %i/%i; Event 0x%08LX (%Lu) Wrong received datatype: %s - Skipping",
 452                     iBlock, evtData.fBlockCnt-1,
 453                     evtData.fEventID, evtData.fEventID,
 454                     DataType2Text(inputDataType).c_str());
 455           continue;
 456         }
 457       else {
 458         HLTDebug("We get the right data type: Block # %i/%i; Event 0x%08LX (%Lu) Received datatype: %s; Block Size: %i",
 459                  iBlock, evtData.fBlockCnt-1,
 460                  evtData.fEventID, evtData.fEventID,
 461                  DataType2Text(inputDataType).c_str(),
 462                  block.fSize);
 463       }
 464
 465
 466       TClonesArray* clusterArray = new TClonesArray("AliTRDcluster"); // would be nice to allocate memory for all clusters here.
 467       AliHLTTRDUtils::ReadClusters(clusterArray, block.fPtr, block.fSize);
 468       HLTDebug("TClonesArray of clusters: nbEntries = %i", clusterArray->GetEntriesFast());
 469       fTracker->LoadClusters(clusterArray);
 470
 471       // maybe it is not so smart to create it each event? clear is enough ?
 472       AliESDEvent *esd = new AliESDEvent();
 473       esd->CreateStdContent();
 474       fTracker->Clusters2Tracks(esd);
 475
 476       Int_t nTracks = esd->GetNumberOfTracks();
 477       HLTInfo("Number of tracks  == %d ==", nTracks);
 478
 479       TClonesArray* trdTracks = 0x0;
 480       //trdTracks = fTracker->GetListOfTracks();
 481
 482       if(nTracks>0){
 483         HLTDebug("We have an output ESDEvent: 0x%x with %i tracks", esd, nTracks);
 484         AliHLTUInt32_t addedSize = AliHLTTRDUtils::AddESDToOutput(esd, (AliHLTUInt8_t*)(outputPtr+offset));
 485         totalSize += addedSize;
 486
 487         // Fill block
 488         AliHLTComponentBlockData bd;
 489         FillBlockData( bd );
 490         //bd.fPtr = outputPtr;
 491         bd.fOffset = offset;
 492         bd.fSize = addedSize;
 493         bd.fSpecification = dBlockSpecification;
 494         bd.fDataType = kAliHLTDataTypeTrack | kAliHLTDataOriginTRD;
 495         outputBlocks.push_back( bd );
 496         HLTDebug("BD fPtr 0x%x, fOffset %i, fSize %i, fSpec 0x%x", bd.fPtr, bd.fOffset, bd.fSize, bd.fSpecification);
 497         offset = totalSize;
 498
 499         if (trdTracks){
 500           //Int_t nbTracks=trdTracks->GetEntriesFast();
 501           //if (nbTracks>0){
 502           HLTDebug("We have an output array: pointer to trdTracks = 0x%x, nbEntries = %i", trdTracks, trdTracks->GetEntriesFast());
 503
 504           AliHLTUInt32_t addedSize = AliHLTTRDUtils::AddTracksToOutput(trdTracks, (AliHLTUInt8_t*)(outputPtr+offset));
 505           totalSize += addedSize;
 506
 507           // Fill block
 508           AliHLTComponentBlockData bd;
 509           FillBlockData( bd );
 510           //bd.fPtr = outputPtr;
 511           bd.fOffset = offset;
 512           bd.fSize = addedSize;
 513           bd.fSpecification = dBlockSpecification;
 514           bd.fDataType = AliHLTTRDDefinitions::fgkTRDSATracksDataType;
 515           outputBlocks.push_back( bd );
 516           HLTDebug("BD fPtr 0x%x, fOffset %i, fSize %i, fSpec 0x%x", bd.fPtr, bd.fOffset, bd.fSize, bd.fSpecification);
 517           offset = totalSize;
 518         }
 519       }
 520
 521
 522       // if (trdTracks)
 523       //        totalSize += TransportTracks(trdTracks, outputPtr, outputBlocks, offset, dBlockSpecification);
 524       // else {
 525       //          HLTDebug("Bad array trdTracks = 0x%x", trdTracks);
 526       // }
 527
 528       HLTDebug("totalSize: %i", totalSize);
 529
 530 //       if ( totalSize > allocSize )
 531 //      {
 532 //        HLTError("Too much data; Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
 533 //        totalSize, size );
 534 //        return EMSGSIZE;
 535 //      }
 536
 537       //here we are deleting clusters (but not the TClonesArray itself)
 538       fTracker->UnloadClusters();
 539       AliTRDReconstructor::SetClusters(0x0);
 540       delete esd;
 541       clusterArray->Delete();
 542       delete clusterArray;
 543
 544       }
 545
 546   size = totalSize;
 547   HLTDebug("Event is done. size written to the output is %i", size);
 548   return 0;
 549 }