x#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 *histogramChi2rough = new TH1D("Chi2", "Chi2", 10000, 0, 100000);
  TH1D *histogramChi2 = new TH1D("Chi2", "Chi2", 10000, 0, 1000);
  TH1D *histogramNDF = new TH1D("NDF", "NDF", 8, 0, 8);


  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 *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();

  TFile* rootOutfile = new TFile("tracks.root","RECREATE");
  cout<<"after rootfile creation"<<endl;
  TTree* outTree = new TTree("at2","testBench analysis tree");
  cout<<"after TTree creation"<<endl;
  outTree->Branch("event","TCevent",&event,32000,99);
  cout<<"after TTree setBranch"<<endl;
  
  TCalign* a = TCalign::getInstance(alignmentFile);
  a->clear();
  a->read(alignmentFile);

  int n_tracks = 0;
  int track_tries=0;

  /*
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  bool skipEvent=false;

  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);
  TCanvas * c = new TCanvas("occupancy","occupancy",10,10,800,600);
  c->Divide(2,1);
  */

  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;
    }
    double chiStore=1000000;
    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());
     }
     {
       const std::list<CsGEMCluster*> clusterList = SI01Y1->GetClusters();
       ampDiffCut(clusterList, clSI1Y, 7,4,alignmentFile);
       histogramSI1Y->Fill(clSI1Y.size());    
     }
     { 
       const std::list<CsGEMCluster*> clusterList = SI02X1->GetClusters();
       ampDiffCut(clusterList, clSI2X, 7,5,alignmentFile);
       histogramSI2X->Fill(clSI2X.size());
     }
     {
       const std::list<CsGEMCluster*> clusterList =SI02Y1->GetClusters();
       ampDiffCut(clusterList, clSI2Y, 7,6,alignmentFile);
       histogramSI2Y->Fill(clSI2Y.size());    
     }

     { 
       const std::list<CsGEMCluster*> clusterList = GM01X1->GetClusters();
       ampRatioCut(clusterList, clGM1X, 1.1,1.2,1,alignmentFile);
       histogramGM2X->Fill(clGM1X.size());
     }

     {
       const std::list<CsGEMCluster*> clusterList = GM01Y1->GetClusters();
       ampRatioCut(clusterList, clGM1Y, 1.1,1.2,2,alignmentFile);
       histogramGM2Y->Fill(clGM1Y.size());    
     }
     { 
       const std::list<CsGEMCluster*> clusterList = GM02X1->GetClusters();
       ampRatioCut(clusterList, clGM2X, 1,0.8,7,alignmentFile);
       histogramGM1X->Fill(clGM2X.size());
     }
     {
       const std::list<CsGEMCluster*> clusterList = GM02Y1->GetClusters();
       ampRatioCut(clusterList, clGM2Y, 1,0.8,8,alignmentFile);
       histogramGM1Y->Fill(clGM2Y.size());    
     }
     
     histogramEvent->Fill(event->nClusters());
     /*
     if(clSI1X.size()>3||
        clSI1Y.size()>3||
        clSI2X.size()>3||
        clSI2Y.size()>3){ 
       continue;
     }
     */

         
     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);
           if(pos.x()>4.5&&pos.x()<7.5&&pos2.x()>4.5&&pos2.x()<7.5){
             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||clSI1X.size()==0||clSI1Y.size()==0||clSI2X.size()==0||clSI2Y.size()==0){
       continue;
     }



     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());
     }
     for(unsigned int i=0;i<clSI1X.size();++i){

       TVector3 spacePoint=clSI1X.at(i).posXYZ();
       histogramSI1Xhitpoint->Fill(spacePoint.x());
     }
     for(unsigned int i=0;i<clSI2X.size();++i){
       TVector3 spacePoint=clSI2X.at(i).posXYZ();
       histogramSI2Xhitpoint->Fill(spacePoint.x());
     }
     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());
     }

     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());
     }
     for(unsigned int i=0;i<clSI1Y.size();++i){
       TVector3 spacePoint=clSI1Y.at(i).posXYZ();
       histogramSI1Yhitpoint->Fill(spacePoint.y());
     }
     for(unsigned int i=0;i<clSI2Y.size();++i){
       TVector3 spacePoint=clSI2Y.at(i).posXYZ();
       histogramSI2Yhitpoint->Fill(spacePoint.y());
     }
     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());
     }




     histogramStart->Fill(startClusters.size());
     histogramEnd->Fill(endClusters.size());
     
     //   cout<<"before loop startClusters.size "<<startClusters.size()<<endl;
     //looping over starting clusters 
     TCtrack* track_store=NULL;
     for(unsigned int i=0;i<startClusters.size()/2;++i){
       
       //looping over ending clusters---------------------------
       // for(unsigned int j=0;j<endClusters.size()/2;++j){

         //looping over SI1X clusters----------------------------
         for(unsigned int k=0;k<clSI1X.size();++k){
           //looping over SI1Y cluster---------------------------
           for(unsigned int l=0;l<clSI1Y.size();++l){
             //looping over SI2X clusters-------------------------
             for(unsigned int m=0;m<clSI2X.size();++m){
               //looping over SI2Y clusters--------------------------
               for(unsigned int n=0;n<clSI2Y.size();++n){
                 
                 track_tries++;
                 TCtrack* track = new TCtrack;
                 vector<TCcluster> trackCl;
                 trackCl.push_back(startClusters.at(2*i));
                 trackCl.push_back(startClusters.at(2*i+1));
                 trackCl.push_back(clSI1X.at(k));
                 trackCl.push_back(clSI1Y.at(l));
                 trackCl.push_back(clSI2X.at(m));
                 trackCl.push_back(clSI2Y.at(n));
                 //                 trackCl.push_back(endClusters.at(2*j));
                 //                 trackCl.push_back(endClusters.at(2*j+1));
                 track->addClusters(trackCl);

// #if 0
//                  TVector3 first;
//                  TVector3 last;
//                  first = startClusters.at(2*i).posXYZ();
//                  //                 last = endClusters.at(2*j+1).posXYZ();
//                  double startAX,startAY,startBX,startBY;
//                  TVector3 startDir = last-first;
//                  startAX=startDir.X()/startDir.Z();
//                  startAY=startDir.Y()/startDir.Z();
//                  double startLambda=-1.*first.Z()/startDir.Z();
//                  startBX=first.X()+startLambda*startDir.X();
//                  startBY=first.Y()+startLambda*startDir.Y();
         

//                  TVector3 first, last, dir;
//                  first = startClusters.at(2*i).posXYZ();
//                  //                 last = endClusters.at(2*j).posXYZ();
//                  dir = last-first;
//                  double startAX=dir.X()/dir.Z();
//                  double startLambda=-1.*first.Z()/dir.Z();
//                  double startBX=first.X()+startLambda*dir.X();

//                  first = startClusters.at(2*i+1).posXYZ();
//                  //                 last = endClusters.at(2*j+1).posXYZ();
//                  dir = last-first;
//                  double startAY=dir.Y()/dir.Z();
//                  startLambda=-1.*first.Z()/dir.Z();
//                  double startBY=first.Y()+startLambda*dir.Y();
                 


//                  Double_t vstart[4] = {startAX,startBX,startAY,startBY};
                 
//                  for(int a=0;a<track->nCl();a++){
                   
//                    histogramChi2->Fill(track->getCl(a).getChi2(vstart)/4);
//                    histogramChi2rough->Fill(track->getCl(a).getChi2(vstart)/4);
//                    TCcluster tmpcl = track->getCl(a);
//                    TVector3 resid=tmpcl.getResid(vstart);
//                    double uCl=tmpcl.posUVW().x();
//                    if(tmpcl.getFit()){
//                      switch(tmpcl.getId()){
//                      case 1:
//                        histogramGM1Xresi->Fill(resid.x());
//                        histogramGM1XresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 2:
//                        histogramGM1Yresi->Fill(resid.x());
//                        histogramGM1YresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 3:
//                        histogramSI1Xresi->Fill(resid.x());
//                        histogramSI1XresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 4:                         
//                        histogramSI1Yresi->Fill(resid.x());
//                        histogramSI1YresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 5:
//                        histogramSI2Xresi->Fill(resid.x());
//                        histogramSI2XresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 6:
//                        histogramSI2Yresi->Fill(resid.x());
//                        histogramSI2YresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 7:
//                        histogramGM2Xresi->Fill(resid.x());
//                        histogramGM2XresiVu->Fill(uCl,resid.x());
//                        break;
//                      case 8:
//                        histogramGM2Yresi->Fill(resid.x());
//                        histogramGM2YresiVu->Fill(uCl,resid.x());
//                        break;
//                      default:
//                       cout<<"unknown id "<<tmpcl.getId()<<endl;
//                      }
//                    }
//                  }
// #endif   
                 //cout<<"nclFit"<<track->nClFit()<<endl;
                 
                 if(track->fit(1,2,3,4,5,6,7,8)){
                   
                   //cout<<"track fit converged with chi2 "<<track->getChi2()/track->getNDF()<<endl;

                   if(track->getChi2()/track->getNDF() < chiStore){
                     track_store = track;
                     chiStore=track->getChi2()/track->getNDF();
                   }
                   else{
                     continue; 
                   }
                   
                   n_tracks++;
                   //cout<<"track nr. "<<n_tracks;
                   //                   cout<<"chi2/ndf "<<track->getChi2()/track->getNDF()<<endl;
                   //cout<<"ax "<<track->getAx()<<" bx "<<track->getBx()<<endl;
                   //cout<<"ay "<<track->getAy()<<" by "<<track->getBy()<<endl;
                   // double par[4]={track->getAX(),track->getBX(),track->getAY(),track->getBY()};

                
                   

 
                   
                 }
                                  //                 cout<<"bla"<<endl;
               }//end looping over SI2Y clusters
             }//end looping over SI2X clusters
           }//end looping over SI1Y clusters
         }//end looping over SI1X clusters
         // }//end looping over END clusters
     }//end looping over START clusters
     
     if(track_store!=NULL){
       histogramChi2->Fill(track_store->getChi2()/track_store->getNDF());
       histogramChi2rough->Fill(track_store->getChi2()/track_store->getNDF());
       histogramNDF->Fill(track_store->getNDF());
       for(unsigned int cl=0;cl<track_store->nCl();cl++){
         
         TCcluster tmpcl = track_store->getCl(cl);
         double uCl=tmpcl.posUVW().x();                   
         TVector3 resid=tmpcl.getRes();
         if(tmpcl.getFit()){
           switch(tmpcl.getId()){
           case 1:
             histogramGM1Xresi->Fill(resid.x());
             histogramGM1XresiVu->Fill(uCl,resid.x());
             break;
           case 2:
             histogramGM1Yresi->Fill(resid.x());
             histogramGM1YresiVu->Fill(uCl,resid.x());
             break;
           case 3:
             histogramSI1Xresi->Fill(resid.x());
             histogramSI1XresiVu->Fill(uCl,resid.x());
             break;
           case 4:                         
             histogramSI1Yresi->Fill(resid.x());
             histogramSI1YresiVu->Fill(uCl,resid.x());
             break;
           case 5:
             histogramSI2Xresi->Fill(resid.x());
             histogramSI2XresiVu->Fill(uCl,resid.x());
             break;
           case 6:
             histogramSI2Yresi->Fill(resid.x());
             histogramSI2YresiVu->Fill(uCl,resid.x());
             break;
           case 7:
             histogramGM2Xresi->Fill(resid.x());
             histogramGM2XresiVu->Fill(uCl,resid.x());
             break;
           case 8:
             histogramGM2Yresi->Fill(resid.x());
             histogramGM2YresiVu->Fill(uCl,resid.x());
             break;
           default:
             cout<<"unknown id "<<tmpcl.getId()<<endl;
           }
         }
                   
       }//end filling residual histograms
       event->addTrack(track_store);
     }else{
       continue;
     }
     
     /*
     If(event->nTracks()==0){
       continue;new 
     }
     */
     histogramTrack->Fill(event->nTracks());
     /*
     cout<<"press enter for next event"<<endl;
     string tmp;
     cin>>tmp;
     cout<<"next event coming up"<<endl;
     */    
  
     outTree->Fill();
     //if(i_ev>2000) break;
  }
  cout<<"Tried trackfits"<<track_tries<<endl;
  cout<<"n tracks "<<n_tracks<<endl;
  outTree->Write();
  rootOutfile->Close();

  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();
  histogramChi2rough->Write();
  histogramNDF->Write();

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

  histogramGM1Xhitpoint->Write();
  histogramGM1Yhitpoint->Write();
  histogramGM2Xhitpoint->Write();
  histogramGM2Yhitpoint->Write();


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


  histogramSI1XresiVu->Write();
  histogramSI1YresiVu->Write();
  histogramSI2XresiVu->Write();
  histogramSI2YresiVu->Write();
  
  histogramGM1XresiVu->Write();
  histogramGM1YresiVu->Write();
  histogramGM2XresiVu->Write();
  histogramGM2YresiVu->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;
}