#include "../common/MyStyle.h"
#include "../common/PitCBNT.h"

#include "TH1F.h"
#include "TChain.h"
#include "TFile.h"

#include <iostream>
#include <cmath>


#include "header.hh"

using namespace std;

//
// anaMask.C
//  mask Occu > 1% strip


int main()
{
  MyFunctions::SetMyStyle();
  gStyle->SetPalette(1);
  gStyle->SetCanvasDefW(1200);
  gStyle->SetCanvasDefH(2400);
  
  //====== make Canvas();
  
  //TCanvas *c1 = new TCanvas();
  
  //====== setting ======
  
  Char_t n_nt[30];
  Int_t i_nt;

  const Int_t n_strip = 768; //strip
  
  std::cout << "In ntule 4keta " << endl;
  std::cin >> i_nt; 
  
  sprintf(n_nt,"~/PitCBNT/ntuplePit%d.root",i_nt);
  
  TChain * tChain = new TChain("Cosmic/tree", "");
  tChain->Add(n_nt);
  
  std::cout << tChain->GetEntries() << std::endl;
  PitCBNT * br = new PitCBNT((TTree*)tChain);
  
  // ======= Histgram setting ======
  
  Int_t i_layer;
  Int_t i_phi;
  Int_t i_eta;
  Int_t i_z;
  Int_t i_side;
 
  char s[30];
  
  Int_t n_layer=4;
  Int_t n_phi=56;
  Int_t n_z=2;
  Int_t n_eta=7;
  Int_t n_side=2;
  Int_t n_phimax[4]={32,40,48,56};

  TH1F *h[n_layer][n_phi][n_z][n_eta][n_side];
  
  for(i_layer=0;i_layer<=3;i_layer++){
    for(i_phi=0;i_phi< n_phimax[i_layer] ;i_phi++){
      for(i_eta=-6;i_eta<=6;i_eta++){
	if(i_eta != 0){
	  if(i_eta > 0)i_z = 1;else i_z =0;
           for(i_side=0;i_side<=1;i_side++){
	  sprintf(s,"Layer%d_Phi%d_Eta%d_Side%d",i_layer,i_phi,i_eta,i_side);
	  h[i_layer][i_phi][i_z][abs(i_eta)][i_side] = new TH1F(s,s,n_strip,0,n_strip);
         }
	}
      }
    }
   }
  
  //===== To start Loop =====
  Int_t nEvents=0;
  Int_t nentries=0;
  Int_t true_ent=0;
  
  if(nEvents==0) nentries = (Int_t)tChain->GetEntries();
  else nentries=nEvents;
  
  for(Int_t jentry=0;jentry<nentries;jentry++){
    if(jentry%1000==0)cout<<"# "<<jentry/1000<<"K th  event start"<<endl;
    Int_t ientry = br->LoadTree(jentry);
    if(ientry<0)break;
    tChain->GetEntry(jentry);

    true_ent = jentry;
     
    if(br->sct_rdo_nhits < 10){break;}
 
    for(Int_t i=0;i<br->sct_rdo_nhits;i++){
      //tbin XXX
      //if(br->sct_rdo_tbin[i] == 4 || br->sct_rdo_tbin[i] == 5 || br->sct_rdo_tbin[i] == 6 || br->sct_rdo_tbin[i] == 7){ //1XX
      if(br->sct_rdo_tbin[i] == 2 || br->sct_rdo_tbin[i] == 3 || br->sct_rdo_tbin[i] == 6 || br->sct_rdo_tbin[i] == 7){ //X1X
	//if(br->sct_rdo_tbin[i] == 1 || br->sct_rdo_tbin[i] == 3 || br->sct_rdo_tbin[i] == 5 || br->sct_rdo_tbin[i] == 7){ //XX1
	//if(br->sct_rdo_tbin[i] == 1 || br->sct_rdo_tbin[i] == 2 || br->sct_rdo_tbin[i] == 4){
	//if(br->sct_rdo_tbin[i] != 0 ){
	
	i_layer = (Int_t)br->sct_rdo_layer[i];
	i_phi   = (Int_t)br->sct_rdo_phi[i];
	i_eta   = (Int_t)br->sct_rdo_eta[i];
	i_side  = (Int_t)br->sct_rdo_side[i];


	if(i_eta > 0)i_z = 1;else i_z =0;
	h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->Fill(br->sct_rdo_strips[i]);
      }//end if tbin
    } //end hit
  } //end event

  TFile fsave("cosmicMask.root", "RECREATE");
  
  TH1F *f[2][4];
       f[0][0] = new TH1F("Occ00","Occ00",500,0,500);
       f[0][1] = new TH1F("Occ01","Occ01",500,0,500);
       f[0][2] = new TH1F("Occ02","Occ02",500,0,500);
       f[0][3] = new TH1F("Occ03","Occ03",500,0,500);
       f[1][0] = new TH1F("Occ10","Occ10",500,0,500);
       f[1][1] = new TH1F("Occ11","Occ12",500,0,500);
       f[1][2] = new TH1F("Occ12","Occ13",500,0,500);
       f[1][3] = new TH1F("Occ13","Occ14",500,0,500);
  TH2F *occuMap[2][4];
  for (i_layer = 0; i_layer < 4; i_layer++){
     for( i_side = 0; i_side <2 ; i_side++){
       sprintf(s,"OccuMap_layer%d_side%d_phi_vs_eta",i_layer,i_side);
       occuMap[i_side][i_layer] = new TH2F(s,s,13,-6.,6.,n_phimax[i_layer],0.,(Double_t)n_phimax[i_layer]);
     }
  }

  Double_t sum = 0;
  Double_t ent = nentries;
  Int_t index[2] = {0};
  Double_t all[2] = {0};
  Double_t allocc[2][4] = {{0,0,0,0},{0,0,0,0}};
  
  Double_t hitLimit = 0.01*(Double_t)nentries;//suppress hot channel
  fsave.mkdir("BadSideHits");
  fsave.cd("BadSideHits");

  for(i_layer=0;i_layer<=3;i_layer++){
    index[0]=0;
    index[1]=0;
    for(i_phi=0;i_phi< n_phimax[i_layer] ;i_phi++){
      for(i_eta=-6;i_eta<=6;i_eta++){
	if(i_eta != 0){if(i_eta > 0)i_z = 1;else i_z =0;
	 for(i_side=0;i_side<=1;i_side++){

        Double_t moduOccu = 0;
        sum = 0;

	for(Float_t i_strip=0 ;i_strip <= n_strip;i_strip++){
	  if(i_eta > 0)i_z = 1;else i_z =0;
          Double_t hitCount = (Double_t) h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->GetBinContent((int)i_strip);

	  if(hitCount < hitLimit ){
		 sum += hitCount;
	  }
	  else{
		 h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->SetBinContent((int)i_strip,0.);// remove hot channel but Entries unchanged
	  }
	} //end i_strip
	
	sprintf(s,"Layer%d_Phi%d_Eta%d",i_layer,i_phi,i_eta);
	if(sum != 0){
        ent = true_ent;
  //      moduOccu = h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->GetEntries() / (ent * n_strip);
        moduOccu = sum / (ent * n_strip);
	if (moduOccu > 50e-5){
            std::cout << s << " Occupancy =" << moduOccu << std::endl;
            h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->Write();
        }
        all[i_side] += moduOccu;
        allocc[i_side][i_layer] += moduOccu;
        index[i_side]++;
        f[i_side][i_layer]->SetBinContent(index[i_side],moduOccu);
        occuMap[i_side][i_layer]->SetBinContent(i_eta+7,i_phi+1,moduOccu);
        } // bin 0:under flow 1 throu nbin: contents, nbin+1 overflow

          } //end side	
	}} //end eta
    } //end phi
  } //end layer
  
  std::cout << "side0 =" << all[0] << std::endl;
  std::cout << "side1 =" << all[1] << std::endl;

  std::cout << "true_ent =" << true_ent << endl;
  fsave.cd();
  fsave.mkdir("sideHits");
  fsave.cd("sideHits");
    for(i_layer=0;i_layer<=3;i_layer++){
    for(i_phi=0;i_phi< n_phimax[i_layer] ;i_phi++){
      for(i_eta=-6;i_eta<=6;i_eta++){
        if(i_eta != 0){
          if(i_eta > 0)i_z = 1;else i_z =0;
           for(i_side=0;i_side<=1;i_side++){
          h[i_layer][i_phi][i_z][abs(i_eta)][i_side]->Write();
         }
        }
      }
    }
   }
                                                                                
  fsave.cd();
  fsave.mkdir("moduleOccu");  
  fsave.cd("moduleOccu");

  Double_t ave[2][4] ={{0,0,0,0},{0,0,0,0}};
  Double_t layerModules[4]={268,339,396,448};

  for(i_layer=0;i_layer<=3;i_layer++){
    for(i_side=0;i_side<=1;i_side++){
  f[i_side][i_layer]->Write();
  occuMap[i_side][i_layer]->Write();
  ave[i_side][i_layer] = allocc[i_side][i_layer] / layerModules[i_layer]; 
  std::cout << "ave of " << i_layer << " =" <<   ave[i_side][i_layer] << std::endl;
   } //end side   
  } //end layer
  fsave.Close();

  return 0;
  
}

