///////////////////////////////////////////////////////////////////////////////
// Dec 27 2002 P.Murat: tau ID cuts
//
// if tau candidate passes all the cuts, ID word is 0, a bit set to 1
// represents failed cut
// angles: da(trk-seed) < max(fMaxTrkAngle[0], fMaxTrkAngle[1]/E(vis))
// angles: da(pi0-seed) < max(fMaxPi0Angle[0], fMaxPi0Angle[1]/E(vis))
// ID histograms: [0] : no ID cuts
//                [1] : all ID cuts except cut on this parameter
//                [2] : before applying this cut in sequence
//                [3] : after applying this cut insequence
//                [4] : all the cuts applied
///////////////////////////////////////////////////////////////////////////////
#include <cmath>
#include "TH1F.h"
#include "Stntuple/alg/TStnTauID.hh"
#include "Stntuple/obj/TStnTau.hh"
#include "Stntuple/obj/TStnTrack.hh"
#include "Stntuple/alg/TStntuple.hh"

ClassImp(TStnTauID)

//_____________________________________________________________________________
void TStnTauID::ReadV1(TBuffer& R__b) {

  struct TStnTauID_V1_t {
    TObject    fObject;
    Float_t    fMinEt;
    Float_t    fMinSeedTrackPt;
    Float_t    fMinSeedTowerEt;
    Float_t    fMaxCalMass;
    Float_t    fMaxTrkMass;
    Float_t    fMaxVisMass;
    Float_t    fMaxTrkAngle[2];
    Float_t    fMaxPi0Angle[2];
    
    Float_t    fMaxEmfr;
    Float_t    fKsi;

    Float_t    fMaxDetEta;
    Float_t    fMinZCes;
    Float_t    fMaxZCes;
    
    Int_t      fMinNAxHits;	// min number of hits on a segment
    Int_t      fMinNStHits;
    Int_t      fMinNAxSeg;	// min number of segments on a track
    Int_t      fMinNStSeg;
    Float_t    fMaxSeedTrackD0;
    Float_t    fMaxSeedTrackZ0;
    Float_t    fMaxSeedTrackDz;
    Float_t    fMaxCalIso;
    Float_t    fMaxCalIso1;
    Float_t    fMaxTrkIso;
    Float_t    fMaxTrkIso1;
    Float_t    fMaxNTrk1030;
    Float_t    fMaxNPi01030;
    Float_t    fMaxNMuStubs;
//-----------------------------------------------------------------------------
//  modes
//-----------------------------------------------------------------------------
    Int_t      fSelectElectrons;
    Int_t      fSelectCosmics;
  } id ;


  id.fObject.Streamer(R__b);
  R__b >> id.fMinEt;
  R__b >> id.fMinSeedTrackPt;
  R__b >> id.fMinSeedTowerEt;
  R__b >> id.fMaxCalMass;
  R__b >> id.fMaxTrkMass;
  R__b >> id.fMaxVisMass;
  R__b.ReadStaticArray(id.fMaxTrkAngle);
  R__b.ReadStaticArray(id.fMaxPi0Angle);
  R__b >> id.fMaxEmfr;
  R__b >> id.fKsi;
  R__b >> id.fMaxDetEta;
  R__b >> id.fMinZCes;
  R__b >> id.fMaxZCes;
  R__b >> id.fMinNAxHits;
  R__b >> id.fMinNStHits;
  R__b >> id.fMinNAxSeg;
  R__b >> id.fMinNStSeg;
  R__b >> id.fMaxSeedTrackD0;
  R__b >> id.fMaxSeedTrackZ0;
  R__b >> id.fMaxSeedTrackDz;
  R__b >> id.fMaxCalIso;
  R__b >> id.fMaxCalIso1;
  R__b >> id.fMaxTrkIso;
  R__b >> id.fMaxTrkIso1;
  R__b >> id.fMaxNTrk1030;
  R__b >> id.fMaxNPi01030;
  R__b >> id.fMaxNMuStubs;
  R__b >> id.fSelectElectrons;
  R__b >> id.fSelectCosmics;
  fUseMask = 0xFFFFFFFF;
//-----------------------------------------------------------------------------
// in principle need to copy back - later !
//-----------------------------------------------------------------------------
}

