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