#define Nada_look_cxx
#include "Nada_look.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"

void Nada_look::Loop()
{
//   In a Root session, you can do:
//      Root > .L Nada_look.C
//      Root > Nada_look t
//      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
//       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(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch
   if (fChain == 0) return;

// Define 1D Histogramms
   TH1F *ecal[2];
   TH1F *eta[2];
   TH1F *phi[2];
   TH1F *lyr[2];
   TH1F *nnad[2];
   TH1F *pt[2];
   TH1F *ptjet[5];
   TH1F *etm[5];
   
// Name them ... 

   char hcyes[500];
   char hcno[500];
   char leset[300];
   sprintf(hcyes,"due to HC");
   sprintf(hcno,"HC suppressed");

   for (Int_t j=0; j<5;j++) {
     char ptjn[500];
     char ptjnb[500];
     char etmn[500];
     char etmnb[500];
     char lo[500]; 
     sprintf(ptjnb,"Pt jets");
     sprintf(lo,"leading object");
     sprintf(etmnb,"Et Miss");
     if (j == 0) {
       sprintf(ptjn,"%s ",ptjnb);
       sprintf(etmn,"%s ",etmnb);}

     if (j == 1) {
       sprintf(ptjn,"%s, %s",ptjnb,lo);
       sprintf(etmn,"%s  %s ,set 1 ",etmnb,hcno);}

     if (j == 2) {
       sprintf(ptjn,"%s matching a HC set 1 ",ptjnb);
       sprintf(etmn,"%s  %s ,set 2 ",etmnb,hcno);}

     if (j == 3) {
       sprintf(ptjn,"%s %s, set 1",ptjnb,hcno);
       sprintf(etmn,"%s  %s ,set 1 ",etmnb,hcyes);}

     if (j == 4) {
       sprintf(etmn,"%s  %s ,set 2 ",etmnb,hcyes);
       sprintf(ptjn,"Nothing 2");}

     ptjet[j]=new TH1F(ptjn,ptjn,50,-1,250);
     etm[j]=new TH1F(etmn,etmn,50,0,250);
   }
   
   for (Int_t j=0; j<2;j++) {
     char nolayer[500];
     char ecaln[300];
     char etan[300];
     char phin[300];
     char lyrn[300];
     char ptn[300];
     char nnadn[300];
     
     Int_t les;
     les=j+1;

     sprintf(leset,"set %d",les);
     sprintf(ecaln,"E Cal ; %s",leset);
     sprintf(etan,"Eta ; %s",leset);
     sprintf(phin,"Phi ; %s",leset);
     sprintf(lyrn,"Layer ; %s",leset);
     sprintf(ptn,"Pt ; %s",leset);
     sprintf(nnadn,"N Nada cell ; %s",leset);
   

     ecal[j]=new TH1F(ecaln,ecaln,50,0,250);

     eta[j]=new TH1F(etan,etan,101,-50,50);

     phi[j]=new TH1F(phin,phin,65,0,64);
     
     pt[j]=new TH1F(ptn,ptn,50,0,250);
     
     lyr[j]=new TH1F(lyrn,lyrn,18,-0.5,17.5);
     
     nnad[j]=new TH1F(nnadn,nnadn,21,0.,20);

   }

 //  Variable of minimal approach to match a jet and a HC 
   Float_t drcone = 0.2;

// 2D Histo :

   TH2F *etaphij = new TH2F("etaphij"," phi vs eta- jets and HC (red)" ,40,-4,4,35,0,7);
   TH2F *etaphinad = new TH2F("etaphinad"," phi vs eta- HC set 1 " ,40,-4,4,35,0,7);

   TH2F *etaphijdr = new TH2F("etaphijdr"," phi vs eta- jet match HC set 1" ,40,-4,4,35,0,7);
   TH2F *etaphinaddr = new TH2F("etaphinaddr"," phi vs eta- HC matched with jets" ,40,-4,4,35,0,7);

   TH2F *pphijdr = new TH2F("pphijdr"," Pt vs eta- jets match HC set 1" ,40,-4,4,50,0,250);
//////////////////////////////

// Will loop on events :

   Int_t nentries = Int_t(fChain->GetEntries());
   Int_t nbytes = 0, nb = 0;
   cout << " nentries " << nentries << endl;

///////////////

   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;


 ///////////////////////////////
//  Et miss calculation  : Initialization
//  3 variables : i = x,y, Et ; j = #set ; k = 1 : all  2 : mEt HC suppressed ;mEt  due to HC ;
      Float_t etmiss[3][2][3];
      for (Int_t i=0; i<3;i++) {
	for (Int_t j=0; j<2;j++) {
	  for (Int_t k=0; k<3;k++) {
	    etmiss[i][j][k]=0.;
	  }
	}
      }
      Float_t letheta;
      for (Int_t calnum=0; calnum<CAL_ncal;calnum++) {
	if (CAL_e_cal[calnum] > 0.1) {
	// limit at 100 Mev on ecal
	  letheta = 2*(atan(exp(-(CAL_ieta_cal[calnum]/10.)))) ;
	  etmiss[0][0][0]=etmiss[0][0][0]+CAL_e_cal[calnum]*sin(letheta)*cos(CAL_iphi_cal[calnum]/10.);
	  etmiss[1][0][0]=etmiss[1][0][0]+CAL_e_cal[calnum]*sin(letheta)*sin(CAL_iphi_cal[calnum]/10.);
	}
      }
// Same Et for both sets
      etmiss[0][1][0]=etmiss[0][0][0];
      etmiss[1][1][0]=etmiss[1][0][0];
// At the beginning, same Et for all and met due to HC
      etmiss[0][0][1]=etmiss[0][0][0];
      etmiss[1][0][1]=etmiss[1][0][0];

      etmiss[0][1][1]=etmiss[0][0][0];
      etmiss[1][1][1]=etmiss[1][0][0];

// End of Et initialization
///////////////////////////////


///// Initialize jets variables
      Float_t pt_cc[10];
      for (Int_t ijet =0; ijet<10 ; ijet++) {
	pt_cc[ijet]=0;
      }
      Float_t leadjet=0.001;
      Int_t ileadjet=100;
      for (Int_t ijet =0; ijet<JETS_c_njets_c ; ijet++) {
	pt_cc[ijet]=JETS_c_pt_c[ijet];
	if (JETS_c_pt_c[ijet] > leadjet) {
	  ileadjet=ijet;
	  leadjet=JETS_c_pt_c[ijet];
	}
	ptjet[0]->Fill(JETS_c_pt_c[ijet]);
	etaphij->Fill(JETS_c_eta_c[ijet],JETS_c_phi_c[ijet]);
      //////////////////////////
//  See if jets matches a HC of set 1
	Int_t matchetil = 0;
//
	for (Int_t i=0; i<NADA_nnada[0];i++) {
	  Float_t  Dr = Ledr(NADA_ieta_nad[i][0]/10.,NADA_iphi_nad[i][0]/10.,JETS_c_eta_c[ijet],JETS_c_phi_c[ijet]);
	  if(Dr < drcone) {
//		 cout << "no of entry "<< jentry << endl;
//		 cout << " Dr  = " << Dr << endl;
	    matchetil = 1;
	    etaphinaddr->Fill(NADA_ieta_nad[i][0]/10.,NADA_iphi_nad[i][0]/10.);
	  // Suppress HC Pt : 
	    pt_cc[ijet]=pt_cc[ijet]-NADA_pt_nad[i][0];
//	    cout << "no of entry "<< jentry << endl;
//	    cout << "pt_cc[ijet]" << pt_cc[ijet] << endl;
	  }
	}
	if (matchetil == 1) { 
	  ptjet[2]->Fill(JETS_c_pt_c[ijet]);
	  etaphijdr->Fill(JETS_c_eta_c[ijet],JETS_c_phi_c[ijet]);
	  ptjet[3]->Fill(pt_cc[ijet]);
	  pphijdr->Fill(JETS_c_eta_c[ijet],JETS_c_pt_c[ijet]);
	}
      }
      if (ileadjet != 100) ptjet[1]->Fill(JETS_c_pt_c[ileadjet]);
//
///////////////////// 
    //  Turn on nada :
//////////////////////
//loop on set : 
      for (Int_t j=0; j<2;j++) { 
//
      //loop on hot cell
	for (Int_t i=0; i<NADA_nnada[j];i++) {

/*   Check matching between calo and nada (to see if it already suppressed in the code). 	
     Int_t cestok=0;
     for (Int_t calnum=0; calnum<CAL_ncal;calnum++) {
     if (CAL_ieta_cal[calnum] == NADA_ieta_nad[i][j] && CAL_iphi_cal[calnum] == NADA_iphi_nad[i][j] && CAL_ilyr_cal[calnum] == NADA_ilyr_nad[i][j] ) {
     if (CAL_e_cal[calnum] == NADA_e_nad[i][j]) {
     cestok=1;
     cout << "calnum " << calnum << endl; 
     cout << " i " << i << " j " << j << endl;
     }
     else {
     cout << "e cal " << CAL_e_cal[calnum] << endl;
     cout << "e nad " << NADA_e_nad[i][j] << endl;
     }
     }
     }
     if (cestok == 1) cout << "OK it is matched" << endl;
     if (cestok == 0) cout << "It is not matched" << endl;
*/
	  nnad[j]->Fill(NADA_nnada[j]);

	  ecal[j]->Fill(NADA_e_nad[i][j]);

	  eta[j]->Fill(NADA_ieta_nad[i][j]);
	  
	  phi[j]->Fill(NADA_iphi_nad[i][j]);
	  
	  pt[j]->Fill(NADA_pt_nad[i][j]);
	  
	  lyr[j]->Fill(NADA_ilyr_nad[i][j]);
	  if (j == 0) etaphinad->Fill(NADA_ieta_nad[i][j]/10.,NADA_iphi_nad[i][j]/10.);
//
/////////////////////////
// Calculate Etmiss without hot cell contribution :
// => calculate Et miss for x and y and subtracte it from all Etmiss 

	  letheta = 2*(atan(exp(-(NADA_ieta_nad[i][j]/10.)))) ;
	//  Et miss on x nada subtracted
	  etmiss[0][j][1]=etmiss[0][j][1]-NADA_e_nad[i][j]*sin(letheta)*cos(NADA_iphi_nad[i][j]/10.);
	//  Et miss on y nada subtracted
	  etmiss[1][j][1]=etmiss[1][j][1]-NADA_e_nad[i][j]*sin(letheta)*sin(NADA_iphi_nad[i][j]/10.);
	//  Et miss on x nada only
	  etmiss[0][j][2]=etmiss[0][j][2] + NADA_e_nad[i][j]*sin(letheta)*cos(NADA_iphi_nad[i][j]/10.);
	//  Et miss on y nada only
	  etmiss[1][j][2]=etmiss[1][j][2]+NADA_e_nad[i][j]*sin(letheta)*sin(NADA_iphi_nad[i][j]/10.);
	}
      
	etmiss[2][j][0] = sqrt((etmiss[0][j][0])**2+(etmiss[1][j][0])**2);
	etmiss[2][j][1] = sqrt((etmiss[0][j][1])**2+(etmiss[1][j][1])**2);
	etmiss[2][j][2] = sqrt((etmiss[0][j][2])**2+(etmiss[1][j][2])**2);
	
      }
      
      if (etmiss[2][0][0] > 0.01) etm[0]->Fill(etmiss[2][0][0]);
      if (etmiss[2][0][1] > 0.01) etm[1]->Fill(etmiss[2][0][1]);
      if (etmiss[2][1][1] > 0.01) etm[2]->Fill(etmiss[2][1][1]);
      if (etmiss[2][0][2] > 0.01) etm[3]->Fill(etmiss[2][0][2]);
      if (etmiss[2][1][2] > 0.01) etm[4]->Fill(etmiss[2][1][2]);
   }   
