ROOT logo
///////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                           //
//                                                                                           //
//     DIGParticle                                                                           //
//                                                                                           //
//     particle class                                                                        //
//       contains: -entry and exit point of the particle into the plane                      //
//                 -segment list (Charge, position) of the track                             //
//                 -pixel list (number, charge) where charge has been collected              //
//                                                                                           //
//                                                                                           //
///////////////////////////////////////////////////////////////////////////////////////////////
#include <digparticle.h>

#include <TROOT.h> // for gROOT object
#include <TMath.h>
#include <TMatrixD.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TAxis.h>
#include <TRandom3.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TClonesArray.h>
#include <TF2.h>
#include <TProfile2D.h>
#include <TH2.h>


#include "digplane.h"
#include "digadc.h"
#include "digtransport.h"

using namespace std;

extern Int_t GlobalSeed;

//_______________________________________________________________________________________
//
  Double_t Lorentz2D(Double_t *x, Double_t *par){ 
    //x[0] = x
    //x[1] = y
    // par[0] = Gamma
    // par[1] = x0
    // par[2] = y0
    // par[3] = norm
  if(par[0]>0){
    Double_t Pi = 3.141592653;
    return par[3]*par[0]/Pi/((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[2])*(x[1]-par[2])+par[0]*par[0]) ; 
  }else{
    return 0;
  }
}
//_______________________________________________________________________________________
//
Double_t SumGaus2D(Double_t *x, Double_t *par){
  //par[0] Norm_1
  //par[1] x0_1
  //par[2] sigma_x_1
  //par[3] y0_1
  //par[4] sigma_y_1
  //par[5] Norm_2
  //par[6] x0_2
  //par[7] sigma_x_2
  //par[8] y0_2
  //par[9] sigma_y_2
  if((par[2]!=0.0) && (par[4]!=0.0) && (par[7]!=0.0) && (par[9]!=0.0)){
    double rx = (x[0]-par[1])/par[2];
    double ry = (x[1]-par[3])/par[4];
    double rx2 = (x[0]-par[6])/par[7];
    double ry2 = (x[1]-par[8])/par[9];
    return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)+par[5]*TMath::Exp(-(rx2*rx2+ry2*ry2)/2.0) ) ;
  }else{
    return par[0]+par[1]+par[2]+par[3]+par[4]+par[5]+par[6]+par[7]+par[8]+par[9];
    //return 0;
  }
}
//_______________________________________________________________________________________
//
Double_t SumGausLorentz2D(Double_t *x, Double_t *par){
  //par[0] Norm_1
  //par[1] x0_1
  //par[2] sigma_x_1
  //par[3] y0_1
  //par[4] sigma_y_1
  // par[5] = Gamma
  // par[6] = x0
  // par[7] = y0
  // par[8] = norm
  //0 0.146571+-0 0+-0 7.55129+-0 0+-0 5.94133+-0 3.2357+-0 0+-0 0+-0 27.3458+-0 chi2 / NDF = 3126.16 / 10009 = 0.312335
  //1 0.1459+-0 0+-0 14.7913+-0 0+-0 11.7393+-0 6.86661+-0 0+-0 0+-0 58.7235+-0 chi2 / NDF = 7458.34 / 10009 = 0.745163
  //2 0.149045+-0 0+-0 22.2327+-0 0+-0 17.6162+-0 9.20632+-0 0+-0 0+-0 90.0768+-0 chi2 / NDF = 12700.5 / 10009 = 1.26891
  //3 0.168693+-0 0+-0 27.9107+-0 0+-0 21.8062+-0 12.9946+-0 0+-0 0+-0 102.019+-0 chi2 / NDF = 8275.98 / 10009 = 0.826854
  Double_t Pi = 3.141592653;
  if((par[2]!=0.0) && (par[4]!=0.0) ){
    double rx = (x[0]-par[1])/par[2];
    double ry = (x[1]-par[3])/par[4];
    return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)
    		    +par[8]*par[5]/Pi/((x[0]-par[6])*(x[0]-par[6])+(x[1]-par[7])*(x[1]-par[7])+par[5]*par[5])
    		    );
  }else{
    return 0;  
  }
}
//==============================================================================
ClassImp(DIGParticle)
DIGParticle::DIGParticle()  
{
  //
  // default constructor
  //

}
//______________________________________________________________________________
//  
DIGParticle::DIGParticle(Float_t EntryX, Float_t EntryY, Float_t EntryZ, 
			 Float_t ExitX, Float_t ExitY, Float_t ExitZ, Float_t Energy_deposited)
{
  fEntryX = EntryX;
  fEntryY = EntryY;
  fEntryZ = EntryZ;
  fExitX = ExitX ;
  fExitY = ExitY ;
  fExitZ = ExitZ ;
  fEnergy_deposited = Energy_deposited;
  fNSegment =0;
  fSegmentX.clear();
  fSegmentY.clear();
  fSegmentZ.clear();
  fSegmentCharge.clear();

  fNpixels=0;
  fPixelMap.clear();
  fAnalogChargeMap.clear();
  fDigitalChargeMap.clear();//not used ?
 
}
//______________________________________________________________________________
//  
DIGParticle::DIGParticle(DIGParticle & adigparticle)  : TObject()
{
  fEntryX = adigparticle.GetEntryX();
  fEntryY = adigparticle.GetEntryY();
  fEntryZ = adigparticle.GetEntryZ();
  fExitX = adigparticle.GetExitX();
  fExitY = adigparticle.GetExitY();
  fExitZ = adigparticle.GetExitZ();
  fEnergy_deposited = adigparticle.GetEnergy_deposited();
  fNSegment =adigparticle.GetNSegment();
  fSegmentX.resize(fNSegment);
  fSegmentY.resize(fNSegment);
  fSegmentZ.resize(fNSegment);
  fSegmentCharge.resize(fNSegment);
  // a modifier:
  for (Int_t i=0 ; i<fNSegment ; i++){
    fSegmentX[i] = adigparticle.GetSegmentX()[i];
    fSegmentY[i] = adigparticle.GetSegmentY()[i];
    fSegmentZ[i] = adigparticle.GetSegmentZ()[i];
    fSegmentCharge[i] = adigparticle.GetSegmentCharge()[i];
  }
  fNpixels = adigparticle.GetNpixels();
  fPixelMap.resize(fNpixels);
  fAnalogChargeMap.resize(fNpixels);
  fDigitalChargeMap.resize(fNpixels);
  for (Int_t i=0 ; i<fNpixels ; i++){
    fPixelMap[i] = adigparticle.GetPixelMap()[i];
    fAnalogChargeMap[i] = adigparticle.GetAnalogCharge()[i];
    fDigitalChargeMap[i] = adigparticle.GetDigitalCharge()[i];
  }
}
//______________________________________________________________________________
//  


