//  DragonSelector.cxx.
/// \class   DragonSelector DragonSelector.h.
/// This class is an example DragonSelector sort code for use with
/// ROOT.
/// The class definition in DragonSelector.h has been generated automatically
/// by the ROOT utility TTree::MakeSelector(). This class is derived
/// from the ROOT class TSelector. For more information on the TSelector
/// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
/// The following methods are defined in this file:
///    Begin():        called everytime a loop on the tree starts,
///                    a convenient place to create your histograms.
///    SlaveBegin():   called after Begin(), when on PROOF called only on the
///                    slave servers.
///    Process():      called for each event, in this function you decide what
///                    to read and fill your histograms.
///    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
///                    called only on the slave servers.
///    Terminate():    called at the end of the loop on the tree,
///                    a convenient place to draw/fit your histograms.
///
/// To use this file, try the following session on your Tree T:
///
/// Root > T->Process("DragonSelector.C")
/// Root > T->Process("DragonSelector.C","some options")
/// Root > T->Process("DragonSelector.C+")
/// \verbinclude example/README
/// \author  JEP.
/// \brief   DragonSelector example sorting class.
/// \version 1.1.
/// \date    Nov 2004.
///
#define DragonSelector_cxx

#include <sstream>
#include <functional>

#include <TH2.h>
#include <TStyle.h>
#include <TCutG.h>
#include <TBrowser.h>

#include "DragonSelector.h"
#include "TDataItem.h"
#include "TBgoDetector.h"
#include "TSbDetector.h"
#include "TIcDetector.h"
#include "TMcpDetector.h"

using std::stringstream;

void DragonSelector::Begin(TTree *)
{
  /// The Begin() function is called at the start of the query.
  /// When running with PROOF Begin() is only called on the client.
  /// The tree argument is deprecated (on PROOF 0 is passed).
  ///
  SetHistograms();
  SetCuts();
  SetDetectors();
  cout << "Beginning sort..." << endl;
  fProcessNum = 0;
  nCoincident = 0;
  nMCP = 0;
  nIC = 0;
  nGamma = 0;
}

void DragonSelector::SlaveBegin(TTree *tree)
{
  /// The SlaveBegin() function is called after the Begin() function.
  /// When running with PROOF SlaveBegin() is called on each slave server.
  /// The tree argument is deprecated (on PROOF 0 is passed).
  ///
  Init(tree);
  TString option = GetOption();
}

