#define QCAnal_cxx
#include "QCAnal.h"
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>

void QCAnal::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L QCAnal.C
//      Root > QCAnal t(runNumber, file_directory, file_prefix, sigType )
//      i.e. QCAnal t(0, "/tmp" , "SbtNtpTest", "MC" )
//      if you want to analyze MC file /tmp/SbtNtpTest_Run0.root
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   _nEntries = int(fChain->GetEntries ());


   std::cout << "Number of entries = " << _nEntries << std::endl;
   std::cout << "Number of entries = " << nentries << std::endl;

   Long64_t nbytes = 0, nb = 0;   

   for (Long64_t jentry=0; jentry<_nEntries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) { 
        std::cout << "jentry == " << jentry << "\n";
        std::cout << "ientry == " << ientry << "\n";
        break;
      }
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      if (Cut(ientry) < 0) continue;


      //reset cluster counters
      for (int i=0; i<_nDet; i++){
        _nClustU[i]=0;
        _nClustV[i]=0;
      }
      _nClustPxl=0;


      // fill beam profile histograms
      // using cluster informations

      for (int iClust= 0; iClust< ncluster ; iClust++){

	if (ClustSide[iClust]==0) _nClustU[ClustLayer[iClust]]++;
        if (ClustSide[iClust]==1) _nClustV[ClustLayer[iClust]]++;
        if (ClustSide[iClust]==-1) _nClustPxl++;


	FillBeamProfileHisto(ClustPos[iClust], ClustLayer[iClust], 
			     ClustSide[iClust]);
	
	// here there is a problem with number of cluster per layer
	// there is not this information in the ntuple need to extract it from the ntuple

	FillClusterHisto(ClustPH[iClust], ClustSize[iClust],
 			 ClustLayer[iClust], ClustSide[iClust] );


	if ( ClustPxlUPos[iClust] > -999.0 && //to be improved
	     ClustPxlVPos[iClust] > -999.0)	
	  _BeamProfPxl->Fill(ClustPxlUPos[iClust],ClustPxlVPos[iClust]);

      // fill beam profile layer-layer correlations
	if ( ClustPos[iClust] > -999.0 && //to be improved
	     ClustLayer[iClust] == 0) {
	  FillBeamLayerLayerCorrHisto(ClustPos[iClust], ClustSide[iClust]);
	}	
	

      } // end loop on clusters

      FillClusterMul();


      // loop on Digi

      for (int iDigi= 0; iDigi< ndigi ; iDigi++){

	FillOccHisto(DigiChan[iDigi], DigiRow[iDigi],DigiColumn[iDigi],
		     DigiLayer[iDigi],DigiSide[iDigi]);

	FillAdcHisto(DigiAdc[iDigi], DigiLayer[iDigi], DigiSide[iDigi]);

      } // end loop on Digi


      //fill trk informations
      _ntrk->Fill(ntrk);
      for (int iTrk= 0; iTrk< ntrk ; iTrk++){

	FillTrkHisto(chi2[iTrk], ndof[iTrk], 
		     fitStatus[iTrk], ItpX[iTrk], ItpY[iTrk],
		     SlpX[iTrk], SlpY[iTrk]);      
      }

      FillResidualHisto();

   } //loop on events
   
      // normalize occupancy histograms
   for (int iDet=0; iDet<_nDet; iDet++){
     Normalize1DHisto(_OccU[iDet]);
     Normalize1DHisto(_OccV[iDet]);
   }
   Normalize2DHisto(_OccPxl);
      

}

void QCAnal::FillBeamProfileHisto(double ClustPos, int ClustLayer, 
				  int ClustSide){

	for (int iDet= 0; iDet< _nDet ; iDet++){
	  //fill beam profile histo with cluster info
	  
	  //check for strip and striplets detectors, U side clusters
	  if ( ClustPos > -999.0 && //to be improved
	       ClustSide == 0 && ClustLayer == iDet) {
	    _BeamProfU[iDet]->Fill(ClustPos);
	  } 
	  //check for strip and striplets detectors, V side clusters
	  else if ( ClustPos > -999.0 && //to be improved
		    ClustSide == 1 && ClustLayer == iDet) {
	    _BeamProfV[iDet]->Fill(ClustPos);
	  } 
	} // loop on detectorElem

}

