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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Control class for Alice C++                                              //
  21 //  Only one single instance of this class exists.                           //
  22 //  The object is created in main program aliroot                            //
  23 //  and is pointed by the global gAlice.                                     //
  24 //                                                                           //
  25 //   -Supports the list of all Alice Detectors (fModules).                 //
  26 //   -Supports the list of particles (fParticles).                           //
  27 //   -Supports the Trees.                                                    //
  28 //   -Supports the geometry.                                                 //
  29 //   -Supports the event display.                                            //
  30 //Begin_Html
  31 /*
  32 <img src="picts/AliRunClass.gif">
  33 */
  34 //End_Html
  35 //Begin_Html
  36 /*
  37 <img src="picts/alirun.gif">
  38 */
  39 //End_Html
  40 //                                                                           //
  41 ///////////////////////////////////////////////////////////////////////////////
  42
  43 #include <Riostream.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 #include <string.h>
  47
  48 #include <TBRIK.h>
  49 #include <TBrowser.h>
  50 #include <TCint.h>
  51 #include <TFile.h>
  52 #include <TFolder.h>
  53 #include <TGeometry.h>
  54 #include <TKey.h>
  55 #include <TNode.h>
  56 #include <TObjectTable.h>
  57 #include <TParticle.h>
  58 #include <TROOT.h>
  59 #include <TRandom.h>
  60 #include <TRandom3.h>
  61 #include <TSystem.h>
  62 #include <TTree.h>
  63 #include <TVirtualMC.h>
  64
  65 #include "AliConfig.h"
  66 #include "AliDetector.h"
  67 #include "AliDisplay.h"
  68 #include "AliGenEventHeader.h"
  69 #include "AliGenerator.h"
  70 #include "AliHeader.h"
  71 #include "AliHit.h"
  72 #include "AliLego.h"
  73 #include "AliLegoGenerator.h"
  74 #include "AliLoader.h"
  75 #include "AliMCQA.h"
  76 #include "AliMagFC.h"
  77 #include "AliMagFCM.h"
  78 #include "AliMagFDM.h"
  79 #include "AliPDG.h"
  80 #include "AliRun.h"
  81 #include "AliRunLoader.h"
  82 #include "AliStack.h"
  83 #include "AliTrackReference.h"
  84
  85 AliRun *gAlice;
  86
  87 ClassImp(AliRun)
  88
  89 //_______________________________________________________________________
  90 AliRun::AliRun():
  91   fRun(0),
  92   fEvent(0),
  93   fEventNrInRun(0),
  94   fEventsPerRun(0),
  95   fDebug(0),
  96   fModules(0),
  97   fGeometry(0),
  98   fDisplay(0),
  99   fTimer(),
 100   fField(0),
 101   fMC(0),
 102   fImedia(0),
 103   fNdets(0),
 104   fTrRmax(1.e10),
 105   fTrZmax(1.e10),
 106   fGenerator(0),
 107   fInitDone(kFALSE),
 108   fLego(0),
 109   fPDGDB(0),  //Particle factory object
 110   fHitLists(0),
 111   fEventEnergy(0),
 112   fSummEnergy(0),
 113   fSum2Energy(0),
 114   fConfigFunction("\0"),
 115   fRandom(0),
 116   fMCQA(0),
 117   fTransParName("\0"),
 118   fRunLoader(0x0),
 119   fTrackReferences(0)
 120 {
 121   //
 122   // Default constructor for AliRun
 123   //
 124   AliConfig::Instance();//skowron 29 Feb 2002
 125                         //ensures that the folder structure is build
 126 }
 127
 128 //_______________________________________________________________________
 129 AliRun::AliRun(const AliRun& arun):
 130   TVirtualMCApplication(arun),
 131   fRun(0),
 132   fEvent(0),
 133   fEventNrInRun(0),
 134   fEventsPerRun(0),
 135   fDebug(0),
 136   fModules(0),
 137   fGeometry(0),
 138   fDisplay(0),
 139   fTimer(),
 140   fField(0),
 141   fMC(0),
 142   fImedia(0),
 143   fNdets(0),
 144   fTrRmax(1.e10),
 145   fTrZmax(1.e10),
 146   fGenerator(0),
 147   fInitDone(kFALSE),
 148   fLego(0),
 149   fPDGDB(0),  //Particle factory object
 150   fHitLists(0),
 151   fEventEnergy(0),
 152   fSummEnergy(0),
 153   fSum2Energy(0),
 154   fConfigFunction("\0"),
 155   fRandom(0),
 156   fMCQA(0),
 157   fTransParName("\0"),
 158   fTrackReferences(new TClonesArray("AliTrackReference", 100)),
 159   fRunLoader(0x0)
 160 {
 161   //
 162   // Copy constructor for AliRun
 163   //
 164   arun.Copy(*this);
 165 }
 166
 167 //_____________________________________________________________________________
 168 AliRun::AliRun(const char *name, const char *title):
 169   TVirtualMCApplication(name,title),
 170   fRun(0),
 171   fEvent(0),
 172   fEventNrInRun(0),
 173   fEventsPerRun(0),
 174   fDebug(0),
 175   fModules(new TObjArray(77)), // Support list for the Detectors
 176   fGeometry(0),
 177   fDisplay(0),
 178   fTimer(),
 179   fField(0),
 180   fMC(gMC),
 181   fImedia(new TArrayI(1000)),
 182   fNdets(0),
 183   fTrRmax(1.e10),
 184   fTrZmax(1.e10),
 185   fGenerator(0),
 186   fInitDone(kFALSE),
 187   fLego(0),
 188   fPDGDB(TDatabasePDG::Instance()),        //Particle factory object!
 189   fHitLists(new TList()),                  // Create HitLists list
 190   fEventEnergy(0),
 191   fSummEnergy(0),
 192   fSum2Energy(0),
 193   fConfigFunction("Config();"),
 194   fRandom(new TRandom3()),
 195   fMCQA(0),
 196   fTransParName("\0"),
 197   fRunLoader(0x0)
 198 {
 199   //
 200   //  Constructor for the main processor.
 201   //  Creates the geometry
 202   //  Creates the list of Detectors.
 203   //  Creates the list of particles.
 204   //
 205
 206   gAlice     = this;
 207
 208   // Set random number generator
 209   gRandom = fRandom;
 210
 211   if (gSystem->Getenv("CONFIG_SEED")) {
 212      gRandom->SetSeed(static_cast<UInt_t>(atoi(gSystem->Getenv("CONFIG_SEED"))));
 213   }
 214
 215   // Add to list of browsable
 216   gROOT->GetListOfBrowsables()->Add(this,name);
 217   // Create the TNode geometry for the event display
 218   BuildSimpleGeometry();
 219
 220   // Create default mag field
 221   SetField();
 222
 223   // Prepare the tracking medium lists
 224   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
 225
 226   // Add particle list to configuration
 227   AliConfig::Instance()->Add(fPDGDB);
 228
 229   // Set transport parameters
 230   SetTransPar();
 231 }
 232
 233
 234 //_______________________________________________________________________
 235 AliRun::~AliRun()
 236 {
 237   //
 238   // Default AliRun destructor
 239   //
 240   gROOT->GetListOfBrowsables()->Remove(this);
 241
 242   if (fRunLoader)
 243    {
 244     TFolder* evfold = fRunLoader->GetEventFolder();
 245     TFolder* modfold = dynamic_cast<TFolder*>(evfold->FindObjectAny(AliConfig::GetModulesFolderName()));
 246     TIter next(fModules);
 247     AliModule *mod;
 248     while((mod = (AliModule*)next()))
 249      {
 250        modfold->Remove(mod);
 251      }
 252    }
 253
 254   delete fImedia;
 255   delete fField;
 256   // delete fMC;
 257   delete gMC; gMC=0;
 258   delete fGeometry;
 259   delete fDisplay;
 260   delete fGenerator;
 261   delete fLego;
 262   if (fModules) {
 263     fModules->Delete();
 264     delete fModules;
 265   }
 266
 267   delete fHitLists;
 268   delete fPDGDB;
 269   delete fMCQA;
 270   // Delete track references
 271   if (fTrackReferences) {
 272     fTrackReferences->Delete();
 273     delete fTrackReferences;
 274     fTrackReferences     = 0;
 275   }
 276
 277 }
 278
 279 //_______________________________________________________________________
 280 void AliRun::Copy(AliRun &) const
 281 {
 282   Fatal("Copy","Not implemented!\n");
 283 }
 284
 285 //_______________________________________________________________________
 286 void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 287 {
 288   //
 289   //  Add a hit to detector id
 290   //
 291   TObjArray &dets = *fModules;
 292   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 293 }
 294
 295 //_______________________________________________________________________
 296 void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 297 {
 298   //
 299   // Add digit to detector id
 300   //
 301   TObjArray &dets = *fModules;
 302   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 303 }
 304
 305 //_______________________________________________________________________
 306 void AliRun::Browse(TBrowser *b)
 307 {
 308   //
 309   // Called when the item "Run" is clicked on the left pane
 310   // of the Root browser.
 311   // It displays the Root Trees and all detectors.
 312   //
 313   //detectors are in folders anyway
 314   b->Add(fMCQA,"AliMCQA");
 315 }
 316
 317 //_______________________________________________________________________
 318 void AliRun::Build()
 319 {
 320   //
 321   // Initialize Alice geometry
 322   // Dummy routine
 323   //
 324 }
 325
 326 //_______________________________________________________________________
 327 void AliRun::BuildSimpleGeometry()
 328 {
 329   //
 330   // Create a simple TNode geometry used by Root display engine
 331   //
 332   // Initialise geometry
 333   //
 334   fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
 335   new TMaterial("void","Vacuum",0,0,0);  //Everything is void
 336   TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
 337   brik->SetVisibility(0);
 338   new TNode("alice","alice","S_alice");
 339 }
 340
 341 //_______________________________________________________________________
 342 void AliRun::CleanDetectors()
 343 {
 344   //
 345   // Clean Detectors at the end of event
 346   //
 347    fRunLoader->CleanDetectors();
 348 }
 349
 350 //_______________________________________________________________________
 351 Int_t AliRun::DistancetoPrimitive(Int_t, Int_t) const
 352 {
 353   //
 354   // Return the distance from the mouse to the AliRun object
 355   // Dummy routine
 356   //
 357   return 9999;
 358 }
 359
 360 //_______________________________________________________________________
 361 void AliRun::DumpPart (Int_t i) const
 362 {
 363   //
 364   // Dumps particle i in the stack
 365   //
 366    if (fRunLoader->Stack())
 367     fRunLoader->Stack()->DumpPart(i);
 368 }
 369
 370 //_______________________________________________________________________
 371 void AliRun::DumpPStack () const
 372 {
 373   //
 374   // Dumps the particle stack
 375   //
 376    if (fRunLoader->Stack())
 377     fRunLoader->Stack()->DumpPStack();
 378 }
 379
 380 //_______________________________________________________________________
 381 void  AliRun::SetField(AliMagF* magField)
 382 {
 383     // Set Magnetic Field Map
 384     fField = magField;
 385     fField->ReadField();
 386 }
 387
 388 //_______________________________________________________________________
 389 void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
 390                       Float_t maxField, char* filename)
 391 {
 392   //
 393   //  Set magnetic field parameters
 394   //  type      Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
 395   //  version   Magnetic field map version (only 1 active now)
 396   //  scale     Scale factor for the magnetic field
 397   //  maxField  Maximum value for the magnetic field
 398
 399   //
 400   // --- Sanity check on mag field flags
 401   if(fField) delete fField;
 402   if(version==1) {
 403     fField = new AliMagFC("Map1"," ",type,scale,maxField);
 404   } else if(version<=2) {
 405     fField = new AliMagFCM("Map2-3",filename,type,scale,maxField);
 406     fField->ReadField();
 407   } else if(version==3) {
 408     fField = new AliMagFDM("Map4",filename,type,scale,maxField);
 409     fField->ReadField();
 410   } else {
 411     Warning("SetField","Invalid map %d\n",version);
 412   }
 413 }
 414
 415 //_____________________________________________________________________________
 416
 417 void AliRun::InitLoaders()
 418 {
 419   //creates list of getters
 420   if (GetDebug()) Info("InitLoaders","");
 421   TIter next(fModules);
 422   AliModule *mod;
 423   while((mod = (AliModule*)next()))
 424    {
 425      AliDetector *det = dynamic_cast<AliDetector*>(mod);
 426      if (det)
 427       {
 428         if (GetDebug()) Info("InitLoaders"," Adding %s ",det->GetName());
 429         fRunLoader->AddLoader(det);
 430       }
 431    }
 432   if (GetDebug()) Info("InitLoaders","Done");
 433 }
 434 //_____________________________________________________________________________
 435
 436 void AliRun::FinishRun()
 437 {
 438   //
 439   // Called at the end of the run.
 440   //
 441
 442   if(fLego)
 443    {
 444     if (GetDebug()) Info("FinishRun"," Finish Lego");
 445     fRunLoader->CdGAFile();
 446     fLego->FinishRun();
 447    }
 448
 449   // Clean detector information
 450   TIter next(fModules);
 451   AliModule *detector;
 452   while((detector = dynamic_cast<AliModule*>(next()))) {
 453     if (GetDebug()) Info("FinishRun"," %s->FinishRun()",detector->GetName());
 454     detector->FinishRun();
 455   }
 456
 457   //Output energy summary tables
 458   if (GetDebug()) Info("FinishRun"," EnergySummary()");
 459   EnergySummary();
 460
 461   if (GetDebug()) Info("FinishRun"," fRunLoader->WriteHeader(OVERWRITE)");
 462   fRunLoader->WriteHeader("OVERWRITE");
 463
 464   // Write AliRun info and all detectors parameters
 465   fRunLoader->CdGAFile();
 466   Write(0,TObject::kOverwrite);//write AliRun
 467   fRunLoader->Write(0,TObject::kOverwrite);//write RunLoader itself
 468
 469   // Clean tree information
 470   if (GetDebug()) Info("FinishRun"," fRunLoader->Stack()->FinishRun()");
 471   fRunLoader->Stack()->FinishRun();
 472
 473   // Clean detector information
 474   if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()");
 475   fGenerator->FinishRun();
 476   fRunLoader->Synchronize();
 477 }
 478
 479 //_______________________________________________________________________
 480 void AliRun::FlagTrack(Int_t track)
 481 {
 482   // Delegate to stack
 483   //
 484     fRunLoader->Stack()->FlagTrack(track);
 485 }
 486
 487 //_______________________________________________________________________
 488 void AliRun::EnergySummary()
 489 {
 490   //
 491   // Print summary of deposited energy
 492   //
 493
 494   Int_t ndep=0;
 495   Float_t edtot=0;
 496   Float_t ed, ed2;
 497   Int_t kn, i, left, j, id;
 498   const Float_t kzero=0;
 499   Int_t ievent=fRunLoader->GetHeader()->GetEvent()+1;
 500   //
 501   // Energy loss information
 502   if(ievent) {
 503     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 504     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 505       ed=fSummEnergy[kn];
 506       if(ed>0) {
 507         fEventEnergy[ndep]=kn;
 508         if(ievent>1) {
 509           ed=ed/ievent;
 510           ed2=fSum2Energy[kn];
 511           ed2=ed2/ievent;
 512           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
 513         } else
 514           ed2=99;
 515         fSummEnergy[ndep]=ed;
 516         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
 517         edtot+=ed;
 518         ndep++;
 519       }
 520     }
 521     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 522       left=ndep-kn*3;
 523       for(i=0;i<(3<left?3:left);i++) {
 524         j=kn*3+i;
 525         id=Int_t (fEventEnergy[j]+0.1);
 526         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 527       }
 528       printf("\n");
 529     }
 530     //
 531     // Relative energy loss in different detectors
 532     printf("******************** Relative Energy Loss per event ********************\n");
 533     printf("Total energy loss per event %10.3f GeV\n",edtot);
 534     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 535       left=ndep-kn*5;
 536       for(i=0;i<(5<left?5:left);i++) {
 537         j=kn*5+i;
 538         id=Int_t (fEventEnergy[j]+0.1);
 539         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 540       }
 541       printf("\n");
 542     }
 543     for(kn=0;kn<75;kn++) printf("*");
 544     printf("\n");
 545   }
 546   //
 547   // Reset the TArray's
 548   //  fEventEnergy.Set(0);
 549   //  fSummEnergy.Set(0);
 550   //  fSum2Energy.Set(0);
 551 }
 552
 553 //_______________________________________________________________________
 554 void AliRun::Announce() const
 555 {
 556   //
 557   // Announce the current version of AliRoot
 558   //
 559   printf("%70s",
 560          "****************************************************************\n");
 561   printf("%6s","*");printf("%64s","*\n");
 562
 563   printf("%6s","*");
 564   printf("    You are running AliRoot version NewIO\n");
 565
 566   printf("%6s","*");
 567   printf("    The cvs tag for the current program is $Name$\n");
 568
 569   printf("%6s","*");printf("%64s","*\n");
 570   printf("%70s",
 571          "****************************************************************\n");
 572 }
 573
 574 //_______________________________________________________________________
 575 AliModule *AliRun::GetModule(const char *name) const
 576 {
 577   //
 578   // Return pointer to detector from name
 579   //
 580   return dynamic_cast<AliModule*>(fModules->FindObject(name));
 581 }
 582
 583 //_______________________________________________________________________
 584 AliDetector *AliRun::GetDetector(const char *name) const
 585 {
 586   //
 587   // Return pointer to detector from name
 588   //
 589   return dynamic_cast<AliDetector*>(fModules->FindObject(name));
 590 }
 591
 592 //_______________________________________________________________________
 593 Int_t AliRun::GetModuleID(const char *name) const
 594 {
 595   //
 596   // Return galice internal detector identifier from name
 597   //
 598   Int_t i=-1;
 599   TObject *mod=fModules->FindObject(name);
 600   if(mod) i=fModules->IndexOf(mod);
 601   return i;
 602 }
 603
 604 //_______________________________________________________________________
 605 Int_t AliRun::GetEvent(Int_t event)
 606 {
 607 //
 608 // Reloads data containers in folders # event
 609 // Set branch addresses
 610 //
 611   if (fRunLoader == 0x0)
 612    {
 613      Error("GetEvent","RunLoader is not set. Can not load data.");
 614      return -1;
 615    }
 616 /*****************************************/
 617 /****   P R E    R E L O A D I N G    ****/
 618 /*****************************************/
 619 // Reset existing structures
 620   ResetHits();
 621   ResetTrackReferences();
 622   ResetDigits();
 623   ResetSDigits();
 624
 625 /*****************************************/
 626 /****       R  E  L  O  A  D          ****/
 627 /*****************************************/
 628
 629   fRunLoader->GetEvent(event);
 630
 631 /*****************************************/
 632 /****  P O S T    R E L O A D I N G   ****/
 633 /*****************************************/
 634
 635   // Set Trees branch addresses
 636   TIter next(fModules);
 637   AliModule *detector;
 638   while((detector = dynamic_cast<AliModule*>(next())))
 639    {
 640      detector->SetTreeAddress();
 641    }
 642
 643   return fRunLoader->GetHeader()->GetNtrack();
 644 }
 645
 646 //_______________________________________________________________________
 647 TGeometry *AliRun::GetGeometry()
 648 {
 649   //
 650   // Import Alice geometry from current file
 651   // Return pointer to geometry object
 652   //
 653   if (!fGeometry) fGeometry = dynamic_cast<TGeometry*>(gDirectory->Get("AliceGeom"));
 654   //
 655   // Unlink and relink nodes in detectors
 656   // This is bad and there must be a better way...
 657   //
 658
 659   TIter next(fModules);
 660   AliModule *detector;
 661   while((detector = dynamic_cast<AliModule*>(next()))) {
 662     TList *dnodes=detector->Nodes();
 663     Int_t j;
 664     TNode *node, *node1;
 665     for ( j=0; j<dnodes->GetSize(); j++) {
 666       node = dynamic_cast<TNode*>(dnodes->At(j));
 667       node1 = fGeometry->GetNode(node->GetName());
 668       dnodes->Remove(node);
 669       dnodes->AddAt(node1,j);
 670     }
 671   }
 672   return fGeometry;
 673 }
 674
 675 //_______________________________________________________________________
 676 Int_t AliRun::GetPrimary(Int_t track) const
 677 {
 678   //
 679   // return number of primary that has generated track
 680   //
 681     return fRunLoader->Stack()->GetPrimary(track);
 682 }
 683
 684 //_______________________________________________________________________
 685 void AliRun::MediaTable()
 686 {
 687   //
 688   // Built media table to get from the media number to
 689   // the detector id
 690   //
 691
 692   Int_t kz, nz, idt, lz, i, k, ind;
 693   //  Int_t ibeg;
 694   TObjArray &dets = *gAlice->Detectors();
 695   AliModule *det;
 696   //
 697   // For all detectors
 698   for (kz=0;kz<fNdets;kz++) {
 699     // If detector is defined
 700     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
 701         TArrayI &idtmed = *(det->GetIdtmed());
 702         for(nz=0;nz<100;nz++) {
 703         // Find max and min material number
 704         if((idt=idtmed[nz])) {
 705           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 706           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 707         }
 708       }
 709       if(det->LoMedium() > det->HiMedium()) {
 710         det->LoMedium() = 0;
 711         det->HiMedium() = 0;
 712       } else {
 713         if(det->HiMedium() > fImedia->GetSize()) {
 714           Error("MediaTable","Increase fImedia from %d to %d",
 715                 fImedia->GetSize(),det->HiMedium());
 716           return;
 717         }
 718         // Tag all materials in rage as belonging to detector kz
 719         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 720           (*fImedia)[lz]=kz;
 721         }
 722       }
 723     }
 724   }
 725   //
 726   // Print summary table
 727   printf(" Traking media ranges:\n");
 728   for(i=0;i<(fNdets-1)/6+1;i++) {
 729     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
 730       ind=i*6+k;
 731       det=dynamic_cast<AliModule*>(dets[ind]);
 732       if(det)
 733         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 734                det->HiMedium());
 735       else
 736         printf(" %6s: %3d -> %3d;","NULL",0,0);
 737     }
 738     printf("\n");
 739   }
 740 }
 741
 742 //_______________________________________________________________________
 743 void AliRun::SetGenerator(AliGenerator *generator)
 744 {
 745   //
 746   // Load the event generator
 747   //
 748   if(!fGenerator) fGenerator = generator;
 749 }
 750
 751 //_______________________________________________________________________
 752 void AliRun::ResetGenerator(AliGenerator *generator)
 753 {
 754   //
 755   // Load the event generator
 756   //
 757   if(fGenerator)
 758     if(generator)
 759       Warning("ResetGenerator","Replacing generator %s with %s\n",
 760               fGenerator->GetName(),generator->GetName());
 761     else
 762       Warning("ResetGenerator","Replacing generator %s with NULL\n",
 763               fGenerator->GetName());
 764   fGenerator = generator;
 765 }
 766
 767 //_______________________________________________________________________
 768 void AliRun::SetTransPar(const char *filename)
 769 {
 770   //
 771   // Sets the file name for transport parameters
 772   //
 773   fTransParName = filename;
 774 }
 775
 776 //_______________________________________________________________________
 777 void AliRun::SetBaseFile(const char *filename)
 778 {
 779   fBaseFileName = filename;
 780 }
 781
 782 //_______________________________________________________________________
 783 void AliRun::ReadTransPar()
 784 {
 785   //
 786   // Read filename to set the transport parameters
 787   //
 788
 789
 790   const Int_t kncuts=10;
 791   const Int_t knflags=11;
 792   const Int_t knpars=kncuts+knflags;
 793   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 794                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 795                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 796                                "MULS","PAIR","PHOT","RAYL"};
 797   char line[256];
 798   char detName[7];
 799   char* filtmp;
 800   Float_t cut[kncuts];
 801   Int_t flag[knflags];
 802   Int_t i, itmed, iret, ktmed, kz;
 803   FILE *lun;
 804   //
 805   // See whether the file is there
 806   filtmp=gSystem->ExpandPathName(fTransParName.Data());
 807   lun=fopen(filtmp,"r");
 808   delete [] filtmp;
 809   if(!lun) {
 810     Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
 811     return;
 812   }
 813   //
 814   if(fDebug) {
 815     printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 816     printf(" *%59s\n","*");
 817     printf(" *       Please check carefully what you are doing!%10s\n","*");
 818     printf(" *%59s\n","*");
 819   }
 820   //
 821   while(1) {
 822     // Initialise cuts and flags
 823     for(i=0;i<kncuts;i++) cut[i]=-99;
 824     for(i=0;i<knflags;i++) flag[i]=-99;
 825     itmed=0;
 826     for(i=0;i<256;i++) line[i]='\0';
 827     // Read up to the end of line excluded
 828     iret=fscanf(lun,"%[^\n]",line);
 829     if(iret<0) {
 830       //End of file
 831       fclose(lun);
 832       if(fDebug){
 833         printf(" *%59s\n","*");
 834         printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 835       }
 836       return;
 837     }
 838     // Read the end of line
 839     fscanf(lun,"%*c");
 840     if(!iret) continue;
 841     if(line[0]=='*') continue;
 842     // Read the numbers
 843     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",
 844                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
 845                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
 846                 &flag[8],&flag[9],&flag[10]);
 847     if(!iret) continue;
 848     if(iret<0) {
 849       //reading error
 850       Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
 851       continue;
 852     }
 853     // Check that the module exist
 854     AliModule *mod = GetModule(detName);
 855     if(mod) {
 856       // Get the array of media numbers
 857       TArrayI &idtmed = *mod->GetIdtmed();
 858       // Check that the tracking medium code is valid
 859       if(0<=itmed && itmed < 100) {
 860         ktmed=idtmed[itmed];
 861         if(!ktmed) {
 862           Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
 863           continue;
 864         }
 865         // Set energy thresholds
 866         for(kz=0;kz<kncuts;kz++) {
 867           if(cut[kz]>=0) {
 868             if(fDebug) printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
 869                    kpars[kz],cut[kz],itmed,mod->GetName());
 870             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
 871           }
 872         }
 873         // Set transport mechanisms
 874         for(kz=0;kz<knflags;kz++) {
 875           if(flag[kz]>=0) {
 876             if(fDebug) printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
 877                    kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
 878             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
 879           }
 880         }
 881       } else {
 882         Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
 883         continue;
 884       }
 885     } else {
 886       if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
 887       continue;
 888     }
 889   }
 890 }
 891 //_____________________________________________________________________________
 892
 893 void AliRun::BeginEvent()
 894 {
 895   //
 896   // Clean-up previous event
 897   // Energy scores
 898   if (GetDebug())
 899    {
 900      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 901      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 902      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 903      Info("BeginEvent","          BEGINNING EVENT               ");
 904      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 905      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 906      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 907    }
 908
 909   /*******************************/
 910   /*   Clean after eventual      */
 911   /*   previous event            */
 912   /*******************************/
 913
 914
 915   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
 916   fRunLoader->SetEventNumber(++fEventNrInRun);// sets new files, cleans the previous event stuff, if necessary, etc.,
 917   if (GetDebug()) Info("BeginEvent","EventNr is %d",fEventNrInRun);
 918
 919   fEventEnergy.Reset();
 920     // Clean detector information
 921
 922   if (fRunLoader->Stack())
 923     fRunLoader->Stack()->Reset();//clean stack -> tree is unloaded
 924   else
 925     fRunLoader->MakeStack();//or make a new one
 926
 927   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(K)");
 928   fRunLoader->MakeTree("K");
 929   if (GetDebug()) Info("BeginEvent","  gMC->SetStack(fRunLoader->Stack())");
 930   gMC->SetStack(fRunLoader->Stack());//Was in InitMC - but was moved here
 931                                      //because we don't have guarantee that
 932                                      //stack pointer is not going to change from event to event
 933                          //since it bellobgs to header and is obtained via RunLoader
 934   //
 935   //  Reset all Detectors & kinematics & make/reset trees
 936   //
 937
 938   fRunLoader->GetHeader()->Reset(fRun,fEvent,fEventNrInRun);
 939 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 940
 941   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTrackRefsContainer()");
 942   fRunLoader->MakeTrackRefsContainer();//for insurance
 943
 944   if (GetDebug()) Info("BeginEvent","  ResetHits()");
 945   ResetHits();
 946   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(H)");
 947   fRunLoader->MakeTree("H");
 948
 949   //
 950   if(fLego)
 951    {
 952     fLego->BeginEvent();
 953     return;
 954    }
 955
 956   //create new branches and SetAdresses
 957   TIter next(fModules);
 958   AliModule *detector;
 959   while((detector = (AliModule*)next()))
 960    {
 961     if (GetDebug()) Info("BeginEvent","  %s->MakeBranch(H)",detector->GetName());
 962     detector->MakeBranch("H");
 963     if (GetDebug()) Info("BeginEvent","  %s->MakeBranchTR()",detector->GetName());
 964     detector->MakeBranchTR();
 965     if (GetDebug()) Info("BeginEvent","  %s->SetTreeAddress()",detector->GetName());
 966     detector->SetTreeAddress();
 967    }
 968   // make branch for AliRun track References
 969   TTree * treeTR = fRunLoader->TreeTR();
 970   if (treeTR){
 971     // make branch for central track references
 972     if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
 973     TBranch *branch;
 974     branch = treeTR->Branch("AliRun",&fTrackReferences);
 975     branch->SetAddress(&fTrackReferences);
 976   }
 977   //
 978 }
 979
 980 //_______________________________________________________________________
 981 TParticle* AliRun::Particle(Int_t i) const
 982 {
 983   if (fRunLoader)
 984    if (fRunLoader->Stack())
 985     return fRunLoader->Stack()->Particle(i);
 986   return 0x0;
 987 }
 988
 989 //_______________________________________________________________________
 990 void AliRun::ResetDigits()
 991 {
 992   //
 993   //  Reset all Detectors digits
 994   //
 995   TIter next(fModules);
 996   AliModule *detector;
 997   while((detector = dynamic_cast<AliModule*>(next()))) {
 998      detector->ResetDigits();
 999   }
1000 }
1001
1002 //_______________________________________________________________________
1003 void AliRun::ResetSDigits()
1004 {
1005   //
1006   //  Reset all Detectors digits
1007   //
1008   TIter next(fModules);
1009   AliModule *detector;
1010   while((detector = dynamic_cast<AliModule*>(next()))) {
1011      detector->ResetSDigits();
1012   }
1013 }
1014
1015 //_______________________________________________________________________
1016 void AliRun::ResetHits()
1017 {
1018   //
1019   //  Reset all Detectors hits
1020   //
1021   TIter next(fModules);
1022   AliModule *detector;
1023   while((detector = dynamic_cast<AliModule*>(next()))) {
1024      detector->ResetHits();
1025   }
1026 }
1027 //_______________________________________________________________________
1028
1029 void  AliRun::AddTrackReference(Int_t label){
1030   //
1031   // add a trackrefernce to the list
1032   if (!fTrackReferences) {
1033     cerr<<"Container trackrefernce not active\n";
1034     return;
1035   }
1036   Int_t nref = fTrackReferences->GetEntriesFast();
1037   TClonesArray &lref = *fTrackReferences;
1038   new(lref[nref]) AliTrackReference(label);
1039 }
1040
1041
1042
1043 void AliRun::ResetTrackReferences()
1044 {
1045   //
1046   //  Reset all  references
1047   //
1048   if (fTrackReferences)   fTrackReferences->Clear();
1049
1050   TIter next(fModules);
1051   AliModule *detector;
1052   while((detector = dynamic_cast<AliModule*>(next()))) {
1053      detector->ResetTrackReferences();
1054   }
1055 }
1056
1057 void AliRun::RemapTrackReferencesIDs(Int_t *map)
1058 {
1059   //
1060   // Remapping track reference
1061   // Called at finish primary
1062   //
1063   if (!fTrackReferences) return;
1064   for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
1065     AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
1066     if (ref) {
1067       Int_t newID = map[ref->GetTrack()];
1068       if (newID>=0) ref->SetTrack(newID);
1069       else {
1070         //ref->SetTrack(-1);
1071         ref->SetBit(kNotDeleted,kFALSE);
1072         fTrackReferences->RemoveAt(i);
1073       }
1074     }
1075   }
1076   fTrackReferences->Compress();
1077 }
1078
1079
1080
1081 //_______________________________________________________________________
1082
1083 void AliRun::ResetPoints()
1084 {
1085   //
1086   // Reset all Detectors points
1087   //
1088   TIter next(fModules);
1089   AliModule *detector;
1090   while((detector = dynamic_cast<AliModule*>(next()))) {
1091      detector->ResetPoints();
1092   }
1093 }
1094 //_______________________________________________________________________
1095
1096 void AliRun::InitMC(const char *setup)
1097 {
1098   //
1099   // Initialize the Alice setup
1100   //
1101   Announce();
1102
1103   if(fInitDone) {
1104     Warning("Init","Cannot initialise AliRun twice!\n");
1105     return;
1106   }
1107
1108   gROOT->LoadMacro(setup);
1109   gInterpreter->ProcessLine(fConfigFunction.Data());
1110
1111   // Register MC in configuration
1112   AliConfig::Instance()->Add(gMC);
1113
1114   InitLoaders();
1115
1116   fRunLoader->MakeTree("E");
1117   fRunLoader->LoadKinematics("RECREATE");
1118   fRunLoader->LoadTrackRefs("RECREATE");
1119   fRunLoader->LoadHits("all","RECREATE");
1120
1121
1122   fRunLoader->CdGAFile();
1123
1124   gMC->DefineParticles();  //Create standard MC particles
1125   AliPDG::AddParticlesToPdgDataBase();
1126
1127   TObject *objfirst, *objlast;
1128
1129   fNdets = fModules->GetLast()+1;
1130
1131   //
1132   //=================Create Materials and geometry
1133   gMC->Init();
1134
1135   // Added also after in case of interactive initialisation of modules
1136   fNdets = fModules->GetLast()+1;
1137
1138   TIter next(fModules);
1139   AliModule *detector;
1140   while((detector = dynamic_cast<AliModule*>(next())))
1141    {
1142      objlast = gDirectory->GetList()->Last();
1143
1144      // Add Detector histograms in Detector list of histograms
1145      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
1146      else         objfirst = gDirectory->GetList()->First();
1147      while (objfirst)
1148       {
1149         detector->Histograms()->Add(objfirst);
1150         objfirst = gDirectory->GetList()->After(objfirst);
1151       }
1152    }
1153    ReadTransPar(); //Read the cuts for all materials
1154
1155    MediaTable(); //Build the special IMEDIA table
1156
1157    //Initialise geometry deposition table
1158    fEventEnergy.Set(gMC->NofVolumes()+1);
1159    fSummEnergy.Set(gMC->NofVolumes()+1);
1160    fSum2Energy.Set(gMC->NofVolumes()+1);
1161
1162    //Compute cross-sections
1163    gMC->BuildPhysics();
1164
1165    //Write Geometry object to current file.
1166    fRunLoader->WriteGeometry();
1167
1168    fInitDone = kTRUE;
1169
1170    fMCQA = new AliMCQA(fNdets);
1171
1172    //
1173    // Save stuff at the beginning of the file to avoid file corruption
1174    Write();
1175    fEventNrInRun = -1; //important - we start Begin event from increasing current number in run
1176 }
1177
1178 //_______________________________________________________________________
1179
1180 void AliRun::RunMC(Int_t nevent, const char *setup)
1181 {
1182   //
1183   // Main function to be called to process a galice run
1184   // example
1185   //    Root > gAlice.Run();
1186   // a positive number of events will cause the finish routine
1187   // to be called
1188   //
1189   fEventsPerRun = nevent;
1190   // check if initialisation has been done
1191   if (!fInitDone) InitMC(setup);
1192
1193   // Create the Root Tree with one branch per detector
1194   //Hits moved to begin event -> now we are crating separate tree for each event
1195
1196   gMC->ProcessRun(nevent);
1197
1198   // End of this run, close files
1199   if(nevent>0) FinishRun();
1200 }
1201
1202 //_______________________________________________________________________
1203 void AliRun::RunReco(const char *selected, Int_t first, Int_t last)
1204 {
1205   //
1206   // Main function to be called to reconstruct Alice event
1207   //
1208    Int_t nev = fRunLoader->GetNumberOfEvents();
1209    if (GetDebug()) Info("RunReco","Found %d events",nev);
1210    Int_t nFirst = first;
1211    Int_t nLast  = (last < 0)? nev : last;
1212
1213    for (Int_t nevent = nFirst; nevent <= nLast; nevent++) {
1214      if (GetDebug()) Info("RunReco","Processing event %d",nevent);
1215      GetEvent(nevent);
1216      Digits2Reco(selected);
1217    }
1218 }
1219
1220 //_______________________________________________________________________
1221
1222 void AliRun::Hits2Digits(const char *selected)
1223 {
1224
1225    // Convert Hits to sumable digits
1226    //
1227    for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
1228      GetEvent(nevent);
1229      Hits2SDigits(selected);
1230      SDigits2Digits(selected);
1231    }
1232 }
1233
1234
1235 //_______________________________________________________________________
1236
1237 void AliRun::Tree2Tree(Option_t *option, const char *selected)
1238 {
1239   //
1240   // Function to transform the content of
1241   //
1242   // - TreeH to TreeS (option "S")
1243   // - TreeS to TreeD (option "D")
1244   // - TreeD to TreeR (option "R")
1245   //
1246   // If multiple options are specified ("SDR"), transformation will be done in sequence for
1247   // selected detector and for all detectors if none is selected (detector string
1248   // can contain blank separated list of detector names).
1249
1250
1251    const char *oS = strstr(option,"S");
1252    const char *oD = strstr(option,"D");
1253    const char *oR = strstr(option,"R");
1254
1255    TObjArray *detectors = Detectors();
1256
1257    TIter next(detectors);
1258
1259    AliDetector *detector = 0;
1260
1261    while((detector = dynamic_cast<AliDetector*>(next()))) {
1262      if (selected)
1263        if (strcmp(detector->GetName(),selected)) continue;
1264      if (detector->IsActive())
1265       {
1266
1267        AliLoader* loader = detector->GetLoader();
1268        if (loader == 0x0) continue;
1269
1270        if (oS)
1271         {
1272           if (GetDebug()) Info("Tree2Tree","Processing Hits2SDigits for %s ...",detector->GetName());
1273           loader->LoadHits("read");
1274           if (loader->TreeS() == 0x0) loader->MakeTree("S");
1275           detector->MakeBranch(option);
1276           detector->SetTreeAddress();
1277           detector->Hits2SDigits();
1278           loader->UnloadHits();
1279           loader->UnloadSDigits();
1280         }
1281        if (oD)
1282         {
1283           if (GetDebug()) Info("Tree2Tree","Processing SDigits2Digits for %s ...",detector->GetName());
1284           loader->LoadSDigits("read");
1285           if (loader->TreeD() == 0x0) loader->MakeTree("D");
1286           detector->MakeBranch(option);
1287           detector->SetTreeAddress();
1288           detector->SDigits2Digits();
1289           loader->UnloadSDigits();
1290           loader->UnloadDigits();
1291         }
1292        if (oR)
1293         {
1294           if (GetDebug()) Info("Tree2Tree","Processing Digits2Reco for %s ...",detector->GetName());
1295           loader->LoadDigits("read");
1296           if (loader->TreeR() == 0x0) loader->MakeTree("R");
1297           detector->MakeBranch(option);
1298           detector->SetTreeAddress();
1299           detector->Digits2Reco();
1300           loader->UnloadDigits();
1301           loader->UnloadRecPoints();
1302
1303         }
1304      }
1305    }
1306 }
1307
1308 //_______________________________________________________________________
1309 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1310                      Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1311                      Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
1312 {
1313   //
1314   // Generates lego plots of:
1315   //    - radiation length map phi vs theta
1316   //    - radiation length map phi vs eta
1317   //    - interaction length map
1318   //    - g/cm2 length map
1319   //
1320   //  ntheta    bins in theta, eta
1321   //  themin    minimum angle in theta (degrees)
1322   //  themax    maximum angle in theta (degrees)
1323   //  nphi      bins in phi
1324   //  phimin    minimum angle in phi (degrees)
1325   //  phimax    maximum angle in phi (degrees)
1326   //  rmin      minimum radius
1327   //  rmax      maximum radius
1328   //
1329   //
1330   //  The number of events generated = ntheta*nphi
1331   //  run input parameters in macro setup (default="Config.C")
1332   //
1333   //  Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1334   //Begin_Html
1335   /*
1336     <img src="picts/AliRunLego1.gif">
1337   */
1338   //End_Html
1339   //Begin_Html
1340   /*
1341     <img src="picts/AliRunLego2.gif">
1342   */
1343   //End_Html
1344   //Begin_Html
1345   /*
1346     <img src="picts/AliRunLego3.gif">
1347   */
1348   //End_Html
1349   //
1350
1351   // check if initialisation has been done
1352   if (!fInitDone) InitMC(setup);
1353   // Save current generator
1354   AliGenerator *gen=Generator();
1355   // If runloader has been initialized, set the number of events per file to nc1 * nc2
1356   if (fRunLoader) fRunLoader->SetNumberOfEventsPerFile(nc1 * nc2);
1357   // Set new generator
1358   if (!gener) gener  = new AliLegoGenerator();
1359   ResetGenerator(gener);
1360   //
1361   // Configure Generator
1362   gener->SetRadiusRange(rmin, rmax);
1363   gener->SetZMax(zmax);
1364   gener->SetCoor1Range(nc1, c1min, c1max);
1365   gener->SetCoor2Range(nc2, c2min, c2max);
1366
1367
1368   //Create Lego object
1369   fLego = new AliLego("lego",gener);
1370
1371   //Prepare MC for Lego Run
1372   gMC->InitLego();
1373
1374   //Run Lego Object
1375
1376   //gMC->ProcessRun(nc1*nc2+1);
1377   gMC->ProcessRun(nc1*nc2);
1378
1379   // Create only the Root event Tree
1380   fRunLoader->MakeTree("E");
1381
1382   // End of this run, close files
1383   FinishRun();
1384   // Restore current generator
1385   ResetGenerator(gen);
1386   // Delete Lego Object
1387   delete fLego; fLego=0;
1388 }
1389
1390 //_______________________________________________________________________
1391 void AliRun::SetConfigFunction(const char * config)
1392 {
1393   //
1394   // Set the signature of the function contained in Config.C to configure
1395   // the run
1396   //
1397   fConfigFunction=config;
1398 }
1399
1400 //_______________________________________________________________________
1401 void AliRun::SetCurrentTrack(Int_t track)
1402 {
1403   //
1404   // Set current track number
1405   //
1406     fRunLoader->Stack()->SetCurrentTrack(track);
1407 }
1408
1409 //_______________________________________________________________________
1410 void AliRun::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1411                       Float_t *vpos, Float_t *polar, Float_t tof,
1412                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
1413 {
1414 // Delegate to stack
1415 //
1416     fRunLoader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1417                      mech, ntr, weight, is);
1418 }
1419
1420 //_______________________________________________________________________
1421 void AliRun::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1422                       Double_t px, Double_t py, Double_t pz, Double_t e,
1423                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1424                       Double_t polx, Double_t poly, Double_t polz,
1425                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
1426 {
1427   // Delegate to stack
1428   //
1429   fRunLoader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1430                                 polx, poly, polz, mech, ntr, weight, is);
1431 }
1432
1433 //_______________________________________________________________________
1434 void AliRun::SetHighWaterMark(const Int_t nt)
1435 {
1436     //
1437     // Set high water mark for last track in event
1438     fRunLoader->Stack()->SetHighWaterMark(nt);
1439 }
1440
1441 //_______________________________________________________________________
1442 void AliRun::KeepTrack(const Int_t track)
1443 {
1444   //
1445   // Delegate to stack
1446   //
1447     fRunLoader->Stack()->KeepTrack(track);
1448 }
1449
1450 //
1451 // MC Application
1452 //
1453
1454 //_______________________________________________________________________
1455 void  AliRun::ConstructGeometry()
1456 {
1457   //
1458   // Create modules, materials, geometry
1459   //
1460
1461     TStopwatch stw;
1462     TIter next(fModules);
1463     AliModule *detector;
1464     if (GetDebug()) Info("ConstructGeometry","Geometry creation:");
1465     while((detector = dynamic_cast<AliModule*>(next()))) {
1466       stw.Start();
1467       // Initialise detector materials and geometry
1468       detector->CreateMaterials();
1469       detector->CreateGeometry();
1470       printf("%10s R:%.2fs C:%.2fs\n",
1471              detector->GetName(),stw.RealTime(),stw.CpuTime());
1472     }
1473 }
1474
1475 //_______________________________________________________________________
1476 void  AliRun::InitGeometry()
1477 {
1478   //
1479   // Initialize detectors and display geometry
1480   //
1481
1482    printf("Initialisation:\n");
1483     TStopwatch stw;
1484     TIter next(fModules);
1485     AliModule *detector;
1486     while((detector = dynamic_cast<AliModule*>(next()))) {
1487       stw.Start();
1488       // Initialise detector and display geometry
1489       detector->Init();
1490       detector->BuildGeometry();
1491       printf("%10s R:%.2fs C:%.2fs\n",
1492              detector->GetName(),stw.RealTime(),stw.CpuTime());
1493     }
1494
1495 }
1496 //_______________________________________________________________________
1497
1498 void  AliRun::GeneratePrimaries()
1499 {
1500   //
1501   // Generate primary particles and fill them in the stack.
1502   //
1503
1504   Generator()->Generate();
1505 }
1506 //_______________________________________________________________________
1507
1508 void AliRun::BeginPrimary()
1509 {
1510   //
1511   // Called  at the beginning of each primary track
1512   //
1513
1514   // Reset Hits info
1515   gAlice->ResetHits();
1516   gAlice->ResetTrackReferences();
1517
1518 }
1519
1520 //_______________________________________________________________________
1521 void AliRun::PreTrack()
1522 {
1523      TObjArray &dets = *fModules;
1524      AliModule *module;
1525
1526      for(Int_t i=0; i<=fNdets; i++)
1527        if((module = dynamic_cast<AliModule*>(dets[i])))
1528          module->PreTrack();
1529
1530      fMCQA->PreTrack();
1531 }
1532
1533 //_______________________________________________________________________
1534 void AliRun::Stepping()
1535 {
1536   //
1537   // Called at every step during transport
1538   //
1539
1540   Int_t id = DetFromMate(gMC->GetMedium());
1541
1542   if (id < 0) return;
1543
1544   //
1545   // --- If lego option, do it and leave
1546   if (fLego)
1547     fLego->StepManager();
1548   else {
1549     Int_t copy;
1550     //Update energy deposition tables
1551     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
1552     //
1553     // write tracke reference for track which is dissapearing - MI
1554     if (gMC->IsTrackDisappeared()) {
1555       if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber());
1556     }
1557
1558     //Call the appropriate stepping routine;
1559     AliModule *det = dynamic_cast<AliModule*>(fModules->At(id));
1560     if(det && det->StepManagerIsEnabled()) {
1561       fMCQA->StepManager(id);
1562       det->StepManager();
1563     }
1564   }
1565 }
1566
1567 //_______________________________________________________________________
1568 void AliRun::PostTrack()
1569 {
1570      TObjArray &dets = *fModules;
1571      AliModule *module;
1572
1573      for(Int_t i=0; i<=fNdets; i++)
1574        if((module = dynamic_cast<AliModule*>(dets[i])))
1575          module->PostTrack();
1576 }
1577
1578 //_______________________________________________________________________
1579 void AliRun::FinishPrimary()
1580 {
1581   //
1582   // Called  at the end of each primary track
1583   //
1584
1585   //  static Int_t count=0;
1586   //  const Int_t times=10;
1587   // This primary is finished, purify stack
1588   fRunLoader->Stack()->PurifyKine();
1589
1590   TIter next(fModules);
1591   AliModule *detector;
1592   while((detector = dynamic_cast<AliModule*>(next()))) {
1593     detector->FinishPrimary();
1594     if(detector->GetLoader())
1595      {
1596        detector->GetLoader()->TreeH()->Fill();
1597      }
1598   }
1599
1600   // Write out track references if any
1601   if (fRunLoader->TreeTR())
1602    {
1603     fRunLoader->TreeTR()->Fill();
1604    }
1605 }
1606
1607 //_______________________________________________________________________
1608 void AliRun::FinishEvent()
1609 {
1610   //
1611   // Called at the end of the event.
1612   //
1613
1614   //
1615   if(fLego) fLego->FinishEvent();
1616
1617   TIter next(fModules);
1618   AliModule *detector;
1619   while((detector = dynamic_cast<AliModule*>(next()))) {
1620     detector->FinishEvent();
1621   }
1622
1623   //Update the energy deposit tables
1624   Int_t i;
1625   for(i=0;i<fEventEnergy.GetSize();i++)
1626    {
1627     fSummEnergy[i]+=fEventEnergy[i];
1628     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
1629    }
1630
1631   AliHeader* header = fRunLoader->GetHeader();
1632   AliStack* stack = fRunLoader->Stack();
1633   if ( (header == 0x0) || (stack == 0x0) )
1634    {//check if we got header and stack. If not cry and exit aliroot
1635     Fatal("AliRun","Can not get the stack or header from LOADER");
1636     return;//never reached
1637    }
1638   // Update Header information
1639   header->SetNprimary(stack->GetNprimary());
1640   header->SetNtrack(stack->GetNtrack());
1641
1642
1643   // Write out the kinematics
1644   stack->FinishEvent();
1645
1646   // Write out the event Header information
1647   TTree* treeE = fRunLoader->TreeE();
1648   if (treeE)
1649    {
1650       header->SetStack(stack);
1651       treeE->Fill();
1652    }
1653   else
1654    {
1655     Error("FinishEvent","Can not get TreeE from RL");
1656    }
1657
1658   fRunLoader->WriteKinematics("OVERWRITE");
1659   fRunLoader->WriteTrackRefs("OVERWRITE");
1660   fRunLoader->WriteHits("OVERWRITE");
1661
1662   if (GetDebug())
1663    {
1664      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1665      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1666      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1667      Info("FinishEvent","          FINISHING EVENT               ");
1668      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1669      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1670      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1671    }
1672 }
1673
1674 //_______________________________________________________________________
1675 void AliRun::Field(const Double_t* x, Double_t *b) const
1676 {
1677    Float_t xfloat[3];
1678    for (Int_t i=0; i<3; i++) xfloat[i] = x[i];
1679
1680    if (Field()) {
1681          Float_t bfloat[3];
1682          Field()->Field(xfloat,bfloat);
1683          for (Int_t j=0; j<3; j++) b[j] = bfloat[j];
1684    }
1685    else {
1686          printf("No mag field defined!\n");
1687          b[0]=b[1]=b[2]=0.;
1688    }
1689
1690 }
1691
1692 //
1693 // End of MC Application
1694 //
1695
1696 //_______________________________________________________________________
1697 void AliRun::Streamer(TBuffer &R__b)
1698 {
1699   // Stream an object of class AliRun.
1700
1701   if (R__b.IsReading()) {
1702     if (!gAlice) gAlice = this;
1703     AliRun::Class()->ReadBuffer(R__b, this);
1704     gROOT->GetListOfBrowsables()->Add(this,"Run");
1705
1706     gRandom = fRandom;
1707   } else {
1708     AliRun::Class()->WriteBuffer(R__b, this);
1709   }
1710 }
1711
1712
1713 //_______________________________________________________________________
1714 Int_t AliRun::GetCurrentTrackNumber() const {
1715   //
1716   // Returns current track
1717   //
1718   return fRunLoader->Stack()->GetCurrentTrackNumber();
1719 }
1720
1721 //_______________________________________________________________________
1722 Int_t AliRun::GetNtrack() const {
1723   //
1724   // Returns number of tracks in stack
1725   //
1726   return fRunLoader->Stack()->GetNtrack();
1727 }
1728 //_______________________________________________________________________
1729
1730 //_______________________________________________________________________
1731 TObjArray* AliRun::Particles() const {
1732   //
1733   // Returns pointer to Particles array
1734   //
1735   if (fRunLoader)
1736    if (fRunLoader->Stack())
1737     return fRunLoader->Stack()->Particles();
1738   return 0x0;
1739 }
1740
1741 //___________________________________________________________________________
1742
1743 //_______________________________________________________________________
1744 void AliRun::SetGenEventHeader(AliGenEventHeader* header)
1745 {
1746   fRunLoader->GetHeader()->SetGenEventHeader(header);
1747 }
1748
1749 //___________________________________________________________________________
1750
1751 Int_t AliRun::GetEvNumber() const
1752 {
1753 //Returns number of current event
1754   if (fRunLoader == 0x0)
1755    {
1756      Error("GetEvent","RunLoader is not set. Can not load data.");
1757      return -1;
1758    }
1759
1760   return fRunLoader->GetEventNumber();
1761 }
1762
1763 void AliRun::SetRunLoader(AliRunLoader* rloader)
1764 {
1765   fRunLoader = rloader;
1766   if (fRunLoader == 0x0) return;
1767
1768   TString evfoldname;
1769   TFolder* evfold = fRunLoader->GetEventFolder();
1770   if (evfold) evfoldname = evfold->GetName();
1771   else Warning("SetRunLoader","Did not get Event Folder from Run Loader");
1772
1773   if ( fRunLoader->GetAliRun() )
1774    {//if alrun already exists in folder
1775     if (fRunLoader->GetAliRun() != this )
1776      {//and is different than this - crash
1777        Fatal("AliRun","AliRun is already in Folder and it is not this object");
1778        return;//pro forma
1779      }//else do nothing
1780    }
1781   else
1782    {
1783      evfold->Add(this);//Post this AliRun to Folder
1784    }
1785
1786   TIter next(fModules);
1787   AliModule *module;
1788   while((module = (AliModule*)next()))
1789    {
1790      if (evfold) AliConfig::Instance()->Add(module,evfoldname);
1791      AliDetector* detector = dynamic_cast<AliDetector*>(module);
1792      if (detector)
1793       {
1794         AliLoader* loader = fRunLoader->GetLoader(detector);
1795         if (loader == 0x0)
1796          {
1797            Error("SetRunLoader","Can not get loader for detector %s",detector->GetName());
1798          }
1799         else
1800          {
1801            if (GetDebug()) Info("SetRunLoader","Setting loader for detector %s",detector->GetName());
1802            detector->SetLoader(loader);
1803          }
1804       }
1805    }
1806 }
1807
1808 void AliRun::AddModule(AliModule* mod)
1809 {
1810   if (mod == 0x0) return;
1811   if (strlen(mod->GetName()) == 0) return;
1812   if (GetModuleID(mod->GetName()) >= 0) return;
1813
1814   if (GetDebug()) Info("AddModule","%s",mod->GetName());
1815   if (fRunLoader == 0x0) AliConfig::Instance()->Add(mod);
1816   else AliConfig::Instance()->Add(mod,fRunLoader->GetEventFolder()->GetName());
1817
1818   Modules()->Add(mod);
1819 }