/* $Id$ */
+// This class is extracted from the AliRun class
+// and contains all the MC-related functionality
+// The number of dependencies has to be reduced...
+// Author: F.Carminati
+// Federico.Carminati@cern.ch
+
+#include <RVersion.h>
#include <TBrowser.h>
#include <TStopwatch.h>
#include <TSystem.h>
#include <TVirtualMC.h>
+#include "TGeant3.h"
+
+#include "AliLog.h"
#include "AliDetector.h"
#include "AliGenerator.h"
#include "AliHeader.h"
#include "AliMCQA.h"
#include "AliRun.h"
#include "AliStack.h"
+#include "AliMagF.h"
#include "AliTrackReference.h"
fTrackReferences(0)
{
+ //default constructor
+ DecayLimits();
}
//_______________________________________________________________________
fHitLists(new TList()),
fTrackReferences(new TClonesArray("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()
{
+ //destructor
delete fGenerator;
delete fImedia;
delete fMCQA;
}
//_______________________________________________________________________
-void AliMC::Copy(AliMC &) const
+void AliMC::Copy(TObject &) const
{
- Fatal("Copy","Not implemented!\n");
+ //dummy Copy function
+ AliFatal("Not implemented!");
}
//_______________________________________________________________________
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()));
}
}
// Initialize detectors and display geometry
//
- printf("Initialisation:\n");
+ AliInfo("Initialisation:");
TStopwatch stw;
TIter next(gAlice->Modules());
AliModule *detector;
// Initialise detector and display geometry
detector->Init();
detector->BuildGeometry();
- 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()));
}
}
//_______________________________________________________________________
-void AliMC::GeneratePrimaries()
+void AliMC::GeneratePrimaries()
{
//
// Generate primary particles and fill them in the stack.
//
if(fGenerator)
if(generator)
- Warning("ResetGenerator","Replacing generator %s with %s\n",
- fGenerator->GetName(),generator->GetName());
+ AliWarning(Form("Replacing generator %s with %s",
+ fGenerator->GetName(),generator->GetName()))
else
- Warning("ResetGenerator","Replacing generator %s with NULL\n",
- fGenerator->GetName());
+ 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());
}
//_______________________________________________________________________
//_______________________________________________________________________
void AliMC::PreTrack()
{
+ // Actions before the track's transport
TObjArray &dets = *gAlice->Modules();
AliModule *module;
//
// Called at every step during transport
//
-
+
Int_t id = DetFromMate(gMC->GetMedium());
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())
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()) {
+ if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber());
+ }
//Call the appropriate stepping routine;
AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
//
// 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();
//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
runloader->Stack()->Reset();//clean stack -> tree is unloaded
else
runloader->MakeStack();//or make a new one
-
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(K)");
- runloader->MakeTree("K");
- if (GetDebug()) Info("BeginEvent"," gMC->SetStack(fRunLoader->Stack())");
+
+ if(gAlice->Lego() == 0x0)
+ {
+ AliDebug(1, "fRunLoader->MakeTree(K)");
+ runloader->MakeTree("K");
+ }
+
+ AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
gMC->SetStack(gAlice->GetRunLoader()->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
gAlice->GetEventNrInRun());
// fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTrackRefsContainer()");
- runloader->MakeTrackRefsContainer();//for insurance
-
- if (GetDebug()) Info("BeginEvent"," ResetHits()");
- ResetHits();
- if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(H)");
- runloader->MakeTree("H");
-
- //
if(gAlice->Lego())
{
gAlice->Lego()->BeginEvent();
return;
}
+
+
+ AliDebug(1, "ResetHits()");
+ ResetHits();
+
+ AliDebug(1, "fRunLoader->MakeTree(H)");
+ runloader->MakeTree("H");
+
+ AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
+ runloader->MakeTrackRefsContainer();//for insurance
+
//create new branches and SetAdresses
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = (AliModule*)next()))
{
- if (GetDebug()) Info("BeginEvent"," %s->MakeBranch(H)",detector->GetName());
+ AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
detector->MakeBranch("H");
- if (GetDebug()) Info("BeginEvent"," %s->MakeBranchTR()",detector->GetName());
+ AliDebug(2, Form("%s->MakeBranchTR()",detector->GetName()));
detector->MakeBranchTR();
- if (GetDebug()) Info("BeginEvent"," %s->SetTreeAddress()",detector->GetName());
+ AliDebug(2, Form("%s->SetTreeAddress()",detector->GetName()));
detector->SetTreeAddress();
}
+ // make branch for AliRun track References
+ TTree * treeTR = runloader->TreeTR();
+ if (treeTR){
+ // make branch for central track references
+ if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
+ TBranch *branch;
+ branch = treeTR->Branch("AliRun",&fTrackReferences);
+ branch->SetAddress(&fTrackReferences);
+ }
+ //
}
//_______________________________________________________________________
//_______________________________________________________________________
void AliMC::PostTrack()
{
+ // Posts tracks for each module
TObjArray &dets = *gAlice->Modules();
AliModule *module;
// static Int_t count=0;
// const Int_t times=10;
// This primary is finished, purify stack
+#if ROOT_VERSION_CODE > 262152
+ if (!(gMC->SecondariesAreOrdered()))
+ runloader->Stack()->ReorderKine();
+#endif
runloader->Stack()->PurifyKine();
-
+
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->FinishPrimary();
- if(detector->GetLoader())
- detector->GetLoader()->TreeH()->Fill();
+ AliLoader* loader = detector->GetLoader();
+ if(loader)
+ {
+ TTree* treeH = loader->TreeH();
+ if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
+ }
}
// Write out track references if any
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
}
else
{
- Error("FinishEvent","Can not get TreeE from RL");
+ AliError("Can not get TreeE from RL");
}
- 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","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
+ if(gAlice->Lego() == 0x0)
+ {
+ runloader->WriteKinematics("OVERWRITE");
+ runloader->WriteTrackRefs("OVERWRITE");
+ runloader->WriteHits("OVERWRITE");
}
+
+ 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()
{
+ // MC initialization
//=================Create Materials and geometry
gMC->Init();
- gMC->DefineParticles(); //Create standard MC particles
-
//Read the cuts for all materials
ReadTransPar();
//Build the special IMEDIA table
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;
}
}
//_______________________________________________________________________
TParticle* AliMC::Particle(Int_t i) const
{
+ // Returns the i-th particle from the stack taking into account
+ // the remaping done by PurifyKine
AliRunLoader * runloader = gAlice->GetRunLoader();
if (runloader)
if (runloader->Stack())
//_______________________________________________________________________
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)
+ TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
{
// Delegate to stack
//
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
//
}
//_______________________________________________________________________
-void AliMC::SetHighWaterMark(const Int_t nt)
+void AliMC::SetHighWaterMark(Int_t nt) const
{
//
// Set high water mark for last track in event
}
//_______________________________________________________________________
-void AliMC::KeepTrack(const Int_t track)
+void AliMC::KeepTrack(Int_t track) const
{
//
// Delegate to stack
}
//_______________________________________________________________________
-void AliMC::FlagTrack(Int_t track)
+void AliMC::FlagTrack(Int_t track) const
{
// Delegate to stack
//
}
//_______________________________________________________________________
-void AliMC::SetCurrentTrack(Int_t track)
+void AliMC::SetCurrentTrack(Int_t track) const
{
//
// Set current track number
//
// add a trackrefernce to the list
if (!fTrackReferences) {
- cerr<<"Container trackrefernce not active\n";
+ AliError("Container trackrefernce not active");
return;
}
Int_t nref = fTrackReferences->GetEntriesFast();
}
fTrackReferences->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 = gAlice->Field()->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
+ //
+ TGeant3 * geant = (TGeant3*) gMC;
+ geant->Gcphys()->sumlif = t / p.Beta() / p.Gamma();
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff