CbmMvdSensorDigitizerTask.cxx

Go to the documentation of this file.

// -------------------------------------------------------------------------
// -----                  CbmMvdSensorDigitizerTask source file              -----
// -----                  Created 02.02.2012 by M. Deveaux            -----
// -------------------------------------------------------------------------
#include "CbmMvdSensorDigitizerTask.h"
 
#include "TClonesArray.h"
#include "TObjArray.h"
 
#include "CbmMvdPileupManager.h"
#include "CbmMvdPoint.h"
 
#include "FairRuntimeDb.h"
 
// Includes from FairRoot
#include "FairEventHeader.h"
#include "FairLogger.h"
#include "FairMCEventHeader.h"
#include "FairRunAna.h"
#include "FairRunSim.h"
 
 
// Includes from C++
#include <iomanip>
#include <iostream>
#include <map>
#include <vector>
 
using std::cout;
using std::endl;
using std::fixed;
using std::ios_base;
using std::left;
using std::map;
using std::pair;
using std::right;
using std::setprecision;
using std::setw;
using std::vector;
 
 
// -----   Default constructor   ------------------------------------------
CbmMvdSensorDigitizerTask::CbmMvdSensorDigitizerTask()
  : CbmMvdSensorTask()
  , fcurrentFrameNumber(0)
  , fEpiTh(0.0014)
  , fSegmentLength(0.0001)
  , fDiffusionCoefficient(0.0055)
  , fElectronsPerKeV(276.)
  , fWidthOfCluster(3.5)
  , fPixelSizeX(0.0030)
  , fPixelSizeY(0.0030)
  , fCutOnDeltaRays(0.00169720)
  , fChargeThreshold(100.)
  , fFanoSilicium(0.115)
  , fEsum(0.)
  , fSegmentDepth(0.)
  , fCurrentTotalCharge(0.)
  , fCurrentParticleMass(0.)
  , fCurrentParticleMomentum(0.)
  , fCurrentParticlePdg(0)
  , fRandomGeneratorTestHisto(NULL)
  , fPosXY(NULL)
  , fpZ(NULL)
  , fPosXinIOut(NULL)
  , fAngle(NULL)
  , fSegResolutionHistoX(NULL)
  , fSegResolutionHistoY(NULL)
  , fSegResolutionHistoZ(NULL)
  , fTotalChargeHisto(NULL)
  , fTotalSegmentChargeHisto(NULL)
  , fLorentzY0(-6.1)
  , fLorentzXc(0.)
  , fLorentzW(1.03)
  , fLorentzA(477.2)
  , fLorentzNorm(1.)
  , fLandauMPV(877.4)
  , fLandauSigma(204.93)
  , fLandauGain(3.3)
  , fLandauRandom(new TRandom3())
  , fPixelSize(0.0025)
  , fPar0(520.)
  , fPar1(0.34)
  , fPar2(-1.2)
  , fCompression(1.)
  , fResolutionHistoX(NULL)
  , fResolutionHistoY(NULL)
  , fNumberOfSegments(0)
  , fCurrentLayer(0)
  , fEvent(0)
  , fVolumeId(0)
  , fNPixelsX(0)
  , fNPixelsY(0)
  , fPixelCharge(new TClonesArray("CbmMvdPixelCharge"))
  , fDigis(NULL)
  , fDigiMatch(NULL)
  , frand(nullptr)
  , fproduceNoise(kFALSE)
  , fPixelChargeShort()
  , fPixelScanAccelerator(NULL)
  , fChargeMap()
  , fChargeMapIt()
  , fsensorDataSheet(NULL)
  , fMode(0)
  , fSigmaX(0.0005)
  , fSigmaY(0.0005)
  , fReadoutTime(0.00005)
  , fEfficiency(-1.)
  , fMergeDist(-1.)
  , fFakeRate(-1.)
  , fNPileup(0)
  , fNDeltaElect(0)
  , fDeltaBufferSize(0)
  , fBgBufferSize(0)
  , fBranchName("")
  , fBgFileName("")
  , fDeltaFileName("")
  , fInputPoints(NULL)
  , fPoints(new TRefArray())
  , fRandGen()
  , fTimer()
  , fPileupManager(NULL)
  , fDeltaManager(NULL)
  , fNEvents(0)
  , fNPoints(0.)
  , fNReal(0.)
  , fNBg(0.)
  , fNFake(0.)
  , fNLost(0.)
  , fNMerged(0.)
  , fTime(0.)
  , fSignalPoints()
  , h_trackLength(NULL)
  , h_numSegments(NULL)
  , h_LengthVsAngle(NULL)
  , h_LengthVsEloss(NULL)
  , h_ElossVsMomIn(NULL) {
  fRandGen.SetSeed(2736);
  frand         = new TRandom3(0);
  fproduceNoise = kFALSE;
}
// -------------------------------------------------------------------------
 
 
// -----   Standard constructor   ------------------------------------------
CbmMvdSensorDigitizerTask::CbmMvdSensorDigitizerTask(Int_t iMode)
  : CbmMvdSensorTask()
  , fcurrentFrameNumber(0)
  , fEpiTh(0.0014)
  , fSegmentLength(0.0001)
  , fDiffusionCoefficient(0.0055)
  , fElectronsPerKeV(276.)
  , fWidthOfCluster(3.5)
  , fPixelSizeX(0.0030)
  , fPixelSizeY(0.0030)
  , fCutOnDeltaRays(0.00169720)
  , fChargeThreshold(100.)
  , fFanoSilicium(0.115)
  , fEsum(0.)
  , fSegmentDepth(0.)
  , fCurrentTotalCharge(0.)
  , fCurrentParticleMass(0.)
  , fCurrentParticleMomentum(0.)
  , fCurrentParticlePdg(0)
  , fRandomGeneratorTestHisto(NULL)
  , fPosXY(NULL)
  , fpZ(NULL)
  , fPosXinIOut(NULL)
  , fAngle(NULL)
  , fSegResolutionHistoX(NULL)
  , fSegResolutionHistoY(NULL)
  , fSegResolutionHistoZ(NULL)
  , fTotalChargeHisto(NULL)
  , fTotalSegmentChargeHisto(NULL)
  , fLorentzY0(-6.1)
  , fLorentzXc(0.)
  , fLorentzW(1.03)
  , fLorentzA(477.2)
  , fLorentzNorm(1.)
  , fLandauMPV(877.4)
  , fLandauSigma(204.93)
  , fLandauGain(3.3)
  , fLandauRandom(new TRandom3())
  , fPixelSize(0.0025)
  , fPar0(520.)
  , fPar1(0.34)
  , fPar2(-1.2)
  , fCompression(1.)
  , fResolutionHistoX(NULL)
  , fResolutionHistoY(NULL)
  , fNumberOfSegments(0)
  , fCurrentLayer(0)
  , fEvent(0)
  , fVolumeId(0)
  , fNPixelsX(0)
  , fNPixelsY(0)
  , fPixelCharge(new TClonesArray("CbmMvdPixelCharge", 100000))
  , fDigis(NULL)
  , fDigiMatch(NULL)
  , frand(nullptr)
  , fproduceNoise(kFALSE)
  , fPixelChargeShort()
  , fPixelScanAccelerator(NULL)
  , fChargeMap()
  , fChargeMapIt()
  , fsensorDataSheet(NULL)
  , fMode(iMode)
  , fSigmaX(0.0005)
  , fSigmaY(0.0005)
  , fReadoutTime(0.00005)
  , fEfficiency(-1.)
  , fMergeDist(-1.)
  , fFakeRate(-1.)
  , fNPileup(0)
  , fNDeltaElect(0)
  , fDeltaBufferSize(0)
  , fBgBufferSize(0)
  , fBranchName("MvdPoint")
  , fBgFileName("")
  , fDeltaFileName("")
  , fInputPoints(NULL)
  , fPoints(new TRefArray())
  , fRandGen()
  , fTimer()
  , fPileupManager(NULL)
  , fDeltaManager(NULL)
  , fNEvents(0)
  , fNPoints(0.)
  , fNReal(0.)
  , fNBg(0.)
  , fNFake(0.)
  , fNLost(0.)
  , fNMerged(0.)
  , fTime(0.)
  , fSignalPoints()
  , h_trackLength(NULL)
  , h_numSegments(NULL)
  , h_LengthVsAngle(NULL)
  , h_LengthVsEloss(NULL)
  , h_ElossVsMomIn(NULL) {
  cout << "Starting CbmMvdSensorDigitizerTask::CbmMvdSensorDigitizerTask() "
       << endl;
 
  fRandGen.SetSeed(2736);
  fEvent       = 0;
  fTime        = 0.;
  fSigmaX      = 0.0005;
  fSigmaY      = 0.0005;
  fReadoutTime = 0.00005;
 
 
  fEpiTh         = 0.0014;
  fSegmentLength = 0.0001;
  fDiffusionCoefficient =
    0.0055;  // correspondes to the sigma of the gauss with the max drift length
  fElectronsPerKeV         = 276;     //3.62 eV for e-h creation
  fWidthOfCluster          = 3.5;     // in sigmas
  fPixelSizeX              = 0.0025;  // in cm
  fPixelSizeY              = 0.0025;
  fCutOnDeltaRays          = 0.00169720;  //MeV
  fChargeThreshold         = 100.;        //electrons change 1 to 10
  fFanoSilicium            = 0.115;
  fEsum                    = 0;
  fSegmentDepth            = 0;
  fCurrentTotalCharge      = 0;
  fCurrentParticleMass     = 0;
  fCurrentParticleMomentum = 0;
  fPixelScanAccelerator    = 0;
 
 
  fPixelSize   = 0.0025;
  fPar0        = 520.;
  fPar1        = 0.34;
  fPar2        = -1.2;
  fLandauMPV   = 877.4;
  fLandauSigma = 204.93;
  fLandauGain  = 3.3;
 
  fLorentzY0 = -6.1;
  fLorentzXc = 0.;
  fLorentzW  = 1.03;
  fLorentzA  = 477.2;
 
 
  fCompression = 1.;
 
 
  //fLorentzNorm=0.00013010281679422413;
  fLorentzNorm = 1;
 
 
  fcurrentFrameNumber = 0;
 
  frand         = new TRandom3(0);
  fproduceNoise = kFALSE;
}
 
 
// -------------------------------------------------------------------------
 
 
// -----   Destructor   ----------------------------------------------------
CbmMvdSensorDigitizerTask::~CbmMvdSensorDigitizerTask() {
 
  if (fPileupManager) delete fPileupManager;
  if (fDeltaManager) delete fDeltaManager;
  if (fInputPoints) {
    fInputPoints->Delete();
    delete fInputPoints;
  }
  if (fOutputBuffer) {
    fOutputBuffer->Delete();
    delete fOutputBuffer;
  }
}
// ------------------------------------------------------------------------
 
 
// -----    Virtual private method ReadSensorInformation   ----------------
InitStatus CbmMvdSensorDigitizerTask::ReadSensorInformation() {
 
  CbmMvdSensorDataSheet* sensorData;
  sensorData = fSensor->GetDataSheet();
  if (!sensorData) { return kERROR; }
 
  fPixelSizeX = sensorData->GetPixelPitchX();
  fPixelSizeY = sensorData->GetPixelPitchY();
  fNPixelsX   = sensorData->GetNPixelsX();
  fNPixelsY   = sensorData->GetNPixelsY();
 
  fChargeThreshold = sensorData->GetChargeThreshold();
 
  fPar0 = sensorData->GetLorentzPar0();
  fPar1 = sensorData->GetLorentzPar1();
  fPar2 =
    sensorData
      ->GetLorentzPar2();  //cout<< endl << " LorentzPar2 is now set to " << fPar2 << endl;
 
  fLandauMPV =
    sensorData
      ->GetLandauMPV();  //cout << endl << " Landau MPV is now set to " << fLandauMPV << endl;
  fLandauSigma =
    sensorData
      ->GetLandauSigma();  //cout << endl << " Landau Sigma is now set to " << fLandauSigma << endl;
  fLandauGain =
    sensorData
      ->GetLandauGain();  //cout << endl << " Landau Gain is now set to " << fLandauGain << endl;
  fEpiTh =
    sensorData
      ->GetEpiThickness();  //cout << endl << " Epitaxial thickness is now set to " << fEpiTh << endl;
  fCompression = fPixelSizeX / fPixelSizeY;
 
  if (fCompression != 1)
    //cout << endl << "working with non uniform pixels" << endl;
    if (fCompression <= 0) fCompression = 1;
  return kSUCCESS;
}
 
 
// -----------------------------------------------------------------------------
 