Bool_t DragonSelector::Process(Long64_t entry)
{
  /// The Process() function is called for each entry in the tree (or possibly
  /// keyed object in the case of PROOF) to be processed. The entry argument
  /// specifies which entry in the currently loaded tree is to be processed.
  /// It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
  /// to read either all or the required parts of the data. When processing
  /// keyed objects with PROOF, the object is already loaded and is available
  /// via the fObject pointer.
  ///
  /// This function should contain the "body" of the analysis. It can contain
  /// simple or elaborate selection criteria, run algorithms on the data
  /// of the event and typically fill histograms.
  ///
  /// WARNING when a selector is used with a TChain, you must use
  ///  the pointer to the current TTree to call GetEntry(entry).
  ///  The entry is always the local entry number in the current tree.
  ///  Assuming that fChain is the pointer to the TChain being processed,
  ///  use fChain->GetTree()->GetEntry(entry).
  ///
  fChain->GetTree()->GetEntry(entry);
  const vector<TSubEvent> &fGmemSubEvents = fDragonEvent->GetGmemSubEvents();
  const vector<TSubEvent> &fHmemSubEvents = fDragonEvent->GetHmemSubEvents();
  const vector<TSubEvent> &fCgmeSubEvents = fDragonEvent->GetCgmeSubEvents();
  const vector<TSubEvent> &fChmeSubEvents = fDragonEvent->GetChmeSubEvents();

  for(Int_t k=0; k<30;k++)
  { 
    SBGOEnergies[k] = 0;
    for(Int_t l=0; l<3;l++) {SBGOVector[k][l] = 0;}
  }
  // Gamma-ray singles.
  Int_t schan = 0;
  Float_t sval = 0;
  for (UInt_t i = 0; i < fGmemSubEvents.size(); i++)
  {
    SrfG = 0;
    //vector<TDataItem<UShort_t> > rf= fGmemSubEvents[i].GetTdcValues(62);
    if (fGmemSubEvents[i].GetTdcValues(62).size()) SrfG = fGmemSubEvents[i].GetTdcValues(62)[0].GetValue();
    for (UInt_t j = 0; j < fDetsG.size(); j++) 
    {
      fDetsG[j]->SetSignals(fGmemSubEvents[i]);
      fDetsG[j]->Calibrate();
      if (fDetsG[j]->IsHit())
      {
        if (dynamic_cast<TBgoDetector *>(fDetsG[j])) fDetsG[j]->PlotEnergies(hsBgoE);
        TH1F *temp2 = new TH1F("temp2","histogram to extract fucking channels",30,0,30);
        for(UInt_t k = 0; k<30;k++)
        {
          temp2->SetBinContent(temp2->FindBin(k),SchannelHist->GetBinContent(SchannelHist->FindBin(k)));
        }
        fDetsG[j]->PlotChannels(SchannelHist);
        schan = 0;
        for(UInt_t k = 0; k<30;k++)
        {
          if(temp2->GetBinContent(temp2->FindBin(k)) != SchannelHist->GetBinContent(SchannelHist->FindBin(k))) schan = k;
        }
        delete temp2;
        sval = 0;
        for (UInt_t k = 0; k < fDetsG[j]->GetEnergies().size(); k++)
        {
          sval = fDetsG[j]->GetEnergies()[k]->GetValue();
        }
        for(Int_t l=0; l<30;l++)
        {
          if(l == schan)
          {
            SBGOEnergies[l] = sval;
            for(Int_t m=0; m<3;m++)
            {
              SBGOVector[l][0] =  fDetsG[j]->GetCoordinates().X();
              SBGOVector[l][1] =  fDetsG[j]->GetCoordinates().Y();
              SBGOVector[l][2] =  fDetsG[j]->GetCoordinates().Z();
            }
          }
         } 
         GammaSinglesTree->Fill();
         for(Int_t k=0; k<30;k++)
         { 
           SBGOEnergies[k] = 0;
           for(Int_t l=0; l<3;l++) {SBGOVector[k][l] = 0;}
         }
      }
    }
    // This puts the detectors in highest->lowest order based on energy.
    sort(fDetsG.begin(), fBgoLast, GreaterDetectorE());
    if (dynamic_cast<TBgoDetector *>(fDetsG[0])) {
      fDetsG[0]->PlotEnergies(hsBgoE0);
      fDetsG[0]->PlotChannels(hsBgo0);
    }
    if (fDetsG[0]->IsHit() && fDetsG[1]->IsHit())
      hsBgoE0BgoE1->Fill(fDetsG[0]->GetEnergies()[0]->GetValue(),
                         fDetsG[1]->GetEnergies()[0]->GetValue());
    // We just plot the first rf value here, but any up to n is good.
    if (fGmemSubEvents[i].GetTdcValues(62).size())
      hsBgoRf0->Fill(fGmemSubEvents[i].GetTdcValues(62)[0].GetValue());
    for (UInt_t j = 0; j < fDetsG.size(); j++) fDetsG[j]->Clear();
    SrfG = 0;
  }

  // Heavy-ion singles.
  Float_t s2val = 0;
  Float_t s2chan = 0;
  for (UInt_t i = 0; i < fHmemSubEvents.size(); i++) 
  {
    SeventTof = 0;
    Stac = 0;
    SrfH = 0;
    TIcDetector *ic = 0;
    TMcpDetector *mcp  = 0;
    //vector<TDataItem<UShort_t> > rf = fHmemSubEvents[i].GetTdcValues(62);
    if (fHmemSubEvents[i].GetTdcValues(62).size()) SrfH = fHmemSubEvents[i].GetTdcValues(62)[0].GetValue();
    UShort_t tacMcp = fHmemSubEvents[i].GetAdcValue(34);
    for (UShort_t j = 0; j < fDetsH.size(); j++) 
    {
      fDetsH[j]->SetSignals(fHmemSubEvents[i]);
      fDetsH[j]->Calibrate();
      if (fDetsH[j]->IsHit()) 
      {
        if (dynamic_cast<TSbDetector *>(fDetsH[j])) 
        {
          fDetsH[j]->PlotEnergies(hsSbE);
          if(fScalerBinX< 60) fDetsH[j]->PlotEnergies(hsSbE2min);  
      	  if (fHmemSubEvents[i].GetTdcValues(62).size()) fDetsH[j]->PlotEvT(hsSbERf, 
              fHmemSubEvents[i].GetTdcValues(62)[0].GetValue());
        }
        else if (dynamic_cast<TMcpDetector *>(fDetsH[j])) 
        {
          mcp = dynamic_cast<TMcpDetector *>(fDetsH[j]);
          hsMcpT->Fill(tacMcp);
        }
        else if (dynamic_cast<TIcDetector *>(fDetsH[j])) 
        {
	  ic = dynamic_cast<TIcDetector *>(fDetsH[j]);
          SicPlate1E = 0;
          SicPlate2E = 0;
          ic->PlotAnodeSum(hsIcSum);
          // const vector< const TDataItem<Float_t > * > & s2fEnergies = fDetsH[j]->GetEnergies();
          s2val = 0;
          s2chan = 0;
          for (UInt_t l = 0; l < fDetsH[j]->GetEnergies().size(); l++)
          {
            s2val = fDetsH[j]->GetEnergies()[l]->GetValue();
            s2chan = fDetsH[j]->GetEnergies()[l]->GetChannel();
            if( s2chan == 0) SicPlate1E = s2val; 
            if( s2chan == 1) SicPlate2E = s2val; 
          }
          ic->PlotEdE(hsIcEdE, 1, 0);
          ic->PlotEnergies(hsIcE);
        }
        Stac = tacMcp;
        //cout << "Stac " << Stac << " Plate 0 " << SicPlate1E <<" Plate 1 " << SicPlate2E << endl;
        IcSinglesTree->Fill();
        SicPlate1E = 0;
        SicPlate2E = 0;
        Stac = 0;
      }
    }
    if (ic && mcp)
    { 
      ic->PlotEvT(hsIcEMcpT, tacMcp);
    }
    for (UShort_t k = 0; k < fDetsH.size(); k++) fDetsH[k]->Clear();
    SrfH = 0;
  }
  
  for(Int_t k=0; k<30;k++)
  { 
    BGOEnergies[k] = 0;
    for(Int_t l=0; l<3;l++) {BGOVector[k][l] = 0;}
  }
    
  // Gamma-ray and heavy-ion coincidence events.
  Int_t chan = 0;
  Float_t val = 0;
  if (fCgmeSubEvents.size() == fChmeSubEvents.size()) 
  { 
    nCoincident = nCoincident+1;
    for (UInt_t i = 0; i < fCgmeSubEvents.size(); i++) 
    { 
      nGamma = nGamma+1;
      //vector<TDataItem<UShort_t> > rfG = fCgmeSubEvents[i].GetTdcValues(62);
      //vector<TDataItem<UShort_t> > rfH = fChmeSubEvents[i].GetTdcValues(62);
      UShort_t tdcGh  = fCgmeSubEvents[i].GetTdcValue(63);
      UShort_t tdcHg  = fChmeSubEvents[i].GetTdcValue(63);
      UShort_t tacMcp = fChmeSubEvents[i].GetAdcValue(34);
      //Reset, re-load and then store the event tree values
      eventTof = 0;
      eventTof = tdcGh;
      icPlate1E = 0;
      icPlate2E = 0;
      icPlate1E = fChmeSubEvents[i].GetAdcValue(35);
      icPlate2E = fChmeSubEvents[i].GetAdcValue(36);
      tac = 0;
      // Sort the gamma-ray coincident subevents
      for (UInt_t j = 0; j < fDetsG.size(); j++) 
      {
        fDetsG[j]->SetSignals(fCgmeSubEvents[i]);
        fDetsG[j]->Calibrate();
        if (fDetsG[j]->IsHit())
        {
          if (dynamic_cast<TBgoDetector *>(fDetsG[j])) fDetsG[j]->PlotEnergies(hcBgoE);
          TH1F *temp = new TH1F("temp","histogram to extract fucking channels",30,0,30);
          for(UInt_t k = 0; k<30;k++) 
          {
          temp->SetBinContent(temp->FindBin(k),channelHist->GetBinContent(channelHist->FindBin(k)));
          }
          fDetsG[j]->PlotChannels(channelHist);
          chan = 0;
          for(UInt_t k = 0; k<30;k++) 
          {
          if(temp->GetBinContent(temp->FindBin(k)) != channelHist->GetBinContent(channelHist->FindBin(k))) chan = k;
          }
          delete temp;
          //const vector< const TDataItem<Float_t > * > & fEnergies = fDetsG[j]->GetEnergies(); 
          val = 0;
          for (UInt_t k = 0; k < fDetsG[j]->GetEnergies().size(); k++) 
          { 
	    val = fDetsG[j]->GetEnergies()[k]->GetValue();
          }
          for(Int_t l=0; l<30;l++)
          { 
            if(l == chan)
            { 
              BGOEnergies[l] = val;
              for(Int_t m=0; m<3;m++) 
              {
                BGOVector[l][0] =  fDetsG[j]->GetCoordinates().X();
                BGOVector[l][1] =  fDetsG[j]->GetCoordinates().Y();
                BGOVector[l][2] =  fDetsG[j]->GetCoordinates().Z();
              }
            }
          }
        }
      }  
      // This puts the BGO dets in highest->lowest order based on energy.
      sort(fDetsG.begin(), fBgoLast, GreaterDetectorE());
      if (dynamic_cast<TBgoDetector *>(fDetsG[0]))
      {
          fDetsG[0]->PlotEnergies(hcBgoE0);
          fDetsG[0]->PlotChannels(hcBgo0);
      }
      if (fDetsG[0]->IsHit() && fDetsG[1]->IsHit())
      {
      hcBgoE0BgoE1->Fill(fDetsG[0]->GetEnergies()[0]->GetValue(),fDetsG[1]->GetEnergies()[0]->GetValue());
      }
        
      // Sort the heavy-ion coincident subevents
      TIcDetector *ic = 0;
      TMcpDetector *mcp = 0;
      for (UInt_t j = 0; j < fDetsH.size(); j++) 
      {
        fDetsH[j]->SetSignals(fChmeSubEvents[i]);
        fDetsH[j]->Calibrate();
        if (fDetsH[j]->IsHit()) 
        {
          if (dynamic_cast<TMcpDetector *>(fDetsH[j])) 
          { 
            mcp = dynamic_cast<TMcpDetector *>(fDetsH[j]);
            nMCP = nMCP+1;
          }
          else if ((ic = dynamic_cast<TIcDetector *>(fDetsH[j]))) 
          {
            nIC = nIC +1;
            ic->PlotAnodeSum(hcIcSum);
            ic->PlotEnergies(hcIcE);
            ic->PlotEdE(hcIcEdE, 1, 0);
            if (tdcGh && tdcHg)
              hcGrTofHiTof->Fill(tdcGh, tdcHg);
            if (tdcGh < 8500 && tdcGh > 6500) ic->PlotEdE(hcIcEdEtof,1,0);
            //if (tdcGh < 11000 && tdcGh > 9500) ic->PlotEdE(hcIcEdEtof,1,0);
	  }          
          if (fDetsG[0]->IsHit())
	  {
          hcTargZ0->Fill(0.1 * fDetsG[0]->GetCoordinates().Z());
          }
        }
      }
    //cout <<"Before Filling channel0 value is " << BGOchannel0 << endl;
    if (ic && mcp) 
    {
      ic->PlotEvT(hcIcEMcpT, tacMcp);
      tac = fChmeSubEvents[i].GetAdcValue(34);
    }
    eventTree->Fill();
    for(Int_t k=0; k<30;k++)
    { 
      BGOEnergies[k] = 0;
      for(Int_t l=0; l<3;l++) {BGOVector[k][l] = 0;}
    }
    for (UInt_t j = 0; j < fDetsG.size(); j++) fDetsG[j]->Clear();
    for (UInt_t j = 0; j < fDetsH.size(); j++) fDetsH[j]->Clear();
    }
  }

  // Deal with the scaler event.
  const vector<TDataItem<UInt_t> > &s = fScalerEvent->GetScalers();
  for (UInt_t i = 0; i < s.size(); i++)
  {
  if(i != 23) hScaler[s[i].GetChannel()]->Fill(fScalerBinX, s[i].GetValue());
  if(i == 23 &&  fScalerBinX < 60 ) { hScaler[s[23].GetChannel()]->Fill(fScalerBinX,s[15].GetValue());}
  }
  if (s.size())
    {
    fScalerBinX++;
    }
 
  const vector<TDataItem<Float_t> > &epcs = fEpicsEvent->GetEpics();
  for (UInt_t i = 0; i < epcs.size(); i++)
  {
  hEpics[epcs[i].GetChannel()]->Fill(epcs[i].GetValue());
  }

return kTRUE;
}

