//-----------------------------------------------------------------------------------//
//                                                                                   //
//   Filter module for the B -> pi/pi0/eta/eta'/rho/rho0/omega/gamma l nu analyses   //
//                                                                                   //
//     Sylvie Brunet  2003 Universite de Montreal                                    //
//     David Cote     2003 Universite de Montreal                                    //
//                                                                                   //
//   Documentation: BAD 740                                                          //
//                                                                                   //
//-----------------------------------------------------------------------------------//


#include "FilterTools/XSLBtoXulnuFilter.hh"

///////////////
#include "CLHEP/Alist/AList.h"
#include "CLHEP/Alist/AIterator.h"
#include "AbsEnv/AbsEnv.hh"
#include "AbsEvent/AbsEvent.hh"
#include "AbsEvent/getTmpAList.hh"
#include "ErrLogger/ErrLog.hh"
#include "AbsEventTag/AbsEventTag.hh"
#include "BetaCoreTools/BtaOpMakeTree.hh"
#include "Beta/EventInfo.hh"
#include "BetaCoreTools/BtaBooster.hh"
#include "PDT/Pdt.hh"

#include <iomanip.h>

//For testing
#include "HepTuple/Tuple.h"
#include "HepTuple/TupleManager.h"
#include "GenEnv/GenEnv.hh"


XSLBtoXulnuFilter::XSLBtoXulnuFilter( const char* const theName, 
				      const char* const theDescription )
  : TagFilterModule( theName, theDescription )
  , mB0(5.2794), mB(5.2790)
  
  //UsrEventBlock
  , _eventData("eventData",this,"XSLBtoXulnuEventData")
  , _decayMode("decayMode")

  //I/O lists
  , _eventInfoList("eventInfoList", this, "Default" )
  , _InputLeptonList("inputLeptonList",this,"eLHBremGTL")
  , _InputPiList("inputPiList",this,"piLHLooseGTL")
  , _InputPi0List("inputPi0List",this,"pi0AllVeryLoose") 
  , _InputEtaList("inputEtaList",this,"etaAllSpecial3Pi")
  , _InputEtapList("inputEtapList",this,"etaPAllPID")
  , _InputRho0List("inputRho0List",this,"rho0VeryLooseMassAndPid")
  , _InputRhoCList("inputRhoCList",this,"rhoCVeryLooseMassAndPid")  
  , _InputOmegaList("inputOmegaList",this,"omegaVeryLooseAndPid")
  , _InputGammaList("inputGammaList",this,"GoodPhotonDefault")
  , _OutputYList("outputYList",this,"XSLBtoXulnuSkimmedYlist") 

  //event parameters
  , _r2allMAX("r2allMAX",this,0.5)//no cut:>=1
  , _ishMIN("ishMIN",this,1)//no cut:0, cut:1
  , _trigMIN("trigMIN",this,0)//no cut:0, cut:1
  //The one lepton tight requirement in the evt is fulfilled automatically by the taking tight electron and muon lists

  //list building parameters
  , _pLABmin_electron("pLABmin_electron",this,0.5)//no cut:0
  , _pLABmin_muon("pLABmin_muon",this,1)//no cut:0
  , _thetaLABmin_electron("thetaLABmin_electron",this,0.41)//no cut:0
  , _thetaLABmin_muon("thetaLABmin_muon",this,0.41)//no cut:0
  , _thetaLABmax_electron("thetaLABmax_electron",this, 2.54)//no cut: 3.14
  , _thetaLABmax_muon("thetaLABmax_muon",this, 2.54)//no cut: 3.14
  , _pLABmax_pi0("pLABmax_pi0",this,10)//no cut:100
  , _pLABmax_eta("pLABmax_eta",this,10)//no cut:100
  , _pLABmax_etap("pLABmax_etap",this,10)//no cut:100
  , _pLABmax_rho0("pLABmax_rho0",this,10)//no cut:100
  , _pLABmax_rhoC("pLABmax_rhoC",this,10)//no cut:100
  , _pLABmax_omega("pLABmax_omega",this,10)//no cut:100
  , _pLABmax_gamma("pLABmax_gamma",this,10)//no cut:100

  //Y selection
  //The cosBY cut is the same for all hadrons so we don't duplicate parameters
  //The cosBY cut is also the same for all leptons, since we're using BremRecovered electrons, so we don't duplicate parameters
  , _cosBYmin("cosBYmin",this,-1.1)//no cut:-5000
  , _cosBYmax("cosBYmax",this,1.3)//no cut:+1000
  , _pStar2Dslope_pi("pStar2Dslope_pi",this,1)
  , _pStar2Dslope_pi0("pStar2Dslope_pi0",this,1)
  , _pStar2Dslope_eta("pStar2Dslope_eta",this,1)
  , _pStar2Dslope_etap("pStar2Dslope_etap",this,0.69)
  , _pStar2Dslope_rhoC("pStar2Dslope_rhoC",this,1)
  , _pStar2Dslope_rho0("pStar2Dslope_rho0",this,1)
  , _pStar2Dslope_omega("pStar2Dslope_omega",this,1)
  , _pStar2Dsum_eta("pStar2Dsum_eta",this,2.8) //no cut: 0
  , _pStar2Dsum_pi0("pStar2Dsum_pi0",this,2.8)//no cut: 0
  , _pStar2Dsum_pi("pStar2Dsum_pi",this,2.8)//no cut: 0
  , _pStar2Dsum_etap("pStar2Dsum_etap",this,2.4)//no cut: 0
  , _pStar2Dsum_rhoC("pStar2Dsum_rhoC",this,2.65)//no cut: 0
  , _pStar2Dsum_rho0("pStar2Dsum_rho0",this,2.65)//no cut: 0
  , _pStar2Dsum_omega("pStar2Dsum_omega",this,2.65)//no cut: 0
  , _pXuUpsMin_pi("pXuUpsMin_pi",this,1.3)//no cut:0
  , _pXuUpsMin_pi0("pXuUpsMin_pi0",this,1.3)//no cut:0
  , _pXuUpsMin_eta("pXuUpsMin_eta",this,1.3)//no cut:0
  , _pXuUpsMin_etap("pXuUpsMin_etap",this,1.65)//no cut:0
  , _pXuUpsMin_rhoC("pXuUpsMin_rhoC",this,1.3)//no cut:0
  , _pXuUpsMin_rho0("pXuUpsMin_rho0",this,1.3)//no cut:0
  , _pXuUpsMin_omega("pXuUpsMin_omega",this,1.3)//no cut:0
  , _pLepUpsMin_pi("pLepUpsMin_pi",this,2.2)//no cut:0
  , _pLepUpsMin_pi0("pLepUpsMin_pi0",this,2.2)//no cut:0
  , _pLepUpsMin_eta("pLepUpsMin_eta",this,2.1)//no cut:0
  , _pLepUpsMin_etap("pLepUpsMin_etap",this,2)//no cut:0
  , _pLepUpsMin_rhoC("pLepUpsMin_rhoC",this,2)//no cut:0
  , _pLepUpsMin_rho0("pLepUpsMin_rho0",this,2)//no cut:0
  , _pLepUpsMin_omega("pLepUpsMin_omega",this,2)//no cut:0

  //modes to be studied
  , _analyzeChargedPi("AnalyzeChargedPi",this,true)
  , _analyzeNeutralPi("AnalyzeNeutralPi",this,true)
  , _analyzeEta("AnalyzeEta",this,true)
  , _analyzeEtap("AnalyzeEtap",this,true)
  , _analyzeChargedRho("AnalyzeChargedRho",this,true)
  , _analyzeNeutralRho("AnalyzeNeutralRho",this,true)
  , _analyzeOmega("AnalyzeOmega",this,true)
  , _analyzeGamma("AnalyzeGamma",this,true)

  //bool
  , _MyVerbose("MyVerbose",this,false)
  , _doPionPidYourself("doPionPidYourself",this,false)
{
  //UsrEventBlock
  commands()->append( &_eventData );
  _eventBlock.addVariable( _decayMode );

  //I/O lists
  commands()->append( & _eventInfoList);
  commands()->append( & _InputLeptonList);
  commands()->append( & _InputPiList);
  commands()->append( & _InputPi0List);
  commands()->append( & _InputEtaList);
  commands()->append( & _InputEtapList);
  commands()->append( & _InputRho0List);
  commands()->append( & _InputRhoCList);
  commands()->append( & _InputOmegaList);
  commands()->append( & _InputGammaList);
  commands()->append( & _OutputYList);

  //evt cuts
  commands()->append( & _r2allMAX );
  commands()->append( & _ishMIN );
  commands()->append( & _trigMIN );

  //trk building cuts
  commands()->append( & _pLABmin_electron );
  commands()->append( & _pLABmin_muon );
  commands()->append( & _thetaLABmin_electron );
  commands()->append( & _thetaLABmin_muon );
  commands()->append( & _thetaLABmax_electron );
  commands()->append( & _thetaLABmax_muon );
  commands()->append( & _pLABmax_pi0 );
  commands()->append( & _pLABmax_eta );
  commands()->append( & _pLABmax_etap );
  commands()->append( & _pLABmax_rhoC );
  commands()->append( & _pLABmax_rho0 );
  commands()->append( & _pLABmax_omega );
  commands()->append( & _pLABmax_gamma );

  //Y selection
  commands()->append( & _cosBYmin );
  commands()->append( & _cosBYmax );

  commands()->append( & _pStar2Dslope_pi );
  commands()->append( & _pStar2Dslope_pi0 );
  commands()->append( & _pStar2Dslope_eta );
  commands()->append( & _pStar2Dslope_etap );
  commands()->append( & _pStar2Dslope_rhoC );
  commands()->append( & _pStar2Dslope_rho0 );
  commands()->append( & _pStar2Dslope_omega );
  commands()->append( & _pStar2Dsum_eta );
  commands()->append( & _pStar2Dsum_pi0 );
  commands()->append( & _pStar2Dsum_pi );
  commands()->append( & _pStar2Dsum_etap );
  commands()->append( & _pStar2Dsum_rhoC );
  commands()->append( & _pStar2Dsum_rho0 );
  commands()->append( & _pStar2Dsum_omega );
  commands()->append( & _pXuUpsMin_pi );
  commands()->append( & _pXuUpsMin_pi0 );
  commands()->append( & _pXuUpsMin_eta );
  commands()->append( & _pXuUpsMin_etap );
  commands()->append( & _pXuUpsMin_rho0 );
  commands()->append( & _pXuUpsMin_rhoC );
  commands()->append( & _pXuUpsMin_omega );
  commands()->append( & _pLepUpsMin_pi );
  commands()->append( & _pLepUpsMin_pi0 );
  commands()->append( & _pLepUpsMin_eta );
  commands()->append( & _pLepUpsMin_etap );
  commands()->append( & _pLepUpsMin_rho0 );
  commands()->append( & _pLepUpsMin_rhoC );
  commands()->append( & _pLepUpsMin_omega );

  //bool
  commands()->append( & _analyzeChargedPi );
  commands()->append( & _analyzeNeutralPi );
  commands()->append( & _analyzeEta );
  commands()->append( & _analyzeEtap );
  commands()->append( & _analyzeChargedRho );
  commands()->append( & _analyzeNeutralRho );
  commands()->append( & _analyzeOmega );
  commands()->append( & _analyzeGamma );
  commands()->append( & _MyVerbose );
  commands()->append( & _doPionPidYourself );
}

