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