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