DIGParticle::~DIGParticle() { // 
  // virtual destructor
  //
  //  delete fLayers;
}
//______________________________________________________________________________
//  
void DIGParticle::Clear(const Option_t *) 
{

}
//______________________________________________________________________________
//  
void DIGParticle::ComputeChargeDeposition(Float_t StartingSegmentSize, Float_t MaximumSegmentSize,
					  Float_t MaximumChargePerSegment) 
{
  Float_t SegmentSize = StartingSegmentSize;
  Float_t TotalLength = TMath::Sqrt((fExitX-fEntryX)*(fExitX-fEntryX)
				    +(fExitY-fEntryY)*(fExitY-fEntryY)
				    +(fExitZ-fEntryZ)*(fExitZ-fEntryZ));

  Float_t ChargePerSegment = 0.0;
  if(SegmentSize<0.0){
    SegmentSize=0.001;
  }
  fNSegment = int(TotalLength*1.000001/SegmentSize) ;
  if(fNSegment<1){
    fNSegment=1;
  }
  SegmentSize = TotalLength/float(fNSegment);
  ChargePerSegment = fEnergy_deposited / float(fNSegment);

  while((SegmentSize>MaximumSegmentSize)&&(ChargePerSegment > MaximumChargePerSegment)){
    Int_t newNSegment = int(fNSegment *1.1);
    if(newNSegment==fNSegment){fNSegment++;}
    SegmentSize = TotalLength/float(fNSegment);
    ChargePerSegment = fEnergy_deposited / float(fNSegment);
  }
  Float_t xstep = fExitX-fEntryX;
  Float_t ystep = fExitY-fEntryY;
  Float_t zstep = fExitZ-fEntryZ;


  for (Int_t i=0 ; i<fNSegment ; i++){
    fSegmentX.push_back(fEntryX + (float(i+0.5)* xstep/float(fNSegment)) );
    fSegmentY.push_back(fEntryY + (float(i+0.5)* ystep/float(fNSegment)) );
    fSegmentZ.push_back(fEntryZ + (float(i+0.5)* zstep/float(fNSegment)) );
    fSegmentCharge.push_back(ChargePerSegment );
  }
    

}
//______________________________________________________________________________
//  