void QCAnal::FillBeamLayerLayerCorrHisto(double L0ClustPos, int L0ClustSide){

  for (int jClust= 0; jClust< ncluster ; jClust++){
    
      
      //check for strip and striplets detectors, U side clusters
    if ( L0ClustSide == 0 && ClustPos[jClust] > -999.0 && //to be improved
	   ClustSide[jClust] == 0  && ClustLayer[jClust] != 0) {
	_BeamLayerLayerCorrU[ClustLayer[jClust]-1]->Fill(L0ClustPos, ClustPos[jClust]);
	_DiffLayerCorrU[ClustLayer[jClust]-1]->Fill(ClustPos[jClust]-L0ClustPos);
      } 
      //check for strip and striplets detectors, V side clusters
      else if ( L0ClustSide == 1 && ClustPos[jClust] > -999.0 && 
		ClustSide[jClust] == 1 && ClustLayer[jClust] != 0 ) {
	_BeamLayerLayerCorrV[ClustLayer[jClust]-1]->Fill(L0ClustPos, ClustPos[jClust]);
	_DiffLayerCorrV[ClustLayer[jClust]-1]->Fill(ClustPos[jClust]-L0ClustPos);
      } 
      //check for pixel detectors
      else if ( L0ClustSide == 0 && ClustPxlUPos[jClust] > -999.0 && ClustLayer[jClust] != 0){
	_BeamLayerLayerCorrU[ClustLayer[jClust]-1]->Fill(L0ClustPos, ClustPxlUPos[jClust]);
      } else if ( L0ClustSide == 1 && ClustPxlVPos[jClust] > -999.0 && ClustLayer[jClust] != 0){
	_BeamLayerLayerCorrV[ClustLayer[jClust]-1]->Fill(L0ClustPos, ClustPxlVPos[jClust]);
      }

    
    
  }
}


void QCAnal::FillOccHisto(int DigiChan, int DigiRow, int DigiColumn, int DigiLayer, 
				  int DigiSide){

  if ( DigiChan > -1 && DigiSide == 0){ //strip detector check to be improved
       
    _OccU[DigiLayer]->Fill(DigiChan);
    
  }else if ( DigiChan > -1 && DigiSide == 1){ //strip detector check to be improved
       
    _OccV[DigiLayer]->Fill(DigiChan);
    
  } else if (DigiRow > -1 && DigiColumn >-1) {

    _OccPxl->Fill(DigiColumn, DigiRow);

  }

}


void QCAnal::FillAdcHisto(int DigiAdc,  int DigiLayer, int DigiSide){

  if ( DigiAdc > -1 && DigiSide == 0){ //strip detector check to be improved
       
    _AdcU[DigiLayer]->Fill(DigiAdc);
    
  }else if ( DigiAdc > -1 && DigiSide == 1){ //strip detector check to be improved
       
    _AdcV[DigiLayer]->Fill(DigiAdc);
    
  } 

}


void QCAnal::FillClusterHisto(int clustPH, int clustMul,
                              int clustLayer, int clustSide){

  if ( clustSide == 0 ) {
    _clustPhU[clustLayer]->Fill(clustPH);
    _clustMulU[clustLayer]->Fill(clustMul);
  }
  else if ( clustSide == 1 ) {
    _clustPhV[clustLayer]->Fill(clustPH);
    _clustMulV[clustLayer]->Fill(clustMul);
  }
  else if ( clustSide == -1 ){
    _clustMulPxl->Fill(clustMul);
  }

}

void QCAnal::FillClusterMul(){

   for (int iDet=0; iDet<_nDet; iDet++){
   _nclustU[iDet]->Fill(_nClustU[iDet]);
   _nclustV[iDet]->Fill(_nClustV[iDet]);
   }
  _nclustPxl->Fill(_nClustPxl);

 }

void QCAnal::FillTrkHisto(double Chi2, double Ndof, int FitStatus, 
			  double itpX, double itpY, double slopeX, double slopeY ){

   _trkChi2->Fill(Chi2);
   _trkProbChi2->Fill(TMath::Prob(Chi2,int(Ndof)));
   _trkSlpX->Fill(slopeX);
   _trkSlpY->Fill(slopeY);
   _trkItpX->Fill(itpX);
   _trkItpY->Fill(itpY);
   _trkFitStatus->Fill(FitStatus);
 }

