Add fast gaussian fitting routine (Marian)

[u/mrichter/AliRoot.git] / STEER / AliMC.cxx

diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx
index 099237fa341ed32ecde29cabbf830ca7456cc308..8d9be17aca1355070030e307198f5b9e8907d365 100644 (file)
--- a/STEER/AliMC.cxx
+++ b/STEER/AliMC.cxx
@@ -23,12 +23,12 @@
 
 #include <RVersion.h>
 #include <TBrowser.h>
+#include <TClonesArray.h>
+#include <TGeoManager.h>
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
-#include <TGeoManager.h>
-#include "TGeant3.h" 
-
+ 
 #include "AliLog.h"
 #include "AliDetector.h"
 #include "AliGenerator.h"
@@ -53,7 +53,7 @@ AliMC::AliMC() :
   fTrRmax(1.e10),
   fTrZmax(1.e10),
   fRDecayMax(1.e10),
-  fRDecayMin(0),
+  fRDecayMin(-1.),
   fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
@@ -76,7 +76,7 @@ AliMC::AliMC(const char *name, const char *title) :
   fTrRmax(1.e10),
   fTrZmax(1.e10),
   fRDecayMax(1.e10),
-  fRDecayMin(0),
+  fRDecayMin(-1.),
   fDecayPdg(0),
   fImedia(new TArrayI(1000)),
   fTransParName("\0"),
@@ -102,7 +102,7 @@ AliMC::AliMC(const AliMC &mc) :
   fTrRmax(1.e10),
   fTrZmax(1.e10),
   fRDecayMax(1.e10),
-  fRDecayMin(0),
+  fRDecayMin(-1.),
   fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
@@ -151,7 +151,7 @@ void  AliMC::ConstructGeometry()
   //
 
   if(gAlice->IsRootGeometry()){
-    // Load geometry 
+    // Load geometry
     const char *geomfilename = gAlice->GetGeometryFileName();
     if(gSystem->ExpandPathName(geomfilename)){
       AliInfo(Form("Loading geometry from file:\n %40s\n\n",geomfilename));
@@ -160,9 +160,7 @@ void  AliMC::ConstructGeometry()
       AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
       return;
     }
-
   }else{
-
     // Create modules, materials, geometry
     TStopwatch stw;
     TIter next(gAlice->Modules());
@@ -176,32 +174,61 @@ void  AliMC::ConstructGeometry()
       AliInfo(Form("%10s R:%.2fs C:%.2fs",
                   detector->GetName(),stw.RealTime(),stw.CpuTime()));
     }
-
   }
   
 }
 
+//_______________________________________________________________________
+void  AliMC::ConstructOpGeometry() 
+{
+  //
+  // Loop all detector modules and call DefineOpticalProperties() method 
+  //
+
+  TIter next(gAlice->Modules());
+  AliModule *detector;
+  AliInfo("Optical properties definition");
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    // Initialise detector optical properties
+    detector->DefineOpticalProperties();
+  }  
+}
+
 //_______________________________________________________________________
 void  AliMC::InitGeometry()
 { 
   //
-  // Initialize detectors and display geometry
+  // Initialize detectors
   //
 
-   AliInfo("Initialisation:");
-    TStopwatch stw;
-    TIter next(gAlice->Modules());
-    AliModule *detector;
-    while((detector = dynamic_cast<AliModule*>(next()))) {
-      stw.Start();
-      // Initialise detector and display geometry
-      if (gAlice->IsRootGeometry()) detector->CreateMaterials();
-      detector->Init();
-      detector->BuildGeometry();
-      AliInfo(Form("%10s R:%.2fs C:%.2fs",
-                  detector->GetName(),stw.RealTime(),stw.CpuTime()));
-    }
- 
+  AliInfo("Initialisation:");
+  TStopwatch stw;
+  TIter next(gAlice->Modules());
+  AliModule *detector;
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    stw.Start();
+    // Initialise detector geometry
+    if(gAlice->IsRootGeometry()) detector->CreateMaterials();
+    detector->Init();
+    AliInfo(Form("%10s R:%.2fs C:%.2fs",
+                detector->GetName(),stw.RealTime(),stw.CpuTime()));
+  }
+}
+
+//_______________________________________________________________________
+void  AliMC::SetAllAlignableVolumes()
+{ 
+  //
+  // Add alignable volumes (TGeoPNEntries) looping on all
+  // active modules
+  //
+
+  AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
+  AliModule *detector;
+  TIter next(gAlice->Modules());
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    detector->AddAlignableVolumes();
+  }
 }
 
 //_______________________________________________________________________
