#define mytest_cxx
#include "mytest.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"

void mytest::FileLoop(char* fname,char* outname)
{

  cout << "output will go to " << outname << endl;

  MakeHist(outname);  //  ReadHotCells();

  char filename[80], iFileName[80];
  strcpy(filename, fname);
  cout << "open input file " << filename << endl;
  ifstream RootListFile(filename, ios::in);
  while (RootListFile >> iFileName) 
    {
      if (!strcmp(iFileName,"//"))
        {
          RootListFile >> iFileName;
          continue;
        }
      TTree* tree;
      cout << "open file " << iFileName << endl;
      TFile *newfile = TFile::Open(iFileName,"READ");
      tree = (TTree*) gDirectory->Get("Global");
      Init(tree);
      //      HotCalCells();
      Loop();
      newfile->Close();
    } // while (RootListFile)

  CleanUpHis();

} // FileLoop()

void mytest::Loop()
{
//   In a Root session, you can do:
//      Root > .L mytest.C
//      Root > mytest 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:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch
   if (fChain == 0) return;
   int first = 0;
   Int_t nentries = Int_t(fChain->GetEntries());
   cout << " number of entries " << nentries << endl;
   Int_t nbytes = 0, nb = 0;
   Int_t cvtx = 0;
   for (Int_t jentry=0; jentry<nentries;jentry++) {
      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      Int_t nchvtx = CHVTX_nchvtx;
      if(first == 0){
        Int_t lum = EVENT_lumblk;
        first = 1;
      }
      Int_t nch = CH_nch;
      
      if(nch < 2) continue;

      Int_t lumblk = EVENT_lumblk - lum;
      Float_t flumblk = lumblk;
      H_nvtx->Fill(nchvtx);
      if(nchvtx == 1){
        cvtx++;
        Float_t vx = CHVTX_chxvtx[0];
        Float_t vy = CHVTX_chyvtx[0];
        Float_t vz = CHVTX_chzvtx[0];
        H_vx->Fill(vx);
        H_vy->Fill(vy);
        H_vz->Fill(vz);
        H_vx_vy->Fill(vx,vy);
        H_vx_vz->Fill(vx,vz);
        H_lum_vx->Fill(flumblk,vx);
        
      }
      else continue;
      Float_t vx0 = .4;
      Float_t vy0 = .5;
      Float_t dvx = vx - vx0;
      Float_t dvy = vy - vy0;

      Float_t r = sqrt(dvx*dvx+dvy*dvy);
      H_r->Fill(r);
      
      if (r > 0.2) continue:

      
      H_nch->Fill(nch);
      H_nch_vx->Fill(nch,vx);
      for(Int_t ich=0; ich < nch; ++ich){
        Float_t phip = CH_chphi[ich]/3.14159/2.;
        H_phip->Fill(phip);
        Float_t eta = CH_cheta;
        H_eta->Fill(eta);
        if(eta == 0.0) continue;
        Float_t chpt = CH_chpt[ich];
        H_chpt->Fill(CH_chpt[ich]);
        if (CH_chpt[ich] < 1.)continue;
        Float_t svx = CH_ch_svx[ich]-vx;
        Float_t svy = CH_ch_svy[ich]-vy;
        Float_t svz = CH_ch_svz[ich]-vz;
        H_svx->Fill(svx);
        H_svy->Fill(svy);
        H_svz->Fill(svz);
        H_svx_svy->Fill(svx,svy);
        H_phip_svx->Fill(phip,svx);
      }

      Int_t ngtrk = GTRACK_ngtrk;
      H_ngtrk->Fill(ngtrk);
      for(Int_t ig = 0; ig < ngtrk; ig++){
        Float_t chi_2 = GTRACK_chi_2[ig];
        
        if(chi_2 > 40.) continue;
        Float_t qoverpt = GTRACK_qoverpt[ig];
        if(fabs(qoverpt)>.5)continue;
        Int_t nsmt = GTRACK_nsmt[ig];
        Float_t phi = GTRACK_phidir[ig];
        Float_t phip = phi/3.14159/2.;
        if(phip < 0) phip +=1.;
        Float_t rsig = GTRACK_rsig[ig];
        Float_t drsig = GTRACK_drsig[ig];
        H_nsmt->Fill(nsmt);
        if(nsmt < 4) continue;
        H_rsig->Fill(rsig);
        H_drsig->Fill(drsig);
        if(abs(rsig)<.01) H_gphip->Fill(phip);
        H_phip_rsig->Fill(phip,rsig);
        H_rc->Fill(GTRACK_rc[ig]);
        if(drsig < 0.02){
          H_rsig_cosphi->Fill(2.*rsig*cos(phi));
          H_rsig_sinphi->Fill(2.*rsig*sin(phi));
	  H_rsig_phi->Fill(2.*rsig*sin(phi),2.*rsig*cos(phi));
        }
        if(qoverpt != 0){
          Float_t pt = 1./qoverpt;
          H_pt->Fill(pt);
          H_pt_rsig->Fill(pt,rsig);
          H_qoverpt_rsig->Fill(qoverpt,rsig);
        }
      }


      Int_t nem = EMPART_S_S_n_EM;
      Int_t ngoodem=0;
      Bool_t isW=0;
      Bool_t isZ=0;
      Int_t emidx[20];


      // loop over all em candidates and look for electrons
      for (Int_t iem=0; iem < nem; ++iem){
        Float_t emfrac = EMPART_S_S_EMfrac[iem];
        Float_t eiso   = EMPART_S_S_iso[iem];
        Float_t ehmx9  = EMPART_S_S_HMx9[iem];
        Float_t ept    = EMPART_S_S_pT[iem];
        Float_t ephi   = EMPART_S_S_phi[iem];
        Float_t eeta   = EMPART_S_S_eta[iem];
        Float_t eE     = EMPART_S_S_E[iem];
        Float_t ePx    = EMPART_S_S_pX[iem];
        Float_t ePy    = EMPART_S_S_pY[iem];
        Float_t ePz    = EMPART_S_S_pZ[iem];
        H_em_frac->Fill(emfrac);
        H_em_iso->Fill(eiso);
        H_em_etaphi->Fill(eeta, ephi);
      }

      // if (Cut(ientry) < 0) continue;
   }
   cout << "number of events with vertex = " << cvtx << endl;

}

