]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCClusterParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Correct assignment of clusters to the catgeory ALL, which was in fact wrongly filled...
[u/mrichter/AliRoot.git] / TPC / AliTPCClusterParam.cxx
diff --git a/TPC/AliTPCClusterParam.cxx b/TPC/AliTPCClusterParam.cxx
index 703ffb0d05b7f5da8fa8f167ae259a78baefcf81..ab29b113138ea56e09ff2e6ae433227c4a52b1f7 100644 (file)
--- a/TPC/AliTPCClusterParam.cxx
+++ b/TPC/AliTPCClusterParam.cxx
@@ -164,6 +164,7 @@ AliTPCClusterParam::AliTPCClusterParam():
   //
   fPosYcor[0]   = 0;   fPosYcor[1]   = 0;   fPosYcor[2]   = 0; 
   fPosZcor[0]   = 0;   fPosZcor[1]   = 0;   fPosZcor[2]   = 0; 
+  fErrorRMSSys[0]=0;   fErrorRMSSys[1]=0; 
 }
 
 AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
@@ -271,24 +272,24 @@ void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t
   // Fit z - angular dependence of resolution 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"Resol:AngleM:Zm");
-  char varErr[100];
-  sprintf(varErr,"Sigma:AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
-  //
-  Int_t entries = tree->Draw(varVal,varCut);
+  TString varVal;
+  varVal="Resol:AngleM:Zm";
+  TString varErr;
+  varErr="Sigma:AngleS:Zs";
+  TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
+  //
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[10000], py[10000], pz[10000];
   Float_t ex[10000], ey[10000], ez[10000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -324,24 +325,24 @@ void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float
   // Fit z - angular dependence of resolution 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"Resol:AngleM:Zm");
-  char varErr[100];
-  sprintf(varErr,"Sigma:AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
-  //
-  Int_t entries = tree->Draw(varVal,varCut);
+ TString varVal;
+  varVal="Resol:AngleM:Zm";
+ TString varErr;
+  varErr="Sigma:AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
+  //
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[10000], py[10000], pz[10000];
   Float_t ex[10000], ey[10000], ez[10000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -389,24 +390,24 @@ void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Flo
   // Fit z - angular dependence of resolution - pad length scaling 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length");
-  char varErr[100];
-  sprintf(varErr,"Sigma:AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&QMean<0",dim);
-  //
-  Int_t entries = tree->Draw(varVal,varCut);
+ TString varVal;
+  varVal="Resol:AngleM*sqrt(Length):Zm/Length";
+ TString varErr;
+  varErr="Sigma:AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&QMean<0",dim);
+  //
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[10000], py[10000], pz[10000];
   Float_t ex[10000], ey[10000], ez[10000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -441,27 +442,27 @@ void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t
   // Fit z - angular dependence of resolution - Q scaling 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
+ TString varVal;
+  varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
   char varVal0[100];
-  sprintf(varVal0,"Resol:AngleM:Zm");
+  snprintf(varVal0,100,"Resol:AngleM:Zm");
   //
-  char varErr[100];
-  sprintf(varErr,"Sigma:AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
+ TString varErr;
+  varErr="Sigma:AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
   //
-  Int_t entries = tree->Draw(varVal,varCut);
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
   Float_t ex[20000], ey[20000], ez[20000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -508,27 +509,27 @@ void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float
   // Fit z - angular dependence of resolution - Q scaling  - parabolic correction
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean");
+ TString varVal;
+  varVal="Resol:AngleM/sqrt(QMean):Zm/QMean";
   char varVal0[100];
-  sprintf(varVal0,"Resol:AngleM:Zm");
+  snprintf(varVal0,100,"Resol:AngleM:Zm");
   //
-  char varErr[100];
-  sprintf(varErr,"Sigma:AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
+ TString varErr;
+  varErr="Sigma:AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
   //
-  Int_t entries = tree->Draw(varVal,varCut);
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
   Float_t ex[20000], ey[20000], ez[20000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -588,24 +589,24 @@ void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *p
   // Fit z - angular dependence of resolution 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"RMSm:AngleM:Zm");
-  char varErr[100];
-  sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
-  //
-  Int_t entries = tree->Draw(varVal,varCut);
+ TString varVal;
+  varVal="RMSm:AngleM:Zm";
+ TString varErr;
+  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
+  //
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[10000], py[10000], pz[10000];
   Float_t ex[10000], ey[10000], ez[10000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -640,24 +641,24 @@ void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float
   // Fit z - angular dependence of resolution - pad length scaling 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"RMSm:AngleM*Length:Zm");
-  char varErr[100];
-  sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&QMean<0",dim);
-  //
-  Int_t entries = tree->Draw(varVal,varCut);
+ TString varVal;
+  varVal="RMSm:AngleM*Length:Zm";
+ TString varErr;
+  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad";
+ TString varCut;
+  varCut=Form("Dim==%d&&QMean<0",dim);
+  //
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[10000], py[10000], pz[10000];
   Float_t type[10000], ey[10000], ez[10000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     type[ipoint] = tree->GetV3()[ipoint];
     ey[ipoint]   = tree->GetV2()[ipoint];
     ez[ipoint]   = tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -695,27 +696,27 @@ void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *p
   // Fit z - angular dependence of resolution - Q scaling 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean");
+ TString varVal;
+  varVal="RMSm:AngleM/sqrt(QMean):Zm/QMean";
   char varVal0[100];
-  sprintf(varVal0,"RMSm:AngleM:Zm");
+  snprintf(varVal0,100,"RMSm:AngleM:Zm");
   //
-  char varErr[100];
-  sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs");
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type);
+ TString varErr;
+  varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs";
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type);
   //
-  Int_t entries = tree->Draw(varVal,varCut);
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000];
   Float_t ex[20000], ey[20000], ez[20000];
   //
-  tree->Draw(varErr,varCut);  
+  tree->Draw(varErr.Data(),varCut);  
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     ex[ipoint]= tree->GetV3()[ipoint];
     ey[ipoint]= tree->GetV2()[ipoint];
     ez[ipoint]= tree->GetV1()[ipoint];
   } 
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV3()[ipoint];
     py[ipoint]= tree->GetV2()[ipoint];
@@ -763,16 +764,16 @@ void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_
   // Fit z - angular dependence of resolution - Q scaling 
   //
   // Int_t dim=0, type=0;
-  char varVal[100];
-  sprintf(varVal,"RMSs:RMSm");
+  TString varVal;
+  varVal="RMSs:RMSm";
   //
-  char varCut[100];
-  sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type);
+ TString varCut;
+  varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type);
   //
-  Int_t entries = tree->Draw(varVal,varCut);
+  Int_t entries = tree->Draw(varVal.Data(),varCut);
   Float_t px[20000], py[20000];
   //
-  tree->Draw(varVal,varCut);
+  tree->Draw(varVal.Data(),varCut);
   for (Int_t ipoint=0; ipoint<entries; ipoint++){
     px[ipoint]= tree->GetV2()[ipoint];
     py[ipoint]= tree->GetV1()[ipoint];
@@ -1155,19 +1156,19 @@ void AliTPCClusterParam::Test(TTree * tree, const char *output){
       char hcut1[300];
       char hexp1[300];
       //
-      sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad);
-      sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
-      sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
+      snprintf(hname1,300,"Delta0 Dir %d Pad %d",idim,ipad);
+      snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
+      snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
       TH1F  his1DRel0(hname1, hname1, 100,-0.2, 0.2);
-      sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
+      snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
       tree->Draw(hexp1,hcut1,"");
       his1DRel0.Write();
       //
-      sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad);
-      sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
-      sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
+      snprintf(hname1,300,"Delta0Par Dir %d Pad %d",idim,ipad);
+      snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
+      snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1);
       TH1F  his1DRel0Par(hname1, hname1, 100,-0.2, 0.2);
-      sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
+      snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
       tree->Draw(hexp1,hcut1,"");
       his1DRel0Par.Write();
       //
@@ -1182,19 +1183,19 @@ void AliTPCClusterParam::Test(TTree * tree, const char *output){
       char hcut1[300];
       char hexp1[300];
       //
-      sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad);
-      sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
-      sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
+      snprintf(hname1,300,"2DDelta0 Dir %d Pad %d",idim,ipad);
+      snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
+      snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
       TProfile2D  profDRel0(hname1, hname1, 6,0,250,6,0,1);
-      sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
+      snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
       tree->Draw(hexp1,hcut1,"");
       profDRel0.Write();
       //
-      sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad);
-      sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
-      sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
+      snprintf(hname1,300,"2DDelta0Par Dir %d Pad %d",idim,ipad);
+      snprintf(hcut1,300,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad);
+      snprintf(hexp1,300,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1);
       TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1);
-      sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
+      snprintf(hname1,300,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad);
       tree->Draw(hexp1,hcut1,"");
       profDRel0Par.Write();
       //
@@ -1597,20 +1598,26 @@ Double_t  AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_
   //TF1 f1("f1","AliTPCClusterParam::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
   //TF2 f2("f2","AliTPCClusterParam::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
   //
-  const Float_t kEpsilon = 0.0001;
+  const Double_t kEpsilon = 0.0001;
+  const Double_t twoPi = TMath::TwoPi();
+  const Double_t hnorm = 0.5/TMath::Sqrt(twoPi);
+  const Double_t sqtwo = TMath::Sqrt(2.);
+
   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
     // small angular effect
-    Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
+    Double_t val = TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1)/(s0*s1*twoPi);
     return val;
   }
   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
-  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
-  //
-  Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);                       
-  Double_t erf0 = AliMathBase::ErfFast( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
-  Double_t erf1 = AliMathBase::ErfFast( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
-  Double_t norm = 1./TMath::Sqrt(sigma2);
-  norm/=2.*TMath::Sqrt(2.*TMath::Pi());
+  Double_t sigma = TMath::Sqrt(sigma2);
+  Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2.*sigma2));
+  //
+  Double_t sigmaErf =  1./(2.*s0*s1*sqtwo*sigma);
+  Double_t k0s1s1 = 2.*k0*s1*s1;
+  Double_t k1s0s0 = 2.*k1*s0*s0;
+  Double_t erf0 = AliMathBase::ErfFast((sigma2-k0s1s1*x0-k1s0s0*x1)*sigmaErf);
+  Double_t erf1 = AliMathBase::ErfFast((sigma2+k0s1s1*x0+k1s0s0*x1)*sigmaErf);
+  Double_t norm = hnorm/sigma;
   Double_t val  = norm*exp0*(erf0+erf1);
   return val;
 }