/* ************************************************************************* */
/*  $Id: TauEmbedModule.cc,v 1.1 2003/03/27 18:44:28 murat Exp $         */
/* Hadronic Tau reconstruction parameters                                    */
/* Author:    P.Murat(FNAL)                                                  */
/* Date:      08/21/2002                                                     */
/* ************************************************************************* */

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
#include <math.h>
#include "inc/bcs.h"
#include "BaBar/Cdf.hh"
#include "ErrorLogger_i/gERRLOG.hh"

#include "RawDataBanks/COTD_StorableBank.hh"

#include "TFolder.h"
#include "TauEmbedModule.hh"

#include "Edm/Handle.hh"

#include "Trybos/TRY_Bank_Number.hh"
#include "ElectronObjects/CdfEmObject.hh"

#include "CalorObjects/CalData.hh"
#include "CalorObjects/CesCluster.hh"
#include "CalorObjects/CprCluster.hh"
#include "TrackingObjects/Storable/CT_HitSet.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/COTQ.hh"
#include "RawDataBanks/CESD_StorableBank.hh"
#include "RawDataBanks/CPRD_StorableBank.hh"
#include "RawDataBanks/PESD_StorableBank.hh"
#include "RawDataBanks/TOFD_StorableBank.hh"
#include "StorableBanks/IMUE_StorableBank.hh"

#include "SimulationObjects/HEPG_StorableBank.hh"
#include "SimulationObjects/MVTX_StorableBank.hh"
#include "TrackingObjects/Storable/StorableRun2SiStripSet.hh"

#include "Stntuple/obj/TStnNode.hh"
#include "Stntuple/obj/TStnEvent.hh"
#include "Stntuple/alg/TStnElectronID.hh"
#include "Stntuple/obj/TStnElectronBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/obj/TCesDataBlock.hh"
#include "Stntuple/obj/TCprDataBlock.hh"
#include "Stntuple/mod/StntupleMakerModule.hh"
#include "Stntuple/mod/StntupleUtilities.hh"
#include "Stntuple/mod/StntupleBeamManager.hh"
#include "Alignment/CotBeam.hh"
// #include "MyCalDataMaker.hh"

namespace {
  static const char rcsid[] = "$Id: ";
}

//_____________________________________________________________________________
TauEmbedModule::TauEmbedModule (const char* fName, const char* fDescription) 
  : StntupleModule (fName, fDescription)
 , fProcessName     ("processName"     ,this,"PROD")
  // , fPdgCode         ("pdgCode"         ,this,15)
  //
  // quick and dirty fix: test with the electrons
  //
 , fPdgCode         ("pdgCode"         ,this,11)
 , fEmbeddingMode   ("embeddingMode"   ,this,2)
{

  commands()->append( &fProcessName);
  commands()->append( &fPdgCode);
  commands()->append( &fEmbeddingMode);
  
  TStnNode* node;

  std::cout << "Using " << rcsid << " instance " << name() << std::endl;
  fEvent = new TStnEvent();

  fEvent->SetErrorLogger(StntupleModule::fgErrorLogger);

  fEvent->AddDataBlock("ElectronBlock","TStnElectronBlock",node);
  fElectronBlock = (TStnElectronBlock*) node->GetDataBlock();
  fElectronBlock->SetCollName("PROD:CalDriven");

  fEvent->AddDataBlock("TrackBlock","TStnTrackBlock",node);
  fTrackBlock = (TStnTrackBlock*) node->GetDataBlock();
  fTrackBlock->SetCollName("PROD:default");


  fListOfTightElectrons    = new TObjArray();
  fListOfLooseElectrons    = new TObjArray();
  fListOfUnbiasedElectrons = new TObjArray();
  fListOfBiasedElectrons   = new TObjArray();
//-----------------------------------------------------------------------------
// define electron cuts
//-----------------------------------------------------------------------------
  fEleID = new TStnElectronID();
  fRandom = new TRandom();

  //  fCalDataMaker = new MyCalDataMaker();
}

//_____________________________________________________________________________
TauEmbedModule::~TauEmbedModule() {
  // later will delete everything

  delete fRandom;
  delete fEleID;
  //  delete fCalDataMaker;
}
  
