///////////////////////////////////////////////////////////////////////////////
//
// Summary of fPrintLevel values:
// ----------------------
// fPrintLevel =   1: print Muon block for all the events
// fPrintLevel = 101: print line per event with good muon
// fPrintLevel = 102: found matching non-trigger jet
///////////////////////////////////////////////////////////////////////////////
#include "TF1.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TText.h"

#include "Stntuple/obj/TStnHeaderBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/alg/TStntuple.hh"

#include "Stntuple/loop/TStnAna.hh"

#include "Stntuple/geom/TTrajectoryPoint.hh"
#include "Stntuple/geom/TSimpleExtrapolator.hh"

#include "TMuoAnaModule.hh"
//#include "TTauUtils.hh"

ClassImp(TMuoAnaModule)
//_____________________________________________________________________________
TMuoAnaModule::TMuoAnaModule(const char* name, const char* title): 
  TStnModule(name,title)
{

  //  int unused;

  fMinNGoodMuons     = 1;
  fMinDrMetJet       = 0.5;
  
  fListOfGoodMuons = new TObjArray(20);
//-----------------------------------------------------------------------------
//  ID bits
//-----------------------------------------------------------------------------
  fDetectorBit       = (1 <<  0);
  fPtBit             = (1 <<  1);
  fCmuDelXBit        = (1 <<  2);
  fCmuDelZBit        = (1 <<  3);
  fCmpDelXBit        = (1 <<  4);
  fCmxDelXBit        = (1 <<  5);
  fTrackZ0Bit        = (1 <<  6);
  fTrackD0Bit        = (1 <<  7);
  fTIsoBit           = (1 <<  8);
  fEIsoBit           = (1 <<  9);
  fEHadBit           = (1 << 10);
  fEEmBit            = (1 << 11);
  fCotNAxBit         = (1 << 12);
  fCotNStBit         = (1 << 13);
  fSeedIEtaBit       = (1 << 14);
  fWedgeNumberBit    = (1 << 15);
  //  unused             = (1 << 16);
  fIso1Bit           = (1 << 17);

//-----------------------------------------------------------------------------
// want al least CMU or CMP stub
//-----------------------------------------------------------------------------
  fDetectorMask      = 0x3;
  fMaxEEm            =  2.;
  fMaxEHad           =  6.;
  fMinGoodMuonPt     = 20.;
  fMaxCmuDelX        =  3.;
  fMaxCmuDelZ        = 1.e10;
  fMaxCmpDelX        =  6.;
  fMaxCmxDelX        = 100.;
  fMaxEIso           = 2.;
  fMaxTIso           = 1000.;
  fMaxZ0             = 60;
  fMaxD0             = 1.0;
  fMinNCotAxHits     = 25;
  fMinNCotStHits     = 25;
  fExtrapolator      = new TSimpleExtrapolator();

//-----------------------------------------------------------------------------
// modes:
//-----------------------------------------------------------------------------
  fDoQCDFakes        = 0;
  fCutMode           = 0; 
//-----------------------------------------------------------------------------
// fudging [0]: data; [1]: MC
//-----------------------------------------------------------------------------
  fFudgeScaleFactor[0] =  0.012;
  fFudgeScaleFactor[1] = -0.006;
  fFudgeSigma      [0] = 0.;
  fFudgeSigma      [1] = 0.02;
  fDoFudging           = 0;
}

//_____________________________________________________________________________
TMuoAnaModule::~TMuoAnaModule() {

  fListOfGoodMuons->Clear();
  delete fListOfGoodMuons;

  fListOfNonTrigJets->Clear();
  delete fListOfNonTrigJets;
  
  fListOfTriggerJets->Clear();
  delete fListOfTriggerJets;
}


//_____________________________________________________________________________
int TMuoAnaModule::BeginJob() {
					// register the data blocks

  RegisterDataBlock("MuonBlock"    , "TStnMuonBlock"   , &fMuonBlock);
  RegisterDataBlock("TauBlock"     , "TStnTauBlock"    , &fTauBlock);
  RegisterDataBlock("TriggerBlock" , "TStnTriggerBlock", &fTriggerBlock);
  RegisterDataBlock("TrackBlock"   , "TStnTrackBlock"  , &fTrackBlock);
  RegisterDataBlock("MetBlock"     , "TStnMetBlock"    , &fMetBlock);
  RegisterDataBlock("JetBlock"     , "TStnJetBlock"    , &fJetBlock);

					// book histograms
  BookHistograms();

  return 0;
}