XSLBtoXulnuFilter::~XSLBtoXulnuFilter()
{
}


AppResult
XSLBtoXulnuFilter::beginRun( AbsEvent* anEvent ) {    
  return AppResult::OK;
}


AppResult
XSLBtoXulnuFilter::beginJob( AbsEvent* anEvent )
{
  ErrMsg(routine)<<"begin filter Job"<<endmsg; 

  _TotalNbEvents=0;
  _NbAfterTrigger=0;
  _NbAfterISHCuts=0;
  _NbAfterR2AllCuts=0;
  _NbAfterOneLeptonTight=0;
  _NbAfterEventCuts=0;
  _NbAllModesAfterCoupleCuts=0;
  _NbPiAfterCoupleCuts=0;
  _NbPiAfterLeptonSelCuts=0;
  _NbPiAfterCosBY=0;
  _NbPiAfterpStar2D=0;
  _NbPi0AfterCoupleCuts=0;
  _NbPi0AfterLeptonSelCuts=0;
  _NbPi0AfterCosBY=0;
  _NbPi0AfterpStar2D=0;
  _NbEtaAfterCoupleCuts=0;
  _NbEtaAfterLeptonSelCuts=0;
  _NbEtaAfterCosBY=0;
  _NbEtaAfterpStar2D=0;
  _NbEtapAfterCoupleCuts=0;
  _NbEtapAfterLeptonSelCuts=0;
  _NbEtapAfterCosBY=0;
  _NbEtapAfterpStar2D=0;
  _NbRhoCAfterCoupleCuts=0;
  _NbRhoCAfterLeptonSelCuts=0;
  _NbRhoCAfterCosBY=0;
  _NbRhoCAfterpStar2D=0;
  _NbRho0AfterCoupleCuts=0;
  _NbRho0AfterLeptonSelCuts=0;
  _NbRho0AfterCosBY=0;
  _NbRho0AfterpStar2D=0;
  _NbOmeAfterCoupleCuts=0;
  _NbOmeAfterLeptonSelCuts=0;
  _NbOmeAfterCosBY=0;
  _NbOmeAfterpStar2D=0;
  _NbGammaAfterCoupleCuts=0;
  _NbGammaAfterLeptonSelCuts=0;
  _NbGammaAfterCosBY=0;
  _NbCoupleAllModesAfterAllCuts=0;
  _NbCouplePiAfterAllCuts=0;
  _NbCouplePi0AfterAllCuts=0;
  _NbCoupleEtaAfterAllCuts=0;
  _NbCoupleRhoCAfterAllCuts=0;
  _NbCoupleRho0AfterAllCuts=0;
  _NbCoupleOmegaAfterAllCuts=0;
  _NbCoupleGammaAfterAllCuts=0;

  //For testing:
  HepTupleManager* manager = gblEnv->getGen()->ntupleManager();
  assert(manager != 0);
  
  //ntuples
  _TestingNtuple = manager->ntuple("BananaJoe",666);
  assert( NULL != _TestingNtuple);


  return AppResult::OK;
}


