// Macro for running Hades  with Geant3  or Geant4 (M. Al-Turany)

{
  TStopwatch timer;
  timer.Start();
  //  gDebug=1;

   // Parameter file name
  TString parFile = "params.root";

  // Load basic libraries
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C"); 
  basiclibs();
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libMCStack");
  gSystem->Load("libGen");
//  gSystem->Load("libPassive");
  gSystem->Load("libHades");
 
  
  FairRunSim *fRun = new FairRunSim();
  
  // set the MC version used 
  fRun->SetName("TGeant3"); 
//  fRun->SetGeoModel("G3Native");
  fRun->SetOutputFile("hadsim.root");

  ////// Set Material file Name  /////////////////////////////////////
  
 // fRun->SetMaterials("media.geo");  
  fRun->SetMaterials("mediaAll.geo");  
  ///////////////////Create and add detectors/////////////////////////

  FairModule *Cave= new HadesCave("cave");
  Cave->SetGeometryFileName("cave_280207162003.geo");
  fRun->AddModule(Cave);

  FairModule *sec= new HadesSectors("sector");
  sec->SetGeometryFileName("sect_280207162003.geo");
  fRun->AddModule(sec);
  
 
  FairModule *Magnet= new HadesCoils("Coils");
  Magnet->SetGeometryFileName("coils_280207162003.geo");
  fRun->AddModule(Magnet);

  FairModule *frames= new HadesFrames("FRAMES");
  frames->SetGeometryFileName("frames_280207162003.geo");
  fRun->AddModule(frames);
 
 

  FairModule *rich= new HadesRich("Rich", kTRUE);
  rich->SetGeometryFileName("rich_280207162003.geo");
  fRun->AddModule(rich);

 FairModule *target= new HadesTarget("Target");
  target->SetGeometryFileName("target_280207162003.geo");
  fRun->AddModule(target);


  FairModule *Mdc= new HadesMdc("Mdc", kTRUE);
  Mdc->SetGeometryFileName("mdc_280207162003.geo");
  fRun->AddModule(Mdc);

  FairModule *tof= new HadesTof("Tof", kFALSE);
  tof->SetGeometryFileName("tof_280207162003.geo");
  fRun->AddModule(tof);

  FairModule *shower= new HadesShower("Shower", kTRUE);
  shower->SetGeometryFileName("shower_280207162003.geo");
  fRun->AddModule(shower);


///////////////////Create and Set Event Generator ///////////////////////
/********************   Define input  **********************************/

  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  fRun->SetGenerator(primGen);

/********************   Define input  **********************************/

// -- Using  Particle Generator

  FairParticleGenerator *fGen2= new FairParticleGenerator();  
  
  /*
  Double_t pp = (gRandom->Rndm()*500+100) * 1./1000; // [100,500 ] Mev/c  
  Double_t ctheta = cos ( 3.1415 * 1./2. ) ; 
  Double_t cphi   = cos ( 3.1415* 55./180 ); 
  Double_t sphi   = sin ( 3.1415* 55./180 );
  cout << " px " << pp *ctheta * cphi << " py " <<  pp * ctheta * sphi << "pz" <<  pp * ctheta << endl;
  
  fGen2->SetMomentum(pp *ctheta * cphi , pp * ctheta * sphi, pp * ctheta);
  fGen2->SetVertex(0,0,0);
  fGen2->SetPDGType(2212);
  
  
  primGen->AddGenerator(fGen2);

 */
  
  FairBoxGenerator* boxGen = new FairBoxGenerator(11, 1); // 13 = muon; 1 = multipl.
  //FairBoxGenerator* boxGen = new FairBoxGenerator(2212, 1); // 13 = muon; 1 = multipl.
    boxGen->SetPRange(.05,.7); // GeV/c
  //boxGen->SetPtRange(1.,1.1); // GeV/c
  boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
//  boxGen->SetThetaRange(35., 36.); // Polar angle in lab system range [degre
 boxGen->SetThetaRange(20., 85.); // Polar angle in lab system range [degree]e]
  boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
  boxGen->Init();
  primGen->AddGenerator(boxGen);

   


// setting a field map
    HadesField *fMagField= new HadesField();
    
// setting a constant field

//  CbmConstField *fMagField=new CbmConstField();
//  fMagField->SetFieldXYZ(0, 10 ,0 ); // values are in kG
// MinX=-75, MinY=-40,MinZ=-12 ,MaxX=75, MaxY=40 ,MaxZ=124 );  // values are in cm
//  fMagField->SetFieldRegions(-74, -39 ,-22 , 74, 39 , 160 );   

    fRun->SetField(fMagField);

 // -----   Fill parameter containers   ------------------------------------
/*  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  Bool_t kParameterMerged = kTRUE;
  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
  parOut->open(parFile.Data());



  rtdb->setOutput(parOut);
  rtdb->saveOutput();
  rtdb->print();

*/
   fRun->SetStoreTraj(kTRUE);
   
  // Init
  fRun->Init();

 // -Trajectories Visualization (TGeoManager Only )
 // -----------------------------------------------


 // Set cuts for storing the trajectpries
   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);
  

  // Propagate events
   fRun->Run(1000);

 timer.Stop();
 Double_t rtime = timer.RealTime();
 Double_t ctime = timer.CpuTime();
 
 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
 
}