STEER/AliRun.cxx

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