void DragonSelector::SlaveTerminate()
{
  /// The SlaveTerminate() function is called after all entries or objects
  /// have been processed. When running with PROOF SlaveTerminate() is called
  /// on each slave server.
  ///
}

void DragonSelector::Terminate()
{
  /// The Terminate() function is the last function to be called during
  /// a query. It always runs on the client, it can be used to present
  /// the results graphically or save the results to file.
  ///
  cout << "Finishing sort..." << endl;
  // new TBrowser();
  cout << "number of coincident events is " << nCoincident << endl;
  cout << "number of coincident IC events is " << nIC << endl;
  cout << "number of coincident MCP events is " << nMCP << endl;
  cout << "number of coincident GAMMA events is " << nGamma << endl;
}

void DragonSelector::SetHistograms()
{
  /// Set up the histograms for the sort.  Add any folders or other
  /// histograms you want in the same fashion as below.
  ///
  cout << "Setting histograms..." << endl;
  
  eventTree = new TTree("eventTree","a simplistic coincidence event tree");
  eventTree->Branch("icPlate1E",&icPlate1E,"icPlate1E/F");  
  eventTree->Branch("icPlate2E",&icPlate2E,"icPlate2E/F");
  eventTree->Branch("eventTof",&eventTof,"eventTof/F");  
  eventTree->Branch("tac",&tac,"tac/F");  
  eventTree->Branch("BGOEnergies",BGOEnergies,"BGOEnergies[30]/F");
  eventTree->Branch("BGOVector",BGOVector,"BGOVector[30][3]/F");

  IcSinglesTree = new TTree("IcSinglesTree","a simplistic gamma singles event tree");
  IcSinglesTree->Branch("SicPlate1E",&SicPlate1E,"SicPlate1E/F");  
  IcSinglesTree->Branch("SicPlate2E",&SicPlate2E,"SicPlate2E/F");
  IcSinglesTree->Branch("SeventTof",&SeventTof,"SeventTof/F");  
  IcSinglesTree->Branch("Stac",&Stac,"Stac/F");  
  IcSinglesTree->Branch("SrfH",&SrfH,"SrfH/F");  
  
  GammaSinglesTree = new TTree("GammaSinglesTree","a simplistic gamma singles event tree");
  GammaSinglesTree->Branch("SBGOEnergies",SBGOEnergies,"SBGOEnergies[30]/F");
  GammaSinglesTree->Branch("SBGOVector",SBGOVector,"SBGOVector[30][3]/F");
  GammaSinglesTree->Branch("SrfG",&SrfG,"SrfG/F");  

  // Make some scaler names for the scaler histograms below.
  string names[scalerN];
  for (Int_t i = 0; i < scalerN; i++) {
    stringstream index;
    index << i;
    names[i] = "hScaler";
    names[i].append(index.str());
  }

  string enames[epicsN];
  for (Int_t i = 0; i < epicsN; i++) {
    stringstream index2;
    index2 << i;
    enames[i] = "hEpics";
    enames[i].append(index2.str());
  }


  // Titles for the scaler spectra.
  const Char_t *titles[scalerN] = {
    "BGO triggers presented",
    "BGO triggers acquired",
    "Disconnected G2",
    "Disconnected G3",
    "NaI triggers",
    "NaI triggers random",
    "Disconnected G6",
    "Beta-monitor",
    "cBeta-monitor",
    "Disconnected G9",
    "Disconnected G10",
    "Disconnected G11",
    "Tail triggers presented",
    "Tail triggers acquired",
    "SB0 triggers",
    "SB1 triggers",
    "DSSSD triggers",
    "MCP triggers",
    "IC0 triggers",
    "Disconnected H7",
    "Disconnected H8",
    "Disconnected H9",
    "Disconnected H10",
    "SB1 triggers for 2 min"
  };
  // Titles for the epics spectra. 
  string etitles[epicsN];
  for (Int_t i=0; i<epicsN;i++)
  {
    stringstream index3;
    index3 << i;
    etitles[i] = "Epics";
    etitles[i].append(index3.str());
  }

  // Tidy up before a re-run and setup Histogram folder.
  gDirectory->Delete("h*");
  TList *l = gROOT->GetListOfBrowsables();
  if (TFolder *f = dynamic_cast<TFolder *>
      (l->Remove(l->FindObject("DragonRoot"))))
    f->Delete();
  fHistFolder = new TFolder("DragonRoot", "Offline");
  l->Add(fHistFolder, "E1024");

  // Make folders to put the spectra in.
  TFolder *fSingle = fHistFolder->AddFolder("Singles", "Singles");
  TFolder *fsBgo = fSingle->AddFolder("Bgo", "BGO detector");
  TFolder *fsIc = fSingle->AddFolder("Ic", "IC detector");
  TFolder *fsMcp = fSingle->AddFolder("Mcp", "MCP detector");
  TFolder *fsSb = fSingle->AddFolder("Sb", "SB detector");
  TFolder *fCoincidence = fHistFolder->AddFolder("Coincidences",
                                                 "Coincidences");
  TFolder *fcBgo = fCoincidence->AddFolder("Bgo", "BGO detector");
  TFolder *fcIc = fCoincidence->AddFolder("Ic", "IC detector");
  TFolder *fAnalysis = fHistFolder->AddFolder("Analysis", "Analysis");
  TFolder *fScaler = fHistFolder->AddFolder("Scalers", "Scalers");
  TFolder *fScalerH = fScaler->AddFolder("HeavyIon", "HI scaler");
  TFolder *fScalerG = fScaler->AddFolder("GammaRay", "GR scaler");
  TFolder *fEpicsEvents = fScaler->AddFolder("EpicsEvents","Epics");  

  // Make some scaler spectra, put them in some browsable folders.
  for (Int_t i = 0; i < scalerN; i++) {
    hScaler[i] = new TH1F(names[i].c_str(), titles[i], 5000, 0, 5000);
    (i <= 12) ? fScalerG->Add(hScaler[i]) : fScalerH->Add(hScaler[i]);
  }
  // Make some epics spectra, put them in a browsable folder.
  for (Int_t i = 0; i< epicsN; i++){
  hEpics[i] = new TH1F(enames[i].c_str(), etitles[i].c_str(), 7000,0,7000);
  fEpicsEvents->Add(hEpics[i]);
  }

  // Now add detector folders and spectra.
  channelHist  = new TH1F("channelHist", "Wtf was jonty thinking",30,0,30);
  SchannelHist  = new TH1F("SchannelHist", "Wtf was jonty thinking",30,0,30);
  S2channelHist  = new TH2F("S2channelHist", "Wtf was jonty thinking 2",4096,0,4096,3,0,3);
  hsBgoE  = new TH2F("hsBgoE", "BGO chan v E", 1024, 0, 1024, 30, 0, 30);
  hsBgo0  = new TH1F("hsBgo0", "BGO_{0} chan", 30, 0, 30);
  hsBgoE0 = new TH2F("hsBgoE0", "BGO_{0} chan v E", 1024, 0, 1024, 30, 0, 30);
  hsBgoE0BgoE1 = new TH2F("hsBgoE0BgoE1", "BGO_{1} v BGO_{0}",
                          256, 0, 1024, 256, 0, 1024);
  hsBgoRf0 = new TH1F("hsBgoRf0", "BGO_{0} rf tdc", 512, 0, 512);
  hcBgoE  = new TH2F("hcBgoE", "BGO adc chan v val", 1024, 0, 1024, 30, 0, 30);
  hcBgo0  = new TH1F("hcBgo0", "BGO_{0} chan", 30, 0, 30);
  hcBgoE0 = new TH2F("hcBgoE0", "BGO_{0} adc chan v val",
                     1024, 0, 1024, 30, 0, 30);
  hcBgoE0BgoE1 = new TH2F("hcBgoE0BgoE1", "BGO_{1} v BGO_{0}",
                          256, 0, 1024, 256, 0, 1024);
  hsIcE   = new TH2F("hsIcE", "IC adc chan v val", 4096, 0, 4096, 3, 0, 3);
  hsIcEdE = new TH2F("hsIcEdE", "IC dE v E", 128, 0, 2048, 128, 0, 2048);
  hsIcSum = new TH1F("hsIcSum", "IC anode sum", 4096, 0, 4096);
  hsIcEMcpT = new TH2F("hsIcEMcpT", "IC0 E v TAC", 512, 0, 2048, 512, 0, 2048);
  hcIcE   = new TH2F("hcIcE", "IC adc chan v val", 4096, 0, 4096, 3, 0, 3);
  hcIcSum = new TH1F("hcIcSum", "IC anode sum", 4096, 0, 4096);
  hcIcEdE = new TH2F("hcIcEdE", "IC dE0 v dE1", 512, 0, 2048, 512, 0, 2048);
  hcIcEdEtof = new TH2F("hcIcEdEtof","IC dE0 v dE1",512,0,2048,512,0,2048);
  hcIcEMcpT = new TH2F("hcIcEMcpT", "IC0 E v TAC", 512, 0, 2048, 512, 0, 2048);
  hsMcpT = new TH1F("hsMcpT", "MCP TAC", 4096, 0, 4096);
  hsSbE   = new TH2F("hsSbE", "SB adc chan v val", 4096, 0, 4096, 2, 0, 2);
  hsSbE2min   = new TH2F("hsSbE2min", "SB1 adc chan v val, 2 min", 4096, 0, 4096, 2, 0, 2);
  hsSbERf = new TH2F("hsSbERf", "SB adc val v rf", 128, 128, 384, 512, 0, 4096);
  hcTargZ0 = new TH1F("hcTargZ0", "GR_{c0} target z-coordinate", 10, -20, 20);
  hcGrTofHiTof = new TH2F("hcGrTofHiTof", "GR_{HI} TOF v HI_{GR} TOF",
                          128, 0, 24576, 128, 0, 24576);
  fsBgo->Add(hsBgoE);
  fsBgo->Add(hsBgo0);
  fsBgo->Add(hsBgoE0);
  fsBgo->Add(hsBgoE0BgoE1);
  fsBgo->Add(hsBgoRf0);
  fcBgo->Add(hcBgoE);
  fcBgo->Add(hcBgo0);
  fcBgo->Add(hcBgoE0);
  fcBgo->Add(hcBgoE0BgoE1);
  fsIc->Add(hsIcE);
  fsIc->Add(hsIcEdE);
  fsIc->Add(hsIcSum);
  fsIc->Add(hsIcEMcpT);
  fcIc->Add(hcIcE);
  fcIc->Add(hcIcSum);
  fcIc->Add(hcIcEdE);
  fcIc->Add(hcIcEMcpT);
  fsMcp->Add(hsMcpT);
  fsSb->Add(hsSbE);
  fsSb->Add(hsSbE2min);
  fsSb->Add(hsSbERf);
  fAnalysis->Add(hcTargZ0);
  fAnalysis->Add(hcGrTofHiTof);
  fAnalysis->Add(hcIcEdEtof);
  fAnalysis->Add(eventTree);
  fAnalysis->Add(GammaSinglesTree);
  fAnalysis->Add(IcSinglesTree);
  fAnalysis->Add(channelHist);
  fAnalysis->Add(SchannelHist);
  fAnalysis->Add(S2channelHist);
}

void DragonSelector::SetDetectors()
{
  /// Set up detectors here.  This will change from experiment to
  /// experiment.  Just add the detectors that were used in you
  /// experiment by following the examples, make sure you include the
  /// correct header.
  ///
  cout << "Setting detectors..." << endl;

  // Use this database.
  const char *database = "mysql://isdaq04:3306/E1024";

  // Add some bgo detectors.
  const Int_t bgoDetsN = 30;
  const Int_t coords[bgoDetsN][3] = {{ 118, -48, -153}, {  0, -101, -122},
                                     { 118,  50, -122}, {  0,   99,  -92},
                                     {   0,  80,  -31}, {  0,   80,   31},
                                     {   0,  99,   92}, {  0, -101,  122},
                                     {   0,  50,  122}, {  0,  -48,  153},
                                     {-116, -26,  -92}, {118,  -26,  -92},
                                     {-116, -79,  -61}, {118,  -79,  -61},
                                     {-116,  27,  -61}, {118,   27,  -61},
                                     {-116, -26,  -31}, {118,  -26,  -31},
                                     {-116, -79,    0}, {118,  -79,    0},
                                     {-116,  27,    0}, {118,   27,    0},
                                     {-116, -26,   31}, {118,  -26,   31},
                                     {-116, -79,   61}, {118,  -79,   61},
                                     {-116,  27,   61}, {118,   27,   61},
                                     {-116, -26,   92}, {118,  -26,   92}};

  for (Int_t i = 0; i < bgoDetsN; i++) {
    stringstream index;
    index << "BGO" << i;
    string name = index.str();
    TBgoDetector *bgo = new TBgoDetector(name.c_str(), 1);
    bgo->SetDatabase(database);
    bgo->SetCalibration();
    bgo->SetXYZ(coords[i][0], coords[i][1], coords[i][2]);
    bgo->Print();
    fDetsG.push_back(bgo);
  }

  // Now the heavy-ion detectors.
  const Int_t sbDetsN = 2;
  for (Int_t i = 0; i < sbDetsN; i++) {
    stringstream index;
    index << "SB" << i;
    string name = index.str();
    TSbDetector *sb = new TSbDetector(name.c_str(), 1);
    sb->SetDatabase(database);
    sb->SetCalibration();
    sb->Print();
    fDetsH.push_back(sb);
  }
  TIcDetector *ic = new TIcDetector("IC0", 3);
  ic->SetDatabase(database);
  ic->SetCalibration();
  ic->Print();
  fDetsH.push_back(ic);
  TMcpDetector *mcp = new TMcpDetector("MCP0", 1);
  mcp->SetDatabase(database);
  mcp->SetCalibration();
  mcp->Print();
  fDetsH.push_back(mcp);

  // Partition the gamma-ray detector vector so that the BGO's are first.
  // If you don't understand what this does then leave it as is.
  fBgoLast = partition(fDetsG.begin(), fDetsG.end(),
                       EqualDetector<TBgoDetector>());
}