//_____________________________________________________________________________
int TMuoAnaModule::BeginRun() {

  int rn = GetHeaderBlock()->RunNumber();
  TStntuple::Init(rn);

  return 0;
}


//_____________________________________________________________________________
void TMuoAnaModule::FillMuonHistograms(MuoHist_t& Hist, TStnMuon* Muo) {

  const TLorentzVector* mom = Muo->Momentum();
  double  pt = mom->Pt();

  Hist.fPt->Fill(pt);
  Hist.fPhi->Fill(TVector2::Phi_0_2pi(mom->Phi()));
  Hist.fEIso->Fill(Muo->Iso());
  Hist.fHadEnergy->Fill(Muo->HadEnergy());
  Hist.fEmEnergy->Fill(Muo->EmEnergy());

  Hist.fDetectorVsEta->Fill(Muo->Momentum()->Eta(),Muo->Detector(),1);

//-----------------------------------------------------------------------------
//  fill track histograms
//-----------------------------------------------------------------------------
  //  TLorentzVector* tmom;
  int it = Muo->TrackNumber(); 
  if (it >= 0) {
    TStnTrack* trk = fTrackBlock->Track(it);
    //    tmom  = trk->Momentum();
    Hist.fTIso->Fill(trk->Iso4());
    Hist.fTrackD0->Fill(trk->D0());
    Hist.fTrackD0Corr->Fill(TStntuple::CorrectedD0(trk));
    Hist.fTrackZ0->Fill(trk->Z0());
    Hist.fTrackChi2->Fill(trk->Chi2()/(trk->NCotHitsTot()+trk->NSvxHits()-4.9999));
    Hist.fTrackCotChi2->Fill(trk->Chi2Cot()/(trk->NCotHitsTot()-4.9999));
    Hist.fTrackNCotHitsTot->Fill(trk->NCotHitsTot());
    Hist.fTrackNCotHitsAx->Fill(trk->NCotHitsAx());
    Hist.fTrackNCotHitsSt->Fill(trk->NCotHitsSt());
    Hist.fTrackAlgorithm->Fill(trk->Algorithm());
  }
}