//_____________________________________________________________________________
void TStnTauID::ReadV2(TBuffer& R__b) {

  struct TStnTauID_V2_t {
    TNamed     fObject;
//-----------------------------------------------------------------------------
//  modes
//-----------------------------------------------------------------------------
    Int_t      fSelectElectrons;
    Int_t      fSelectCosmics;
    Int_t      fMinNAxHits;	// min number of hits on a segment
    Int_t      fMinNStHits;
    Int_t      fMinNAxSeg;	// min number of segments on a track
    Int_t      fMinNStSeg;
//-----------------------------------------------------------------------------
// cut values
//-----------------------------------------------------------------------------
    Float_t    fMinEt;
    Float_t    fMinSeedTrackPt;
    Float_t    fMinSeedTowerEt;
    Float_t    fMaxCalMass;
    Float_t    fMaxTrkMass;
    Float_t    fMaxVisMass;
    Float_t    fMaxTrkAngle[2];
    Float_t    fMaxPi0Angle[2];

    Float_t    fMaxEmfr;
    Float_t    fKsi[2];

    Float_t    fMaxDetEta;
    Float_t    fMinZCes;
    Float_t    fMaxZCes;

    Float_t    fMaxSeedTrackD0;
    Float_t    fMaxSeedTrackZ0;
    Float_t    fMaxSeedTrackDz;
    Float_t    fMaxCalIso;
    Float_t    fMaxCalIso1;
    Float_t    fMaxTrkIso;
    Float_t    fMaxTrkIso1;
    Float_t    fMaxNTrk1030;
    Float_t    fMaxNPi01030;
    Float_t    fMaxNMuStubs;

    void*      fEOR;		// ! end-of-record
  } id;

  int nwi = ((Int_t*  )&id.fMinEt)-&id.fSelectElectrons;
  int nwf = ((Float_t*)&id.fEOR  )-&id.fMinEt;

  id.fObject.Streamer(R__b);
  R__b.ReadFastArray(&id.fSelectElectrons,nwi);
  R__b.ReadFastArray(&id.fMinEt          ,nwf);
  fUseMask = 0xFFFFFFFF;
//-----------------------------------------------------------------------------
//  and copy it into the object if necessary and when necessary
//-----------------------------------------------------------------------------
}


//______________________________________________________________________________
void TStnTauID::Streamer(TBuffer &R__b) {
  // Stream an object of class TStnTauID, current version 2

  //  UInt_t R__s, R__c;

  int nwi = ((Int_t*  )&fMinVisEt)-&fSelectElectrons;
  int nwf = ((Float_t*)&fEOR     )-&fMinVisEt;

  if (R__b.IsReading()) {
    Version_t R__v = R__b.ReadVersion(); 

    if      (R__v == 1) ReadV1(R__b);
    else if (R__v == 2) ReadV2(R__b);
    else if (R__v >= 3) { 
//-----------------------------------------------------------------------------
//  version 3,4
//-----------------------------------------------------------------------------
      TNamed::Streamer(R__b);
      R__b.ReadFastArray(&fSelectElectrons,nwi);
      R__b.ReadFastArray(&fMinVisEt       ,nwf);
    }
  } 
  else {
    R__b.WriteVersion(TStnTauID::IsA());
    TNamed::Streamer(R__b);
    R__b.WriteFastArray(&fSelectElectrons,nwi);
    R__b.WriteFastArray(&fMinVisEt       ,nwf);
  }
}

//_____________________________________________________________________________
TStnTauID::TStnTauID(const char* Name, const char* Title): TNamed(Name,Title) {
  // kinematics cuts: mostly, but not exactly, correspond to L3 TAU_MET filter

  fMinSeedTrackPt =  4.5;
  fMinSeedTowerEt = 10. ;
  fMinVisEt       = 20. ;
  fMinCalEt       =  0. ;	// do not cut by default, L3 does at 20
  fMaxDetEta      = 1.0 ;	// no such cut really at L3, need for analysis
  fMaxTrkMass     = 2.0 ;
  fMaxCalMass     = 4.0 ;	// this is the 4.3.0 L3 TAU_MET cut 
  fMaxVisMass     = 5.0 ;	// offline
  fKsi[0]         = 1.00;
  fKsi[1]         = 0.15;
					// fiducial cuts
  fMinZCes        = 9;
  fMaxZCes        = 230;
					// isolation cuts (track iso is the 
					// relative isolation from the very 
					// beginning)
  fMaxNTrk1030    = 0;
  fMaxNPi01030    = 10;
  fMaxCalIso      = 100.;
  fMaxCalIso1     = 0.1;
  fMaxTrkIso      = 1.e6;		// turn track iso cut off by default

  fMinNAxHits     = 5;			// these are offline cuts 
  fMinNStHits     = 5;			// for 4.8.4 tracking
  fMinNAxSeg      = 3;			// for 5.x tracking should become
  fMinNStSeg      = 3;			// 5,5, 3,3

  fMaxSeedTrackD0 = 0.2;	        // same as for high-Pt electrons and muons
  fMaxSeedTrackZ0 = 60.;
  fMaxSeedTrackDz = 60.;		// what is this? 

  fMaxTrkAngle[0] = 100.;		// John wants 10 and 5 (?)
  fMaxTrkAngle[1] = 100.;
					// don't cut on pi0 for the moment
  fMaxPi0Angle[0] = 100.;
  fMaxPi0Angle[1] = 100.;

  fMaxTrkEta      = 1.e6;	        // unfortunately used at L3 in 2002
  fMaxNTowers     = 100;

  fSelectElectrons = 0;
  fSelectCosmics   = 0;

  for (int i=0; i<4; i++) fInteger[0] = 0;
  for (int i=0; i<5; i++) fFloat  [0] = 0;
  fUseMask = 0xFFFFFFFF;
}


//_____________________________________________________________________________
TStnTauID::~TStnTauID() {
}


