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