#include <string.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 <TMCVerbose.h>
#include <TTree.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 "AliMagF.h"
#include "AliRun.h"
#include "AliSimulation.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliTransportMonitor.h"
using std::endl;
using std::cout;
ClassImp(AliMC)
AliMC::AliMC() :
fGenerator(0),
fSaveRndmStatus(kFALSE),
fSaveRndmEventStatus(kFALSE),
fReadRndmStatus(kFALSE),
fUseMonitoring(kFALSE),
fRndmFileName("random.root"),
fEventEnergy(0),
fSummEnergy(0),
fSum2Energy(0),
fTrRmax(1.e10),
fTrZmax(1.e10),
fRDecayMax(1.e10),
fRDecayMin(-1.),
fDecayPdg(0),
fImedia(0),
fTransParName("\0"),
fMonitor(0),
fHitLists(0),
fTmpTreeTR(0),
fTmpFileTR(0),
fTrackReferences(),
fTmpTrackReferences()
{
DecayLimits();
}
AliMC::AliMC(const char *name, const char *title) :
TVirtualMCApplication(name, title),
fGenerator(0),
fSaveRndmStatus(kFALSE),
fSaveRndmEventStatus(kFALSE),
fReadRndmStatus(kFALSE),
fUseMonitoring(kFALSE),
fRndmFileName("random.root"),
fEventEnergy(0),
fSummEnergy(0),
fSum2Energy(0),
fTrRmax(1.e10),
fTrZmax(1.e10),
fRDecayMax(1.e10),
fRDecayMin(-1.),
fDecayPdg(0),
fImedia(new TArrayI(1000)),
fTransParName("\0"),
fMonitor(0),
fHitLists(new TList()),
fTmpTreeTR(0),
fTmpFileTR(0),
fTrackReferences("AliTrackReference", 100),
fTmpTrackReferences("AliTrackReference", 100)
{
SetTransPar();
DecayLimits();
for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
}
AliMC::~AliMC()
{
delete fGenerator;
delete fImedia;
delete fHitLists;
delete fMonitor;
}
void AliMC::ConstructGeometry()
{
if(AliSimulation::Instance()->IsGeometryFromFile()){
if(IsGeometryFromCDB()){
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{
const char *geomfilename = AliSimulation::Instance()->GetGeometryFile();
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;
}
}
TVirtualMC::GetMC()->SetRootGeometry();
}else{
if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
TStopwatch stw;
TIter next(gAlice->Modules());
AliModule *detector;
AliDebug(1, "Geometry creation:");
while((detector = dynamic_cast<AliModule*>(next()))) {
stw.Start();
detector->CreateMaterials();
detector->CreateGeometry();
AliInfo(Form("%10s R:%.2fs C:%.2fs",
detector->GetName(),stw.RealTime(),stw.CpuTime()));
}
}
}
Bool_t AliMC::MisalignGeometry()
{
if(!AliSimulation::Instance()->IsGeometryFromFile()){
SetAllAlignableVolumes();
}
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()
{
TIter next(gAlice->Modules());
AliModule *detector;
AliInfo("Optical properties definition");
while((detector = dynamic_cast<AliModule*>(next()))) {
if(AliSimulation::Instance()->IsGeometryFromFile()) detector->CreateMaterials();
detector->DefineOpticalProperties();
}
}
#include <TPDGCode.h>
void AliMC::AddParticles()
{
cout << "########## AliMC::AddParticles" << endl;
TVirtualMC::GetMC()->DefineParticle(1010010030, "HyperTriton", kPTHadron, 2.99131 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1010010030, "AntiHyperTriton", kPTHadron, 2.99131 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1010010040, "Hyperhydrog4", kPTHadron, 3.931 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1010010040, "AntiHyperhydrog4", kPTHadron, 3.931 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1010020040, "Hyperhelium4", kPTHadron, 3.929 , 2.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1010020040, "AntiHyperhelium4", kPTHadron, 3.929 , 2.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1010000020, "LambdaNeutron", kPTNeutron, 2.054 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1010000020, "AntiLambdaNeutron", kPTNeutron, 2.054 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1020000020, "Hdibaryon", kPTNeutron, 2.23, 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1020000020, "AntiHdibaryon", kPTNeutron, 2.23, 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1030000020, "Xi0Proton", kPTHadron, 2.248 , 1.0, 1.333e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1030000020, "AntiXi0Proton", kPTHadron, 2.248 , 1.0, 1.333e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(1010000030, "LambdaNeutronNeutron", kPTNeutron, 2.982 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(-1010000030, "AntiLambdaNeutronNeutron", kPTNeutron, 2.982 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
TVirtualMC::GetMC()->DefineParticle(9010221, "f0_980", kPTNeutron, 0.98 , 0.0, 9.403e-24,"Hadron", 7e-2, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
TVirtualMC::GetMC()->DefineParticle(225, "f2_1270", kPTNeutron, 1.275 , 0.0, 3.558e-24,"Hadron", 0.185, 0, 0, 0, 0, 0, 0, 0, 0, kTRUE);
Int_t mode[6][3];
Float_t bratio[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio[kz] = 0.;
mode[kz][0] = 0;
mode[kz][1] = 0;
mode[kz][2] = 0;
}
bratio[0] = 50.;
mode[0][0] = 1000020030;
mode[0][1] = -211;
bratio[1] = 50.;
mode[1][0] = 1000010020;
mode[1][1] = 2212;
mode[1][2] = -211;
TVirtualMC::GetMC()->SetDecayMode(1010010030,bratio,mode);
Int_t amode[6][3];
Float_t abratio[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio[kz] = 0.;
amode[kz][0] = 0;
amode[kz][1] = 0;
amode[kz][2] = 0;
}
abratio[0] = 50.;
amode[0][0] = -1000020030;
amode[0][1] = 211;
abratio[1] = 50.;
amode[1][0] = -1000010020;
amode[1][1] = -2212;
amode[1][2] = 211;
TVirtualMC::GetMC()->SetDecayMode(-1010010030,abratio,amode);
Int_t mode3[6][3];
Float_t bratio3[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio3[kz] = 0.;
mode3[kz][0] = 0;
mode3[kz][1] = 0;
mode3[kz][2] = 0;
}
bratio3[0] = 50.;
mode3[0][0] = 1000020040;
mode3[0][1] = -211;
bratio3[1] = 50.;
mode3[1][0] = 1000010030;
mode3[1][1] = 2212;
mode3[1][2] = -211;
TVirtualMC::GetMC()->SetDecayMode(1010010040,bratio3,mode3);
Int_t amode3[6][3];
Float_t abratio3[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio3[kz] = 0.;
amode3[kz][0] = 0;
amode3[kz][1] = 0;
amode3[kz][2] = 0;
}
abratio3[0] = 50.;
amode3[0][0] = -1000020040;
amode3[0][1] = 211;
abratio3[1] = 50.;
amode3[1][0] = -1000010030;
amode3[1][1] = -2212;
amode3[1][2] = 211;
TVirtualMC::GetMC()->SetDecayMode(-1010010040,abratio3,amode3);
Int_t mode4[6][3];
Float_t bratio4[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio4[kz] = 0.;
mode4[kz][0] = 0;
mode4[kz][1] = 0;
mode4[kz][2] = 0;
}
bratio4[0] = 100.;
mode4[0][0] = 1000020030;
mode4[0][1] = -211;
mode4[0][2] = 2212;
TVirtualMC::GetMC()->SetDecayMode(1010020040,bratio4,mode4);
Int_t amode4[6][3];
Float_t abratio4[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio4[kz] = 0.;
amode4[kz][0] = 0;
amode4[kz][1] = 0;
amode4[kz][2] = 0;
}
abratio4[0] = 100.;
amode4[0][0] = -1000020030;
amode4[0][1] = 211;
amode4[0][2] = -2212;
TVirtualMC::GetMC()->SetDecayMode(-1010020040,abratio4,amode4);
Int_t mode1[6][3];
Float_t bratio1[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio1[kz] = 0.;
mode1[kz][0] = 0;
mode1[kz][1] = 0;
mode1[kz][2] = 0;
}
bratio1[0] = 100.;
mode1[0][0] = 1000010020;
mode1[0][1] = -211;
TVirtualMC::GetMC()->SetDecayMode(1010000020,bratio1,mode1);
Int_t amode1[6][3];
Float_t abratio1[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio1[kz] = 0.;
amode1[kz][0] = 0;
amode1[kz][1] = 0;
amode1[kz][2] = 0;
}
abratio1[0] = 100.;
amode1[0][0] = -1000010020;
amode1[0][1] = 211;
TVirtualMC::GetMC()->SetDecayMode(-1010000020,abratio1,amode1);
Int_t mode2[6][3];
Float_t bratio2[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio2[kz] = 0.;
mode2[kz][0] = 0;
mode2[kz][1] = 0;
mode2[kz][2] = 0;
}
bratio2[0] = 100.;
mode2[0][0] = 3122;
mode2[0][1] = 2212;
mode2[0][2] = -211;
TVirtualMC::GetMC()->SetDecayMode(1020000020,bratio2,mode2);
Int_t amode2[6][3];
Float_t abratio2[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio2[kz] = 0.;
amode2[kz][0] = 0;
amode2[kz][1] = 0;
amode2[kz][2] = 0;
}
abratio2[0] = 100.;
amode2[0][0] = -3122;
amode2[0][1] = -2212;
amode2[0][2] = 211;
TVirtualMC::GetMC()->SetDecayMode(-1020000020,abratio2,amode2);
Int_t mode5[6][3];
Float_t bratio5[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio5[kz] = 0.;
mode5[kz][0] = 0;
mode5[kz][1] = 0;
mode5[kz][2] = 0;
}
bratio5[0] = 100.;
mode5[0][0] = 3122;
mode5[0][1] = 2212;
TVirtualMC::GetMC()->SetDecayMode(1030000020,bratio5,mode5);
Int_t amode5[6][3];
Float_t abratio5[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio5[kz] = 0.;
amode5[kz][0] = 0;
amode5[kz][1] = 0;
amode5[kz][2] = 0;
}
abratio5[0] = 100.;
amode5[0][0] = -3122;
amode5[0][1] = -2212;
TVirtualMC::GetMC()->SetDecayMode(-1030000020,abratio5,amode5);
Int_t mode6[6][3];
Float_t bratio6[6];
for (Int_t kz = 0; kz < 6; kz++) {
bratio6[kz] = 0.;
mode6[kz][0] = 0;
mode6[kz][1] = 0;
mode6[kz][2] = 0;
}
bratio6[0] = 100.;
mode6[0][0] = 1000010030;
mode6[0][1] = -211;
TVirtualMC::GetMC()->SetDecayMode(1010000030,bratio6,mode6);
Int_t amode6[6][3];
Float_t abratio6[6];
for (Int_t kz = 0; kz < 6; kz++) {
abratio6[kz] = 0.;
amode6[kz][0] = 0;
amode6[kz][1] = 0;
amode6[kz][2] = 0;
}
abratio6[0] = 100.;
amode6[0][0] = -1000010030;
amode6[0][1] = 211;
TVirtualMC::GetMC()->SetDecayMode(-1010000030,abratio6,amode6);
for (Int_t kz = 0; kz < 6; kz++) {
bratio[kz] = 0.;
mode[kz][0] = 0;
mode[kz][1] = 0;
mode[kz][2] = 0;
}
bratio[0] = 100.;
mode[0][0] = 211;
mode[0][1] = -211;
TVirtualMC::GetMC()->SetDecayMode(9010221,bratio,mode);
for (Int_t kz = 0; kz < 6; kz++) {
bratio[kz] = 0.;
mode[kz][0] = 0;
mode[kz][1] = 0;
mode[kz][2] = 0;
}
bratio[0] = 100.;
mode[0][0] = 211;
mode[0][1] = -211;
TVirtualMC::GetMC()->SetDecayMode(225,bratio,mode);
}
void AliMC::InitGeometry()
{
AliInfo("Initialisation:");
TStopwatch stw;
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
stw.Start();
detector->Init();
AliInfo(Form("%10s R:%.2fs C:%.2fs",
detector->GetName(),stw.RealTime(),stw.CpuTime()));
}
}
void AliMC::SetGeometryFromCDB()
{
if(AliCDBManager::Instance()->IsDefaultStorageSet() &&
AliCDBManager::Instance()->GetRun() >= 0)
AliSimulation::Instance()->SetGeometryFile("*OCDB*");
else
AliError("Loading of geometry from CDB ignored. First set a default CDB storage!");
}
Bool_t AliMC::IsGeometryFromCDB() const
{
return (strcmp(AliSimulation::Instance()->GetGeometryFile(),"*OCDB*")==0);
}
void AliMC::SetAllAlignableVolumes()
{
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()
{
Generator()->Generate();
}
void AliMC::SetGenerator(AliGenerator *generator)
{
if(!fGenerator) fGenerator = generator;
}
void AliMC::ResetGenerator(AliGenerator *generator)
{
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()
{
AliDebug(1, "fGenerator->FinishRun()");
fGenerator->FinishRun();
if (fMonitor) {
fMonitor->Print();
fMonitor->Export("timing.root");
}
AliDebug(1, "EnergySummary()");
ToAliDebug(1, EnergySummary());
}
void AliMC::BeginPrimary()
{
ResetHits();
ResetTrackReferences();
}
void AliMC::PreTrack()
{
TObjArray &dets = *gAlice->Modules();
AliModule *module;
for(Int_t i=0; i<=gAlice->GetNdets(); i++)
if((module = dynamic_cast<AliModule*>(dets[i])))
module->PreTrack();
}
void AliMC::Stepping()
{
Int_t id = DetFromMate(TVirtualMC::GetMC()->CurrentMedium());
if (id < 0) return;
if ( TVirtualMC::GetMC()->IsNewTrack() &&
TVirtualMC::GetMC()->TrackTime() == 0. &&
fRDecayMin >= 0. &&
fRDecayMax > fRDecayMin &&
TVirtualMC::GetMC()->TrackPid() == fDecayPdg )
{
FixParticleDecaytime();
}
if (fUseMonitoring) {
if (!fMonitor) {
fMonitor = new AliTransportMonitor(TVirtualMC::GetMC()->NofVolumes()+1);
fMonitor->Start();
}
if (TVirtualMC::GetMC()->IsNewTrack() || TVirtualMC::GetMC()->TrackTime() == 0. || TVirtualMC::GetMC()->TrackStep()<1.1E-10) {
fMonitor->DummyStep();
} else {
Int_t copy;
Int_t volId = TVirtualMC::GetMC()->CurrentVolID(copy);
Int_t pdg = TVirtualMC::GetMC()->TrackPid();
TLorentzVector xyz, pxpypz;
TVirtualMC::GetMC()->TrackPosition(xyz);
TVirtualMC::GetMC()->TrackMomentum(pxpypz);
fMonitor->StepInfo(volId, pdg, pxpypz.E(), xyz.X(), xyz.Y(), xyz.Z());
}
}
if (AliSimulation::Instance()->Lego())
AliSimulation::Instance()->Lego()->StepManager();
else {
Int_t copy;
AddEnergyDeposit(TVirtualMC::GetMC()->CurrentVolID(copy),TVirtualMC::GetMC()->Edep());
if (TVirtualMC::GetMC()->IsTrackDisappeared() && !(TVirtualMC::GetMC()->IsTrackAlive())) {
if (TVirtualMC::GetMC()->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(),
AliTrackReference::kDisappeared);
}
AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
if(det && det->StepManagerIsEnabled()) {
det->StepManager();
}
}
}
void AliMC::EnergySummary()
{
Int_t ndep=0;
Float_t edtot=0;
Float_t ed, ed2;
Int_t kn, i, left, j, id;
const Float_t kzero=0;
Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1;
if(ievent) {
printf("***************** Energy Loss Information per event (GEV) *****************\n");
for(kn=1;kn<fEventEnergy.GetSize();kn++) {
ed=fSummEnergy[kn];
if(ed>0) {
fEventEnergy[ndep]=kn;
if(ievent>1) {
ed=ed/ievent;
ed2=fSum2Energy[kn];
ed2=ed2/ievent;
ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
} else
ed2=99;
fSummEnergy[ndep]=ed;
fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
edtot+=ed;
ndep++;
}
}
for(kn=0;kn<(ndep-1)/3+1;kn++) {
left=ndep-kn*3;
for(i=0;i<(3<left?3:left);i++) {
j=kn*3+i;
id=Int_t (fEventEnergy[j]+0.1);
printf(" %s %10.3f +- %10.3f%%;",TVirtualMC::GetMC()->VolName(id),fSummEnergy[j],fSum2Energy[j]);
}
printf("\n");
}
printf("******************** Relative Energy Loss per event ********************\n");
printf("Total energy loss per event %10.3f GeV\n",edtot);
for(kn=0;kn<(ndep-1)/5+1;kn++) {
left=ndep-kn*5;
for(i=0;i<(5<left?5:left);i++) {
j=kn*5+i;
id=Int_t (fEventEnergy[j]+0.1);
printf(" %s %10.3f%%;",TVirtualMC::GetMC()->VolName(id),100*fSummEnergy[j]/edtot);
}
printf("\n");
}
for(kn=0;kn<75;kn++) printf("*");
printf("\n");
}
}
#include <TFile.h>
void AliMC::BeginEvent()
{
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliDebug(1, " BEGINNING EVENT ");
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
AliRunLoader *runloader=AliRunLoader::Instance();
gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
runloader->SetEventNumber(gAlice->GetEventNrInRun());
AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
fEventEnergy.Reset();
if (runloader->Stack())
runloader->Stack()->Reset();
else
runloader->MakeStack();
if ( fSaveRndmStatus || fSaveRndmEventStatus) {
TString fileName="random";
if ( fSaveRndmEventStatus ) {
fileName += "Evt";
fileName += gAlice->GetEventNrInRun();
}
fileName += ".root";
cout << "Saving random engine status in " << fileName.Data() << endl;
TFile f(fileName.Data(),"RECREATE");
gRandom->Write(fileName.Data());
}
if ( fReadRndmStatus ) {
cout << "Reading random engine status from " << fRndmFileName.Data() << endl;
TFile f(fRndmFileName.Data());
gRandom->Read(fRndmFileName.Data());
}
if(AliSimulation::Instance()->Lego() == 0x0)
{
AliDebug(1, "fRunLoader->MakeTree(K)");
runloader->MakeTree("K");
}
AliDebug(1, "TVirtualMC::GetMC()->SetStack(fRunLoader->Stack())");
TVirtualMC::GetMC()->SetStack(runloader->Stack());
runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(),
gAlice->GetEventNrInRun());
if(AliSimulation::Instance()->Lego())
{
AliSimulation::Instance()->Lego()->BeginEvent();
return;
}
AliDebug(1, "ResetHits()");
ResetHits();
AliDebug(1, "fRunLoader->MakeTree(H)");
runloader->MakeTree("H");
MakeTmpTrackRefsTree();
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = (AliModule*)next()))
{
AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
detector->MakeBranch("H");
}
}
void AliMC::ResetHits()
{
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->ResetHits();
}
}
void AliMC::ResetDigits()
{
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->ResetDigits();
}
}
void AliMC::ResetSDigits()
{
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->ResetSDigits();
}
}
void AliMC::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();
}
void AliMC::FinishPrimary()
{
AliRunLoader *runloader=AliRunLoader::Instance();
#if ROOT_VERSION_CODE > 262152
if (!(TVirtualMC::GetMC()->SecondariesAreOrdered())) {
if (runloader->Stack()->ReorderKine()) RemapHits();
}
#endif
if (runloader->Stack()->PurifyKine()) RemapHits();
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->FinishPrimary();
AliLoader* loader = detector->GetLoader();
if(loader)
{
TTree* treeH = loader->TreeH();
if (treeH) treeH->Fill();
}
}
if (fTmpTreeTR) fTmpTreeTR->Fill();
}
void AliMC::RemapHits()
{
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()));
}
}
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());
}
void AliMC::FinishEvent()
{
if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent();
TIter next(gAlice->Modules());
AliModule *detector;
while((detector = dynamic_cast<AliModule*>(next()))) {
detector->FinishEvent();
}
Int_t i;
for(i=0;i<fEventEnergy.GetSize();i++)
{
fSummEnergy[i]+=fEventEnergy[i];
fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
}
AliRunLoader *runloader=AliRunLoader::Instance();
AliHeader* header = runloader->GetHeader();
AliStack* stack = runloader->Stack();
if ( (header == 0x0) || (stack == 0x0) )
{
AliFatal("Can not get the stack or header from LOADER");
return;
}
header->SetNprimary(stack->GetNprimary());
header->SetNtrack(stack->GetNtrack());
header->SetTimeStamp(AliSimulation::Instance()->GenerateTimeStamp());
if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
if (fTmpTreeTR) ReorderAndExpandTreeTR();
TTree* treeE = runloader->TreeE();
if (treeE)
{
header->SetStack(stack);
treeE->Fill();
}
else
{
AliError("Can not get TreeE from RL");
}
if(AliSimulation::Instance()->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::Init()
{
TVirtualMC::GetMC()->Init();
#if ROOT_VERSION_CODE < 331527
SetAllAlignableVolumes();
#endif
ReadTransPar();
MediaTable();
TVirtualMC::GetMC()->BuildPhysics();
fEventEnergy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
fSummEnergy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
fSum2Energy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
AliConfig::Instance()->Add(TVirtualMC::GetMC());
}
void AliMC::MediaTable()
{
Int_t kz, nz, idt, lz, i, k, ind;
TObjArray &dets = *gAlice->Detectors();
AliModule *det;
Int_t ndets=gAlice->GetNdets();
for (kz=0;kz<ndets;kz++) {
if((det=dynamic_cast<AliModule*>(dets[kz]))) {
TArrayI &idtmed = *(det->GetIdtmed());
for(nz=0;nz<100;nz++) {
if((idt=idtmed[nz])) {
det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
}
}
if(det->LoMedium() > det->HiMedium()) {
det->LoMedium() = 0;
det->HiMedium() = 0;
} else {
if(det->HiMedium() > fImedia->GetSize()) {
AliError(Form("Increase fImedia from %d to %d",
fImedia->GetSize(),det->HiMedium()));
return;
}
for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
(*fImedia)[lz]=kz;
}
}
}
}
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;
det=dynamic_cast<AliModule*>(dets[ind]);
if(det)
printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
det->HiMedium());
else
printf(" %6s: %3d -> %3d;","NULL",0,0);
}
printf("\n");
}
);
}
void AliMC::ReadTransPar()
{
const Int_t kncuts=10;
const Int_t knflags=12;
const Int_t knpars=kncuts+knflags;
const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
"BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
"BREM","COMP","DCAY","DRAY","HADR","LOSS",
"MULS","PAIR","PHOT","RAYL","STRA"};
char line[256];
char detName[7];
char* filtmp;
Float_t cut[kncuts];
Int_t flag[knflags];
Int_t i, itmed, iret, jret, ktmed, kz;
FILE *lun;
filtmp=gSystem->ExpandPathName(fTransParName.Data());
lun=fopen(filtmp,"r");
delete [] filtmp;
if(!lun) {
AliWarning(Form("File %s does not exist!",fTransParName.Data()));
return;
}
while(1) {
for(i=0;i<kncuts;i++) cut[i]=-99;
for(i=0;i<knflags;i++) flag[i]=-99;
itmed=0;
memset(line,0,256);
iret=fscanf(lun,"%255[^\n]",line);
if(iret<0) {
fclose(lun);
return;
}
jret = fscanf(lun,"%*c");
if(!iret) continue;
if(line[0]=='*') continue;
iret=sscanf(line,"%6s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d %d",
detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
&cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
&flag[8],&flag[9],&flag[10],&flag[11]);
if(!iret) continue;
if(iret<0) {
AliWarning(Form("Error reading file %s",fTransParName.Data()));
continue;
}
AliModule *mod = gAlice->GetModule(detName);
if(mod) {
TArrayI &idtmed = *mod->GetIdtmed();
if(0<=itmed && itmed < 100) {
ktmed=idtmed[itmed];
if(!ktmed) {
AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
continue;
}
for(kz=0;kz<kncuts;kz++) {
if(cut[kz]>=0) {
AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
kpars[kz],cut[kz],itmed,mod->GetName()));
TVirtualMC::GetMC()->Gstpar(ktmed,kpars[kz],cut[kz]);
}
}
for(kz=0;kz<knflags;kz++) {
if(flag[kz]>=0) {
AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
TVirtualMC::GetMC()->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
}
}
} else {
AliWarning(Form("Invalid medium code %d",itmed));
continue;
}
} else {
AliDebug(1, Form("%s not present",detName));
continue;
}
}
}
void AliMC::SetTransPar(const char *filename)
{
fTransParName = filename;
}
void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
{
TObjArray &dets = *gAlice->Modules();
if(dets[id]) static_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
}
void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
{
TObjArray &dets = *gAlice->Modules();
if(dets[id]) static_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
}
Int_t AliMC::GetCurrentTrackNumber() const {
return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
}
void AliMC::DumpPart (Int_t i) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
runloader->Stack()->DumpPart(i);
}
void AliMC::DumpPStack () const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
runloader->Stack()->DumpPStack();
}
Int_t AliMC::GetNtrack() const {
Int_t ntracks = -1;
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
ntracks = runloader->Stack()->GetNtrack();
return ntracks;
}
Int_t AliMC::GetPrimary(Int_t track) const
{
Int_t nprimary = -999;
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader->Stack())
nprimary = runloader->Stack()->GetPrimary(track);
return nprimary;
}
TParticle* AliMC::Particle(Int_t i) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
return runloader->Stack()->Particle(i);
return 0x0;
}
const TObjArray* AliMC::Particles() const {
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
return runloader->Stack()->Particles();
return 0x0;
}
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
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
mech, ntr, weight, is);
}
void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
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) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
polx, poly, polz, mech, ntr, weight, is);
}
void AliMC::SetHighWaterMark(Int_t nt) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->SetHighWaterMark(nt);
}
void AliMC::KeepTrack(Int_t track) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->KeepTrack(track);
}
void AliMC::FlagTrack(Int_t track) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->FlagTrack(track);
}
void AliMC::SetCurrentTrack(Int_t track) const
{
AliRunLoader * runloader = AliRunLoader::Instance();
if (runloader)
if (runloader->Stack())
runloader->Stack()->SetCurrentTrack(track);
}
AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id)
{
Int_t primary = GetPrimary(label);
Particle(primary)->SetBit(kKeepBit);
Int_t nref = fTmpTrackReferences.GetEntriesFast();
return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
}
void AliMC::ResetTrackReferences()
{
fTmpTrackReferences.Clear();
}
void AliMC::RemapTrackReferencesIDs(const Int_t *map)
{
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);
}
}
}
fTmpTrackReferences.Compress();
}
void AliMC::FixParticleDecaytime()
{
TLorentzVector p;
TVirtualMC::GetMC()->TrackMomentum(p);
Double_t tmin, tmax;
Double_t b;
Double_t vt = p.Pt() / p.E();
if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) {
Double_t rho = p.Pt() / 0.0003 / b;
Double_t omega = vt / rho;
const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.));
if (fRDecayMax * kOvRhoSqr2 > 1.) return;
tmax = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega;
tmin = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega;
} else {
tmax = fRDecayMax / vt;
tmin = fRDecayMin / vt;
}
Double_t t = tmin + (tmax - tmin) * gRandom->Rndm();
TVirtualMC::GetMC()->ForceDecayTime(t / 2.99792458e10);
}
void AliMC::MakeTmpTrackRefsTree()
{
fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
fTmpTreeTR = new TTree("TreeTR", "Track References");
TClonesArray* pRef = &fTmpTrackReferences;
fTmpTreeTR->Branch("TrackReferences", &pRef, 4000);
}
void AliMC::ReorderAndExpandTreeTR()
{
AliRunLoader *rl = AliRunLoader::Instance();
AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
rl->MakeTrackRefsContainer();
TTree * treeTR = rl->TreeTR();
TClonesArray* pRef = &fTrackReferences;
treeTR->Branch("TrackReferences", &pRef);
AliStack* stack = rl->Stack();
Int_t np = stack->GetNprimary();
Int_t nt = fTmpTreeTR->GetEntries();
Int_t ifills = 0;
Int_t it = 0;
for (Int_t ip = np - 1; ip > -1; ip--) {
TParticle *part = stack->Particle(ip);
Int_t dau1 = part->GetFirstDaughter();
Int_t dau2 = -1;
if (!part->TestBit(kTransportBit)) continue;
fTmpTreeTR->GetEntry(it++);
Int_t nh = fTmpTrackReferences.GetEntries();
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) {
dau2 = -1;
} else {
dau2--;
}
} else {
dau2 = stack->GetNtrack() - 1;
}
inext--;
}
}
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();
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) {
Int_t nref = fTrackReferences.GetEntriesFast();
new(fTrackReferences[nref]) AliTrackReference(*tr);
}
}
treeTR->Fill();
fTrackReferences.Clear();
ifills++;
}
}
it = nt - 1;
for (Int_t ip = 0; ip < np; ip++) {
TParticle* part = stack->Particle(ip);
if (part->TestBit(kTransportBit))
{
fTmpTreeTR->GetEntry(it--);
Int_t nh = fTmpTrackReferences.GetEntries();
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++;
}
if (ifills != stack->GetNtrack())
AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
delete fTmpTreeTR;
fTmpFileTR->Close();
delete fTmpFileTR;
fTmpTrackReferences.Clear();
gSystem->Exec("rm -rf TrackRefsTmp.root");
}