//_____________________________________________________________________________
AppResult TauEmbedModule::beginJob (AbsEvent* aJob) {

  fgStntupleFolder->Add(new TNamed("ProcessName"     ,fProcessName.value().c_str()));
  fgStntupleFolder->Add(new TNamed("TauCollName"     ,"CdfTauCollection"));
  fgStntupleFolder->Add(new TNamed("MuonCollName"    ,"ProductionMuons" ));
  fgStntupleFolder->Add(new TNamed("ElectronCollName","CalDriven"));

  fElectronBlock->SetCollName(fProcessName.value().data(),"CalDriven");
  fTrackBlock->SetCollName(fProcessName.value().data(),"default");
//-----------------------------------------------------------------------------
//  set beam history code to use
//-----------------------------------------------------------------------------
  //  fCalDataMaker->initJob();
//-----------------------------------------------------------------------------
//  set beam history code to use
//-----------------------------------------------------------------------------
  StntupleBeamManager* bm = StntupleBeamManager::Instance();
  bm->SetHistoryCode(-1);

  return AppResult::OK; 
}

//_____________________________________________________________________________
AppResult TauEmbedModule::beginRun (AbsEvent* aRun) {
  AbsEnv* env = AbsEnv::instance();

  int  run_num, mc_flag;

  run_num = env->runNumber();
  mc_flag = env->monteFlag();
//-----------------------------------------------------------------------------
// verbose, useDB, useOffSCL : 4.8.4 production didn't use offline SCL, do not
// do it either
//-----------------------------------------------------------------------------
  //  fCalDataMaker->initRun(false,true,false);  

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

  return AppResult::OK; 
}

//_____________________________________________________________________________
AppResult TauEmbedModule::endJob (AbsEvent* aJob) {
  return AppResult::OK; 
}

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

//_____________________________________________________________________________
int TauEmbedModule::MakeListOfElectrons(TStnElectronID* EleID) 
{
  //
  
  TStnElectron* ele;
  int           tight_id_word, loose_id_word;

  StntupleMakerModule::InitTrackBlock   (fTrackBlock);
  StntupleMakerModule::InitElectronBlock(fElectronBlock);

  // this updates the list of high-Pt electrons, where it is cleared?

  StntupleMakerModule::InitElectronBlockLinks(fElectronBlock);

  fListOfTightElectrons->Clear();
  fListOfLooseElectrons->Clear();

  int nele = fElectronBlock->NElectrons();
  for (int i=0; i<nele; i++) {
    ele     = fElectronBlock->Electron(i);
    tight_id_word = EleID->TightIDWord(ele);
    if (tight_id_word == 0) {
      fListOfTightElectrons->Add((TObject*)ele);
    }
    loose_id_word = EleID->LooseIDWord(ele);
    if (loose_id_word == 0) {
      fListOfLooseElectrons->Add((TObject*)ele);
    }
    
  }
  return 0;
}


