/***************************************************************************
 *
 * $Id: StHiMicroMaker.cxx,v 1.4 2002/05/31 21:50:14 jklay Exp $                                      
 *
 * Author: Bum Choi, UT Austin, Apr 2002
 *
 ***************************************************************************
 *
 * Description:  This Maker creates the highpt uDST's from StEvent
 *               for highpt Analysis.
 *
 ***************************************************************************
 *
 * $Log: StHiMicroMaker.cxx,v $
 * Revision 1.4  2002/05/31 21:50:14  jklay
 * Fixed the way centrality is calculated, see README
 *
 * Revision 1.3  2002/04/03 00:37:41  jklay
 * Fixed some bugs, added new version of dcaz
 *
 * Revision 1.2  2002/04/02 23:35:14  jklay
 * Added L3RichTrigger information
 *
 * Revision 1.1  2002/04/02 20:00:41  jklay
 * Bums highpt uDST Maker
 *
 *
 **************************************************************************/
#include <iostream>

using namespace std;

#include "StHiMicroMaker.h"

#include "TMath.h"
#include "TRandom.h"
#include "TChain.h"

#include "../StHiMicroEvent/StHiMicroEvent.h"
#include "StCtbMatcher.hh"
#include "Helper.h"

#include "StIOMaker/StIOMaker.h"
#include "StThreeVectorF.hh"
#include "StThreeVectorD.hh"
#include "StEventTypes.h"
#include "StMessMgr.h"
#include "SystemOfUnits.h" // tesla...
#include "StTpcDedxPidAlgorithm.h"
#include "StuProbabilityPidAlgorithm.h"
#include "StPhysicalHelixD.hh"
#include "StuRefMult.hh"
#include "StTrackDetectorInfo.h"
#include "tables/St_vertexSeed_Table.h" 

#include "StMuDSTMaker/COMMON/StMuDstMaker.h"
#include "StMuDSTMaker/COMMON/StMuDst.h"
#include "StMuDSTMaker/COMMON/StMuEvent.h"
#include "StMuDSTMaker/COMMON/StMuTrack.h"

#define PR(x) cout << "##### StHiMicroMaker: " << (#x) << " = " << (x) << endl;


// increasing order
inline bool sortCmp(StTpcHit* p1, StTpcHit* p2){
  return p1->position().perp()<p2->position().perp();
}


//    
 //  The constructor. Initialize you data members here.
//    
StHiMicroMaker::StHiMicroMaker(const Char_t *makerName)
  : StMaker(makerName), mDSTFile(0), mDSTTree(0), mIOMaker(0), mMudstMaker(0),
    mNEvent(0), mNAcceptedEvent(0),
    mNAcceptedTrack(0),
    mHitLoop(0), mSaveLowPt(0)
{

}

//
//  Usually ok to leave this as it is.
//
StHiMicroMaker::~StHiMicroMaker() 
{ 
/* noop */ 
}

//
//  Called once at the beginning.
//  This is a good place to book histos and tuples.
//

Int_t StHiMicroMaker::InitRun(int runID) {
  if (mDebug) gMessMgr->Info() << "StHiMicroMaker: InitRun()" << endm;
  PR(runID);

  return kStOK;
}
//---------------------------
Int_t
StHiMicroMaker::Init()
{

    cout << "###StHiMicroMaker::Init():\n";

    const char* set = (mDebug) ? "ON" : "OFF";
    cout << "\tDebug is " << set << endl;

    //Pretty much moved everything that was here to InitRun(int)
    // Moved everything back to be able to read more tahn one run (Mvl jun 04)

    if(mSaveAllEvents) { cout << "\t<I> Saving event info for events without a primary vertex!!" << endl; }

    if(mDebug) cout << "##Creating StHiMicroEvent..." << endl;
    mHiMicroEvent = new StHiMicroEvent;

    mIOMaker = (StIOMaker*)GetMaker("IO");    
    mMudstMaker = (StMuDstMaker*)GetMaker("MuDst");    
    /* Probably won't work, now in Make()
    if(mIOMaker) mInFileName = strrchr(mIOMaker->GetFile(),'/')+1;
    Int_t stat = openFile();
    */

    // beam line from jon gans
    getBeamLine();

    return StMaker::Init();
}

//
//  Called every event after Make(). Usually you do not
//  need to do anything here. Leave it as it is.
//
void
StHiMicroMaker::Clear(Option_t *opt)
{
    StMaker::Clear();
}

//
//  Called once at the end.
//
Int_t
StHiMicroMaker::Finish()
{

  cout << "###StHiMicroMaker::Finish()\n";
  cout << "\tTotal " << mNEvent << " events." << endl;
  cout << "\tProcessed " << mNAcceptedEvent << " events." << endl;
  cout << "\tTracks : " << mNAcceptedTrack << endl << endl;
  closeFile();

  return StMaker::Finish();
}

//
//  This method is called every event. That's the
//  right place to plug in your analysis. 
//
Int_t
StHiMicroMaker::Make()
{
  mNEvent++; 
	
  //
  // it it's a new file, close the old one and open a new one
  //
  TString curFileName;
  if(mIOMaker) curFileName = strrchr(mIOMaker->GetFile(),'/')+1;
  if(mMudstMaker) curFileName = strrchr(mMudstMaker->chain()->GetFile()->GetName(),'/')+1;
  
  if(mInFileName!=curFileName){
    if(mDebug) {
      cout << "##New file found : " << curFileName << endl;
      cout << "##Replacing " << mInFileName << endl;
    }
    closeFile();
    mInFileName = curFileName;
    openFile();
  }

  //
  //   First try MuDst
  //
  StMuDst* mudst=0;
  mudst = (StMuDst *) GetInputDS("MuDst");
  if (mudst) {
    processMuDst(mudst);
    return kStOk;
  }
    
  //
  //  If not, use StEvent
  //
  StEvent* event=0;
  event = (StEvent *) GetInputDS("StEvent");

  if (event) 
    processStEvent(event);
  
  return kStOK;        // if no event, we're done
}

void StHiMicroMaker::processStEvent(StEvent *event) {
  // always require trigger
  if(!event->triggerDetectorCollection()) return;

  //
  if (mSaveAllEvents || accept(event)) {   

    mNAcceptedEvent++;  
    // 
    // fill tracks and hits
    //
    //30.May.2002 nGoodEta is no longer used, as there is a common defintion
    //from StEvent now...
    fillTracks(event);

    //Save fpd info
    fillFpdCollection(event);
    
    computeQuadrant(event);

    //
    // fill StHiMicroEvent
    //
    //30.May.2002 nGoodEta is no longer used, as there is a common defintion
    //from StEvent now...
    fillEvent(event);

    //
    // fill the tree
    //
    mDSTTree->Fill();

  }
  /*
  else {

    if(mSaveAllEvents) {  //If we really want them all...
      fillEvent(event);
      mDSTTree->Fill();
    }
  }
  */
  mHiMicroEvent->Clear(); // clears all the clonesarrays
}
//____________________

void StHiMicroMaker::processMuDst(StMuDst *mudst) {
  StMuEvent *event=mudst->event();
  if (!event)
    return;

  //
  if (mSaveAllEvents || accept(event)) {   

    mNAcceptedEvent++;  
    // 
    // fill tracks and hits
    //
    //30.May.2002 nGoodEta is no longer used, as there is a common defintion
    //from StEvent now...
    fillTracks(mudst);

    // Save fpd info // Not for now...
    // fillFpdCollection(event);
    
    computeQuadrant(mudst);

    //
    // fill StHiMicroEvent
    //
    //30.May.2002 nGoodEta is no longer used, as there is a common defintion
    //from StEvent now...
    fillEvent(mudst);

    //
    // fill the tree
    //
    if (mDebug)
      cout << "mDSTTree at : " << mDSTTree << endl;
    mDSTTree->Fill();

  }
  /*
  else {

    if(mSaveAllEvents) {  //If we really want them all...
      fillEvent(event);
      mDSTTree->Fill();
    }
  }
  */
  mHiMicroEvent->Clear(); // clears all the clonesarrays
}

void StHiMicroMaker::fillCtbZdcBbc(StCtbTriggerDetector &CTB, StZdcTriggerDetector &ZDC, StBbcTriggerDetector &BBC) {
  Float_t ctb  = -1., zdce = -1, zdcw = -1, zdcVertexZ = 999.;
  // get CTB
  for (UInt_t slat=0; slat<CTB.numberOfSlats(); slat++) {
    for (UInt_t tray=0; tray<CTB.numberOfTrays();tray++) {
      ctb += CTB.mips(tray,slat,0);
    }
  }
  //get ZDCe and ZDCw        
  zdce = ZDC.adcSum(east);
  zdcw = ZDC.adcSum(west);
  zdcVertexZ = ZDC.vertexZ();

  // BBC - 01/27/03 BUM
  mHiMicroEvent->SetBBCNHitE(BBC.nHitEast()); 
  mHiMicroEvent->SetBBCNHitW(BBC.nHitWest());
  mHiMicroEvent->SetBBCNHitEL(BBC.nHitEastLarge()); 
  mHiMicroEvent->SetBBCNHitWL(BBC.nHitWestLarge());
  mHiMicroEvent->SetBBCAdcSumE(BBC.adcSumEast()); 
  mHiMicroEvent->SetBBCAdcSumW(BBC.adcSumWest());
  mHiMicroEvent->SetBBCAdcSumEL(BBC.adcSumEastLarge()); 
  mHiMicroEvent->SetBBCAdcSumWL(BBC.adcSumWestLarge());
  mHiMicroEvent->SetBBCTdcEarliestE(BBC.tdcEarliestEast()); 
  mHiMicroEvent->SetBBCTdcEarliestW(BBC.tdcEarliestWest()); 
  mHiMicroEvent->SetBBCVertexZ(BBC.zVertex());
   
  
  mHiMicroEvent->SetCTB(ctb);
  mHiMicroEvent->SetZDCe(zdce);
  mHiMicroEvent->SetZDCw(zdcw);
  mHiMicroEvent->SetZDCVertexZ(zdcVertexZ);
}

