//==================================================================
//
//	TJetCand.cc
//
// This code TEMPLATE plots quantities from the stntuple TStnJetBlock
// for three block descriptions: "JetBlock", "TopTJetBlock", and "TopLJetBlock"
//
//	Version 0.0:      P. Koehn    05/05/01
//	Version 0.1:      P. Koehn    03/14/02
// 		Tagging information is coming soon
//
//==================================================================
#include <iostream>
#include <iomanip>
#include <string>
#include <algorithm>
#include "Stntuple/loop/TStnAna.hh"
#include "TJetCand.hh"

//==================================================================
//TJetCand(TFile* file, const char* name, const char* title):
//==================================================================
TJetCand::TJetCand(TFile* file, const char* name, const char* title):TStnModule(name,title){

// create a subdirectory "dir" in this file
// 	(assumed created prior to this call)
   histFile = file;
   histDir = name;
   dir = histFile->mkdir(histDir);
	
//
// Initialize cuts
	init();
}

//==================================================================
// TJetCand(const char* name, const char* title):
//==================================================================
TJetCand::TJetCand(const char* name, const char* title):TStnModule(name,title){

// create a subdirectory "dir" in this file
// 	(assumed created prior to this call)
   histFile = new TFile("TJetCand.root","RECREATE","test file");
   histDir = name;
   dir = histFile->mkdir(histDir);
	
//
// Initialize cuts
	init();
}


//==================================
// ~TJetCand() 
//==================================
TJetCand::~TJetCand() {}

//==================================
//       init()
//==================================
void TJetCand::init() {

  // set the "description name" of the stntuple
  // data blocks that you need to access
  desc[0] = "JetBlock";
  desc[1] = "TopTJetBlock";
  desc[2] = "TopLJetBlock";
	
  NJetMax = 20;   
  for(Int_t i=0; i < 3; i++){
   NgoodJet[i]=0;   
   totEvents[i]=0;   
   jetEvents[i]=0; 
  }
  
  enableJetHist(false);  
  enableJetEtCut(true); setJetEtCut(10.0);
  enableJetEmfCut(true); setJetEmfCut(0.8);
  enableJetDtetaCut(true); setJetDtetaCut(2.0);
  
}

//==================================
// BeginJob()
//==================================
int TJetCand::BeginJob() {

  printf("----- Begin job: ---- %s\n",GetName());

//
// register the data blocks and enable histograms:
//
// jetblock[i] is the boolean to indicate which TStnJetBlock jet blocks are available:
//  i = 0 -> JetBlock
//  i = 1 -> TopTJetBlock
//  i = 2 -> TopLJetBlock
//
//  STNTUPLE POINTER TO DATA BLOCK SYNTAX: 
//  pointer = ("data block class"*) fAna->RegisterDataBlock("description name","data block class");
//
// Default Jets
  pJetBlock = (TStnJetBlock*) fAna->RegisterDataBlock(desc[0],"TStnJetBlock");
  jetblock[0] = true;
  if (! pJetBlock) {
		printf(" >>> branch *** %s *** doesn't exist \n",desc[0]);
		enableJetHist(false);  // default jets
		jetblock[0] = false;
  }
//
// Tight TopEvent Jets
  pTopTJetBlock = (TStnJetBlock*) fAna->RegisterDataBlock(desc[1],"TStnJetBlock");
  jetblock[1] = true;
  if (! pTopTJetBlock) {
		printf(" >>> branch *** %s *** doesn't exist \n",desc[1]);
		enableTJetHist(false); // tight top jets
		jetblock[1] = false;
  }
//
// Loose Top event Jets
  pTopLJetBlock = (TStnJetBlock*) fAna->RegisterDataBlock(desc[2],"TStnJetBlock");
  jetblock[2] = true;
  if (! pTopLJetBlock) {
		printf(" >>> branch *** %s *** doesn't exist \n",desc[2]);
		enableLJetHist(false); // loose top jets	 
		jetblock[2] = false;
  }

	BookHistograms();
  
	printcuts();

	return 0;
}

