//#====================================================#
//#        Object Oriented Analysis Class              #
//# read the tutorial ooanalysis.html in TopMods/doc   #
//#====================================================#

#include <iostream>
#include "TChain.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"

using std::cout;
using std::endl;

// defines objects like jet, electron, offline track...
#include "HighLevelObjects.hh"

// needed stuff in my analysis
#include "anautil.c"
float zn[3]={0.0, 0.0, 1.0};  // unit vector of z-direction

class ana {
public:

  // chain = list of root files containing the same tree
  TChain* chain;
  
  // The output file with histograms
  TFile* histf;
    
  // The tree branches
  //
  // Arrays of objects per event
  // Make sure the names don't coinside with the names used
  // in HighLevelObjects.cc
  TClonesArray* evtA;    //event
  TClonesArray* smryA;   //summary
  TClonesArray* pvtxA;   //primary vertex 
  TClonesArray* svtxtrA;  //secondary vertex track
  TClonesArray* jpbtrA;  //JetProb tracks
  TClonesArray* jpbassocA;  //JetProb track-jet association
  TClonesArray* trinfA;  //trigger info
  TClonesArray* trnmA;   //trigger name
  TClonesArray* metA;    //missing Et
  TClonesArray* hepgA;   //hepg bank
  TClonesArray* obspA;   //obsp bank
  TClonesArray* eleA;    //electrons
  TClonesArray* muoA;    //muons
  TClonesArray* jetA;    //jets
  //    TClonesArray* jetcA;   //corrected jets
  TClonesArray* otrkA;   //offline tracks
  TClonesArray* chitA;  //cot hits
  TClonesArray* xhitA;   //xft hits
  TClonesArray* xpxlA;   //xft pixels
  TClonesArray* xtrkA;   //xft tracks
  //
  // number of objects in each event (and each TClonesArray)
  int nsvtr;
  int njpbtr;
  int njpbassoc;
  int nmet;
  int nhepg;
  int nobsp;  
  int nele;
  int nmuo;
  int njet;
  int njetc; 
  int notrk;
  int nchit;
  int nxhit;
  int nxpxl;
  int nxtrk;
  //
  // Collection/View descriptions
  //char metCollDescr[40];
  // 
  char otrkCollDescr1[40];
  char otrkCollDescr2[40];
  char otrkCollDescr3[40];
  char otrkCollDescr4[40];
  char otrkCollDescr5[40]; 
  //
  char muoViewDescr1[40];
  char muoViewDescr2[40];
  char muoViewDescr3[40];
  char muoViewDescr4[40];
  char muoViewDescr5[40]; 
  //
  char eleViewDescr1[40];
  char eleViewDescr2[40];
  char eleViewDescr3[40];
  char eleViewDescr4[40];
  char eleViewDescr5[40]; 
  //
  char jetCollDescr1[40];
  char jetCollDescr2[40];
  char jetCollDescr3[40];
  char jetCollDescr4[40];
  char jetCollDescr5[40];
  //
  char svtrkCollDescr[40];  
          
  ana();
  ~ana(){};
    
  void LoadFile(); // To load file
  void OpenHistFile(const char*);
  void CloseHistFile();
  void EventLoop();// the main analysis 
    
  int CheckBranches(int entry);// checks if the branches exist
  int currtreenum;// current tree number in the chain
    
  // Do switches - use them instead of commenting out big chunks of code.
  // 1 - process that branch, 0 - don't process it.
  // A switch will be set to 0 by CheckBranches() if the corresponding
  // branch doesn't exist in the current file!!!
  bool      doevt,
    dosmry ,
    dopvtx ,
    dosvtxtr,
    dojpbtr,
    dojpbassoc,
    dotrinf,
    dotrnm ,
    domet  ,
    dohepg ,    
    doobsp ,
    doele  ,
    domuo  ,     
    dojet  ,     
    dojetc ,
    dootrk ,
    dochit , 
    doxhit ,    
    doxpxl ,    
    doxtrk ;  
  void SetSwitches();             
    
};

void ana::OpenHistFile(const char* filename)
{
  // Open/create the analysis histogram file
  histf = new TFile(filename,"RECREATE");
}
void ana::CloseHistFile()
{
  // Close the histogram file
  histf -> Write();
  histf -> Close();
}

