Update handling of different inputs
[u/mrichter/AliRoot.git] / 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 #include "TDatabasePDG.h"
39
40 #include "AliHLTJETReader.h"
41
42 #include "AliHLTJETConeJetCandidate.h"
43
44 /** ROOT macro for the implementation of ROOT specific class methods */
45 ClassImp(AliHLTJETReader)
46
47 /*
48  * ---------------------------------------------------------------------------------
49  *                            Constructor / Destructor
50  * ---------------------------------------------------------------------------------
51  */
52   
53 // #################################################################################
54 AliHLTJETReader::AliHLTJETReader()
55   : 
56   AliJetReader(),
57   fESD(NULL), 
58   fMC(NULL),
59   fHLTMC(NULL),
60   fAOD(NULL),
61 #ifdef HAVE_FASTJET
62   fMomentumVector(NULL),
63 #endif
64   fGrid(NULL),
65   fNJetCandidates(0),
66   fJetCandidates(NULL),
67   fSeedCuts(NULL),
68   fTrackCuts(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 AliHLTJETReader::~AliHLTJETReader() {
79   // see header file for class documentation
80
81 #ifdef HAVE_FASTJET
82   if ( fMomentumVector )
83     delete fMomentumVector;
84   fMomentumVector = NULL;
85 #endif
86
87   if ( fGrid )
88     delete fGrid;
89   fGrid = NULL;
90
91   if ( fJetCandidates ) {
92     fJetCandidates->Clear();
93     delete fJetCandidates;
94   }
95   fJetCandidates = NULL;
96
97 }
98
99 // #################################################################################
100 Int_t AliHLTJETReader::Initialize() {
101   // see header file for class documentation
102
103   Int_t iResult = 0;
104   AliHLTJETReaderHeader* readerHeader = NULL;
105
106   HLTInfo(" -= AliHLTJETReader =- " );
107
108   // -- Initialize reader header
109   // -----------------------------
110   if ( fReaderHeader ) {
111     readerHeader = GetReaderHeader();
112
113     iResult = readerHeader->Initialize();
114     if ( iResult )
115       HLTError("Error initializing HLT jet reader header");
116   }
117   else {
118     HLTError("Reader Header not present");
119     iResult = -EINPROGRESS;
120   }
121   
122   // -- Initialize Algorithms
123   // --------------------------
124   if ( readerHeader->GetJetAlgorithm() >= AliHLTJETBase::kFFSCSquareCell )
125     iResult = InitializeFFSC();
126   else {
127 #ifdef HAVE_FASTJET
128     iResult = InitializeFastjet();
129 #else  
130     HLTError("Error FastJet not present.");
131     iResult = -EINPROGRESS;
132 #endif
133   }
134  
135   // -- Get ptr to cuts from reader
136   // --------------------------------
137
138   // -- Track cuts
139   if ( ! iResult ) {
140     fTrackCuts = readerHeader->GetTrackCuts();
141     if ( ! fTrackCuts ) {
142       HLTError("Error getting ptr to track cuts.");
143       iResult = -EINPROGRESS;
144     }
145   }
146
147   // -- Final check
148   // ----------------
149   if ( iResult )
150     HLTError("Error initializing HLT jet reader");
151  
152   return iResult;
153 }
154
155 //#################################################################################
156 void AliHLTJETReader::ResetEvent() {
157   // see header file for class documentation
158
159   // -- Reset FFSC algorithms
160   // --------------------------
161   if (  GetReaderHeader()->GetJetAlgorithm() >= AliHLTJETBase::kFFSCSquareCell ) {
162     // -- clear grid
163     fGrid->Reset();
164     
165     // -- clear jet candidates
166     fJetCandidates->Clear();
167     
168     fNJetCandidates = 0;
169   }
170
171   // -- Reset for FastJet algorithms
172   // ---------------------------------
173   else {
174 #ifdef HAVE_FASTJET
175     // -- Clear input vector
176     if ( fMomentumVector )
177       fMomentumVector->clear();
178 #endif
179   }
180
181   return;  
182 }
183
184 /*
185  * ---------------------------------------------------------------------------------
186  *                                     Setter
187  * ---------------------------------------------------------------------------------
188  */
189
190 // #################################################################################
191 void AliHLTJETReader::SetInputEvent(const TObject* esd, const TObject* aod, const TObject* mc) {
192   // see header file for class documentation
193
194   //  Needs "useMC" flag for running in analysis task only
195
196   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
197
198   // -- Fill ESD
199   if ( esd && !readerHeader->GetUseMC() )
200     fESD = dynamic_cast<AliESDEvent*> (const_cast<TObject*>(esd));
201
202   // -- Fill AOD
203   else if ( aod && !readerHeader->GetUseMC() )
204     fAOD = dynamic_cast<AliAODEvent*> (const_cast<TObject*>(aod));
205
206   // -- Fill MC
207   else if ( mc && readerHeader->GetUseMC() ) {
208
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
213     else 
214       fHLTMC = dynamic_cast<AliHLTMCEvent*> (const_cast<TObject*>(mc));
215   }
216   
217   return;
218 }
219
220 /*
221  * ---------------------------------------------------------------------------------
222  *                            Fastjet Reader functionality
223  * ---------------------------------------------------------------------------------
224  */
225 #ifdef HAVE_FASTJET
226 // #################################################################################
227 Bool_t AliHLTJETReader::FillVector() {
228   // see header file for class documentation
229
230   Bool_t bResult = kFALSE;
231
232   if ( fESD )
233     bResult = FillVectorESD();
234   else if ( fMC )
235     bResult = FillVectorMC();
236   else if ( fHLTMC )
237     bResult = FillVectorHLTMC();
238   else if ( fAOD )
239     bResult = FillVectorAOD();
240   
241   return bResult;
242 }
243
244 // #################################################################################
245 Bool_t AliHLTJETReader::FillVectorMC() {
246   // see header file for class documentation
247
248   Bool_t bResult = kTRUE; 
249
250   if ( ! fMC ) {
251     HLTError( "No MC Event present!" );
252     return kFALSE;
253   }
254   
255   // -- Reset Event
256   ResetEvent();
257    
258   Int_t nTracks = 0;
259
260   AliStack* stack = fMC->Stack();
261
262   for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && bResult; iterStack++) {
263
264     TParticle *particle = stack->Particle(iterStack);
265     if ( !particle) {
266       HLTError( "Error reading particle %i out of %i", iterStack, stack->GetNtrack() );
267       bResult = kFALSE;
268       continue;
269     }
270
271     // ------------------------------
272     // -- Basic cuts on MC particle      --> To be done better XXX
273     // ------------------------------
274
275     // -- primary
276     if ( !(stack->IsPhysicalPrimary(iterStack)) )
277       continue;
278     
279     // -- final state
280     if ( particle->GetNDaughters() != 0 )
281       continue;
282
283     // -- particle in DB
284     TParticlePDG * particlePDG = particle->GetPDG();
285     if ( ! particlePDG ) {
286       particlePDG = TDatabasePDG::Instance()->GetParticle( particle->GetPdgCode() );
287
288       if ( ! particlePDG ) {
289         HLTError("Particle %i not in PDG database", particle->GetPdgCode() );
290         bResult = kFALSE;
291         continue;
292       }
293     }
294
295     // ------------------------
296     // -- Standard track cuts
297     // ------------------------
298
299     // -- Apply track cuts     
300     if ( ! fTrackCuts->IsSelected(particle) )
301       continue;
302     
303     // -- Create PseudoJet object      
304     fastjet::PseudoJet part( particle->Px(), particle->Py(), 
305                              particle->Pz(), particle->Energy() ); 
306     
307     // -- label the particle into Fastjet algortihm
308     part.set_user_index( iterStack ); 
309
310 #ifdef HAVE_FASTJET
311     // -- Add to input_particles vector  
312     fMomentumVector->push_back(part);  
313 #endif
314
315     nTracks++;    
316   } // for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && !bResult; iterStack++) {
317
318   HLTDebug(" Number of selected tracks %d", nTracks);
319
320   return kTRUE;
321 }
322
323 // #################################################################################
324 Bool_t AliHLTJETReader::FillVectorHLTMC() {
325   // see header file for class documentation
326
327   if ( ! fHLTMC ) {
328     HLTError( "No HLT MC Event present!" );
329     return kFALSE;
330   }
331      
332   // -- Reset Event
333   ResetEvent();
334
335   Int_t nTracks = 0;
336
337   TParticle* particle = NULL;
338
339   // -- Loop over particles
340   // ------------------------
341   while ( (particle = fHLTMC->NextParticle() ) ) {
342     
343     // -- Apply track cuts 
344     if ( ! fTrackCuts->IsSelected(particle) )
345       continue;
346     
347     // -- Create PseudoJet object      
348     fastjet::PseudoJet part( particle->Px(), particle->Py(), 
349                              particle->Pz(), particle->Energy() ); 
350     
351     // -- label the particle into Fastjet algortihm
352     part.set_user_index( fHLTMC->GetIndex() ); 
353
354 #ifdef HAVE_FASTJET
355     // -- Add to input_particles vector  
356     fMomentumVector->push_back(part);  
357 #endif
358
359     nTracks++;    
360     
361   } // while ( (particle = fHLTMC->NextParticle() ) ) {
362   
363   HLTInfo(" Number of selected tracks %d", nTracks);
364
365   return kTRUE;
366 }
367
368 // #################################################################################
369 Bool_t AliHLTJETReader::FillVectorESD() {
370   // see header file for class documentation
371
372   Bool_t bResult = kTRUE;
373
374   if ( ! fESD ) {
375     HLTError( "No ESD Event present!" );
376     return kFALSE;
377   }
378
379   // -- Reset Event
380   ResetEvent();
381
382   Int_t nTracks = 0;
383
384   // -- Loop over particles
385   // ------------------------
386   for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && bResult; iter++ ) {
387
388     AliESDtrack* esdTrack = fESD->GetTrack(iter);
389     if ( ! esdTrack ) {
390       HLTError("Could not read ESD track %d from %d", iter, fESD->GetNumberOfTracks() );
391       bResult = kFALSE;
392       continue;
393     }
394
395     // -- Apply track cuts 
396     if ( ! fTrackCuts->IsSelected(esdTrack) )
397       continue;
398
399     // -- Create PseudoJet object      
400     fastjet::PseudoJet part( esdTrack->Px(), esdTrack->Py(), 
401                              esdTrack->Pz(), esdTrack->E() ); 
402     
403     // -- label the particle into Fastjet algortihm
404     part.set_user_index( iter ); 
405
406 #ifdef HAVE_FASTJET
407     // -- Add to input_particles vector  
408     fMomentumVector->push_back(part);  
409 #endif
410
411     nTracks++; 
412
413   } // for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !iResult; iter++ ) {
414   
415   HLTInfo(" Number of selected tracks %d", nTracks);
416
417   return bResult;
418 }
419
420 // #################################################################################
421 Bool_t AliHLTJETReader::FillVectorAOD() {
422   // see header file for class documentation
423
424   return kFALSE;
425 }
426 #endif
427
428 /*
429  * ---------------------------------------------------------------------------------
430  *                               Grid functionality
431  * ---------------------------------------------------------------------------------
432  */
433
434 // #################################################################################
435 Bool_t AliHLTJETReader::FillGrid() {
436   // see header file for class documentation
437
438   Bool_t bResult = kFALSE;
439
440   if ( fESD )
441     bResult = FillGridESD();
442   else if ( fMC )
443     bResult = FillGridMC();
444   else if ( fHLTMC )
445     bResult = FillGridHLTMC();
446   else if ( fAOD )
447     bResult = FillGridAOD();
448
449   return bResult;
450 }
451
452 // #################################################################################
453 Bool_t AliHLTJETReader::FillGridMC() {
454   // see header file for class documentation
455
456   Bool_t bResult = kTRUE; 
457
458   if ( ! fMC ) {
459     HLTError( "No MC Event present!" );
460     return kFALSE;
461   }
462
463   // -- Reset Event
464   ResetEvent();
465
466   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
467
468   Int_t nTracks = 0;
469
470   AliStack* stack = fMC->Stack();
471
472   for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && bResult; iterStack++) {
473
474     TParticle *particle = stack->Particle(iterStack);
475     if ( !particle) {
476       HLTError( "Error reading particle %i out of %i", iterStack, stack->GetNtrack() );
477       bResult = kFALSE;
478       continue;
479     }
480
481     // ------------------------------
482     // -- Basic cuts on MC particle      --> To be done better XXX
483     // ------------------------------
484
485     // -- primary
486     if ( !(stack->IsPhysicalPrimary(iterStack)) )
487       continue;
488
489     // -- final state
490     if ( particle->GetNDaughters() != 0 )
491       continue;
492
493     // -- particle in DB
494     TParticlePDG * particlePDG = particle->GetPDG();
495     if ( ! particlePDG ) {
496       particlePDG = TDatabasePDG::Instance()->GetParticle( particle->GetPdgCode() );
497
498       if ( ! particlePDG ) {
499         HLTError("Particle %i not in PDG database", particle->GetPdgCode() );
500         bResult = kFALSE;
501         continue;
502       }
503     }
504
505     // ------------------------
506     // -- Standard track cuts
507     // ------------------------
508
509     // -- Apply track cuts     
510     if ( ! fTrackCuts->IsSelected(particle) )
511       continue;
512
513     const Float_t aEtaPhi[]  = { particle->Eta(), particle->Phi(), particle->Pt() }; 
514           Int_t   aGridIdx[] = { -1, -1, -1, -1, -1 };
515
516     fGrid->FillTrack(particle, aEtaPhi, aGridIdx);
517           
518     nTracks++;    
519
520     // -- Apply seed cuts 
521     if ( ! fSeedCuts->IsSelected(particle) )
522       continue;  
523
524     // -- Add Seed
525     AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx), readerHeader->GetConeRadius());
526     
527   } // for (Int_t iterStack = 0; iterStack < stack->GetNtrack() && !bResult; iterStack++) {
528
529   HLTInfo(" Number of selected tracks %d", nTracks);
530   HLTInfo(" Number of seeds %d", fNJetCandidates);
531
532   return kTRUE;
533 }
534
535 // #################################################################################
536 Bool_t AliHLTJETReader::FillGridHLTMC() {
537   // see header file for class documentation
538
539   if ( ! fHLTMC ) {
540     HLTError( "No HLT MC Event present!" );
541     return kFALSE;
542   }
543
544   // -- Reset Event
545   ResetEvent();
546
547   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
548
549   Int_t nTracks = 0;
550   TParticle* particle = NULL;
551
552   // -- Loop over particles
553   // ------------------------
554   while ( ( particle = fHLTMC->NextParticle() ) ) {
555     
556     //    HLTError("=== nTracks %d ===",nTracks);
557
558     // -- Apply track cuts 
559     if ( ! fTrackCuts->IsSelected(particle) )
560       continue;
561
562     const Float_t aEtaPhi[]  = { particle->Eta(), particle->Phi(), particle->Pt() }; 
563           Int_t   aGridIdx[] = { -1, -1, -1, -1, -1 };
564
565     // -- Fill grid
566     fGrid->FillTrack(particle, aEtaPhi, aGridIdx);
567
568     nTracks++;    
569
570     // -- Apply seed cuts 
571     if ( ! fSeedCuts->IsSelected(particle) )
572       continue;  
573
574     // -- Add Seed
575     AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
576             readerHeader->GetConeRadius());
577   
578   } // while ( (particle = fHLTMC->NextParticle() ) ) {
579   
580   HLTInfo(" Number of selected tracks %d", nTracks);
581   HLTInfo(" Number of seeds %d", fNJetCandidates);
582
583   return kTRUE;
584 }
585
586 // #################################################################################
587 Bool_t AliHLTJETReader::FillGridESD() {
588   // see header file for class documentation
589
590   Bool_t bResult = kTRUE;
591
592   if ( ! fESD ) {
593     HLTError( "No ESD Event present!" );
594     return kFALSE;
595   }
596
597   // -- Reset Event
598   ResetEvent();
599
600   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
601
602   Int_t nTracks = 0;
603
604   // -- Loop over particles
605   // ------------------------
606   for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !bResult; iter++ ) {
607
608     AliESDtrack* esdTrack = fESD->GetTrack(iter);
609     if ( ! esdTrack ) {
610       HLTError("Could not read ESD track %d from %d", iter, fESD->GetNumberOfTracks() );
611       bResult = kFALSE;
612       continue;
613     }
614
615     // -- Apply track cuts 
616     if ( ! fTrackCuts->IsSelected(esdTrack) )
617       continue;
618
619     const Float_t aEtaPhi[]  = { esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt() }; 
620           Int_t   aGridIdx[] = { -1, -1, -1, -1, -1 };
621
622     // -- Fill grid
623     fGrid->FillTrack(esdTrack, aEtaPhi, aGridIdx);
624
625     nTracks++;    
626
627     // -- Apply seed cuts 
628     if ( ! fSeedCuts->IsSelected(esdTrack) )
629       continue;
630
631     // -- Add Seed
632     AddSeed(aEtaPhi, const_cast<const Int_t*> (aGridIdx),
633             readerHeader->GetConeRadius());
634   
635   } // for ( Int_t iter = 0; iter < fESD->GetNumberOfTracks() && !iResult; iter++ ) {
636   
637   HLTInfo(" Number of selected tracks %d", nTracks);
638   HLTInfo(" Number of seeds %d", fNJetCandidates);
639
640   return bResult;
641 }
642
643 // #################################################################################
644 Bool_t AliHLTJETReader::FillGridAOD() {
645   // see header file for class documentation
646
647   return kFALSE;
648 }
649
650 /*
651  * ---------------------------------------------------------------------------------
652  *                                     Seeds
653  * ---------------------------------------------------------------------------------
654  */
655
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
660
661   Bool_t useWholeCell = kTRUE;
662
663   if ( GetReaderHeader()->GetJetAlgorithm() == AliHLTJETBase::kFFSCRadiusCell )
664     useWholeCell = kFALSE ;
665
666   // -- Add track / particle
667   new( (*fJetCandidates) [fNJetCandidates] ) AliHLTJETConeJetCandidate( aEtaPhi, 
668                                                                         aGridIdx,
669                                                                         coneRadius,
670                                                                         useWholeCell );
671   fNJetCandidates++;
672
673   HLTDebug("Added Seed Pt=%f, Eta=%f, Phi=%f", aEtaPhi[kIdxPt], 
674            aEtaPhi[kIdxEta], aEtaPhi[kIdxPhi] );
675
676   return;
677 }
678
679 /*
680  * ---------------------------------------------------------------------------------
681  *                         Initialize - private
682  * ---------------------------------------------------------------------------------
683  */
684
685 // #################################################################################
686 Int_t AliHLTJETReader::InitializeFFSC() {
687   // see header file for class documentation
688
689   Int_t iResult = 0;
690   AliHLTJETReaderHeader* readerHeader = GetReaderHeader();
691
692   // -- Initialize grid
693   // --------------------
694   if ( fGrid )
695     delete fGrid;
696   
697   if ( ! (fGrid = new AliHLTJETConeGrid()) ) {
698     HLTError("Error instanciating grid.");
699     iResult = -EINPROGRESS;
700   }
701   
702   if ( ! iResult ) {
703     fGrid->SetEtaRange(   readerHeader->GetFiducialEtaMin(),
704                           readerHeader->GetFiducialEtaMax(),
705                           readerHeader->GetGridEtaRange() );
706     
707     fGrid->SetPhiRange(   readerHeader->GetFiducialPhiMin(),
708                           readerHeader->GetFiducialPhiMax(),
709                           readerHeader->GetGridPhiRange() );
710     
711     fGrid->SetBinning(    readerHeader->GetGridEtaBinning(),
712                           readerHeader->GetGridEtaBinning() );
713     
714     fGrid->SetConeRadius( readerHeader->GetConeRadius() );
715     
716     iResult = fGrid->Initialize();
717   }
718   
719   // -- Initialize jet candidates
720   // ------------------------------
721   if ( ! iResult ) {
722     fJetCandidates = new TClonesArray("AliHLTJETConeJetCandidate", 30);
723     if ( ! fJetCandidates) {
724       HLTError("Error instanciating jet candidates.");
725       iResult = -EINPROGRESS;
726     }
727   }
728
729   // -- Get ptr to cuts from reader
730   // --------------------------------
731  
732   // -- Seed cuts
733   if ( ! iResult ) {
734     fSeedCuts = readerHeader->GetSeedCuts();
735     if ( ! fSeedCuts ) {
736       HLTError("Error getting ptr to seed cuts.");
737       iResult = -EINPROGRESS;
738     }
739   }
740  
741   return iResult;
742 }
743
744 #ifdef HAVE_FASTJET
745 // #################################################################################
746 Int_t AliHLTJETReader::InitializeFastjet() {
747   // see header file for class documentation
748
749   Int_t iResult = 0;
750
751   // -- Initialize Vector
752   // ----------------------
753   if ( fMomentumVector )
754     delete fMomentumVector;
755   
756   if ( ! (fMomentumVector = new vector<fastjet::PseudoJet>) ) {
757     HLTError("Error instanciating momentum vector.");
758     iResult = -EINPROGRESS;
759   }
760  
761   return iResult;
762 }
763 #endif