run_CbmIonGenerator.C

void run_CbmIonGenerator(Int_t nEvents = 1)
{
  // ========================================================================
  //          Adjust this part according to your requirements

  // Input file
  //TString inPath = "/d/cbm03/urqmd/auau/25gev/centr/";
  //TString inFile = inPath + "urqmd.auau.25gev.centr.0000.ftn14";
  
  // Output file
  TString outFile = Form("sts.mc.root",nEvents);
  
  // Parameter file
  TString parFile = Form("params.root",nEvents);

  // Cave geometry
  TString caveGeom = "cave.geo";
  
  // Target geometry
  TString targetGeom = "target_au_250mu.geo";

  // Beam pipe geometry
  TString pipeGeom = "pipe_standard.geo";

  // Magnet geometry and field map
  TString magnetGeom = "passive/magnet_v09e.geo";
  TString fieldMap    = "field_v10e";
  Double_t fieldZ     = 50.;     // z position of field centre
  Double_t fieldScale = 1.;      // field scaling factor

  // MVD geometry
  TString mvdGeom    = "mvd/mvd_v07a.geo";

  // STS geometry
  TString stsGeom    = "sts/sts_v11a.geo";
  
  //STS geometry for the same z position of all sensors
  //TString stsGeom = "sts_same_z.geo";
  
  targetGeom = "";
  magnetGeom = "";
  stsGeom  = "";
  mvdGeom  = "";
  
  // In general, the following parts need not be touched
  // ========================================================================

  
  

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  
  
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  
  
  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libCbmBase");
  gSystem->Load("libCbmData");
  gSystem->Load("libField");
  gSystem->Load("libGen");
  gSystem->Load("libPassive");
  gSystem->Load("libMvd");
  gSystem->Load("libSts");
  gSystem->Load("libCbmGenerators"); // for CbmIonGenerator
  // ------------------------------------------------------------------------

   
 
  // -----   Create simulation run   ----------------------------------------
  FairRunSim* run = new FairRunSim();
  run->SetName("TGeant3");              // Transport engine
  run->SetOutputFile(outFile);          // Output file
  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  // ------------------------------------------------------------------------
  
  

  // -----   Create media   -------------------------------------------------
  run->SetMaterials("media.geo");       // Materials
  // ------------------------------------------------------------------------

    
  
  // -----   Create geometry   ----------------------------------------------
  FairModule* cave= new CbmCave("CAVE");
  cave->SetGeometryFileName(caveGeom);
  run->AddModule(cave);

  FairModule* pipe= new CbmPipe("PIPE");
  pipe->SetGeometryFileName(pipeGeom);
  run->AddModule(pipe);
  
  if( targetGeom != "") {
    FairModule* target= new CbmTarget("Target");
    target->SetGeometryFileName(targetGeom);
    run->AddModule(target);
  }
  if( magnetGeom != "") {
    FairModule* magnet= new CbmMagnet("MAGNET");
    magnet->SetGeometryFileName(magnetGeom);
    run->AddModule(magnet);
  }
  if( mvdGeom != ""){
    FairDetector* mvd= new CbmMvd("MVD", kTRUE);
    mvd->SetGeometryFileName(mvdGeom); 
    run->AddModule(mvd);
  }
  if( stsGeom != "") {
    FairDetector* sts= new CbmSts("STS", kTRUE);
    sts->SetGeometryFileName(stsGeom); 
    run->AddModule(sts);
  }
  // ------------------------------------------------------------------------

  
  
  // -----   Create magnetic field   ----------------------------------------
  if(magnetGeom!="") {
    CbmFieldMap* magField = new CbmFieldMapSym2(fieldMap);
    magField->SetPosition(0., 0., fieldZ);
    magField->SetScale(fieldScale);
    run->SetField(magField);
  }
  // ------------------------------------------------------------------------

  
  
  // -----   Create PrimaryGenerator   --------------------------------------
  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  //FairUrqmdGenerator*  urqmdGen = new FairUrqmdGenerator(inFile);
  //primGen->AddGenerator(urqmdGen);
  
  // CbmIonGenerator - check in Load libraries: gSystem->Load("libCbmGenerators");
  Int_t nions=10; // number of ions/event
  cout<<"@@@@@> CbmIonGenerator is ON with "<<nions<<" ion(s)/event !!!"<<endl;
  Int_t z=79, a=197, q=79; // Au-ion
  Double_t p=8.; // in AGeV/c
  Double_t vz=-0.0126; // in cm - z-pos. of vertex
  Double_t meanX=0.3;  // in cm - spatial distr. (XOY)
  Double_t meanY=-0.2; // in cm
  Double_t sigmaX = 0.078; // in cm - spatial distr. (XOY)
  Double_t sigmaY = 0.032; // in cm
  Double_t sigmatX = 0.00117885; // in rad (Px/P) - angular distr. (mean=0)
  Double_t sigmatY = 0.00094955; // in rad (Py/P)
  // Parameters of the trapezoid are set with respect to mean of Gaussian, NOT necesseraly in absolute coordinates.
  // x1 < x2 < 0 < x3 < x4
  Double_t x1=-0.10452, x2=-0.06942, x3=0.06942, x4=0.10452; // in cm - trapezoid distr.
  Double_t y1=-0.04448, y2=-0.02688, y3=0.02688, y4=0.04448; // in cm
  Double_t tX1=-0.00157966, tX2=-0.001049177, tX3=0.001049177, tX4=0.00157966; // in rad
  Double_t tY1=-0.00131987, tY2=-0.000797622, tY3=0.000797622, tY4=0.00131987; // in rad

//  CbmIonGenerator *IonGen = new CbmIonGenerator(z, a, q, nions, p, sigmaX, sigmaY, sigmatX, sigmatY);
  
//  CbmIonGenerator *IonGen = new CbmIonGenerator(z, a, q, nions, p, sigmaX, sigmaY, sigmatX, sigmatY, meanX, meanY, vz);
  
  CbmIonGenerator *IonGen = new CbmIonGenerator(z, a, q, nions, p, sigmaX, sigmaY, sigmatX, sigmatY, meanX, meanY, vz,
                                                x1, x2, x3, x4, y1, y2, y3, y4, tX1, tX2, tX3, tX4, tY1, tY2, tY3, tY4);
  
  primGen->AddGenerator(IonGen);
  
  run->SetGenerator(primGen);
  // ------------------------------------------------------------------------

//   run->SetStoreTraj(kTRUE);
  
  // -----   Initialize simulation run   ------------------------------------
  run->Init();
  // ------------------------------------------------------------------------

  
  
  // -----   Runtime database   ---------------------------------------------
  if(magnetGeom!="") {
    CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
    fieldPar->SetParameters(magField);
    fieldPar->setChanged();
    fieldPar->setInputVersion(run->GetRunId(),1);
  }
  Bool_t kParameterMerged = kTRUE;
  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
  parOut->open(parFile.Data());
  rtdb->setOutput(parOut);
  rtdb->saveOutput();
  rtdb->print();
  // ------------------------------------------------------------------------

  
   
  // -----   Start run   ----------------------------------------------------
  run->Run(nEvents);
  // ------------------------------------------------------------------------
  //run->CreateGeometryFile("data/geofile_full.root");
  
  
  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Output file is "    << outFile << endl;
  cout << "Parameter file is " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime 
       << "s" << endl << endl;
  // ------------------------------------------------------------------------
  cout << " Test passed" << endl;
  cout << " All ok " << endl;

  exit(0);
}