HLT/JET/fastjet/AliHLTJETFastJetFinder.cxx

   1 //-*- Mode: C++ -*-
   2 // $Id: AliHLTJETFastJetFinder.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   AliHLTJETFastJetFinder.cxx
  20     @author Jochen Thaeder
  21     @date
  22     @brief  FastJet finder interface
  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 //#include "fastjet/AreaDefinition.hh"
  32 //#include "fastjet/JetDefinition.hh"
  33 // get info on how fastjet was configured
  34 //#include "fastjet/config.h"
  35
  36 //#include <TLorentzVector.h>
  37 //#include <TArrayF.h>
  38 //#include <TClonesArray.h>
  39
  40 //#include "AliAODJet.h"
  41
  42 #include "AliHLTJETReader.h"
  43 #include "AliHLTJETJetCuts.h"
  44
  45 #include "AliHLTJETFastJetFinder.h"
  46 #include "AliHLTJETFastJetHeader.h"
  47
  48 //#include<sstream>  // needed for internal io
  49 //#include<vector>
  50 //#include <cmath>
  51
  52 using namespace std;
  53
  54 /** ROOT macro for the implementation of ROOT specific class methods */
  55 ClassImp(AliHLTJETFastJetFinder)
  56
  57 /*
  58  * ---------------------------------------------------------------------------------
  59  *                            Constructor / Destructor
  60  * ---------------------------------------------------------------------------------
  61  */
  62
  63 // #################################################################################
  64 AliHLTJETFastJetFinder::AliHLTJETFastJetFinder()
  65   :
  66   AliJetFinder(),
  67   fInputVector(NULL),
  68   fJets(NULL) {
  69   // see header file for class documentation
  70   // or
  71   // refer to README to build package
  72   // or
  73   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  74
  75 }
  76
  77 // #################################################################################
  78 AliHLTJETFastJetFinder::~AliHLTJETFastJetFinder() {
  79   // see header file for class documentation
  80
  81 }
  82
  83 /*
  84  * ---------------------------------------------------------------------------------
  85  *                                    Initialize
  86  * ---------------------------------------------------------------------------------
  87  */
  88
  89 // #################################################################################
  90 Int_t AliHLTJETFastJetFinder::Initialize() {
  91   // see header file for class documentation
  92
  93   Int_t iResult = 0;
  94
  95   if ( !fHeader || !fReader ) {
  96     HLTError("No header or reader set!");
  97     return -EINPROGRESS;
  98   }
  99
 100   // -- Initialize Reader
 101   AliHLTJETReader *reader = dynamic_cast<AliHLTJETReader*> (fReader);
 102
 103   iResult = reader->Initialize();
 104   if ( iResult ) {
 105     HLTError( "Initializing Reader failed!");
 106     return iResult;
 107   }
 108
 109   // -- Initialize Header
 110   AliHLTJETFastJetHeader *header = dynamic_cast<AliHLTJETFastJetHeader*> (fHeader);
 111
 112   iResult = header->Initialize();
 113   if ( iResult ) {
 114     HLTError( "Initializing Header failed!");
 115     return iResult;
 116   }
 117
 118   // -- Set ptr to vector
 119   fInputVector = reader->GetVector();
 120   if ( ! fInputVector ) {
 121     HLTError( "Getting ptr to vector failed!");
 122     return -EINPROGRESS;
 123   }
 124
 125   // -- Check ptr to output container
 126   if ( !fJets ) {
 127     HLTError( "Ptr to output container not set!");
 128     return -EINPROGRESS;
 129   }
 130
 131   return iResult;
 132 }
 133
 134 // #################################################################################
 135 void AliHLTJETFastJetFinder::Reset() {
 136   // see header file for class documentation
 137
 138   // -- Reset output container
 139   if (fJets)
 140     fJets->Reset();
 141
 142   return;
 143 }
 144
 145 /*
 146  * ---------------------------------------------------------------------------------
 147  *                                      Process
 148  * ---------------------------------------------------------------------------------
 149  */
 150
 151 // #################################################################################
 152 Bool_t AliHLTJETFastJetFinder::ProcessEvent() {
 153   // see header file for class documentation
 154
 155   // -- Pick up jet reader
 156   AliHLTJETReader *reader = dynamic_cast<AliHLTJETReader*> (fReader);
 157
 158   // -- Reset
 159   Reset();
 160
 161   // -- Fill Vector
 162   if ( !reader->FillVector() ){
 163     HLTError("Error filling vector.");
 164     return kFALSE;
 165   }
 166
 167   // -- Find Jets, fill jets and apply jet cuts
 168   if ( FindFastJets() ) {
 169     HLTError("Error finding jets.");
 170     return kFALSE;
 171   }
 172
 173   return kTRUE;
 174 }
 175
 176 // #################################################################################
 177 Bool_t AliHLTJETFastJetFinder::ProcessHLTEvent() {
 178   // see header file for class documentation
 179
 180   // -- Reset
 181   Reset();
 182
 183   // -- Find Jets, fill jets and apply jet cuts
 184   if ( FindFastJets() ) {
 185     HLTError("Error finding jets.");
 186     return kFALSE;
 187   }
 188
 189   return kTRUE;
 190 }
 191
 192 /*
 193  * ---------------------------------------------------------------------------------
 194  *                                      Helper
 195  * ---------------------------------------------------------------------------------
 196  */
 197
 198 // #################################################################################
 199 void AliHLTJETFastJetFinder::PrintJets( vector<fastjet::PseudoJet> &jets,
 200                                   fastjet::ClusterSequenceArea &clust_seq ) {
 201   // see header file for class documentation
 202
 203   // -- pick up fastjet header
 204   AliHLTJETFastJetHeader *header = dynamic_cast<AliHLTJETFastJetHeader*> (fHeader);
 205
 206   // -- print header info
 207   TString comment = header->GetComment();
 208   comment += TString(clust_seq.strategy_string());
 209
 210   HLTInfo( "--------------------------------------------------------" );
 211   HLTInfo( "%s", comment.Data() );
 212   HLTInfo( "--------------------------------------------------------" );
 213
 214   header->PrintParameters();
 215
 216   // -- print found jets
 217   HLTInfo( "Number of unclustered particles: %i", clust_seq.unclustered_particles().size() );
 218
 219   HLTInfo( "Printing inclusive sub jets with pt > %f GeV", header->GetPtMin() );
 220   HLTInfo( "-------------------------------------------------------" );
 221   HLTInfo( " ijet   rap      phi        Pt         area  +-   err" );
 222
 223   for ( UInt_t jetIter = 0; jetIter < jets.size(); jetIter++ ) {
 224
 225     Double_t area       = clust_seq.area(jets[jetIter]);
 226     Double_t area_error = clust_seq.area_error(jets[jetIter]);
 227
 228     HLTInfo( "%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f",
 229             jetIter,jets[jetIter].rap(),jets[jetIter].phi(),jets[jetIter].perp(), area, area_error );
 230
 231   } // for ( Int_t jetIter = 0; jetIter < jets.size(); jetIter++ ) {
 232
 233   return;
 234 }
 235
 236 /*
 237  * ---------------------------------------------------------------------------------
 238  *                             Process - private
 239  * ---------------------------------------------------------------------------------
 240  */
 241
 242 // #################################################################################
 243 Int_t AliHLTJETFastJetFinder::FindFastJets() {
 244   // see header file for class documentation
 245
 246   // -- pick up fastjet header
 247   AliHLTJETFastJetHeader *header = dynamic_cast<AliHLTJETFastJetHeader*> (fHeader);
 248
 249   // -- Run the jet clustering with the jet definition from the header
 250   fastjet::ClusterSequenceArea clust_seq( (*fInputVector),
 251                                           (*header->GetJetDefinition()),
 252                                           (*header->GetAreaDefinition()) );
 253
 254   // -- Extract the inclusive jets with pt > ptmin, sorted by pt
 255   vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(header->GetPtMin());
 256
 257   // -- Subtract background
 258   vector<fastjet::PseudoJet> sub_jets = clust_seq.subtracted_jets((*header->GetRangeDefinition()),
 259                                                                   header->GetPtMin());
 260
 261   // -- Sort jets into increasing pt
 262   vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets);
 263
 264   // -- Fill jets in output container
 265   // -----------------------------------
 266
 267   // -- Get jet cuts
 268   AliHLTJETJetCuts* jetCuts = header->GetJetCuts();
 269
 270   // -- Loop over jet candidates
 271   for ( UInt_t jetIter = 0; jetIter < jets.size(); jetIter++ ) {
 272
 273     // -- Fill AOD jets
 274     AliAODJet aodjet (jets[jetIter].px(), jets[jetIter].py(),
 275                       jets[jetIter].pz(), jets[jetIter].E());
 276
 277     // -- Apply jet cuts
 278     if ( ! jetCuts->IsSelected(&aodjet) )
 279       continue;
 280
 281     fJets->AddJet(&aodjet);
 282
 283   } // for ( Int_t jetIter = 0; jetIter < jets.size(); jetIter++ ) {
 284
 285 #if 0
 286   // -- Print found jets
 287   PrintJets( jets, clust_seq );
 288 #endif
 289
 290   return 0;
 291 }
 292