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