void StHiMicroMaker::fillQuadrant() {
  // 03/02/03 - quadrant stuff

  mHiMicroEvent->SetNchQLead(mQuad[0].QLead);
  mHiMicroEvent->SetNchQOpp(mQuad[0].QOpp);
  mHiMicroEvent->SetNchQA(mQuad[0].QA);
  mHiMicroEvent->SetNchQB(mQuad[0].QB);

  if (mQuad[0].LeadMom.mag() > 0.001) { 
    mHiMicroEvent->SetLeadPt(mQuad[0].LeadMom.perp());
    mHiMicroEvent->SetLeadPhi(mQuad[0].LeadMom.phi());
    mHiMicroEvent->SetLeadEta(mQuad[0].LeadMom.pseudoRapidity());
  }
  else {
    mHiMicroEvent->SetLeadPt(0);
    mHiMicroEvent->SetLeadPhi(-999);
    mHiMicroEvent->SetLeadEta(-999);
  }

  mHiMicroEvent->SetFTPCeQLead(mQuad[1].QLead);
  mHiMicroEvent->SetFTPCeQOpp(mQuad[1].QOpp);
  mHiMicroEvent->SetFTPCeQA(mQuad[1].QA);
  mHiMicroEvent->SetFTPCeQB(mQuad[1].QB);

  mHiMicroEvent->SetFTPCwQLead(mQuad[2].QLead);
  mHiMicroEvent->SetFTPCwQOpp(mQuad[2].QOpp);
  mHiMicroEvent->SetFTPCwQA(mQuad[2].QA);
  mHiMicroEvent->SetFTPCwQB(mQuad[2].QB);
}  

void
StHiMicroMaker::fillEvent(StEvent* stEvent)
{
  

  //Set it to some really big value, so it will be cut out by relevant analyses
//  const StThreeVectorF& prVtxPos(999.,999.,999.);

  //if (stEvent->primaryVertex()) {
  if (stEvent->primaryVertex() &&
     (fabs(stEvent->primaryVertex()->position().x())>0.001 || 
      fabs(stEvent->primaryVertex()->position().y())>0.001 ||
      fabs(stEvent->primaryVertex()->position().z())>0.001)) { 
    mHiMicroEvent->SetVertexZ(stEvent->primaryVertex()->position().z());
    mHiMicroEvent->SetVertexX(stEvent->primaryVertex()->position().x());
    mHiMicroEvent->SetVertexY(stEvent->primaryVertex()->position().y());
    mHiMicroEvent->SetOriginMult(stEvent->primaryVertex()->numberOfDaughters()); 
   
  } else {
    mHiMicroEvent->SetVertexZ(999.);
    mHiMicroEvent->SetVertexX(999.);
    mHiMicroEvent->SetVertexY(999.);
    mHiMicroEvent->SetOriginMult(0); 
  }
  
  // beam stuff
  if(stEvent->runInfo()){
    StRunInfo* runInfo=stEvent->runInfo();
    mHiMicroEvent->SetCenterOfMassEnergy(runInfo->centerOfMassEnergy());
    mHiMicroEvent->SetMagneticField(runInfo->magneticField());
    mHiMicroEvent->SetBeamMassNumberEast(runInfo->beamMassNumber(east));
    mHiMicroEvent->SetBeamMassNumberWest(runInfo->beamMassNumber(west));

    mHiMicroEvent->SetBackgroundRate(runInfo->backgroundRate());
    mHiMicroEvent->SetZDCeRate(runInfo->zdcEastRate());
    mHiMicroEvent->SetZDCwRate(runInfo->zdcWestRate());
    mHiMicroEvent->SetZDCCoincidenceRate(runInfo->zdcCoincidenceRate());
  }
  else{
    gMessMgr->Info() << "StHiMicroMaker: no Run Info, reverting to year 1 settings "
                     << endm;
    mHiMicroEvent->SetCenterOfMassEnergy(130);
    mHiMicroEvent->SetMagneticField(4.98);
    mHiMicroEvent->SetBeamMassNumberEast(197);
    mHiMicroEvent->SetBeamMassNumberWest(197);

    mHiMicroEvent->SetBackgroundRate(-1);
    mHiMicroEvent->SetZDCeRate(-1);
    mHiMicroEvent->SetZDCwRate(-1);
    mHiMicroEvent->SetZDCCoincidenceRate(-1);

  }

  mHiMicroEvent->SetNUncorrectedNegativePrimaries(uncorrectedNumberOfNegativePrimaries(*stEvent));
  mHiMicroEvent->SetNUncorrectedPrimaries(uncorrectedNumberOfPrimaries(*stEvent));

  //30.May.2002 - fixed this definition to match the rest of STAR
  //same as mNuncorrectedNumberOfPrimaries
  //  mHiMicroEvent->SetCentMult(mHiMicroEvent->NUncorrectedPrimaries());
  mHiMicroEvent->SetCentrality(mHiMicroEvent->NUncorrectedPrimaries());

  mHiMicroEvent->SetRunId((Int_t) stEvent->runId());
  mHiMicroEvent->SetEventId((Int_t) stEvent->id());

//Need CTB Match stuff for good global counts
  // jons CTB stuff
  
  vector<ctbHit> *ctbHits=new vector<ctbHit>;
  StCtbTriggerDetector &ctbDet = stEvent->triggerDetectorCollection()->ctb();

  for (unsigned int slat = 0; slat < ctbDet.numberOfSlats(); slat++) {
    for (unsigned int tray = 0; tray < ctbDet.numberOfTrays(); tray++) {
      ctbHit curHit;
      curHit.adc = ctbDet.mips(tray,slat,0);
      if(curHit.adc > 0){
	ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
	ctbHits->push_back(curHit);
      }
    }
  }
  
//Get the number of Global tracks.  Not sure if this is the best way, but it will work
//JLK

  StTrack *track;
  StSPtrVecTrackNode& theNodes = stEvent->trackNodes();
  long allcnt = 0; long goodcnt = 0; long goodFlagcnt = 0; long goodCTBcnt = 0;

  for (unsigned int i=0; i<theNodes.size(); i++) {
     allcnt += theNodes[i]->entries(global);   
     track = theNodes[i]->track(global);
//goodGlobal checks globals against evr criteria for vertex finding
     bool shouldHitCTB=false;
     double etaInCTBFrame(0),phiInCTBFrame(0);    
     unsigned int matchToCTB = 
       EtaAndPhiToOrriginAtCTB(track->geometry()->helix(),ctbHits,shouldHitCTB,etaInCTBFrame,phiInCTBFrame);

     if (goodGlobal(track)) goodcnt++;
     if (goodGlobalCTBMatch(track,matchToCTB)) goodCTBcnt++;
     if (goodGlobalFlag(track)) goodFlagcnt++;
  }
  mHiMicroEvent->SetNAllGlobals(allcnt);
  mHiMicroEvent->SetNGoodGlobals(goodcnt);
  mHiMicroEvent->SetNFlagGlobals(goodFlagcnt);
  mHiMicroEvent->SetNGoodGlobalsCTBMatch(goodCTBcnt);

  // jon
  if(ctbHits) delete ctbHits;

  // Trigger stuff, the old-fashioned way (for 2001 data)
  //** 01/28/02 - trigger word added
  StL0Trigger* pTrigger = stEvent->l0Trigger();
  if(pTrigger){
    mHiMicroEvent->SetL0TriggerWord(pTrigger->triggerWord());
  }
  StL3Trigger* pL3Trigger = stEvent->l3Trigger();
  if(pL3Trigger){
    mHiMicroEvent->SetL3UnbiasedTrigger(pL3Trigger->l3EventSummary()->unbiasedTrigger());
    mHiMicroEvent->SetL3RichTrigger(false);

    //Now check for L3 Rich Trigger
    const StPtrVecL3AlgorithmInfo& algoInfo = pL3Trigger->l3EventSummary()->algorithmsAcceptingEvent();
    for (UInt_t i = 0; i < algoInfo.size(); i++) {
       if(algoInfo[i]->id() == 4) mHiMicroEvent->SetL3RichTrigger(true);
    }
  }

  //
  // 02/21/03 - BUM -trigger id stuff
  //
  // default to 0
  mHiMicroEvent->ClearTriggerIds();

  StTriggerIdCollection* pId=stEvent->triggerIdCollection();
  if(pId && pId->nominal()){
    
    vector<unsigned int> vv = pId->nominal()->triggerIds();
    for(unsigned int i=0;i<vv.size();i++){
      //cout << i << " : " << vv[i] << endl;
      mHiMicroEvent->AddTriggerId(i,vv[i]);
    }

    // this is a reality check
    if( (mHiMicroEvent->IsTrigger(2001) || mHiMicroEvent->IsTrigger(2003)) 
	!= (pId->nominal()->isTrigger(2001) ||pId->nominal()->isTrigger(2003))){
	  cout << ">>>> WHATS UP WITH THE TRIGGER?" << endl;
    }

  }

  // reality check
  //
  /*
  cout << "-----" << endl;
  for(int i=0;i<32;i++){
    if(mHiMicroEvent->TriggerId(i))
       cout << i << " : " << mHiMicroEvent->TriggerId(i) << endl;
  }
  */

  // taken from flow maker
  //

  StTriggerDetectorCollection *triggers 
    = stEvent->triggerDetectorCollection();
  // 03/15/03 - now require triggers!=0 at Make()
  //            so this should always be true.
  if (triggers) {
    StCtbTriggerDetector &CTB = triggers->ctb();
    StZdcTriggerDetector &ZDC = triggers->zdc();
    StBbcTriggerDetector &BBC = triggers->bbc();
    fillCtbZdcBbc(CTB,ZDC,BBC);
  }
  else
    cout << "WARNING: No trigger detectors found in " << __PRETTY_FUNCTION__ 
	 << " will cause problems without proper initialisation." <<  endl;
  
  // count FTPC
  //if(stEvent->primaryVertex()){
  if(stEvent->primaryVertex() &&
     (fabs(stEvent->primaryVertex()->position().x())>0.001 || 
      fabs(stEvent->primaryVertex()->position().y())>0.001 ||
      fabs(stEvent->primaryVertex()->position().z())>0.001)) { 
    mHiMicroEvent->SetFTPCe(mFTPCe); mHiMicroEvent->SetFTPCw(mFTPCw);
  }
  else{
    mHiMicroEvent->SetFTPCe(-1); mHiMicroEvent->SetFTPCw(-1);
  }

  fillQuadrant();

  if (mDebug) {
    cout << "##StHiMicroMaker::fillEvent()" << endl;
    cout << *mHiMicroEvent << endl; 
  }
}


