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 kTRUE;
 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 kTRUE;
 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
 184     fDirK ->GetObject("TreeK", fTreeK);
 185     // Connect TreeK to MCEvent
 186     fMCEvent->ConnectTreeK(fTreeK);
 187     //Tree TR
 188     if (fFileTR) {
 189         // Check which format has been read
 190         fFileTR->GetObject(folder, fDirTR);
 191         fDirTR->GetObject("TreeTR", fTreeTR);
 192         //
 193         // Connect TR to MCEvent
 194         fMCEvent->ConnectTreeTR(fTreeTR);
 195     }
 196     //
 197     return kTRUE;
 198 }
 199
 200 Bool_t AliMCEventHandler::OpenFile(Int_t i)
 201 {
 202     // Open file i
 203     if (i > 0) {
 204         fExtension = Form("%d", i);
 205     } else {
 206         fExtension = "";
 207     }
 208
 209
 210     delete fFileK;
 211     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
 212     if (!fFileK) {
 213         AliError(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 214         fInitOk = kFALSE;
 215         return kFALSE;
 216     }
 217
 218     if (fReadTR) {
 219         delete fFileTR;
 220         fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
 221         if (!fFileTR) {
 222             AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
 223             fInitOk = kFALSE;
 224             return kFALSE;
 225         }
 226     }
 227
 228     fInitOk = kTRUE;
 229     return kTRUE;
 230 }
 231
 232 Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
 233 {
 234     fParticleSelected.Delete();
 235     fLabelMap.Delete();
 236     // Read the next event
 237     if (entry == -1) {
 238         fEvent++;
 239         entry = fEvent;
 240     } else {
 241         fEvent = entry;
 242     }
 243
 244     if (entry >= fNEvent) {
 245         AliWarning(Form("AliMCEventHandler: Event number out of range %5d %5d\n", entry,fNEvent));
 246         return kFALSE;
 247     }
 248     return GetEvent(entry);
 249 }
 250
 251 void AliMCEventHandler::SelectParticle(Int_t i){
 252   // taking the absolute values here, need to take care
 253   // of negative daughter and mother
 254   // IDs when setting!
 255   if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
 256 }
 257
 258 Bool_t AliMCEventHandler::IsParticleSelected(Int_t i)  {
 259   // taking the absolute values here, need to take
 260   // care with negative daughter and mother
 261   // IDs when setting!
 262   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
 263 }
 264
 265
 266 void AliMCEventHandler::CreateLabelMap(){
 267
 268   //
 269   // this should be called once all selections where done
 270   //
 271
 272   fLabelMap.Delete();
 273   if(!fMCEvent){
 274     fParticleSelected.Delete();
 275     return;
 276   }
 277
 278   VerifySelectedParticles();
 279   AliStack *pStack = fMCEvent->Stack();
 280
 281   Int_t iNew = 0;
 282   for(int i = 0;i < pStack->GetNtrack();++i){
 283     if(IsParticleSelected(i)){
 284       fLabelMap.Add(i,iNew);
 285       iNew++;
 286     }
 287   }
 288 }
 289
 290 Int_t AliMCEventHandler::GetNewLabel(Int_t i) {
 291   // Gets the labe from the new created Map
 292   // Call CreatLabelMap before
 293   // otherwise only 0 returned
 294   return fLabelMap.GetValue(TMath::Abs(i));
 295 }
 296
 297 void  AliMCEventHandler::VerifySelectedParticles(){
 298
 299   //
 300   // Make sure that each particle has at least it's predecessors
 301   // selected so that we have the complete ancestry tree
 302   // Private, should be only called by CreateLabelMap
 303
 304   if(!fMCEvent){
 305     fParticleSelected.Delete();
 306     return;
 307   }
 308   AliStack *pStack = fMCEvent->Stack();
 309
 310   Int_t nprim = pStack->GetNprimary();
 311
 312   for(int i = 0;i < pStack->GetNtrack();++i){
 313     if(i<nprim){
 314       SelectParticle(i);// take all primaries
 315       continue;
 316     }
 317     if(!IsParticleSelected(i))continue;
 318     TParticle *part = pStack->Particle(i);
 319     Int_t imo = part->GetFirstMother();
 320     while((imo >= nprim)&&!IsParticleSelected(imo)){
 321       // Mother not yet selected
 322       SelectParticle(imo);
 323       TParticle *mother = pStack->Particle(imo);
 324       imo = mother->GetFirstMother();
 325     }
 326     // after last step we may have a unselected primary
 327     // mother
 328     if(imo>=0){
 329       if(!IsParticleSelected(imo))
 330         SelectParticle(imo);
 331     }
 332   }// loop over all tracks
 333 }
 334
 335 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 336 {
 337     // Retrieve entry i
 338     if (!fInitOk) {
 339         return 0;
 340     } else {
 341         return (fMCEvent->GetParticleAndTR(i, particle, trefs));
 342     }
 343 }
 344
 345 void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
 346 {
 347     // Retrieve entry i and draw momentum vector and hits
 348     fMCEvent->DrawCheck(i, search);
 349 }
 350
 351 Bool_t AliMCEventHandler::Notify(const char *path)
 352 {
 353   // Notify about directory change
 354   // The directory is taken from the 'path' argument
 355   // Reconnect trees
 356     TString fileName(path);
 357     if(fileName.Contains("AliESDs.root")){
 358         fileName.ReplaceAll("AliESDs.root", "");
 359     }
 360     else if(fileName.Contains("AliAOD.root")){
 361         fileName.ReplaceAll("AliAOD.root", "");
 362     }
 363     else if(fileName.Contains("galice.root")){
 364         // for running with galice and kinematics alone...
 365         fileName.ReplaceAll("galice.root", "");
 366     }
 367     else if (fileName.BeginsWith("root:")) {
 368       fileName.Append("?ZIP=");
 369     }
 370     *fPathName = fileName;
 371     AliInfo(Form("Notify() Path: %s\n", fPathName->Data()));
 372
 373     ResetIO();
 374     InitIO("");
 375
 376     return kTRUE;
 377 }
 378
 379 void AliMCEventHandler::ResetIO()
 380 {
 381 //  Clear header and stack
 382
 383     if (fInitOk) fMCEvent->Clean();
 384
 385 // Delete Tree E
 386     delete fTreeE; fTreeE = 0;
 387
 388 // Reset files
 389     if (fFileE)  {delete fFileE;  fFileE  = 0;}
 390     if (fFileK)  {delete fFileK;  fFileK  = 0;}
 391     if (fFileTR) {delete fFileTR; fFileTR = 0;}
 392     fExtension="";
 393     fInitOk = kFALSE;
 394 }
 395
 396
 397 Bool_t AliMCEventHandler::FinishEvent()
 398 {
 399     // Clean-up after each event
 400     delete fDirTR;  fDirTR = 0;
 401     delete fDirK;   fDirK  = 0;
 402     if (fInitOk) fMCEvent->FinishEvent();
 403     return kTRUE;
 404 }
 405
 406 Bool_t AliMCEventHandler::Terminate()
 407 {
 408     // Dummy
 409     return kTRUE;
 410 }
 411
 412 Bool_t AliMCEventHandler::TerminateIO()
 413 {
 414     // Dummy
 415     return kTRUE;
 416 }
 417
 418
 419 void AliMCEventHandler::SetInputPath(const char* fname)
 420 {
 421     // Set the input path name
 422     delete fPathName;
 423     fPathName = new TString(fname);
 424 }