///////////////////////////////////////////////////////////////////////////////
//  Initialization of STNTUPLE muon block
//  P.Murat
///////////////////////////////////////////////////////////////////////////////
#include "TROOT.h"
#include "TFolder.h"

#include "CalorGeometry/CalorKey.hh"
#include "CalorObjects/CalTower.hh"
#include <CalorObjects/CalData.hh>
#include <TrackingObjects/Storable/CdfTrackColl.hh>
#include <TrackingObjects/Storable/CdfTrackView.hh>
#include <MuonObjects/CdfMuonView.hh>
#include <MuonObjects/CdfMuonColl.hh>
#include <MuonObjects/MuonStubColl.hh>
#include <MuonObjects/MuonStub.hh>
#include "Muon/MuonFiducialTool.hh"
#include "MuonUser/MuonAnalysisTools.hh"
#include <AbsEnv/AbsEnv.hh>

#include <iostream>

#include <Stntuple/obj/TStnEvent.hh>
#include <Stntuple/obj/TStnErrorLogger.hh>
#include "Stntuple/obj/TStnMuonBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"

#include "Stntuple/alg/TStntuple.hh"

#include "Stntuple/mod/InitStntupleDataBlocks.hh"
#include "Stntuple/mod/StntupleUtilities.hh"

namespace {
  TFolder*          fgFolder    = NULL;
  TStnErrorLogger*  fgLogger    = NULL;
}

