STEER/STEERBase/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 <TTreeCache.h>
  39 #include <TFile.h>
  40 #include <TList.h>
  41 #include <TParticle.h>
  42 #include <TString.h>
  43 #include <TClonesArray.h>
  44 #include <TDirectoryFile.h>
  45
  46 ClassImp(AliMCEventHandler)
  47
  48 AliMCEventHandler::AliMCEventHandler() :
  49     AliInputEventHandler(),
  50     fMCEvent(new AliMCEvent()),
  51     fFileE(0),
  52     fFileK(0),
  53     fFileTR(0),
  54     fTreeE(0),
  55     fTreeK(0),
  56     fTreeTR(0),
  57     fDirK(0),
  58     fDirTR(0),
  59     fParticleSelected(0),
  60     fLabelMap(0),
  61     fNEvent(-1),
  62     fEvent(-1),
  63     fPathName(new TString("./")),
  64     fkExtension(""),
  65     fFileNumber(0),
  66     fEventsPerFile(0),
  67     fReadTR(kTRUE),
  68     fInitOk(kFALSE),
  69     fSubsidiaryHandlers(0),
  70     fEventsInContainer(0),
  71     fPreReadMode(kNoPreRead),
  72     fCacheSize(0),
  73     fCacheTK(0),
  74     fCacheTR(0)
  75 {
  76   //
  77   // Default constructor
  78   //
  79   // Be sure to add all particles to the PDG database
  80   AliPDG::AddParticlesToPdgDataBase();
  81 }
  82
  83 AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
  84     AliInputEventHandler(name, title),
  85     fMCEvent(new AliMCEvent()),
  86     fFileE(0),
  87     fFileK(0),
  88     fFileTR(0),
  89     fTreeE(0),
  90     fTreeK(0),
  91     fTreeTR(0),
  92     fDirK(0),
  93     fDirTR(0),
  94     fParticleSelected(0),
  95     fLabelMap(0),
  96     fNEvent(-1),
  97     fEvent(-1),
  98     fPathName(new TString("./")),
  99     fkExtension(""),
 100     fFileNumber(0),
 101     fEventsPerFile(0),
 102     fReadTR(kTRUE),
 103     fInitOk(kFALSE),
 104     fSubsidiaryHandlers(0),
 105     fEventsInContainer(0),
 106     fPreReadMode(kNoPreRead),
 107     fCacheSize(0),
 108     fCacheTK(0),
 109     fCacheTR(0)
 110 {
 111   //
 112   // Constructor
 113   //
 114   // Be sure to add all particles to the PDG database
 115   AliPDG::AddParticlesToPdgDataBase();
 116 }
 117 AliMCEventHandler::~AliMCEventHandler()
 118 {
 119     // Destructor
 120   delete fPathName;
 121     delete fMCEvent;
 122     delete fFileE;
 123     delete fFileK;
 124     delete fFileTR;
 125     delete fCacheTK;
 126     delete fCacheTR;
 127 }
 128
 129 Bool_t AliMCEventHandler::Init(Option_t* opt)
 130 {
 131     // Initialize input
 132     //
 133     if (!(strcmp(opt, "proof")) || !(strcmp(opt, "local"))) return kTRUE;
 134     //
 135     fFileE = TFile::Open(Form("%sgalice.root", fPathName->Data()));
 136     if (!fFileE) {
 137         AliError(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));
 138         fInitOk = kFALSE;
 139         return kFALSE;
 140     }
 141
 142     //
 143     // Tree E
 144     fFileE->GetObject("TE", fTreeE);
 145     // Connect Tree E to the MCEvent
 146     fMCEvent->ConnectTreeE(fTreeE);
 147     fNEvent = fTreeE->GetEntries();
 148     //
 149     // Tree K
 150     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fkExtension));
 151     if (!fFileK) {
 152         AliError(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName->Data()));
 153         fInitOk = kFALSE;
 154         return kTRUE;
 155     }
 156
 157     fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
 158     //
 159     // Tree TR
 160     if (fReadTR) {
 161         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fkExtension));
 162         if (!fFileTR) {
 163             AliError(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName->Data()));
 164             fInitOk = kFALSE;
 165             return kTRUE;
 166         }
 167     }
 168     //
 169     // Reset the event number
 170     fEvent      = -1;
 171     fFileNumber =  0;
 172     AliInfo(Form("Number of events in this directory %5d \n", fNEvent));
 173     fInitOk = kTRUE;
 174
 175
 176     if (fSubsidiaryHandlers) {
 177         TIter next(fSubsidiaryHandlers);
 178         AliMCEventHandler *handler;
 179         while((handler = (AliMCEventHandler*)next())) {
 180             handler->Init(opt);
 181             handler->SetNumberOfEventsInContainer(fNEvent);
 182         }
 183     }
 184
 185     return kTRUE;
 186 }
 187
 188 Bool_t AliMCEventHandler::LoadEvent(Int_t iev)
 189 {
 190     // Load the event number iev
 191     //
 192     // Calculate the file number
 193   if (!fInitOk) return kFALSE;
 194
 195   Int_t inew  = iev / fEventsPerFile;
 196   Bool_t firsttree = (fTreeK==0) ? kTRUE : kFALSE;
 197 //  Bool_t newtree = firsttree;
 198   if (inew != fFileNumber) {
 199 //    newtree = kTRUE;
 200     fFileNumber = inew;
 201     if (!OpenFile(fFileNumber)){
 202       return kFALSE;
 203     }
 204   }
 205   // Folder name
 206   char folder[20];
 207   snprintf(folder, 20, "Event%d", iev);
 208   // TreeE
 209   fTreeE->GetEntry(iev);
 210   // Tree K
 211   fFileK->GetObject(folder, fDirK);
 212   if (!fDirK) {
 213     AliWarning(Form("AliMCEventHandler: Event #%5d - Cannot get kinematics\n", iev));
 214     return kFALSE;
 215   }
 216
 217   fDirK ->GetObject("TreeK", fTreeK);
 218   if (!fTreeK) {
 219     AliError(Form("AliMCEventHandler: Event #%5d - Cannot get TreeK\n",iev));
 220     return kFALSE;
 221   }
 222   // Connect TreeK to MCEvent
 223   fMCEvent->ConnectTreeK(fTreeK);
 224   //Tree TR
 225   if (fFileTR) {
 226     // Check which format has been read
 227     fFileTR->GetObject(folder, fDirTR);
 228     if (!fDirTR) {
 229       AliError(Form("AliMCEventHandler: Event #%5d - Cannot get track references\n",iev));
 230       return kFALSE;
 231     }
 232
 233     fDirTR->GetObject("TreeTR", fTreeTR);
 234     //
 235     if (!fTreeTR) {
 236       AliError(Form("AliMCEventHandler: Event #%5d - Cannot get TreeTR\n",iev));
 237       return kFALSE;
 238     }
 239     // Connect TR to MCEvent
 240     fMCEvent->ConnectTreeTR(fTreeTR);
 241   }
 242   // Now setup the caches if not yet done
 243   if (fCacheSize) {
 244     fTreeK->SetCacheSize(fCacheSize);
 245     fCacheTK = (TTreeCache*) fFileK->GetCacheRead(fTreeK);
 246     TTreeCache::SetLearnEntries(1);
 247     fTreeK->AddBranchToCache("*",kTRUE);
 248     if (firsttree) Info("LoadEvent","Read cache enabled %lld bytes for TreeK",fCacheSize);
 249     if (fTreeTR) {
 250       fTreeTR->SetCacheSize(fCacheSize);
 251       fCacheTR = (TTreeCache*) fFileTR->GetCacheRead(fTreeTR);
 252       TTreeCache::SetLearnEntries(1);
 253       fTreeTR->AddBranchToCache("*",kTRUE);
 254       if (firsttree) Info("LoadEvent","Read cache enabled %lld bytes for TreeTR",fCacheSize);
 255     }
 256 //    } else {
 257       // We need to reuse the previous caches and every new event is a new tree
 258 //      if (fCacheTK) {
 259 //         fCacheTK->ResetCache();
 260 //         if (fFileK) fFileK->SetCacheRead(fCacheTK, fTreeK);
 261 //         fCacheTK->UpdateBranches(fTreeK);
 262 //      }
 263 //      if (fCacheTR) {
 264 //         fCacheTR->ResetCache();
 265 //         if (fFileTR) fFileTR->SetCacheRead(fCacheTR, fTreeTR);
 266 //         fCacheTR->UpdateBranches(fTreeTR);
 267 //      }
 268   }
 269   return kTRUE;
 270 }
 271
 272 Bool_t AliMCEventHandler::OpenFile(Int_t i)
 273 {
 274     // Open file i
 275     if (i > 0) {
 276         fkExtension = Form("%d", i);
 277     } else {
 278         fkExtension = "";
 279     }
 280
 281     if (fFileK && fCacheTK) fFileK->SetCacheRead(0, fTreeK);
 282     delete fFileK;
 283     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fkExtension));
 284     if (!fFileK) {
 285         AliError(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fkExtension, fPathName->Data()));
 286         fInitOk = kFALSE;
 287         return kFALSE;
 288     }
 289
 290     if (fReadTR) {
 291       if (fFileTR && fCacheTR) fFileTR->SetCacheRead(0, fTreeTR);
 292       delete fFileTR;
 293       fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fkExtension));
 294       if (!fFileTR) {
 295         AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fkExtension, fPathName->Data()));
 296         fInitOk = kFALSE;
 297         return kFALSE;
 298       }
 299     }
 300
 301     fInitOk = kTRUE;
 302
 303     return kTRUE;
 304 }
 305
 306 Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
 307 {
 308     // Begin event
 309     fParticleSelected.Delete();
 310     fLabelMap.Delete();
 311     // Read the next event
 312
 313     if (fEventsInContainer != 0) {
 314         entry = (Long64_t) ( entry * Float_t(fNEvent) / Float_t (fEventsInContainer));
 315     }
 316
 317
 318     if (entry == -1) {
 319         fEvent++;
 320         entry = fEvent;
 321     } else {
 322         fEvent = entry;
 323     }
 324
 325     if (entry >= fNEvent) {
 326         AliWarning(Form("AliMCEventHandler: Event number out of range %5lld %5d\n", entry, fNEvent));
 327         return kFALSE;
 328     }
 329
 330     Bool_t result = LoadEvent(entry);
 331
 332     if (fSubsidiaryHandlers) {
 333         TIter next(fSubsidiaryHandlers);
 334         AliMCEventHandler *handler;
 335         while((handler = (AliMCEventHandler*)next())) {
 336             handler->BeginEvent(entry);
 337         }
 338         next.Reset();
 339         while((handler = (AliMCEventHandler*)next())) {
 340             fMCEvent->AddSubsidiaryEvent(handler->MCEvent());
 341         }
 342         fMCEvent->InitEvent();
 343     }
 344
 345     if (fPreReadMode == kLmPreRead) {
 346         fMCEvent->PreReadAll();
 347     }
 348
 349     return result;
 350
 351 }
 352
 353 void AliMCEventHandler::SelectParticle(Int_t i){
 354   // taking the absolute values here, need to take care
 355   // of negative daughter and mother
 356   // IDs when setting!
 357     if (TMath::Abs(i) >= AliMCEvent::BgLabelOffset()) i =  fMCEvent->BgLabelToIndex(TMath::Abs(i));
 358     if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
 359 }
 360
 361 Bool_t AliMCEventHandler::IsParticleSelected(Int_t i)  {
 362   // taking the absolute values here, need to take
 363   // care with negative daughter and mother
 364   // IDs when setting!
 365   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
 366 }
 367
 368
 369 void AliMCEventHandler::CreateLabelMap(){
 370
 371   //
 372   // this should be called once all selections where done
 373   //
 374
 375   fLabelMap.Delete();
 376   if(!fMCEvent){
 377     fParticleSelected.Delete();
 378     return;
 379   }
 380
 381   VerifySelectedParticles();
 382
 383   Int_t iNew = 0;
 384   for(int i = 0;i < fMCEvent->GetNumberOfTracks();++i){
 385     if(IsParticleSelected(i)){
 386       fLabelMap.Add(i,iNew);
 387       iNew++;
 388     }
 389   }
 390 }
 391
 392 Int_t AliMCEventHandler::GetNewLabel(Int_t i) {
 393   // Gets the label from the new created Map
 394   // Call CreatLabelMap before
 395   // otherwise only 0 returned
 396   return fLabelMap.GetValue(TMath::Abs(i));
 397 }
 398
 399 void  AliMCEventHandler::VerifySelectedParticles(){
 400
 401   //
 402   // Make sure that each particle has at least it's predecessors
 403   // selected so that we have the complete ancestry tree
 404   // Private, should be only called by CreateLabelMap
 405
 406   if(!fMCEvent){
 407       fParticleSelected.Delete();
 408       return;
 409   }
 410
 411   Int_t nprim = fMCEvent->GetNumberOfPrimaries();
 412
 413   for(int i = 0;i < fMCEvent->GetNumberOfTracks(); ++i){
 414       if(i < nprim){
 415           SelectParticle(i);// take all primaries
 416           continue;
 417       }
 418
 419       if(!IsParticleSelected(i))continue;
 420
 421       AliMCParticle* mcpart = (AliMCParticle*) fMCEvent->GetTrack(i);
 422       Int_t imo = mcpart->GetMother();
 423       while((imo >= nprim)&&!IsParticleSelected(imo)){
 424           // Mother not yet selected
 425           SelectParticle(imo);
 426           mcpart = (AliMCParticle*) fMCEvent->GetTrack(imo);
 427           imo = mcpart->GetMother();
 428       }
 429     // after last step we may have an unselected primary
 430     // mother
 431     if(imo>=0){
 432       if(!IsParticleSelected(imo))
 433         SelectParticle(imo);
 434     }
 435   }// loop over all tracks
 436 }
 437
 438 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 439 {
 440     // Retrieve entry i
 441     if (!fInitOk) {
 442         return 0;
 443     } else {
 444         return (fMCEvent->GetParticleAndTR(i, particle, trefs));
 445     }
 446 }
 447
 448 void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
 449 {
 450     // Retrieve entry i and draw momentum vector and hits
 451     fMCEvent->DrawCheck(i, search);
 452 }
 453
 454 Bool_t AliMCEventHandler::Notify(const char *path)
 455 {
 456   // Notify about directory change
 457   // The directory is taken from the 'path' argument
 458   // Reconnect trees
 459     TString fileName(path);
 460     if(fileName.Contains("AliESDs.root")){
 461         fileName.ReplaceAll("AliESDs.root", "");
 462     }
 463     else if(fileName.Contains("AliESDs_wSDD.root")){
 464         fileName.ReplaceAll("AliESDs_wSDD.root", "");
 465     }
 466     else if(fileName.Contains("AliESDs_nob.root")){
 467         fileName.ReplaceAll("AliESDs_nob.root", "");
 468     }
 469     else if(fileName.Contains("AliAOD.root")){
 470         fileName.ReplaceAll("AliAOD.root", "");
 471     }
 472     else if(fileName.Contains("galice.root")){
 473         // for running with galice and kinematics alone...
 474         fileName.ReplaceAll("galice.root", "");
 475     }
 476     else if (fileName.BeginsWith("root:")) {
 477       fileName.Append("?ZIP=");
 478     }
 479
 480     *fPathName = fileName;
 481     AliInfo(Form("Path: -%s-\n", fPathName->Data()));
 482
 483     ResetIO();
 484     InitIO("");
 485
 486 // Handle subsidiary handlers
 487     if (fSubsidiaryHandlers) {
 488         TIter next(fSubsidiaryHandlers);
 489         AliMCEventHandler *handler;
 490         while((handler = (AliMCEventHandler*) next())) {
 491             TString* spath = handler->GetInputPath();
 492             if (spath->Contains("merged")) {
 493                 if (! fPathName->IsNull()) {
 494                     handler->Notify(Form("%s/../.", fPathName->Data()));
 495                 } else {
 496                     handler->Notify("../");
 497                 }
 498             }
 499         }
 500     }
 501
 502     return kTRUE;
 503 }
 504
 505 void AliMCEventHandler::ResetIO()
 506 {
 507 //  Clear header and stack
 508
 509     if (fInitOk) fMCEvent->Clean();
 510
 511 // Delete Tree E
 512     delete fTreeE; fTreeE = 0;
 513
 514 // Reset files
 515     if (fFileE)  {delete fFileE;  fFileE  = 0;}
 516     if (fFileK)  {delete fFileK;  fFileK  = 0;}
 517     if (fFileTR) {delete fFileTR; fFileTR = 0;}
 518     fkExtension="";
 519     fInitOk = kFALSE;
 520
 521     if (fSubsidiaryHandlers) {
 522         TIter next(fSubsidiaryHandlers);
 523         AliMCEventHandler *handler;
 524         while((handler = (AliMCEventHandler*)next())) {
 525             handler->ResetIO();
 526         }
 527     }
 528
 529 }
 530
 531
 532 Bool_t AliMCEventHandler::FinishEvent()
 533 {
 534     // Clean-up after each event
 535    if (fFileK && fCacheTK) {
 536       fTreeK->SetCacheSize(0);
 537       fCacheTK = 0;
 538       fFileK->SetCacheRead(0, fTreeK);
 539    }
 540    if (fFileTR && fCacheTR) {
 541       fTreeTR->SetCacheSize(0);
 542       fCacheTR = 0;
 543       fFileTR->SetCacheRead(0, fTreeTR);
 544    }
 545    delete fDirTR;  fDirTR = 0;
 546    delete fDirK;   fDirK  = 0;
 547     if (fInitOk) fMCEvent->FinishEvent();
 548
 549     if (fSubsidiaryHandlers) {
 550         TIter next(fSubsidiaryHandlers);
 551         AliMCEventHandler *handler;
 552         while((handler = (AliMCEventHandler*)next())) {
 553             handler->FinishEvent();
 554         }
 555     }
 556
 557     return kTRUE;
 558 }
 559
 560 Bool_t AliMCEventHandler::Terminate()
 561 {
 562     // Dummy
 563     return kTRUE;
 564 }
 565
 566 Bool_t AliMCEventHandler::TerminateIO()
 567 {
 568     // Dummy
 569     return kTRUE;
 570 }
 571
 572
 573 void AliMCEventHandler::SetInputPath(const char* fname)
 574 {
 575     // Set the input path name
 576     delete fPathName;
 577     fPathName = new TString(fname);
 578 }
 579
 580 void AliMCEventHandler::AddSubsidiaryHandler(AliMCEventHandler* handler)
 581 {
 582     // Add a subsidiary handler. For example for background events
 583
 584     if (!fSubsidiaryHandlers) fSubsidiaryHandlers = new TList();
 585     fSubsidiaryHandlers->Add(handler);
 586 }