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