// --------------------------------------------------------------------------
//
// CBM Computing School, Kolkata, February 2018
// Macro for Exercise 1: Transport Simulation
// V. Friese   15/02/2018
//
// --------------------------------------------------------------------------

void transport() {

  // First, we create a FairRunSim and tell it to use GEANT3 as transport engine.
  FairRunSim* run = new FairRunSim();
  run->SetName("TGeant3");
  
    // --- Logger settings ----------------------------------------------------
  TString logLevel     = "INFO";
  TString logVerbosity = "MEDIUM";
  // ------------------------------------------------------------------------
  FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
  FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());

  
  // Let us call the output file exercise.mc.root.
  run->SetOutputFile("exercise.mc.root");
  
  // Now, let us define the input. We first have to define a FairPrimaryGenerator
  // and register it to the run:
  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  run->SetGenerator(primGen);
    
  // The event vertex shall be at (0., 0., 0.):
  primGen->SetTarget(0., 0.);
  primGen->SetBeam(0., 0., 0., 0.);
  
  // Now, we define the specific inputs. There are two of them: the first one
  // should be a FairBoxGenerator with 10 negative muons per event.
  // The PDG code of mu- is 13.
  FairBoxGenerator* gen1 = new FairBoxGenerator(13, 10);
  
  // The momentum should be randomly distributed from 1 to 10 GeV.
  gen1->SetPRange(1., 10.);
  
  // The polar angle shall be 10 degrees
  gen1->SetThetaRange(10., 10.);
  
  // The azimut angle shall be randomly distributed 
  gen1->SetPhiRange(0., 360.);
  
  // The second input is a UrQMD file, which we find in the input directory.
  // It is unigen format, so we use CbmUnigenGenerator to read it:
  TString srcDir = gSystem->Getenv("VMCWORKDIR");
  TString inputFile = srcDir + "/input/urqmd.auau.10gev.centr.root";
  CbmUnigenGenerator* gen2 = new CbmUnigenGenerator(inputFile);
  
  // Let us not forget to register the two generators to the PrimaryGenerator:
  primGen->AddGenerator(gen1);
  primGen->AddGenerator(gen2);
  
  // Good, we have set the proper input. Now comes the geometry. We are supposed
  // to start from the setup_sis100_muon_jpsi, so let us load it.
  TString setupFileName = srcDir + "/geometry/setup/setup_sis100_muon_jpsi.C";
  gROOT->LoadMacro(setupFileName);
  gROOT->ProcessLine("setup_sis100_muon_jpsi()");
  
  // The setup macro has created a CbmSetup singleton object. 
  // In order to modify the setup, we first get the CbmSetup singleton.
  CbmSetup* setup = CbmSetup::Instance();
  
  // Now we remove TOF ...
  setup->RemoveModule(kTof);
  
  // ... change the magnet geometry to v15a...
  setup->SetModule(kMagnet, "v15a");
  
  // ... and set the field map version to v12b.
  setup->SetField("v12b");
  
  // Good, the setup is defined. We need now create all detectors and modules
  // with their respective geometries. Nice that we have a macro doing that
  // for us:
  TString macroName = srcDir + "/macro/run/modules/registerSetup.C";
  gROOT->LoadMacro(macroName);
  gROOT->ProcessLine("registerSetup()");
  
  // Remember: the target has to be constructed seperately. 
  // Material is Gold, thickness is 2.5 mm, diameter is 2 cm:
  CbmTarget* target = new CbmTarget("Gold", 0.25, 2.);
  
  // The target position is the origin of the coordinate system:
  target->SetPosition(0., 0., 0.);
  
  // And the target must be registered to the run in order to be
  // included in the geometry:
  run->AddModule(target);
  
  // Now we have to take care of the magnetic field. This is
  // created through the CbmSetup object; it must then be registered to
  // the run.
  CbmFieldMap* magField = CbmSetup::Instance()->CreateFieldMap();
  run->SetField(magField);
  
  // -Trajectories Visualization (TGeoManager Only )
  // Switch this on if you want to visualize tracks in the
  // eventdisplay.
  // This is normally switch off, because of the huge files created
  // when it is switched on. 
  run->SetStoreTraj(kTRUE);

  // The run has to be given a list of materials to create the geometry:
  run->SetMaterials("media.geo");

  // Wow! We have defined our simulation. Now the run must be initialised.
  // This creates e.g. the full geometry.
  run->Init();
  
  // Now comes some stuff with parameters. We did not touch this in the
  // lecture, so we just copy it from run_mc.C:
  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  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("exercise.par.root");
  rtdb->setOutput(parOut);
  rtdb->saveOutput();

  // Set cuts for storing the trajectories.
  // Switch this on only if trajectories are stored.
  // Choose this cuts according to your needs, but be aware
  // that the file size of the output file depends on these cuts
  FairTrajFilter* trajFilter = FairTrajFilter::Instance();
  trajFilter->SetStepSizeCut(0.01); // 1 cm
  trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
  trajFilter->SetMomentumCutP(10e-3); // p_lab > 10 MeV
  trajFilter->SetEnergyCut(0., 1.02); // 0 < Etot < 1.04 GeV
  trajFilter->SetStorePrimaries(kTRUE);
  trajFilter->SetStoreSecondaries(kTRUE);

  // Now the run is started. Let the fun begin! Will it really work?
  run->Run(10);
  
  // After finishing, we want to write the geometry to a separate file:
  run->CreateGeometryFile("transport.geo.root");
  
  // Florian can tell why this is necessary.
  RemoveGeoManager();
  
  // Hooray!
  std::cout << "Macro finished successfully." << std::endl;
  	
}