//_____________________________________________________________________________
// StntupleFilterModule.cc (by P.Murat)
// Functionality:  provides various event filters which can be used for
//                 creating secondary datasets
//-----------------------------------------------------------------------------
#include "Edm/EventRecord.hh"
#include "Edm/Handle.hh"

#include <iostream>
#include <iomanip>
#include <math.h>

#include "TFolder.h"
#include <AbsEnv/AbsEnv.hh>


#include "CalorObjects/PhysicsTowerParams.hh"
#include "CalorObjects/PhysicsTowerData.hh"
#include "Calor/PhysicsTowerDataMaker.hh"
#include "Calor/PhysicsTowerDataMaker.hh"

#include "MetObjects/CdfMet.hh"
#include "Met/NaiveEtParams.hh"
#include "Met/NaiveEtCalculator.hh"


#include <Stntuple/obj/TClcDataBlock.hh>
#include <Stntuple/obj/TCmuDataBlock.hh>
#include <Stntuple/obj/TCmpDataBlock.hh>
#include <Stntuple/obj/TCmxDataBlock.hh>
#include <Stntuple/obj/TCalDataBlock.hh>
#include <Stntuple/obj/TStnElectronBlock.hh>
#include <Stntuple/obj/TStnErrorLogger.hh>
#include <Stntuple/obj/TStnJetBlock.hh>
#include <Stntuple/obj/TStnMuonBlock.hh>
#include <Stntuple/obj/TStnMetBlock.hh>
#include <Stntuple/obj/TStnTauBlock.hh>
#include <Stntuple/obj/TStnTrackBlock.hh>
#include <Stntuple/obj/TStnVertexBlock.hh>
#include <Stntuple/obj/TStnTriggerBlock.hh>
#include <Stntuple/obj/TGenpBlock.hh>
#include <Stntuple/obj/TObspBlock.hh>
#include <Stntuple/obj/TFwdDetDataBlock.hh>

#include "Stntuple/alg/TStntuple.hh"
#include "Stntuple/alg/TStnMuonID.hh"

#include <Stntuple/mod/StntupleBeamManager.hh>
#include <Stntuple/mod/InitStntupleDataBlocks.hh>
#include <Stntuple/mod/StntupleUtilities.hh>

#include "Stntuple/mod/StntupleModule.hh"
#include "Stntuple/mod/StntupleEleFilter.hh"
#include "Stntuple/mod/StntupleMuonFilter.hh"
#include "Stntuple/mod/StntupleTrackFilter.hh"
#include "Stntuple/mod/StntupleZmumuFilter.hh"

#include "Stntuple/mod/StntupleFilterModule.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/obj/TStnNode.hh"

ClassImp(StntupleFilterModule)

//_____________________________________________________________________________
StntupleFilterModule::StntupleFilterModule( const char* const theName,
					    const char* const theDescription )
  : StntupleModule( theName, theDescription )

  ,_filteringMode    ("filteringMode"   ,this,"and")
  , fDebug           ("debug"           ,this, 0)
  , fProcessName     ("processName"     ,this,"PROD")
  , fBeamHistoryCode ("beamHistoryCode" ,this, -1, -1, 1000)

  , fJetCollName     ("jetCollName"     ,this,"JetCluModule")
  , fElectronCollName("electronCollName",this,"CalDriven")
  , fMuonCollName    ("muonCollName"    ,this,"ProductionMuons")
  , fTrackCollName   ("trackCollName"   ,this,"default")
  , fVertexCollName  ("vertexCollName"  ,this,"DEFAULT_VXPRIM_VERTICES")
  , fMetName         ("metName"         ,this,"MetDefault")

  ,_clcZ0Cut       ("clcZ0Cut"       ,this,10000.)
  ,_clcTMin        ("clcTMin"        ,this,    0.)
  ,_clcTMax        ("clcTMax"        ,this, 2000.)
  ,_clcAdcThreshold("clcAdcThreshold",this,  300.)
  
  ,_jetEtMin       ("jetEtMin"       ,this,   10.)
  ,_jetEtaMin      ("jetEtaMin"      ,this,  -10.,-10.,10)
  ,_jetEtaMax      ("jetEtaMax"      ,this,   10.,-10.,10)
  ,_jetNMin        ("jetNMin"        ,this,    1 )
  
  ,_monojetEt1     ("monojetEt1"     ,this,   20.)
  ,_monojetEt2     ("monojetEt2"     ,this,   10.)
  ,_monojetNTracks ("monojetNTracks" ,this,   0  )
  ,_monojetTrackPt ("monojetTrackPt" ,this,   2. )
  ,_monojetEtaMin  ("monojetEtaMin"  ,this,  -1.,-10.,10)
  ,_monojetEtaMax  ("monojetEtaMax"  ,this,   1.,-10.,10)
  ,_monojetCone    ("monojetCone"    ,this,   0.4,0.,10)
  ,_monojetDPhi    ("monojetDPhi"    ,this,   -1.,-1.,10)
  
  ,_metMin         ("metMin"         ,this,   20.)
  ,_metMode        ("metMode"        ,this,    0 )

  ,_tauPtMin       ("tauPtMin"       ,this,    5.)
  ,_tauPtMax       ("tauPtMax"       ,this, 1.e5 )
  ,_tauEmfrMin     ("tauEmfrMin"     ,this,    0.)
  ,_tauEmfrMax     ("tauEmfrMax"     ,this,    1.)
  ,_tauNMin        ("tauNMin"        ,this,    1 )
  
  ,_cmuNDHitsMin   ("cmuNDHitsMin"   ,this,    0 )
  ,_cmuNEHitsMin   ("cmuNEHitsMin"   ,this,    3 )
  ,_cmuRequireTrack("cmuRequireTrack",this,    0 )
  ,_cmuTrackPtMin  ("cmuTrackPtMin"  ,this,    false)
  ,_cmuDeltaPhiMax ("cmuDeltaPhiMax" ,this,    0.4)
  
  ,_cmpNCspHitsMin ("cmpNCspHitsMin" ,this,    0 )
  
  , _L1Bits         ("L1Bits"         ,this, 0,100,0)
  , _MyronBucket    ("myronBucket"    ,this, -1)
  , _L3Paths        ("L3Paths"        ,this, 0,100,0)

  , _genpNMin       ("genpNMin"       ,this,   0 )
  , _obspNMin       ("obspNMin"       ,this,   0 )
  