AppResult
XSLBtoXulnuFilter::event( AbsEvent* anEvent )
{
  // Build the tag accessor using the base class
  TagFilterModule::event( anEvent );

  //First thing: create the output list
  //The getTmpAList function makes the list available for all AppModule's children within an event
  //A zero length YList will tell to the reader to reject the event
  HepAList< BtaCandidate >* outYList;
  getTmpAList( anEvent, outYList, _OutputYList.value() );
  assert( outYList!=0 );

  _TotalNbEvents+=1;

  //Global event Cuts
  //The _CMS booster object is not created at this point so we cannot use QuitEvent()!
  if (SurvivedTrigger( anEvent ))  _NbAfterTrigger+=1; 
  else { setPassed(false); return AppResult::OK; } 
  if (SurvivedISH( anEvent ))  _NbAfterISHCuts+=1; 
  else { setPassed(false); return AppResult::OK; } 
  if(SurvivedR2All( anEvent ))  _NbAfterR2AllCuts+=1; 
  else { setPassed(false); return AppResult::OK; }

  // 4-Vector P(e+)+P(e-) CM Upsilon(4S)
  HepAList< EventInfo >* infoList = Ifd<HepAList< EventInfo > >::get(anEvent, _eventInfoList.value());
  if (infoList == 0){ ErrMsg(fatal) << "Could not locate event info list of name ";}
  EventInfo* eventInfo = infoList->first();  
  HepLorentzVector UpsP4 = eventInfo->cmFrame();
  double s = UpsP4.mag2(); 

  //Create the Ups(4S) frame booster
  _CMS = new BtaBooster( UpsP4 );

  //Making lists
  MakeHadronsAndLeptonsList( anEvent );

  if(NewLepList->length()!=0) 
    { 
      _NbAfterEventCuts+=1;  
      _NbAfterOneLeptonTight+=1;    
    }
  else { setPassed(false); return QuitEvent(); }

  //////////////////////////////////////
  //  We now build the YList made out of all Y couples passing their mode-specific criteria
  /////////////////////////////////////
  BtaCandidate* lepton;
  BtaCandidate* hadron;

  HepAListIterator<BtaCandidate> iter_lep(*NewLepList);  
  HepAListIterator<BtaCandidate> iter_pi(NewPiList);  
  HepAListIterator<BtaCandidate> iter_pi0(NewPi0List);  
  HepAListIterator<BtaCandidate> iter_eta(NewEtaList);  
  HepAListIterator<BtaCandidate> iter_etap(NewEtapList);  
  HepAListIterator<BtaCandidate> iter_rho0(NewRho0List);  
  HepAListIterator<BtaCandidate> iter_rhoC(NewRhoCList);  
  HepAListIterator<BtaCandidate> iter_ome(NewOmegaList);  
  HepAListIterator<BtaCandidate> iter_gam(NewGammaList);  


  int piLevel=0,pi0Level=0,etaLevel=0,etapLevel=0,rhoCLevel=0,rho0Level=0,omeLevel=0,gamLevel=0;
  
  while ( 0 != ( lepton = iter_lep() ) ) {
   
    ///////////////////////////////////////////////////////////////////////    
    if( _analyzeChargedPi.value()==true ) {
      iter_pi.rewind();

      while ( 0 != ( hadron = iter_pi() ) ) {	
	

	if( !(hadron->overlaps(*lepton)) && lepton->charge()!=hadron->charge() ) 
	  { 
	    if( piLevel==0 ) { piLevel=1;  _NbPiAfterLeptonSelCuts+=1; }
	    
	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );

	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( piLevel==1 ) { piLevel=2;  _NbPiAfterCosBY+=1; }
		
		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_pi.value()*pXu;

		if( sum >= _pStar2Dsum_pi.value() || pXu>= _pXuUpsMin_pi.value() || pLep>= _pLepUpsMin_pi.value() )		
		  {
		    if( piLevel==2 ) { piLevel=3;  _NbPiAfterpStar2D+=1;  _NbPiAfterCoupleCuts+=1; }
		    
		    BtaOpMakeTree comb;
		    BtaCandidate* Y = comb.create(*hadron,*lepton);
		    if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B0));
		    else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::anti_B0));
		    else cout<<"Houston, we have a problem with the pi mode!"<<endl;
		    outYList->append( new BtaCandidate(*Y) );
		    _NbCouplePiAfterAllCuts+=1;
		    delete Y;
		  }//pStar2D
	      }//cosBY
	  }//opposite charge and no overlap
      }//while ( 0 != ( hadron = iter_pi() ) )
    }//if(_chargedPi)

    ///////////////////////////////////////////////////////////////////////
    if( _analyzeNeutralPi.value()==true ) {
      iter_pi0.rewind();

      while ( 0 != ( hadron = iter_pi0() ) ) {
		
	if( true )
	  {
	    if( pi0Level==0 ) { pi0Level=1;  _NbPi0AfterLeptonSelCuts+=1; }
	    
	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );
	    
	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;
	    
	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( pi0Level==1 ) { pi0Level=2;  _NbPi0AfterCosBY+=1; }
		
		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_pi0.value()*pXu;

		if( sum >= _pStar2Dsum_pi0.value() || pXu>= _pXuUpsMin_pi0.value() || pLep>= _pLepUpsMin_pi0.value() )		
		  {
		    if( pi0Level==2 ) { pi0Level=3;  _NbPi0AfterpStar2D+=1;   _NbPi0AfterCoupleCuts+=1; }
			
		    BtaOpMakeTree comb;
		    BtaCandidate* Y = comb.create(*hadron,*lepton);
		    if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
		    else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
		    else cout<<"Houston, we have a problem with the pi0 mode!"<<endl;
		    outYList->append( new BtaCandidate(*Y) );
		    _NbCouplePi0AfterAllCuts+=1;
		    delete Y;
		  }//pStar2D
	      }//cosBY
	  }//lepton selection
      }//while ( 0 != ( hadron = iter_pi0() ) )
    }//if(_analyzeNeutralPi)
    
    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeEta.value()==true ) {
      iter_eta.rewind();

      while ( 0 != ( hadron = iter_eta() ) ) {	
	
	//We don't want Eta's daughters to be the lepton!
	bool HadNotMixedToLepton = true;
	BtaCandidate* temp;
	HepAListIterator<BtaCandidate> fille = hadron->daughterIterator();
	while (  0 != ( temp = fille() ) )
	  {  if( temp->uid()==lepton->uid() ) HadNotMixedToLepton = false; }
	
	
	if( HadNotMixedToLepton )
	  {
	    if( etaLevel==0 ) { etaLevel=1;   _NbEtaAfterLeptonSelCuts+=1; }
	    
	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );
	    
	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( etaLevel==1 ) { etaLevel=2;  _NbEtaAfterCosBY+=1; }

		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_eta.value()*pXu;

		if( sum >= _pStar2Dsum_eta.value() || pXu>= _pXuUpsMin_eta.value() || pLep>= _pLepUpsMin_eta.value() )		
		  {
		    if( etaLevel==2 ) { etaLevel=3;  _NbEtaAfterpStar2D+=1;  _NbEtaAfterCoupleCuts+=1; }
		    
			BtaOpMakeTree comb;
			BtaCandidate* Y = comb.create(*hadron,*lepton);
			if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
			else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
			else cout<<"Houston, we have a problem with the eta mode!"<<endl;
			outYList->append( new BtaCandidate(*Y) );
			_NbCoupleEtaAfterAllCuts+=1;
			delete Y;
		  }//pStar2D
	      }//cosBY
	  }//if HadNotMixedToLepton
      }//while ( 0 != ( hadron = iter_eta() ) )
    }//if(_analyzeEta)


    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeEtap.value()==true ) {
      iter_etap.rewind();

      while ( 0 != ( hadron = iter_etap() ) ) {	
	
	//We don't want Etap's daughters to be the lepton!
	bool HadNotMixedToLepton = true;
	BtaCandidate* temp;
	HepAListIterator<BtaCandidate> fille = hadron->daughterIterator();
	while (  0 != ( temp = fille() ) )
	  {  if( temp->uid()==lepton->uid() ) HadNotMixedToLepton = false; }
	
	
	if( HadNotMixedToLepton )
	  {
	    if( etapLevel==0 ) { etapLevel=1;   _NbEtapAfterLeptonSelCuts+=1; }
	    
	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );
	    
	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( etapLevel==1 ) { etapLevel=2;  _NbEtapAfterCosBY+=1; }

		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_etap.value()*pXu;

		if( sum >= _pStar2Dsum_etap.value() || pXu>= _pXuUpsMin_etap.value() || pLep>= _pLepUpsMin_etap.value() )		
		  {
		    if( etapLevel==2 ) { etapLevel=3;  _NbEtapAfterpStar2D+=1;  _NbEtapAfterCoupleCuts+=1; }
		    
			BtaOpMakeTree comb;
			BtaCandidate* Y = comb.create(*hadron,*lepton);
			if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
			else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
			else cout<<"Houston, we have a problem with the etap mode!"<<endl;
			outYList->append( new BtaCandidate(*Y) );
			_NbCoupleEtapAfterAllCuts+=1;
			delete Y;
		  }//pStar2D
	      }//cosBY
	  }//if HadNotMixedToLepton
      }//while ( 0 != ( hadron = iter_etap() ) )
    }//if(_analyzeEtap)


    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeChargedRho.value()==true ) {
      iter_rhoC.rewind();

      while ( 0 != ( hadron = iter_rhoC() ) ) {	
	
	//We don't want Rho's daughters to be the lepton!
	bool HadNotMixedToLepton = true;
	BtaCandidate* temp;
	HepAListIterator<BtaCandidate> fille = hadron->daughterIterator();
	while (  0 != ( temp = fille() ) )
	  {  if( temp->uid()==lepton->uid() ) HadNotMixedToLepton = false; }

		
	if( lepton->charge()!=hadron->charge() && HadNotMixedToLepton )
	  {
	    if( rhoCLevel==0 ) { rhoCLevel=1;  _NbRhoCAfterLeptonSelCuts+=1; }

	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );

	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( rhoCLevel==1 ) { rhoCLevel=2;  _NbRhoCAfterCosBY+=1; }

		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_rhoC.value()*pXu;

		if( sum >= _pStar2Dsum_rhoC.value() || pXu>= _pXuUpsMin_rhoC.value() || pLep>= _pLepUpsMin_rhoC.value() )		
		  {
		    if( rhoCLevel==2 ) { rhoCLevel=3;  _NbRhoCAfterpStar2D+=1;   _NbRhoCAfterCoupleCuts+=1; }
		    
		    BtaOpMakeTree comb;
		    BtaCandidate* Y = comb.create(*hadron,*lepton);
		    if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B0));
		    else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::anti_B0));
		    else cout<<"Houston, we have a problem with the rhoC mode!"<<endl;
		    outYList->append( new BtaCandidate(*Y) );
		    _NbCoupleRhoCAfterAllCuts+=1;
		    delete Y;
		  }//pStar2D
	      }//cosBY
	  }//if opposite charges and HadNotMixedToLepton
      }//while ( 0 != ( hadron = iter_rhoC() ) )
    }//if(_analyzeChargedRho)

    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeNeutralRho.value()==true ) {
      iter_rho0.rewind();

      while ( 0 != ( hadron = iter_rho0() ) ) {	
	
	//We don't want Rho's daughters to be the lepton!
	bool HadNotMixedToLepton = true;
	BtaCandidate* temp;
	HepAListIterator<BtaCandidate> fille = hadron->daughterIterator();
	while (  0 != ( temp = fille() ) )
	  {  if( temp->uid()==lepton->uid() ) HadNotMixedToLepton = false; }

	
	if( HadNotMixedToLepton )
	  {
	    if( rho0Level==0 ) { rho0Level=1;  _NbRho0AfterLeptonSelCuts+=1; }

	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );

	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( rho0Level==1 ) { rho0Level=2;  _NbRho0AfterCosBY+=1; }

		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_rho0.value()*pXu;

		if( sum >= _pStar2Dsum_rho0.value() || pXu>= _pXuUpsMin_rho0.value() || pLep>= _pLepUpsMin_rho0.value() )		
		  {
		    if( rho0Level==2 ) { rho0Level=3;  _NbRho0AfterpStar2D+=1;  _NbRho0AfterCoupleCuts+=1; }
		    
		    BtaOpMakeTree comb;
		    BtaCandidate* Y = comb.create(*hadron,*lepton);
		    if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
		    else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
		    else cout<<"Houston, we have a problem with the rho0 mode!"<<endl;
		    outYList->append( new BtaCandidate(*Y) );
		    _NbCoupleRho0AfterAllCuts+=1;
		    delete Y;
		  }//pStar2D
	      }//cosBY
	  }//if HadNotMixedToLepton
      }//while ( 0 != ( hadron = iter_rho0() ) )
    }//if(_analyzeNeutralRho)

    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeOmega.value()==true ) {
      iter_ome.rewind();

      while ( 0 != ( hadron = iter_ome() ) ) {	
	
	//We don't want Omega's daughters to be the lepton!
	bool HadNotMixedToLepton = true;
	BtaCandidate* temp;
	HepAListIterator<BtaCandidate> fille = hadron->daughterIterator();
	while (  0 != ( temp = fille() ) )
	  {  if( temp->uid()==lepton->uid() ) HadNotMixedToLepton = false; }

	
	if( HadNotMixedToLepton )
	  {
	    if( omeLevel==0 ) { omeLevel=1;  _NbOmeAfterLeptonSelCuts+=1; }

	    BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	    BtaCandidate boostedHad = _CMS->boostTo( *hadron );

	    HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	    HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	    
	    double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	    double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	    //the abs is for the offres case. It has no effect otherwise!
	    double cosBY = numerat/denomin;

	    if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	      {
		if( omeLevel==1 ) { omeLevel=2;  _NbOmeAfterCosBY+=1; }

		double pLep = boostedLep.p();
		double pXu = boostedHad.p();
		double sum = pLep + _pStar2Dslope_omega.value()*pXu;

		if( sum >= _pStar2Dsum_omega.value() || pXu>= _pXuUpsMin_omega.value() || pLep>= _pLepUpsMin_omega.value() )		
		  {
		    if( omeLevel==2 ) { omeLevel=3;  _NbOmeAfterpStar2D+=1;   _NbOmeAfterCoupleCuts+=1; }
		    
		    BtaOpMakeTree comb;
		    BtaCandidate* Y = comb.create(*hadron,*lepton);
		    if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
		    else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
		    else cout<<"Houston, we have a problem with the omega l nu mode!"<<endl;
		    outYList->append( new BtaCandidate(*Y) );
		    _NbCoupleOmegaAfterAllCuts+=1;
		    delete Y;
		  }//pStar2D
	      }//cosBY
	  }//if HadNotMixedToLepton
      }//while ( 0 != ( hadron = iter_ome() ) )
    }//if(_analyzeOmega)

    ///////////////////////////////////////////////////////////////////////////////
    if( _analyzeGamma.value()==true ) {
      iter_gam.rewind();
      
      while ( 0 != ( hadron = iter_gam() ) ) {	
	
	  if( gamLevel==0 ) { gamLevel=1;  _NbGammaAfterLeptonSelCuts+=1; }

	  BtaCandidate boostedLep = _CMS->boostTo( *lepton );
	  BtaCandidate boostedHad = _CMS->boostTo( *hadron );

	  HepLorentzVector boostedP4Y = boostedLep.p4() + boostedHad.p4();
	  HepLorentzVector P4Y = lepton->p4() + hadron->p4();
	  
	  double numerat = sqrt(s)*boostedP4Y.e()-(mB0*mB0)- P4Y.mag2();  
	  double denomin = 2*sqrt( abs(s/4-(mB0*mB0)) )*boostedP4Y.vect().mag(); 
	  //the abs is for the offres case. It has no effect otherwise!
	  double cosBY = numerat/denomin;

	  if( cosBY>=_cosBYmin.value() &&  cosBY<=_cosBYmax.value() ) 
	    {
	      if( gamLevel==1 ) { gamLevel=2;  _NbGammaAfterCosBY+=1;  _NbGammaAfterCoupleCuts+=1; }
	      
	      BtaOpMakeTree comb;
	      BtaCandidate* Y = comb.create(*hadron,*lepton);
	      if(lepton->charge()==1) Y->setType(Pdt::lookup(PdtLund::B_plus));
	      else if(lepton->charge()==-1) Y->setType(Pdt::lookup(PdtLund::B_minus));
	      else cout<<"Houston, we have a problem with the gamma l nu mode!"<<endl;
	      outYList->append( new BtaCandidate(*Y) );
	      _NbCoupleGammaAfterAllCuts+=1;
	      delete Y;
	    }//cosBY
      }//while ( 0 != ( hadron = iter_gam() ) )
    }//if(_analyzeGamma)
    
    
  }//while ( 0 != ( lepton = iter_lep() ) )

  if( outYList->length()!=0 ) { _NbAllModesAfterCoupleCuts+=1; _NbCoupleAllModesAfterAllCuts+=outYList->length(); }


  //Testing (voir BruCo.cc pour faire un bloc "rholnu", ce qui serait surement pertinent! ;-) ): 
  _TestingNtuple->column("NbCoupleOK",outYList->length());
  _TestingNtuple->dumpData();


  ///////////////////
  //  Final result
  ///////////////////
  bool isPi=false,isPi0=false,isEta=false,isEtap=false,isRhoC=false,isRho0=false,isOmega=false,isGam=false;
  if(piLevel==3) isPi=true;
  if(pi0Level==3) isPi0=true;
  if(etaLevel==3) isEta=true;
  if(etapLevel==3) isEtap=true;
  if(rhoCLevel==3) isRhoC=true;
  if(rho0Level==3) isRho0=true;
  if(omeLevel==3) isOmega=true;
  if(gamLevel==2) isGam=true;

  setPassed(isPi||isPi0||isEta||isEtap||isRhoC||isRho0||isOmega||isGam); 

  //////////////////////////////////////
  //  user data part (using bit mask):
  //////////////////////////////////////
  int i=0;
  if ( isPi )    i = i | 1; 
  if ( isPi0 )   i = i | 2; 
  if ( isEta )   i = i | 4; 
  if ( isEtap )  i = i | 8; 
  if ( isRhoC )  i = i | 16; 
  if ( isRho0 )  i = i | 32; 
  if ( isOmega ) i = i | 64; 
  if ( isGam )   i = i | 128; 
  _decayMode=i;


  UsrIfd<UsrEventBlock>::put( anEvent, &_eventBlock, _eventData.value() );
  bool ok = true;
  ok &= _eventBlock.put( _decayMode );
  assert( ok );

  return QuitEvent(); 

}


