HLT/BASE/util/AliHLTMCEvent.cxx

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