void StHiMicroMaker::fillEvent(StMuDst* mudst) {

  StMuEvent *event=mudst->event();
  
  if (fabs(event->primaryVertexPosition().x())>0.001 || 
      fabs(event->primaryVertexPosition().y())>0.001 ||
      fabs(event->primaryVertexPosition().z())>0.001) { 
    mHiMicroEvent->SetVertexZ(event->primaryVertexPosition().z());
    mHiMicroEvent->SetVertexX(event->primaryVertexPosition().x());
    mHiMicroEvent->SetVertexY(event->primaryVertexPosition().y());
    mHiMicroEvent->SetOriginMult(0); // Don't set for the moment, we do not have a raw track count in MuEvent 
   
  } else {
    mHiMicroEvent->SetVertexZ(999.);
    mHiMicroEvent->SetVertexX(999.);
    mHiMicroEvent->SetVertexY(999.);
    mHiMicroEvent->SetOriginMult(0); 
  }

  // beam stuff
  StRunInfo &runInfo=event->runInfo();
  mHiMicroEvent->SetCenterOfMassEnergy(runInfo.centerOfMassEnergy());
  mHiMicroEvent->SetMagneticField(runInfo.magneticField());
  mHiMicroEvent->SetBeamMassNumberEast(runInfo.beamMassNumber(east));
  mHiMicroEvent->SetBeamMassNumberWest(runInfo.beamMassNumber(west));

  mHiMicroEvent->SetBackgroundRate(runInfo.backgroundRate());
  mHiMicroEvent->SetZDCeRate(runInfo.zdcEastRate());
  mHiMicroEvent->SetZDCwRate(runInfo.zdcWestRate());
  mHiMicroEvent->SetZDCCoincidenceRate(runInfo.zdcCoincidenceRate());

  mHiMicroEvent->SetNUncorrectedNegativePrimaries(event->refMultNeg());
  //mHiMicroEvent->SetNUncorrectedPositivePrimaries(event->refMultPos());
  mHiMicroEvent->SetNUncorrectedPrimaries(event->refMult());

  //30.May.2002 - fixed this definition to match the rest of STAR
  //same as mNuncorrectedNumberOfPrimaries
  //  mHiMicroEvent->SetCentMult(mHiMicroEvent->NUncorrectedPrimaries());
  mHiMicroEvent->SetCentrality(event->refMult());

  mHiMicroEvent->SetRunId((Int_t) event->runId());
  mHiMicroEvent->SetEventId((Int_t) event->eventId());

//Need CTB Match stuff for good global counts
  // jons CTB stuff
  
  vector<ctbHit> *ctbHits=new vector<ctbHit>;
  StCtbTriggerDetector &ctbDet = event->ctbTriggerDetector();

  for (unsigned int slat = 0; slat < ctbDet.numberOfSlats(); slat++) {
    for (unsigned int tray = 0; tray < ctbDet.numberOfTrays(); tray++) {
      ctbHit curHit;
      curHit.adc = ctbDet.mips(tray,slat,0);
      if(curHit.adc > 0){
	ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
	ctbHits->push_back(curHit);
      }
    }
  }
  
//Get the number of Global tracks.  Not sure if this is the best way, but it will work
//JLK

  StMuTrack *track;
  long goodcnt = 0; long goodFlagcnt = 0; long goodCTBcnt = 0;

  long allcnt=mudst->numberOfGlobalTracks();
  for (unsigned int i=0; i<(unsigned int)allcnt; i++) {
     track = mudst->globalTracks(i);
//goodGlobal checks globals against evr criteria for vertex finding
     bool shouldHitCTB=false;
     double etaInCTBFrame(0),phiInCTBFrame(0);    
     unsigned int matchToCTB = 
       EtaAndPhiToOrriginAtCTB(track->helix(),ctbHits,shouldHitCTB,etaInCTBFrame,phiInCTBFrame);

     if (goodGlobal(track)) goodcnt++;
     if (goodGlobalCTBMatch(track,matchToCTB)) goodCTBcnt++;
     if (goodGlobalFlag(track)) goodFlagcnt++;
  }
  mHiMicroEvent->SetNAllGlobals(allcnt);
  mHiMicroEvent->SetNGoodGlobals(goodcnt);
  mHiMicroEvent->SetNFlagGlobals(goodFlagcnt);
  mHiMicroEvent->SetNGoodGlobalsCTBMatch(goodCTBcnt);

  // jon
  if(ctbHits) delete ctbHits;

  // Trigger stuff, the old-fashioned way (for 2001 data)
  //** 01/28/02 - trigger word added
  StL0Trigger &pTrigger = event->l0Trigger();
  mHiMicroEvent->SetL0TriggerWord(pTrigger.triggerWord());

  mHiMicroEvent->SetL3UnbiasedTrigger(event->l3EventSummary().unbiasedTrigger());
  mHiMicroEvent->SetL3RichTrigger(false);

  // Cannot check for L3 Rich Trigger on MuDst (?) 
  /*
  const StPtrVecL3AlgorithmInfo& algoInfo = pL3Trigger->l3EventSummary()->algorithmsAcceptingEvent();
    for (UInt_t i = 0; i < algoInfo.size(); i++) {
       if(algoInfo[i]->id() == 4) mHiMicroEvent->SetL3RichTrigger(true);
    }
  }
  */
  //
  // 02/21/03 - BUM -trigger id stuff
  //
  // default to 0
  mHiMicroEvent->ClearTriggerIds();

  const StTriggerId &nom_trig=event->triggerIdCollection().nominal();
  vector<unsigned int> vv = nom_trig.triggerIds();
  for(unsigned int i=0;i<vv.size();i++){
    //cout << i << " : " << vv[i] << endl;
    mHiMicroEvent->AddTriggerId(i,vv[i]);
  }
  
  // this is a reality check
  if( (mHiMicroEvent->IsTrigger(2001) || mHiMicroEvent->IsTrigger(2003)) 
      != (nom_trig.isTrigger(2001) || nom_trig.isTrigger(2003))){
    cout << ">>>> WHATS UP WITH THE TRIGGER?" << endl;
  }

  // reality check
  //
  /*
  cout << "-----" << endl;
  for(int i=0;i<32;i++){
    if(mHiMicroEvent->TriggerId(i))
       cout << i << " : " << mHiMicroEvent->TriggerId(i) << endl;
  }
  */

  // taken from flow maker
  //

  StCtbTriggerDetector &CTB = event->ctbTriggerDetector();
  StZdcTriggerDetector &ZDC = event->zdcTriggerDetector();
  StBbcTriggerDetector &BBC = event->bbcTriggerDetector();
  fillCtbZdcBbc(CTB,ZDC,BBC);
  
  // count FTPC
  //if(stEvent->primaryVertex()){
  if(fabs(event->primaryVertexPosition().x())>0.001 || 
     fabs(event->primaryVertexPosition().y())>0.001 ||
     fabs(event->primaryVertexPosition().z())>0.001) { 
    mHiMicroEvent->SetFTPCe(mFTPCe); mHiMicroEvent->SetFTPCw(mFTPCw);
  }
  else{
    mHiMicroEvent->SetFTPCe(-1); mHiMicroEvent->SetFTPCw(-1);
  }

  fillQuadrant();
  
  if (mDebug) {
    cout << "##StHiMicroMaker::fillEvent()" << endl;
    cout << *mHiMicroEvent << endl;
  }
}