void DIGParticle::PrintInfo() {
  std::cout<<"---------Particle properties------------- "<<endl;
  std::cout<<"fEntryX fEntryY fEntryZ"<<endl;
  std::cout<<fEntryX<<" "<< fEntryY<<" "<<fEntryZ<<endl;
  std::cout<<"fExitX fExitY fExitZ fEnergy_deposited"<<endl;
  std::cout<<fExitX <<" "<<fExitY <<" "<<fExitZ <<" "<<fEnergy_deposited <<endl;
  std::cout<<fNSegment<<" Segments X Y Z Charge"<<endl;
  /*  for (Int_t i=0 ; i<fNSegment ; i++){
    std::cout<<i<<" "<< fSegmentX[i]<<" "<<fSegmentY[i] <<" "<<fSegmentZ[i] <<" "<< fSegmentCharge[i]<<endl;
    }*/
  cout<<" size vectors "<<  fPixelMap.size()<<" "<<fAnalogChargeMap.size()<<" "<<fDigitalChargeMap.size()<<endl;
  std::cout<<fNpixels<<" fNpixels map analog digital "<<endl;
  for (Int_t i=0 ; i<fNpixels ; i++){
    std::cout<<i<<" "<<fPixelMap[i]<<" "<<fAnalogChargeMap[i]<<" "<<fDigitalChargeMap[i]<<endl;
  }
}
//______________________________________________________________________________
//   
void DIGParticle::ComputeChargeTransport(DIGPlane *aDIGPlane,DIGTransport *aDIGTransport){

  Float_t pitchX =     aDIGPlane->GetPitchX();
  Float_t pitchY =     aDIGPlane->GetPitchY();

  Float_t lorentz2Dmodel_Cp0 = aDIGTransport->GetLorentz2DModel_Cp0();
  Float_t lorentz2Dmodel_Cp1 = aDIGTransport->GetLorentz2DModel_Cp1();
  Float_t rangelimit_inpitchunit = aDIGTransport->GetRangeLimit_InPitchUnit();

  Int_t ChargeModel = aDIGTransport->GetChargeModel();

  Float_t Gauss2DModel_sigma1_Cp0 = aDIGTransport->GetGauss2DModel_sigma1_Cp0();
  Float_t Gauss2DModel_sigma1_Cp1 = aDIGTransport->GetGauss2DModel_sigma1_Cp1();
  Float_t Gauss2DModel_sigma2_Cp0 = aDIGTransport->GetGauss2DModel_sigma2_Cp0();
  Float_t Gauss2DModel_sigma2_Cp1 = aDIGTransport->GetGauss2DModel_sigma2_Cp1();
  Float_t Gauss2DModel_weight = aDIGTransport->GetGauss2DModel_weight();

  Float_t LorGaussModel_Norm1_Cp0  = aDIGTransport->GetLorGaussModel_Norm1_Cp0();
  Float_t LorGaussModel_Norm1_Cp1 = aDIGTransport->GetLorGaussModel_Norm1_Cp1();
  Float_t LorGaussModel_Norm1_Cp2 = aDIGTransport->GetLorGaussModel_Norm1_Cp2();
  Float_t LorGaussModel_sigma_Cp0 = aDIGTransport->GetLorGaussModel_sigma_Cp0();
  Float_t LorGaussModel_sigma_Cp1 = aDIGTransport->GetLorGaussModel_sigma_Cp1();
  Float_t LorGaussModel_C_Cp0 = aDIGTransport->GetLorGaussModel_C_Cp0();
  Float_t LorGaussModel_C_Cp1 = aDIGTransport->GetLorGaussModel_C_Cp1();
  Float_t LorGaussModel_Norm_Cp0 = aDIGTransport->GetLorGaussModel_Norm_Cp0();
  Float_t LorGaussModel_Norm_Cp1 = aDIGTransport->GetLorGaussModel_Norm_Cp1();


  GlobalSeed++;
  TRandom3 *r3 = new TRandom3(GlobalSeed);
  r3->SetSeed(GlobalSeed);

  TF2 *mymodel2D;
  //-----------------------------------------------------
  // chose model
  //-----------------------------------------------------
  if(ChargeModel==1){
    //-----------------------------------------------------
    // model Lorentz
    //-----------------------------------------------------

    Double_t Cvalue=  lorentz2Dmodel_Cp0 + pitchX * lorentz2Dmodel_Cp1 ;

    mymodel2D = new TF2("funlor2d",Lorentz2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
			  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,4);
    mymodel2D->SetParNames("C","X_{0}","Y_{0}","Norm");
    Double_t params1[] = {Cvalue,0.,0.,1.0};
    mymodel2D->SetParameters(params1);

 }else if(ChargeModel==2){
    //-----------------------------------------------------
    // 2xGaussian model 
    //-----------------------------------------------------
    //Double_t SumGaus2D(Double_t *x, Double_t *par){
    //par[0] Norm_1
    //par[1] x0_1
    //par[2] sigma_x_1
    //par[3] y0_1
    //par[4] sigma_y_1
    //par[5] Norm_2
    //par[6] x0_2
    //par[7] sigma_x_2
    //par[8] y0_2
    //par[9] sigma_y_2
    Double_t Gsigma1 = Gauss2DModel_sigma1_Cp0  + pitchX * Gauss2DModel_sigma1_Cp1 ;
    Double_t Gsigma2 = Gauss2DModel_sigma2_Cp0  + pitchX * Gauss2DModel_sigma2_Cp1 ;

    mymodel2D = new TF2("funlor2d",SumGaus2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
			  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,10);
    mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Norm_2","x0_2","sigma_x_2","y0_2","sigma_y_2");
    Double_t params1[] = {1.0,0.0,Gsigma1,0.0,Gsigma1,Gauss2DModel_weight,0.0,Gsigma2,0.0,Gsigma2,0.0};
    mymodel2D->SetParameters(params1);
  }else if(ChargeModel==3){
    //-----------------------------------------------------
    // Lorentz + Gauss model 
    //-----------------------------------------------------
    //par[0] Norm_1
    //par[1] x0_1
    //par[2] sigma_x_1
    //par[3] y0_1
    //par[4] sigma_y_1
    // par[5] = Gamma
    // par[6] = x0
    // par[7] = y0
    // par[8] = norm    
    Double_t Norm1Value =  LorGaussModel_Norm1_Cp0 + pitchX * LorGaussModel_Norm1_Cp1 + pitchX *pitchX * LorGaussModel_Norm1_Cp2;
    Double_t sigmaValue = LorGaussModel_sigma_Cp0 + pitchX *LorGaussModel_sigma_Cp1;
    Double_t Cvalue= LorGaussModel_C_Cp0 + pitchX *LorGaussModel_C_Cp1;
    Double_t NormValue =LorGaussModel_Norm_Cp0 + pitchX * LorGaussModel_Norm_Cp1;
    mymodel2D = new TF2("funlor2d",SumGausLorentz2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
			  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,9);
    mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Gamma","x0","y0","norm");
    Double_t params1[] = {Norm1Value,0.0,sigmaValue,0.0,sigmaValue,Cvalue,0.0,0.0,NormValue};
    mymodel2D->SetParameters(params1);


  }else{
    Double_t Gsigma1 = Gauss2DModel_sigma1_Cp0  + pitchX * Gauss2DModel_sigma1_Cp1 ;
    Double_t Gsigma2 = Gauss2DModel_sigma2_Cp0  + pitchX * Gauss2DModel_sigma2_Cp1 ;

    mymodel2D = new TF2("funlor2d",SumGaus2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
			  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,10);
    mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Norm_2","x0_2","sigma_x_2","y0_2","sigma_y_2");
    Double_t params1[] = {1.0,0.0,Gsigma1,0.0,Gsigma1,Gauss2DModel_weight,0.0,Gsigma2,0.0,Gsigma2,0.0};
    mymodel2D->SetParameters(params1);
 
  }


    //---------------------------------------------------------------------
    // method A: lorentz function define a probability density function which gives where the charge is collected.
    //---------------------------------------------------------------------
    /*
      Float_t Xdimension =     aDIGPlane->GetXdimension();
      Float_t Ydimension =     aDIGPlane->GetYdimension();
  
      for (Int_t i=0 ; i<fNSegment ; i++){
      //std::cout<<" i="<<i<<" x="<< fSegmentX[i]<<" y="<<fSegmentY[i] <<" z="<<fSegmentZ[i] <<" c="<< fSegmentCharge[i]<<endl;
      Double_t xrandom;
      Double_t yrandom;
      gRandom = new TRandom3;
      GlobalSeed++;
      gRandom->SetSeed(GlobalSeed);
      mylorentz2D->GetRandom2(xrandom,yrandom);
    
      Double_t xpos_aftertransport = fSegmentX[i] + xrandom;
      Double_t ypos_aftertransport = fSegmentY[i] + yrandom;
      //if charge reaches a x,y position outside the size of the detector, consider that the charge is lost.
      if((xpos_aftertransport>=0.0)
      &&(xpos_aftertransport<=Xdimension)
      &&(ypos_aftertransport>=0.0)
      &&(ypos_aftertransport<=Ydimension)){
      Int_t PixelReached =  GetPixelNumber(aDIGPlane, xpos_aftertransport ,  ypos_aftertransport );
      UpdatePixel(fSegmentCharge[i],PixelReached);	  
      }else{
      //  Int_t PixelReached =  GetPixelNumber(aDIGPlane, xpos_aftertransport ,  ypos_aftertransport );
      //     Int_t xnpos, ynpos;
      //      GetXYPixelNumber(xnpos,ynpos, aDIGPlane, PixelReached);
      //      Float_t  xdpos, ydpos;    
      //      GetXYPixelCenter(xdpos,ydpos, aDIGPlane, PixelReached);
      //      std::cout<<"charge outside DIMENSIONS "<<Xdimension <<" x "<<Ydimension <<endl;
      //      std::cout<<" after ranx="<< xrandom<<" rany="<<yrandom<<" x="<<xpos_aftertransport<<" y="<<ypos_aftertransport
      //      <<" n="<<PixelReached<<" xn= "<<xnpos 
      //      <<" yn="<<ynpos <<" xd="<<xdpos<<" yd="<<ydpos<<endl;
      
      }
      }
    */	       
    //---------------------------------------------------------------------
    // method B: more realistic model. The lorentz function allows to compute 25 probabilities (5x5 cluster)
    // the  random number is generated to see where the segment charge is deposited.
    //---------------------------------------------------------------------

    for (Int_t i=0 ; i<fNSegment ; i++){      
    
      Int_t PixelReached =  GetPixelNumber(aDIGPlane, fSegmentX[i], fSegmentY[i]);
      Float_t  xdpos, ydpos;
      GetXYPixelCenter(xdpos, ydpos, aDIGPlane, PixelReached);
      Int_t xpixn,ypixn;
      GetXYPixelNumber( xpixn,ypixn, aDIGPlane, PixelReached);

      Float_t TotalProba=0.0;
      const Int_t Npix=25;
      Float_t PixelposX[Npix];
      Float_t PixelposY[Npix];
      Float_t Pixelproba[Npix];
      Int_t Pixelnumber[Npix];
      for (Int_t j=0 ; j<Npix ; j++){
	PixelposX[j]=0.0;
	PixelposY[j]=0.0;
	Pixelproba[j]=0.0;
	Pixelnumber[j]=0;
      }
    
      for (Int_t j=0 ; j<Npix ; j++){
	Int_t xi = j%5;
	Int_t yi = j/5;     
	Float_t xeval = -(fSegmentX[i]-(xpixn+0.5)*pitchX) + (float(xi)-2.0)*pitchX ;
	Float_t yeval = -(fSegmentY[i]-(ypixn+0.5)*pitchY) + (float(yi)-2.0)*pitchY ;

	// force Lorentz function to be centered in the seed pixel for tests:
	//Float_t xeval = (float(xi)-2.0)*pitchX ;
	//Float_t yeval = (float(yi)-2.0)*pitchY ;


	Pixelproba[j]=mymodel2D->Eval(xeval,yeval);
	TotalProba+=Pixelproba[j];
	PixelposX[j]= xdpos + (float(xi)-2.0)*pitchX;
	PixelposY[j]= ydpos + (float(yi)-2.0)*pitchY;
	if(  (PixelposX[j]<=0.0)||
	     (PixelposX[j]>=aDIGPlane->GetXdimension())||
	     (PixelposY[j]<=0.0)||
	     (PixelposY[j]>=aDIGPlane->GetYdimension())){
	  Pixelnumber[j] = -1;
	}else{
	  Pixelnumber[j] = GetPixelNumber(aDIGPlane, PixelposX[j], PixelposY[j]);
	}
      }
    
      Double_t rando = r3->Rndm()*TotalProba;
      Bool_t reached = false;
      Double_t incrprob = 0.0;
      Int_t icounter=0;
      while((!reached)&&(icounter<Npix)){
	incrprob+=Pixelproba[icounter];
	if(incrprob>rando){
	  reached=true;
	}else{
	  icounter++;
	}
      }

      if(Pixelnumber[icounter]>=0){
	UpdatePixel(fSegmentCharge[i],Pixelnumber[icounter]);
      }
    }
  delete r3;
}
//______________________________________________________________________________
//   

