Line data Source code
1 : //-*- Mode: C++ -*-
2 : // $Id: AliHLTMCEvent.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 AliHLTMCEvent.cxx
20 : @author Jochen Thaeder
21 : @date
22 : @brief Container class for an AliMCEvent
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 "TList.h"
32 : #include "TParticlePDG.h"
33 : #include "TDatabasePDG.h"
34 :
35 : #include "AliGenCocktailEventHeader.h"
36 :
37 : #include "AliHLTMCEvent.h"
38 :
39 : using namespace std;
40 :
41 : /** ROOT macro for the implementation of ROOT specific class methods */
42 8 : ClassImp(AliHLTMCEvent)
43 :
44 : /*
45 : * ---------------------------------------------------------------------------------
46 : * Constructor / Destructor
47 : * ---------------------------------------------------------------------------------
48 : */
49 :
50 : // #################################################################################
51 0 : AliHLTMCEvent::AliHLTMCEvent() :
52 0 : fCurrentParticleIndex(-1),
53 0 : fNParticles(0),
54 0 : fStack( NULL ),
55 0 : fCurrentGenJetIndex(-1),
56 0 : fNGenJets(0),
57 0 : fGenJets(NULL),
58 0 : fHasParticleCuts(kFALSE) {
59 : // see header file for class documentation
60 : // or
61 : // refer to README to build package
62 : // or
63 : // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
64 0 : }
65 :
66 : // #################################################################################
67 0 : AliHLTMCEvent::AliHLTMCEvent( Bool_t applyParticleCuts ) :
68 0 : fCurrentParticleIndex(-1),
69 0 : fNParticles(0),
70 0 : fStack( NULL ),
71 0 : fCurrentGenJetIndex(-1),
72 0 : fNGenJets(0),
73 0 : fGenJets(NULL),
74 0 : fHasParticleCuts(applyParticleCuts) {
75 : // see header file for class documentation
76 0 : }
77 :
78 : // #################################################################################
79 0 : AliHLTMCEvent::~AliHLTMCEvent() {
80 : // see header file for class documentation
81 :
82 0 : if ( fStack ) {
83 0 : fStack->Delete();
84 0 : delete fStack;
85 : }
86 0 : fStack = NULL;
87 :
88 0 : if ( fGenJets ) {
89 0 : fGenJets->Delete();
90 0 : delete fGenJets;
91 : }
92 0 : fGenJets = NULL;
93 0 : }
94 :
95 : /*
96 : * ---------------------------------------------------------------------------------
97 : * Setter - public
98 : * ---------------------------------------------------------------------------------
99 : */
100 :
101 : // #################################################################################
102 : Int_t AliHLTMCEvent::FillMCEvent( AliStack *stack, AliHeader *header ) {
103 : // see header file for class documentation
104 :
105 : Int_t iResult = 0;
106 :
107 0 : if ( stack ) {
108 0 : if ( (iResult = FillMCTracks(stack)) )
109 0 : HLTError("Error filling particles" );
110 : }
111 : else {
112 0 : HLTError("Error reading stack" );
113 : iResult = -2;
114 : }
115 :
116 0 : if ( header ) {
117 0 : AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*> (header->GenEventHeader());
118 0 : if ( pythiaGenHeader ) {
119 0 : if ( (iResult = FillMCJets(pythiaGenHeader)) )
120 0 : HLTError("Error filling jets" );
121 : }
122 : else {
123 0 : HLTError("Error reading PythiaHeader" );
124 : iResult = -2;
125 : }
126 0 : }
127 : else {
128 0 : HLTError("Error reading header" );
129 : iResult = -2;
130 : }
131 :
132 0 : Compress();
133 :
134 0 : return iResult;;
135 : }
136 :
137 : // #################################################################################
138 : Int_t AliHLTMCEvent::FillMCEvent( AliMCEvent *pMCEvent ) {
139 : // see header file for class documentation
140 :
141 : Int_t iResult = 0;
142 :
143 0 : if ( (iResult = FillMCTracks(pMCEvent->Stack())) )
144 0 : HLTError("Error filling particles" );
145 :
146 0 : if ( (iResult = FillMCJets(GetPythiaEventHeader(pMCEvent))) )
147 0 : HLTError("Error filling jets" );
148 :
149 0 : Compress();
150 :
151 0 : return iResult;
152 : }
153 :
154 : /*
155 : * ---------------------------------------------------------------------------------
156 : * Getter
157 : * ---------------------------------------------------------------------------------
158 : */
159 :
160 : // #################################################################################
161 : TParticle* AliHLTMCEvent::Particle( Int_t iParticle ) {
162 : // see header file for class documentation
163 :
164 0 : if ( iParticle >= fNParticles || !fStack )
165 0 : return NULL;
166 :
167 0 : return (TParticle*) (*fStack)[iParticle];
168 0 : }
169 :
170 : // #################################################################################
171 : TParticle* AliHLTMCEvent::NextParticle() {
172 : // see header file for class documentation
173 :
174 0 : fCurrentParticleIndex++;
175 0 : return Particle( fCurrentParticleIndex );
176 : }
177 :
178 : // #################################################################################
179 : AliAODJet* AliHLTMCEvent::GenJet( Int_t iJet ) const {
180 : // see header file for class documentation
181 :
182 0 : if ( iJet >= fNGenJets || !fGenJets )
183 0 : return NULL;
184 :
185 0 : return reinterpret_cast<AliAODJet*> ((*fGenJets)[iJet]);
186 0 : }
187 :
188 : // #################################################################################
189 : AliAODJet* AliHLTMCEvent:: NextGenJet() {
190 : // see header file for class documentation
191 :
192 0 : fCurrentGenJetIndex++;
193 0 : return GenJet( fCurrentGenJetIndex );
194 : }
195 :
196 : /*
197 : * ---------------------------------------------------------------------------------
198 : * Helper
199 : * ---------------------------------------------------------------------------------
200 : */
201 :
202 : // #################################################################################
203 : void AliHLTMCEvent::Compress() {
204 : // see header file for class documentation
205 :
206 0 : if (fStack)
207 0 : fStack->Compress();
208 0 : if (fGenJets)
209 0 : fGenJets->Compress();
210 0 : }
211 :
212 : // #################################################################################
213 : void AliHLTMCEvent::Reset() {
214 : // see header file for class documentation
215 :
216 0 : fCurrentParticleIndex = -1;
217 0 : fCurrentGenJetIndex = -1;
218 0 : }
219 :
220 : /*
221 : * ---------------------------------------------------------------------------------
222 : * Setter - private
223 : * ---------------------------------------------------------------------------------
224 : */
225 :
226 : // #################################################################################
227 : Int_t AliHLTMCEvent::FillMCTracks( AliStack* stack ) {
228 : // see header file for class documentation
229 :
230 : Int_t iResult = 0;
231 :
232 : // -- Create local stack
233 0 : if ( stack && stack->GetNtrack() > 0 )
234 0 : fStack = new TClonesArray("TParticle", stack->GetNtrack() );
235 : else {
236 0 : HLTError( "Error creating local stack" );
237 : iResult = -EINPROGRESS;
238 : }
239 :
240 : // -- Loop over off-line stack and fill local stack
241 0 : for (Int_t iterStack = 0; !iResult &&iterStack < stack->GetNtrack(); iterStack++) {
242 :
243 0 : TParticle *particle = stack->Particle(iterStack);
244 0 : if ( !particle) {
245 0 : HLTError( "Error reading particle %i out of %i", iterStack, stack->GetNtrack() );
246 : iResult = -EINPROGRESS;
247 0 : continue;
248 : }
249 :
250 : // ----------------
251 : // -- Apply cuts
252 : // ----------------
253 :
254 0 : if ( fHasParticleCuts ) {
255 :
256 : // -- primary
257 0 : if ( !(stack->IsPhysicalPrimary(iterStack)) )
258 0 : continue;
259 :
260 : // -- final state
261 0 : if ( particle->GetNDaughters() != 0 )
262 0 : continue;
263 :
264 : // -- particle in DB
265 0 : TParticlePDG * particlePDG = particle->GetPDG();
266 0 : if ( ! particlePDG ) {
267 0 : particlePDG = TDatabasePDG::Instance()->GetParticle( particle->GetPdgCode() );
268 :
269 0 : if ( ! particlePDG ) {
270 0 : HLTError("Particle %i not in PDG database", particle->GetPdgCode() );
271 : iResult = -EINPROGRESS;
272 0 : continue;
273 : }
274 : }
275 0 : }
276 :
277 : // -- Add particle after cuts
278 0 : AddParticle ( particle );
279 0 : }
280 :
281 0 : return iResult;
282 0 : }
283 :
284 : // #################################################################################
285 : void AliHLTMCEvent::AddParticle( const TParticle* particle ) {
286 : // see header file for class documentation
287 :
288 0 : new( (*fStack) [fNParticles] ) TParticle( *particle );
289 0 : fNParticles++;
290 :
291 0 : return;
292 0 : }
293 :
294 : // #################################################################################
295 : Int_t AliHLTMCEvent::FillMCJets( AliGenPythiaEventHeader* header ) {
296 : // see header file for class documentation
297 :
298 : Int_t iResult = 0;
299 :
300 : // -- Check if Header is present
301 0 : if ( !header ) {
302 0 : HLTError( "Error no PythiaHeader present" );
303 : iResult = -EINPROGRESS;
304 0 : return iResult;
305 : }
306 :
307 : // -- Create jet array
308 0 : if ( header->NTriggerJets() > 0 )
309 0 : fGenJets = new TClonesArray("AliAODJet", header->NTriggerJets());
310 : else
311 0 : return iResult;
312 :
313 : // -- Loop over jets in header and fill local array
314 0 : for (Int_t iterJet = 0; iterJet < header->NTriggerJets() && !iResult; iterJet++) {
315 :
316 : // -- Add jet
317 0 : AddGenJet(header, iterJet);
318 :
319 : } // for (Int_t iterJet = 0; iterJet < header->NTriggerJets() && !iResult; iterJet++) {
320 :
321 : HLTDebug("Pythia Jets found: %d", fNGenJets );
322 :
323 0 : return iResult;
324 0 : }
325 :
326 : /*
327 : * ---------------------------------------------------------------------------------
328 : * Pythia jets - private
329 : * ---------------------------------------------------------------------------------
330 : */
331 :
332 : //##################################################################################
333 : AliGenPythiaEventHeader* AliHLTMCEvent::GetPythiaEventHeader(AliMCEvent *mcEvent) {
334 : // see header file for class documentation
335 :
336 : Int_t iResult = 0;
337 :
338 0 : AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
339 :
340 0 : if ( !pythiaGenHeader ) {
341 :
342 : // -- is cocktail header
343 0 : AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(mcEvent->GenEventHeader());
344 0 : if ( !genCocktailHeader ) {
345 0 : HLTError("Error: Unknown header type (not Pythia or Cocktail)");
346 : iResult = -1;
347 0 : }
348 :
349 0 : if ( !iResult ) {
350 : // -- Get Header
351 0 : TList* headerList = genCocktailHeader->GetHeaders();
352 :
353 0 : for (Int_t iter = 0; iter < headerList->GetEntries(); iter++ ) {
354 0 : pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(iter));
355 0 : if ( pythiaGenHeader )
356 0 : break;
357 : }
358 0 : }
359 0 : }
360 :
361 0 : if ( pythiaGenHeader && !iResult )
362 0 : return pythiaGenHeader;
363 : else {
364 0 : HLTError("PythiaHeader not found!");
365 0 : return NULL;
366 : }
367 0 : }
368 :
369 : //##################################################################################
370 : void AliHLTMCEvent::AddGenJet( AliGenPythiaEventHeader* header, Int_t iterJet ) {
371 : // see header file for class documentation
372 :
373 0 : Float_t pJet[] = {0., 0., 0., 0.}; // jet 4-vector ( x y z )
374 :
375 : // -- Get jet
376 0 : header->TriggerJet(iterJet, pJet);
377 :
378 : // -- Create TLorentzVector
379 0 : TLorentzVector v(pJet);
380 :
381 : // -- Add AliAODJet
382 0 : new( (*fGenJets) [fNGenJets] ) AliAODJet(v);
383 0 : fNGenJets++;
384 :
385 : return;
386 0 : }
|