void CbmMvdSensorDigitizerTask::SetInputArray(TClonesArray* inputStream) {
 
  Int_t i       = 0;
  Int_t nInputs = inputStream->GetEntriesFast();
  while (nInputs > i) {
    new ((*fInputPoints)[fInputPoints->GetEntriesFast()])
      CbmMvdPoint(*((CbmMvdPoint*) inputStream->At(i)));
 
    //cout << endl << "New Input registered" << endl;
    i++;
  }
}
// -----------------------------------------------------------------------------
void CbmMvdSensorDigitizerTask::SetInput(CbmMvdPoint* point) {
 
  new ((*fInputPoints)[fInputPoints->GetEntriesFast()])
    CbmMvdPoint(*((CbmMvdPoint*) point));
}
 
 
// -------------- public method ExecChain   ------------------------------------
void CbmMvdSensorDigitizerTask::ExecChain() {
 
  //   cout << endl << "start Digitizer on sensor " << fSensor->GetName() << endl;
 
 
  Exec();
}
 
// -----   Virtual public method Exec   ------------------------------------
void CbmMvdSensorDigitizerTask::Exec() {
 
  if (fPreviousPlugin) {
    fInputPoints->Delete();
    fInputPoints->AbsorbObjects(fPreviousPlugin->GetOutputArray());
  }
  fOutputBuffer->Clear("C");
  fDigis->Clear("C");
  fDigiMatch->Clear("C");
  if (fInputPoints->GetEntriesFast() > 0) {
 
    for (Int_t iPoint = 0; iPoint < fInputPoints->GetEntriesFast(); iPoint++) {
 
      CbmMvdPoint* point = (CbmMvdPoint*) fInputPoints->At(iPoint);
 
      if (!point) {
        cout << "-W-" << GetName() << ":: Exec:" << endl;
        cout << "    -received bad MC-Point. Ignored." << endl;
        continue;
      }
      if (point->GetStationNr() != fSensor->GetSensorNr()) {
        cout << "-W-" << GetName() << ":: Exec:" << endl;
        cout << "    -received bad MC-Point which doesn't belong here. Ignored."
             << endl;
        continue;
      }
      fcurrentFrameNumber = point->GetFrame();
      //The digitizer acts only on particles, which crossed the station.
      //Particles generated in the sensor or being absorbed in this sensor are ignored
      if (TMath::Abs(point->GetZOut() - point->GetZ()) < 0.9 * fEpiTh) {
        LOG(debug) << "hit not on chip with thickness "
                   << 0.9 * 2 * fSensor->GetDZ();
        LOG(debug) << "hit not on chip with length "
                   << TMath::Abs(point->GetZOut() - point->GetZ());
        continue;
      }
      // Reject for the time being light nuclei (no digitization modell yet)
      if (point->GetPdgCode() > 100000) { continue; }
      // Digitize the point
      //cout << endl << " found point make digi" << endl;
      ProduceIonisationPoints(point);
      ProducePixelCharge(point);
    }  //loop on MCpoints
 
 
    Int_t inputNr      = 0;
    Int_t eventNr      = 0;
    Double_t eventTime = .0;
    Int_t nDigis       = 0;
    GetEventInfo(inputNr, eventNr, eventTime);
 
    if (fproduceNoise) ProduceNoise();
 
    for (Int_t i = 0; i < fPixelCharge->GetEntriesFast(); i++) {
      CbmMvdPixelCharge* pixel = (CbmMvdPixelCharge*) fPixelCharge->At(i);
 
      if (pixel->GetCharge() > fChargeThreshold) {
        nDigis = fDigis->GetEntriesFast();
 
        new ((*fDigis)[nDigis]) CbmMvdDigi(fSensor->GetSensorNr(),
                                           pixel->GetX(),
                                           pixel->GetY(),
                                           pixel->GetCharge(),
                                           fPixelSizeX,
                                           fPixelSizeY,
                                           pixel->GetTime(),
                                           pixel->GetFrame());
 
 
        new ((*fOutputBuffer)[nDigis]) CbmMvdDigi(fSensor->GetSensorNr(),
                                                  pixel->GetX(),
                                                  pixel->GetY(),
                                                  pixel->GetCharge(),
                                                  fPixelSizeX,
                                                  fPixelSizeY,
                                                  pixel->GetTime(),
                                                  pixel->GetFrame());
 
        new ((*fDigiMatch)[nDigis]) CbmMatch();
        CbmMatch* match = (CbmMatch*) fDigiMatch->At(nDigis);
        for (Int_t iLink = 0; iLink < pixel->GetNContributors(); iLink++) {
          if (pixel->GetTrackID()[iLink] > -1)
            match->AddLink((Double_t) pixel->GetPointWeight()[iLink],
                           pixel->GetPointID()[iLink],
                           eventNr,
                           inputNr);
          else
            match->AddLink((Double_t) pixel->GetPointWeight()[iLink],
                           pixel->GetPointID()[iLink]);
        }
      } else {
        //cout << endl << "charge under threshold, digi rejected" << endl;
      }
    }
    //LOG(debug) << nDigis <<" new Digis at on Sensor " << fSensor->GetSensorNr() << " from " << fInputPoints->GetEntriesFast()<< " McPoints";
 
  }
 
  else {
    //LOG(debug)<< "No input found on Sensor " << fSensor->GetSensorNr() << " from " << fInputPoints->GetEntriesFast()<< " McPoints";
  }
 
  fPixelCharge->Delete();
  fChargeMap.clear();
  fInputPoints->Delete();
  fSignalPoints.clear();
 
}  // end of exec
 
