Split: HLT/JET -> HLTANALYSIS/JET
[u/mrichter/AliRoot.git] / HLTANALYSIS / 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