//_____________________________________________________________________________
int TauEmbedModule::CreateMcBanks(TObjArray* ListOfElectrons) {

  const double ELE_MASS = 0.511e-3;
  const double TAU_MASS = 1.778;
//-----------------------------------------------------------------------------
//  mark old HEPG banks (if any)
//-----------------------------------------------------------------------------
  for (AbsEvent::ConstIterator ib(fCdfEvent,"HEPG_StorableBank"); 
                                                         ib.is_valid(); ib++) {
    ConstHandle<HEPG_StorableBank> h(ib);
    HEPG_StorableBank* hepg = (HEPG_StorableBank*) &(*h);
    hepg->set_description("UNUSED");
  }
  
  int nele = ListOfElectrons->GetEntriesFast();
  if (nele == 0) return 0;
//-----------------------------------------------------------------------------
//  all the electrons are "PRIMARY"
//-----------------------------------------------------------------------------
  HEPG_StorableBank* hepg = new HEPG_StorableBank(1,nele,0);
  hepg->set_firstDisplacedParticleOffset(nele+1);

  for (int i=0; i<nele; i++) {
    TStnElectron* ele = (TStnElectron*) ListOfElectrons->At(i);

    if (ele->Charge() < 0) hepg->set_idhep  (i, fPdgCode.value());
    else                   hepg->set_idhep  (i,-fPdgCode.value());
//-----------------------------------------------------------------------------
//  MC electrons are "stable", taus are not
//-----------------------------------------------------------------------------
    if (fPdgCode.value() == 11) {
      hepg->set_istdhep(i,1);
      hepg->set_Mass(i,ELE_MASS);
    }
    else if (fPdgCode.value() == 15) {
      hepg->set_istdhep(i,2);
      hepg->set_Mass(i,TAU_MASS);
    }

    hepg->set_jmo1hep(i,0);
    hepg->set_jmo2hep(i,0);
    hepg->set_jda1hep(i,0);
    hepg->set_jda2hep(i,0);
    hepg->set_Px(i,ele->Momentum()->Px());
    hepg->set_Py(i,ele->Momentum()->Py());
    hepg->set_Pz(i,ele->Momentum()->Pz());
    hepg->set_E (i,ele->Momentum()->E());
  }
  hepg->set_description("TauEmbed");
  Handle<HEPG_StorableBank> h(hepg);
  fCdfEvent->append(h);
//-----------------------------------------------------------------------------
//  generate bank with the primary vertex position
//  for starters assume that beam axis is at (0,0)
//  don't forget to drop OBSP and OBSV banks
//  also need to "mask" one of the HEPG banks...
//-----------------------------------------------------------------------------
  TStnElectron* e0 = (TStnElectron*) ListOfElectrons->At(0);
  MVTX_StorableBank* mvtx = new MVTX_StorableBank(fCdfEvent,1);
  mvtx->set_description("TauEmbed");
					// retrieve COT beam position

  StntupleBeamManager* bm = StntupleBeamManager::Instance();
  
  CotBeam* cot_beam = bm->GetCotBeam();

  HepPoint3D pos = cot_beam->position(e0->Z0());

  mvtx->rw(1) = pos.x();
  mvtx->rw(2) = pos.y();
  mvtx->rw(3) = e0->Z0();
  mvtx->rw(4) = 0;

  Handle<MVTX_StorableBank> mvtx_handle(mvtx);
  fCdfEvent->append(mvtx_handle);

  return 0;
}

//_____________________________________________________________________________
AppResult TauEmbedModule::event (AbsEvent* Event) {
//-----------------------------------------------------------------------------
//  this module finds high-Pt electrons in W and Z events 
//  removes their track hits and energy depositions from the event 
//  and adds to an event a HEPG bank with the tau-leptons, which parameters
//  correspond to those of the electron
//-----------------------------------------------------------------------------
  int           n_loose_ele, n_tight_ele, irem;
  bool          passed;
  TStnElectron  *e1, *e2, *te1, *e_u, *e_b;
  TLorentzVector mom;

  fCdfEvent = Event;
  fListOfUnbiasedElectrons->Clear();
  fListOfBiasedElectrons->Clear();
//-----------------------------------------------------------------------------
// step 1: create lists of high-Pt electrons - should we be using STNTUPLE
//         infrastructure for that? - may not be a bad idea...
//-----------------------------------------------------------------------------
  MakeListOfElectrons(fEleID);

  n_tight_ele = fListOfTightElectrons->GetEntriesFast();
  n_loose_ele = fListOfLooseElectrons->GetEntriesFast();

  passed      = false;

  if ((n_tight_ele > 0) && (n_loose_ele == 2)) {

				// 2 EM objects, make sure they make Z

    e1 = (TStnElectron*) fListOfLooseElectrons->UncheckedAt(0);
    e2 = (TStnElectron*) fListOfLooseElectrons->UncheckedAt(1);
    if (e1->Charge()+e2->Charge() == 0) {
      mom  = *e1->Momentum();
      mom += *e2->Momentum();

      double mass = mom.M();
      if ( fabs(mass-91.) < 15.) {
//-----------------------------------------------------------------------------
// OK, Z is found, choose the electron to replace
//-----------------------------------------------------------------------------
	passed = true;
	if (n_tight_ele == 1) {
//-----------------------------------------------------------------------------
// 1 tight , 1 loose electron
//-----------------------------------------------------------------------------
	  te1 = (TStnElectron*) fListOfTightElectrons->UncheckedAt(0);
	  if (te1 == e1) {
	    fListOfUnbiasedElectrons->Add(e2);
	    fListOfBiasedElectrons->Add(e1);
	  }
	  else {
	    fListOfUnbiasedElectrons->Add(e1);
	    fListOfBiasedElectrons->Add(e2);
	  }
	}
	else {
//-----------------------------------------------------------------------------
// both electrons are tight, choose randomly
//-----------------------------------------------------------------------------
	  irem = fRandom->Integer(2);
	  e_u = (TStnElectron*) fListOfTightElectrons->UncheckedAt(irem);
	  fListOfUnbiasedElectrons->Add(e_u);
	  e_b = (TStnElectron*) fListOfTightElectrons->UncheckedAt(1-irem);
	  fListOfBiasedElectrons->Add(e_b);
	}
      }
    }
  }
//-----------------------------------------------------------------------------
// make decision
//-----------------------------------------------------------------------------
  setPassed(passed);
  if (! passed)                                             goto END;
//-----------------------------------------------------------------------------
// this ia a Z->ee event, it passes the cuts
// step 2: loop over the electrons and create HEPG bank, 
// replacing e- with tau- and e+ with tau+
// I want to have such a bank in every passed event, even if the bank is empty
//-----------------------------------------------------------------------------
  TObjArray* list;
  if (fEmbeddingMode.value() == 1) {
//-----------------------------------------------------------------------------
// just take the "biased" tight electron out, create W' -> e nu
//-----------------------------------------------------------------------------
    list = fListOfBiasedElectrons;
  }
  else if (fEmbeddingMode.value() == 2) {
//-----------------------------------------------------------------------------
//  replace the unbiased electrons (there is only one!) with the tau's
//  remove both electrons - create W'->tau nu
//-----------------------------------------------------------------------------
    CreateMcBanks(fListOfUnbiasedElectrons);
    list = fListOfLooseElectrons;
  }
  else {
    strcpy((char*)0,"emoe"); // should never happen, bomb if come here
  }
//-----------------------------------------------------------------------------
//  unpack the corresponding data banks into the appropriate structures
//  assume that the track-to-hit links do exist and I can get from the hits
//  back to the D-bank level information
// loop over the electrons again and remove them from the event record
//
//  tracking:
//-----------------------------------------------------------------------------
  RemoveCotHits(list);
  RemoveSvxHits(list);
//-----------------------------------------------------------------------------
// CES: 
//-----------------------------------------------------------------------------
  RemoveCesHits(list);
//-----------------------------------------------------------------------------
// CPR
//-----------------------------------------------------------------------------
  RemoveCprHits(list);
//-----------------------------------------------------------------------------
// calorimeter
//-----------------------------------------------------------------------------
  RemoveCalHits(list);
//-----------------------------------------------------------------------------
// plug strips
//-----------------------------------------------------------------------------
  RemovePesHits(list);
//-----------------------------------------------------------------------------
// remove TOF hits
//-----------------------------------------------------------------------------
  RemoveTofHits(list);

  //  fCdfEvent->listObjects(std::cout);
 END:;
  return AppResult::OK;
}