ana::ana(){
  cout<<endl;
   
  // Create the chain and add the root file(s) to it
  chain = new TChain("TopTree");
    
  // Initialize the clones arrays
  evtA = new TClonesArray("evt");
  smryA = new TClonesArray("summary");   
  pvtxA = new TClonesArray("privertex");
  svtxtrA = new TClonesArray("secvtxtrack");  
  jpbtrA = new TClonesArray("jetprobtrack");  
  jpbassocA = new TClonesArray("jetprobTJassoc");  
  trinfA = new TClonesArray("trigInfo");  
  trnmA = new TClonesArray("trigName");   
  metA = new TClonesArray("missingEt");   
  hepgA = new TClonesArray("hepg");   
  obspA = new TClonesArray("obsp");   
  eleA = new TClonesArray("electron");    
  muoA = new TClonesArray("muon");    
  jetA = new TClonesArray("jet");    
  //    jetcA = new TClonesArray("jetCorr");  
  otrkA = new TClonesArray("offltrack"); 
  chitA = new TClonesArray("cothit");  
  xhitA = new TClonesArray("xfthit");   
  xpxlA = new TClonesArray("xftpixel");   
  xtrkA = new TClonesArray("xfttrack"); 
    
  // Connect the branches
  chain->SetBranchAddress("evt",&evtA);
  chain->SetBranchAddress("summary",&smryA);
  chain->SetBranchAddress("privertex",&pvtxA);
  chain->SetBranchAddress("secvtxtrack",&svtxtrA);
  chain->SetBranchAddress("jetprobtrack",&jpbtrA);
  chain->SetBranchAddress("jetprobTJassoc",&jpbassocA);
  chain->SetBranchAddress("trigInfo",&trinfA);
  chain->SetBranchAddress("trigName",&trnmA);
  chain->SetBranchAddress("met",&metA);
  chain->SetBranchAddress("hepg",&hepgA);
  chain->SetBranchAddress("obsp",&obspA);  
  chain->SetBranchAddress("electron",&eleA);
  chain->SetBranchAddress("muon",&muoA);
  chain->SetBranchAddress("jet",&jetA);
  //    chain->SetBranchAddress("jetCorr",&jetcA);
  chain->SetBranchAddress("offltrack",&otrkA);
  chain->SetBranchAddress("cothit",&chitA);
  chain->SetBranchAddress("xfthit",&xhitA);
  chain->SetBranchAddress("xftpixel",&xpxlA);
  chain->SetBranchAddress("xfttrack",&xtrkA);
  //
  chain->SetBranchAddress("numSecvtxTracks",&nsvtr);
  chain->SetBranchAddress("numJetProbTracks",&njpbtr);
  chain->SetBranchAddress("numJetProbTJAssoc",&njpbassoc);
  chain->SetBranchAddress("numMetObjs",&nmet);
  chain->SetBranchAddress("numHepgObjs",&nhepg);
  chain->SetBranchAddress("numObspObjs",&nobsp);
  chain->SetBranchAddress("numEmObjs",&nele);  
  chain->SetBranchAddress("numMuObjs",&nmuo);
  chain->SetBranchAddress("numJetObjs",&njet);
  chain->SetBranchAddress("numJetCorrObjs",&njetc);
  chain->SetBranchAddress("numOfflTracks",&notrk); 
  chain->SetBranchAddress("numCotHits",&nchit);    
  chain->SetBranchAddress("numXFTHits",&nxhit);
  chain->SetBranchAddress("numXFTPixels",&nxpxl);
  chain->SetBranchAddress("numXFTTracks",&nxtrk);
  //
  //chain->SetBranchAddress("metCollDescr",metCollDescr); 
  chain->SetBranchAddress("otrkCollDescr1",otrkCollDescr1);
  chain->SetBranchAddress("otrkCollDescr2",otrkCollDescr2);
  chain->SetBranchAddress("otrkCollDescr3",otrkCollDescr3);
  chain->SetBranchAddress("otrkCollDescr4",otrkCollDescr4); 
  chain->SetBranchAddress("otrkCollDescr5",otrkCollDescr5);
  chain->SetBranchAddress("muoViewDescr1",muoViewDescr1);
  chain->SetBranchAddress("muoViewDescr2",muoViewDescr2);
  chain->SetBranchAddress("muoViewDescr3",muoViewDescr3);       
  chain->SetBranchAddress("muoViewDescr4",muoViewDescr4);
  chain->SetBranchAddress("muoViewDescr5",muoViewDescr5);
  chain->SetBranchAddress("eleViewDescr1",eleViewDescr1);
  chain->SetBranchAddress("eleViewDescr2",eleViewDescr2); 
  chain->SetBranchAddress("eleViewDescr3",eleViewDescr3);
  chain->SetBranchAddress("eleViewDescr4",eleViewDescr4);
  chain->SetBranchAddress("eleViewDescr5",eleViewDescr5);
  chain->SetBranchAddress("jetCollDescr1",jetCollDescr1); 
  chain->SetBranchAddress("jetCollDescr2",jetCollDescr2);
  chain->SetBranchAddress("jetCollDescr3",jetCollDescr3);
  chain->SetBranchAddress("jetCollDescr4",jetCollDescr4);
  chain->SetBranchAddress("jetCollDescr5",jetCollDescr5);
  chain->SetBranchAddress("svtrkCollDescr",svtrkCollDescr);
}