//____________________
void
StHiMicroMaker::fillTracks(StEvent* stEvent)
{

  Int_t nHi(0);
  // points to the trigger track for quadrant stuff
  for(int i=0;i<3;i++){ mQuad[i].LeadMom=StThreeVectorF(0,0,0);}

  // count FTPC tracks
  mFTPCe=mFTPCw=0;

  //
  // create the track and hit object to temporarily store the info.
  // these will then be added to the clonesarray
  //

  StHiMicroTrack* hiMicroTrack = new StHiMicroTrack;
  StHiMicroHit*   hiMicroHit   = new StHiMicroHit;

  //
  // pid stuff for contamination
  //
  //StuProbabilityPidAlgorithm uPid(*stEvent);
  StTpcDedxPidAlgorithm tpcDedxAlgo;

  //
  // 03/15/03 - check for primary vertex
  //
  StPrimaryVertex* prVtx=0;
  bool foundVertex=false;
  //if(stEvent->primaryVertex()) { 
  if(stEvent->primaryVertex() &&
     (fabs(stEvent->primaryVertex()->position().x())>0.001 || 
      fabs(stEvent->primaryVertex()->position().y())>0.001 ||
      fabs(stEvent->primaryVertex()->position().z())>0.001)) { 
    foundVertex=true; prVtx=stEvent->primaryVertex(0);
  }

  // jons CTB stuff
  
  vector<ctbHit> *ctbHits=new vector<ctbHit>;
  StCtbTriggerDetector &ctbDet = stEvent->triggerDetectorCollection()->ctb();

  for (unsigned int slat = 0; slat < ctbDet.numberOfSlats(); slat++) {
    for (unsigned int tray = 0; tray < ctbDet.numberOfTrays(); tray++) {
      ctbHit curHit;
      curHit.adc = ctbDet.mips(tray,slat,0);
      if(curHit.adc > 0){
	ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
	ctbHits->push_back(curHit);
      }
    }
  }
  
  //
  // 03/15/03 -loop over nodes - 
  // loop over primary tracks if the vertex is found,
  // all global tracks otherwise
  //
  StSPtrVecTrackNode& nodes = stEvent->trackNodes();  

  for (unsigned int i=0; i<nodes.size(); i++) {
    StTrack* glTrack = nodes[i]->track(global);
    if(!glTrack) { cout << "HUH? " << endl; continue; }
    StTrack* prTrack = nodes[i]->track(primary); // should return 0 if not pr

    // 03/15/03 - if mSaveAllEvents and not a primary track
    //            then glTrack=prTrack
    //          - mSaveAllEvents saves all track nodes
    
    bool isPrimary=false; if(prTrack) isPrimary=true;

    if(!mSaveAllEvents && !prTrack) continue;
    else if(mSaveAllEvents && !prTrack) prTrack=glTrack;

    //
    // find track with highest pt for quadrant stuff
    // only primary tracks for now
    //
    
    if(isPrimary){
      int indx=type(prTrack,stEvent);
      if(indx>=0){
	bool (StHiMicroMaker::* fp)(StTrack*,StEvent*);
	switch(indx){
	case 0: fp=&StHiMicroMaker::acceptLead; break;
	case 1: mFTPCe++; break; 
	case 2: mFTPCw++; break; // count FTPCw
	default: cout << "Unknown indx="<<indx << endl; exit(-1);	
	}
	// find leading TPC track
	if(indx==0 && (this->*fp)(prTrack,stEvent)){
	  if(prTrack->geometry()->momentum().perp()>
	     mQuad[indx].LeadMom.perp())
	    mQuad[indx].LeadMom=prTrack->geometry()->momentum();
	}
      }
    }

    //
    // dont save ftpc tracks for now
    // unless save low pt.  then we save everything
    if(fabs(glTrack->geometry()->momentum().pseudoRapidity())>2.5 &&
       !mSaveLowPt)
      continue;

    //
    // accept the track?
    //
    if(!mSaveLowPt &&
       !accept(glTrack) && !accept(prTrack)) continue;
    

    if(accept(glTrack) || accept(prTrack)) nHi++;
    if(mDebug && (accept(prTrack)||accept(glTrack))){
      cout << "Accepted track" ;
      if(mSaveAllEvents) cout << ": Is primary? " 
			      << (isPrimary ? "yes." : "no.");
      cout << endl;
      dump(prTrack,glTrack);
    }


    // 
    // sort the hits
    //
    
    const StPtrVecHit& hhits = glTrack->detectorInfo()->hits(kTpcId);
    Float_t firstZ(0),lastZ(0);
    Short_t innerPadList(0);
    Int_t firstPadrow(0), lastPadrow(0), outerPadList(0);
    UInt_t firstSector(99), lastSector(99);
    Double_t crossAngle(999);

    vector<StTpcHit*> vec;
    int countHits=0;
    // fill the vector
    for(UInt_t i=0; i<hhits.size(); i++){
      StTpcHit* hit = (StTpcHit*)hhits[i];
      if(!hit) continue;
      countHits++;
      vec.push_back(hit);
    }
    sort(vec.begin(),vec.end(),sortCmp);
    if(vec.size()){
      firstZ=vec[0]->position().z(); lastZ=vec[vec.size()-1]->position().z();
      firstPadrow=vec[0]->padrow(); lastPadrow=vec[vec.size()-1]->padrow();

//JLK  Do bit packing here.  Basically set a mask for each bit to turn on or off according to the padrow numbers
//     in the vec

      Short_t inner = 0;
      Int_t outer = 0;
      Int_t SET_ON;   
      for (UInt_t i=0; i < vec.size(); i++) {
        //printf("vec[%d]->padrow()=%d",i,vec[i]->padrow());
        Int_t val = vec[i]->padrow();
        if (val < 14) { //Inner sector  
          SET_ON = 1 << val;  //Moves "1" over by val places
          inner = inner | SET_ON;
        } else { 
         SET_ON = 1 << val-14;  //Moves "1" over by val places
         outer = outer | SET_ON;
        }
        innerPadList=inner; outerPadList=outer;
      }
      firstSector=vec[0]->sector(); lastSector=vec[vec.size()-1]->sector();
      //crossAngle=crossingAngle(glTrack->helix(),vec[0],stEvent->runInfo()->magneticField());
    }
    else{ 
      if(mDebug)
	cout << "\tno hits" << endl;
      //continue; 
      cout << "No hits on track " << endl;
      cout << "global track: " 
	   << glTrack->detectorInfo()->numberOfPoints() << " hits, ( " 
	   << glTrack->topologyMap().numberOfHits(kSvtId) << " svt, "
	   << glTrack->topologyMap().numberOfHits(kSsdId) << " ssd, " 
	   << glTrack->topologyMap().numberOfHits(kTpcId) << " tpc, " 
	   << glTrack->topologyMap().numberOfHits(kFtpcEastId)+
	glTrack->topologyMap().numberOfHits(kFtpcWestId) << " ftpc )" << endl;
      cout << "fit points: " << glTrack->fitTraits().numberOfFitPoints() << endl;
      cout << "first point: " << glTrack->detectorInfo()->firstPoint() << endl;
      cout << "last point: " << glTrack->detectorInfo()->lastPoint() << endl;
      cout << "primary track: " 
	   << prTrack->detectorInfo()->numberOfPoints() << " hits, ( " 
	   << prTrack->topologyMap().numberOfHits(kSvtId) << " svt, "
	   << prTrack->topologyMap().numberOfHits(kSsdId) << " ssd, " 
	   << prTrack->topologyMap().numberOfHits(kTpcId) << " tpc, " 
	   << prTrack->topologyMap().numberOfHits(kFtpcEastId)+
	prTrack->topologyMap().numberOfHits(kFtpcWestId) << " ftpc )" << endl;
      cout << "fit points: " << prTrack->fitTraits().numberOfFitPoints() << endl;
      cout << "first point: " << prTrack->detectorInfo()->firstPoint() << endl;
      cout << "last point: " << prTrack->detectorInfo()->lastPoint() << endl;
      cout << "pt " << prTrack->geometry()->momentum().perp()
	   << ", eta " << prTrack->geometry()->momentum().pseudoRapidity()
	   << ", phi " << prTrack->geometry()->momentum().phi() << endl;
      cout << "Hit vector global: " << glTrack->detectorInfo()->hits().size()
	   << ", primary : " << prTrack->detectorInfo()->hits().size() << endl;
      cout << "Tpc hit vector global: " << glTrack->detectorInfo()->hits(kTpcId).size()
	   << ", primary : " << prTrack->detectorInfo()->hits(kTpcId).size() << endl;
    }

    //
    // alias to the momentum and helix
    //
    const StThreeVectorF& prMom = prTrack->geometry()->momentum();
    const StThreeVectorF& glMom = glTrack->geometry()->momentum();
    const StPhysicalHelixD& prHelix = prTrack->geometry()->helix();
    const StPhysicalHelixD& glHelix = glTrack->geometry()->helix();

    //
    //***** fill some track info *****
    //
        
    hiMicroTrack->SetPtPr(prMom.perp());
    hiMicroTrack->SetEtaPr(prMom.pseudoRapidity());
    hiMicroTrack->SetPhiPr(prMom.phi());
    
    hiMicroTrack->SetPtGl(glMom.perp());  
    hiMicroTrack->SetEtaGl(glMom.pseudoRapidity());
    hiMicroTrack->SetPhiGl(glMom.phi()); 
   
    //
    // find the sign of the dca	--Could just get this from StPhysicalHelix now - Jamie added it
    //
    if(foundVertex){
      hiMicroTrack->SetDcaGl(glHelix.distance(prVtx->position())); 
      hiMicroTrack->SetDcaXYGl(computeXY(prVtx->position(),glHelix));
      hiMicroTrack->SetDcaZGl(dcaz(glHelix,prVtx->position()));
      
      //cout << glHelix.distance(prVtx->position())
      //	   << " vs " << glTrack->impactParameter() << endl;
     
    }
    else{
      hiMicroTrack->SetDcaGl(999); 
      hiMicroTrack->SetDcaXYGl(999);
      hiMicroTrack->SetDcaZGl(999);
    }

    hiMicroTrack->SetChi2(prTrack->fitTraits().chi2());
    hiMicroTrack->SetIsPrimary(isPrimary);

    //
    // match to CTB from jon gans
    //
    
    bool shouldHitCTB=false;
    double etaInCTBFrame(0),phiInCTBFrame(0);
    
    unsigned int matchToCTB = 
      EtaAndPhiToOrriginAtCTB(glHelix,ctbHits,shouldHitCTB,etaInCTBFrame,phiInCTBFrame);
    
    hiMicroTrack->SetMatchToCTB(matchToCTB);
    
    //
    // dca to beam line from jon
    //
    
    if (mBeamLineHelix) {
      pairD pathLengths = glHelix.pathLengths(*mBeamLineHelix);
      double dcaToBeamLine = (glHelix.at(pathLengths.first) - mBeamLineHelix->at(pathLengths.second)).mag();
      hiMicroTrack->SetDcaToBeamLine((Float_t)dcaToBeamLine);
    }
    else
      hiMicroTrack->SetDcaToBeamLine(-99);
    //
    // dedx
    //
    StDedxPidTraits* pid=0;
    StPtrVecTrackPidTraits traits = prTrack->pidTraits(kTpcId);
    
    for (UInt_t i = 0; i < traits.size(); i++) {
      pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
      if (pid && pid->method() == kTruncatedMeanId) break;
    }
    
    hiMicroTrack->SetDedx((pid) ? pid->mean() : 0);
    hiMicroTrack->SetDedxPts((pid) ? pid->numberOfPoints() : 0);

    hiMicroTrack->SetFitPts(glTrack->fitTraits().numberOfFitPoints(kTpcId));
    hiMicroTrack->SetAllPts(glTrack->detectorInfo()->numberOfPoints(kTpcId));
    hiMicroTrack->SetMaxPossPts(glTrack->numberOfPossiblePoints());
    hiMicroTrack->SetFlag(glTrack->flag());
    hiMicroTrack->SetCharge(glTrack->geometry()->charge());

    // hit stuff
    //
    hiMicroTrack->SetFirstZ(firstZ);
    hiMicroTrack->SetLastZ(lastZ);
    hiMicroTrack->SetFirstPadrow(firstPadrow);
    hiMicroTrack->SetInnerPadList(innerPadList);
    hiMicroTrack->SetOuterPadList(outerPadList);
    hiMicroTrack->SetLastPadrow(lastPadrow);
    hiMicroTrack->SetFirstSector(firstSector);
    hiMicroTrack->SetLastSector(lastSector);

    // 01/27/03 - note, at first hit for global
    hiMicroTrack->SetCrossingAngle(crossAngle);  

    Int_t curTrackIndex = mHiMicroEvent->NTrack(); // 1st track is 0
    Int_t firstHitIndex = mHiMicroEvent->NHit(); // first hit is 0
    Int_t nHit(0); // just in case
    
    // ****** find some hit info ******
    //
       
    if (mHitLoop ) {
    
      const StPtrVecHit& hits = prTrack->detectorInfo()->hits(kTpcId);
      for(UInt_t iHit=0 ; iHit<hits.size(); iHit++){
	const StTpcHit* hit = dynamic_cast<StTpcHit*> (hits[iHit]);
	if(!hit){cout << "Not a tpc hit?" << endl; continue;
	}
	nHit++;
	// 01/27/03 include all hits
	//if (!hit->usedInFit()) continue;
	
	//
	// repeat some info in the hits
	// 01/27/03 BUM killed redundant info

	//
	// general hit info
	//
	hiMicroHit->SetUsedInFit((Bool_t)hit->usedInFit());
	hiMicroHit->SetFlag(hit->flag());
	hiMicroHit->SetR(hit->position().perp());
	hiMicroHit->SetZ(hit->position().z());   
	hiMicroHit->SetPadRow( (Int_t) hit->padrow());
	hiMicroHit->SetSector( (Int_t) hit->sector());
	
	// double check
	hiMicroHit->SetTrackIndex(curTrackIndex);

	// 01/27/03 crossing angle
	// assume MagneticField is in kG
	float bfield = mHiMicroEvent->MagneticField()*0.1*tesla;
	hiMicroHit->SetCrossingAngleGl(crossingAngle(glHelix,hit,bfield));
	hiMicroHit->SetCrossingAnglePr(crossingAngle(prHelix,hit,bfield));
	
	//
	// residuals
	//
	//StThreeVectorF glHitPos(glHelix.at(sGl));
	//StThreeVectorF prHitPos(prHelix.at(sPr));
	
	hiMicroHit->SetZResGl(dcaz(glHelix,hit->position()));
	hiMicroHit->SetZResPr(dcaz(prHelix,hit->position()));
	
	Float_t resXYGl = computeXY(hit->position(),glTrack->geometry()->helix());
	Float_t resXYPr = computeXY(hit->position(),prTrack->geometry()->helix());
	hiMicroHit->SetXYResGl(resXYGl);
	hiMicroHit->SetXYResPr(resXYPr);
	//
	// add the hit!
	//
	
	mHiMicroEvent->AddHit(hiMicroHit);
      } // hit loop
    } // if do hit loop

    hiMicroTrack->SetFirstHitIndex(firstHitIndex);
    hiMicroTrack->SetLastHitIndex(firstHitIndex+nHit-1);
    hiMicroTrack->SetNHit(nHit);

    //
    // add the track!
    //
    mHiMicroEvent->AddTrack(hiMicroTrack);


  } // track loop

  if(hiMicroTrack) delete hiMicroTrack;
  if(hiMicroHit) delete hiMicroHit;
  
  // jon
  if(ctbHits) delete ctbHits;
  cout << "\thi pt tracks : " << nHi << endl;
  //  cout << "\tFTPC tracks : " << endl;
  mNAcceptedTrack += nHi;

  //return nGoodTrackEta;
    
}

