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
  31 #include "TFile.h"
  32 #include "TChain.h"
  33
  34 #include "AliGeomManager.h"
  35 #include "AliCDBManager.h"
  36 #include "AliESDEvent.h"
  37 #include "AliMagFMaps.h"
  38 #include "AliESDfriend.h"
  39
  40 #include "AliTRDcalibDB.h"
  41 #include "AliTRDReconstructor.h"
  42 #include "AliTRDtrackerV1.h"
  43 #include "AliTRDrecoParam.h"
  44
  45 #include <cstdlib>
  46 #include <cerrno>
  47 #include <string>
  48
  49 // this is a global object used for automatic component registration, do not use this
  50 AliHLTTRDTrackerV1Component gAliHLTTRDTrackerV1Component;
  51
  52 ClassImp(AliHLTTRDTrackerV1Component);
  53
  54 AliHLTTRDTrackerV1Component::AliHLTTRDTrackerV1Component()
  55   : AliHLTProcessor()
  56   , fOutputPercentage(100) // By default we copy to the output exactly what we got as input
  57   , fStrorageDBpath("local://$ALICE_ROOT")
  58   , fCDB(NULL)
  59   , fField(NULL)
  60   , fGeometryFileName("")
  61   , fGeometryFile(NULL)
  62   , fGeoManager(NULL)
  63   , fTracker(NULL)
  64   , fRecoParam(NULL)
  65   , fReconstructor(NULL)
  66 {
  67   // Default constructor
  68
  69   fGeometryFileName = getenv("ALICE_ROOT");
  70   fGeometryFileName += "/HLT/TRD/geometry.root";
  71 }
  72
  73 AliHLTTRDTrackerV1Component::~AliHLTTRDTrackerV1Component()
  74 {
  75   // Destructor
  76 }
  77
  78 const char* AliHLTTRDTrackerV1Component::GetComponentID()
  79 {
  80   // Return the component ID const char *
  81   return "TRDTrackerV1"; // The ID of this component
  82 }
  83
  84 void AliHLTTRDTrackerV1Component::GetInputDataTypes( vector<AliHLTComponent_DataType>& list)
  85 {
  86   // Get the list of input data
  87   list.clear(); // We do not have any requirements for our input data type(s).
  88   list.push_back( AliHLTTRDDefinitions::fgkClusterDataType );
  89 }
  90
  91 AliHLTComponent_DataType AliHLTTRDTrackerV1Component::GetOutputDataType()
  92 {
  93   // Get the output data type
  94   return AliHLTTRDDefinitions::fgkClusterDataType;
  95 }
  96
  97 void AliHLTTRDTrackerV1Component::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
  98 {
  99   // Get the output data size
 100   constBase = 0;
 101   inputMultiplier = ((double)fOutputPercentage)/100.0;
 102 }
 103
 104 // Spawn function, return new instance of this class
 105 AliHLTComponent* AliHLTTRDTrackerV1Component::Spawn()
 106 {
 107   return new AliHLTTRDTrackerV1Component;
 108 };
 109
 110 int AliHLTTRDTrackerV1Component::DoInit( int argc, const char** argv )
 111 {
 112   // perform initialization. We check whether our relative output size is specified in the arguments.
 113   fOutputPercentage = 100;
 114   int i = 0;
 115   char* cpErr;
 116
 117
 118   Int_t iRecoParamType = -1; // default will be the low flux
 119   Int_t iNtimeBins = -1;     // number of time bins for the tracker to use
 120   Int_t iMagneticField = -1; // magnetic field: 0==OFF and 1==ON
 121   Bool_t bHLTMode = kTRUE, bWriteClusters = kFALSE;
 122
 123   while ( i < argc )
 124     {
 125       HLTDebug("argv[%d] == %s", i, argv[i] );
 126       if ( !strcmp( argv[i], "output_percentage" ) )
 127         {
 128           if ( i+1>=argc )
 129             {
 130               HLTError("Missing output_percentage parameter");
 131               return ENOTSUP;
 132             }
 133           HLTDebug("argv[%d+1] == %s", i, argv[i+1] );
 134           fOutputPercentage = strtoul( argv[i+1], &cpErr, 0 );
 135           if ( *cpErr )
 136             {
 137               HLTError("Cannot convert output_percentage parameter '%s'", argv[i+1] );
 138               return EINVAL;
 139             }
 140           HLTInfo("Output percentage set to %lu %%", fOutputPercentage );
 141           i += 2;
 142
 143         }
 144
 145       else if ( !strcmp( argv[i], "-NTimeBins" ) )
 146         {
 147           if ( i+1>=argc )
 148             {
 149               HLTError("Missing -NTimeBins parameter");
 150               return ENOTSUP;
 151             }
 152           HLTDebug("Arguments", "argv[%d+1] == %s", i, argv[i+1] );
 153           iNtimeBins = strtoul( argv[i+1], &cpErr, 0 );
 154           if ( *cpErr )
 155             {
 156               HLTError("Wrong Argument. Cannot convert -NTimeBins parameter '%s'", argv[i+1] );
 157               return EINVAL;
 158             }
 159           i += 2;
 160
 161         }
 162
 163       else if ( strcmp( argv[i], "-cdb" ) == 0)
 164         {
 165           if ( i+1 >= argc )
 166             {
 167               HLTError( "Missing -cdb argument");
 168               return ENOTSUP;
 169             }
 170           fStrorageDBpath = argv[i+1];
 171           HLTInfo("DB storage is %s", fStrorageDBpath.c_str() );
 172           i += 2;
 173
 174         }
 175
 176       else if ( strcmp( argv[i], "-geometry" ) == 0)
 177         {
 178           if ( i+1 >= argc )
 179             {
 180               HLTError("Missing -geometry argument");
 181               return ENOTSUP;
 182             }
 183           fGeometryFileName = argv[i+1];
 184           HLTInfo("GeomFile storage is %s",
 185                   fGeometryFileName.c_str() );
 186           i += 2;
 187
 188         }
 189
 190       // the flux parametrizations
 191       else if ( strcmp( argv[i], "-lowflux" ) == 0)
 192         {
 193           iRecoParamType = 0;
 194           HLTDebug("Low flux reco selected.");
 195           i++;
 196
 197         }
 198
 199       else if ( strcmp( argv[i], "-highflux" ) == 0)
 200         {
 201           iRecoParamType = 1;
 202           HLTDebug("Low flux reco selected.");
 203           i++;
 204         }
 205
 206       else if ( strcmp( argv[i], "-cosmics" ) == 0)
 207         {
 208           iRecoParamType = 2;
 209           HLTDebug("Cosmic test reco selected.");
 210           i++;
 211         }
 212
 213       else if ( strcmp( argv[i], "-magnetic_field_ON" ) == 0)
 214         {
 215           iMagneticField = 1;
 216           i++;
 217         }
 218       else if ( strcmp( argv[i], "-magnetic_field_OFF" ) == 0)
 219         {
 220           iMagneticField = 0;
 221           i++;
 222         }
 223       else if ( strcmp( argv[i], "-writeClusters" ) == 0)
 224         {
 225           bWriteClusters = kTRUE;
 226           HLTDebug("input clusters are expected to be in a TTree.");
 227           i++;
 228         }
 229       else if ( strcmp( argv[i], "-offlineMode" ) == 0)
 230         {
 231           bHLTMode=kFALSE;
 232           HLTDebug("Using standard offline tracking.");
 233           i++;
 234         }
 235       else {
 236         HLTError("Unknown option '%s'", argv[i] );
 237         return EINVAL;
 238       }
 239
 240     }
 241
 242   // THE "REAL" INIT COMES HERE
 243   // offline condition data base
 244   fCDB = AliCDBManager::Instance();
 245   if (!fCDB)
 246     {
 247       HLTError("Could not get CDB instance", "fCDB 0x%x", fCDB);
 248       return -1;
 249     }
 250   else
 251     {
 252       fCDB->SetRun(0); // THIS HAS TO BE RETRIEVED !!!
 253       fCDB->SetDefaultStorage(fStrorageDBpath.c_str());
 254       HLTDebug("CDB instance", "fCDB 0x%x", fCDB);
 255     }
 256
 257   // check if the N of time bins make sense
 258   if (iNtimeBins <= 0)
 259     {
 260       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.");
 261       return -1;
 262     }
 263
 264   if (iNtimeBins < 24 || iNtimeBins > 30)
 265     {
 266       HLTWarning("The number of time bins seems to be strange = %d. But okay. Let's try it...", iNtimeBins);
 267     }
 268
 269   HLTDebug("The number of time bins = %d.", iNtimeBins);
 270   AliTRDtrackerV1::SetNTimeBins(iNtimeBins);
 271
 272   // !!!! THIS IS IMPORTANT
 273   // init alifield map - temporarly via parameter - should come from a DB or DCS ?
 274   // !!!!
 275   if (iMagneticField < 0)
 276     {
 277       iMagneticField = 0;
 278       HLTWarning("No magnetic field switch stated. Use -magnetic_field_ON or -magnetic_field_OFF flag. Defaulting to OFF = NO MAGNETIC FIELD");
 279     }
 280
 281   if (iMagneticField == 0)
 282     {
 283       // magnetic field OFF
 284       fField = new AliMagFMaps("Maps","Maps", 2, 0., 10., 1);
 285       HLTDebug("Magnetic field is OFF.");
 286     }
 287
 288   if (iMagneticField == 1)
 289     {
 290       // magnetic field ON
 291       fField = new AliMagFMaps("Maps","Maps", 2, 1., 10., 1);
 292       HLTDebug("Magnetic field is ON.");
 293     }
 294
 295   if (fField == 0)
 296     {
 297       HLTError("Unable to init the field. Trouble at this point.");
 298       return -1;
 299     }
 300
 301   // kTRUE sets the map uniform
 302   AliTracker::SetFieldMap(fField,kTRUE);
 303
 304   // reconstruction parameters
 305   if (iRecoParamType < 0 || iRecoParamType > 2)
 306     {
 307       HLTWarning("No reco param selected. Use -lowflux -highflux -cosmics flags. Defaulting to low flux.");
 308       iRecoParamType = 0;
 309     }
 310
 311   if (iRecoParamType == 0)
 312     {
 313       fRecoParam = AliTRDrecoParam::GetLowFluxParam();
 314       HLTDebug("Low flux params init.");
 315     }
 316
 317   if (iRecoParamType == 1)
 318     {
 319       fRecoParam = AliTRDrecoParam::GetHighFluxParam();
 320       HLTDebug("High flux params init.");
 321     }
 322
 323   if (iRecoParamType == 2)
 324     {
 325       fRecoParam = AliTRDrecoParam::GetCosmicTestParam();
 326       HLTDebug("Cosmic Test params init.");
 327     }
 328
 329   if (fRecoParam == 0)
 330     {
 331       HLTError("No reco params initialized. Sniffing big trouble!");
 332       return -1;
 333     }
 334
 335   fReconstructor = new AliTRDReconstructor();
 336 //   fRecoParam->SetChi2Y(.1);
 337 //   fRecoParam->SetChi2Z(5.);
 338   fReconstructor->SetRecoParam(fRecoParam);
 339   // write clusters [cw] = true
 340   // track seeding (stand alone tracking) [sa] = true
 341   // PID method in reconstruction (NN) [nn] = true
 342   // write online tracklets [tw] = false
 343   // drift gas [ar] = false
 344   // sl_tr_0 = StreamLevel_task_Level
 345   //  fReconstructor->SetOption("sa,!cw,hlt,sl_tr_0");
 346   TString recoOptions="sa,sl_cf_0";
 347
 348   if (bWriteClusters)
 349     {
 350       recoOptions += ",cw";
 351     }
 352   else
 353     {
 354       recoOptions += ",!cw";
 355     }
 356   if (bHLTMode)
 357     recoOptions += ",hlt";
 358
 359   fReconstructor->SetOption(recoOptions.Data());
 360   HLTDebug("Reconstructor options are: %s",recoOptions.Data());
 361
 362   fGeometryFile = 0;
 363   fGeometryFile = TFile::Open(fGeometryFileName.c_str());
 364   if (fGeometryFile)
 365     {
 366       AliGeomManager::LoadGeometry(fGeometryFileName.c_str());
 367     }
 368   else
 369     {
 370       HLTError("Unable to open file. FATAL!");
 371       return -1;
 372     }
 373
 374   // create the tracker
 375   fTracker = new AliTRDtrackerV1();
 376   fTracker->SetReconstructor(fReconstructor);
 377   HLTDebug("TRDTracker at 0x%x", fTracker);
 378
 379   if (fTracker == 0)
 380     {
 381       HLTError("Unable to create the tracker!");
 382       return -1;
 383     }
 384
 385   return 0;
 386 }
 387
 388 int AliHLTTRDTrackerV1Component::DoDeinit()
 389 {
 390   // Deinitialization of the component
 391
 392   delete fField;
 393   fField = 0x0;
 394
 395   // fTracker->SetClustersOwner(kFALSE);
 396   delete fTracker;
 397   fTracker = 0x0;
 398
 399   // We need to set clusters in Reconstructor to null to prevent from
 400   // double deleting, since we delete TClonesArray by ourself in DoEvent.
 401   fReconstructor->SetClusters(0x0);
 402   delete fReconstructor;
 403   fReconstructor = 0x0;
 404
 405   AliTRDcalibDB::Terminate();
 406
 407   if (fGeometryFile)
 408     {
 409       fGeometryFile->Close();
 410       delete fGeometryFile;
 411       fGeometryFile = 0x0;
 412     }
 413
 414   return 0;
 415 }
 416
 417 int AliHLTTRDTrackerV1Component::DoEvent( const AliHLTComponentEventData & /*evtData*/,
 418                                           AliHLTComponentTriggerData & /*trigData*/ )
 419 {
 420   // Process an event
 421   Bool_t bWriteClusters = fReconstructor->IsWritingClusters();
 422
 423   HLTDebug("Output percentage set to %lu %%", fOutputPercentage );
 424   HLTDebug("NofBlocks %lu", GetNumberOfInputBlocks() );
 425
 426   AliHLTUInt32_t dBlockSpecification = 0;
 427
 428
 429   //implement a usage of the following
 430   //   AliHLTUInt32_t triggerDataStructSize = trigData.fStructSize;
 431   //   AliHLTUInt32_t triggerDataSize = trigData.fDataSize;
 432   //   void *triggerData = trigData.fData;
 433   //HLTDebug("Struct size %d Data size %d Data location 0x%x", trigData.fStructSize, trigData.fDataSize, (UInt_t*)trigData.fData);
 434
 435   AliHLTComponentBlockData *dblock = (AliHLTComponentBlockData *)GetFirstInputBlock( AliHLTTRDDefinitions::fgkClusterDataType );
 436   if (dblock != 0)
 437     {
 438       AliHLTComponentDataType inputDataType = dblock->fDataType;
 439       HLTDebug( "Event 0x%08LX (%Lu) received datatype: %s",
 440                 GetEventId(), GetEventId(),
 441                 DataType2Text(inputDataType).c_str());
 442       dBlockSpecification = dblock->fSpecification;
 443     }
 444   else
 445     {
 446       HLTWarning("First Input Block not found! 0x%x", dblock);
 447       return 0;
 448     }
 449
 450   int ibForce = 0;
 451
 452   TObject *tobjin = 0x0;
 453
 454   TTree *clusterTree = 0x0;
 455   TClonesArray *clusterArray = 0x0;
 456   if (bWriteClusters){
 457     tobjin = (TObject *)GetFirstInputObject( AliHLTTRDDefinitions::fgkClusterDataType, "TTree", ibForce);
 458     HLTDebug("1stBLOCK; Pointer = 0x%x", tobjin);
 459     clusterTree = (TTree*)tobjin;
 460
 461     //     while (tobjin != 0)
 462     //         {
 463     if (clusterTree)
 464       {
 465         HLTDebug("CLUSTERS; Pointer to TTree = 0x%x Name = %s", clusterTree, clusterTree->GetName());
 466         HLTDebug("TTree of clusters: nbEntries = %i", clusterTree->GetEntriesFast());
 467         fTracker->LoadClusters(clusterTree);
 468       }
 469     else
 470       {
 471         HLTError("First Input Block not a TTree 0x%x", tobjin);
 472       }
 473   }
 474   else{
 475     tobjin = (TObject *)GetFirstInputObject( AliHLTTRDDefinitions::fgkClusterDataType, "TClonesArray", ibForce);
 476     HLTDebug("1stBLOCK; Pointer = 0x%x", tobjin);
 477     clusterArray = (TClonesArray*)tobjin;
 478     if (clusterArray)
 479       {
 480         HLTDebug("CLUSTERS; Pointer to TClonesArray = 0x%x Name = %s", clusterArray, clusterArray->GetName());
 481         HLTDebug("TClonesArray of clusters: nbEntries = %i", clusterArray->GetEntriesFast());
 482         Int_t nb = clusterArray->GetEntriesFast();
 483         for (Int_t i=0; i<nb; i++){
 484           //AliTRDcluster * cluster = (AliTRDcluster* ) clusterArray->At(i);
 485           //HLTDebug("Cluster[%i]: detector %i", i, cluster->GetDetector());
 486         }
 487
 488         fTracker->LoadClusters(clusterArray);
 489       }
 490     else
 491       {
 492         HLTError("First Input Block not a TClonesArray 0x%x", tobjin);
 493       }
 494   }
 495
 496   // maybe it is not so smart to create it each event? clear is enough ?
 497   AliESDEvent *esd = new AliESDEvent();
 498   esd->CreateStdContent();
 499   fTracker->Clusters2Tracks(esd);
 500
 501   //here transport the esd tracks further
 502   Int_t nTracks = esd->GetNumberOfTracks();
 503   Int_t nTRDTracks = esd->GetNumberOfTrdTracks();
 504   HLTInfo("Number of tracks  == %d == Number of TRD tracks %d", nTracks, nTRDTracks);
 505   //esd->Print();
 506   //PushBack(esd, AliHLTTRDDefinitions::fgkTRDSAEsdDataType);
 507
 508   // extract the friend ?
 509   //   AliESDfriend *esdFriend = new AliESDfriend();
 510   //   esd->GetESDfriend(esdFriend);
 511   //   PushBack(esdFriend, AliHLTTRDDefinitions::fgkTRDSAEsdDataType);
 512   //   delete esdFriend;
 513
 514    TClonesArray* trdTracks = fTracker->GetListOfTracks();
 515
 516    if (trdTracks)
 517     {
 518       nTracks=trdTracks->GetEntriesFast();
 519       if (nTracks>0)
 520         {
 521           HLTDebug("We have an output array: pointer to trdTracks = 0x%x, nbEntries = %i", trdTracks, nTracks);
 522       HLTDebug("We have an output array: pointer to trdTracks = 0x%x, nbEntries = %i", trdTracks, nTracks);
 523
 524       trdTracks->BypassStreamer(kFALSE);
 525       PushBack(trdTracks, AliHLTTRDDefinitions::fgkTRDSATracksDataType);
 526       //trdTracks->Delete();
 527       //delete trdTracks;
 528         }
 529     }
 530   else
 531     HLTDebug("Bad array trdTracks = 0x%x", trdTracks);
 532
 533
 534
 535   //here we are deleting clusters (but not the TClonesArray itself)
 536   fTracker->UnloadClusters();
 537   delete esd;
 538   if (bWriteClusters)
 539     delete clusterTree;
 540   else
 541     delete clusterArray;
 542
 543   HLTDebug("Event done.");
 544   return 0;
 545 }
 546
 547 ///////////////////////////////
 548 /*
 549   consider transporting TRD tracks only as they might have signigicantly smaller size... on the other hand you will need to prodece ESDs at some point...
 550
 551   // this is for ESDtrack
 552   //   for (Int_t it = 0; it < nTracks; it++)
 553   //     {
 554   //       AliESDtrack* track = esd->GetTrack(it);
 555   //       HLTInfo("Track %d 0x%x Pt %1.2f", it, track, track->Pt());
 556   //       //PushBack(track, AliHLTTRDDefinitions::fgkTRDSATracksDataType, ++dBlockSpecification);
 557   //       PushBack(track, AliHLTTRDDefinitions::fgkTRDSATracksDataType);
 558   //     }
 559
 560   // one can do similar things with the TRDtrack
 561   esd->GetNumberOfTrdTracks();
 562   and then
 563   for (i;;)
 564   AliESDTrdTrack *trdtrack = esd->GetTrdTrack(i)
 565 */
 566