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. *
7 * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de> *
8 * for The ALICE HLT Project. *
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 **************************************************************************/
19 /** @file AliHLTJETReader.cxx
20 @author Jochen Thaeder
22 @brief Reader for jet finder
25 // see header file for class documentation
27 // refer to README to build package
29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
35 #include "TLorentzVector.h"
36 #include "TParticle.h"
37 #include "TParticlePDG.h"
38 #include "TDatabasePDG.h"
40 #include "AliHLTJETReader.h"
42 #include "AliHLTJETConeJetCandidate.h"
44 /** ROOT macro for the implementation of ROOT specific class methods */
45 ClassImp(AliHLTJETReader)
48 * ---------------------------------------------------------------------------------
49 * Constructor / Destructor
50 * ---------------------------------------------------------------------------------
53 // #################################################################################
54 AliHLTJETReader::AliHLTJETReader()
62 fMomentumVector(NULL),
69 // see header file for class documentation
71 // refer to README to build package
73 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
77 // #################################################################################
78 AliHLTJETReader::~AliHLTJETReader() {
79 // see header file for class documentation
82 if ( fMomentumVector )
83 delete fMomentumVector;
84 fMomentumVector = NULL;
91 if ( fJetCandidates ) {
92 fJetCandidates->Clear();
93 delete fJetCandidates;
95 fJetCandidates = NULL;
99 // #################################################################################
100 Int_t AliHLTJETReader::Initialize() {
101 // see header file for class documentation
104 AliHLTJETReaderHeader* readerHeader = NULL;
106 HLTInfo(" -= AliHLTJETReader =- " );
108 // -- Initialize reader header
109 // -----------------------------
110 if ( fReaderHeader ) {
111 readerHeader = GetReaderHeader();
113 iResult = readerHeader->Initialize();
115 HLTError("Error initializing HLT jet reader header");
118 HLTError("Reader Header not present");
119 iResult = -EINPROGRESS;
122 // -- Initialize Algorithms
123 // --------------------------
124 if ( readerHeader->GetJetAlgorithm() >= AliHLTJETBase::kFFSCSquareCell )
125 iResult = InitializeFFSC();
128 iResult = InitializeFastjet();
130 HLTError("Error FastJet not present.");
131 iResult = -EINPROGRESS;
135 // -- Get ptr to cuts from reader
136 // --------------------------------
140 fTrackCuts = readerHeader->GetTrackCuts();
141 if ( ! fTrackCuts ) {
142 HLTError("Error getting ptr to track cuts.");
143 iResult = -EINPROGRESS;
150 HLTError("Error initializing HLT jet reader");
155 //#################################################################################
156 void AliHLTJETReader::ResetEvent() {
157 // see header file for class documentation
159 // -- Reset FFSC algorithms
160 // --------------------------
161 if ( GetReaderHeader()->GetJetAlgorithm() >= AliHLTJETBase::kFFSCSquareCell ) {
165 // -- clear jet candidates
166 fJetCandidates->Clear();
171 // -- Reset for FastJet algorithms
172 // ---------------------------------
175 // -- Clear input vector
176 if ( fMomentumVector )
177 fMomentumVector->clear();
185 * ---------------------------------------------------------------------------------
187 * ---------------------------------------------------------------------------------
190 // #################################################################################
191 void AliHLTJETReader::SetInputEvent(const TObject* esd, const TObject* aod, const TObject* mc) {
192 // see header file for class documentation
194 // Needs "useMC" flag for running in analysis task only
196 AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
199 if ( esd && !readerHeader->GetUseMC() )
200 fESD = dynamic_cast<AliESDEvent*> (const_cast<TObject*>(esd));
203 else if ( aod && !readerHeader->GetUseMC() )
204 fAOD = dynamic_cast<AliAODEvent*> (const_cast<TObject*>(aod));
207 else if ( mc && readerHeader->GetUseMC() ) {
209 // -- if off-line MC event,
210 if ( !strcmp (mc->ClassName(),"AliMCEvent") )
211 fMC = dynamic_cast<AliMCEvent*> (const_cast<TObject*>(mc));
212 // -- if on-line MC event
214 fHLTMC = dynamic_cast<AliHLTMCEvent*> (const_cast<TObject*>(mc));
221 * ---------------------------------------------------------------------------------
222 * Fastjet Reader functionality
223 * ---------------------------------------------------------------------------------
226 // #################################################################################
227 Bool_t AliHLTJETReader::FillVector() {
228 // see header file for class documentation
230 Bool_t bResult = kFALSE;
233 bResult = FillVectorESD();
235 bResult = FillVectorMC();
237 bResult = FillVectorHLTMC();
239 bResult = FillVectorAOD();
244 // #################################################################################
245 Bool_t AliHLTJETReader::FillVectorMC() {
246 // see header file for class documentation
248 Bool_t bResult = kTRUE;
251 HLTError( "No MC Event present!" );
260 AliStack* stack = fMC->Stack();
262 for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && bResult; iterStack++) {
264 TParticle *particle = stack->Particle(iterStack);
266 HLTError( "Error reading particle %i out of %i", iterStack, stack->GetNtrack() );
271 // ------------------------------
272 // -- Basic cuts on MC particle --> To be done better XXX
273 // ------------------------------
276 if ( !(stack->IsPhysicalPrimary(iterStack)) )
280 if ( particle->GetNDaughters() != 0 )
284 TParticlePDG * particlePDG = particle->GetPDG();
285 if ( ! particlePDG ) {
286 particlePDG = TDatabasePDG::Instance()->GetParticle( particle->GetPdgCode() );
288 if ( ! particlePDG ) {
289 HLTError("Particle %i not in PDG database", particle->GetPdgCode() );
295 // ------------------------
296 // -- Standard track cuts
297 // ------------------------
299 // -- Apply track cuts
300 if ( ! fTrackCuts->IsSelected(particle) )
303 // -- Create PseudoJet object
304 fastjet::PseudoJet part( particle->Px(), particle->Py(),
305 particle->Pz(), particle->Energy() );
307 // -- label the particle into Fastjet algortihm
308 part.set_user_index( iterStack );
311 // -- Add to input_particles vector
312 fMomentumVector->push_back(part);
316 } // for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && !bResult; iterStack++) {
318 HLTDebug(" Number of selected tracks %d", nTracks);
323 // #################################################################################
324 Bool_t AliHLTJETReader::FillVectorHLTMC() {
325 // see header file for class documentation
328 HLTError( "No HLT MC Event present!" );
337 TParticle* particle = NULL;
339 // -- Loop over particles
340 // ------------------------
341 while ( (particle = fHLTMC->NextParticle() ) ) {
343 // -- Apply track cuts
344 if ( ! fTrackCuts->IsSelected(particle) )
347 // -- Create PseudoJet object
348 fastjet::PseudoJet part( particle->Px(), particle->Py(),
349 particle->Pz(), particle->Energy() );
351 // -- label the particle into Fastjet algortihm
352 part.set_user_index( fHLTMC->GetIndex() );
355 // -- Add to input_particles vector
356 fMomentumVector->push_back(part);
361 } // while ( (particle = fHLTMC->NextParticle() ) ) {
363 HLTInfo(" Number of selected tracks %d", nTracks);
368 // #################################################################################
369 Bool_t AliHLTJETReader::FillVectorESD() {
370 // see header file for class documentation
372 Bool_t bResult = kTRUE;
375 HLTError( "No ESD Event present!" );
384 // -- Loop over particles
385 // ------------------------
386 for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && bResult; iter++ ) {
388 AliESDtrack* esdTrack = fESD->GetTrack(iter);
390 HLTError("Could not read ESD track %d from %d", iter, fESD->GetNumberOfTracks() );
395 // -- Apply track cuts
396 if ( ! fTrackCuts->IsSelected(esdTrack) )
399 // -- Create PseudoJet object
400 fastjet::PseudoJet part( esdTrack->Px(), esdTrack->Py(),
401 esdTrack->Pz(), esdTrack->E() );
403 // -- label the particle into Fastjet algortihm
404 part.set_user_index( iter );
407 // -- Add to input_particles vector
408 fMomentumVector->push_back(part);
413 } // for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !iResult; iter++ ) {
415 HLTInfo(" Number of selected tracks %d", nTracks);
420 // #################################################################################
421 Bool_t AliHLTJETReader::FillVectorAOD() {
422 // see header file for class documentation
429 * ---------------------------------------------------------------------------------
431 * ---------------------------------------------------------------------------------
434 // #################################################################################
435 Bool_t AliHLTJETReader::FillGrid() {
436 // see header file for class documentation
438 Bool_t bResult = kFALSE;
441 bResult = FillGridESD();
443 bResult = FillGridMC();
445 bResult = FillGridHLTMC();
447 bResult = FillGridAOD();
452 // #################################################################################
453 Bool_t AliHLTJETReader::FillGridMC() {
454 // see header file for class documentation
456 Bool_t bResult = kTRUE;
459 HLTError( "No MC Event present!" );
466 AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
470 AliStack* stack = fMC->Stack();
472 for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && bResult; iterStack++) {
474 TParticle *particle = stack->Particle(iterStack);
476 HLTError( "Error reading particle %i out of %i", iterStack, stack->GetNtrack() );
481 // ------------------------------
482 // -- Basic cuts on MC particle --> To be done better XXX
483 // ------------------------------
486 if ( !(stack->IsPhysicalPrimary(iterStack)) )
490 if ( particle->GetNDaughters() != 0 )
494 TParticlePDG * particlePDG = particle->GetPDG();
495 if ( ! particlePDG ) {
496 particlePDG = TDatabasePDG::Instance()->GetParticle( particle->GetPdgCode() );
498 if ( ! particlePDG ) {
499 HLTError("Particle %i not in PDG database", particle->GetPdgCode() );
505 // ------------------------
506 // -- Standard track cuts
507 // ------------------------
509 // -- Apply track cuts
510 if ( ! fTrackCuts->IsSelected(particle) )
513 const Float_t aEtaPhi[] = { particle->Eta(), particle->Phi(), particle->Pt() };
514 Int_t aGridIdx[] = { -1, -1, -1, -1, -1 };
516 fGrid->FillTrack(particle, aEtaPhi, aGridIdx);
520 // -- Apply seed cuts
521 if ( ! fSeedCuts->IsSelected(particle) )
525 AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx), readerHeader->GetConeRadius());
527 } // for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && !bResult; iterStack++) {
529 HLTInfo(" Number of selected tracks %d", nTracks);
530 HLTInfo(" Number of seeds %d", fNJetCandidates);
535 // #################################################################################
536 Bool_t AliHLTJETReader::FillGridHLTMC() {
537 // see header file for class documentation
540 HLTError( "No HLT MC Event present!" );
547 AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
550 TParticle* particle = NULL;
552 // -- Loop over particles
553 // ------------------------
554 while ( ( particle = fHLTMC->NextParticle() ) ) {
556 // HLTError("=== nTracks %d ===",nTracks);
558 // -- Apply track cuts
559 if ( ! fTrackCuts->IsSelected(particle) )
562 const Float_t aEtaPhi[] = { particle->Eta(), particle->Phi(), particle->Pt() };
563 Int_t aGridIdx[] = { -1, -1, -1, -1, -1 };
566 fGrid->FillTrack(particle, aEtaPhi, aGridIdx);
570 // -- Apply seed cuts
571 if ( ! fSeedCuts->IsSelected(particle) )
575 AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
576 readerHeader->GetConeRadius());
578 } // while ( (particle = fHLTMC->NextParticle() ) ) {
580 HLTInfo(" Number of selected tracks %d", nTracks);
581 HLTInfo(" Number of seeds %d", fNJetCandidates);
586 // #################################################################################
587 Bool_t AliHLTJETReader::FillGridESD() {
588 // see header file for class documentation
590 Bool_t bResult = kTRUE;
593 HLTError( "No ESD Event present!" );
600 AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
604 // -- Loop over particles
605 // ------------------------
606 for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !bResult; iter++ ) {
608 AliESDtrack* esdTrack = fESD->GetTrack(iter);
610 HLTError("Could not read ESD track %d from %d", iter, fESD->GetNumberOfTracks() );
615 // -- Apply track cuts
616 if ( ! fTrackCuts->IsSelected(esdTrack) )
619 const Float_t aEtaPhi[] = { esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt() };
620 Int_t aGridIdx[] = { -1, -1, -1, -1, -1 };
623 fGrid->FillTrack(esdTrack, aEtaPhi, aGridIdx);
627 // -- Apply seed cuts
628 if ( ! fSeedCuts->IsSelected(esdTrack) )
632 AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
633 readerHeader->GetConeRadius());
635 } // for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !iResult; iter++ ) {
637 HLTInfo(" Number of selected tracks %d", nTracks);
638 HLTInfo(" Number of seeds %d", fNJetCandidates);
643 // #################################################################################
644 Bool_t AliHLTJETReader::FillGridAOD() {
645 // see header file for class documentation
651 * ---------------------------------------------------------------------------------
653 * ---------------------------------------------------------------------------------
656 //#################################################################################
657 void AliHLTJETReader::AddSeed( const Float_t* aEtaPhi, const Int_t* aGridIdx,
658 const Float_t coneRadius ) {
659 // see header file for class documentation
661 Bool_t useWholeCell = kTRUE;
663 if ( GetReaderHeader()->GetJetAlgorithm() == AliHLTJETBase::kFFSCRadiusCell )
664 useWholeCell = kFALSE ;
666 // -- Add track / particle
667 new( (*fJetCandidates) [fNJetCandidates] ) AliHLTJETConeJetCandidate( aEtaPhi,
673 HLTDebug("Added Seed Pt=%f, Eta=%f, Phi=%f", aEtaPhi[kIdxPt],
674 aEtaPhi[kIdxEta], aEtaPhi[kIdxPhi] );
680 * ---------------------------------------------------------------------------------
681 * Initialize - private
682 * ---------------------------------------------------------------------------------
685 // #################################################################################
686 Int_t AliHLTJETReader::InitializeFFSC() {
687 // see header file for class documentation
690 AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
692 // -- Initialize grid
693 // --------------------
697 if ( ! (fGrid = new AliHLTJETConeGrid()) ) {
698 HLTError("Error instanciating grid.");
699 iResult = -EINPROGRESS;
703 fGrid->SetEtaRange( readerHeader->GetFiducialEtaMin(),
704 readerHeader->GetFiducialEtaMax(),
705 readerHeader->GetGridEtaRange() );
707 fGrid->SetPhiRange( readerHeader->GetFiducialPhiMin(),
708 readerHeader->GetFiducialPhiMax(),
709 readerHeader->GetGridPhiRange() );
711 fGrid->SetBinning( readerHeader->GetGridEtaBinning(),
712 readerHeader->GetGridEtaBinning() );
714 fGrid->SetConeRadius( readerHeader->GetConeRadius() );
716 iResult = fGrid->Initialize();
719 // -- Initialize jet candidates
720 // ------------------------------
722 fJetCandidates = new TClonesArray("AliHLTJETConeJetCandidate", 30);
723 if ( ! fJetCandidates) {
724 HLTError("Error instanciating jet candidates.");
725 iResult = -EINPROGRESS;
729 // -- Get ptr to cuts from reader
730 // --------------------------------
734 fSeedCuts = readerHeader->GetSeedCuts();
736 HLTError("Error getting ptr to seed cuts.");
737 iResult = -EINPROGRESS;
745 // #################################################################################
746 Int_t AliHLTJETReader::InitializeFastjet() {
747 // see header file for class documentation
751 // -- Initialize Vector
752 // ----------------------
753 if ( fMomentumVector )
754 delete fMomentumVector;
756 if ( ! (fMomentumVector = new vector<fastjet::PseudoJet>) ) {
757 HLTError("Error instanciating momentum vector.");
758 iResult = -EINPROGRESS;