CCDDigitizer.h

/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */

#ifndef CCDDigitizer_h
#define CCDDigitizer_h 1

#include "marlin/Processor.h"
#include "lcio.h"
#include "EVENT/SimTrackerHit.h"
#include "IMPL/TrackerHitImpl.h"
#include "IMPL/SimTrackerHitImpl.h"
#include <string>
#include <vector>
#include "MyG4UniversalFluctuationForSi.h"
#include "EVENT/LCIO.h"
#include <IMPL/LCCollectionVec.h>



#ifdef CCD_diagnostics
#include <AIDA/IHistogramFactory.h>
#include <AIDA/IHistogram1D.h>
#include <AIDA/IHistogram2D.h>
#include <AIDA/IAxis.h>
#include <AIDA/ITree.h>

#include <marlin/Global.h>
#include <marlin/VerbosityLevels.h>
#include <marlin/AIDAProcessor.h>
#endif

using namespace lcio ;
using namespace marlin ;

struct IonisationPoint {
  double x;
  double y;
  double z;
  double eloss;
};


typedef std::vector<IonisationPoint> IonisationPointVec;
typedef std::vector<TrackerHitImpl*> TrackerHitImplVec;
typedef std::vector<SimTrackerHitImpl*> SimTrackerHitImplVec;



#define maxpixx 9// grid size in local(ladder) x axis
#define maxpixy 9// grid size in local(ladder) y axis
//grid size in which diffusion is computed, the charge, which diffuses outside this grid, is lost;value may be redefined again after sigmacoefficient or thickness of active layer is changed

//we approach the continual distribution of the diffusion with discrete points:
#define Numstepx 10// Number of points at which amplitude of diffusion is calculated within one pixel in x direction
#define Numstepy 10// Number of points at which amplitude of diffusion is calculated within one pixel in y direction
//according to first tests increasing the number of steps effects the performance of the processor only slightly and non- systematically

//if using a table, numhitstep must be adjusted when changing Numstep
//table
//#define numhitstepx 6//must be =(Numstepx/2)+1
//#define numhitstepy 6//must be =(Numstepy/2)+1
//#define numsigstep 20
//table//


class CCDDigitizer : public Processor {
  
 public:

  CCDDigitizer(const CCDDigitizer&) = delete;
  CCDDigitizer& operator=(const CCDDigitizer&) = delete;
  
  virtual Processor*  newProcessor() { return new CCDDigitizer ; }
  
  
  CCDDigitizer() ;
  
  virtual void init() ;
  
  virtual void processRunHeader( LCRunHeader* run ) ;
  
  virtual void processEvent( LCEvent * evt ) ; 
  
  virtual void check( LCEvent * evt ) ; 
  
  virtual void end() ;
  
  
 protected:

  std::string _colName{};
  std::string _outputCollectionName{};
  std::string _colVTXRelation{};

  int _nRun{};
  
  int _nEvt{};
  double _cutOnDeltaRays{};
  double _diffusionCoefficient{}; 



  int _numberOfLayers{};
  double _pixelSizeX{};
  double _pixelSizeY{};
  double _electronsPerKeV{};
  double _segmentDepth{};
  double _currentTotalCharge{};

  std::vector<int> _laddersInLayer{};
  std::vector<float> _layerRadius{};
  std::vector<float> _layerThickness{};
  std::vector<float> _layerHalfThickness{};
  std::vector<float> _layerLadderLength{};
  std::vector<float> _layerLadderHalfWidth{};
  std::vector<float> _layerPhiOffset{};
  std::vector<float> _layerActiveSiOffset{};
  std::vector<float> _layerHalfPhi{};
  std::vector<float> _layerLadderGap{};
  std::vector<float> _bkgdHitsInLayer{};
  std::vector<float> _layerLadderWidth{};

  int _currentLayer{};
  int _currentModule{};
  int _generateBackground{};
  double _currentParticleMomentum{};
  double _currentParticleEnergy{};
  double _currentParticleMass{};
  double _currentPhi{};

  double PI{},TWOPI{},PI2{};

  int _produceFullPattern{};
  int _numberOfSegments{};
  int _debug{};
  int _PoissonSmearing{};
  int _electronicEffects{};
  int _useMCPMomentum{};
 
  int _recmethod{};
  double _threshold{};
  double _saturation{};
  int _framesize{};
  int maxnionpoint{};
  
  double _currentLocalPosition[3]{};
  double _electronicNoise{};
  double _segmentLength{};

  IonisationPointVec _ionisationPoints{};
 

  MyG4UniversalFluctuationForSi * _fluctuate{};

   
  //Finds coordinates    
  void FindLocalPosition(SimTrackerHit * hit,
                         double * localPosition,
                         double * localDirection);

  void TransformToLab(double * xLoc, double * xLab);
  void ProduceIonisationPoints( SimTrackerHit * hit);
  void diffusion(double xdif,double ydif, double sigma);

  void ProduceHits(SimTrackerHitImplVec & simTrkVec);
  void TransformXYToCellID(double x, double y, 
                           int & ix, 
                           int & iy,double & xdif,double & ydif);  
  void TransformCellIDToXY(int ix, int iy,
                           double & x, double & y);
  void PoissonSmearer( SimTrackerHitImplVec & simTrkVec );
  void GainSmearer( SimTrackerHitImplVec & simTrkVec );
  void PrintInfo( SimTrackerHit * simTrkHit, TrackerHitImpl * recoHit);
  TrackerHitImpl * ReconstructTrackerHit(SimTrackerHitImplVec & simTrkVec );
  void TrackerHitToLab( TrackerHitImpl * recoHit );

  void generateBackground(LCCollectionVec * col);
  void settanlorentzangle(double B, double E, double mu, double T);
  void settanlorentzangleb(double B, double E, double mu, double T);

//  double _currentLocalPosition[3];
//   double _currentEntryPoint[3];
//   double _currentExitPoint[3];
 
  double _energyLoss{};
  std::vector <SimTrackerHitImplVec> _hitsInLayer{};

  double depdep{};
  double undep{};
  double epitaxdep{};

  double pxl[maxpixx][maxpixy]{};
  int midpixx{};
  int midpixy{};
  double stepx{};
  double stepy{};
  double xobsoffset{};
  double yobsoffset{};  
 
  double sigmacoefficient{};
  double sigmin{};
  double _difcoef{};
  double _efield{};
  double _biasvolt{};
  double _bfield{};
  double _T{};
  double _mu{};
  double TanLorentzAngle{};


#ifdef CCD_diagnostics
  double dirraw[3]{};
  double posraw[3]{};
  AIDA::IHistogram1D* histdist{};
  AIDA::IHistogram1D* histcluster{};
  AIDA::IHistogram2D* histclustxy{};
  AIDA::IHistogram1D* histcharge{};
  AIDA::IHistogram2D* histdistxy{};
  AIDA::IHistogram1D* histNionpoint{};
  AIDA::IHistogram1D* histzcoord{};
  AIDA::IHistogram1D* histenergy{};
  AIDA::IHistogram1D* histsignal{};
  AIDA::IHistogram1D* histsignalframe{};
  AIDA::IHistogram1D* histenergycentre{};
  int Nionpoint{};
#endif

 // table
 //  double sigstep;
//   double table [numsigstep][numhitstepx][numhitstepy][maxpixx][maxpixy];
//   void diffusiontable(double xdif,double ydif, double sigma);
//   void settable();
// table//

};



#endif