EMCAL/AliEMCALUnfolding.cxx

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