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