//==================================================================
//
//	TTopHepg.cc
//
// This code TEMPLATE reads the stntuple "TGenpBlock" (filled from HEPG bank)
// and searches for and plots simple t - W - b particle properties from
// ttbar decay. 
//
//	Version 0.0:      P. Koehn    03/04/02
//
//==================================================================
#include "Stntuple/loop/TStnAna.hh"
#include "TTopHepg.hh"

//=============================================
// TTopHepg(const char* name, const char* title)
//=============================================
TTopHepg::TTopHepg(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("TTopHepg.root","RECREATE","test file");
   histDir = name;
   dir = histFile->mkdir(histDir);
	
}
//==================================================================
//TTopHepg(TFile* file, const char* name, const char* title):
//==================================================================
TTopHepg::TTopHepg(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);
	
	pdg_code = 0;
	im = 0;
	bm = 0;
	wm = 0;
	totev = 0;
	
}

//=============================================
// ~TTopHepg()
//=============================================
TTopHepg::~TTopHepg() {
  fHistogramList->Delete();
  delete fHistogramList;
}

//=============================================
// BookHistograms() 
//=============================================
void TTopHepg::BookHistograms() {
  char name [200];
  char title[200];
  GetHistogramList()->Delete();

  sprintf(name, "%s_nparticles",GetName());
  sprintf(title,"%s: N(Particles)",GetName());
  fHist.fNParticles = new TH1F(name, title,100 ,0,1000);
  fHist.fNParticles->SetXTitle("particle type"); fHist.fNParticles->SetYTitle("number of particles");
  AddHistogram(fHist.fNParticles);

  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: Particle Pt",GetName());
  fHist.fPt = new TH1F(name, title,150 ,0,300);
  fHist.fPt->SetXTitle("Pt (GeV/c)"); fHist.fPt ->SetYTitle("number of particles");
  AddHistogram(fHist.fPt);

  sprintf(name, "%s_parent",GetName());
  sprintf(title,"%s: Particle Parent (abs value)",GetName());
  fHist.fParent = new TH1F(name, title,1000 ,0,1000);
  fHist.fParent->SetXTitle("parent particle type"); fHist.fParent->SetYTitle("number of particles");
  AddHistogram(fHist.fParent);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: Particle Eta",GetName());
  fHist.fEta = new TH1F(name, title,120 ,-6,6);
  fHist.fEta->SetXTitle(" eta "); fHist.fEta ->SetYTitle("number of particles");
  AddHistogram(fHist.fEta);
  
  // top
  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: top Pt",GetName());
  fHist.tPt = new TH1F(name, title,150 ,0,300);
  fHist.tPt->SetXTitle("Pt (GeV/c)"); fHist.tPt ->SetYTitle("number of particles");
  AddHistogram(fHist.tPt);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: top Eta",GetName());
  fHist.tEta = new TH1F(name, title,120 ,-6,6);
  fHist.tEta->SetXTitle(" eta "); fHist.tEta ->SetYTitle("number of particles");
  AddHistogram(fHist.tEta);
  
  // bottom
  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: bottom Pt",GetName());
  fHist.bPt = new TH1F(name, title,150 ,0,300);
  fHist.bPt->SetXTitle("Pt (GeV/c)"); fHist.bPt ->SetYTitle("number of particles");
  AddHistogram(fHist.bPt);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: bottom Eta",GetName());
  fHist.bEta = new TH1F(name, title,120 ,-6,6);
  fHist.bEta->SetXTitle(" eta "); fHist.bEta ->SetYTitle("number of particles");
  AddHistogram(fHist.bEta);
  
  // W boson
  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: W Pt",GetName());
  fHist.wPt = new TH1F(name, title,150 ,0,300);
  fHist.wPt->SetXTitle("Pt (GeV/c)"); fHist.wPt ->SetYTitle("number of particles");
  AddHistogram(fHist.wPt);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: W Eta",GetName());
  fHist.wEta = new TH1F(name, title,120 ,-6,6);
  fHist.wEta->SetXTitle(" eta "); fHist.wEta ->SetYTitle("number of particles");
  AddHistogram(fHist.wEta);
  
  // lepton from W
  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: Lepton Pt",GetName());
  fHist.lPt = new TH1F(name, title,150 ,0,300);
  fHist.lPt->SetXTitle("Pt (GeV/c)"); fHist.lPt ->SetYTitle("number of particles");
  AddHistogram(fHist.lPt);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: lepton Eta",GetName());
  fHist.lEta = new TH1F(name, title,120 ,-6,6);
  fHist.lEta->SetXTitle(" eta "); fHist.lEta ->SetYTitle("number of particles");
  AddHistogram(fHist.lEta);
  
}


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