//_____________________________________________________________________________
void TMuoAnaModule::FillMuonIDHistograms(MuoIDHist_t&  Hist, TStnMuon* Muo) 
{
//-----------------------------------------------------------------------------
// distributions for ID variables (with n-1 ID cuts applied). 
// Plot 3 sets of histograms:
// --------------------------
// set 1: for all the muons
// set 2: for the parameter given all the rest cuts were successfull
// set 3&4: for the cuts in sequence they are applied
//-----------------------------------------------------------------------------
  double  iso1;
  int     id_word, imuo;
 
  // TStnTrack* trk;

  const TLorentzVector* mom = Muo->Momentum();
  double  pt = mom->Pt();

  imuo    = Muo->Number();
  id_word = fMuoIDWord[imuo];
//-----------------------------------------------------------------------------
//  1. muon detector
//-----------------------------------------------------------------------------
  Hist.fDetector[0]->Fill(Muo->Detector());
  if ((id_word & ~fDetectorBit) == 0) Hist.fDetector[1]->Fill(Muo->Detector());
  if (id_word == 0) Hist.fDetector[4]->Fill(Muo->Detector());

  Hist.fTrackPt[0]->Fill(pt);
  if ((id_word & ~fPtBit) == 0) Hist.fTrackPt[1]->Fill(pt);
  if (id_word == 0) Hist.fTrackPt[4]->Fill(pt);
//-----------------------------------------------------------------------------
//  1. calorimetry isolation
//-----------------------------------------------------------------------------
  Hist.fEIso[0]->Fill(Muo->Iso());
  if ((id_word & ~fEIsoBit) == 0) Hist.fEIso[1]->Fill(Muo->Iso());
  if (id_word == 0) Hist.fEIso[4]->Fill(Muo->Iso());

  iso1 = Muo->Iso()/pt;

  Hist.fIso1[0]->Fill(iso1);
  if ((id_word & ~fIso1Bit) == 0) Hist.fIso1[1]->Fill(iso1);
  if (id_word == 0) Hist.fIso1[4]->Fill(iso1);
//-----------------------------------------------------------------------------
//  2. calorimetry
//-----------------------------------------------------------------------------
  Hist.fHadEnergy[0]->Fill(Muo->HadEnergy());
  if ((id_word & ~fEHadBit) == 0) Hist.fHadEnergy[1]->Fill(Muo->HadEnergy());
  if (id_word == 0) Hist.fHadEnergy[4]->Fill(Muo->HadEnergy());

  Hist.fEmEnergy[0]->Fill(Muo->EmEnergy());
  if ((id_word & ~fEEmBit) == 0) Hist.fEmEnergy[1]->Fill(Muo->EmEnergy());
  if (id_word == 0) Hist.fEmEnergy[4]->Fill(Muo->EmEnergy());
//-----------------------------------------------------------------------------
//  7. CMU residuals
//-----------------------------------------------------------------------------
  if (Muo->HasCmuStub()) {
    Hist.fCmuDelX[0]->Fill(Muo->CmuDelX());
    if ((id_word & ~fCmuDelXBit) == 0) Hist.fCmuDelX[1]->Fill(Muo->CmuDelX());
    if (id_word == 0) Hist.fCmuDelX[4]->Fill(Muo->CmuDelX());

    Hist.fCmuDelZ[0]->Fill(Muo->CmuDelZ());
    if ((id_word & ~fCmuDelZBit) == 0) Hist.fCmuDelZ[1]->Fill(Muo->CmuDelZ());
    if (id_word == 0) Hist.fCmuDelZ[4]->Fill(Muo->CmuDelZ());
  }
//-----------------------------------------------------------------------------
//  8. CMP residuals
//-----------------------------------------------------------------------------
  if (Muo->HasCmpStub()) {
    Hist.fCmpDelX[0]->Fill(Muo->CmpDelX());
    if ((id_word & ~fCmpDelXBit) == 0) Hist.fCmpDelX[1]->Fill(Muo->CmpDelX());
    if (id_word == 0) Hist.fCmpDelX[4]->Fill(Muo->CmpDelX());
  }
//-----------------------------------------------------------------------------
//  9. CMX residuals
//-----------------------------------------------------------------------------
  if (Muo->HasCmxStub()) {
    Hist.fCmxDelX[0]->Fill(Muo->CmxDelX());
    if ((id_word & ~fCmxDelXBit) == 0) Hist.fCmxDelX[1]->Fill(Muo->CmxDelX());
    if (id_word == 0) Hist.fCmxDelX[4]->Fill(Muo->CmxDelX());
  }
//-----------------------------------------------------------------------------
// 10. fill track histograms - cut on corrected D0
//-----------------------------------------------------------------------------
  Hist.fTIso[0]->Fill(Muo->TIso());
  if ((id_word & ~fTIsoBit) == 0) Hist.fTIso[1]->Fill(Muo->TIso());
  if (id_word == 0) Hist.fTIso[4]->Fill(Muo->TIso());

  double d0_corr;
  int    n_ax_hits, n_st_hits;

  int it = Muo->TrackNumber(); 
  if (it < 0) {
//-----------------------------------------------------------------------------
// muon track is not defined
//-----------------------------------------------------------------------------
    d0_corr   = 1.e6;
    n_ax_hits = -1;
    n_st_hits = -1;
  }
  else {
//-----------------------------------------------------------------------------
// muon has a track
//-----------------------------------------------------------------------------
    TStnTrack* trk     = fTrackBlock->Track(it);
    d0_corr = TStntuple::CorrectedD0(trk);
    n_ax_hits = trk->NCotHitsAx();
    n_st_hits = trk->NCotHitsSt();
  }

  Hist.fTrackD0[0]->Fill(d0_corr);
  if ((id_word & ~fTrackD0Bit) == 0) Hist.fTrackD0[1]->Fill(d0_corr);
  if (id_word == 0) Hist.fTrackD0[4]->Fill(d0_corr);

  Hist.fTrackZ0[0]->Fill(Muo->Z0());
  if ((id_word & ~fTrackZ0Bit) == 0) Hist.fTrackZ0[1]->Fill(Muo->Z0());
  if (id_word == 0) Hist.fTrackZ0[4]->Fill(Muo->Z0());

  Hist.fCotNAxHits[0]->Fill(n_ax_hits);
  if ((id_word & ~fCotNAxBit) == 0) Hist.fCotNAxHits[1]->Fill(n_ax_hits);
  if (id_word == 0) Hist.fCotNAxHits[4]->Fill(n_ax_hits);

  Hist.fCotNStHits[0]->Fill(n_st_hits);
  if ((id_word & ~fCotNStBit) == 0) Hist.fCotNStHits[1]->Fill(n_st_hits);
  if (id_word == 0) Hist.fCotNStHits[4]->Fill(n_st_hits);
//-----------------------------------------------------------------------------
//  single histogram showing how often every particular cut failed
//-----------------------------------------------------------------------------
  for (int bit=0; bit<32; bit++) {
    if (((id_word >> bit) & 0x1) == 1) {
      Hist.fFailedBits->Fill(bit);
    }
  }
  Hist.fPassed->Fill(id_word == 0);
//-----------------------------------------------------------------------------
//  ***** now histogram effect of cuts in the order they are applied
//-----------------------------------------------------------------------------
  Hist.fDetector[2]->Fill(Muo->Detector());
  if ((id_word & fDetectorBit) != 0)                     goto END;
  Hist.fDetector[3]->Fill(Muo->Detector());

  Hist.fTrackPt[2]->Fill(pt);
  if ((id_word & fPtBit) != 0)                        goto END;
  Hist.fTrackPt[3]->Fill(pt);

  Hist.fEIso[2]->Fill(Muo->Iso());
  if ((id_word & fEIsoBit) != 0)                     goto END;
  Hist.fEIso[3]->Fill(Muo->Iso());

  Hist.fIso1[2]->Fill(Muo->TIso());
  if ((id_word & fIso1Bit) != 0)                      goto END;
  Hist.fIso1[3]->Fill(Muo->TIso());

  Hist.fHadEnergy[2]->Fill(Muo->HadEnergy());
  if ((id_word & fEHadBit) != 0)                     goto END;
  Hist.fHadEnergy[3]->Fill(Muo->HadEnergy());

  Hist.fEmEnergy[2]->Fill(Muo->EmEnergy());
  if ((id_word & fEEmBit) != 0)                     goto END;
  Hist.fEmEnergy[3]->Fill(Muo->EmEnergy());

  if (Muo->HasCmuStub()) {
    Hist.fCmuDelX[2]->Fill(Muo->CmuDelX());
    if ((id_word & fCmuDelXBit) != 0)                      goto END;
    Hist.fCmuDelX[3]->Fill(Muo->CmuDelX());

    Hist.fCmuDelZ[2]->Fill(Muo->CmuDelZ());
    if ((id_word & fCmuDelZBit) != 0)                      goto END;
    Hist.fCmuDelZ[3]->Fill(Muo->CmuDelZ());
  }

  Hist.fTrackD0[2]->Fill(d0_corr);
  if ((id_word & fTrackD0Bit) != 0)                   goto END;
  Hist.fTrackD0[3]->Fill(d0_corr);

  Hist.fTrackZ0[2]->Fill(Muo->Z0());
  if ((id_word & fTrackZ0Bit) != 0)                   goto END;
  Hist.fTrackZ0[3]->Fill(Muo->Z0());

  Hist.fCotNAxHits[2]->Fill(n_ax_hits);
  if ((id_word & fCotNAxBit) != 0)                    goto END;
  Hist.fCotNAxHits[3]->Fill(n_ax_hits);

  Hist.fCotNStHits[2]->Fill(n_st_hits);
  if ((id_word & fCotNStBit) != 0)                    goto END;
  Hist.fCotNStHits[3]->Fill(n_st_hits);

 END:;
}