bool
XSLBtoXulnuFilter::SurvivedTrigger(AbsEvent* anEvent)
{
  bool evtOK=false;
  float trigOK=-1;

  AbsEventTag* tag = Ifd<AbsEventTag>::get( anEvent );
  if(tag==0) 
    { 
      cout<<"NOTHING IN AbsEvent!"<<endl;
    }
  else
    {
      bool L3EMC,L3DCH;
      if(tag->getBool(L3EMC,"L3OutEmc")
	 &&tag->getBool(L3DCH,"L3OutDch"))	
	{
	  if(L3EMC||L3DCH) trigOK=1;
	  else trigOK=0;
	}
      else ErrMsg(fatal) << "Could not find L3 Trigger bits!";
    }
  
  if( trigOK>= _trigMIN.value() ) evtOK=true;
  
  return evtOK;
}

bool
XSLBtoXulnuFilter::SurvivedISH(AbsEvent* anEvent)
{
  bool evtOK=false;
  float ish=-1;

  AbsEventTag* tag = Ifd<AbsEventTag>::get( anEvent );
  if(tag==0) 
    { 
      cout<<"NOTHING IN AbsEvent!"<<endl;
    }
  else
    {
      bool ismultihad;
      if(tag->getBool(ismultihad,"BGFMultiHadron"))	{
	if(ismultihad) ish=1;
	else ish=0;
      }
      else ErrMsg(fatal) << "Could not find BGFMultiHadron";
    }
  
  if( ish>= _ishMIN.value() ) evtOK=true;
  
  return evtOK;
}

