STEER/AliMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // This class is extracted from the AliRun class
  19 // and contains all the MC-related functionality
  20 // The number of dependencies has to be reduced...
  21 // Author: F.Carminati
  22 //         Federico.Carminati@cern.ch
  23
  24 #include <RVersion.h>
  25 #include <TBrowser.h>
  26 #include <TClonesArray.h>
  27 #include <TGeoManager.h>
  28 #include <TStopwatch.h>
  29 #include <TSystem.h>
  30 #include <TVirtualMC.h>
  31 #include <TParticle.h>
  32 #include <TROOT.h>
  33
  34 #include "AliLog.h"
  35 #include "AliDetector.h"
  36 #include "AliGenerator.h"
  37 #include "AliHeader.h"
  38 #include "AliLego.h"
  39 #include "AliMC.h"
  40 #include "AliMCQA.h"
  41 #include "AliRun.h"
  42 #include "AliHit.h"
  43 #include "AliStack.h"
  44 #include "AliMagF.h"
  45 #include "AliTrackReference.h"
  46 #include "AliSimulation.h"
  47 #include "AliGeomManager.h"
  48 #include "AliCDBManager.h"
  49 #include "AliCDBStorage.h"
  50 #include "AliCDBEntry.h"
  51
  52
  53 ClassImp(AliMC)
  54
  55 //_______________________________________________________________________
  56 AliMC::AliMC() :
  57   fGenerator(0),
  58   fEventEnergy(0),
  59   fSummEnergy(0),
  60   fSum2Energy(0),
  61   fTrRmax(1.e10),
  62   fTrZmax(1.e10),
  63   fRDecayMax(1.e10),
  64   fRDecayMin(-1.),
  65   fDecayPdg(0),
  66   fImedia(0),
  67   fTransParName("\0"),
  68   fMCQA(0),
  69   fHitLists(0),
  70   fTrackReferences(0)
  71
  72 {
  73   //default constructor
  74   DecayLimits();
  75 }
  76
  77 //_______________________________________________________________________
  78 AliMC::AliMC(const char *name, const char *title) :
  79   TVirtualMCApplication(name, title),
  80   fGenerator(0),
  81   fEventEnergy(0),
  82   fSummEnergy(0),
  83   fSum2Energy(0),
  84   fTrRmax(1.e10),
  85   fTrZmax(1.e10),
  86   fRDecayMax(1.e10),
  87   fRDecayMin(-1.),
  88   fDecayPdg(0),
  89   fImedia(new TArrayI(1000)),
  90   fTransParName("\0"),
  91   fMCQA(0),
  92   fHitLists(new TList()),
  93   fTrackReferences(new TClonesArray("AliTrackReference", 100))
  94 {
  95   //constructor
  96   // Set transport parameters
  97   SetTransPar();
  98   DecayLimits();
  99   // Prepare the tracking medium lists
 100   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
 101 }
 102
 103 //_______________________________________________________________________
 104 AliMC::AliMC(const AliMC &mc) :
 105   TVirtualMCApplication(mc),
 106   fGenerator(0),
 107   fEventEnergy(0),
 108   fSummEnergy(0),
 109   fSum2Energy(0),
 110   fTrRmax(1.e10),
 111   fTrZmax(1.e10),
 112   fRDecayMax(1.e10),
 113   fRDecayMin(-1.),
 114   fDecayPdg(0),
 115   fImedia(0),
 116   fTransParName("\0"),
 117   fMCQA(0),
 118   fHitLists(0),
 119   fTrackReferences(0)
 120 {
 121   //
 122   // Copy constructor for AliMC
 123   //
 124   mc.Copy(*this);
 125 }
 126
 127 //_______________________________________________________________________
 128 AliMC::~AliMC()
 129 {
 130   //destructor
 131   delete fGenerator;
 132   delete fImedia;
 133   delete fMCQA;
 134   delete fHitLists;
 135   // Delete track references
 136   if (fTrackReferences) {
 137     fTrackReferences->Delete();
 138     delete fTrackReferences;
 139     fTrackReferences     = 0;
 140   }
 141
 142 }
 143
 144 //_______________________________________________________________________
 145 void AliMC::Copy(TObject &) const
 146 {
 147   //dummy Copy function
 148   AliFatal("Not implemented!");
 149 }
 150
 151 //_______________________________________________________________________
 152 void  AliMC::ConstructGeometry()
 153 {
 154   //
 155   // Either load geometry from file or create it through usual
 156   // loop on detectors. In the first case the method
 157   // AliModule::CreateMaterials() only builds fIdtmed and is postponed
 158   // at InitGeometry().
 159   //
 160
 161   if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file
 162     if(gAlice->IsGeomFromCDB()){
 163       AliInfo("Loading geometry from CDB default storage");
 164       AliCDBPath path("GRP","Geometry","Data");
 165       AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
 166       if(!entry) AliFatal("Unable to load geometry from CDB!");
 167       entry->SetOwner(0);
 168       gGeoManager = (TGeoManager*) entry->GetObject();
 169       if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
 170     }else{
 171       // Load geometry
 172       const char *geomfilename = gAlice->GetGeometryFileName();
 173       if(gSystem->ExpandPathName(geomfilename)){
 174         AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
 175         TGeoManager::Import(geomfilename);
 176       }else{
 177         AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
 178         return;
 179       }
 180     }
 181   }else{
 182     // Create modules, materials, geometry
 183     TStopwatch stw;
 184     TIter next(gAlice->Modules());
 185     AliModule *detector;
 186     AliDebug(1, "Geometry creation:");
 187     while((detector = dynamic_cast<AliModule*>(next()))) {
 188       stw.Start();
 189       // Initialise detector materials and geometry
 190       detector->CreateMaterials();
 191       detector->CreateGeometry();
 192       AliInfo(Form("%10s R:%.2fs C:%.2fs",
 193                    detector->GetName(),stw.RealTime(),stw.CpuTime()));
 194     }
 195   }
 196
 197 }
 198
 199 //_______________________________________________________________________
 200 Bool_t  AliMC::MisalignGeometry()
 201 {
 202 // Call misalignment code if AliSimulation object was defined.
 203
 204    if(!gAlice->IsRootGeometry()){
 205      //Set alignable volumes for the whole geometry
 206      SetAllAlignableVolumes();
 207    }
 208    // Misalign geometry via AliSimulation instance
 209    if (!AliSimulation::GetInstance()) return kFALSE;
 210    AliGeomManager::SetGeometry(gGeoManager);
 211    return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader());
 212 }
 213
 214 //_______________________________________________________________________
 215 void  AliMC::ConstructOpGeometry()
 216 {
 217   //
 218   // Loop all detector modules and call DefineOpticalProperties() method
 219   //
 220
 221   TIter next(gAlice->Modules());
 222   AliModule *detector;
 223   AliInfo("Optical properties definition");
 224   while((detector = dynamic_cast<AliModule*>(next()))) {
 225     // Initialise detector optical properties
 226     detector->DefineOpticalProperties();
 227   }
 228 }
 229
 230 //_______________________________________________________________________
 231 void  AliMC::InitGeometry()
 232 {
 233   //
 234   // Initialize detectors
 235   //
 236
 237   AliInfo("Initialisation:");
 238   TStopwatch stw;
 239   TIter next(gAlice->Modules());
 240   AliModule *detector;
 241   while((detector = dynamic_cast<AliModule*>(next()))) {
 242     stw.Start();
 243     // Initialise detector geometry
 244     if(gAlice->IsRootGeometry()) detector->CreateMaterials();
 245     detector->Init();
 246     AliInfo(Form("%10s R:%.2fs C:%.2fs",
 247                  detector->GetName(),stw.RealTime(),stw.CpuTime()));
 248   }
 249 }
 250
 251 //_______________________________________________________________________
 252 void  AliMC::SetAllAlignableVolumes()
 253 {
 254   //
 255   // Add alignable volumes (TGeoPNEntries) looping on all
 256   // active modules
 257   //
 258
 259   AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
 260   AliModule *detector;
 261   TIter next(gAlice->Modules());
 262   while((detector = dynamic_cast<AliModule*>(next()))) {
 263     detector->AddAlignableVolumes();
 264   }
 265 }
 266
 267 //_______________________________________________________________________
 268 void  AliMC::GeneratePrimaries()
 269 {
 270   //
 271   // Generate primary particles and fill them in the stack.
 272   //
 273
 274   Generator()->Generate();
 275 }
 276
 277 //_______________________________________________________________________
 278 void AliMC::SetGenerator(AliGenerator *generator)
 279 {
 280   //
 281   // Load the event generator
 282   //
 283   if(!fGenerator) fGenerator = generator;
 284 }
 285
 286 //_______________________________________________________________________
 287 void AliMC::ResetGenerator(AliGenerator *generator)
 288 {
 289   //
 290   // Load the event generator
 291   //
 292   if(fGenerator) {
 293     if(generator) {
 294       AliWarning(Form("Replacing generator %s with %s",
 295                       fGenerator->GetName(),generator->GetName()));
 296     }
 297     else {
 298       AliWarning(Form("Replacing generator %s with NULL",
 299                       fGenerator->GetName()));
 300     }
 301   }
 302   fGenerator = generator;
 303 }
 304
 305 //_______________________________________________________________________
 306 void AliMC::FinishRun()
 307 {
 308   // Clean generator information
 309   AliDebug(1, "fGenerator->FinishRun()");
 310   fGenerator->FinishRun();
 311
 312   //Output energy summary tables
 313   AliDebug(1, "EnergySummary()");
 314   ToAliDebug(1, EnergySummary());
 315 }
 316
 317 //_______________________________________________________________________
 318 void AliMC::BeginPrimary()
 319 {
 320   //
 321   // Called  at the beginning of each primary track
 322   //
 323
 324   // Reset Hits info
 325   ResetHits();
 326   ResetTrackReferences();
 327
 328 }
 329
 330 //_______________________________________________________________________
 331 void AliMC::PreTrack()
 332 {
 333   // Actions before the track's transport
 334      TObjArray &dets = *gAlice->Modules();
 335      AliModule *module;
 336
 337      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 338        if((module = dynamic_cast<AliModule*>(dets[i])))
 339          module->PreTrack();
 340
 341      fMCQA->PreTrack();
 342 }
 343
 344 //_______________________________________________________________________
 345 void AliMC::Stepping()
 346 {
 347   //
 348   // Called at every step during transport
 349   //
 350
 351   Int_t id = DetFromMate(gMC->CurrentMedium());
 352   if (id < 0) return;
 353
 354
 355   if ( gMC->IsNewTrack()            &&
 356        gMC->TrackTime() == 0.       &&
 357        fRDecayMin >= 0.             &&
 358        fRDecayMax > fRDecayMin      &&
 359        gMC->TrackPid() == fDecayPdg )
 360   {
 361       FixParticleDecaytime();
 362   }
 363
 364
 365
 366   //
 367   // --- If lego option, do it and leave
 368   if (gAlice->Lego())
 369     gAlice->Lego()->StepManager();
 370   else {
 371     Int_t copy;
 372     //Update energy deposition tables
 373     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
 374     //
 375     // write tracke reference for track which is dissapearing - MI
 376     if (gMC->IsTrackDisappeared()) {
 377         if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber(),
 378                                                 AliTrackReference::kDisappeared);
 379     }
 380
 381     //Call the appropriate stepping routine;
 382     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
 383     if(det && det->StepManagerIsEnabled()) {
 384       if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id);
 385       det->StepManager();
 386     }
 387   }
 388 }
 389
 390 //_______________________________________________________________________
 391 void AliMC::EnergySummary()
 392 {
 393   //e
 394   // Print summary of deposited energy
 395   //
 396
 397   Int_t ndep=0;
 398   Float_t edtot=0;
 399   Float_t ed, ed2;
 400   Int_t kn, i, left, j, id;
 401   const Float_t kzero=0;
 402   Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
 403   //
 404   // Energy loss information
 405   if(ievent) {
 406     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 407     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 408       ed=fSummEnergy[kn];
 409       if(ed>0) {
 410         fEventEnergy[ndep]=kn;
 411         if(ievent>1) {
 412           ed=ed/ievent;
 413           ed2=fSum2Energy[kn];
 414           ed2=ed2/ievent;
 415           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
 416         } else
 417           ed2=99;
 418         fSummEnergy[ndep]=ed;
 419         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
 420         edtot+=ed;
 421         ndep++;
 422       }
 423     }
 424     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 425       left=ndep-kn*3;
 426       for(i=0;i<(3<left?3:left);i++) {
 427         j=kn*3+i;
 428         id=Int_t (fEventEnergy[j]+0.1);
 429         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 430       }
 431       printf("\n");
 432     }
 433     //
 434     // Relative energy loss in different detectors
 435     printf("******************** Relative Energy Loss per event ********************\n");
 436     printf("Total energy loss per event %10.3f GeV\n",edtot);
 437     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 438       left=ndep-kn*5;
 439       for(i=0;i<(5<left?5:left);i++) {
 440         j=kn*5+i;
 441         id=Int_t (fEventEnergy[j]+0.1);
 442         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 443       }
 444       printf("\n");
 445     }
 446     for(kn=0;kn<75;kn++) printf("*");
 447     printf("\n");
 448   }
 449   //
 450   // Reset the TArray's
 451   //  fEventEnergy.Set(0);
 452   //  fSummEnergy.Set(0);
 453   //  fSum2Energy.Set(0);
 454 }
 455
 456 //_____________________________________________________________________________
 457 void AliMC::BeginEvent()
 458 {
 459   //
 460   // Clean-up previous event
 461   // Energy scores
 462   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 463   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 464   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 465   AliDebug(1, "          BEGINNING EVENT               ");
 466   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 467   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 468   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 469
 470     AliRunLoader *runloader=gAlice->GetRunLoader();
 471
 472   /*******************************/
 473   /*   Clean after eventual      */
 474   /*   previous event            */
 475   /*******************************/
 476
 477
 478   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
 479   gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
 480   runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
 481   AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
 482
 483   fEventEnergy.Reset();
 484     // Clean detector information
 485
 486   if (runloader->Stack())
 487     runloader->Stack()->Reset();//clean stack -> tree is unloaded
 488   else
 489     runloader->MakeStack();//or make a new one
 490
 491 //  runloader->Stack()->BeginEvent();
 492
 493 if(gAlice->Lego() == 0x0)
 494    {
 495      AliDebug(1, "fRunLoader->MakeTree(K)");
 496      runloader->MakeTree("K");
 497    }
 498
 499   AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
 500   gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here
 501                                      //because we don't have guarantee that
 502                                      //stack pointer is not going to change from event to event
 503                          //since it bellobgs to header and is obtained via RunLoader
 504   //
 505   //  Reset all Detectors & kinematics & make/reset trees
 506   //
 507
 508   runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
 509     gAlice->GetEventNrInRun());
 510 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 511
 512   if(gAlice->Lego())
 513    {
 514     gAlice->Lego()->BeginEvent();
 515     return;
 516    }
 517
 518
 519   AliDebug(1, "ResetHits()");
 520   ResetHits();
 521
 522   AliDebug(1, "fRunLoader->MakeTree(H)");
 523   runloader->MakeTree("H");
 524
 525   AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
 526   runloader->MakeTrackRefsContainer();//for insurance
 527
 528
 529   //create new branches and SetAdresses
 530   TIter next(gAlice->Modules());
 531   AliModule *detector;
 532   while((detector = (AliModule*)next()))
 533    {
 534     AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
 535     detector->MakeBranch("H");
 536    }
 537   // make branch for AliRun track References
 538   TTree * treeTR = runloader->TreeTR();
 539   if (treeTR){
 540     // make branch for central track references
 541     if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
 542     TBranch *branch;
 543     branch = treeTR->Branch("TrackReferences",&fTrackReferences);
 544     branch->SetAddress(&fTrackReferences);
 545   }
 546   //
 547 }
 548
 549 //_______________________________________________________________________
 550 void AliMC::ResetHits()
 551 {
 552   //
 553   //  Reset all Detectors hits
 554   //
 555   TIter next(gAlice->Modules());
 556   AliModule *detector;
 557   while((detector = dynamic_cast<AliModule*>(next()))) {
 558      detector->ResetHits();
 559   }
 560 }
 561
 562 //_______________________________________________________________________
 563 void AliMC::PostTrack()
 564 {
 565   // Posts tracks for each module
 566   TObjArray &dets = *gAlice->Modules();
 567   AliModule *module;
 568
 569   for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 570     if((module = dynamic_cast<AliModule*>(dets[i])))
 571       module->PostTrack();
 572 }
 573
 574 //_______________________________________________________________________
 575 void AliMC::FinishPrimary()
 576 {
 577   //
 578   // Called  at the end of each primary track
 579   //
 580   AliRunLoader *runloader=gAlice->GetRunLoader();
 581   //  static Int_t count=0;
 582   //  const Int_t times=10;
 583   // This primary is finished, purify stack
 584 #if ROOT_VERSION_CODE > 262152
 585   if (!(gMC->SecondariesAreOrdered())) {
 586       runloader->Stack()->ReorderKine();
 587       RemapHits();
 588   }
 589 #endif
 590   runloader->Stack()->PurifyKine();
 591   RemapHits();
 592
 593   TIter next(gAlice->Modules());
 594   AliModule *detector;
 595   while((detector = dynamic_cast<AliModule*>(next()))) {
 596     detector->FinishPrimary();
 597     AliLoader* loader = detector->GetLoader();
 598     if(loader)
 599      {
 600        TTree* treeH = loader->TreeH();
 601        if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
 602      }
 603   }
 604
 605   // Write out track references if any
 606   if (runloader->TreeTR()) runloader->TreeTR()->Fill();
 607 }
 608
 609 void AliMC::RemapHits()
 610 {
 611 //
 612 // Remaps the track labels of the hits
 613     AliRunLoader *runloader=gAlice->GetRunLoader();
 614     AliStack* stack = runloader->Stack();
 615     TList* hitLists = GetHitLists();
 616     TIter next(hitLists);
 617     TCollection *hitList;
 618
 619     while((hitList = dynamic_cast<TCollection*>(next()))) {
 620         TIter nexthit(hitList);
 621         AliHit *hit;
 622         while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 623             hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
 624         }
 625     }
 626
 627     //
 628     // This for detectors which have a special mapping mechanism
 629     // for hits, such as TPC and TRD
 630     //
 631
 632     TObjArray* modules = gAlice->Modules();
 633     TIter nextmod(modules);
 634     AliDetector *detector;
 635     while((detector = dynamic_cast<AliDetector*>(nextmod()))) {
 636         detector->RemapTrackHitIDs(stack->TrackLabelMap());
 637     }
 638     //
 639     RemapTrackReferencesIDs(stack->TrackLabelMap());
 640 }
 641
 642 //_______________________________________________________________________
 643 void AliMC::FinishEvent()
 644 {
 645   //
 646   // Called at the end of the event.
 647   //
 648
 649   //
 650
 651   if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
 652
 653   TIter next(gAlice->Modules());
 654   AliModule *detector;
 655   while((detector = dynamic_cast<AliModule*>(next()))) {
 656     detector->FinishEvent();
 657   }
 658
 659   //Update the energy deposit tables
 660   Int_t i;
 661   for(i=0;i<fEventEnergy.GetSize();i++)
 662    {
 663     fSummEnergy[i]+=fEventEnergy[i];
 664     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
 665    }
 666
 667   AliRunLoader *runloader=gAlice->GetRunLoader();
 668
 669   AliHeader* header = runloader->GetHeader();
 670   AliStack* stack = runloader->Stack();
 671   if ( (header == 0x0) || (stack == 0x0) )
 672    {//check if we got header and stack. If not cry and exit aliroot
 673     AliFatal("Can not get the stack or header from LOADER");
 674     return;//never reached
 675    }
 676   // Update Header information
 677   header->SetNprimary(stack->GetNprimary());
 678   header->SetNtrack(stack->GetNtrack());
 679
 680
 681   // Write out the kinematics
 682   if (!gAlice->Lego()) stack->FinishEvent();
 683
 684   // Write out the event Header information
 685   TTree* treeE = runloader->TreeE();
 686   if (treeE)
 687    {
 688       header->SetStack(stack);
 689       treeE->Fill();
 690    }
 691   else
 692    {
 693     AliError("Can not get TreeE from RL");
 694    }
 695
 696   if(gAlice->Lego() == 0x0)
 697    {
 698      runloader->WriteKinematics("OVERWRITE");
 699      runloader->WriteTrackRefs("OVERWRITE");
 700      runloader->WriteHits("OVERWRITE");
 701    }
 702
 703   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 704   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 705   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 706   AliDebug(1, "          FINISHING EVENT               ");
 707   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 708   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 709   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 710 }
 711
 712 //_______________________________________________________________________
 713 void AliMC::Field(const Double_t* x, Double_t* b) const
 714 {
 715   // Calculates field "b" at point "x"
 716     gAlice->Field(x,b);
 717 }
 718
 719 //_______________________________________________________________________
 720 void AliMC::Init()
 721 {
 722   // MC initialization
 723
 724    //=================Create Materials and geometry
 725    gMC->Init();
 726   // Set alignable volumes for the whole geometry (with old root)
 727 #if ROOT_VERSION_CODE < 331527
 728   SetAllAlignableVolumes();
 729 #endif
 730    //Read the cuts for all materials
 731    ReadTransPar();
 732    //Build the special IMEDIA table
 733    MediaTable();
 734
 735    //Compute cross-sections
 736    gMC->BuildPhysics();
 737
 738    //Initialise geometry deposition table
 739    fEventEnergy.Set(gMC->NofVolumes()+1);
 740    fSummEnergy.Set(gMC->NofVolumes()+1);
 741    fSum2Energy.Set(gMC->NofVolumes()+1);
 742
 743    //
 744    fMCQA = new AliMCQA(gAlice->GetNdets());
 745
 746    // Register MC in configuration
 747    AliConfig::Instance()->Add(gMC);
 748
 749 }
 750
 751 //_______________________________________________________________________
 752 void AliMC::MediaTable()
 753 {
 754   //
 755   // Built media table to get from the media number to
 756   // the detector id
 757   //
 758
 759   Int_t kz, nz, idt, lz, i, k, ind;
 760   //  Int_t ibeg;
 761   TObjArray &dets = *gAlice->Detectors();
 762   AliModule *det;
 763   Int_t ndets=gAlice->GetNdets();
 764   //
 765   // For all detectors
 766   for (kz=0;kz<ndets;kz++) {
 767     // If detector is defined
 768     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
 769         TArrayI &idtmed = *(det->GetIdtmed());
 770         for(nz=0;nz<100;nz++) {
 771
 772         // Find max and min material number
 773         if((idt=idtmed[nz])) {
 774           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 775           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 776         }
 777       }
 778       if(det->LoMedium() > det->HiMedium()) {
 779         det->LoMedium() = 0;
 780         det->HiMedium() = 0;
 781       } else {
 782         if(det->HiMedium() > fImedia->GetSize()) {
 783           AliError(Form("Increase fImedia from %d to %d",
 784                         fImedia->GetSize(),det->HiMedium()));
 785           return;
 786         }
 787         // Tag all materials in rage as belonging to detector kz
 788         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 789           (*fImedia)[lz]=kz;
 790         }
 791       }
 792     }
 793   }
 794   //
 795   // Print summary table
 796   AliInfo("Tracking media ranges:");
 797   ToAliInfo(
 798   for(i=0;i<(ndets-1)/6+1;i++) {
 799     for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
 800       ind=i*6+k;
 801       det=dynamic_cast<AliModule*>(dets[ind]);
 802       if(det)
 803         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 804                det->HiMedium());
 805       else
 806         printf(" %6s: %3d -> %3d;","NULL",0,0);
 807     }
 808     printf("\n");
 809   }
 810   )
 811 }
 812
 813 //_______________________________________________________________________
 814 void AliMC::ReadTransPar()
 815 {
 816   //
 817   // Read filename to set the transport parameters
 818   //
 819
 820
 821   const Int_t kncuts=10;
 822   const Int_t knflags=11;
 823   const Int_t knpars=kncuts+knflags;
 824   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 825                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 826                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 827                                "MULS","PAIR","PHOT","RAYL"};
 828   char line[256];
 829   char detName[7];
 830   char* filtmp;
 831   Float_t cut[kncuts];
 832   Int_t flag[knflags];
 833   Int_t i, itmed, iret, ktmed, kz;
 834   FILE *lun;
 835   //
 836   // See whether the file is there
 837   filtmp=gSystem->ExpandPathName(fTransParName.Data());
 838   lun=fopen(filtmp,"r");
 839   delete [] filtmp;
 840   if(!lun) {
 841     AliWarning(Form("File %s does not exist!",fTransParName.Data()));
 842     return;
 843   }
 844   //
 845   while(1) {
 846     // Initialise cuts and flags
 847     for(i=0;i<kncuts;i++) cut[i]=-99;
 848     for(i=0;i<knflags;i++) flag[i]=-99;
 849     itmed=0;
 850     for(i=0;i<256;i++) line[i]='\0';
 851     // Read up to the end of line excluded
 852     iret=fscanf(lun,"%[^\n]",line);
 853     if(iret<0) {
 854       //End of file
 855       fclose(lun);
 856       return;
 857     }
 858     // Read the end of line
 859     fscanf(lun,"%*c");
 860     if(!iret) continue;
 861     if(line[0]=='*') continue;
 862     // Read the numbers
 863     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",
 864                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
 865                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
 866                 &flag[8],&flag[9],&flag[10]);
 867     if(!iret) continue;
 868     if(iret<0) {
 869       //reading error
 870       AliWarning(Form("Error reading file %s",fTransParName.Data()));
 871       continue;
 872     }
 873     // Check that the module exist
 874     AliModule *mod = gAlice->GetModule(detName);
 875     if(mod) {
 876       // Get the array of media numbers
 877       TArrayI &idtmed = *mod->GetIdtmed();
 878       // Check that the tracking medium code is valid
 879       if(0<=itmed && itmed < 100) {
 880         ktmed=idtmed[itmed];
 881         if(!ktmed) {
 882           AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
 883           continue;
 884         }
 885         // Set energy thresholds
 886         for(kz=0;kz<kncuts;kz++) {
 887           if(cut[kz]>=0) {
 888             AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
 889                              kpars[kz],cut[kz],itmed,mod->GetName()));
 890             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
 891           }
 892         }
 893         // Set transport mechanisms
 894         for(kz=0;kz<knflags;kz++) {
 895           if(flag[kz]>=0) {
 896             AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
 897                              kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
 898             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
 899           }
 900         }
 901       } else {
 902         AliWarning(Form("Invalid medium code %d",itmed));
 903         continue;
 904       }
 905     } else {
 906       AliDebug(1, Form("%s not present",detName));
 907       continue;
 908     }
 909   }
 910 }
 911
 912 //_______________________________________________________________________
 913 void AliMC::SetTransPar(const char *filename)
 914 {
 915   //
 916   // Sets the file name for transport parameters
 917   //
 918   fTransParName = filename;
 919 }
 920
 921 //_______________________________________________________________________
 922 void AliMC::Browse(TBrowser *b)
 923 {
 924   //
 925   // Called when the item "Run" is clicked on the left pane
 926   // of the Root browser.
 927   // It displays the Root Trees and all detectors.
 928   //
 929   //detectors are in folders anyway
 930   b->Add(fMCQA,"AliMCQA");
 931 }
 932
 933 //_______________________________________________________________________
 934 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 935 {
 936   //
 937   //  Add a hit to detector id
 938   //
 939   TObjArray &dets = *gAlice->Modules();
 940   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 941 }
 942
 943 //_______________________________________________________________________
 944 void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 945 {
 946   //
 947   // Add digit to detector id
 948   //
 949   TObjArray &dets = *gAlice->Modules();
 950   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 951 }
 952
 953 //_______________________________________________________________________
 954 Int_t AliMC::GetCurrentTrackNumber() const {
 955   //
 956   // Returns current track
 957   //
 958   return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
 959 }
 960
 961 //_______________________________________________________________________
 962 void AliMC::DumpPart (Int_t i) const
 963 {
 964   //
 965   // Dumps particle i in the stack
 966   //
 967   AliRunLoader * runloader = gAlice->GetRunLoader();
 968    if (runloader->Stack())
 969     runloader->Stack()->DumpPart(i);
 970 }
 971
 972 //_______________________________________________________________________
 973 void AliMC::DumpPStack () const
 974 {
 975   //
 976   // Dumps the particle stack
 977   //
 978   AliRunLoader * runloader = gAlice->GetRunLoader();
 979    if (runloader->Stack())
 980     runloader->Stack()->DumpPStack();
 981 }
 982
 983 //_______________________________________________________________________
 984 Int_t AliMC::GetNtrack() const {
 985   //
 986   // Returns number of tracks in stack
 987   //
 988   Int_t ntracks = -1;
 989   AliRunLoader * runloader = gAlice->GetRunLoader();
 990    if (runloader->Stack())
 991      ntracks = runloader->Stack()->GetNtrack();
 992    return ntracks;
 993 }
 994
 995 //_______________________________________________________________________
 996 Int_t AliMC::GetPrimary(Int_t track) const
 997 {
 998   //
 999   // return number of primary that has generated track
1000   //
1001   Int_t nprimary = -999;
1002   AliRunLoader * runloader = gAlice->GetRunLoader();
1003   if (runloader->Stack())
1004     nprimary = runloader->Stack()->GetPrimary(track);
1005   return nprimary;
1006 }
1007
1008 //_______________________________________________________________________
1009 TParticle* AliMC::Particle(Int_t i) const
1010 {
1011   // Returns the i-th particle from the stack taking into account
1012   // the remaping done by PurifyKine
1013   AliRunLoader * runloader = gAlice->GetRunLoader();
1014   if (runloader)
1015    if (runloader->Stack())
1016     return runloader->Stack()->Particle(i);
1017   return 0x0;
1018 }
1019
1020 //_______________________________________________________________________
1021 TObjArray* AliMC::Particles() const {
1022   //
1023   // Returns pointer to Particles array
1024   //
1025   AliRunLoader * runloader = gAlice->GetRunLoader();
1026   if (runloader)
1027    if (runloader->Stack())
1028     return runloader->Stack()->Particles();
1029   return 0x0;
1030 }
1031
1032 //_______________________________________________________________________
1033 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1034                       Float_t *vpos, Float_t *polar, Float_t tof,
1035                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1036 {
1037 // Delegate to stack
1038 //
1039   AliRunLoader * runloader = gAlice->GetRunLoader();
1040   if (runloader)
1041     if (runloader->Stack())
1042       runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1043                                     mech, ntr, weight, is);
1044 }
1045
1046 //_______________________________________________________________________
1047 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1048                       Double_t px, Double_t py, Double_t pz, Double_t e,
1049                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1050                       Double_t polx, Double_t poly, Double_t polz,
1051                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1052 {
1053   // Delegate to stack
1054   //
1055   AliRunLoader * runloader = gAlice->GetRunLoader();
1056   if (runloader)
1057     if (runloader->Stack())
1058       runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1059                                     polx, poly, polz, mech, ntr, weight, is);
1060 }
1061
1062 //_______________________________________________________________________
1063 void AliMC::SetHighWaterMark(Int_t nt) const
1064 {
1065     //
1066     // Set high water mark for last track in event
1067   AliRunLoader * runloader = gAlice->GetRunLoader();
1068   if (runloader)
1069     if (runloader->Stack())
1070       runloader->Stack()->SetHighWaterMark(nt);
1071 }
1072
1073 //_______________________________________________________________________
1074 void AliMC::KeepTrack(Int_t track) const
1075 {
1076   //
1077   // Delegate to stack
1078   //
1079   AliRunLoader * runloader = gAlice->GetRunLoader();
1080   if (runloader)
1081     if (runloader->Stack())
1082       runloader->Stack()->KeepTrack(track);
1083 }
1084
1085 //_______________________________________________________________________
1086 void AliMC::FlagTrack(Int_t track) const
1087 {
1088   // Delegate to stack
1089   //
1090   AliRunLoader * runloader = gAlice->GetRunLoader();
1091   if (runloader)
1092     if (runloader->Stack())
1093       runloader->Stack()->FlagTrack(track);
1094 }
1095
1096 //_______________________________________________________________________
1097 void AliMC::SetCurrentTrack(Int_t track) const
1098 {
1099   //
1100   // Set current track number
1101   //
1102   AliRunLoader * runloader = gAlice->GetRunLoader();
1103   if (runloader)
1104     if (runloader->Stack())
1105       runloader->Stack()->SetCurrentTrack(track);
1106 }
1107
1108 //_______________________________________________________________________
1109 AliTrackReference*  AliMC::AddTrackReference(Int_t label, Int_t id)
1110 {
1111   //
1112   // add a trackrefernce to the list
1113   if (!fTrackReferences) {
1114     AliError("Container trackrefernce not active");
1115     return NULL;
1116   }
1117   Int_t primary = GetPrimary(label);
1118   Particle(primary)->SetBit(kKeepBit);
1119
1120   Int_t nref = fTrackReferences->GetEntriesFast();
1121   TClonesArray &lref = *fTrackReferences;
1122   return new(lref[nref]) AliTrackReference(label, id);
1123 }
1124
1125
1126
1127 //_______________________________________________________________________
1128 void AliMC::ResetTrackReferences()
1129 {
1130   //
1131   //  Reset all  references
1132   //
1133   if (fTrackReferences)   fTrackReferences->Clear();
1134 }
1135
1136 void AliMC::RemapTrackReferencesIDs(Int_t *map)
1137 {
1138   //
1139   // Remapping track reference
1140   // Called at finish primary
1141   //
1142   if (!fTrackReferences) return;
1143   Int_t nEntries = fTrackReferences->GetEntries();
1144   for (Int_t i=0; i < nEntries; i++){
1145       AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
1146       if (ref) {
1147           Int_t newID = map[ref->GetTrack()];
1148           if (newID>=0) ref->SetTrack(newID);
1149           else {
1150               ref->SetBit(kNotDeleted,kFALSE);
1151               fTrackReferences->RemoveAt(i);
1152           }
1153       } // if ref
1154   }
1155   fTrackReferences->Compress();
1156 }
1157
1158 void AliMC::FixParticleDecaytime()
1159 {
1160     //
1161     // Fix the particle decay time according to rmin and rmax for decays
1162     //
1163
1164     TLorentzVector p;
1165     gMC->TrackMomentum(p);
1166     Double_t tmin, tmax;
1167     Double_t b;
1168
1169     // Transverse velocity
1170     Double_t vt    = p.Pt() / p.E();
1171
1172     if ((b = gAlice->Field()->SolenoidField()) > 0.) {     // [kG]
1173
1174         // Radius of helix
1175
1176         Double_t rho   = p.Pt() / 0.0003 / b; // [cm]
1177
1178         // Revolution frequency
1179
1180         Double_t omega = vt / rho;
1181
1182         // Maximum and minimum decay time
1183         //
1184         // Check for curlers first
1185         if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
1186
1187         //
1188
1189         tmax  = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega;   // [ct]
1190         tmin  = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega;   // [ct]
1191     } else {
1192         tmax =  fRDecayMax / vt;                                                      // [ct]
1193         tmin =  fRDecayMin / vt;                                                      // [ct]
1194     }
1195
1196     //
1197     // Dial t using the two limits
1198     Double_t t = tmin + (tmax - tmin) * gRandom->Rndm();                              // [ct]
1199     //
1200     //
1201     // Force decay time in transport code
1202     //
1203     gMC->ForceDecayTime(t / 2.99792458e10);
1204 }