//_____________________________________________________________________________
int TauEmbedModule::RemoveSvxHits(TObjArray* ListOfElectrons) {
//-----------------------------------------------------------------------------
// step 3a: silicon hits. Have to unpack SVX banks into the structures, where
//          the hits can be deleted from. Find bob's track, mark strips
//          corresponding to the clusters, associated with the electron track
//-----------------------------------------------------------------------------
  TStnElectron* ele;
  int nele = ListOfElectrons->GetEntriesFast();
  for (int i=0; i<nele; i++) {
    ele = (TStnElectron*) ListOfElectrons->UncheckedAt(i);
    CdfEmObject*  bob = ele->GetCdfEmObject();
    const CdfTrack* trk = bob->maxPtTrack().operator->();

    if (trk) {

      std::vector<const SiHit*>::const_iterator it_si = trk->beginSIHits();

      if (AppModule::_verbose.value()) printf(" ------------ electron strips\n");

      while (it_si != trk->endSIHits()) {
	const SiHit* hit = *it_si;
	const SiCluster* cl = hit->getCluster();
    
	SiStrip* strip = (SiStrip*) cl->firstStrip();
	while (strip != cl->lastStrip()) {
	  if (AppModule::_verbose.value() > 0) {
	    printf(" strip: 0x%08x number: %4i pulse = %5i adc = %i5 noise = %5i \n",
		   strip,
		   (*strip).getStripNumber(),
		   (*strip).getPulseHeight(),
		   (*strip).getScaledAdc  (),
		   (*strip).getNoise      ()
		   );
	  }
	  strip->setPulseHeight(-100);
	  if (AppModule::_verbose.value() > 0) {
	    printf(" strip: 0x%08x number: %4i pulse = %5i adc = %i5 noise = %5i \n",
		   strip,
		   (*strip).getStripNumber(),
		   (*strip).getPulseHeight(),
		   (*strip).getScaledAdc  (),
		   (*strip).getNoise      ()
		   );
	  }
	  strip++;
	};
      
	it_si++;
      }
    }
  }

  if (AppModule::_verbose.value()) printf(" ------------ electron strips\n");
//-----------------------------------------------------------------------------
// mark old SIXD and ISLD banks, also stripset...
//-----------------------------------------------------------------------------
  EventRecord::ConstIterator it;
  it = StntupleGetIterator(fCdfEvent,"SIXD_StorableBank","","");
  if (it.is_valid()) {
    ConstHandle<SIXD_StorableBank> h_sixd(it);
    SIXD_StorableBank* sixd = (SIXD_StorableBank*) h_sixd.operator->();
    sixd->set_description("UNUSED");
  }

  it = StntupleGetIterator(fCdfEvent,"ISLD_StorableBank","","");
  if (it.is_valid()) {
    ConstHandle<ISLD_StorableBank> h_isld(it);
    ISLD_StorableBank* isld = (ISLD_StorableBank*) h_isld.operator->();
    isld->set_description("UNUSED");
  }

  it = StntupleGetIterator(fCdfEvent,"StorableRun2SiStripSet","",
			   "GlobalSI_StripSet");

  if (! it.is_valid()) return -1;
  StorableRun2SiStripSet_ch ch(it);
  StorableRun2SiStripSet* ss = (StorableRun2SiStripSet*) ch.operator->();
  ss->set_description("UNUSED");

//    int geom_printed;
//    for (SiStripSet::const_digi_iterator di=(*ch).contents().begin(); 
//                                         di!=(*ch).contents().end(); di++) {
//      SiDigiCode digi=(*di).first;
//      SiHalfLadderCode higi(digi);

//      geom_printed = 0;
//      for (SiStripSet::const_iterator si=(*di).second.begin(); 
//  	 //	 si!=(*di).second.end() && (*di).second.size()>1;
//  	 si!=(*di).second.end(); si++) {

//        if (! geom_printed) {
//  	geom_printed = 1;
//  	printf(" ------ barrel = %2i hladder = %2i lseg = %2i wedge = %2i ",
//  	       digi.getBarrel(),
//  	       digi.getHalfLadder(),
//  	       digi.getLadderSeg(),
//  	       digi.getPhiWedge() );

//  	printf(" layer = %2i ladder = %2i pnside = %2i \n",
//  	       digi.getLayer(),
//  	       digi.getLadder(),
//  	       digi.getPnSide());
//        }
//  //-----------------------------------------------------------------------------
//  // suspect we can print strip info here
//  //-----------------------------------------------------------------------------
//        printf(" strip: 0x%08x number: %4i pulse = %5i adc = %i5 noise = %5i \n",
//               &(*si),
//  	     (*si).getStripNumber(),
//  	     (*si).getPulseHeight(),
//  	     (*si).getScaledAdc  (),
//  	     (*si).getNoise      ()
//  	     );
//      } 
//    }


//-----------------------------------------------------------------------------
// step 3a: create new stripset, copy old one into the new one w/o marked 
//          strips and create SIXD /ISLD bank from the new stripset
//          how do I know which stripset
//-----------------------------------------------------------------------------

				// first step: just print the stripset
  //  ch->print();

  StorableRun2SiStripSet* new_ss = new StorableRun2SiStripSet();
  new_ss->set_description("TauEmbed");

  for (SiStripSet::const_digi_iterator di=(*ch).contents().begin(); 
                                            di!=(*ch).contents().end(); di++) {
    SiDigiCode digi=(*di).first;
    SiHalfLadderCode higi(digi);

    for (SiStripSet::const_iterator si=(*di).second.begin(); 
	                                        si!=(*di).second.end(); si++) {

      if (si->getPulseHeight() != -100) {
	new_ss->accumulate(digi,(*si));
      }
    }
  }

  Handle<StorableRun2SiStripSet> new_ss_h(new_ss);
  new_ss_h->writeBanks(fCdfEvent);
  fCdfEvent->append(new_ss_h);

  return 0;
}


