//_____________________________________________________________________________
// CLC analysis module: below are the examples of what one can do in 
// the interactive session
// 
// note that because of the present limitations of the script compiler this
// module need to reside in the current directory
//_____________________________________________________________________________

/*

  ****** to run:

//  TStnAna x("dataval_ar0192ea_clc.root");
  TStnAna x("stnmaker.root");
  gInterpreter->ProcessLine(".L /data12/murat/anyes/./TMcAnaModule.cc+");
  x.AddModule(new TMcAnaModule());
					// using version 1 of the calib const
					// default is 0

  TMcAnaModule* m = (TMcAnaModule*) x.GetModule("McAna");
  x.Run();

..........

  x.DeleteModule("McAna");
  gInterpreter->ProcessLine(".U /data12/murat/anyes/./TMcAnaModule.so");


					// to fit residual miscalibration:
  TMcAnaModule* m = (TMcAnaModule*) x.GetModule("McAna"); 

  gStyle->SetStatW(0.13);
  gStyle->SetStatH(0.12);
  gStyle->SetStatX(0.4);
  gStyle->SetStatY(0.5);

  ****** to display a histogram

  TMcAnaModule* m = (TMcAnaModule*) x.GetModule("McAna"); 
  TMcAnaModule::Hist_t* h;
  h = m->GetHist();
  h->fPt->Draw();

  ****** to process one event (#15 in the ntuple) and to print the CLC data"

  x.ProcessEvent(15)
  TMcAnaModule* m = (TMcAnaModule*) x.GetModule("ClcAna2");
  TClcData* data = m->GetClcBlock();
  data->Print();

*/

#include "TF1.h"
#include "TPad.h"
#include "TCanvas.h"

#include "TMcAnaModule.hh"

//_____________________________________________________________________________
TMcAnaModule::TMcAnaModule(const char* name, const char* title):
  TStnModule(name,title)
{
}

//_____________________________________________________________________________
TMcAnaModule::~TMcAnaModule() {
}

//_____________________________________________________________________________
void TMcAnaModule::BookHistograms() {
  char name [200];
  char title[200];
				// book histograms
  Delete("hist");

  sprintf(name, "%s_nparticles",GetName());
  sprintf(title,"%s: N(Particles)",GetName());
  fHist.fNParticles = new TH1F(name, title,100 ,0,1000);
  AddHistogram(fHist.fNParticles);

  sprintf(name, "%s_pt",GetName());
  sprintf(title,"%s: Particle Pt",GetName());
  fHist.fPt = new TH1F(name, title,100 ,0,100);
  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);
  AddHistogram(fHist.fParent);

  sprintf(name, "%s_eta",GetName());
  sprintf(title,"%s: Particle Eta",GetName());
  fHist.fEta = new TH1F(name, title,100 ,-5,5);
  AddHistogram(fHist.fEta);
}


//_____________________________________________________________________________
int TMcAnaModule::BeginJob() {

				// register the data block

  fGenpBlock = (TGenpBlock*) RegisterDataBlock("GenpBlock","TGenpBlock");

					// book histograms
  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);
  return 0;
}

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

//_____________________________________________________________________________
int TMcAnaModule::Event(int ientry) {
					// in case of a TChain, ientry is the 
					// entry number in the current file
					// read the trigger block to look at
					// the myron flag
  fGenpBlock->GetEntry(ientry);
					// do whatever you want
  if (PrintLevel()) {
    fGenpBlock->Print();
  }

  fHist.fNParticles->Fill(fGenpBlock->NParticles());

  for (int i=0; i<fGenpBlock->NParticles(); i++) {
    TGenParticle* p = fGenpBlock->Particle(i);
    int pdg_code = p->GetPdgCode();
					// fill the final state electrons and
					// muons only
    if (p->GetStatusCode() == 1) {
      if (p->IsElectron() || p->IsMuon()) {

	fHist.fPt->Fill(p->Pt());
	fHist.fEta->Fill(p->Eta());
					// `im' is a C++ index
	int im = p->GetFirstMother();
	TGenParticle* mom = fGenpBlock->Particle(im);

					// treat PYTHIA case when a particle
					// sometimes goes to itself - we want
					// to see the real parent
	if (mom) {
	  while (mom->GetPdgCode() == pdg_code) {
	    im = mom->GetFirstMother();
	    mom = fGenpBlock->Particle(im);
	  }
	  fHist.fParent->Fill(TMath::Abs(mom->GetPdgCode()));
	}
      }
    }
  }
  return 0;			       
}
//_____________________________________________________________________________
int TMcAnaModule::EndRun() {
  return 0;
}

//_____________________________________________________________________________
int TMcAnaModule::EndJob() {
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}

//_____________________________________________________________________________
void TMcAnaModule::Print(Option_t* opt) const {
}


//_____________________________________________________________________________
void TMcAnaModule::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);

  if (slide == 1) {
//-----------------------------------------------------------------------------
//  basic MC plots
//-----------------------------------------------------------------------------
    //    TPostScript ps("clc_adc.ps",-111);

    TCanvas* c = NewSlide(name,title,2,3);
    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();
  }
  else if (slide == 2) {
				// plot your own histograms ...
  }
}