#define omega2pi_selector_cxx
// The class definition in omega2pi_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("omega2pi_selector.C")
// Root > T->Process("omega2pi_selector.C","some options")
// Root > T->Process("omega2pi_selector.C+")
//

#include "omega2pi_selector.h"
#include <TH2.h>
#include <TStyle.h>


void omega2pi_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).

   nSignal=0;
   user_BdtCut = -999.;

   TString option = GetOption();
   FriendOutput = option.Contains("Signal");
   FriendOutputBackground = option.Contains("Background");
   bggenData = option.Contains("bggenData");
   EvalTMVA = option.Contains("EvalTMVA");

   if(EvalTMVA) {
     Int_t firstEqual = option.First("=");
     user_BdtCut = atof(option(firstEqual+1,4).Data());
   }

   if(EvalTMVA) { //only used after BDT training completed

     // setup for TMVA evaluation
     tmvaReader = new TMVA::Reader( "!Color:!Silent" );

     // 3pi PID variables
     tmvaReader->AddVariable("PiPlus1__timingFOM := (PiPlus1__ChiSq_Timing_KinFit==PiPlus1__ChiSq_Timing_KinFit) ? TMath::Prob(PiPlus1__ChiSq_Timing_KinFit,PiPlus1__NDF_Timing) : 0.",&PiPlus1__Timing_FOM);
     tmvaReader->AddVariable("PiPlus1__DCdEdx_FOM := TMath::Prob(PiPlus1__ChiSq_DCdEdx,PiPlus1__NDF_DCdEdx)",&PiPlus1__DCdEdx_FOM);
     tmvaReader->AddVariable("PiPlus1__NDF_Tracking",&PiPlus1__NDF_Tracking_Alt);
     tmvaReader->AddVariable("PiPlus2__Timing_FOM := (PiPlus2__ChiSq_Timing_KinFit==PiPlus2__ChiSq_Timing_KinFit) ? TMath::Prob(PiPlus2__ChiSq_Timing_KinFit,PiPlus2__NDF_Timing) : 0.",&PiPlus2__Timing_FOM);
     tmvaReader->AddVariable("PiPlus2__DCdEdx_FOM := TMath::Prob(PiPlus2__ChiSq_DCdEdx,PiPlus2__NDF_DCdEdx)",&PiPlus2__DCdEdx_FOM);
     tmvaReader->AddVariable("PiPlus2__NDF_Tracking",&PiPlus2__NDF_Tracking_Alt);
     tmvaReader->AddVariable("PiMinus__Timing_FOM := (PiMinus__ChiSq_Timing_KinFit==PiMinus__ChiSq_Timing_KinFit) ? TMath::Prob(PiMinus__ChiSq_Timing_KinFit,PiMinus__NDF_Timing) : 0.",&PiMinus__Timing_FOM);
     tmvaReader->AddVariable("PiMinus__DCdEdx_FOM := TMath::Prob(PiMinus__ChiSq_DCdEdx,PiMinus__NDF_DCdEdx)",&PiMinus__DCdEdx_FOM);
     tmvaReader->AddVariable("PiMinus__NDF_Tracking",&PiMinus__NDF_Tracking_Alt);

     // Computed variables (friend tree)
     tmvaReader->AddVariable("Unused__Max_Proton_FOM",&Unused__Max_Proton_FOM);
     tmvaReader->AddVariable("Unused__Max_KPlus_FOM",&Unused__Max_KPlus_FOM);
     tmvaReader->AddVariable("Unused__Max_KMinus_FOM" ,&Unused__Max_KMinus_FOM);

     // book the MVA methods
     TString methodName = "BDT method";
     TString weightFile = "weights/n3pi_BDT.weights.xml";
     tmvaReader->BookMVA( methodName, weightFile );
   }

}

void omega2pi_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 omega2pi_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 omega2pi_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);
   if(entry==0) cout<<"--- Processing: "<<fChain->GetTree()->GetEntries()<<" total combos"<<endl;
   if(entry%10000 == 0) cout<<"--- ... Processing combo: "<<entry<<endl;

   // Find the actual true signal by cycling through the events to match the appropriate pions to the omega
