STEER/AliMCEventHandler.cxx

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