@@ -229,13 +256,16 @@ void AliMC::ResetGenerator(AliGenerator *generator)
   //
   // Load the event generator
   //
-  if(fGenerator)
-    if(generator)
+  if(fGenerator) {
+    if(generator) {
       AliWarning(Form("Replacing generator %s with %s",
-                     fGenerator->GetName(),generator->GetName()))
-    else
+                     fGenerator->GetName(),generator->GetName()));
+    }
+    else {
       AliWarning(Form("Replacing generator %s with NULL",
                      fGenerator->GetName()));
+    }
+  }
   fGenerator = generator;
 }
 
@@ -285,13 +315,13 @@ void AliMC::Stepping()
   // Called at every step during transport
   //
     
-  Int_t id = DetFromMate(gMC->GetMedium());
+  Int_t id = DetFromMate(gMC->CurrentMedium());
   if (id < 0) return;
 
 
   if ( gMC->IsNewTrack()            && 
        gMC->TrackTime() == 0.       &&
-       fRDecayMin > 0.              &&  
+       fRDecayMin >= 0.             &&  
        fRDecayMax > fRDecayMin      &&
        gMC->TrackPid() == fDecayPdg ) 
   {
@@ -317,7 +347,7 @@ void AliMC::Stepping()
     //Call the appropriate stepping routine;
     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
     if(det && det->StepManagerIsEnabled()) {
-      fMCQA->StepManager(id);
+      if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id);
       det->StepManager();
     }
   }
@@ -549,6 +579,7 @@ void AliMC::FinishEvent()
   //
   
   //
+    
   if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
 
   TIter next(gAlice->Modules());
@@ -580,7 +611,7 @@ void AliMC::FinishEvent()
 
   
   // Write out the kinematics
-  stack->FinishEvent();
+  if (!gAlice->Lego()) stack->FinishEvent();
    
   // Write out the event Header information
   TTree* treeE = runloader->TreeE();
@@ -624,6 +655,8 @@ void AliMC::Init()
 
    //=================Create Materials and geometry
    gMC->Init();
+   //Set alignable volumes for the whole geometry
+   SetAllAlignableVolumes();
    //Read the cuts for all materials
    ReadTransPar();
    //Build the special IMEDIA table
@@ -632,9 +665,6 @@ void AliMC::Init()
    //Compute cross-sections
    gMC->BuildPhysics();
    
-   //Write Geometry object to current file.
-   gAlice->GetRunLoader()->WriteGeometry();
-
    //Initialise geometry deposition table
    fEventEnergy.Set(gMC->NofVolumes()+1);
    fSummEnergy.Set(gMC->NofVolumes()+1);
@@ -830,7 +860,6 @@ void AliMC::Browse(TBrowser *b)
   b->Add(fMCQA,"AliMCQA");
 }
 
-//PH
 //_______________________________________________________________________
 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 {
@@ -1043,17 +1072,17 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map)
   // Called at finish primary
   //
   if (!fTrackReferences) return;
-  for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
-    AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
-    if (ref) {
-      Int_t newID = map[ref->GetTrack()];
-      if (newID>=0) ref->SetTrack(newID);
-      else {
-        //ref->SetTrack(-1);
-        ref->SetBit(kNotDeleted,kFALSE);
-        fTrackReferences->RemoveAt(i);  
-      }      
-    }
+  Int_t nEntries = fTrackReferences->GetEntries();
+  for (Int_t i=0; i < nEntries; i++){
+      AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
+      if (ref) {
+         Int_t newID = map[ref->GetTrack()];
+         if (newID>=0) ref->SetTrack(newID);
+         else {
+             ref->SetBit(kNotDeleted,kFALSE);
+             fTrackReferences->RemoveAt(i);  
+         }      
+      } // if ref
   }
   fTrackReferences->Compress();
 }
@@ -1103,6 +1132,5 @@ void AliMC::FixParticleDecaytime()
     //
     // Force decay time in transport code
     //
-    TGeant3 * geant = (TGeant3*) gMC;
-    geant->Gcphys()->sumlif = t / p.Beta() / p.Gamma();
+    gMC->ForceDecayTime(t / 2.99792458e10);
 }