//_____________________________________________________________________________
int TauEmbedModule::RemoveCesHits(TObjArray* ListOfElectrons) {
// remove all the strips belonging to CES clusters associated to the electron

  TCesDataBlock ces_block;

  const vector<ScStrip> strip;
  const vector<ScWire>  wire;

  EventRecord::ConstIterator bank_iter(fCdfEvent,"CESD_StorableBank");
  if(! bank_iter.is_valid() )   return -1;

  ConstHandle<CESD_StorableBank> cesd_bank(bank_iter);

  CESD_StorableBank* cesd = (CESD_StorableBank*) cesd_bank.operator->();
  cesd->set_description("UNUSED");

  StntupleInitCesDataBlock(&ces_block,fCdfEvent,0);
//-----------------------------------------------------------------------------
//  loop over the electrons and mark electron strips as unused
//-----------------------------------------------------------------------------
  int nele = ListOfElectrons->GetEntriesFast();
  for (int iele=0; iele<nele; iele++) {
    TStnElectron* ele = (TStnElectron*) ListOfElectrons->UncheckedAt(iele);
    CdfEmObject*  bob = ele->GetCdfEmObject();

    const CdfEmObject::Ces_CollType& ces_coll = bob->matchingCesClusters("TRACK");
    CdfEmObject::Ces_CollType::const_iterator ces_it;
    CdfEmObject::Ces_CollType::const_iterator last;
    ces_it = ces_coll.contents().begin();
    last   = ces_coll.contents().end();
    
    for ( ; ces_it != last; ++ces_it) {
      CesCluster* ces_cluster = (CesCluster*) (*ces_it).operator->();
      if (_verbose.value()) ces_cluster->print(cout);

      if (ces_cluster->view() == CesCluster::WIRE) {
//-----------------------------------------------------------------------------
//  wire cluster
//-----------------------------------------------------------------------------
	int nwires = ces_cluster->wires().size();
	for (int iw=0; iw<nwires; iw++) {
	  ScWire w = ces_cluster->wires()[iw];
	  ces_block.SetWireData(w.getSide(),w.getModule(),w.getWireNo(),0);
	}
      }
      else {
//-----------------------------------------------------------------------------
//  strip cluster
//-----------------------------------------------------------------------------
	int nstrips = ces_cluster->strips().size();
	for (int is=0; is<nstrips; is++) {
	  ScStrip s = ces_cluster->strips()[is];
	  ces_block.SetStripData(s.getSide(),s.getModule(),s.getStripNo(),0);
	}
      }
    }
  }
//-----------------------------------------------------------------------------
//  create and fill new CESD bank
//-----------------------------------------------------------------------------
  CESD_StorableBank* new_cesd;
  StntupleFillCesd(&ces_block,new_cesd);
  new_cesd->set_description("TauEmbed");

  Handle<CESD_StorableBank> new_handle(new_cesd);
  fCdfEvent->append(new_handle);

  return 0;
}

