/**
* Solution_2.C
* @author J. Wilson <j.wilson2@physics.ox.ac.uk>
* Carry out the exercises in tutorial 2.
**/
void Solutions(std::string filename)
{
  // Read in file
  RAT::DU::DSReader dsreader(filename);

  //Get pmtInfo
  const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
  
  // Create the histograms (work with double precision throughout)
  // Radius of simulated events - expect them all inside the scintillator, but because of the neck this
  // goes up to R = 10000mm.
  TH1D *h_rsim = new TH1D("h_rsim","Simulate event radius",100,0,10000);
  h_rsim->GetXaxis()->SetTitle("R /mm");
  // radius of hit PMTs (should all be around 8400mm but with some spread)
  TH1D *h_rpmt = new TH1D("h_rpmt","Radius of hit PMTs",100,8000,9000);
  h_rpmt->GetXaxis()->SetTitle("R_{pmt}/mm)");
  // 2D histogram of Nhits vs simulated radius
  TH2D *h_nhits_r = new TH2D("h_nhits_r","Nhits versus simulated radius", 100,0,10000,100,0,1000);
  h_nhits_r->GetXaxis()->SetTitle("R/mm");
  h_nhits_r->GetYaxis()->SetTitle("Nhits");
  // raw PMT time
  TH1D *h_traw = new TH1D("h_traw","PMT raw time", 1000, -500,500);
  h_traw->GetXaxis()->SetTitle("raw time (ns)");
  // PMT time residual
  TH1D *h_tres = new TH1D("h_tres","PMT residual time", 500, -100,400);
  h_tres->GetXaxis()->SetTitle("residual time (ns)");
  // to calculate the time residual we need a couple of constants
  double c = 300; // speed of light mm/ns
  double RI = 1.505; // take the average RI for scintillator as an approx (cos obviously some of the path is in H2O and AV)
  double ftriggerdelay = 250.; // looked this up in DAQ.ratdb

  // Loop through events
  for(int i=0; i<dsreader.GetEntryCount();i++){
    RAT::DS::Entry& rds = dsreader.GetEntry(i);
    // Get the MC information
    RAT::DS::MC& rmc = rds.GetMC();
    // we know there is only one particle - the electron - get it
    RAT::DS::MCParticle& rmcparticle = rmc.GetMCParticle(0);
    // get the position
    TVector3 mcpos = rmcparticle->GetPosition();
    // The quick way to get the radius
    // (or you could do pow((mcpos.X()*mcpos.X()+mcpos.Y()*mcpos.Y()+mcpos.Z()*mcpos.Z()+),0.5))
    double mc_r = mcpos.Mag();
    // put this in the first histogram
    h_rsim->Fill(mc_r);

    // now access the events - we'll loop through them in case there is more than one triggered event per
    // simulated event
    int nevC = rds.GetEVCount();
    for(int iev=0;iev<nevC; iev++){
      // Get the Event information
      RAT::DS::EV& rev = rds.GetEV(iev);
      // We want the Nhits for the 2D histogram
      int nHits = rev.GetNhits();
      // now fill the 2D histogram
      h_nhits_r->Fill(mc_r,nHits);
      // Now we need to loop through all the PMTs
      RAT::DS::CalPMTs& calpmts = rev.GetCalPMTs();
      for(int ipmt=0;ipmt<calpmts.GetCount();ipmt++){
        TVector3 pmtpos = pmtInfo.GetPosition(calpmts.GetPMT(ipmt).GetID());
        double pmt_r = pmtpos.Mag();
        // put this in the pmt radius histogram
        h_rpmt->Fill(pmt_r);
        // To get at the time we need to access the PMT sample (assume only one sample per PMT)
        double timeRaw = (calpmts.GetPMT(ipmt)).GetTime();
        // put this in a histogram
        h_traw->Fill(timeRaw);
        // Lastly we want to calculate the time residual
        TVector3 path = pmtpos-mcpos;
        double calcTime = path.Mag()*RI/c;
        double timeRes = timeRaw - ftriggerdelay - calcTime;
        // and histogram this
        h_tres->Fill(timeRes);
      }
    }
  }
  
  // Done - now all we need to do is plot the histograms
  TCanvas *C1 = new TCanvas("C1");
  C1->Divide(2,2);
  C1->cd(1);
  h_rsim->Draw();
  C1->cd(2);
  h_rpmt->Draw();
  C1->cd(3);
  gPad->SetLogy();		// choose a log scale in y
  h_traw->Draw();
  C1->cd(4);
  gPad->SetLogy();		// choose a log scale in y
  h_tres->Draw();
  TCanvas *C2 = new TCanvas("C2");
  C2->cd();
  gStyle->SetPalette(1);
  h_nhits_r->Draw("COLZ1");
}