// -------------------------------------------------------------------------
 
// -------------------------------------------------------------------------
void CbmMvdSensorDigitizerTask::GetEventInfo(Int_t& inputNr,
                                             Int_t& eventNr,
                                             Double_t& eventTime) {
 
  // --- The event number is taken from the FairRootManager
  eventNr = FairRootManager::Instance()->GetEntryNr();
 
  // --- In a FairRunAna, take the information from FairEventHeader
  if (FairRunAna::Instance()) {
    FairEventHeader* event = FairRunAna::Instance()->GetEventHeader();
    inputNr                = event->GetInputFileId();
    eventTime              = event->GetEventTime();
  }
 
  // --- In a FairRunSim, the input number and event time are always zero.
  else {
    if (!FairRunSim::Instance())
      LOG(fatal) << GetName() << ": neither SIM nor ANA run.";
    inputNr   = 0;
    eventTime = 0.;
  }
}
// -------------------------------------------------------------------------
 
// -------------------------------------------------------------------------
void CbmMvdSensorDigitizerTask::ProduceIonisationPoints(CbmMvdPoint* point) {
  //Option_t* opt1;
  //cout << endl << "Computing Point "   << endl;
  //point->Print(opt1);
 
  // Int_t pdgCode = point->GetPdgCode();
 
  //Transform coordinates of the point into sensor frame
 
  Double_t globalPositionIn[3]  = {point->GetX(), point->GetY(), point->GetZ()};
  Double_t globalPositionOut[3] = {
    point->GetXOut(), point->GetYOut(), point->GetZOut()};
 
  Double_t localPositionIn[3] = {0, 0, 0};
 
  Double_t localPositionOut[3] = {0, 0, 0};
 
  fSensor->TopToLocal(globalPositionIn, localPositionIn);
  fSensor->TopToLocal(globalPositionOut, localPositionOut);
 
  Int_t pixelX, pixelY;
  fSensor->LocalToPixel(&localPositionIn[0], pixelX, pixelY);
 
  // Copy results into variables used by earlier versions
 
  Double_t entryX = localPositionIn[0];
  Double_t exitX  = localPositionOut[0];
  Double_t entryY = localPositionIn[1];
  Double_t exitY  = localPositionOut[1];
  Double_t entryZ = localPositionIn[2];
  Double_t exitZ  = localPositionOut[2];
 
  // entry and exit from the epi layer ( det ref frame ) :
  Double_t entryZepi = -fEpiTh / 2;
  Double_t exitZepi  = fEpiTh / 2;
 
 
  TVector3 a(entryX, entryY, entryZ);  // entry in the detector
  TVector3 b(exitX, exitY, exitZ);     // exit from the detector
  TVector3 c;
 
  c = b - a;  // AB in schema
 
  TVector3 d;
  TVector3 e;
 
  Double_t scale1 = (entryZepi - entryZ) / (exitZ - entryZ);
  Double_t scale2 = (exitZepi - entryZ) / (exitZ - entryZ);
 
 
  d = c * scale1;
  e = c * scale2;
 
  TVector3 entryEpiCoord;
  TVector3 exitEpiCoord;
 
  entryEpiCoord = d + a;
  exitEpiCoord  = e + a;
 
 
  //Get x and y coordinates at the ENTRY of the epi layer
  Double_t entryXepi = entryEpiCoord.X();
  Double_t entryYepi = entryEpiCoord.Y();
  entryZepi          = entryEpiCoord.Z();
 
  //Get x and y coordinates at the EXIT of the epi layer
  Double_t exitXepi = exitEpiCoord.X();
  Double_t exitYepi = exitEpiCoord.Y();
  exitZepi          = exitEpiCoord.Z();
 
 
  Double_t lx = -(entryXepi - exitXepi);  //length of segment x-direction
  Double_t ly = -(entryYepi - exitYepi);
  Double_t lz = -(entryZepi - exitZepi);
 
 
  //-----------------------------------------------------------
 
 
  Double_t rawLength =
    sqrt(lx * lx + ly * ly
         + lz * lz);  //length of the track inside the epi-layer, in cm
  Double_t trackLength = 0;
 
  if (rawLength < 1.0e+3) {
    trackLength = rawLength;
  }
 
  else {
    cout << "-W- " << GetName() << " : rawlength > 1.0e+3 : " << rawLength
         << endl;
    trackLength = 1.0e+3;
  }
 
  //Smear the energy on each track segment
  Double_t charge =
    fLandauRandom->Landau(fLandauGain, fLandauSigma / fLandauMPV);
 
  //if (fShowDebugHistos )  cout << endl << "charge " << charge << endl;
 
  if (charge > (12000 / fLandauMPV)) {
    charge = 12000 / fLandauMPV;
  }  //limit Random generator behaviour
 
  if (fShowDebugHistos) {
    fRandomGeneratorTestHisto->Fill(charge * fLandauMPV);
  }
  //Translate the charge to normalized energy
 
  //     cout << endl << "charge after random generator " << charge << endl;
  Double_t dEmean = charge / (fElectronsPerKeV * 1e6);
  //   cout << endl << "dEmean " << dEmean << endl;
  fNumberOfSegments = int(trackLength / fSegmentLength) + 1;
 
  dEmean =
    dEmean
    * ((Double_t) trackLength / fEpiTh);  //scale the energy to the track length
 
  dEmean =
    dEmean
    / ((Double_t)
         fNumberOfSegments);  // From this point dEmean corresponds to the E lost per segment.
 
 
  fSignalPoints.resize(fNumberOfSegments);
 
  fEsum = 0.0;
 
  //Double_t segmentLength_update = trackLength/((Double_t)fNumberOfSegments);
 
  if (lz != 0) {
    fSegmentDepth = fEpiTh / ((Double_t) fNumberOfSegments);
  } else {  //condition added 05/08/08
    fSegmentDepth = 0;
    cout << "-W- " << GetName()
         << " Length of track in detector (z-direction) is 0!!!" << endl;
  }
 
 
  Double_t x = 0, y = 0, z = 0;
 
  Double_t xDebug = 0, yDebug = 0, zDebug = 0;
  Float_t totalSegmentCharge = 0;
 
  for (int i = 0; i < fNumberOfSegments; ++i) {
 
    z = -fEpiTh / 2
        + ((double) (i) + 0.5)
            * fSegmentDepth;  //middle position of the segment; zdirection
    x =
      entryXepi
      + ((double) (i) + 0.5)
          * (lx
             / ((Double_t)
                  fNumberOfSegments));  //middle position of the segment; xdirection
    y =
      entryYepi
      + ((double) (i) + 0.5)
          * (ly
             / ((Double_t)
                  fNumberOfSegments));  //middle position of the segment; ydirection
 
    if (fShowDebugHistos) {
      xDebug = xDebug + x;
      yDebug = yDebug + y;
      zDebug = zDebug + z;
    };
 
    SignalPoint* sPoint = &fSignalPoints[i];
 
    fEsum         = fEsum + dEmean;
    sPoint->eloss = dEmean;
    sPoint->x =
      x;  //here the coordinates x,y,z are given in the sensor reference frame.
    sPoint->y = y;
    sPoint->z = z;
    charge    = 1.0e+6 * dEmean * fElectronsPerKeV;
    //cout << endl << "charge " << charge << endl;
    sPoint->sigmaX     = fPixelSize;
    sPoint->sigmaY     = fPixelSize;
    sPoint->charge     = charge;
    totalSegmentCharge = totalSegmentCharge + charge;
  }
  //if (fShowDebugHistos  )cout << endl << "totalSegmentCharge " << totalSegmentCharge << endl;
  if (fShowDebugHistos) {
    fTotalSegmentChargeHisto->Fill(totalSegmentCharge * fLandauMPV);
    fSegResolutionHistoX->Fill(xDebug / fNumberOfSegments
                               - (point->GetX() + point->GetXOut()) / 2
                               - fSensor->GetX());
    fSegResolutionHistoY->Fill(yDebug / fNumberOfSegments
                               - (point->GetY() + point->GetYOut()) / 2
                               - fSensor->GetY());
    fSegResolutionHistoZ->Fill(zDebug / fNumberOfSegments
                               - (point->GetZ() + point->GetZOut()) / 2
                               - fSensor->GetZ());
  }
}
 
 
// -------------------------------------------------------------------------
 
