// Data Utility Examples

// The DSReader class makes it easy to get the run information and loop over all of the events in a rootfile.
void CountEvents(std::string filename){
  RAT::DU::DSReader dsreader(filename);

  RAT::DS::Run& run = dsreader.GetRun();
  for(int iEv=0; iEv<dsreader.GetEntryCount();iEv++)
    std::cout << iEv << std::endl;
}

// Once a file has been loaded, the DU classes contain the DB information. E.g. PMT information
void PrintPMTPositions(std::string filename){
  RAT::DU::DSReader dsreader(filename);
  // This static class member always exists. Just get a reference through utility.
  RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
  for(unsigned int iPMT = 0; iPMT < pmtInfo.GetCount(); iPMT++){
    TVector3 pmtPos = pmtInfo.GetPosition(iPMT);
    pmtPos.Print(); // ROOT TVector3 class has a built in print method
  }
}

// Data Structure Examples
//Print the date of the run, and the Nhits for each event
void PrintNhits(std::string filename){
  RAT::DU::DSReader dsreader(filename);
  RAT::DS::Run& run = dsreader.GetRun();
  std::cout << "Run Date:   " << run.GetDate();

  //Loop over entries
  for(int iEntry=0; iEntry<dsreader.GetEntryCount(); iEntry++){
    RAT::DS::Entry& ds = dsreader.GetEntry(iEntry);
    //Loop over triggered events in each entry
    for(int iEv=0; iEv<ds.GetEVCount(); iEv++){
      RAT::DS::EV& ev = ds.GetEV(iEv);
      std::cout << ev.GetNhits() << std::endl;
    }
  }
}
//A more complicated example: plot the QHS (integrated charge) on the PMTs in a set of events
TH1D* PlotQHS(std::string filename){
  RAT::DU::DSReader dsreader(filename);
  TH1D* qHist = new TH1D("","",10,0,10);
  // Loop over entries
  for(int iEntry=0; iEntry<dsreader.GetEntryCount(); iEntry++){
    RAT::DS::Entry& ds = dsreader.GetEntry(iEntry);
    // Loop over triggered events in each entry
    for(int iEv=0; iEv<ds.GetEVCount(); iEv++){
      RAT::DS::EV& ev = ds.GetEV(iEv);
      RAT::DS::CalPMTs& calPMTs = ev.GetCalPMTs(); //Get the set of calibrated hit PMTs
      //Loop over PMTs
      for(int iPMT=0; iPMT<calPMTs.GetCount(); iPMT++){
	RAT::DS::PMTCal& pmt = calPMTs.GetPMT(iPMT);
	qHist->Fill(pmt.GetQHS());
      } //CalPMTs
    } //Triggered events
  } //rootfile entries
  return qHist;
}

// Database Example: Direct Access. Print the refractive index of mineral oil as a function of wavelength

void PrintRefractiveIndex(std::string filename){

  // Database will be loaded automatically when ROOT files are opened
  RAT::DU::DSReader reader(filename);
  // Get Pointer to global database, like Utility there is just one instance that always exists
  RAT::DB *db = RAT::DB::Get();
  // Get a pointer to OPTICS table entry for the desired material (see rat/data/OPTICS.ratdb)
  // "OPTICS" refers to the table, "mineral_oil" is the entry name.
  RAT::DBLinkPtr optics = db->GetLink("OPTICS","mineral_oil");
  // Use this to get references to the data, one array for the wavelength and another for the RI
  // The arguments of GetDArray match the names in OPTICS.ratdb
  // Note that the Get methods are typed.
  const vector<double>& wavelengths     = optics->GetDArray("RINDEX_value1");
  const vector<double>& refractiveIndex = optics->GetDArray("RINDEX_value2");

  //Print out
  for(unsigned int i=0; i< wavelengths.size(); i++){
    std::cout << "Wavelength/nm: \t"     << wavelengths[i]     << "\t"
	      << "Refractive Index: \t"  << refractiveIndex[i]
	      << std::endl;
  }
}


// A larger script to test your understanding
void Test(std::string filename){
  // Open the file
  RAT::DU::DSReader dsreader(filename);

  // How many events are there
  int nEvents = dsreader.GetEntryCount();
  // Loop through 10 of them
  for(int i=0; i<10;++i){
    RAT::DS::Entry& ds = dsreader.GetEntry(i);

    // Get the MC information
    RAT::DS::MC& rmc = ds.GetMC();
    std::cout << "Event " << rmc.GetMCID() << std::endl;

    std::cout << rmc.GetMCParticleCount() << " particles, "
	      << rmc.GetMCPECount() << " photoelectrons, "
	      << rmc.GetScintEnergyDeposit() << " MeV deposited in scintillator "
	      << std::endl;
    // Now get some info on the first generated particle
    RAT::DS::MCParticle& rmcparticle = rmc.GetMCParticle(0);
    std::cout  << "PDG code = "<< rmcparticle.GetPDGCode()
	       << ", x-position = "<< rmcparticle.GetPosition().x()
	       << "mm, y-momentum ="<< rmcparticle.GetMomentum().y() << " MeV/c "
	       << std::endl;

    //Get the pmtInfo data utility
    RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();

    // Get the event info for any events triggered
    int nevC = ds.GetEVCount();
    std::cout << "Triggered " << nevC << " events " <<std::endl;
    if(nevC>0){
      RAT::DS::EV& rev = ds.GetEV(0);
      RAT::DS::CalPMTs& calpmts = rev.GetCalPMTs();
      std::cout << "#Normal Calibrated PMTs = "
		<< calpmts.GetCount()
		<< ", Nhits = "
		<< rev.GetNhits()
		<< std::endl;

      // Loop through the PMT hits and output charge and time and position
      for(int iPMT=0;iPMT<calpmts.GetCount();iPMT++){
	RAT::DS::PMTCal& pmt = calpmts.GetPMT(iPMT);

	//Output
	std::cout << "PMT " << pmt.GetID()
		  << " q = " << pmt.GetQHS()
		  << " t = " << pmt.GetTime()
		  << " x = " << pmtInfo.GetPosition(pmt.GetID()).x()
		  << " y = " << pmtInfo.GetPosition(pmt.GetID()).y()
		  << " z = " << pmtInfo.GetPosition(pmt.GetID()).z()
		  << std::endl;
      }
    }
  }
}