// Author: F.Carminati
// Federico.Carminati@cern.ch
-#include <TBrowser.h>
+#include <RVersion.h>
+#include <TArrayI.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TGeoGlobalMagField.h>
+#include <TGeoManager.h>
+#include <TParticle.h>
+#include <TROOT.h>
#include <TStopwatch.h>
#include <TSystem.h>
#include <TVirtualMC.h>
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
#include "AliDetector.h"
#include "AliGenerator.h"
+#include "AliGeomManager.h"
#include "AliHeader.h"
+#include "AliHit.h"
#include "AliLego.h"
+#include "AliLog.h"
#include "AliMC.h"
-#include "AliMCQA.h"
+#include "AliMagF.h"
#include "AliRun.h"
+#include "AliSimulation.h"
#include "AliStack.h"
#include "AliTrackReference.h"
-
ClassImp(AliMC)
//_______________________________________________________________________
fSum2Energy(0),
fTrRmax(1.e10),
fTrZmax(1.e10),
+ fRDecayMax(1.e10),
+ fRDecayMin(-1.),
+ fDecayPdg(0),
fImedia(0),
fTransParName("\0"),
- fMCQA(0),
fHitLists(0),
- fTrackReferences(0)
+ fTmpTreeTR(0),
+ fTmpFileTR(0),
+ fTrackReferences(),
+ fTmpTrackReferences()
{
//default constructor
+ DecayLimits();
}
//_______________________________________________________________________
fSum2Energy(0),
fTrRmax(1.e10),
fTrZmax(1.e10),
+ fRDecayMax(1.e10),
+ fRDecayMin(-1.),
+ fDecayPdg(0),
fImedia(new TArrayI(1000)),
fTransParName("\0"),
- fMCQA(0),
fHitLists(new TList()),
- fTrackReferences(new TClonesArray("AliTrackReference", 100))
+ fTmpTreeTR(0),
+ fTmpFileTR(0),
+ fTrackReferences("AliTrackReference", 100),
+ fTmpTrackReferences("AliTrackReference", 100)
{
//constructor
// Set transport parameters
SetTransPar();
-
+ DecayLimits();
// Prepare the tracking medium lists
for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
-
-}
-
-//_______________________________________________________________________
-AliMC::AliMC(const AliMC &mc) :
- TVirtualMCApplication(mc),
- fGenerator(0),
- fEventEnergy(0),
- fSummEnergy(0),
- fSum2Energy(0),
- fTrRmax(1.e10),
- fTrZmax(1.e10),
- fImedia(0),
- fTransParName("\0"),
- fMCQA(0),
- fHitLists(0),
- fTrackReferences(0)
-{
- //
- // Copy constructor for AliRun
- //
- mc.Copy(*this);
}
//_______________________________________________________________________
//destructor
delete fGenerator;
delete fImedia;
- delete fMCQA;
delete fHitLists;
// Delete track references
- if (fTrackReferences) {
- fTrackReferences->Delete();
- delete fTrackReferences;
- fTrackReferences = 0;
- }
-
-}
-
-//_______________________________________________________________________
-void AliMC::Copy(AliMC &) const
-{
- //dummy Copy function
- Fatal("Copy","Not implemented!\n");
}
//_______________________________________________________________________
void AliMC::ConstructGeometry()
{
//
- // Create modules, materials, geometry
- //
-
+ // Either load geometry from file or create it through usual
+ // loop on detectors. In the first case the method
+ // AliModule::CreateMaterials() only builds fIdtmed and is postponed
+ // at InitGeometry().
+ //
+
+ if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file
+ if(gAlice->IsGeomFromCDB()){
+ AliInfo("Loading geometry from CDB default storage");
+ AliCDBPath path("GRP","Geometry","Data");
+ AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
+ if(!entry) AliFatal("Unable to load geometry from CDB!");
+ entry->SetOwner(0);
+ gGeoManager = (TGeoManager*) entry->GetObject();
+ if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
+ }else{
+ // Load geometry
+ const char *geomfilename = gAlice->GetGeometryFileName();
+ if(gSystem->ExpandPathName(geomfilename)){
+ AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
+ TGeoManager::Import(geomfilename);
+ }else{
+ AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
+ return;
+ }
+ }
+ }else{
+ // Create modules, materials, geometry
+ if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
TStopwatch stw;
TIter next(gAlice->Modules());
AliModule *detector;
- if (GetDebug()) Info("ConstructGeometry","Geometry creation:");
+ AliDebug(1, "Geometry creation:");
while((detector = dynamic_cast<AliModule*>(next()))) {
stw.Start();
// Initialise detector materials and geometry
detector->CreateMaterials();
detector->CreateGeometry();
- printf("%10s R:%.2fs C:%.2fs\n",
- detector->GetName(),stw.RealTime(),stw.CpuTime());
+ AliInfo(Form("%10s R:%.2fs C:%.2fs",
+ detector->GetName(),stw.RealTime(),stw.CpuTime()));
}
+ }
+
+}
+
+//_______________________________________________________________________
+Bool_t AliMC::MisalignGeometry()
+{
+// Call misalignment code if AliSimulation object was defined.
+
+ if(!gAlice->IsRootGeometry()){
+ //Set alignable volumes for the whole geometry
+ SetAllAlignableVolumes();
+ }
+ // Misalign geometry via AliSimulation instance
+ if (!AliSimulation::Instance()) return kFALSE;
+ AliGeomManager::SetGeometry(gGeoManager);
+ if(!AliGeomManager::CheckSymNamesLUT("ALL"))
+ AliFatal("Current loaded geometry differs in the definition of symbolic names!");
+
+ return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance());
+}
+
+//_______________________________________________________________________
+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
//
- printf("Initialisation:\n");
- TStopwatch stw;
- TIter next(gAlice->Modules());
- AliModule *detector;
- while((detector = dynamic_cast<AliModule*>(next()))) {
- stw.Start();
- // Initialise detector and display geometry
- detector->Init();
- detector->BuildGeometry();
- printf("%10s R:%.2fs C:%.2fs\n",
- 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::GeneratePrimaries()
+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();
+ }
+}
+
+//_______________________________________________________________________
+void AliMC::GeneratePrimaries()
{
//
// Generate primary particles and fill them in the stack.
//
// Load the event generator
//
- if(fGenerator)
- if(generator)
- Warning("ResetGenerator","Replacing generator %s with %s\n",
- fGenerator->GetName(),generator->GetName());
- else
- Warning("ResetGenerator","Replacing generator %s with NULL\n",
- fGenerator->GetName());
+ if(fGenerator) {
+ if(generator) {
+ AliWarning(Form("Replacing generator %s with %s",
+ fGenerator->GetName(),generator->GetName()));
+ }
+ else {
+ AliWarning(Form("Replacing generator %s with NULL",
+ fGenerator->GetName()));
+ }
+ }
fGenerator = generator;
}
void AliMC::FinishRun()
{
// Clean generator information
- if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()");
+ AliDebug(1, "fGenerator->FinishRun()");
fGenerator->FinishRun();
//Output energy summary tables
- if (GetDebug()) Info("FinishRun"," EnergySummary()");
- EnergySummary();
-
+ AliDebug(1, "EnergySummary()");
+ ToAliDebug(1, EnergySummary());
}
//_______________________________________________________________________
// Reset Hits info
ResetHits();
ResetTrackReferences();
-
}
//_______________________________________________________________________
if((module = dynamic_cast<AliModule*>(dets[i])))
module->PreTrack();
- fMCQA->PreTrack();
}
//_______________________________________________________________________
//
// 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. &&
+ fRDecayMax > fRDecayMin &&
+ gMC->TrackPid() == fDecayPdg )
+ {
+ FixParticleDecaytime();
+ }
+
+
+
//
// --- If lego option, do it and leave
- if (gAlice->Lego())
- gAlice->Lego()->StepManager();
+ if (AliSimulation::Instance()->Lego())
+ AliSimulation::Instance()->Lego()->StepManager();
else {
Int_t copy;
//Update energy deposition tables
AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
-
+ //
+ // write tracke reference for track which is dissapearing - MI
+
+ if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) {
+ if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(),
+ AliTrackReference::kDisappeared);
+
+
+ }
+
//Call the appropriate stepping routine;
AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
if(det && det->StepManagerIsEnabled()) {
- fMCQA->StepManager(id);
det->StepManager();
}
}
//_______________________________________________________________________
void AliMC::EnergySummary()
{
- //
+ //e
// Print summary of deposited energy
//
Float_t ed, ed2;
Int_t kn, i, left, j, id;
const Float_t kzero=0;
- Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
+ Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1;
//
// Energy loss information
if(ievent) {
//
// Clean-up previous event
// Energy scores
- if (GetDebug())
- {
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- Info("BeginEvent"," BEGINNING EVENT ");
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- }
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
+ AliDebug(1, " BEGINNING EVENT ");
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
+ AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
- AliRunLoader *runloader=gAlice->GetRunLoader();
+ AliRunLoader *runloader=AliRunLoader::Instance();
/*******************************/
/* Clean after eventual */
//Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
- if (GetDebug()) Info("BeginEvent","EventNr is %d",gAlice->GetEventNrInRun());
+ AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
fEventEnergy.Reset();
// Clean detector information
if (runloader->Stack())
- runloader->Stack()->Reset();//clean stack -> tree is unloaded
+ runloader->Stack()->Reset();//clean stack -> tree is unloaded
else
- runloader->MakeStack();//or make a new one
+ runloader->MakeStack();//or make a new one
- if(gAlice->Lego() == 0x0)
- {
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(K)");
- runloader->MakeTree("K");
- }
-
- if (GetDebug()) Info("BeginEvent"," gMC->SetStack(fRunLoader->Stack())");
- gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here
+
+ if(AliSimulation::Instance()->Lego() == 0x0)
+ {
+ AliDebug(1, "fRunLoader->MakeTree(K)");
+ runloader->MakeTree("K");
+ }
+
+ AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
+ gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here
//because we don't have guarantee that
//stack pointer is not going to change from event to event
//since it bellobgs to header and is obtained via RunLoader
// Reset all Detectors & kinematics & make/reset trees
//
- runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
- gAlice->GetEventNrInRun());
+ runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(),
+ gAlice->GetEventNrInRun());
// fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so
- if(gAlice->Lego())
- {
- gAlice->Lego()->BeginEvent();
- return;
- }
+ if(AliSimulation::Instance()->Lego())
+ {
+ AliSimulation::Instance()->Lego()->BeginEvent();
+ return;
+ }
- if (GetDebug()) Info("BeginEvent"," ResetHits()");
+ AliDebug(1, "ResetHits()");
ResetHits();
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(H)");
+ AliDebug(1, "fRunLoader->MakeTree(H)");
runloader->MakeTree("H");
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTrackRefsContainer()");
- runloader->MakeTrackRefsContainer();//for insurance
+ MakeTmpTrackRefsTree();
//create new branches and SetAdresses
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = (AliModule*)next()))
{
- if (GetDebug()) Info("BeginEvent"," %s->MakeBranch(H)",detector->GetName());
- detector->MakeBranch("H");
- if (GetDebug()) Info("BeginEvent"," %s->MakeBranchTR()",detector->GetName());
- detector->MakeBranchTR();
- if (GetDebug()) Info("BeginEvent"," %s->SetTreeAddress()",detector->GetName());
- detector->SetTreeAddress();
+ AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
+ detector->MakeBranch("H");
}
}
}
}
+//_______________________________________________________________________
+void AliMC::ResetDigits()
+{
+ //
+ // Reset all Detectors digits
+ //
+ TIter next(gAlice->Modules());
+ AliModule *detector;
+ while((detector = dynamic_cast<AliModule*>(next()))) {
+ detector->ResetDigits();
+ }
+}
+
+//_______________________________________________________________________
+void AliMC::ResetSDigits()
+{
+ //
+ // Reset all Detectors digits
+ //
+ TIter next(gAlice->Modules());
+ AliModule *detector;
+ while((detector = dynamic_cast<AliModule*>(next()))) {
+ detector->ResetSDigits();
+ }
+}
+
//_______________________________________________________________________
void AliMC::PostTrack()
{
// Posts tracks for each module
- TObjArray &dets = *gAlice->Modules();
- AliModule *module;
-
- for(Int_t i=0; i<=gAlice->GetNdets(); i++)
-
if((module = dynamic_cast<AliModule*>(dets[i])))
- module->PostTrack();
+ TObjArray &dets = *gAlice->Modules();
+ AliModule *module;
+
+ for(Int_t i=0; i<=gAlice->GetNdets(); i++)
+ if((module = dynamic_cast<AliModule*>(dets[i])))
+ module->PostTrack();
}
//_______________________________________________________________________
//
// Called at the end of each primary track
//
- AliRunLoader *runloader=gAlice->GetRunLoader();
+ AliRunLoader *runloader=AliRunLoader::Instance();
// static Int_t count=0;
// const Int_t times=10;
// This primary is finished, purify stack
- runloader->Stack()->PurifyKine();
+#if ROOT_VERSION_CODE > 262152
+ if (!(gMC->SecondariesAreOrdered())) {
+ if (runloader->Stack()->ReorderKine()) RemapHits();
+ }
+#endif
+ if (runloader->Stack()->PurifyKine()) RemapHits();
TIter next(gAlice->Modules());
AliModule *detector;
}
// Write out track references if any
- if (runloader->TreeTR()) runloader->TreeTR()->Fill();
+ if (fTmpTreeTR) fTmpTreeTR->Fill();
+}
+
+void AliMC::RemapHits()
+{
+//
+// Remaps the track labels of the hits
+ AliRunLoader *runloader=AliRunLoader::Instance();
+ AliStack* stack = runloader->Stack();
+ TList* hitLists = GetHitLists();
+ TIter next(hitLists);
+ TCollection *hitList;
+
+ while((hitList = dynamic_cast<TCollection*>(next()))) {
+ TIter nexthit(hitList);
+ AliHit *hit;
+ while((hit = dynamic_cast<AliHit*>(nexthit()))) {
+ hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
+ }
+ }
+
+ //
+ // This for detectors which have a special mapping mechanism
+ // for hits, such as TPC and TRD
+ //
+
+
+ TObjArray* modules = gAlice->Modules();
+ TIter nextmod(modules);
+ AliModule *module;
+ while((module = (AliModule*) nextmod())) {
+ AliDetector* det = dynamic_cast<AliDetector*> (module);
+ if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
+ }
+ //
+ RemapTrackReferencesIDs(stack->TrackLabelMap());
}
//_______________________________________________________________________
//
//
- if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
+
+ if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent();
TIter next(gAlice->Modules());
AliModule *detector;
fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
}
- AliRunLoader *runloader=gAlice->GetRunLoader();
+ AliRunLoader *runloader=AliRunLoader::Instance();
AliHeader* header = runloader->GetHeader();
AliStack* stack = runloader->Stack();
if ( (header == 0x0) || (stack == 0x0) )
{//check if we got header and stack. If not cry and exit aliroot
- Fatal("AliRun","Can not get the stack or header from LOADER");
+ AliFatal("Can not get the stack or header from LOADER");
return;//never reached
}
// Update Header information
header->SetNprimary(stack->GetNprimary());
header->SetNtrack(stack->GetNtrack());
-
// Write out the kinematics
- stack->FinishEvent();
+ if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
+
+ // Synchronize the TreeTR with TreeK
+ if (fTmpTreeTR) ReorderAndExpandTreeTR();
// Write out the event Header information
TTree* treeE = runloader->TreeE();
}
else
{
- Error("FinishEvent","Can not get TreeE from RL");
+ AliError("Can not get TreeE from RL");
}
- if(gAlice->Lego() == 0x0)
+ if(AliSimulation::Instance()->Lego() == 0x0)
{
runloader->WriteKinematics("OVERWRITE");
runloader->WriteTrackRefs("OVERWRITE");
runloader->WriteHits("OVERWRITE");
}
- if (GetDebug())
- {
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- Info("FinishEvent"," FINISHING EVENT ");
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
- }
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, " FINISHING EVENT ");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
}
-//_______________________________________________________________________
-void AliMC::Field(const Double_t* x, Double_t* b) const
-{
- // Calculates field "b" at point "x"
- gAlice->Field(x,b);
-}
-
//_______________________________________________________________________
void AliMC::Init()
{
//=================Create Materials and geometry
gMC->Init();
- gMC->DefineParticles(); //Create standard MC particles
-
+ // Set alignable volumes for the whole geometry (with old root)
+#if ROOT_VERSION_CODE < 331527
+ SetAllAlignableVolumes();
+#endif
//Read the cuts for all materials
ReadTransPar();
//Build the special IMEDIA table
//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);
fSum2Energy.Set(gMC->NofVolumes()+1);
- //
- fMCQA = new AliMCQA(gAlice->GetNdets());
-
// Register MC in configuration
AliConfig::Instance()->Add(gMC);
if((det=dynamic_cast<AliModule*>(dets[kz]))) {
TArrayI &idtmed = *(det->GetIdtmed());
for(nz=0;nz<100;nz++) {
+
// Find max and min material number
if((idt=idtmed[nz])) {
det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
det->HiMedium() = 0;
} else {
if(det->HiMedium() > fImedia->GetSize()) {
- Error("MediaTable","Increase fImedia from %d to %d",
- fImedia->GetSize(),det->HiMedium());
+ AliError(Form("Increase fImedia from %d to %d",
+ fImedia->GetSize(),det->HiMedium()));
return;
}
// Tag all materials in rage as belonging to detector kz
}
//
// Print summary table
- printf(" Traking media ranges:\n");
+ AliInfo("Tracking media ranges:");
+ ToAliInfo(
for(i=0;i<(ndets-1)/6+1;i++) {
for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
ind=i*6+k;
}
printf("\n");
}
+ )
}
//_______________________________________________________________________
lun=fopen(filtmp,"r");
delete [] filtmp;
if(!lun) {
- Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
+ AliWarning(Form("File %s does not exist!",fTransParName.Data()));
return;
}
//
- if(fDebug) {
- printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
- printf(" *%59s\n","*");
- printf(" * Please check carefully what you are doing!%10s\n","*");
- printf(" *%59s\n","*");
- }
- //
while(1) {
// Initialise cuts and flags
for(i=0;i<kncuts;i++) cut[i]=-99;
if(iret<0) {
//End of file
fclose(lun);
- if(fDebug){
- printf(" *%59s\n","*");
- printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
- }
return;
}
// Read the end of line
if(!iret) continue;
if(iret<0) {
//reading error
- Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
+ AliWarning(Form("Error reading file %s",fTransParName.Data()));
continue;
}
// Check that the module exist
if(0<=itmed && itmed < 100) {
ktmed=idtmed[itmed];
if(!ktmed) {
- Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
+ AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
continue;
}
// Set energy thresholds
for(kz=0;kz<kncuts;kz++) {
if(cut[kz]>=0) {
- if(fDebug) printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
- kpars[kz],cut[kz],itmed,mod->GetName());
+ AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
+ kpars[kz],cut[kz],itmed,mod->GetName()));
gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
}
}
// Set transport mechanisms
for(kz=0;kz<knflags;kz++) {
if(flag[kz]>=0) {
- if(fDebug) printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
- kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
+ AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
+ kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
}
}
} else {
- Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
+ AliWarning(Form("Invalid medium code %d",itmed));
continue;
}
} else {
- if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
+ AliDebug(1, Form("%s not present",detName));
continue;
}
}
fTransParName = filename;
}
-//_______________________________________________________________________
-void AliMC::Browse(TBrowser *b)
-{
- //
- // Called when the item "Run" is clicked on the left pane
- // of the Root browser.
- // It displays the Root Trees and all detectors.
- //
- //detectors are in folders anyway
- b->Add(fMCQA,"AliMCQA");
-}
-
-//PH
//_______________________________________________________________________
void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
{
//
// Returns current track
//
- return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
+ return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
}
//_______________________________________________________________________
//
// Dumps particle i in the stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
runloader->Stack()->DumpPart(i);
}
//
// Dumps the particle stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
runloader->Stack()->DumpPStack();
}
// Returns number of tracks in stack
//
Int_t ntracks = -1;
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
ntracks = runloader->Stack()->GetNtrack();
return ntracks;
// return number of primary that has generated track
//
Int_t nprimary = -999;
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
nprimary = runloader->Stack()->GetPrimary(track);
return nprimary;
{
// Returns the i-th particle from the stack taking into account
// the remaping done by PurifyKine
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
return runloader->Stack()->Particle(i);
}
//_______________________________________________________________________
-TObjArray* AliMC::Particles() const {
+const TObjArray* AliMC::Particles() const {
//
// Returns pointer to Particles array
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
return runloader->Stack()->Particles();
}
//_______________________________________________________________________
-void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
- Float_t *vpos, Float_t *polar, Float_t tof,
- TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
+void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
+ const Float_t *vpos, const Float_t *polar, Float_t tof,
+ TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
{
// Delegate to stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
Double_t px, Double_t py, Double_t pz, Double_t e,
Double_t vx, Double_t vy, Double_t vz, Double_t tof,
Double_t polx, Double_t poly, Double_t polz,
- TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
+ TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
{
// Delegate to stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
}
//_______________________________________________________________________
-void AliMC::SetHighWaterMark(const Int_t nt)
+void AliMC::SetHighWaterMark(Int_t nt) const
{
//
// Set high water mark for last track in event
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->SetHighWaterMark(nt);
}
//_______________________________________________________________________
-void AliMC::KeepTrack(const Int_t track)
+void AliMC::KeepTrack(Int_t track) const
{
//
// Delegate to stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->KeepTrack(track);
}
//_______________________________________________________________________
-void AliMC::FlagTrack(Int_t track)
+void AliMC::FlagTrack(Int_t track) const
{
// Delegate to stack
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->FlagTrack(track);
}
//_______________________________________________________________________
-void AliMC::SetCurrentTrack(Int_t track)
+void AliMC::SetCurrentTrack(Int_t track) const
{
//
// Set current track number
//
- AliRunLoader * runloader = gAlice->GetRunLoader();
+ AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->SetCurrentTrack(track);
}
//_______________________________________________________________________
-void AliMC::AddTrackReference(Int_t label){
+AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id)
+{
//
// add a trackrefernce to the list
- if (!fTrackReferences) {
- cerr<<"Container trackrefernce not active\n";
- return;
- }
- Int_t nref = fTrackReferences->GetEntriesFast();
- TClonesArray &lref = *fTrackReferences;
- new(lref[nref]) AliTrackReference(label);
+ Int_t primary = GetPrimary(label);
+ Particle(primary)->SetBit(kKeepBit);
+
+ Int_t nref = fTmpTrackReferences.GetEntriesFast();
+ return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
}
//
// Reset all references
//
- if (fTrackReferences) fTrackReferences->Clear();
-
- TIter next(gAlice->Modules());
- AliModule *detector;
- while((detector = dynamic_cast<AliModule*>(next()))) {
- detector->ResetTrackReferences();
- }
+ fTmpTrackReferences.Clear();
}
void AliMC::RemapTrackReferencesIDs(Int_t *map)
// Remapping track reference
// 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 = fTmpTrackReferences.GetEntries();
+ for (Int_t i=0; i < nEntries; i++){
+ AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences.UncheckedAt(i));
+ if (ref) {
+ Int_t newID = map[ref->GetTrack()];
+ if (newID>=0) ref->SetTrack(newID);
+ else {
+ ref->SetBit(kNotDeleted,kFALSE);
+ fTmpTrackReferences.RemoveAt(i);
+ }
+ } // if ref
}
- fTrackReferences->Compress();
+ fTmpTrackReferences.Compress();
+}
+
+void AliMC::FixParticleDecaytime()
+{
+ //
+ // Fix the particle decay time according to rmin and rmax for decays
+ //
+
+ TLorentzVector p;
+ gMC->TrackMomentum(p);
+ Double_t tmin, tmax;
+ Double_t b;
+
+ // Transverse velocity
+ Double_t vt = p.Pt() / p.E();
+
+ if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) { // [kG]
+
+ // Radius of helix
+
+ Double_t rho = p.Pt() / 0.0003 / b; // [cm]
+
+ // Revolution frequency
+
+ Double_t omega = vt / rho;
+
+ // Maximum and minimum decay time
+ //
+ // Check for curlers first
+ if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
+
+ //
+
+ tmax = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega; // [ct]
+ tmin = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega; // [ct]
+ } else {
+ tmax = fRDecayMax / vt; // [ct]
+ tmin = fRDecayMin / vt; // [ct]
+ }
+
+ //
+ // Dial t using the two limits
+ Double_t t = tmin + (tmax - tmin) * gRandom->Rndm(); // [ct]
+ //
+ //
+ // Force decay time in transport code
+ //
+ gMC->ForceDecayTime(t / 2.99792458e10);
+}
+
+void AliMC::MakeTmpTrackRefsTree()
+{
+ // Make the temporary track reference tree
+ fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
+ fTmpTreeTR = new TTree("TreeTR", "Track References");
+ TClonesArray* pRef = &fTmpTrackReferences;
+ fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &pRef, 4000);
}
+
+void AliMC::ReorderAndExpandTreeTR()
+{
+//
+// Reorder and expand the temporary track reference tree in order to match the kinematics tree
+//
+
+ AliRunLoader *rl = AliRunLoader::Instance();
+//
+// TreeTR
+ AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
+ rl->MakeTrackRefsContainer();
+ TTree * treeTR = rl->TreeTR();
+ if (treeTR){
+ // make branch for central track references
+ TBranch *branch;
+ TClonesArray* pRef = &fTrackReferences;
+ branch = treeTR->Branch("TrackReferences", &pRef);
+ branch->SetAddress(&pRef);
+ }
+
+ AliStack* stack = rl->Stack();
+ Int_t np = stack->GetNprimary();
+ Int_t nt = fTmpTreeTR->GetEntries();
+ //
+ // Loop over tracks and find the secondaries with the help of the kine tree
+ Int_t ifills = 0;
+ Int_t it = 0;
+ for (Int_t ip = np - 1; ip > -1; ip--) {
+ TParticle *part = stack->Particle(ip);
+ //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
+
+ // Skip primaries that have not been transported
+ Int_t dau1 = part->GetFirstDaughter();
+ Int_t dau2 = -1;
+ if (!part->TestBit(kTransportBit)) continue;
+ //
+ fTmpTreeTR->GetEntry(it++);
+ Int_t nh = fTmpTrackReferences.GetEntries();
+ // Determine range of secondaries produced by this primary
+ if (dau1 > -1) {
+ Int_t inext = ip - 1;
+ while (dau2 < 0) {
+ if (inext >= 0) {
+ part = stack->Particle(inext);
+ dau2 = part->GetFirstDaughter();
+ if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
+// if (dau2 == -1 || dau2 < np) {
+ dau2 = -1;
+ } else {
+ dau2--;
+ }
+ } else {
+ dau2 = stack->GetNtrack() - 1;
+ }
+ inext--;
+ } // find upper bound
+ } // dau2 < 0
+// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
+ //
+ // Loop over reference hits and find secondary label
+ for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
+ for (Int_t ih = 0; ih < nh; ih++) {
+ AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
+ Int_t label = tr->Label();
+ // Skip primaries
+ if (label == ip) continue;
+ if (label > dau2 || label < dau1)
+ AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
+ if (label == id) {
+ // secondary found
+ Int_t nref = fTrackReferences.GetEntriesFast();
+ new(fTrackReferences[nref]) AliTrackReference(*tr);
+ }
+ } // hits
+ treeTR->Fill();
+ fTrackReferences.Clear();
+ ifills++;
+ } // daughters
+ } // tracks
+ //
+ // Now loop again and write the primaries
+ it = nt - 1;
+ for (Int_t ip = 0; ip < np; ip++) {
+ TParticle* part = stack->Particle(ip);
+// if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np)
+ if (part->TestBit(kTransportBit))
+ {
+ // Skip particles that have not been transported
+ fTmpTreeTR->GetEntry(it--);
+ Int_t nh = fTmpTrackReferences.GetEntries();
+ //
+ // Loop over reference hits and find primary labels
+ for (Int_t ih = 0; ih < nh; ih++) {
+ AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
+ Int_t label = tr->Label();
+ if (label == ip) {
+ Int_t nref = fTrackReferences.GetEntriesFast();
+ new(fTrackReferences[nref]) AliTrackReference(*tr);
+ }
+ }
+ }
+ treeTR->Fill();
+ fTrackReferences.Clear();
+ ifills++;
+ } // tracks
+ // Check
+ if (ifills != stack->GetNtrack())
+ AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
+//
+// Clean-up
+ delete fTmpTreeTR;
+ fTmpFileTR->Close();
+ delete fTmpFileTR;
+ fTmpTrackReferences.Clear();
+ gSystem->Exec("rm -rf TrackRefsTmp.root");
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff