// #include "TFile.h"
// #include "TTree.h"
// #include "TH2F.h"
// #include "TObjArray.h"
// #include "TEllipse.h"
// #include "TCanvas.h"
// #include "TPad.h"
// #include "TMath.h"
// #include "TGraph.h"
// #include "TVector3.h"

// #include <iostream>
// using namespace std;


// // global: list of cones
// TObjArray* cones;
// //global last mouse click for drag-ang-drop
// Int_t downx, downy, upx,upy;

// //global histogram used for axis 
// TH2F *axis;

// enum bmActionType {none, stretching};
// bmActionType thisaction;

TPad* paddy;
TCanvas* buttoncanvas;  
TEllipse* te;
TObjArray* cones;


void hand_analyze3() {

  //	gROOT->SetStyle("Plain");

  // Open ROOT file
  cout << "Enter the filename you want to analyze: " << endl;
  TString infilename;
  cin >> infilename;
  
    
  TFile *f = new TFile(infilename);
  TString outfilename = infilename.Tokenize(".")->First()->GetName();
  outfilename+="_cones.root";
  cout << "Data will be written to " << outfilename << ", press return to accept or type new name " << endl;
  TString userinput = Getline("");
 
  if (userinput.Length() != 1) //user has asked for a different filename
    { outfilename = userinput; }

  outfilename.Remove(outfilename.Length()-1); //strip carriage return
  TTree* tree = (TTree*)f->Get("tree");
  int firstevent=0;

  if (tree==NULL) {cout << "The file you specified does not contain Cerenkov data.\n"; exit(); }
  
  cout << "This file has a tree with " << tree->GetEntries() <<  " events.\n";
  firstevent = atoi(Getline("Which event (0-N) do you want to start at?"));

  // Variables
  Int_t nHits, j, k, l, m, n, p, q;
  Double_t x,y,z, phi, r, s;
  char buffer[1000];
  Int_t cCnt = 0; // Cerenkov photon counter for histo. filling
  Double_t angle; // angle between betas
  Double_t dE; // Beta KE energy difference
  Double_t hc = 1241.172; // hc in eV*nm
  
  // Tree branch variables
  Int_t nHitsTracker1; // num hits on tracker 1 (cylindrical tracker)
  Int_t nHitsTracker2; // num hits on tracker 2 (top tracker plate)
  Int_t nHitsTracker3; // num hits on tracker 3 (bottom tracker plate)
  Int_t nHitsTracker4; // not used for cyl. geom.
  Int_t nHitsTracker5; // not used for cyl. geom.
  Int_t nHitsTracker6; // not used for cyl. geom.
  Double_t PhotonHitsX[100000]; // Photon X hit positions [cm]
  Double_t PhotonHitsY[100000]; // Photon Y hit positions [cm]
  Double_t PhotonHitsZ[100000]; // Photon Z hit positions [cm]
  Double_t PhotonTimes_ps[100000]; // Photon hit times [ps]
  Int_t TrackerNumber[100000]; // Tracker number hit by photon (1 = cylinder, 2 = top, 3 = bottom)
  Double_t TargetEdep_keV; // Edep in target [keV]
  TVector3* PrimaryVertexBeta1_cm = new TVector3(); // Beta 1 vertex vector [cm]
  TVector3* PrimaryVertexBeta2_cm = new TVector3(); // Beta 2 vertex vector [cm]
  TVector3* PrimaryMomentumBeta1 = new TVector3(); // Beta 1 momentum vector
  TVector3* PrimaryMomentumBeta2 = new TVector3(); // Beta 2 momentum vector
  Double_t EnergyBeta1_keV; // Beta 1 KE [keV]
  Double_t EnergyBeta2_keV; // Beta 2 KE [keV]
  Double_t PhotonEnergy_eV[100000]; // Photon energy [eV]
  Int_t HitType[100000]; // Photon hit type (0 = Cerenkov, 1 = Scintillation, -1 = Other)
  Int_t nHitsTotal; // Total photon hits in this event
	
  // Histograms and graphs
  TH1F *hAngle = new TH1F("hAngle","Angle between 0nbb b's",100,0,360);
  TH1F *hAsym = new TH1F("hAsym","Energy symmetry between 0nbb b's",100,-5,5);
  TH1F *hCerenkov = new TH1F("hCerenkov","Total number of Cerenkov photons produced",100,0,1000);
  TGraph *gAngle_dE = new TGraph();
  TGraph *event_graph; // graph of a single event
	
  // Set branch addresses
  tree->SetBranchAddress("nHitsTracker1",&nHitsTracker1);
  tree->SetBranchAddress("nHitsTracker2",&nHitsTracker2);
  tree->SetBranchAddress("nHitsTracker3",&nHitsTracker3);
  tree->SetBranchAddress("nHitsTracker4",&nHitsTracker4);
  tree->SetBranchAddress("nHitsTracker5",&nHitsTracker5);
  tree->SetBranchAddress("nHitsTracker6",&nHitsTracker6);
  tree->SetBranchAddress("nHitsTotal",&nHitsTotal);
  tree->SetBranchAddress("TargetEdep_keV",&TargetEdep_keV);
  tree->SetBranchAddress("PhotonHitsX",PhotonHitsX);
  tree->SetBranchAddress("PhotonHitsY",PhotonHitsY);
  tree->SetBranchAddress("PhotonHitsZ",PhotonHitsZ);
  tree->SetBranchAddress("PhotonTimes_ps",PhotonTimes_ps);
  tree->SetBranchAddress("PrimaryVertexBeta1_cm",&PrimaryVertexBeta1_cm);
  tree->SetBranchAddress("PrimaryVertexBeta2_cm",&PrimaryVertexBeta2_cm);
  tree->SetBranchAddress("PrimaryMomentumBeta1",&PrimaryMomentumBeta1);
  tree->SetBranchAddress("PrimaryMomentumBeta2",&PrimaryMomentumBeta2);
  tree->SetBranchAddress("EnergyBeta1_keV",&EnergyBeta1_keV);
  tree->SetBranchAddress("EnergyBeta2_keV",&EnergyBeta2_keV);
  tree->SetBranchAddress("PhotonEnergy_eV",PhotonEnergy_eV);
  tree->SetBranchAddress("HitType",HitType);
  tree->SetBranchAddress("TrackerNumber",TrackerNumber);
	
  /////////////////
  //create tree to save analyzed events
  //open (or create new) file for output
  TFile *fout = new TFile(outfilename,"recreate");
  TTree* ntcones;
  if (fout->Get("ntcones")==NULL)
    {
      ntcones = new TTree("ntcones","cerenkov cones list");
    }
  else
    {
      ntcones = (TTree*)fout->Get("ntcones");
    }

  const int nmaxellipse = 100;
  Double_t r1[nmaxellipse];
  Double_t r2[nmaxellipse];
  Double_t x1[nmaxellipse];
  Double_t y1[nmaxellipse];
  Int_t nphotons[nmaxellipse];
  Int_t ncones;
  Int_t eventno;
  ntcones->Branch("ncones",&ncones,"ncones/I");
  ntcones->Branch("eventno",&eventno,"eventno/I");
  ntcones->Branch("r1",r1,"r1[ncones]/D");
  ntcones->Branch("r2",r2,"r2[ncones]/D");
  ntcones->Branch("x1",x1,"x1[ncones]/D");
  ntcones->Branch("y1",y1,"y1[ncones]/D");
  ntcones->Branch("nphotons",nphotons,"nphotons[ncones]/I");
		

  // Set up canvas and pads
  for (Int_t i=0; i<tree->GetEntries(); i++) {

    tree->GetEntry(i);
    // calculate angle between betas
    angle = (180./TMath::Pi())*(PrimaryMomentumBeta1->Angle(*PrimaryMomentumBeta2));
    // beta KE energy difference
    dE = EnergyBeta1_keV - EnergyBeta2_keV;
    // count Cerenkov photons (same as nHitsTotal if no scintillation)
    for (Int_t j=0; j<nHitsTotal; j++) if (HitType[j]==0) cCnt++;
    // fill histograms
    hAngle->Fill(angle);
    gAngle_dE->SetPoint(i,dE,angle);
    hCerenkov->Fill(cCnt);
    cCnt = 0; // reset Cerenkov photon counter for next event
		
    // Uncomment to print out data from each event
    /*
      cout <<  "Nhits " <<  "fTarEdep " <<  "fn1 " <<  "fn2 "
      <<  "fn3 " <<  "fn4 " <<  "fn5 " <<  "fn6 " <<  "NCerPhot "
      <<  "B1energy " <<  "B2energy " <<  "SumEnergy" << endl;
      cout << nHitsTotal << " " <<  TargetEdep_keV << " "
      <<  nHitsTracker1 << " " <<  nHitsTracker2 << " "
      <<  nHitsTracker2 << " " <<  nHitsTracker4 << " "
      <<  nHitsTracker5 << " " <<  nHitsTracker6 << " "
      <<  cCnt << " " <<  EnergyBeta1_keV << " " <<  EnergyBeta2_keV
      << " " <<  EnergyBeta1_keV + EnergyBeta2_keV << endl << endl;
    */
  }
		
  // 	// separate canvas with a button on it
  // 	TCanvas* c4 = new TCanvas("next","next",100,100);
  // 	buttoncanvas = c4; 
  // 	TBox* tb = new TBox(0.1,0.1,0.9,0.9);
  // 	tb->SetFillColor(kRed);
  // 	tt = new TText();
  // 	tt->DrawText(0.2,0.2,"NEXT");
	
  // Plot single events in y-phi space (y is the central axis of the cylinder,
  // phi is the polar coordinate)
  TCanvas *c3 = new TCanvas();
  c3->AddExec("ex1",".x draw_circle.C"); 
  //c3->AddExec("ex1",".x test_exec.C"); 
	


  // Construct a standard axis for plotting
  TH2F* axis = new TH2F("axis","Cylindrical tracker",10,0,360,10,-30,30);
  axis->GetXaxis()->SetTitle("Angle [deg]");
  axis->GetYaxis()->SetTitle("Y [cm]");
  axis->SetStats(kFALSE);
  event_graph = new TGraph();
  // iterate over 200 events, plot, print to PS file
  for (Long_t evt=firstevent; evt<tree->GetEntries(); evt++) {
    eventno=evt;
    cout << "event number " << evt  << endl;
    c3->cd();
    event_graph->Set(0); // clear event graph
    k = 0;
    tree->GetEntry((Int_t) evt);
    for (Int_t j=0; j<nHitsTotal; j++) {
      if (TrackerNumber[j]==1) { // only grab photons that hit cylindrical tracker
	x = PhotonHitsX[j];
	y = PhotonHitsY[j];
	z = PhotonHitsZ[j];
	// unwrap {x,y,z} -> {y,phi}
	if (x>0 && z>0) phi = (180./TMath::Pi())*atan(x/z);
	if (x>0 && z<=0) phi = 90.+(180./TMath::Pi())*atan(TMath::Abs(z/x));
	if (x<=0 && z<=0) phi = 180.+(180./TMath::Pi())*atan(TMath::Abs(x/z));
	if (x<=0 && z>0) phi = 270.+(180./TMath::Pi())*atan(TMath::Abs(z/x));
	event_graph->SetPoint(k,phi,y); // populate a point on event graph
	k++;
      }
    }
    axis->Draw();
    gPad->SetGridx();
    gPad->SetGridy();
    event_graph->SetMarkerStyle(7);
    event_graph->SetMarkerColor(kBlue);
    event_graph->SetEditable(false);
    event_graph->Draw("PSAME");

    gPad->Update();
    paddy = gPad;// some sort of hack!  Need gpad in own pointer otherwise draw_circle can't find it.
		
		
    cones = new TObjArray();
		
		
    ////////////////
    //The following block is waiting for user keyboard input.
    //This keeps the pad alive while the user is clicking 
    //and activating the AddExec that we called earlier.
    //////////////

    TTimer  *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE);
    char    *input;
    Bool_t done = kFALSE;
		
    do {
      timer->TurnOn();
      timer->Reset();
      // Now let's read the input, we can use here any
      // stdio or iostream reading methods. like std::cin >> myinputl;
      input = Getline("<return> for next event, <q> to finish, <d> to delete last ellipse: "); 
      timer->TurnOff();
		  
      // Now useful stuff with the input!
      cout << input << endl;

      if (input[0]=='q')
	{
	  break;
	}
      if (input[0]=='d') 
	{
	  TEllipse* killthiste = cones->Last();
	  cones->RemoveLast();
	  delete killthiste;
	  continue;
	}

      if (input) done = kTRUE;

		  
    } while (!done);  //end of timer-wait-loop
		
	
		
    ////////////////
    //The following block runs when the user is done fiddling 
    //with one event screen and has typed "q"
    //////////////
    cout << "begin saving/analysis of ellipses" << endl;


    TIter nextring(cones);
    TEllipse* thisellipse;
    ncones = 0;
    while (thisellipse = (TEllipse*)nextring())
      {
	r1[ncones] = thisellipse->GetR1();
	r2[ncones] = thisellipse->GetR2();
	x1[ncones] = thisellipse->GetX1();	
	y1[ncones] = thisellipse->GetY1();	
	ncones++;
	if (ncones==nmaxellipse) continue;
      }
    ntcones->Fill();


    if (input[0]=='q') 
      {
	break;// if user requested to exit event loop
      }
		

  }
  ntcones->Write();
  fout->Close();
  c3->Close();
}