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