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
 477   fRunLoader->Synchronize();
 478 }
 479
 480 //_______________________________________________________________________
 481 void AliRun::FlagTrack(Int_t track)
 482 {
 483   // Delegate to stack
 484   //
 485     fRunLoader->Stack()->FlagTrack(track);
 486 }
 487
 488 //_______________________________________________________________________
 489 void AliRun::EnergySummary()
 490 {
 491   //
 492   // Print summary of deposited energy
 493   //
 494
 495   Int_t ndep=0;
 496   Float_t edtot=0;
 497   Float_t ed, ed2;
 498   Int_t kn, i, left, j, id;
 499   const Float_t kzero=0;
 500   Int_t ievent=fRunLoader->GetHeader()->GetEvent()+1;
 501   //
 502   // Energy loss information
 503   if(ievent) {
 504     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 505     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 506       ed=fSummEnergy[kn];
 507       if(ed>0) {
 508         fEventEnergy[ndep]=kn;
 509         if(ievent>1) {
 510           ed=ed/ievent;
 511           ed2=fSum2Energy[kn];
 512           ed2=ed2/ievent;
 513           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
 514         } else
 515           ed2=99;
 516         fSummEnergy[ndep]=ed;
 517         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
 518         edtot+=ed;
 519         ndep++;
 520       }
 521     }
 522     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 523       left=ndep-kn*3;
 524       for(i=0;i<(3<left?3:left);i++) {
 525         j=kn*3+i;
 526         id=Int_t (fEventEnergy[j]+0.1);
 527         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 528       }
 529       printf("\n");
 530     }
 531     //
 532     // Relative energy loss in different detectors
 533     printf("******************** Relative Energy Loss per event ********************\n");
 534     printf("Total energy loss per event %10.3f GeV\n",edtot);
 535     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 536       left=ndep-kn*5;
 537       for(i=0;i<(5<left?5:left);i++) {
 538         j=kn*5+i;
 539         id=Int_t (fEventEnergy[j]+0.1);
 540         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 541       }
 542       printf("\n");
 543     }
 544     for(kn=0;kn<75;kn++) printf("*");
 545     printf("\n");
 546   }
 547   //
 548   // Reset the TArray's
 549   //  fEventEnergy.Set(0);
 550   //  fSummEnergy.Set(0);
 551   //  fSum2Energy.Set(0);
 552 }
 553
 554 //_______________________________________________________________________
 555 void AliRun::Announce() const
 556 {
 557   //
 558   // Announce the current version of AliRoot
 559   //
 560   printf("%70s",
 561          "****************************************************************\n");
 562   printf("%6s","*");printf("%64s","*\n");
 563
 564   printf("%6s","*");
 565   printf("    You are running AliRoot version NewIO\n");
 566
 567   printf("%6s","*");
 568   printf("    The cvs tag for the current program is $Name$\n");
 569
 570   printf("%6s","*");printf("%64s","*\n");
 571   printf("%70s",
 572          "****************************************************************\n");
 573 }
 574
 575 //_______________________________________________________________________
 576 AliModule *AliRun::GetModule(const char *name) const
 577 {
 578   //
 579   // Return pointer to detector from name
 580   //
 581   return dynamic_cast<AliModule*>(fModules->FindObject(name));
 582 }
 583
 584 //_______________________________________________________________________
 585 AliDetector *AliRun::GetDetector(const char *name) const
 586 {
 587   //
 588   // Return pointer to detector from name
 589   //
 590   return dynamic_cast<AliDetector*>(fModules->FindObject(name));
 591 }
 592
 593 //_______________________________________________________________________
 594 Int_t AliRun::GetModuleID(const char *name) const
 595 {
 596   //
 597   // Return galice internal detector identifier from name
 598   //
 599   Int_t i=-1;
 600   TObject *mod=fModules->FindObject(name);
 601   if(mod) i=fModules->IndexOf(mod);
 602   return i;
 603 }
 604
 605 //_______________________________________________________________________
 606 Int_t AliRun::GetEvent(Int_t event)
 607 {
 608 //
 609 // Reloads data containers in folders # event
 610 // Set branch addresses
 611 //
 612   if (fRunLoader == 0x0)
 613    {
 614      Error("GetEvent","RunLoader is not set. Can not load data.");
 615      return -1;
 616    }
 617 /*****************************************/
 618 /****   P R E    R E L O A D I N G    ****/
 619 /*****************************************/
 620 // Reset existing structures
 621   ResetHits();
 622   ResetTrackReferences();
 623   ResetDigits();
 624   ResetSDigits();
 625
 626 /*****************************************/
 627 /****       R  E  L  O  A  D          ****/
 628 /*****************************************/
 629
 630   fRunLoader->GetEvent(event);
 631
 632 /*****************************************/
 633 /****  P O S T    R E L O A D I N G   ****/
 634 /*****************************************/
 635
 636   // Set Trees branch addresses
 637   TIter next(fModules);
 638   AliModule *detector;
 639   while((detector = dynamic_cast<AliModule*>(next())))
 640    {
 641      detector->SetTreeAddress();
 642    }
 643
 644   return fRunLoader->GetHeader()->GetNtrack();
 645 }
 646
 647 //_______________________________________________________________________
 648 TGeometry *AliRun::GetGeometry()
 649 {
 650   //
 651   // Import Alice geometry from current file
 652   // Return pointer to geometry object
 653   //
 654   if (!fGeometry) fGeometry = dynamic_cast<TGeometry*>(gDirectory->Get("AliceGeom"));
 655   //
 656   // Unlink and relink nodes in detectors
 657   // This is bad and there must be a better way...
 658   //
 659
 660   TIter next(fModules);
 661   AliModule *detector;
 662   while((detector = dynamic_cast<AliModule*>(next()))) {
 663     TList *dnodes=detector->Nodes();
 664     Int_t j;
 665     TNode *node, *node1;
 666     for ( j=0; j<dnodes->GetSize(); j++) {
 667       node = dynamic_cast<TNode*>(dnodes->At(j));
 668       node1 = fGeometry->GetNode(node->GetName());
 669       dnodes->Remove(node);
 670       dnodes->AddAt(node1,j);
 671     }
 672   }
 673   return fGeometry;
 674 }
 675
 676 //_______________________________________________________________________
 677 Int_t AliRun::GetPrimary(Int_t track) const
 678 {
 679   //
 680   // return number of primary that has generated track
 681   //
 682     return fRunLoader->Stack()->GetPrimary(track);
 683 }
 684
 685 //_______________________________________________________________________
 686 void AliRun::MediaTable()
 687 {
 688   //
 689   // Built media table to get from the media number to
 690   // the detector id
 691   //
 692
 693   Int_t kz, nz, idt, lz, i, k, ind;
 694   //  Int_t ibeg;
 695   TObjArray &dets = *gAlice->Detectors();
 696   AliModule *det;
 697   //
 698   // For all detectors
 699   for (kz=0;kz<fNdets;kz++) {
 700     // If detector is defined
 701     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
 702         TArrayI &idtmed = *(det->GetIdtmed());
 703         for(nz=0;nz<100;nz++) {
 704         // Find max and min material number
 705         if((idt=idtmed[nz])) {
 706           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 707           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 708         }
 709       }
 710       if(det->LoMedium() > det->HiMedium()) {
 711         det->LoMedium() = 0;
 712         det->HiMedium() = 0;
 713       } else {
 714         if(det->HiMedium() > fImedia->GetSize()) {
 715           Error("MediaTable","Increase fImedia from %d to %d",
 716                 fImedia->GetSize(),det->HiMedium());
 717           return;
 718         }
 719         // Tag all materials in rage as belonging to detector kz
 720         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 721           (*fImedia)[lz]=kz;
 722         }
 723       }
 724     }
 725   }
 726   //
 727   // Print summary table
 728   printf(" Traking media ranges:\n");
 729   for(i=0;i<(fNdets-1)/6+1;i++) {
 730     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
 731       ind=i*6+k;
 732       det=dynamic_cast<AliModule*>(dets[ind]);
 733       if(det)
 734         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 735                det->HiMedium());
 736       else
 737         printf(" %6s: %3d -> %3d;","NULL",0,0);
 738     }
 739     printf("\n");
 740   }
 741 }
 742
 743 //_______________________________________________________________________
 744 void AliRun::SetGenerator(AliGenerator *generator)
 745 {
 746   //
 747   // Load the event generator
 748   //
 749   if(!fGenerator) fGenerator = generator;
 750 }
 751
 752 //_______________________________________________________________________
 753 void AliRun::ResetGenerator(AliGenerator *generator)
 754 {
 755   //
 756   // Load the event generator
 757   //
 758   if(fGenerator)
 759     if(generator)
 760       Warning("ResetGenerator","Replacing generator %s with %s\n",
 761               fGenerator->GetName(),generator->GetName());
 762     else
 763       Warning("ResetGenerator","Replacing generator %s with NULL\n",
 764               fGenerator->GetName());
 765   fGenerator = generator;
 766 }
 767
 768 //_______________________________________________________________________
 769 void AliRun::SetTransPar(const char *filename)
 770 {
 771   //
 772   // Sets the file name for transport parameters
 773   //
 774   fTransParName = filename;
 775 }
 776
 777 //_______________________________________________________________________
 778 void AliRun::SetBaseFile(const char *filename)
 779 {
 780   fBaseFileName = filename;
 781 }
 782
 783 //_______________________________________________________________________
 784 void AliRun::ReadTransPar()
 785 {
 786   //
 787   // Read filename to set the transport parameters
 788   //
 789
 790
 791   const Int_t kncuts=10;
 792   const Int_t knflags=11;
 793   const Int_t knpars=kncuts+knflags;
 794   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 795                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 796                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 797                                "MULS","PAIR","PHOT","RAYL"};
 798   char line[256];
 799   char detName[7];
 800   char* filtmp;
 801   Float_t cut[kncuts];
 802   Int_t flag[knflags];
 803   Int_t i, itmed, iret, ktmed, kz;
 804   FILE *lun;
 805   //
 806   // See whether the file is there
 807   filtmp=gSystem->ExpandPathName(fTransParName.Data());
 808   lun=fopen(filtmp,"r");
 809   delete [] filtmp;
 810   if(!lun) {
 811     Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
 812     return;
 813   }
 814   //
 815   if(fDebug) {
 816     printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 817     printf(" *%59s\n","*");
 818     printf(" *       Please check carefully what you are doing!%10s\n","*");
 819     printf(" *%59s\n","*");
 820   }
 821   //
 822   while(1) {
 823     // Initialise cuts and flags
 824     for(i=0;i<kncuts;i++) cut[i]=-99;
 825     for(i=0;i<knflags;i++) flag[i]=-99;
 826     itmed=0;
 827     for(i=0;i<256;i++) line[i]='\0';
 828     // Read up to the end of line excluded
 829     iret=fscanf(lun,"%[^\n]",line);
 830     if(iret<0) {
 831       //End of file
 832       fclose(lun);
 833       if(fDebug){
 834         printf(" *%59s\n","*");
 835         printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 836       }
 837       return;
 838     }
 839     // Read the end of line
 840     fscanf(lun,"%*c");
 841     if(!iret) continue;
 842     if(line[0]=='*') continue;
 843     // Read the numbers
 844     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",
 845                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
 846                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
 847                 &flag[8],&flag[9],&flag[10]);
 848     if(!iret) continue;
 849     if(iret<0) {
 850       //reading error
 851       Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
 852       continue;
 853     }
 854     // Check that the module exist
 855     AliModule *mod = GetModule(detName);
 856     if(mod) {
 857       // Get the array of media numbers
 858       TArrayI &idtmed = *mod->GetIdtmed();
 859       // Check that the tracking medium code is valid
 860       if(0<=itmed && itmed < 100) {
 861         ktmed=idtmed[itmed];
 862         if(!ktmed) {
 863           Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
 864           continue;
 865         }
 866         // Set energy thresholds
 867         for(kz=0;kz<kncuts;kz++) {
 868           if(cut[kz]>=0) {
 869             if(fDebug) printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
 870                    kpars[kz],cut[kz],itmed,mod->GetName());
 871             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
 872           }
 873         }
 874         // Set transport mechanisms
 875         for(kz=0;kz<knflags;kz++) {
 876           if(flag[kz]>=0) {
 877             if(fDebug) printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
 878                    kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
 879             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
 880           }
 881         }
 882       } else {
 883         Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
 884         continue;
 885       }
 886     } else {
 887       if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
 888       continue;
 889     }
 890   }
 891 }
 892 //_____________________________________________________________________________
 893
 894 void AliRun::BeginEvent()
 895 {
 896   //
 897   // Clean-up previous event
 898   // Energy scores
 899   if (GetDebug())
 900    {
 901      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 902      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 903      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 904      Info("BeginEvent","          BEGINNING EVENT               ");
 905      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 906      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 907      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 908    }
 909
 910   /*******************************/
 911   /*   Clean after eventual      */
 912   /*   previous event            */
 913   /*******************************/
 914
 915
 916   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
 917   fRunLoader->SetEventNumber(++fEventNrInRun);// sets new files, cleans the previous event stuff, if necessary, etc.,
 918   if (GetDebug()) Info("BeginEvent","EventNr is %d",fEventNrInRun);
 919
 920   fEventEnergy.Reset();
 921     // Clean detector information
 922
 923   if (fRunLoader->Stack())
 924     fRunLoader->Stack()->Reset();//clean stack -> tree is unloaded
 925   else
 926     fRunLoader->MakeStack();//or make a new one
 927
 928   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(K)");
 929   fRunLoader->MakeTree("K");
 930   if (GetDebug()) Info("BeginEvent","  gMC->SetStack(fRunLoader->Stack())");
 931   gMC->SetStack(fRunLoader->Stack());//Was in InitMC - but was moved here
 932                                      //because we don't have guarantee that
 933                                      //stack pointer is not going to change from event to event
 934                          //since it bellobgs to header and is obtained via RunLoader
 935   //
 936   //  Reset all Detectors & kinematics & make/reset trees
 937   //
 938
 939   fRunLoader->GetHeader()->Reset(fRun,fEvent,fEventNrInRun);
 940 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 941
 942   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTrackRefsContainer()");
 943   fRunLoader->MakeTrackRefsContainer();//for insurance
 944
 945   if (GetDebug()) Info("BeginEvent","  ResetHits()");
 946   ResetHits();
 947   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(H)");
 948   fRunLoader->MakeTree("H");
 949
 950   //
 951   if(fLego)
 952    {
 953     fLego->BeginEvent();
 954     return;
 955    }
 956
 957   //create new branches and SetAdresses
 958   TIter next(fModules);
 959   AliModule *detector;
 960   while((detector = (AliModule*)next()))
 961    {
 962     if (GetDebug()) Info("BeginEvent","  %s->MakeBranch(H)",detector->GetName());
 963     detector->MakeBranch("H");
 964     if (GetDebug()) Info("BeginEvent","  %s->MakeBranchTR()",detector->GetName());
 965     detector->MakeBranchTR();
 966     if (GetDebug()) Info("BeginEvent","  %s->SetTreeAddress()",detector->GetName());
 967     detector->SetTreeAddress();
 968    }
 969   // make branch for AliRun track References
 970   TTree * treeTR = fRunLoader->TreeTR();
 971   if (treeTR){
 972     // make branch for central track references
 973     if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
 974     TBranch *branch;
 975     branch = treeTR->Branch("AliRun",&fTrackReferences);
 976     branch->SetAddress(&fTrackReferences);
 977   }
 978   //
 979 }
 980
 981 //_______________________________________________________________________
 982 TParticle* AliRun::Particle(Int_t i) const
 983 {
 984   if (fRunLoader)
 985    if (fRunLoader->Stack())
 986     return fRunLoader->Stack()->Particle(i);
 987   return 0x0;
 988 }
 989
 990 //_______________________________________________________________________
 991 void AliRun::ResetDigits()
 992 {
 993   //
 994   //  Reset all Detectors digits
 995   //
 996   TIter next(fModules);
 997   AliModule *detector;
 998   while((detector = dynamic_cast<AliModule*>(next()))) {
 999      detector->ResetDigits();
1000   }
1001 }
1002
1003 //_______________________________________________________________________
1004 void AliRun::ResetSDigits()
1005 {
1006   //
1007   //  Reset all Detectors digits
1008   //
1009   TIter next(fModules);
1010   AliModule *detector;
1011   while((detector = dynamic_cast<AliModule*>(next()))) {
1012      detector->ResetSDigits();
1013   }
1014 }
1015
1016 //_______________________________________________________________________
1017 void AliRun::ResetHits()
1018 {
1019   //
1020   //  Reset all Detectors hits
1021   //
1022   TIter next(fModules);
1023   AliModule *detector;
1024   while((detector = dynamic_cast<AliModule*>(next()))) {
1025      detector->ResetHits();
1026   }
1027 }
1028 //_______________________________________________________________________
1029
1030 void  AliRun::AddTrackReference(Int_t label){
1031   //
1032   // add a trackrefernce to the list
1033   if (!fTrackReferences) {
1034     cerr<<"Container trackrefernce not active\n";
1035     return;
1036   }
1037   Int_t nref = fTrackReferences->GetEntriesFast();
1038   TClonesArray &lref = *fTrackReferences;
1039   new(lref[nref]) AliTrackReference(label);
1040 }
1041
1042
1043
1044 void AliRun::ResetTrackReferences()
1045 {
1046   //
1047   //  Reset all  references
1048   //
1049   if (fTrackReferences)   fTrackReferences->Clear();
1050
1051   TIter next(fModules);
1052   AliModule *detector;
1053   while((detector = dynamic_cast<AliModule*>(next()))) {
1054      detector->ResetTrackReferences();
1055   }
1056 }
1057
1058 void AliRun::RemapTrackReferencesIDs(Int_t *map)
1059 {
1060   //
1061   // Remapping track reference
1062   // Called at finish primary
1063   //
1064   if (!fTrackReferences) return;
1065   for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
1066     AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
1067     if (ref) {
1068       Int_t newID = map[ref->GetTrack()];
1069       if (newID>=0) ref->SetTrack(newID);
1070       else {
1071         //ref->SetTrack(-1);
1072         ref->SetBit(kNotDeleted,kFALSE);
1073         fTrackReferences->RemoveAt(i);
1074       }
1075     }
1076   }
1077   fTrackReferences->Compress();
1078 }
1079
1080
1081
1082 //_______________________________________________________________________
1083
1084 void AliRun::ResetPoints()
1085 {
1086   //
1087   // Reset all Detectors points
1088   //
1089   TIter next(fModules);
1090   AliModule *detector;
1091   while((detector = dynamic_cast<AliModule*>(next()))) {
1092      detector->ResetPoints();
1093   }
1094 }
1095 //_______________________________________________________________________
1096
1097 void AliRun::InitMC(const char *setup)
1098 {
1099   //
1100   // Initialize the Alice setup
1101   //
1102   Announce();
1103
1104   if(fInitDone) {
1105     Warning("Init","Cannot initialise AliRun twice!\n");
1106     return;
1107   }
1108
1109   gROOT->LoadMacro(setup);
1110   gInterpreter->ProcessLine(fConfigFunction.Data());
1111
1112   // Register MC in configuration
1113   AliConfig::Instance()->Add(gMC);
1114
1115   InitLoaders();
1116
1117   fRunLoader->MakeTree("E");
1118   fRunLoader->LoadKinematics("RECREATE");
1119   fRunLoader->LoadTrackRefs("RECREATE");
1120   fRunLoader->LoadHits("all","RECREATE");
1121
1122
1123   fRunLoader->CdGAFile();
1124
1125   gMC->DefineParticles();  //Create standard MC particles
1126   AliPDG::AddParticlesToPdgDataBase();
1127
1128   TObject *objfirst, *objlast;
1129
1130   fNdets = fModules->GetLast()+1;
1131
1132   //
1133   //=================Create Materials and geometry
1134   gMC->Init();
1135
1136   // Added also after in case of interactive initialisation of modules
1137   fNdets = fModules->GetLast()+1;
1138
1139   TIter next(fModules);
1140   AliModule *detector;
1141   while((detector = dynamic_cast<AliModule*>(next())))
1142    {
1143      objlast = gDirectory->GetList()->Last();
1144
1145      // Add Detector histograms in Detector list of histograms
1146      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
1147      else         objfirst = gDirectory->GetList()->First();
1148      while (objfirst)
1149       {
1150         detector->Histograms()->Add(objfirst);
1151         objfirst = gDirectory->GetList()->After(objfirst);
1152       }
1153    }
1154    ReadTransPar(); //Read the cuts for all materials
1155
1156    MediaTable(); //Build the special IMEDIA table
1157
1158    //Initialise geometry deposition table
1159    fEventEnergy.Set(gMC->NofVolumes()+1);
1160    fSummEnergy.Set(gMC->NofVolumes()+1);
1161    fSum2Energy.Set(gMC->NofVolumes()+1);
1162
1163    //Compute cross-sections
1164    gMC->BuildPhysics();
1165
1166    //Write Geometry object to current file.
1167    fRunLoader->WriteGeometry();
1168
1169    fInitDone = kTRUE;
1170
1171    fMCQA = new AliMCQA(fNdets);
1172
1173    //
1174    // Save stuff at the beginning of the file to avoid file corruption
1175    Write();
1176    fEventNrInRun = -1; //important - we start Begin event from increasing current number in run
1177 }
1178
1179 //_______________________________________________________________________
1180
1181 void AliRun::RunMC(Int_t nevent, const char *setup)
1182 {
1183   //
1184   // Main function to be called to process a galice run
1185   // example
1186   //    Root > gAlice.Run();
1187   // a positive number of events will cause the finish routine
1188   // to be called
1189   //
1190   fEventsPerRun = nevent;
1191   // check if initialisation has been done
1192   if (!fInitDone) InitMC(setup);
1193
1194   // Create the Root Tree with one branch per detector
1195   //Hits moved to begin event -> now we are crating separate tree for each event
1196
1197   gMC->ProcessRun(nevent);
1198
1199   // End of this run, close files
1200   if(nevent>0) FinishRun();
1201 }
1202
1203 //_______________________________________________________________________
1204 void AliRun::RunReco(const char *selected, Int_t first, Int_t last)
1205 {
1206   //
1207   // Main function to be called to reconstruct Alice event
1208   //
1209    Int_t nev = fRunLoader->GetNumberOfEvents();
1210    if (GetDebug()) Info("RunReco","Found %d events",nev);
1211    Int_t nFirst = first;
1212    Int_t nLast  = (last < 0)? nev : last;
1213
1214    for (Int_t nevent = nFirst; nevent <= nLast; nevent++) {
1215      if (GetDebug()) Info("RunReco","Processing event %d",nevent);
1216      GetEvent(nevent);
1217      Digits2Reco(selected);
1218    }
1219 }
1220
1221 //_______________________________________________________________________
1222
1223 void AliRun::Hits2Digits(const char *selected)
1224 {
1225
1226    // Convert Hits to sumable digits
1227    //
1228    for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
1229      GetEvent(nevent);
1230      Hits2SDigits(selected);
1231      SDigits2Digits(selected);
1232    }
1233 }
1234
1235
1236 //_______________________________________________________________________
1237
1238 void AliRun::Tree2Tree(Option_t *option, const char *selected)
1239 {
1240   //
1241   // Function to transform the content of
1242   //
1243   // - TreeH to TreeS (option "S")
1244   // - TreeS to TreeD (option "D")
1245   // - TreeD to TreeR (option "R")
1246   //
1247   // If multiple options are specified ("SDR"), transformation will be done in sequence for
1248   // selected detector and for all detectors if none is selected (detector string
1249   // can contain blank separated list of detector names).
1250
1251
1252    const char *oS = strstr(option,"S");
1253    const char *oD = strstr(option,"D");
1254    const char *oR = strstr(option,"R");
1255
1256    TObjArray *detectors = Detectors();
1257
1258    TIter next(detectors);
1259
1260    AliDetector *detector = 0;
1261
1262    while((detector = dynamic_cast<AliDetector*>(next()))) {
1263      if (selected)
1264        if (strcmp(detector->GetName(),selected)) continue;
1265      if (detector->IsActive())
1266       {
1267
1268        AliLoader* loader = detector->GetLoader();
1269        if (loader == 0x0) continue;
1270
1271        if (oS)
1272         {
1273           if (GetDebug()) Info("Tree2Tree","Processing Hits2SDigits for %s ...",detector->GetName());
1274           loader->LoadHits("read");
1275           if (loader->TreeS() == 0x0) loader->MakeTree("S");
1276           detector->MakeBranch(option);
1277           detector->SetTreeAddress();
1278           detector->Hits2SDigits();
1279           loader->UnloadHits();
1280           loader->UnloadSDigits();
1281         }
1282        if (oD)
1283         {
1284           if (GetDebug()) Info("Tree2Tree","Processing SDigits2Digits for %s ...",detector->GetName());
1285           loader->LoadSDigits("read");
1286           if (loader->TreeD() == 0x0) loader->MakeTree("D");
1287           detector->MakeBranch(option);
1288           detector->SetTreeAddress();
1289           detector->SDigits2Digits();
1290           loader->UnloadSDigits();
1291           loader->UnloadDigits();
1292         }
1293        if (oR)
1294         {
1295           if (GetDebug()) Info("Tree2Tree","Processing Digits2Reco for %s ...",detector->GetName());
1296           loader->LoadDigits("read");
1297           if (loader->TreeR() == 0x0) loader->MakeTree("R");
1298           detector->MakeBranch(option);
1299           detector->SetTreeAddress();
1300           detector->Digits2Reco();
1301           loader->UnloadDigits();
1302           loader->UnloadRecPoints();
1303
1304         }
1305      }
1306    }
1307 }
1308
1309 //_______________________________________________________________________
1310 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1311                      Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1312                      Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
1313 {
1314   //
1315   // Generates lego plots of:
1316   //    - radiation length map phi vs theta
1317   //    - radiation length map phi vs eta
1318   //    - interaction length map
1319   //    - g/cm2 length map
1320   //
1321   //  ntheta    bins in theta, eta
1322   //  themin    minimum angle in theta (degrees)
1323   //  themax    maximum angle in theta (degrees)
1324   //  nphi      bins in phi
1325   //  phimin    minimum angle in phi (degrees)
1326   //  phimax    maximum angle in phi (degrees)
1327   //  rmin      minimum radius
1328   //  rmax      maximum radius
1329   //
1330   //
1331   //  The number of events generated = ntheta*nphi
1332   //  run input parameters in macro setup (default="Config.C")
1333   //
1334   //  Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1335   //Begin_Html
1336   /*
1337     <img src="picts/AliRunLego1.gif">
1338   */
1339   //End_Html
1340   //Begin_Html
1341   /*
1342     <img src="picts/AliRunLego2.gif">
1343   */
1344   //End_Html
1345   //Begin_Html
1346   /*
1347     <img src="picts/AliRunLego3.gif">
1348   */
1349   //End_Html
1350   //
1351
1352   // check if initialisation has been done
1353   if (!fInitDone) InitMC(setup);
1354   //Save current generator
1355   AliGenerator *gen=Generator();
1356
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   if (id < 0) return;
1542
1543   //
1544   // --- If lego option, do it and leave
1545   if (fLego)
1546     fLego->StepManager();
1547   else {
1548     Int_t copy;
1549     //Update energy deposition tables
1550     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
1551     //
1552     // write tracke reference for track which is dissapearing - MI
1553     if (gMC->IsTrackDisappeared()) {
1554       if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber());
1555     }
1556
1557     //Call the appropriate stepping routine;
1558     AliModule *det = dynamic_cast<AliModule*>(fModules->At(id));
1559     if(det && det->StepManagerIsEnabled()) {
1560       fMCQA->StepManager(id);
1561       det->StepManager();
1562     }
1563   }
1564 }
1565
1566 //_______________________________________________________________________
1567 void AliRun::PostTrack()
1568 {
1569      TObjArray &dets = *fModules;
1570      AliModule *module;
1571
1572      for(Int_t i=0; i<=fNdets; i++)
1573        if((module = dynamic_cast<AliModule*>(dets[i])))
1574          module->PostTrack();
1575 }
1576
1577 //_______________________________________________________________________
1578 void AliRun::FinishPrimary()
1579 {
1580   //
1581   // Called  at the end of each primary track
1582   //
1583
1584   //  static Int_t count=0;
1585   //  const Int_t times=10;
1586   // This primary is finished, purify stack
1587   fRunLoader->Stack()->PurifyKine();
1588
1589   TIter next(fModules);
1590   AliModule *detector;
1591   while((detector = dynamic_cast<AliModule*>(next()))) {
1592     detector->FinishPrimary();
1593     if(detector->GetLoader())
1594      {
1595        detector->GetLoader()->TreeH()->Fill();
1596      }
1597   }
1598
1599   // Write out track references if any
1600   if (fRunLoader->TreeTR())
1601    {
1602     fRunLoader->TreeTR()->Fill();
1603    }
1604 }
1605
1606 //_______________________________________________________________________
1607 void AliRun::FinishEvent()
1608 {
1609   //
1610   // Called at the end of the event.
1611   //
1612
1613   //
1614   if(fLego) fLego->FinishEvent();
1615
1616   TIter next(fModules);
1617   AliModule *detector;
1618   while((detector = dynamic_cast<AliModule*>(next()))) {
1619     detector->FinishEvent();
1620   }
1621
1622   //Update the energy deposit tables
1623   Int_t i;
1624   for(i=0;i<fEventEnergy.GetSize();i++)
1625    {
1626     fSummEnergy[i]+=fEventEnergy[i];
1627     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
1628    }
1629
1630   AliHeader* header = fRunLoader->GetHeader();
1631   AliStack* stack = fRunLoader->Stack();
1632   if ( (header == 0x0) || (stack == 0x0) )
1633    {//check if we got header and stack. If not cry and exit aliroot
1634     Fatal("AliRun","Can not get the stack or header from LOADER");
1635     return;//never reached
1636    }
1637   // Update Header information
1638   header->SetNprimary(stack->GetNprimary());
1639   header->SetNtrack(stack->GetNtrack());
1640
1641
1642   // Write out the kinematics
1643   stack->FinishEvent();
1644
1645   // Write out the event Header information
1646   TTree* treeE = fRunLoader->TreeE();
1647   if (treeE)
1648    {
1649       header->SetStack(stack);
1650       treeE->Fill();
1651    }
1652   else
1653    {
1654     Error("FinishEvent","Can not get TreeE from RL");
1655    }
1656
1657   fRunLoader->WriteKinematics("OVERWRITE");
1658   fRunLoader->WriteTrackRefs("OVERWRITE");
1659   fRunLoader->WriteHits("OVERWRITE");
1660
1661   if (GetDebug())
1662    {
1663      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1664      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1665      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1666      Info("FinishEvent","          FINISHING EVENT               ");
1667      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1668      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1669      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1670    }
1671 }
1672
1673 //_______________________________________________________________________
1674 void AliRun::Field(const Double_t* x, Double_t *b) const
1675 {
1676    Float_t xfloat[3];
1677    for (Int_t i=0; i<3; i++) xfloat[i] = x[i];
1678
1679    if (Field()) {
1680          Float_t bfloat[3];
1681          Field()->Field(xfloat,bfloat);
1682          for (Int_t j=0; j<3; j++) b[j] = bfloat[j];
1683    }
1684    else {
1685          printf("No mag field defined!\n");
1686          b[0]=b[1]=b[2]=0.;
1687    }
1688
1689 }
1690
1691 //
1692 // End of MC Application
1693 //
1694
1695 //_______________________________________________________________________
1696 void AliRun::Streamer(TBuffer &R__b)
1697 {
1698   // Stream an object of class AliRun.
1699
1700   if (R__b.IsReading()) {
1701     if (!gAlice) gAlice = this;
1702     AliRun::Class()->ReadBuffer(R__b, this);
1703     gROOT->GetListOfBrowsables()->Add(this,"Run");
1704
1705     gRandom = fRandom;
1706   } else {
1707     AliRun::Class()->WriteBuffer(R__b, this);
1708   }
1709 }
1710
1711
1712 //_______________________________________________________________________
1713 Int_t AliRun::GetCurrentTrackNumber() const {
1714   //
1715   // Returns current track
1716   //
1717   return fRunLoader->Stack()->GetCurrentTrackNumber();
1718 }
1719
1720 //_______________________________________________________________________
1721 Int_t AliRun::GetNtrack() const {
1722   //
1723   // Returns number of tracks in stack
1724   //
1725   return fRunLoader->Stack()->GetNtrack();
1726 }
1727 //_______________________________________________________________________
1728
1729 //_______________________________________________________________________
1730 TObjArray* AliRun::Particles() const {
1731   //
1732   // Returns pointer to Particles array
1733   //
1734   if (fRunLoader)
1735    if (fRunLoader->Stack())
1736     return fRunLoader->Stack()->Particles();
1737   return 0x0;
1738 }
1739
1740 //___________________________________________________________________________
1741
1742 //_______________________________________________________________________
1743 void AliRun::SetGenEventHeader(AliGenEventHeader* header)
1744 {
1745   fRunLoader->GetHeader()->SetGenEventHeader(header);
1746 }
1747
1748 //___________________________________________________________________________
1749
1750 Int_t AliRun::GetEvNumber() const
1751 {
1752 //Returns number of current event
1753   if (fRunLoader == 0x0)
1754    {
1755      Error("GetEvent","RunLoader is not set. Can not load data.");
1756      return -1;
1757    }
1758
1759   return fRunLoader->GetEventNumber();
1760 }
1761
1762 void AliRun::SetRunLoader(AliRunLoader* rloader)
1763 {
1764   fRunLoader = rloader;
1765   if (fRunLoader == 0x0) return;
1766
1767   TString evfoldname;
1768   TFolder* evfold = fRunLoader->GetEventFolder();
1769   if (evfold) evfoldname = evfold->GetName();
1770   else Warning("SetRunLoader","Did not get Event Folder from Run Loader");
1771
1772   if ( fRunLoader->GetAliRun() )
1773    {//if alrun already exists in folder
1774     if (fRunLoader->GetAliRun() != this )
1775      {//and is different than this - crash
1776        Fatal("AliRun","AliRun is already in Folder and it is not this object");
1777        return;//pro forma
1778      }//else do nothing
1779    }
1780   else
1781    {
1782      evfold->Add(this);//Post this AliRun to Folder
1783    }
1784
1785   TIter next(fModules);
1786   AliModule *module;
1787   while((module = (AliModule*)next()))
1788    {
1789      if (evfold) AliConfig::Instance()->Add(module,evfoldname);
1790      AliDetector* detector = dynamic_cast<AliDetector*>(module);
1791      if (detector)
1792       {
1793         AliLoader* loader = fRunLoader->GetLoader(detector);
1794         if (loader == 0x0)
1795          {
1796            Error("SetRunLoader","Can not get loader for detector %s",detector->GetName());
1797          }
1798         else
1799          {
1800            if (GetDebug()) Info("SetRunLoader","Setting loader for detector %s",detector->GetName());
1801            detector->SetLoader(loader);
1802          }
1803       }
1804    }
1805 }
1806
1807 void AliRun::AddModule(AliModule* mod)
1808 {
1809   if (mod == 0x0) return;
1810   if (strlen(mod->GetName()) == 0) return;
1811   if (GetModuleID(mod->GetName()) >= 0) return;
1812
1813   if (GetDebug()) Info("AddModule","%s",mod->GetName());
1814   if (fRunLoader == 0x0) AliConfig::Instance()->Add(mod);
1815   else AliConfig::Instance()->Add(mod,fRunLoader->GetEventFolder()->GetName());
1816
1817   Modules()->Add(mod);
1818 }