{
//-----------------------------------------------------------------------------
// talk-to parameters 
//-----------------------------------------------------------------------------
  _filteringMode.addDescription(
" \t and or or");
  commands()->append(&_filteringMode);

  fDebug.addDescription(
" \t integer level of debugging, default: 0");
  commands()->append(&fDebug);

  fProcessName.addDescription(
" \t process name, default: PROD");
  commands()->append(&fProcessName);

  commands()->append(&fBeamHistoryCode);
  fBeamHistoryCode.addDescription("\tPreferred beamline history code");
//-----------------------------------------------------------------------------
//  defaults for collection names
//-----------------------------------------------------------------------------
  commands()->append(&fJetCollName);
  commands()->append(&fMuonCollName);
  commands()->append(&fTrackCollName);
  commands()->append(&fElectronCollName);
  commands()->append(&fMetName);
//-----------------------------------------------------------------------------
//  CLC cuts
//-----------------------------------------------------------------------------
  _clcZ0Cut.addDescription(
" \t Z(max) : events with |Z| >= Z0Cut are rejected");
  commands()->append(&_clcZ0Cut);

  _clcTMin.addDescription(
" \t time averaging procedure uses times tMin < T < tMax");
  commands()->append(&_clcTMin);
 
  _clcTMax.addDescription(
" \t time averaging procedure uses times tMin < T < tMax");
  commands()->append(&_clcTMax);

  _clcAdcThreshold.addDescription(
" \t time averaging procedure uses channels with A(Adc) > adcThreshold");
  commands()->append(&_clcAdcThreshold);
//-----------------------------------------------------------------------------
//  jet cuts
//-----------------------------------------------------------------------------
  _jetEtMin.addDescription(" \t min Et of the jet, default = 10 GeV");
  commands()->append(&_jetEtMin);

  _jetEtaMin.addDescription(" \t min Eta of the jet, default = -10");
  commands()->append(&_jetEtaMin);

  _jetEtaMax.addDescription(" \t max Eta of the jet, default = 10");
  commands()->append(&_jetEtaMax);

  _jetNMin.addDescription(" \t min N(jets), default = 0");
  commands()->append(&_jetNMin);
//-----------------------------------------------------------------------------
//  monojet cuts
//-----------------------------------------------------------------------------
  _monojetEt1.addDescription(" \t min Et of the monojet, default = 20 GeV");
  commands()->append(&_monojetEt1);

  _monojetEt2.addDescription(" \t max Et of the other jet, default = 10 GeV");
  commands()->append(&_monojetEt2);

  _monojetEtaMin.addDescription(" \t min Eta of the monojet, default = -1");
  commands()->append(&_monojetEtaMin);

  _monojetEtaMax.addDescription(" \t max Eta of the monojet, default = 1");
  commands()->append(&_monojetEtaMax);

  _monojetNTracks.addDescription(" \t min number of monojet tracks, default = 1");
  commands()->append(&_monojetNTracks);

  _monojetTrackPt.addDescription(" \t min pt of monojet tracks, default = 2 GeV");
  commands()->append(&_monojetTrackPt);

  _monojetCone.addDescription(" \t cone within the tracks are searched, default = 0.4");
  commands()->append(&_monojetCone);

  _monojetDPhi.addDescription(" \t min dphi) between MET and the other jet, dflt = -1");
  commands()->append(&_monojetDPhi);
//-----------------------------------------------------------------------------
//  met cuts
//-----------------------------------------------------------------------------
  _metMin.addDescription(" \t min MET, default = 20 GeV");
  commands()->append(&_metMin);

  _metMode.addDescription(" \t way of calculating the MET, default = 0");
  commands()->append(&_metMode);
//-----------------------------------------------------------------------------
//  tau cuts
//-----------------------------------------------------------------------------
  _tauPtMin.addDescription(" \t min Pt of the tau candidate, default = 5 GeV");
  commands()->append(&_tauPtMin);

  _tauPtMax.addDescription(" \t max Pt of the tau candidate, default = 1.e5 GeV");
  commands()->append(&_tauPtMax);

  _tauNMin.addDescription(" \t min N(tau candidates) with Pt > tauPtMin, default = 1");
  commands()->append(&_tauNMin);

  _tauEmfrMin.addDescription(" \t min EM(frackion of the tau cluster, default = 0");
  commands()->append(&_tauEmfrMin);

  _tauEmfrMax.addDescription(" \t max EM(frackion of the tau cluster, default = 1");
  commands()->append(&_tauEmfrMax);
//-----------------------------------------------------------------------------
//  CMU cuts
//-----------------------------------------------------------------------------
  _cmuNDHitsMin.addDescription(" \t min N(CMU D hits), default = 0");
  commands()->append(&_cmuNDHitsMin);

  _cmuNEHitsMin.addDescription(" \t min N(CMU E hits), default = 0");
  commands()->append(&_cmuNEHitsMin);

  _cmuRequireTrack.addDescription(" \t if true, require track (default = false)");
  commands()->append(&_cmuRequireTrack);

  _cmuTrackPtMin.addDescription(" \t min Pt of a muon track(default = 3 GeV)");
  commands()->append(&_cmuTrackPtMin);

  _cmuDeltaPhiMax.addDescription(
"         max delta(phi) in radians between the track momentum in the vertex \n\
          and the center of the wedge (default = 0.4)");
  commands()->append(&_cmuDeltaPhiMax);
//-----------------------------------------------------------------------------
//  GENP/OBSP cuts
//-----------------------------------------------------------------------------
  _genpNMin.addDescription(" \t min N(GENP part) passing the cuts, dflt = 0");
  commands()->append(&_genpNMin);

  _obspNMin.addDescription(" \t min N(OBSP part) passing the cuts, dflt = 0");
  commands()->append(&_obspNMin);
//-----------------------------------------------------------------------------
//  create local data blocks
//-----------------------------------------------------------------------------
  TStnNode*  node;

  fEvent            = new TStnEvent();
  fEvent->SetErrorLogger(ErrorLogger());
//-----------------------------------------------------------------------------
//  it is just not possible to have filters completely independent w/o 
//  huge overhead: hence ways for communication
//  remember, that TStnEvent owns everything in its list of objects
//-----------------------------------------------------------------------------
  fListOfMetVertices = new TClonesArray("TStnVertex",10);
  fListOfMetVertices->SetName("ListOfMetVertices");

  fEvent->AddObject(fListOfMetVertices);

  fEvent->AddDataBlock("CalDataBlock","TCalDataBlock",node);
  _calDataBlock     = (TCalDataBlock*) node->GetDataBlock();

  fEvent->AddDataBlock("ClcDataBlock","TClcDataBlock",node);
  _clcDataBlock     = (TClcDataBlock*) node->GetDataBlock();

  fEvent->AddDataBlock("CmuDataBlock","TCmuDataBlock",node);
  _cmuDataBlock     = (TCmuDataBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("CmpDataBlock","TCmpDataBlock",node);
  _cmpDataBlock     = (TCmpDataBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("CmxDataBlock","TCmxDataBlock",node);
  _cmxDataBlock     = (TCmxDataBlock*) node->GetDataBlock();

  fEvent->AddDataBlock("JetBlock","TStnJetBlock",node);
  _jetBlock         = (TStnJetBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("MuonBlock","TStnMuonBlock",node);
  _muoBlock         = (TStnMuonBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("MetBlock","TStnMetBlock",node);
  _metBlock         = (TStnMetBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("ElectronBlock","TStnElectronBlock",node);
  _eleBlock         = (TStnElectronBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("TauBlock","TStnTauBlock",node);
  _tauBlock         = (TStnTauBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("TrackBlock","TStnTrackBlock",node);
  _trkBlock         = (TStnTrackBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("TriggerBlock","TStnTriggerBlock",node);
  _trgBlock         = (TStnTriggerBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("GenpBlock","TGenpBlock",node);
  _genpBlock        = (TGenpBlock*) node->GetDataBlock();
  fEvent->AddDataBlock("ObspBlock","TObspBlock",node);
  _obspBlock        = (TObspBlock*) node->GetDataBlock();

  fEvent->AddDataBlock("FwdDetDataBlock","TFwdDetDataBlock",node);
  _fwddetDataBlock  = (TFwdDetDataBlock*) node->GetDataBlock();
//-----------------------------------------------------------------------------
//  list of filters
//-----------------------------------------------------------------------------
  fListOfFilters = new TObjArray(10);
  fListOfFilters->SetName("ListOfFilters");
  TModule::fFolder->Add(fListOfFilters);
//-----------------------------------------------------------------------------
//  CMP cuts
//-----------------------------------------------------------------------------
  _cmpNCspHitsMin.addDescription(" \t min N(CSP hits), default = 0");
  commands()->append(&_cmpNCspHitsMin);

//-----------------------------------------------------------------------------
//  L1 cuts
//-----------------------------------------------------------------------------
  _L1Bits.addDescription(" \t fired L1 bits, default = none");
  commands()->append(&_L1Bits);

  _MyronBucket.addDescription(
"         if > 0: Myron bucket to accept, default = -1 (all accepted)");
  commands()->append(&_MyronBucket);

//-----------------------------------------------------------------------------
//  L3 cuts
//-----------------------------------------------------------------------------
  _L3Paths.addDescription(" \t passed L3 paths, default = none");
  commands()->append(&_L3Paths);
//-----------------------------------------------------------------------------
//  GENP and OBSP cuts
//-----------------------------------------------------------------------------
  AddNewCommand(new StntupleFilterModule::Command("useCalFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useClcFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useCmuFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useCmpFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useCmxFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useEleFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useJetFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useL1Filter",this));
  AddNewCommand(new StntupleFilterModule::Command("useL3Filter",this));
  AddNewCommand(new StntupleFilterModule::Command("useMetFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useMuonFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useZmumuFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useMonojetFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useTauFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useTrackFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useGenpFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useObspFilter",this));
  AddNewCommand(new StntupleFilterModule::Command("useFwdDetFilter",this));

  AddNewCommand(new StntupleFilterModule::Command("genpParticle",this));
  AddNewCommand(new StntupleFilterModule::Command("obspParticle",this));
}

//_____________________________________________________________________________
StntupleFilterModule::~StntupleFilterModule() {
  fListOfFilters->Delete();
  delete fListOfFilters;
  delete fEvent;
}

//_____________________________________________________________________________
AppResult StntupleFilterModule::beginJob(AbsEvent* aJob) {

					// define names of the collections used
					// for filtering

  char proc[200], desc[200];

  StntupleGetProcessName(fJetCollName.value().data(),proc,desc);
  if (proc[0] == 0) strcpy(proc,fProcessName.value().data());
  _jetBlock->SetCollName(proc,desc);

  StntupleGetProcessName(fTrackCollName.value().data(),proc,desc);
  if (proc[0] == 0) strcpy(proc,fProcessName.value().data());
  _trkBlock->SetCollName(proc,desc);

  StntupleGetProcessName(fMuonCollName.value().data(),proc,desc);
  if (proc[0] == 0) strcpy(proc,fProcessName.value().data());
  _muoBlock->SetCollName(proc,desc);

  _eleBlock->SetCollName(fProcessName.value().data(),"CalDriven");
  _tauBlock->SetCollName(fProcessName.value().data(),"CdfTauCollection");

				// loop over the filters and pass locations
				// of the data blocks to them

  int nfilters = NFilters();
  for (int i=0; i<nfilters; i++) {
    StntupleFilter* filter = Filter(i);
    filter->Init(fEvent);
  }
//-----------------------------------------------------------------------------
//  set beam history code to use
//-----------------------------------------------------------------------------
  StntupleBeamManager* bm = StntupleBeamManager::Instance();
  bm->SetHistoryCode(fBeamHistoryCode.value());

  return AppResult::OK;
}

//_____________________________________________________________________________
AppResult StntupleFilterModule::beginRun( AbsEvent* aRun ) {
  // read calibrations

//   if (fClcFilter.Used()) {
//     _clcDataBlock->ReadCalibrations(gblEnv->runNumber());
//   }

  AbsEnv* env  = AbsEnv::instance();
  int  run_num = env->runNumber();
  int  mc_flag = env->monteFlag();

  StntupleBeamManager* bm = StntupleBeamManager::Instance();
  bm->BeginRun(run_num,mc_flag);

  TStntuple::Init(run_num);

  return AppResult::OK;
}

//_____________________________________________________________________________
AppResult StntupleFilterModule::event( AbsEvent* event ) {

  bool pass = true;

  int and_mode;

  TStnErrorLogger* logger = fEvent->GetErrorLogger();
  logger->Connect("Report(Int_t, const char*)",
		  "StntupleModule",
		  this,
		  "LogError(const char*)");

  if (strcmp(_filteringMode.value().data(),"and") == 0) and_mode = 1;
  else                                                  and_mode = 0;

//-----------------------------------------------------------------------------
//  clear communications
//-----------------------------------------------------------------------------
  fListOfMetVertices->Clear();
//-----------------------------------------------------------------------------
//  list includes only used filters, not all of them
//-----------------------------------------------------------------------------
  int nfilters = NFilters();
  for (int i=0; i<nfilters; i++) {
    StntupleFilter* filter = Filter(i);
    filter->ProcessEvent(event);
    if (filter->Passed()) {
				// filter passed, in OR mode it is enough
				// to mark event as passed
      if (! and_mode) {
	pass = true;
	goto SET_PASSED;
      }
    }
    else {
				// filter not passed, in AND more this
				// means that event is rejected
      if (and_mode) {
	pass = false;
	goto SET_PASSED;
      }
    }
  }
  
  if (! and_mode) pass = false;

SET_PASSED: 
  if (fDebug.value() > 0) {
    printf(" <%s> : run = %7i event = %7i passed = %i\n",
	   name(), 
	   AbsEnv::instance()->runNumber(), 
	   AbsEnv::instance()->trigNumber(), 
	   pass);
  }
  this->setPassed(pass);
  
//-----------------------------------------------------------------------------
// disconnect from the error reporting signal and return back to AC++
//-----------------------------------------------------------------------------
  logger->Disconnect("Report(Int_t,const char*)",
		     this,"LogError(Int_t,const char*)");

  return AppResult::OK;
}

//_____________________________________________________________________________
AppResult StntupleFilterModule::endRun( AbsEvent* aRun ) {
   return AppResult::OK;
}

//_____________________________________________________________________________
AppResult StntupleFilterModule::endJob( AbsEvent* aJob ) {

  printf("------------------------------------------------------------\n");
  printf("        %s ran in  <%s> mode\n",name(),
	 _filteringMode.value().c_str());

  int nfilters = NFilters();

  double passed_fraction;

  printf("    filter         N(input)  N(passed)  fraction \n");
  printf("------------------------------------------------------------\n");

  for (int i=0; i<nfilters; i++) {
    StntupleFilter* filter = Filter(i);
    if (filter->NInput() != 0) {
      passed_fraction = filter->NPassed()/(filter->NInput()+0.);
    }
    else {
      passed_fraction = 0;
    }

    printf(" %-20s %2i %10i %6.3f\n",
	   filter->GetName(),
	   filter->NInput(),
	   filter->NPassed(),
	   passed_fraction);
    filter->Print();
    printf("------------------------------------------------------------\n");
  }
  return AppResult::OK;
}


//_____________________________________________________________________________
AppResult StntupleFilterModule::abortJob( AbsEvent* aJob ){
  return AppResult::OK;
}


//_____________________________________________________________________________
AppModule* StntupleFilterModule::clone(const char* cloneName) {
  return new StntupleFilterModule(cloneName,
			    "this module is a clone of StntupleFilterModule");
}


//_____________________________________________________________________________
bool StntupleFilterModule::CalFilter(AbsEvent* event){
  StntupleInitCalDataBlock(_calDataBlock,event,0);
  return true;
}

//_____________________________________________________________________________
bool StntupleFilterModule::ClcFilter(AbsEvent* event){

					// unpack CLC banks

  StntupleInitClcDataBlock(_clcDataBlock,event,0);

					// calculate timing
  float time[2];
  
  for (int im=0; im<2; im++) {
    TClcModule* m = _clcDataBlock->Module(im);

					// calculate mean time for the module
					// (and also fo the layers)

    m->CalculateMeanTime(_clcTMin.value(),
			 _clcTMax.value(),
			 _clcAdcThreshold.value());

    time[im] = m->MeanTime();
  }
					// calculate vertex position in cm;
  float z;

  if ( (time[0] > _clcTMin.value()) && (time[0] <= _clcTMax.value()) ) {
    if ( (time[1] > _clcTMin.value()) && (time[1] <= _clcTMax.value()) ) {
      z = (time[0]-time[1])/2.*29.997;
    }
    else {
      z = -700.;
    }
  }
  else {
    if ((time[1] > _clcTMin.value()) && (time[1] <= _clcTMax.value()) ) {
      z = -800;
    }
    else {
					// both times are not defined
      z = -900;
    }
  }
					// now it is time to make sure that 
					// the mean times make sense
  if (fabs(z) < _clcZ0Cut.value()) {
    return true;
  }
  else { 
    return false;
  }
}


//_____________________________________________________________________________
bool StntupleFilterModule::CmuFilter(AbsEvent* event) {
  // to be use  CMU filter one needs to do track reconstruction

  int ntracks;

  StntupleInitCmuDataBlock(_cmuDataBlock,event,0);
  StntupleInitTrackBlock  (_trkBlock    ,event,0);
  ntracks = _trkBlock->NTracks();

  bool rc = false;
  if (_cmuDataBlock->NDHits() >= _cmuNDHitsMin.value()) {

				// quick check if total number of hits in CMU
				// is above the threshold

    int nhits = _cmuDataBlock->NEHits();
    if ( nhits >= _cmuNEHitsMin.value()) {

				// number of hits is large enough, check if
				// there is enough hits in one wedge and also
				// check if there is a high-Pt track pointing
				// to this wedge

      int nhits_per_wedge[2][24];
      int nhits_per_layer[2][24][4];
      int n_hit_layers   [2][24];

      memset(nhits_per_wedge,0,2*24*sizeof(int));
      memset(nhits_per_layer,0,2*24*4*sizeof(int));
      memset(n_hit_layers   ,0,2*24*sizeof(int));

      for (int i=0; i<nhits; i++) {
	TCmueHit* hit = _cmuDataBlock->EHit(i);
	int kb, kw, kl;

	kb = hit->SideNumber();
	kw = hit->WedgeNumber();
	kl = hit->LayerNumber();
	nhits_per_layer[kb][kw][kl]++;
	nhits_per_wedge[kb][kw]++;
	if (nhits_per_layer[kb][kw][kl] == 1) {
	  n_hit_layers[kb][kw]++;
	}
      }

      for (int ib=0; ib<2; ib++) {
	for (int iw=0; iw<24; iw++) {
	  if (n_hit_layers[ib][iw] >= _cmuNEHitsMin.value()) {
	    if (_cmuRequireTrack.value() != 0) {

	      double phi_wedge = (7.5+iw*15.)*M_PI/180.;

				// number of hit layers is OK, check for the
				// presence of high-Pt track

	      for (int i=0; i<ntracks; i++) {
		TStnTrack* trk = _trkBlock->Track(i);
		double pt  = trk->Momentum()->Pt();
		if (pt > _cmuTrackPtMin.value()) {
		  double phi_trk = trk->Momentum()->Phi();
		  double dphi = TVector2::Phi_mpi_pi(phi_trk-phi_wedge);
		  if (fabs(dphi) < _cmuDeltaPhiMax.value()) {
		    rc = true;
		    goto END;
		  }
		}
	      }
	    }
	    else {
				// tracks are not required
	      rc = true;
	      goto END;
	    }
	  }
	}
      }
    }
  }
 END: return rc;
}

//_____________________________________________________________________________
bool StntupleFilterModule::CmpFilter(AbsEvent* event) {
  // for the moment cuts only on the number of CSP hits - there is a 
  // noise in CMP, so cutting on the number of hits in CMPD doesnt' make
  // much sense

  StntupleInitCmpDataBlock(_cmpDataBlock,event,0);

  bool rc = false;
  if (_cmpDataBlock->NCspHits() >= _cmpNCspHitsMin.value()) {
    rc = true;
  }
  return rc;
}

//_____________________________________________________________________________
bool StntupleFilterModule::CmxFilter(AbsEvent* event){
  StntupleInitCmxDataBlock(_cmxDataBlock,event,0);
  return true;
}

//_____________________________________________________________________________
bool StntupleFilterModule::JetFilter(AbsEvent* event){
  // returns true, if the number of jets with 
  // Et >= jetEtMin and jetEtaMin <= eta < jetEtaMax
  // is >= jetNMin

  StntupleInitJetBlock(_jetBlock,event,0);

  int njets = _jetBlock->NJets();

  int nj    = 0;
  for (int i=0; i<njets; i++) {
    TStnJet* jet = _jetBlock->Jet(i);
    double et  = jet->Et();
    double eta = jet->Momentum()->Eta();
    if ( (et >= _jetEtMin.value()) && 
	 (eta >= _jetEtaMin.value()) && (eta < _jetEtaMax.value())) {
      nj++;
    }
  }
  return (nj >= _jetNMin.value());
}

//_____________________________________________________________________________
bool StntupleFilterModule::MonojetFilter(AbsEvent* event) {
  // returns true, if there is 1 jet with Et >= monojetEt1
  // and monojetEtaMin <= eta < monojetEtaMax
  // and no other jet with Et > monojetEt2 within metJetDPhi

  StntupleInitJetBlock  (_jetBlock,event,0);
  StntupleInitTrackBlock(_trkBlock,event,0);
  StntupleInitMetBlock  (_metBlock,event,0);

  int njets = _jetBlock->NJets();

  TStnJet*  monojet;

  int nother = 0;
  int nj     = 0;

  for (int i=0; i<njets; i++) {
    TStnJet* jet = _jetBlock->Jet(i);
    double et  = jet->Et();
    double eta = jet->Momentum()->Eta();
    if ( (et >= _monojetEt1.value()) && 
	 (eta >= _monojetEtaMin.value()) && (eta < _monojetEtaMax.value())) {
      nj++;
      monojet = jet;
    }
  }
  
  bool passed = false;

  if (nj == 1) {
    for (int i=0; i<njets; i++) {
					// OK, only 1 jet of interest
      TStnJet* jet = _jetBlock->Jet(i);
      if (jet != monojet) {
	double et      = jet->Et();
	double jet_phi = jet->Momentum()->Phi();
	if (et >= _monojetEt2.value()) {
	  double dphi = TVector2::Phi_mpi_pi(_metBlock->MetPhi(0)-jet_phi);
	  if (fabs(dphi) < _monojetDPhi.value()) {
	    nother++;
	  }
	}
      }
    }
    if (nother == 0) passed = true;
  }

  if (passed) {
//-----------------------------------------------------------------------------
//  have one jet, check if there are tracks around - this is a very efficient
//  and cheap way to get rid of cosmics - but this also kills photons....
//-----------------------------------------------------------------------------

    int nt    = 0;
    int ntrk  = _trkBlock->NTracks();
    double pt, dphi, deta, dr;

    TLorentzVector* jet_mom = monojet->Momentum();
    double mj_phi= jet_mom->Phi();
    double mj_eta= jet_mom->Eta();

    for (int i=0; i<ntrk; i++) {
      TStnTrack* ti = _trkBlock->Track(i);
      TLorentzVector* mi = ti->Momentum();
      pt  = mi->Pt();
      if (pt > _monojetTrackPt.value()) {

				// high-Pt track, check delta(eta-phi)

	dphi = TVector2::Phi_mpi_pi(mj_phi-ti->Momentum()->Phi());
	deta = mj_eta-ti->Momentum()->Eta();
	dr   = sqrt(dphi*dphi+deta*deta);
	if (dr < _monojetCone.value()) {
	  nt++;
	}
      }
    }
    passed = (nt >= _monojetNTracks.value());
  }

  return passed;
}

//_____________________________________________________________________________
bool StntupleFilterModule::MetFilter(AbsEvent* Event) {
  // returns true, if MET > metMin

  static PhysicsTowerDataMaker ptd_maker;
  static NaiveEtCalculator     met_calculator;
  static TStnMuonID            muonID;

  static int first_call(1);

  if (first_call) {
    first_call = 0;
    muonID.SetMinTrackBcPt   ( 0.);
    muonID.SetMinTrackPt     (15.);
    muonID.SetMinLooseTrackPt(10.);
  }

  CdfMet_h                     met_h;
  PhysicsTowerData_ch          towers_ch;

  double                       met;
  TVector2                     met_v2, corr_met_v2;

  TStnVertex*                  v;
  bool                         status;
  bool                         passed = false;

  StntupleInitMetBlock  (_metBlock  ,Event,0);

  if (_metMode.value() == 0) {
//-----------------------------------------------------------------------------
// cut on non-corrected MET, calculated at Z=0
//-----------------------------------------------------------------------------
    met    = _metBlock->Met(0);
  }
  else if (_metMode.value() == 1) {
//-----------------------------------------------------------------------------
// cut on MET calculated in the highest sumPt Z-vertex corrected for muons
// tight muons only...
//-----------------------------------------------------------------------------
    StntupleInitTrackBlock(_trkBlock,Event,0);
    StntupleInitMuonBlock (_muoBlock,Event,0);
    StntupleMuonBlockLinks(_muoBlock,Event,0);

    met_v2.Set(_metBlock->MetX(3),_metBlock->MetY(3));

    TStntuple::CorrectMetForMuons(_muoBlock,&muonID,&met_v2,&corr_met_v2,0);
    met = corr_met_v2.Mod();
  }
  else if (_metMode.value() == 2) {
//-----------------------------------------------------------------------------
// loop over the list of vertices in fListOfMetVertices, calculate MET for
// each of the vertices in that list and see if for any of them the calculated
// value of MET is above the threshold
//-----------------------------------------------------------------------------
    int vtx_strategy = 1;
    NaiveEtParams met_parms(0.1,vtx_strategy);
    int nv = fListOfMetVertices->GetEntriesFast();
    for (int iv=0; iv<nv; iv++) {
      v     = (TStnVertex*) fListOfMetVertices->UncheckedAt(iv);
      met_h = new CdfMet(v->Z(),vtx_strategy);
//-----------------------------------------------------------------------------
// second parameter =2: MetCalculator
//-----------------------------------------------------------------------------
      PhysicsTowerParams pdt_parms(0,2,0.1,v->Z());
      towers_ch = ptd_maker.create(Event,pdt_parms,"StntupleFilterModule");
      status    = met_calculator.calculateMet(met_h,towers_ch);
      met       = met_h->met();
      if (! status) {
	ERRLOG(ELerror,"StntupleMetFilter: ") << "failed to fill CdfMet" 
					      << endmsg;
      }
    }
  }
//-----------------------------------------------------------------------------
// clean up the memory and return
//-----------------------------------------------------------------------------
  passed = (met >= _metMin.value());
  return passed;
}

//_____________________________________________________________________________
bool StntupleFilterModule::TauFilter(AbsEvent* event){
  // returns true, if the number of tau candidates with 
  // Pt >= tauPtMin is >= tauNMin

  StntupleInitTauBlock(_tauBlock,event,0);
  int ntaus = _tauBlock->NTaus();

  int nt = 0;
  for (int i=0; i<ntaus; i++) {
    TStnTau* tau = _tauBlock->Tau(i);
    double pt   = tau->Momentum()->Pt();
    double emfr = tau->Emfr();
    if ((pt   >= _tauPtMin.value  ()) && 
	(pt   <  _tauPtMax.value  ()) &&
	(emfr >= _tauEmfrMin.value()) &&
	(emfr <  _tauEmfrMax.value())    ) {
      nt++;
    }
  }
  return (nt >= _tauNMin.value());
}

//_____________________________________________________________________________
bool StntupleFilterModule::GenpFilter(AbsEvent* event) {
  //
  int nfound [100]; 

  memset(nfound,0,4*sizeof(int));

  StntupleInitGenpBlock(_genpBlock,event,0);
  int npart = _genpBlock->NParticles();

  //  int np = 0;
  Double_t eta;

  int nfound_tot = 0;

  int ncuts = _genpCuts.size();

  for (int i=0; i<npart; i++) {
    TGenParticle* p = _genpBlock->Particle(i);
    for (int icut=0; icut<ncuts; icut++) {
      ParticleCuts_t& genp = _genpCuts[icut];
      if (p->GetPdgCode() == genp.fPdgCode) {
	double pt  = p->Pt();
	if (pt >= genp.fPtMin) {
	  eta = p->Eta();
	  if ((eta >= genp.fEtaMin) && (eta < genp.fEtaMax)) {
	    nfound[icut]++;
	    nfound_tot++;
	  }
	}
      }
    }
  }

  int passed = 1;

  if (_genpNMin.value() == 0) {
				// check if each cut is satisfied

    for (int icut=0; icut<ncuts; icut++) {
      ParticleCuts_t& genp = _genpCuts[icut];
      if (nfound[icut] < genp.fN) {
	passed = 0;
	break;
      }
    }
  }
  else {
				// check if total number of particles
				// passing the cuts is above the threshold

    if (nfound_tot < _genpNMin.value()) {
      passed = 0;
    }
  }

  if (passed && verbose()) {
    _genpBlock->Print();
  }
  return passed;
}

//_____________________________________________________________________________
bool StntupleFilterModule::ObspFilter(AbsEvent* event) {
  //
  int nfound [100]; 

  memset(nfound,0,4*sizeof(int));

  StntupleInitObspBlock(_obspBlock,event,0);
  int npart = _obspBlock->NParticles();

  int nfound_tot = 0;

  //  int np = 0;
  Double_t eta;

  int ncuts = _obspCuts.size();

  for (int i=0; i<npart; i++) {
    TObspParticle* p = _obspBlock->Particle(i);
    for (int icut=0; icut<ncuts; icut++) {
      ParticleCuts_t& obsp = _obspCuts[icut];
      if (p->PdgCode() == obsp.fPdgCode) {
	double pt  = p->Momentum()->Pt();
	if (pt >= obsp.fPtMin) {
	  eta = p->Momentum()->Eta();
	  if ((eta >= obsp.fEtaMin) && (eta < obsp.fEtaMax)) {
	    nfound[icut]++;
	    nfound_tot++;
	  }
	}
      }
    }
  }

  int passed = 1;

  if (_obspNMin.value() == 0) {
				// see if each of cuts is satisfied

    for (int icut=0; icut<ncuts; icut++) {
      ParticleCuts_t& obsp = _obspCuts[icut];
      if (nfound[icut] < obsp.fN) {
	passed = 0;
	break;
      }
    }
  }
  else {
				// check if total number of particles
				// passing the cuts is above the threshold

    if (nfound_tot < _obspNMin.value()) {
      passed = 0;
    }
  }

  if (passed && verbose()) {
    _obspBlock->Print();
  }
  return passed;
}

//_____________________________________________________________________________
bool StntupleFilterModule::L1Filter(AbsEvent* event){
  // returns true, if at leas one from the L1 trigger bits defined by the user
  // is fired

  bool passed = true;

  StntupleInitTriggerBlock(_trgBlock,event,0);

  TTfrd* fred = _trgBlock->Tfrd();

  int nbits = _L1Bits.size();
  if (nbits > 0) {
				// need to check L1 bits
    passed = false;
    for (AbsParmList<int>::ConstIterator bit = _L1Bits.begin();
	 bit != _L1Bits.end(); ++bit) {
      if (fred->PrescaledL1Bit(*bit) == 1) {
	passed = true;
	break;
      }
    }
    if (! passed) goto EXIT;
  }
				// check number of Myron bucket

  if (_MyronBucket.value() >= 0) {
    TTsid* tsid = _trgBlock->Tsid();
    if (tsid->MyronFlag() != _MyronBucket.value()) {
      passed = false;
    }
    if (! passed) goto EXIT;
  }


EXIT:  return passed;
}


//_____________________________________________________________________________
bool StntupleFilterModule::L3Filter(AbsEvent* event){
  // returns true, if at least one from the requested L3 trigger 
  // paths is passed

  bool passed = true;

  StntupleInitTriggerBlock(_trgBlock,event,0);

  TTl3d* tl3d = _trgBlock->Tl3d();

  int npaths = _L3Paths.size();
  if (npaths > 0) {
				// need to check L3 paths
    passed = false;
    for (AbsParmList<int>::ConstIterator bit = _L3Paths.begin();
	 bit != _L3Paths.end(); ++bit) {
      if ((*bit                   >= 0            ) && 
	  (*bit                   < tl3d->NPaths()) && 
	  (tl3d->PathPassed(*bit) == 1            )   ) {
	passed = true;
	break;
      }
    }
    if (! passed) goto EXIT;
  }

EXIT:  return passed;
}

//_____________________________________________________________________________
bool StntupleFilterModule::FwdDetFilter(AbsEvent* event){
  StntupleInitFwdDetDataBlock(_fwddetDataBlock,event,0);
  return true;
}

//_____________________________________________________________________________
StntupleFilterModule::Command::Command() {
}

//_____________________________________________________________________________
StntupleFilterModule::Command::~Command() {
}


//_____________________________________________________________________________
StntupleFilterModule::Command::Command(const char* const name, AppModule* module) : 
    APPCommand(name,module)
{
}


//_____________________________________________________________________________
int StntupleFilterModule::Command::handle(int argc, char* argv[]) {
  //

  int rc        = 0;
  ParticleCuts_t part;

  StntupleFilterModule* module = (StntupleFilterModule*) target();

  TObjArray* list = module->GetListOfFilters();

  if (strncmp(command(),"genpParticle",12)==0) {
    if (argc == 6) {
      part.fPdgCode  = atoi(argv[1]);
      part.fN        = atoi(argv[2]);
      part.fPtMin    = atof(argv[3]);
      part.fEtaMin   = atof(argv[4]);
      part.fEtaMax   = atof(argv[5]);
      module->_genpCuts.push_back(part);
    }
    else {
      rc = -1;
    }
  }
  else if (strncmp(command(),"obspParticle",12)==0) {
    if (argc == 6) {
      part.fPdgCode  = atoi(argv[1]);
      part.fN        = atoi(argv[2]);
      part.fPtMin    = atof(argv[3]);
      part.fEtaMin   = atof(argv[4]);
      part.fEtaMax   = atof(argv[5]);
      module->_obspCuts.push_back(part);
    }
    else {
      rc = -1;
    }
  }
  else if (strcmp(command(),"useCalFilter") == 0) {
    list->Add(new StntupleFilter("CalFilter",module,StntupleCalFilter));
  }
  else if (strcmp(command(),"useClcFilter") == 0) {
    list->Add(new StntupleFilter("ClcFilter",module,StntupleClcFilter));
  }
  else if (strcmp(command(),"useCmuFilter") == 0) {
    list->Add(new StntupleFilter("CmuFilter",module,StntupleCmuFilter));
  }
  else if (strcmp(command(),"useCmpFilter") == 0) {
    list->Add(new StntupleFilter("CmpFilter",module,StntupleCmpFilter));
  }
  else if (strcmp(command(),"useCmxFilter") == 0) {
    list->Add(new StntupleFilter("CmxFilter",module,StntupleCmxFilter));
  }
  else if (strcmp(command(),"useEleFilter") == 0) {
    list->Add(new StntupleEleFilter("EleFilter",module));
  }
  else if (strcmp(command(),"useJetFilter") == 0) {
    list->Add(new StntupleFilter("JetFilter",module,StntupleJetFilter));
  }
  else if (strcmp(command(),"useMetFilter") == 0) {
    list->Add(new StntupleFilter("MetFilter",module,StntupleMetFilter));
  }
  else if (strcmp(command(),"useMuonFilter") == 0) {
    list->Add(new StntupleMuonFilter("MuonFilter",module));
  }
  else if (strcmp(command(),"useZmumuFilter") == 0) {
    list->Add(new StntupleZmumuFilter("ZmumuFilter",module));
  }
  else if (strcmp(command(),"useMonojetFilter") == 0) {
    list->Add(new StntupleFilter("MonojetFilter",module,StntupleMonojetFilter));
  }
  else if (strcmp(command(),"useTauFilter") == 0) {
    list->Add(new StntupleFilter("TauFilter",module,StntupleTauFilter));
  }
  else if (strcmp(command(),"useTrackFilter") == 0) {
    list->Add(new StntupleTrackFilter("TrackFilter",module));
  }
  else if (strcmp(command(),"useL1Filter") == 0) {
    list->Add(new StntupleFilter("L1Filter",module,StntupleL1Filter));
  }
  else if (strcmp(command(),"useL3Filter") == 0) {
    list->Add(new StntupleFilter("L3Filter",module,StntupleL3Filter));
  }
  else if (strcmp(command(),"useGenpFilter") == 0) {
    list->Add(new StntupleFilter("GenpFilter",module,StntupleGenpFilter));
  }
  else if (strcmp(command(),"useObspFilter") == 0) {
    list->Add(new StntupleFilter("ObspFilter",module,StntupleObspFilter));
  }
  else if (strcmp(command(),"useFwdDetFilter") == 0) {
    list->Add(new StntupleFilter("FwdDetFilter",module,StntupleFwdDetFilter));
  }
//-----------------------------------------------------------------------------
//  electron commands
//-----------------------------------------------------------------------------
  else if (strncmp(command(),"ele",3) == 0) {
    StntupleFilter* filter = module->Filter("EleFilter");
    if (filter) {
      return filter->Command(command(),argc,argv);
    }
    else {
      return -1;
    }
  }
//-----------------------------------------------------------------------------
//  muon commands
//-----------------------------------------------------------------------------
  else if (strncmp(command(),"muo",3) == 0) {
    StntupleFilter* filter = module->Filter("MuonFilter");
    if (filter) {
      return filter->Command(command(),argc,argv);
    }
    else {
      return -1;
    }
  }
//-----------------------------------------------------------------------------
//  track commands
//-----------------------------------------------------------------------------
  else if (strncmp(command(),"trk",3) == 0) {
    StntupleFilter* filter = module->Filter("TrackFilter");
    if (filter) {
      return filter->Command(command(),argc,argv);
    }
    else {
      return -1;
    }
  }

  return rc;
}

//_____________________________________________________________________________
void StntupleFilterModule::Command::show() const {
  if(strncmp(command(),"DETECTORS",9)==0) {
  }
}

//_____________________________________________________________________________
bool StntupleFilterModule::Command::isShowable() const {
  return true;
}

//_____________________________________________________________________________
string StntupleFilterModule::Command::description() const {
  string retval;
  if (strcmp(command(),"useEleFilter") == 0) {
    retval = "Use Electron Filter and read electron cuts menu";
  }
  if (strcmp(command(),"useMuonFilter") == 0) {
    retval = "Use Muon Filter";
  }
  if (strcmp(command(),"useTrackFilter") == 0) {
    retval = "Use Track Filter";
  }
  return retval;
}


//_____________________________________________________________________________
APPCommand*  StntupleFilterModule::NewCommand(const char* Name) {
  return new StntupleFilterModule::Command(Name,this);
}