//_____________________________________________________________________________
int TStnTauID::IsElectron(TStnTau* Tau) {
//-----------------------------------------------------------------------------
// electron rejection
// fSelectElectrons = 0       : standard "John-type" electron rejection cut
// fSelectElectrons = 2       : use E/P(seed) rather than E/sum(P)
//                  = 11      : do nothing
// fSelectElectrons = 100+mode: inverse of 'mode'

  TLorentzVector*   cmom;
  float             cal_e, cal_et, ep, esump;
  int               mode;
  int               is_electron = 0;

  mode     = fSelectElectrons % 100;

  cmom     = Tau->ClusterMomentum();
  cal_e    = cmom->E();
  cal_et   = cmom->Pt()*cal_e/cmom->P();
  esump    = cal_e/Tau->SumTrackP();
  ep       = cal_et/Tau->SeedTrackPt();

  if      (mode == 0) {
    if (Tau->Emfr()      > fKsi[0] - fKsi[1]/esump ) is_electron = 1;
  }
  else if (mode == 2) {
    if (Tau->Emfr()      > fKsi[0] - fKsi[1]/ep    ) is_electron = 1;
  }
  else if (mode == 11) {
  }

  if (fSelectElectrons > 100) is_electron = 1-is_electron;

  return (is_electron & fUseMask);
}

//_____________________________________________________________________________
int TStnTauID::IDWord(TStnTau* Tau) {

  TLorentzVector  *vmom, *tmom;
  TStnTrack*      seed_track;

  Float_t         d0, z0, abs_deta, max_angle, trk_eta, cal_et, abs_teta;
  Float_t         vis_e, vis_et, xces, zces, z_cot, chi2_cot, z_cem;
  Int_t           n_ax_seg, n_st_seg;
  Int_t           id_word = 0;
//-----------------------------------------------------------------------------
// precalculate some numbers, use visible momentum to cut on
//-----------------------------------------------------------------------------
  vmom       = Tau->VisMomentum();
  tmom       = Tau->TrackMomentum();

  trk_eta    = tmom->Eta();
  cal_et     = Tau->CalEt();

  vis_e      = vmom->E();
  vis_et     = vmom->Pt()*vis_e/vmom->P();
  abs_deta   = TMath::Abs(Tau->DetEta());
  abs_teta   = TMath::Abs(trk_eta);
  seed_track = Tau->SeedTrack();
//-----------------------------------------------------------------------------
// calorimetry
//-----------------------------------------------------------------------------
  if (abs_deta           > fMaxDetEta     )           id_word |= kDetEtaBit;
  if (abs_teta           > fMaxTrkEta     )           id_word |= kTrkEtaBit;
  if (cal_et             < fMinCalEt      ) {
    id_word |= kCalEtBit ;
  }
  if (vis_et             < fMinVisEt      )           id_word |= kVisEtBit ;
  if (Tau->TowerEt(0)    < fMinSeedTowerEt)           id_word |= kSeedTowerEtBit;
  if (Tau->SeedTrackPt() < fMinSeedTrackPt)           id_word |= kSeedTrackPtBit;
//-----------------------------------------------------------------------------
// check EM fraction - it may depend on both - cal and trk 
// for electron rejection cuts we have whole spectrum of choices, 
// for the moment use John's thing
//-----------------------------------------------------------------------------
  if (IsElectron(Tau))                                id_word |= kEmfrBit;
//-----------------------------------------------------------------------------
// tau-specific isolation - check NT1030 first because this is done at L3, do
// not cut on the number of pi0's for the moment
//-----------------------------------------------------------------------------
  if (Tau->NTracks1030()           > fMaxNTrk1030   ) id_word |= kNTrk1030Bit;
  if (Tau->NPi01030()              > fMaxNPi01030   ) id_word |= kNPi01030Bit;
//-----------------------------------------------------------------------------
// "normal" isolation cuts
//-----------------------------------------------------------------------------
  double cal_iso1 = Tau->CalIso04()/Tau->Et();

  if (Tau->TrackIso()              > fMaxTrkIso     ) id_word |= kTrkIsoBit;
  if (Tau->CalIso04()              > fMaxCalIso     ) id_word |= kCalIsoBit;
  if (cal_iso1                     > fMaxCalIso1    ) id_word |= kCalIso1Bit;
//-----------------------------------------------------------------------------
// seed track quality cuts
//-----------------------------------------------------------------------------
  n_ax_seg = -1;
  n_st_seg = -1;
  d0       = 1.e5;
  z0       = TMath::Abs(Tau->Zv());
//-----------------------------------------------------------------------------
//  count the segments "by hand": number of hits issue
//-----------------------------------------------------------------------------
  if (seed_track) { 
    d0       = TStntuple::CorrectedD0(seed_track);
    n_ax_seg = seed_track->NCotAxSeg(fMinNAxHits);
    n_st_seg = seed_track->NCotStSeg(fMinNStHits);
  }

  z_cot     = seed_track->Z0()+seed_track->Lam0()*112.;
  chi2_cot  = seed_track->Chi2Cot()/(seed_track->NCotHitsTot()-4.9999);

  if (n_ax_seg           < fMinNAxSeg ) id_word |= kNAxSegBit;
  if (n_st_seg           < fMinNStSeg ) id_word |= kNStSegBit;
  if (fabs(z_cot)        > 155.       ) id_word |= kCotFidBit;
  if (chi2_cot           > 3.         ) id_word |= kCotChi2Bit;

  if (fSelectCosmics == 0) {
    if (TMath::Abs(d0)   > fMaxSeedTrackD0) id_word |= kSeedTrackD0Bit;
  }
  else {
    if (TMath::Abs(d0)   < fMaxSeedTrackD0) id_word |= kSeedTrackD0Bit;
  }

  if (z0                 > fMaxSeedTrackZ0) id_word |= kSeedTrackZ0Bit;
//-----------------------------------------------------------------------------
// now that we know the seed track is OK, check Z(CES)
// for 1-prong tau candidates require also the seed to be CES fiducial
//-----------------------------------------------------------------------------
  zces     = Tau->SeedTrackZCes();
  if (( zces < fMinZCes) || (zces > fMaxZCes)) id_word |= kZCesBit;
//-----------------------------------------------------------------------------
//  check that at R_CEM_MIN Z coordinate of the tau candidate 
//  is also within the CES fiducial .. .assume that ZCes is always local
//-----------------------------------------------------------------------------
  if (seed_track->Side() == 0) {
    z_cem = zces+seed_track->Lam0()*(184.15-167.64);
  }
  else {
    z_cem = zces-seed_track->Lam0()*(184.15-167.64);
  }
  
  if (fabs(z_cem) < fMinZCes)                  id_word |= kZCesBit;

  if (Tau->NTracks10() == 1) {
    xces     = seed_track->XCes();
    if (TMath::Abs(xces) > 21.5)               id_word |= kXCesBit;
  }
//-----------------------------------------------------------------------------
// mass cuts
//-----------------------------------------------------------------------------
  if (Tau->Mass     ()   > fMaxCalMass    )  id_word |= kCalMassBit;
  if (Tau->TrackMass()   > fMaxTrkMass    )  id_word |= kTrkMassBit;
  if (Tau->VisMass  ()   > fMaxVisMass    )  id_word |= kVisMassBit;
//-----------------------------------------------------------------------------
// angles - John's variables
//-----------------------------------------------------------------------------
  max_angle = TMath::Min(fMaxTrkAngle[0],fMaxTrkAngle[1]/vis_e);
  if (Tau->MaxTrkAngle() > max_angle      )  id_word |= kTrkAngleBit;

  max_angle = TMath::Min(fMaxPi0Angle[0],fMaxPi0Angle[1]/vis_e);
  if (Tau->MaxPi0Angle() > max_angle      )  id_word |= kPi0AngleBit;
//-----------------------------------------------------------------------------
// cut on N(towers)
//-----------------------------------------------------------------------------
  if (Tau->NTowers  ()   > fMaxNTowers    )  id_word |= kNTowersBit;

  return (id_word & fUseMask);
}

