EMCAL/EMCALrec/AliEMCALUnfolding.cxx

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