#include <TH1F.h>
#include <TH2.h>
#include <TF1.h>
#include <TObjArray.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <cmath>
#include <TApplication.h>
#include <TMath.h>
#include <TGraph2D.h>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <TCutG.h>
#include <TRoot.h>
#include <TFile.h>


#define MyClass_cxx

#include "MyClass2015_17.h"

using namespace std;

// Globals
TCanvas *c1 = 0;
TCanvas *c2 = 0;
TObject *obj=0;


TH2D *dwc0 = new TH2D("DWC0", "DWC0", 200, -0.05, 0.05, 200, -0.05, 0.05);
TH2D *dwc1 = new TH2D("DWC1", "DWC1", 100, -0.01, 0.01, 100, -0.01, 0.01);
TH2D *dwc2 = new TH2D("DWC2", "DWC2", 200, -0.05, 0.05, 200, -0.05, 0.05);
TH2D *Target_posn = new TH2D("Target_posn", "Target_posn", 100, -0.005, 0.005, 100, -0.005, 0.005);
TH2D *thetas = new TH2D("Theta_open", "ThetaPlot open", 100, -10, 10, 50, -5, 5);
TH2D *thetaB = new TH2D("Theta_cut", "ThetaPlot cuts", 100, -10, 10, 50, -5, 5);

TH1F *bgo = new TH1F("BGO0", "BGO1", 400, 0, 1.3);
TH1F *bgoa = new TH1F("BGO1a", "BGO1", 400, 0, 1.3);

TLegend *leg = new TLegend(0.7,0.7,0.95,0.95);
TGraph2D *cham = new TGraph2D();

float thetax, thetay;  // Incident angle phase space
float T_posx, T_posy;  // Extrapolation of trajectory onto target



void MyClass::Loop()
{
    if (fChain == 0) return;
    Long64_t nentries = fChain->GetEntriesFast();    
    Long64_t nbytes = 0, nb = 0;
 
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        
        thetax = ((dwc1x_abs - dwc0x_abs)/2.77)*1000;  // now in mrad !
        thetay = ((dwc1y_abs - dwc0y_abs)/2.77)*1000;  // now in mrad !

		// using the geometry of the experiment .... 
		// DWC0 <----  2.77-0.22-0.5 ---> Diamond <---- 0.22 + 0.05 ---> DWC1
        T_posx = 0.27*((dwc0x_abs - dwc1x_abs)/2.77)+ dwc1x_abs; // x = x_L * m_x + c_x, considering dwc1 at origin
        T_posy = 0.27*((dwc0y_abs - dwc1y_abs)/2.77)+ dwc1y_abs; // y = y_L * m_y + c_y, considering dwc1 at origin
		
		
		  if (dwc0valid == true  && dwc1valid == true 
			  && bgo1energy > 0.22 && bgo1energy < .8  
				  //&& abs(thetax)<.5 && abs(thetay)<.5 
					// && ((dwc0x_abs - (-0.00222))/ 0.02) * ((dwc0x_abs - (-0.00222))/0.02)  +
					// ((dwc0y_abs - (0))/0.007)  * ((dwc0y_abs - (0))/0.007) < 1 &&
					//&& ((dwc1x_abs - (-0.0))/ 0.002) * ((dwc1x_abs - (-0.00))/0.002)  +
					// ((dwc1y_abs - (0))/0.002)  * ((dwc1y_abs - (0))/0.002) < 1
              )
				  
		  {
					  
				
			dwc0->Fill(dwc0x_abs, dwc0y_abs);
            dwc1->Fill(dwc1x_abs, dwc1y_abs);
            Target_posn->Fill(T_posx, T_posy);
            //dwc2->Fill(dwc2x_abs, dwc2y_abs);
            bgo->Fill(bgo1energy, bgo1energy);
           
			//bgob->Fill(bgo1energy, bgo1energy); //POWAspectrum
			//thetas->Fill(trackThetaX, trackThetaY);
			thetas->Fill(thetax, thetay);
			
            
				}
        }
            
	
		c1->cd(1);
		bgo->SetLineColor(kGreen);
        bgo->Draw();
		
		c1->cd(2);
		thetas->Draw("COLZ");
			
		c1->cd(3);
		dwc0->Draw();

        c1->cd(4);
		dwc1->Draw("COLZ");
	
		c1->cd(5);
		Target_posn->Draw("COLZ");
		
    	//snipped; this code has the computer wait for you to draw the cut.
        cout << "\n ho ho ho - draw the cut ..... \n";
	    c2 = new TCanvas("c2","c2",100,100,1000,700); //canvas instantiated
		c2->Update();
		c2->cd();
		Target_posn->Draw("COLZ");
		c2->Update();
	    TCutG* cutg=(TCutG*)gPad->WaitPrimitive("CUTG","CutG"); // making cut, store to CUTG
		delete c2;

        /*
        cham->SetPoint(0,dwc0x_abs,0,dwc0y_abs);
        cham->SetPoint(222,dwc1x_abs,222,dwc1y_abs);
        cham->SetPoint(225,dwc2x_abs,225,dwc2y_abs);
         
		leg->AddEntry(bgob,"Incremented Reading (Original)");
        leg->AddEntry(bgoa,"Clear Signal");
        leg->AddEntry(bgo,"Noise Factor");
        leg->Draw("SAME");
        */
	

   
		