//_____________________________________________________________________________
int TauEmbedModule::RemoveCprHits(TObjArray* ListOfElectrons) {
  // take out of the event ALL the CPR clusters associated with the electron

  TCprDataBlock cpr_block;

  EventRecord::ConstIterator bank_iter(fCdfEvent,"CPRD_StorableBank");
  if (! bank_iter.is_valid()) return -1;

  ConstHandle<CPRD_StorableBank> cprd_bank(bank_iter);
  
  CPRD_StorableBank* cprd = (CPRD_StorableBank*) cprd_bank.operator->();
  cprd->set_description("UNUSED");

  StntupleInitCprDataBlock(&cpr_block,fCdfEvent,0);
//-----------------------------------------------------------------------------
//  loop over the electrons and set electron wires to zero
//-----------------------------------------------------------------------------
  vector<CprWire> wire;

  int nele = ListOfElectrons->GetEntriesFast();
  for (int i=0; i<nele; i++) {
    TStnElectron* ele = (TStnElectron*) ListOfElectrons->UncheckedAt(i);
    CdfEmObject*  bob = ele->GetCdfEmObject();
    //    const CdfTrack* trk = bob->maxPtTrack().operator->();

    const IndexViewVector<CprClusterColl>& cpr_list=bob->matchingCprClusters();
    vector< ValueVectorLink<CprClusterColl> >::const_iterator cprCluster;

    for (cprCluster  = cpr_list.contents().begin();
	 cprCluster != cpr_list.contents().end  (); ++cprCluster) {
//-----------------------------------------------------------------------------
// loop over the wires included into electron cluster and set their pulse
// hights to 0 in the CPR data block
//-----------------------------------------------------------------------------
      if (_verbose.value()) {
	(*cprCluster)->print(cout);
      }
      wire = (*cprCluster)->wires();
      int nwires = wire.size();
      for (int iw=0; iw<nwires; iw++) {
	CprWire w = wire[iw];
	cpr_block.SetWireData(w.getSide(),w.getModule(),w.getWireNo(),0);
      }
    }
  }
//-----------------------------------------------------------------------------
//  create new CPRD bank
//-----------------------------------------------------------------------------
  CPRD_StorableBank* new_bank;
  StntupleFillCprd(&cpr_block,new_bank);
  new_bank->set_description("TauEmbed");

  Handle<CPRD_StorableBank> new_handle(new_bank);
  fCdfEvent->append(new_handle);

  return 0;
}

