EMCAL/AliEMCALUnfolding.cxx

   1 /**************************************************************************\r\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r\r
   3  *                                                                        *\r\r
   4  * Author: The ALICE Off-line Project.                                    *\r\r
   5  * Contributors are mentioned in the code where appropriate.              *\r\r
   6  *                                                                        *\r\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r\r
   9  * without fee, provided that the above copyright notice appears in all   *\r\r
  10  * copies and that both the copyright notice and this permission notice   *\r\r
  11  * appear in the supporting documentation. The authors make no claims     *\r\r
  12  * about the suitability of this software for any purpose. It is          *\r\r
  13  * provided "as is" without express or implied warranty.                  *\r\r
  14  **************************************************************************/\r\r
  15 \r\r
  16 //_________________________________________________________________________\r\r
  17 //  Base class for the cluster unfolding algorithm \r\r
  18 //*-- Author: Adam Matyja (SUBATECH)\r\r
  19 //  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis\r\r
  20 //-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1\r\r
  21 //-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)\r\r
  22 //\r\r
  23 //  unfolds the clusters having several local maxima. \r\r
  24 //////////////////////////////////////////////////////////////////////////////\r\r
  25 \r\r
  26 // --- ROOT system ---\r\r
  27 #include "TClonesArray.h"\r\r
  28 #include <TMath.h> \r\r
  29 #include <TMinuit.h>\r\r
  30 \r\r
  31 // --- Standard library ---\r\r
  32 #include <cassert>\r\r
  33 \r\r
  34 // --- AliRoot header files ---\r\r
  35 #include "AliEMCALUnfolding.h"\r\r
  36 #include "AliEMCALGeometry.h"\r\r
  37 #include "AliRunLoader.h"\r\r
  38 #include "AliRun.h"\r\r
  39 #include "AliEMCAL.h"\r\r
  40 #include "AliEMCALRecParam.h"\r\r
  41 #include "AliEMCALRecPoint.h"\r\r
  42 #include "AliEMCALDigit.h"\r\r
  43 #include "AliEMCALReconstructor.h"\r\r
  44 \r\r
  45 #include "AliLog.h"\r\r
  46 #include "AliCDBManager.h"\r\r
  47 class AliCDBStorage;\r\r
  48 #include "AliCDBEntry.h"\r\r
  49 \r\r
  50 Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r\r
  51 Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};\r\r
  52 Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};\r\r
  53 \r\r
  54 ClassImp(AliEMCALUnfolding)\r\r
  55   \r\r
  56 //____________________________________________________________________________\r\r
  57 AliEMCALUnfolding::AliEMCALUnfolding():\r\r
  58   fNumberOfECAClusters(0),\r\r
  59   fECALocMaxCut(0),\r\r
  60   fThreshold(0.01),//10 MeV\r\r
  61   fGeom(NULL),\r\r
  62   fRecPoints(NULL),\r\r
  63   fDigitsArr(NULL)\r\r
  64 {\r\r
  65   // ctor with the indication of the file where header Tree and digits Tree are stored\r\r
  66   Init() ;\r\r
  67 }\r\r
  68 \r\r
  69 //____________________________________________________________________________\r\r
  70 AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):\r\r
  71   fNumberOfECAClusters(0),\r\r
  72   fECALocMaxCut(0),\r\r
  73   fThreshold(0.01),//10 MeV\r\r
  74   fGeom(geometry),\r\r
  75   fRecPoints(NULL),\r\r
  76   fDigitsArr(NULL)\r\r
  77 {\r\r
  78   // ctor with the indication of the file where header Tree and digits Tree are stored\r\r
  79   // use this contructor to avoid usage of Init() which uses runloader\r\r
  80   // change needed by HLT - MP\r\r
  81   if (!fGeom)\r\r
  82   {\r\r
  83     AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r\r
  84   }\r\r
  85 \r\r
  86 }\r\r
  87 \r\r
  88 //____________________________________________________________________________\r\r
  89 AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):\r\r
  90   fNumberOfECAClusters(0),\r\r
  91   fECALocMaxCut(ECALocMaxCut),\r\r
  92   fThreshold(0.01),//10 MeV\r\r
  93   fGeom(geometry),\r\r
  94   fRecPoints(NULL),\r\r
  95   fDigitsArr(NULL)\r\r
  96 {\r\r
  97   // ctor with the indication of the file where header Tree and digits Tree are stored\r\r
  98   // use this contructor to avoid usage of Init() which uses runloader\r\r
  99   // change needed by HLT - MP\r\r
 100   if (!fGeom)\r\r
 101   {\r\r
 102     AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r\r
 103   }\r\r
 104   Int_t i=0;\r\r
 105   for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];\r\r
 106   for (i = 0; i < 3; i++) {\r\r
 107     fgPar5[i] = Par5[i];\r\r
 108     fgPar6[i] = Par6[i];\r\r
 109   }\r\r
 110 \r\r
 111 }\r\r
 112 \r\r
 113 //____________________________________________________________________________\r\r
 114 void AliEMCALUnfolding::Init()\r\r
 115 {\r\r
 116   // Make all memory allocations which can not be done in default constructor.\r\r
 117   // Attach the Clusterizer task to the list of EMCAL tasks\r\r
 118 \r\r
 119   AliRunLoader *rl = AliRunLoader::Instance();\r\r
 120   if (rl && rl->GetAliRun()){\r\r
 121     AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));\r\r
 122     if(emcal)fGeom = emcal->GetGeometry();\r\r
 123   }\r\r
 124   \r\r
 125   if(!fGeom)\r\r
 126     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());\r\r
 127   \r\r
 128   AliDebug(1,Form("geom %p",fGeom));\r\r
 129   \r\r
 130   if(!gMinuit) \r\r
 131     //    gMinuit = new TMinuit(100) ;//the same is in FindFitV2\r\r
 132     gMinuit = new TMinuit(30) ;//the same is in FindFitV2\r\r
 133   \r\r
 134 }\r\r
 135 \r\r
 136 //____________________________________________________________________________\r\r
 137   AliEMCALUnfolding::~AliEMCALUnfolding()\r\r
 138 {\r\r
 139   // dtor\r\r
 140 }\r\r
 141 \r\r
 142 //____________________________________________________________________________\r\r
 143 void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)\r\r
 144 {\r\r
 145   //\r\r
 146   //Set input for unfolding purposes\r\r
 147   //\r\r
 148   SetNumberOfECAClusters(numberOfECAClusters);\r\r
 149   SetRecPoints(recPoints);\r\r
 150   SetDigitsArr(digitsArr);\r\r
 151 }\r\r
 152 \r\r
 153 //____________________________________________________________________________\r\r
 154 void AliEMCALUnfolding::MakeUnfolding()\r\r
 155 {\r\r
 156   // Unfolds clusters using the shape of an ElectroMagnetic shower\r\r
 157   // Performs unfolding of all clusters\r\r
 158   \r\r
 159   if(fNumberOfECAClusters > 0){\r\r
 160     if (fGeom==0)\r\r
 161       AliFatal("Did not get geometry from EMCALLoader") ;\r\r
 162     //Int_t nModulesToUnfold = fGeom->GetNCells();\r\r
 163     \r\r
 164     Int_t numberOfClustersToUnfold=fNumberOfECAClusters;\r\r
 165     //we unfold only clusters present in the array untill now\r\r
 166     //fNumberOfECAClusters may change due to unfilded clusters\r\r
 167     //so 0 to numberOfClustersToUnfold-1: clusters before unfolding\r\r
 168     //numberOfClustersToUnfold to the end: new clusters from unfolding\r\r
 169     //of course numberOfClustersToUnfold also is decreased but we don't loop over clusters added in UF \r\r
 170     Int_t index ;\r\r
 171     for(index = 0 ; index < numberOfClustersToUnfold ; index++){\r\r
 172       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r\r
 173       if(recPoint){\r\r
 174         Int_t nMultipl = recPoint->GetMultiplicity() ;\r\r
 175         AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r\r
 176         Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r\r
 177         Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r\r
 178         if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0\r\r
 179           if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r\r
 180             //if unfolding correct remove old recPoint\r\r
 181             fRecPoints->Remove(recPoint);\r\r
 182             fRecPoints->Compress() ;//is it really needed\r\r
 183             index-- ;\r\r
 184             fNumberOfECAClusters-- ;\r\r
 185             numberOfClustersToUnfold--;\r\r
 186           }\r\r
 187         } else{\r\r
 188           recPoint->SetNExMax(1) ; //Only one local maximum\r\r
 189         }\r\r
 190         \r\r
 191         delete[] maxAt ;\r\r
 192         delete[] maxAtEnergy ;\r\r
 193       } else {\r\r
 194         //AliError("RecPoint NULL"); //end of check if recPoint exist\r\r
 195         Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ;\r\r
 196       }\r\r
 197     } // rec point loop\r\r
 198   }//end of check fNumberOfECAClusters\r\r
 199   // End of Unfolding of clusters\r\r
 200 }\r\r
 201 \r\r
 202 //____________________________________________________________________________\r\r
 203 Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, \r\r
 204                                           Int_t nMax, \r\r
 205                                           AliEMCALDigit ** maxAt, \r\r
 206                                           Float_t * maxAtEnergy)\r\r
 207 {\r\r
 208   // Extended to whole EMCAL \r\r
 209   \r\r
 210   //**************************** part 1 *******************************************\r\r
 211   // Performs the unfolding of a cluster with nMax overlapping showers \r\r
 212   \r\r
 213   //cout<<"unfolding check here part 1"<<endl;\r\r
 214   //printf("Original cluster E %f\n",iniTower->GetEnergy());\r\r
 215   \r\r
 216   Int_t nPar = 3 * nMax ;\r\r
 217   Float_t * fitparameters = new Float_t[nPar] ;\r\r
 218   \r\r
 219   if (fGeom==0)\r\r
 220     AliFatal("Did not get geometry from EMCALLoader") ;\r\r
 221   \r\r
 222   Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r\r
 223   if( !rv ) \r\r
 224   {\r\r
 225     // Fit failed, return (and remove cluster? - why? I leave the cluster)\r\r
 226     iniTower->SetNExMax(-1) ;\r\r
 227     delete[] fitparameters ;\r\r
 228     return kFALSE;\r\r
 229   }\r\r
 230   \r\r
 231   //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold\r\r
 232   if(nMax==2){\r\r
 233     if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){\r\r
 234       AliDebug(1,"One of fitted energy below threshold");\r\r
 235       iniTower->SetNExMax(1) ;\r\r
 236       delete[] fitparameters ;\r\r
 237       return kFALSE;\r\r
 238     }\r\r
 239   }\r\r
 240 \r\r
 241   //**************************** part 2 *******************************************\r\r
 242   // create unfolded rec points and fill them with new energy lists\r\r
 243   // First calculate energy deposited in each sell in accordance with\r\r
 244   // fit (without fluctuations): efit[]\r\r
 245   // and later correct this number in acordance with actual energy\r\r
 246   // deposition\r\r
 247   \r\r
 248   //  cout<<"unfolding check here part 2"<<endl;\r\r
 249   Int_t nDigits = iniTower->GetMultiplicity() ;\r\r
 250   Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r\r
 251   Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r\r
 252   \r\r
 253   AliEMCALDigit * digit = 0 ;\r\r
 254   Int_t * digitsList = iniTower->GetDigitsList() ;\r\r
 255   \r\r
 256   Int_t iSupMod =  0 ;\r\r
 257   Int_t iTower  =  0 ;\r\r
 258   Int_t iIphi   =  0 ;\r\r
 259   Int_t iIeta   =  0 ;\r\r
 260   Int_t iphi    =  0 ;//x direction\r\r
 261   Int_t ieta    =  0 ;//z direstion\r\r
 262   \r\r
 263   Int_t iparam = 0 ;\r\r
 264   Int_t iDigit = 0 ;\r\r
 265   \r\r
 266   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r\r
 267   {\r\r
 268     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r\r
 269     if(digit)\r\r
 270     {\r\r
 271       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 272       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 273                                          iIphi, iIeta,iphi,ieta);\r\r
 274       EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
 275       \r\r
 276       efit[iDigit] = 0.;\r\r
 277       iparam = 0;\r\r
 278       while(iparam < nPar )\r\r
 279       {\r\r
 280         xpar = fitparameters[iparam] ;\r\r
 281         zpar = fitparameters[iparam+1] ;\r\r
 282         epar = fitparameters[iparam+2] ;\r\r
 283 \r\r
 284         efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
 285         iparam += 3 ;\r\r
 286       }\r\r
 287 \r\r
 288     } else AliDebug(1,"Digit NULL part 2!");\r\r
 289     \r\r
 290   }//digit loop\r\r
 291   \r\r
 292   //**************************** part 3 *******************************************\r\r
 293   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r\r
 294   // so that energy deposited in each cell is distributed between new clusters proportionally\r\r
 295   // to its contribution to efit\r\r
 296   \r\r
 297   Float_t * energiesList = iniTower->GetEnergiesList() ;\r\r
 298   Float_t ratio = 0. ;//0 -> 0. changed\r\r
 299   Float_t eDigit = 0. ;\r\r
 300   Int_t nSplittedClusters=(Int_t)nPar/3;\r\r
 301   \r\r
 302   Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters];\r\r
 303   //above - temporary table with energies after unfolding.\r\r
 304   //the order is following: \r\r
 305   //first cluster <first cell - last cell>, \r\r
 306   //second cluster <first cell - last cell>, etc.\r\r
 307   \r\r
 308   //**************************** sub-part 3.1 *************************************\r\r
 309   //here we check if energy of the cell in the cluster after unfolding is above threshold. \r\r
 310   //If not the energy from a given cell in the cluster is divided in correct proportions \r\r
 311   //in accordance to the other clusters and added to them and set to 0.\r\r
 312   \r\r
 313   //  cout<<"unfolding check here part 3.1"<<endl;\r\r
 314 \r\r
 315   iparam = 0 ;\r\r
 316   while(iparam < nPar )\r\r
 317   {\r\r
 318     xpar = fitparameters[iparam] ;\r\r
 319     zpar = fitparameters[iparam+1] ;\r\r
 320     epar = fitparameters[iparam+2] ;\r\r
 321     \r\r
 322     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r\r
 323     {\r\r
 324       digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
 325       if(digit)\r\r
 326       {\r\r
 327         fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 328         fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 329                                            iIphi, iIeta,iphi,ieta);\r\r
 330         \r\r
 331         EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
 332        \r\r
 333         if(efit[iDigit]==0) \r\r
 334         {//just for sure\r\r
 335           correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here\r\r
 336           continue;\r\r
 337         }\r\r
 338         \r\r
 339         ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r\r
 340         eDigit = energiesList[iDigit] * ratio ;\r\r
 341         \r\r
 342         //add energy to temporary matrix\r\r
 343         correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;\r\r
 344         \r\r
 345       } else AliDebug(1,"NULL digit part 3");\r\r
 346     }//digit loop \r\r
 347     iparam += 3 ;\r\r
 348   }//while\r\r
 349   \r\r
 350   //**************************** sub-part 3.2 *************************************\r\r
 351   //here we correct energy for each cell and cluster\r\r
 352   //  cout<<"unfolding check here part 3.2"<<endl;\r\r
 353 \r\r
 354   Float_t maximumEne=0.;//0 -> 0. changed\r\r
 355   Int_t maximumIndex=0;\r\r
 356   Bool_t isAnyBelowThreshold=kFALSE;\r\r
 357   //  Float_t Threshold=0.01;\r\r
 358   Float_t * energyFraction = new Float_t[nSplittedClusters];\r\r
 359   Int_t iparam2 = 0 ;\r\r
 360   for(iDigit = 0 ; iDigit < nDigits ; iDigit++)\r\r
 361   {\r\r
 362     isAnyBelowThreshold=kFALSE;\r\r
 363     maximumEne=0.;//0 -> 0. changed\r\r
 364     for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
 365     {\r\r
 366 \r\r
 367       if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;\r\r
 368       if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne) \r\r
 369       {\r\r
 370         maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
 371         maximumIndex = iparam;\r\r
 372       }\r\r
 373     }//end of loop over clusters after unfolding\r\r
 374     \r\r
 375     if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold \r\r
 376 \r\r
 377 \r\r
 378     if(maximumEne < fThreshold) \r\r
 379     {//add all cluster cells and put energy into max index, other set to 0\r\r
 380       maximumEne=0.;\r\r
 381       for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
 382       {\r\r
 383         maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
 384         correctedEnergyList[iparam/3*nDigits+iDigit]=0;\r\r
 385       }\r\r
 386       correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;\r\r
 387       continue;\r\r
 388     }//end if\r\r
 389     \r\r
 390     //divide energy of cell below threshold in the correct proportion and add to other cells\r\r
 391     maximumEne=0.;//not used any more so use it for the energy sum \r\r
 392     for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
 393     {//calculate energy sum\r\r
 394       if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r\r
 395       else \r\r
 396       {\r\r
 397         energyFraction[iparam/3]=1.;\r\r
 398         maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
 399       }\r\r
 400     }//end of loop over clusters after unfolding\r\r
 401     if(maximumEne>0.)\r\r
 402     {\r\r
 403       for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r\r
 404         energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;\r\r
 405       }\r\r
 406       \r\r
 407       for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
 408       {//add energy from cells below threshold to others\r\r
 409         if(energyFraction[iparam/3]>0.) continue;\r\r
 410         else\r\r
 411         {\r\r
 412           for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)\r\r
 413           {\r\r
 414             correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * \r\r
 415                                                       correctedEnergyList[iparam/3*nDigits+iDigit]) ;\r\r
 416           }//inner loop\r\r
 417           correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;\r\r
 418         }\r\r
 419       }\r\r
 420     }\r\r
 421     else\r\r
 422     {\r\r
 423       //digit energy to be set to 0\r\r
 424       for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
 425       {\r\r
 426         correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;\r\r
 427       }\r\r
 428     }//correction for: is energy>0\r\r
 429     \r\r
 430   }//end of loop over digits\r\r
 431   delete[] energyFraction;\r\r
 432   \r\r
 433   //**************************** sub-part 3.3 *************************************\r\r
 434   //here we add digits to recpoints with corrected energy\r\r
 435   //  cout<<"unfolding check here part 3.3"<<endl;\r\r
 436 \r\r
 437   iparam = 0 ;\r\r
 438   while(iparam < nPar )\r\r
 439   {\r\r
 440     AliEMCALRecPoint * recPoint = 0 ;\r\r
 441     \r\r
 442     if(fNumberOfECAClusters >= fRecPoints->GetSize())\r\r
 443       fRecPoints->Expand(2*fNumberOfECAClusters) ;\r\r
 444     \r\r
 445     //add recpoint\r\r
 446     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;//fNumberOfECAClusters-1 is old cluster before unfolding\r\r
 447     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r\r
 448     \r\r
 449     if(recPoint){//recPoint present -> good\r\r
 450       recPoint->SetNExMax(nSplittedClusters) ;\r\r
 451       \r\r
 452       for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) {\r\r
 453         digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
 454         if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. ){\r\r
 455           //if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold);\r\r
 456           recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case\r\r
 457         } else {\r\r
 458           AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));\r\r
 459         }\r\r
 460       }//digit loop\r\r
 461 \r\r
 462       if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list\r\r
 463         delete (*fRecPoints)[fNumberOfECAClusters];\r\r
 464         fRecPoints->RemoveAt(fNumberOfECAClusters);\r\r
 465         //fNumberOfECAClusters--; //it should be like that we don't want to decrease cluster number\r\r
 466         nSplittedClusters--;\r\r
 467       } else {//recPoint exists and has digits associated -> very good increase number of clusters \r\r
 468         fNumberOfECAClusters++ ; \r\r
 469       }\r\r
 470       \r\r
 471     } else {//recPoint empty -> remove from list\r\r
 472       AliError("NULL RecPoint");\r\r
 473       //protection from recpoint with no digits\r\r
 474 //      if(recPoint->GetMultiplicity()==0)\r\r
 475 //      {\r\r
 476       delete (*fRecPoints)[fNumberOfECAClusters];\r\r
 477       fRecPoints->RemoveAt(fNumberOfECAClusters);\r\r
 478       //  fNumberOfECAClusters--;//it should be like that we don't want to decrease cluster number\r\r
 479       nSplittedClusters--;\r\r
 480 //      }\r\r
 481     \r\r
 482     }\r\r
 483 \r\r
 484     iparam += 3 ;\r\r
 485   }//while\r\r
 486   \r\r
 487   delete[] fitparameters ;\r\r
 488   delete[] efit ;\r\r
 489   delete[] correctedEnergyList ;\r\r
 490   \r\r
 491   //  cout<<"end of unfolding check part 3.3"<<endl;\r\r
 492   return kTRUE;\r\r
 493 }\r\r
 494 \r\r
 495 \r\r
 496 //____________________________________________________________________________\r\r
 497 Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, \r\r
 498                                           Int_t nMax, \r\r
 499                                           AliEMCALDigit ** maxAt, \r\r
 500                                           Float_t * maxAtEnergy)\r\r
 501 {\r\r
 502   // Extended to whole EMCAL \r\r
 503   // Performs the unfolding of a cluster with nMax overlapping showers \r\r
 504   \r\r
 505   Int_t nPar = 3 * nMax ;\r\r
 506   Float_t * fitparameters = new Float_t[nPar] ;\r\r
 507   \r\r
 508   if (fGeom==0)\r\r
 509     AliFatal("Did not get geometry from EMCALLoader") ;\r\r
 510   \r\r
 511   Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r\r
 512   if( !rv ) {\r\r
 513     // Fit failed, return (and remove cluster? - why? I leave the cluster)\r\r
 514     iniTower->SetNExMax(-1) ;\r\r
 515     delete[] fitparameters ;\r\r
 516     return kFALSE;\r\r
 517   }\r\r
 518   \r\r
 519   // create unfolded rec points and fill them with new energy lists\r\r
 520   // First calculate energy deposited in each sell in accordance with\r\r
 521   // fit (without fluctuations): efit[]\r\r
 522   // and later correct this number in acordance with actual energy\r\r
 523   // deposition\r\r
 524   \r\r
 525   Int_t nDigits = iniTower->GetMultiplicity() ;\r\r
 526   Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r\r
 527   Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r\r
 528   \r\r
 529   AliEMCALDigit * digit = 0 ;\r\r
 530   Int_t * digitsList = iniTower->GetDigitsList() ;\r\r
 531   \r\r
 532   Int_t iSupMod =  0 ;\r\r
 533   Int_t iTower  =  0 ;\r\r
 534   Int_t iIphi   =  0 ;\r\r
 535   Int_t iIeta   =  0 ;\r\r
 536   Int_t iphi    =  0 ;//x direction\r\r
 537   Int_t ieta    =  0 ;//z direstion\r\r
 538   \r\r
 539   Int_t iparam = 0 ;\r\r
 540   Int_t iDigit = 0 ;\r\r
 541   \r\r
 542   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
 543     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r\r
 544     if(digit){\r\r
 545       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 546       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 547                                          iIphi, iIeta,iphi,ieta);\r\r
 548       EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
 549       \r\r
 550       efit[iDigit] = 0.;\r\r
 551       iparam = 0;\r\r
 552       while(iparam < nPar ){\r\r
 553         xpar = fitparameters[iparam] ;\r\r
 554         zpar = fitparameters[iparam+1] ;\r\r
 555         epar = fitparameters[iparam+2] ;\r\r
 556         iparam += 3 ;\r\r
 557         \r\r
 558         efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
 559       }\r\r
 560     } else AliError("Digit NULL!");\r\r
 561     \r\r
 562   }//digit loop\r\r
 563   \r\r
 564   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r\r
 565   // so that energy deposited in each cell is distributed between new clusters proportionally\r\r
 566   // to its contribution to efit\r\r
 567   \r\r
 568   Float_t * energiesList = iniTower->GetEnergiesList() ;\r\r
 569   Float_t ratio = 0 ;\r\r
 570   \r\r
 571   iparam = 0 ;\r\r
 572   while(iparam < nPar ){\r\r
 573     xpar = fitparameters[iparam] ;\r\r
 574     zpar = fitparameters[iparam+1] ;\r\r
 575     epar = fitparameters[iparam+2] ;\r\r
 576     iparam += 3 ;\r\r
 577     \r\r
 578     AliEMCALRecPoint * recPoint = 0 ;\r\r
 579     \r\r
 580     if(fNumberOfECAClusters >= fRecPoints->GetSize())\r\r
 581       fRecPoints->Expand(2*fNumberOfECAClusters) ;\r\r
 582     \r\r
 583     //add recpoint\r\r
 584     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r\r
 585     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r\r
 586     \r\r
 587     if(recPoint){\r\r
 588       \r\r
 589       fNumberOfECAClusters++ ;\r\r
 590       recPoint->SetNExMax((Int_t)nPar/3) ;\r\r
 591       \r\r
 592       Float_t eDigit = 0. ;\r\r
 593       for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
 594         digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
 595         if(digit){\r\r
 596           fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 597           fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 598                                              iIphi, iIeta,iphi,ieta);\r\r
 599           EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
 600           if(efit[iDigit]==0) continue;//just for sure\r\r
 601           ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r\r
 602           eDigit = energiesList[iDigit] * ratio ;\r\r
 603           recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r\r
 604         } else AliError("NULL digit");\r\r
 605       }//digit loop \r\r
 606     } else AliError("NULL RecPoint");\r\r
 607   }//while\r\r
 608   \r\r
 609   delete[] fitparameters ;\r\r
 610   delete[] efit ;\r\r
 611   \r\r
 612   return kTRUE;\r\r
 613 }\r\r
 614 \r\r
 615 \r\r
 616 //____________________________________________________________________________\r\r
 617 Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, \r\r
 618                                         const Float_t* maxAtEnergy,\r\r
 619                                         Int_t nPar, Float_t * fitparameters) const\r\r
 620 {\r\r
 621   // Calls TMinuit to fit the energy distribution of a cluster with several maxima\r\r
 622   // The initial values for fitting procedure are set equal to the\r\r
 623   // positions of local maxima.       \r\r
 624   // Cluster will be fitted as a superposition of nPar/3\r\r
 625   // electromagnetic showers\r\r
 626 \r\r
 627   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r\r
 628         \r\r
 629   if(!gMinuit){\r\r
 630     //    gMinuit = new TMinuit(100) ;//max 100 parameters\r\r
 631     if(nPar<30) gMinuit = new TMinuit(30);\r\r
 632     else gMinuit = new TMinuit(nPar) ;//max nPar parameters\r\r
 633     //\r\r
 634   } else {\r\r
 635     if(gMinuit->fMaxpar < nPar) {\r\r
 636       delete gMinuit;\r\r
 637       gMinuit = new TMinuit(nPar);\r\r
 638     }\r\r
 639   }\r\r
 640 \r\r
 641   gMinuit->mncler();                     // Reset Minuit's list of paramters\r\r
 642   gMinuit->SetPrintLevel(-1) ;           // No Printout\r\r
 643   gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;\r\r
 644   // To set the address of the minimization function\r\r
 645   TList * toMinuit = new TList();\r\r
 646   toMinuit->AddAt(recPoint,0) ;\r\r
 647   toMinuit->AddAt(fDigitsArr,1) ;\r\r
 648   toMinuit->AddAt(fGeom,2) ;\r\r
 649 \r\r
 650   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare\r\r
 651 \r\r
 652   // filling initial values for fit parameters\r\r
 653   AliEMCALDigit * digit ;\r\r
 654 \r\r
 655   Int_t ierflg  = 0;\r\r
 656   Int_t index   = 0 ;\r\r
 657   Int_t nDigits = (Int_t) nPar / 3 ;\r\r
 658 \r\r
 659   Int_t iDigit ;\r\r
 660 \r\r
 661   Int_t iSupMod =  0 ;\r\r
 662   Int_t iTower  =  0 ;\r\r
 663   Int_t iIphi   =  0 ;\r\r
 664   Int_t iIeta   =  0 ;\r\r
 665   Int_t iphi    =  0 ;//x direction\r\r
 666   Int_t ieta    =  0 ;//z direstion\r\r
 667 \r\r
 668   for(iDigit = 0; iDigit < nDigits; iDigit++){\r\r
 669     digit = maxAt[iDigit];\r\r
 670     if(digit==0) AliError("energy of digit = 0!");\r\r
 671     fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 672     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 673                                        iIphi, iIeta,iphi,ieta);\r\r
 674 \r\r
 675     Float_t energy = maxAtEnergy[iDigit] ;\r\r
 676 \r\r
 677     //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original\r\r
 678     gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;\r\r
 679     index++ ;\r\r
 680     if(ierflg != 0){\r\r
 681       Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: x=%d, param.id=%d, nMaxima=%d",iphi,index-1,nPar/3 ) ;\r\r
 682       toMinuit->Clear();\r\r
 683       delete toMinuit ;\r\r
 684       return kFALSE;\r\r
 685     }\r\r
 686     //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original\r\r
 687     gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;\r\r
 688     index++ ;\r\r
 689     if(ierflg != 0){\r\r
 690       Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: z=%d, param.id=%d, nMaxima=%d", ieta, index-1,nPar/3) ;\r\r
 691       toMinuit->Clear();\r\r
 692       delete toMinuit ;\r\r
 693       return kFALSE;\r\r
 694     }\r\r
 695     //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original\r\r
 696     gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r\r
 697     index++ ;\r\r
 698     if(ierflg != 0){\r\r
 699       Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: energy = %f, param.id=%d, nMaxima=%d", energy, index-1, nPar/3) ;\r\r
 700       toMinuit->Clear();\r\r
 701       delete toMinuit ;\r\r
 702       return kFALSE;\r\r
 703     }\r\r
 704   }\r\r
 705 \r\r
 706   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; \r\r
 707                       // The number of function call slightly depends on it.\r\r
 708   //  Double_t p1 = 1.0 ;// par to gradient \r\r
 709   Double_t p2 = 0.0 ;\r\r
 710   //  Double_t p3 = 3.0 ;\r\r
 711   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls\r\r
 712   //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient\r\r
 713   gMinuit->SetMaxIterations(5);//was 5\r\r
 714   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings\r\r
 715   //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts\r\r
 716 \r\r
 717   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize\r\r
 718   //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize\r\r
 719   if(ierflg == 4){  // Minimum not found\r\r
 720     AliDebug(1,"EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;\r\r
 721     toMinuit->Clear();\r\r
 722     delete toMinuit ;\r\r
 723     return kFALSE ;\r\r
 724   }\r\r
 725   for(index = 0; index < nPar; index++){\r\r
 726     Double_t err = 0. ;\r\r
 727     Double_t val = 0. ;\r\r
 728     gMinuit->GetParameter(index, val, err) ;    // Returns value and error ofOA parameter index\r\r
 729     fitparameters[index] = val ;\r\r
 730   }\r\r
 731 \r\r
 732   toMinuit->Clear();\r\r
 733   delete toMinuit ;\r\r
 734 \r\r
 735   if(gMinuit->fMaxpar>30) delete gMinuit;\r\r
 736 \r\r
 737   return kTRUE;\r\r
 738 \r\r
 739 }\r\r
 740 \r\r
 741 //____________________________________________________________________________\r\r
 742 Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)\r\r
 743 { \r\r
 744   // extended to whole EMCAL \r\r
 745   // Shape of the shower\r\r
 746   // If you change this function, change also the gradient evaluation in ChiSquare()\r\r
 747 \r\r
 748   Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);\r\r
 749   Double_t rp1  = TMath::Power(r, fgSSPars[1]) ;\r\r
 750   Double_t rp5  = TMath::Power(r, fgSSPars[5]) ;\r\r
 751   Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;\r\r
 752   return shape ;\r\r
 753 }\r\r
 754 \r\r
 755 //____________________________________________________________________________\r\r
 756 void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r\r
 757                                              Double_t & fret,\r\r
 758                                              Double_t * x, Int_t iflag)\r\r
 759 {\r\r
 760   // Calculates the Chi square for the cluster unfolding minimization\r\r
 761   // Number of parameters, Gradient, Chi squared, parameters, what to do\r\r
 762   \r\r
 763   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r\r
 764   if(toMinuit){\r\r
 765     AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;\r\r
 766     TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;\r\r
 767     // A bit buggy way to get an access to the geometry\r\r
 768     // To be revised!\r\r
 769     AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r\r
 770     \r\r
 771     if(recPoint && digits && geom){\r\r
 772       \r\r
 773       Int_t * digitsList     = recPoint->GetDigitsList() ;\r\r
 774       \r\r
 775       Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r\r
 776       \r\r
 777       Float_t * energiesList = recPoint->GetEnergiesList() ;\r\r
 778       \r\r
 779       fret = 0. ;\r\r
 780       Int_t iparam = 0 ;\r\r
 781       \r\r
 782       if(iflag == 2)\r\r
 783         for(iparam = 0 ; iparam < nPar ; iparam++)\r\r
 784           Grad[iparam] = 0 ; // Will evaluate gradient\r\r
 785       \r\r
 786       Double_t efit = 0. ;\r\r
 787       \r\r
 788       AliEMCALDigit * digit ;\r\r
 789       Int_t iDigit ;\r\r
 790       \r\r
 791       Int_t iSupMod =  0 ;\r\r
 792       Int_t iTower  =  0 ;\r\r
 793       Int_t iIphi   =  0 ;\r\r
 794       Int_t iIeta   =  0 ;\r\r
 795       Int_t iphi    =  0 ;//x direction\r\r
 796       Int_t ieta    =  0 ;//z direstion\r\r
 797       \r\r
 798       \r\r
 799       for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r\r
 800         if(energiesList[iDigit]==0) continue;\r\r
 801         \r\r
 802         digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r\r
 803         \r\r
 804         if(digit){\r\r
 805         geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
 806         geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
 807                                           iIphi, iIeta,iphi,ieta);\r\r
 808         EvalParsPhiDependence(digit->GetId(),geom);\r\r
 809         \r\r
 810         if(iflag == 2){  // calculate gradient\r\r
 811           Int_t iParam = 0 ;\r\r
 812           efit = 0. ;\r\r
 813           while(iParam < nPar ){\r\r
 814             Double_t dx = ((Float_t)iphi - x[iParam]) ;\r\r
 815             iParam++ ;\r\r
 816             Double_t dz = ((Float_t)ieta - x[iParam]) ;\r\r
 817             iParam++ ;\r\r
 818             efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r\r
 819             iParam++ ;\r\r
 820           }\r\r
 821           \r\r
 822           Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r\r
 823           iParam = 0 ;\r\r
 824           while(iParam < nPar ){\r\r
 825             Double_t xpar = x[iParam] ;\r\r
 826             Double_t zpar = x[iParam+1] ;\r\r
 827             Double_t epar = x[iParam+2] ;\r\r
 828             \r\r
 829             Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r\r
 830             Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
 831             Double_t rp1  = TMath::Power(dr, fgSSPars[1]) ;\r\r
 832             Double_t rp5  = TMath::Power(dr, fgSSPars[5]) ;\r\r
 833             \r\r
 834             Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * \r\r
 835             (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - \r\r
 836              (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + \r\r
 837               fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );\r\r
 838             \r\r
 839             //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )\r\r
 840             //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r\r
 841             \r\r
 842             Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x\r\r
 843             iParam++ ;\r\r
 844             Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z\r\r
 845             iParam++ ;\r\r
 846             Grad[iParam] += shape ;                                  // Derivative over energy\r\r
 847             iParam++ ;\r\r
 848           }\r\r
 849         }\r\r
 850         efit = 0;\r\r
 851         iparam = 0 ;\r\r
 852         \r\r
 853         while(iparam < nPar ){\r\r
 854           Double_t xpar = x[iparam] ;\r\r
 855           Double_t zpar = x[iparam+1] ;\r\r
 856           Double_t epar = x[iparam+2] ;\r\r
 857           iparam += 3 ;\r\r
 858           efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
 859         }\r\r
 860         \r\r
 861         fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r\r
 862         // Here we assume, that sigma = sqrt(E) \r\r
 863         } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker\r\r
 864       } // digit loop\r\r
 865     } // recpoint, digits and geom not NULL\r\r
 866   }// List is not NULL\r\r
 867   \r\r
 868 }\r\r
 869 \r\r
 870 \r\r
 871 //____________________________________________________________________________\r\r
 872 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r\r
 873   for(UInt_t i=0;i<7;++i)\r\r
 874     fgSSPars[i]=pars[i];\r\r
 875   if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0\r\r
 876 }\r\r
 877 \r\r
 878 //____________________________________________________________________________\r\r
 879 void AliEMCALUnfolding::SetPar5(Double_t *pars){\r\r
 880   for(UInt_t i=0;i<3;++i)\r\r
 881     fgPar5[i]=pars[i];\r\r
 882 }\r\r
 883 \r\r
 884 //____________________________________________________________________________\r\r
 885 void AliEMCALUnfolding::SetPar6(Double_t *pars){\r\r
 886   for(UInt_t i=0;i<3;++i)\r\r
 887     fgPar6[i]=pars[i];\r\r
 888 }\r\r
 889 \r\r
 890 //____________________________________________________________________________\r\r
 891 void AliEMCALUnfolding::EvalPar5(Double_t phi){\r\r
 892   //\r\r
 893   //Evaluate the 5th parameter of the shower shape function\r\r
 894   //phi in degrees range (-10,10)\r\r
 895   //\r\r
 896   //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r\r
 897   fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];\r\r
 898 }\r\r
 899 \r\r
 900 //____________________________________________________________________________\r\r
 901 void AliEMCALUnfolding::EvalPar6(Double_t phi){\r\r
 902   //\r\r
 903   //Evaluate the 6th parameter of the shower shape function\r\r
 904   //phi in degrees range (-10,10)\r\r
 905   //\r\r
 906   //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r\r
 907   fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];\r\r
 908 }\r\r
 909 \r\r
 910 //____________________________________________________________________________\r\r
 911 void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){\r\r
 912   //\r\r
 913   // calculate params p5 and p6 depending on the phi angle in global coordinate\r\r
 914   // for the cell with given absId index\r\r
 915   //\r\r
 916   Double_t etaGlob = 0.;//eta in global c.s. - unused\r\r
 917   Double_t phiGlob = 0.;//phi in global c.s. in radians\r\r
 918   geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);\r\r
 919   phiGlob*=180./TMath::Pi();\r\r
 920   phiGlob-=90.;\r\r
 921   phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);\r\r
 922 \r\r
 923   EvalPar5(phiGlob);\r\r
 924   EvalPar6(phiGlob);\r\r
 925 }\r\r
 926 \r\r