/***************************************************************************
                         InDetSoftwareTutorial.cxx  -  Description
                             -------------------
    begin   : 10-05-2005
    authors : Edward Moyse(CERN PH-ATC)
    email   : Edward.Moyse@cern.ch
 ***************************************************************************/

#include "GaudiKernel/MsgStream.h"

#include "StoreGate/StoreGateSvc.h"
#include "DataModel/DataVector.h"
#include "GaudiKernel/IHistogramSvc.h"
//#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/SmartDataPtr.h"  //NTupleFilePtr
                                                                                
// Athena Histogramming Stuff
#include "AIDA/IHistogram1D.h"
//#include "AIDA/IHistogram2D.h"

#include "CLHEP/Matrix/Vector.h"

#include "TrkEventPrimitives/ParamDefs.h"
#include "TrkTrack/Track.h"
#include "TrkTrack/TrackCollection.h"
//#include "TrkParameters/Perigee.h"
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "TrkSurfaces/StraightLineSurface.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkParameters/MeasuredTrackParameters.h"

#include "InDetPrepRawData/SCT_Cluster.h"
#include "InDetPrepRawData/SiCluster.h"
#include "TrkEventUtils/SiWidth.h"
// Geometry Stuff
#include "Identifier/Identifier.h"
#include "AtlasDetDescr/AtlasDetectorID.h"
#include "IdDictDetDescr/IdDictManager.h"
//#include "InDetReadoutGeometry/PixelDetectorManager.h"
#include "InDetReadoutGeometry/SCT_DetectorManager.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "InDetReadoutGeometry/SiDetectorElementCollection.h"
//#include "DetDescrConditions/AlignableTransform.h"

#include "InDetSoftwareTutorial/InDetSoftwareTutorial.h"

