HLT/JET/AliHLTJETReader.cxx

   1 //-*- Mode: C++ -*-
   2 // $Id: AliHLTJETReaderHeader.cxx  $
   3 /**************************************************************************
   4  * This file is property of and copyright by the ALICE HLT Project        *
   5  * ALICE Experiment at CERN, All rights reserved.                         *
   6  *                                                                        *
   7  * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
   8  *                  for The ALICE HLT 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   AliHLTJETReader.cxx
  20     @author Jochen Thaeder
  21     @date
  22     @brief  Reader for jet finder
  23 */
  24
  25 // see header file for class documentation
  26 // or
  27 // refer to README to build package
  28 // or
  29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  30
  31 #if __GNUC__ >= 3
  32 using namespace std;
  33 #endif
  34
  35 #include "TLorentzVector.h"
  36 #include "TParticle.h"
  37 #include "TParticlePDG.h"
  38
  39 #include "AliHLTJETReader.h"
  40
  41 #include "AliHLTJETConeJetCandidate.h"
  42
  43 /** ROOT macro for the implementation of ROOT specific class methods */
  44 ClassImp(AliHLTJETReader)
  45
  46 /*
  47  * ---------------------------------------------------------------------------------
  48  *                            Constructor / Destructor
  49  * ---------------------------------------------------------------------------------
  50  */
  51
  52 // #################################################################################
  53 AliHLTJETReader::AliHLTJETReader()
  54   :
  55   AliJetReader(),
  56   fESD(NULL),
  57   fMC(NULL),
  58   fAOD(NULL),
  59 #ifdef HAVE_FASTJET
  60   fMomentumVector( new vector<fastjet::PseudoJet> ),
  61 #endif
  62   fGrid(NULL),
  63   fNJetCandidates(0),
  64   fJetCandidates(NULL),
  65   fSeedCuts(NULL),
  66   fTrackCuts(NULL) {
  67   // see header file for class documentation
  68   // or
  69   // refer to README to build package
  70   // or
  71   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  72
  73 }
  74
  75 // #################################################################################
  76 AliHLTJETReader::~AliHLTJETReader() {
  77   // see header file for class documentation
  78
  79 #ifdef HAVE_FASTJET
  80   if ( fMomentumVector )
  81     delete fMomentumVector;
  82   fMomentumVector = NULL;
  83 #endif
  84
  85   if ( fGrid )
  86     delete fGrid;
  87   fGrid = NULL;
  88
  89   if ( fJetCandidates ) {
  90     fJetCandidates->Clear();
  91     delete fJetCandidates;
  92   }
  93   fJetCandidates = NULL;
  94
  95 }
  96
  97 // #################################################################################
  98 Int_t AliHLTJETReader::Initialize() {
  99   // see header file for class documentation
 100
 101   Int_t iResult = 0;
 102   AliHLTJETReaderHeader* readerHeader = NULL;
 103
 104   HLTInfo(" -= AliHLTJETReader =- " );
 105
 106   // -- Initialize reader header
 107   // -----------------------------
 108   if ( fReaderHeader ) {
 109     readerHeader = GetReaderHeader();
 110
 111     iResult = readerHeader->Initialize();
 112     if ( iResult )
 113       HLTError("Error initializing HLT jet reader header");
 114   }
 115   else {
 116     HLTError("Reader Header not present");
 117     iResult = -EINPROGRESS;
 118   }
 119
 120   // -- Initialize grid
 121   // --------------------
 122   if ( ! iResult ) {
 123
 124     if ( fGrid )
 125       delete fGrid;
 126
 127     if ( ! (fGrid = new AliHLTJETConeGrid()) ) {
 128       HLTError("Error instanciating grid.");
 129       iResult = -EINPROGRESS;
 130     }
 131   }
 132
 133   if ( ! iResult ) {
 134       fGrid->SetEtaRange(   readerHeader->GetFiducialEtaMin(),
 135                             readerHeader->GetFiducialEtaMax(),
 136                             readerHeader->GetGridEtaRange() );
 137
 138       fGrid->SetPhiRange(   readerHeader->GetFiducialPhiMin(),
 139                             readerHeader->GetFiducialPhiMax(),
 140                             readerHeader->GetGridPhiRange() );
 141
 142       fGrid->SetBinning(    readerHeader->GetGridEtaBinning(),
 143                             readerHeader->GetGridEtaBinning() );
 144
 145       fGrid->SetConeRadius( readerHeader->GetConeRadius() );
 146
 147       iResult = fGrid->Initialize();
 148   }
 149
 150   // -- Initialize jet candidates
 151   // ------------------------------
 152   if ( ! iResult ) {
 153     fJetCandidates = new TClonesArray("AliHLTJETConeJetCandidate", 30);
 154     if ( ! fJetCandidates) {
 155       HLTError("Error instanciating jet candidates.");
 156       iResult = -EINPROGRESS;
 157     }
 158   }
 159
 160   // -- Initialize cuts
 161   // --------------------
 162
 163   // -- Seed cuts
 164   if ( ! iResult ) {
 165     fSeedCuts = readerHeader->GetSeedCuts();
 166     if ( ! fSeedCuts ) {
 167       HLTError("Error getting ptr to seed cuts.");
 168       iResult = -EINPROGRESS;
 169     }
 170     else {
 171       HLTInfo(" -= SeedCuts =- " );
 172     }
 173   }
 174
 175   // -- Track cuts
 176   if ( ! iResult ) {
 177     fTrackCuts = readerHeader->GetTrackCuts();
 178     if ( ! fTrackCuts ) {
 179       HLTError("Error getting ptr to track cuts.");
 180       iResult = -EINPROGRESS;
 181     }
 182   }
 183
 184   // -- Final check
 185   // ----------------
 186   if ( iResult )
 187     HLTError("Error initializing HLT jet reader");
 188
 189   return iResult;
 190 }
 191
 192 //#################################################################################
 193 void AliHLTJETReader::ResetEvent() {
 194   // see header file for class documentation
 195
 196   // -- clear grid
 197   fGrid->Reset();
 198
 199   // -- clear jet candidates
 200   fJetCandidates->Clear();
 201
 202   fNJetCandidates = 0;
 203
 204   return;
 205 }
 206
 207 /*
 208  * ---------------------------------------------------------------------------------
 209  *                            Fastjet Reader functionality
 210  * ---------------------------------------------------------------------------------
 211  */
 212 #ifdef HAVE_FASTJET
 213 // #################################################################################
 214 Bool_t AliHLTJETReader::FillMomentumArrayFast() {
 215   // see header file for class documentation
 216
 217   Bool_t bResult = kFALSE;
 218
 219   if ( fESD )
 220     bResult = FillMomentumArrayFastESD();
 221   else if ( fMC )
 222     bResult = FillMomentumArrayFastMC();
 223   else if ( fAOD )
 224     bResult = FillMomentumArrayFastAOD();
 225
 226   return bResult;
 227 }
 228
 229 // #################################################################################
 230 Bool_t AliHLTJETReader::FillMomentumArrayFastMC() {
 231   // see header file for class documentation
 232
 233   if ( ! fMC ) {
 234     HLTError( "No MC Event present!" );
 235     return kFALSE;
 236   }
 237
 238   // -- Clear input vector
 239   if ( fMomentumVector )
 240     fMomentumVector->clear();
 241
 242   Int_t nTracks = 0;
 243
 244   TParticle* particle = NULL;
 245
 246   // -- Loop over particles
 247   // ------------------------
 248   while ( (particle = fMC->NextParticle() ) ) {
 249
 250     // -- Apply cuts
 251     if ( ! fTrackCuts->IsSelected(particle) )
 252       continue;
 253
 254     // -- Create PseudoJet object
 255     fastjet::PseudoJet part( particle->Px(), particle->Py(),
 256                              particle->Pz(), particle->Energy() );
 257
 258     // -- label the particle into Fastjet algortihm
 259     part.set_user_index( fMC->GetIndex() );
 260
 261     // -- Add to input_particles vector
 262     fMomentumVector->push_back(part);
 263
 264     nTracks++;
 265
 266   } // while ( (particle = fMC->NextParticle() ) ) {
 267
 268   HLTInfo(" Number of selected tracks %d \n", nTracks);
 269
 270   return kTRUE;
 271 }
 272
 273 // #################################################################################
 274 Bool_t AliHLTJETReader::FillMomentumArrayFastESD() {
 275   // see header file for class documentation
 276
 277   return kTRUE;
 278 }
 279
 280 // #################################################################################
 281 Bool_t AliHLTJETReader::FillMomentumArrayFastAOD() {
 282   // see header file for class documentation
 283
 284   return kFALSE;
 285 }
 286 #endif
 287
 288 /*
 289  * ---------------------------------------------------------------------------------
 290  *                               Grid functionality
 291  * ---------------------------------------------------------------------------------
 292  */
 293
 294 // #################################################################################
 295 Bool_t AliHLTJETReader::FillGrid() {
 296   // see header file for class documentation
 297
 298   Bool_t bResult = kFALSE;
 299
 300   if ( fESD )
 301     bResult = FillGridESD();
 302   else if ( fMC )
 303     bResult = FillGridMC();
 304   else if ( fAOD )
 305     bResult = FillGridAOD();
 306
 307   return bResult;
 308 }
 309
 310 // #################################################################################
 311 Bool_t AliHLTJETReader::FillGridMC() {
 312   // see header file for class documentation
 313
 314   if ( ! fMC ) {
 315     HLTError( "No MC Event present!" );
 316     return kFALSE;
 317   }
 318
 319   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
 320
 321   Int_t nTracks = 0;
 322   TParticle* particle = NULL;
 323
 324   // -- Loop over particles
 325   // ------------------------
 326   while ( ( particle = fMC->NextParticle() ) ) {
 327
 328     // -- Apply track cuts
 329     if ( ! fTrackCuts->IsSelected(particle) )
 330       continue;
 331
 332     const Float_t aEtaPhi[]  = { particle->Eta(), particle->Phi(), particle->Pt() };
 333           Int_t   aGridIdx[] = { -1, -1, -1, -1, -1 };
 334
 335     fGrid->FillTrack(particle, aEtaPhi, aGridIdx);
 336
 337     nTracks++;
 338
 339     // -- Apply seed cuts
 340     if ( ! fSeedCuts->IsSelected(particle) )
 341       continue;
 342
 343     // -- Add Seed
 344     AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
 345             readerHeader->GetConeRadius());
 346
 347   } // while ( (particle = fMC->NextParticle() ) ) {
 348
 349   HLTDebug(" Number of selected tracks %d", nTracks);
 350   HLTDebug(" Number of seeds %d", fNJetCandidates);
 351
 352   return kTRUE;
 353 }
 354
 355 // #################################################################################
 356 Bool_t AliHLTJETReader::FillGridESD() {
 357   // see header file for class documentation
 358
 359   Bool_t bResult = kTRUE;
 360
 361   if ( ! fESD ) {
 362     HLTError( "No ESD Event present!" );
 363     return kFALSE;
 364   }
 365
 366   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
 367
 368   Int_t nTracks = 0;
 369
 370   // -- Loop over particles
 371   // ------------------------
 372   for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !bResult; iter++ ) {
 373
 374     AliESDtrack* esdTrack = fESD->GetTrack(iter);
 375     if ( ! esdTrack ) {
 376       HLTError("Could not read ESD track %d from %d\n", iter, fESD->GetNumberOfTracks() );
 377       bResult = kFALSE;
 378       continue;
 379     }
 380
 381     // -- Apply track cuts
 382     if ( ! fTrackCuts->IsSelected(esdTrack) )
 383       continue;
 384
 385     const Float_t aEtaPhi[]  = { esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt() };
 386           Int_t   aGridIdx[] = { -1, -1, -1, -1, -1 };
 387
 388     // -- Fill grid
 389           fGrid->FillTrack(esdTrack, aEtaPhi, aGridIdx);
 390
 391     nTracks++;
 392
 393     // -- Apply seed cuts
 394     if ( ! fSeedCuts->IsSelected(esdTrack) )
 395       continue;
 396
 397     // -- Add Seed
 398     AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
 399             readerHeader->GetConeRadius());
 400
 401   } // for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !iResult; iter++ ) {
 402
 403   HLTDebug(" Number of selected tracks %d", nTracks);
 404   HLTDebug(" Number of seeds %d", fNJetCandidates);
 405
 406   return bResult;
 407 }
 408
 409 // #################################################################################
 410 Bool_t AliHLTJETReader::FillGridAOD() {
 411   // see header file for class documentation
 412
 413   return kFALSE;
 414 }
 415
 416 /*
 417  * ---------------------------------------------------------------------------------
 418  *                                     Setter
 419  * ---------------------------------------------------------------------------------
 420  */
 421
 422 // #################################################################################
 423 void AliHLTJETReader::SetInputEvent(TObject* esd, TObject* aod, TObject* mc) {
 424   // see header file for class documentation
 425
 426   if ( esd )
 427     fESD = dynamic_cast<AliESDEvent*> (esd);
 428   else if ( aod )
 429     fAOD = dynamic_cast<AliAODEvent*> (aod);
 430   else if ( mc )
 431     fMC = dynamic_cast<AliHLTMCEvent*> (mc);
 432
 433   return;
 434 }
 435
 436 /*
 437  * ---------------------------------------------------------------------------------
 438  *                                     Seeds
 439  * ---------------------------------------------------------------------------------
 440  */
 441
 442 //#################################################################################
 443 void AliHLTJETReader::AddSeed( const Float_t* aEtaPhi, const Int_t* aGridIdx,
 444                                const Float_t coneRadius ) {
 445   // see header file for class documentation
 446
 447   Bool_t useWholeCell = kTRUE ; // XXXXXXXXXXXXXXXXx get reader header finder type balhh
 448   useWholeCell = kFALSE ;
 449   // -- Add track / particle
 450
 451   new( (*fJetCandidates) [fNJetCandidates] ) AliHLTJETConeJetCandidate( aEtaPhi,
 452                                                                         aGridIdx,
 453                                                                         coneRadius,
 454                                                                         useWholeCell );
 455   fNJetCandidates++;
 456
 457   HLTDebug("Added Seed Pt=%f, Eta=%f, Phi=%f", aEtaPhi[kIdxPt],
 458           aEtaPhi[kIdxEta], aEtaPhi[kIdxPhi] );
 459
 460   return;
 461 }