// Register data block
  fGenpBlock = (TGenpBlock*) RegisterDataBlock("GenpBlock","TGenpBlock");

  BookHistograms();
  
// add particles presently missing in the  database
  TDatabasePDG* db = TDatabasePDG::Instance();

  db->AddParticle("K*0(1410)"    ,"K*(0) 1410"     ,1.410, kTRUE, -1, 0.,"K-meson", 100313);
  db->AddParticle("K*0_bar(1410)","anti K*(0) 1410",1.410, kTRUE, -1, 0.,"K-meson",-100313);

  db->AddParticle("K*+(1410)","K*(+) 1410",1.410, kTRUE, -1, 0.,"K-meson", 100323);
  db->AddParticle("K*-(1410)","K*(-) 1410",1.410, kTRUE, -1, 0.,"K-meson",-100323);

	pdg_code = 0;
	im = 0;
	bm = 0;
	wm = 0;
	totev = 0;

  return 0;
}

//=============================================
// BeginRun() 
//=============================================
int TTopHepg::BeginRun() {
  // initialize run-dependent calibration constants
  
  return 0;
}

//=============================================
// Event(int ientry)
//=============================================
int TTopHepg::Event(int ientry) {
	fGenpBlock->GetEntry(ientry);
	
	totev++;
	
	if (fDebugLevel)fGenpBlock->Print();

	fHist.fNParticles->Fill(fGenpBlock->NParticles());
	
	for (int i=0; i<fGenpBlock->NParticles(); i++) {
		TGenParticle* p = fGenpBlock->Particle(i);
		pdg_code = p->GetPdgCode();
				
		//
		// look at all stable particles
		//
		if (p->GetStatusCode() == 1) {  // status code == 1, final state particles

			fHist.fPt->Fill(p->Pt());
			fHist.fEta->Fill(p->Eta());

			im = p->GetFirstMother();
			TGenParticle* mom = fGenpBlock->Particle(im);
			
			//
 			// This while-loop accounts for some confusing listings in the Hepg bank as filled by pythia:
			// 		1) a particle referring to itself as a mother
			// 		2) initial state colliding p and pbar listed as stable particles, 
			// 			and referring to themselves as a mother
			//
			if ( mom ) {
				while ( (mom->GetPdgCode() == pdg_code) && (im !=0 ) && (im != 1) ) { 
					im = mom->GetFirstMother();
					mom = fGenpBlock->Particle(im);
				//printf("\n INSIDE while: ith part = %d, TGenParticle* mom = %d , pdg_code = %d, mom pdg code = %d, mom index = %d \n", i,mom, pdg_code, mom->GetPdgCode(), im);
				}
				fHist.fParent->Fill(abs(mom->GetPdgCode()));
			}
			 				
		} // look at all stable particles
		
		//
		// find the bottom quarks from top decay
		//
		if (p->GetStatusCode() == 3) {  // status code == 3, particles before fragmentation		
			if(abs(pdg_code)==5){        // find the bottom quarks
				bm = p->GetFirstMother();
				TGenParticle* bmom = fGenpBlock->Particle(bm);
				// treat PYTHIA case when a particle sometimes goes to itself - we want
		   	// to see the real parent - top
				if (bmom) {              // find the bottom quark's mother
					while (bmom->GetPdgCode() == pdg_code) {
						bm = bmom->GetFirstMother();
						bmom = fGenpBlock->Particle(bm);
					}
					// top properties
					if(abs(bmom->GetPdgCode())==6){ //  is bottom quark's mother top ?
				   	//printf("\n --->b's mother pdg_code = %d, index = %d, status = %d",bmom->GetPdgCode(),bm, bmom->GetStatusCode());
						//printf("\n     |--->b particle is pdg_code = %d, index = %d, and status = %d",pdg_code, i,p->GetStatusCode());
						fHist.bPt->Fill(p->Pt());
						fHist.bEta->Fill(p->Eta());
						fHist.tPt->Fill(bmom->Pt());
						fHist.tEta->Fill(bmom->Eta());
					} // top
					
				} // bmother
				
			} // bottom
			
		} // before fragmentation
		
		//
		// find the leptons from the W's from top decay
		//
		if (p->GetStatusCode() == 1) {  // status code == 1, final state particles
			fHist.fPt->Fill(p->Pt());
			fHist.fEta->Fill(p->Eta());
			if ((abs(pdg_code) == 11) || (abs(pdg_code) == 13) ) { // get e or mu
				im = p->GetFirstMother();
				TGenParticle* mom = fGenpBlock->Particle(im);
				if (mom) {                         // find lepton mother
					while (mom->GetPdgCode() == pdg_code) {
						im = mom->GetFirstMother();
						mom = fGenpBlock->Particle(im);
					}
				   if(abs(mom->GetPdgCode())==24){ // if it is a DubYa...					
						wm = mom->GetFirstMother();
						TGenParticle* wmom = fGenpBlock->Particle(wm);
						if (wmom) {                  // find DubYa's mom
							while (wmom->GetPdgCode() == pdg_code) {
								wm = wmom->GetFirstMother();
								wmom = fGenpBlock->Particle(wm);
							}
							if(abs(wmom->GetPdgCode())==6){ // if DubYa's mom is top
								//printf("\n --->LEPTON's mother pdg_code = %d, index = %d , and status = %d",mom->GetPdgCode(),im,mom->GetStatusCode());
				      		//printf("\n        |--->LEPTON:  pdg_code = %d, index = %d, and status = %d",pdg_code,i,p->GetStatusCode());
								fHist.wPt->Fill(mom->Pt());
								fHist.wEta->Fill(mom->Eta());
								fHist.lPt->Fill(p->Pt());
								fHist.lEta->Fill(p->Eta());
							} // if DubYa's mom is top					
						} // find DubYa's mom					
					} // if it is a DubYa
				} // find lepton mother				
			} // get stable leptons						
		} // status code == 1, final state particles		
				
	} // loop over all particles
	
	
	//if (totev%500==0) printf("\n Events processed so far = %d\n",totev);
	
	return 0;			       
} // loop over Events

