void sim(Int_t seed = 9999, Int_t nEvents = 1000000){
  //---------------------Files-----------------------------------------------
  TString outFile;
  outFile.Form("sim%d.root",seed);
  TString parFile;
  parFile.Form("par%d.root",seed);
  // ------------------------------------------------------------------------
  gRandom->SetSeed(seed);
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------
 
  // -----   Create simulation run   ----------------------------------------
  FairRunSim* run = new FairRunSim();
  /** Select transport engine
  * TGeant3
  * TGeant4
  **/
  run->SetName("TGeant4");              // Transport engine
  run->SetOutputFile(outFile.Data());          // Output file
  // ------------------------------------------------------------------------
  
  // -----   Runtime database   ---------------------------------------------
  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  // ------------------------------------------------------------------------
  
  // -----   Create media   -------------------------------------------------
  run->SetMaterials("media.geo");       // Materials
  // ------------------------------------------------------------------------

  //-------- Set MC event header --------------------------------------------
  ERMCEventHeader* header = new ERMCEventHeader();
  run->SetMCEventHeader(header);
  //-------------------------------------------------------------------------

  // -----   Create detectors  ---------------------------------------------- 
  FairModule* cave= new ERCave("CAVE");
  cave->SetGeometryFileName("cave.geo");
  run->AddModule(cave);
  
  ERCollimator* collimator = new ERCollimator();
  collimator->SetGeometryFileName("collimator.geo.root");
  run->AddModule(collimator);

  // ER NeuRad definition
  /* Select verbosity level
   * 1 - only standard logs
   * 2 - Print points after each event
   * 3 - - GEANT Step information
  */
  Int_t verbose = 1;
  ERNeuRad* neuRad= new ERNeuRad("ERNeuRad", kTRUE,verbose);
  neuRad->SetGeometryFileName("NeuRad.v4.geo.root");
  /* Select storing stepss
   * not store steps
   * SetStorePrimarySteps() - store only primary particle step
   * SetStoreAllSteps() - store all steps. WARNING - very slow
  */
  //neuRad->SetStorePrimarySteps();
  run->AddModule(neuRad);
  // ------------------------------------------------------------------------
  
  // -----   Create PrimaryGenerator   --------------------------------------
  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  // Int_t pdgId = 2112; // neutron  beam
  // Double32_t theta1 = 0.;  // polar angle distribution
  // Double32_t theta2 = 7.;
  // Double32_t kin_energy = .500; //GeV
  Int_t pdgId = 22; // neutron  beam
  Double32_t theta1 = 0.;  // polar angle distribution
  Double32_t theta2 = 180.;
  Double32_t kin_energy = 0.001333; //GeV

  Double_t mass = TDatabasePDG::Instance()->GetParticle(pdgId)->Mass();
  Double32_t momentum = TMath::Sqrt(kin_energy*kin_energy + 2.*kin_energy*mass); //GeV
  FairBoxGenerator* boxGen = new FairBoxGenerator(pdgId, 1);
  boxGen->SetThetaRange(theta1, theta1);
  boxGen->SetPRange(momentum, momentum);
  // boxGen->SetPhiRange(90, 90);
  boxGen->SetPhiRange(0, 0);
  // boxGen->SetBoxXYZ(0.,0,0.3,0.3,0.);
  boxGen->SetBoxXYZ(-0.149,-0.149,0.149,0.149,19.4);
  primGen->AddGenerator(boxGen);
  run->SetGenerator(primGen);
  // ------------------------------------------------------------------------
  
  //-------Set visualisation flag to true------------------------------------
  //run->SetStoreTraj(kTRUE);
  
  //-------Set LOG verbosity  ----------------------------------------------- 
  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
  
  // -----   Initialize simulation run   ------------------------------------
  run->Init();
  Int_t nSteps = -15000;
  //gMC->SetMaxNStep(nSteps);
  
  // -----   Runtime database   ---------------------------------------------
  Bool_t kParameterMerged = kTRUE;
  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
  parOut->open(parFile.Data());
  rtdb->setOutput(parOut);
  rtdb->saveOutput();
  rtdb->print();
  // ---------------------------------------------------------
  
  // -----   Run simulation  ------------------------------------------------
  run->Run(nEvents);
  
  // -----   Finish   -------------------------------------------------------
  neuRad->WriteHistos();
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Output file is sim.root" << endl;
  cout << "Parameter file is par.root" << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime
      << "s" << endl << endl;
}