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