STAGE_POOL and STAGE_HOST defined
[u/mrichter/AliRoot.git] / STEER / AliRun.cxx
CommitLineData
99d554c8 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/*
17$Log$
6781433e 18Revision 1.57 2001/02/23 17:40:23 buncic
19All trees needed for simulation created in RunMC(). TreeR and its branches
20are now created in new RunReco() method.
21
d47c658f 22Revision 1.56 2001/02/14 15:45:20 hristov
23Algorithmic way of getting entry index in fParticleMap. Protection of fParticleFileMap (I.Hrivnacova)
24
a41d61a9 25Revision 1.55 2001/02/12 15:52:54 buncic
26Removed OpenBaseFile().
27
39de14fb 28Revision 1.54 2001/02/07 10:39:05 hristov
29Remove default value for argument
30
1d994b80 31Revision 1.53 2001/02/06 11:02:26 hristov
32New SetTrack interface added, added check for unfilled particles in FinishEvent (I.Hrivnacova)
33
89bbad6f 34Revision 1.52 2001/02/05 16:22:25 buncic
35Added TreeS to GetEvent().
36
82711e7a 37Revision 1.51 2001/02/02 15:16:20 morsch
38SetHighWaterMark method added to mark last particle in event.
39
4d69d91e 40Revision 1.50 2001/01/27 10:32:00 hristov
41Leave the loop when primaries are filled (I.Hrivnacova)
42
cb94ac2a 43Revision 1.49 2001/01/26 19:58:48 hristov
44Major upgrade of AliRoot code
45
2ab0c725 46Revision 1.48 2001/01/17 10:50:50 hristov
47Corrections to destructors
48
e460afec 49Revision 1.47 2000/12/18 10:44:01 morsch
50Possibility to set field map by passing pointer to objet of type AliMagF via
51SetField().
52Example:
53gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
54
d8408e76 55Revision 1.46 2000/12/14 19:29:27 fca
56galice.cuts was not read any more
57
327136d2 58Revision 1.45 2000/11/30 07:12:49 alibrary
59Introducing new Rndm and QA classes
60
65fb704d 61Revision 1.44 2000/10/26 13:58:59 morsch
62Add possibility to choose the lego generator (of type AliGeneratorLego or derived) when running
63RunLego(). Default is the base class AliGeneratorLego.
64
0a520a66 65Revision 1.43 2000/10/09 09:43:17 fca
66Special remapping of hits for TPC and TRD. End-of-primary action introduced
67
24de2263 68Revision 1.42 2000/10/02 21:28:14 fca
69Removal of useless dependecies via forward declarations
70
94de3818 71Revision 1.41 2000/07/13 16:19:09 fca
72Mainly coding conventions + some small bug fixes
73
ef42d733 74Revision 1.40 2000/07/12 08:56:25 fca
75Coding convention correction and warning removal
76
8918e700 77Revision 1.39 2000/07/11 18:24:59 fca
78Coding convention corrections + few minor bug fixes
79
aee8290b 80Revision 1.38 2000/06/20 13:05:45 fca
81Writing down the TREE headers before job starts
82
51e0e89d 83Revision 1.37 2000/06/09 20:05:11 morsch
84Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat
85
f1b9d7c3 86Revision 1.36 2000/06/08 14:03:58 hristov
87Only one initializer for a default argument
88
d33c0226 89Revision 1.35 2000/06/07 10:13:14 hristov
90Delete only existent objects.
91
d2ecd553 92Revision 1.34 2000/05/18 10:45:38 fca
93Delete Particle Factory properly
94
c222d2b0 95Revision 1.33 2000/05/16 13:10:40 fca
96New method IsNewTrack and fix for a problem in Father-Daughter relations
97
a01a8b12 98Revision 1.32 2000/04/27 10:38:21 fca
99Correct termination of Lego Run and introduce Lego getter in AliRun
100
838edcaf 101Revision 1.31 2000/04/26 10:17:32 fca
102Changes in Lego for G4 compatibility
103
dffd31ef 104Revision 1.30 2000/04/18 19:11:40 fca
105Introduce variable Config.C function signature
106
45189757 107Revision 1.29 2000/04/07 11:12:34 fca
108G4 compatibility changes
109
875c717b 110Revision 1.28 2000/04/05 06:51:06 fca
111Workaround for an HP compiler problem
112
5eb58812 113Revision 1.27 2000/03/22 18:08:07 fca
114Rationalisation of the virtual MC interfaces
115
80762cb1 116Revision 1.26 2000/03/22 13:42:26 fca
117SetGenerator does not replace an existing generator, ResetGenerator does
118
ee1dd322 119Revision 1.25 2000/02/23 16:25:22 fca
120AliVMC and AliGeant3 classes introduced
121ReadEuclid moved from AliRun to AliModule
122
b13db077 123Revision 1.24 2000/01/19 17:17:20 fca
124Introducing a list of lists of hits -- more hits allowed for detector now
125
1cedd08a 126Revision 1.23 1999/12/03 11:14:31 fca
127Fixing previous wrong checking
128
00719c1b 129Revision 1.21 1999/11/25 10:40:08 fca
130Fixing daughters information also in primary tracks
131
ae23d366 132Revision 1.20 1999/10/04 18:08:49 fca
133Adding protection against inconsistent Euclid files
134
3fcc96a1 135Revision 1.19 1999/09/29 07:50:40 fca
136Introduction of the Copyright and cvs Log
137
99d554c8 138*/
139
fe4da5cc 140///////////////////////////////////////////////////////////////////////////////
141// //
142// Control class for Alice C++ //
143// Only one single instance of this class exists. //
144// The object is created in main program aliroot //
145// and is pointed by the global gAlice. //
146// //
8494b010 147// -Supports the list of all Alice Detectors (fModules). //
fe4da5cc 148// -Supports the list of particles (fParticles). //
149// -Supports the Trees. //
150// -Supports the geometry. //
151// -Supports the event display. //
152//Begin_Html
153/*
1439f98e 154<img src="picts/AliRunClass.gif">
fe4da5cc 155*/
156//End_Html
157//Begin_Html
158/*
1439f98e 159<img src="picts/alirun.gif">
fe4da5cc 160*/
161//End_Html
162// //
163///////////////////////////////////////////////////////////////////////////////
164
94de3818 165#include <stdlib.h>
166#include <stdio.h>
167#include <string.h>
2ab0c725 168#include <iostream.h>
94de3818 169
fe4da5cc 170#include <TFile.h>
171#include <TRandom.h>
172#include <TBRIK.h>
173#include <TNode.h>
fe4da5cc 174#include <TCint.h>
175#include <TSystem.h>
a8f1fb7c 176#include <TObjectTable.h>
94de3818 177#include <TTree.h>
178#include <TGeometry.h>
179#include <TROOT.h>
180#include "TBrowser.h"
fe4da5cc 181
1578254f 182#include "TParticle.h"
fe4da5cc 183#include "AliRun.h"
fe4da5cc 184#include "AliDisplay.h"
875c717b 185#include "AliMC.h"
dffd31ef 186#include "AliLego.h"
aee8290b 187#include "AliMagFC.h"
188#include "AliMagFCM.h"
189#include "AliMagFDM.h"
94de3818 190#include "AliHit.h"
65fb704d 191#include "TRandom3.h"
192#include "AliMCQA.h"
193#include "AliGenerator.h"
194#include "AliLegoGenerator.h"
94de3818 195
196#include "AliDetector.h"
fe4da5cc 197
fe4da5cc 198AliRun *gAlice;
199
aee8290b 200static AliHeader *gAliHeader;
fe4da5cc 201
fe4da5cc 202ClassImp(AliRun)
203
204//_____________________________________________________________________________
205AliRun::AliRun()
206{
207 //
208 // Default constructor for AliRun
209 //
aee8290b 210 gAliHeader=&fHeader;
fe4da5cc 211 fRun = 0;
212 fEvent = 0;
213 fCurrent = -1;
8494b010 214 fModules = 0;
fe4da5cc 215 fGenerator = 0;
216 fTreeD = 0;
217 fTreeK = 0;
218 fTreeH = 0;
219 fTreeE = 0;
220 fTreeR = 0;
2ab0c725 221 fTreeS = 0;
fe4da5cc 222 fParticles = 0;
223 fGeometry = 0;
224 fDisplay = 0;
225 fField = 0;
875c717b 226 fMC = 0;
fe4da5cc 227 fNdets = 0;
228 fImedia = 0;
229 fTrRmax = 1.e10;
230 fTrZmax = 1.e10;
fe4da5cc 231 fInitDone = kFALSE;
232 fLego = 0;
1578254f 233 fPDGDB = 0; //Particle factory object!
1cedd08a 234 fHitLists = 0;
45189757 235 fConfigFunction = "\0";
65fb704d 236 fRandom = 0;
237 fMCQA = 0;
238 fTransParName = "\0";
39de14fb 239 fBaseFileName = ".\0";
2ab0c725 240 fParticleBuffer = 0;
241 fParticleMap = new TObjArray(10000);
fe4da5cc 242}
243
244//_____________________________________________________________________________
245AliRun::AliRun(const char *name, const char *title)
246 : TNamed(name,title)
247{
248 //
249 // Constructor for the main processor.
250 // Creates the geometry
251 // Creates the list of Detectors.
252 // Creates the list of particles.
253 //
254 Int_t i;
255
256 gAlice = this;
257 fTreeD = 0;
258 fTreeK = 0;
259 fTreeH = 0;
260 fTreeE = 0;
261 fTreeR = 0;
2ab0c725 262 fTreeS = 0;
fe4da5cc 263 fTrRmax = 1.e10;
264 fTrZmax = 1.e10;
1141f8e4 265 fGenerator = 0;
fe4da5cc 266 fInitDone = kFALSE;
267 fLego = 0;
268 fField = 0;
45189757 269 fConfigFunction = "Config();";
65fb704d 270
271 // Set random number generator
272 gRandom = fRandom = new TRandom3();
2ab0c725 273
274 if (gSystem->Getenv("CONFIG_SEED")) {
275 gRandom->SetSeed((UInt_t)atoi(gSystem->Getenv("CONFIG_SEED")));
276 }
fe4da5cc 277
278 gROOT->GetListOfBrowsables()->Add(this,name);
279 //
280 // create the support list for the various Detectors
8494b010 281 fModules = new TObjArray(77);
fe4da5cc 282 //
283 // Create the TNode geometry for the event display
284
285 BuildSimpleGeometry();
286
287
288 fNtrack=0;
289 fHgwmk=0;
290 fCurrent=-1;
aee8290b 291 gAliHeader=&fHeader;
fe4da5cc 292 fRun = 0;
293 fEvent = 0;
294 //
295 // Create the particle stack
2ab0c725 296 fParticles = new TClonesArray("TParticle",1000);
fe4da5cc 297
298 fDisplay = 0;
299 //
300 // Create default mag field
301 SetField();
302 //
875c717b 303 fMC = gMC;
fe4da5cc 304 //
305 // Prepare the tracking medium lists
306 fImedia = new TArrayI(1000);
307 for(i=0;i<1000;i++) (*fImedia)[i]=-99;
1578254f 308 //
309 // Make particles
310 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
1cedd08a 311 //
312 // Create HitLists list
313 fHitLists = new TList();
65fb704d 314 //
315 SetTransPar();
39de14fb 316 fBaseFileName = ".\0";
2ab0c725 317 fParticleBuffer = 0;
318 fParticleMap = new TObjArray(10000);
fe4da5cc 319}
320
aee8290b 321
322//_____________________________________________________________________________
fe4da5cc 323AliRun::~AliRun()
324{
325 //
2ab0c725 326 // Default AliRun destructor
fe4da5cc 327 //
fe4da5cc 328 delete fImedia;
329 delete fField;
875c717b 330 delete fMC;
fe4da5cc 331 delete fGeometry;
332 delete fDisplay;
333 delete fGenerator;
334 delete fLego;
335 delete fTreeD;
336 delete fTreeK;
337 delete fTreeH;
338 delete fTreeE;
339 delete fTreeR;
2ab0c725 340 delete fTreeS;
8494b010 341 if (fModules) {
342 fModules->Delete();
343 delete fModules;
fe4da5cc 344 }
345 if (fParticles) {
346 fParticles->Delete();
347 delete fParticles;
348 }
1cedd08a 349 delete fHitLists;
c222d2b0 350 delete fPDGDB;
e460afec 351 delete fMCQA;
fe4da5cc 352}
353
354//_____________________________________________________________________________
355void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
356{
357 //
358 // Add a hit to detector id
359 //
8494b010 360 TObjArray &dets = *fModules;
361 if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
fe4da5cc 362}
363
364//_____________________________________________________________________________
365void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
366{
367 //
368 // Add digit to detector id
369 //
8494b010 370 TObjArray &dets = *fModules;
371 if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
fe4da5cc 372}
373
374//_____________________________________________________________________________
375void AliRun::Browse(TBrowser *b)
376{
377 //
378 // Called when the item "Run" is clicked on the left pane
379 // of the Root browser.
380 // It displays the Root Trees and all detectors.
381 //
382 if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
383 if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
384 if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
385 if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
386 if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
2ab0c725 387 if (fTreeS) b->Add(fTreeS,fTreeS->GetName());
fe4da5cc 388
8494b010 389 TIter next(fModules);
390 AliModule *detector;
391 while((detector = (AliModule*)next())) {
fe4da5cc 392 b->Add(detector,detector->GetName());
393 }
65fb704d 394 b->Add(fMCQA,"AliMCQA");
fe4da5cc 395}
396
397//_____________________________________________________________________________
398void AliRun::Build()
399{
400 //
401 // Initialize Alice geometry
402 // Dummy routine
403 //
404}
405
406//_____________________________________________________________________________
407void AliRun::BuildSimpleGeometry()
408{
409 //
410 // Create a simple TNode geometry used by Root display engine
411 //
412 // Initialise geometry
413 //
414 fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
415 new TMaterial("void","Vacuum",0,0,0); //Everything is void
416 TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
417 brik->SetVisibility(0);
418 new TNode("alice","alice","S_alice");
419}
420
421//_____________________________________________________________________________
422void AliRun::CleanDetectors()
423{
424 //
425 // Clean Detectors at the end of event
426 //
8494b010 427 TIter next(fModules);
428 AliModule *detector;
429 while((detector = (AliModule*)next())) {
fe4da5cc 430 detector->FinishEvent();
431 }
432}
433
434//_____________________________________________________________________________
435void AliRun::CleanParents()
436{
437 //
438 // Clean Particles stack.
1578254f 439 // Set parent/daughter relations
fe4da5cc 440 //
2ab0c725 441 TObjArray &particles = *fParticleMap;
1578254f 442 TParticle *part;
fe4da5cc 443 int i;
2ab0c725 444 for(i=0; i<fHgwmk+1; i++) {
445 part = (TParticle *)particles.At(i);
446 if(part) if(!part->TestBit(kDaughtersBit)) {
1578254f 447 part->SetFirstDaughter(-1);
448 part->SetLastDaughter(-1);
fe4da5cc 449 }
450 }
451}
452
453//_____________________________________________________________________________
454Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
455{
456 //
457 // Return the distance from the mouse to the AliRun object
458 // Dummy routine
459 //
460 return 9999;
461}
462
463//_____________________________________________________________________________
94de3818 464void AliRun::DumpPart (Int_t i) const
fe4da5cc 465{
466 //
467 // Dumps particle i in the stack
468 //
2ab0c725 469 ((TParticle*) (*fParticleMap)[i])->Print();
fe4da5cc 470}
471
472//_____________________________________________________________________________
94de3818 473void AliRun::DumpPStack () const
fe4da5cc 474{
475 //
476 // Dumps the particle stack
477 //
2ab0c725 478 TObjArray &particles = *fParticleMap;
fe4da5cc 479 printf(
480 "\n\n=======================================================================\n");
481 for (Int_t i=0;i<fNtrack;i++)
482 {
1578254f 483 printf("-> %d ",i); ((TParticle*) particles[i])->Print();
fe4da5cc 484 printf("--------------------------------------------------------------\n");
485 }
486 printf(
487 "\n=======================================================================\n\n");
488}
489
d8408e76 490void AliRun::SetField(AliMagF* magField)
491{
492 // Set Magnetic Field Map
493 fField = magField;
494 fField->ReadField();
495}
496
fe4da5cc 497//_____________________________________________________________________________
498void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
499 Float_t maxField, char* filename)
500{
501 //
502 // Set magnetic field parameters
503 // type Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
504 // version Magnetic field map version (only 1 active now)
505 // scale Scale factor for the magnetic field
506 // maxField Maximum value for the magnetic field
507
508 //
509 // --- Sanity check on mag field flags
fe4da5cc 510 if(fField) delete fField;
511 if(version==1) {
d8408e76 512 fField = new AliMagFC("Map1"," ",type,scale,maxField);
f1b9d7c3 513 } else if(version<=2) {
d8408e76 514 fField = new AliMagFCM("Map2-3",filename,type,scale,maxField);
fe4da5cc 515 fField->ReadField();
f1b9d7c3 516 } else if(version==3) {
d8408e76 517 fField = new AliMagFDM("Map4",filename,type,scale,maxField);
f1b9d7c3 518 fField->ReadField();
fe4da5cc 519 } else {
23370b7a 520 Warning("SetField","Invalid map %d\n",version);
fe4da5cc 521 }
522}
fe4da5cc 523
524//_____________________________________________________________________________
65fb704d 525void AliRun::PreTrack()
526{
527 TObjArray &dets = *fModules;
528 AliModule *module;
529
530 for(Int_t i=0; i<=fNdets; i++)
531 if((module = (AliModule*)dets[i]))
532 module->PreTrack();
533
534 fMCQA->PreTrack();
535}
536
537//_____________________________________________________________________________
538void AliRun::PostTrack()
539{
540 TObjArray &dets = *fModules;
541 AliModule *module;
542
543 for(Int_t i=0; i<=fNdets; i++)
544 if((module = (AliModule*)dets[i]))
545 module->PostTrack();
546}
547
548//_____________________________________________________________________________
fe4da5cc 549void AliRun::FinishPrimary()
550{
551 //
552 // Called at the end of each primary track
553 //
554
6c9704e6 555 // static Int_t count=0;
556 // const Int_t times=10;
fe4da5cc 557 // This primary is finished, purify stack
80762cb1 558 PurifyKine();
fe4da5cc 559
24de2263 560 TIter next(fModules);
561 AliModule *detector;
562 while((detector = (AliModule*)next())) {
563 detector->FinishPrimary();
564 }
565
fe4da5cc 566 // Write out hits if any
567 if (gAlice->TreeH()) {
568 gAlice->TreeH()->Fill();
569 }
570
571 // Reset Hits info
572 gAlice->ResetHits();
a8f1fb7c 573
574 //
575 // if(++count%times==1) gObjectTable->Print();
fe4da5cc 576}
577
578//_____________________________________________________________________________
579void AliRun::FinishEvent()
580{
581 //
582 // Called at the end of the event.
583 //
7fb01480 584
dffd31ef 585 //
586 if(fLego) fLego->FinishEvent();
587
fe4da5cc 588 //Update the energy deposit tables
589 Int_t i;
875c717b 590 for(i=0;i<fEventEnergy.GetSize();i++) {
591 fSummEnergy[i]+=fEventEnergy[i];
592 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
fe4da5cc 593 }
875c717b 594 fEventEnergy.Reset();
fe4da5cc 595
596 // Clean detector information
597 CleanDetectors();
598
599 // Write out the kinematics
600 if (fTreeK) {
601 CleanParents();
a41d61a9 602 if(fTreeK->GetEntries() ==0) {
603 // set the fParticleFileMap size for the first time
6781433e 604 if (fHgwmk+1 > fParticleFileMap.GetSize())
a41d61a9 605 fParticleFileMap.Set(fHgwmk+1);
606 }
607 // fTreeK->Fill();
89bbad6f 608 Bool_t allFilled = kFALSE;
2ab0c725 609 TObject *part;
610 for(i=0; i<fHgwmk+1; ++i) if((part=fParticleMap->At(i))) {
611 fParticleBuffer = (TParticle*) part;
612 fParticleFileMap[i]= (Int_t) fTreeK->GetEntries();
613 fTreeK->Fill();
89bbad6f 614 (*fParticleMap)[i]=0;
615
616 // When all primaries were filled no particle!=0
617 // should be left => to be removed later.
618 if (allFilled) printf("Why != 0 part # %d?\n",i);
619 }
620 else {
621 // // printf("Why = 0 part # %d?\n",i); => We know.
622 // break;
623 // we don't break now in order to be sure there is no
624 // particle !=0 left.
625 // To be removed later and replaced with break.
626 if(!allFilled) allFilled = kTRUE;
627 }
fe4da5cc 628 }
629
a41d61a9 630 // Set number of tracks to event header
631 fHeader.SetNtrack(fNtrack);
632
fe4da5cc 633 // Write out the digits
634 if (fTreeD) {
635 fTreeD->Fill();
636 ResetDigits();
637 }
2ab0c725 638
639 if (fTreeS) {
640 fTreeS->Fill();
641 ResetSDigits();
642 }
643
fe4da5cc 644 // Write out reconstructed clusters
645 if (fTreeR) {
646 fTreeR->Fill();
647 }
648
649 // Write out the event Header information
650 if (fTreeE) fTreeE->Fill();
651
652 // Reset stack info
653 ResetStack();
654
655 // Write Tree headers
51e0e89d 656 if (fTreeK) fTreeK->Write(0,TObject::kOverwrite);
51e0e89d 657 if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
51e0e89d 658 if (fTreeD) fTreeD->Write(0,TObject::kOverwrite);
51e0e89d 659 if (fTreeR) fTreeR->Write(0,TObject::kOverwrite);
2ab0c725 660 if (fTreeS) fTreeS->Write(0,TObject::kOverwrite);
661
875c717b 662 ++fEvent;
fe4da5cc 663}
664
665//_____________________________________________________________________________
666void AliRun::FinishRun()
667{
668 //
669 // Called at the end of the run.
670 //
671
dffd31ef 672 //
673 if(fLego) fLego->FinishRun();
674
fe4da5cc 675 // Clean detector information
8494b010 676 TIter next(fModules);
677 AliModule *detector;
678 while((detector = (AliModule*)next())) {
fe4da5cc 679 detector->FinishRun();
680 }
681
682 //Output energy summary tables
683 EnergySummary();
2ab0c725 684
685 TFile *file = fTreeE->GetCurrentFile();
686
aee8290b 687 file->cd();
2ab0c725 688
51e0e89d 689 fTreeE->Write(0,TObject::kOverwrite);
fe4da5cc 690
2ab0c725 691 // Write AliRun info and all detectors parameters
d47c658f 692 Write(0,TObject::kOverwrite);
2ab0c725 693
fe4da5cc 694 // Clean tree information
d2ecd553 695 if (fTreeK) {
696 delete fTreeK; fTreeK = 0;
697 }
698 if (fTreeH) {
699 delete fTreeH; fTreeH = 0;
700 }
701 if (fTreeD) {
702 delete fTreeD; fTreeD = 0;
703 }
704 if (fTreeR) {
705 delete fTreeR; fTreeR = 0;
706 }
707 if (fTreeE) {
708 delete fTreeE; fTreeE = 0;
709 }
2ab0c725 710
fe4da5cc 711 // Close output file
aee8290b 712 file->Write();
fe4da5cc 713}
714
715//_____________________________________________________________________________
716void AliRun::FlagTrack(Int_t track)
717{
718 //
719 // Flags a track and all its family tree to be kept
720 //
721 int curr;
1578254f 722 TParticle *particle;
fe4da5cc 723
724 curr=track;
725 while(1) {
2ab0c725 726 particle=(TParticle*)fParticleMap->At(curr);
fe4da5cc 727
728 // If the particle is flagged the three from here upward is saved already
aee8290b 729 if(particle->TestBit(kKeepBit)) return;
fe4da5cc 730
731 // Save this particle
aee8290b 732 particle->SetBit(kKeepBit);
fe4da5cc 733
734 // Move to father if any
1578254f 735 if((curr=particle->GetFirstMother())==-1) return;
fe4da5cc 736 }
737}
738
739//_____________________________________________________________________________
740void AliRun::EnergySummary()
741{
742 //
743 // Print summary of deposited energy
744 //
745
fe4da5cc 746 Int_t ndep=0;
747 Float_t edtot=0;
748 Float_t ed, ed2;
749 Int_t kn, i, left, j, id;
aee8290b 750 const Float_t kzero=0;
fe4da5cc 751 Int_t ievent=fHeader.GetEvent()+1;
752 //
753 // Energy loss information
754 if(ievent) {
755 printf("***************** Energy Loss Information per event (GEV) *****************\n");
875c717b 756 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
757 ed=fSummEnergy[kn];
fe4da5cc 758 if(ed>0) {
875c717b 759 fEventEnergy[ndep]=kn;
fe4da5cc 760 if(ievent>1) {
761 ed=ed/ievent;
875c717b 762 ed2=fSum2Energy[kn];
fe4da5cc 763 ed2=ed2/ievent;
aee8290b 764 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
fe4da5cc 765 } else
766 ed2=99;
875c717b 767 fSummEnergy[ndep]=ed;
aee8290b 768 fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
fe4da5cc 769 edtot+=ed;
770 ndep++;
771 }
772 }
773 for(kn=0;kn<(ndep-1)/3+1;kn++) {
774 left=ndep-kn*3;
775 for(i=0;i<(3<left?3:left);i++) {
776 j=kn*3+i;
875c717b 777 id=Int_t (fEventEnergy[j]+0.1);
778 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
fe4da5cc 779 }
780 printf("\n");
781 }
782 //
783 // Relative energy loss in different detectors
784 printf("******************** Relative Energy Loss per event ********************\n");
785 printf("Total energy loss per event %10.3f GeV\n",edtot);
786 for(kn=0;kn<(ndep-1)/5+1;kn++) {
787 left=ndep-kn*5;
788 for(i=0;i<(5<left?5:left);i++) {
789 j=kn*5+i;
875c717b 790 id=Int_t (fEventEnergy[j]+0.1);
791 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
fe4da5cc 792 }
793 printf("\n");
794 }
795 for(kn=0;kn<75;kn++) printf("*");
796 printf("\n");
797 }
798 //
799 // Reset the TArray's
875c717b 800 // fEventEnergy.Set(0);
801 // fSummEnergy.Set(0);
802 // fSum2Energy.Set(0);
fe4da5cc 803}
804
805//_____________________________________________________________________________
94de3818 806AliModule *AliRun::GetModule(const char *name) const
fe4da5cc 807{
808 //
809 // Return pointer to detector from name
810 //
8494b010 811 return (AliModule*)fModules->FindObject(name);
fe4da5cc 812}
813
814//_____________________________________________________________________________
94de3818 815AliDetector *AliRun::GetDetector(const char *name) const
a68348e9 816{
817 //
818 // Return pointer to detector from name
819 //
820 return (AliDetector*)fModules->FindObject(name);
821}
822
823//_____________________________________________________________________________
94de3818 824Int_t AliRun::GetModuleID(const char *name) const
fe4da5cc 825{
826 //
827 // Return galice internal detector identifier from name
828 //
23370b7a 829 Int_t i=-1;
830 TObject *mod=fModules->FindObject(name);
831 if(mod) i=fModules->IndexOf(mod);
832 return i;
fe4da5cc 833}
834
835//_____________________________________________________________________________
836Int_t AliRun::GetEvent(Int_t event)
837{
838 //
839 // Connect the Trees Kinematics and Hits for event # event
840 // Set branch addresses
841 //
fe4da5cc 842
843 // Reset existing structures
2ab0c725 844 // ResetStack();
fe4da5cc 845 ResetHits();
846 ResetDigits();
2ab0c725 847 ResetSDigits();
fe4da5cc 848
849 // Delete Trees already connected
850 if (fTreeK) delete fTreeK;
851 if (fTreeH) delete fTreeH;
852 if (fTreeD) delete fTreeD;
853 if (fTreeR) delete fTreeR;
2ab0c725 854 if (fTreeS) delete fTreeS;
59fe9bd2 855
856 // Get header from file
857 if(fTreeE) fTreeE->GetEntry(event);
858 else Error("GetEvent","Cannot file Header Tree\n");
2ab0c725 859 TFile *file = fTreeE->GetCurrentFile();
860
861 file->cd();
fe4da5cc 862
863 // Get Kine Tree from file
864 char treeName[20];
865 sprintf(treeName,"TreeK%d",event);
866 fTreeK = (TTree*)gDirectory->Get(treeName);
2ab0c725 867 if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticleBuffer);
23370b7a 868 else Error("GetEvent","cannot find Kine Tree for event:%d\n",event);
2ab0c725 869 // Create the particle stack
870 if(!fParticles) fParticles = new TClonesArray("TParticle",1000);
871 // Build the pointer list
872 if(fParticleMap) {
873 fParticleMap->Clear();
874 fParticleMap->Expand(fTreeK->GetEntries());
875 } else
876 fParticleMap = new TObjArray(fTreeK->GetEntries());
fe4da5cc 877
2ab0c725 878 file->cd();
879
fe4da5cc 880 // Get Hits Tree header from file
881 sprintf(treeName,"TreeH%d",event);
882 fTreeH = (TTree*)gDirectory->Get(treeName);
883 if (!fTreeH) {
23370b7a 884 Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
fe4da5cc 885 }
886
2ab0c725 887 file->cd();
888
fe4da5cc 889 // Get Digits Tree header from file
890 sprintf(treeName,"TreeD%d",event);
891 fTreeD = (TTree*)gDirectory->Get(treeName);
892 if (!fTreeD) {
2ab0c725 893 // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
fe4da5cc 894 }
82711e7a 895
896 file->cd();
897
898 // Get SDigits Tree header from file
899 sprintf(treeName,"TreeS%d",event);
900 fTreeS = (TTree*)gDirectory->Get(treeName);
901 if (!fTreeS) {
902 // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
903 }
904
2ab0c725 905 file->cd();
fe4da5cc 906
907 // Get Reconstruct Tree header from file
908 sprintf(treeName,"TreeR%d",event);
909 fTreeR = (TTree*)gDirectory->Get(treeName);
910 if (!fTreeR) {
911 // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
912 }
2ab0c725 913
914 file->cd();
fe4da5cc 915
916 // Set Trees branch addresses
8494b010 917 TIter next(fModules);
918 AliModule *detector;
919 while((detector = (AliModule*)next())) {
fe4da5cc 920 detector->SetTreeAddress();
921 }
922
2ab0c725 923 fNtrack = Int_t (fTreeK->GetEntries());
fe4da5cc 924 return fNtrack;
925}
926
927//_____________________________________________________________________________
928TGeometry *AliRun::GetGeometry()
929{
930 //
931 // Import Alice geometry from current file
932 // Return pointer to geometry object
933 //
934 if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
935 //
936 // Unlink and relink nodes in detectors
937 // This is bad and there must be a better way...
938 //
fe4da5cc 939
8494b010 940 TIter next(fModules);
941 AliModule *detector;
942 while((detector = (AliModule*)next())) {
fe4da5cc 943 TList *dnodes=detector->Nodes();
944 Int_t j;
945 TNode *node, *node1;
946 for ( j=0; j<dnodes->GetSize(); j++) {
947 node = (TNode*) dnodes->At(j);
52d0ab00 948 node1 = fGeometry->GetNode(node->GetName());
fe4da5cc 949 dnodes->Remove(node);
950 dnodes->AddAt(node1,j);
951 }
952 }
953 return fGeometry;
954}
955
956//_____________________________________________________________________________
957void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
958 Float_t &e, Float_t *vpos, Float_t *polar,
959 Float_t &tof)
960{
961 //
962 // Return next track from stack of particles
963 //
a8f1fb7c 964 TVector3 pol;
fe4da5cc 965 fCurrent=-1;
1578254f 966 TParticle *track;
fe4da5cc 967 for(Int_t i=fNtrack-1; i>=0; i--) {
2ab0c725 968 track=(TParticle*) fParticleMap->At(i);
969 if(track) if(!track->TestBit(kDoneBit)) {
fe4da5cc 970 //
2ab0c725 971 // The track exists and has not yet been processed
fe4da5cc 972 fCurrent=i;
1578254f 973 ipart=track->GetPdgCode();
974 pmom[0]=track->Px();
975 pmom[1]=track->Py();
976 pmom[2]=track->Pz();
977 e =track->Energy();
978 vpos[0]=track->Vx();
979 vpos[1]=track->Vy();
980 vpos[2]=track->Vz();
a8f1fb7c 981 track->GetPolarisation(pol);
982 polar[0]=pol.X();
983 polar[1]=pol.Y();
984 polar[2]=pol.Z();
1578254f 985 tof=track->T();
aee8290b 986 track->SetBit(kDoneBit);
fe4da5cc 987 break;
988 }
989 }
990 mtrack=fCurrent;
991 //
992 // stop and start timer when we start a primary track
993 Int_t nprimaries = fHeader.GetNprimary();
994 if (fCurrent >= nprimaries) return;
995 if (fCurrent < nprimaries-1) {
996 fTimer.Stop();
2ab0c725 997 track=(TParticle*) fParticleMap->At(fCurrent+1);
1578254f 998 // track->SetProcessTime(fTimer.CpuTime());
fe4da5cc 999 }
1000 fTimer.Start();
1001}
1002
1003//_____________________________________________________________________________
94de3818 1004Int_t AliRun::GetPrimary(Int_t track) const
fe4da5cc 1005{
1006 //
1007 // return number of primary that has generated track
1008 //
1009 int current, parent;
1578254f 1010 TParticle *part;
fe4da5cc 1011 //
1012 parent=track;
1013 while (1) {
1014 current=parent;
2ab0c725 1015 part = (TParticle *)fParticleMap->At(current);
1578254f 1016 parent=part->GetFirstMother();
fe4da5cc 1017 if(parent<0) return current;
1018 }
1019}
1020
1021//_____________________________________________________________________________
875c717b 1022void AliRun::InitMC(const char *setup)
fe4da5cc 1023{
1024 //
1025 // Initialize the Alice setup
1026 //
1027
51e0e89d 1028 if(fInitDone) {
1029 Warning("Init","Cannot initialise AliRun twice!\n");
1030 return;
1031 }
2ab0c725 1032
fe4da5cc 1033 gROOT->LoadMacro(setup);
45189757 1034 gInterpreter->ProcessLine(fConfigFunction.Data());
fe4da5cc 1035
2ab0c725 1036
cfce8870 1037 gMC->DefineParticles(); //Create standard MC particles
fe4da5cc 1038
1039 TObject *objfirst, *objlast;
1040
23370b7a 1041 fNdets = fModules->GetLast()+1;
1042
fe4da5cc 1043 //
875c717b 1044 //=================Create Materials and geometry
1045 gMC->Init();
1046
65fb704d 1047 // Added also after in case of interactive initialisation of modules
1048 fNdets = fModules->GetLast()+1;
1049
8494b010 1050 TIter next(fModules);
1051 AliModule *detector;
1052 while((detector = (AliModule*)next())) {
fe4da5cc 1053 detector->SetTreeAddress();
1054 objlast = gDirectory->GetList()->Last();
1055
fe4da5cc 1056 // Add Detector histograms in Detector list of histograms
1057 if (objlast) objfirst = gDirectory->GetList()->After(objlast);
1058 else objfirst = gDirectory->GetList()->First();
1059 while (objfirst) {
1060 detector->Histograms()->Add(objfirst);
1061 objfirst = gDirectory->GetList()->After(objfirst);
1062 }
1063 }
327136d2 1064 ReadTransPar(); //Read the cuts for all materials
fe4da5cc 1065
1066 MediaTable(); //Build the special IMEDIA table
1067
fe4da5cc 1068 //Initialise geometry deposition table
875c717b 1069 fEventEnergy.Set(gMC->NofVolumes()+1);
1070 fSummEnergy.Set(gMC->NofVolumes()+1);
1071 fSum2Energy.Set(gMC->NofVolumes()+1);
fe4da5cc 1072
fe4da5cc 1073 //Compute cross-sections
875c717b 1074 gMC->BuildPhysics();
fe4da5cc 1075
1076 //Write Geometry object to current file.
1077 fGeometry->Write();
1078
1079 fInitDone = kTRUE;
51e0e89d 1080
65fb704d 1081 fMCQA = new AliMCQA(fNdets);
1082
51e0e89d 1083 //
1084 // Save stuff at the beginning of the file to avoid file corruption
1085 Write();
fe4da5cc 1086}
1087
1088//_____________________________________________________________________________
1089void AliRun::MediaTable()
1090{
1091 //
1092 // Built media table to get from the media number to
1093 // the detector id
1094 //
ad51aeb0 1095 Int_t kz, nz, idt, lz, i, k, ind;
1096 // Int_t ibeg;
fe4da5cc 1097 TObjArray &dets = *gAlice->Detectors();
8494b010 1098 AliModule *det;
fe4da5cc 1099 //
1100 // For all detectors
1101 for (kz=0;kz<fNdets;kz++) {
1102 // If detector is defined
8494b010 1103 if((det=(AliModule*) dets[kz])) {
ad51aeb0 1104 TArrayI &idtmed = *(det->GetIdtmed());
1105 for(nz=0;nz<100;nz++) {
fe4da5cc 1106 // Find max and min material number
ad51aeb0 1107 if((idt=idtmed[nz])) {
fe4da5cc 1108 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
1109 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
1110 }
1111 }
1112 if(det->LoMedium() > det->HiMedium()) {
1113 det->LoMedium() = 0;
1114 det->HiMedium() = 0;
1115 } else {
1116 if(det->HiMedium() > fImedia->GetSize()) {
ad51aeb0 1117 Error("MediaTable","Increase fImedia from %d to %d",
1118 fImedia->GetSize(),det->HiMedium());
fe4da5cc 1119 return;
1120 }
1121 // Tag all materials in rage as belonging to detector kz
1122 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
1123 (*fImedia)[lz]=kz;
1124 }
1125 }
1126 }
1127 }
1128 //
1129 // Print summary table
1130 printf(" Traking media ranges:\n");
1131 for(i=0;i<(fNdets-1)/6+1;i++) {
1132 for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
1133 ind=i*6+k;
8494b010 1134 det=(AliModule*)dets[ind];
fe4da5cc 1135 if(det)
1136 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
1137 det->HiMedium());
1138 else
1139 printf(" %6s: %3d -> %3d;","NULL",0,0);
1140 }
1141 printf("\n");
1142 }
1143}
1144
1145//____________________________________________________________________________
1146void AliRun::SetGenerator(AliGenerator *generator)
1147{
1148 //
1149 // Load the event generator
1150 //
ee1dd322 1151 if(!fGenerator) fGenerator = generator;
1152}
1153
1154//____________________________________________________________________________
1155void AliRun::ResetGenerator(AliGenerator *generator)
1156{
1157 //
1158 // Load the event generator
1159 //
b13db077 1160 if(fGenerator)
838edcaf 1161 if(generator)
1162 Warning("ResetGenerator","Replacing generator %s with %s\n",
1163 fGenerator->GetName(),generator->GetName());
1164 else
1165 Warning("ResetGenerator","Replacing generator %s with NULL\n",
1166 fGenerator->GetName());
b13db077 1167 fGenerator = generator;
fe4da5cc 1168}
1169
1170//____________________________________________________________________________
65fb704d 1171void AliRun::SetTransPar(char *filename)
1172{
1173 fTransParName = filename;
1174}
1175
1176//____________________________________________________________________________
2ab0c725 1177void AliRun::SetBaseFile(char *filename)
1178{
39de14fb 1179 fBaseFileName = filename;
2ab0c725 1180}
1181
1182//____________________________________________________________________________
65fb704d 1183void AliRun::ReadTransPar()
fe4da5cc 1184{
1185 //
1186 // Read filename to set the transport parameters
1187 //
1188
fe4da5cc 1189
aee8290b 1190 const Int_t kncuts=10;
1191 const Int_t knflags=11;
1192 const Int_t knpars=kncuts+knflags;
1193 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
fe4da5cc 1194 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
1195 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
1196 "MULS","PAIR","PHOT","RAYL"};
1197 char line[256];
ad51aeb0 1198 char detName[7];
fe4da5cc 1199 char* filtmp;
aee8290b 1200 Float_t cut[kncuts];
1201 Int_t flag[knflags];
fe4da5cc 1202 Int_t i, itmed, iret, ktmed, kz;
1203 FILE *lun;
1204 //
1205 // See whether the file is there
65fb704d 1206 filtmp=gSystem->ExpandPathName(fTransParName.Data());
fe4da5cc 1207 lun=fopen(filtmp,"r");
1208 delete [] filtmp;
1209 if(!lun) {
65fb704d 1210 Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
fe4da5cc 1211 return;
1212 }
1213 //
1214 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1215 printf(" *%59s\n","*");
1216 printf(" * Please check carefully what you are doing!%10s\n","*");
1217 printf(" *%59s\n","*");
1218 //
1219 while(1) {
1220 // Initialise cuts and flags
aee8290b 1221 for(i=0;i<kncuts;i++) cut[i]=-99;
1222 for(i=0;i<knflags;i++) flag[i]=-99;
fe4da5cc 1223 itmed=0;
1224 for(i=0;i<256;i++) line[i]='\0';
1225 // Read up to the end of line excluded
1226 iret=fscanf(lun,"%[^\n]",line);
1227 if(iret<0) {
1228 //End of file
1229 fclose(lun);
1230 printf(" *%59s\n","*");
1231 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1232 return;
1233 }
1234 // Read the end of line
1235 fscanf(lun,"%*c");
1236 if(!iret) continue;
1237 if(line[0]=='*') continue;
1238 // Read the numbers
ad51aeb0 1239 iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
1240 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1241 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1242 &flag[8],&flag[9],&flag[10]);
fe4da5cc 1243 if(!iret) continue;
1244 if(iret<0) {
1245 //reading error
65fb704d 1246 Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
fe4da5cc 1247 continue;
1248 }
ad51aeb0 1249 // Check that the module exist
1250 AliModule *mod = GetModule(detName);
1251 if(mod) {
1252 // Get the array of media numbers
1253 TArrayI &idtmed = *mod->GetIdtmed();
1254 // Check that the tracking medium code is valid
1255 if(0<=itmed && itmed < 100) {
1256 ktmed=idtmed[itmed];
1257 if(!ktmed) {
65fb704d 1258 Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
ad51aeb0 1259 continue;
fe4da5cc 1260 }
ad51aeb0 1261 // Set energy thresholds
aee8290b 1262 for(kz=0;kz<kncuts;kz++) {
ad51aeb0 1263 if(cut[kz]>=0) {
23370b7a 1264 printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
aee8290b 1265 kpars[kz],cut[kz],itmed,mod->GetName());
1266 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
ad51aeb0 1267 }
fe4da5cc 1268 }
ad51aeb0 1269 // Set transport mechanisms
aee8290b 1270 for(kz=0;kz<knflags;kz++) {
ad51aeb0 1271 if(flag[kz]>=0) {
1272 printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
aee8290b 1273 kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
1274 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
ad51aeb0 1275 }
1276 }
1277 } else {
65fb704d 1278 Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
ad51aeb0 1279 continue;
fe4da5cc 1280 }
1281 } else {
65fb704d 1282 Warning("ReadTransPar","Module %s not present\n",detName);
fe4da5cc 1283 continue;
1284 }
1285 }
1286}
1287
1288//_____________________________________________________________________________
2ab0c725 1289void AliRun::MakeBranchInTree(TTree *tree, const char* name, void* address, Int_t size, char *file)
1290{
d47c658f 1291 //
1292 // Makes branch in given tree and diverts them to a separate file
1293 //
2ab0c725 1294 if (GetDebug()>1)
1295 printf("* MakeBranch * Making Branch %s \n",name);
1296
1297 TBranch *branch = tree->Branch(name,address,size);
1298
1299 if (file) {
39de14fb 1300 char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
1301 sprintf(outFile,"%s/%s",gAlice->GetBaseFile(),file);
2ab0c725 1302 TDirectory *cwd = gDirectory;
39de14fb 1303 branch->SetFile(outFile);
2ab0c725 1304 TIter next( branch->GetListOfBranches());
1305 while ((branch=(TBranch*)next())) {
39de14fb 1306 branch->SetFile(outFile);
2ab0c725 1307 }
1308 if (GetDebug()>1)
1309 printf("* MakeBranch * Diverting Branch %s to file %s\n",name,file);
1310 cwd->cd();
39de14fb 1311 delete outFile;
2ab0c725 1312 }
1313}
1314
1315//_____________________________________________________________________________
1316void AliRun::MakeBranchInTree(TTree *tree, const char* name, const char *classname, void* address, Int_t size, Int_t splitlevel, char *file)
1317{
d47c658f 1318 //
1319 // Makes branch in given tree and diverts them to a separate file
1320 //
2ab0c725 1321 TDirectory *cwd = gDirectory;
1322 TBranch *branch = tree->Branch(name,classname,address,size,splitlevel);
1323 if (GetDebug()>1)
1324 printf("* MakeBranch * Making Branch %s \n",name);
1325 if (file) {
39de14fb 1326 char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
1327 sprintf(outFile,"%s/%s",gAlice->GetBaseFile(),file);
1328 branch->SetFile(outFile);
2ab0c725 1329 TIter next( branch->GetListOfBranches());
1330 while ((branch=(TBranch*)next())) {
39de14fb 1331 branch->SetFile(outFile);
2ab0c725 1332 }
1333 if (GetDebug()>1)
1334 printf("* MakeBranch * Diverting Branch %s to file %s\n",name,file);
1335 cwd->cd();
39de14fb 1336 delete outFile;
2ab0c725 1337 }
1338}
1339//_____________________________________________________________________________
1340void AliRun::MakeTree(Option_t *option, char *file)
fe4da5cc 1341{
1342 //
1343 // Create the ROOT trees
1344 // Loop on all detectors to create the Root branch (if any)
1345 //
1346
b13db077 1347 char hname[30];
fe4da5cc 1348 //
1349 // Analyse options
aee8290b 1350 char *oK = strstr(option,"K");
1351 char *oH = strstr(option,"H");
1352 char *oE = strstr(option,"E");
1353 char *oD = strstr(option,"D");
1354 char *oR = strstr(option,"R");
2ab0c725 1355 char *oS = strstr(option,"S");
fe4da5cc 1356 //
2ab0c725 1357
aee8290b 1358 if (oK && !fTreeK) {
b13db077 1359 sprintf(hname,"TreeK%d",fEvent);
1360 fTreeK = new TTree(hname,"Kinematics");
1361 // Create a branch for particles
2ab0c725 1362 MakeBranchInTree(fTreeK,
1363 "Particles", "TParticle", &fParticleBuffer, 4000, 1, file) ;
51e0e89d 1364 fTreeK->Write();
b13db077 1365 }
aee8290b 1366 if (oH && !fTreeH) {
b13db077 1367 sprintf(hname,"TreeH%d",fEvent);
1368 fTreeH = new TTree(hname,"Hits");
1369 fTreeH->SetAutoSave(1000000000); //no autosave
51e0e89d 1370 fTreeH->Write();
b13db077 1371 }
aee8290b 1372 if (oD && !fTreeD) {
b13db077 1373 sprintf(hname,"TreeD%d",fEvent);
1374 fTreeD = new TTree(hname,"Digits");
51e0e89d 1375 fTreeD->Write();
b13db077 1376 }
2ab0c725 1377 if (oS && !fTreeS) {
1378 sprintf(hname,"TreeS%d",fEvent);
1379 fTreeS = new TTree(hname,"SDigits");
1380 fTreeS->Write();
1381 }
aee8290b 1382 if (oR && !fTreeR) {
b13db077 1383 sprintf(hname,"TreeR%d",fEvent);
1384 fTreeR = new TTree(hname,"Reconstruction");
51e0e89d 1385 fTreeR->Write();
b13db077 1386 }
aee8290b 1387 if (oE && !fTreeE) {
b13db077 1388 fTreeE = new TTree("TE","Header");
1389 // Create a branch for Header
2ab0c725 1390 MakeBranchInTree(fTreeE,
1391 "Header", "AliHeader", &gAliHeader, 4000, 1, file) ;
51e0e89d 1392 fTreeE->Write();
b13db077 1393 }
2ab0c725 1394
fe4da5cc 1395 //
1396 // Create a branch for hits/digits for each detector
1397 // Each branch is a TClonesArray. Each data member of the Hits classes
1398 // will be in turn a subbranch of the detector master branch
8494b010 1399 TIter next(fModules);
1400 AliModule *detector;
1401 while((detector = (AliModule*)next())) {
d47c658f 1402 if (oH) detector->MakeBranch(option,file);
fe4da5cc 1403 }
fe4da5cc 1404}
1405
1406//_____________________________________________________________________________
1407Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks)
1408{
1409 //
1410 // PurifyKine with external parameters
1411 //
1412 fHgwmk = lastSavedTrack;
1413 fNtrack = nofTracks;
1414 PurifyKine();
1415 return fHgwmk;
1416}
1417
1418//_____________________________________________________________________________
2ab0c725 1419TParticle* AliRun::Particle(Int_t i)
1420{
1421 if(!(*fParticleMap)[i]) {
1422 Int_t nentries = fParticles->GetEntries();
a41d61a9 1423
1424 // algorithmic way of getting entry index
1425 // (primary particles are filled after secondaries)
1426 Int_t entry;
1427 if (i<fHeader.GetNprimary())
1428 entry = i+fHeader.GetNsecondary();
1429 else
1430 entry = i-fHeader.GetNprimary();
1431
1432 // only check the algorithmic way and give
1433 // the fatal error if it is wrong
1434 if (entry != fParticleFileMap[i]) {
1435 Fatal("Particle",
1436 "!!!! The algorithmic way is WRONG: !!!\n entry: %d map: %d",
1437 entry, fParticleFileMap[i]);
1438 }
1439
2ab0c725 1440 fTreeK->GetEntry(fParticleFileMap[i]);
1441 new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
1442 fParticleMap->AddAt((*fParticles)[nentries],i);
1443 }
1444 return (TParticle *) (*fParticleMap)[i];
1445}
1446
1447//_____________________________________________________________________________
fe4da5cc 1448void AliRun::PurifyKine()
1449{
1450 //
1451 // Compress kinematic tree keeping only flagged particles
1452 // and renaming the particle id's in all the hits
1453 //
2ab0c725 1454 // TClonesArray &particles = *fParticles;
1455 TObjArray &particles = *fParticleMap;
fe4da5cc 1456 int nkeep=fHgwmk+1, parent, i;
2ab0c725 1457 TParticle *part, *father;
1458 TArrayI map(particles.GetLast()+1);
fe4da5cc 1459
1460 // Save in Header total number of tracks before compression
1461 fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk);
1462
2ab0c725 1463 // If no tracks generated return now
1464 if(fHgwmk+1 == fNtrack) return;
1465
1466 Int_t toshrink = fNtrack-fHgwmk-1;
1467
a01a8b12 1468 // First pass, invalid Daughter information
fe4da5cc 1469 for(i=0; i<fNtrack; i++) {
a01a8b12 1470 // Preset map, to be removed later
2ab0c725 1471 if(i<=fHgwmk) map[i]=i ;
1472 else {
1473 map[i] = -99;
1474 // particles.UncheckedAt(i)->ResetBit(kDaughtersBit);
1475 if((part=(TParticle*) particles.At(i))) part->ResetBit(kDaughtersBit);
1476 }
a01a8b12 1477 }
2ab0c725 1478 // Invalid daughter information for the parent of the first particle
1479 // generated. This may or may not be the current primary according to
1480 // whether decays have been recorded among the primaries
1481 part = (TParticle *)particles.At(fHgwmk+1);
1482 particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
fe4da5cc 1483 // Second pass, build map between old and new numbering
1484 for(i=fHgwmk+1; i<fNtrack; i++) {
2ab0c725 1485 if(particles.At(i)->TestBit(kKeepBit)) {
fe4da5cc 1486
1487 // This particle has to be kept
1488 map[i]=nkeep;
2ab0c725 1489 // If old and new are different, have to move the pointer
1490 if(i!=nkeep) particles[nkeep]=particles.At(i);
1491 part = (TParticle*) particles.At(nkeep);
fe4da5cc 1492
1493 // as the parent is always *before*, it must be already
1494 // in place. This is what we are checking anyway!
2ab0c725 1495 if((parent=part->GetFirstMother())>fHgwmk)
1496 if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
1497 else part->SetFirstMother(map[parent]);
1498
fe4da5cc 1499 nkeep++;
1500 }
1501 }
fe4da5cc 1502
1578254f 1503 // Fix daughters information
2ab0c725 1504 for (i=fHgwmk+1; i<nkeep; i++) {
1505 part = (TParticle *)particles.At(i);
1578254f 1506 parent = part->GetFirstMother();
ae23d366 1507 if(parent>=0) {
2ab0c725 1508 father = (TParticle *)particles.At(parent);
aee8290b 1509 if(father->TestBit(kDaughtersBit)) {
fe4da5cc 1510
ae23d366 1511 if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
1512 if(i>father->GetLastDaughter()) father->SetLastDaughter(i);
1513 } else {
2ab0c725 1514 // Initialise daughters info for first pass
ae23d366 1515 father->SetFirstDaughter(i);
1516 father->SetLastDaughter(i);
aee8290b 1517 father->SetBit(kDaughtersBit);
ae23d366 1518 }
fe4da5cc 1519 }
1520 }
1521
1cedd08a 1522 // Now loop on all registered hit lists
1523 TIter next(fHitLists);
1524 TCollection *hitList;
1525 while((hitList = (TCollection*)next())) {
1526 TIter nexthit(hitList);
1527 AliHit *hit;
1528 while((hit = (AliHit*)nexthit())) {
1529 hit->SetTrack(map[hit->GetTrack()]);
1530 }
1531 }
fe4da5cc 1532
24de2263 1533 //
1534 // This for detectors which have a special mapping mechanism
1535 // for hits, such as TPC and TRD
1536 //
1537
1538 TIter nextmod(fModules);
1539 AliModule *detector;
1540 while((detector = (AliModule*)nextmod())) {
2ab0c725 1541 detector->RemapTrackHitIDs(map.GetArray());
24de2263 1542 }
1543
2ab0c725 1544 // Now the output bit, from fHgwmk to nkeep we write everything and we erase
1545 if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
1546
24de2263 1547
2ab0c725 1548 for (i=fHgwmk+1; i<nkeep; ++i) {
1549 fParticleBuffer = (TParticle*) particles.At(i);
1550 fParticleFileMap[i]=(Int_t) fTreeK->GetEntries();
1551 fTreeK->Fill();
1552 particles[i]=0;
1553 }
1554
1555 for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
1556
1557 fLoadPoint-=toshrink;
1558 for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
1559
1560 fNtrack=nkeep;
1561 fHgwmk=nkeep-1;
1562 // delete [] map;
fe4da5cc 1563}
1564
1565//_____________________________________________________________________________
dffd31ef 1566void AliRun::BeginEvent()
fe4da5cc 1567{
1568 //
1569 // Reset all Detectors & kinematics & trees
1570 //
59fe9bd2 1571 char hname[30];
1572 //
dffd31ef 1573
1574 //
1575 if(fLego) {
1576 fLego->BeginEvent();
1577 return;
1578 }
1579
1580 //
fe4da5cc 1581 ResetStack();
1582 ResetHits();
1583 ResetDigits();
2ab0c725 1584 ResetSDigits();
fe4da5cc 1585
1586 // Initialise event header
875c717b 1587 fHeader.Reset(fRun,fEvent);
fe4da5cc 1588
59fe9bd2 1589 if(fTreeK) {
1590 fTreeK->Reset();
875c717b 1591 sprintf(hname,"TreeK%d",fEvent);
59fe9bd2 1592 fTreeK->SetName(hname);
1593 }
1594 if(fTreeH) {
1595 fTreeH->Reset();
875c717b 1596 sprintf(hname,"TreeH%d",fEvent);
59fe9bd2 1597 fTreeH->SetName(hname);
1598 }
1599 if(fTreeD) {
1600 fTreeD->Reset();
875c717b 1601 sprintf(hname,"TreeD%d",fEvent);
59fe9bd2 1602 fTreeD->SetName(hname);
1603 }
2ab0c725 1604 if(fTreeS) {
1605 fTreeS->Reset();
1606 sprintf(hname,"TreeS%d",fEvent);
1607 fTreeS->SetName(hname);
1608 }
59fe9bd2 1609 if(fTreeR) {
1610 fTreeR->Reset();
875c717b 1611 sprintf(hname,"TreeR%d",fEvent);
59fe9bd2 1612 fTreeR->SetName(hname);
1613 }
fe4da5cc 1614}
fe4da5cc 1615//_____________________________________________________________________________
1616void AliRun::ResetDigits()
1617{
1618 //
1619 // Reset all Detectors digits
1620 //
8494b010 1621 TIter next(fModules);
1622 AliModule *detector;
1623 while((detector = (AliModule*)next())) {
fe4da5cc 1624 detector->ResetDigits();
1625 }
1626}
1627
1628//_____________________________________________________________________________
2ab0c725 1629void AliRun::ResetSDigits()
1630{
1631 //
1632 // Reset all Detectors digits
1633 //
1634 TIter next(fModules);
1635 AliModule *detector;
1636 while((detector = (AliModule*)next())) {
1637 detector->ResetSDigits();
1638 }
1639}
1640
1641//_____________________________________________________________________________
fe4da5cc 1642void AliRun::ResetHits()
1643{
1644 //
1645 // Reset all Detectors hits
1646 //
8494b010 1647 TIter next(fModules);
1648 AliModule *detector;
1649 while((detector = (AliModule*)next())) {
fe4da5cc 1650 detector->ResetHits();
1651 }
1652}
1653
1654//_____________________________________________________________________________
1655void AliRun::ResetPoints()
1656{
1657 //
1658 // Reset all Detectors points
1659 //
8494b010 1660 TIter next(fModules);
1661 AliModule *detector;
1662 while((detector = (AliModule*)next())) {
fe4da5cc 1663 detector->ResetPoints();
1664 }
1665}
1666
1667//_____________________________________________________________________________
875c717b 1668void AliRun::RunMC(Int_t nevent, const char *setup)
fe4da5cc 1669{
1670 //
1671 // Main function to be called to process a galice run
1672 // example
1673 // Root > gAlice.Run();
1674 // a positive number of events will cause the finish routine
1675 // to be called
1676 //
1677
fe4da5cc 1678 // check if initialisation has been done
875c717b 1679 if (!fInitDone) InitMC(setup);
fe4da5cc 1680
fe4da5cc 1681 // Create the Root Tree with one branch per detector
2ab0c725 1682
d47c658f 1683 MakeTree("ESD");
1684
2ab0c725 1685 if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
39de14fb 1686 MakeTree("K","Kine.root");
1687 MakeTree("H","Hits.root");
2ab0c725 1688 } else {
d47c658f 1689 MakeTree("KH");
2ab0c725 1690 }
fe4da5cc 1691
875c717b 1692 gMC->ProcessRun(nevent);
1693
fe4da5cc 1694 // End of this run, close files
80762cb1 1695 if(nevent>0) FinishRun();
fe4da5cc 1696}
1697
1698//_____________________________________________________________________________
d47c658f 1699void AliRun::RunReco(const char *detector)
1700{
1701 //
1702 // Main function to be called to reconstruct Alice event
1703 //
1704
1705 MakeTree("R");
1706 Digits2Reco(detector);
1707}
1708
1709//_____________________________________________________________________________
2ab0c725 1710
1711void AliRun::Hits2Digits(const char *selected)
1712{
d47c658f 1713 // Convert Hits to sumable digits
1714 //
2ab0c725 1715 Hits2SDigits(selected);
1716 SDigits2Digits(selected);
1717}
1718
d47c658f 1719
2ab0c725 1720//_____________________________________________________________________________
1721
d47c658f 1722void AliRun::Tree2Tree(Option_t *option, const char *selected)
2ab0c725 1723{
1724 //
d47c658f 1725 // Function to transform the content of
1726 //
1727 // - TreeH to TreeS (option "S")
1728 // - TreeS to TreeD (option "D")
1729 // - TreeD to TreeR (option "R")
1730 //
1731 // If multiple options are specified ("SDR"), transformation will be done in sequence for
1732 // selected detector and for all detectors if none is selected (detector string
1733 // can contain blank separated list of detector names).
2ab0c725 1734
2ab0c725 1735
d47c658f 1736 char *oS = strstr(option,"S");
1737 char *oD = strstr(option,"D");
1738 char *oR = strstr(option,"R");
2ab0c725 1739
2ab0c725 1740 gAlice->GetEvent(0);
1741
1742 TObjArray *detectors = gAlice->Detectors();
1743
1744 TIter next(detectors);
1745
d47c658f 1746 AliDetector *detector = 0;
2ab0c725 1747
1748 TDirectory *cwd = gDirectory;
1749
d47c658f 1750 char outFile[32];
2ab0c725 1751
1752 while((detector = (AliDetector*)next())) {
d47c658f 1753 if (selected)
2ab0c725 1754 if (strcmp(detector->GetName(),selected)) continue;
2ab0c725 1755 if (detector->IsActive()){
d47c658f 1756 if (GetDebug()>0)
2ab0c725 1757 cout << "Processing " << detector->GetName() << "..." << endl;
d47c658f 1758 if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
1759 if (oS) {
1760 sprintf(outFile,"SDigits.%s.root",detector->GetName());
1761 detector->MakeBranch("S",outFile);
1762 }
1763 if (oD) {
1764 sprintf(outFile,"Digits.%s.root",detector->GetName());
1765 detector->MakeBranch("D",outFile);
1766 }
1767 if (oR) {
1768 sprintf(outFile,"Reco.%s.root",detector->GetName());
1769 detector->MakeBranch("R",outFile);
1770 }
2ab0c725 1771 } else {
d47c658f 1772 detector->MakeBranch(option);
1773 }
1774
1775 cwd->cd();
1776
1777 if (oS)
1778 detector->Hits2SDigits();
1779 if (oD)
1780 detector->SDigits2Digits();
1781 if (oR)
1782 detector->Digits2Reco();
1783
1784 cwd->cd();
1785
2ab0c725 1786 }
1787 }
1788}
1789
d47c658f 1790
2ab0c725 1791//_____________________________________________________________________________
0a520a66 1792void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1793 Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1794 Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
fe4da5cc 1795{
1796 //
1797 // Generates lego plots of:
1798 // - radiation length map phi vs theta
1799 // - radiation length map phi vs eta
1800 // - interaction length map
1801 // - g/cm2 length map
1802 //
1803 // ntheta bins in theta, eta
1804 // themin minimum angle in theta (degrees)
1805 // themax maximum angle in theta (degrees)
1806 // nphi bins in phi
1807 // phimin minimum angle in phi (degrees)
1808 // phimax maximum angle in phi (degrees)
1809 // rmin minimum radius
1810 // rmax maximum radius
1811 //
1812 //
1813 // The number of events generated = ntheta*nphi
1814 // run input parameters in macro setup (default="Config.C")
1815 //
1816 // Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1817 //Begin_Html
1818 /*
1439f98e 1819 <img src="picts/AliRunLego1.gif">
fe4da5cc 1820 */
1821 //End_Html
1822 //Begin_Html
1823 /*
1439f98e 1824 <img src="picts/AliRunLego2.gif">
fe4da5cc 1825 */
1826 //End_Html
1827 //Begin_Html
1828 /*
1439f98e 1829 <img src="picts/AliRunLego3.gif">
fe4da5cc 1830 */
1831 //End_Html
1832 //
1833
1834 // check if initialisation has been done
875c717b 1835 if (!fInitDone) InitMC(setup);
838edcaf 1836 //Save current generator
1837 AliGenerator *gen=Generator();
1838
0a520a66 1839 // Set new generator
1840 if (!gener) gener = new AliLegoGenerator();
1841 ResetGenerator(gener);
1842 //
1843 // Configure Generator
1844 gener->SetRadiusRange(rmin, rmax);
1845 gener->SetZMax(zmax);
1846 gener->SetCoor1Range(nc1, c1min, c1max);
1847 gener->SetCoor2Range(nc2, c2min, c2max);
1848
1849
b13db077 1850 //Create Lego object
0a520a66 1851 fLego = new AliLego("lego",gener);
b13db077 1852
dffd31ef 1853 //Prepare MC for Lego Run
1854 gMC->InitLego();
1855
b13db077 1856 //Run Lego Object
0a520a66 1857
1858 gMC->ProcessRun(nc1*nc2+1);
fe4da5cc 1859
1860 // Create only the Root event Tree
80762cb1 1861 MakeTree("E");
fe4da5cc 1862
1863 // End of this run, close files
80762cb1 1864 FinishRun();
0a520a66 1865 // Restore current generator
1866 ResetGenerator(gen);
838edcaf 1867 // Delete Lego Object
1868 delete fLego; fLego=0;
fe4da5cc 1869}
1870
1871//_____________________________________________________________________________
94de3818 1872void AliRun::SetConfigFunction(const char * config)
1873{
1874 //
1875 // Set the signature of the function contained in Config.C to configure
1876 // the run
1877 //
1878 fConfigFunction=config;
1879}
1880
1881//_____________________________________________________________________________
fe4da5cc 1882void AliRun::SetCurrentTrack(Int_t track)
1883{
1884 //
1885 // Set current track number
1886 //
1887 fCurrent = track;
1888}
1889
1890//_____________________________________________________________________________
1578254f 1891void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
fe4da5cc 1892 Float_t *vpos, Float_t *polar, Float_t tof,
65fb704d 1893 AliMCProcess mech, Int_t &ntr, Float_t weight)
fe4da5cc 1894{
1895 //
1896 // Load a track on the stack
1897 //
1898 // done 0 if the track has to be transported
1899 // 1 if not
1900 // parent identifier of the parent track. -1 for a primary
1578254f 1901 // pdg particle code
fe4da5cc 1902 // pmom momentum GeV/c
1903 // vpos position
1904 // polar polarisation
1905 // tof time of flight in seconds
1906 // mecha production mechanism
1907 // ntr on output the number of the track stored
1908 //
1909 TClonesArray &particles = *fParticles;
1578254f 1910 TParticle *particle;
fe4da5cc 1911 Float_t mass;
aee8290b 1912 const Int_t kfirstdaughter=-1;
1913 const Int_t klastdaughter=-1;
1914 const Int_t kS=0;
1578254f 1915 // const Float_t tlife=0;
fe4da5cc 1916
1578254f 1917 //
1918 // Here we get the static mass
1919 // For MC is ok, but a more sophisticated method could be necessary
1920 // if the calculated mass is required
1921 // also, this method is potentially dangerous if the mass
1922 // used in the MC is not the same of the PDG database
1923 //
1924 mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
fe4da5cc 1925 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1926 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1927
aee8290b 1928 //printf("Loading particle %s mass %f ene %f No %d ip %d pos %f %f %f mom %f %f %f kS %d m %s\n",
1929 //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS,mecha);
2ab0c725 1930
1931 particle=new(particles[fLoadPoint++]) TParticle(pdg,kS,parent,-1,kfirstdaughter,
aee8290b 1932 klastdaughter,pmom[0],pmom[1],pmom[2],
1578254f 1933 e,vpos[0],vpos[1],vpos[2],tof);
2ab0c725 1934 particle->SetPolarisation(TVector3(polar[0],polar[1],polar[2]));
1935 particle->SetWeight(weight);
65fb704d 1936 particle->SetUniqueID(mech);
aee8290b 1937 if(!done) particle->SetBit(kDoneBit);
2ab0c725 1938 // Declare that the daughter information is valid
1939 particle->SetBit(kDaughtersBit);
1940 // Add the particle to the stack
1941 fParticleMap->AddAtAndExpand(particle,fNtrack);
1942
fe4da5cc 1943 if(parent>=0) {
2ab0c725 1944 particle=(TParticle*) fParticleMap->At(parent);
1578254f 1945 particle->SetLastDaughter(fNtrack);
1946 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
fe4da5cc 1947 } else {
1948 //
1949 // This is a primary track. Set high water mark for this event
1950 fHgwmk=fNtrack;
1951 //
1952 // Set also number if primary tracks
1953 fHeader.SetNprimary(fHgwmk+1);
1954 fHeader.SetNtrack(fHgwmk+1);
1955 }
1956 ntr = fNtrack++;
89bbad6f 1957
1958/*
1959 //
1960 // Here we get the static mass
1961 // For MC is ok, but a more sophisticated method could be necessary
1962 // if the calculated mass is required
1963 // also, this method is potentially dangerous if the mass
1964 // used in the MC is not the same of the PDG database
1965 //
1966 Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1967 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1968 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1969
1970 SetTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
1971 vpos[0], vpos[1], vpos[2], tof, polar[0],polar[1],polar[2],
1972 mech, ntr, weight);
1973*/
fe4da5cc 1974}
1975
89bbad6f 1976//_____________________________________________________________________________
1977void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg,
1978 Double_t px, Double_t py, Double_t pz, Double_t e,
1979 Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1980 Double_t polx, Double_t poly, Double_t polz,
1d994b80 1981 AliMCProcess mech, Int_t &ntr, Float_t weight)
89bbad6f 1982{
1983 //
1984 // Load a track on the stack
1985 //
1986 // done 0 if the track has to be transported
1987 // 1 if not
1988 // parent identifier of the parent track. -1 for a primary
1989 // pdg particle code
1990 // kS generation status code
1991 // px, py, pz momentum GeV/c
1992 // vx, vy, vz position
1993 // polar polarisation
1994 // tof time of flight in seconds
1995 // mech production mechanism
1996 // ntr on output the number of the track stored
1997 //
1998 // New method interface:
1999 // arguments were changed to be in correspondence with TParticle
2000 // constructor.
2001 // Note: the energy is not calculated from the static mass but
2002 // it is passed by argument e.
2003
2004 TClonesArray &particles = *fParticles;
2005
2006 const Int_t kS=0;
2007 const Int_t kFirstDaughter=-1;
2008 const Int_t kLastDaughter=-1;
2009
2010 TParticle* particle
2011 = new(particles[fLoadPoint++]) TParticle(pdg, kS, parent, -1,
2012 kFirstDaughter, kLastDaughter,
2013 px, py, pz, e, vx, vy, vz, tof);
2014
2015 particle->SetPolarisation(polx, poly, polz);
2016 particle->SetWeight(weight);
2017 particle->SetUniqueID(mech);
2018
2019 if(!done) particle->SetBit(kDoneBit);
2020
2021 // Declare that the daughter information is valid
2022 particle->SetBit(kDaughtersBit);
2023 // Add the particle to the stack
2024 fParticleMap->AddAtAndExpand(particle,fNtrack);
2025
2026 if(parent>=0) {
2027 particle=(TParticle*) fParticleMap->At(parent);
2028 particle->SetLastDaughter(fNtrack);
2029 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
2030 } else {
2031 //
2032 // This is a primary track. Set high water mark for this event
2033 fHgwmk=fNtrack;
2034 //
2035 // Set also number if primary tracks
2036 fHeader.SetNprimary(fHgwmk+1);
2037 fHeader.SetNtrack(fHgwmk+1);
2038 }
2039 ntr = fNtrack++;
2040}
2041
2042//_____________________________________________________________________________
4d69d91e 2043void AliRun::SetHighWaterMark(const Int_t nt)
2044{
2045 //
2046 // Set high water mark for last track in event
2047 fHgwmk=fNtrack-1;
2048 //
2049 // Set also number if primary tracks
2050 fHeader.SetNprimary(fHgwmk+1);
2051 fHeader.SetNtrack(fHgwmk+1);
2052}
2053
fe4da5cc 2054//_____________________________________________________________________________
2055void AliRun::KeepTrack(const Int_t track)
2056{
2057 //
2058 // flags a track to be kept
2059 //
2ab0c725 2060 fParticleMap->At(track)->SetBit(kKeepBit);
fe4da5cc 2061}
2062
2063//_____________________________________________________________________________
875c717b 2064void AliRun::StepManager(Int_t id)
fe4da5cc 2065{
2066 //
2067 // Called at every step during transport
2068 //
2069
fe4da5cc 2070 //
2071 // --- If lego option, do it and leave
b13db077 2072 if (fLego)
fe4da5cc 2073 fLego->StepManager();
b13db077 2074 else {
2075 Int_t copy;
2076 //Update energy deposition tables
875c717b 2077 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
fe4da5cc 2078
b13db077 2079 //Call the appropriate stepping routine;
2080 AliModule *det = (AliModule*)fModules->At(id);
65fb704d 2081 if(det) {
2082 fMCQA->StepManager(id);
2083 det->StepManager();
2084 }
fe4da5cc 2085 }
fe4da5cc 2086}
2087
2088//_____________________________________________________________________________
fe4da5cc 2089void AliRun::Streamer(TBuffer &R__b)
2090{
2ab0c725 2091 // Stream an object of class AliRun.
2092
2093 if (R__b.IsReading()) {
2094 if (!gAlice) gAlice = this;
2095
2096 AliRun::Class()->ReadBuffer(R__b, this);
2097 //
2098 gROOT->GetListOfBrowsables()->Add(this,"Run");
2099
2100 fTreeE = (TTree*)gDirectory->Get("TE");
2101 if (fTreeE) fTreeE->SetBranchAddress("Header", &gAliHeader);
2102 else Error("Streamer","cannot find Header Tree\n");
59fe9bd2 2103 fTreeE->GetEntry(0);
2ab0c725 2104
65fb704d 2105 gRandom = fRandom;
2ab0c725 2106 } else {
2107 AliRun::Class()->WriteBuffer(R__b, this);
2108 }
2109}
2110