void DIGParticle::AddPixel(Float_t AnalogCharge, Int_t PixelNumber){
  fNpixels++;
  fPixelMap.push_back(PixelNumber);
  fAnalogChargeMap.push_back(AnalogCharge);
  fDigitalChargeMap.push_back(0);
}

//______________________________________________________________________________
//   
void DIGParticle::UpdatePixel(Float_t AnalogCharge, Int_t PixelNumber){
  Bool_t found = false;
  Int_t i=0;
  while((!found)&&(i<fNpixels)){
    if(PixelNumber==fPixelMap[i]){     
      found=true;
      fAnalogChargeMap[i]+=AnalogCharge;
    }
    i++;
  }
  if(!found){
    AddPixel(AnalogCharge, PixelNumber);
  } 

}
//______________________________________________________________________________
//   
void DIGParticle::AddRandomNoise(DIGPlane *myDIGPlane){
  Double_t Noisefix = myDIGPlane->GetNoiseElectrons();
  Double_t Noise;
  for (Int_t i = 0; i < fNpixels ; i++){
    Noise =   GaussianLaw(0.0, Noisefix);
    fAnalogChargeMap[i]+=Noise;
  }
}
//______________________________________________________________________________
//   
void DIGParticle::AnalogToDigitalconversion(DIGADC *myDIGADC, DIGPlane *myDIGPlane ){
  Float_t Noisefix = myDIGPlane->GetNoiseElectrons();
  if(Noisefix<=0.0){
    std::cout <<"<---- DIGParticle::AnalogToDigitalconversion --->"<<endl;
    std::cout<<"WARNING negative or null Noise is not physical, please correct the input file"<<endl; 
    Noisefix = 1.0;
  }
  Float_t *myADCthreshold = 0;
  myADCthreshold = myDIGADC->GetADC_thresholds();
  for (Int_t i = 0; i < fNpixels ; i++){
    if (fAnalogChargeMap[i]<=0.0){
      fDigitalChargeMap[i]=0;
    }else{
      Bool_t thresholdfound = false;
      Int_t ithres = 0;      
      fDigitalChargeMap[i]=0;
      while((thresholdfound==false)&&(ithres< (myDIGADC->GetNThresholds()) )){
	if( (fAnalogChargeMap[i]/Noisefix) < myADCthreshold[ithres] ){
	  thresholdfound = true;
	}else{
	  fDigitalChargeMap[i]++;
	  ithres++;
	}
      }
    }
  }
}
//______________________________________________________________________________
//   
Int_t DIGParticle::GetPixelNumber(DIGPlane *myDIGPlane,  Float_t Xpos, Float_t Ypos){
  Int_t XPixelnumber = int(Xpos/(myDIGPlane->GetPitchX()));
  Int_t YPixelnumber = int(Ypos/(myDIGPlane->GetPitchY()));
  Int_t PixelNumber = XPixelnumber + (myDIGPlane->GetNpixelsX())*YPixelnumber;
 
  if(
     (Xpos<0.0)||
     (Xpos>myDIGPlane->GetXdimension())||
     (Ypos<0.0)||
     (Ypos>myDIGPlane->GetYdimension())){
    cout<<" DIGParticle::GetPixelNumber WARNING  charge is going outside the plane limits"<<endl;
    return 0;
  }else{
    return PixelNumber;
  }
}
//______________________________________________________________________________
//   
void DIGParticle::GetXYPixelNumber(Int_t &Xpix, Int_t &Ypix, DIGPlane *myDIGPlane, Int_t PixelNumber){

  Xpix = PixelNumber%(myDIGPlane->GetNpixelsX());
  Ypix = PixelNumber/(myDIGPlane->GetNpixelsX());

}
//______________________________________________________________________________
//   
void DIGParticle::GetXYPixelCenter(Float_t &Xpix, Float_t &Ypix, DIGPlane *myDIGPlane, Int_t PixelNumber){

  Xpix = (myDIGPlane->GetPitchY()) * (0.5+PixelNumber%(myDIGPlane->GetNpixelsX()));
  Ypix = (myDIGPlane->GetPitchX()) * (0.5+PixelNumber/(myDIGPlane->GetNpixelsX()));

}
//______________________________________________________________________________
//  
Double_t  DIGParticle::GaussianLaw(Double_t mean, Double_t sigma) 
{ 
  Double_t x;
  TRandom3 *r3 = new TRandom3(0);
  GlobalSeed++;
  r3->SetSeed(GlobalSeed);

  x = r3->Gaus(mean,sigma);

  delete r3;
  return x;
}
//______________________________________________________________________________
//  
Float_t DIGParticle::GetTotalAnalogCharge(){
  Float_t TotalCharge = 0;
  for (Int_t i=0 ; i<fNpixels ; i++){
    TotalCharge+=fAnalogChargeMap[i];
  }
  return TotalCharge;
}
//______________________________________________________________________________
//  
Int_t DIGParticle::GetTotalDigitalCharge(){
  Int_t TotalCharge = 0;
  for (Int_t i=0 ; i<fNpixels ; i++){
    TotalCharge+=fDigitalChargeMap[i];
  }
  return TotalCharge;
}
//==============================================================================
 digparticle.cxx:1
 digparticle.cxx:2
 digparticle.cxx:3
 digparticle.cxx:4
 digparticle.cxx:5
 digparticle.cxx:6
 digparticle.cxx:7
 digparticle.cxx:8
 digparticle.cxx:9
 digparticle.cxx:10
 digparticle.cxx:11
 digparticle.cxx:12
 digparticle.cxx:13
 digparticle.cxx:14
 digparticle.cxx:15
 digparticle.cxx:16
 digparticle.cxx:17
 digparticle.cxx:18
 digparticle.cxx:19
 digparticle.cxx:20
 digparticle.cxx:21
 digparticle.cxx:22
 digparticle.cxx:23
 digparticle.cxx:24
 digparticle.cxx:25
 digparticle.cxx:26
 digparticle.cxx:27
 digparticle.cxx:28
 digparticle.cxx:29
 digparticle.cxx:30
 digparticle.cxx:31
 digparticle.cxx:32
 digparticle.cxx:33
 digparticle.cxx:34
 digparticle.cxx:35
 digparticle.cxx:36
 digparticle.cxx:37
 digparticle.cxx:38
 digparticle.cxx:39
 digparticle.cxx:40
 digparticle.cxx:41
 digparticle.cxx:42
 digparticle.cxx:43
 digparticle.cxx:44
 digparticle.cxx:45
 digparticle.cxx:46
 digparticle.cxx:47
 digparticle.cxx:48
 digparticle.cxx:49
 digparticle.cxx:50
 digparticle.cxx:51
 digparticle.cxx:52
 digparticle.cxx:53
 digparticle.cxx:54
 digparticle.cxx:55
 digparticle.cxx:56
 digparticle.cxx:57
 digparticle.cxx:58
 digparticle.cxx:59
 digparticle.cxx:60
 digparticle.cxx:61
 digparticle.cxx:62
 digparticle.cxx:63
 digparticle.cxx:64
 digparticle.cxx:65
 digparticle.cxx:66
 digparticle.cxx:67
 digparticle.cxx:68
 digparticle.cxx:69
 digparticle.cxx:70
 digparticle.cxx:71
 digparticle.cxx:72
 digparticle.cxx:73
 digparticle.cxx:74
 digparticle.cxx:75
 digparticle.cxx:76
 digparticle.cxx:77
 digparticle.cxx:78
 digparticle.cxx:79
 digparticle.cxx:80
 digparticle.cxx:81
 digparticle.cxx:82
 digparticle.cxx:83
 digparticle.cxx:84
 digparticle.cxx:85
 digparticle.cxx:86
 digparticle.cxx:87
 digparticle.cxx:88
 digparticle.cxx:89
 digparticle.cxx:90
 digparticle.cxx:91
 digparticle.cxx:92
 digparticle.cxx:93
 digparticle.cxx:94
 digparticle.cxx:95
 digparticle.cxx:96
 digparticle.cxx:97
 digparticle.cxx:98
 digparticle.cxx:99
 digparticle.cxx:100
 digparticle.cxx:101
 digparticle.cxx:102
 digparticle.cxx:103
 digparticle.cxx:104
 digparticle.cxx:105
 digparticle.cxx:106
 digparticle.cxx:107
 digparticle.cxx:108
 digparticle.cxx:109
 digparticle.cxx:110
 digparticle.cxx:111
 digparticle.cxx:112
 digparticle.cxx:113
 digparticle.cxx:114
 digparticle.cxx:115
 digparticle.cxx:116
 digparticle.cxx:117
 digparticle.cxx:118
 digparticle.cxx:119
 digparticle.cxx:120
 digparticle.cxx:121
 digparticle.cxx:122
 digparticle.cxx:123
 digparticle.cxx:124
 digparticle.cxx:125
 digparticle.cxx:126
 digparticle.cxx:127
 digparticle.cxx:128
 digparticle.cxx:129
 digparticle.cxx:130
 digparticle.cxx:131
 digparticle.cxx:132
 digparticle.cxx:133
 digparticle.cxx:134
 digparticle.cxx:135
 digparticle.cxx:136
 digparticle.cxx:137
 digparticle.cxx:138
 digparticle.cxx:139
 digparticle.cxx:140
 digparticle.cxx:141
 digparticle.cxx:142
 digparticle.cxx:143
 digparticle.cxx:144
 digparticle.cxx:145
 digparticle.cxx:146
 digparticle.cxx:147
 digparticle.cxx:148
 digparticle.cxx:149
 digparticle.cxx:150
 digparticle.cxx:151
 digparticle.cxx:152
 digparticle.cxx:153
 digparticle.cxx:154
 digparticle.cxx:155
 digparticle.cxx:156
 digparticle.cxx:157
 digparticle.cxx:158
 digparticle.cxx:159
 digparticle.cxx:160
 digparticle.cxx:161
 digparticle.cxx:162
 digparticle.cxx:163
 digparticle.cxx:164
 digparticle.cxx:165
 digparticle.cxx:166
 digparticle.cxx:167
 digparticle.cxx:168
 digparticle.cxx:169
 digparticle.cxx:170
 digparticle.cxx:171
 digparticle.cxx:172
 digparticle.cxx:173
 digparticle.cxx:174
 digparticle.cxx:175
 digparticle.cxx:176
 digparticle.cxx:177
 digparticle.cxx:178
 digparticle.cxx:179
 digparticle.cxx:180
 digparticle.cxx:181
 digparticle.cxx:182
 digparticle.cxx:183
 digparticle.cxx:184
 digparticle.cxx:185
 digparticle.cxx:186
 digparticle.cxx:187
 digparticle.cxx:188
 digparticle.cxx:189
 digparticle.cxx:190
 digparticle.cxx:191
 digparticle.cxx:192
 digparticle.cxx:193
 digparticle.cxx:194
 digparticle.cxx:195
 digparticle.cxx:196
 digparticle.cxx:197
 digparticle.cxx:198
 digparticle.cxx:199
 digparticle.cxx:200
 digparticle.cxx:201
 digparticle.cxx:202
 digparticle.cxx:203
 digparticle.cxx:204
 digparticle.cxx:205
 digparticle.cxx:206
 digparticle.cxx:207
 digparticle.cxx:208
 digparticle.cxx:209
 digparticle.cxx:210
 digparticle.cxx:211
 digparticle.cxx:212
 digparticle.cxx:213
 digparticle.cxx:214
 digparticle.cxx:215
 digparticle.cxx:216
 digparticle.cxx:217
 digparticle.cxx:218
 digparticle.cxx:219
 digparticle.cxx:220
 digparticle.cxx:221
 digparticle.cxx:222
 digparticle.cxx:223
 digparticle.cxx:224
 digparticle.cxx:225
 digparticle.cxx:226
 digparticle.cxx:227
 digparticle.cxx:228
 digparticle.cxx:229
 digparticle.cxx:230
 digparticle.cxx:231
 digparticle.cxx:232
 digparticle.cxx:233
 digparticle.cxx:234
 digparticle.cxx:235
 digparticle.cxx:236
 digparticle.cxx:237
 digparticle.cxx:238
 digparticle.cxx:239
 digparticle.cxx:240
 digparticle.cxx:241
 digparticle.cxx:242
 digparticle.cxx:243
 digparticle.cxx:244
 digparticle.cxx:245
 digparticle.cxx:246
 digparticle.cxx:247
 digparticle.cxx:248
 digparticle.cxx:249
 digparticle.cxx:250
 digparticle.cxx:251
 digparticle.cxx:252
 digparticle.cxx:253
 digparticle.cxx:254
 digparticle.cxx:255
 digparticle.cxx:256
 digparticle.cxx:257
 digparticle.cxx:258
 digparticle.cxx:259
 digparticle.cxx:260
 digparticle.cxx:261
 digparticle.cxx:262
 digparticle.cxx:263
 digparticle.cxx:264
 digparticle.cxx:265
 digparticle.cxx:266
 digparticle.cxx:267
 digparticle.cxx:268
 digparticle.cxx:269
 digparticle.cxx:270
 digparticle.cxx:271
 digparticle.cxx:272
 digparticle.cxx:273
 digparticle.cxx:274
 digparticle.cxx:275
 digparticle.cxx:276
 digparticle.cxx:277
 digparticle.cxx:278
 digparticle.cxx:279
 digparticle.cxx:280
 digparticle.cxx:281
 digparticle.cxx:282
 digparticle.cxx:283
 digparticle.cxx:284
 digparticle.cxx:285
 digparticle.cxx:286
 digparticle.cxx:287
 digparticle.cxx:288
 digparticle.cxx:289
 digparticle.cxx:290
 digparticle.cxx:291
 digparticle.cxx:292
 digparticle.cxx:293
 digparticle.cxx:294
 digparticle.cxx:295
 digparticle.cxx:296
 digparticle.cxx:297
 digparticle.cxx:298
 digparticle.cxx:299
 digparticle.cxx:300
 digparticle.cxx:301
 digparticle.cxx:302
 digparticle.cxx:303
 digparticle.cxx:304
 digparticle.cxx:305
 digparticle.cxx:306
 digparticle.cxx:307
 digparticle.cxx:308
 digparticle.cxx:309
 digparticle.cxx:310
 digparticle.cxx:311
 digparticle.cxx:312
 digparticle.cxx:313
 digparticle.cxx:314
 digparticle.cxx:315
 digparticle.cxx:316
 digparticle.cxx:317
 digparticle.cxx:318
 digparticle.cxx:319
 digparticle.cxx:320
 digparticle.cxx:321
 digparticle.cxx:322
 digparticle.cxx:323
 digparticle.cxx:324
 digparticle.cxx:325
 digparticle.cxx:326
 digparticle.cxx:327
 digparticle.cxx:328
 digparticle.cxx:329
 digparticle.cxx:330
 digparticle.cxx:331
 digparticle.cxx:332
 digparticle.cxx:333
 digparticle.cxx:334
 digparticle.cxx:335
 digparticle.cxx:336
 digparticle.cxx:337
 digparticle.cxx:338
 digparticle.cxx:339
 digparticle.cxx:340
 digparticle.cxx:341
 digparticle.cxx:342
 digparticle.cxx:343
 digparticle.cxx:344
 digparticle.cxx:345
 digparticle.cxx:346
 digparticle.cxx:347
 digparticle.cxx:348
 digparticle.cxx:349
 digparticle.cxx:350
 digparticle.cxx:351
 digparticle.cxx:352
 digparticle.cxx:353
 digparticle.cxx:354
 digparticle.cxx:355
 digparticle.cxx:356
 digparticle.cxx:357
 digparticle.cxx:358
 digparticle.cxx:359
 digparticle.cxx:360
 digparticle.cxx:361
 digparticle.cxx:362
 digparticle.cxx:363
 digparticle.cxx:364
 digparticle.cxx:365
 digparticle.cxx:366
 digparticle.cxx:367
 digparticle.cxx:368
 digparticle.cxx:369
 digparticle.cxx:370
 digparticle.cxx:371
 digparticle.cxx:372
 digparticle.cxx:373
 digparticle.cxx:374
 digparticle.cxx:375
 digparticle.cxx:376
 digparticle.cxx:377
 digparticle.cxx:378
 digparticle.cxx:379
 digparticle.cxx:380
 digparticle.cxx:381
 digparticle.cxx:382
 digparticle.cxx:383
 digparticle.cxx:384
 digparticle.cxx:385
 digparticle.cxx:386
 digparticle.cxx:387
 digparticle.cxx:388
 digparticle.cxx:389
 digparticle.cxx:390
 digparticle.cxx:391
 digparticle.cxx:392
 digparticle.cxx:393
 digparticle.cxx:394
 digparticle.cxx:395
 digparticle.cxx:396
 digparticle.cxx:397
 digparticle.cxx:398
 digparticle.cxx:399
 digparticle.cxx:400
 digparticle.cxx:401
 digparticle.cxx:402
 digparticle.cxx:403
 digparticle.cxx:404
 digparticle.cxx:405
 digparticle.cxx:406
 digparticle.cxx:407
 digparticle.cxx:408
 digparticle.cxx:409
 digparticle.cxx:410
 digparticle.cxx:411
 digparticle.cxx:412
 digparticle.cxx:413
 digparticle.cxx:414
 digparticle.cxx:415
 digparticle.cxx:416
 digparticle.cxx:417
 digparticle.cxx:418
 digparticle.cxx:419
 digparticle.cxx:420
 digparticle.cxx:421
 digparticle.cxx:422
 digparticle.cxx:423
 digparticle.cxx:424
 digparticle.cxx:425
 digparticle.cxx:426
 digparticle.cxx:427
 digparticle.cxx:428
 digparticle.cxx:429
 digparticle.cxx:430
 digparticle.cxx:431
 digparticle.cxx:432
 digparticle.cxx:433
 digparticle.cxx:434
 digparticle.cxx:435
 digparticle.cxx:436
 digparticle.cxx:437
 digparticle.cxx:438
 digparticle.cxx:439
 digparticle.cxx:440
 digparticle.cxx:441
 digparticle.cxx:442
 digparticle.cxx:443
 digparticle.cxx:444
 digparticle.cxx:445
 digparticle.cxx:446
 digparticle.cxx:447
 digparticle.cxx:448
 digparticle.cxx:449
 digparticle.cxx:450
 digparticle.cxx:451
 digparticle.cxx:452
 digparticle.cxx:453
 digparticle.cxx:454
 digparticle.cxx:455
 digparticle.cxx:456
 digparticle.cxx:457
 digparticle.cxx:458
 digparticle.cxx:459
 digparticle.cxx:460
 digparticle.cxx:461
 digparticle.cxx:462
 digparticle.cxx:463
 digparticle.cxx:464
 digparticle.cxx:465
 digparticle.cxx:466
 digparticle.cxx:467
 digparticle.cxx:468
 digparticle.cxx:469
 digparticle.cxx:470
 digparticle.cxx:471
 digparticle.cxx:472
 digparticle.cxx:473
 digparticle.cxx:474
 digparticle.cxx:475
 digparticle.cxx:476
 digparticle.cxx:477
 digparticle.cxx:478
 digparticle.cxx:479
 digparticle.cxx:480
 digparticle.cxx:481
 digparticle.cxx:482
 digparticle.cxx:483
 digparticle.cxx:484
 digparticle.cxx:485
 digparticle.cxx:486
 digparticle.cxx:487
 digparticle.cxx:488
 digparticle.cxx:489
 digparticle.cxx:490
 digparticle.cxx:491
 digparticle.cxx:492
 digparticle.cxx:493
 digparticle.cxx:494
 digparticle.cxx:495
 digparticle.cxx:496
 digparticle.cxx:497
 digparticle.cxx:498
 digparticle.cxx:499
 digparticle.cxx:500
 digparticle.cxx:501
 digparticle.cxx:502
 digparticle.cxx:503
 digparticle.cxx:504
 digparticle.cxx:505
 digparticle.cxx:506
 digparticle.cxx:507
 digparticle.cxx:508
 digparticle.cxx:509
 digparticle.cxx:510
 digparticle.cxx:511
 digparticle.cxx:512
 digparticle.cxx:513
 digparticle.cxx:514
 digparticle.cxx:515
 digparticle.cxx:516
 digparticle.cxx:517
 digparticle.cxx:518
 digparticle.cxx:519
 digparticle.cxx:520
 digparticle.cxx:521
 digparticle.cxx:522
 digparticle.cxx:523
 digparticle.cxx:524
 digparticle.cxx:525
 digparticle.cxx:526
 digparticle.cxx:527
 digparticle.cxx:528
 digparticle.cxx:529
 digparticle.cxx:530
 digparticle.cxx:531
 digparticle.cxx:532
 digparticle.cxx:533
 digparticle.cxx:534
 digparticle.cxx:535
 digparticle.cxx:536
 digparticle.cxx:537
 digparticle.cxx:538
 digparticle.cxx:539
 digparticle.cxx:540
 digparticle.cxx:541
 digparticle.cxx:542
 digparticle.cxx:543
 digparticle.cxx:544
 digparticle.cxx:545
 digparticle.cxx:546
 digparticle.cxx:547
 digparticle.cxx:548
 digparticle.cxx:549
 digparticle.cxx:550
 digparticle.cxx:551
 digparticle.cxx:552
 digparticle.cxx:553
 digparticle.cxx:554
 digparticle.cxx:555
 digparticle.cxx:556
 digparticle.cxx:557
 digparticle.cxx:558
 digparticle.cxx:559
 digparticle.cxx:560
 digparticle.cxx:561
 digparticle.cxx:562
 digparticle.cxx:563
 digparticle.cxx:564
 digparticle.cxx:565
 digparticle.cxx:566
 digparticle.cxx:567
 digparticle.cxx:568
 digparticle.cxx:569
 digparticle.cxx:570
 digparticle.cxx:571
 digparticle.cxx:572
 digparticle.cxx:573
 digparticle.cxx:574
 digparticle.cxx:575
 digparticle.cxx:576
 digparticle.cxx:577
 digparticle.cxx:578
 digparticle.cxx:579
 digparticle.cxx:580
 digparticle.cxx:581
 digparticle.cxx:582
 digparticle.cxx:583
 digparticle.cxx:584
 digparticle.cxx:585
 digparticle.cxx:586
 digparticle.cxx:587
 digparticle.cxx:588
 digparticle.cxx:589
 digparticle.cxx:590
 digparticle.cxx:591
 digparticle.cxx:592
 digparticle.cxx:593
 digparticle.cxx:594
 digparticle.cxx:595
 digparticle.cxx:596
 digparticle.cxx:597
 digparticle.cxx:598
 digparticle.cxx:599