void run(Int_t nEvents = 1000)
{
  TStopwatch timer;
  timer.Start();
  gDebug = 0;
  
  // Use non default gconfig and geometry directories
  TString dir = getenv("VMCWORKDIR");
  TString tutdir = dir + "/macro/torinoDetector";
  
  TString tut_geomdir = tutdir + "/geometry";
  gSystem -> Setenv("GEOMPATH", tut_geomdir.Data());
  
  TString tut_configdir = tutdir + "/gconfig";
  gSystem -> Setenv("CONFIG_DIR", tut_configdir.Data());
  
  TString rootLogonFile = tut_configdir + "/rootlogon.C";
  
  // Load basic libraries
  gROOT -> LoadMacro(rootLogonFile);
  rootlogon();
  gSystem -> Load("libPndTorinoDetector.so");
  
  // create Instance of Run manager class
  FairRunSim *fRun = new FairRunSim();
  
  
  // set the MC version used
  //--------------------------------------------------
  
  fRun -> SetName("TGeant3");
  // Choose the Geant Navigation System
  // fRun -> SetGeoModel("G3Native");
  
  fRun -> SetOutputFile("data/testrun.root");
  
  //----------------Magnetic field-------------------
  //======== Constant Field ====================
  PndConstField *fMagField = new PndConstField();
  fMagField -> SetField(0., 10., 0.); // values are in kG
  fMagField -> SetFieldRegion(-50, 50, -50, 50, 350, 450);// values are in cm (xmin, xmax, ymin, ymax, zmin, zmax)
  fRun      -> SetField(fMagField);
  //======== End for magnetic field setting =======
  
  //------------ Materials setting ------------------
  fRun -> SetMaterials("media.geo");
  
  //-------------------------------------------------
  // Create and add detectors
  FairModule *Cave = new PndCave("CAVE");
  Cave -> SetGeometryFileName("cave.geo");
  fRun -> AddModule(Cave);
  
  FairModule *Magnet = new PndMagnet("MAGNET");
  Magnet -> SetGeometryFileName("magnet.geo");
  fRun -> AddModule(Magnet);
  
  FairDetector *Torino = new PndTorinoDetector("TORINO", kTRUE);
  Torino -> SetGeometryFileName("torino.geo");
  fRun -> AddModule(Torino);
  
  //---------------------------------------------------
  // Create and set Event Generator
  
  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  fRun -> SetGenerator(primGen);
  
  // Box Generator
  FairBoxGenerator* boxGen = new FairBoxGenerator(2212, 3); //2212 = proton 3 = multipl.
  boxGen -> SetPRange(2., 2.); // GeV/c setPrange vs setPtRange
  boxGen -> SetPhiRange(0, 360); // Azimuth angle range [degree]
  boxGen -> SetThetaRange(3, 10); //Polar angle in lab system range [degree]
  boxGen -> SetCosTheta(); //uniform generation on all the solid angle(default)
  
  // boxGen -> SetXYZ(0., 0.37, 0.);
  primGen -> AddGenerator(boxGen);
  
  fRun -> SetStoreTraj(kTRUE);
  
  fRun -> Init();
  
  // -Trajectories Visualization(TGeoManager Only)
  //-----------------------------------------------------
  
  // Set cuts for storing the trajectories
  /* 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);
   */
  
  // Fill the Parameter containers for this run
  //-----------------------------------------------------
  
  FairRuntimeDb *rtdb = fRun -> GetRuntimeDb();
  Bool_t kParamterMerged = kTRUE;
  FairParRootFileIo* output = new FairParRootFileIo(kParamterMerged);
  output -> open("data/testparams.root");
  rtdb -> setOutput(output);
  
  rtdb -> saveOutput();
  rtdb -> print();
  
  // Transport nEvents
  //---------------------------------------------------------
  
  // Int_t nEvents = 1;
  fRun -> Run(nEvents);
  
  fRun -> CreateGeometryFile("data/geofile_full.root");
  
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
  
}