-
Eoin Clerkin authored
Decision to not use doxygen for licence headers. Removes doxygen formatting and file tag.
Eoin Clerkin authoredDecision to not use doxygen for licence headers. Removes doxygen formatting and file tag.
sim.C 6.16 KiB
/* Copyright (C) 2010 Petersburg Nuclear Physics Institute named by B.P.Konstantinov of National Research Centre "Kurchatov Institute", Gatchina
SPDX-License-Identifier: GPL-3.0-only
Authors: Evgeny Kryshen [committer] */
void sim(Int_t file_nr = 0, Int_t nEvents = 1, Int_t s = 0, Int_t seed = 1)
{
Char_t filenr[4];
sprintf(filenr, "%04d", file_nr);
printf("Filenr: %s\n", filenr);
TString signal;
if (s != 3312 && s != 3334) signal = "la";
else if (s == 3312)
signal = "xi";
else if (s == 3334)
signal = "om";
// ----- Paths and file names --------------------------------------------
TString inDir = gSystem->Getenv("URQMD_INPUT_PATH");
TString inFile = inDir + "/urqmd.auau.25gev.centr." + filenr + ".ftn14";
TString outDir = TString(filenr);
outDir += "/";
gSystem->mkdir(outDir.Data());
TString outFile = outDir + signal + ".mc.root";
TString parFile = outDir + signal + ".par.root";
// ----- Geometries -----------------------------------------------------
TString caveGeom = "cave.geo";
TString pipeGeom = "pipe_standard.geo";
TString targetGeom = "target_au_250mu.geo";
TString magnetGeom = "magnet_standard.geo";
TString stsGeom = "sts_Standard_s3055AAFK5.SecD.geo";
// ----- Magnetic field -----------------------------------------------
TString fieldMap = "FieldMuonMagnet"; // name of field map
Double_t fieldZ = 50.; // field centre z position
Double_t fieldScale = 1.; // field scaling factor
gDebug = 0;
gRandom->SetSeed(seed);
// ---- Load libraries -------------------------------------------------
cout << endl << "=== much_sim.C : Loading libraries ..." << endl;
gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
basiclibs();
gROOT->LoadMacro("$VMCWORKDIR/analysis/hyperon/analysislibs.C");
analysislibs();
FairRunSim* fRun = new FairRunSim();
fRun->SetName("TGeant3");
fRun->SetOutputFile(outFile.Data());
FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
fRun->SetMaterials("media.geo");
// ----- Create detectors and passive volumes -------------------------
cout << endl << "=== much_sim.C : Create geeometry ..." << endl;
if (caveGeom != "") {
FairModule* cave = new CbmCave("CAVE");
cave->SetGeometryFileName(caveGeom);
fRun->AddModule(cave);
cout << " --- " << caveGeom << endl;
}
if (pipeGeom != "") {
FairModule* pipe = new CbmPipe("PIPE");
pipe->SetGeometryFileName(pipeGeom);
fRun->AddModule(pipe);
cout << " --- " << pipeGeom << endl;
}
if (targetGeom != "") {
FairModule* target = new CbmTarget("Target");
target->SetGeometryFileName(targetGeom);
fRun->AddModule(target);
cout << " --- " << targetGeom << endl;
}
if (magnetGeom != "") {
FairModule* magnet = new CbmMagnet("MAGNET");
magnet->SetGeometryFileName(magnetGeom);
fRun->AddModule(magnet);
cout << " --- " << magnetGeom << endl;
}
if (stsGeom != "") {
FairDetector* sts = new CbmSts("STS", kTRUE);
sts->SetGeometryFileName(stsGeom);
cout << " --- " << stsGeom << endl;
fRun->AddModule(sts);
}
// ----- Create magnetic field ----------------------------------------
cout << endl << "=== much_sim.C : Create magnetic field ..." << endl;
CbmFieldMap* magField = NULL;
if (fieldMap == "FieldActive" || fieldMap == "FieldIron") magField = new CbmFieldMapSym3(fieldMap);
else if (fieldMap == "FieldAlligator")
magField = new CbmFieldMapSym2(fieldMap);
else if (fieldMap == "FieldMuonMagnet")
magField = new CbmFieldMapSym3(fieldMap);
else {
cout << "===> ERROR: Field map " << fieldMap << " unknown! " << endl;
exit;
}
magField->SetPosition(0., 0., fieldZ);
magField->SetScale(fieldScale);
fRun->SetField(magField);
// ----- Create PrimaryGenerator --------------------------------------
cout << endl << "=== much_sim.C : Create generators ..." << endl;
FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
if (signal == "om") {
delete gRandom;
gRandom = new TRandom3();
CbmAnaHypYPtGenerator* gen = new CbmAnaHypYPtGenerator(3334);
gen->SetDistributionPt(0.1603); // 6 GeV
gen->SetDistributionY(1.277, 0.412); // 6 GeV
//gen->SetDistributionPt(0.149808); // 25 GeV
//gen->SetDistributionY(1.9875,0.546669); // 25 GeV
gen->Init();
primGen->AddGenerator(gen);
}
else if (signal == "xi") {
delete gRandom;
gRandom = new TRandom3();
CbmAnaHypYPtGenerator* gen = new CbmAnaHypYPtGenerator(3312);
gen->SetDistributionPt(0.1603); // 6 GeV
gen->SetDistributionY(1.277, 0.412); // 6 GeV
//gen->SetDistributionPt(0.154319); // 25 GeV
//gen->SetDistributionY(1.98604,0.617173); // 25 GeV
gen->Init();
primGen->AddGenerator(gen);
}
FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(inFile.Data());
primGen->AddGenerator(urqmdGen);
fRun->SetGenerator(primGen);
// ------------------------------------------------------------------------
fRun->Init();
// ------------------------------------------------------------------------
CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
fieldPar->SetParameters(magField);
fieldPar->setChanged();
fieldPar->setInputVersion(fRun->GetRunId(), 1);
Bool_t kParameterMerged = kTRUE;
FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
parOut->open(parFile);
rtdb->setOutput(parOut);
rtdb->saveOutput();
rtdb->print();
// ------------------------------------------------------------------------
TStopwatch timer;
timer.Start();
fRun->Run(nEvents);
// ----- Finish -------------------------------------------------------
timer.Stop();
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
cout << endl << endl;
cout << "=== sim.C : Macro finished successfully." << endl;
cout << "=== sim.C : Output file is " << outFile << endl;
cout << "=== sim.C : Real time used: " << rtime << "s " << endl;
cout << "=== sim.C : CPU time used : " << ctime << "s " << endl;
cout << endl << endl;
// ------------------------------------------------------------------------
}