void DragonSelector::SetCuts()
{
  /// Set up any cuts that you want in here.
  ///
  cout << "Setting cuts..." << endl;

  // Set up a cut for the BGO hits.
  Float_t x[5] = {100, 200, 200, 100, 100};
  Float_t y[5] = {5, 5, 10, 10, 5};
  fBgoCut = new TCutG("bgoCut", 5, x, y);
}

/*
vector<TDataItem<UShort_t> >
DragonSelector::GetChannelHits(const vector<TDataItem<UShort_t> > &module,
                           const UShort_t channel) const
{
  /// Get the all hits from one channel in a vector.  This uses
  /// reverse iterators, as often hardware readout of each channel is
  /// in a largest to smallest fashion (LIFO-like).
  /// \param [in] module The hardware vector to be searched.
  /// \param [in] channel The channel to be searched for in the
  /// vector.
  /// \returns A vector containing hits with only channel channel.
  /// \note This is quite complex use of STL if you don't understand
  /// it then just leave as is.
  ///
  ptrdiff_t n = count_if(module.begin(), module.end(),
                         EqualChannel<UShort_t>(channel));
  vector<TDataItem<UShort_t> > items(n);
  remove_copy_if(module.rbegin(), module.rend(),
                 items.begin(), not1(EqualChannel<UShort_t>(channel)));
  return items;
}

UShort_t
DragonSelector::GetChannelHit(const vector<TDataItem<UShort_t> > &module,
                          const UShort_t channel) const
{
  /// Get a hit from a vector.  This uses reverse iterators, as often
  /// hardware readout of each channel is in a largest to smallest
  /// fashion (LIFO-like).
  /// \param [in] module The hardware vector to be searched.
  /// \param [in] channel The channel to be searched for in the
  /// vector.
  /// \returns The value of the hit with channel \a channel.
  /// \note This is quite complex use of STL if you don't understand
  /// it then just leave as is.
  ///
  UShort_t value = 0;
  vector<TDataItem<UShort_t> >::const_reverse_iterator item =
    find_if(module.rbegin(), module.rend(), EqualChannel<UShort_t>(channel));
  if (item != module.rend())
    value = item->GetValue();
  return value;
}
*/