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