//_____________________________________________________________________________
int TStnTauID::LooseIDWord(TStnTau* Tau) {
  // "loose" tau id: tau object above the thresholds with electrons removed

  TLorentzVector  *vmom, *cmom;
  TStnTrack*      seed_track;

  Float_t         cal_et, vis_et, vis_e, abs_deta, d0, zces;
//    Int_t           n_ax_seg, n_st_seg;
  Int_t           id_word = 0;
//-----------------------------------------------------------------------------
// precalculate some numbers, use visible momentum to cut on
//-----------------------------------------------------------------------------
  vmom      = Tau->VisMomentum();
  vis_e     = vmom->E();
  vis_et    = vmom->Pt()*vis_e/vmom->P();
  cmom      = Tau->ClusterMomentum();
  cal_et    = cmom->Pt()*cmom->E()/cmom->P();

  abs_deta  = TMath::Abs(Tau->DetEta());
//    esump    = e/Tau->SumTrackP();

//    //  itau  = Tau->Number();

  seed_track = Tau->SeedTrack();
//-----------------------------------------------------------------------------
// calorimetry
//-----------------------------------------------------------------------------
  if (abs_deta           > fMaxDetEta     )  id_word |= kDetEtaBit;
  if (cal_et             < fMinCalEt      )  id_word |= kCalEtBit ;
  if (vis_et             < fMinVisEt      )  id_word |= kVisEtBit ;
  if (Tau->TowerEt(0)    < fMinSeedTowerEt)  id_word |= kSeedTowerEtBit;
  if (Tau->SeedTrackPt() < fMinSeedTrackPt)  id_word |= kSeedTrackPtBit;
  d0       = 1.e5;
//-----------------------------------------------------------------------------
//  always cut on impact parameter to reject cosmics
//-----------------------------------------------------------------------------
  if (seed_track) { 
    d0       = TStntuple::CorrectedD0(seed_track);
  }

  if (fSelectCosmics == 0) {
    if (TMath::Abs(d0)   > fMaxSeedTrackD0) id_word |= kSeedTrackD0Bit;
  }
  else {
    if (TMath::Abs(d0)   < fMaxSeedTrackD0) id_word |= kSeedTrackD0Bit;
  }
//-----------------------------------------------------------------------------
// check ZCES to reject electron background
//-----------------------------------------------------------------------------
  zces     = Tau->SeedTrackZCes();
  if (( zces < fMinZCes) || (zces > fMaxZCes)) id_word |= kZCesBit;
//-----------------------------------------------------------------------------
// check EM fraction - it may depend on both - cal and trk 
// for electron rejection cuts we have whole spectrum of choices, 
// for the moment use John's thing
// do not reject electrons for the purpose of loose tau ID
//-----------------------------------------------------------------------------
  //  if (IsElectron(Tau))                       id_word |= kEmfrBit;

  return (id_word & fUseMask);
}