void CbmMvdSensorDigitizerTask::ProducePixelCharge(CbmMvdPoint* point) {
  fCurrentTotalCharge = 0.0;
 
 
  // MDx - Variables needed in order to compute a "Monte Carlo Center of Gravity" of the cluster
 
  Float_t xCharge = 0., yCharge = 0., totClusterCharge = 0.;
  CbmMvdPixelCharge* pixel;
 
  pair<Int_t, Int_t> thispoint;
 
  Double_t xCentre, yCentre, sigmaX, sigmaY, xLo, xUp, yLo, yUp;
 
  SignalPoint* sPoint;
  sPoint = &fSignalPoints[0];
 
  xCentre = sPoint->x;  //of segment
  yCentre = sPoint->y;  
  sigmaX  = sPoint->sigmaX;
  sigmaY  = sPoint->sigmaY;
 
  xLo = sPoint->x - fWidthOfCluster * sigmaX;
  xUp = sPoint->x + fWidthOfCluster * sigmaX;
  yLo = sPoint->y - fWidthOfCluster * sigmaY;
  yUp = sPoint->y + fWidthOfCluster * sigmaY;
 
  if (fNumberOfSegments < 2) {
    Fatal("-E- CbmMvdDigitizer: ",
          "fNumberOfSegments < 2, this makes no sense, check parameters.");
  }
 
  Int_t* lowerXArray = new Int_t[fNumberOfSegments];
  Int_t* upperXArray = new Int_t[fNumberOfSegments];
  Int_t* lowerYArray = new Int_t[fNumberOfSegments];
  Int_t* upperYArray = new Int_t[fNumberOfSegments];
 
  Int_t ixLo, ixUp, iyLo, iyUp;
 
 
  Double_t minCoord[] = {xLo, yLo};
  Double_t maxCoord[] = {xUp, yUp};
 
 
  fSensor->LocalToPixel(minCoord, lowerXArray[0], lowerYArray[0]);
  fSensor->LocalToPixel(maxCoord, upperXArray[0], upperYArray[0]);
 
 
  if (lowerXArray[0] < 0) lowerXArray[0] = 0;
  if (lowerYArray[0] < 0) lowerYArray[0] = 0;
  if (upperXArray[0] > fNPixelsX) upperXArray[0] = fNPixelsX;
  if (upperYArray[0] > fNPixelsY) upperYArray[0] = fNPixelsY;
 
  ixLo = lowerXArray[0];
  iyLo = lowerYArray[0];
  ixUp = upperXArray[0];
  iyUp = upperYArray[0];
 
 
  for (Int_t i = 1; i < fNumberOfSegments; i++) {
 
    sPoint = &fSignalPoints[i];
    sigmaX = sPoint->sigmaX;
    sigmaY = sPoint->sigmaY;
 
    minCoord[0] = sPoint->x - fWidthOfCluster * sigmaX;
    minCoord[1] = sPoint->y - fWidthOfCluster * sigmaY;
    maxCoord[0] = sPoint->x + fWidthOfCluster * sigmaX;
    maxCoord[1] = sPoint->y + fWidthOfCluster * sigmaY;
 
 
    fSensor->LocalToPixel(minCoord, lowerXArray[i], lowerYArray[i]);
    fSensor->LocalToPixel(maxCoord, upperXArray[i], upperYArray[i]);
 
    if (lowerXArray[i] < 0) lowerXArray[i] = 0;
    if (lowerYArray[i] < 0) lowerYArray[i] = 0;
 
    if (upperXArray[i] > fNPixelsX) upperXArray[i] = fNPixelsX;
    if (upperYArray[i] > fNPixelsY) upperYArray[i] = fNPixelsY;
 
    if (ixLo > lowerXArray[i]) { ixLo = lowerXArray[i]; }
    if (ixUp < upperXArray[i]) { ixUp = upperXArray[i]; }
    if (iyLo > lowerYArray[i]) { iyLo = lowerYArray[i]; }
    if (iyUp < upperYArray[i]) { iyUp = upperYArray[i]; }
  }
 
 
  //cout << "Scanning from x= " << ixLo << " to " <<ixUp <<" and  y= "<<iyLo<< " to " << iyUp << endl;
 
  // loop over all pads of interest.
  fPixelChargeShort.clear();
  Int_t ix, iy;
 
 
  for (ix = ixLo; ix < ixUp + 1; ix++) {
 
    for (iy = iyLo; iy < iyUp + 1; iy++) {
 
      //calculate the position of the current pixel in the lab-system
 
      Double_t Current[3];
      fSensor->PixelToLocal(ix, iy, Current);
      pixel = 0;  //decouple pixel-pointer from previous pixel
                  //loop over segments, check if the pad received some charge
 
      for (Int_t i = 0; i < fNumberOfSegments; ++i) {
        //                      cout << endl << "check segment nr. " << i << " from " << fNumberOfSegments << endl;
        // ignore pads, which are out of reach for this segments
        if (ix < lowerXArray[i]) { continue; }
        if (iy < lowerYArray[i]) { continue; }
        if (ix > upperXArray[i]) { continue; }
        if (iy > upperYArray[i]) { continue; }
 
        // cout << endl << "found vallied pad " << i << endl;
        sPoint = &fSignalPoints[i];
 
        xCentre = sPoint->x;  //of segment
        yCentre = sPoint->y;  // idem
        sigmaX  = sPoint->sigmaX;
        sigmaY  = sPoint->sigmaY;
 
        fCurrentTotalCharge += sPoint->charge;
 
        //compute the charge distributed to this pixel by this segment
        Float_t totCharge =
          (sPoint->charge * fLorentzNorm * (0.5 * fPar0 * fPar1 / TMath::Pi())
             / TMath::Max(1.e-10,
                          (((Current[0] - xCentre) * (Current[0] - xCentre))
                           + ((Current[1] - yCentre) * (Current[1] - yCentre)))
                              / fPixelSize / fPixelSize
                            + 0.25 * fPar1 * fPar1)
           + fPar2);
 
        if (totCharge < 1) {
 
          //                     cout << endl << "charge is " << totCharge << " < 1 electron thus charge is negligible" << endl;
          continue;
        }  //ignore negligible charge (< 1 electron)
        if (!pixel) {
          // cout << endl << "charge is " << totCharge << " > 1 electron thus pixel is firred at "<< ix << " " << iy << endl;
          // Look for pixel in charge map if not yet linked.
          thispoint = std::make_pair(ix, iy);
          //                            cout << endl << "creat pair at "<< ix << " " << iy << endl;
          fChargeMapIt = fChargeMap.find(thispoint);
          //                            cout << endl << "found pair at "<< ix << " " << iy << endl;
          //                            cout << endl << "Map size is now " << fChargeMap.size() << endl;
          // Pixel not yet in map -> Add new pixel
          if (fChargeMapIt == fChargeMap.end()) {
            pixel = new ((*fPixelCharge)[fPixelCharge->GetEntriesFast()])
              CbmMvdPixelCharge(totCharge,
                                ix,
                                iy,
                                point->GetPointId(),
                                point->GetTrackID(),
                                (point->GetX() + point->GetXOut()) / 2,
                                (point->GetY() + point->GetXOut()) / 2,
                                point->GetTime(),
                                point->GetFrame());
            //cout << endl << "new charched pixel with charge " << totCharge << " at " << ix << " " << iy << endl;
            //                                            fPixelChargeShort.push_back(pixel);
            //                          cout << endl << "added pixel to ChargeShort vector " << endl;
 
            fChargeMap[thispoint] = pixel;
 
          }
 
          // Pixel already in map -> Add charge
          else {
            pixel = fChargeMapIt->second;
            //if ( ! pixel ) Fatal("AddChargeToPixel", "Zero pointer in charge map!");
            pixel->AddCharge(totCharge);
            //                                          if(pixel->GetCharge()>150)
            //                                  {cout << endl << "added charge to pixel summing up to "<< pixel->GetCharge() << endl;}
          }
          fPixelChargeShort.push_back(pixel);
        } else {  //pixel already linked => add charge only
          pixel->AddCharge(totCharge);
          //                            cout<<"put charge" << endl;
        }
 
 
        if (fShowDebugHistos) {
          xCharge          = xCharge + Current[0] * totCharge;
          yCharge          = yCharge + Current[1] * totCharge;
          totClusterCharge = totClusterCharge + totCharge;
        }  // end if
      }    // end for (track segments)
 
 
    }  //for y
  }    // for x
       //    cout << endl << "End of loops " << endl;
  std::vector<CbmMvdPixelCharge*>::size_type vectorSize =
    fPixelChargeShort.size();
 
  for (ULong64_t f = 0; f < vectorSize; f++) {
    CbmMvdPixelCharge* pixelCharge = fPixelChargeShort.at(f);
    if (pixelCharge) {
      pixelCharge->DigestCharge(
        ((float) (point->GetX() + point->GetXOut()) / 2),
        ((float) (point->GetY() + point->GetYOut()) / 2),
        point->GetPointId(),
        point->GetTrackID());
    } else {
      cout << endl << "Warning working on broken pixel " << endl;
    }
  }
 
 
  if (fShowDebugHistos) {
    //cout << endl << "produced " << fPixelChargeShort.size() << " Digis with total charge of " << totClusterCharge << endl;
    TVector3 momentum, position;
    if (totClusterCharge > 0) fTotalChargeHisto->Fill(totClusterCharge);
    point->Position(position);
    point->Momentum(momentum);
    fPosXY->Fill(position.X(), position.Y());
    fpZ->Fill(momentum.Z());
    fPosXinIOut->Fill(point->GetZ() - point->GetZOut());
    fAngle->Fill(TMath::ATan(momentum.Pt() / momentum.Pz()) * 180
                 / TMath::Pi());
  }
 
  delete[] lowerXArray;
  delete[] upperXArray;
  delete[] lowerYArray;
  delete[] upperYArray;
 
 
}  //end of function
 
 
// -------------------------------------------------------------------------
 