//_____________________________________________________________________________
// fill muon ID word: this is a bit-packed word, where each bit tells
// whether an muon passed certain cut
// if muon passes the cut, the corresposing bit is set to 0, otherwise 
// the bit is set to 1
//-----------------------------------------------------------------------------
int TMuoAnaModule::MuoIDWord(const TStnMuon* Muo) 
{
  TStnTrack* trk;
  int        it;
  Int_t      CutWord = 0;

  double d0;
  double pt   = Muo->Momentum()->Pt();
  double iso1 = Muo->TIso()/pt;
//-----------------------------------------------------------------------------
// muon cuts
//-----------------------------------------------------------------------------
  if ((Muo->Detector() & fDetectorMask) == 0     ) CutWord |= fDetectorBit;
  if (pt                   < fMinGoodMuonPt      ) CutWord |= fPtBit;
  if (Muo->EmEnergy()      > fMaxEEm             ) CutWord |= fEEmBit;
  if (Muo->HadEnergy()     > fMaxEHad            ) CutWord |= fEHadBit;

  if (Muo->HasCmuStub()) {
    if (fabs(Muo->CmuDelX()) > fMaxCmuDelX)   CutWord |= fCmuDelXBit;
    if (fabs(Muo->CmuDelZ()) > fMaxCmuDelZ)   CutWord |= fCmuDelZBit;
  }

  if (Muo->HasCmpStub()) {
    if (fabs(Muo->CmpDelX()) > fMaxCmpDelX)   CutWord |= fCmpDelXBit;
  }

  if (Muo->HasCmxStub()) {
    if (fabs(Muo->CmxDelX()) > fMaxCmxDelX)   CutWord |= fCmxDelXBit;
  }

//-----------------------------------------------------------------------------
// isolation
//-----------------------------------------------------------------------------
  if (Muo->Iso()           > fMaxEIso            ) CutWord |= fEIsoBit;

  it = Muo->TrackNumber(); 
  if (it < 0) {
//-----------------------------------------------------------------------------
// muon track is not defined
//-----------------------------------------------------------------------------
    CutWord |= fTIsoBit; 
    CutWord |= fTrackZ0Bit;
    CutWord |= fTrackD0Bit;
    CutWord |= fCotNAxBit;
    CutWord |= fCotNStBit;
  }
  else {
//-----------------------------------------------------------------------------
// muon has a track
//-----------------------------------------------------------------------------
    trk = fTrackBlock->Track(it);
    if (trk->Iso4()       > fMaxTIso)      CutWord |= fTIsoBit;
    if (fabs(trk->Z0())   > fMaxZ0  )      CutWord |= fTrackZ0Bit;

    d0 = TStntuple::CorrectedD0(trk);
    if (TMath::Abs(d0)    > fMaxD0  )      CutWord |= fTrackD0Bit;

    if (trk->NCotHitsAx() < fMinNCotAxHits) CutWord |= fCotNAxBit;
    if (trk->NCotHitsSt() < fMinNCotStHits) CutWord |= fCotNStBit;
  }

  return CutWord;
}
  
