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