bool
XSLBtoXulnuFilter::SurvivedR2All(AbsEvent* anEvent)
{
  bool evtOK=false;
  float r2all=2;

  AbsEventTag* tag = Ifd<AbsEventTag>::get( anEvent );
  if(tag==0) 
    { 
      cout<<"NOTHING IN AbsEvent!"<<endl;
    }
  else
    {
      if(tag->getFloat(r2all,"R2All"));
      else ErrMsg(fatal) << "Could not find R2All bit!";
    }  
  
  if( r2all<=_r2allMAX.value() ) evtOK=true;

  return evtOK;
}


void
XSLBtoXulnuFilter::MakeHadronsAndLeptonsList(AbsEvent* anEvent)
{
  //We apply here additional kinematic and quality cuts to the various lists
  //Note that the base list are already somewhat customized in various CompositionSequences modules
  BtaCandidate* candid;  


  //PID Lists   //PidLHElectrons, piLHLoose, cutBased muSelector
  getTmpAList (anEvent, _eTightLH,  IfdStrKey(HepString("PidLHElectrons")) );
  if (_eTightLH == 0)  ErrMsg(fatal) << "Could not locate list eTightLH" << endmsg;

  getTmpAList (anEvent, _piLooseLH,  IfdStrKey(HepString("piLHLoose")) );
  if (_piLooseLH == 0)  ErrMsg(fatal) << "Could not locate list piLHLoose" << endmsg;

  getTmpAList (anEvent, _muTight,  IfdStrKey(HepString("muMicroTight")) );
  if (_muTight == 0)  ErrMsg(fatal) << "Could not locate list muMicroTight" << endmsg;
  getTmpAList (anEvent, _muVeryTight,  IfdStrKey(HepString("muMicroVeryTight")) );
  if (_muVeryTight == 0)  ErrMsg(fatal) << "Could not locate list muMicroVeryTight" << endmsg;


  ///////////////
  //  LEPTONS  //
  ///////////////
  getTmpAList(anEvent,InputLeptonList,_InputLeptonList.value());
  if (InputLeptonList == 0)   ErrMsg(fatal)<<"No leptonList!!"<<endmsg;
  HepAListIterator<BtaCandidate> iter_InputLepton(*InputLeptonList);  

  getTmpAList(anEvent,NewLepList,IfdStrKey(HepString("NewLepList")));
  

  while ( 0 != ( candid = iter_InputLepton()) ) 
    {
      bool ePID = electronPID( candid );
      bool muPID = muonPID( candid );

      bool eCut=false;
      bool muCut=false;
      if(candid->p() >= _pLABmin_electron.value() && ePID ) eCut=true;
      if(candid->p() >= _pLABmin_muon.value() && muPID ) muCut=true;
      
      if( ( eCut || muCut)	
	  && candid->p()<10  
	  && candid->p3().theta() >= _thetaLABmin_electron.value()
	  && candid->p3().theta() <= _thetaLABmax_electron.value() )  
	{
	  double ch = candid->charge();
	  BtaCandidate* newLep = new BtaCandidate(*candid);
	  if(eCut && ch>0) newLep->setType(Pdt::lookup(PdtLund::e_plus));
	  else if(eCut && ch<0) newLep->setType(Pdt::lookup(PdtLund::e_minus));
	  else if(muCut && ch>0) newLep->setType(Pdt::lookup(PdtLund::mu_plus));
	  else if(muCut && ch<0) newLep->setType(Pdt::lookup(PdtLund::mu_minus));
	  else  ErrMsg(fatal)<<"Something's wrong when trying to build NewLep..."<<endmsg;

	  NewLepList->append( newLep );
	 }
    }


  ///////////////
  //  HADRONS  //
  ///////////////

  //Put Hadron Selection cuts here! 

  if(_analyzeChargedPi.value()==true)
    {
      getTmpAList(anEvent,InputPiList,_InputPiList.value());
      if (InputPiList == 0)   ErrMsg(fatal)<<"No piList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputPi(*InputPiList);  
      while ( 0 != ( candid = iter_InputPi()) ) 
	{
	  bool daughtersOK=false;
	  if( _doPionPidYourself.value()==false) daughtersOK=true;
	  else daughtersOK = pionPID(candid);

	  if( daughtersOK && candid->nDaughters()==0 ) NewPiList.append( candid ); 
	}
    }

  if(_analyzeNeutralPi.value()==true)
    {
      getTmpAList(anEvent,InputPi0List,_InputPi0List.value());
      if (InputPi0List == 0)   ErrMsg(fatal)<<"No pi0List!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputPi0(*InputPi0List);  
      while ( 0 != ( candid = iter_InputPi0()) ) 
	{
	  if( candid->p() <= _pLABmax_pi0.value() ) NewPi0List.append( candid );
	}
    }

  if(_analyzeEta.value()==true)
    {
      getTmpAList(anEvent,InputEtaList,_InputEtaList.value());
      if (InputEtaList == 0)   ErrMsg(fatal)<<"No etaList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputEta(*InputEtaList);  
      while ( 0 != ( candid = iter_InputEta()) ) 
	{
	  bool daughtersOK=false;
	  if(candid->nDaughters()==2 || _doPionPidYourself.value()==false ) daughtersOK=true; //eta-> gam gam mode
	  else if( candid->nDaughters()==3 && _doPionPidYourself.value()==true ) //eta->pi+ pi- pi0 mode
	    {
	      BtaCandidate* cand;
	      int g=0;
	      HepAListIterator<BtaCandidate>iter_fille = candid->daughterIterator();
	      //We ask the two charged pions daughters to pass the piLH selector
	      while ( 0 != (cand =iter_fille()) && g>=0 ) 
		{ if(cand->charge()!=0 && !pionPID(cand) ) g=-1; }
	      if(g==0) daughtersOK=true;
	    }

	  if( daughtersOK && candid->p()<=_pLABmax_eta.value() ) NewEtaList.append( candid );
	}
    }

  if(_analyzeEtap.value()==true)
    {
      getTmpAList(anEvent,InputEtapList,_InputEtapList.value());
      if (InputEtapList == 0)   ErrMsg(fatal)<<"No etapList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputEtap(*InputEtapList);  
      while ( 0 != ( candid = iter_InputEtap()) ) 
	{
	  if( candid->p()<=_pLABmax_etap.value() ) NewEtapList.append( candid );
	}
    }

  if(_analyzeNeutralRho.value()==true)
    {
      getTmpAList(anEvent,InputRho0List,_InputRho0List.value());
      if (InputRho0List == 0)   ErrMsg(fatal)<<"No rho0List!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputRho0(*InputRho0List);  
      while ( 0 != ( candid = iter_InputRho0()) ) 
	{
	  bool daughtersOK=false;

	  if(_doPionPidYourself.value()==false) daughtersOK=true;
	  else 
	    {
	      BtaCandidate* cand;
	      int g=0;
	      HepAListIterator<BtaCandidate>iter_fille = candid->daughterIterator();
	      //We ask the two charged pions daughters to pass the piLH selector
	      while ( 0 != (cand =iter_fille()) && g>=0 ) 
		{ if(cand->charge()!=0 && !pionPID(cand) ) g=-1; }
	      if(g==0) daughtersOK=true;
	    }

	  if( daughtersOK && candid->p() <= _pLABmax_rho0.value() ) NewRho0List.append( candid );
	}
    }

  if(_analyzeChargedRho.value()==true)
    {
      getTmpAList(anEvent,InputRhoCList,_InputRhoCList.value());
      if (InputRhoCList == 0)   ErrMsg(fatal)<<"No rhoCList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputRhoC(*InputRhoCList);  
      while ( 0 != ( candid = iter_InputRhoC()) ) 
	{
	  bool daughtersOK=false;

	  if(_doPionPidYourself.value()==false) daughtersOK=true;
	  else 
	    {
	      BtaCandidate* cand;
	      int g=0;
	      HepAListIterator<BtaCandidate>iter_fille = candid->daughterIterator();
	      //We ask the two charged pions daughters to pass the piLH selector
	      while ( 0 != (cand =iter_fille()) && g>=0 ) 
		{ if(cand->charge()!=0 && !pionPID(cand) ) g=-1; }
	      if(g==0) daughtersOK=true;
	    }

	  if( daughtersOK && candid->p()<=_pLABmax_rhoC.value() ) NewRhoCList.append( candid );
	}
    }

  if(_analyzeOmega.value()==true)
    {
      getTmpAList(anEvent,InputOmegaList,_InputOmegaList.value());
      if (InputOmegaList == 0)   ErrMsg(fatal)<<"No omegaList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputOmega(*InputOmegaList);  
      while ( 0 != ( candid = iter_InputOmega()) ) 
	{
	  bool daughtersOK=false;

	  if(_doPionPidYourself.value()==false) daughtersOK=true;
	  else 
	    {
	      BtaCandidate* cand;
	      int g=0;
	      HepAListIterator<BtaCandidate>iter_fille = candid->daughterIterator();
	      //We ask the two charged pions daughters to pass the piLH selector
	      while ( 0 != (cand =iter_fille()) && g>=0 ) 
		{ if(cand->charge()!=0 && !pionPID(cand) ) g=-1; }
	      if(g==0) daughtersOK=true;
	    }

	  if( daughtersOK && candid->p()<=_pLABmax_omega.value() ) NewOmegaList.append( candid );
	}
    }
  
  if(_analyzeGamma.value()==true)
    {
      getTmpAList(anEvent,InputGammaList,_InputGammaList.value());
      if (InputGammaList == 0)   ErrMsg(fatal)<<"No gammaList!!"<<endmsg;
      HepAListIterator<BtaCandidate> iter_InputGamma(*InputGammaList);  
      while ( 0 != ( candid = iter_InputGamma()) ) 
	{
	  if( candid->p() <= _pLABmax_gamma.value() ) NewGammaList.append( candid );
	}
    } // if _analyzeGamma
}



