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