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 fMCEvent;
 114     delete fFileE;
 115     delete fFileK;
 116     delete fFileTR;
 117 }
 118
 119 Bool_t AliMCEventHandler::Init(Option_t* opt)
 120 {
 121     // Initialize input
 122     //
 123     if (!(strcmp(opt, "proof")) || !(strcmp(opt, "local"))) return kTRUE;
 124     //
 125     fFileE = TFile::Open(Form("%sgalice.root", fPathName->Data()));
 126     if (!fFileE) {
 127         AliError(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));
 128         fInitOk = kFALSE;
 129         return kFALSE;
 130     }
 131
 132     //
 133     // Tree E
 134     fFileE->GetObject("TE", fTreeE);
 135     // Connect Tree E to the MCEvent
 136     fMCEvent->ConnectTreeE(fTreeE);
 137     fNEvent = fTreeE->GetEntries();
 138     //
 139     // Tree K
 140     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
 141     if (!fFileK) {
 142         AliError(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName));
 143         fInitOk = kFALSE;
 144         return kTRUE;
 145     }
 146
 147     fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
 148     //
 149     // Tree TR
 150     if (fReadTR) {
 151         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
 152         if (!fFileTR) {
 153             AliError(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName->Data()));
 154             fInitOk = kFALSE;
 155             return kTRUE;
 156         }
 157     }
 158     //
 159     // Reset the event number
 160     fEvent      = -1;
 161     fFileNumber =  0;
 162     AliInfo(Form("Number of events in this directory %5d \n", fNEvent));
 163     fInitOk = kTRUE;
 164
 165
 166     if (fSubsidiaryHandlers) {
 167         TIter next(fSubsidiaryHandlers);
 168         AliMCEventHandler *handler;
 169         while((handler = (AliMCEventHandler*)next())) {
 170             handler->Init(opt);
 171             handler->SetNumberOfEventsInContainer(fNEvent);
 172         }
 173     }
 174
 175     return kTRUE;
 176 }
 177
 178 Bool_t AliMCEventHandler::GetEvent(Int_t iev)
 179 {
 180     // Load the event number iev
 181     //
 182     // Calculate the file number
 183   if (!fInitOk) return kFALSE;
 184
 185   Int_t inew  = iev / fEventsPerFile;
 186   if (inew != fFileNumber) {
 187     fFileNumber = inew;
 188     if (!OpenFile(fFileNumber)){
 189       return kFALSE;
 190     }
 191   }
 192   // Folder name
 193   char folder[20];
 194   sprintf(folder, "Event%d", iev);
 195   // TreeE
 196   fTreeE->GetEntry(iev);
 197   // Tree K
 198   fFileK->GetObject(folder, fDirK);
 199   if (!fDirK) {
 200     AliWarning(Form("AliMCEventHandler: Event #%5d - Cannot get kinematics\n", iev));
 201     return kFALSE;
 202   }
 203
 204   fDirK ->GetObject("TreeK", fTreeK);
 205   if (!fTreeK) {
 206     AliError(Form("AliMCEventHandler: Event #%5d - Cannot get TreeK\n",iev));
 207     return kFALSE;
 208   }
 209   // Connect TreeK to MCEvent
 210   fMCEvent->ConnectTreeK(fTreeK);
 211   //Tree TR
 212   if (fFileTR) {
 213     // Check which format has been read
 214     fFileTR->GetObject(folder, fDirTR);
 215     if (!fDirTR) {
 216       AliError(Form("AliMCEventHandler: Event #%5d - Cannot get track references\n",iev));
 217       return kFALSE;
 218     }
 219
 220     fDirTR->GetObject("TreeTR", fTreeTR);
 221     //
 222     if (!fTreeTR) {
 223       AliError(Form("AliMCEventHandler: Event #%5d - Cannot get TreeTR\n",iev));
 224       return kFALSE;
 225     }
 226     // Connect TR to MCEvent
 227     fMCEvent->ConnectTreeTR(fTreeTR);
 228   }
 229
 230   //
 231   return kTRUE;
 232 }
 233
 234 Bool_t AliMCEventHandler::OpenFile(Int_t i)
 235 {
 236     // Open file i
 237     if (i > 0) {
 238         fExtension = Form("%d", i);
 239     } else {
 240         fExtension = "";
 241     }
 242
 243
 244     delete fFileK;
 245     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
 246     if (!fFileK) {
 247         AliError(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 248         fInitOk = kFALSE;
 249         return kFALSE;
 250     }
 251
 252     if (fReadTR) {
 253         delete fFileTR;
 254         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
 255         if (!fFileTR) {
 256             AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 257             fInitOk = kFALSE;
 258             return kFALSE;
 259         }
 260     }
 261
 262     fInitOk = kTRUE;
 263
 264     return kTRUE;
 265 }
 266
 267 Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
 268 {
 269     fParticleSelected.Delete();
 270     fLabelMap.Delete();
 271     // Read the next event
 272
 273     if (fEventsInContainer != 0) {
 274         entry = (Long64_t) ( entry * Float_t(fNEvent) / Float_t (fEventsInContainer));
 275     }
 276
 277
 278     if (entry == -1) {
 279         fEvent++;
 280         entry = fEvent;
 281     } else {
 282         fEvent = entry;
 283     }
 284
 285     if (entry >= fNEvent) {
 286         AliWarning(Form("AliMCEventHandler: Event number out of range %5d %5d\n", entry, fNEvent));
 287         return kFALSE;
 288     }
 289
 290     Bool_t result = GetEvent(entry);
 291
 292     if (fSubsidiaryHandlers) {
 293         TIter next(fSubsidiaryHandlers);
 294         AliMCEventHandler *handler;
 295         while((handler = (AliMCEventHandler*)next())) {
 296             handler->BeginEvent(entry);
 297         }
 298         next.Reset();
 299         while((handler = (AliMCEventHandler*)next())) {
 300             fMCEvent->AddSubsidiaryEvent(handler->MCEvent());
 301         }
 302         fMCEvent->InitEvent();
 303     }
 304
 305     if (fPreReadMode == kLmPreRead) {
 306         fMCEvent->PreReadAll();
 307     }
 308
 309     return result;
 310
 311 }
 312
 313 void AliMCEventHandler::SelectParticle(Int_t i){
 314   // taking the absolute values here, need to take care
 315   // of negative daughter and mother
 316   // IDs when setting!
 317     if (TMath::Abs(i) >= AliMCEvent::BgLabelOffset()) i =  fMCEvent->BgLabelToIndex(TMath::Abs(i));
 318     if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
 319 }
 320
 321 Bool_t AliMCEventHandler::IsParticleSelected(Int_t i)  {
 322   // taking the absolute values here, need to take
 323   // care with negative daughter and mother
 324   // IDs when setting!
 325   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
 326 }
 327
 328
 329 void AliMCEventHandler::CreateLabelMap(){
 330
 331   //
 332   // this should be called once all selections where done
 333   //
 334
 335   fLabelMap.Delete();
 336   if(!fMCEvent){
 337     fParticleSelected.Delete();
 338     return;
 339   }
 340
 341   VerifySelectedParticles();
 342
 343   Int_t iNew = 0;
 344   for(int i = 0;i < fMCEvent->GetNumberOfTracks();++i){
 345     if(IsParticleSelected(i)){
 346       fLabelMap.Add(i,iNew);
 347       iNew++;
 348     }
 349   }
 350 }
 351
 352 Int_t AliMCEventHandler::GetNewLabel(Int_t i) {
 353   // Gets the label from the new created Map
 354   // Call CreatLabelMap before
 355   // otherwise only 0 returned
 356   return fLabelMap.GetValue(TMath::Abs(i));
 357 }
 358
 359 void  AliMCEventHandler::VerifySelectedParticles(){
 360
 361   //
 362   // Make sure that each particle has at least it's predecessors
 363   // selected so that we have the complete ancestry tree
 364   // Private, should be only called by CreateLabelMap
 365
 366   if(!fMCEvent){
 367       fParticleSelected.Delete();
 368       return;
 369   }
 370
 371   Int_t nprim = fMCEvent->GetNumberOfPrimaries();
 372
 373   for(int i = 0;i < fMCEvent->GetNumberOfTracks(); ++i){
 374       if(i < nprim){
 375           SelectParticle(i);// take all primaries
 376           continue;
 377       }
 378
 379       if(!IsParticleSelected(i))continue;
 380
 381       AliMCParticle* mcpart = (AliMCParticle*) fMCEvent->GetTrack(i);
 382       Int_t imo = mcpart->GetMother();
 383       while((imo >= nprim)&&!IsParticleSelected(imo)){
 384           // Mother not yet selected
 385           SelectParticle(imo);
 386           mcpart = (AliMCParticle*) fMCEvent->GetTrack(imo);
 387           imo = mcpart->GetMother();
 388       }
 389     // after last step we may have an unselected primary
 390     // mother
 391     if(imo>=0){
 392       if(!IsParticleSelected(imo))
 393         SelectParticle(imo);
 394     }
 395   }// loop over all tracks
 396 }
 397
 398 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 399 {
 400     // Retrieve entry i
 401     if (!fInitOk) {
 402         return 0;
 403     } else {
 404         return (fMCEvent->GetParticleAndTR(i, particle, trefs));
 405     }
 406 }
 407
 408 void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
 409 {
 410     // Retrieve entry i and draw momentum vector and hits
 411     fMCEvent->DrawCheck(i, search);
 412 }
 413
 414 Bool_t AliMCEventHandler::Notify(const char *path)
 415 {
 416   // Notify about directory change
 417   // The directory is taken from the 'path' argument
 418   // Reconnect trees
 419     TString fileName(path);
 420     if(fileName.Contains("AliESDs.root")){
 421         fileName.ReplaceAll("AliESDs.root", "");
 422     }
 423     else if(fileName.Contains("AliAOD.root")){
 424         fileName.ReplaceAll("AliAOD.root", "");
 425     }
 426     else if(fileName.Contains("galice.root")){
 427         // for running with galice and kinematics alone...
 428         fileName.ReplaceAll("galice.root", "");
 429     }
 430     else if (fileName.BeginsWith("root:")) {
 431       fileName.Append("?ZIP=");
 432     }
 433
 434     *fPathName = fileName;
 435     AliInfo(Form("Path: -%s-\n", fPathName->Data()));
 436
 437     ResetIO();
 438     InitIO("");
 439
 440 // Handle subsidiary handlers
 441     if (fSubsidiaryHandlers) {
 442         TIter next(fSubsidiaryHandlers);
 443         AliMCEventHandler *handler;
 444         while((handler = (AliMCEventHandler*) next())) {
 445             TString* spath = handler->GetInputPath();
 446             if (spath->Contains("merged")) {
 447                 if (! fPathName->IsNull()) {
 448                     handler->Notify(Form("%s/../.", fPathName->Data()));
 449                 } else {
 450                     handler->Notify("../");
 451                 }
 452             }
 453         }
 454     }
 455
 456     return kTRUE;
 457 }
 458
 459 void AliMCEventHandler::ResetIO()
 460 {
 461 //  Clear header and stack
 462
 463     if (fInitOk) fMCEvent->Clean();
 464
 465 // Delete Tree E
 466     delete fTreeE; fTreeE = 0;
 467
 468 // Reset files
 469     if (fFileE)  {delete fFileE;  fFileE  = 0;}
 470     if (fFileK)  {delete fFileK;  fFileK  = 0;}
 471     if (fFileTR) {delete fFileTR; fFileTR = 0;}
 472     fExtension="";
 473     fInitOk = kFALSE;
 474
 475     if (fSubsidiaryHandlers) {
 476         TIter next(fSubsidiaryHandlers);
 477         AliMCEventHandler *handler;
 478         while((handler = (AliMCEventHandler*)next())) {
 479             handler->ResetIO();
 480         }
 481     }
 482
 483 }
 484
 485
 486 Bool_t AliMCEventHandler::FinishEvent()
 487 {
 488     // Clean-up after each event
 489     delete fDirTR;  fDirTR = 0;
 490     delete fDirK;   fDirK  = 0;
 491     if (fInitOk) fMCEvent->FinishEvent();
 492
 493     if (fSubsidiaryHandlers) {
 494         TIter next(fSubsidiaryHandlers);
 495         AliMCEventHandler *handler;
 496         while((handler = (AliMCEventHandler*)next())) {
 497             handler->FinishEvent();
 498         }
 499     }
 500
 501     return kTRUE;
 502 }
 503
 504 Bool_t AliMCEventHandler::Terminate()
 505 {
 506     // Dummy
 507     return kTRUE;
 508 }
 509
 510 Bool_t AliMCEventHandler::TerminateIO()
 511 {
 512     // Dummy
 513     return kTRUE;
 514 }
 515
 516
 517 void AliMCEventHandler::SetInputPath(const char* fname)
 518 {
 519     // Set the input path name
 520     delete fPathName;
 521     fPathName = new TString(fname);
 522 }
 523
 524 void AliMCEventHandler::AddSubsidiaryHandler(AliMCEventHandler* handler)
 525 {
 526     // Add a subsidiary handler. For example for background events
 527
 528     if (!fSubsidiaryHandlers) fSubsidiaryHandlers = new TList();
 529     fSubsidiaryHandlers->Add(handler);
 530 }