PHOS/AliPHOSTrackSegmentMakerv1.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 // Algorithm class to construct track segments connection RecPoints in
  18 // EMCA and Ppsd. Unfolds also the clusters in EMCA.
  19 //*-- Author : D. Peressounko  SUBATECH
  20 //////////////////////////////////////////////////////////////////////////////
  21
  22 // --- ROOT system ---
  23
  24 #include "TObjArray.h"
  25 #include "TClonesArray.h"
  26 #include "TObjectTable.h"
  27
  28 // --- Standard library ---
  29
  30 #include <iostream>
  31 #include <cassert>
  32
  33 // --- AliRoot header files ---
  34
  35 #include "AliPHOSTrackSegmentMakerv1.h"
  36 #include "AliPHOSTrackSegment.h"
  37 #include "AliPHOSLink.h"
  38 #include "AliPHOSv0.h"
  39 #include "AliRun.h"
  40
  41 extern void UnfoldingChiSquare(Int_t &nPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) ;
  42
  43 ClassImp( AliPHOSTrackSegmentMakerv1)
  44
  45
  46 //____________________________________________________________________________
  47  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1()
  48 {
  49   // ctor
  50   fR0 = 4. ;
  51   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
  52   //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0),
  53   fDelta = fR0 + geom->GetCrystalSize(0) ;
  54   fMinuit = new TMinuit(100) ;
  55 }
  56
  57 //____________________________________________________________________________
  58  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
  59 {
  60   // dtor
  61   delete fMinuit ;
  62 }
  63 //____________________________________________________________________________
  64 Bool_t  AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
  65                                     Int_t nPar, Float_t * fitparameters)
  66 {
  67   // Calls TMinuit for fitting cluster with several maxima
  68
  69   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
  70
  71   gMinuit->SetPrintLevel(-1) ;           // No Printout
  72   gMinuit->SetFCN(UnfoldingChiSquare) ;  // To set the address of the minimization function
  73   gMinuit->SetObjectFit(emcRP) ;         // To tranfer pointer to UnfoldingChiSquare
  74
  75   // filling initial values for fit parameters
  76   AliPHOSDigit * digit ;
  77
  78   Int_t ierflg  = 0;
  79   Int_t index   = 0 ;
  80   Int_t nDigits = (Int_t) nPar / 3 ;
  81
  82   Int_t iDigit ;
  83
  84
  85   for(iDigit = 0; iDigit < nDigits; iDigit++){
  86     digit = (AliPHOSDigit *) maxAt[iDigit];
  87
  88     Int_t relid[4] ;
  89     Float_t x ;
  90     Float_t z ;
  91     geom->AbsToRelNumbering(digit->GetId(), relid) ;
  92     geom->RelPosInModule(relid, x, z) ;
  93
  94     Float_t energy = maxAtEnergy[iDigit] ;
  95
  96     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
  97     index++ ;
  98     if(ierflg != 0){
  99       cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
 100       return kFALSE;
 101     }
 102     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 103     index++ ;
 104     if(ierflg != 0){
 105       cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
 106       return kFALSE;
 107     }
 108     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 109     index++ ;
 110     if(ierflg != 0){
 111       cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : energy = " << energy << endl ;
 112       return kFALSE;
 113     }
 114   }
 115
 116   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 117                       //  depends on it.
 118   Double_t p1 = 1.0 ;
 119   Double_t p2 = 0.0 ;
 120
 121   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TgMinuit to reduce function calls
 122   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 123   gMinuit->SetMaxIterations(5);
 124   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 125   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 126   if(ierflg == 4){  // Minimum not found
 127     cout << "PHOS Unfolding>  Fit not converged, cluster abandoned "<< endl ;
 128     return kFALSE ;
 129   }
 130   for(index = 0; index < nPar; index++){
 131     Double_t err ;
 132     Double_t val ;
 133     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 134     fitparameters[index] = val ;
 135    }
 136   return kTRUE;
 137
 138 }
 139
 140 //____________________________________________________________________________
 141 void  AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * Dl, RecPointsList * emcIn, TObjArray * emcOut,
 142                                         RecPointsList * ppsdIn, TObjArray * ppsdOutUp,
 143                                         TObjArray * ppsdOutLow, Int_t & phosmod, Int_t & emcStopedAt,
 144                                         Int_t & ppsdStopedAt)
 145 {
 146   // Unfold clusters and fill xxxOut arrays with clusters from one PHOS module
 147
 148   AliPHOSEmcRecPoint *  emcRecPoint  ;
 149   AliPHOSPpsdRecPoint * ppsdRecPoint ;
 150   Int_t index ;
 151
 152   Int_t nEmcUnfolded = emcIn->GetEntries() ;
 153   for(index = emcStopedAt; index < nEmcUnfolded; index++){
 154     emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ;
 155
 156     if(emcRecPoint->GetPHOSMod() != phosmod )
 157        break ;
 158
 159     Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 160     Int_t maxAt[nMultipl] ;
 161     Float_t maxAtEnergy[nMultipl] ;
 162     Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ;
 163
 164     if(nMax <= 1)     // if cluster is very flat (no pronounced maximum) then nMax = 0
 165       emcOut->Add(emcRecPoint) ;
 166     else {
 167       UnfoldClusters(Dl, emcIn, emcRecPoint, nMax, maxAt, maxAtEnergy, emcOut) ;
 168       emcIn->Remove(emcRecPoint);
 169       emcIn->Compress() ;
 170       nEmcUnfolded-- ;
 171       index-- ;
 172     }
 173   }
 174   emcStopedAt = index ;
 175
 176   for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){
 177     ppsdRecPoint = (AliPHOSPpsdRecPoint *) (*ppsdIn)[index] ;
 178     if(ppsdRecPoint->GetPHOSMod() != phosmod )
 179       break ;
 180     if(ppsdRecPoint->GetUp() )
 181       ppsdOutUp->Add(ppsdRecPoint) ;
 182     else
 183       ppsdOutLow->Add(ppsdRecPoint) ;
 184   }
 185   ppsdStopedAt = index ;
 186
 187   emcOut->Sort() ;
 188   ppsdOutUp->Sort() ;
 189   ppsdOutLow->Sort() ;
 190 }
 191 //____________________________________________________________________________
 192 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar)
 193 {
 194   Float_t r = fR0 ;
 195
 196   TVector3 vecEmc ;
 197   TVector3 vecPpsd ;
 198
 199   emcclu->GetLocalPosition(vecEmc) ;
 200   PpsdClu->GetLocalPosition(vecPpsd)  ;
 201   if(emcclu->GetPHOSMod() == PpsdClu->GetPHOSMod()){
 202     if(vecPpsd.X() >= vecEmc.X() - fDelta ){
 203       if(vecPpsd.Z() >= vecEmc.Z() - fDelta ){
 204         AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 205         // Correct to difference in CPV and EMC position due to different distance to center.
 206         // we assume, that particle moves from center
 207         Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 208         Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 209         dEMC         = dEMC / dCPV ;
 210         vecPpsd = dEMC * vecPpsd  - vecEmc ;
 211         r = vecPpsd.Mag() ;
 212       } // if  zPpsd >= zEmc - fDelta
 213       toofar = kFALSE ;
 214     } // if  xPpsd >= xEmc - fDelta
 215     else
 216       toofar = kTRUE ;
 217   }
 218   else
 219     toofar = kTRUE ;
 220
 221   return r ;
 222 }
 223
 224 //____________________________________________________________________________
 225 void  AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp,
 226                                      TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray,
 227                                      TClonesArray *linkupArray)
 228 {
 229   //Finds distanses (links) between all EMC and PPSD clusters, which are not further from each other than fR0
 230
 231   TIter nextEmc(emcRecPoints) ;
 232   Int_t iEmcClu = 0 ;
 233
 234   AliPHOSPpsdRecPoint * ppsdlow ;
 235   AliPHOSPpsdRecPoint * ppsdup ;
 236   AliPHOSEmcRecPoint * emcclu ;
 237
 238   Int_t iLinkLow = 0 ;
 239   Int_t iLinkUp  = 0 ;
 240
 241   while( (emcclu = (AliPHOSEmcRecPoint*)nextEmc() ) ) {
 242     Bool_t toofar ;
 243     TIter nextPpsdLow(ppsdRecPointsLow ) ;
 244     Int_t iPpsdLow = 0 ;
 245
 246     while( (ppsdlow = (AliPHOSPpsdRecPoint*)nextPpsdLow() ) ) {
 247       Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdlow, toofar) ;
 248
 249       if(toofar)
 250         break ;
 251       if(r < fR0)
 252         new( (*linklowArray)[iLinkLow++]) AliPHOSLink(r, iEmcClu, iPpsdLow) ;
 253
 254       iPpsdLow++ ;
 255
 256     }
 257
 258     TIter nextPpsdUp(ppsdRecPointsUp ) ;
 259     Int_t iPpsdUp = 0 ;
 260
 261     while( (ppsdup = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) {
 262       Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdup, toofar) ;
 263
 264       if(toofar)
 265         break ;
 266       if(r < fR0)
 267         new( (*linkupArray)[iLinkUp++]) AliPHOSLink(r, iEmcClu, iPpsdUp) ;
 268
 269       iPpsdUp++ ;
 270
 271     }
 272
 273     iEmcClu++ ;
 274
 275   } // while nextEmC
 276
 277   linklowArray->Sort() ; //first links with smallest distances
 278   linkupArray->Sort() ;
 279 }
 280
 281 //____________________________________________________________________________
 282 void  AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp,
 283                                     TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray,
 284                                     TClonesArray * linkupArray, TrackSegmentsList * trsl)
 285 {
 286
 287   // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance
 288
 289   TIter nextLow(linklowArray) ;
 290   TIter nextUp(linkupArray) ;
 291
 292   AliPHOSLink * linkLow ;
 293   AliPHOSLink * linkUp ;
 294
 295   AliPHOSEmcRecPoint * emc ;
 296   AliPHOSPpsdRecPoint * ppsdLow ;
 297   AliPHOSPpsdRecPoint * ppsdUp ;
 298
 299   AliPHOSRecPoint * nullpointer = 0 ;
 300
 301   while ( (linkLow =  (AliPHOSLink *)nextLow() ) ){
 302     emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(linkLow->GetEmc()) ;
 303     ppsdLow = (AliPHOSPpsdRecPoint *) ppsdRecPointsLow->At(linkLow->GetPpsd()) ;
 304     if( (emc) && (ppsdLow) ){ // RecPoints not removed yet
 305          ppsdUp = 0 ;
 306
 307          while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){
 308            if(linkLow->GetEmc() == linkUp->GetEmc() ){
 309              ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ;
 310              break ;
 311            }
 312
 313          } // while nextUp
 314
 315          nextUp.Reset();
 316          AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
 317          trsl->Add(subtr) ;
 318          emcRecPoints->AddAt(nullpointer,linkLow->GetEmc()) ;
 319          ppsdRecPointsLow->AddAt(nullpointer,linkLow->GetPpsd()) ;
 320
 321          if(ppsdUp)
 322            ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
 323
 324     }
 325   }
 326
 327   TIter nextEmc(emcRecPoints) ;
 328   nextEmc.Reset() ;
 329
 330   while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no ppsdlow
 331     ppsdLow = 0 ;
 332     ppsdUp  = 0 ;
 333
 334     while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){
 335
 336       if(emcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){
 337         ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ;
 338         break ;
 339       }
 340
 341     }
 342     AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
 343     trsl->Add(subtr) ;
 344
 345     if(ppsdUp)
 346       ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
 347   }
 348
 349 }
 350
 351 //____________________________________________________________________________
 352 void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl,
 353                                         RecPointsList * ppsdl, TrackSegmentsList * trsl)
 354 {
 355   // main function, does the job
 356
 357   Int_t phosmod      = 1 ;
 358   Int_t emcStopedAt  = 0 ;
 359   Int_t ppsdStopedAt = 0 ;
 360
 361   TObjArray * emcRecPoints     = new TObjArray(100) ;  // these arrays keep pointers
 362   TObjArray * ppsdRecPointsUp  = new TObjArray(100) ;  // to RecPoints, which are
 363   TObjArray * ppsdRecPointsLow = new TObjArray(100) ;  // kept in TClonesArray's emcl and ppsdl
 364
 365
 366   TClonesArray * linklowArray = new TClonesArray("AliPHOSLink", 100);
 367   TClonesArray * linkupArray  = new TClonesArray("AliPHOSLink", 100);
 368
 369   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 370
 371   while(phosmod <= geom->GetNModules() ){
 372
 373     FillOneModule(DL, emcl, emcRecPoints, ppsdl, ppsdRecPointsUp, ppsdRecPointsLow, phosmod, emcStopedAt, ppsdStopedAt) ;
 374
 375     MakeLinks(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray) ;
 376
 377     MakePairs(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray, trsl) ;
 378
 379     emcRecPoints->Clear() ;
 380
 381     ppsdRecPointsUp->Clear() ;
 382
 383     ppsdRecPointsLow->Clear() ;
 384
 385     linkupArray->Clear() ;
 386
 387     linklowArray->Clear() ;
 388
 389     phosmod++ ;
 390   }
 391   delete emcRecPoints ;
 392   emcRecPoints = 0 ;
 393
 394   delete ppsdRecPointsUp ;
 395   ppsdRecPointsUp = 0 ;
 396
 397   delete ppsdRecPointsLow ;
 398   ppsdRecPointsLow = 0 ;
 399
 400   delete linkupArray ;
 401   linkupArray = 0  ;
 402
 403   delete linklowArray ;
 404   linklowArray = 0 ;
 405 }
 406
 407 //____________________________________________________________________________
 408 Double_t  AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r)
 409 {
 410 // If you change this function, change also gradiend evaluation  in ChiSquare()
 411   Double_t r4    = r*r*r*r ;
 412   Double_t r295  = TMath::Power(r, 2.95) ;
 413   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 414   return shape ;
 415 }
 416
 417 //____________________________________________________________________________
 418 void  AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * DL, RecPointsList * emcIn,  AliPHOSEmcRecPoint * iniEmc,
 419                                          Int_t nMax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList)
 420 {
 421   // fits cluster with nMax overlapping showers
 422
 423   Int_t nPar = 3 * nMax ;
 424   Float_t fitparameters[nPar] ;
 425   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 426
 427   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 428   if( !rv )  // Fit failed, return and remove cluster
 429     return ;
 430
 431   Float_t xDigit ;
 432   Float_t zDigit ;
 433   Int_t relid[4] ;
 434
 435   Int_t nDigits = iniEmc->GetMultiplicity() ;
 436   Float_t xpar  ;
 437   Float_t zpar  ;
 438   Float_t epar  ;
 439   Float_t distance ;
 440   Float_t ratio ;
 441   Float_t efit[nDigits] ;
 442   Int_t iparam ;
 443   Int_t iDigit ;
 444
 445   AliPHOSDigit * digit ;
 446   AliPHOSEmcRecPoint * emcRP ;
 447   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 448   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 449
 450   Int_t iRecPoint = emcIn->GetEntries() ;
 451
 452   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 453     digit = (AliPHOSDigit *) emcDigits[iDigit];
 454     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 455     geom->RelPosInModule(relid, xDigit, zDigit) ;
 456     efit[iDigit] = 0;
 457     iparam = 0 ;
 458
 459     while(iparam < nPar ){
 460       xpar = fitparameters[iparam] ;
 461       zpar = fitparameters[iparam+1] ;
 462       epar = fitparameters[iparam+2] ;
 463       iparam += 3 ;
 464       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 465       distance =  TMath::Sqrt(distance) ;
 466       efit[iDigit] += epar * ShowerShape(distance) ;
 467     }
 468   }
 469
 470   iparam = 0 ;
 471   Float_t eDigit ;
 472
 473   while(iparam < nPar ){
 474     xpar = fitparameters[iparam] ;
 475     zpar = fitparameters[iparam+1] ;
 476     epar = fitparameters[iparam+2] ;
 477     iparam += 3 ;
 478     new ((*emcIn)[iRecPoint]) AliPHOSEmcRecPoint( iniEmc->GetLogWeightCut(), iniEmc->GetLocMaxCut() ) ;
 479     emcRP = (AliPHOSEmcRecPoint *) emcIn->At(iRecPoint++);
 480
 481     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 482       digit = (AliPHOSDigit *) emcDigits[iDigit];
 483       geom->AbsToRelNumbering(digit->GetId(), relid) ;
 484       geom->RelPosInModule(relid, xDigit, zDigit) ;
 485       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 486       distance =  TMath::Sqrt(distance) ;
 487       ratio = epar * ShowerShape(distance) / efit[iDigit] ;
 488       eDigit = emcEnergies[iDigit] * ratio ;
 489       emcRP->AddDigit( *digit, eDigit ) ;
 490     }
 491
 492     emcList->Add(emcRP) ;
 493
 494   }
 495 }
 496
 497 //______________________________________________________________________________
 498 void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 499 {
 500
 501 // Number of parameters, Gradient, Chi squared, parameters, what to do
 502
 503   AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 504
 505   AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
 506   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 507   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 508   fret = 0. ;
 509   Int_t iparam ;
 510
 511   if(iflag == 2)
 512     for(iparam = 0 ; iparam < nPar ; iparam++)
 513       Grad[iparam] = 0 ; // Will evaluate gradient
 514
 515   Double_t efit ;
 516
 517   AliPHOSDigit * digit ;
 518   Int_t iDigit = 0 ;
 519
 520   while ( (digit = (AliPHOSDigit *)emcDigits[iDigit] )){
 521     Int_t relid[4] ;
 522     Float_t xDigit ;
 523     Float_t zDigit ;
 524     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 525     geom->RelPosInModule(relid, xDigit, zDigit) ;
 526
 527      if(iflag == 2){  // calculate gradient
 528        Int_t iParam = 0 ;
 529        efit = 0 ;
 530        while(iParam < nPar ){
 531          Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
 532          iParam++ ;
 533          distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
 534          distance = TMath::Sqrt( distance ) ;
 535          iParam++ ;
 536          efit += x[iParam] * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
 537          iParam++ ;
 538        }
 539        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 540        iParam = 0 ;
 541        while(iParam < nPar ){
 542          Double_t xpar = x[iParam] ;
 543          Double_t zpar = x[iParam+1] ;
 544          Double_t epar = x[iParam+2] ;
 545          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 546          Double_t shape = sum * AliPHOSTrackSegmentMakerv1::ShowerShape(dr) ;
 547          Double_t r4 = dr*dr*dr*dr ;
 548          Double_t r295 = TMath::Power(dr,2.95) ;
 549          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 550                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 551
 552          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 553          iParam++ ;
 554          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 555          iParam++ ;
 556          Grad[iParam] += shape ;                                  // Derivative over energy
 557          iParam++ ;
 558        }
 559      }
 560      efit = 0;
 561      iparam = 0 ;
 562      while(iparam < nPar ){
 563        Double_t xpar = x[iparam] ;
 564        Double_t zpar = x[iparam+1] ;
 565        Double_t epar = x[iparam+2] ;
 566        iparam += 3 ;
 567        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 568        distance =  TMath::Sqrt(distance) ;
 569        efit += epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
 570      }
 571      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 572      // Here we assume, that sigma = sqrt(E)
 573      iDigit++ ;
 574   }
 575 }