//_____________________________________________________________________________
Int_t StntupleInitMuon(TStnMuon* Muon, const CdfMuon* Ken) 
{
  // init muon object from a CdfMuon (Run II only)

  double phi, dr;

  HepLorentzVector mom = Ken->fourMomentum();
  Muon->fMomentum.SetXYZT(mom.px(),mom.py(),mom.pz(),mom.e());

  const CdfTrack* trk;

  trk                 = Ken->bestTrack().operator->();
  Muon->fCdfMuon      = (CdfMuon*) Ken;

  Muon->fTrind        = -1;
  Muon->fIdwrd        = 0;
  Muon->fIswrd        = 0;
  Muon->fFid          = 0;

  Muon->fCharge       = trk->charge();
  Muon->fHitMask      = 0;
  Muon->fHitMask      = (Ken->muonType() & 0xff) << 24;

  Muon->fCray         = TStnDataBlock::kUndefined;
  Muon->fCmuStub      = NULL;
  Muon->fCmpStub      = NULL;
  Muon->fCmxStub      = NULL;
  Muon->fBmuStub      = NULL;
//-----------------------------------------------------------------------------
// unused variables - future reserve
//-----------------------------------------------------------------------------
  for (int i=0; i<3; i++) {
    Muon->fCalEta[i] = TStnDataBlock::kUndefined; 
    Muon->fCalPhi[i] = TStnDataBlock::kUndefined; 
  }
//-----------------------------------------------------------------------------
//  Fiducial
//-----------------------------------------------------------------------------
  Muon->fCmuFidX        = TStnDataBlock::kUndefined;
  Muon->fCmuFidZ        = TStnDataBlock::kUndefined;
  Muon->fCmpFidX        = TStnDataBlock::kUndefined;
  Muon->fCmpFidZ        = TStnDataBlock::kUndefined;
  Muon->fCmxFidX        = TStnDataBlock::kUndefined;
  Muon->fCmxFidZ        = TStnDataBlock::kUndefined;
  Muon->fCmxFidX2       = TStnDataBlock::kUndefined;
  Muon->fCmxFidZ2       = TStnDataBlock::kUndefined;
  Muon->fBmuFidX        = TStnDataBlock::kUndefined;
  Muon->fBmuFidZ        = TStnDataBlock::kUndefined;
  MuonFiducialTool fidTool;
  UChar_t CMUFid = 0, CMUFid2 = 0;
  UChar_t CMPFid = 0, CMPFid2 = 0;
  UChar_t CMXFid = 0, CMXFid2 = 0;
  UChar_t CMXFidt= 0, CMXFid2t= 0;
  UChar_t BMUFid = 0, BMUFid2 = 0;
  double d=0,x=0,z=0;
  if(fidTool.isFiducial(*trk,MuonDigiCode::CMU,d,x,z)) {
    CMUFid=1;
    if(x<0.0 && z<-3.0) CMUFid2=1;
    Muon->fCmuFidX        = x;
    Muon->fCmuFidZ        = z;
  }
  if(fidTool.isFiducial(*trk,MuonDigiCode::CMP,d,x,z)) {
    CMPFid=1;
    if(x<0.0 && z<-3.0) CMPFid2=1;
    Muon->fCmpFidX        = x;
    Muon->fCmpFidZ        = z;
  }
  if(fidTool.isFiducial(*trk,MuonDigiCode::CMX,d,x,z)) {
    CMXFid=1;
    if(x<0.0 && z<-3.0) CMXFid2=1;
    Muon->fCmxFidX        = x;
    Muon->fCmxFidZ        = z;
  }
  if(fidTool.isFiducial(*trk,MuonDigiCode::BMU,d,x,z)) {
    BMUFid=1;
    if(x<0.0 && z<-3.0) BMUFid2=1;
    Muon->fBmuFidX        = x;
    Muon->fBmuFidZ        = z;
  }   

  // these bits are to be the same as the track block
  Muon->fFid |= ( (CMUFid << 0) |
		  (CMPFid << 1) |
		  (CMXFid << 2) |
		  (BMUFid << 3) );
  Muon->fFid |= ( (CMUFid2 << 4) |
		  (CMPFid2 << 5) |
		  (CMXFid2 << 6) |
		  (BMUFid2 << 7) );
  

  if (Ken->hasCmp() == 1) {
    if(MuonAnalysisTools::MuonIsInBluebeam(*Ken))  Muon->fFid  |= ( 1 <<  8);
  }
  if (Ken->hasCmx() == 1) {
    if(MuonAnalysisTools::MuonIsInMiniskirt(*Ken)) Muon->fFid  |= ( 1 <<  9);
    if(MuonAnalysisTools::MuonIsInKeystone(*Ken))  Muon->fFid  |= ( 1 << 10);
  }

  // now do the fiducial ignoring the keystone and miniskirt
  fidTool.setMiniskirtLive(false);
  fidTool.setKeystoneLive (false);
  if (fidTool.isFiducial(*trk,MuonDigiCode::CMX,d,x,z)) {
    CMXFidt = 1;
    if ((x < 0.0) && (z < -3.0)) CMXFid2t = 1;
    Muon->fCmxFidX2        = x;
    Muon->fCmxFidZ2        = z;
  }
  Muon->fFid |= ( (CMXFidt  << 11) |
		  (CMXFid2t << 12) );
//-----------------------------------------------------------------------------
//  calorimetry
//-----------------------------------------------------------------------------
  Muon->fEHad         = Ken->hadronEnergy();
  Muon->fEEm          = Ken->emEnergy();
  Muon->fEHadNeighb   = TStnDataBlock::kUndefined;
  Muon->fEEmNeighb    = TStnDataBlock::kUndefined;

  Muon->fConeR4Et     = Ken->coneR4Et();
  Muon->fConeR7Et     = Ken->coneR7Et();
  Muon->fConeR4EmFrac = Ken->coneR4EmFraction();
  Muon->fConeR7EmFrac = Ken->coneR7EmFraction();
  Muon->fHadTdcTime   = TStnDataBlock::kUndefined;
  Muon->fSigmaT       = TStnDataBlock::kUndefined;

  Muon->fNHitTowers   = Ken->numEmTowers();
  Muon->fNHitTowers  += (Ken->numHadTowers() << 16);
//-----------------------------------------------------------------------------
//  cached tracking info
//-----------------------------------------------------------------------------
  Muon->fZv           = trk->z0();       
  Muon->fD0           = trk->d0();
  Muon->fTIso         = TStnDataBlock::kUndefined;

  Muon->fNStubs       = 0;
  Muon->fNHits[0]     = 0;
  Muon->fNHits[1]     = 0;
  Muon->fDet          = 0;
//-----------------------------------------------------------------------------
//  CMU stub
//-----------------------------------------------------------------------------
  Muon->fCmuStubX       = TStnDataBlock::kUndefined;
  Muon->fCmuStubDxDy    = TStnDataBlock::kUndefined;
  Muon->fCmuStubZ       = TStnDataBlock::kUndefined;
  Muon->fCmuStubDzDy    = TStnDataBlock::kUndefined;

  Muon->fCmuDelX        = TStnDataBlock::kUndefined;
  Muon->fCmuDelZ        = TStnDataBlock::kUndefined;  
  Muon->fCmuDelDxDy     = TStnDataBlock::kUndefined;
  Muon->fCmuDelDzDy     = TStnDataBlock::kUndefined;
  Muon->fCmuStubChi2X   = TStnDataBlock::kUndefined;
  Muon->fCmuStubChi2Z   = TStnDataBlock::kUndefined;
  Muon->fCmuStubPos [0] = TStnDataBlock::kUndefined;
  Muon->fCmuStubPos [1] = TStnDataBlock::kUndefined;
  Muon->fCmuStubPos [2] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackPos[0] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackPos[1] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackPos[2] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackDir[0] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackDir[1] = TStnDataBlock::kUndefined;
  Muon->fCmuTrackDir[2] = TStnDataBlock::kUndefined;
  Muon->fCmuChi2Link    = TStnDataBlock::kUndefined;

  if (Ken->hasCmu() == 1) {
				// however the stub pointer can be 0
    Muon->fNStubs++;
    const MuonStub* stub = Ken->stub(CdfMuon::CMU);
    Muon->fCmuStub       = (MuonStub*) stub;

    Muon->fCmuDelX      = Ken->cmu().drphi();
    Muon->fCmuDelZ      = Ken->cmu().dz();
    Muon->fCmuDelDxDy   = Ken->cmu().dphi();
    Muon->fCmuDelDzDy   = Ken->cmu().dtheta();
    Muon->fCmuStubChi2Z = Ken->cmu().chsqZ();
//-----------------------------------------------------------------------------
// for the time being use Tom's chi2
//-----------------------------------------------------------------------------
    Muon->fCmuChi2Link  = Ken->cmu().chsqXPosition();

    if (stub) {
      MuonStub::AssHitCIter first = stub->beginAssHits();
      MuonStub::AssHitCIter last  = stub->endAssHits();

      for (MuonStub::AssHitCIter hit=first; hit!=last; hit++) {
	int lay = hit->getHit()->getCode().getLayer();
	Muon->fHitMask |= (1 << lay);
      }

      Muon->fDet          |= first->getHit()->getCode().getWedge();
      Muon->fDet          |= first->getHit()->getCode().getSide() << 5;

      Muon->fCmuStubChi2X  = stub->chiSquared();
      Muon->fCmuStubPos[0] = stub->spacePoint().x();
      Muon->fCmuStubPos[1] = stub->spacePoint().y();
      Muon->fCmuStubPos[2] = stub->spacePoint().z();

      phi = stub->spacePoint().phi();
      dr  = Ken->cmu().drphi();
					// sign can be wrong here, this has 
					// to be validated

      Muon->fCmuTrackPos[0] = stub->spacePoint().x()-dr*cos(phi);
      Muon->fCmuTrackPos[1] = stub->spacePoint().y()-dr*sin(phi);
      Muon->fCmuTrackPos[2] = stub->spacePoint().z()-Ken->cmu().dz();
      Muon->fNHits[0]      += stub->numAssHits();
    }
  }
//-----------------------------------------------------------------------------
//  CMP stub
//-----------------------------------------------------------------------------
  Muon->fCmpStubX       = TStnDataBlock::kUndefined;
  Muon->fCmpStubDxDy    = TStnDataBlock::kUndefined;
  Muon->fCmpDelX        = TStnDataBlock::kUndefined;
  Muon->fCmpDelDxDy     = TStnDataBlock::kUndefined;
  Muon->fCmpStubChi2X   = TStnDataBlock::kUndefined;
  Muon->fCmpStubPos [0] = TStnDataBlock::kUndefined;
  Muon->fCmpStubPos [1] = TStnDataBlock::kUndefined;
  Muon->fCmpStubPos [2] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackPos[0] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackPos[1] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackPos[2] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackDir[0] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackDir[1] = TStnDataBlock::kUndefined;
  Muon->fCmpTrackDir[2] = TStnDataBlock::kUndefined;
  Muon->fCmpChi2Link    = TStnDataBlock::kUndefined;

  if (Ken->hasCmp() == 1) {
    Muon->fNStubs++;
    const MuonStub* stub = Ken->stub(CdfMuon::CMP);
    Muon->fCmpStub       = (MuonStub*) stub;

    Muon->fCmpDelX       = Ken->cmp().drphi();
    Muon->fCmpDelDxDy    = Ken->cmp().dphi ();
    Muon->fCmpChi2Link   = Ken->cmp().chsqXPosition();

    if (stub) {
      MuonStub::AssHitCIter first = stub->beginAssHits();
      MuonStub::AssHitCIter last  = stub->endAssHits();

      for (MuonStub::AssHitCIter hit=first; hit!=last; hit++) {
	int lay = hit->getHit()->getCode().getLayer();

					// CMP hits - bits 4-7

	Muon->fHitMask |= (1 << (lay+4));
      }

      Muon->fDet          |= (*first).getHit()->getCode().getCMPPart()  << 8;
      Muon->fDet          |= (*first).getHit()->getCode().getRegion() << 11;

      Muon->fCmpStubChi2X  = stub->chiSquared();
      Muon->fCmpStubPos[0] = stub->spacePoint().x();
      Muon->fCmpStubPos[1] = stub->spacePoint().y();
      Muon->fCmpStubPos[2] = stub->spacePoint().z();

      phi = stub->spacePoint().phi();
      dr  = Ken->cmp().drphi();
					// sign can be wrong here, this has 
					// to be validated

      Muon->fCmpTrackPos[0] = stub->spacePoint().x()-dr*cos(phi);
      Muon->fCmpTrackPos[1] = stub->spacePoint().y()-dr*sin(phi);
      Muon->fCmpTrackPos[2] = stub->spacePoint().z()-Ken->cmp().dz();

      Muon->fNHits[0]      += (stub->numAssHits() << 8);
    }
  }
//-----------------------------------------------------------------------------
//  CSP: not defined yet
//-----------------------------------------------------------------------------
  Muon->fCspX[0]        = TStnDataBlock::kUndefined;
  Muon->fCspX[1]        = TStnDataBlock::kUndefined;
  Muon->fCspZ[0]        = TStnDataBlock::kUndefined;
  Muon->fCspZ[1]        = TStnDataBlock::kUndefined;
//-----------------------------------------------------------------------------
//  CMX stub
//-----------------------------------------------------------------------------
  Muon->fCmxStubX       = TStnDataBlock::kUndefined;
  Muon->fCmxStubZ       = TStnDataBlock::kUndefined;
  Muon->fCmxStubDxDy    = TStnDataBlock::kUndefined;
  Muon->fCmxStubDzDy    = TStnDataBlock::kUndefined;
  Muon->fCmxDelX        = TStnDataBlock::kUndefined;
  Muon->fCmxDelZ        = TStnDataBlock::kUndefined;  
  Muon->fCmxDelDxDy     = TStnDataBlock::kUndefined;
  Muon->fCmxDelDzDy     = TStnDataBlock::kUndefined;
  Muon->fCmxStubChi2X   = TStnDataBlock::kUndefined;
  Muon->fCmxStubChi2Z   = TStnDataBlock::kUndefined;
  Muon->fCmxStubPos [0] = TStnDataBlock::kUndefined;
  Muon->fCmxStubPos [1] = TStnDataBlock::kUndefined;
  Muon->fCmxStubPos [2] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackPos[0] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackPos[1] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackPos[2] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackDir[0] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackDir[1] = TStnDataBlock::kUndefined;
  Muon->fCmxTrackDir[2] = TStnDataBlock::kUndefined;
  Muon->fCmxChi2Link    = TStnDataBlock::kUndefined;
  
  if (Ken->hasCmx() == 1) {
    Muon->fNStubs++;
    const MuonStub* stub = Ken->stub(CdfMuon::CMX);
    Muon->fCmpStub       = (MuonStub*) stub;

    Muon->fCmxDelX       = Ken->cmx().drphi();
    Muon->fCmxDelZ       = Ken->cmx().dz();
    Muon->fCmxDelDxDy    = Ken->cmx().dphi();
    Muon->fCmxDelDzDy    = Ken->cmx().dtheta();

					// match ch**2 Z is not well defined,
					// stub is not fit in Z

    Muon->fCmxStubChi2Z  = Ken->cmx().chsqZ(); // TO BE FIXED
    Muon->fCmxChi2Link   = Ken->cmx().chsqXPosition();

    if (stub) {
      MuonStub::AssHitCIter first = stub->beginAssHits();
      MuonStub::AssHitCIter last  = stub->endAssHits();

      for (MuonStub::AssHitCIter hit=first; hit!=last; hit++) {
	int lay = hit->getHit()->getCode().getLayer();

					// CMX hits - bits 8-13

	Muon->fHitMask |= (1 << (lay+8));
      }
      Muon->fDet |= ((*first).getHit()->getCode().getWedge() << 16);
      Muon->fDet |= ((*first).getHit()->getCode().getSide () << 21);

      Muon->fCmxStubChi2X  = stub->chiSquared();
      Muon->fCmxStubPos[0] = stub->spacePoint().x();
      Muon->fCmxStubPos[1] = stub->spacePoint().y();
      Muon->fCmxStubPos[2] = stub->spacePoint().z();

      //Anyes add track position
      
      

      Muon->fNHits[0]      = (stub->numAssHits() << 16);
    }
  }
//-----------------------------------------------------------------------------
//  BMU stub
//-----------------------------------------------------------------------------
  Muon->fBmuStubX       = TStnDataBlock::kUndefined;
  Muon->fBmuStubDxDy    = TStnDataBlock::kUndefined;
  Muon->fBmuStubZ       = TStnDataBlock::kUndefined;
  Muon->fBmuStubDzDy    = TStnDataBlock::kUndefined;

  Muon->fBmuDelX        = TStnDataBlock::kUndefined;
  Muon->fBmuDelZ        = TStnDataBlock::kUndefined;  
  Muon->fBmuDelDxDy     = TStnDataBlock::kUndefined;
  Muon->fBmuDelDzDy     = TStnDataBlock::kUndefined;
  Muon->fBmuStubChi2X   = TStnDataBlock::kUndefined;
  Muon->fBmuStubChi2Z   = TStnDataBlock::kUndefined;
  Muon->fBmuStubPos [0] = TStnDataBlock::kUndefined;
  Muon->fBmuStubPos [1] = TStnDataBlock::kUndefined;
  Muon->fBmuStubPos [2] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackPos[0] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackPos[1] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackPos[2] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackDir[0] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackDir[1] = TStnDataBlock::kUndefined;
  Muon->fBmuTrackDir[2] = TStnDataBlock::kUndefined;
  Muon->fBmuChi2Link    = TStnDataBlock::kUndefined;
  
  if (Ken->hasBmu() == 1) {
    Muon->fNStubs++;
    const MuonStub* stub = Ken->stub(CdfMuon::BMU);
    Muon->fCmpStub       = (MuonStub*) stub;

    Muon->fBmuDelX       = Ken->bmu().drphi();
    Muon->fBmuDelZ       = Ken->bmu().dz();
    Muon->fBmuDelDxDy    = Ken->bmu().dphi();
    Muon->fBmuDelDzDy    = Ken->bmu().dtheta();

					// match ch**2 Z is not well defined,
					// stub is not fit in Z

    Muon->fBmuStubChi2Z  = Ken->bmu().chsqZ(); // TO BE FIXED
    Muon->fBmuChi2Link   = Ken->bmu().chsqXPosition();

    if (stub) {
      MuonStub::AssHitCIter first = stub->beginAssHits();
      MuonStub::AssHitCIter last  = stub->endAssHits();

      for (MuonStub::AssHitCIter hit=first; hit!=last; hit++) {
	int lay = hit->getHit()->getCode().getLayer();
					// BMU hits - bits 14-17
	Muon->fHitMask |= (1 << (lay+14));
      }

      if(first!=last){ //security against previous version which do not have hit links
	Muon->fDet |= ((*first).getHit()->getCode().getWedge() << 24);
	Muon->fDet |= ((*first).getHit()->getCode().getSide () << 29);
      }      

      Muon->fBmuStubChi2X  = stub->chiSquared();
      Muon->fBmuStubPos[0] = stub->spacePoint().x();
      Muon->fBmuStubPos[1] = stub->spacePoint().y();
      Muon->fBmuStubPos[2] = stub->spacePoint().z();
      Muon->fNHits[0]      = (stub->numAssHits() << 24);

    }
  }
//-----------------------------------------------------------------------------
// energy in the towers neighbouring to cracks yet to be defined
//-----------------------------------------------------------------------------
  Muon->fBte      = TStnDataBlock::kUndefined;
  //anyes comment out not in TStnMuon.hh
  //  Muon->fNeighbEt = Ken->neighborEt();

  Muon->fTrind    = -1;
					// muon status code:
					// we about to learn what it is

  Muon->fStat       = TStnDataBlock::kUndefined;
  Muon->fHadTdcTime = TStnDataBlock::kUndefined;
  Muon->fSigmaT     = TStnDataBlock::kUndefined;
  Muon->fCdfMuon    = (CdfMuon*) Ken;

  return 0;
}