//==================================
// BeginRun() 
//==================================
int TJetCand::BeginRun() {
  return 0;
}

//==================================
// BookHistograms()
//==================================
void TJetCand::BookHistograms() {

  printf("----- BookHistograms ---- %s\n",GetName());

   histFile->cd(histDir);    // make the "dir" directory the current directory

// loop for TStnJetBlock jet blocks available:
//  i = 0 -> JetBlock		   // default jets 
//  i = 1 -> TopTJetBlock		// tight jets from topeventmodule
//  i = 2 -> TopLJetBlock		// loose jets from topeventmodule
	for (Int_t i = 0; i < 3 && jetblock[i]; i++){ 

	  	NJetHist[i] = new TH1F("NJetHist","Njet",11,-0.5,10.5);
		NJetHist[i]->SetXTitle("Njet"); NJetHist[i]->SetYTitle("events");
		
		JetNtracksHist[i] = new TH1F("JetNtracksHist","Ntracks",31,-0.5,30.5);
		JetNtracksHist[i]->SetXTitle("number o' tracks");JetNtracksHist[i]->SetYTitle("jets");
	
		JetEtHist[i] = new TH1F("JetEtHist","Jet Et",150,0,300);
		JetEtHist[i]->SetXTitle("Et (GeV)");JetEtHist[i]->SetYTitle("jets");
	
		JetEmfHist[i] = new TH1F("JetEmfHist","Jet EMF",120,0.0,1.2);
		JetEmfHist[i]->SetXTitle("Em fraction");JetEmfHist[i]->SetYTitle("jets");
	
		JetDtetaHist[i] = new TH1F("JetDtetaHist","Detector Eta",100,-4,4);
		JetDtetaHist[i]->SetXTitle("detector eta");JetDtetaHist[i]->SetYTitle("jets");
	
		JetSeedEtaHist[i] = new TH1F("JetSeedEtaHist","Seed Eta",100,-4,4);
		JetSeedEtaHist[i]->SetXTitle("Seed eta");JetSeedEtaHist[i]->SetYTitle("jets");
	
		JetSeedPhiHist[i] = new TH1F("JetSeedPhiHist","Seed Phi",126,0,6.3);
		JetSeedPhiHist[i]->SetXTitle("Seed Phi (radians)");JetSeedPhiHist[i]->SetYTitle("jets");
	
		JetPhiHist[i] = new TH1F("JetPhiHist","jet Phi",126,0,6.3);
		JetPhiHist[i]->SetXTitle("jet Phi (radians)");JetPhiHist[i]->SetYTitle("jets");
	
		JetRapHist[i] = new TH1F("JetRapHist","jet rapidity",100,-4,4);
		JetRapHist[i]->SetXTitle("jet pseudorapidity");JetRapHist[i]->SetYTitle("jets");
	
		JetZvtxHist[i] = new TH1F("JetZvtxHist","z vtx",200,-100,100);
		JetZvtxHist[i]->SetXTitle("z vertex (cm)");JetZvtxHist[i]->SetYTitle("jets");
	
		JetMassHist[i] = new TH1F("JetMassHist","mass",100,0.0,50);
		JetMassHist[i]->SetXTitle("jet mass");JetMassHist[i]->SetYTitle("jets");
	
		JetCorfdHist[i] = new TH1F("JetCorfdHist","JetCand cor FD",50, 0.0, 5.0);
		JetCorfdHist[i]->SetXTitle("correction factor FD");JetCorfdHist[i]->SetYTitle("jets");
	
		JetCorfmHist[i] = new TH1F("JetCorfmHist","JetCand cor FM",50, 0.0, 5.0);
		JetCorfmHist[i]->SetXTitle("correction factor FM");JetCorfmHist[i]->SetYTitle("jets");
	
//----- with cuts --------
		NJetCutHist[i] = new TH1F("NJetCutHist","Njet JetCand",11,-0.5,10.5);
		NJetCutHist[i]->SetXTitle("jets");NJetCutHist[i]->SetYTitle("events");

		JetEtCutHist[i] = new TH1F("JetEtCutHist","JetCutCand Et",150,0,300);
		JetEtCutHist[i]->SetXTitle("Et (GeV)");JetEtCutHist[i]->SetYTitle("jets");
	
		JetDtetaCutHist[i] = new TH1F("JetDtetaCutHist","JetCutCand Det Eta",120,-6,6);
		JetDtetaCutHist[i]->SetXTitle("det eta");JetDtetaCutHist[i]->SetYTitle("jets");
	
		JetEmfCutHist[i] = new TH1F("JetEmfCutHist","JetCutCand Emf",150,0.0,1.5);
		JetEmfCutHist[i]->SetXTitle("Em fraction");JetEmfCutHist[i]->SetYTitle("jets");
 		
	   cout << "\n Histograms for jet block " << desc[i] << " are booked." << endl;
	}
	
} //BookHistograms()