/*
   Int_t PiPlusMatchToOmega = 0;
   Int_t PiMinusMatchToOmega = 0;

   for(UInt_t i=0; i<NumThrown; i++){
     // select unused tracks with some criteria for max PID FOM
     //if(Unused__NDF_Tracking[i] > 5) {
          //Int_t           Thrown__ParentID[NumThrownMax];   //[NumThrown]
   	  //Int_t           Thrown__PID[NumThrownMax];   //[NumThrown]
   	  //Int_t           Thrown__MatchID[NumThrownMax];
	cout << "Index = " << i <<" Thrown PID = " << Thrown__PID[i] << " Thrown Parent ID = " << Thrown__ParentID[i] <<endl;
	if(Thrown__PID[i] == 211 && Thrown__ParentID[i] == 223){
		PiPlusMatchToOmega = Thrown__MatchID[i];
		cout << "PiPlusMatch = " << PiPlusMatchToOmega << endl;
	}
	if(Thrown__PID[i] == -211 && Thrown__ParentID[i] == 223){
                PiMinusMatchToOmega = Thrown__MatchID[i];
		cout << "PiMinusMatch = " << PiMinusMatchToOmega << endl;
        }
	// select highest CL for proton and kaon hypothesis
       //Float_t PID_FOM = TMath::Prob(Unused__ChiSq_DCdEdx[i]+Unused__ChiSq_Timing[i],Unused__NDF_DCdEdx[i]+Unused__NDF_Timing[i]);
       //if(Unused__PID[i] == 321  && PID_FOM > Unused__Max_KPlus_FOM) Unused__Max_KPlus_FOM=PID_FOM;
       //if(Unused__PID[i] == -321 && PID_FOM > Unused__Max_KMinus_FOM) Unused__Max_KMinus_FOM=PID_FOM;
     //}	
   }
*/ 
   TrueSignal=kFALSE;

