#define TMBTree_bu_cxx
#include "TMBTree_bu.h"
#include "tmb_tree/TMBBCJet.hpp"
#include "tmb_tree/TMBBCMuTagger.hpp"
#include "TLorentzVector.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TRefArray.h"
#include "TString.h"
#include <iostream>

#ifdef __MAKECINT__
#pragma link C++ class TMBTree_bu;
#pragma link C++ class std::vector<TH1*>;
#endif

using namespace std;

void TMBTree_bu::Walk( TMBMCpart * mcpart, Bool_t debug )
{
   Int_t pdgid = mcpart->pdgid();
   if (debug) cout << " Next particle's PDG ID " << pdgid << endl; 
   // Track down the particle
   TMBMCvtx  * mcvtx = mcpart->getDMCvtx();
   if (!mcvtx) {
      if (debug) cout << "   particle is stable " << endl; 
      return;
   }
   Int_t ndaughters = mcvtx->ndaughters();
   if (debug) cout << "    particle has " << ndaughters << " daughters" << endl;
   if (ndaughters == 0) return;
   for (Int_t ii=0; ii<ndaughters; ii++) {
      TMBMCpart * mcpartd = mcvtx->getDaughter(ii);
      pdgid = mcpartd->pdgid();
      if (debug) 
	 cout << "      PDG ID of " << ii << "th daughter " << pdgid << endl;
      Walk (mcpartd, debug);
   }
}


