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