//==================================
// PlotHistograms()
//==================================
void TJetCand::PlotHistograms() {

//
//  Plot for all jets in all 3 blocks
//
// 3 possible TStnJetBlock jet blocks :
	char title[3][50]  = {{"All Jets"},{"TopEventModule tight jet candidates"},
	                                   {"TopEventModule loose jet candidates"}};
	
	
	TCanvas* can[3];
	can[0] = new TCanvas(title[0], title[0], 80, 80,700,700);
	can[1] = new TCanvas(title[1], title[1], 90, 90,700,700);
	can[2] = new TCanvas(title[2], title[2],100,100,700,700);
	can[0]->Divide(3,4);
	can[1]->Divide(3,4);
	can[2]->Divide(3,4);
	
// loop for available TStnJetBlock jet blocks :
//  i = 0 -> JetBlock		   // default jets 
//  i = 1 -> TopTJetBlock		// tight jets from topeventmodule
//  i = 2 -> TopLJetBlock		// loose jets from topeventmodule

	for (Int_t i = 0; i < 3 && jetblock[i]; i++){

		printf("----- PlotHistograms(): ---- %s\n",GetName());	
		can[i]->cd(1);
		NJetHist[i]->Draw();
		can[i]->cd(2);
		JetEtHist[i]->Draw();
		can[i]->cd(3);
		JetEmfHist[i]->Draw();
		can[i]->cd(4);
		JetDtetaHist[i]->Draw();
		can[i]->cd(5);
		JetRapHist[i]->Draw();
		can[i]->cd(6);
		JetZvtxHist[i]->Draw();
		can[i]->cd(7);
		JetMassHist[i]->Draw();
		can[i]->cd(8);
		JetPhiHist[i]->Draw();
		can[i]->cd(9);
		JetNtracksHist[i]->Draw();
		can[i]->cd(10);
		JetCorfdHist[i]->Draw();
		can[i]->cd(11);
		JetCorfmHist[i]->Draw();

		can[i]->Modified();
		can[i]->Update();
		
		cout << "\n Histograms for jet block " << desc[i] << " are being plotted." << endl;
	}
//
//  Plot cut histograms for all jets in all 3 blocks
//
	char cuttitle[3][50]  = {{"Cut histos:All Jets"},{"Cut histos:TopEventModule tight jet candidates"},
	                                   {"Cut histos:TopEventModule loose jet candidates"}};
	
	TCanvas* cutcan[3];
	cutcan[0] = new TCanvas(cuttitle[0], cuttitle[0], 85, 80,500,500);
	cutcan[1] = new TCanvas(cuttitle[1], cuttitle[1], 95, 90,500,500);
	cutcan[2] = new TCanvas(cuttitle[2], cuttitle[2],105,100,500,500);
	cutcan[0]->Divide(2,2);
	cutcan[1]->Divide(2,2);
	cutcan[2]->Divide(2,2);
	
// loop for available TStnJetBlock jet blocks :
//  i = 0 -> JetBlock		   // default jets 
//  i = 1 -> TopTJetBlock		// tight jets from topeventmodule
//  i = 2 -> TopLJetBlock		// loose jets from topeventmodule

	for (Int_t i = 0; i < 3 && jetblock[i]; i++){

		printf("----- Plot Cut Histograms(): ---- %s\n",GetName());	
		cutcan[i]->cd(1);
		NJetCutHist[i]->Draw();
		
		cutcan[i]->cd(2);
		JetEtCutHist[i]->Draw();
		
		cutcan[i]->cd(3);
		JetEmfCutHist[i]->Draw();
		
		cutcan[i]->cd(4);
		JetDtetaCutHist[i]->Draw();

		cutcan[i]->Modified();
		cutcan[i]->Update();

		cout << "\n cut Histograms for jet block " << desc[i] << " are being plotted." << endl;
	}
	
} // PlotHistograms()