AppResult
XSLBtoXulnuFilter::endJob( AbsEvent* anEvent )
{

  if(_MyVerbose.value())
    {	
      //Setting the precision of the numbers to be printed to a higher value...
      ostream& os = ErrMsg(routine);
      int p=os.precision(); 
      os<<setprecision(8);

      cout<<"Results from B->Xu l nu XSL skim"<<endl;
      cout<<endl;
      cout<<"SUMMARY:"<<endl;
      cout<<"Total number of events before any cut:              "<<_TotalNbEvents<<endl;
      cout<<"Events passing Event cuts:                          "<<_NbAfterEventCuts<<" --> "<<_NbAfterEventCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (All Modes):             "<<_NbAllModesAfterCoupleCuts<<" --> "<<_NbAllModesAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Pi Mode):               "<<_NbPiAfterCoupleCuts<<" --> "<<_NbPiAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Pi0 Mode):              "<<_NbPi0AfterCoupleCuts<<" --> "<<_NbPi0AfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Eta Mode):              "<<_NbEtaAfterCoupleCuts<<" --> "<<_NbEtaAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Eta' Mode):             "<<_NbEtapAfterCoupleCuts<<" --> "<<_NbEtapAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Rho0 Mode):             "<<_NbRho0AfterCoupleCuts<<" --> "<<_NbRho0AfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (RhoC Mode):             "<<_NbRhoCAfterCoupleCuts<<" --> "<<_NbRhoCAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Omega Mode):            "<<_NbOmeAfterCoupleCuts<<" --> "<<_NbOmeAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cuts (Gamma Mode):            "<<_NbGammaAfterCoupleCuts<<" --> "<<_NbGammaAfterCoupleCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"BtaCandidates INFO:"<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (All Modes):  "<<_NbCoupleAllModesAfterAllCuts<<" --> "<<_NbCoupleAllModesAfterAllCuts/_NbAllModesAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (pi Mode):    "<<_NbCouplePiAfterAllCuts<<" --> "<<_NbCouplePiAfterAllCuts/_NbPiAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (pi0 Mode):   "<<_NbCouplePi0AfterAllCuts<<" --> "<<_NbCouplePi0AfterAllCuts/_NbPi0AfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (eta Mode):   "<<_NbCoupleEtaAfterAllCuts<<" --> "<<_NbCoupleEtaAfterAllCuts/_NbEtaAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (eta' Mode):   "<<_NbCoupleEtapAfterAllCuts<<" --> "<<_NbCoupleEtapAfterAllCuts/_NbEtapAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (rhoC Mode):  "<<_NbCoupleRhoCAfterAllCuts<<" --> "<<_NbCoupleRhoCAfterAllCuts/_NbRhoCAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (rho0 Mode):  "<<_NbCoupleRho0AfterAllCuts<<" --> "<<_NbCoupleRho0AfterAllCuts/_NbRho0AfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (omega Mode): "<<_NbCoupleOmegaAfterAllCuts<<" --> "<<_NbCoupleOmegaAfterAllCuts/_NbOmeAfterCoupleCuts<<endl;
      cout<<"Nb of BtaCandidates/accepted event after all cuts (gamma Mode): "<<_NbCoupleGammaAfterAllCuts<<" --> "<<_NbCoupleGammaAfterAllCuts/_NbGammaAfterCoupleCuts<<endl;
      cout<<endl;
      cout<<"CUT BY CUT RESULTS:"<<endl;
      cout<<"Events passing Trigger:                             "<<_NbAfterTrigger<<" --> "<<_NbAfterTrigger/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing ISH:                                 "<<_NbAfterISHCuts<<" --> "<<_NbAfterISHCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing R2All:                               "<<_NbAfterR2AllCuts<<" --> "<<_NbAfterR2AllCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing >=1 lepton tight:                    "<<_NbAfterOneLeptonTight<<" --> "<<_NbAfterOneLeptonTight/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Pi Mode):    "<<_NbPiAfterLeptonSelCuts<<" --> "<<_NbPiAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Pi Mode):          "<<_NbPiAfterCosBY<<" --> "<<_NbPiAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Pi Mode):        "<<_NbPiAfterpStar2D<<" --> "<<_NbPiAfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Pi0 Mode):   "<<_NbPi0AfterLeptonSelCuts<<" --> "<<_NbPi0AfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Pi0 Mode):         "<<_NbPi0AfterCosBY<<" --> "<<_NbPi0AfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Pi0 Mode):       "<<_NbPi0AfterpStar2D<<" --> "<<_NbPi0AfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Eta Mode):   "<<_NbEtaAfterLeptonSelCuts<<" --> "<<_NbEtaAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Eta Mode):         "<<_NbEtaAfterCosBY<<" --> "<<_NbEtaAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Eta Mode):       "<<_NbEtaAfterpStar2D<<" --> "<<_NbEtaAfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Eta' Mode):  "<<_NbEtapAfterLeptonSelCuts<<" --> "<<_NbEtapAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Eta' Mode):        "<<_NbEtapAfterCosBY<<" --> "<<_NbEtapAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Eta' Mode):      "<<_NbEtapAfterpStar2D<<" --> "<<_NbEtapAfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Rho0 Mode):  "<<_NbRho0AfterLeptonSelCuts<<" --> "<<_NbRho0AfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Rho0 Mode):        "<<_NbRho0AfterCosBY<<" --> "<<_NbRho0AfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Rho0 Mode):      "<<_NbRho0AfterpStar2D<<" --> "<<_NbRho0AfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (RhoC Mode):  "<<_NbRhoCAfterLeptonSelCuts<<" --> "<<_NbRhoCAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (RhoC Mode):        "<<_NbRhoCAfterCosBY<<" --> "<<_NbRhoCAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (RhoC Mode):      "<<_NbRhoCAfterpStar2D<<" --> "<<_NbRhoCAfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Omega Mode): "<<_NbOmeAfterLeptonSelCuts<<" --> "<<_NbOmeAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Omega Mode):       "<<_NbOmeAfterCosBY<<" --> "<<_NbOmeAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple pStar2D cut (Omega Mode):     "<<_NbOmeAfterpStar2D<<" --> "<<_NbOmeAfterpStar2D/_TotalNbEvents*100<<"%"<<endl;
      cout<<endl;
      cout<<"Events passing Couple lepton Sel cuts (Gamma Mode): "<<_NbGammaAfterLeptonSelCuts<<" --> "<<_NbGammaAfterLeptonSelCuts/_TotalNbEvents*100<<"%"<<endl;
      cout<<"Events passing Couple cosBY cut (Gamma Mode):       "<<_NbGammaAfterCosBY<<" --> "<<_NbGammaAfterCosBY/_TotalNbEvents*100<<"%"<<endl;
      cout<<"No pStar2D cut for the Gamma Mode"<<endl;
      cout<<endl;
      cout<<"C't'a peu pres ca pour a soir... Byebye... Byebye..."<<endl;

      //...setting it back to its original precision.
      os<<setprecision(p)<<endmsg;
    }
  
  return AppResult::OK;
  
}

