Last minute changes (C.Blume)
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.10  2001/05/07 08:06:44  cblume
19 Speedup of the code. Create only AliTRDcluster
20
21 Revision 1.9  2000/11/01 14:53:20  cblume
22 Merge with TRD-develop
23
24 Revision 1.1.4.5  2000/10/15 23:40:01  cblume
25 Remove AliTRDconst
26
27 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
28 Made Getters const
29
30 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
31 Replace include files by forward declarations
32
33 Revision 1.1.4.2  2000/09/22 14:49:49  cblume
34 Adapted to tracking code
35
36 Revision 1.8  2000/10/02 21:28:19  fca
37 Removal of useless dependecies via forward declarations
38
39 Revision 1.7  2000/06/27 13:08:50  cblume
40 Changed to Copy(TObject &A) to appease the HP-compiler
41
42 Revision 1.6  2000/06/09 11:10:07  cblume
43 Compiler warnings and coding conventions, next round
44
45 Revision 1.5  2000/06/08 18:32:58  cblume
46 Make code compliant to coding conventions
47
48 Revision 1.4  2000/06/07 16:27:01  cblume
49 Try to remove compiler warnings on Sun and HP
50
51 Revision 1.3  2000/05/08 16:17:27  cblume
52 Merge TRD-develop
53
54 Revision 1.1.4.1  2000/05/08 15:09:01  cblume
55 Introduce AliTRDdigitsManager
56
57 Revision 1.1  2000/02/28 18:58:54  cblume
58 Add new TRD classes
59
60 */
61
62 ///////////////////////////////////////////////////////////////////////////////
63 //                                                                           //
64 // TRD cluster finder for the slow simulator. 
65 //                                                                           //
66 ///////////////////////////////////////////////////////////////////////////////
67
68 #include <TF1.h>
69 #include <TTree.h>
70 #include <TH1.h>
71
72 #include "AliRun.h"
73
74 #include "AliTRD.h"
75 #include "AliTRDclusterizerV1.h"
76 #include "AliTRDmatrix.h"
77 #include "AliTRDgeometry.h"
78 #include "AliTRDdigitizer.h"
79 #include "AliTRDdataArrayF.h"
80 #include "AliTRDdataArrayI.h"
81 #include "AliTRDdigitsManager.h"
82
83 ClassImp(AliTRDclusterizerV1)
84
85 //_____________________________________________________________________________
86 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
87 {
88   //
89   // AliTRDclusterizerV1 default constructor
90   //
91
92   fDigitsManager = NULL;
93
94   fClusMaxThresh = 0;
95   fClusSigThresh = 0;
96
97   fUseLUT        = kFALSE;
98
99 }
100
101 //_____________________________________________________________________________
102 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
103                     :AliTRDclusterizer(name,title)
104 {
105   //
106   // AliTRDclusterizerV1 default constructor
107   //
108
109   fDigitsManager = new AliTRDdigitsManager();
110
111   Init();
112
113 }
114
115 //_____________________________________________________________________________
116 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
117 {
118   //
119   // AliTRDclusterizerV1 copy constructor
120   //
121
122   ((AliTRDclusterizerV1 &) c).Copy(*this);
123
124 }
125
126 //_____________________________________________________________________________
127 AliTRDclusterizerV1::~AliTRDclusterizerV1()
128 {
129   //
130   // AliTRDclusterizerV1 destructor
131   //
132
133   if (fDigitsManager) {
134     delete fDigitsManager;
135   }
136
137 }
138
139 //_____________________________________________________________________________
140 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
141 {
142   //
143   // Assignment operator
144   //
145
146   if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
147   return *this;
148
149 }
150
151 //_____________________________________________________________________________
152 void AliTRDclusterizerV1::Copy(TObject &c)
153 {
154   //
155   // Copy function
156   //
157
158   ((AliTRDclusterizerV1 &) c).fUseLUT        = fUseLUT;
159   ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
160   ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
161   ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
162   for (Int_t ilut = 0; ilut < kNlut; ilut++) {
163     ((AliTRDclusterizerV1 &) c).fLUT[ilut] = fLUT[ilut];
164   }
165
166   AliTRDclusterizer::Copy(c);
167
168 }
169
170 //_____________________________________________________________________________
171 void AliTRDclusterizerV1::Init()
172 {
173   //
174   // Initializes the cluster finder
175   //
176
177   // The default parameter for the clustering
178   fClusMaxThresh = 3;
179   fClusSigThresh = 1;
180
181   // Use the lookup table for the position determination
182   fUseLUT        = kTRUE;
183
184   // The lookup table from Bogdan
185   Float_t lut[128] = {  
186     0.0068,  0.0198,  0.0318,  0.0432,  0.0538,  0.0642,  0.0742,  0.0838,
187     0.0932,  0.1023,  0.1107,  0.1187,  0.1268,  0.1347,  0.1423,  0.1493,  
188     0.1562,  0.1632,  0.1698,  0.1762,  0.1828,  0.1887,  0.1947,  0.2002,  
189     0.2062,  0.2118,  0.2173,  0.2222,  0.2278,  0.2327,  0.2377,  0.2428,  
190     0.2473,  0.2522,  0.2567,  0.2612,  0.2657,  0.2697,  0.2743,  0.2783,  
191     0.2822,  0.2862,  0.2903,  0.2943,  0.2982,  0.3018,  0.3058,  0.3092,  
192     0.3128,  0.3167,  0.3203,  0.3237,  0.3268,  0.3302,  0.3338,  0.3368,  
193     0.3402,  0.3433,  0.3462,  0.3492,  0.3528,  0.3557,  0.3587,  0.3613,  
194     0.3643,  0.3672,  0.3702,  0.3728,  0.3758,  0.3783,  0.3812,  0.3837,  
195     0.3862,  0.3887,  0.3918,  0.3943,  0.3968,  0.3993,  0.4017,  0.4042,  
196     0.4067,  0.4087,  0.4112,  0.4137,  0.4157,  0.4182,  0.4207,  0.4227,  
197     0.4252,  0.4272,  0.4293,  0.4317,  0.4338,  0.4358,  0.4383,  0.4403,  
198     0.4423,  0.4442,  0.4462,  0.4482,  0.4502,  0.4523,  0.4543,  0.4563,  
199     0.4582,  0.4602,  0.4622,  0.4638,  0.4658,  0.4678,  0.4697,  0.4712,  
200     0.4733,  0.4753,  0.4767,  0.4787,  0.4803,  0.4823,  0.4837,  0.4857,  
201     0.4873,  0.4888,  0.4908,  0.4922,  0.4942,  0.4958,  0.4972,  0.4988  
202   }; 
203   for (Int_t ilut = 0; ilut < kNlut; ilut++) {
204     fLUT[ilut] = lut[ilut];
205   }
206
207 }
208
209 //_____________________________________________________________________________
210 Bool_t AliTRDclusterizerV1::ReadDigits()
211 {
212   //
213   // Reads the digits arrays from the input aliroot file
214   //
215
216   if (!fInputFile) {
217     printf("AliTRDclusterizerV1::ReadDigits -- ");
218     printf("No input file open\n");
219     return kFALSE;
220   }
221
222   // Read in the digit arrays
223   return (fDigitsManager->ReadDigits());  
224
225 }
226
227 //_____________________________________________________________________________
228 Bool_t AliTRDclusterizerV1::MakeClusters()
229 {
230   //
231   // Generates the cluster.
232   //
233
234   Int_t row, col, time;
235
236   if (fTRD->IsVersion() != 1) {
237     printf("AliTRDclusterizerV1::MakeCluster -- ");
238     printf("TRD must be version 1 (slow simulator).\n");
239     return kFALSE; 
240   }
241
242   // Get the geometry
243   AliTRDgeometry *geo = fTRD->GetGeometry();
244
245   printf("AliTRDclusterizerV1::MakeCluster -- ");
246   printf("Start creating clusters.\n");
247
248   AliTRDdataArrayI *digits;
249   AliTRDdataArrayI *track0;
250   AliTRDdataArrayI *track1;
251   AliTRDdataArrayI *track2; 
252
253   // Threshold value for the maximum
254   Int_t maxThresh = fClusMaxThresh;   
255   // Threshold value for the digit signal
256   Int_t sigThresh = fClusSigThresh;   
257
258   // Iteration limit for unfolding procedure
259   const Float_t kEpsilon = 0.01;             
260
261   const Int_t   kNclus   = 3;  
262   const Int_t   kNsig    = 5;
263   const Int_t   kNtrack  = 3 * kNclus;
264
265   // For the LUT
266   const Float_t kLUTmin  = 0.106113;
267   const Float_t kLUTmax  = 0.995415;
268
269   Int_t   iType          = 0;
270   Int_t   iUnfold        = 0;
271
272   Float_t ratioLeft      = 1.0;
273   Float_t ratioRight     = 1.0;
274
275   Float_t padSignal[kNsig];   
276   Float_t clusterSignal[kNclus];
277   Float_t clusterPads[kNclus];   
278   Int_t   clusterDigit[kNclus];
279   Int_t   clusterTracks[kNtrack];   
280
281   Int_t chamBeg = 0;
282   Int_t chamEnd = AliTRDgeometry::Ncham();
283   if (fTRD->GetSensChamber()  >= 0) {
284     chamBeg = fTRD->GetSensChamber();
285     chamEnd = chamBeg + 1;
286   }
287   Int_t planBeg = 0;
288   Int_t planEnd = AliTRDgeometry::Nplan();
289   if (fTRD->GetSensPlane()    >= 0) {
290     planBeg = fTRD->GetSensPlane();
291     planEnd = planBeg + 1;
292   }
293   Int_t sectBeg = 0;
294   Int_t sectEnd = AliTRDgeometry::Nsect();
295
296   // Start clustering in every chamber
297   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
298     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
299       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
300
301         if (fTRD->GetSensSector() >= 0) {
302           Int_t sens1 = fTRD->GetSensSector();
303           Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
304           sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
305                  * AliTRDgeometry::Nsect();
306           if (sens1 < sens2) {
307             if ((isect < sens1) || (isect >= sens2)) continue;
308           }
309           else {
310             if ((isect < sens1) && (isect >= sens2)) continue;
311           }
312         }
313
314         Int_t idet = geo->GetDetector(iplan,icham,isect);
315
316         Int_t nClusters      = 0;
317         Int_t nClusters2pad  = 0;
318         Int_t nClusters3pad  = 0;
319         Int_t nClusters4pad  = 0;
320         Int_t nClusters5pad  = 0;
321         Int_t nClustersLarge = 0;
322
323         printf("AliTRDclusterizerV1::MakeCluster -- ");
324         printf("Analyzing chamber %d, plane %d, sector %d.\n"
325               ,icham,iplan,isect);
326
327         Int_t nRowMax     = geo->GetRowMax(iplan,icham,isect);
328         Int_t nColMax     = geo->GetColMax(iplan);
329         Int_t nTimeBefore = geo->GetTimeBefore();
330         Int_t nTimeTotal  = geo->GetTimeTotal();  
331
332         // Get the digits
333         digits = fDigitsManager->GetDigits(idet);
334         digits->Expand();
335         track0 = fDigitsManager->GetDictionary(idet,0);
336         track0->Expand();
337         track1 = fDigitsManager->GetDictionary(idet,1);
338         track1->Expand();
339         track2 = fDigitsManager->GetDictionary(idet,2); 
340         track2->Expand();
341
342         // Loop through the chamber and find the maxima 
343         for ( row = 0;  row <  nRowMax;    row++) {
344           for ( col = 2;  col <  nColMax;    col++) {
345             for (time = 0; time < nTimeTotal; time++) {
346
347               Int_t signalL = digits->GetDataUnchecked(row,col  ,time);
348               Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
349               Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
350  
351               // Look for the maximum
352               if (signalM >= maxThresh) {
353                 if (((signalL >= sigThresh) &&
354                      (signalL <  signalM))  ||
355                     ((signalR >= sigThresh) &&
356                      (signalR <  signalM))) {
357                   // Maximum found, mark the position by a negative signal
358                   digits->SetDataUnchecked(row,col-1,time,-signalM);
359                 }
360               }
361
362             }  
363           }    
364         }      
365
366         // Now check the maxima and calculate the cluster position
367         for ( row = 0;  row <  nRowMax  ;  row++) {
368           for (time = 0; time < nTimeTotal; time++) {
369             for ( col = 1;  col <  nColMax-1;  col++) {
370
371               // Maximum found ?             
372               if (digits->GetDataUnchecked(row,col,time) < 0) {
373
374                 Int_t iPad;
375                 for (iPad = 0; iPad < kNclus; iPad++) {
376                   Int_t iPadCol = col - 1 + iPad;
377                   clusterSignal[iPad]     = TMath::Abs(digits->GetDataUnchecked(row
378                                                                                ,iPadCol
379                                                                                ,time));
380                   clusterDigit[iPad]      = digits->GetIndexUnchecked(row,iPadCol,time);
381                   clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
382                   clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
383                   clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
384                 }
385
386                 // Count the number of pads in the cluster
387                 Int_t nPadCount = 0;
388                 Int_t ii        = 0;
389                 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii  ,time))
390                                                                   >= sigThresh) {
391                   nPadCount++;
392                   ii++;
393                   if (col-ii   <        0) break;
394                 }
395                 ii = 0;
396                 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
397                                                                   >= sigThresh) {
398                   nPadCount++;
399                   ii++;
400                   if (col+ii+1 >= nColMax) break;
401                 }
402
403                 nClusters++;
404                 switch (nPadCount) {
405                 case 2:
406                   iType = 0;
407                   nClusters2pad++;
408                   break;
409                 case 3:
410                   iType = 1;
411                   nClusters3pad++;
412                   break;
413                 case 4:
414                   iType = 2;
415                   nClusters4pad++;
416                   break;
417                 case 5:
418                   iType = 3;
419                   nClusters5pad++;
420                   break;
421                 default:
422                   iType = 4;
423                   nClustersLarge++;
424                   break;
425                 };
426
427                 // Don't analyze large clusters
428                 //if (iType == 4) continue;
429
430                 // Look for 5 pad cluster with minimum in the middle
431                 Bool_t fivePadCluster = kFALSE;
432                 if (col < nColMax-3) {
433                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
434                     fivePadCluster = kTRUE;
435                   }
436                   if ((fivePadCluster) && (col < nColMax-5)) {
437                     if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
438                       fivePadCluster = kFALSE;
439                     }
440                   }
441                   if ((fivePadCluster) && (col >         1)) {
442                     if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
443                       fivePadCluster = kFALSE;
444                     }
445                   }
446                 }
447
448                 // 5 pad cluster
449                 // Modify the signal of the overlapping pad for the left part 
450                 // of the cluster which remains from a previous unfolding
451                 if (iUnfold) {
452                   clusterSignal[0] *= ratioLeft;
453                   iType   = 3;
454                   iUnfold = 0;
455                 }
456
457                 // Unfold the 5 pad cluster
458                 if (fivePadCluster) {
459                   for (iPad = 0; iPad < kNsig; iPad++) {
460                     padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
461                                                                          ,col-1+iPad
462                                                                          ,time));
463                   }
464                   // Unfold the two maxima and set the signal on 
465                   // the overlapping pad to the ratio
466                   ratioRight        = Unfold(kEpsilon,padSignal);
467                   ratioLeft         = 1.0 - ratioRight; 
468                   clusterSignal[2] *= ratioRight;
469                   iType   = 3;
470                   iUnfold = 1;
471                 }
472
473                 Float_t clusterCharge = clusterSignal[0]
474                                       + clusterSignal[1]
475                                       + clusterSignal[2];
476                 
477                 // The position of the cluster
478                 clusterPads[0] = row + 0.5;
479                 // Take the shift of the additional time bins into account
480                 clusterPads[2] = time - nTimeBefore + 0.5;
481
482                 if (fUseLUT) {
483
484                   // Calculate the position of the cluster by using the
485                   // lookup table method
486                   Float_t ratioLUT;
487                   Float_t signLUT;
488                   Float_t lut = 0.0;
489                   if (clusterSignal[0] > clusterSignal[2]) {
490                     ratioLUT = clusterSignal[0] / clusterSignal[1];
491                     signLUT  = -1.0;
492                   }
493                   else {
494                     ratioLUT = clusterSignal[2] / clusterSignal[1];
495                     signLUT  =  1.0;
496                   }
497                   if      (ratioLUT < kLUTmin) {
498                     lut = 0.0;
499                   }
500                   else if (ratioLUT > kLUTmax) {
501                     lut = 0.5;
502                   }
503                   else {
504                     Int_t indexLUT = TMath::Nint ((kNlut-1) * (ratioLUT - kLUTmin)  
505                                                             / (kLUTmax  - kLUTmin)); 
506                     lut = fLUT[indexLUT];
507                   }
508                   clusterPads[1] = col + 0.5 + signLUT * lut;
509
510                 }
511                 else {
512
513                   // Calculate the position of the cluster by using the
514                   // center of gravity method
515                   clusterPads[1] = col + 0.5 
516                                  + (clusterSignal[2] - clusterSignal[0]) 
517                                  / clusterCharge;
518
519                 }
520
521                 if (fVerbose) {
522                   printf("-----------------------------------------------------------\n");
523                   printf("Create cluster no. %d\n",nClusters);
524                   printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
525                                                                     ,clusterPads[1]
526                                                                     ,clusterPads[2]);
527                   printf("Indices: %d, %d, %d\n",clusterDigit[0]
528                                                 ,clusterDigit[1]
529                                                 ,clusterDigit[2]);
530                   printf("Total charge = %f\n",clusterCharge);
531                   printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
532                                                     ,clusterTracks[1]
533                                                     ,clusterTracks[2]);
534                   printf("        pad1 %d, %d, %d\n",clusterTracks[3]
535                                                     ,clusterTracks[4]
536                                                     ,clusterTracks[5]);
537                   printf("        pad2 %d, %d, %d\n",clusterTracks[6]
538                                                     ,clusterTracks[7]
539                                                     ,clusterTracks[8]);
540                   printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
541                 }
542
543                 // Add the cluster to the output array 
544                 fTRD->AddCluster(clusterPads
545                                 ,clusterDigit
546                                 ,idet
547                                 ,clusterCharge
548                                 ,clusterTracks
549                                 ,iType);
550
551               }
552             } 
553           }   
554         }     
555
556         // Compress the arrays
557         digits->Compress(1,0);
558         track0->Compress(1,0);
559         track1->Compress(1,0);
560         track2->Compress(1,0);
561
562         // Write the cluster and reset the array
563         WriteClusters(idet);
564         fTRD->ResetRecPoints();
565
566         printf("AliTRDclusterizerV1::MakeCluster -- ");
567         printf("Found %d clusters in total.\n"
568               ,nClusters);
569         printf("                                    2pad:  %d\n",nClusters2pad);
570         printf("                                    3pad:  %d\n",nClusters3pad);
571         printf("                                    4pad:  %d\n",nClusters4pad);
572         printf("                                    5pad:  %d\n",nClusters5pad);
573         printf("                                    Large: %d\n",nClustersLarge);
574
575       }    
576     }      
577   }        
578
579   printf("AliTRDclusterizerV1::MakeCluster -- ");
580   printf("Done.\n");
581
582   return kTRUE;
583
584 }
585
586 //_____________________________________________________________________________
587 Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
588 {
589   //
590   // Method to unfold neighbouring maxima.
591   // The charge ratio on the overlapping pad is calculated
592   // until there is no more change within the range given by eps.
593   // The resulting ratio is then returned to the calling method.
594   //
595
596   Int_t   itStep            = 0;      // Count iteration steps
597
598   Float_t ratio             = 0.5;    // Start value for ratio
599   Float_t prevRatio         = 0;      // Store previous ratio
600
601   Float_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
602   Float_t newRightSignal[3] = {0};    // Array to store right cluster signal
603
604   // Start the iteration
605   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
606
607     itStep++;
608     prevRatio = ratio;
609
610     // Cluster position according to charge ratio
611     Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
612                      / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
613     Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
614                      / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
615
616     // Set cluster charge ratio
617     Float_t ampLeft  = padSignal[1];
618     Float_t ampRight = padSignal[3];
619
620     // Apply pad response to parameters
621     newLeftSignal[0]  = ampLeft  * PadResponse(-1 - maxLeft);
622     newLeftSignal[1]  = ampLeft  * PadResponse( 0 - maxLeft);
623     newLeftSignal[2]  = ampLeft  * PadResponse( 1 - maxLeft);
624
625     newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
626     newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
627     newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
628
629     // Calculate new overlapping ratio
630     ratio = TMath::Min(1.0,newLeftSignal[2] / 
631                           (newLeftSignal[2] + newRightSignal[0]));
632
633   }
634
635   return ratio;
636
637 }
638
639 //_____________________________________________________________________________
640 Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
641 {
642   //
643   // The pad response for the chevron pads. 
644   // We use a simple Gaussian approximation which should be good
645   // enough for our purpose.
646   // Updated for new PRF 1/5/01.
647   //
648
649   // The parameters for the response function
650   const Float_t kA  =  0.8303; 
651   const Float_t kB  = -0.00392;
652   const Float_t kC  =  0.472 * 0.472;
653   const Float_t kD  =  2.19;
654
655   Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
656
657   return (pr);
658
659 }