if(FriendOutput){
   Int_t PiMinusParent = -1;
   Int_t PiPlusParent = -1;
   Int_t Photon1Parent = -1;
   Int_t Photon2Parent = -1;
   if (PiMinus1__MatchID != -1 && PiMinus2__MatchID != -1 && PiPlus1__MatchID != -1 && PiPlus2__MatchID != -1 && Photon1__MatchID != -1 && Photon2__MatchID != -1 && Proton__MatchID != -1){ 
   	PiMinusParent = Thrown__ParentID[PiMinus1__MatchID]; 
   	PiPlusParent = Thrown__ParentID[PiPlus1__MatchID];
	Photon1Parent = Thrown__ParentID[Photon1__MatchID];
	Photon2Parent = Thrown__ParentID[Photon2__MatchID];
        //cout << "Entering main check" << endl; 
   	if ( PiMinusParent != -1 &&  PiPlusParent != -1 && Photon1Parent != -1 && Photon2Parent != -1 
		&& Thrown__PID[PiMinusParent] == 223 && Thrown__PID[PiPlusParent] == 223 && Thrown__PID[Photon1Parent] == 111 
		&& Thrown__PID[Photon2Parent] == 111 && Thrown__PID[PiMinus1__MatchID] == -211 && Thrown__PID[PiMinus2__MatchID] == -211
		&& Thrown__PID[PiPlus1__MatchID] == 211 && Thrown__PID[PiPlus2__MatchID] == 211 && Thrown__PID[Proton__MatchID] == 2212){
   	//if (PiPlusMatchToOmega == PiPlus1__MatchID && PiMinusMatchToOmega == PiMinus1__MatchID){
		if(Thrown__ParentID[Photon1Parent] != -1 && Thrown__ParentID[Photon2Parent] != -1 && Thrown__PID[Thrown__ParentID[Photon1Parent]] == 223 && Thrown__PID[Thrown__ParentID[Photon2Parent]] == 223){
   			//if (Thrown__ParentID[PiMinus2__MatchID] == -1 && Thrown__ParentID[PiPlus2__MatchID] == -1){
				if (FriendOutput){
				TrueSignal=kTRUE;
				//cout << "The parent PIDs are " << (Thrown__ParentID[PiMinus2__MatchID] == -1 ? -1 : Thrown__PID[ Thrown__ParentID[PiMinus2__MatchID] ]) << " and " << (Thrown__ParentID[PiPlus2__MatchID] == -1 ? -1 : Thrown__PID[ Thrown__ParentID[PiPlus2__MatchID] ]) << " Event Number " << EventNumber << " entry " << entry << endl;
			}
/*
			else{;
				cout << "The parent PIDs are " << (Thrown__ParentID[PiMinus2__MatchID] == -1 ? -1 : Thrown__PID[ Thrown__ParentID[PiMinus2__MatchID] ]) << " and " << (Thrown__ParentID[PiPlus2__MatchID] == -1 ? -1 : Thrown__PID[ Thrown__ParentID[PiPlus2__MatchID] ]) << " Event Number " << EventNumber << " entry " << entry << endl;
			}
*/ 
		} 
   	}

   	else { TrueSignal = kFALSE;}
   }	 
}
   // compute sums from unused calorimeter energies and maximum PID FOMs
   Unused__Max_KPlus_FOM=0.;
   Unused__Max_KMinus_FOM=0.;
   Unused__TotalEnergy_BCAL=0.;
   Unused__TotalEnergy_FCAL=0.;
   // loop over unused hypotheses
   for(UInt_t i=0; i<NumUnused; i++){
     Unused__TotalEnergy_FCAL += Unused__Energy_FCAL[i];
     Unused__TotalEnergy_BCAL += Unused__Energy_BCAL[i];
     // select unused tracks with some criteria for max PID FOM
     if(Unused__NDF_Tracking[i] > 5) {
       // select highest CL for proton and kaon hypothesis
       Float_t PID_FOM = TMath::Prob(Unused__ChiSq_DCdEdx[i]+Unused__ChiSq_Timing[i],Unused__NDF_DCdEdx[i]+Unused__NDF_Timing[i]);
       if(Unused__PID[i] == 321  && PID_FOM > Unused__Max_KPlus_FOM) Unused__Max_KPlus_FOM=PID_FOM;
       if(Unused__PID[i] == -321 && PID_FOM > Unused__Max_KMinus_FOM) Unused__Max_KMinus_FOM=PID_FOM;
     }
   }
   // construct measured missing mass
   TLorentzVector measuredMissingP4(0,0,BeamPhoton__P4_Measured->Energy(),BeamPhoton__P4_Measured->Energy()+0.938);
   Measured__MissingMass = measuredMissingP4.M();
   measuredMissingP4 = measuredMissingP4 - *Proton__P4_Measured - *PiPlus1__P4_Measured - *PiPlus2__P4_Measured - *PiMinus1__P4_Measured - *PiMinus2__P4_Measured - *Photon1__P4_Measured - *Photon2__P4_Measured;
   Measured__MissingMass = measuredMissingP4.M();
   //cout << "Passed calculations" <<endl;
   TLorentzVector pi0__P4(0,0,0,0);
   pi0__P4 = *Photon1__P4_Measured + *Photon2__P4_Measured; 
   MeasuredPi0MassMinusPi0Mass = pi0__P4.M();

   TLorentzVector omega__P4(0,0,0,0);
   omega__P4 = *PiPlus1__P4_Measured + *PiMinus1__P4_Measured + *Photon1__P4_Measured + *Photon2__P4_Measured;
   Omega_MeasuredMass = omega__P4.M();

   if(EvalTMVA) ProcessTMVA();
   if(FriendOutput || FriendOutputBackground) outTree->Fill();
   //cout << "Written" << endl;
   return kTRUE;
}

void omega2pi_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 omega2pi_selector::Terminate()
{
	cout << "Entering Terminate()" << endl;
   if(FriendOutput || FriendOutputBackground){
     //cout<<"found "<<nSignal<<" true 3pi n events"<<endl;
     outFile->Write();
     outFile->Close();
   }
   // 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.

}