void
StHiMicroMaker::fillTracks(StMuDst* mudst)
{
  StMuEvent *event=mudst->event();
  Int_t nHi(0);
  // points to the trigger track for quadrant stuff
  for(int i=0;i<3;i++){ mQuad[i].LeadMom=StThreeVectorF(0,0,0);}

  // count FTPC tracks
  mFTPCe=mFTPCw=0;

  //
  // create the track and hit object to temporarily store the info.
  // these will then be added to the clonesarray
  //

  StHiMicroTrack* hiMicroTrack = new StHiMicroTrack;

  //
  // pid stuff for contamination
  //
  //StuProbabilityPidAlgorithm uPid(*stEvent);
  StTpcDedxPidAlgorithm tpcDedxAlgo;

  //
  // 03/15/03 - check for primary vertex
  //
  StThreeVectorF prVtx;
  bool foundVertex=false;
  //if(stEvent->primaryVertex()) { 
  if(fabs(event->primaryVertexPosition().x())>0.001 || 
     fabs(event->primaryVertexPosition().y())>0.001 ||
     fabs(event->primaryVertexPosition().z())>0.001) { 
    foundVertex=true; 
    prVtx=event->primaryVertexPosition();
  }

  // jons CTB stuff
  
  vector<ctbHit> *ctbHits=new vector<ctbHit>;
  StCtbTriggerDetector &ctbDet = event->ctbTriggerDetector();

  for (unsigned int slat = 0; slat < ctbDet.numberOfSlats(); slat++) {
    for (unsigned int tray = 0; tray < ctbDet.numberOfTrays(); tray++) {
      ctbHit curHit;
      curHit.adc = ctbDet.mips(tray,slat,0);
      if(curHit.adc > 0){
	ctb_get_slat_from_data(slat, tray, curHit.phi, curHit.eta);
	ctbHits->push_back(curHit);
      }
    }
  }
  
  //
  // 03/15/03 -loop over nodes - 
  // loop over primary tracks if the vertex is found,
  // all global tracks otherwise
  //
  // 09/02/05 - Suppport for copying all global tracks not implemented
  // fo MuDst. Need StEvent for that

  unsigned int n_prim=mudst->numberOfPrimaryTracks();
  
  //cout<<n_prim<<endl;

  for (unsigned int i=0; i < n_prim; i++) {
    StMuTrack* prTrack = mudst->primaryTracks(i);
    StMuTrack *glTrack = prTrack->globalTrack();

    // 03/15/03 - if mSaveAllEvents and not a primary track
    //            then glTrack=prTrack
    //          - mSaveAllEvents saves all track nodes
    
    bool isPrimary=true;

    //cout<<n_prim<<" "<<i<<endl;

    //
    // find track with highest pt for quadrant stuff
    // only primary tracks for now
    //
    
    if(isPrimary){
      int indx=type(prTrack,event);
      if(indx>=0){
	bool (StHiMicroMaker::* fp)(StMuTrack*,StMuEvent*);
	switch(indx){
	case 0: fp=&StHiMicroMaker::acceptLead; break;
	  //case 1: mFTPCe++; break; 
	  //case 2: mFTPCw++; break; // count FTPCw
	default: cout << "Unknown indx="<<indx << endl; exit(-1);	
	}
	// find leading TPC track
	if(indx==0 && (this->*fp)(prTrack,event)){
	  if(prTrack->pt() > mQuad[indx].LeadMom.perp())
	    mQuad[indx].LeadMom=prTrack->momentum();
	  
	}
      }
    }

    //
    // dont save ftpc tracks for now
    // unless save low pt.  then we save everything
    // JOERN
    //if(fabs(prTrack->eta())>2.5 &&
    // !mSaveLowPt)
    //continue;

    //
    // accept the track?
    //

    // if (fabs(prTrack->eta())<2.0)
//       {
// 	if(!mSaveLowPt &&
// 	   !accept(glTrack) && !accept(prTrack)) continue;
	
// 	if(accept(glTrack) || accept(prTrack)) nHi++;
	
	
// 	if(mDebug && (accept(prTrack)||accept(glTrack))){
// 	  cout << "Accepted track" ;
// 	  if(mSaveAllEvents) cout << ": Is primary? " 
// 				  << (isPrimary ? "yes." : "no.");
// 	  cout << endl;
// 	  dump(prTrack,glTrack);
// 	}
//       }
//     else
//       {

// 	if (prTrack->eta()>0)
// 	  mFTPCw++;
// 	else
// 	  mFTPCe++;

// 	if (prTrack->pt()>=1.5)
// 	  {
// 	    nHi++;
// 	    if(mDebug)
// 	      {
// 		if (glTrack) {
// 		  dump(prTrack,glTrack);
// 		}
// 		else
// 		  {
// 		    cout<<"Global missing print pr:"<<endl;
// 		    dump(prTrack,prTrack);
// 		  }
// 		cout<<"-----"<<endl;
// 	      }
// 	  }
// 	else continue;
//       }

    if(!mSaveLowPt &&
       !accept(glTrack) && !accept(prTrack)) continue;


    if(accept(glTrack) || accept(prTrack)) nHi++;
    if(mDebug && (accept(prTrack)||accept(glTrack))){
      cout << "Accepted track" ;
      if(mSaveAllEvents) cout << ": Is primary? "
                              << (isPrimary ? "yes." : "no.");
      cout << endl;
      dump(prTrack,glTrack);
    }


    //
    //***** fill some track info *****
    //
        
    hiMicroTrack->SetPtPr(prTrack->pt());
    hiMicroTrack->SetEtaPr(prTrack->eta());
    hiMicroTrack->SetPhiPr(prTrack->phi());

    // Fill Pid infos

    hiMicroTrack->SetnSigmaElectron(prTrack->nSigmaElectron());
    hiMicroTrack->SetnSigmaPion(prTrack->nSigmaPion());
    hiMicroTrack->SetnSigmaKaon(prTrack->nSigmaKaon());
    hiMicroTrack->SetnSigmaProton(prTrack->nSigmaProton());

    // Debug
    //cout<<prTrack->nSigmaPion()<<" "<<hiMicroTrack->nSigmaPion()<<" "<<nHi<<endl;

    if (glTrack) {
    hiMicroTrack->SetPtGl(glTrack->pt());  
    hiMicroTrack->SetEtaGl(glTrack->eta());
    hiMicroTrack->SetPhiGl(glTrack->phi()); 
     }
    else {
      // This is a rare case
      hiMicroTrack->SetPtGl(-99);
      hiMicroTrack->SetEtaGl(-99);
      hiMicroTrack->SetPhiGl(-99); 
    }
    //
    // find the sign of the dca	--Could just get this from StPhysicalHelix now - Jamie added it
    //
    hiMicroTrack->SetDcaGl(prTrack->dcaGlobal().mag()); 
    
    if (glTrack) {
    hiMicroTrack->SetDcaXYGl(computeXY(prVtx,glTrack->helix()));
    hiMicroTrack->SetDcaZGl(dcaz(glTrack->helix(),prVtx));      
    }
    else {
      // This is a rare case
      hiMicroTrack->SetDcaXYGl(-99);
      hiMicroTrack->SetDcaZGl(-99);
    }
     
    hiMicroTrack->SetChi2(prTrack->chi2());
    hiMicroTrack->SetIsPrimary(isPrimary);

    //
    // match to CTB from jon gans
    //
    
    // JOERN

    bool shouldHitCTB=false;
    double etaInCTBFrame(0),phiInCTBFrame(0);
    unsigned int matchToCTB=0;

     if (glTrack && fabs(prTrack->eta())<1.5) {
       matchToCTB = 
	 EtaAndPhiToOrriginAtCTB(glTrack->helix(),ctbHits,shouldHitCTB,etaInCTBFrame,phiInCTBFrame);
     }
     else matchToCTB=0;

    hiMicroTrack->SetMatchToCTB(matchToCTB);
    
    //
    // dca to beam line from jon
    //

    if (mBeamLineHelix && glTrack && fabs(prTrack->eta())<1.5) {
      pairD pathLengths = glTrack->helix().pathLengths(*mBeamLineHelix);
      double dcaToBeamLine = (glTrack->helix().at(pathLengths.first) - mBeamLineHelix->at(pathLengths.second)).mag();
      hiMicroTrack->SetDcaToBeamLine((Float_t)dcaToBeamLine);
    }
    else
      hiMicroTrack->SetDcaToBeamLine(-99);

    //
    // dedx
    //
    hiMicroTrack->SetDedx(prTrack->dEdx());
    hiMicroTrack->SetDedxPts(prTrack->nHitsDedx());

    if (glTrack) {
      // So, normally we take these thinsg from the global track 
      // (for historical reasons)
      hiMicroTrack->SetFitPts(glTrack->nHitsFit(kTpcId));
      hiMicroTrack->SetAllPts(glTrack->nHits());
      hiMicroTrack->SetMaxPossPts(glTrack->nHitsPoss());
      hiMicroTrack->SetFlag(glTrack->flag());
      hiMicroTrack->SetCharge(glTrack->charge());
    }
    else {
      // But in the rare case that there is no global track, 
      // we use the primary instead
      hiMicroTrack->SetFitPts(prTrack->nHitsFit(kTpcId));
      hiMicroTrack->SetAllPts(prTrack->nHits());
      hiMicroTrack->SetMaxPossPts(prTrack->nHitsPoss());
      hiMicroTrack->SetFlag(prTrack->flag());
      hiMicroTrack->SetCharge(prTrack->charge());
    }

    // hit stuff
    //
    hiMicroTrack->SetFirstZ(prTrack->firstPoint().z());
    hiMicroTrack->SetLastZ(prTrack->lastPoint().z());
    /*
    hiMicroTrack->SetFirstPadrow(firstPadrow);
    hiMicroTrack->SetInnerPadList(innerPadList);
    hiMicroTrack->SetOuterPadList(outerPadList);
    hiMicroTrack->SetLastPadrow(lastPadrow);
    hiMicroTrack->SetFirstSector(firstSector);
    hiMicroTrack->SetLastSector(lastSector);
    */
    hiMicroTrack->SetFirstPadrow(-1);
    hiMicroTrack->SetInnerPadList(-1);
    hiMicroTrack->SetOuterPadList(-1);
    hiMicroTrack->SetLastPadrow(-1);
    hiMicroTrack->SetFirstSector(-1);
    hiMicroTrack->SetLastSector(-1);

    // 01/27/03 - note, at first hit for global
    //hiMicroTrack->SetCrossingAngle(crossAngle);  
    hiMicroTrack->SetCrossingAngle(-99);
    
    // No Hits in MuDst
    hiMicroTrack->SetFirstHitIndex(-1);
    hiMicroTrack->SetLastHitIndex(-1);
    hiMicroTrack->SetNHit(0);

    //
    // add the track!
    //
    //cout<<"Before Add track ..."<<endl;
    mHiMicroEvent->AddTrack(hiMicroTrack);
    //cout<<"After Add track ..."<<endl;


  } // track loop

  if(hiMicroTrack) delete hiMicroTrack;
  
  // jon
  //if(ctbHits) delete ctbHits;
  if (mDebug)
    cout << "\thi pt tracks : " << nHi << endl;
  //  cout << "\tFTPC tracks : " << endl;
  mNAcceptedTrack += nHi;

  //return nGoodTrackEta;
}