void ana::SetSwitches(){
  // Do switches - use them instead of commenting out big chunks of code.
  // 1 - process that branch, 0 - don't process it.
  // A switch will be set to 0 by CheckBranches() if the corresponding
  // branch doesn't exist in the current file!!!
    
  doevt=   1;
  dosmry = 1;
  dopvtx = 1;
  dosvtxtr=1;
  dojpbtr =1;
  dojpbassoc=1;
  dotrinf= 0;
  dotrnm = 0;
  domet  = 0;
  dohepg = 1;    
  doobsp = 0;
  doele  = 0;
  domuo  = 0;     
  dojet  = 1;     
  dojetc = 0;
  dootrk = 1;
  dochit = 1; 
  doxhit = 0;    
  doxpxl = 0;    
  doxtrk = 1; 
}   
  
void ana::EventLoop(){
  // Initialize the histograms
  TH1F *evt_evnum   = new TH1F("evt_evnum","evt_evnum",200,0,8000);
  TH1F *smry_nTightJet   = new TH1F("smry_nTightJet","smry_nTightJet",
				    100, 0, 100);
  TH1F *pvtx_x   = new TH1F("pvtx_x","pvtx_x",200,-1050,50);
  TH1F *jet_Pt   = new TH1F("jet_Pt","jet_Pt",200,0,20);            
  TH1F *otrk_Pt   = new TH1F("otrk_Pt","otrk_Pt",200,0,20);
  TH1F *xtrk_Pt   = new TH1F("xtrk_Pt","xtrk_Pt",200,0,20);  
  TH2F *angdist = new TH2F("angdist","Angular dist.",
			   6, -1.0, 1.0, 6, -1.0, 1.0);  
  
  // Get the number of events/entries in the file chain
  int nEvents = chain -> GetEntries(); 
  cout << "Entry=" << nEvents << endl;
  
  for(int ev=0; ev<nEvents; ev++) { // loop over the events
    //for(int ev=0; ev<10; ev++) { // Just ONE time
    
    if ((ev+1)%100==0) { cout << "index=" << ev+1 << endl;}
    //cout << "index=" << ev+1 << endl;

    CheckBranches(ev);//disable processing of branches that don't exist
    chain -> GetEntry(ev);
    
    if (doevt) {// event
      evt* event = (evt*) evtA->At(0); 
      evt_evnum -> Fill(event->eventNumber);
    }// event
    
    if (dosmry) {// summary
      summary* smry = (summary*) smryA->At(0);
      smry_nTightJet -> Fill(smry->fnTightJet);
    }// summary
    
    if (dopvtx) {// primary vertex
      privertex* pvtx = (privertex*) pvtxA->At(0);
      pvtx_x -> Fill(pvtx->x);
    }// primary vertex
    
    //////////////////
    // hepg bank
    ////////////////////
    if (dohepg) {// hepg

      // define lorentz vectors for parton level
      float t1l[4], b1l[4], l1l[4], n1l[4];
      float t2l[4], b2l[4], l2l[4], n2l[4];
      float b1t[4], l1t[4], n1t[4];
      float b2t[4], l2t[4], n2t[4];

      // main loop in a event
      for(int indx = 0; indx<nhepg; indx++) {
	hepg* hg = (hepg*) hepgA->At(indx);
	bool cond=false;
	if (hg->ID==6 && hg->Stdhep==3) { // top
	  t1l[0]=hg->Px;
	  t1l[1]=hg->Py;
	  t1l[2]=hg->Pz;
	  t1l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==-6 && hg->Stdhep==3) { // anti-top
	  t2l[0]=hg->Px;
	  t2l[1]=hg->Py;
	  t2l[2]=hg->Pz;
	  t2l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==24 && hg->IDparent==6) { // W+ from top
	  cond=true;
	}
	if (hg->ID==-24 && hg->IDparent==-6) { // W- from anti-top
	  cond=true;
	}
	if (hg->ID==5 && hg->IDparent==6) { // b from top
	  b1l[0]=hg->Px;
	  b1l[1]=hg->Py;
	  b1l[2]=hg->Pz;
	  b1l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==-5 && hg->IDparent==-6) { // bbar from anti-top
	  b2l[0]=hg->Px;
	  b2l[1]=hg->Py;
	  b2l[2]=hg->Pz;
	  b2l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==-11 && hg->IDparent==24) { // e+ from W+
	  l1l[0]=hg->Px;
	  l1l[1]=hg->Py;
	  l1l[2]=hg->Pz;
	  l1l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==11 && hg->IDparent==-24) { // e- from W-
	  l2l[0]=hg->Px;
	  l2l[1]=hg->Py;
	  l2l[2]=hg->Pz;
	  l2l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==12 && hg->IDparent==24) { // nu from W+
	  n1l[0]=hg->Px;
	  n1l[1]=hg->Py;
	  n1l[2]=hg->Pz;
	  n1l[3]=hg->E;
	  cond=true;
	}
	if (hg->ID==-12 && hg->IDparent==-24) { // nu- from W-
	  n2l[0]=hg->Px;
	  n2l[1]=hg->Py;
	  n2l[2]=hg->Pz;
	  n2l[3]=hg->E;
	  cond=true;
	}
	if (0 && cond) {
	  cout << hg->Ind << ": " << hg->ID << ": " << hg->IDparent << ": "
	       << hg->Stdhep << ": " << hg->Mo1 << ": ("
	       << hg->Px << ", " << hg->Py << ", " << hg->Pz << "; "
	       << hg->E << ") " << hg->Mass <<endl;
	}
      } // END of loop
      loren4(t1l, l1l, l1t);
      loren4(t2l, l2l, l2t);
      angdist->Fill(cos12(l1t, zn), cos12(l2t, zn));
    }// jet

    if (dojet) {// jet
      for(int indx = 0; indx<njet; indx++) {
	jet* jt = (jet*) jetA->At(indx);
	jet_Pt -> Fill(jt->Pt);
      }
    }// jet
    
    if (dootrk) {// offline tracks
      for(int itr=0; itr<notrk; itr++) {
	offltrack* otrk = (offltrack*) otrkA->At(itr);
	otrk_Pt -> Fill(otrk->Pt);
      } 
    }// offline tracks
    
    if (doxtrk) {// xft tracks
      for(int itr=0; itr<nxtrk; itr++) {
	xfttrack* xtrk = (xfttrack*) xtrkA->At(itr);
	xtrk_Pt -> Fill(xtrk->ptVal);
      } 
    }// xft tracks
    
  }// loop over the events
}