//////////////////////////////   
//   Make the plots !!!
/////////////////////////////
   TCanvas *can = new TCanvas("can","test",600,800);
   can->SetFillColor(0);
   can->SetFillStyle(4000);
// Define Postscript
   Int_t type = 111;
   TPostScript ps("nada_look.ps",type);

   gStyle->SetOptStat(111111);
   gStyle->SetStatW(0.3);
   gStyle->SetTextSize(1.2);
//
// Make the plots  
//// The number of HC for each set
   ps->NewPage();
   can->Divide(1,2);
   for (Int_t ecran=0; ecran<2;ecran++) {
     Int_t ecranb=ecran+1; 
     can->cd(ecranb);
     nnad[ecran]->SetLineWidth(2); 
     nnad[ecran]->Draw(); 
   }
   can.Update();
// HC parameters for  set 1 and 2
   for (Int_t set=0; set<2; set++) {
     ps->NewPage();
     can->Clear();
     can->Divide(2,3);
     for (Int_t ecran=0; ecran<5;ecran++) {
       ecranb=ecran+1;
       can->cd(ecranb);
       gStyle->SetOptStat(111111);
       gStyle->SetStatW(0.5);
       gStyle->SetTextSize(1.2);
       if (ecran == 0) {
	 ecal[set]->SetLineWidth(2); 
	 ecal[set]->Draw(); 
       }
       if (ecran == 1) {
	 pt[set]->SetLineWidth(2); 
	 pt[set]->Draw(); 
       }
       if (ecran == 2) {
	 eta[set]->SetLineWidth(2); 
	 eta[set]->Draw(); 
       }
       if (ecran == 3) {
	 phi[set]->SetLineWidth(2); 
	 phi[set]->Draw(); 
       }
       if (ecran == 4) {
	 lyr[set]->SetLineWidth(2); 
	 lyr[set]->Draw(); 
       }
     }
     can.Update();
   }