//____________________
void
StHiMicroMaker::fillFpdCollection(StEvent* stEvent)
{

  //Get the fpd collection from StEvent
  StFpdCollection* stFpd = stEvent->fpdCollection();

  mHiMicroEvent->fpd()->setToken(stFpd->token());

  for (unsigned int i=0;i<mHiMicroEvent->fpd()->numberOfADC();i++){
    unsigned short* adc = stFpd->adc();
    mHiMicroEvent->fpd()->setAdc(i,adc[i]);
  }
  for (unsigned int i=0;i<mHiMicroEvent->fpd()->numberOfTDC();i++){
    unsigned short* tdc = stFpd->tdc();
    mHiMicroEvent->fpd()->setTdc(i,tdc[i]);
  }
  for (unsigned int i=0;i<mHiMicroEvent->fpd()->numberOfRegisters();i++){
    mHiMicroEvent->fpd()->setRegister(i,stFpd->registers(i));
  }
  for (unsigned int i=0;i<mHiMicroEvent->fpd()->numberOfPedestal();i++){
    unsigned short* ped = stFpd->pedestal();
    mHiMicroEvent->fpd()->setPedestal(i,ped[i]);
  }
  for (unsigned int i=0;i<mHiMicroEvent->fpd()->numberOfScalers();i++){
    mHiMicroEvent->fpd()->setScaler(i,stFpd->scaler(i));
  }

}