//_____________________________________________________________________________
int TauEmbedModule::RemovePesHits(TObjArray* ListOfElectrons) {
  // 

  EventRecord::ConstIterator bank_iter(fCdfEvent,"PESD_StorableBank");
  if (! bank_iter.is_valid()) return -1;
  ConstHandle<PESD_StorableBank> pesd_bank(bank_iter);
  
  PESD_StorableBank* pesd = (PESD_StorableBank*) pesd_bank.operator->();
  pesd->set_description("PRIMARY");

  return 0;
}

//_____________________________________________________________________________
int TauEmbedModule::RemoveCalHits(TObjArray* ListOfElectrons) {
  // create new version of CalData not containing electron hits...

  int   ieta, iphi;
  CalData::TowerData*            tower;

  CalData_ch caldata_ch;
  CalData::Error result = CalData::find(caldata_ch);
  if (result != CalData::OK) return -1;
//-----------------------------------------------------------------------------
//  mark old CalData as unused
//-----------------------------------------------------------------------------
  CalData* caldata = (CalData*) caldata_ch.operator->();
  caldata->set_description("UNUSED");


  if (AppModule::_verbose.value()) {
    caldata->print(cout);
  }
//-----------------------------------------------------------------------------
// create new CalData object
//-----------------------------------------------------------------------------
  CalData* new_caldata = new CalData(*caldata);
  new_caldata->set_description("TauEmbed");
  if (AppModule::_verbose.value()) {
    new_caldata->print(cout);
  }
//-----------------------------------------------------------------------------
//  loop over the electrons and delete their towers from new CalData
//-----------------------------------------------------------------------------
  tower = new_caldata->getData();

  int nele = ListOfElectrons->GetEntriesFast();
  for (int i=0; i<nele; i++) {
    TStnElectron* ele = (TStnElectron*) ListOfElectrons->UncheckedAt(i);
    CdfEmObject*  bob = ele->GetCdfEmObject();
    EmCluster*    cl  = (EmCluster*) bob->getEmCluster();

    const std::vector<TowerKey>& towlist = cl->towerIds();
    int ntowers = towlist.size();
    for (int itower=0; itower<ntowers; itower++) {
      const TowerKey tkey =  towlist[itower];

      ieta = tkey.iEta();
      iphi = tkey.iPhi();

      delete (*tower)[ieta][iphi];
      (*tower)[ieta][iphi] = 0;
    }
  }

  if (AppModule::_verbose.value()) {
    new_caldata->print(std::cout);
  }

  Handle<CalData> h_caldata(new_caldata);
  fCdfEvent->append(h_caldata);

  return 0;
}


//_____________________________________________________________________________
int TauEmbedModule::RemoveTofHits(TObjArray* ListOfElectrons) {

  EventRecord::ConstIterator bank_iter(fCdfEvent,"TOFD_StorableBank");
  if (! bank_iter.is_valid()) return -1;
  ConstHandle<TOFD_StorableBank> tofd_bank(bank_iter);
  
  TOFD_StorableBank* tofd = (TOFD_StorableBank*) tofd_bank.operator->();
  tofd->set_description("PRIMARY");

  return 0;
}

//_____________________________________________________________________________
int TauEmbedModule::RemoveImuHits(TObjArray* ListOfElectrons) {

  EventRecord::ConstIterator bank_iter(fCdfEvent,"IMUE_StorableBank");
  if (! bank_iter.is_valid()) return -1;

  ConstHandle<IMUE_StorableBank> imue_bank(bank_iter);
  
  IMUE_StorableBank* imue = (IMUE_StorableBank*) imue_bank.operator->();
  imue->set_description("PRIMARY");

  return 0;
}

