TPC/fastSimul/AliTPCclusterFast.cxx

   1 /*
   2   gSystem->Load("libSTAT.so");
   3   .x ~/NimStyle.C
   4
   5  .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
   6  //
   7  AliTPCclusterFast::fPRF = new TF1("fprf","gausn",-5,5);
   8  AliTPCclusterFast::fTRF = new TF1("ftrf","gausn",-5,5);
   9  AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
  10  AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
  11
  12  //
  13  AliTPCtrackFast::Simul("trackerSimul.root",100);
  14 // AliTPCclusterFast::Simul("cluterSimul.root",20000);
  15 */
  16
  17 #include "TObject.h"
  18 #include "TF1.h"
  19 #include "TMath.h"
  20 #include "TRandom.h"
  21 #include "TVectorD.h"
  22 #include "TMatrixD.h"
  23 #include "TH1.h"
  24 #include "TClonesArray.h"
  25 #include "TTreeStream.h"
  26
  27 class AliTPCclusterFast: public TObject {
  28 public:
  29   AliTPCclusterFast();
  30   virtual ~AliTPCclusterFast();
  31   void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
  32   void GenerElectrons();
  33   void Digitize();
  34   Double_t GetQtot(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
  35   Double_t GetQmax(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
  36   Double_t GetQmaxCorr(Float_t rmsy0, Float_t rmsz0);
  37   Double_t GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr);
  38
  39   Double_t GetNsec();
  40   static void Simul(const char* simul, Int_t npoints);
  41   static Double_t GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1);
  42   static Double_t GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1);
  43   static Double_t GaussGamma4(Double_t x, Double_t s0, Double_t p1);
  44   static Double_t Gamma4(Double_t x, Double_t p0, Double_t p1);
  45 public:
  46   Float_t fMNprim;     // mean number of primary electrons
  47   //                   //electrons part input
  48   Int_t   fNprim;      // mean number of primary electrons
  49   Int_t   fNtot;       // total number of  electrons
  50   Float_t fQtot;       // total charge - Gas gain flucuation taken into account
  51   //
  52   Float_t fDiff;       // diffusion sigma
  53   Float_t fY;          // y position
  54   Float_t fZ;          // z postion
  55   Float_t fAngleY;     // y angle - tan(y)
  56   Float_t fAngleZ;     // z angle - tan z
  57   //
  58   //
  59   //                   // electron part simul
  60   TVectorD fSec;       //! number of secondary electrons
  61   TVectorD fPosY;      //! position y for each electron
  62   TVectorD fPosZ;      //! position z for each electron
  63   TVectorD fGain;      //! gg for each electron
  64   //
  65   TVectorD fStatY;     //!stat Y
  66   TVectorD fStatZ;     //!stat Y
  67   //
  68   // digitization part
  69   //
  70   TMatrixD fDigits;    // response matrix
  71   static TF1* fPRF;    // Pad response
  72   static TF1* fTRF;    // Time response function
  73   ClassDef(AliTPCclusterFast,1)  // container for
  74 };
  75
  76
  77 class AliTPCtrackFast: public TObject {
  78 public:
  79   AliTPCtrackFast();
  80   void Add(AliTPCtrackFast &track2);
  81   void MakeTrack();
  82   static void Simul(const char* simul, Int_t ntracks);
  83   Double_t  CookdEdxNtot(Double_t f0,Float_t f1);
  84   Double_t  CookdEdxQtot(Double_t f0,Float_t f1);
  85   //
  86   Double_t  CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t corr = kTRUE);
  87   Double_t  CookdEdxDmax(Double_t f0,Float_t f1,Float_t gain,Float_t thr, Float_t noise, Bool_t corr=kTRUE);
  88   //
  89   Double_t  CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1);
  90   //
  91   Float_t fMNprim;     // mean number of primary electrons
  92   Float_t fAngleY;     // y angle - tan(y)
  93   Float_t fAngleZ;     // z angle - tan z
  94   Float_t fDiff;       // diffusion
  95   Int_t   fN;          // number of clusters
  96   TClonesArray *fCl;   // array of clusters
  97   //
  98   Bool_t   fInit;      // initialization flag
  99   //
 100   //
 101   ClassDef(AliTPCtrackFast,2)  // container for
 102 };
 103
 104
 105
 106 ClassImp(AliTPCclusterFast)
 107 ClassImp(AliTPCtrackFast)
 108
 109
 110
 111
 112
 113 TF1 *AliTPCclusterFast::fPRF=0;
 114 TF1 *AliTPCclusterFast::fTRF=0;
 115
 116
 117 AliTPCtrackFast::AliTPCtrackFast():
 118   TObject(),
 119   fMNprim(0),
 120   fAngleY(0),
 121   fAngleZ(0),
 122   fN(0),
 123   fCl(0),
 124   fInit(kFALSE)
 125 {
 126   //
 127   //
 128   //
 129 }
 130
 131 void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
 132   if (!track2.fInit) return;
 133
 134 }
 135
 136
 137
 138 void AliTPCtrackFast::MakeTrack(){
 139   //
 140   //
 141   //
 142   if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
 143   for (Int_t i=0;i<fN;i++){
 144     Double_t tY = i*fAngleY;
 145     Double_t tZ = i*fAngleZ;
 146     AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
 147     if (!cluster) cluster =   new ((*fCl)[i]) AliTPCclusterFast;
 148     //
 149     Double_t posY = tY-TMath::Nint(tY);
 150     Double_t posZ = tZ-TMath::Nint(tZ);
 151     cluster->SetParam(fMNprim,fDiff,posY,posZ,fAngleY,fAngleZ);
 152     //
 153     cluster->GenerElectrons();
 154     cluster->Digitize();
 155   }
 156 }
 157
 158 Double_t  AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
 159   //
 160   Double_t amp[160];
 161   for (Int_t i=0;i<fN;i++){
 162     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 163     amp[i]=cluster->fNtot;
 164   }
 165   return CookdEdx(fN,amp,f0,f1);
 166 }
 167
 168 Double_t  AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
 169   //
 170   Double_t amp[160];
 171   for (Int_t i=0;i<fN;i++){
 172     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 173     amp[i]=cluster->fQtot;
 174   }
 175   return CookdEdx(fN,amp,f0,f1);
 176 }
 177
 178 Double_t   AliTPCtrackFast::CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
 179   //
 180   //
 181   //
 182   Double_t amp[160];
 183   Int_t over=0;
 184   for (Int_t i=0;i<fN;i++){
 185     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 186     Float_t camp = cluster->GetQtot(gain,thr,noise);
 187     if (camp==0) continue;
 188     Float_t corr =  1;
 189     if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
 190     amp[over]=camp/corr;
 191     over++;
 192   }
 193   return CookdEdx(over,amp,f0,f1);
 194
 195 }
 196
 197 Double_t   AliTPCtrackFast::CookdEdxDmax(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
 198   //
 199   //
 200   //
 201   Double_t amp[160];
 202   Int_t over=0;
 203   for (Int_t i=0;i<fN;i++){
 204     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 205     Float_t camp = cluster->GetQmax(gain,thr,noise);
 206     if (camp==0) continue;
 207     Float_t corr =  1;
 208     if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
 209     amp[over]=camp/corr;
 210     over++;
 211   }
 212   return CookdEdx(over,amp,f0,f1);
 213
 214 }
 215
 216
 217 Double_t  AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1){
 218   //
 219   //
 220   //
 221   Int_t index[160];
 222   TMath::Sort(npoints,amp,index,kFALSE);
 223   Float_t sum0=0, sum1=0,sum2=0;
 224   for (Int_t i=0;i<npoints;i++){
 225     if (i<npoints*f0) continue;
 226     if (i>npoints*f1) continue;
 227     sum0++;
 228     sum1+= amp[index[i]];
 229     sum2+= amp[index[i]];
 230   }
 231   if (sum0<=0) return 0;
 232   return sum1/sum0;
 233 }
 234
 235 void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks){
 236   //
 237   //
 238   //
 239   AliTPCtrackFast fast;
 240   TTreeSRedirector cstream(fname,"recreate");
 241   for (Int_t itr=0; itr<ntracks; itr++){
 242     //
 243     fast.fMNprim=(5+50*gRandom->Rndm());
 244     fast.fDiff =0.01 +0.35*gRandom->Rndm();
 245     //
 246     fast.fAngleY   = 4.0*(gRandom->Rndm()-0.5);
 247     fast.fAngleZ   = 4.0*(gRandom->Rndm()-0.5);
 248     fast.fN  = TMath::Nint(80.+gRandom->Rndm()*80.);
 249     fast.MakeTrack();
 250     if (itr%100==0) printf("%d\n",itr);
 251     cstream<<"simulTrack"<<
 252       "tr.="<<&fast<<
 253       "\n";
 254   }
 255   fast.Write("track");
 256 }
 257
 258
 259
 260 AliTPCclusterFast::AliTPCclusterFast(){
 261   //
 262   //
 263   fDigits.ResizeTo(5,7);
 264 }
 265
 266 AliTPCclusterFast::~AliTPCclusterFast(){
 267 }
 268
 269
 270 void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
 271   //
 272   //
 273   fMNprim = mnprim; fDiff = diff;
 274   fY=y; fZ=z;
 275   fAngleY=ky; fAngleZ=kz;
 276 }
 277 Double_t AliTPCclusterFast::GetNsec(){
 278   //
 279   // Generate number of secondary electrons
 280   // copy of procedure implemented in geant
 281   //
 282   const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
 283   const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
 284   const Double_t W=20.77E-9;
 285   Float_t RAN = gRandom->Rndm();
 286   //Double_t edep = TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO);
 287   //edep = TMath::Min(edep, EEND);
 288   //return TMath::Nint(edep/W);
 289   return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
 290 }
 291
 292 void AliTPCclusterFast::GenerElectrons(){
 293   //
 294   //
 295   //
 296   //
 297   const Int_t knMax=1000;
 298   if (fPosY.GetNrows()<knMax){
 299     fPosY.ResizeTo(knMax);
 300     fPosZ.ResizeTo(knMax);
 301     fGain.ResizeTo(knMax);
 302     fSec.ResizeTo(knMax);
 303     fStatY.ResizeTo(3);
 304     fStatZ.ResizeTo(3);
 305   }
 306   fNprim = gRandom->Poisson(fMNprim);  //number of primary electrons
 307   fNtot=0; //total number of electrons
 308   fQtot=0; //total number of electrons after gain multiplification
 309   //
 310   Double_t sumQ=0;
 311   Double_t sumYQ=0;
 312   Double_t sumZQ=0;
 313   Double_t sumY2Q=0;
 314   Double_t sumZ2Q=0;
 315   for (Int_t i=0;i<knMax;i++){
 316     fSec[i]=0;
 317   }
 318   for (Int_t iprim=0; iprim<fNprim;iprim++){
 319     Float_t dN   =  GetNsec();
 320     fSec[iprim]=dN;
 321     Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
 322     Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
 323     for (Int_t isec=0;isec<=dN;isec++){
 324       //
 325       //
 326       Double_t y = gRandom->Gaus(0,fDiff)+yc;
 327       Double_t z = gRandom->Gaus(0,fDiff)+zc;
 328       Double_t gg = -TMath::Log(gRandom->Rndm());
 329       fPosY[fNtot]=y;
 330       fPosZ[fNtot]=z;
 331       fGain[fNtot]=gg;
 332       fQtot+=gg;
 333       fNtot++;
 334       sumQ+=gg;
 335       sumYQ+=gg*y;
 336       sumY2Q+=gg*y*y;
 337       sumZQ+=gg*z;
 338       sumZ2Q+=gg*z*z;
 339       if (fNtot>=knMax) break;
 340     }
 341     if (fNtot>=knMax) break;
 342   }
 343   if (sumQ>0){
 344     fStatY[0]=sumQ;
 345     fStatY[1]=sumYQ/sumQ;
 346     fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
 347     fStatZ[0]=sumQ;
 348     fStatZ[1]=sumZQ/sumQ;
 349     fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
 350   }
 351 }
 352
 353 void AliTPCclusterFast::Digitize(){
 354   //
 355   //
 356   //
 357   // 1. Clear digits
 358   for (Int_t i=0; i<5;i++)
 359     for (Int_t j=0; j<7;j++){
 360       fDigits(i,j)=0;
 361     }
 362   //
 363   // Fill digits
 364   for (Int_t iel = 0; iel<fNtot; iel++){
 365     for (Int_t di=-2; di<=2;di++)
 366       for (Int_t dj=-3; dj<=3;dj++){
 367         Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
 368         fac*=fGain[iel];
 369         fDigits(2+di,3+dj)+=fac;
 370       }
 371   }
 372
 373 }
 374
 375
 376
 377 void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
 378   //
 379   // Calc rms
 380   //
 381   AliTPCclusterFast fast;
 382   TTreeSRedirector cstream(fname);
 383   for (Int_t icl=0; icl<npoints; icl++){
 384     Float_t nprim=(10+20*gRandom->Rndm());
 385     Float_t diff =0.01 +0.35*gRandom->Rndm();
 386     Float_t posY = gRandom->Rndm()-0.5;
 387     Float_t posZ = gRandom->Rndm()-0.5;
 388     //
 389     Float_t ky   = 4.0*(gRandom->Rndm()-0.5);
 390     Float_t kz   = 4.0*(gRandom->Rndm()-0.5);
 391     fast.SetParam(nprim,diff,posY,posZ,ky,kz);
 392     fast.GenerElectrons();
 393     fast.Digitize();
 394     if (icl%10000==0) printf("%d\n",icl);
 395     cstream<<"simul"<<
 396       "s.="<<&fast<<
 397       "\n";
 398   }
 399 }
 400
 401
 402 Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 403   //
 404   //
 405   //
 406   Float_t sum =0;
 407   for (Int_t ip=0;ip<5;ip++){
 408     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 409     for (Int_t it=0;it<7;it++){
 410       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 411       if (baddPedestal) amp+=pedestal;
 412       if (brounding) amp=TMath::Nint(amp);
 413       if (amp>thr) sum+=amp;
 414     }
 415   }
 416   return sum;
 417 }
 418
 419 Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 420   //
 421   //
 422   //
 423   Float_t max =0;
 424   for (Int_t ip=0;ip<5;ip++){
 425     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 426     for (Int_t it=0;it<7;it++){
 427       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 428       if (baddPedestal) amp+=pedestal;
 429       if (brounding) amp=TMath::Nint(amp);
 430       if (amp>max && amp>thr) max=amp;
 431     }
 432   }
 433   return max;
 434 }
 435
 436
 437
 438 Double_t  AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
 439   //
 440   // Gaus distribution convolueted with rectangular
 441   // Gaus width sy and sz is determined by RF width and diffusion
 442   // Integral of Q is equal 1
 443   // Q max is calculated at position fY,fX
 444   //
 445   //
 446   //
 447   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 448   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 449   return GaussConvolution(fY,fZ, fAngleY,fAngleZ,sy,sz);
 450 }
 451
 452
 453 Double_t  AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
 454   //
 455   //  Calculates the fraction of the charge over threshol to total charge
 456   //  The response function
 457   //
 458   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 459   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 460   Double_t sumAll=0,sumThr=0;
 461   Double_t qtot = GetQtot(gain,thr,0); // sum of signal over threshold
 462   //
 463   Double_t corr =1;
 464   Double_t qnorm=qtot;
 465   for (Int_t iter=0;iter<2;iter++){
 466     for (Int_t iy=-2;iy<=2;iy++)
 467       for (Int_t iz=-2;iz<=2;iz++){
 468         Double_t val = GaussConvolution(fY-iy,fZ-iz, fAngleY,fAngleZ,sy,sz);
 469         Double_t qlocal =TMath::Nint(qnorm*val);
 470         if (qlocal>thr) sumThr+=qlocal;
 471         sumAll+=qlocal;
 472       }
 473     if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
 474     //corr = sumThr;
 475     if (corr>0) qnorm=qtot/corr;
 476
 477   }
 478   return corr;
 479 }
 480
 481
 482
 483
 484
 485 Double_t  AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
 486   //
 487   // 2 D gaus convoluted with angular effect
 488   // See in mathematica:
 489   //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
 490   //
 491   //TF1 f1("f1","AliTPCclusterFast::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
 492   //TF2 f2("f2","AliTPCclusterFast::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
 493   //
 494   const Float_t kEpsilon = 0.0001;
 495   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
 496     // small angular effect
 497     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
 498     return val;
 499   }
 500
 501   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
 502   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
 503   //
 504   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);
 505   Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
 506   Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
 507   Double_t norm = 1./TMath::Sqrt(sigma2);
 508   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
 509   Double_t val  = norm*exp0*(erf0+erf1);
 510   return val;
 511
 512 }
 513
 514
 515 Double_t  AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
 516  //
 517   // 2 D gaus convoluted with exponential
 518   // Integral nomalized to 1
 519   // See in mathematica:
 520   //Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-x1*t1],{x1,0,Infinity}]]
 521   // TF1 fgexp("fgexp","AliTPCclusterFast::GaussExpConvolution(x,0.5,1)",-2,2)
 522   Double_t exp1 = (s0*s0*t1-2*x0)*t1/2.;
 523   exp1 = TMath::Exp(exp1);
 524   Double_t erf = 1+TMath::Erf((-s0*s0*t1+x0)/(s0*TMath::Sqrt(2.)));
 525   Double_t val = exp1*erf;
 526   val *=t1/(2.);
 527   return val;
 528
 529 }
 530
 531
 532 Double_t  AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
 533   //
 534   // Gamma 4 Time response function of ALTRO
 535   //
 536   if (x<0) return 0;
 537   Double_t g1 = TMath::Exp(-4.*x/p1);
 538   Double_t g2 = TMath::Power(x/p1,4);
 539   return p0*g1*g2;
 540 }
 541
 542
 543
 544 Double_t  AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
 545   //
 546   // Gamma 4 Time response function of ALTRO convoluted with Gauss
 547   // Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-4*x1/p1]*(x/p1)^4/s0,{x1,0,Infinity}]]
 548   //TF1 fgg4("fgg4","AliTPCclusterFast::GaussGamma4(x,0.5,0.5)",-2,2)
 549
 550   Double_t exp1 = (8*s0*s0-4.*p1*x)/(p1*p1);
 551   exp1 = TMath::Exp(exp1);
 552   Double_t erf1 = 1+TMath::Erf((-4*s0/p1+x/s0)/TMath::Sqrt(2));
 553   //  Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
 554   return exp1*erf1;
 555
 556
 557 }
 558
 559 // Analytical sollution only in 1D - too long expression
 560 // Simplify[Integrate[Exp[-(x0-(x1-k*x2))*(x0-(x1-k*x2))/(2*s0*s0)]*Exp[-(x1*t1-k*x2)],{x2,-1,1}]]
 561 //
 562 //
 563 // No analytical solution
 564 //
 565 //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-xt-k1*xd)*(x1-xt-k1*xd)/(2*s1*s1)]*Exp[-kt*xt]/(s0*s1),{xd,-1/2,1/2},{xt,0,Infinity}]]