void mytest::MakeHist(char* outputname){
  _resultFile = TFile::Open(outputname,"recreate");
  TDirectory *emtest = _resultFile->mkdir("em test");
  emtest->cd();

  H_nvtx = new TH1F("nvtx","nvtx",11,-1.,10.);
  H_lum_vx = new TH2F("lum_vx","lum_vx",100,-5.,100.,100,-2.,2.);
  H_nch = new TH1F("nch", "nch", 51,-1.,50.);
  H_vx = new TH1F("vx", "vx", 100,-2., 2. );
  H_vy = new TH1F("vy", "vy", 100,-2., 2. );
  H_vz = new TH1F("vz", "vz", 100,-100., 100 );
  H_vx_vy = new TH2F("vx_vy","vx_vy",100,-2.,2.,100,-2.,2.);
  H_vx_vz = new TH2F("vx_vz","vx_vz",100,-2.,2.,100,-100.,100.);
  H_phip_svx = new TH2F("phip_svx","phip_svx",100,0.,1.,100,-2.,2.);
  H_nch_vx = new TH2F("nch_vx","nch_vx",20,0.,20.,100,-2.,2.);
  H_phip = new TH1F("phip","phip",100,0.,1.);
  H_gphip = new TH1F("gphip","gphip",200,0.,1.);
  H_eta = new TH1F("eta","eta",100,-4.,4.);
  H_svx = new TH1F("svx", "svx", 100,-2., 2. );
  H_svy = new TH1F("svy", "svy", 100,-2., 2. );
  H_svz = new TH1F("svz", "svz", 100,-2., 2. );
  H_svx_svy = new TH2F("svx_svx","svx_svy",100,-2.,2.,100,-2.,2.);

  H_em_frac = new TH1F("emfrac", "emfrac",100,-.1,1.4);
  H_em_iso  = new TH1F("emiso ", "emiso ",100,-2., 2.);
  H_em_etaphi  = new TH2F("eta-phi","eta-phi",40,-4., 4., 70, 0., 7.);
  H_r = new TH1F("r","r",100,0.,2.);
  H_ngtrk = new TH1F("ngtrk","ngtrk",40,0.,40.);
  H_nsmt = new TH1F("nsmt","nsmt",10,0.,10.);
  H_rsig = new TH1F("rsig","rsig",100,-2.,2.);
  H_drsig = new TH1F("drsig","drsig",100,0.,.1);
  H_rsig_cosphi = new TH1F("rsigcosphi","rsigcosphi",100,-2.,2.);
  H_rsig_sinphi = new TH1F("rsigsinphi","rsigsinphi",100,-2.,2.);
  H_phip_rsig = new TH2F("phip_rsig","phip_rsig",100,0.,1.,100,-2.,2.);
  H_pt = new TH1F("pt","pt",100,-50.,50.);
  H_pt_rsig = new TH2F("pt_rsig","pt_rsig",100,-5.,5.,100,-2.,2.);
  H_qoverpt_rsig = new TH2F("qoverpt_rsig","qoverpt_rsig",100,-20.,20.,100,-2.,2.);
  H_rc = new TH1F("rc","rc",100,-100.,100.);
  H_chpt = new TH1F("chpt","chpt",100,-20.,20.);
  H_rsig_phi = new TH2F("rsigphi","rsigphi",100,-4.,4.,100,-4.,4.);
}

void mytest::~mytest(){
  delete H_nvtx;
  delete H_lum_vx;
  delete H_nch;
}

void mytest::CleanUpHis()
{
  // write histograms to output file
  _resultFile->Write();

} // ClanUpHis()