//___________________
// second loop to do quadrant stuff

void
StHiMicroMaker::computeQuadrant(StEvent* stEvent)
{
  if (mDebug)
    cout << "Computing quadrants..." << endl;

  float leadPhi[3]={9999};
  for(int i=0; i<3; i++){
    mQuad[i].QLead=0; mQuad[i].QOpp=0; mQuad[i].QA=0; mQuad[i].QB=0;
    
  }
  if(mQuad[0].LeadMom.x() != 0 || mQuad[0].LeadMom.y() != 0) {
    leadPhi[0]=mQuad[0].LeadMom.phi();
  } else {
    cout << "No leading track? " << endl; return;
  }


  // get out if no primary vertex
  if(!stEvent->primaryVertex()) return;

  StSPtrVecPrimaryTrack& prTracks = stEvent->primaryVertex(0)->daughters();
  
  for(UInt_t i=0; i<prTracks.size(); i++){
    StPrimaryTrack* prTrack = prTracks[i];
    if(!prTrack) { cout << "No primary?" << endl; continue; }

    //StTrack* glTrack=prTrack->node()->track(global);

    // nch, ftcpe, or ftpcw, or none of above
    int j=type(prTrack,stEvent);
    if(j<0) continue;

    float phi = prTrack->geometry()->momentum().phi();
    float dphi= phi-leadPhi[0]; // w/ respect to the TPC track
    
    dphi = (fabs(dphi) < M_PI ) ? dphi :
      ((dphi<0) ? 2.0*M_PI + dphi : 2*M_PI - dphi);
    if(dphi<0) dphi += M_PI*2.0;
    
    if(dphi>=7.0*M_PI/4.0 || dphi<M_PI/4.0) mQuad[j].QLead++;
    else if(dphi >= M_PI/4.0 && dphi < 3.0*M_PI/4) mQuad[j].QA++;
    else if(dphi >= 3.0*M_PI/4.0 && dphi < 5.0*M_PI/4.0) mQuad[j].QOpp++;
    else mQuad[j].QB++;
  
  }
  if (mDebug)
    cout << "...done" << endl;

  //if(mNchQLead==0) { cout << "HUH? None in Lead Q?" << endl; }

}

//___________________
// second loop to do quadrant stuff

void
StHiMicroMaker::computeQuadrant(StMuDst* mudst)
{
  if (mDebug)
    cout << "Computing quadrants..." << endl;
  StMuEvent *event=mudst->event();

  float leadPhi[3]={9999};
  for(int i=0; i<3; i++){
    mQuad[i].QLead=0; mQuad[i].QOpp=0; mQuad[i].QA=0; mQuad[i].QB=0;
    
  }
  if( mQuad[0].LeadMom.x()!=0 || mQuad[0].LeadMom.y()!=0) {
      leadPhi[0]=mQuad[0].LeadMom.phi();
  } else {
    cout << "No leading track? " << endl; return;
  }


  // get out if no primary vertex
  UInt_t n_prim=mudst->numberOfPrimaryTracks();
  for(UInt_t i=0; i < n_prim; i++){
    StMuTrack* prTrack = mudst->primaryTracks(i);

    // nch, ftcpe, or ftpcw, or none of above
    int j=type(prTrack,event);
    if(j<0) continue;

    float phi = prTrack->phi();
    float dphi= phi-leadPhi[0]; // w/ respect to the TPC track
    
    dphi = (fabs(dphi) < M_PI ) ? dphi :
      ((dphi<0) ? 2.0*M_PI + dphi : 2*M_PI - dphi);
    if(dphi<0) dphi += M_PI*2.0;
    
    if(dphi>=7.0*M_PI/4.0 || dphi<M_PI/4.0) mQuad[j].QLead++;
    else if(dphi >= M_PI/4.0 && dphi < 3.0*M_PI/4) mQuad[j].QA++;
    else if(dphi >= 3.0*M_PI/4.0 && dphi < 5.0*M_PI/4.0) mQuad[j].QOpp++;
    else mQuad[j].QB++;
  
  }
  if (mDebug)
    cout << "...done" << endl;

  //if(mNchQLead==0) { cout << "HUH? None in Lead Q?" << endl; }

}

//____________________
void
StHiMicroMaker::dump(StTrack* prTrack,StTrack* glTrack)
{
  cout << "\tptGl=" 
       << glTrack->geometry()->momentum().perp() 
       << ",etaGl=" 
       << glTrack->geometry()->momentum().pseudoRapidity() 
       << ",fitpts="
       << glTrack->fitTraits().numberOfFitPoints(kTpcId)
       << ",allpts=" 
       <<  glTrack->detectorInfo()->numberOfPoints(kTpcId)<<endl
       << "\tptPr=" <<prTrack->geometry()->momentum().perp()
       << ",etaPr="
       << prTrack->geometry()->momentum().pseudoRapidity() <<endl;
}
void
StHiMicroMaker::dump(StMuTrack* prTrack,StMuTrack* glTrack)
{
  cout << "\tptGl=" 
       << glTrack->pt() 
       << ",etaGl=" 
       << glTrack->eta()
       << ",fitpts="
       << glTrack->nHitsFit(kTpcId)
       << ",allpts=" 
       <<  glTrack->nHits()<<endl
       << "\tptPr=" << prTrack->pt()
       << ",etaPr="
       << prTrack->eta() <<endl;
}
//____________________

Int_t
StHiMicroMaker::openFile()
{
  cout << "###StHiMicroMaker::openFile()" << endl;
  
  //
  // for the output root file, replace event.root with himicro.root
  //
  TString outFileName(mInFileName);

  //
  // assume the extension is of the form ( {\..*?\.root} (perl regexp)) 
  //

  char cTemp[100];
  strcpy(cTemp,outFileName.Data());
  
  cout << "We think the filename is " << outFileName.Data() << endl;
  TString replace;
  if(strstr(cTemp,"dst.root")) replace=".dst.root";
  else if(strstr(cTemp,"event.root")) replace=".event.root";
  else if(strstr(cTemp,"MuDst.root")) replace=".MuDst.root";
  else { cout << "unknown extension " << endl; exit(1); }
  
  //  outFileName.Replace(outFileName.Length(),
  //		      replace.Length(),".himicro.root");
  cout << replace << endl;

  outFileName.ReplaceAll(replace.Data(),".himicro.root");

  cout << outFileName << endl;

  outFileName.Prepend(mOutDir + "/");

  // use default compression 1
  mDSTFile = new TFile(outFileName.Data(),"RECREATE");

  if(!mDSTFile){ 
    gMessMgr->Error() 
      << "Cannot create = " << outFileName << endm;
    return kStErr;
  }
  //mDSTFile->SetFormat(1); 
  cout << "\toutfile = " << outFileName << endl;

  //
  // top level tree
  //
  if(mDebug) 
    cout << "##Creating the top level tree..." << endl;
  mDSTTree = new TTree("StHiMicroTree","StHiMicroTree");
  cout << "mDSTTree created at " << mDSTTree << endl; 
  if(!mDSTTree){
    gMessMgr->Error() 
      << "Cannot create mDSTTree" << endm;
    return kStErr;
  }
  if(mDebug)cout << "##...done" << endl;
  //#if ROOT_VERSION_CODE >= ROOT_VERSION(3,01,05)
  //  mDSTTree->SetBranchStyle(0);
  //#endif 
  mDSTTree->SetAutoSave(10000000); // 10 MB
  
  // set the event branch
  Int_t bufSZ = 64000;
  mDSTTree->Branch("StHiMicroEvent","StHiMicroEvent",
		   &mHiMicroEvent, bufSZ,1);
  return kStOk;
}
//__________________________

Int_t
StHiMicroMaker::closeFile()
{
  cout << "###StHiMicroMaker::closeFile()" << endl;
  
  cout << "##Writing " << mInFileName << endl;

  if(mDSTFile && mDSTFile->IsOpen()){
    mDSTFile->Write();
    mDSTFile->Close();
  }

  cout << "##...done\n";

  return kStOk;
}
  
//____________________________
// from jon gans

void
StHiMicroMaker::getBeamLine()
{
  if(mDebug) cout << "StHiMicroMaker::getBeamLine()..." << endl;

  TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic");
  if(!dbDataSet) {
    cout << "Couldn't get the database, not using beam line info" << endl;
    mBeamLineHelix = 0;
    return;
  }

  vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
  double x0 = vSeed->x0;
  double y0 = vSeed->y0;
  double dxdz = vSeed->dxdz;
  double dydz = vSeed->dydz;
  cout << "BeamLine Constraint: " << endl;
  cout << "x(z) = " << x0 << " + " << dxdz << " * z" << endl;
  cout << "y(z) = " << y0 << " + " << dydz << " * z" << endl << endl;
  double pt  = 88889999;   
  double nxy=sqrt(dxdz*dxdz + dydz*dydz);
  if(nxy<1.e-5){ // beam line _MUST_ be tilted
    nxy=dxdz=1.e-5;
  } 
  double p0=pt/nxy; 
  double px   = p0*dxdz;
  double py   = p0*dydz;
  double pz   = p0; // approximation: nx,ny<<0
  StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
  StThreeVectorD origin(x0, y0, 0.);

  //
  // NOTE THE HARD CODED B FIELD
  //
  mBeamLineHelix = new StPhysicalHelixD (MomFstPt, origin, 0.5*tesla, 1.);

  cout << *mBeamLineHelix << endl;

  if(mDebug) cout << "...done" << endl;
}


