FPCCDSiliconTracking_MarlinTrk.h

#ifndef FPCCDSILICONTRACKING_MarlinTrk_H
#define FPCCDSILICONTRACKING_MarlinTrk_H 1

#include "marlin/Processor.h"
#include <marlin/Global.h>
#include "lcio.h"
#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <cmath>
#include <IMPL/TrackImpl.h>
#include "ClusterExtended.h"
#include "TrackExtended.h"
#include "TrackerHitExtended.h"
#include "HelixClass_double.h"
#include <EVENT/SimTrackerHit.h>
#include <EVENT/MCParticle.h>
#include <IMPL/TrackImpl.h>

#include "MarlinTrk/IMarlinTrack.h"

#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>

#include "TTree.h"
#include "TNtupleD.h"
#include "TFile.h"
#include "TVector3.h"
#include "TStopwatch.h"
#include "./MCPMap.h"
#include "../src/GetPurity.h" //by Mori

namespace DiagnosticsHistograms {
  class Histograms ;
}

namespace DiagnosticsHistograms2D {
  class Histograms ;
}

using namespace lcio ;
using namespace marlin ;

namespace  gear{
  class GearMgr ;
}

namespace MarlinTrk {
  class HelixFit;
  class IMarlinTrkSystem ;
}

namespace UTIL{
  class LCRelationNavigator ;
}


class FPCCDSiliconTracking_MarlinTrk : public Processor {
  
public:
  
  virtual Processor*  newProcessor() { return new FPCCDSiliconTracking_MarlinTrk ; }
  
  
  FPCCDSiliconTracking_MarlinTrk() ;
  
  virtual void init() ;
  
  virtual void processRunHeader( LCRunHeader* run ) ;
  
  virtual void processEvent( LCEvent * evt ) ; 
  
  
  virtual void check( LCEvent * evt ) ; 
  
  
  virtual void end() ;
  
  
protected:
  
  int _nRun ;
  int _nEvt ;
  EVENT::LCEvent* _current_event;
  
  int _nDivisionsInPhi;
  int _nDivisionsInTheta;
  int _nLayers;
  
  MarlinTrk::HelixFit* _fastfitter;
  
  MarlinTrk::IMarlinTrkSystem* _trksystem ;
  bool _runMarlinTrkDiagnostics;
  std::string _MarlinTrkDiagnosticsName;
  
  bool _MSOn, _ElossOn, _SmoothOn ;
  
  float _initialTrackError_d0;
  float _initialTrackError_phi0;
  float _initialTrackError_omega;
  float _initialTrackError_z0;
  float _initialTrackError_tanL;
  
  double _maxChi2PerHit;  
  double _maxChi2PerHit2nd;  
  
  bool  _UseEventDisplay;
  std::vector<int> _colours;  
  
  void drawEvent();
  
  
  // histogram member variables
  
  bool  _createDiagnosticsHistograms;
  DiagnosticsHistograms::Histograms* _histos ;

  
  int _ntriplets, _ntriplets_good, _ntriplets_2MCP, _ntriplets_3MCP, _ntriplets_1MCP_Bad, _ntriplets_bad;
  
  
  LCCollection* GetCollection(  LCEvent * evt, std::string colName ) ;
  
  LCRelationNavigator* GetRelations( LCEvent * evt, std::string RelName ) ;
  
  std::string  _colNameMCParticles;
  
  
  struct compare_TrackExtended{
    // n.b.: a and b should be TrackExtended const *, but the getters are not const :-(
    bool operator()(TrackExtended *a, TrackExtended *b) const {
      if ( a == b ) return false;
      return (a->getChi2()/a->getNDF() < b->getChi2()/b->getNDF() );
    }
  };
  
  
  std::string _VTXHitCollection;
  std::string _FTDPixelHitCollection;
  std::string _FTDSpacePointCollection;
  std::string _SITHitCollection;
  std::string _siTrkCollection;
  
  std::vector< LCCollection* > _colTrackerHits;
  std::map< LCCollection*, std::string > _colNamesTrackerHits;
  
  std::vector<TrackerHitExtendedVec> _sectors;
  std::vector<TrackerHitExtendedVec> _sectorsFTD;
  
  class TracksWithNHitsContainer {
  public:
    void clear();
    
    inline void resize(size_t maxHits) {
      _tracksNHits.resize(maxHits-2);
      _maxIndex=(maxHits-3);
    }
    