int ana::CheckBranches(int entry) {
  // Must be called before each GetEntry.
  // Checks if the branches exist if a new root 
  // file from the chain is loaded.
  // If a branch doesn't exist it's not processed
  // regardless of the corresponding switch in SetSwitches()!

  if (!chain) return -5;
  // entry is the event number in the whole chain
  // centry is the local event number in the current file
  int centry = chain->LoadTree(entry);
  if (centry < 0) return centry;
  if (chain->IsA() != TChain::Class()) return centry;
  if (chain->GetTreeNumber() != currtreenum) { // the next file has been loaded

    currtreenum = chain->GetTreeNumber(); 	
    SetSwitches(); // set the branches the user requested for processing
	
    // Check if the branches exist in the new file
    // and disable the processing of the ones that don't exist.
    //
    TBranch *pbr;	
    cout<<"/\\                                                        /\\"<<endl;
    cout<<"||                                                        ||"<<endl;
    cout<< "processing file " << currtreenum<<
      ";if missing branches-you'll see 'unknown branch' above."<<endl<<endl<<endl;
    //
    pbr = chain->GetBranch("evt");              if (pbr==0) doevt =0;  
    pbr = chain->GetBranch("summary");          if (pbr==0) dosmry =0; 
    pbr = chain->GetBranch("privertex");        if (pbr==0) dopvtx =0;
    pbr = chain->GetBranch("secvtxtrack");      if (pbr==0) dosvtxtr =0;  
    pbr = chain->GetBranch("trigInfo");         if (pbr==0) dotrinf =0;  
    pbr = chain->GetBranch("trigName");         if (pbr==0) dotrnm =0;  
    pbr = chain->GetBranch("met");              if (pbr==0) domet =0; 
    pbr = chain->GetBranch("hepg");             if (pbr==0) dohepg =0; 
    pbr = chain->GetBranch("obsp");             if (pbr==0) doobsp =0;  
    pbr = chain->GetBranch("numEmObjs");        if (pbr==0) doele =0;  
    pbr = chain->GetBranch("numMuObjs");        if (pbr==0) domuo =0;  
    pbr = chain->GetBranch("numJetObjs");       if (pbr==0) dojet =0;  
    pbr = chain->GetBranch("numJetCorrObjs");   if (pbr==0) dojetc =0; 
    pbr = chain->GetBranch("numOfflTracks");    if (pbr==0) dootrk =0;
    pbr = chain->GetBranch("numCotHits");       if (pbr==0) dochit =0;  
    pbr = chain->GetBranch("numXFTHits");       if (pbr==0) doxhit =0;  
    pbr = chain->GetBranch("numXFTPixels");     if (pbr==0) doxpxl =0;  
    pbr = chain->GetBranch("numXFTTracks");     if (pbr==0) doxtrk =0;  
  }
  return centry;
}

#include "filedef.cc"
// EOF ////