Float_t
StHiMicroMaker::computeXY(const StThreeVectorF& pos, const StPhysicalHelixD &helix)
{
  //
  // find the distance between the center of the circle and pos.
  // if this greater than the radius of curvature, then call 
  // it negative.
  //
  double xCenter = helix.xcenter();
  double yCenter = helix.ycenter();
  double radius  = 1.0/helix.curvature();
  
  double dPosCenter 
    = TMath::Sqrt( (pos.x()-xCenter) * (pos.x()-xCenter) +
		   (pos.y()-yCenter) * (pos.y()-yCenter));

  return (Float_t) radius - dPosCenter;
}


// 
// i like this one better. dip angle in s-z plane
// should be exact assuming a circle in x-y
// (redundant StTrack* for debugging)
double StHiMicroMaker::dcaz(const StPhysicalHelixD& helix, 
			    const StThreeVectorF& point)  
{  
  double z0       = helix.origin().z();
  double phi      = atan2(point.y()-helix.ycenter(),
                          point.x()-helix.xcenter());
  int    h        = helix.h(); // -sign(q*B) (+ := ccw looking down from +z)
  double dphi     = h*(phi-helix.phase());

  // half circle assumption
  dphi           = (fabs(dphi) < M_PI ) ? dphi :
                    ((dphi<0) ? 2*M_PI + dphi : 2*M_PI - dphi);
   
  double arclength= (1./helix.curvature()) * dphi;
  
  double dcaZ =  (point.z() - (z0 + arclength*tan(helix.dipAngle())));
  /*
  if(gRandom->Rndm(1)<0.1){
    cout << "--------" << endl;
    cout << "zpoint=" << point.z() << endl;
    if(track)
      cout << "pt=" << track->geometry()->momentum().perp()
           << ",fit pts=" << track->fitTraits().numberOfFitPoints(kTpcId)
           << endl;
  
    cout << ">>>zdca=" << dcaZ
         << ", dphi=" << dphi*180./M_PI
         << ", z0=" << z0 << ",arclength=" << arclength << endl
         << ", dip=" << helix.dipAngle()*180./M_PI
         << ", arc/tan(dip)="<< arclength/tan(helix.dipAngle()) << endl
         << ">>>dca.z= " << point.z()-helix.at(helix.pathLength(point)).z()
         << endl;
    cout << "--------" << endl;
  }
  */
  return dcaZ;
         
}

//____________________
//
// primary vertex exists 
//
bool StHiMicroMaker::accept(StEvent* stEvent)
{

  //return (stEvent->primaryVertex());
  return (stEvent->primaryVertex() &&
     (fabs(stEvent->primaryVertex()->position().x())>0.001 || 
      fabs(stEvent->primaryVertex()->position().y())>0.001 ||
      fabs(stEvent->primaryVertex()->position().z())>0.001)); 
}

//
// primary vertex exists 
//
bool StHiMicroMaker::accept(StMuEvent* event)
{

  //return (stEvent->primaryVertex());
  return (fabs(event->primaryVertexPosition().x())>0.001 || 
	  fabs(event->primaryVertexPosition().y())>0.001 ||
	  fabs(event->primaryVertexPosition().z())>0.001); 
}


//____________________
//
bool StHiMicroMaker::accept(StTrack* track)
{

  return (track && track->flag() > 0 &&
	  track->geometry()->momentum().perp()>=1.5);
}

bool StHiMicroMaker::accept(StMuTrack* track)
{

  return (track && track->flag() > 0 &&
	  track->pt()>=1.5);
}
//____________________
//
// only works if the input is a primary track

bool StHiMicroMaker::acceptNch(StTrack* track,StEvent* stEvent){
  return (track && track->flag() &&  
	  track->fitTraits().numberOfFitPoints(kTpcId)>10 &&
	  fabs(track->geometry()->momentum().pseudoRapidity()) < 0.5&&
	  track->geometry()->helix().distance(stEvent->primaryVertex()->position())<3);
}

bool StHiMicroMaker::acceptNch(StMuTrack* track,StMuEvent* event){
  return (track && track->flag() &&  
	  track->nHitsFit(kTpcId)>10 &&
	  fabs(track->eta()) < 0.5 &&
	  track->dcaGlobal().mag() < 3);  
  // DCA cut above seems to be primary DCA, took global here....
}
//_____________________
//
// 
bool StHiMicroMaker::acceptLead(StTrack* track,StEvent* stEvent){
  return (track && track->flag() && 
	  track->fitTraits().numberOfFitPoints(kTpcId)>=20 &&
	  fabs(track->geometry()->momentum().pseudoRapidity()) < 0.5 &&
	  track->geometry()->helix().distance(stEvent->primaryVertex()->position())<1 &&
	  track->geometry()->momentum().perp()<12);
}
bool StHiMicroMaker::acceptLead(StMuTrack* track,StMuEvent* event){
  return (track && track->flag() && 
	  track->nHitsFit(kTpcId)>=20 &&
	  fabs(track->eta()) < 0.5 &&
	  track->dcaGlobal().mag() <1 &&
	  track->pt() < 12);
}

//_____________________
//

bool StHiMicroMaker::acceptFTPC(StTrack* prTrack,StEvent* stEvent){
  if(!prTrack) return false;
  StTrack* glTrack=prTrack->node()->track(global);
  
  return (prTrack && stEvent->primaryVertex() &&
	  glTrack->geometry()->helix().distance(stEvent->primaryVertex()->position())<3&&
	  // glTrack->geometry()->momentum().perp() > 0 &&
	  prTrack->geometry()->momentum().perp() < 3 &&
	  glTrack->fitTraits().numberOfFitPoints(kTpcId) >=5 &&
	  fabs(prTrack->geometry()->momentum().pseudoRapidity())>=2.8 &&
	  fabs(prTrack->geometry()->momentum().pseudoRapidity())<3.8
	  );
}

bool StHiMicroMaker::acceptFTPC(StMuTrack* prTrack,StMuEvent* event){
  if(!prTrack) return false;
  StMuTrack* glTrack=prTrack->globalTrack();
  
  if (glTrack) {
    return (prTrack && prTrack->dcaGlobal().mag() <3&&
	    // glTrack->geometry()->momentum().perp() > 0 &&
	    prTrack->pt() < 3 &&
	    glTrack->nHitsFit(kTpcId) >=5 &&
	    fabs(prTrack->eta())>=2.8 &&
	    fabs(prTrack->eta())<3.8);
  }
  else {
    // if no global track use fitpoints (-1) from primary track
    return (prTrack && prTrack->dcaGlobal().mag() <3&&
	    // glTrack->geometry()->momentum().perp() > 0 &&
	    prTrack->pt() < 3 &&
	    prTrack->nHitsFit(kTpcId)-1 >=5 &&
	    fabs(prTrack->eta())>=2.8 &&
	    fabs(prTrack->eta())<3.8);
  }
}

//_____________________

bool StHiMicroMaker::goodGlobal(StTrack *track) 
{
  return (track && track->flag() > 0 && track->geometry()->charge() != 0 &&
	   track->fitTraits().numberOfFitPoints(kTpcId) >= 15);
  
}

bool StHiMicroMaker::goodGlobal(StMuTrack *track) 
{
  return (track && track->flag() > 0 && track->charge() != 0 &&
	   track->nHitsFit(kTpcId) >= 15);
  
}

bool StHiMicroMaker::goodGlobalCTBMatch(StTrack *track,unsigned int matchToCTB) 
{

  return (track && track->geometry()->charge() != 0 &&
	   track->fitTraits().numberOfFitPoints(kTpcId) >= 15 && matchToCTB);

}

bool StHiMicroMaker::goodGlobalCTBMatch(StMuTrack *track,unsigned int matchToCTB) 
{

  return (track && track->charge() != 0 &&
	   track->nHitsFit(kTpcId) >= 15 && matchToCTB);

}

bool StHiMicroMaker::goodGlobalFlag(StTrack *track) 
{
  return (track && track->flag() > 0 );

}

bool StHiMicroMaker::goodGlobalFlag(StMuTrack *track) 
{
  return (track && track->flag() > 0 );

}
//_________________

Int_t 
StHiMicroMaker::type(StTrack* track,StEvent* event){
  if(acceptNch(track,event)) return 0;
  else if(acceptFTPC(track,event)&&track->geometry()->momentum().pseudoRapidity()<0) return 1;
  else if(acceptFTPC(track,event)&&track->geometry()->momentum().pseudoRapidity()>0) return 2;
  else return -1;
}

Int_t 
StHiMicroMaker::type(StMuTrack* track,StMuEvent* event){
  if(acceptNch(track,event)) return 0;
  else if(acceptFTPC(track,event)&&track->eta()<0) return 1;
  else if(acceptFTPC(track,event)&&track->eta()>0) return 2;
  else return -1;
}

//
ClassImp(StHiMicroMaker)