{
//
// Script for only optics analysis
// Make sure $DB_DIR is set correctly !!!
//

// configuration files
const char* outdef = "./optics.odef";
const char* cutdef = "./optics.cuts";

// Set up the standard analyzer, if not already done
THaAnalyzer* analyzer = THaAnalyzer::GetInstance();
if( analyzer ) {
  analyzer->Close();
} else {
  analyzer = new THaAnalyzer;
}

// Get the run to analyze
int nrun;
bool found = false;
const char* path = 0;
while( !found ) {
  cout << "Run number?"<<endl;
  cin >> nrun;
  if( nrun<0 ) break;
  path = Form("./raw/happexsp_%d.dat.0",nrun);
  if( path ) {
    found = true;
    cout << "Will analyze file " << path << endl;
  }
  else
    cout << "File not found. Try again, or -1 to quit." << endl;
}
if(nrun<0) return;
THaRun run(path);

// Get number of events
int nmin = 1;
int nmax = -1;
cout << "Number of events to analyze? (-1=all, 0=quit)"<<endl;
cin >> nmax;
//nmax = 50000;
if(nmax == 0) return;

// Make sure we start with empty lists
gHaApps->Delete();
gHaPhysics->Delete();

// add the two spectrometers with the "standard" configuration
// (VDC planes, S1, and S2)
THaApparatus* HRSL = new THaHRS("L","Left arm HRS");
gHaApps->Add( HRSL );
THaApparatus* HRSR = new THaHRS("R","Right arm HRS");
gHaApps->Add( HRSR );

// beam
//THaApparatus* BEAM = new THaIdealBeam("Beam", "Ideal beam");
//gHaApps->Add( BEAM );
// Rastered beam
gHaApps->Add(new THaRasteredBeam("rb","Rastered Beamline"));
// Unrastered beam
THaApparatus* BEAM = new THaUnRasteredBeam("urb","Beamline");
gHaApps->Add( BEAM );

gSystem->Load("libDynCheck.so");
gHaApps->Add(new THaDynData("D","Misc. Decoder Data")) ;

// Prepare the physics modules to use
// Watch the electron kinematics
Double_t amu = 931.494*1.e-3;  // amu to GeV
Double_t mass_he3 = 2.809413; // in GeV
Double_t mass_C12 = 12.0*amu; // AMU to GeV
Double_t mass_H1 = 938.3*1.e-3; 
Double_t mass_O16 = 16.0*amu;
Double_t mass_tg  = mass_C12;  // chosen for Waterfall O16 elastic calibration runs

// Set up physics modules
gHaPhysics->Add( new THaReactionPoint("rpr","Reaction vertex for Right","R","urb"));
gHaPhysics->Add( new THaReactionPoint("rpl","Reaction vertex for Left", "L","urb"));
gHaPhysics->Add( new THaGoldenTrack( "R.gold", "HRS-R Golden Track", "R" ));
gHaPhysics->Add( new THaGoldenTrack( "L.gold", "HRS-L Golden Track", "L" ));
//gHaPhysics->Add( new THaExtTarCor("CorR","HRS-R xtg correction","R","rpr"));
//gHaPhysics->Add( new THaExtTarCor("CorL","HRS-L xtg correction","L","rpl"));

analyzer->SetCompressionLevel(1);
analyzer->SetOdefFile(outdef);
analyzer->SetCutFile(cutdef);
analyzer->EnableBenchmarks();

// Set output file names
analyzer->SetOutFile(Form("./rootfiles/optics_%d.root",nrun));
 
// Set number of events to analyze
run.SetFirstEvent(nmin);
run.SetLastEvent(nmax);

// Analyze selected run
analyzer->Process(run);

}