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