EMCAL/AliEMCALUnfolding.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 //_________________________________________________________________________\r
  17 //  Base class for the cluster unfolding algorithm \r
  18 //*-- Author: Adam Matyja (SUBATECH)\r
  19 //  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis\r
  20 //-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1\r
  21 //-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)\r
  22 //\r
  23 //  unfolds the clusters having several local maxima. \r
  24 //////////////////////////////////////////////////////////////////////////////\r
  25 \r
  26 // --- ROOT system ---\r
  27 #include "TClonesArray.h"\r
  28 //#include "TTree.h"\r
  29 //#include <TFile.h> \r
  30 //class TFolder;\r
  31 #include <TMath.h> \r
  32 #include <TMinuit.h>\r
  33 //#include <TTree.h> \r
  34 //class TSystem; \r
  35 //#include <TBenchmark.h>\r
  36 //#include <TBrowser.h>\r
  37 //#include <TROOT.h>\r
  38 \r
  39 // --- Standard library ---\r
  40 #include <cassert>\r
  41 \r
  42 // --- AliRoot header files ---\r
  43 #include "AliEMCALUnfolding.h"\r
  44 #include "AliEMCALGeometry.h"\r
  45 #include "AliRunLoader.h"\r
  46 #include "AliRun.h"\r
  47 #include "AliEMCAL.h"\r
  48 #include "AliEMCALRecParam.h"\r
  49 #include "AliEMCALRecPoint.h"\r
  50 #include "AliEMCALDigit.h"\r
  51 #include "AliEMCALReconstructor.h"\r
  52 //#include "AliEMCALClusterizer.h"\r
  53 \r
  54 \r
  55 \r
  56 #include "AliLog.h"\r
  57 \r
  58 #include "AliCDBManager.h"\r
  59 //#include "AliCaloCalibPedestal.h"\r
  60 //#include "AliEMCALCalibData.h"\r
  61 class AliCDBStorage;\r
  62 #include "AliCDBEntry.h"\r
  63 \r
  64 Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r
  65 Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r
  66 Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};\r
  67 \r
  68 ClassImp(AliEMCALUnfolding)\r
  69   \r
  70 //____________________________________________________________________________\r
  71 AliEMCALUnfolding::AliEMCALUnfolding():\r
  72   fNumberOfECAClusters(0),\r
  73   fECALocMaxCut(0),\r
  74   fGeom(NULL),\r
  75   fRecPoints(NULL),\r
  76   fDigitsArr(NULL)\r
  77 {\r
  78   // ctor with the indication of the file where header Tree and digits Tree are stored\r
  79  \r
  80   Init() ;\r
  81 }\r
  82 \r
  83 //____________________________________________________________________________\r
  84 AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):\r
  85   fNumberOfECAClusters(0),\r
  86   fECALocMaxCut(0),\r
  87   fGeom(geometry),\r
  88   fRecPoints(NULL),\r
  89   fDigitsArr(NULL)\r
  90 {\r
  91   // ctor with the indication of the file where header Tree and digits Tree are stored\r
  92   // use this contructor to avoid usage of Init() which uses runloader\r
  93   // change needed by HLT - MP\r
  94   if (!fGeom)\r
  95   {\r
  96     AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
  97   }\r
  98   \r
  99 }\r
 100 \r
 101 //____________________________________________________________________________\r
 102 AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):\r
 103   fNumberOfECAClusters(0),\r
 104   fECALocMaxCut(ECALocMaxCut),\r
 105   fGeom(geometry),\r
 106   fRecPoints(NULL),\r
 107   fDigitsArr(NULL)\r
 108 {\r
 109   // ctor with the indication of the file where header Tree and digits Tree are stored\r
 110   // use this contructor to avoid usage of Init() which uses runloader\r
 111   // change needed by HLT - MP\r
 112   if (!fGeom)\r
 113   {\r
 114     AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
 115   }\r
 116   Int_t i=0;\r
 117   for (i = 0; i < 8; i++) fSSPars[i] = SSPars[i];\r
 118   for (i = 0; i < 3; i++) {\r
 119     fPar5[i] = Par5[i];\r
 120     fPar6[i] = Par6[i];\r
 121   }\r
 122 \r
 123 }\r
 124 \r
 125 //____________________________________________________________________________\r
 126 void AliEMCALUnfolding::Init()\r
 127 {\r
 128   // Make all memory allocations which can not be done in default constructor.\r
 129   // Attach the Clusterizer task to the list of EMCAL tasks\r
 130 \r
 131   AliRunLoader *rl = AliRunLoader::Instance();\r
 132   if (rl && rl->GetAliRun()){\r
 133     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));\r
 134     if(emcal)fGeom = emcal->GetGeometry();\r
 135   }\r
 136   \r
 137   if(!fGeom)\r
 138     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());\r
 139   \r
 140   AliDebug(1,Form("geom %p",fGeom));\r
 141   \r
 142   if(!gMinuit) \r
 143     gMinuit = new TMinuit(100) ;\r
 144   \r
 145 }\r
 146 \r
 147 //____________________________________________________________________________\r
 148   AliEMCALUnfolding::~AliEMCALUnfolding()\r
 149 {\r
 150   // dtor\r
 151 }\r
 152 \r
 153 //____________________________________________________________________________\r
 154 void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)\r
 155 {\r
 156   //\r
 157   //Set input for unfolding purposes\r
 158   SetNumberOfECAClusters(numberOfECAClusters);\r
 159   SetRecPoints(recPoints);\r
 160   SetDigitsArr(digitsArr);\r
 161 }\r
 162 \r
 163 //____________________________________________________________________________\r
 164 void AliEMCALUnfolding::MakeUnfolding()\r
 165 {\r
 166   // Unfolds clusters using the shape of an ElectroMagnetic shower\r
 167   // Performs unfolding of all clusters\r
 168   \r
 169   if(fNumberOfECAClusters > 0){\r
 170     if (fGeom==0)\r
 171       AliFatal("Did not get geometry from EMCALLoader") ;\r
 172     //Int_t nModulesToUnfold = fGeom->GetNCells();\r
 173     \r
 174     Int_t numberofNotUnfolded = fNumberOfECAClusters ;\r
 175     Int_t index ;\r
 176     for(index = 0 ; index < numberofNotUnfolded ; index++){\r
 177       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r
 178       if(recPoint){\r
 179         //do we really need it?\r
 180         //      TVector3 gpos;\r
 181         //      Int_t absId = -1;\r
 182         //      recPoint->GetGlobalPosition(gpos);\r
 183         //      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r
 184         //      if(absId > nModulesToUnfold)\r
 185         //        break ;\r
 186         \r
 187         Int_t nMultipl = recPoint->GetMultiplicity() ;\r
 188         AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
 189         Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
 190         Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
 191         \r
 192         if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0\r
 193           if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
 194             fRecPoints->Remove(recPoint);\r
 195             fRecPoints->Compress() ;//is it really needed\r
 196             index-- ;\r
 197             fNumberOfECAClusters-- ;\r
 198             numberofNotUnfolded-- ;\r
 199           }\r
 200         }\r
 201         else{\r
 202           recPoint->SetNExMax(1) ; //Only one local maximum\r
 203         }\r
 204         \r
 205         delete[] maxAt ;\r
 206         delete[] maxAtEnergy ;\r
 207       } else AliError("RecPoint NULL");\r
 208     } // rec point loop\r
 209   }\r
 210   // End of Unfolding of clusters\r
 211 }\r
 212 \r
 213 //____________________________________________________________________________\r
 214 Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, \r
 215                                           Int_t nMax, \r
 216                                           AliEMCALDigit ** maxAt, \r
 217                                           Float_t * maxAtEnergy)\r
 218 {\r
 219   // Extended to whole EMCAL \r
 220   // Performs the unfolding of a cluster with nMax overlapping showers \r
 221   \r
 222   Int_t nPar = 3 * nMax ;\r
 223   Float_t * fitparameters = new Float_t[nPar] ;\r
 224   \r
 225   if (fGeom==0)\r
 226     AliFatal("Did not get geometry from EMCALLoader") ;\r
 227   \r
 228   Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
 229   if( !rv ) {\r
 230     // Fit failed, return (and remove cluster? - why? I leave the cluster)\r
 231     iniTower->SetNExMax(-1) ;\r
 232     delete[] fitparameters ;\r
 233     return kFALSE;\r
 234   }\r
 235   \r
 236   // create unfolded rec points and fill them with new energy lists\r
 237   // First calculate energy deposited in each sell in accordance with\r
 238   // fit (without fluctuations): efit[]\r
 239   // and later correct this number in acordance with actual energy\r
 240   // deposition\r
 241   \r
 242   Int_t nDigits = iniTower->GetMultiplicity() ;\r
 243   Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
 244   Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r
 245   \r
 246   AliEMCALDigit * digit = 0 ;\r
 247   Int_t * digitsList = iniTower->GetDigitsList() ;\r
 248   \r
 249   Int_t iSupMod =  0 ;\r
 250   Int_t iTower  =  0 ;\r
 251   Int_t iIphi   =  0 ;\r
 252   Int_t iIeta   =  0 ;\r
 253   Int_t iphi    =  0 ;//x direction\r
 254   Int_t ieta    =  0 ;//z direstion\r
 255   \r
 256   Int_t iparam = 0 ;\r
 257   Int_t iDigit = 0 ;\r
 258   \r
 259   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
 260     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
 261     if(digit){\r
 262       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
 263       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
 264                                          iIphi, iIeta,iphi,ieta);\r
 265       EvalParsPhiDependence(digit->GetId(),fGeom);\r
 266       \r
 267       efit[iDigit] = 0.;\r
 268       iparam = 0;\r
 269       while(iparam < nPar ){\r
 270         xpar = fitparameters[iparam] ;\r
 271         zpar = fitparameters[iparam+1] ;\r
 272         epar = fitparameters[iparam+2] ;\r
 273         iparam += 3 ;\r
 274         \r
 275         efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
 276       }\r
 277     } else AliError("Digit NULL!");\r
 278     \r
 279   }//digit loop\r
 280   \r
 281   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
 282   // so that energy deposited in each cell is distributed between new clusters proportionally\r
 283   // to its contribution to efit\r
 284   \r
 285   Float_t * energiesList = iniTower->GetEnergiesList() ;\r
 286   Float_t ratio = 0 ;\r
 287   \r
 288   iparam = 0 ;\r
 289   while(iparam < nPar ){\r
 290     xpar = fitparameters[iparam] ;\r
 291     zpar = fitparameters[iparam+1] ;\r
 292     epar = fitparameters[iparam+2] ;\r
 293     iparam += 3 ;\r
 294     \r
 295     AliEMCALRecPoint * recPoint = 0 ;\r
 296     \r
 297     if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
 298       fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
 299     \r
 300     //add recpoint\r
 301     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
 302     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
 303     \r
 304     if(recPoint){\r
 305       \r
 306       fNumberOfECAClusters++ ;\r
 307       recPoint->SetNExMax((Int_t)nPar/3) ;\r
 308       \r
 309       Float_t eDigit = 0. ;\r
 310       for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
 311         digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
 312         if(digit){\r
 313           fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
 314           fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
 315                                              iIphi, iIeta,iphi,ieta);\r
 316           EvalParsPhiDependence(digit->GetId(),fGeom);\r
 317           if(efit[iDigit]==0) continue;//just for sure\r
 318           ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
 319           eDigit = energiesList[iDigit] * ratio ;\r
 320           recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
 321         } else AliError("NULL digit");\r
 322       }//digit loop \r
 323     } else AliError("NULL RecPoint");\r
 324   }//while\r
 325   \r
 326   delete[] fitparameters ;\r
 327   delete[] efit ;\r
 328   \r
 329   return kTRUE;\r
 330 }\r
 331 \r
 332 //____________________________________________________________________________\r
 333 Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, \r
 334                                         const Float_t* maxAtEnergy,\r
 335                                         Int_t nPar, Float_t * fitparameters) const\r
 336 {\r
 337   // Calls TMinuit to fit the energy distribution of a cluster with several maxima\r
 338   // The initial values for fitting procedure are set equal to the\r
 339   // positions of local maxima.       \r
 340   // Cluster will be fitted as a superposition of nPar/3\r
 341   // electromagnetic showers\r
 342 \r
 343   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r
 344         \r
 345   if(!gMinuit)\r
 346     gMinuit = new TMinuit(100) ;//max 100 parameters\r
 347 \r
 348   gMinuit->mncler();                     // Reset Minuit's list of paramters\r
 349   gMinuit->SetPrintLevel(-1) ;           // No Printout\r
 350   gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;\r
 351   // To set the address of the minimization function\r
 352   TList * toMinuit = new TList();\r
 353   toMinuit->AddAt(recPoint,0) ;\r
 354   toMinuit->AddAt(fDigitsArr,1) ;\r
 355   toMinuit->AddAt(fGeom,2) ;\r
 356 \r
 357   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare\r
 358 \r
 359   // filling initial values for fit parameters\r
 360   AliEMCALDigit * digit ;\r
 361 \r
 362   Int_t ierflg  = 0;\r
 363   Int_t index   = 0 ;\r
 364   Int_t nDigits = (Int_t) nPar / 3 ;\r
 365 \r
 366   Int_t iDigit ;\r
 367 \r
 368   Int_t iSupMod =  0 ;\r
 369   Int_t iTower  =  0 ;\r
 370   Int_t iIphi   =  0 ;\r
 371   Int_t iIeta   =  0 ;\r
 372   Int_t iphi    =  0 ;//x direction\r
 373   Int_t ieta    =  0 ;//z direstion\r
 374 \r
 375   for(iDigit = 0; iDigit < nDigits; iDigit++){\r
 376     digit = maxAt[iDigit];\r
 377     fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
 378     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
 379                                        iIphi, iIeta,iphi,ieta);\r
 380 \r
 381     Float_t energy = maxAtEnergy[iDigit] ;\r
 382 \r
 383     //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original\r
 384     gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;\r
 385     index++ ;\r
 386     if(ierflg != 0){\r
 387       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;\r
 388       toMinuit->Clear();\r
 389       delete toMinuit ;\r
 390       return kFALSE;\r
 391     }\r
 392     //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original\r
 393     gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;\r
 394     index++ ;\r
 395     if(ierflg != 0){\r
 396       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;\r
 397       toMinuit->Clear();\r
 398       delete toMinuit ;\r
 399       return kFALSE;\r
 400     }\r
 401     //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original\r
 402     gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r
 403     index++ ;\r
 404     if(ierflg != 0){\r
 405       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;\r
 406       toMinuit->Clear();\r
 407       delete toMinuit ;\r
 408       return kFALSE;\r
 409     }\r
 410   }\r
 411 \r
 412   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; \r
 413                       // The number of function call slightly depends on it.\r
 414   //  Double_t p1 = 1.0 ;// par to gradient \r
 415   Double_t p2 = 0.0 ;\r
 416   //  Double_t p3 = 3.0 ;\r
 417   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls\r
 418   //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient\r
 419   gMinuit->SetMaxIterations(5);//was 5\r
 420   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings\r
 421   //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts\r
 422 \r
 423   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize\r
 424   //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize\r
 425   if(ierflg == 4){  // Minimum not found\r
 426     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;\r
 427     toMinuit->Clear();\r
 428     delete toMinuit ;\r
 429     return kFALSE ;\r
 430   }\r
 431   for(index = 0; index < nPar; index++){\r
 432     Double_t err = 0. ;\r
 433     Double_t val = 0. ;\r
 434     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index\r
 435     fitparameters[index] = val ;\r
 436   }\r
 437 \r
 438   toMinuit->Clear();\r
 439   delete toMinuit ;\r
 440   return kTRUE;\r
 441 \r
 442 }\r
 443 \r
 444 //____________________________________________________________________________\r
 445 Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)\r
 446 { \r
 447   // extended to whole EMCAL \r
 448   // Shape of the shower\r
 449   // If you change this function, change also the gradient evaluation in ChiSquare()\r
 450 \r
 451   Double_t r = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r
 452   Double_t rp1  = TMath::Power(r, fSSPars[1]) ;\r
 453   Double_t rp5  = TMath::Power(r, fSSPars[5]) ;\r
 454   Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[6] * rp5) ) ) ;\r
 455   return shape ;\r
 456 }\r
 457 \r
 458 //____________________________________________________________________________\r
 459 void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r
 460                                                  Double_t & fret,\r
 461                                                  Double_t * x, Int_t iflag)\r
 462 {\r
 463   // Calculates the Chi square for the cluster unfolding minimization\r
 464   // Number of parameters, Gradient, Chi squared, parameters, what to do\r
 465   \r
 466   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
 467   if(toMinuit){\r
 468     AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;\r
 469     TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;\r
 470     // A bit buggy way to get an access to the geometry\r
 471     // To be revised!\r
 472     AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
 473     \r
 474     if(recPoint && digits && geom){\r
 475       \r
 476       Int_t * digitsList     = recPoint->GetDigitsList() ;\r
 477       \r
 478       Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
 479       \r
 480       Float_t * energiesList = recPoint->GetEnergiesList() ;\r
 481       \r
 482       fret = 0. ;\r
 483       Int_t iparam = 0 ;\r
 484       \r
 485       if(iflag == 2)\r
 486         for(iparam = 0 ; iparam < nPar ; iparam++)\r
 487           Grad[iparam] = 0 ; // Will evaluate gradient\r
 488       \r
 489       Double_t efit = 0. ;\r
 490       \r
 491       AliEMCALDigit * digit ;\r
 492       Int_t iDigit ;\r
 493       \r
 494       Int_t iSupMod =  0 ;\r
 495       Int_t iTower  =  0 ;\r
 496       Int_t iIphi   =  0 ;\r
 497       Int_t iIeta   =  0 ;\r
 498       Int_t iphi    =  0 ;//x direction\r
 499       Int_t ieta    =  0 ;//z direstion\r
 500       \r
 501       \r
 502       for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
 503         if(energiesList[iDigit]==0) continue;\r
 504         \r
 505         digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
 506         \r
 507         if(digit){\r
 508         geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
 509         geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
 510                                           iIphi, iIeta,iphi,ieta);\r
 511         EvalParsPhiDependence(digit->GetId(),geom);\r
 512         \r
 513         if(iflag == 2){  // calculate gradient\r
 514           Int_t iParam = 0 ;\r
 515           efit = 0. ;\r
 516           while(iParam < nPar ){\r
 517             Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
 518             iParam++ ;\r
 519             Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
 520             iParam++ ;\r
 521             efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
 522             iParam++ ;\r
 523           }\r
 524           \r
 525           Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
 526           iParam = 0 ;\r
 527           while(iParam < nPar ){\r
 528             Double_t xpar = x[iParam] ;\r
 529             Double_t zpar = x[iParam+1] ;\r
 530             Double_t epar = x[iParam+2] ;\r
 531             \r
 532             Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r
 533             Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
 534             Double_t rp1  = TMath::Power(dr, fSSPars[1]) ;\r
 535             Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;\r
 536             \r
 537             Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
 538             (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
 539              (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
 540               fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\r
 541             \r
 542             //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )\r
 543             //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
 544             \r
 545             Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x\r
 546             iParam++ ;\r
 547             Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z\r
 548             iParam++ ;\r
 549             Grad[iParam] += shape ;                                  // Derivative over energy\r
 550             iParam++ ;\r
 551           }\r
 552         }\r
 553         efit = 0;\r
 554         iparam = 0 ;\r
 555         \r
 556         while(iparam < nPar ){\r
 557           Double_t xpar = x[iparam] ;\r
 558           Double_t zpar = x[iparam+1] ;\r
 559           Double_t epar = x[iparam+2] ;\r
 560           iparam += 3 ;\r
 561           efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
 562         }\r
 563         \r
 564         fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
 565         // Here we assume, that sigma = sqrt(E) \r
 566         } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n");\r
 567       } // digit loop\r
 568     } // recpoint, digits and geom not NULL\r
 569   }// List is not NULL\r
 570   \r
 571 }\r
 572 \r
 573 \r
 574 //____________________________________________________________________________\r
 575 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r
 576   for(UInt_t i=0;i<7;++i)\r
 577     fSSPars[i]=pars[i];\r
 578   if(pars[2]==0. && pars[3]==0.) fSSPars[2]=1.;//to avoid dividing by 0\r
 579 }\r
 580 \r
 581 //____________________________________________________________________________\r
 582 void AliEMCALUnfolding::SetPar5(Double_t *pars){\r
 583   for(UInt_t i=0;i<3;++i)\r
 584     fPar5[i]=pars[i];\r
 585 }\r
 586 \r
 587 //____________________________________________________________________________\r
 588 void AliEMCALUnfolding::SetPar6(Double_t *pars){\r
 589   for(UInt_t i=0;i<3;++i)\r
 590     fPar6[i]=pars[i];\r
 591 }\r
 592 \r
 593 //____________________________________________________________________________\r
 594 void AliEMCALUnfolding::EvalPar5(Double_t phi){\r
 595   //\r
 596   //Evaluate the 5th parameter of the shower shape function\r
 597   //phi in degrees range (-10,10)\r
 598   //\r
 599   //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r
 600   fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[2];\r
 601 }\r
 602 \r
 603 //____________________________________________________________________________\r
 604 void AliEMCALUnfolding::EvalPar6(Double_t phi){\r
 605   //\r
 606   //Evaluate the 6th parameter of the shower shape function\r
 607   //phi in degrees range (-10,10)\r
 608   //\r
 609   //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r
 610   fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r
 611 }\r
 612 \r
 613 //____________________________________________________________________________\r
 614 void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom){\r
 615   //\r
 616   // calculate params p5 and p6 depending on the phi angle in global coordinate\r
 617   // for the cell with given absId index\r
 618   //\r
 619   Double_t etaGlob = 0.;//eta in global c.s. - unused\r
 620   Double_t phiGlob = 0.;//phi in global c.s. in radians\r
 621   geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);\r
 622   phiGlob*=180./TMath::Pi();\r
 623   phiGlob-=90.;\r
 624   phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);\r
 625 \r
 626   EvalPar5(phiGlob);\r
 627   EvalPar6(phiGlob);\r
 628 }\r
 629 \r