1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.56 2001/02/14 15:45:20 hristov
19 Algorithmic way of getting entry index in fParticleMap. Protection of fParticleFileMap (I.Hrivnacova)
21 Revision 1.55 2001/02/12 15:52:54 buncic
22 Removed OpenBaseFile().
24 Revision 1.54 2001/02/07 10:39:05 hristov
25 Remove default value for argument
27 Revision 1.53 2001/02/06 11:02:26 hristov
28 New SetTrack interface added, added check for unfilled particles in FinishEvent (I.Hrivnacova)
30 Revision 1.52 2001/02/05 16:22:25 buncic
31 Added TreeS to GetEvent().
33 Revision 1.51 2001/02/02 15:16:20 morsch
34 SetHighWaterMark method added to mark last particle in event.
36 Revision 1.50 2001/01/27 10:32:00 hristov
37 Leave the loop when primaries are filled (I.Hrivnacova)
39 Revision 1.49 2001/01/26 19:58:48 hristov
40 Major upgrade of AliRoot code
42 Revision 1.48 2001/01/17 10:50:50 hristov
43 Corrections to destructors
45 Revision 1.47 2000/12/18 10:44:01 morsch
46 Possibility to set field map by passing pointer to objet of type AliMagF via
49 gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
51 Revision 1.46 2000/12/14 19:29:27 fca
52 galice.cuts was not read any more
54 Revision 1.45 2000/11/30 07:12:49 alibrary
55 Introducing new Rndm and QA classes
57 Revision 1.44 2000/10/26 13:58:59 morsch
58 Add possibility to choose the lego generator (of type AliGeneratorLego or derived) when running
59 RunLego(). Default is the base class AliGeneratorLego.
61 Revision 1.43 2000/10/09 09:43:17 fca
62 Special remapping of hits for TPC and TRD. End-of-primary action introduced
64 Revision 1.42 2000/10/02 21:28:14 fca
65 Removal of useless dependecies via forward declarations
67 Revision 1.41 2000/07/13 16:19:09 fca
68 Mainly coding conventions + some small bug fixes
70 Revision 1.40 2000/07/12 08:56:25 fca
71 Coding convention correction and warning removal
73 Revision 1.39 2000/07/11 18:24:59 fca
74 Coding convention corrections + few minor bug fixes
76 Revision 1.38 2000/06/20 13:05:45 fca
77 Writing down the TREE headers before job starts
79 Revision 1.37 2000/06/09 20:05:11 morsch
80 Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat
82 Revision 1.36 2000/06/08 14:03:58 hristov
83 Only one initializer for a default argument
85 Revision 1.35 2000/06/07 10:13:14 hristov
86 Delete only existent objects.
88 Revision 1.34 2000/05/18 10:45:38 fca
89 Delete Particle Factory properly
91 Revision 1.33 2000/05/16 13:10:40 fca
92 New method IsNewTrack and fix for a problem in Father-Daughter relations
94 Revision 1.32 2000/04/27 10:38:21 fca
95 Correct termination of Lego Run and introduce Lego getter in AliRun
97 Revision 1.31 2000/04/26 10:17:32 fca
98 Changes in Lego for G4 compatibility
100 Revision 1.30 2000/04/18 19:11:40 fca
101 Introduce variable Config.C function signature
103 Revision 1.29 2000/04/07 11:12:34 fca
104 G4 compatibility changes
106 Revision 1.28 2000/04/05 06:51:06 fca
107 Workaround for an HP compiler problem
109 Revision 1.27 2000/03/22 18:08:07 fca
110 Rationalisation of the virtual MC interfaces
112 Revision 1.26 2000/03/22 13:42:26 fca
113 SetGenerator does not replace an existing generator, ResetGenerator does
115 Revision 1.25 2000/02/23 16:25:22 fca
116 AliVMC and AliGeant3 classes introduced
117 ReadEuclid moved from AliRun to AliModule
119 Revision 1.24 2000/01/19 17:17:20 fca
120 Introducing a list of lists of hits -- more hits allowed for detector now
122 Revision 1.23 1999/12/03 11:14:31 fca
123 Fixing previous wrong checking
125 Revision 1.21 1999/11/25 10:40:08 fca
126 Fixing daughters information also in primary tracks
128 Revision 1.20 1999/10/04 18:08:49 fca
129 Adding protection against inconsistent Euclid files
131 Revision 1.19 1999/09/29 07:50:40 fca
132 Introduction of the Copyright and cvs Log
136 ///////////////////////////////////////////////////////////////////////////////
138 // Control class for Alice C++ //
139 // Only one single instance of this class exists. //
140 // The object is created in main program aliroot //
141 // and is pointed by the global gAlice. //
143 // -Supports the list of all Alice Detectors (fModules). //
144 // -Supports the list of particles (fParticles). //
145 // -Supports the Trees. //
146 // -Supports the geometry. //
147 // -Supports the event display. //
150 <img src="picts/AliRunClass.gif">
155 <img src="picts/alirun.gif">
159 ///////////////////////////////////////////////////////////////////////////////
164 #include <iostream.h>
172 #include <TObjectTable.h>
174 #include <TGeometry.h>
176 #include "TBrowser.h"
178 #include "TParticle.h"
180 #include "AliDisplay.h"
183 #include "AliMagFC.h"
184 #include "AliMagFCM.h"
185 #include "AliMagFDM.h"
187 #include "TRandom3.h"
189 #include "AliGenerator.h"
190 #include "AliLegoGenerator.h"
192 #include "AliDetector.h"
196 static AliHeader *gAliHeader;
200 //_____________________________________________________________________________
204 // Default constructor for AliRun
229 fPDGDB = 0; //Particle factory object!
231 fConfigFunction = "\0";
234 fTransParName = "\0";
235 fBaseFileName = ".\0";
237 fParticleMap = new TObjArray(10000);
240 //_____________________________________________________________________________
241 AliRun::AliRun(const char *name, const char *title)
245 // Constructor for the main processor.
246 // Creates the geometry
247 // Creates the list of Detectors.
248 // Creates the list of particles.
265 fConfigFunction = "Config();";
267 // Set random number generator
268 gRandom = fRandom = new TRandom3();
270 if (gSystem->Getenv("CONFIG_SEED")) {
271 gRandom->SetSeed((UInt_t)atoi(gSystem->Getenv("CONFIG_SEED")));
274 gROOT->GetListOfBrowsables()->Add(this,name);
276 // create the support list for the various Detectors
277 fModules = new TObjArray(77);
279 // Create the TNode geometry for the event display
281 BuildSimpleGeometry();
291 // Create the particle stack
292 fParticles = new TClonesArray("TParticle",1000);
296 // Create default mag field
301 // Prepare the tracking medium lists
302 fImedia = new TArrayI(1000);
303 for(i=0;i<1000;i++) (*fImedia)[i]=-99;
306 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
308 // Create HitLists list
309 fHitLists = new TList();
312 fBaseFileName = ".\0";
314 fParticleMap = new TObjArray(10000);
318 //_____________________________________________________________________________
322 // Default AliRun destructor
342 fParticles->Delete();
350 //_____________________________________________________________________________
351 void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
354 // Add a hit to detector id
356 TObjArray &dets = *fModules;
357 if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
360 //_____________________________________________________________________________
361 void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
364 // Add digit to detector id
366 TObjArray &dets = *fModules;
367 if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
370 //_____________________________________________________________________________
371 void AliRun::Browse(TBrowser *b)
374 // Called when the item "Run" is clicked on the left pane
375 // of the Root browser.
376 // It displays the Root Trees and all detectors.
378 if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
379 if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
380 if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
381 if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
382 if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
383 if (fTreeS) b->Add(fTreeS,fTreeS->GetName());
385 TIter next(fModules);
387 while((detector = (AliModule*)next())) {
388 b->Add(detector,detector->GetName());
390 b->Add(fMCQA,"AliMCQA");
393 //_____________________________________________________________________________
397 // Initialize Alice geometry
402 //_____________________________________________________________________________
403 void AliRun::BuildSimpleGeometry()
406 // Create a simple TNode geometry used by Root display engine
408 // Initialise geometry
410 fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
411 new TMaterial("void","Vacuum",0,0,0); //Everything is void
412 TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
413 brik->SetVisibility(0);
414 new TNode("alice","alice","S_alice");
417 //_____________________________________________________________________________
418 void AliRun::CleanDetectors()
421 // Clean Detectors at the end of event
423 TIter next(fModules);
425 while((detector = (AliModule*)next())) {
426 detector->FinishEvent();
430 //_____________________________________________________________________________
431 void AliRun::CleanParents()
434 // Clean Particles stack.
435 // Set parent/daughter relations
437 TObjArray &particles = *fParticleMap;
440 for(i=0; i<fHgwmk+1; i++) {
441 part = (TParticle *)particles.At(i);
442 if(part) if(!part->TestBit(kDaughtersBit)) {
443 part->SetFirstDaughter(-1);
444 part->SetLastDaughter(-1);
449 //_____________________________________________________________________________
450 Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
453 // Return the distance from the mouse to the AliRun object
459 //_____________________________________________________________________________
460 void AliRun::DumpPart (Int_t i) const
463 // Dumps particle i in the stack
465 ((TParticle*) (*fParticleMap)[i])->Print();
468 //_____________________________________________________________________________
469 void AliRun::DumpPStack () const
472 // Dumps the particle stack
474 TObjArray &particles = *fParticleMap;
476 "\n\n=======================================================================\n");
477 for (Int_t i=0;i<fNtrack;i++)
479 printf("-> %d ",i); ((TParticle*) particles[i])->Print();
480 printf("--------------------------------------------------------------\n");
483 "\n=======================================================================\n\n");
486 void AliRun::SetField(AliMagF* magField)
488 // Set Magnetic Field Map
493 //_____________________________________________________________________________
494 void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
495 Float_t maxField, char* filename)
498 // Set magnetic field parameters
499 // type Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
500 // version Magnetic field map version (only 1 active now)
501 // scale Scale factor for the magnetic field
502 // maxField Maximum value for the magnetic field
505 // --- Sanity check on mag field flags
506 if(fField) delete fField;
508 fField = new AliMagFC("Map1"," ",type,scale,maxField);
509 } else if(version<=2) {
510 fField = new AliMagFCM("Map2-3",filename,type,scale,maxField);
512 } else if(version==3) {
513 fField = new AliMagFDM("Map4",filename,type,scale,maxField);
516 Warning("SetField","Invalid map %d\n",version);
520 //_____________________________________________________________________________
521 void AliRun::PreTrack()
523 TObjArray &dets = *fModules;
526 for(Int_t i=0; i<=fNdets; i++)
527 if((module = (AliModule*)dets[i]))
533 //_____________________________________________________________________________
534 void AliRun::PostTrack()
536 TObjArray &dets = *fModules;
539 for(Int_t i=0; i<=fNdets; i++)
540 if((module = (AliModule*)dets[i]))
544 //_____________________________________________________________________________
545 void AliRun::FinishPrimary()
548 // Called at the end of each primary track
551 // static Int_t count=0;
552 // const Int_t times=10;
553 // This primary is finished, purify stack
556 TIter next(fModules);
558 while((detector = (AliModule*)next())) {
559 detector->FinishPrimary();
562 // Write out hits if any
563 if (gAlice->TreeH()) {
564 gAlice->TreeH()->Fill();
571 // if(++count%times==1) gObjectTable->Print();
574 //_____________________________________________________________________________
575 void AliRun::FinishEvent()
578 // Called at the end of the event.
582 if(fLego) fLego->FinishEvent();
584 //Update the energy deposit tables
586 for(i=0;i<fEventEnergy.GetSize();i++) {
587 fSummEnergy[i]+=fEventEnergy[i];
588 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
590 fEventEnergy.Reset();
592 // Clean detector information
595 // Write out the kinematics
598 if(fTreeK->GetEntries() ==0) {
599 // set the fParticleFileMap size for the first time
600 fParticleFileMap.Set(fHgwmk+1);
603 Bool_t allFilled = kFALSE;
605 for(i=0; i<fHgwmk+1; ++i) if((part=fParticleMap->At(i))) {
606 fParticleBuffer = (TParticle*) part;
607 fParticleFileMap[i]= (Int_t) fTreeK->GetEntries();
609 (*fParticleMap)[i]=0;
611 // When all primaries were filled no particle!=0
612 // should be left => to be removed later.
613 if (allFilled) printf("Why != 0 part # %d?\n",i);
616 // // printf("Why = 0 part # %d?\n",i); => We know.
618 // we don't break now in order to be sure there is no
619 // particle !=0 left.
620 // To be removed later and replaced with break.
621 if(!allFilled) allFilled = kTRUE;
625 // Set number of tracks to event header
626 fHeader.SetNtrack(fNtrack);
628 // Write out the digits
639 // Write out reconstructed clusters
644 // Write out the event Header information
645 if (fTreeE) fTreeE->Fill();
650 // Write Tree headers
651 if (fTreeK) fTreeK->Write(0,TObject::kOverwrite);
652 if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
653 if (fTreeD) fTreeD->Write(0,TObject::kOverwrite);
654 if (fTreeR) fTreeR->Write(0,TObject::kOverwrite);
655 if (fTreeS) fTreeS->Write(0,TObject::kOverwrite);
660 //_____________________________________________________________________________
661 void AliRun::FinishRun()
664 // Called at the end of the run.
668 if(fLego) fLego->FinishRun();
670 // Clean detector information
671 TIter next(fModules);
673 while((detector = (AliModule*)next())) {
674 detector->FinishRun();
677 //Output energy summary tables
680 TFile *file = fTreeE->GetCurrentFile();
684 fTreeE->Write(0,TObject::kOverwrite);
686 // Write AliRun info and all detectors parameters
687 Write(0,TObject::kOverwrite);
689 // Clean tree information
691 delete fTreeK; fTreeK = 0;
694 delete fTreeH; fTreeH = 0;
697 delete fTreeD; fTreeD = 0;
700 delete fTreeR; fTreeR = 0;
703 delete fTreeE; fTreeE = 0;
710 //_____________________________________________________________________________
711 void AliRun::FlagTrack(Int_t track)
714 // Flags a track and all its family tree to be kept
721 particle=(TParticle*)fParticleMap->At(curr);
723 // If the particle is flagged the three from here upward is saved already
724 if(particle->TestBit(kKeepBit)) return;
726 // Save this particle
727 particle->SetBit(kKeepBit);
729 // Move to father if any
730 if((curr=particle->GetFirstMother())==-1) return;
734 //_____________________________________________________________________________
735 void AliRun::EnergySummary()
738 // Print summary of deposited energy
744 Int_t kn, i, left, j, id;
745 const Float_t kzero=0;
746 Int_t ievent=fHeader.GetEvent()+1;
748 // Energy loss information
750 printf("***************** Energy Loss Information per event (GEV) *****************\n");
751 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
754 fEventEnergy[ndep]=kn;
759 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
762 fSummEnergy[ndep]=ed;
763 fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
768 for(kn=0;kn<(ndep-1)/3+1;kn++) {
770 for(i=0;i<(3<left?3:left);i++) {
772 id=Int_t (fEventEnergy[j]+0.1);
773 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
778 // Relative energy loss in different detectors
779 printf("******************** Relative Energy Loss per event ********************\n");
780 printf("Total energy loss per event %10.3f GeV\n",edtot);
781 for(kn=0;kn<(ndep-1)/5+1;kn++) {
783 for(i=0;i<(5<left?5:left);i++) {
785 id=Int_t (fEventEnergy[j]+0.1);
786 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
790 for(kn=0;kn<75;kn++) printf("*");
794 // Reset the TArray's
795 // fEventEnergy.Set(0);
796 // fSummEnergy.Set(0);
797 // fSum2Energy.Set(0);
800 //_____________________________________________________________________________
801 AliModule *AliRun::GetModule(const char *name) const
804 // Return pointer to detector from name
806 return (AliModule*)fModules->FindObject(name);
809 //_____________________________________________________________________________
810 AliDetector *AliRun::GetDetector(const char *name) const
813 // Return pointer to detector from name
815 return (AliDetector*)fModules->FindObject(name);
818 //_____________________________________________________________________________
819 Int_t AliRun::GetModuleID(const char *name) const
822 // Return galice internal detector identifier from name
825 TObject *mod=fModules->FindObject(name);
826 if(mod) i=fModules->IndexOf(mod);
830 //_____________________________________________________________________________
831 Int_t AliRun::GetEvent(Int_t event)
834 // Connect the Trees Kinematics and Hits for event # event
835 // Set branch addresses
838 // Reset existing structures
844 // Delete Trees already connected
845 if (fTreeK) delete fTreeK;
846 if (fTreeH) delete fTreeH;
847 if (fTreeD) delete fTreeD;
848 if (fTreeR) delete fTreeR;
849 if (fTreeS) delete fTreeS;
851 // Get header from file
852 if(fTreeE) fTreeE->GetEntry(event);
853 else Error("GetEvent","Cannot file Header Tree\n");
854 TFile *file = fTreeE->GetCurrentFile();
858 // Get Kine Tree from file
860 sprintf(treeName,"TreeK%d",event);
861 fTreeK = (TTree*)gDirectory->Get(treeName);
862 if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticleBuffer);
863 else Error("GetEvent","cannot find Kine Tree for event:%d\n",event);
864 // Create the particle stack
865 if(!fParticles) fParticles = new TClonesArray("TParticle",1000);
866 // Build the pointer list
868 fParticleMap->Clear();
869 fParticleMap->Expand(fTreeK->GetEntries());
871 fParticleMap = new TObjArray(fTreeK->GetEntries());
875 // Get Hits Tree header from file
876 sprintf(treeName,"TreeH%d",event);
877 fTreeH = (TTree*)gDirectory->Get(treeName);
879 Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
884 // Get Digits Tree header from file
885 sprintf(treeName,"TreeD%d",event);
886 fTreeD = (TTree*)gDirectory->Get(treeName);
888 // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
893 // Get SDigits Tree header from file
894 sprintf(treeName,"TreeS%d",event);
895 fTreeS = (TTree*)gDirectory->Get(treeName);
897 // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
902 // Get Reconstruct Tree header from file
903 sprintf(treeName,"TreeR%d",event);
904 fTreeR = (TTree*)gDirectory->Get(treeName);
906 // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
911 // Set Trees branch addresses
912 TIter next(fModules);
914 while((detector = (AliModule*)next())) {
915 detector->SetTreeAddress();
918 fNtrack = Int_t (fTreeK->GetEntries());
922 //_____________________________________________________________________________
923 TGeometry *AliRun::GetGeometry()
926 // Import Alice geometry from current file
927 // Return pointer to geometry object
929 if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
931 // Unlink and relink nodes in detectors
932 // This is bad and there must be a better way...
935 TIter next(fModules);
937 while((detector = (AliModule*)next())) {
938 TList *dnodes=detector->Nodes();
941 for ( j=0; j<dnodes->GetSize(); j++) {
942 node = (TNode*) dnodes->At(j);
943 node1 = fGeometry->GetNode(node->GetName());
944 dnodes->Remove(node);
945 dnodes->AddAt(node1,j);
951 //_____________________________________________________________________________
952 void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
953 Float_t &e, Float_t *vpos, Float_t *polar,
957 // Return next track from stack of particles
962 for(Int_t i=fNtrack-1; i>=0; i--) {
963 track=(TParticle*) fParticleMap->At(i);
964 if(track) if(!track->TestBit(kDoneBit)) {
966 // The track exists and has not yet been processed
968 ipart=track->GetPdgCode();
976 track->GetPolarisation(pol);
981 track->SetBit(kDoneBit);
987 // stop and start timer when we start a primary track
988 Int_t nprimaries = fHeader.GetNprimary();
989 if (fCurrent >= nprimaries) return;
990 if (fCurrent < nprimaries-1) {
992 track=(TParticle*) fParticleMap->At(fCurrent+1);
993 // track->SetProcessTime(fTimer.CpuTime());
998 //_____________________________________________________________________________
999 Int_t AliRun::GetPrimary(Int_t track) const
1002 // return number of primary that has generated track
1004 int current, parent;
1010 part = (TParticle *)fParticleMap->At(current);
1011 parent=part->GetFirstMother();
1012 if(parent<0) return current;
1016 //_____________________________________________________________________________
1017 void AliRun::InitMC(const char *setup)
1020 // Initialize the Alice setup
1024 Warning("Init","Cannot initialise AliRun twice!\n");
1028 gROOT->LoadMacro(setup);
1029 gInterpreter->ProcessLine(fConfigFunction.Data());
1032 gMC->DefineParticles(); //Create standard MC particles
1034 TObject *objfirst, *objlast;
1036 fNdets = fModules->GetLast()+1;
1039 //=================Create Materials and geometry
1042 // Added also after in case of interactive initialisation of modules
1043 fNdets = fModules->GetLast()+1;
1045 TIter next(fModules);
1046 AliModule *detector;
1047 while((detector = (AliModule*)next())) {
1048 detector->SetTreeAddress();
1049 objlast = gDirectory->GetList()->Last();
1051 // Add Detector histograms in Detector list of histograms
1052 if (objlast) objfirst = gDirectory->GetList()->After(objlast);
1053 else objfirst = gDirectory->GetList()->First();
1055 detector->Histograms()->Add(objfirst);
1056 objfirst = gDirectory->GetList()->After(objfirst);
1059 ReadTransPar(); //Read the cuts for all materials
1061 MediaTable(); //Build the special IMEDIA table
1063 //Initialise geometry deposition table
1064 fEventEnergy.Set(gMC->NofVolumes()+1);
1065 fSummEnergy.Set(gMC->NofVolumes()+1);
1066 fSum2Energy.Set(gMC->NofVolumes()+1);
1068 //Compute cross-sections
1069 gMC->BuildPhysics();
1071 //Write Geometry object to current file.
1076 fMCQA = new AliMCQA(fNdets);
1079 // Save stuff at the beginning of the file to avoid file corruption
1083 //_____________________________________________________________________________
1084 void AliRun::MediaTable()
1087 // Built media table to get from the media number to
1090 Int_t kz, nz, idt, lz, i, k, ind;
1092 TObjArray &dets = *gAlice->Detectors();
1095 // For all detectors
1096 for (kz=0;kz<fNdets;kz++) {
1097 // If detector is defined
1098 if((det=(AliModule*) dets[kz])) {
1099 TArrayI &idtmed = *(det->GetIdtmed());
1100 for(nz=0;nz<100;nz++) {
1101 // Find max and min material number
1102 if((idt=idtmed[nz])) {
1103 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
1104 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
1107 if(det->LoMedium() > det->HiMedium()) {
1108 det->LoMedium() = 0;
1109 det->HiMedium() = 0;
1111 if(det->HiMedium() > fImedia->GetSize()) {
1112 Error("MediaTable","Increase fImedia from %d to %d",
1113 fImedia->GetSize(),det->HiMedium());
1116 // Tag all materials in rage as belonging to detector kz
1117 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
1124 // Print summary table
1125 printf(" Traking media ranges:\n");
1126 for(i=0;i<(fNdets-1)/6+1;i++) {
1127 for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
1129 det=(AliModule*)dets[ind];
1131 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
1134 printf(" %6s: %3d -> %3d;","NULL",0,0);
1140 //____________________________________________________________________________
1141 void AliRun::SetGenerator(AliGenerator *generator)
1144 // Load the event generator
1146 if(!fGenerator) fGenerator = generator;
1149 //____________________________________________________________________________
1150 void AliRun::ResetGenerator(AliGenerator *generator)
1153 // Load the event generator
1157 Warning("ResetGenerator","Replacing generator %s with %s\n",
1158 fGenerator->GetName(),generator->GetName());
1160 Warning("ResetGenerator","Replacing generator %s with NULL\n",
1161 fGenerator->GetName());
1162 fGenerator = generator;
1165 //____________________________________________________________________________
1166 void AliRun::SetTransPar(char *filename)
1168 fTransParName = filename;
1171 //____________________________________________________________________________
1172 void AliRun::SetBaseFile(char *filename)
1174 fBaseFileName = filename;
1177 //____________________________________________________________________________
1178 void AliRun::ReadTransPar()
1181 // Read filename to set the transport parameters
1185 const Int_t kncuts=10;
1186 const Int_t knflags=11;
1187 const Int_t knpars=kncuts+knflags;
1188 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
1189 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
1190 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
1191 "MULS","PAIR","PHOT","RAYL"};
1195 Float_t cut[kncuts];
1196 Int_t flag[knflags];
1197 Int_t i, itmed, iret, ktmed, kz;
1200 // See whether the file is there
1201 filtmp=gSystem->ExpandPathName(fTransParName.Data());
1202 lun=fopen(filtmp,"r");
1205 Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
1209 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1210 printf(" *%59s\n","*");
1211 printf(" * Please check carefully what you are doing!%10s\n","*");
1212 printf(" *%59s\n","*");
1215 // Initialise cuts and flags
1216 for(i=0;i<kncuts;i++) cut[i]=-99;
1217 for(i=0;i<knflags;i++) flag[i]=-99;
1219 for(i=0;i<256;i++) line[i]='\0';
1220 // Read up to the end of line excluded
1221 iret=fscanf(lun,"%[^\n]",line);
1225 printf(" *%59s\n","*");
1226 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1229 // Read the end of line
1232 if(line[0]=='*') continue;
1234 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",
1235 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1236 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1237 &flag[8],&flag[9],&flag[10]);
1241 Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
1244 // Check that the module exist
1245 AliModule *mod = GetModule(detName);
1247 // Get the array of media numbers
1248 TArrayI &idtmed = *mod->GetIdtmed();
1249 // Check that the tracking medium code is valid
1250 if(0<=itmed && itmed < 100) {
1251 ktmed=idtmed[itmed];
1253 Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
1256 // Set energy thresholds
1257 for(kz=0;kz<kncuts;kz++) {
1259 printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
1260 kpars[kz],cut[kz],itmed,mod->GetName());
1261 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
1264 // Set transport mechanisms
1265 for(kz=0;kz<knflags;kz++) {
1267 printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
1268 kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
1269 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
1273 Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
1277 Warning("ReadTransPar","Module %s not present\n",detName);
1283 //_____________________________________________________________________________
1284 void AliRun::MakeBranchInTree(TTree *tree, const char* name, void* address, Int_t size, char *file)
1287 // Makes branch in given tree and diverts them to a separate file
1290 printf("* MakeBranch * Making Branch %s \n",name);
1292 TBranch *branch = tree->Branch(name,address,size);
1295 char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
1296 sprintf(outFile,"%s/%s",gAlice->GetBaseFile(),file);
1297 TDirectory *cwd = gDirectory;
1298 branch->SetFile(outFile);
1299 TIter next( branch->GetListOfBranches());
1300 while ((branch=(TBranch*)next())) {
1301 branch->SetFile(outFile);
1304 printf("* MakeBranch * Diverting Branch %s to file %s\n",name,file);
1310 //_____________________________________________________________________________
1311 void AliRun::MakeBranchInTree(TTree *tree, const char* name, const char *classname, void* address, Int_t size, Int_t splitlevel, char *file)
1314 // Makes branch in given tree and diverts them to a separate file
1316 TDirectory *cwd = gDirectory;
1317 TBranch *branch = tree->Branch(name,classname,address,size,splitlevel);
1319 printf("* MakeBranch * Making Branch %s \n",name);
1321 char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
1322 sprintf(outFile,"%s/%s",gAlice->GetBaseFile(),file);
1323 branch->SetFile(outFile);
1324 TIter next( branch->GetListOfBranches());
1325 while ((branch=(TBranch*)next())) {
1326 branch->SetFile(outFile);
1329 printf("* MakeBranch * Diverting Branch %s to file %s\n",name,file);
1334 //_____________________________________________________________________________
1335 void AliRun::MakeTree(Option_t *option, char *file)
1338 // Create the ROOT trees
1339 // Loop on all detectors to create the Root branch (if any)
1345 char *oK = strstr(option,"K");
1346 char *oH = strstr(option,"H");
1347 char *oE = strstr(option,"E");
1348 char *oD = strstr(option,"D");
1349 char *oR = strstr(option,"R");
1350 char *oS = strstr(option,"S");
1353 if (oK && !fTreeK) {
1354 sprintf(hname,"TreeK%d",fEvent);
1355 fTreeK = new TTree(hname,"Kinematics");
1356 // Create a branch for particles
1357 MakeBranchInTree(fTreeK,
1358 "Particles", "TParticle", &fParticleBuffer, 4000, 1, file) ;
1361 if (oH && !fTreeH) {
1362 sprintf(hname,"TreeH%d",fEvent);
1363 fTreeH = new TTree(hname,"Hits");
1364 fTreeH->SetAutoSave(1000000000); //no autosave
1367 if (oD && !fTreeD) {
1368 sprintf(hname,"TreeD%d",fEvent);
1369 fTreeD = new TTree(hname,"Digits");
1372 if (oS && !fTreeS) {
1373 sprintf(hname,"TreeS%d",fEvent);
1374 fTreeS = new TTree(hname,"SDigits");
1377 if (oR && !fTreeR) {
1378 sprintf(hname,"TreeR%d",fEvent);
1379 fTreeR = new TTree(hname,"Reconstruction");
1382 if (oE && !fTreeE) {
1383 fTreeE = new TTree("TE","Header");
1384 // Create a branch for Header
1385 MakeBranchInTree(fTreeE,
1386 "Header", "AliHeader", &gAliHeader, 4000, 1, file) ;
1391 // Create a branch for hits/digits for each detector
1392 // Each branch is a TClonesArray. Each data member of the Hits classes
1393 // will be in turn a subbranch of the detector master branch
1394 TIter next(fModules);
1395 AliModule *detector;
1396 while((detector = (AliModule*)next())) {
1397 if (oH) detector->MakeBranch(option,file);
1401 //_____________________________________________________________________________
1402 Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks)
1405 // PurifyKine with external parameters
1407 fHgwmk = lastSavedTrack;
1408 fNtrack = nofTracks;
1413 //_____________________________________________________________________________
1414 TParticle* AliRun::Particle(Int_t i)
1416 if(!(*fParticleMap)[i]) {
1417 Int_t nentries = fParticles->GetEntries();
1419 // algorithmic way of getting entry index
1420 // (primary particles are filled after secondaries)
1422 if (i<fHeader.GetNprimary())
1423 entry = i+fHeader.GetNsecondary();
1425 entry = i-fHeader.GetNprimary();
1427 // only check the algorithmic way and give
1428 // the fatal error if it is wrong
1429 if (entry != fParticleFileMap[i]) {
1431 "!!!! The algorithmic way is WRONG: !!!\n entry: %d map: %d",
1432 entry, fParticleFileMap[i]);
1435 fTreeK->GetEntry(fParticleFileMap[i]);
1436 new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
1437 fParticleMap->AddAt((*fParticles)[nentries],i);
1439 return (TParticle *) (*fParticleMap)[i];
1442 //_____________________________________________________________________________
1443 void AliRun::PurifyKine()
1446 // Compress kinematic tree keeping only flagged particles
1447 // and renaming the particle id's in all the hits
1449 // TClonesArray &particles = *fParticles;
1450 TObjArray &particles = *fParticleMap;
1451 int nkeep=fHgwmk+1, parent, i;
1452 TParticle *part, *father;
1453 TArrayI map(particles.GetLast()+1);
1455 // Save in Header total number of tracks before compression
1456 fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk);
1458 // If no tracks generated return now
1459 if(fHgwmk+1 == fNtrack) return;
1461 Int_t toshrink = fNtrack-fHgwmk-1;
1463 // First pass, invalid Daughter information
1464 for(i=0; i<fNtrack; i++) {
1465 // Preset map, to be removed later
1466 if(i<=fHgwmk) map[i]=i ;
1469 // particles.UncheckedAt(i)->ResetBit(kDaughtersBit);
1470 if((part=(TParticle*) particles.At(i))) part->ResetBit(kDaughtersBit);
1473 // Invalid daughter information for the parent of the first particle
1474 // generated. This may or may not be the current primary according to
1475 // whether decays have been recorded among the primaries
1476 part = (TParticle *)particles.At(fHgwmk+1);
1477 particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
1478 // Second pass, build map between old and new numbering
1479 for(i=fHgwmk+1; i<fNtrack; i++) {
1480 if(particles.At(i)->TestBit(kKeepBit)) {
1482 // This particle has to be kept
1484 // If old and new are different, have to move the pointer
1485 if(i!=nkeep) particles[nkeep]=particles.At(i);
1486 part = (TParticle*) particles.At(nkeep);
1488 // as the parent is always *before*, it must be already
1489 // in place. This is what we are checking anyway!
1490 if((parent=part->GetFirstMother())>fHgwmk)
1491 if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
1492 else part->SetFirstMother(map[parent]);
1498 // Fix daughters information
1499 for (i=fHgwmk+1; i<nkeep; i++) {
1500 part = (TParticle *)particles.At(i);
1501 parent = part->GetFirstMother();
1503 father = (TParticle *)particles.At(parent);
1504 if(father->TestBit(kDaughtersBit)) {
1506 if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
1507 if(i>father->GetLastDaughter()) father->SetLastDaughter(i);
1509 // Initialise daughters info for first pass
1510 father->SetFirstDaughter(i);
1511 father->SetLastDaughter(i);
1512 father->SetBit(kDaughtersBit);
1517 // Now loop on all registered hit lists
1518 TIter next(fHitLists);
1519 TCollection *hitList;
1520 while((hitList = (TCollection*)next())) {
1521 TIter nexthit(hitList);
1523 while((hit = (AliHit*)nexthit())) {
1524 hit->SetTrack(map[hit->GetTrack()]);
1529 // This for detectors which have a special mapping mechanism
1530 // for hits, such as TPC and TRD
1533 TIter nextmod(fModules);
1534 AliModule *detector;
1535 while((detector = (AliModule*)nextmod())) {
1536 detector->RemapTrackHitIDs(map.GetArray());
1539 // Now the output bit, from fHgwmk to nkeep we write everything and we erase
1540 if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
1543 for (i=fHgwmk+1; i<nkeep; ++i) {
1544 fParticleBuffer = (TParticle*) particles.At(i);
1545 fParticleFileMap[i]=(Int_t) fTreeK->GetEntries();
1550 for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
1552 fLoadPoint-=toshrink;
1553 for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
1560 //_____________________________________________________________________________
1561 void AliRun::BeginEvent()
1564 // Reset all Detectors & kinematics & trees
1571 fLego->BeginEvent();
1581 // Initialise event header
1582 fHeader.Reset(fRun,fEvent);
1586 sprintf(hname,"TreeK%d",fEvent);
1587 fTreeK->SetName(hname);
1591 sprintf(hname,"TreeH%d",fEvent);
1592 fTreeH->SetName(hname);
1596 sprintf(hname,"TreeD%d",fEvent);
1597 fTreeD->SetName(hname);
1601 sprintf(hname,"TreeS%d",fEvent);
1602 fTreeS->SetName(hname);
1606 sprintf(hname,"TreeR%d",fEvent);
1607 fTreeR->SetName(hname);
1610 //_____________________________________________________________________________
1611 void AliRun::ResetDigits()
1614 // Reset all Detectors digits
1616 TIter next(fModules);
1617 AliModule *detector;
1618 while((detector = (AliModule*)next())) {
1619 detector->ResetDigits();
1623 //_____________________________________________________________________________
1624 void AliRun::ResetSDigits()
1627 // Reset all Detectors digits
1629 TIter next(fModules);
1630 AliModule *detector;
1631 while((detector = (AliModule*)next())) {
1632 detector->ResetSDigits();
1636 //_____________________________________________________________________________
1637 void AliRun::ResetHits()
1640 // Reset all Detectors hits
1642 TIter next(fModules);
1643 AliModule *detector;
1644 while((detector = (AliModule*)next())) {
1645 detector->ResetHits();
1649 //_____________________________________________________________________________
1650 void AliRun::ResetPoints()
1653 // Reset all Detectors points
1655 TIter next(fModules);
1656 AliModule *detector;
1657 while((detector = (AliModule*)next())) {
1658 detector->ResetPoints();
1662 //_____________________________________________________________________________
1663 void AliRun::RunMC(Int_t nevent, const char *setup)
1666 // Main function to be called to process a galice run
1668 // Root > gAlice.Run();
1669 // a positive number of events will cause the finish routine
1673 // check if initialisation has been done
1674 if (!fInitDone) InitMC(setup);
1676 // Create the Root Tree with one branch per detector
1680 if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
1681 MakeTree("K","Kine.root");
1682 MakeTree("H","Hits.root");
1687 gMC->ProcessRun(nevent);
1689 // End of this run, close files
1690 if(nevent>0) FinishRun();
1693 //_____________________________________________________________________________
1694 void AliRun::RunReco(const char *detector)
1697 // Main function to be called to reconstruct Alice event
1701 Digits2Reco(detector);
1704 //_____________________________________________________________________________
1706 void AliRun::Hits2Digits(const char *selected)
1708 // Convert Hits to sumable digits
1710 Hits2SDigits(selected);
1711 SDigits2Digits(selected);
1715 //_____________________________________________________________________________
1717 void AliRun::Tree2Tree(Option_t *option, const char *selected)
1720 // Function to transform the content of
1722 // - TreeH to TreeS (option "S")
1723 // - TreeS to TreeD (option "D")
1724 // - TreeD to TreeR (option "R")
1726 // If multiple options are specified ("SDR"), transformation will be done in sequence for
1727 // selected detector and for all detectors if none is selected (detector string
1728 // can contain blank separated list of detector names).
1731 char *oS = strstr(option,"S");
1732 char *oD = strstr(option,"D");
1733 char *oR = strstr(option,"R");
1735 gAlice->GetEvent(0);
1737 TObjArray *detectors = gAlice->Detectors();
1739 TIter next(detectors);
1741 AliDetector *detector = 0;
1743 TDirectory *cwd = gDirectory;
1747 while((detector = (AliDetector*)next())) {
1749 if (strcmp(detector->GetName(),selected)) continue;
1750 if (detector->IsActive()){
1752 cout << "Processing " << detector->GetName() << "..." << endl;
1753 if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
1755 sprintf(outFile,"SDigits.%s.root",detector->GetName());
1756 detector->MakeBranch("S",outFile);
1759 sprintf(outFile,"Digits.%s.root",detector->GetName());
1760 detector->MakeBranch("D",outFile);
1763 sprintf(outFile,"Reco.%s.root",detector->GetName());
1764 detector->MakeBranch("R",outFile);
1767 detector->MakeBranch(option);
1773 detector->Hits2SDigits();
1775 detector->SDigits2Digits();
1777 detector->Digits2Reco();
1786 //_____________________________________________________________________________
1787 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1788 Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1789 Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
1792 // Generates lego plots of:
1793 // - radiation length map phi vs theta
1794 // - radiation length map phi vs eta
1795 // - interaction length map
1796 // - g/cm2 length map
1798 // ntheta bins in theta, eta
1799 // themin minimum angle in theta (degrees)
1800 // themax maximum angle in theta (degrees)
1802 // phimin minimum angle in phi (degrees)
1803 // phimax maximum angle in phi (degrees)
1804 // rmin minimum radius
1805 // rmax maximum radius
1808 // The number of events generated = ntheta*nphi
1809 // run input parameters in macro setup (default="Config.C")
1811 // Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1814 <img src="picts/AliRunLego1.gif">
1819 <img src="picts/AliRunLego2.gif">
1824 <img src="picts/AliRunLego3.gif">
1829 // check if initialisation has been done
1830 if (!fInitDone) InitMC(setup);
1831 //Save current generator
1832 AliGenerator *gen=Generator();
1834 // Set new generator
1835 if (!gener) gener = new AliLegoGenerator();
1836 ResetGenerator(gener);
1838 // Configure Generator
1839 gener->SetRadiusRange(rmin, rmax);
1840 gener->SetZMax(zmax);
1841 gener->SetCoor1Range(nc1, c1min, c1max);
1842 gener->SetCoor2Range(nc2, c2min, c2max);
1845 //Create Lego object
1846 fLego = new AliLego("lego",gener);
1848 //Prepare MC for Lego Run
1853 gMC->ProcessRun(nc1*nc2+1);
1855 // Create only the Root event Tree
1858 // End of this run, close files
1860 // Restore current generator
1861 ResetGenerator(gen);
1862 // Delete Lego Object
1863 delete fLego; fLego=0;
1866 //_____________________________________________________________________________
1867 void AliRun::SetConfigFunction(const char * config)
1870 // Set the signature of the function contained in Config.C to configure
1873 fConfigFunction=config;
1876 //_____________________________________________________________________________
1877 void AliRun::SetCurrentTrack(Int_t track)
1880 // Set current track number
1885 //_____________________________________________________________________________
1886 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1887 Float_t *vpos, Float_t *polar, Float_t tof,
1888 AliMCProcess mech, Int_t &ntr, Float_t weight)
1891 // Load a track on the stack
1893 // done 0 if the track has to be transported
1895 // parent identifier of the parent track. -1 for a primary
1896 // pdg particle code
1897 // pmom momentum GeV/c
1899 // polar polarisation
1900 // tof time of flight in seconds
1901 // mecha production mechanism
1902 // ntr on output the number of the track stored
1904 TClonesArray &particles = *fParticles;
1905 TParticle *particle;
1907 const Int_t kfirstdaughter=-1;
1908 const Int_t klastdaughter=-1;
1910 // const Float_t tlife=0;
1913 // Here we get the static mass
1914 // For MC is ok, but a more sophisticated method could be necessary
1915 // if the calculated mass is required
1916 // also, this method is potentially dangerous if the mass
1917 // used in the MC is not the same of the PDG database
1919 mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1920 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1921 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1923 //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",
1924 //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS,mecha);
1926 particle=new(particles[fLoadPoint++]) TParticle(pdg,kS,parent,-1,kfirstdaughter,
1927 klastdaughter,pmom[0],pmom[1],pmom[2],
1928 e,vpos[0],vpos[1],vpos[2],tof);
1929 particle->SetPolarisation(TVector3(polar[0],polar[1],polar[2]));
1930 particle->SetWeight(weight);
1931 particle->SetUniqueID(mech);
1932 if(!done) particle->SetBit(kDoneBit);
1933 // Declare that the daughter information is valid
1934 particle->SetBit(kDaughtersBit);
1935 // Add the particle to the stack
1936 fParticleMap->AddAtAndExpand(particle,fNtrack);
1939 particle=(TParticle*) fParticleMap->At(parent);
1940 particle->SetLastDaughter(fNtrack);
1941 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
1944 // This is a primary track. Set high water mark for this event
1947 // Set also number if primary tracks
1948 fHeader.SetNprimary(fHgwmk+1);
1949 fHeader.SetNtrack(fHgwmk+1);
1955 // Here we get the static mass
1956 // For MC is ok, but a more sophisticated method could be necessary
1957 // if the calculated mass is required
1958 // also, this method is potentially dangerous if the mass
1959 // used in the MC is not the same of the PDG database
1961 Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1962 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1963 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1965 SetTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
1966 vpos[0], vpos[1], vpos[2], tof, polar[0],polar[1],polar[2],
1971 //_____________________________________________________________________________
1972 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg,
1973 Double_t px, Double_t py, Double_t pz, Double_t e,
1974 Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1975 Double_t polx, Double_t poly, Double_t polz,
1976 AliMCProcess mech, Int_t &ntr, Float_t weight)
1979 // Load a track on the stack
1981 // done 0 if the track has to be transported
1983 // parent identifier of the parent track. -1 for a primary
1984 // pdg particle code
1985 // kS generation status code
1986 // px, py, pz momentum GeV/c
1987 // vx, vy, vz position
1988 // polar polarisation
1989 // tof time of flight in seconds
1990 // mech production mechanism
1991 // ntr on output the number of the track stored
1993 // New method interface:
1994 // arguments were changed to be in correspondence with TParticle
1996 // Note: the energy is not calculated from the static mass but
1997 // it is passed by argument e.
1999 TClonesArray &particles = *fParticles;
2002 const Int_t kFirstDaughter=-1;
2003 const Int_t kLastDaughter=-1;
2006 = new(particles[fLoadPoint++]) TParticle(pdg, kS, parent, -1,
2007 kFirstDaughter, kLastDaughter,
2008 px, py, pz, e, vx, vy, vz, tof);
2010 particle->SetPolarisation(polx, poly, polz);
2011 particle->SetWeight(weight);
2012 particle->SetUniqueID(mech);
2014 if(!done) particle->SetBit(kDoneBit);
2016 // Declare that the daughter information is valid
2017 particle->SetBit(kDaughtersBit);
2018 // Add the particle to the stack
2019 fParticleMap->AddAtAndExpand(particle,fNtrack);
2022 particle=(TParticle*) fParticleMap->At(parent);
2023 particle->SetLastDaughter(fNtrack);
2024 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
2027 // This is a primary track. Set high water mark for this event
2030 // Set also number if primary tracks
2031 fHeader.SetNprimary(fHgwmk+1);
2032 fHeader.SetNtrack(fHgwmk+1);
2037 //_____________________________________________________________________________
2038 void AliRun::SetHighWaterMark(const Int_t nt)
2041 // Set high water mark for last track in event
2044 // Set also number if primary tracks
2045 fHeader.SetNprimary(fHgwmk+1);
2046 fHeader.SetNtrack(fHgwmk+1);
2049 //_____________________________________________________________________________
2050 void AliRun::KeepTrack(const Int_t track)
2053 // flags a track to be kept
2055 fParticleMap->At(track)->SetBit(kKeepBit);
2058 //_____________________________________________________________________________
2059 void AliRun::StepManager(Int_t id)
2062 // Called at every step during transport
2066 // --- If lego option, do it and leave
2068 fLego->StepManager();
2071 //Update energy deposition tables
2072 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
2074 //Call the appropriate stepping routine;
2075 AliModule *det = (AliModule*)fModules->At(id);
2077 fMCQA->StepManager(id);
2083 //_____________________________________________________________________________
2084 void AliRun::Streamer(TBuffer &R__b)
2086 // Stream an object of class AliRun.
2088 if (R__b.IsReading()) {
2089 if (!gAlice) gAlice = this;
2091 AliRun::Class()->ReadBuffer(R__b, this);
2093 gROOT->GetListOfBrowsables()->Add(this,"Run");
2095 fTreeE = (TTree*)gDirectory->Get("TE");
2096 if (fTreeE) fTreeE->SetBranchAddress("Header", &gAliHeader);
2097 else Error("Streamer","cannot find Header Tree\n");
2098 fTreeE->GetEntry(0);
2102 AliRun::Class()->WriteBuffer(R__b, this);