    // Sort all track vectors according to chi2/nDof
    //      void sort();
    
    inline TrackExtendedVec & getTracksWithNHitsVec( size_t nHits ) {
      //return _tracksNHits[ std::min(nHits-3, _maxIndex) ];
      // for debugging: with boundary check
      return _tracksNHits.at(std::min(nHits-3, _maxIndex));
    }
    
  protected:
    std::vector< TrackExtendedVec > _tracksNHits;
    size_t _maxIndex; 
  };
  
  TracksWithNHitsContainer _tracksWithNHitsContainer;
  
  int InitialiseVTX(LCEvent * evt);
  int InitialiseFTD(LCEvent * evt);
  void ProcessOneSector(int iSectorPhi, int iSectorTheta);
  void ProcessOneSectorVer2(int iSectorPhi, int iSectorTheta);
  void CleanUp();
  TrackExtended * TestTriplet(TrackerHitExtended * outerHit, 
                              TrackerHitExtended * middleHit,
                              TrackerHitExtended * innerHit,
                              HelixClass_double & helix,
                              int omegamode);
  
  int BuildTrack_KalFit(TrackerHitExtended * outerHit, 
                 TrackerHitExtended * middleHit,
                 TrackerHitExtended * innerHit,
                 HelixClass_double & helix, 
                 int innerlayer,
                 TrackExtended * trackAR);

  int _useBuildTrackForHighPt;
  double _cosThetaRangeForBuildTrackForHighPt;
  double _phiRangeForBuildTrackForHighPt;
  void getPhiThetaRegionForHighPt(int* boundaries,TrackExtended* trackAR);

  void Sorting( TrackExtendedVec & trackVec);
  void CreateTrack(TrackExtended * trackAR );
  void AttachRemainingVTXHitsSlow();
  void AttachRemainingFTDHitsSlow();
  void AttachRemainingVTXHitsFast();
  void AttachRemainingFTDHitsFast();
  void AttachRemainingVTXHitsVeryFast();
  void TrackingInFTD();
  int BuildTrackFTD(TrackExtended* trackAR, int* nLR, int iS);
  int AttachHitToTrack(TrackExtended * trackAR, TrackerHitExtended * hit, int iopt);
  int AttachHitToTrack_KalFit(TrackExtended * trackAR, TrackerHitExtended * hit);
  
  void FinalRefit(LCCollectionVec* trk_col, LCCollectionVec* rel_col);
  
  float _bField;
  float _chi2WRPhiTriplet;
  float _chi2WRPhiQuartet;
  float _chi2WRPhiSeptet;
  float _chi2WZTriplet;
  float _chi2WZQuartet;
  float _chi2WZSeptet;
  float _minDistCutAttachForFTD;
  float _minDistCutAttachForVXD;
  int _minimalLayerToAttach;
  int _minMissAddition;
  
  // two pi is not a constant in cmath. Calculate it, once!
  static const double TWOPI;
  
  double _dPhi;
  double _dTheta;
  double _dPhiFTD;
  

  
  std::vector<int> _Combinations;
  std::vector<int> _CombinationsFTD;
  
  float _resolutionRPhiVTX;
  float _resolutionZVTX;
  
  float _resolutionRPhiFTD;
  float _resolutionZFTD;
  
  float _resolutionRPhiSIT;
  float _resolutionZSIT;
  
  float _phiCutForMerging;
  float _tanlambdaCutForMerging;
  float _angleCutForMerging;
  //float _angleCutForMerging_highPt;
  //float _angleCutForMerging_lowPt;
  
  int _print;
  int _checkForDelta;
  float _minDistToDelta;
  
  float _distRPhi;
  float _distZ;
  float _chi2FitCut;
  float _chi2FitCut_lowPt;
  
  
  TrackExtendedVec _trackImplVec;
  
  
  float _cutOnD0, _cutOnZ0, _cutOnOmegaVXD, _cutOnOmegaFTD;
  double _cutOnPtVXD,_cutOnPtFTD;
  
  int _minimalHits;
  int _nHitsChi2;
  int _attachVXD;
  int _attachFTD;
  
  int _max_hits_per_sector;
  
  int _nTotalVTXHits,_nTotalFTDHits,_nTotalSITHits;
  int _useSIT;
  int _useFTD;
  
  
  //  int _createMap;
  