//_____________________________________________________________________________
void TMuoAnaModule::FillEventHistograms(EventHist_t&  Hist) {
  Hist.fNMuons->Fill(fNRecoMuons);
  Hist.fNGoodMuons->Fill(fNGoodMuons);

  Hist.fNCmuStubs->Fill(fMuonBlock->NCmuStubs());
  Hist.fNCmpStubs->Fill(fMuonBlock->NCmpStubs());
  Hist.fNCmxStubs->Fill(fMuonBlock->NCmxStubs());

//    Hist.fMet0      [0]->Fill(fMetBlock->Met(0));
//    Hist.fMet1      [0]->Fill(fMetBlock->Met(1));
//    if (fNGoodMuons) {
//      Hist.fMet0      [1]->Fill(fMetBlock->Met(0));
//      Hist.fMet1      [1]->Fill(fMetBlock->Met(1));
//    }

//-----------------------------------------------------------------------------
//  transverse mass histograms
//-----------------------------------------------------------------------------
  double metm[2], metx[2], mety[2], q[2];
  //  double metm[2], metx[2], mety[2], mt[2], q[2];

  for (int i=0; i<2; i++) {
    metm[i] = fMetBlock->Met(i);
    metx[i] = metm[i]*cos(fMetBlock->MetPhi(i));
    mety[i] = metm[i]*sin(fMetBlock->MetPhi(i));
  }

  if (fNGoodMuons > 0) {

    for (int i1=0; i1<fNRecoMuons; i1++) {
      TStnMuon*   e1 = fMuonBlock->Muon(i1);
      const TLorentzVector* m1 = e1->Momentum();
      if (fMuoIDWord[i1] == 0) {
	
//         	mt[0] = sqrt(2*(e1->Et()*metm[0] -m1->Px()*metx[0]-m1->Py()*mety[0]));
//         	mt[1] = sqrt(2*(e1->Et()*metm[1] -m1->Px()*metx[1]-m1->Py()*mety[1]));
	
//  	Hist.fMt[0]->Fill(mt[0]);
//  	Hist.fMt[1]->Fill(mt[1]);

//-----------------------------------------------------------------------------
// distributions for isolation vs MET
//-----------------------------------------------------------------------------
//  	fHist.fIsoVsMet->Fill(fMetBlock->Met(1),e1->Iso());
//  	fHist.fIso1VsMet->Fill(fMetBlock->Met(1),e1->Iso()/e1->Et());
//-----------------------------------------------------------------------------
// distributions for the component of MET perpendicular to the muon 
// momentum
//-----------------------------------------------------------------------------
	q[0] = (metx[0]*m1->Py()-mety[0]*m1->Px())/m1->Pt();
	q[1] = (metx[1]*m1->Py()-mety[1]*m1->Px())/m1->Pt();

	Hist.fMetPerp[0]->Fill(q[0]);
	Hist.fMetPerp[1]->Fill(q[1]);

      }
    }
  }

//   TLorentzVector m12;
//   double  mass;
//   for (int i1=0; i1<nemo-1; i1++) {
//     TStnMuon* e1 = fMuonBlock->Muon(i1);
//     const TLorentzVector* m1 = e1->Momentum();
//     for (int i2=i1+1; i2<nemo; i2++) {
//       TStnMuon* e2 = fMuonBlock->Muon(i2);
//       const TLorentzVector* m2 = e2->Momentum();
//       m12 =  *m1;
//       m12 += *m2;
//       mass = m12.M();
//       Hist.fMass->Fill(mass);
//     }
//   }

//   Hist.fNEmo10->Fill(ne10);
}


