#include <list>
#include <iostream>
#include "../Hits.h"
#include <cassert>
#include"TApplication.h"

#include <TROOT.h>
#include <TObject.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1I.h>
#include "../../../src/TCtrack.h"
#include "../../../src/TCevent.h"
#include <TCanvas.h>
#include <TGraph.h>
#include "helpers.h"
#include <TStyle.h>
#include <TF1.h>
#include <Math/Polynomial.h>
#include <TH2D.h>
#include <TProfile.h>
int main(int argc,char **argv){
  assert(argc==3);

  TApplication theApp("theApp",NULL,NULL);

  TFile::Open(argv[1]);
  TTree *tree = (TTree*)gROOT->FindObject("Hits");
  std::string alignmentFile = argv[2];
  CsGEMPlane* GM01X1=new CsGEMPlane();
  CsGEMPlane* GM01Y1=new CsGEMPlane();
  CsGEMPlane* GM02X1=new CsGEMPlane();
  CsGEMPlane* GM02Y1=new CsGEMPlane();

  CsGEMPlane* SI01X1=new CsGEMPlane();
  CsGEMPlane* SI01Y1=new CsGEMPlane();
  CsGEMPlane* SI02X1=new CsGEMPlane();
  CsGEMPlane* SI02Y1=new CsGEMPlane();
  
  std::list<CsGEMCluster*> GM01X1___fClusters;
  std::list<CsGEMHit*> GM01X1___fHits;

  TBranch *branchGM01X1 =tree->GetBranch("GM01X1__.");
  TBranch *branchGM01Y1 =tree->GetBranch("GM01Y1__.");
  TBranch *branchGM02X1 =tree->GetBranch("GM02X1__.");
  TBranch *branchGM02Y1 =tree->GetBranch("GM02Y1__.");

  TBranch *branchSI01X1 =tree->GetBranch("SI01X1__.");
  TBranch *branchSI01Y1 =tree->GetBranch("SI01Y1__.");
  TBranch *branchSI02X1 =tree->GetBranch("SI02X1__.");
  TBranch *branchSI02Y1 =tree->GetBranch("SI02Y1__.");
  
  
  branchGM01X1->SetAddress(&GM01X1);
  branchGM01Y1->SetAddress(&GM01Y1);
  branchGM02X1->SetAddress(&GM02X1);
  branchGM02Y1->SetAddress(&GM02Y1);
  
  branchSI01X1->SetAddress(&SI01X1);  
  branchSI01Y1->SetAddress(&SI01Y1);  
  branchSI02X1->SetAddress(&SI02X1);  
  branchSI02Y1->SetAddress(&SI02Y1);  

  TH1I *histogramSI1X = new TH1I("Multiplicity  SI01X", "Multiplicity  SI01X", 20, 0, 20);
  TH1I *histogramSI1Y = new TH1I("Multiplicity  SI01Y", "Multiplicity  SI01Y", 20, 0, 20);
  TH1I *histogramSI2X = new TH1I("Multiplicity  SI02X", "Multiplicity  SI02X", 20, 0, 20);
  TH1I *histogramSI2Y = new TH1I("Multiplicity  SI02Y", "Multiplicity  SI02Y", 20, 0, 20);

  TH1I *histogramGM1X = new TH1I("Multiplicity  GM01X", "Multiplicity  GM01X", 20, 0, 20);
  TH1I *histogramGM1Y = new TH1I("Multiplicity  GM01Y", "Multiplicity  GM01Y", 20, 0, 20);
  TH1I *histogramGM2X = new TH1I("Multiplicity  GM02X", "Multiplicity  GM02X", 20, 0, 20);
  TH1I *histogramGM2Y = new TH1I("Multiplicity  GM02Y", "Multiplicity  GM02Y", 20, 0, 20);

  TH1I *histogramEvent = new TH1I("Event multiplicity", "Event multiplicity", 40, 0, 40);
  TH1I *histogramStart = new TH1I("Start cluster multiplicity", "Start cluster multiplicity", 40, 0, 40);
  TH1I *histogramEnd = new TH1I("End cluster multiplicity", "End cluster multiplicity", 40, 0, 40);
  TH1I *histogramTrack = new TH1I("Track multiplicity", "Track multiplicity", 40, 0, 40);
  TH1D *histogramChi2 = new TH1D("Chi2", "Chi2", 100000, 0, 100000);


  TH1D *histogramSI1Xhitpoint = new TH1D("SI1Xhitpoint", "SI1Xhitpoint", 200, 5, 25);
  TH1D *histogramSI1Yhitpoint = new TH1D("SI1Yhitpoint", "SI1Yhitpoint", 200, 15, 30);
  TH1D *histogramSI2Xhitpoint = new TH1D("SI2Xhitpoint", "SI2Xhitpoint", 200, 5, 25);
  TH1D *histogramSI2Yhitpoint = new TH1D("SI2Yhitpoint", "SI2Yhitpoint", 200, 15, 30);
  
  TH1D *histogramGM1Xhitpoint = new TH1D("GM1Xhitpoint", "GM1Xhitpoint", 200, 5, 25);
  TH1D *histogramGM1Yhitpoint = new TH1D("GM1Yhitpoint", "GM1Yhitpoint", 200, 15, 30);
  TH1D *histogramGM2Xhitpoint = new TH1D("GM2Xhitpoint", "GM2Xhitpoint", 200, 5, 25);
  TH1D *histogramGM2Yhitpoint = new TH1D("GM2Yhitpoint", "GM2Yhitpoint", 200, 15, 30);


  TH1D *histogramSI1XhitStrip = new TH1D("SI1Xhitstrip", "SI1Xhitstrip", 384, 0, 384);
  TH1D *histogramSI1YhitStrip = new TH1D("SI1Yhitstrip", "SI1Yhitstrip", 384, 0, 384);
  TH1D *histogramSI2XhitStrip = new TH1D("SI2Xhitstrip", "SI2Xhitstrip", 384, 0, 384);
  TH1D *histogramSI2YhitStrip = new TH1D("SI2Yhitstrip", "SI2Yhitstrip", 384, 0, 384);
  
  TH1D *histogramGM1XhitStrip = new TH1D("GM1Xhitstrip", "GM1Xhitstrip", 256, 0, 256);
  TH1D *histogramGM1YhitStrip = new TH1D("GM1Yhitstrip", "GM1Yhitstrip", 256, 0, 256);
  TH1D *histogramGM2XhitStrip = new TH1D("GM2Xhitstrip", "GM2Xhitstrip", 256, 0, 256);
  TH1D *histogramGM2YhitStrip = new TH1D("GM2Yhitstrip", "GM2Yhitstrip", 256, 0, 256);



  TH1D *histogramSI1XclSize = new TH1D("SI1XclSize", "SI1XclSize", 10, 0, 10);
  TH1D *histogramSI1YclSize = new TH1D("SI1YclSize", "SI1YclSize", 10, 0, 10);
  TH1D *histogramSI2XclSize = new TH1D("SI2XclSize", "SI2XclSize", 10, 0, 10);
  TH1D *histogramSI2YclSize = new TH1D("SI2YclSize", "SI2YclSize", 10, 0, 10);
  
  TH1D *histogramGM1XclSize = new TH1D("GM1XclSize", "GM1XclSize", 10, 0, 10);
  TH1D *histogramGM1YclSize = new TH1D("GM1YclSize", "GM1YclSize", 10, 0, 10);
  TH1D *histogramGM2XclSize = new TH1D("GM2XclSize", "GM2XclSize", 10, 0, 10);
  TH1D *histogramGM2YclSize = new TH1D("GM2YclSize", "GM2YclSize", 10, 0, 10);

  
  TH1D *histogramSI1Xresi = new TH1D("resSI1X", "resSI1X", 200, -3, 3);
  TH1D *histogramSI1Yresi = new TH1D("resSI1Y", "resSI1Y", 200, -3, 3);
  TH1D *histogramSI2Xresi = new TH1D("resSI2X", "resSI2X", 200, -3, 3);
  TH1D *histogramSI2Yresi = new TH1D("resSI2Y", "resSI2Y", 200, -3, 3);
  
  TH1D *histogramGM1Xresi = new TH1D("resGM1X", "resGM1X", 200, -3, 3);
  TH1D *histogramGM1Yresi = new TH1D("resGM1Y", "resGM1Y", 200, -3, 3);
  TH1D *histogramGM2Xresi = new TH1D("resGM2X", "resGM2X", 200, -3, 3);
  TH1D *histogramGM2Yresi = new TH1D("resGM2Y", "resGM2Y", 200, -3, 3);

  TProfile *histogramSI1XresiVu = new TProfile("resVuSI1X", "resVuSI1X", 200, 0,2,-10,10);
  TProfile *histogramSI1YresiVu = new TProfile("resVuSI1Y", "resVuSI1Y", 200, 0,2,-10,10);
  TProfile *histogramSI2XresiVu = new TProfile("resVuSI2X", "resVuSI2X", 200, 0,2,-10,10);
  TProfile *histogramSI2YresiVu = new TProfile("resVuSI2Y", "resVuSI2Y", 200, 0,2,-10,10);
  
  TProfile *histogramGM1XresiVu = new TProfile("resVuGM1X", "resVuGM1X", 1000, 0,10,-10,10);
  TProfile *histogramGM1YresiVu = new TProfile("resVuGM1Y", "resVuGM1Y", 1000, 0,10,-10,10);
  TProfile *histogramGM2XresiVu = new TProfile("resVuGM2X", "resVuGM2X", 1000, 0,10,-10,10);
  TProfile *histogramGM2YresiVu = new TProfile("resVuGM2Y", "resVuGM2Y", 1000, 0,10,-10,10);



  int nEvents=tree->GetEntries();
  /*
  looking at plots 
  cutting at 
  x<1.1 and y<1.1 for GM01X1 and GM01Y1
  x<0.8 and y<0.4 for GM02
  si still unknown
  */  
  using namespace std;
  TCevent* event = new TCevent();

  
  TCalign* a = TCalign::getInstance(alignmentFile);
  a->clear();
  a->read(alignmentFile);
  int n_tracks = 0;
  int track_tries=0;
  bool skipEvent=false;

  for(int i_ev=0;i_ev<nEvents;i_ev++) {
   
    if(i_ev%100==0){
      std::cout<<i_ev<<" "<<track_tries<<" "<<n_tracks<<"*******************************************************************************************************"<<std::endl;
    }
 

    tree->GetEntry(i_ev);
    std::vector<TCcluster> clSI1X;
    std::vector<TCcluster> clSI1Y;
    std::vector<TCcluster> clSI2X;
    std::vector<TCcluster> clSI2Y;
    std::vector<TCcluster> clGM1X;
    std::vector<TCcluster> clGM1Y;
    std::vector<TCcluster> clGM2X;
    std::vector<TCcluster> clGM2Y;

     event->clear();
     { 
       const std::list<CsGEMCluster*> clusterList = SI01X1->GetClusters();
       ampDiffCut(clusterList, clSI1X, 7,3,alignmentFile);
       histogramSI1X->Fill(clSI1X.size());
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramSI1XclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clSI1X.size();i++){
         histogramSI1XhitStrip->Fill(clSI1X.at(i).posUVW().x()/0.005);
       }
     }
     {
       const std::list<CsGEMCluster*> clusterList = SI01Y1->GetClusters();
       ampDiffCut(clusterList, clSI1Y, 7,4,alignmentFile);
       histogramSI1Y->Fill(clSI1Y.size()); 
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramSI1YclSize->Fill((*it)->GetSize());
       }       
       for(unsigned int i=0;i<clSI1Y.size();i++){
         histogramSI1YhitStrip->Fill(clSI1Y.at(i).posUVW().x()/0.005);
       }
     }
     { 
       const std::list<CsGEMCluster*> clusterList = SI02X1->GetClusters();
       ampDiffCut(clusterList, clSI2X, 7,5,alignmentFile);
       histogramSI2X->Fill(clSI2X.size());
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramSI2XclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clSI2X.size();i++){
         histogramSI2XhitStrip->Fill(clSI2X.at(i).posUVW().x()/0.005);
       }
     }
     {
       const std::list<CsGEMCluster*> clusterList =SI02Y1->GetClusters();
       ampDiffCut(clusterList, clSI2Y, 7,6,alignmentFile);
       histogramSI2Y->Fill(clSI2Y.size());    
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramSI2YclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clSI2Y.size();i++){
         histogramSI2YhitStrip->Fill(clSI2Y.at(i).posUVW().x()/0.005);
       }
     }
     { 
       const std::list<CsGEMCluster*> clusterList = GM01X1->GetClusters();
       ampRatioCut(clusterList, clGM1X, 1.1,1.2,1,alignmentFile);
       histogramGM1X->Fill(clGM1X.size());
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramGM1XclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clGM1X.size();i++){
         histogramGM2XhitStrip->Fill(clGM1X.at(i).posUVW().x()/0.04);
       }
     }
     {
       const std::list<CsGEMCluster*> clusterList = GM01Y1->GetClusters();
       ampRatioCut(clusterList, clGM1Y, 1.1,1.2,2,alignmentFile);
       histogramGM1Y->Fill(clGM1Y.size());    
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramGM1YclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clGM2Y.size();i++){
         histogramGM1YhitStrip->Fill(clGM1Y.at(i).posUVW().x()/0.04);
       }
     }
     { 
       const std::list<CsGEMCluster*> clusterList = GM02X1->GetClusters();
       ampRatioCut(clusterList, clGM2X,1,0.8,7,alignmentFile);
       histogramGM2X->Fill(clGM2X.size());
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramGM2XclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clGM2X.size();i++){
         histogramGM2XhitStrip->Fill(clGM2X.at(i).posUVW().x()/0.04);
       }
     }
     {
       const std::list<CsGEMCluster*> clusterList = GM02Y1->GetClusters();
       ampRatioCut(clusterList, clGM2Y, 1,0.8,8,alignmentFile);
       histogramGM2Y->Fill(clGM2Y.size());  
       for(std::list<CsGEMCluster*>::const_iterator it =clusterList.begin();it!=clusterList.end();it++ ){
         histogramGM2YclSize->Fill((*it)->GetSize());
       }
       for(unsigned int i=0;i<clGM2Y.size();i++){
         histogramGM2YhitStrip->Fill(clGM2Y.at(i).posUVW().x()/0.04);
       }
     }
     

     
         
     vector<TCcluster> startClusters;
     vector<TCcluster> endClusters;
     for(unsigned int i=0;i<clGM1X.size();++i){
       for(unsigned int j=0;j<clGM1Y.size();++j){
         double ratio=clGM1X.at(i).getAmp()/clGM1Y.at(j).getAmp();
         if(ratio>0.9&&ratio<1.1){
           TVector3 pos(clGM1X.at(i).posUVW().x(),clGM1Y.at(j).posUVW().x()-5,0);
           TVector3 err(clGM1X.at(i).getErr().x(),0.2,0.2);
           double amp(clGM1X.at(i).getAmp());
           int detID(clGM1X.at(i).getId()); 
           TCcluster _c(pos,err,amp,detID); 
           _c.setFit(true);

           TVector3 pos2(clGM1Y.at(j).posUVW().x(),clGM1X.at(i).posUVW().x()+5,0);
           TVector3 err2(clGM1Y.at(j).getErr().x(),0.2,0.2);;
           double amp2(clGM1Y.at(j).getAmp());
           int detID2(clGM1Y.at(j).getId()); 
           TCcluster _c2(pos2,err2,amp2,detID2); 
           _c2.setFit(true);

           startClusters.push_back(_c);
           startClusters.push_back(_c2);
           //cout<<"gm1 pulse cor"<<endl;
         }
       } //end looping over GM1Y clusters
     } //end looping over GM1X clusters
         
     
     for(unsigned int i=0;i<clGM2X.size();++i){
       for(unsigned int j=0;j<clGM2Y.size();++j){
         double ratio=clGM2X.at(i).getAmp()/clGM2Y.at(j).getAmp();
         if(ratio>0.9&&ratio<1.1){
           TVector3 pos(clGM2X.at(i).posUVW().x(),clGM2Y.at(j).posUVW().x()-5,0);
           TVector3 err(clGM2X.at(i).getErr().x(),0.2,0.2);
           double amp(clGM2X.at(i).getAmp());      
           int detID(clGM2X.at(i).getId()); 
           TCcluster _c(pos,err,amp,detID); 
           _c.setFit(true);

           TVector3 pos2(clGM2Y.at(j).posUVW().x(),clGM2X.at(i).posUVW().x()+5,0);
           TVector3 err2(clGM2Y.at(j).getErr().x(),0.2,0.2);;
           double amp2(clGM2Y.at(j).getAmp());
           int detID2(clGM2Y.at(j).getId()); 
           TCcluster _c2(pos2,err2,amp2,detID2); 
           _c2.setFit(true);

           endClusters.push_back(_c);
           endClusters.push_back(_c2);
           //cout<<"gm1 pulse cor"<<endl;
         }
       } //end looping over GM1Y clusters
     } //end looping over GM1X clusters
     

     
     if(startClusters.size()==0||endClusters.size()==0){
       continue;
     }

     int n_points = 0;

     for(unsigned int i=0;i<startClusters.size()/2;i++){
       int id =2*i;
       //cout<<"start x "<<startClusters.size()<<" "<<id<<endl;
       TVector3 spacePoint=startClusters.at(id).posXYZ();
       histogramGM1Xhitpoint->Fill(spacePoint.x());
       n_points++;
     }
     for(unsigned int i=0;i<clSI1X.size();++i){

       TVector3 spacePoint=clSI1X.at(i).posXYZ();
       histogramSI1Xhitpoint->Fill(spacePoint.x());
       n_points++;
     }
     for(unsigned int i=0;i<clSI2X.size();++i){
       TVector3 spacePoint=clSI2X.at(i).posXYZ();
       histogramSI2Xhitpoint->Fill(spacePoint.x());
       n_points++;
     }
     for(unsigned int i=0;i<endClusters.size()/2;i++){
       int id =2*i;
       //cout<<"end x "<<endClusters.size()<<" "<<id<<endl;
       TVector3 spacePoint=endClusters.at(id).posXYZ();
       histogramGM2Xhitpoint->Fill(spacePoint.x());
       n_points++;
     }

     n_points=0;
     
     for(unsigned int i=0;i<startClusters.size()/2;i++){
       int id =2*i+1;
       //       cout<<"start y "<<startClusters.size()<<" "<<id<<endl;
       TVector3 spacePoint=startClusters.at(id).posXYZ();
       histogramGM1Yhitpoint->Fill(spacePoint.y());
       n_points++;
     }
     for(unsigned int i=0;i<clSI1Y.size();++i){
       TVector3 spacePoint=clSI1Y.at(i).posXYZ();
       histogramSI1Yhitpoint->Fill(spacePoint.y());
       n_points++;
     }
     for(unsigned int i=0;i<clSI2Y.size();++i){
       TVector3 spacePoint=clSI2Y.at(i).posXYZ();
       histogramSI2Yhitpoint->Fill(spacePoint.y());
       n_points++;
     }
     for(unsigned int i=0;i<endClusters.size()/2;i++){
       int id =2*i+1;
       //       cout<<"end y "<<endClusters.size()<<" "<<id<<endl;
       TVector3 spacePoint=endClusters.at(id).posXYZ();
       histogramGM2Yhitpoint->Fill(spacePoint.y());
       n_points++;
     }
     
     histogramEvent->Fill(event->nClusters());

     



  }
  cout<<"Tried trackfits"<<track_tries<<endl;
  cout<<"n tracks "<<n_tracks<<endl;



  TFile* file = new TFile("clusterMultipCut.root","RECREATE");
  histogramGM1X->Write();
  histogramGM1Y->Write();
  histogramGM2X->Write();
  histogramGM2Y->Write();

  histogramSI1X->Write();
  histogramSI1Y->Write();
  histogramSI2X->Write(); 
  histogramSI2Y->Write();
     
  histogramEvent->Write();
  histogramStart->Write();
  histogramEnd->Write();
  histogramTrack->Write();
  histogramChi2->Write();

  histogramSI1XhitStrip->Write();
  histogramSI1YhitStrip->Write();
  histogramSI2XhitStrip->Write();
  histogramSI2YhitStrip->Write();

  histogramGM1XhitStrip->Write();
  histogramGM1YhitStrip->Write();
  histogramGM2XhitStrip->Write();
  histogramGM2YhitStrip->Write();


  histogramSI1Xresi->Write();
  histogramSI1Yresi->Write();
  histogramSI2Xresi->Write();
  histogramSI2Yresi->Write();
  
  histogramGM1Xresi->Write();
  histogramGM1Yresi->Write();
  histogramGM2Xresi->Write();
  histogramGM2Yresi->Write();

  histogramSI1Xhitpoint->Write();
  histogramSI1Yhitpoint->Write();
  histogramSI2Xhitpoint->Write();
  histogramSI2Yhitpoint->Write();

  histogramGM1Xhitpoint->Write();
  histogramGM1Yhitpoint->Write();
  histogramGM2Xhitpoint->Write();
  histogramGM2Yhitpoint->Write();
  file->Close();
  delete histogramGM1X;
  
  delete  histogramGM1Y;
  delete  histogramGM2X;
  delete  histogramGM2Y;

  delete histogramSI1X;
  delete  histogramSI1Y;
  delete  histogramSI2X;
  delete  histogramSI2Y;

  delete GM01X1;
  delete GM01Y1;
  delete GM02X1;
  delete GM02Y1;
  delete SI01X1;
  delete SI01Y1;
  delete SI02X1;
  delete SI02Y1;

  
return 0;
}