#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){
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){
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];
}
}
Double_t SumGausLorentz2D(Double_t *x, Double_t *par){
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()
{
}
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();
}
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);
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() {
}
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;
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;
if(ChargeModel==1){
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){
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){
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);
}
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 ;
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;
}