Target_posn->Reset();	
nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    
    thetax = ((dwc1x_abs - dwc0x_abs)/2.77)*1000;
    thetay = ((dwc1y_abs - dwc0y_abs)/2.77)*1000;

    T_posx = 0.27*((dwc0x_abs - dwc1x_abs)/2.77)+ dwc1x_abs; // x = x_L * m_x + c_x, considering dwc1 at origin
    T_posy = 0.27*((dwc0y_abs - dwc1y_abs)/2.77)+ dwc1y_abs; // y = y_L * m_y + c_y, considering dwc1 at origin
    
		  if (dwc0valid == true  && dwc1valid == true
			&& bgo1energy > 0.005 && bgo1energy < 0.6
            //&& abs(thetax)<0.5 && abs(thetay)<0.5
			&& ((dwc0x_abs - (-0.00222))/ 0.02) * ((dwc0x_abs - (-0.00222))/0.02)  +
			((dwc0y_abs - (0))/0.007)  * ((dwc0y_abs - (0))/0.007) < 1 
			//&& ((dwc1x_abs - (-0.0))/ 0.002) * ((dwc1x_abs - (-0.00))/0.002)  +
			//((dwc1y_abs - (0))/0.002)  * ((dwc1y_abs - (0))/0.002) < 1  
			// && ((thetax - (-4))/ 2) * ((thetax - (-4))/2)  +
			// ((thetay - (0))/1)  * ((thetay - (0))/1) < 1  
			&& cutg->IsInside(T_posx, T_posy)
              ){
              
              //dwc0->Fill(dwc0x_abs, dwc0y_abs);
              //dwc2->Fill(dwc2x_abs, dwc2y_abs);
			  Target_posn->Fill(T_posx, T_posy);
              bgoa->Fill(bgo1energy, bgo1energy); //PS
              thetaB->Fill(thetax, thetay);
              //thetaB->Fill(thetax, thetay);              
              
                    }  // end of the IF statement
              }  // end of the second for loop

     
	c1->cd(5);
    Target_posn->Draw("COLZ");

    c1->cd(1);
    bgoa->SetLineColor(kBlue);
    bgoa->Draw("SAME");
   
	c1->cd(6);
	thetaB->Draw("COLZ");
	c1->Update();


	
}  // end of loop function




int main(int argc, char **argv) {
    //Reset ROOT and connect tree file
    char filenameDef[]="hello.root";
    char* filename=filenameDef;
    if (argc>1) filename=argv[1];
    gROOT->Reset();
    TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
    if (!f) {
        f = new TFile(filename);
    }
    TTree *our_tree = (TTree*)gDirectory->Get("tree");
	
	if (!our_tree) {
        fprintf(stderr,"tree ntuple not found in file, exiting.\n");
        exit(1);
    }
    new TApplication("DAQanalisys",&argc,argv);
    c1 = new TCanvas("c1","c1",1200,800); //canvas instantiated
	c1->Divide(3,2);
    
    //instantiate an object of type MyClass call ev
    MyClass *ev = new MyClass(our_tree);
    
    ev->Loop();
    
    c1->Update();
    
    printf("\n\n\n Done  %s.  Now waiting...\n",filename);
    
    std::cout << "\n= Double click mouse button in graphics window to end program.\n";
    
	
    int n=1;
    while (n>0) {
        cout << n << ", ";
        n++;
//		obj= c1->WaitPrimitive();
		obj=(TObject*)gPad->WaitPrimitive();
        printf("Loop i=%d, found objIsA=%s, name=%s\n", n,obj->ClassName(),obj->GetName());
        if (obj==0) break;
        // Waits until you click the middle mouse button in the graphics window 
    }

	
	
    delete c1;
    delete f;
    delete ev;
	
    delete dwc0;
    delete dwc1;
    delete dwc2;
    delete bgo;
    delete bgoa;
    delete cham;
    delete leg;
    delete thetas;
    delete thetaB;
	
}