//_____________________________________________________________________________
Int_t StntupleInitMuonBlock(TStnDataBlock* block, 
			    AbsEvent*      event, 
			    int            mode) 
{
  // fill muon data block with the event data

  int   ev_number, rn_number, inum[100], stub_number;
  char  process[100], description[100], coll_type[100];

  ev_number = AbsEnv::instance()->trigNumber();
  rn_number = AbsEnv::instance()->runNumber();

  if (! fgFolder) {
//-----------------------------------------------------------------------------
//  initialize local static variables
//-----------------------------------------------------------------------------
    fgFolder = (TFolder*) gROOT->GetRootFolder()->FindObject("StntupleFolder");
    fgLogger = (TStnErrorLogger* ) fgFolder->FindObject("StntupleErrorLogger");
  }

  if (block->Initialized(ev_number,rn_number)) return 0;

  TStnMuonBlock* data = (TStnMuonBlock*) block;
  data->Clear();

  StntupleGetProcessName(data->CollName()->Data(),
			 process,
			 description,
			 coll_type);

  EventRecord::ConstIterator ci;

  int               nmuons = 0;
  vector<CdfMuon*>  mulist;
  
  if (strcmp(coll_type,"VIEW") == 0) {
    ci = StntupleGetIterator(event,"CdfMuonView",process,description);
    if ( ci.is_valid()){
      CdfMuonView_ch handle(ci);
      const vector<ValueVectorLink<CdfMuonColl> >& kens_list = handle->contents();
      nmuons     = kens_list.size();
      for (int i=0; i<nmuons; i++) {
	const CdfMuon* ken = kens_list[i].operator->();
	mulist.push_back((CdfMuon*) ken);
      }
    }
  }else {
    ci = StntupleGetIterator(event,"CdfMuonColl",process,description);
    if ( ci.is_valid()){
      CdfMuonColl_ch handle(ci);
      const vector<CdfMuon>& kens_list = handle->contents();
      nmuons     = kens_list.size();
      for (int i=0; i<nmuons; i++) {
	mulist.push_back((CdfMuon*) &kens_list[i]);
      }
    }
  }
  
  const CdfMuon *ken,  *mu1, *mu2;;
  TStnMuon*      muon;
  // order muons in Pt  
  for (int i=0; i<nmuons; i++) inum[i] = i;
  
  for (int i1=0; i1<nmuons-1; i1++) {
    mu1 = mulist[i1];
    for (int i2=i1+1; i2<nmuons; i2++) {
      mu2 = mulist[i2];
      if (mu1->bestTrack()->pt() < mu2->bestTrack()->pt()) {

					// swap the muons
	int itmp = inum[i1];
	inum[i1] = inum[i2];
	inum[i2] = itmp;
	mu1      = mu2;
      }
    }
  }

  for (int i=0; i<nmuons; i++) {
    ken = mulist[inum[i]];
    muon  = data->NewMuon();
    StntupleInitMuon(muon,ken);
    
    if(muon->NStubs()==0) data->fNCmio++;
  }
//-----------------------------------------------------------------------------
// process CMU stubs
//-----------------------------------------------------------------------------
  TBaseMuonStub* stub;
  const MuonHit* muon_hit;

  ConstEventIter cmu_stub_event_iter = StntupleGetIterator(event,
							   "MuonStubColl",
							   process,
							   "CMUStubs");

  if (cmu_stub_event_iter.is_valid()) {
    ConstHandle<MuonStubColl> cmu_coll(cmu_stub_event_iter);
    const MuonStubColl::container_type& cmu = cmu_coll->contents();
    for (MuonStubColl::container_type::const_iterator cmu_iter=cmu.begin(); 
                                        cmu_iter != cmu.end(); cmu_iter++) {

      muon_hit = cmu_iter->beginAssHits()->getHit();

      stub = data->NewCmuStub();
      stub_number = data->NCmuStubs()-1;

      stub->SetNHits(cmu_iter->numAssHits());
      stub->SetNcHits(cmu_iter->numCounterHits());
      stub->SetChi2XY(cmu_iter->chiSquaredX());
      stub->SetChi2XZ(cmu_iter->chiSquaredZ());
      stub->SetPhi(cmu_iter->spacePoint().phi());
      stub->SetTheta(cmu_iter->spacePoint().theta());
      if (muon_hit) {
	stub->SetDigiCode(muon_hit->getCode().getKey());
      }
      else {
	stub->SetDigiCode(0);
      }
				        // muon corresponding to this stub
      for (int i=0; i<nmuons; i++) {
	TStnMuon* muon = data->Muon(i);
	MuonStub* cmus = muon->GetCmuStub();
	if (cmus == &(*cmu_iter)) {
	  muon->SetCmuStubNumber(stub_number);
	}
      }
    }
  }
//-----------------------------------------------------------------------------
// process CMP stubs
//-----------------------------------------------------------------------------
  ConstEventIter cmp_stub_event_iter = StntupleGetIterator(event,
							   "MuonStubColl",
							   process,
							   "CMPStubs");
  if (cmp_stub_event_iter.is_valid()) {
    ConstHandle<MuonStubColl> cmp_coll(cmp_stub_event_iter);
    const MuonStubColl::container_type& cmp = cmp_coll->contents();
    for (MuonStubColl::container_type::const_iterator cmp_iter=cmp.begin(); 
                                        cmp_iter != cmp.end(); cmp_iter++) {

      muon_hit    = cmp_iter->beginAssHits()->getHit();
      stub        = data->NewCmpStub();
      stub_number = data->NCmpStubs()-1;

      stub->SetNHits(cmp_iter->numAssHits());
      stub->SetNcHits(cmp_iter->numCounterHits());
      stub->SetChi2XY(cmp_iter->chiSquaredX());
      stub->SetPhi(cmp_iter->spacePoint().phi());
      stub->SetTheta(cmp_iter->spacePoint().theta());
      if (muon_hit) {
	stub->SetDigiCode(muon_hit->getCode().getKey());
      }
      else {
	stub->SetDigiCode(0);
      }
				        // muon corresponding to this stub
      for (int i=0; i<nmuons; i++) {
	TStnMuon* muon = data->Muon(i);
	MuonStub* cmps = muon->GetCmpStub();
	if (cmps == &(*cmp_iter)) {
	  muon->SetCmpStubNumber(stub_number);
	}
      }
    }
  }
//-----------------------------------------------------------------------------
// process CMX stubs
//-----------------------------------------------------------------------------
  ConstEventIter cmx_stub_iter = StntupleGetIterator(event,
						     "MuonStubColl",
						     process,
						     "CMXStubs");
  if (cmx_stub_iter.is_valid()) {
    ConstHandle<MuonStubColl> cmx_coll(cmx_stub_iter);
    const MuonStubColl::container_type& cmx = cmx_coll->contents();
    for (MuonStubColl::container_type::const_iterator dave=cmx.begin(); 
	                                     dave != cmx.end(); dave++) {
      stub        = data->NewCmxStub();
      stub_number = data->NCmxStubs()-1;

      stub->SetNHits(dave->numAssHits());
      stub->SetNcHits(dave->numCounterHits());
      stub->SetChi2XY(dave->chiSquaredX());
      stub->SetPhi(dave->spacePoint().phi());
      stub->SetTheta(dave->spacePoint().theta());
      const MuonHit* daves_hit = dave->beginAssHits()->getHit();
      if (daves_hit) {
	stub->SetDigiCode(daves_hit->getCode().getKey());
      }
      else {
	stub->SetDigiCode(0);
      }
				        // muon corresponding to this stub
      for (int i=0; i<nmuons; i++) {
	TStnMuon* muon = data->Muon(i);
	MuonStub* cmxs = muon->GetCmxStub();
	if (cmxs == &(*dave)) {
	  muon->SetCmxStubNumber(stub_number);
	}
      }
    }
  }

//-----------------------------------------------------------------------------
// process BMU stubs
//-----------------------------------------------------------------------------

  ConstEventIter bmu_stub_iter = StntupleGetIterator(event,
						     "MuonStubColl",
						     process,
						     "BMUStubs");
  if (bmu_stub_iter.is_valid()) {

    ConstHandle<MuonStubColl> bmu_coll(bmu_stub_iter);
    const MuonStubColl::container_type& bmu = bmu_coll->contents();
    for (MuonStubColl::container_type::const_iterator istub=bmu.begin(); 
	 istub != bmu.end(); istub++) {

      stub        = data->NewBmuStub();
      stub_number = data->NBmuStubs()-1;

      stub->SetNHits (istub->numAssHits()        );
      stub->SetNcHits(istub->numCounterHits()    );
      stub->SetChi2XY(istub->chiSquaredX()       );
      stub->SetChi2XZ(istub->chiSquaredZ()       );
      stub->SetPhi   (istub->spacePoint().phi()  );
      stub->SetTheta (istub->spacePoint().theta());
      // this protection is needed against BMU stub bug in 6.1.1
      // where iterators are empty for 2-hit stub occasionally
      if (istub->beginAssHits() != istub->endAssHits()) {
	muon_hit = istub->beginAssHits()->getHit();
	if(muon_hit) {
	  stub->SetDigiCode(muon_hit->getCode().getKey());
	}
	else {
	  stub->SetDigiCode(0);
	}
      }
				        // muon corresponding to this stub
      for (int i=0; i<nmuons; i++) {
	TStnMuon* muon = data->Muon(i);
	MuonStub* bmus = muon->GetBmuStub();
	if (bmus == &(*istub)) {
	  muon->SetBmuStubNumber(stub_number);
	}
      }
    }
  }

					// set event and run numbers to 
					// mark block as initialized
  data->f_RunNumber   = rn_number;
  data->f_EventNumber = ev_number;

  return 0;
}