//=============================================
// EndRun()
//=============================================
//int TTopHepg::EndRun() {
//  return 0;
//}

//=============================================
// EndJob()
//=============================================
int TTopHepg::EndJob() {

  printf("\n Total events processed = %d\n",totev);

  printf("\n----- end job: ---- %s\n",GetName());
  
  return 0;
}

//=============================================
// spew()
//=============================================
void TTopHepg::spew() {
	printf("\n urp...slop...get the mop!");
}


//=============================================
// PlotHistograms(int run_number, int slide)
//=============================================
void TTopHepg::PlotHistograms(int run_number, int slide) {

  // plot slides
  char name[120], canvas_name[120], title[120], pad_name[120];
  sprintf(name,"run_%i_slide_%i",run_number,slide);
  sprintf(title,"run %i slide %i ",run_number, slide);
  sprintf(canvas_name,"%s_%s",GetName(),name);
  sprintf(pad_name,"%s_p1",canvas_name);

  //
  //  all particles
  //
  if (slide == 1) { 

    TCanvas* c = NewSlide(name,title,2,2);
    TPad* p1 = (TPad*) c->GetPrimitive(pad_name);

    p1->cd(1); 
    fHist.fNParticles->Draw();

    p1->cd(2); 
    fHist.fPt->Draw();

    p1->cd(3); 
    fHist.fEta->Draw();

    p1->cd(4); 
    fHist.fParent->Draw();

    gPad->Update();
	 
  //
  //  top and bottom
  //
  } else if (slide == 2) { 
    TCanvas* c2 = NewSlide(name,title,2,2);
    TPad* p2 = (TPad*) c2->GetPrimitive(pad_name);

    p2->cd(1); 
    fHist.tPt->Draw();

    p2->cd(2); 
    fHist.tEta->Draw();

    p2->cd(3); 
    fHist.bPt->Draw();

    p2->cd(4); 
    fHist.bEta->Draw();

    gPad->Update();
				
  //
  //  Dubya and lepton(electron or muon)
  //
  } else if (slide == 3) { 
    TCanvas* c3 = NewSlide(name,title,2,2);
    TPad* p3 = (TPad*) c3->GetPrimitive(pad_name);

    p3->cd(1); 
    fHist.wPt->Draw();

    p3->cd(2); 
    fHist.wEta->Draw();

    p3->cd(3); 
    fHist.lPt->Draw();

    p3->cd(4); 
    fHist.lEta->Draw();

    gPad->Update();
				
  }
  
}