//_____________________________________________________________________________
int TMuoAnaModule::Event(int ientry) {
  
  Int_t           id_word;
//-----------------------------------------------------------------------------
//  read the necessary data
//-----------------------------------------------------------------------------
  fMuonBlock->GetEntry(ientry);
  fTrackBlock->GetEntry(ientry);
  fTriggerBlock->GetEntry(ientry);
  fMetBlock->GetEntry(ientry);
  fJetBlock->GetEntry(ientry);
  fTauBlock->GetEntry(ientry);

  if (PrintLevel() == 1) {
    GetHeaderBlock()->Print("MUON BLOCK for all the events");
    fMuonBlock->Print();
  }
//-----------------------------------------------------------------------------
//  always fill histogram for the run number
//-----------------------------------------------------------------------------
  fHist.fRunNumber->Fill(GetHeaderBlock()->RunNumber());
//-----------------------------------------------------------------------------
//  dont forget to clear the list
//-----------------------------------------------------------------------------
  fNRecoMuons = fMuonBlock->NMuons();
  fListOfGoodMuons->Clear();
  fNGoodMuons = 0;

  //  TLorentzVector* mom;
  for (int i=0; i<fNRecoMuons; i++) {
    TStnMuon* muo = fMuonBlock->Muon(i);
    //    mom = (TLorentzVector*) muo->Momentum();
//-----------------------------------------------------------------------------
// apply ID cuts
//-----------------------------------------------------------------------------
    if      (fCutMode == 0) id_word = MuoIDWord      (muo);
    fMuoIDWord[i] = id_word;
//-----------------------------------------------------------------------------
// and fill muon ID histograms:
//-----------------------------------------------------------------------------
    FillMuonIDHistograms(fHist.fMuoID[0],muo);
//-----------------------------------------------------------------------------
// muon histograms: set 0 = all CdfMuons
//-----------------------------------------------------------------------------
    FillMuonHistograms(fHist.fMuo[0],muo);

    if (fMuoIDWord[i] == 0) { 
//-----------------------------------------------------------------------------
// set 1: GOOD central muons
//-----------------------------------------------------------------------------
      fNGoodMuons++;
      fListOfGoodMuons->Add(muo);
      FillMuonHistograms(fHist.fMuo[1],muo);
//-----------------------------------------------------------------------------
// DEBUG
//-----------------------------------------------------------------------------
      if (PrintLevel() == 101) {
	GetHeaderBlock()->Print("All GOOD MUONS");
      }
    }
    else {
//-----------------------------------------------------------------------------
// set 2: not-GOOD, but central muons
//-----------------------------------------------------------------------------
      FillMuonHistograms(fHist.fMuo[2],muo);
    }
//-----------------------------------------------------------------------------
// set 3: debug, muons with EM energy=0
//-----------------------------------------------------------------------------
    if ( muo->EmEnergy() < 0.01) {
      FillMuonHistograms(fHist.fMuo[3],muo);
    }
//-----------------------------------------------------------------------------
// set 4: debug, muons with EM+HAD energy=0
//-----------------------------------------------------------------------------
    if ( muo->EmEnergy()+muo->HadEnergy() < 0.01) {
      FillMuonHistograms(fHist.fMuo[4],muo);
    }
//-----------------------------------------------------------------------------
// set 5: CMIO's
//-----------------------------------------------------------------------------
    if ( muo->Detector() == 0) {
      FillMuonHistograms(fHist.fMuo[5],muo);
    }
  }
//-----------------------------------------------------------------------------
//  event (MET and other histograms)
//-----------------------------------------------------------------------------
  FillEventHistograms(fHist.fEvent);
//-----------------------------------------------------------------------------
//  filtering part
//-----------------------------------------------------------------------------
  int passed = (fNGoodMuons >= fMinNGoodMuons);

  SetPassed(passed);

  fHist.fFilterResult->Fill(passed);
  return 0;			       
}


