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