]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliRun.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding convention corrections + few minor bug fixes
[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$
aee8290b 18Revision 1.38 2000/06/20 13:05:45 fca
19Writing down the TREE headers before job starts
20
51e0e89d 21Revision 1.37 2000/06/09 20:05:11 morsch
22Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat
23
f1b9d7c3 24Revision 1.36 2000/06/08 14:03:58 hristov
25Only one initializer for a default argument
26
d33c0226 27Revision 1.35 2000/06/07 10:13:14 hristov
28Delete only existent objects.
29
d2ecd553 30Revision 1.34 2000/05/18 10:45:38 fca
31Delete Particle Factory properly
32
c222d2b0 33Revision 1.33 2000/05/16 13:10:40 fca
34New method IsNewTrack and fix for a problem in Father-Daughter relations
35
a01a8b12 36Revision 1.32 2000/04/27 10:38:21 fca
37Correct termination of Lego Run and introduce Lego getter in AliRun
38
838edcaf 39Revision 1.31 2000/04/26 10:17:32 fca
40Changes in Lego for G4 compatibility
41
dffd31ef 42Revision 1.30 2000/04/18 19:11:40 fca
43Introduce variable Config.C function signature
44
45189757 45Revision 1.29 2000/04/07 11:12:34 fca
46G4 compatibility changes
47
875c717b 48Revision 1.28 2000/04/05 06:51:06 fca
49Workaround for an HP compiler problem
50
5eb58812 51Revision 1.27 2000/03/22 18:08:07 fca
52Rationalisation of the virtual MC interfaces
53
80762cb1 54Revision 1.26 2000/03/22 13:42:26 fca
55SetGenerator does not replace an existing generator, ResetGenerator does
56
ee1dd322 57Revision 1.25 2000/02/23 16:25:22 fca
58AliVMC and AliGeant3 classes introduced
59ReadEuclid moved from AliRun to AliModule
60
b13db077 61Revision 1.24 2000/01/19 17:17:20 fca
62Introducing a list of lists of hits -- more hits allowed for detector now
63
1cedd08a 64Revision 1.23 1999/12/03 11:14:31 fca
65Fixing previous wrong checking
66
00719c1b 67Revision 1.21 1999/11/25 10:40:08 fca
68Fixing daughters information also in primary tracks
69
ae23d366 70Revision 1.20 1999/10/04 18:08:49 fca
71Adding protection against inconsistent Euclid files
72
3fcc96a1 73Revision 1.19 1999/09/29 07:50:40 fca
74Introduction of the Copyright and cvs Log
75
99d554c8 76*/
77
fe4da5cc 78///////////////////////////////////////////////////////////////////////////////
79// //
80// Control class for Alice C++ //
81// Only one single instance of this class exists. //
82// The object is created in main program aliroot //
83// and is pointed by the global gAlice. //
84// //
8494b010 85// -Supports the list of all Alice Detectors (fModules). //
fe4da5cc 86// -Supports the list of particles (fParticles). //
87// -Supports the Trees. //
88// -Supports the geometry. //
89// -Supports the event display. //
90//Begin_Html
91/*
1439f98e 92<img src="picts/AliRunClass.gif">
fe4da5cc 93*/
94//End_Html
95//Begin_Html
96/*
1439f98e 97<img src="picts/alirun.gif">
fe4da5cc 98*/
99//End_Html
100// //
101///////////////////////////////////////////////////////////////////////////////
102
103#include <TFile.h>
104#include <TRandom.h>
105#include <TBRIK.h>
106#include <TNode.h>
fe4da5cc 107#include <TCint.h>
108#include <TSystem.h>
a8f1fb7c 109#include <TObjectTable.h>
fe4da5cc 110
1578254f 111#include "TParticle.h"
fe4da5cc 112#include "AliRun.h"
fe4da5cc 113#include "AliDisplay.h"
875c717b 114#include "AliMC.h"
dffd31ef 115#include "AliLego.h"
aee8290b 116#include "AliMagFC.h"
117#include "AliMagFCM.h"
118#include "AliMagFDM.h"
fe4da5cc 119
fe4da5cc 120#include <stdlib.h>
121#include <stdio.h>
122#include <string.h>
123
124AliRun *gAlice;
125
aee8290b 126static AliHeader *gAliHeader;
fe4da5cc 127
fe4da5cc 128ClassImp(AliRun)
129
130//_____________________________________________________________________________
131AliRun::AliRun()
132{
133 //
134 // Default constructor for AliRun
135 //
aee8290b 136 gAliHeader=&fHeader;
fe4da5cc 137 fRun = 0;
138 fEvent = 0;
139 fCurrent = -1;
8494b010 140 fModules = 0;
fe4da5cc 141 fGenerator = 0;
142 fTreeD = 0;
143 fTreeK = 0;
144 fTreeH = 0;
145 fTreeE = 0;
146 fTreeR = 0;
147 fParticles = 0;
148 fGeometry = 0;
149 fDisplay = 0;
150 fField = 0;
875c717b 151 fMC = 0;
fe4da5cc 152 fNdets = 0;
153 fImedia = 0;
154 fTrRmax = 1.e10;
155 fTrZmax = 1.e10;
fe4da5cc 156 fInitDone = kFALSE;
157 fLego = 0;
1578254f 158 fPDGDB = 0; //Particle factory object!
1cedd08a 159 fHitLists = 0;
45189757 160 fConfigFunction = "\0";
fe4da5cc 161}
162
163//_____________________________________________________________________________
164AliRun::AliRun(const char *name, const char *title)
165 : TNamed(name,title)
166{
167 //
168 // Constructor for the main processor.
169 // Creates the geometry
170 // Creates the list of Detectors.
171 // Creates the list of particles.
172 //
173 Int_t i;
174
175 gAlice = this;
176 fTreeD = 0;
177 fTreeK = 0;
178 fTreeH = 0;
179 fTreeE = 0;
180 fTreeR = 0;
181 fTrRmax = 1.e10;
182 fTrZmax = 1.e10;
1141f8e4 183 fGenerator = 0;
fe4da5cc 184 fInitDone = kFALSE;
185 fLego = 0;
186 fField = 0;
45189757 187 fConfigFunction = "Config();";
fe4da5cc 188
189 gROOT->GetListOfBrowsables()->Add(this,name);
190 //
191 // create the support list for the various Detectors
8494b010 192 fModules = new TObjArray(77);
fe4da5cc 193 //
194 // Create the TNode geometry for the event display
195
196 BuildSimpleGeometry();
197
198
199 fNtrack=0;
200 fHgwmk=0;
201 fCurrent=-1;
aee8290b 202 gAliHeader=&fHeader;
fe4da5cc 203 fRun = 0;
204 fEvent = 0;
205 //
206 // Create the particle stack
1578254f 207 fParticles = new TClonesArray("TParticle",100);
fe4da5cc 208
209 fDisplay = 0;
210 //
211 // Create default mag field
212 SetField();
213 //
875c717b 214 fMC = gMC;
fe4da5cc 215 //
216 // Prepare the tracking medium lists
217 fImedia = new TArrayI(1000);
218 for(i=0;i<1000;i++) (*fImedia)[i]=-99;
1578254f 219 //
220 // Make particles
221 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
1cedd08a 222 //
223 // Create HitLists list
224 fHitLists = new TList();
fe4da5cc 225}
226
aee8290b 227//_____________________________________________________________________________
228AliRun::AliRun(const AliRun &run)
229{
230 //
231 // Copy constructor
232 //
233 run.Copy(*this);
234}
235
236
fe4da5cc 237//_____________________________________________________________________________
238AliRun::~AliRun()
239{
240 //
241 // Defaullt AliRun destructor
242 //
fe4da5cc 243 delete fImedia;
244 delete fField;
875c717b 245 delete fMC;
fe4da5cc 246 delete fGeometry;
247 delete fDisplay;
248 delete fGenerator;
249 delete fLego;
250 delete fTreeD;
251 delete fTreeK;
252 delete fTreeH;
253 delete fTreeE;
254 delete fTreeR;
8494b010 255 if (fModules) {
256 fModules->Delete();
257 delete fModules;
fe4da5cc 258 }
259 if (fParticles) {
260 fParticles->Delete();
261 delete fParticles;
262 }
1cedd08a 263 delete fHitLists;
c222d2b0 264 delete fPDGDB;
fe4da5cc 265}
266
267//_____________________________________________________________________________
268void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
269{
270 //
271 // Add a hit to detector id
272 //
8494b010 273 TObjArray &dets = *fModules;
274 if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
fe4da5cc 275}
276
277//_____________________________________________________________________________
278void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
279{
280 //
281 // Add digit to detector id
282 //
8494b010 283 TObjArray &dets = *fModules;
284 if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
fe4da5cc 285}
286
287//_____________________________________________________________________________
288void AliRun::Browse(TBrowser *b)
289{
290 //
291 // Called when the item "Run" is clicked on the left pane
292 // of the Root browser.
293 // It displays the Root Trees and all detectors.
294 //
295 if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
296 if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
297 if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
298 if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
299 if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
300
8494b010 301 TIter next(fModules);
302 AliModule *detector;
303 while((detector = (AliModule*)next())) {
fe4da5cc 304 b->Add(detector,detector->GetName());
305 }
306}
307
308//_____________________________________________________________________________
309void AliRun::Build()
310{
311 //
312 // Initialize Alice geometry
313 // Dummy routine
314 //
315}
316
317//_____________________________________________________________________________
318void AliRun::BuildSimpleGeometry()
319{
320 //
321 // Create a simple TNode geometry used by Root display engine
322 //
323 // Initialise geometry
324 //
325 fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
326 new TMaterial("void","Vacuum",0,0,0); //Everything is void
327 TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
328 brik->SetVisibility(0);
329 new TNode("alice","alice","S_alice");
330}
331
332//_____________________________________________________________________________
333void AliRun::CleanDetectors()
334{
335 //
336 // Clean Detectors at the end of event
337 //
8494b010 338 TIter next(fModules);
339 AliModule *detector;
340 while((detector = (AliModule*)next())) {
fe4da5cc 341 detector->FinishEvent();
342 }
343}
344
345//_____________________________________________________________________________
346void AliRun::CleanParents()
347{
348 //
349 // Clean Particles stack.
1578254f 350 // Set parent/daughter relations
fe4da5cc 351 //
352 TClonesArray &particles = *(gAlice->Particles());
1578254f 353 TParticle *part;
fe4da5cc 354 int i;
355 for(i=0; i<fNtrack; i++) {
1578254f 356 part = (TParticle *)particles.UncheckedAt(i);
aee8290b 357 if(!part->TestBit(kDaughtersBit)) {
1578254f 358 part->SetFirstDaughter(-1);
359 part->SetLastDaughter(-1);
fe4da5cc 360 }
361 }
362}
363
364//_____________________________________________________________________________
365Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
366{
367 //
368 // Return the distance from the mouse to the AliRun object
369 // Dummy routine
370 //
371 return 9999;
372}
373
374//_____________________________________________________________________________
375void AliRun::DumpPart (Int_t i)
376{
377 //
378 // Dumps particle i in the stack
379 //
380 TClonesArray &particles = *fParticles;
1578254f 381 ((TParticle*) particles[i])->Print();
fe4da5cc 382}
383
384//_____________________________________________________________________________
385void AliRun::DumpPStack ()
386{
387 //
388 // Dumps the particle stack
389 //
390 TClonesArray &particles = *fParticles;
391 printf(
392 "\n\n=======================================================================\n");
393 for (Int_t i=0;i<fNtrack;i++)
394 {
1578254f 395 printf("-> %d ",i); ((TParticle*) particles[i])->Print();
fe4da5cc 396 printf("--------------------------------------------------------------\n");
397 }
398 printf(
399 "\n=======================================================================\n\n");
400}
401
402//_____________________________________________________________________________
403void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
404 Float_t maxField, char* filename)
405{
406 //
407 // Set magnetic field parameters
408 // type Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
409 // version Magnetic field map version (only 1 active now)
410 // scale Scale factor for the magnetic field
411 // maxField Maximum value for the magnetic field
412
413 //
414 // --- Sanity check on mag field flags
415 if(type<0 || type > 2) {
23370b7a 416 Warning("SetField",
417 "Invalid magnetic field flag: %5d; Helix tracking chosen instead\n"
fe4da5cc 418 ,type);
419 type=2;
420 }
421 if(fField) delete fField;
422 if(version==1) {
423 fField = new AliMagFC("Map1"," ",type,version,scale,maxField);
f1b9d7c3 424 } else if(version<=2) {
fe4da5cc 425 fField = new AliMagFCM("Map2-3",filename,type,version,scale,maxField);
426 fField->ReadField();
f1b9d7c3 427 } else if(version==3) {
428 fField = new AliMagFDM("Map4",filename,type,version,scale,maxField);
429 fField->ReadField();
fe4da5cc 430 } else {
23370b7a 431 Warning("SetField","Invalid map %d\n",version);
fe4da5cc 432 }
433}
434
435//_____________________________________________________________________________
436void AliRun::FillTree()
437{
438 //
439 // Fills all AliRun TTrees
440 //
441 if (fTreeK) fTreeK->Fill();
442 if (fTreeH) fTreeH->Fill();
443 if (fTreeD) fTreeD->Fill();
444 if (fTreeR) fTreeR->Fill();
445}
446
447//_____________________________________________________________________________
448void AliRun::FinishPrimary()
449{
450 //
451 // Called at the end of each primary track
452 //
453
6c9704e6 454 // static Int_t count=0;
455 // const Int_t times=10;
fe4da5cc 456 // This primary is finished, purify stack
80762cb1 457 PurifyKine();
fe4da5cc 458
459 // Write out hits if any
460 if (gAlice->TreeH()) {
461 gAlice->TreeH()->Fill();
462 }
463
464 // Reset Hits info
465 gAlice->ResetHits();
a8f1fb7c 466
467 //
468 // if(++count%times==1) gObjectTable->Print();
fe4da5cc 469}
470
471//_____________________________________________________________________________
472void AliRun::FinishEvent()
473{
474 //
475 // Called at the end of the event.
476 //
7fb01480 477
dffd31ef 478 //
479 if(fLego) fLego->FinishEvent();
480
fe4da5cc 481 //Update the energy deposit tables
482 Int_t i;
875c717b 483 for(i=0;i<fEventEnergy.GetSize();i++) {
484 fSummEnergy[i]+=fEventEnergy[i];
485 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
fe4da5cc 486 }
875c717b 487 fEventEnergy.Reset();
fe4da5cc 488
489 // Clean detector information
490 CleanDetectors();
491
492 // Write out the kinematics
493 if (fTreeK) {
494 CleanParents();
495 fTreeK->Fill();
496 }
497
498 // Write out the digits
499 if (fTreeD) {
500 fTreeD->Fill();
501 ResetDigits();
502 }
503
504 // Write out reconstructed clusters
505 if (fTreeR) {
506 fTreeR->Fill();
507 }
508
509 // Write out the event Header information
510 if (fTreeE) fTreeE->Fill();
511
512 // Reset stack info
513 ResetStack();
514
515 // Write Tree headers
59fe9bd2 516 // Int_t ievent = fHeader.GetEvent();
517 // char hname[30];
518 // sprintf(hname,"TreeK%d",ievent);
51e0e89d 519 if (fTreeK) fTreeK->Write(0,TObject::kOverwrite);
59fe9bd2 520 // sprintf(hname,"TreeH%d",ievent);
51e0e89d 521 if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
59fe9bd2 522 // sprintf(hname,"TreeD%d",ievent);
51e0e89d 523 if (fTreeD) fTreeD->Write(0,TObject::kOverwrite);
59fe9bd2 524 // sprintf(hname,"TreeR%d",ievent);
51e0e89d 525 if (fTreeR) fTreeR->Write(0,TObject::kOverwrite);
875c717b 526
527 ++fEvent;
fe4da5cc 528}
529
530//_____________________________________________________________________________
531void AliRun::FinishRun()
532{
533 //
534 // Called at the end of the run.
535 //
536
dffd31ef 537 //
538 if(fLego) fLego->FinishRun();
539
fe4da5cc 540 // Clean detector information
8494b010 541 TIter next(fModules);
542 AliModule *detector;
543 while((detector = (AliModule*)next())) {
fe4da5cc 544 detector->FinishRun();
545 }
546
547 //Output energy summary tables
548 EnergySummary();
549
550 // file is retrieved from whatever tree
aee8290b 551 TFile *file = 0;
552 if (fTreeK) file = fTreeK->GetCurrentFile();
553 if ((!file) && (fTreeH)) file = fTreeH->GetCurrentFile();
554 if ((!file) && (fTreeD)) file = fTreeD->GetCurrentFile();
555 if ((!file) && (fTreeE)) file = fTreeE->GetCurrentFile();
556 if( NULL==file ) {
fe4da5cc 557 Error("FinishRun","There isn't root file!");
558 exit(1);
559 }
aee8290b 560 file->cd();
51e0e89d 561 fTreeE->Write(0,TObject::kOverwrite);
fe4da5cc 562
563 // Clean tree information
d2ecd553 564 if (fTreeK) {
565 delete fTreeK; fTreeK = 0;
566 }
567 if (fTreeH) {
568 delete fTreeH; fTreeH = 0;
569 }
570 if (fTreeD) {
571 delete fTreeD; fTreeD = 0;
572 }
573 if (fTreeR) {
574 delete fTreeR; fTreeR = 0;
575 }
576 if (fTreeE) {
577 delete fTreeE; fTreeE = 0;
578 }
fe4da5cc 579
580 // Write AliRun info and all detectors parameters
581 Write();
582
583 // Close output file
aee8290b 584 file->Write();
fe4da5cc 585}
586
587//_____________________________________________________________________________
588void AliRun::FlagTrack(Int_t track)
589{
590 //
591 // Flags a track and all its family tree to be kept
592 //
593 int curr;
1578254f 594 TParticle *particle;
fe4da5cc 595
596 curr=track;
597 while(1) {
1578254f 598 particle=(TParticle*)fParticles->UncheckedAt(curr);
fe4da5cc 599
600 // If the particle is flagged the three from here upward is saved already
aee8290b 601 if(particle->TestBit(kKeepBit)) return;
fe4da5cc 602
603 // Save this particle
aee8290b 604 particle->SetBit(kKeepBit);
fe4da5cc 605
606 // Move to father if any
1578254f 607 if((curr=particle->GetFirstMother())==-1) return;
fe4da5cc 608 }
609}
610
611//_____________________________________________________________________________
612void AliRun::EnergySummary()
613{
614 //
615 // Print summary of deposited energy
616 //
617
fe4da5cc 618 Int_t ndep=0;
619 Float_t edtot=0;
620 Float_t ed, ed2;
621 Int_t kn, i, left, j, id;
aee8290b 622 const Float_t kzero=0;
fe4da5cc 623 Int_t ievent=fHeader.GetEvent()+1;
624 //
625 // Energy loss information
626 if(ievent) {
627 printf("***************** Energy Loss Information per event (GEV) *****************\n");
875c717b 628 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
629 ed=fSummEnergy[kn];
fe4da5cc 630 if(ed>0) {
875c717b 631 fEventEnergy[ndep]=kn;
fe4da5cc 632 if(ievent>1) {
633 ed=ed/ievent;
875c717b 634 ed2=fSum2Energy[kn];
fe4da5cc 635 ed2=ed2/ievent;
aee8290b 636 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
fe4da5cc 637 } else
638 ed2=99;
875c717b 639 fSummEnergy[ndep]=ed;
aee8290b 640 fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
fe4da5cc 641 edtot+=ed;
642 ndep++;
643 }
644 }
645 for(kn=0;kn<(ndep-1)/3+1;kn++) {
646 left=ndep-kn*3;
647 for(i=0;i<(3<left?3:left);i++) {
648 j=kn*3+i;
875c717b 649 id=Int_t (fEventEnergy[j]+0.1);
650 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
fe4da5cc 651 }
652 printf("\n");
653 }
654 //
655 // Relative energy loss in different detectors
656 printf("******************** Relative Energy Loss per event ********************\n");
657 printf("Total energy loss per event %10.3f GeV\n",edtot);
658 for(kn=0;kn<(ndep-1)/5+1;kn++) {
659 left=ndep-kn*5;
660 for(i=0;i<(5<left?5:left);i++) {
661 j=kn*5+i;
875c717b 662 id=Int_t (fEventEnergy[j]+0.1);
663 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
fe4da5cc 664 }
665 printf("\n");
666 }
667 for(kn=0;kn<75;kn++) printf("*");
668 printf("\n");
669 }
670 //
671 // Reset the TArray's
875c717b 672 // fEventEnergy.Set(0);
673 // fSummEnergy.Set(0);
674 // fSum2Energy.Set(0);
fe4da5cc 675}
676
677//_____________________________________________________________________________
8494b010 678AliModule *AliRun::GetModule(const char *name)
fe4da5cc 679{
680 //
681 // Return pointer to detector from name
682 //
8494b010 683 return (AliModule*)fModules->FindObject(name);
fe4da5cc 684}
685
a68348e9 686//_____________________________________________________________________________
687AliDetector *AliRun::GetDetector(const char *name)
688{
689 //
690 // Return pointer to detector from name
691 //
692 return (AliDetector*)fModules->FindObject(name);
693}
694
fe4da5cc 695//_____________________________________________________________________________
8494b010 696Int_t AliRun::GetModuleID(const char *name)
fe4da5cc 697{
698 //
699 // Return galice internal detector identifier from name
700 //
23370b7a 701 Int_t i=-1;
702 TObject *mod=fModules->FindObject(name);
703 if(mod) i=fModules->IndexOf(mod);
704 return i;
fe4da5cc 705}
706
707//_____________________________________________________________________________
708Int_t AliRun::GetEvent(Int_t event)
709{
710 //
711 // Connect the Trees Kinematics and Hits for event # event
712 // Set branch addresses
713 //
fe4da5cc 714
715 // Reset existing structures
716 ResetStack();
717 ResetHits();
718 ResetDigits();
719
720 // Delete Trees already connected
721 if (fTreeK) delete fTreeK;
722 if (fTreeH) delete fTreeH;
723 if (fTreeD) delete fTreeD;
724 if (fTreeR) delete fTreeR;
59fe9bd2 725
726 // Get header from file
727 if(fTreeE) fTreeE->GetEntry(event);
728 else Error("GetEvent","Cannot file Header Tree\n");
fe4da5cc 729
730 // Get Kine Tree from file
731 char treeName[20];
732 sprintf(treeName,"TreeK%d",event);
733 fTreeK = (TTree*)gDirectory->Get(treeName);
734 if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticles);
23370b7a 735 else Error("GetEvent","cannot find Kine Tree for event:%d\n",event);
fe4da5cc 736
737 // Get Hits Tree header from file
738 sprintf(treeName,"TreeH%d",event);
739 fTreeH = (TTree*)gDirectory->Get(treeName);
740 if (!fTreeH) {
23370b7a 741 Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
fe4da5cc 742 }
743
744 // Get Digits Tree header from file
745 sprintf(treeName,"TreeD%d",event);
746 fTreeD = (TTree*)gDirectory->Get(treeName);
747 if (!fTreeD) {
07a68c1d 748 Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
fe4da5cc 749 }
750
751
752 // Get Reconstruct Tree header from file
753 sprintf(treeName,"TreeR%d",event);
754 fTreeR = (TTree*)gDirectory->Get(treeName);
755 if (!fTreeR) {
756 // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
757 }
758
759 // Set Trees branch addresses
8494b010 760 TIter next(fModules);
761 AliModule *detector;
762 while((detector = (AliModule*)next())) {
fe4da5cc 763 detector->SetTreeAddress();
764 }
765
766 if (fTreeK) fTreeK->GetEvent(0);
767 fNtrack = Int_t (fParticles->GetEntries());
768 return fNtrack;
769}
770
771//_____________________________________________________________________________
772TGeometry *AliRun::GetGeometry()
773{
774 //
775 // Import Alice geometry from current file
776 // Return pointer to geometry object
777 //
778 if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
779 //
780 // Unlink and relink nodes in detectors
781 // This is bad and there must be a better way...
782 //
fe4da5cc 783
8494b010 784 TIter next(fModules);
785 AliModule *detector;
786 while((detector = (AliModule*)next())) {
fe4da5cc 787 detector->SetTreeAddress();
788 TList *dnodes=detector->Nodes();
789 Int_t j;
790 TNode *node, *node1;
791 for ( j=0; j<dnodes->GetSize(); j++) {
792 node = (TNode*) dnodes->At(j);
52d0ab00 793 node1 = fGeometry->GetNode(node->GetName());
fe4da5cc 794 dnodes->Remove(node);
795 dnodes->AddAt(node1,j);
796 }
797 }
798 return fGeometry;
799}
800
801//_____________________________________________________________________________
802void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
803 Float_t &e, Float_t *vpos, Float_t *polar,
804 Float_t &tof)
805{
806 //
807 // Return next track from stack of particles
808 //
a8f1fb7c 809 TVector3 pol;
fe4da5cc 810 fCurrent=-1;
1578254f 811 TParticle *track;
fe4da5cc 812 for(Int_t i=fNtrack-1; i>=0; i--) {
1578254f 813 track=(TParticle*) fParticles->UncheckedAt(i);
aee8290b 814 if(!track->TestBit(kDoneBit)) {
fe4da5cc 815 //
816 // The track has not yet been processed
817 fCurrent=i;
1578254f 818 ipart=track->GetPdgCode();
819 pmom[0]=track->Px();
820 pmom[1]=track->Py();
821 pmom[2]=track->Pz();
822 e =track->Energy();
823 vpos[0]=track->Vx();
824 vpos[1]=track->Vy();
825 vpos[2]=track->Vz();
a8f1fb7c 826 track->GetPolarisation(pol);
827 polar[0]=pol.X();
828 polar[1]=pol.Y();
829 polar[2]=pol.Z();
1578254f 830 tof=track->T();
aee8290b 831 track->SetBit(kDoneBit);
fe4da5cc 832 break;
833 }
834 }
835 mtrack=fCurrent;
836 //
837 // stop and start timer when we start a primary track
838 Int_t nprimaries = fHeader.GetNprimary();
839 if (fCurrent >= nprimaries) return;
840 if (fCurrent < nprimaries-1) {
841 fTimer.Stop();
1578254f 842 track=(TParticle*) fParticles->UncheckedAt(fCurrent+1);
843 // track->SetProcessTime(fTimer.CpuTime());
fe4da5cc 844 }
845 fTimer.Start();
846}
847
848//_____________________________________________________________________________
849Int_t AliRun::GetPrimary(Int_t track)
850{
851 //
852 // return number of primary that has generated track
853 //
854 int current, parent;
1578254f 855 TParticle *part;
fe4da5cc 856 //
857 parent=track;
858 while (1) {
859 current=parent;
1578254f 860 part = (TParticle *)fParticles->UncheckedAt(current);
861 parent=part->GetFirstMother();
fe4da5cc 862 if(parent<0) return current;
863 }
864}
865
866//_____________________________________________________________________________
875c717b 867void AliRun::InitMC(const char *setup)
fe4da5cc 868{
869 //
870 // Initialize the Alice setup
871 //
872
51e0e89d 873 if(fInitDone) {
874 Warning("Init","Cannot initialise AliRun twice!\n");
875 return;
876 }
877
fe4da5cc 878 gROOT->LoadMacro(setup);
45189757 879 gInterpreter->ProcessLine(fConfigFunction.Data());
fe4da5cc 880
cfce8870 881 gMC->DefineParticles(); //Create standard MC particles
fe4da5cc 882
883 TObject *objfirst, *objlast;
884
23370b7a 885 fNdets = fModules->GetLast()+1;
886
fe4da5cc 887 //
875c717b 888 //=================Create Materials and geometry
889 gMC->Init();
890
8494b010 891 TIter next(fModules);
892 AliModule *detector;
893 while((detector = (AliModule*)next())) {
fe4da5cc 894 detector->SetTreeAddress();
895 objlast = gDirectory->GetList()->Last();
896
fe4da5cc 897 // Add Detector histograms in Detector list of histograms
898 if (objlast) objfirst = gDirectory->GetList()->After(objlast);
899 else objfirst = gDirectory->GetList()->First();
900 while (objfirst) {
901 detector->Histograms()->Add(objfirst);
902 objfirst = gDirectory->GetList()->After(objfirst);
903 }
904 }
905 SetTransPar(); //Read the cuts for all materials
906
907 MediaTable(); //Build the special IMEDIA table
908
fe4da5cc 909 //Initialise geometry deposition table
875c717b 910 fEventEnergy.Set(gMC->NofVolumes()+1);
911 fSummEnergy.Set(gMC->NofVolumes()+1);
912 fSum2Energy.Set(gMC->NofVolumes()+1);
fe4da5cc 913
fe4da5cc 914 //Compute cross-sections
875c717b 915 gMC->BuildPhysics();
fe4da5cc 916
917 //Write Geometry object to current file.
918 fGeometry->Write();
919
920 fInitDone = kTRUE;
51e0e89d 921
922 //
923 // Save stuff at the beginning of the file to avoid file corruption
924 Write();
fe4da5cc 925}
926
927//_____________________________________________________________________________
928void AliRun::MediaTable()
929{
930 //
931 // Built media table to get from the media number to
932 // the detector id
933 //
ad51aeb0 934 Int_t kz, nz, idt, lz, i, k, ind;
935 // Int_t ibeg;
fe4da5cc 936 TObjArray &dets = *gAlice->Detectors();
8494b010 937 AliModule *det;
fe4da5cc 938 //
939 // For all detectors
940 for (kz=0;kz<fNdets;kz++) {
941 // If detector is defined
8494b010 942 if((det=(AliModule*) dets[kz])) {
ad51aeb0 943 TArrayI &idtmed = *(det->GetIdtmed());
944 for(nz=0;nz<100;nz++) {
fe4da5cc 945 // Find max and min material number
ad51aeb0 946 if((idt=idtmed[nz])) {
fe4da5cc 947 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
948 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
949 }
950 }
951 if(det->LoMedium() > det->HiMedium()) {
952 det->LoMedium() = 0;
953 det->HiMedium() = 0;
954 } else {
955 if(det->HiMedium() > fImedia->GetSize()) {
ad51aeb0 956 Error("MediaTable","Increase fImedia from %d to %d",
957 fImedia->GetSize(),det->HiMedium());
fe4da5cc 958 return;
959 }
960 // Tag all materials in rage as belonging to detector kz
961 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
962 (*fImedia)[lz]=kz;
963 }
964 }
965 }
966 }
967 //
968 // Print summary table
969 printf(" Traking media ranges:\n");
970 for(i=0;i<(fNdets-1)/6+1;i++) {
971 for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
972 ind=i*6+k;
8494b010 973 det=(AliModule*)dets[ind];
fe4da5cc 974 if(det)
975 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
976 det->HiMedium());
977 else
978 printf(" %6s: %3d -> %3d;","NULL",0,0);
979 }
980 printf("\n");
981 }
982}
983
984//____________________________________________________________________________
985void AliRun::SetGenerator(AliGenerator *generator)
ee1dd322 986{
987 //
988 // Load the event generator
989 //
990 if(!fGenerator) fGenerator = generator;
991}
992
993//____________________________________________________________________________
994void AliRun::ResetGenerator(AliGenerator *generator)
fe4da5cc 995{
996 //
997 // Load the event generator
998 //
b13db077 999 if(fGenerator)
838edcaf 1000 if(generator)
1001 Warning("ResetGenerator","Replacing generator %s with %s\n",
1002 fGenerator->GetName(),generator->GetName());
1003 else
1004 Warning("ResetGenerator","Replacing generator %s with NULL\n",
1005 fGenerator->GetName());
b13db077 1006 fGenerator = generator;
fe4da5cc 1007}
1008
1009//____________________________________________________________________________
1010void AliRun::SetTransPar(char* filename)
1011{
1012 //
1013 // Read filename to set the transport parameters
1014 //
1015
fe4da5cc 1016
aee8290b 1017 const Int_t kncuts=10;
1018 const Int_t knflags=11;
1019 const Int_t knpars=kncuts+knflags;
1020 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
fe4da5cc 1021 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
1022 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
1023 "MULS","PAIR","PHOT","RAYL"};
1024 char line[256];
ad51aeb0 1025 char detName[7];
fe4da5cc 1026 char* filtmp;
aee8290b 1027 Float_t cut[kncuts];
1028 Int_t flag[knflags];
fe4da5cc 1029 Int_t i, itmed, iret, ktmed, kz;
1030 FILE *lun;
1031 //
1032 // See whether the file is there
1033 filtmp=gSystem->ExpandPathName(filename);
1034 lun=fopen(filtmp,"r");
1035 delete [] filtmp;
1036 if(!lun) {
ad51aeb0 1037 Warning("SetTransPar","File %s does not exist!\n",filename);
fe4da5cc 1038 return;
1039 }
1040 //
1041 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1042 printf(" *%59s\n","*");
1043 printf(" * Please check carefully what you are doing!%10s\n","*");
1044 printf(" *%59s\n","*");
1045 //
1046 while(1) {
1047 // Initialise cuts and flags
aee8290b 1048 for(i=0;i<kncuts;i++) cut[i]=-99;
1049 for(i=0;i<knflags;i++) flag[i]=-99;
fe4da5cc 1050 itmed=0;
1051 for(i=0;i<256;i++) line[i]='\0';
1052 // Read up to the end of line excluded
1053 iret=fscanf(lun,"%[^\n]",line);
1054 if(iret<0) {
1055 //End of file
1056 fclose(lun);
1057 printf(" *%59s\n","*");
1058 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1059 return;
1060 }
1061 // Read the end of line
1062 fscanf(lun,"%*c");
1063 if(!iret) continue;
1064 if(line[0]=='*') continue;
1065 // Read the numbers
ad51aeb0 1066 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",
1067 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1068 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1069 &flag[8],&flag[9],&flag[10]);
fe4da5cc 1070 if(!iret) continue;
1071 if(iret<0) {
1072 //reading error
ad51aeb0 1073 Warning("SetTransPar","Error reading file %s\n",filename);
fe4da5cc 1074 continue;
1075 }
ad51aeb0 1076 // Check that the module exist
1077 AliModule *mod = GetModule(detName);
1078 if(mod) {
1079 // Get the array of media numbers
1080 TArrayI &idtmed = *mod->GetIdtmed();
1081 // Check that the tracking medium code is valid
1082 if(0<=itmed && itmed < 100) {
1083 ktmed=idtmed[itmed];
1084 if(!ktmed) {
1085 Warning("SetTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
1086 continue;
fe4da5cc 1087 }
ad51aeb0 1088 // Set energy thresholds
aee8290b 1089 for(kz=0;kz<kncuts;kz++) {
ad51aeb0 1090 if(cut[kz]>=0) {
23370b7a 1091 printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
aee8290b 1092 kpars[kz],cut[kz],itmed,mod->GetName());
1093 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
ad51aeb0 1094 }
fe4da5cc 1095 }
ad51aeb0 1096 // Set transport mechanisms
aee8290b 1097 for(kz=0;kz<knflags;kz++) {
ad51aeb0 1098 if(flag[kz]>=0) {
1099 printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
aee8290b 1100 kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
1101 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
ad51aeb0 1102 }
1103 }
1104 } else {
1105 Warning("SetTransPar","Invalid medium code %d *\n",itmed);
1106 continue;
fe4da5cc 1107 }
1108 } else {
ad51aeb0 1109 Warning("SetTransPar","Module %s not present\n",detName);
fe4da5cc 1110 continue;
1111 }
1112 }
1113}
1114
1115//_____________________________________________________________________________
1116void AliRun::MakeTree(Option_t *option)
1117{
1118 //
1119 // Create the ROOT trees
1120 // Loop on all detectors to create the Root branch (if any)
1121 //
1122
b13db077 1123 char hname[30];
fe4da5cc 1124 //
1125 // Analyse options
aee8290b 1126 char *oK = strstr(option,"K");
1127 char *oH = strstr(option,"H");
1128 char *oE = strstr(option,"E");
1129 char *oD = strstr(option,"D");
1130 char *oR = strstr(option,"R");
fe4da5cc 1131 //
aee8290b 1132 if (oK && !fTreeK) {
b13db077 1133 sprintf(hname,"TreeK%d",fEvent);
1134 fTreeK = new TTree(hname,"Kinematics");
1135 // Create a branch for particles
1136 fTreeK->Branch("Particles",&fParticles,4000);
51e0e89d 1137 fTreeK->Write();
b13db077 1138 }
aee8290b 1139 if (oH && !fTreeH) {
b13db077 1140 sprintf(hname,"TreeH%d",fEvent);
1141 fTreeH = new TTree(hname,"Hits");
1142 fTreeH->SetAutoSave(1000000000); //no autosave
51e0e89d 1143 fTreeH->Write();
b13db077 1144 }
aee8290b 1145 if (oD && !fTreeD) {
b13db077 1146 sprintf(hname,"TreeD%d",fEvent);
1147 fTreeD = new TTree(hname,"Digits");
51e0e89d 1148 fTreeD->Write();
b13db077 1149 }
aee8290b 1150 if (oR && !fTreeR) {
b13db077 1151 sprintf(hname,"TreeR%d",fEvent);
1152 fTreeR = new TTree(hname,"Reconstruction");
51e0e89d 1153 fTreeR->Write();
b13db077 1154 }
aee8290b 1155 if (oE && !fTreeE) {
b13db077 1156 fTreeE = new TTree("TE","Header");
1157 // Create a branch for Header
aee8290b 1158 fTreeE->Branch("Header","AliHeader",&gAliHeader,4000);
51e0e89d 1159 fTreeE->Write();
b13db077 1160 }
fe4da5cc 1161 //
1162 // Create a branch for hits/digits for each detector
1163 // Each branch is a TClonesArray. Each data member of the Hits classes
1164 // will be in turn a subbranch of the detector master branch
8494b010 1165 TIter next(fModules);
1166 AliModule *detector;
1167 while((detector = (AliModule*)next())) {
aee8290b 1168 if (oH || oD || oR) detector->MakeBranch(option);
fe4da5cc 1169 }
fe4da5cc 1170}
1171
1172//_____________________________________________________________________________
1173Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks)
1174{
1175 //
1176 // PurifyKine with external parameters
1177 //
1178 fHgwmk = lastSavedTrack;
1179 fNtrack = nofTracks;
1180 PurifyKine();
1181 return fHgwmk;
1182}
1183
1184//_____________________________________________________________________________
1185void AliRun::PurifyKine()
1186{
1187 //
1188 // Compress kinematic tree keeping only flagged particles
1189 // and renaming the particle id's in all the hits
1190 //
1191 TClonesArray &particles = *fParticles;
1192 int nkeep=fHgwmk+1, parent, i;
1578254f 1193 TParticle *part, *partnew, *father;
fe4da5cc 1194 int *map = new int[particles.GetEntries()];
1195
1196 // Save in Header total number of tracks before compression
1197 fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk);
1198
a01a8b12 1199 // First pass, invalid Daughter information
fe4da5cc 1200 for(i=0; i<fNtrack; i++) {
a01a8b12 1201 // Preset map, to be removed later
1202 if(i<=fHgwmk) map[i]=i ; else map[i] = -99;
aee8290b 1203 ((TParticle *)particles.UncheckedAt(i))->ResetBit(kDaughtersBit);
a01a8b12 1204 }
fe4da5cc 1205 // Second pass, build map between old and new numbering
1206 for(i=fHgwmk+1; i<fNtrack; i++) {
1578254f 1207 part = (TParticle *)particles.UncheckedAt(i);
aee8290b 1208 if(part->TestBit(kKeepBit)) {
fe4da5cc 1209
1210 // This particle has to be kept
1211 map[i]=nkeep;
1212 if(i!=nkeep) {
1213
1214 // Old and new are different, have to copy
1578254f 1215 partnew = (TParticle *)particles.UncheckedAt(nkeep);
5eb58812 1216 // Change due to a bug in the HP compiler
1217 // *partnew = *part;
1218 memcpy(partnew,part,sizeof(TParticle));
fe4da5cc 1219 } else partnew = part;
1220
1221 // as the parent is always *before*, it must be already
1222 // in place. This is what we are checking anyway!
1578254f 1223 if((parent=partnew->GetFirstMother())>fHgwmk) {
fe4da5cc 1224 if(map[parent]==-99) printf("map[%d] = -99!\n",parent);
1578254f 1225 partnew->SetFirstMother(map[parent]);
fe4da5cc 1226 }
1227 nkeep++;
1228 }
1229 }
1230 fNtrack=nkeep;
1231
1578254f 1232 // Fix daughters information
ae23d366 1233 for (i=0; i<fNtrack; i++) {
1578254f 1234 part = (TParticle *)particles.UncheckedAt(i);
1235 parent = part->GetFirstMother();
ae23d366 1236 if(parent>=0) {
1237 father = (TParticle *)particles.UncheckedAt(parent);
aee8290b 1238 if(father->TestBit(kDaughtersBit)) {
fe4da5cc 1239
ae23d366 1240 if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
1241 if(i>father->GetLastDaughter()) father->SetLastDaughter(i);
1242 } else {
1243 // Iitialise daughters info for first pass
1244 father->SetFirstDaughter(i);
1245 father->SetLastDaughter(i);
aee8290b 1246 father->SetBit(kDaughtersBit);
ae23d366 1247 }
fe4da5cc 1248 }
1249 }
1250
1cedd08a 1251#ifdef old
fe4da5cc 1252 // Now loop on all detectors and reset the hits
b13db077 1253 AliHit *OneHit;
8494b010 1254 TIter next(fModules);
1255 AliModule *detector;
1256 while((detector = (AliModule*)next())) {
fe4da5cc 1257 if (!detector->Hits()) continue;
1258 TClonesArray &vHits=*(detector->Hits());
1259 if(vHits.GetEntries() != detector->GetNhits())
1260 printf("vHits.GetEntries()!=detector->GetNhits(): %d != %d\n",
1261 vHits.GetEntries(),detector->GetNhits());
1262 for (i=0; i<detector->GetNhits(); i++) {
1263 OneHit = (AliHit *)vHits.UncheckedAt(i);
1264 OneHit->SetTrack(map[OneHit->GetTrack()]);
1265 }
1266 }
1cedd08a 1267#else
1268
1269 // Now loop on all registered hit lists
1270 TIter next(fHitLists);
1271 TCollection *hitList;
1272 while((hitList = (TCollection*)next())) {
1273 TIter nexthit(hitList);
1274 AliHit *hit;
1275 while((hit = (AliHit*)nexthit())) {
1276 hit->SetTrack(map[hit->GetTrack()]);
1277 }
1278 }
1279#endif
fe4da5cc 1280
1281 fHgwmk=nkeep-1;
1282 particles.SetLast(fHgwmk);
1283 delete [] map;
1284}
1285
1286//_____________________________________________________________________________
dffd31ef 1287void AliRun::BeginEvent()
fe4da5cc 1288{
1289 //
1290 // Reset all Detectors & kinematics & trees
1291 //
59fe9bd2 1292 char hname[30];
dffd31ef 1293 //
1294
1295 //
1296 if(fLego) {
1297 fLego->BeginEvent();
1298 return;
1299 }
1300
59fe9bd2 1301 //
fe4da5cc 1302 ResetStack();
1303 ResetHits();
1304 ResetDigits();
1305
1306 // Initialise event header
875c717b 1307 fHeader.Reset(fRun,fEvent);
fe4da5cc 1308
59fe9bd2 1309 if(fTreeK) {
1310 fTreeK->Reset();
875c717b 1311 sprintf(hname,"TreeK%d",fEvent);
59fe9bd2 1312 fTreeK->SetName(hname);
1313 }
1314 if(fTreeH) {
1315 fTreeH->Reset();
875c717b 1316 sprintf(hname,"TreeH%d",fEvent);
59fe9bd2 1317 fTreeH->SetName(hname);
1318 }
1319 if(fTreeD) {
1320 fTreeD->Reset();
875c717b 1321 sprintf(hname,"TreeD%d",fEvent);
59fe9bd2 1322 fTreeD->SetName(hname);
1323 }
1324 if(fTreeR) {
1325 fTreeR->Reset();
875c717b 1326 sprintf(hname,"TreeR%d",fEvent);
59fe9bd2 1327 fTreeR->SetName(hname);
1328 }
fe4da5cc 1329}
1330
1331//_____________________________________________________________________________
1332void AliRun::ResetDigits()
1333{
1334 //
1335 // Reset all Detectors digits
1336 //
8494b010 1337 TIter next(fModules);
1338 AliModule *detector;
1339 while((detector = (AliModule*)next())) {
fe4da5cc 1340 detector->ResetDigits();
1341 }
1342}
1343
1344//_____________________________________________________________________________
1345void AliRun::ResetHits()
1346{
1347 //
1348 // Reset all Detectors hits
1349 //
8494b010 1350 TIter next(fModules);
1351 AliModule *detector;
1352 while((detector = (AliModule*)next())) {
fe4da5cc 1353 detector->ResetHits();
1354 }
1355}
1356
1357//_____________________________________________________________________________
1358void AliRun::ResetPoints()
1359{
1360 //
1361 // Reset all Detectors points
1362 //
8494b010 1363 TIter next(fModules);
1364 AliModule *detector;
1365 while((detector = (AliModule*)next())) {
fe4da5cc 1366 detector->ResetPoints();
1367 }
1368}
1369
1370//_____________________________________________________________________________
875c717b 1371void AliRun::RunMC(Int_t nevent, const char *setup)
fe4da5cc 1372{
1373 //
1374 // Main function to be called to process a galice run
1375 // example
1376 // Root > gAlice.Run();
1377 // a positive number of events will cause the finish routine
1378 // to be called
1379 //
1380
fe4da5cc 1381 // check if initialisation has been done
875c717b 1382 if (!fInitDone) InitMC(setup);
fe4da5cc 1383
fe4da5cc 1384 // Create the Root Tree with one branch per detector
80762cb1 1385 MakeTree("KHDER");
fe4da5cc 1386
875c717b 1387 gMC->ProcessRun(nevent);
1388
fe4da5cc 1389 // End of this run, close files
80762cb1 1390 if(nevent>0) FinishRun();
fe4da5cc 1391}
1392
1393//_____________________________________________________________________________
1394void AliRun::RunLego(const char *setup,Int_t ntheta,Float_t themin,
1395 Float_t themax,Int_t nphi,Float_t phimin,Float_t phimax,
1396 Float_t rmin,Float_t rmax,Float_t zmax)
1397{
1398 //
1399 // Generates lego plots of:
1400 // - radiation length map phi vs theta
1401 // - radiation length map phi vs eta
1402 // - interaction length map
1403 // - g/cm2 length map
1404 //
1405 // ntheta bins in theta, eta
1406 // themin minimum angle in theta (degrees)
1407 // themax maximum angle in theta (degrees)
1408 // nphi bins in phi
1409 // phimin minimum angle in phi (degrees)
1410 // phimax maximum angle in phi (degrees)
1411 // rmin minimum radius
1412 // rmax maximum radius
1413 //
1414 //
1415 // The number of events generated = ntheta*nphi
1416 // run input parameters in macro setup (default="Config.C")
1417 //
1418 // Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1419 //Begin_Html
1420 /*
1439f98e 1421 <img src="picts/AliRunLego1.gif">
fe4da5cc 1422 */
1423 //End_Html
1424 //Begin_Html
1425 /*
1439f98e 1426 <img src="picts/AliRunLego2.gif">
fe4da5cc 1427 */
1428 //End_Html
1429 //Begin_Html
1430 /*
1439f98e 1431 <img src="picts/AliRunLego3.gif">
fe4da5cc 1432 */
1433 //End_Html
1434 //
1435
1436 // check if initialisation has been done
875c717b 1437 if (!fInitDone) InitMC(setup);
b13db077 1438
838edcaf 1439 //Save current generator
1440 AliGenerator *gen=Generator();
1441
b13db077 1442 //Create Lego object
1443 fLego = new AliLego("lego",ntheta,themin,themax,nphi,phimin,phimax,rmin,rmax,zmax);
1444
dffd31ef 1445 //Prepare MC for Lego Run
1446 gMC->InitLego();
1447
b13db077 1448 //Run Lego Object
dffd31ef 1449 gMC->ProcessRun(ntheta*nphi+1);
fe4da5cc 1450
1451 // Create only the Root event Tree
80762cb1 1452 MakeTree("E");
fe4da5cc 1453
1454 // End of this run, close files
80762cb1 1455 FinishRun();
838edcaf 1456
1457 // Delete Lego Object
1458 delete fLego; fLego=0;
1459
1460 // Restore current generator
1461 SetGenerator(gen);
fe4da5cc 1462}
1463
1464//_____________________________________________________________________________
1465void AliRun::SetCurrentTrack(Int_t track)
1466{
1467 //
1468 // Set current track number
1469 //
1470 fCurrent = track;
1471}
1472
1473//_____________________________________________________________________________
1578254f 1474void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
fe4da5cc 1475 Float_t *vpos, Float_t *polar, Float_t tof,
d33c0226 1476 const char *mecha, Int_t &ntr, Float_t weight)
fe4da5cc 1477{
1478 //
1479 // Load a track on the stack
1480 //
1481 // done 0 if the track has to be transported
1482 // 1 if not
1483 // parent identifier of the parent track. -1 for a primary
1578254f 1484 // pdg particle code
fe4da5cc 1485 // pmom momentum GeV/c
1486 // vpos position
1487 // polar polarisation
1488 // tof time of flight in seconds
1489 // mecha production mechanism
1490 // ntr on output the number of the track stored
1491 //
1492 TClonesArray &particles = *fParticles;
1578254f 1493 TParticle *particle;
fe4da5cc 1494 Float_t mass;
aee8290b 1495 const Int_t kfirstdaughter=-1;
1496 const Int_t klastdaughter=-1;
1497 const Int_t kS=0;
1578254f 1498 // const Float_t tlife=0;
fe4da5cc 1499
1578254f 1500 //
1501 // Here we get the static mass
1502 // For MC is ok, but a more sophisticated method could be necessary
1503 // if the calculated mass is required
1504 // also, this method is potentially dangerous if the mass
1505 // used in the MC is not the same of the PDG database
1506 //
1507 mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
fe4da5cc 1508 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1509 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1510
aee8290b 1511 //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",
1512 //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS,mecha);
fe4da5cc 1513
aee8290b 1514 particle=new(particles[fNtrack]) TParticle(pdg,kS,parent,-1,kfirstdaughter,
1515 klastdaughter,pmom[0],pmom[1],pmom[2],
1578254f 1516 e,vpos[0],vpos[1],vpos[2],tof);
1517 // polar[0],polar[1],polar[2],tof,
1518 // mecha,weight);
1519 ((TParticle*)particles[fNtrack])->SetPolarisation(TVector3(polar[0],polar[1],polar[2]));
1520 ((TParticle*)particles[fNtrack])->SetWeight(weight);
aee8290b 1521 if(!done) particle->SetBit(kDoneBit);
a01a8b12 1522 //Declare that the daughter information is valid
aee8290b 1523 ((TParticle*)particles[fNtrack])->SetBit(kDaughtersBit);
fe4da5cc 1524
1525 if(parent>=0) {
1578254f 1526 particle=(TParticle*) fParticles->UncheckedAt(parent);
1527 particle->SetLastDaughter(fNtrack);
1528 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
fe4da5cc 1529 } else {
1530 //
1531 // This is a primary track. Set high water mark for this event
1532 fHgwmk=fNtrack;
1533 //
1534 // Set also number if primary tracks
1535 fHeader.SetNprimary(fHgwmk+1);
1536 fHeader.SetNtrack(fHgwmk+1);
1537 }
1538 ntr = fNtrack++;
1539}
1540
1541//_____________________________________________________________________________
1542void AliRun::KeepTrack(const Int_t track)
1543{
1544 //
1545 // flags a track to be kept
1546 //
1547 TClonesArray &particles = *fParticles;
aee8290b 1548 ((TParticle*)particles[track])->SetBit(kKeepBit);
fe4da5cc 1549}
1550
1551//_____________________________________________________________________________
875c717b 1552void AliRun::StepManager(Int_t id)
fe4da5cc 1553{
1554 //
1555 // Called at every step during transport
1556 //
1557
fe4da5cc 1558 //
1559 // --- If lego option, do it and leave
b13db077 1560 if (fLego)
fe4da5cc 1561 fLego->StepManager();
b13db077 1562 else {
1563 Int_t copy;
1564 //Update energy deposition tables
875c717b 1565 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
fe4da5cc 1566
b13db077 1567 //Call the appropriate stepping routine;
1568 AliModule *det = (AliModule*)fModules->At(id);
1569 if(det) det->StepManager();
fe4da5cc 1570 }
fe4da5cc 1571}
1572
aee8290b 1573//_____________________________________________________________________________
1574AliRun& AliRun::operator = (const AliRun &run)
1575{
1576 //
1577 // Assignment operator
1578 //
1579 if(this != &run) run.Copy(*this);
1580 return (*this);
1581}
1582
1583
1584//_____________________________________________________________________________
1585void AliRun::Copy(AliRun &run) const
1586{
1587 //
1588 // Copy *this onto run, not implemented for AliRun
1589 //
1590 Fatal("operator =","Assignment operator not implemented!\n");
1591}
1592
fe4da5cc 1593//_____________________________________________________________________________
1594void AliRun::Streamer(TBuffer &R__b)
1595{
1596 //
1597 // Stream an object of class AliRun.
1598 //
1599 if (R__b.IsReading()) {
1600 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1601 TNamed::Streamer(R__b);
1602 if (!gAlice) gAlice = this;
1603 gROOT->GetListOfBrowsables()->Add(this,"Run");
59fe9bd2 1604 fTreeE = (TTree*)gDirectory->Get("TE");
aee8290b 1605 if (fTreeE) fTreeE->SetBranchAddress("Header", &gAliHeader);
59fe9bd2 1606 else Error("Streamer","cannot find Header Tree\n");
fe4da5cc 1607 R__b >> fNtrack;
1608 R__b >> fHgwmk;
1609 R__b >> fDebug;
1610 fHeader.Streamer(R__b);
8494b010 1611 R__b >> fModules;
fe4da5cc 1612 R__b >> fParticles;
1613 R__b >> fField;
875c717b 1614 // R__b >> fMC;
fe4da5cc 1615 R__b >> fNdets;
1616 R__b >> fTrRmax;
1617 R__b >> fTrZmax;
1618 R__b >> fGenerator;
59fe9bd2 1619 if(R__v>1) {
1620 R__b >> fPDGDB; //Particle factory object!
1621 fTreeE->GetEntry(0);
1622 } else {
1623 fHeader.SetEvent(0);
1624 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
1625 }
45189757 1626 if(R__v>2) {
1627 fConfigFunction.Streamer(R__b);
1628 } else {
1629 fConfigFunction="Config();";
1630 }
fe4da5cc 1631 } else {
1632 R__b.WriteVersion(AliRun::IsA());
1633 TNamed::Streamer(R__b);
1634 R__b << fNtrack;
1635 R__b << fHgwmk;
1636 R__b << fDebug;
1637 fHeader.Streamer(R__b);
8494b010 1638 R__b << fModules;
fe4da5cc 1639 R__b << fParticles;
1640 R__b << fField;
875c717b 1641 // R__b << fMC;
fe4da5cc 1642 R__b << fNdets;
1643 R__b << fTrRmax;
1644 R__b << fTrZmax;
1645 R__b << fGenerator;
1578254f 1646 R__b << fPDGDB; //Particle factory object!
45189757 1647 fConfigFunction.Streamer(R__b);
fe4da5cc 1648 }
1649}