//_____________________________________________________________________________
int TauEmbedModule::RemoveCotHits(TObjArray* ListOfElectrons) {
//-----------------------------------------------------------------------------
// step 3b: mark COT hits, associated with the electron track
//-----------------------------------------------------------------------------
  fCT_HitSet= CT_HitSet::instance();

  CT_HitSet* fcc= CT_HitSet::instance_raw();

  int nele = ListOfElectrons->GetEntriesFast();
  for (int i=0; i<nele; i++) {
    TStnElectron* ele = (TStnElectron*) ListOfElectrons->UncheckedAt(i);
    CdfEmObject*  bob = ele->GetCdfEmObject();

    const CdfTrack* trk = bob->maxPtTrack().operator->();

    if (trk) {
  //  fCT_HitSet->loadHitsFromCOTD();
//-----------------------------------------------------------------------------
// track-to-hit links are saved, however CT_HitSet is recostructed on the fly
//-----------------------------------------------------------------------------
      CT_HitLink* it = (CT_HitLink*) trk->beginCTHitBlock();
      while (it != trk->endCTHitBlock()) {
					// this marks hits in the instance
	CT_Hit* h = (CT_Hit*) &it->hit();
	h->set_packed(-1); 
	it++;
      }
    }
  }

  if (AppModule::_verbose.value()) {
    fCT_HitSet->print();
  }
//-----------------------------------------------------------------------------
// pack all the COT hits except marked ones in the new COTD bank
// need to pack a RAW hitset!
//-----------------------------------------------------------------------------
  TDC_D_BankMaker maker;
  std::vector<int> obsp_numbers;
  std::vector<int4> bankVector;
    
  obsp_numbers.clear();

  int nh = fCT_HitSet->numHits();
  const CT_Configuration& cfg = CT_Configuration::reference();
  

  CT_CellCursor ct_cc;

  for (int sl=0; sl < cfg.numSuperLayers(); ++sl) {
    for (CT_CellCursor c=fCT_HitSet->cellCursor(sl,0); 
	 c.numSteps() < cfg.numCells(sl); ++c) {
	
      ct_cc                   = fcc->cellCursor(sl,c.numSteps());
      CT_HitSet::Iterator ihh = ct_cc.beginHit();

      for (CT_HitSet::Iterator ih=c.beginHit(); ih != c.endHit(); ++ih) {
	if (ih->packed() != -1) {
	  CT_WireCode wc(sl,c.cellNumber(),ih->wire());
	  int channelID = COTD_StorableBank::getChannelID(wc);
	  int moduleID = COTD_StorableBank::getModuleID(wc);
	  TDC_D_BankMaker::TDCFormatAddress address = 
	    COTD_StorableBank::getFormatAddress(wc);
	  
//  	  maker.addHit(address,channelID, moduleID, ih->tdc(),
//  		       obsp_numbers, ih->width());

	  maker.addHit(address,channelID, moduleID, ihh->tdc(),
		       obsp_numbers, ihh->width());
	}
	ihh++;
      }
    }
  }
//-----------------------------------------------------------------------------
// housekeeping in the event record: mark old COT hit banks to deletion
//-----------------------------------------------------------------------------
  COTD_StorableBank* cotd;
  EventRecord::ConstIterator icotd(fCdfEvent,"COTD_StorableBank");
  if (icotd.is_valid()) {
    ConstHandle<COTD_StorableBank> hcotd(icotd);
    cotd = (COTD_StorableBank*) hcotd.operator->();
    cotd->set_description("UNUSED");
  }

  EventRecord::ConstIterator icotq(fCdfEvent,"COTQ");
  if (icotq.is_valid()) {
    ConstHandle<COTQ> hcotq(icotq);
    COTQ* cotq = (COTQ*) hcotq.operator->();
    cotq->set_description("UNUSED");
  }
//-----------------------------------------------------------------------------
// Done adding Hits, now make the bank and collection
//-----------------------------------------------------------------------------
  maker.formVectors(bankVector,obsp_numbers);
  
  cotd = new COTD_StorableBank(TRY_Bank_Number(2),
			       bankVector.size(),
			       &(bankVector[0]));
  cotd->set_description("TauEmbed");

  if (AppModule::_verbose.value()) {
    cotd->print();
  }

  Handle<COTD_StorableBank> COTD_h(cotd);
  fCdfEvent->append(COTD_h);

  return 0;
}