#define KsHunt_selector_cxx
// The class definition in KsHunt_selector.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("KsHunt_selector.C")
// Root > T->Process("KsHunt_selector.C","some options")
// Root > T->Process("KsHunt_selector.C+")
//

#include "KsHunt_selector.h"

using namespace std;

void KsHunt_selector::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

   // Outfile ROOT file name is passed in as option,
   // and is created by main analyze_KsHunt.
   // Open the file and cd into it.
   outfile = new TFile(option.Data(),"UPDATE");
   outfile->cd();

   // kfit info
   hkfitCL = new TH1F("hkfitCL",";kfit CL;",200,0,1);

   // M(pi+,pi-)
   hM2pi = new TH1F("hM2pi",";M(#pi^{+},#pi^{-});",400,0,2.0);
   // vrho vs M(pi+,pi-)
   hvrho_M2pi = new TH2F("hvrho_M2pi",";M(#pi^{+},#pi^{-});v_{#rho}",400,0,2.0,200,0,30);
   // vz vs M(pi+,pi-)
   hvz_M2pi = new TH2F("hvz_M2pi",";M(#pi^{+},#pi^{-});v_{z}",400,0,2.0,400,0,400);
   // M(pi+,pi-) with selection on vrho, vz
   hM2pi_allVertexSelect = new TH1F("hM2pi_allVertexSelect",";M(#pi^{+},#pi^{-});",400,0,2.0);

   // M(pi+,pi-,unused pi+-)
   hM3pi = new TH1F("hM3pi",";M(#pi^{+},#pi^{-},#pi^{#pm}_{unused});",400,0,2.0);
   // M(pi+,pi-,unused pi) with selection on Ks, vertices
   hM3pi_allVertexKsSelect = new TH1F("hM3pi_allVertexKsSelect",";M(#pi^{+},#pi^{-},#pi^{#pm}_{unused})",400,0,2.0);

   // M(unused p,unused pi-)
   hMppi = new TH1F("hMppi",";M(p_{unused},#pi^{-}_{unused});",200,1.0,2.0);
   hMppi_allVertexKsSelect = new TH1F("hMppi_allVertexKsSelect",";M(p_{unused},#pi^{-}_{unused})",400,1.0,2.0);

   // MM from Ks vs M(pi+,pi-)
   hMM_M2pi = new TH2F("hMM_M2pi",";M(#pi^{+},#pi^{-});MM",400,0,2.0,400,0,2.0);
   hMM_M2pi_allVertexKsSelect = new TH2F("hMM_M2pi_allVertexKsSelect",";M(#pi^{+},#pi^{-});MM",400,0,2.0,400,0,2.0);

   // Momentum
   hPpip_Ppim = new TH2F("hPpip_Ppim",";p_{#pi+};p_{#pi-}",400,0,10,400,0,10);
   hPpip_Ppim_allVertexKsSelect = new TH2F("hPpip_Ppim_allVertexKsSelect",";p_{#pi+};p_{#pi-}",400,0,10,400,0,10);
   hPKs = new TH1F("hPKs",";p_{Ks};",400,0,10);
   hPKs_allVertexKsSelect = new TH1F("hPKs_allVertexKsSelect",";p_{Ks};",400,0,10);

   // dE/dx vs |p|
   // all pi+, pi-
   hdEdx_pip_CDC = new TH2F("hdEdx_pip_CDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_pip_FDC = new TH2F("hdEdx_pip_FDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_pip_ST = new TH2F("hdEdx_pip_ST",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   hdEdx_pim_CDC = new TH2F("hdEdx_pim_CDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_pim_FDC = new TH2F("hdEdx_pim_FDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_pim_ST = new TH2F("hdEdx_pim_ST",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   // pi+, pi- in Ks mass range
   hdEdx_pip_CDC_selectKs = new TH2F("hdEdx_pip_CDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_pip_FDC_selectKs = new TH2F("hdEdx_pip_FDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_pip_ST_selectKs = new TH2F("hdEdx_pip_ST_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   hdEdx_pim_CDC_selectKs = new TH2F("hdEdx_pim_CDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_pim_FDC_selectKs = new TH2F("hdEdx_pim_FDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_pim_ST_selectKs = new TH2F("hdEdx_pim_ST_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   // unused K+-
   hdEdx_kp_CDC = new TH2F("hdEdx_kp_CDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_kp_FDC = new TH2F("hdEdx_kp_FDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_kp_ST = new TH2F("hdEdx_kp_ST",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   hdEdx_km_CDC = new TH2F("hdEdx_km_CDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_km_FDC = new TH2F("hdEdx_km_FDC",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_km_ST = new TH2F("hdEdx_km_ST",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   hdEdx_kp_CDC_selectKs = new TH2F("hdEdx_kp_CDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_kp_FDC_selectKs = new TH2F("hdEdx_kp_FDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_kp_ST_selectKs = new TH2F("hdEdx_kp_ST_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   hdEdx_km_CDC_selectKs = new TH2F("hdEdx_km_CDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.01);
   hdEdx_km_FDC_selectKs = new TH2F("hdEdx_km_FDC_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,0,0.001);
   hdEdx_km_ST_selectKs = new TH2F("hdEdx_km_ST_selectKs",";|p| (GeV/c);dE/dx",200,0,5,200,-10,10);

   // vertices
   hvy_vx_pip_measured = new TH2F("hvy_vx_pip_measured",";v_{x} for #pi^{+};v_{y} for #pi^{+}",200,-20,20,200,-20,20);
   hvy_vx_pip = new TH2F("hvy_vx_pip",";v_{x} for #pi^{+};v_{y} for #pi^{+}",200,-20,20,200,-20,20);

   hvz_pip_measured = new TH1F("hvz_pip_measured",";v_{z} for #pi^{-};",400,0,200);
   hvz_pip = new TH1F("hvz_pip",";v_{z} for #pi^{-};",400,0,200);

   hvy_vx_pim_measured = new TH2F("hvy_vx_pim_measured",";v_{x} for #pi^{-};v_{y} for #pi^{-}",200,-20,20,200,-20,20);
   hvy_vx_pim = new TH2F("hvy_vx_pim",";v_{x} for #pi^{-};v_{y} for #pi^{-}",200,-20,20,200,-20,20);

   hvz_pim_measured = new TH1F("hvz_pim_measured",";v_{z} for #pi^{-};",400,0,200);
   hvz_pim = new TH1F("hvz_pim",";v_{z} for #pi^{-};",400,0,200);

   hvz_pip_pim_measured = new TH2F("hvz_pip_pim_measured",";v_{z} for #pi^{-};v_{z} for #pi^{+}",400,0,200,400,0,200);
   hvz_pip_pim = new TH2F("hvz_pip_pim",";v_{z} for #pi^{-};v_{z} for #pi^{+}",400,0,200,400,0,200);
}

void KsHunt_selector::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}

Bool_t KsHunt_selector::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either KsHunt_selector::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

  GetEntry(entry);
  // cout << "working on entry" << setw(8) << entry << " (EventNumber is " << EventNumber << ")" << endl;
  if(entry%10000==0) cout << " --- Processing combo: " << setw(8) << entry << "/" << fChain->GetTree()->GetEntries() << endl;

  // kfit info
  hkfitCL->Fill(TMath::Prob(ChiSq_KinFit,NDF_KinFit));

  // Cut on kfit CL
  if(NDF_KinFit<=0 || TMath::Prob(ChiSq_KinFit,NDF_KinFit) < 0.10) return kTRUE;

  //////////////////////////
  //      Get 4-moms.     //
  //////////////////////////

  // no kfit
  p4init_photon_PreKinFit = *BeamPhoton__P4_Measured;
  p4init_proton_PreKinFit.SetXYZT(0,0,0,ParticleMass(Proton));
  p4pip_PreKinFit = *PiPlus__P4_Measured;
  p4pim_PreKinFit = *PiMinus__P4_Measured;

  x4pip_PreKinFit = *PiPlus__X4_Measured;
  x4pim_PreKinFit = *PiMinus__X4_Measured;

  p4init_photon = *BeamPhoton__P4_Measured;
  p4init_proton.SetXYZT(0,0,0,ParticleMass(Proton));
  p4pip = *PiPlus__P4_KinFit;
  p4pim = *PiMinus__P4_KinFit;

  x4pip = *PiPlus__X4_KinFit;
  x4pim = *PiMinus__X4_KinFit;

  // Set up selection criteria
  Bool_t isGoodVrho   = false;
  Bool_t isGoodVz     = false;
  Bool_t isGoodKsMass = false;

  const Double_t VRHOCUT = 0.0;
  const Double_t VZCUT   = 70.0;
  const Double_t KSCUT   = 0.013 * 3.;

  if(x4pip.Vect().Perp()                              > VRHOCUT) isGoodVrho   = true;
  if(x4pip.Z()                                        > VZCUT  ) isGoodVz     = true;
  if(fabs((p4pip + p4pim).M() - ParticleMass(KShort)) < KSCUT  ) isGoodKsMass = true;

  /*
  cout << setw(8) << x4pip.X() << setw(8) << x4pim.X() << endl;
  cout << setw(8) << x4pip.Y() << setw(8) << x4pim.Y() << endl;
  cout << setw(8) << x4pip.Z() << setw(8) << x4pim.Z() << endl;
  cout << setw(8) << x4pip.T() << setw(8) << x4pim.T() << endl;
  */

  // M(pi+,pi-)
  hM2pi->Fill((p4pip + p4pim).M());
  hvrho_M2pi->Fill((p4pip + p4pim).M(), x4pip.Vect().Perp());
  hvz_M2pi->Fill((p4pip + p4pim).M(), x4pip.Z());

  // momentum
  hPpip_Ppim->Fill(p4pip.P(),p4pim.P());
  hPKs->Fill((p4pip+p4pim).P());

  // MM vs M(pi+,pi-)
  hMM_M2pi->Fill((p4pip + p4pim).M(),(p4init_photon + p4init_proton - p4pip - p4pim).M());

  // Selection on vrho, vz
  if(isGoodVrho && isGoodVz){
    hM2pi_allVertexSelect->Fill((p4pip + p4pim).M());
  }

  // Selection on vrho, vz, Ks mass
  if(isGoodVrho && isGoodVz && isGoodKsMass){
    // p_pi+ vs p_pi-
    hPpip_Ppim_allVertexKsSelect->Fill(p4pip.P(),p4pim.P());
    // p_Ks
    hPKs_allVertexKsSelect->Fill((p4pip+p4pim).P());
    // MM vs M(pi+,pi-)
    hMM_M2pi_allVertexKsSelect->Fill((p4pip + p4pim).M(),(p4init_photon + p4init_proton - p4pip - p4pim).M());
  }

  // dE/dx
  hdEdx_pip_CDC->Fill(p4pip.P(),PiPlus__dEdx_CDC);
  hdEdx_pip_FDC->Fill(p4pip.P(),PiPlus__dEdx_FDC);
  if(PiPlus__dEdx_ST != 0)
    hdEdx_pip_ST->Fill(p4pip.P(),PiPlus__dEdx_ST);

  hdEdx_pim_CDC->Fill(p4pim.P(),PiMinus__dEdx_CDC);
  hdEdx_pim_FDC->Fill(p4pim.P(),PiMinus__dEdx_FDC);
  if(PiMinus__dEdx_ST != 0)
    hdEdx_pim_ST->Fill(p4pim.P(),PiMinus__dEdx_ST);

  // pi+, pi- in Ks mass range
  if(fabs((p4pip + p4pim).M() - ParticleMass(KShort)) < 0.013 * 3.){
    hdEdx_pip_CDC_selectKs->Fill(p4pip.P(),PiPlus__dEdx_CDC);
    hdEdx_pip_FDC_selectKs->Fill(p4pip.P(),PiPlus__dEdx_FDC);
    if(PiPlus__dEdx_ST != 0)
      hdEdx_pip_ST_selectKs->Fill(p4pip.P(),PiPlus__dEdx_ST);

    hdEdx_pim_CDC_selectKs->Fill(p4pim.P(),PiMinus__dEdx_CDC);
    hdEdx_pim_FDC_selectKs->Fill(p4pim.P(),PiMinus__dEdx_FDC);
    if(PiMinus__dEdx_ST != 0)
      hdEdx_pim_ST_selectKs->Fill(p4pim.P(),PiMinus__dEdx_ST);
  }

  // vertices
  hvy_vx_pip_measured->Fill(x4pip_PreKinFit.X(),x4pip_PreKinFit.Y());
  hvy_vx_pip->Fill(x4pip.X(),x4pip.Y());

  hvz_pip_measured->Fill(x4pip_PreKinFit.Z());
  hvz_pip->Fill(x4pip.Z());

  hvy_vx_pim_measured->Fill(x4pim_PreKinFit.X(),x4pim_PreKinFit.Y());
  hvy_vx_pim->Fill(x4pim.X(),x4pim.Y());

  hvz_pim_measured->Fill(x4pim_PreKinFit.Z());
  hvz_pim->Fill(x4pim.Z());

  hvz_pip_pim_measured->Fill(x4pim_PreKinFit.Z(),x4pip_PreKinFit.Z());
  hvz_pip_pim->Fill(x4pim.Z(),x4pip.Z());


  // Look at unused tracks
  assert(NumUnused < NumUnusedMax);

  Int_t nUnusedPiMinus = 0;
  Int_t nUnusedProton  = 0;
  for(Int_t i=0;i<NumUnused;i++){
    if(Unused__NDF_Tracking[i]>0){
      p4unused[i] = (TLorentzVector*)Unused__P4_Measured->ConstructedAt(i);
      p4unused_ = *p4unused[i];
      x4unused[i] = (TLorentzVector*)Unused__X4_Measured->ConstructedAt(i);
      x4unused_ = *x4unused[i];

      // charged pions
      if(fabs(p4unused[i]->M() - ParticleMass(PiPlus)) < 1.e-2){
	hM3pi->Fill((p4pip + p4pim + p4unused_).M());

	// Selection on vrho, vz, Ks mass
	if(isGoodVrho && isGoodVz && isGoodKsMass){
	  hM3pi_allVertexKsSelect->Fill((p4pip + p4pim + p4unused_).M());
	}

	// check charged track PID
	if(Unused__PID[i] == PDGtype(PiMinus)){
	  p4unused_PiMinus[nUnusedPiMinus] = p4unused_;
	  x4unused_PiMinus[nUnusedPiMinus] = x4unused_;
	  objectID_PiMinus[nUnusedPiMinus] = Unused__ObjectID[i];
	  nUnusedPiMinus++;
	}
      }

      // protons
      if(fabs(p4unused[i]->M() - ParticleMass(Proton)) < 1.e-2){
	// check charged track PID
	if(Unused__PID[i] == PDGtype(Proton)){
	  p4unused_Proton[nUnusedProton] = p4unused_;
	  x4unused_Proton[nUnusedProton] = x4unused_;
	  objectID_Proton[nUnusedProton] = Unused__ObjectID[i];
	  nUnusedProton++;
	}
      }
      
      // unused K+
      if(Unused__PID[i] == PDGtype(KPlus)){
	hdEdx_kp_CDC->Fill(p4unused_.P(),Unused__dEdx_CDC[i]);
	hdEdx_kp_FDC->Fill(p4unused_.P(),Unused__dEdx_FDC[i]);
	if(Unused__dEdx_ST[i] != 0)
	  hdEdx_kp_ST->Fill(p4unused_.P(),Unused__dEdx_ST[i]);

	if(fabs((p4pip + p4pim).M() - ParticleMass(KShort)) < 0.013 * 3.){
	  hdEdx_kp_CDC_selectKs->Fill(p4unused_.P(),Unused__dEdx_CDC[i]);
	  hdEdx_kp_FDC_selectKs->Fill(p4unused_.P(),Unused__dEdx_FDC[i]);
	  if(Unused__dEdx_ST[i] != 0)
	    hdEdx_kp_ST_selectKs->Fill(p4unused_.P(),Unused__dEdx_ST[i]);
	}
      }

      // unused K-
      if(Unused__PID[i] == PDGtype(KMinus)){
	hdEdx_km_CDC->Fill(p4unused_.P(),Unused__dEdx_CDC[i]);
	hdEdx_km_FDC->Fill(p4unused_.P(),Unused__dEdx_FDC[i]);
	if(Unused__dEdx_ST[i] != 0)
	  hdEdx_km_ST->Fill(p4unused_.P(),Unused__dEdx_ST[i]);

	if(fabs((p4pip + p4pim).M() - ParticleMass(KShort)) < 0.013 * 3.){
	  hdEdx_km_CDC_selectKs->Fill(p4unused_.P(),Unused__dEdx_CDC[i]);
	  hdEdx_km_FDC_selectKs->Fill(p4unused_.P(),Unused__dEdx_FDC[i]);
	  if(Unused__dEdx_ST[i] != 0)
	    hdEdx_km_ST_selectKs->Fill(p4unused_.P(),Unused__dEdx_ST[i]);
	}
      }
      
    } // end of NDF for tracking of unused
  } // end of loop over unused particles

  // cout << "--------------------------------------------" << endl;
  // cout << "total number of unused: " << NumUnused << endl;
  // for(Int_t i=0;i<nUnusedPiMinus;i++){
  //   cout << "PiMinus Id = " << objectID_PiMinus[i] << "("
  //        << p4unused_PiMinus[i].X() << ", "
  //        << p4unused_PiMinus[i].Y() << ", "
  //        << p4unused_PiMinus[i].Z() << ", "
  //        << p4unused_PiMinus[i].T() << ")" << endl;
  // }

  // for(Int_t i=0;i<nUnusedProton;i++){
  //   cout << "Proton Id = " << objectID_Proton[i] << "("
  //        << p4unused_Proton[i].X() << ", "
  //        << p4unused_Proton[i].Y() << ", "
  //        << p4unused_Proton[i].Z() << ", "
  //        << p4unused_Proton[i].T() << ")" << endl;
  // }

  for(Int_t iPiMinus=0;iPiMinus<nUnusedPiMinus;iPiMinus++){
    for(Int_t iProton=0;iProton<nUnusedProton;iProton++){
      if(objectID_PiMinus[iPiMinus] != objectID_Proton[iProton]){
	hMppi->Fill((p4unused_PiMinus[iPiMinus] + p4unused_Proton[iProton]).M());

	// Selection on vrho, vz, Ks mass
	if(isGoodVrho && isGoodVz && isGoodKsMass){
	  hMppi_allVertexKsSelect->Fill((p4unused_PiMinus[iPiMinus] + p4unused_Proton[iProton]).M());
	}
      }
    }
  }

  return kTRUE;
}

void KsHunt_selector::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void KsHunt_selector::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
  outfile->Write();

}