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