void QCAnal::FillResidualHisto(){


  for  (int iInt=0; iInt < nIntersect; iInt++){

    for (int iDet=0; iDet < _nDet; iDet++){
      
      //      if ( intLayer[iInt] == iDet && intInside[iInt]==1 ){
      if ( intLayer[iInt] == iDet ){
	
 	if (intInside[iInt]==1){
 	  _intInside[iDet]->Fill(intXPos[iInt],intYPos[iInt]);
 	}
	/*
	  for (int iClust=0; iClust<ncluster; iClust++){
	  
          if (ClustLayer[iClust] == iDet) {
	  
	  if (ClustLayerType[iClust] != 2 && ClustSide[iClust] ==0 ){
	  _resU[iDet]->Fill(ClustPos[iClust]-intUPos[iInt]);
	    } 	else if (ClustLayerType[iClust] != 2 && ClustSide[iClust] ==1 ){
	    _resV[iDet]->Fill(ClustPos[iClust]-intVPos[iInt]);	   
	    }
	    else if (ClustLayerType[iClust] == 2){
            _resU[iDet]->Fill(ClustPxlUPos[iClust]-intUPos[iInt]);
            _resV[iDet]->Fill(ClustPxlVPos[iClust]-intVPos[iInt]);
	    }
	    
	    
	    }
	    
	    }
	*/


        for (int iSP=0; iSP<nSP; iSP++){
	  
	  if (spLayer[iSP] == iDet) {
	    
	    _resU[iDet]->Fill(spUPos[iSP]-intUPos[iInt]);
	    _resV[iDet]->Fill(spVPos[iSP]-intVPos[iInt]);	   
	    _resX[iDet]->Fill(spXPos[iSP]-intXPos[iInt]);
	    _resY[iDet]->Fill(spYPos[iSP]-intYPos[iInt]);	   
	    	    
          }

        }

      }

    }
  }

}



void QCAnal::InitBeamProfileHisto(){
  //assume 7 detectors:
  //1 map detector. 4 strip + 2 striplet
  //n. 0-6 strip+striplet detectors (n.4 pixel)

  for (int i=0; i< _nDet ; i++){

    if (4==i) continue;

    char HistoNameU[100];
    sprintf(HistoNameU,"beamProfU_Det%i",i);
    char HistoTitleU[100];
    sprintf(HistoTitleU,"Beam U coord profile for Det%i",i);
    _BeamProfU[i] = new TH1D(HistoNameU,HistoTitleU,200,_stripULow,_stripUHigh);
      
    char HistoNameV[100];
    sprintf(HistoNameV,"beamProfV_Det%i",i);
    char HistoTitleV[100];
    sprintf(HistoTitleV,"Beam V coord profile for Det%i",i);
    _BeamProfV[i] = new TH1D(HistoNameV,HistoTitleV,200,_stripULow,_stripUHigh);
      
  }
  // profile histo for Pxl detector
  _BeamProfPxl = new TH2D("beamProfPxl", "Beam profile for Pxl det", 
			  200, _pxlULow, _pxlUHigh, 200, _pxlVLow, _pxlVHigh);
   
  
  //build here the layer 0-layer n correlation
  for (int iDet=1; iDet<4; iDet++){

    char HistoNameU[100];
    sprintf(HistoNameU,"beamProfUL0L%i",iDet);
    char HistoTitleU[100];
    sprintf(HistoTitleU,"Beam Profile U: Layer0 - Layer%i Corr",iDet);

    _BeamLayerLayerCorrU[iDet-1] = new TH2D(HistoNameU,HistoTitleU,
					 200, _stripULow, _stripUHigh,
					 200, _stripULow, _stripUHigh);
    
    char HistoNameV[100];
    sprintf(HistoNameV,"beamProfVL0L%i",iDet);
    char HistoTitleV[100];
    sprintf(HistoTitleV,"Beam Profile V: Layer0 - Layer%i Corr",iDet);
    
    _BeamLayerLayerCorrV[iDet-1] = new TH2D(HistoNameV,HistoTitleV,
					 200, _stripVLow, _stripVHigh,
					 200, _stripVLow, _stripVHigh);
    
  }
  char HistoNameU[100];
  sprintf(HistoNameU,"beamProfUL0L4");
  char HistoTitleU[100];
  sprintf(HistoTitleU,"Beam Profile U: Layer0 - Layer4 Corr");
  
  _BeamLayerLayerCorrU[3] = new TH2D(HistoNameU,HistoTitleU,
				     200, _stripULow, _stripUHigh,
				     200, _pxlULow, _pxlUHigh);

  char HistoNameV[100];
  sprintf(HistoNameV,"beamProfVL0L4");
  char HistoTitleV[100];
  sprintf(HistoTitleV,"Beam Profile V: Layer0 - Layer4 Corr");
  
  _BeamLayerLayerCorrV[3] = new TH2D(HistoNameV,HistoTitleV,
				     200, _stripVLow, _stripVHigh,
				     200, _pxlVLow, _pxlVHigh);
    
    
  //build here the layer 0-layer n correlation differences
  for (int iDet=1; iDet<4; iDet++){

    char HistoNameU[100];
    sprintf(HistoNameU,"diffUL0L%i",iDet);
    char HistoTitleU[100];
    sprintf(HistoTitleU,"Diff U: Layer0 - Layer%i",iDet);

    _DiffLayerCorrU[iDet-1] = new TH1D(HistoNameU,HistoTitleU,
				       200, _stripULow/2., _stripUHigh/2.);
    
    char HistoNameV[100];
    sprintf(HistoNameV,"diffVL0L%i",iDet);
    char HistoTitleV[100];
    sprintf(HistoTitleV,"Diff V: Layer0 - Layer%i",iDet);
    
    _DiffLayerCorrV[iDet-1] = new TH1D(HistoNameV,HistoTitleV,
				       200, _stripVLow/2., _stripVHigh/2.);
  }
}