void CbmMvdSensorDigitizerTask::ProduceNoise() {
  Int_t fmaxNoise = 100;
  Int_t noiseHits = (Int_t)(frand->Rndm(fmaxNoise) * fmaxNoise);
  Int_t xPix, yPix;
  CbmMvdPixelCharge* pixel;
  pair<Int_t, Int_t> thispoint;
 
  for (Int_t i = 0; i <= noiseHits; i++) {
    xPix = frand->Integer(fNPixelsX);
    yPix = frand->Integer(fNPixelsY);
 
    Double_t Current[3];
    fSensor->PixelToLocal(xPix, yPix, Current);
 
    thispoint = std::make_pair(xPix, yPix);
 
    fChargeMapIt = fChargeMap.find(thispoint);
 
    if (fChargeMapIt == fChargeMap.end()) {
      pixel = new ((*fPixelCharge)[fPixelCharge->GetEntriesFast()])
        CbmMvdPixelCharge(1000, xPix, yPix, 0, -4, Current[0], Current[1]);
      pixel->DigestCharge(Current[0], Current[1], 0, -4);
      fChargeMap[thispoint] = pixel;
    } else {
      pixel = fChargeMapIt->second;
      pixel->AddCharge(1000);
      pixel->DigestCharge(Current[0], Current[1], 0, -4);
    }
  }
}
// -------------------------------------------------------------------------
 