//_____________________________________________________________________________
int TMuoAnaModule::EndJob() {
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}


//_____________________________________________________________________________
void TMuoAnaModule::PlotHistograms(int run_number, int slide) {
  // plot slides

  char name[120], canvas_name[120], title[120], pad_name[120];

  sprintf(name,"run_%i_%i",run_number,slide);

  sprintf(title,"run %i slide %i ",run_number, slide);

  sprintf(canvas_name,"%s_%s",GetName(),name);
  sprintf(pad_name,"%s_p1",canvas_name);

  if (slide == 1) {
//-----------------------------------------------------------------------------
//  
//-----------------------------------------------------------------------------
    //    TPostScript ps("l1ana.ps",-111);

    TCanvas* c = NewSlide(name,title,2,3);
    TPad* p1 = (TPad*) c->GetPrimitive(pad_name);

//      p1->cd(1); 
//      fHist.fMuo[0].fEtVsEta->Draw();

//      p1->cd(3); 
//      fHist.fMuo[0].fEtVsEta->ProjectionX()->Draw();

//      p1->cd(5); 
//      fHist.fMuo[0].fEtVsEta->ProjectionY()->Draw();

//      p1->cd(2); 
//      fHist.fMuo[0].fPhi->Draw();

//      p1->cd(4); 
//      fHist.fMuo[0].fN10->Draw();

//      p1->cd(6); 
//      fHist.fMuo[0].fMass->Draw();

    gPad->Update();
  }
  else if (slide == 2) {
//-----------------------------------------------------------------------------
//  
//-----------------------------------------------------------------------------
    //    TPostScript ps("l1ana.ps",-111);

    TCanvas* c = NewSlide(name,title,2,3);
    TPad* p1 = (TPad*) c->GetPrimitive(pad_name);

//      p1->cd(1); 
//      fHist.fMuo[0].fIso->Draw();

//      p1->cd(2); 
//      fHist.fMuo[0].fDelX->Draw();

//      p1->cd(3); 
//      fHist.fMuo[0].fDelZ->Draw();

//      p1->cd(4); 
//      fHist.fMuo[0].fEOverPPos->Draw();
//      fHist.fMuo[0].fEOverPNeg->Draw("same");

//      p1->cd(5); 
//      fHist.fMuo[0].fHadEm->Draw();

    gPad->Update();
  }
}