void QCAnal::InitDigiHisto(){

  for (int iDet= 0; iDet< _nDet ; iDet++){
    //fill beam profile histo with cluster info
	  

  //define here the histo for occupancy and ADC
  char HistoNameU[100];
  char HistoTitleU[100];
  char HistoNameV[100];
  char HistoTitleV[100];

  sprintf(HistoNameU,"OccU_Det%i",iDet);
  sprintf(HistoTitleU,"Occupancy U side for Det%i",iDet);
  _OccU[iDet] = new TH1D(HistoNameU,HistoTitleU,384,-0.5,383.5);


  sprintf(HistoNameV,"OccV_Det%i",iDet);
  sprintf(HistoTitleV,"Occupancy V side for Det%i",iDet);
  _OccV[iDet] = new TH1D(HistoNameV,HistoTitleV,384,-0.5,383.5);

  sprintf(HistoNameU,"AdcU_Det%i",iDet);
  sprintf(HistoTitleU,"Adc U side for Det%i",iDet);
  _AdcU[iDet] = new TH1D(HistoNameU,HistoTitleU,8,-0.5,7.5);


  sprintf(HistoNameV,"AdcV_Det%i",iDet);
  sprintf(HistoTitleV,"Adc V side for Det%i",iDet);
  _AdcV[iDet] = new TH1D(HistoNameV,HistoTitleV,8,-0.5,7.5);

  }


  // profile histo for Pxl detector
  _OccPxl = new TH2D("OccPxl", "Occ for Pxl det",128, -0.5, 127.5, 
		     32, -0.5, 31.5);

}

void QCAnal::InitClusterHisto(){


  //define here the histo for clusters
  char HistoNameU[100];
  char HistoTitleU[100];
  char HistoNameV[100];
  char HistoTitleV[100];

  for (int iDet= 0; iDet< _nDet ; iDet++){
    //fill  histo with cluster info


  sprintf(HistoNameU,"nclustU_Det%i",iDet);
  sprintf(HistoTitleU,"Number of cluster U side for Det%i",iDet);
  _nclustU[iDet] = new TH1D(HistoNameU,HistoTitleU,50,-0.5, 49.5);

  sprintf(HistoNameV,"nclustV_Det%i",iDet);
  sprintf(HistoTitleV,"Number of cluster V side for Det%i",iDet);
  _nclustV[iDet] = new TH1D(HistoNameV,HistoTitleV,50,-0.5, 49.5);

  sprintf(HistoNameU,"clustMulU_Det%i",iDet);
  sprintf(HistoTitleU,"Cluster multiplicity U side for Det%i",iDet);
 _clustMulU[iDet] = new TH1D(HistoNameU,HistoTitleU,50,-0.5, 49.5);

  sprintf(HistoNameV,"clustMulV_Det%i",iDet);
  sprintf(HistoTitleV,"Cluster multiplicity V side for Det%i",iDet);
  _clustMulV[iDet] = new TH1D(HistoNameV,HistoTitleV,20,-0.5, 19.5);

  sprintf(HistoNameU,"clustPhU_Det%i",iDet);
  sprintf(HistoTitleU,"Cluster Pulse Height U side for Det%i",iDet);
  _clustPhU[iDet] = new TH1D(HistoNameU,HistoTitleU,20,-0.5, 19.5);

  sprintf(HistoNameV,"clustPhV_Det%i",iDet);
  sprintf(HistoTitleV,"Cluster Pulse Height V side for Det%i",iDet);
  _clustPhV[iDet] = new TH1D(HistoNameV,HistoTitleV,20,-0.5, 19.5);

  }

  _nclustPxl = new TH1D("nclustPxl","Number of cluster for pxl detector",50,-0.5, 49.5);
  _clustMulPxl = new TH1D("clustMulPxl","Multiplicity of the cluster for pxl det",50,-0.5, 49.5);


}



