//-----------------------------------------------------------------------------
// June 6 2001 A.Dominguez: initialization of the STNTUPLE Geant/silicon hit block
//-----------------------------------------------------------------------------
#include "AbsEnv/AbsEnv.hh"
#include "Edm/ConstHandle.hh"

#include "Stntuple/obj/TSiGeantIsectBlock.hh"
#include "Stntuple/mod/InitStntupleDataBlocks.hh"
#include "Stntuple/obj/TObspBlock.hh"
#include "Stntuple/obj/TStnEvent.hh"

#include "SimulationObjects/PropagatedSiParticle.hh"
#include "SimulationObjects/PropagatedSiParticleColl.hh"

#include "TLorentzVector.h"

#include "CLHEP/Geometry/Point3D.h"
#include <cmath>

#include "ErrorLogger_i/gERRLOG.hh"

//_____________________________________________________________________________
Int_t StntupleInitSiGeantIsect(TStnSiGeantIsect& isect, 
			       const PropagatedSiParticle& prop) {
  // initialize STNTUPLE Geant/Si intersects starting from CDF Prop.Parts.
  
  isect.SetUniqueID(prop.id()-1); // Silly fortran numbering
  isect.fDigiCode.fDigiCode = prop.value();
  isect.fEntryMomentum = TLorentzVector(prop.entryFourMomentum().x(),
					prop.entryFourMomentum().y(),
					prop.entryFourMomentum().z(),
					prop.entryFourMomentum().t());
  isect.fExitMomentum = TLorentzVector(prop.exitFourMomentum().x(),
				       prop.exitFourMomentum().y(),
				       prop.exitFourMomentum().z(),
				       prop.exitFourMomentum().t());
  isect.fEntry = TVector3(prop.entryPoint().x(),
			  prop.entryPoint().y(),
			  prop.entryPoint().z());
  isect.fExit = TVector3(prop.exitPoint().x(),
			 prop.exitPoint().y(),
			 prop.exitPoint().z());
  isect.fDE = prop.dE();
  isect.fRadLen = prop.radLen();
  for (int i = 0; i<8; i++){
    if ( (prop.flag() >> i) & 0x1) isect.fActiveRegion.SetBit(i);
  }
  isect.fPhiStrip = prop.stripNum()[0];
  isect.fZStrip   = prop.stripNum()[1];
  isect.fInitPhi  = prop.modelInitStrip()[0];
  isect.fFinalPhi = prop.modelFinalStrip()[0];
  isect.fInitZ  = prop.modelInitStrip()[1];
  isect.fFinalZ = prop.modelFinalStrip()[1];

  return 0;
}

//_____________________________________________________________________________
Int_t StntupleInitSiGeantIsectBlock(TStnDataBlock* Block, 
				    AbsEvent*      Event, 
				    int            Mode) {
  // fill silicon Geant intersection block

  int ev_number, rn_number;
					// check if block has already been 
					// initialized

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

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

  TSiGeantIsectBlock* data = (TSiGeantIsectBlock*) Block;
  data->Clear();

  EventRecord::ConstIterator it(Event,"PropagatedSiParticleColl");
  if (! it.is_valid()) { // No hit list.  Go home
    return -1;
  }
  // Assume that we want the first Coll in the event record
  ConstHandle<PropagatedSiParticleColl> ph(it);
  if (ph->contents().empty()) { // Hit list is empty, go home.
    return -1;
  }
  
  // Loop over all the intersections and make TStnSiGeantIsects out of them
  for (int i=0; i<ph->contents().size(); i++){
    
    TStnSiGeantIsect* isect = data->NewSiGeantIsect(); // new space for 1 TStnSiGeantIsect in TClonesArray
    StntupleInitSiGeantIsect( *isect, ph->contents()[i] ); // Fill branch

  }

  data->f_EventNumber = ev_number;
  data->f_RunNumber   = rn_number;

  return 0;
}


// After geant intersection block and obsp blocks have been filled,
// make links between them
Int_t StntupleSiGeantIsectBlockLinks(TStnDataBlock* Block, 
				     AbsEvent*      Event, 
				     int            Mode) {
  TSiGeantIsectBlock* data = (TSiGeantIsectBlock*) Block;
  if (data->NSiGeantIsect()<1) return 0; // no geant intersections
  
  TObspBlock* odata;
  odata = (TObspBlock* ) Block->GetEvent()->GetDataBlock("ObspBlock");
  if (! odata) return 0; // No OBSP block in event, bail out gracefully
  
  if (odata->NParticles()<0) return -1;  // What kind of event is this?

  data->IsectLinkObsp()->Init( odata->NParticles() );
  int nobsp=0;

  // Geant isect block is stored in increasing obsp numbers in
  // contiguous blocks of the same obsp number.  This allows for fast
  // mapping in one loop thru the data.  (This also assumes that the
  // obsp block is a superset of obsp numbers of the obsp numbers in
  // the intersection block.)
  int ig=0;
  TStnSiGeantIsect *g = data->SiGeantIsect(ig);
  for (int i=0; i<odata->NParticles(); i++) {
    int iobsp = odata->Particle(i)->Number();
    data->IsectLinkObsp()->InitLinks(iobsp);
    
    int gobsp = g->GetUniqueID();
    if ( gobsp == iobsp ) { // Loop thru geant isect block
      ++nobsp;
      while (ig<data->NSiGeantIsect() && g->GetUniqueID()==iobsp) {
	data->IsectLinkObsp()->Add( i, ig );
	++ig;
	if (ig<data->NSiGeantIsect()) g = data->SiGeantIsect(ig);
      }
    }
  }

  if (ig != data->NSiGeantIsect()) { // Something went wrong
    data->IsectLinkObsp()->Clear(); // with our assumptions!
    return -1;
  }

  return 0;
}