//_____________________________________________________________________________
void TStnTauID::FillHistograms(Hist_t& Hist, TStnTau* Tau) {
  // this method uses only bit assignments, thus declared static

  float       d0, vis_et, cal_et, trk_eta, cot_z, cot_chi2;
  int         n_ax_seg, n_st_seg, id_word, tight_ele_bit, loose_ele_bit;

  TStnTrack*  seed_track;

  TLorentzVector* tau_mom  = (TLorentzVector*) Tau->VisMomentum();
  //  TLorentzVector* cmom     = (TLorentzVector*) Tau->ClusterMomentum();
  TLorentzVector* tmom     = (TLorentzVector*) Tau->TrackMomentum();

  vis_et  = tau_mom->Pt()*tau_mom->E()/tau_mom->P();
  cal_et  = Tau->CalEt();
  trk_eta = tmom->Eta();

  id_word = Tau->IDWord();
//-----------------------------------------------------------------------------
//  L3 cuts on detector eta ..
//-----------------------------------------------------------------------------
  Hist.fDetEta[0]->Fill(Tau->DetEta());
  if ((id_word & ~kDetEtaBit) == 0) Hist.fDetEta[1]->Fill(Tau->DetEta());
  if (id_word == 0) Hist.fDetEta[4]->Fill(Tau->DetEta());

  Hist.fTrkEta[0]->Fill(trk_eta);
  if ((id_word & ~kTrkEtaBit) == 0) Hist.fTrkEta[1]->Fill(trk_eta);
  if (id_word == 0) Hist.fTrkEta[4]->Fill(trk_eta);
//-----------------------------------------------------------------------------
//  ET/PT cuts
//-----------------------------------------------------------------------------
  Hist.fCalEt[0]->Fill(cal_et);
  if ((id_word & ~kCalEtBit) == 0) Hist.fCalEt[1]->Fill(cal_et);
  if (id_word == 0) Hist.fCalEt[4]->Fill(cal_et);

  Hist.fVisEt[0]->Fill(vis_et);
  if ((id_word & ~kVisEtBit) == 0) Hist.fVisEt[1]->Fill(vis_et);
  if (id_word == 0) Hist.fVisEt[4]->Fill(vis_et);

  Hist.fSeedTowerEt[0]->Fill(Tau->TowerEt(0));
  if ((id_word & ~kSeedTowerEtBit) == 0)
    Hist.fSeedTowerEt[1]->Fill(Tau->TowerEt(0));
  if (id_word == 0) Hist.fSeedTowerEt[4]->Fill(Tau->TowerEt(0));

  Hist.fSeedTrackPt[0]->Fill(Tau->SeedTrackPt());
  if ((id_word & ~kSeedTrackPtBit) == 0) 
    Hist.fSeedTrackPt[1]->Fill(Tau->SeedTrackPt());
  if (id_word == 0) Hist.fSeedTrackPt[4]->Fill(Tau->SeedTrackPt());
//-----------------------------------------------------------------------------
//  electron removal - see CDF-3932, for example... choice of 0.15 is 
//  not critical, could be 0.10 as well
//-----------------------------------------------------------------------------
  Hist.fEmfr[0]->Fill(Tau->Emfr());
  if ((id_word & ~ kEmfrBit) == 0) Hist.fEmfr[1]->Fill(Tau->Emfr());
  if (id_word == 0) Hist.fEmfr[4]->Fill(Tau->Emfr());
//-----------------------------------------------------------------------------
//  N(tracks and pi0's)
//-----------------------------------------------------------------------------
  Hist.fNTrk1030[0]->Fill(Tau->NTracks1030());
  if ((id_word & ~ kNTrk1030Bit) == 0) 
    Hist.fNTrk1030[1]->Fill(Tau->NTracks1030());
  if (id_word == 0) Hist.fNTrk1030[4]->Fill(Tau->NTracks1030());

  Hist.fNPi01030[0]->Fill(Tau->NPi01030());
  if ((id_word & ~ kNPi01030Bit) == 0) 
    Hist.fNPi01030[1]->Fill(Tau->NPi01030());
  if (id_word == 0) Hist.fNPi01030[4]->Fill(Tau->NPi01030());
//-----------------------------------------------------------------------------
// isolation cuts
// cut on cal iso is the most suspicious one, apply it last
//-----------------------------------------------------------------------------
  Hist.fCalorIso04[0]->Fill(Tau->CalIso04());
  if ((id_word & ~ kCalIsoBit) == 0) 
    Hist.fCalorIso04[1]->Fill(Tau->CalIso04());
  if (id_word == 0) Hist.fCalorIso04[4]->Fill(Tau->CalIso04());

  Hist.fCalIso1[0]->Fill(Tau->CalIso1());
  if ((id_word & ~ kCalIso1Bit) == 0) 
    Hist.fCalIso1[1]->Fill(Tau->CalIso1());
  if (id_word == 0) Hist.fCalIso1[4]->Fill(Tau->CalIso1());

  Hist.fTrackIso04[0]->Fill(Tau->TrackIso());
  if ((id_word & ~ kTrkIsoBit) == 0) 
    Hist.fTrackIso04[1]->Fill(Tau->TrackIso());
  if (id_word == 0) Hist.fTrackIso04[4]->Fill(Tau->TrackIso());
//-----------------------------------------------------------------------------
//  require seed track to be good enough
//-----------------------------------------------------------------------------
  seed_track = Tau->SeedTrack();
  if (seed_track) {
    d0        = TStntuple::CorrectedD0(seed_track);
    cot_z     = seed_track->Z0()+seed_track->Lam0()*112.;
    cot_chi2  = seed_track->Chi2Cot()/(seed_track->NCotHitsTot()-4.9999);
  }
  else {
    d0        = 999.99;
    cot_z     = 999.99;
    cot_chi2  = 999.99;
  }

  n_ax_seg = seed_track->NAxSeg();
  n_st_seg = seed_track->NStSeg();

  Hist.fSeedTrackNAxSeg[0]->Fill(n_ax_seg);
  if ((id_word & ~kNAxSegBit) == 0) 
    Hist.fSeedTrackNAxSeg[1]->Fill(n_ax_seg);
  if (id_word == 0) Hist.fSeedTrackNAxSeg[4]->Fill(n_ax_seg);

  Hist.fSeedTrackNStSeg[0]->Fill(n_st_seg);
  if ((id_word & ~kNStSegBit) == 0) 
    Hist.fSeedTrackNStSeg[1]->Fill(n_st_seg);
  if (id_word == 0) Hist.fSeedTrackNStSeg[4]->Fill(n_st_seg);

  Hist.fSeedTrackD0[0]->Fill(d0);
  if ((id_word & ~kSeedTrackD0Bit) == 0) Hist.fSeedTrackD0[1]->Fill(d0);
  if (id_word == 0) Hist.fSeedTrackD0[4]->Fill(d0);

  Hist.fSeedTrackZ0[0]->Fill(Tau->SeedTrackZ());
  if ((id_word & ~kSeedTrackZ0Bit) == 0) 
    Hist.fSeedTrackZ0[1]->Fill(Tau->SeedTrackZ());
  if (id_word == 0) Hist.fSeedTrackZ0[4]->Fill(Tau->SeedTrackZ());

  Hist.fSeedTrackDZ0[0]->Fill(seed_track->DZ0());
  if ((id_word & ~kSeedTrackDzBit) == 0) 
    Hist.fSeedTrackDZ0[1]->Fill(seed_track->DZ0());
  if (id_word == 0) Hist.fSeedTrackDZ0[4]->Fill(seed_track->DZ0());

  Hist.fSeedTrackCotZ[0]->Fill(cot_z);
  if ((id_word & ~kCotFidBit) == 0) 
    Hist.fSeedTrackCotZ[1]->Fill(cot_z);
  if (id_word == 0) Hist.fSeedTrackCotZ[4]->Fill(cot_z);

  Hist.fSeedTrackCotChi2[0]->Fill(cot_chi2);
  if ((id_word & ~kCotChi2Bit) == 0) 
    Hist.fSeedTrackCotChi2[1]->Fill(cot_chi2);
  if (id_word == 0) Hist.fSeedTrackCotChi2[4]->Fill(cot_chi2);
//-----------------------------------------------------------------------------
// CES coordinates
//-----------------------------------------------------------------------------
  Hist.fZCes[0]->Fill(Tau->SeedTrackZCes());
  if ((id_word & ~ kZCesBit) == 0) Hist.fZCes[1]->Fill(Tau->SeedTrackZCes());
  if (id_word == 0) Hist.fZCes[4]->Fill(Tau->SeedTrackZCes());

  Hist.fXCes[0]->Fill(Tau->SeedTrackXCes());
  if ((id_word & ~ kXCesBit) == 0) Hist.fXCes[1]->Fill(Tau->SeedTrackXCes());
  if (id_word == 0) Hist.fXCes[4]->Fill(Tau->SeedTrackXCes());
//-----------------------------------------------------------------------------
//  various masses
//-----------------------------------------------------------------------------
  Hist.fCalorMass[0]->Fill(Tau->Mass());
  if ((id_word & ~ kCalMassBit) == 0) Hist.fCalorMass[1]->Fill(Tau->Mass());
  if (id_word == 0) Hist.fCalorMass[4]->Fill(Tau->Mass());

  Hist.fTrackMass[0]->Fill(Tau->TrackMass());
  if ((id_word & ~ kTrkMassBit) == 0) 
    Hist.fTrackMass[1]->Fill(Tau->TrackMass());
  if (id_word == 0) Hist.fTrackMass[4]->Fill(Tau->TrackMass());

  Hist.fVisMass[0]->Fill(Tau->VisMass());
  if ((id_word & ~ kVisMassBit) == 0) 
    Hist.fVisMass[1]->Fill(Tau->VisMass());
  if (id_word == 0) Hist.fVisMass[4]->Fill(Tau->VisMass());
//-----------------------------------------------------------------------------
// John's max angles
//-----------------------------------------------------------------------------
  Hist.fMaxTrkAngle[0]->Fill(Tau->MaxTrkAngle());
  if ((id_word & ~kTrkAngleBit) == 0)
    Hist.fMaxTrkAngle[1]->Fill(Tau->MaxTrkAngle());
  if (id_word == 0) Hist.fMaxTrkAngle[4]->Fill(Tau->MaxTrkAngle());

  Hist.fMaxPi0Angle[0]->Fill(Tau->MaxPi0Angle());
  if ((id_word & ~kPi0AngleBit) == 0)
    Hist.fMaxPi0Angle[1]->Fill(Tau->MaxPi0Angle());
  if (id_word == 0) Hist.fMaxPi0Angle[4]->Fill(Tau->MaxPi0Angle());
//-----------------------------------------------------------------------------
// N(towers) in the calorimeter
//-----------------------------------------------------------------------------
  Hist.fNTowers[0]->Fill(Tau->NTowers());
  if ((id_word & ~kNTowersBit) == 0)
    Hist.fNTowers[1]->Fill(Tau->NTowers());
  if (id_word == 0) Hist.fNTowers[4]->Fill(Tau->NTowers());
//-----------------------------------------------------------------------------
// tight and loose electron ID bits
//-----------------------------------------------------------------------------
  tight_ele_bit = (id_word & kTightElectronBit) != 0;
  Hist.fTightEleBit[0]->Fill(tight_ele_bit);
  if ((id_word & ~kTightElectronBit) == 0) 
    Hist.fTightEleBit[1]->Fill(tight_ele_bit);
  if (id_word == 0) Hist.fTightEleBit[4]->Fill(tight_ele_bit);

  loose_ele_bit = (id_word & kLooseElectronBit) != 0 ;
  Hist.fLooseEleBit[0]->Fill(loose_ele_bit);
  Hist.fEleChi2    [0]->Fill(Tau->Tmp(37));
  if ((id_word & ~kLooseElectronBit) == 0) {
    Hist.fLooseEleBit[1]->Fill(loose_ele_bit);
    Hist.fEleChi2    [1]->Fill(Tau->Tmp(37));
  }
  if (id_word == 0) {
    Hist.fLooseEleBit[4]->Fill(loose_ele_bit);
    Hist.fEleChi2    [4]->Fill(Tau->Tmp(37));
  }
//-----------------------------------------------------------------------------
//    Hist.fNMuStubs[0]->Fill(Tau->NMuStubs());
//    if ((id_word & ~ kNMuStubsBit) == 0) 
//      Hist.fNMuStubs[1]->Fill(Tau->NMuStubs());
//    if (id_word == 0) Hist.fNMuStubs[4]->Fill(Tau->NMuStubs());
//-----------------------------------------------------------------------------
//  now histogram effect of cuts in the order they are applied
//  cal_eta is what L3 is cutting on
//-----------------------------------------------------------------------------
  Hist.fDetEta[2]->Fill(Tau->DetEta());
  if (id_word & kDetEtaBit)                                 goto END;
  Hist.fDetEta[3]->Fill(Tau->DetEta());

  Hist.fTrkEta[2]->Fill(trk_eta);
  if (id_word & kTrkEtaBit)                                 goto END;
  Hist.fTrkEta[3]->Fill(trk_eta);

  Hist.fCalEt[2]->Fill(cal_et);
  if (id_word & kCalEtBit)                                  goto END;
  Hist.fCalEt[3]->Fill(cal_et);

  Hist.fVisEt[2]->Fill(vis_et);
  if (id_word & kVisEtBit)                                  goto END;
  Hist.fVisEt[3]->Fill(vis_et);

  Hist.fSeedTowerEt[2]->Fill(Tau->TowerEt(0));
  if (id_word & kSeedTowerEtBit)                            goto END;
  Hist.fSeedTowerEt[3]->Fill(Tau->TowerEt(0));

  Hist.fSeedTrackPt[2]->Fill(Tau->SeedTrackPt());
  if (id_word & kSeedTrackPtBit)                            goto END;
  Hist.fSeedTrackPt[3]->Fill(Tau->SeedTrackPt());
//-----------------------------------------------------------------------------
//  electron removal - see CDF-3932, for example... choice of 0.15 is 
//  not critical, could be 0.10 as well
//-----------------------------------------------------------------------------
  Hist.fEmfr[2]->Fill(Tau->Emfr());
  if (id_word & kEmfrBit)                                   goto END;
  Hist.fEmfr[3]->Fill(Tau->Emfr());
//-----------------------------------------------------------------------------
// tau-specific isolation - check NT1030 first because this is done at L3,
// for the moment do not cut on the number of pi0's
//-----------------------------------------------------------------------------
  Hist.fNTrk1030[2]->Fill(Tau->NTracks1030());
  if (id_word & kNTrk1030Bit)                               goto END;
  Hist.fNTrk1030[3]->Fill(Tau->NTracks1030());

  Hist.fNPi01030[2]->Fill(Tau->NPi01030());
  if (id_word & kNPi01030Bit)                               goto END;
  Hist.fNPi01030[3]->Fill(Tau->NPi01030());
//-----------------------------------------------------------------------------
// "normal" calorimetry and track isolations
//-----------------------------------------------------------------------------
  Hist.fCalorIso04[2]->Fill(Tau->CalIso04());
  if (id_word & kCalIsoBit)                                 goto END;
  Hist.fCalorIso04[3]->Fill(Tau->CalIso04());

  Hist.fCalIso1[2]->Fill(Tau->CalIso1());
  if (id_word & kCalIso1Bit)                                goto END;
  Hist.fCalIso1[3]->Fill(Tau->CalIso1());

  Hist.fTrackIso04[2]->Fill(Tau->TrackIso());
  if (id_word & kTrkIsoBit)                                 goto END;
  Hist.fTrackIso04[3]->Fill(Tau->TrackIso());
//-----------------------------------------------------------------------------
//  require seed track to be good enough
//-----------------------------------------------------------------------------
  Hist.fSeedTrackNAxSeg[2]->Fill(n_ax_seg);
  if (id_word & kNAxSegBit)                                 goto END;
  Hist.fSeedTrackNAxSeg[3]->Fill(n_ax_seg);

  Hist.fSeedTrackNStSeg[2]->Fill(n_st_seg);
  if (id_word & kNStSegBit)                                 goto END;
  Hist.fSeedTrackNStSeg[3]->Fill(n_st_seg);
//-----------------------------------------------------------------------------
//  cut on the (corrected) impact parameter, Z0 and DZ of the seed track
//-----------------------------------------------------------------------------
  Hist.fSeedTrackD0[2]->Fill(d0);
  if (id_word & kSeedTrackD0Bit)                            goto END;
  Hist.fSeedTrackD0[3]->Fill(d0);

  Hist.fSeedTrackZ0[2]->Fill(Tau->SeedTrackZ());
  if (id_word & kSeedTrackZ0Bit)                            goto END;
  Hist.fSeedTrackZ0[3]->Fill(Tau->SeedTrackZ());

  Hist.fSeedTrackDZ0[2]->Fill(seed_track->DZ0());
  if (id_word & kSeedTrackDzBit)                            goto END;
  Hist.fSeedTrackDZ0[3]->Fill(seed_track->DZ0());

  Hist.fSeedTrackCotZ[2]->Fill(cot_z);
  if (id_word & kCotFidBit)                                 goto END;
  Hist.fSeedTrackCotZ[3]->Fill(cot_z);

  Hist.fSeedTrackCotChi2[2]->Fill(cot_chi2);
  if (id_word & kCotChi2Bit)                                goto END;
  Hist.fSeedTrackCotChi2[3]->Fill(cot_chi2);

  Hist.fZCes[2]->Fill(Tau->SeedTrackZCes());
  if (id_word & kZCesBit)                                   goto END;
  Hist.fZCes[3]->Fill(Tau->SeedTrackZCes());
//-----------------------------------------------------------------------------
//  now: masses
//-----------------------------------------------------------------------------
  Hist.fCalorMass[2]->Fill(Tau->Mass());
  if (id_word & kCalMassBit)                                goto END;
  Hist.fCalorMass[3]->Fill(Tau->Mass());

  Hist.fTrackMass[2]->Fill(Tau->TrackMass());
  if (id_word & kTrkMassBit)                                goto END;
  Hist.fTrackMass[3]->Fill(Tau->TrackMass());

  Hist.fVisMass[2]->Fill(Tau->VisMass());
  if (id_word & kVisMassBit)                                goto END;
  Hist.fVisMass[3]->Fill(Tau->VisMass());

//    Hist.fNMuStubs[2]->Fill(Tau->NMuStubs());
//    if (id_word & kNMuStubsBit )                goto END;
//    Hist.fNMuStubs[3]->Fill(Tau->NMuStubs());
//-----------------------------------------------------------------------------
//  John's angles
//-----------------------------------------------------------------------------
  Hist.fMaxTrkAngle[2]->Fill(Tau->MaxTrkAngle());
  if (id_word & kTrkAngleBit)                               goto END;
  Hist.fMaxTrkAngle[3]->Fill(Tau->MaxTrkAngle());

  Hist.fMaxPi0Angle[2]->Fill(Tau->MaxPi0Angle());
  if (id_word & kPi0AngleBit)                               goto END;
  Hist.fMaxPi0Angle[3]->Fill(Tau->MaxPi0Angle());
//-----------------------------------------------------------------------------
//  N(towers) cut
//-----------------------------------------------------------------------------
  Hist.fNTowers[2]->Fill(Tau->NTowers());
  if (id_word & kNTowersBit)                               goto END;
  Hist.fNTowers[3]->Fill(Tau->NTowers());
//-----------------------------------------------------------------------------
//  tight and loose electron bits
//-----------------------------------------------------------------------------
  Hist.fTightEleBit[2]->Fill(tight_ele_bit);
  if (id_word & kTightElectronBit)                         goto END;
  Hist.fTightEleBit[3]->Fill(tight_ele_bit);

  Hist.fLooseEleBit[2]->Fill(loose_ele_bit);
  Hist.fEleChi2    [2]->Fill(Tau->Tmp(37));
  if (id_word & kLooseElectronBit)                         goto END;
  Hist.fLooseEleBit[3]->Fill(loose_ele_bit);
  Hist.fEleChi2    [3]->Fill(Tau->Tmp(37));
//-----------------------------------------------------------------------------
// cut on N(muon hits) should not be used - we have too much background from
// the accelerator there
//  if (Tau->NMuHits()          > 0              ) goto END;
//-----------------------------------------------------------------------------
 END:;

//-----------------------------------------------------------------------------
//  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);
}



//_____________________________________________________________________________
void TStnTauID::Print(const char* Opt) const {
  Error("Print","not implemented yet");
}