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