// -------------------------------------------------------------------------
 
 
// -----   Virtual private method SetParContainers   -----------------------
void CbmMvdSensorDigitizerTask::SetParContainers() {
  //    FairRunAna*    ana  = FairRunAna::Instance();
  //    FairRuntimeDb* rtdb = ana->GetRuntimeDb();
}
// -------------------------------------------------------------------------
 
 
// -----    Virtual private method Init   ----------------------------------
void CbmMvdSensorDigitizerTask::InitTask(CbmMvdSensor* mySensor) {
 
  //Read information on the sensor von data base
  fSensor = mySensor;
 
  // cout << "-I- " << GetName() << ": Initialisation of sensor " << fSensor->GetName() << endl;
 
  fDigis     = new TClonesArray("CbmMvdDigi", 10000);
  fDigiMatch = new TClonesArray("CbmMatch", 10000);
 
  fOutputBuffer = new TClonesArray("CbmMvdDigi", 10000);
  fInputPoints  = new TClonesArray("CbmMvdPoint", 10000);
 
  if (!fSensor) {
    Fatal(GetName(),
          "Fatal error: Init(CbmMvdSensor*) called without valid pointer, "
          "don't know how to proceed.");
  };
 
  ReadSensorInformation();
 
  // Initialize histogramms used for debugging
 
  if (fShowDebugHistos) {
    cout << endl << "Show debug histos in this Plugin" << endl;
    fRandomGeneratorTestHisto =
      new TH1F("TestHisto", "TestHisto", 1000, 0, 12000);
    fResolutionHistoX =
      new TH1F("DigiResolutionX", "DigiResolutionX", 1000, -.005, .005);
    fResolutionHistoY =
      new TH1F("DigiResolutionY", "DigiResolutionY", 1000, -.005, .005);
    fPosXY = new TH2F("DigiPointXY", "DigiPointXY", 100, -6, 6, 100, -6, 6);
    fpZ    = new TH1F("DigiPointPZ", "DigiPointPZ", 1000, -0.5, 0.5);
    fPosXinIOut = new TH1F("DigiZInOut", "Digi Z InOut", 1000, -0.04, 0.04);
    fAngle      = new TH1F("DigiAngle", "DigiAngle", 1000, 0, 90);
    fSegResolutionHistoX =
      new TH1F("SegmentResolutionX", "SegmentResolutionX", 1000, -.005, .005);
    fSegResolutionHistoY =
      new TH1F("SegmentResolutionY", "SegmentResolutionY", 1000, -.005, .005);
    fSegResolutionHistoZ =
      new TH1F("SegmentResolutionZ", "SegmentResolutionZ", 1000, -.005, .005);
    fTotalChargeHisto =
      new TH1F("TotalChargeHisto", "TotalChargeHisto", 1000, 0, 12000);
    fTotalSegmentChargeHisto = new TH1F(
      "TotalSegmentChargeHisto", "TotalSegmentChargeHisto", 1000, 0, 12000);
  }
 
  fPreviousPlugin = NULL;
  initialized     = kTRUE;
}
// -------------------------------------------------------------------------
 
 
// -----   Virtual public method Reinit   ----------------------------------
void CbmMvdSensorDigitizerTask::ReInit(CbmMvdSensor* sensor) {
 
  delete fOutputBuffer;
 
  InitTask(sensor);
}
// -------------------------------------------------------------------------
 
 
// -----   Virtual method Finish   -----------------------------------------
void CbmMvdSensorDigitizerTask::Finish() {
 
 
  // PrintParameters();
 
 
  if (fShowDebugHistos) {
    TCanvas* c = new TCanvas("DigiCanvas", "DigiCanvas", 150, 10, 800, 600);
    c->Divide(3, 3);
    c->cd(1);
    fResolutionHistoX->Draw();
    fResolutionHistoX->Write();
    c->cd(2);
    fResolutionHistoY->Draw();
    fResolutionHistoY->Write();
    c->cd(3);
    fPosXY->Draw();
    fPosXY->Write();
    c->cd(4);
    fpZ->Draw();
    fpZ->Write();
    c->cd(5);
    fPosXinIOut->Draw();
    fPosXinIOut->Write();
    c->cd(6);
    fAngle->Draw();
    fAngle->Write();
    c->cd(7);
    //fSegResolutionHistoX->Draw();
    fSegResolutionHistoX->Write();
    fTotalSegmentChargeHisto->Draw();
    fTotalSegmentChargeHisto->Write();
    c->cd(8);
    fRandomGeneratorTestHisto->Draw();
    fRandomGeneratorTestHisto->Write();
 
    fSegResolutionHistoY->Write();
    c->cd(9);
    fTotalChargeHisto->Draw();
    fTotalChargeHisto->Write();
    cout << "-I- CbmMvdDigitizerL::Finish - Fit of the total cluster charge"
         << endl;
    fTotalChargeHisto->Fit("landau");
    cout << "=============================================================="
         << endl;
    // new TCanvas();
    //fTotalChargeHisto->Draw();
  };
}
// -------------------------------------------------------------------------
 
 
// -----   Private method Reset   ------------------------------------------
void CbmMvdSensorDigitizerTask::Reset() {}
// -------------------------------------------------------------------------
 
 
// -----   Private method PrintParameters   --------------------------------
void CbmMvdSensorDigitizerTask::PrintParameters() {
 
  cout.setf(ios_base::fixed, ios_base::floatfield);
  cout << "============================================================"
       << endl;
  cout << "============== Parameters of the Lorentz - Digitizer ======="
       << endl;
  cout << "============================================================"
       << endl;
 
 
  cout << "Pixel Size X               : " << setw(8) << setprecision(2)
       << fPixelSizeX * 10000. << " mum" << endl;
  cout << "Pixel Size Y               : " << setw(8) << setprecision(2)
       << fPixelSizeY * 10000. << " mum" << endl;
  cout << "Epitaxial layer thickness  : " << setw(8) << setprecision(2)
       << fEpiTh * 10000. << " mum" << endl;
  cout << "Segment Length             : " << setw(8) << setprecision(2)
       << fSegmentLength * 10000. << " mum" << endl;
  cout << "Diffusion Coefficient      : " << setw(8) << setprecision(2)
       << fDiffusionCoefficient * 10000. << " mum" << endl;
  cout << "Width of Cluster           : " << setw(8) << setprecision(2)
       << fWidthOfCluster << " * sigma " << endl;
  cout << "ElectronsPerKeV 3.62 eV/eh : " << setw(8) << setprecision(2)
       << fElectronsPerKeV << endl;
  cout << "CutOnDeltaRays             : " << setw(8) << setprecision(8)
       << fCutOnDeltaRays << " MeV " << endl;
  cout << "ChargeThreshold            : " << setw(8) << setprecision(2)
       << fChargeThreshold << endl;
  cout << "Pileup: " << fNPileup << endl;
  cout << "Delta - Pileup: " << fNDeltaElect << endl;
  cout << "=============== End Parameters Digitizer ==================="
       << endl;
}
// -------------------------------------------------------------------------
 
 
ClassImp(CbmMvdSensorDigitizerTask);