{
  //
  // Script for CER cali root gen
  // Make sure $DB_DIR is set correctly !!!
  //
  // Version 1.0    Nov 28, 2003     R. Michaels
  //
  // Modifications: yqiang May 13 2004

  
  // 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
  gHaApps->Add( new THaIdealBeam("Beam", "Ideal beam"));
  gHaApps->Add( new THaUnRasteredBeam("urb","Unrastered beam"));

  // add detectors that are not in the default config
  HRSL->AddDetector( new THaCherenkov("a1", "Aerogel 1" ));
  HRSL->AddDetector( new THaCherenkov("a2", "Aerogel 2" ));
  //  HRSL->AddDetector( new THaScintillator("s0","Scintillator S0"));
  HRSL->AddDetector( new THaShower("prl1", "Preshower counter"));
  HRSL->AddDetector( new THaShower("prl2", "Shower counter"));

  HRSR->AddDetector( new THaCherenkov("cer", "Gas Cherenkov counter" ));
  //  HRSR->AddDetector( new THaShower("ps", "Preshower counter"));
  //  HRSR->AddDetector( new THaShower("sh", "Shower counter"));
  
  // Bob's diagnostic library for online debugging
  //  THaDecData* bobdec = new THaDecData("D","Misc. Decoder Data");
  //  gHaApps->Add( bobdec ) ;

  // Prepare the physics modules to use
  
  // Watch the electron kinematics
  Double_t mass_el   = 0.000511; // electron mass 
  Double_t mass_kaon = 0.493677; // kaon mass in GeV
  Double_t mass_pion = 0.139;    // pion mass in GeV
  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_hidr  = 0.93827;//
  Double_t mass_tg  = mass_hidr;  // chosen for hydrogen runs

  /*
  // now, coincidence time using only single scintillator detectors
  gHaPhysics->Add( new THaS2CoincTime("CTpi","Corrected Coincidence time from only S2","L", "R", mass_pion, mass_el));

  // Cheating for now, and assuming a perfect, unrastered electron beam
  gHaPhysics->Add( new THaReactionPoint("ReactPt_R","Reaction vertex for Right","R","Beam"));
  gHaPhysics->Add( new THaReactionPoint("ReactPt_L","Reaction vertex for Left","L","Beam"));
  gHaPhysics->Add( new THaGoldenTrack( "R.gold", "HRS-R Golden Track", "R" ));
  gHaPhysics->Add( new THaGoldenTrack( "L.gold", "HRS-L Golden Track", "L" ));

  // Correct for using an Extended target
  gHaPhysics->Add( new THaExtTarCor("ExTgtCor_R","Corrected for extended target, HRS-R","R","ReactPt_R"));
  gHaPhysics->Add( new THaExtTarCor("ExTgtCor_L","Corrected for extended target, HRS-L","L","ReactPt_L"));
  */

  // 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
  
  system("rm -f Afile.root");
  int nrun;
  int found = 0;
  const char** path = 0;
  char filename[200];
  FILE *fd;
  const char* paths[] = {  
    "/adaql2/data1",
    "/adaql2/data2",
    "/adaql2/data3",
    "/adaql2/data4",
    "/cache/mss/halla/e04012/raw",
    0
  };

  path=paths;
  while( found==0 ) {
    cout << "Run number?\n";
    cin >> nrun;
    if( nrun<0 ) break;

    int ped;
    cout << "Is this a pedastal run? (0=No, 1=Yes)";
    cin >> ped;
    
    while ( path && *path ) {
      if(!ped) sprintf(filename,"%s/e04012_%d.dat.0",*path,nrun);
      else sprintf(filename,"%s/e04012_ped_%d.dat",*path,nrun);
      fd = fopen(filename,"r");
      if (fd != NULL) {
	found = 1;
	cout << "Will analyze file "<<filename<<'\n';
	fclose(fd);
	break;
      }
      path++;
    }
    if ( (!path || !*path) && !found ) {
      cout << "File not found.  Try again, or -1 to quit.\n";
    }
  }
  if(nrun<0) return;
  /*  
  int detid;
  cout << "Is this for Cherenkov or for Shower calibration (1=CH, 2=SH)?";
  cin >> detid;
  */
  THaRun run(filename);
  
  // Define the analysis parameters
  
  analyzer->SetCompressionLevel(1);
  analyzer->SetOdefFile("cere_root.odef");
  if (!ped) analyzer->SetCutFile("cere_root.cuts");
  //  else analyzer->SetCutFile("cere_root_ped.cuts");
  else analyzer->SetCutFile("cere_root.cuts");

  // Set output file names

  analyzer->SetSummaryFile( Form("./summary/summary_%d.log",nrun) );

  char rootfile[100];
  if (!ped)  sprintf(rootfile,"./rootfiles/cere_%d.root",nrun);
  else sprintf(rootfile,"./rootfiles/ped_%d.root",nrun);
  //  if (!ped) sprintf(rootfile,"/adaql3/work1/e04012/onlroot/e94107_cere_%d.root",nrun);
  //  else sprintf(rootfile,"/adaql3/work1/e04012/onlroot/e04012_ped_%d.root",nrun);
  analyzer->SetOutFile(rootfile);
 
  // make a symbolic link to Afile.root 
  //  char command[100];
  //strcpy(command,"ln -s ");
  //strcat(command,rootfile);
  //strcat(command," Afile.root");
  //system(command);
  
  // Set number of events to analyze
  /*
  if (detid==2){
    run.SetLastEvent(-1);
    cout << "For Shower calibration, will process the whole run"<<endl;
  }
  else{
    run.SetLastEvent(50000);
    cout << "For Cherenkov calibration, will process 50K events"<<endl;
  }
  */
  run.SetLastEvent(500000);

  // Analyze selected run
  analyzer->Process(run);
  
  // Clean up in case we re-run this script; //delete analyzer; //urgh this closes the outfile!
  run->Close();  
  gHaApps->Delete();
  gHaPhysics->Delete();

  cout << "To re-run this script, you must quit. A fix is in the works.\n\n";

}