void TMBTree_bu::Loop()
{
   //   In a ROOT session, you can do:
   //      Root > .L T.C+
   //      Root > T 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:
   // disable all branches:
   //    fChain->SetBranchStatus("*",0);  
   // activate branchname (DON'T FORGET THE "*" AT THE END!)
   //    fChain->SetBranchStatus("branchname*",1);  

   if (fChain == 0) return;
   
   if (!fH.size()) {
      fH.resize(40);
      fH[0] = 0; // whoever set this up was used to FORTRAN...
      fH[1] = new TH1F("h1","Trks_px:",50,-10.,10.);
      fH[2] = new TH1F("h2","Emcl_px:",50,-10.,10.);
      fH[3] = new TH1F("h3","emcl_px-chtr_px - correct pair",
		       50,-10.,10.);
      fH[4] = new TH1F("h4","emcl_px(all)-chtr_px(1) - all pairs",
		       50,-10.,10.);
      fH[5] = new TH1F("h5","Jets_px:",50,-50.,50.);
      fH[6] = new TH1F("h6","Number of Charged Tracks in a jet all",
		       20,0.,10.);
      fH[7] = new TH1F("h7","Cos(jet sum_chp - correct pair) ",
		       50,-1.,1.);
      fH[8] = new TH1F("h8","Cos(1st_jet sum_chp) ",50,-1.,1.);
      fH[9] = new TH1F("h9","Number of Charged Tracks in a jet algo JCCA",
		       20,0.,10.);
      fH[10] = new TH1F("h10","Number of Charged Tracks in a jet algo JCMG",
			20,0.,10.);
      fH[11] = new TH1F("h11","Number of Charged Tracks in a jet algo JKCA",
			20,0.,10.);
      fH[12] = new TH1F("h12","Number of Charged Tracks in a jet algo JNC1",
			20,0.,10.);
      fH[13] = new TH1F("h13","Missing ET",100,0.,100.);
      fH[14] = new TH1F("h14","Missing ET",100,0.,100.);
      fH[15] = new TH1F("h15","Missing ET",100,0.,100.);
      fH[16] = new TH1F("h16","Ring EMx",20,-10.,10.);
      fH[17] = new TH1F("h17","Jet_vertex_z",100,-100.,100.);
      fH[18] = new TH1F("h18","All_vertex_z",100,-100.,100.);
      fH[19] = new TH1F("h19","MCpart multiplicity",1000,0.,1000.);
      fH[20] = new TH1F("h20","Event Numbers",1000,0.,500.);
      fH[21] = new TH1F("h21","pT of Type 1 Taus",20,0.,80.);
      fH[22] = new TH1F("h22","pT of Type 2 Taus",20,0.,80.);
      fH[23] = new TH2F("h23","Tau tracks E vs Tau E", 100, 0., 100., 100, 0., 100.);
      fH[24] = 0;
      fH[25] = new TH1F("h25","emcl pt  scone ",100, 0.,50.);
      fH[26] = new TH1F("h26","emcl pt CellNN ",100, 0.,50.);
      fH[27] = new TH1F("h27","emcl pt of electrons/positrons ",100, 0.,50.);
      fH[28] = new TH1F("h28","emcl pt of photons ",100, 0.,50.);
      fH[29] = new TH1F("h29","Number of hits on a track ",100, 0, 50);
      fH[30] = new TH1F("h30","err(dca) 5 < Nhit < 10 ",100, 0., 0.15);
      fH[31] = new TH1F("h31","err(dca) 10 <= Nhit < 20 ",100, 0., 0.15);
      fH[32] = new TH1F("h32","err(dca) Nhit >= 20 ",100, 0., 0.15);
      fH[33] = new TH1F("h33","Number of cells/emobj ",100, 0.,100.);
      fH[34] = new TH1F("h34","Total Cell Energy ",100, 0.,100.);
      fH[35] = new TH1F("h35","Cell Energy - EMobj Energy ",100, -5.,5.);
      fH[36] = new TH1F("h36","Masses of charged tracks ",300, -1.,2.);
      
      fH[37] = new TH1F("h37","pt of tagged bjet",100, 0., 50.);
      fH[38] = new TH1F("h38","PtRel of muon wrt tagged bjet",100, 0., 10.);
      fH[39] = new TH1F("h39","Charge of the muon in a tagged bjet",10, -2., 2.);
   }

   Int_t entry=0;
   while (GetEntry(entry++)) {
      Bool_t debug = entry<3;//kFALSE;

      //cerr << std::endl << "event: " << entry << std::endl;
      //cerr << "Size of Jets array: " << fJets->GetLast()+1 << std::endl;
      //cerr << "Size of Trks array: " << fTrks->GetLast()+1 << std::endl;
      //cerr << "Size of Emcl array: " << fEmcl->GetLast()+1 << std::endl;

      //Event number
      fH[20]->Fill( fGlob->evtno() );
      if(debug) std::cout << "evt " << fGlob->evtno() << std::endl ;

      //History
      if (debug) {
	 TIter iHist(fHist);
	 TMBHist* hist=0;
	 while ((hist=(TMBHist*)iHist()))
	    std::cout << hist->program() << "  " << hist->version() << std::endl;
      }

      //Vertices:
      TIter iVrts(fVrts);
      TMBVrts* vrts=0;
      while ((vrts=(TMBVrts*)iVrts()))
	 fH[18]->Fill( vrts->vz() );

      //Tracks:
      TIter iTrks(fTrks);
      TMBTrks* trks=0;
      Int_t i=0;
      while ((trks=(TMBTrks*)iTrks())) {
	 Float_t mass = (trks->lorentz_vector()).M();
	 if (debug) std::cout << "mass of particle " << i << " = " << mass << std::endl;
	 fH[36]->Fill( mass );
	 fH[1]->Fill( trks->px() );
	 Int_t Nhit = 0;
	 if (debug) std::cout << "Track: " << i++ << std::endl;
	 for (Int_t mask=0; mask<3; mask++){
	    UInt_t hitmask = *(trks->hitmask(mask));
	    if (debug) std::cout << "mask: " << mask << " hitmask: " << hitmask << std::endl;
	    for (Int_t bit=0; bit<32; bit++) {
	       if (hitmask & (1 << bit) ) {
		  ++Nhit;
		  if (debug) std::cout << "bit: " << bit << " is set" << std::endl;
	       }
	    } // bit
	 } // mask
	 fH[29]->Fill( Nhit );
	 if ( Nhit > 5 ) {
	    if ( Nhit < 10 ) fH[30]->Fill( trks->dca1err() );
	    if ( Nhit >= 10 && Nhit < 20 ) fH[31]->Fill( trks->dca1err() );
	    if ( Nhit >= 20 ) fH[32]->Fill( trks->dca1err() );
	 }
	 const Double_t* trpars = trks->trpars();
	 const Double_t* trerrs = trks->trerrs();
	 if (debug) {
	    std::cout << std::endl;
	    for (Int_t itr=0; itr<5; itr++)
	       std::cout << "track param " << itr << " = " << *(trpars+itr) << std::endl;
	    std::cout << std::endl;
	    for (Int_t itr=0; itr<15; itr++)
	       std::cout << "track error " << itr << " = " << *(trerrs+itr) << std::endl;
	 }
      }// fTrks
      //Emcl:
      TIter iEmcl(fEmcl);
      TMBEmcl* emcl=0;
      while ((emcl=(TMBEmcl*)iEmcl())) {
	 fH[2]->Fill( emcl->px() );
	 //testing algo name
	 if (emcl->algoname() == TString("A")) std::cout <<"=======ERROR===="<<std::endl;
	 if (emcl->algoname() == TString("scone"))
	    fH[25]->Fill( emcl->pT() );
	 if (emcl->algoname() == TString("CellNN"))
	    fH[26]->Fill( emcl->pT() );

	 if ( emcl->charge() < 0. || emcl->charge() > 0.)
	    fH[27]->Fill( emcl->pT() );
	 else fH[28]->Fill( emcl->pT() );

	 //cerr << "Emcl is at " << emcl << " with " << emcl->px() << std::endl;
	 //cerr << "TRef is at " << &(emcl->_chptr) << std::endl;
	 //emcl->_chptr.Dump();

	 TMBTrks* trks = emcl->GetChargedTrack();
	 //         cerr << "Tracks is at " << trks <<  std::endl;
	 if (trks) fH[3]->Fill( emcl->px()-trks->px() );
	 TMBTrks* trks0=(TMBTrks*) fTrks->At(0);
	 if (trks0) fH[4]->Fill( emcl->px()-trks0->px() );

	 fH[33]->Fill(emcl->ncells());
	 Float_t Sum_E = 0.; 
	 for(Int_t j=0; j<emcl->ncells(); j++)
	    Sum_E += emcl->GetEmCell(j)->E();
	 fH[34]->Fill(Sum_E);
	 fH[35]->Fill(Sum_E-emcl->E());
      } //Emcl

      //Jets:
      TIter iJets(fJets);
      TMBJets* jets=0;
      while ((jets=(TMBJets*) iJets())) {
	 TString algo_name = jets->algoname();
	 //std::cout << "algoname: " << jets->algoname() << std::endl;

	 TMBVrts* vrts = jets->GetVertex();
	 if (vrts) fH[17]->Fill( vrts->vz() );

	 TMBJets* jets0 = (TMBJets*) fJets->At(0);
	 fH[5]->Fill(jets->px());
	 fH[6]->Fill(jets->Ntr());

	 if (jets->algoname() == TString("JCCA"))   
	    fH[9]->Fill(jets->Ntr());
	 if (jets->algoname() == TString("JCMG"))       
	    fH[10]->Fill(jets->Ntr());
	 if (jets->algoname() == TString("JKCA")) 
	    fH[11]->Fill(jets->Ntr());
	 if (jets->algoname() == TString("JNC1"))   
	    fH[12]->Fill(jets->Ntr());

	 //cerr << " nmbr of tracks: " << jets->Ntr() << std::endl; 

	 TVector3 Sum_chp;
	 //cerr << "before jet loop " << std::endl; 

	 for (Int_t j=0; j<jets->Ntr(); j++) {
	    TMBTrks* trk=jets->GetChargedTrack(j);
	    if (trk)
	       Sum_chp+=trk->lorentz_vector().Vect();
	 }

	 //cerr << "Jets coordinates: " << Sum_chp.X() << " " << Sum_chp.Y() << " " << Sum_chp.Z()<< std::endl; 

	 if (jets->Ntr() > 0) { 
	    const TVector3 jetv3=jets->lorentz_vector().Vect();
	    Double_t cos_JetSumchp = Sum_chp.Dot(jetv3);
	    cos_JetSumchp /= Sum_chp.Mag();
	    cos_JetSumchp /= jetv3.Mag();
	    fH[7]->Fill( cos_JetSumchp );

	    //cerr << "histo 7 filled " << std::endl; 

	    const TVector3 jetv30=jets0->lorentz_vector().Vect();
	    cos_JetSumchp = Sum_chp.Dot(jetv30);
	    cos_JetSumchp /= Sum_chp.Mag();
	    cos_JetSumchp /= jetv30.Mag();
	    fH[8]->Fill( cos_JetSumchp );

	    //cerr << "histo 8 filled " << std::endl; 
	 }
      }

      // MET:
      //double vMET = fMet->getMET();
      //double vMET_noeta = fMet->getMET_noeta();
      //double vMET_weta = fMet->getMET_weta();
      //fH[13]->Fill(fMet->getMET());
      fH[14]->Fill(fMet->getMET_noeta());
      fH[15]->Fill(fMet->getMET_weta());
      //fH[13]->Fill(vMET);
      //fH[14]->Fill(vMET_noeta);
      //fH[15]->Fill(vMET_weta);

      const Float_t *ringEMx ;
      const Float_t *ringEMy ;
      const Float_t *ringHDx ;
      const Float_t *ringHDy ;
      fMet->getRings(ringEMx,ringEMy,ringHDx,ringHDy);
      Int_t Nrings = fMet->numRings();
      for (Int_t i=0; i<Nrings; i++)
	 {
	    /*
	      std::cout << i << " EMX EMy HDx HDy rings: " << 
	      ringEMx[i] << " " << ringEMy[i] << " " <<
	      ringHDx[i] << " " << ringHDy[i] << std::endl; 
	    */
	    fH[16]->Fill( ringEMx[i] );  
	 }
      //cerr << "end of event: " << entry << std::endl;
      // MC particles
      Int_t mult_mc=0;
      if (debug) std::cout << "Number of MC particle: " << fMCpart->GetLast()+1 << std::endl;

      TIter iMCpart(fMCpart);
      TMBMCpart* mcpart=0;
      while ((mcpart =(TMBMCpart*) iMCpart())) {
	 if (debug) std::cout << " pointer of " << mult_mc << " particle: " << mcpart << std::endl;
	 Int_t pdgid = mcpart->pdgid();
	 TMBMCvtx* mcvtx = mcpart->getPMCvtx();
	 if (mcvtx) {
	    Int_t nparents = mcvtx->nparents();
	    debug = nparents > 1;
	    //         pdgid = abs(pdgid);
	    //         debug = pdgid == 1000005 || pdgid == 2000005;
	    //         debug = debug || pdgid == 1000001 || pdgid == 2000001;
	    //         debug = debug || pdgid == 1000002 || pdgid == 2000002;
	    //         debug = debug || pdgid == 1000023 || pdgid == 2000024;
	    if (debug) std::cout << mult_mc << "th particle's PDG ID " << pdgid << 
			  " number of parents " << nparents << std::endl; 
	 }// if mcvtx

	 // Get daughters for each particle 
	 if (debug) {
	    mcvtx = mcpart->getDMCvtx();
	    if (mcvtx) {
	       Int_t ndaughters = mcvtx->ndaughters();
	       std::cout << "   particle has " << ndaughters << " daughters" << std::endl;
	       for (Int_t ii=0; ii<ndaughters; ii++) 
		  std::cout << "      PDG ID of " << ii << "th daughter " 
			    << mcvtx->getDaughter(ii)->pdgid() << std::endl;
	    }
	    else std::cout << "   particle is stable " << std::endl;
	 }
	 mult_mc++;
      } // fMCpart
      fH[19]->Fill( mult_mc);
      //std::cout << std::endl << std::endl;
      // Track down kth particle
      Int_t k = 7;
      mcpart=(TMBMCpart*) fMCpart->At(k); 
      if (debug) std::cout << "Track down particle " << k << std::endl;
      if (mcpart) Walk(mcpart,(entry==0));
      if (debug) std::cout << "Done tracking down particle " << k << std::endl;

      //Taus
      TIter iTaus(fTaus);
      TMBTaus* taus=0;
      while ((taus=(TMBTaus*) iTaus())) {
	 if(taus->type()==1) fH[21]->Fill( taus->pT() );
	 else if(taus->type()==2) fH[22]->Fill( taus->pT() );

	 for(Int_t j=0; j<taus->ntrk(); j++) 
	    fH[23]->Fill( taus->E(), taus->GetChargedTrack(j)->E() ); 
      }  // fTaus

      // BCJets
      if (fBCJets) {
	 Int_t nbcjets = fBCJets->GetLast()+1;
	 if (debug) std::cout << " number of bc jets: " << nbcjets << std::endl;

	 TIter iBCJets(fBCJets);
	 TMBBCJet* bcjet=0;
	 Int_t ibcjet=0;
	 while ((bcjet=(TMBBCJet*) iBCJets())) {
	    if (debug) {
	       std::cout << " bcjet pointer: " << bcjet << " for bcjet: " << ibcjet++ << std::endl;
	       // check identities filling this histo
	       std::cout << "    jet pointer: " << bcjet->GetJet() << std::endl; 
	    }
         
	    if (bcjet->GetJet()) fH[37]->Fill(bcjet->GetJet()->pT());
	    // get its muon tagger
	    TMBBCMuTagger mugetter;
	    TMBBCMuTagger* mutag = (TMBBCMuTagger*) mugetter.Get(bcjet);
	    if (debug) std::cout << "    mutag pointer: " << mutag << std::endl;
	    if (!mutag) continue;
	    // get the number of available muon tags for this jet
	    // (identical with the number of associated muons)
	    Int_t ntagmuons = mutag->GetNumTags();
	    for (Int_t ibcmuon=0; ibcmuon<ntagmuons; ibcmuon++){
	       // fill the muon tag's PtRel from the muon
	       // tagger into the next histogram
	       fH[38]->Fill(mutag->GetPtRel(ibcmuon));
	       // fill the muon's charge as reported by the muon
	       // itself into the next hiostogram
	       TMBMuon* bcmuon = mutag->GetMuon(ibcmuon);        
	       if (debug) std::cout << "       muon pointer: " << bcmuon << std::endl;
	       if (!bcmuon) continue;
	       fH[39]->Fill(bcmuon->charge());
	    }
	 }
      } // if fBCJets
   }
   TCanvas* canv=new TCanvas("tmb_tree_canv","TMBTree Example Macro Results",
                             0,0,1000,900);
   canv->Divide(4,10);
   for (UInt_t i=0; i<fH.size(); i++) {
      if (!fH[i]) continue;
      canv->cd(i);
      fH[i]->Draw();
   }

} // TMBTree_bu::Loop