void QCAnal::InitTrkHisto(){

  //define here the histo for clusters

  char HistoName[100];
  char HistoTitle[100];

  sprintf(HistoName,"ntrk");
  sprintf(HistoTitle,"Number of reco tracks");
  _ntrk = new TH1D(HistoName,HistoTitle,50,-0.5, 49.5);

  sprintf(HistoName,"trkSlpX");
  sprintf(HistoTitle,"Track Slope X");
  _trkSlpX = new TH1D(HistoName,HistoTitle,50,-.1, .1);

  sprintf(HistoName,"trkSlpY");
  sprintf(HistoTitle,"Track Slope Y");
  _trkSlpY = new TH1D(HistoName,HistoTitle,50,-.1, .1);

  sprintf(HistoName,"trkItpX");
  sprintf(HistoTitle,"Track Intercept X");
  _trkItpX = new TH1D(HistoName,HistoTitle,50,-1.0, 1.0);

  sprintf(HistoName,"trkItpY");
  sprintf(HistoTitle,"Track Intercept Y");
  _trkItpY = new TH1D(HistoName,HistoTitle,50,-1.0, 1.0);

  sprintf(HistoName,"trkChi2");
  sprintf(HistoTitle,"Track Chi2 ");
  _trkChi2 = new TH1D(HistoName,HistoTitle,50,-0.5, 49.5);

  sprintf(HistoName,"trkProbChi2");
  sprintf(HistoTitle,"Track Prob Chi2 ");
  _trkProbChi2 = new TH1D(HistoName,HistoTitle,50,-0.05,1.05 );

  sprintf(HistoName,"trkFitStatus");
  sprintf(HistoTitle,"Track Fit Status ");
  _trkFitStatus = new TH1D(HistoName,HistoTitle,4,-0.5,3.5 );

}

void QCAnal::InitResidualHisto(){


  //define here the histo for local residuals
  char HistoNameU[100];
  char HistoTitleU[100];
  char HistoNameV[100];
  char HistoTitleV[100];

  for (int iDet= 0; iDet< _nDet ; iDet++){
    //fill  histo with cluster info

  sprintf(HistoNameU,"ResU_Det%i",iDet);
  sprintf(HistoTitleU,"Residual (Meas-Extr) U side for Det%i",iDet);
  _resU[iDet] = new TH1D(HistoNameU,HistoTitleU,1000, -0.0100, 0.0100);

  sprintf(HistoNameV,"ResV_Det%i",iDet);
  sprintf(HistoTitleV,"Residual (Meas-Extr) V side for Det%i",iDet);
  _resV[iDet] = new TH1D(HistoNameV,HistoTitleV,1000, -0.0100, 0.0100);


  //define here the histo for local residuals
  char HistoNameX[100];
  char HistoTitleX[100];
  char HistoNameY[100];
  char HistoTitleY[100];

  sprintf(HistoNameX,"ResX_Det%i",iDet);
  sprintf(HistoTitleX,"Residual (Meas-Extr) X side for Det%i",iDet);
  _resX[iDet] = new TH1D(HistoNameX,HistoTitleX,1000, -1, 1);

  sprintf(HistoNameY,"ResY_Det%i",iDet);
  sprintf(HistoTitleY,"Residual (Meas-Extr) Y side for Det%i",iDet);
  _resY[iDet] = new TH1D(HistoNameY,HistoTitleY,1000, -1, 1);

  // define here the histo for local residuals
  char HistoName2D[100];
  char HistoTitle2D[100];
  sprintf(HistoName2D,"intInside_Det%i",iDet);
  sprintf(HistoTitle2D,"Scatter plot of intercept for Det%i",iDet);
  _intInside[iDet] = new TH2D(HistoName2D,HistoTitle2D,1000, -1, 1, 1000, -1, 1);

  }

}


void QCAnal::Normalize1DHisto(TH1D* h1){

  h1->Scale(double(1.0/_nEntries));

}

void QCAnal::Normalize2DHisto(TH2D* h2){

  h2->Scale(double(1.0/_nEntries));

}