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