//_____________________________________________________________________________
Int_t StntupleMuonBlockLinks(TStnDataBlock* block,AbsEvent* event,int mode) 
{
  // set links between the electrons and the objects in the other blocks
  // return 0 if everything is OK

  int                ev_number, rn_number;
  TStnTrackBlock*    track_block;
  TStnMuonBlock*     muo_block = (TStnMuonBlock*) block;

  TStnMuon*     mu;
  TStnTrack*    trk;
  CdfMuon*      ken;
  const CdfTrack*     best_track;
  const CdfTrack*     next_track;
  Int_t         nmu, it, rc;

  ev_number = AbsEnv::instance()->trigNumber();
  rn_number = AbsEnv::instance()->runNumber();

  if (! block->Initialized(ev_number,rn_number)) return -1;

					// do not do initialize links 2nd time

  if (block->LinksInitialized()) return 0;

					// need pointer to the default track 
					// block

  //-----------------------------------------------------------------------------
  // fill tower link block
  //-----------------------------------------------------------------------------

  TStnEvent* ev = muo_block->GetEvent();
  track_block   = (TStnTrackBlock*) ev->GetDataBlock("TrackBlock");

  TObjArray* hptl = ev->GetListOfHptl();

  rc = 0;
  nmu = muo_block->NMuons();
  for (int i=0; i<nmu; i++) {
    mu = muo_block->Muon(i);
    ken = mu->GetCdfMuon();

    int tower_key;
    for (CdfMuon::TowerKeySet::const_iterator itow = ken->emTowers().begin();
	 itow != ken->emTowers().end(); ++itow) {
      tower_key = (((itow->iPhi() & 0x3f  )      ) |
		   ((itow->iEta() & 0x3f  ) << 8 ) |
		   ((itow->isValid()) << 16) );
      muo_block->EmTowerLinkList()->Add(i,tower_key);
    }

    int    nhits = 0;
    int    nt = 0;
    double t0 = 0;
    double t2 = 0;

    CalData_ch cal_data;
    CalData::Error result = CalData::find(cal_data);

    if (result == CalData::OK) {
      for (CdfMuon::TowerKeySet::const_iterator itow = ken->hadTowers().begin();
	   itow != ken->hadTowers().end(); ++itow) {

	const CalTower* tower = NULL;
	if(itow->isValid()) tower = cal_data->tower(*itow);

	tower_key = (((itow->iPhi() & 0x3f  )      ) |
		     ((itow->iEta() & 0x3f  ) << 8 ) |
		     ((itow->isValid()) << 16) );
	muo_block->HadTowerLinkList()->Add(i,tower_key);

	if (tower) {
	  CellKey key(itow->iEta(),itow->iPhi(),CHA);
	  for (int i=0; i<1; i++) {
	    nhits = tower->nTDCHits(key);
	    if ((nhits > 0) && (nhits < 3)) {
	      t0 += tower->time(key);
	      t2 += t0*t0;
	      nt++;
	    }
	  }
	}
      } // end tower loop
//-----------------------------------------------------------------------------
//  set HAD TDC times (protect against round-off errors)
//-----------------------------------------------------------------------------
      if (nt > 0) {
	t0 = t0/nt;
	mu->SetHadTdcTime(t0);
	mu->SetSigmaT    (sqrt(t2/nt-t0*t0+1.e-5));
      }
    }      // end CalData OK
//-----------------------------------------------------------------------------
//  link to tracks (best one and 2nd best one)
//-----------------------------------------------------------------------------
    best_track = ken->bestTrack().operator->();
    it  = track_block->TrackNumber(best_track);
    mu->SetTrackNumber(it);
    trk = track_block->Track(it);
    mu->SetTrack(trk);
    if (it == -1) {
//-----------------------------------------------------------------------------
//  best track should always exist !
//-----------------------------------------------------------------------------
      fgLogger->Report(-112, Form("StntupleMuonBlockLinks: no track for muon %i",
				     i));
      mu->SetGenpNumber(-1);
      mu->SetObspNumber(-1);
    } 
    else {
      mu->SetTIso(trk->Iso4());
      mu->SetBcPt(fabs(TStntuple::TrackPt(trk->BcC0())));
      mu->SetGenpNumber( trk->GenpNumber() );
      mu->SetObspNumber( trk->ObspNumber() );
    }
					// process 2nd best track

    next_track = ken->secondBestTrack().operator->();
    if (next_track == NULL) it = -1;
    else            it  = track_block->TrackNumber(next_track);
    mu->SetSecondBestTrackNumber(it);
//-----------------------------------------------------------------------------
//  update list of good HPT leptons to calculate MET (if any...) - 
//  we do have isolation at this point
//-----------------------------------------------------------------------------
    if (TStntuple::GoodHptMuon(mu)) {
      hptl->Add(mu);
    }
  }
//-----------------------------------------------------------------------------
// mark links as initialized
//-----------------------------------------------------------------------------
  block->SetLinksInitialized();

  return rc;
}