diff --git a/macro/run/run_tra_file_sis19.C b/macro/run/run_tra_file_sis19.C new file mode 100644 index 0000000000000000000000000000000000000000..a4cd6c2e77f5fd68ff116ca2b72cb3821cffab40 --- /dev/null +++ b/macro/run/run_tra_file_sis19.C @@ -0,0 +1,231 @@ +/* Copyright (C) 2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Volker Friese [committer] */ + +/** @file run_tra_file_sis18.C + ** @author Volker Friese <v.friese@gsi.de> + ** @since 3 November 2020 + **/ + +// clang-format off + +// Includes needed for IDE +#if !defined(__CLING__) +#include "CbmTransport.h" + +#include "FairSystemInfo.h" + +#include "TStopwatch.h" +#endif + + +/** @brief Macro for transport simulation of events from file + ** @author Volker Friese <v.friese@gsi.de> + ** @since 3 November 2020 + ** @param input Name of input file + ** @param nEvents Number of events to transport + ** @param output Name of output data set + ** @param setup Name of predefined geometry setup + ** @param engine kGeant3 or kGeant4 + ** + ** This macro performs the transport simulation of externally generated + ** events from an input file. The macro demonstrates the minimal settings + ** needed to perform the transport: input and output file names, + ** geometry setup, number of events and the transport engine. These can be + ** specified through the arguments. Also for the arguments, defaults are + ** defined, such that the macro should execute from scratch with no arguments. + ** + ** First argument [input]: input file name. If left empty, a default input + ** file (urqmd) distributed with cbmroot will be used. If the file name + ** starts with "pluto", CbmPlutoGenerator will be used, otherwise + ** the input is assumed to be in unigen format. Other file formats require + ** intervention into the macro code. + ** + ** Second argument: [nEvents]: Number of events to be processed. + ** + ** Third argument [output]: specifies the names of the output files: + ** - [output].tra.root: Output MC data (MCPoints and MCTracks) + ** - [output].par.root: Parameter file with e.g. the geometry and field in + ** - [output].geo.root: Used geometry in TGeoManager format + ** If left empty, the term "test" will be used. + ** + ** Fourth argument [setup]: specifies the geometry setup to be used in the + ** transport run. This has to be chosen from the available predefined + ** setups (see the folder geometry/setup). By default, the setup + ** sis100_electron will be used. + ** + ** Fifth argument [engine]: transport engine. The user can choose between + ** GEANT3 (default) and GEANT4. + ** + ** Some other required definitions are coded in the macro body: + ** the target properties and the beam profile. The event vertex will be + ** sampled along the intersection of the beam trajectory with the target + ** volume. A beam trajectory is sampled for each event from the specified + ** beam profile. The beam angle is zero. Random rotation of the event + ** plane around the z axis is activated. + ** + ** All other settings are the default ones in the steering class + ** CbmTransport. For the options, consult the documentation of that class: + ** http://computing.gitpages.cbm.gsi.de/cbmroot/classCbmTransport.html + **/ + +void SetTrack(CbmTransport*, Double_t, Int_t, Double_t, Double_t, Double_t); + +void run_tra_file_sis19(const char* input = "", Int_t nEvents = 1, const char* output = "data/test", + const char* setup = "sis19_hadron", ECbmEngine engine = kGeant3, int randomSeed = 0, + Bool_t overwrite = kTRUE) +{ + + // --- Logger settings ---------------------------------------------------- + FairLogger::GetLogger()->SetLogScreenLevel("INFO"); + FairLogger::GetLogger()->SetLogVerbosityLevel("LOW"); + // ------------------------------------------------------------------------ + + // ----- Environment -------------------------------------------------- + TString myName = "run_tra_file_sis18"; // this macro's name for screen output + TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + + // ----- Input file --------------------------------------------------- + // Use default input if not specified by the user + TString defaultInput = srcDir + "/input/urqmd.auau.10gev.centr.root"; + TString inFile = (strcmp(input, "") == 0 ? defaultInput : TString(input)); + std::cout << "-I- " << myName << ": Using input file " << inFile << std::endl; + // ------------------------------------------------------------------------ + + + // ----- Output files ------------------------------------------------- + TString dataset(output); + if (dataset.IsNull()) dataset = "test"; + TString outFile = dataset + ".tra.root"; + TString parFile = dataset + ".par.root"; + TString geoFile = dataset + ".geo.root"; + std::cout << std::endl; + // ------------------------------------------------------------------------ + + + // ----- Target properties -------------------------------------------- + // For flexibility, the target volume is not part of the predefined + // geometry setup. It thus has to be specified in this macro. By default, + // a Gold target of 250 micrometer thickness and 2.5 cm diameter is used. + // The target shifts by 4 cm upstream by decision of the technical board + // in Apr 2020 which is 40 cm from center of magnet. + const char* targetMedium = "Gold"; + Double_t targetThickness = 0.025; // in cm + Double_t targetDiameter = 2.5; // in cm + + Double_t targetZpos = -35.5; // position 5 mm downstream of 1st MVD v20d_tr station + // Double_t targetZpos = -36.0; // position of 1st MVD v20d_tr station + // Double_t targetZpos = 12.5; // SIS18 5 cm upstream of 4th STS station + // Double_t targetZpos = 7.0; // SIS18 location of 3rd STS station + // Double_t targetZpos = -12.5; // SIS18 30 cm upstream of 4th STS station + // Double_t targetZpos = -44.0; + // The target position is at z=0 for the old coordinate system but is intended + // to be moved to -4cm in the DEC21 release. The global coordinate system will + // also shift from the old target position to the center of the magnet which + // is a net displacement of -40 cm. In terms of the new coordainte system + // the target is therefore to be at -44 cm. In order not to cause forgetting + // we will automate the shifting process for a short time, until the full move + // has been completed. + // if (strstr(setup, "_APR21")) targetZpos = 0.0; + // if (strstr(setup, "_DEC21")) targetZpos = -44.0; + std::cout << "Target is at " << targetZpos << "cm from origin" << std::endl; + // ------------------------------------------------------------------------ + + // ----- Beam properties ---------------------------------------------- + // The beam particles themselves are not subject to the transport, + // but the properties of the beam determine where the main collision + // vertex is placed and how the event is rotated. The collision vertex + // is determined from the intersection of the beam through the target. + // The beam defined here is along the z axis and its distributions in + // (x,y) are Gaussian with mean values 0 and sigma 1 mm. + Double_t beamPosX = 0.; + Double_t beamPosY = 0.; + Double_t beamSigmaX = 0.1; // cm + Double_t beamSigmaY = 0.1; // cm + // ------------------------------------------------------------------------ + + + // ----- Event plane rotation ----------------------------------------- + // Some generators, e.g., UrQMD, generate the event plane always at + // zero angle, i.e., the impact parameter is in the x direction. In order + // get a realistic distribution, the event plane has to be randomly + // rotated around the z axis. This changes the momentum direction of + // all primary particles coherently. + Bool_t rotateEvents = kFALSE; // kTRUE; + // ------------------------------------------------------------------------ + + + // ----- Timer -------------------------------------------------------- + TStopwatch timer; + timer.Start(); + // ------------------------------------------------------------------------ + + + // --- Transport run ---------------------------------------------------- + CbmTransport run; + + // run.AddInput(new FairParticleGenerator(-13, 1, -72.0, +72.0, 72.0)); + + // ACC // geometrical acceptance + Double_t beamRotY = 0; + // ACC // + SetTrack(&run, beamRotY,-13, -72.0, +72.0, 72.0); + SetTrack(&run, beamRotY,-13, -72.0, 0.0, 72.0); + SetTrack(&run, beamRotY,-13, -72.0, -72.0, 72.0); + // ACC // + SetTrack(&run, beamRotY, 13, 0.0, +72.0, 72.0); + SetTrack(&run, beamRotY, 13, 0.0, 0.0, 72.0); + SetTrack(&run, beamRotY, 13, 0.0, -72.0, 72.0); + // ACC // + SetTrack(&run, beamRotY,-13, +72.0, +72.0, 72.0); + SetTrack(&run, beamRotY,-13, +72.0, 0.0, 72.0); + SetTrack(&run, beamRotY,-13, +72.0, -72.0, 72.0); + + run.SetEngine(engine); + // comment the following line to remove target interaction +// if (inFile.Contains("pluto", TString::kIgnoreCase)) { run.AddInput(inFile, kPluto); } +// else +// run.AddInput(inFile, kUnigen); + run.SetOutFileName(outFile, overwrite); + run.SetParFileName(parFile); + run.SetGeoFileName(geoFile); + run.LoadSetup(setup); + run.SetTarget(targetMedium, targetThickness, targetDiameter, 0, 0, targetZpos); + run.SetBeamPosition(beamPosX, beamPosY, beamSigmaX, beamSigmaY); + if (rotateEvents) run.SetRandomEventPlane(); + run.SetRandomSeed(randomSeed); + run.StoreTrajectories(); + run.Run(nEvents); + // ------------------------------------------------------------------------ + + + // ----- Finish ------------------------------------------------------- + timer.Stop(); + Double_t rtime = timer.RealTime(); + Double_t ctime = timer.CpuTime(); + std::cout << "Macro finished successfully." << std::endl; + std::cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << std::endl; + FairSystemInfo sysInfo; + Float_t maxMemory = sysInfo.GetMaxMemory(); + std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">"; + std::cout << maxMemory; + std::cout << "</DartMeasurement>" << std::endl; + Float_t cpuUsage = ctime / rtime; + std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">"; + std::cout << cpuUsage; + std::cout << "</DartMeasurement>" << std::endl; + // ------------------------------------------------------------------------ + +} // End of macro + +void SetTrack(CbmTransport* run, Double_t beamRotY, Int_t pdgid, Double_t x, Double_t y, Double_t z) +{ + TVector3 v; + v.SetXYZ(x, y, z); + v.RotateY(-beamRotY * acos(-1.) / 180.); + cout << "X " << v.X() << " Y " << v.Y() << " Z " << v.Z() << endl; + + run->AddInput(new FairParticleGenerator(pdgid, 1, v.X(), v.Y(), v.Z())); // single electron along beam axis +}