  UTIL::BitField64* _encoder;
  int getDetectorID(TrackerHit* hit) { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::subdet]; }
  int getSideID(TrackerHit* hit)     { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::side]; };
  int getLayerID(TrackerHit* hit)    { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::layer]; };
  int getModuleID(TrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::module]; };
  int getSensorID(TrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::sensor]; };

  int getDetectorID(SimTrackerHit* hit) { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::subdet]; }
  int getSideID(SimTrackerHit* hit)     { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::side]; };
  int getLayerID(SimTrackerHit* hit)    { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::layer]; };
  int getModuleID(SimTrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::module]; };
  int getSensorID(SimTrackerHit* hit)   { _encoder->setValue(hit->getCellID0()); return (*_encoder)[lcio::ILDCellID0::sensor]; };
  
  void setupGearGeom( const gear::GearMgr* gearMgr ) ;
  
  
  unsigned int _nLayersVTX;
  
  unsigned int _nLayersSIT;
  
  
  std::vector<float> _zLayerFTD;
  
  unsigned int _nlayersFTD;
  bool _petalBasedFTDWithOverlaps;
  int _nPhiFTD; 

  int _output_track_col_quality;
  static const int _output_track_col_quality_GOOD;
  static const int _output_track_col_quality_FAIR;
  static const int _output_track_col_quality_POOR;

  int _sw_theta;//search window theta
  float _chi2FitCut_kalman;
  bool _useClusterRejection;
  float _minDotOf2Clusters;

  int getIntersectionEasy(HelixClass_double& helix, TrackerHit* curInmos , int layer, double* isec, double* ref); 
  int getIntersectionEasyTest(HelixClass_double& helix, TrackerHit* basisTrkhit , int layer, std::vector<double> &vec ); 
  int CheckTiltOf2Clusters(TrackerHit* A, TrackerHit* B, int level);
  float DotOf2Clusters(TrackerHit* A, TrackerHit* B);
  int KalFit(int& ndf, float& Chi2, TrackerHitVec trkHits,TrackerHitVec& hits_in_fit, TrackerHitVec& outliers, float* par , float* epar, HelixClass_double& helix);
  int getPhiThetaRegion(TrackExtended* trackAR, int layer, int* Boundaries);

  struct GeoData_t {
    int nladder;
    double rmin;  // distance of inner surface of sensitive region from IP
    double dphi;  // azimuthal angle step of each ladder
    double phi0;  // aximuthal angle offset
    std::vector<double> cosphi;  // cos[phi_ladder], cos_phi of each ladder
    std::vector<double> sinphi;  // sin[phi_ladder], sin_phi of each ladder
    std::vector<double> phi;  // phi of each ladder
    std::vector<double> phiAtXiMin;  // phiAtXiMin of each ladder
    std::vector<double> phiAtXiMax;  // phiAtXiMax of each ladder
    std::vector<double> ladder_incline;//the tilt of the line of the ladder expressed by phi
    double sthick;  // sensitive region thickness
    double sximin;  // minimum xi of sensitive region.
    double sximax;  // maximum xi of sensitive region
    double hlength; // ladder's half length in z
    int num_xi_pixel;      // Number of xi pixel in this ladder
    int num_zeta_pixel;    // Number of zeta pixel in this ladder
    double rmes; //distance in R of measurement surface 
  };
  
  struct vxdGeoData_t {
     int nLayer;
     int maxLadder;
     std::vector<GeoData_t> geodata;
  }_vxd,_sit;

  void InitVXDGeometry();
  void InitSITGeometry();
  FloatVec _pixelSizeVec;
  float _pixelheight;
  TVector3 LocalToGlobal(TVector3 local,int layer,int ladder);
  //purityMCP CheckOriginOf2Clusters(TrackerHit* A, TrackerHit* B);
  void calcTrackParameterOfMCP(MCParticle* pmcp, double* par);
  
  class ClusterStatus{
    public :
     int cellid0 ;
     int cellid1;
     unsigned int layer ;
     unsigned int ladder;
     unsigned int   xiwidth ;
     unsigned int zetawidth ;
     unsigned int      nPix ;
     unsigned int      tilt ;
     unsigned int   quality ;

     ClusterStatus(TrackerHit* hit){
        cellid0 = hit->getCellID0(); 
        cellid1 = hit->getCellID1();
        //layer  = ( cellid0 >> 5 ) & 0x000001ff;//maybe bug.
        //ladder = (cellid0 >> 15 ) & 0x000000ff;
        xiwidth =  cellid1 & 0x000001ff ;
        zetawidth = ( cellid1 >> 9 ) & 0x000001ff;
        nPix = ( cellid1 >> 18) & 0x000003ff;
        tilt = ( cellid1 >> 28) & 0x00000003;
        quality = ( cellid1 >> 30) & 0x00000003;
        const int cellId_trk = hit->getCellID0();
        UTIL::BitField64 encoder_trk( lcio::ILDCellID0::encoder_string ) ; 
        encoder_trk.setValue(cellId_trk) ;
        layer    = encoder_trk[lcio::ILDCellID0::layer] ;
        ladder    = encoder_trk[lcio::ILDCellID0::module] ;
     }
     ClusterStatus(){;}
  };

  //Old Ver//////////////////
  typedef std::map< std::pair< int, int >, int > RangeMap;
  RangeMap _phiRangeForTriplet;

  double getNeededPhiSectors(double Pt, int outly , int inly); //Difference of radian is returned

  //New Ver////////////////////==under construction====
  typedef std::map< std::vector<int> , std::vector<int> > RangeMapVer2;
  RangeMapVer2 _phiRangeForTripletVer2;
  void getNeededPhiSectorsVer2(double Pt, std::vector<int> layers, std::vector<double>& phiDiff);

  float _safetyPhiRange_ratio;
  int _safetyPhiRange_fix;
  float _fudgeFactorForSITsr_rphi;
  float _fudgeFactorForSITsr_z;
  int _fudgePhiRange;
  int _fudgeThetaRange;

  std::string _colNameVXDTrackerHitRelations;
  std::string _colNameSITSpacePointRelations;
  std::string _colNameFTDSpacePointRelations;
  std::string _colNameFTDPixelTrackerHitRelations;
  LCRelationNavigator* _navVXD;
  LCRelationNavigator* _navSIT;
  LCRelationNavigator* _navFTDsp;
  LCRelationNavigator* _navFTDpix;
  std::string _colNameVXDSimHit;
  std::string _colNameSITSimHit;
  std::string _colNameFTDspSimHit;
  std::string _colNameFTDpixSimHit;
  LCCollection* _simVXD;
  LCCollection* _simSIT;
  LCCollection* _simFTDsp;
  LCCollection* _simFTDpix;



  double _nSigmaBuild_phi;
  double _nSigmaBuild_theta;

  bool _keepCandidate;//used in AttachRemainingVTXHitsVeryFast <-- under construction for now

  moriUTIL* _moriUtil;
  GetPurityUtil* _purityUtil;

  bool _mydebug;
  bool _mydebugKalFit;
  bool _mydebugIntersection;
  bool _mydebugRoot;
  bool _mydebugRootVXD;
  bool _mydebugRootFTD;
  bool _mydebugRootVXDFTD;
  bool _mydebugVXDHits;
  bool _mydebugSITHits;
  bool _mydebugTriplet;
  int  _mydebugTripletMode;
  bool _mydebugTripletVXD;
  bool _mydebugTripletFTD;
  bool _mydebugTripletVXDFTD;
  bool _mydebugBuildTrack;
  int  _mydebugBuildTrackMode;
  bool _mydebugAttachVXD;
  bool _mydebugPrintMCP;
  bool _stopwatch;
  class Timer{
   public:
    TStopwatch inAnEvt;
    TStopwatch InitialiseVTX;
    TStopwatch InitialiseFTD;
    TStopwatch ProcessOneSector;
    TStopwatch TrackingInFTD;
    TStopwatch Sorting;
    TStopwatch CreateTrack;
    TStopwatch AttachRemainingVXD;
    TStopwatch AttachRemainingFTD;
    TStopwatch FinalRefit;
    void reset(){
      inAnEvt.Reset();
      InitialiseVTX.Reset();
      InitialiseFTD.Reset();
      ProcessOneSector.Reset();
      TrackingInFTD.Reset();
      Sorting.Reset();
      CreateTrack.Reset();
      AttachRemainingVXD.Reset();
      AttachRemainingFTD.Reset();
      FinalRefit.Reset();
    }
    void cout(){
      printf("inAnEvt           : RT=%.3f s, CPU=%.3f s \n",inAnEvt.RealTime(),inAnEvt.CpuTime());
      printf("InitialiseVTX     : RT=%.3f s, CPU=%.3f s \n",InitialiseVTX.RealTime(),InitialiseVTX.CpuTime());
      printf("InitialiseFTD     : RT=%.3f s, CPU=%.3f s \n",InitialiseFTD.RealTime(),InitialiseFTD.CpuTime());
      printf("ProcessOneSector  : RT=%.3f s, CPU=%.3f s \n",ProcessOneSector.RealTime(),ProcessOneSector.CpuTime());
      printf("TrackingInFTD     : RT=%.3f s, CPU=%.3f s \n",TrackingInFTD.RealTime(),TrackingInFTD.CpuTime());
      printf("Sorting           : RT=%.3f s, CPU=%.3f s \n",Sorting.RealTime(),Sorting.CpuTime());
      printf("CreateTrack       : RT=%.3f s, CPU=%.3f s \n",CreateTrack.RealTime(),CreateTrack.CpuTime());
      printf("AttachRemainingVXD: RT=%.3f s, CPU=%.3f s \n",AttachRemainingVXD.RealTime(),AttachRemainingVXD.CpuTime());
      printf("AttachRemainingFTD: RT=%.3f s, CPU=%.3f s \n",AttachRemainingFTD.RealTime(),AttachRemainingFTD.CpuTime());
      printf("FinalRefit        : RT=%.3f s, CPU=%.3f s \n",FinalRefit.RealTime(),FinalRefit.CpuTime());
    }
    Timer(){reset();}
  }_timer;
  TStopwatch _timer2Triplet;
  TStopwatch _timer2Build;
  bool  _mydebugstopwatch2;

  float _currentPurity;
  MCPMap LoadMCPMap();
  std::map< MCParticle*, SimTrackerHitVec > _mcpVXD;
  std::map< MCParticle*, SimTrackerHitVec > _mcpVXDFTD;
  std::map< MCParticle*, SimTrackerHitVec > _mcpVXDSIT;
  std::map< MCParticle*, SimTrackerHitVec > _mcpVXDFTDSIT;
  std::map< MCParticle*, SimTrackerHitVec > _mcpFTD;
  std::map< MCParticle*, SimTrackerHitVec > _mcpFTDSIT;
  std::map< MCParticle*, SimTrackerHitVec > _mcpSIT;


  enum MCPContributions{
     contVXD,
     contFTD,
     contSIT,
     contVXDFTD,
     contVXDSIT,
     contFTDSIT,
     contVXDFTDSIT,
     contBG,
     contSize
  };

  enum MCPContributions getMCPContribution(IntVec nsub);

  class Triplet : public IMPL::TrackImpl{
   public :
     static int numOfProcesses;
     int id;
     int quality_code;
     int detID[3];//内側に近いところから詰める。
     int layer[3];
     float probability;
     float purity;
     SimTrackerHitVec simVecFromTrk;//TrackerHitVecに対応するもの
     std::pair< MCParticle* , SimTrackerHitVec > mcTruth ;//Dominant MCPについての情報

     //new (2013_08_10)
     TrackerHitVec truthTrkHits ;//Dominant MCPのsimthitsに対応するTrackerHitVec.
     enum MCPContributions mcpcont;
     IntVec nsub;//vxd:0,sit:1,ftd:2 <--for Dominant MCP
     

     float mcp_d0;
     float mcp_z0;
     float mcp_omega;
     float mcp_phi0;
     float mcp_tanL;

     //bool availableInBuildedTrack;

     Triplet(bool incrementID = true){
       id = numOfProcesses; 
       if(incrementID){
         numOfProcesses++;
       }
       quality_code = -1;
       detID[0] = detID[1] = detID[2] = -1;
       layer[0] = layer[1] = layer[2] = -1;
       probability = 1e20;
       purity = 1e20;
       simVecFromTrk.clear();
       mcTruth.first = 0; mcTruth.second.clear();
       truthTrkHits.clear();
       mcpcont = contSize;
       nsub.clear();
       mcp_d0 = mcp_z0 = mcp_omega = mcp_phi0 = mcp_tanL = 1e20;
     }
     Triplet(const Triplet& obj){
       *this = obj;
     }
  };
  Triplet* _curtriplet;
  

  enum BuildTrackResult{
    NormalEnd,
    ManyMisAssignments,
    ManyOutliers,
    PhiThetaError,
    BuildTrackResultSize
  };
  class BuildedTrack : public IMPL::TrackImpl{
   public :
     static int numOfProcesses;
     int id;

     BuildTrackResult result;
     int nMisAssign;
     Triplet triplet;
     const MCParticle* truthmcp;
     SimTrackerHitVec simvec;
     float probability;
     float purity;
     SimTrackerHitVec simVecFromTrk;//TrackerHitVecに対応するもの
     BuildedTrack():triplet(false){
        id = numOfProcesses; 
        numOfProcesses++;
        result = BuildTrackResultSize;
        nMisAssign = -1;
        truthmcp = 0;
        simvec.clear();
        probability = 1e20;
        purity = 1e20;
        simVecFromTrk.clear();
     }
     BuildedTrack(const BuildedTrack& obj){
       *this = obj;
     }
  };
  bool _availableInBuildedTrack;

  typedef std::vector< BuildedTrack > BuildedTrackVec;
  BuildedTrackVec _buildedTrackContainer;
  void SetBuildedTrack(TrackExtended* trackAR, BuildedTrackVec& btrackvec);
  void BuildedTrackDebuger1(BuildedTrackVec::iterator begin, BuildedTrackVec::iterator end);
  void BuildedTrackDebuger2(std::vector<const BuildedTrack*>::iterator begin, std::vector<const BuildedTrack*>::iterator end);

  //時間がアレばclass templateにしてみたい。
  class MCP_BuildedTrack{
   public :
     std::pair<MCParticle*,SimTrackerHitVec> pair;
     IntVec nsub;//nhits in sub detectors. 0 VXD, 1 FTD, 2 SIT
     std::vector<const BuildedTrack*> BuildedTrackVec;
     std::vector<const BuildedTrack*> BuildedTrack2ndVec;
     std::vector<const BuildedTrack*> BuildedTrack3rdVec;
  };
  typedef std::map<MCParticle*,MCP_BuildedTrack> BTMap; 
  BTMap _mcpBuildedTrackContainer;

  std::vector< std::map<MCParticle*, MCP_BuildedTrack> > _mcpRemainingBTContainerA;
  //purity 100%のトリプレットが生成されなかった生成されなかったmcp
  std::vector< std::map<MCParticle*, MCP_BuildedTrack> > _mcpRemainingBTContainerB;
  //purity 100%のトリプレットが生成されたものの、quality_code != 0により結局再構成されなかったmcp
  



  std::vector< Triplet > _tripletContainer;

  class MCP_Triplet{
   public :
     std::pair<MCParticle*,SimTrackerHitVec> pair;
     IntVec nsub;//nhits in sub detectors. 0 VXD, 1 FTD, 2 SIT
     std::vector<const Triplet*> tripvec;
     std::vector<const Triplet*> trip2ndvec;
     std::vector<const Triplet*> trip3rdvec;
  };
  typedef std::map<MCParticle*,MCP_Triplet> TripMap; 
  TripMap _mcpTripletContainer;


  std::vector< std::map<MCParticle*, MCP_Triplet> > _mcpRemainingTRContainerA;
  //purity 100%のトリプレットが生成されなかった生成されなかったmcp
  std::vector< std::map<MCParticle*, MCP_Triplet> > _mcpRemainingTRContainerB;
  //purity 100%のトリプレットが生成されたものの、quality_code != 0により結局再構成されなかったmcp
  
  void TripletDebugerWithMCPRemain(int nth, std::vector< std::map<MCParticle*, MCP_Triplet> > A);

  MCPMap _mcpMap;
  std::vector<LCRelationNavigator*> _naviVec;
  IntVec getNHitsInSubDet(SimTrackerHitVec simvec);
  void TripletDebuger1(std::vector<Triplet>::iterator begin, std::vector<Triplet>::iterator end);
  void TripletDebuger2(std::vector<const Triplet*>::iterator begin, std::vector<const Triplet*>::iterator end);

  
};

int FPCCDSiliconTracking_MarlinTrk::Triplet::numOfProcesses = 0;
int FPCCDSiliconTracking_MarlinTrk::BuildedTrack::numOfProcesses = 0;

#endif