namespace InDet
{
  
InDetSoftwareTutorial::InDetSoftwareTutorial(const std::string &n, ISvcLocator *pSvcLoc) 
  : 
  Algorithm(n, pSvcLoc),
  m_fitTool(0),
  m_maxD0(0.25),
  m_storeGate(0),
  m_tracksName("Tracks") // set sensible defaults, so you don't rely on jobOptions
{
  declareProperty("TracksName",m_tracksName);
  declareProperty("maxD0", m_maxD0);
}

InDetSoftwareTutorial::~InDetSoftwareTutorial() {}

StatusCode InDetSoftwareTutorial::initialize()
{
  MsgStream log(msgSvc(), name());
  StatusCode sc;

  // Retrieve StoreGate service 
  sc = service("StoreGateSvc", m_storeGate);
  if (sc.isFailure())
  {
    log << MSG::ERROR << "Could not find StoreGateSvc." << endreq;
    return sc;
  }
    m_idHelper = new AtlasDetectorID;

  // get DetectorStore service
  sc = service("DetectorStore", m_detStore);                                    
  if (sc.isFailure()) {
       log << MSG::FATAL << "Could not find Detector service." << endreq;
  return sc;
  }

    // Get the dictionary manager from the detector store
    const IdDictManager*  idDictMgr = 0;
    sc = m_detStore->retrieve(idDictMgr, "IdDict");
    if (sc.isFailure()) {
      log << MSG::ERROR << "Could not get IdDictManager !" << endreq;
      return sc;
    }
    // Initialize the helper with the dictionary information.
    if (idDictMgr) {
      if (idDictMgr->initializeHelper(*m_idHelper)) { // Returns 1 if there is a problem
        log << MSG::ERROR << "Unable to initialize ID helper." << endreq;
        return StatusCode::FAILURE;
      }
    } else {
      log << MSG::ERROR << "IdDictManager pointer is zero. "
           << "Unable to initialize id helper. " << endreq;
      return StatusCode::FAILURE;
    }
    sc = m_detStore->retrieve(m_SCT_Manager, "SCT");
    if (sc.isFailure()) {
      log << MSG::ERROR << "Could not get SCT_Manager !" << endreq;
      return sc;
    }
    m_histo1D_Barrel_Residual = histoSvc()->book("/stat/residuals/Barrel/","Barrel Residual","SCT Barrel Residual", 300, -0.15, 0.15);
    m_histo1D_Barrel_0_Residual = histoSvc()->book("/stat/residuals/Layer0/","Barrel Layer 0 Residual","SCT Barrel Layer 0 Residual", 300, -0.15, 0.15);
    m_histo1D_Barrel_1_Residual = histoSvc()->book("/stat/residuals/Layer1/","Barrel Layer 1 Residual","SCT Barrel Layer 1 Residual", 300, -0.15, 0.15);
    m_histo1D_Barrel_2_Residual = histoSvc()->book("/stat/residuals/Layer2/","Barrel Layer 2 Residual","SCT Barrel Layer 2 Residual", 300, -0.15, 0.15);
    m_histo1D_Barrel_3_Residual = histoSvc()->book("/stat/residuals/Layer3/","Barrel Layer 3 Residual","SCT Barrel Layer 3 Residual", 300, -0.15, 0.15);                                                                                
 
  log << MSG::INFO << "Initialization successful" << endreq;

  return StatusCode::SUCCESS;
}

StatusCode InDetSoftwareTutorial::finalize()
{
  return StatusCode::SUCCESS;
}

StatusCode InDetSoftwareTutorial::execute()
{
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "Execution track retrieve" << endreq;
  // Retrieve tracks from StoreGate 
  const TrackCollection* tracks;
  StatusCode sc = m_storeGate->retrieve( tracks, m_tracksName);
  if (sc.isFailure())
  {
    log << MSG::DEBUG << "Could not find DataVector<Trk::Track> " << m_tracksName << " in StoreGate." << endreq;
    return StatusCode::SUCCESS; 
    // it may seem strange to return SUCCESS when in fact we failed to find any tracks
    // but the point is that returning FAILURE can cause Athena to stop 
    // (which should only happen with catasprohic errors)
  }
  log << MSG::DEBUG << "Retrieved " << tracks->size() << " tracks " << endreq;
//    std::vector<const Trk::Track*> selectedTracks;
    DataVector<Trk::Track>::const_iterator trackIt    = tracks->begin();
    DataVector<Trk::Track>::const_iterator trackItEnd = tracks->end();    
    for ( ; trackIt!= trackItEnd; trackIt++)
    {
      // loop over tracks
      //DetectorDescription/GeoModel/HitDisplay/src/TrackSystem.cxx
//InDetAlignGenAlgs/src/AlignARes.cxx but error inebitable
//taka      const Trk::MeasuredPerigee *perigee = dynamic_cast<const Trk::MeasuredPerigee *>((**trackIt).perigeeParameters());
//      const Trk::Perigee *perigee=**trackIt->perigeeParameters();
//taka      const HepVector& perigeeParms = perigee->parameters();
//taka      double d0 = perigeeParms[Trk::d0];
 
//taka       log << MSG::DEBUG << "d0= " << d0 <<  endreq;

       //log << MSG::DEBUG << "Begin HitLoop over " << trackIt->trackStateOnSurfaces()->size() << " TrkSurfaces" << endreq;
       HitLoop(*trackIt);
    } // end of Track Loop

      
  return StatusCode::SUCCESS;
}
void InDetSoftwareTutorial::HitLoop(const Trk::Track* trackIt) {
      MsgStream log(msgSvc(),name());

      int nSurface = 0;
      int nHitsEndCap =0;
      int nHitsBarrel = 0;
       for (std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsos=trackIt->trackStateOnSurfaces()->begin();tsos!=trackIt->trackStateOnSurfaces()->end();
++tsos) {
      ++nSurface;
      const Trk::RIO_OnTrack* hitOnTrack = dynamic_cast <const Trk::RIO_OnTrack*>((*tsos)->measurementOnTrack());
//10.0      const Trk::RIO_OnTrack* hitOnTrack = ((*tsos)->rioOnTrack());
      //log << MSG::DEBUG << "RIO_OnTrack read out" << endreq;
      // TSOS must contain an RIO_OnTrack
      if (hitOnTrack != 0) {

        const Identifier& surfaceID = hitOnTrack->identify();
      //  log << MSG::DEBUG << "Hit on " << m_idHelper->show_to_string(surfaceID,0,'/') << endreq;
        const Identifier& SCT_ModuleID = m_SCT_Manager->getIdHelper()->module_id(surfaceID); //from wafer id to module id

        // SCT-Hits only
        if (m_idHelper->is_sct(surfaceID)) {

          const Trk::TrackParameters* trackParams = ((*tsos)->trackParameters());

          // TSOS must also contain TrackParameters
          if (trackParams!=0) {

  //          log << MSG::DEBUG << "TrackParameter:" << endreq;
  //          log << MSG::DEBUG << (*trackParams) << endreq;
  //          log << MSG::DEBUG << "Rio_OnTrack:" << endreq;
  //          log << MSG::DEBUG << (*hitOnTrack) << endreq;
  //          log << MSG::DEBUG << "Local x Position of strip " << m_idHelper->show_to_string(surfaceID,0,'/') << ": " <<
  //                                m_SCT_Manager->getDetectorElement(surfaceID)->localPositionOfCell(surfaceID).x() << endreq;
  //          log << MSG::DEBUG << "Local y Position of strip " << m_idHelper->show_to_string(surfaceID,0,'/') << ": " <<
 //                                 m_SCT_Manager->getDetectorElement(surfaceID)->localPositionOfCell(surfaceID).y() << endreq;
// need cluster info
// InDetRecEvent/InDetPrepRawData/SCT_Cluster.h:Si_Cluster.h
     /** return width class reference */
//         virtual const Trk::SiWidth&  width()  const
// example : iPatTrackFollower/LayerPrediction.h
//class SiCluster :   public Trk::PrepRawData{ why ?
             const SiCluster* siCluster=dynamic_cast<const SiCluster *>(hitOnTrack->prepRawData());
             if(siCluster){
             const SCT_Cluster* sctCluster=dynamic_cast<const SCT_Cluster *>(siCluster);
             Trk::SiWidth clusterWidth=sctCluster->width();
             Hep2Vector vec=clusterWidth.colRow();
            log << MSG::DEBUG << "clusterWidth.colRow().x():" << vec.x()<<"y():"<< vec.y() << endreq;
            }
	     //wrong double residual=hitOnTrack->localParameters()[0]-trackParams->parameters()[0];
//             double residualFirstStrip=m_SCT_Manager->getDetectorElement(surfaceID)->localPositionOfCell(surfaceID).x()-trackParams->parameters()[0];
  //           log << MSG::DEBUG << residual2 << endreq;
//             log << MSG::DEBUG << residual << endreq;
              double residual=m_SCT_Manager->getDetectorElement(surfaceID)->localPositionOfCell(surfaceID).x()-trackParams->parameters()[0];
     // cutted  and saved in file cutted.txt
            if (m_SCT_Manager->getIdHelper()->is_barrel(surfaceID)) {
              ++nHitsBarrel;
              m_histo1D_Barrel_Residual->fill(residual);
              if (m_SCT_Manager->getIdHelper()->layer_disk(surfaceID) == 0) { // Barrel Layer 0
                m_histo1D_Barrel_0_Residual->fill(residual);
              } else if (m_SCT_Manager->getIdHelper()->layer_disk(surfaceID) == 1) { // Barrel Layer 1
                m_histo1D_Barrel_1_Residual->fill(residual);
              } else if (m_SCT_Manager->getIdHelper()->layer_disk(surfaceID) == 2) { // Barrel Layer 2
                m_histo1D_Barrel_2_Residual->fill(residual);
              } else if (m_SCT_Manager->getIdHelper()->layer_disk(surfaceID) == 3) { // Barrel Layer 3
                m_histo1D_Barrel_3_Residual->fill(residual);
              }
            } // end of "is Barrel"
                // EndCap Loop
               else {
                ++nHitsEndCap;
            }
          } else{
            log << MSG::INFO << "No track parameter information found for Surface " << nSurface << endreq;
          } // end of do trackParameters exist?

        } // end of "is_SCT"

      } // end of "is RIO"
     } // Hit Loop

    log << MSG::DEBUG << nHitsBarrel << " Barrel SCT Hits. " << nHitsEndCap << " Forward SCT Hits." << endreq;
}

}//end of namespace