//
// Draw etmiss
   ps->NewPage();
   can->Clear();
   can->Divide(2,2);
   for (Int_t ecran=0; ecran<4; ecran++) {
     ecranb=ecran+1;
     can->cd(ecranb);
     gStyle->SetOptStat(111111);
     gStyle->SetStatW(0.5);
     gStyle->SetTextSize(1.2);
     if (ecran > 0 ) Int_t ecrant = ecran-1;
     if (ecran == 0 ) Int_t ecrant = ecran;
     etm[ecrant]->SetLineWidth(2); 
     etm[ecrant]->Draw(); 
   }
   can.Update();
// Draw jets

// First 2 D eta-phi
   gStyle->SetOptStat(111111);
   gStyle->SetStatW(0.3);
   gStyle->SetTextSize(1.2);
   ps->NewPage();
   can->Clear();
   can->Divide(1,2);
   can->cd(1);
   etaphij->Draw("box");
   etaphinad->SetFillColor(2);
   etaphinad->Draw("samebox");
   can->cd(2);
   etaphinad->SetFillColor(2);
   etaphinad->Draw("cont1");
   can.Update();
//  
   ps->NewPage();
   can->Divide(1,2);
   can->cd(1);
   etaphijdr->Draw("box");
   etaphinaddr->SetFillColor(2);
   etaphinaddr->Draw("samebox");
   can->cd(2);
   pphijdr->Draw("box");
   can.Update();  
//
   ps->NewPage();
   can->Clear();
   can->Divide(2,2);
   for (Int_t ecran=0; ecran<4;ecran++) {
     Int_t ecranb=ecran+1;
     can->cd(ecranb);
     gStyle->SetOptStat(111111);
     gStyle->SetStatW(0.5);
     gStyle->SetTextSize(1.2);
     ptjet[ecran]->SetLineWidth(2); 
     ptjet[ecran]->Draw();
   }  
   can.Update();
   ps.Close();
}
// Define function
Float_t 
Nada_look::Ledr(Float_t& eh, Float_t& ph, Float_t& ej, Float_t& pj)
//Float_t Nada_look::Ledr(Float_t *eh, Float_t *ph, Float_t *ej, Float_t *pj)
     {
       Float_t pi = 3.1415926;
       Float_t dp=fabs(pj-ph);
       if (dp > pi) dp=(2*pi-dp);
       Float_t dr = sqrt(dp*dp + (ej-eh)*(ej-eh));
       return dr;
     }