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