//==================================
//       Event(int ientry)
//==================================
int TJetCand::Event(int ientry) {

	if (fDebugLevel > 0)printf(" ----TJetCand ientry = %7i\n",ientry);
		
// loop for TStnJetBlock jet blocks available:
//  i = 0 -> JetBlock		   // default jets 
//  i = 1 -> TopTJetBlock		// tight jets from topeventmodule
//  i = 2 -> TopLJetBlock		// loose jets from topeventmodule
	for (Int_t i = 0; i < 3 && jetblock[i]; i++){
	
// How many events have we processed ?
	   totEvents[i]++;

		foundJet[i] = false;
		TLorentzVector * Mo;
		NgoodJet[i] = 0;
		njets[i] = 0;
		ncutjets[i] = 0;
		
		// default
		if (i==0 ){
			pJetBlock->GetEntry(ientry);
			njets[i] = pJetBlock->NJets();
			if (pJetBlock->NJets() > NJetMax) {
				printf("Jet number (%d) too high! set to %d\n",pJetBlock->NJets(),NJetMax);
				njets[i] = NJetMax;
			}
		}
		// topeventmodule tight jets	
		if (i==1 ){
			pTopTJetBlock->GetEntry(ientry);
			njets[i] = pTopTJetBlock->NJets();
			if (pTopTJetBlock->NJets() > NJetMax) {
				printf("Jet number (%d) too high! set to %d\n",pTopTJetBlock->NJets(),NJetMax);
				njets[i] = NJetMax;
			}
		}
		// topeventmodule loose jets	
		if (i==2 ){
			pTopLJetBlock->GetEntry(ientry);
			njets[i] = pTopLJetBlock->NJets();
			if (pTopLJetBlock->NJets() > NJetMax) {
				printf("Jet number (%d) too high! set to %d\n",pTopLJetBlock->NJets(),NJetMax);
				njets[i] = NJetMax;
			}
		}

		NJetHist[i]->Fill(njets[i]);		

		//
		// loop over jets 
		for (Int_t j=0; j < njets[i]; j++) {

			if (i==0)JET = pJetBlock->Jet(j);
			if (i==1)JET = pTopTJetBlock->Jet(j);
			if (i==2)JET = pTopLJetBlock->Jet(j);
			
			Mo = JET->Momentum();
			//printf("jet et %f det eta %f emf %f \n", JET->Et(), JET->DetEta(),JET->EmFract());

//
// Fill the histograms if requested
			if (isEnableJetHist()) {

				JetEtHist[i]->Fill(JET->Et());
				JetDtetaHist[i]->Fill(JET->DetEta());
				JetEmfHist[i]->Fill(JET->EmFract());
				JetRapHist[i]->Fill(Mo->PseudoRapidity());
				JetCorfdHist[i]->Fill(JET->Corfd());
				JetCorfmHist[i]->Fill(JET->Corfm());
				JetZvtxHist[i]->Fill(JET->Vz());
				JetMassHist[i]->Fill(JET->M());
				JetPhiHist[i]->Fill( Phi_0_2pi(Mo->Phi()) );
				JetSeedEtaHist[i]->Fill(JET->SeedIEta());
				JetSeedPhiHist[i]->Fill(JET->SeedIPhi());
				JetNtracksHist[i]->Fill(JET->NTracks());
			
//
// These histos are made for each cut var, requiring the 
// enabled cuts of other vars
				if(
					((!JetDtetaCutEnabled) || (fabs(JET->DetEta()) < JetDtetaCut)) &&
					((!JetEmfCutEnabled)   || (JET->EmFract() < JetEmfCut))
					) JetEtCutHist[i]->Fill(JET->Et());
		
				if(
					((!JetEtCutEnabled) || (JET->Et() > JetEtCut)) &&
					((!JetEmfCutEnabled) || (JET->EmFract() < JetEmfCut))
					) JetDtetaCutHist[i]->Fill(JET->DetEta());
	
				if(
					((!JetEtCutEnabled) || (JET->Et() > JetEtCut)) &&
					((!JetDtetaCutEnabled) || (fabs(JET->DetEta()) < JetDtetaCut))
					) JetEmfCutHist[i]->Fill(JET->EmFract());
					
				if(
					((!JetEtCutEnabled) || (JET->Et() > JetEtCut)) &&
					((!JetEmfCutEnabled) || (JET->EmFract() < JetEmfCut)) &&
					((!JetDtetaCutEnabled) || (fabs(JET->DetEta()) < JetDtetaCut))
					) ncutjets[i]++;
			
			} // isEnableHist()

//
// Make the cuts
			if(
				((!JetEtCutEnabled) || (JET->Et() > JetEtCut)) &&
				((!JetDtetaCutEnabled) || (fabs(JET->DetEta()) < JetDtetaCut)) &&
				((!JetEmfCutEnabled) || (JET->EmFract() < JetEmfCut))
			){
				//printf("jet passed\n");
				NgoodJetList[j] = true;
				foundJet[i] = true;
				if (NgoodJet[i] < NJetMax) {
					goodJet[NgoodJet[i]] = JET;
					NgoodJet[i]++;
				}
			
			} else {
				//printf("jet failed\n");
				NgoodJetList[j] = false;
			
			} // good jets
	  
		} // loop over jets
		
		NJetCutHist[i]->Fill(ncutjets[i]);
		
		if (foundJet[i]) jetEvents[i]++;
		if (ientry%1000 == 0)
		printf(" In jet block %s, Events processed = %d with good jets = %d \n",desc[i],totEvents[i],jetEvents[i]);

	} // loop over available jet blocks

return 0;		       

} // Event(int ientry)

//==================================
//       EndJob()
//==================================
int TJetCand::EndJob() {
	printf("----- end job: ---- %s\n",GetName());
	for( int i=0; i<3; i++ )
	if(jetblock[i])printf(" In jet block %s, Events processed = %d with %d jets that pass current selection.\n",
	                      desc[i],totEvents[i],jetEvents[i]);
	printcuts();
return 0;
}

//==================================
//       printcuts()
//==================================
void TJetCand::printcuts() {
	printf("\n Current Jet Selection Cuts:\n");
	printf("  Jet Et cut enabled     %d value %f\n", isEnableEtCut(),getJetEtCut());
	printf("  Jet DetEta cut enabled %d value %f\n", isEnableDtetaCut(),getJetDtetaCut());
	printf("  Jet Emf cut enabled    %d value %f\n\n", isEnableEmfCut(),getJetEmfCut());
}