//////////////

AppResult
XSLBtoXulnuFilter::endRun( AbsEvent* anEvent ) {    
  return AppResult::OK;
}
 
bool XSLBtoXulnuFilter::pionPID( BtaCandidate* candid )
{
  BtaCandidate* PIDcand;
  
  HepAListIterator<BtaCandidate> piL(*_piLooseLH);
  while( PIDcand = piL() ){
    if( PIDcand->overlaps(*candid) )
      { 
	return true;
      }
  }
  
  return false;  
}

bool XSLBtoXulnuFilter::electronPID( BtaCandidate* candid )
{
  BtaCandidate* PIDcand;
  
  HepAListIterator<BtaCandidate> eT(*_eTightLH);
  while( PIDcand = eT() ){
    if( PIDcand->overlaps(*candid) )
      { 
	return true;
      }
  }
  
  return false;  
}

bool XSLBtoXulnuFilter::muonPID( BtaCandidate* candid )
{
  BtaCandidate* PIDcand;
  
  HepAListIterator<BtaCandidate> muT(*_muTight);
  while( PIDcand = muT() ){
    if( PIDcand->overlaps(*candid) )
      { 
	return true;
      }
  }

  HepAListIterator<BtaCandidate> muVT(*_muVeryTight);
  while( PIDcand = muVT() ){
    if( PIDcand->overlaps(*candid) )
      { 
	return true;
      }
  }

    
  return false;
}

// clone
//
AppModule* 
XSLBtoXulnuFilter::clone(const char* cloneName)
{
  return new XSLBtoXulnuFilter(cloneName, "clone of XSLBtoXulnuFilter");
}


AppResult
XSLBtoXulnuFilter::QuitEvent()
{
  //CLEANING the list before quiting the event...
   
  //NewLepList.removeAll();
  NewPiList.removeAll();
  NewPi0List.removeAll();
  NewEtaList.removeAll();
  NewEtapList.removeAll();
  NewOmegaList.removeAll();
  NewRho0List.removeAll();
  NewRhoCList.removeAll();
  NewGammaList.removeAll();
  
  //HepAListDeleteAll( NewLepList );
  HepAListDeleteAll( NewPiList );
  HepAListDeleteAll( NewPi0List );
  HepAListDeleteAll( NewEtaList );
  HepAListDeleteAll( NewEtapList );
  HepAListDeleteAll( NewOmegaList );
  HepAListDeleteAll( NewRho0List );
  HepAListDeleteAll( NewRhoCList );
  HepAListDeleteAll( NewGammaList );
  
  delete _CMS;

  return AppResult::OK;

}