STEER/AliMCEventHandler.cxx

   1 /************************************************************************* * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved *
   2  *                                                                        *
   3  * Author: The ALICE Off-line Project.                                    *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 /* $Id$ */
  16 //---------------------------------------------------------------------------------
  17 //                          Class AliMCEventHandler
  18 // This class gives access to MC truth during the analysis.
  19 // Monte Carlo truth is containe in the kinematics tree (produced particles) and
  20 // the tree of reference hits.
  21 //
  22 // Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch
  23 //---------------------------------------------------------------------------------
  24
  25
  26
  27 #include "AliMCEventHandler.h"
  28 #include "AliMCEvent.h"
  29 #include "AliPDG.h"
  30 #include "AliTrackReference.h"
  31 #include "AliHeader.h"
  32 #include "AliStack.h"
  33 #include "AliLog.h"
  34
  35 #include <TTree.h>
  36 #include <TFile.h>
  37 #include <TParticle.h>
  38 #include <TString.h>
  39 #include <TClonesArray.h>
  40 #include <TDirectoryFile.h>
  41
  42 ClassImp(AliMCEventHandler)
  43
  44 AliMCEventHandler::AliMCEventHandler() :
  45     AliVEventHandler(),
  46     fMCEvent(new AliMCEvent()),
  47     fFileE(0),
  48     fFileK(0),
  49     fFileTR(0),
  50     fTreeE(0),
  51     fTreeK(0),
  52     fTreeTR(0),
  53     fDirK(0),
  54     fDirTR(0),
  55     fParticleSelected(0),
  56     fLabelMap(0),
  57     fNEvent(-1),
  58     fEvent(-1),
  59     fPathName(new TString("./")),
  60     fExtension(""),
  61     fFileNumber(0),
  62     fEventsPerFile(0),
  63     fReadTR(kTRUE),
  64     fInitOk(kFALSE)
  65 {
  66   //
  67   // Default constructor
  68   //
  69   // Be sure to add all particles to the PDG database
  70   AliPDG::AddParticlesToPdgDataBase();
  71 }
  72
  73 AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
  74     AliVEventHandler(name, title),
  75     fMCEvent(new AliMCEvent()),
  76     fFileE(0),
  77     fFileK(0),
  78     fFileTR(0),
  79     fTreeE(0),
  80     fTreeK(0),
  81     fTreeTR(0),
  82     fDirK(0),
  83     fDirTR(0),
  84     fParticleSelected(0),
  85     fLabelMap(0),
  86     fNEvent(-1),
  87     fEvent(-1),
  88     fPathName(new TString("./")),
  89     fExtension(""),
  90     fFileNumber(0),
  91     fEventsPerFile(0),
  92     fReadTR(kTRUE),
  93     fInitOk(kFALSE)
  94 {
  95   //
  96   // Constructor
  97   //
  98   // Be sure to add all particles to the PDG database
  99   AliPDG::AddParticlesToPdgDataBase();
 100 }
 101 AliMCEventHandler::~AliMCEventHandler()
 102 {
 103     // Destructor
 104     delete fMCEvent;
 105     delete fFileE;
 106     delete fFileK;
 107     delete fFileTR;
 108 }
 109
 110 Bool_t AliMCEventHandler::Init(Option_t* opt)
 111 {
 112     // Initialize input
 113     //
 114     if (!(strcmp(opt, "proof")) || !(strcmp(opt, "local"))) return kTRUE;
 115     //
 116     fFileE = TFile::Open(Form("%sgalice.root", fPathName->Data()));
 117     if (!fFileE) {
 118         AliError(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));
 119         fInitOk = kFALSE;
 120         return kFALSE;
 121     }
 122
 123     //
 124     // Tree E
 125     fFileE->GetObject("TE", fTreeE);
 126     // Connect Tree E to the MCEvent
 127     fMCEvent->ConnectTreeE(fTreeE);
 128     fNEvent = fTreeE->GetEntries();
 129     //
 130     // Tree K
 131     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
 132     if (!fFileK) {
 133         AliError(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName));
 134         fInitOk = kFALSE;
 135         return kFALSE;
 136     }
 137
 138     fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
 139     //
 140     // Tree TR
 141     if (fReadTR) {
 142         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
 143         if (!fFileTR) {
 144             AliError(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName->Data()));
 145             fInitOk = kFALSE;
 146             return kFALSE;
 147         }
 148     }
 149     //
 150     // Reset the event number
 151     fEvent      = -1;
 152     fFileNumber =  0;
 153     AliInfo(Form("Number of events in this directory %5d \n", fNEvent));
 154     fInitOk = kTRUE;
 155     return kTRUE;
 156 }
 157
 158 Bool_t AliMCEventHandler::GetEvent(Int_t iev)
 159 {
 160     // Load the event number iev
 161     //
 162     // Calculate the file number
 163     if (!fInitOk) return kFALSE;
 164
 165     Int_t inew  = iev / fEventsPerFile;
 166     if (inew != fFileNumber) {
 167         fFileNumber = inew;
 168         if (!OpenFile(fFileNumber)){
 169             return kFALSE;
 170         }
 171     }
 172     // Folder name
 173     char folder[20];
 174     sprintf(folder, "Event%d", iev);
 175     // TreeE
 176     fTreeE->GetEntry(iev);
 177     // Tree K
 178     fFileK->GetObject(folder, fDirK);
 179     if (!fDirK) {
 180         AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
 181         return kFALSE;
 182     }
 183     fDirK ->GetObject("TreeK", fTreeK);
 184     // Connect TreeK to MCEvent
 185     fMCEvent->ConnectTreeK(fTreeK);
 186     //Tree TR
 187     if (fFileTR) {
 188         // Check which format has been read
 189         fFileTR->GetObject(folder, fDirTR);
 190         fDirTR->GetObject("TreeTR", fTreeTR);
 191         //
 192         // Connect TR to MCEvent
 193         fMCEvent->ConnectTreeTR(fTreeTR);
 194     }
 195     //
 196     return kTRUE;
 197 }
 198
 199 Bool_t AliMCEventHandler::OpenFile(Int_t i)
 200 {
 201     // Open file i
 202     if (i > 0) {
 203         fExtension = Form("%d", i);
 204     } else {
 205         fExtension = "";
 206     }
 207
 208
 209     delete fFileK;
 210     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
 211     if (!fFileK) {
 212         AliError(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 213         fInitOk = kFALSE;
 214         return kFALSE;
 215     }
 216
 217     if (fReadTR) {
 218         delete fFileTR;
 219         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
 220         if (!fFileTR) {
 221             AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 222             fInitOk = kFALSE;
 223             return kFALSE;
 224         }
 225     }
 226
 227     fInitOk = kTRUE;
 228     return kTRUE;
 229 }
 230
 231 Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
 232 {
 233     fParticleSelected.Delete();
 234     fLabelMap.Delete();
 235     // Read the next event
 236     if (entry == -1) {
 237         fEvent++;
 238         entry = fEvent;
 239     } else {
 240         fEvent = entry;
 241     }
 242
 243     if (entry >= fNEvent) {
 244         AliWarning(Form("AliMCEventHandler: Event number out of range %5d %5d\n", entry,fNEvent));
 245         return kFALSE;
 246     }
 247     return GetEvent(entry);
 248 }
 249
 250 void AliMCEventHandler::SelectParticle(Int_t i){
 251   // taking the absolute values here, need to take care
 252   // of negative daughter and mother
 253   // IDs when setting!
 254   if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
 255 }
 256
 257 Bool_t AliMCEventHandler::IsParticleSelected(Int_t i)  {
 258   // taking the absolute values here, need to take
 259   // care with negative daughter and mother
 260   // IDs when setting!
 261   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
 262 }
 263
 264
 265 void AliMCEventHandler::CreateLabelMap(){
 266
 267   //
 268   // this should be called once all selections where done
 269   //
 270
 271   fLabelMap.Delete();
 272   if(!fMCEvent){
 273     fParticleSelected.Delete();
 274     return;
 275   }
 276
 277   VerifySelectedParticles();
 278   AliStack *pStack = fMCEvent->Stack();
 279
 280   Int_t iNew = 0;
 281   for(int i = 0;i < pStack->GetNtrack();++i){
 282     if(IsParticleSelected(i)){
 283       fLabelMap.Add(i,iNew);
 284       iNew++;
 285     }
 286   }
 287 }
 288
 289 Int_t AliMCEventHandler::GetNewLabel(Int_t i) {
 290   // Gets the labe from the new created Map
 291   // Call CreatLabelMap before
 292   // otherwise only 0 returned
 293   return fLabelMap.GetValue(TMath::Abs(i));
 294 }
 295
 296 void  AliMCEventHandler::VerifySelectedParticles(){
 297
 298   //
 299   // Make sure that each particle has at least it's predecessors
 300   // selected so that we have the complete ancestry tree
 301   // Private, should be only called by CreateLabelMap
 302
 303   if(!fMCEvent){
 304     fParticleSelected.Delete();
 305     return;
 306   }
 307   AliStack *pStack = fMCEvent->Stack();
 308
 309   Int_t nprim = pStack->GetNprimary();
 310
 311   for(int i = 0;i < pStack->GetNtrack();++i){
 312     if(i<nprim){
 313       SelectParticle(i);// take all primaries
 314       continue;
 315     }
 316     if(!IsParticleSelected(i))continue;
 317     TParticle *part = pStack->Particle(i);
 318     Int_t imo = part->GetFirstMother();
 319     while((imo >= nprim)&&!IsParticleSelected(imo)){
 320       // Mother not yet selected
 321       SelectParticle(imo);
 322       TParticle *mother = pStack->Particle(imo);
 323       imo = mother->GetFirstMother();
 324     }
 325     // after last step we may have a unselected primary
 326     // mother
 327     if(imo>=0){
 328       if(!IsParticleSelected(imo))
 329         SelectParticle(imo);
 330     }
 331   }// loop over all tracks
 332 }
 333
 334 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 335 {
 336     // Retrieve entry i
 337     if (!fInitOk) {
 338         return 0;
 339     } else {
 340         return (fMCEvent->GetParticleAndTR(i, particle, trefs));
 341     }
 342 }
 343
 344 void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
 345 {
 346     // Retrieve entry i and draw momentum vector and hits
 347     fMCEvent->DrawCheck(i, search);
 348 }
 349
 350 Bool_t AliMCEventHandler::Notify(const char *path)
 351 {
 352   // Notify about directory change
 353   // The directory is taken from the 'path' argument
 354   // Reconnect trees
 355     TString fileName(path);
 356     if(fileName.Contains("AliESDs.root")){
 357         fileName.ReplaceAll("AliESDs.root", "");
 358     }
 359     else if(fileName.Contains("AliAOD.root")){
 360         fileName.ReplaceAll("AliAOD.root", "");
 361     }
 362     else if(fileName.Contains("galice.root")){
 363         // for running with galice and kinematics alone...
 364         fileName.ReplaceAll("galice.root", "");
 365     }
 366
 367     *fPathName = fileName;
 368     AliInfo(Form("Notify() Path: %s\n", fPathName->Data()));
 369
 370     ResetIO();
 371     InitIO("");
 372
 373     return kTRUE;
 374 }
 375
 376 void AliMCEventHandler::ResetIO()
 377 {
 378 //  Clear header and stack
 379     if (fInitOk) fMCEvent->Clean();
 380
 381 // Delete Tree E
 382     delete fTreeE; fTreeE = 0;
 383
 384 // Reset files
 385     if (fFileE)  {delete fFileE;  fFileE  = 0;}
 386     if (fFileK)  {delete fFileK;  fFileK  = 0;}
 387     if (fFileTR) {delete fFileTR; fFileTR = 0;}
 388     fExtension="";
 389     fInitOk = kFALSE;
 390 }
 391
 392
 393 Bool_t AliMCEventHandler::FinishEvent()
 394 {
 395     // Clean-up after each event
 396     delete fDirTR;  fDirTR = 0;
 397     delete fDirK;   fDirK  = 0;
 398     if (fInitOk) fMCEvent->FinishEvent();
 399     return kTRUE;
 400 }
 401
 402 Bool_t AliMCEventHandler::Terminate()
 403 {
 404     // Dummy
 405     return kTRUE;
 406 }
 407
 408 Bool_t AliMCEventHandler::TerminateIO()
 409 {
 410     // Dummy
 411     return kTRUE;
 412 }
 413
 414
 415 void AliMCEventHandler::SetInputPath(const char* fname)
 416 {
 417     // Set the input path name
 418     delete fPathName;
 419     fPathName = new TString(fname);
 420 }