STEER/AliRun.cxx

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