set size of geometry data member fParSM to 3, change corresponding method to get...
[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 && 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       if(recPoint){\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       } else AliError("RecPoint NULL");\r
208     } // rec point loop\r
209   }\r
210   // End of Unfolding of clusters\r
211 }\r
212 \r
213 //____________________________________________________________________________\r
214 Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, \r
215                                           Int_t nMax, \r
216                                           AliEMCALDigit ** maxAt, \r
217                                           Float_t * maxAtEnergy)\r
218 {\r
219   // Extended to whole EMCAL \r
220   // Performs the unfolding of a cluster with nMax overlapping showers \r
221   \r
222   Int_t nPar = 3 * nMax ;\r
223   Float_t * fitparameters = new Float_t[nPar] ;\r
224   \r
225   if (fGeom==0)\r
226     AliFatal("Did not get geometry from EMCALLoader") ;\r
227   \r
228   Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
229   if( !rv ) {\r
230     // Fit failed, return (and remove cluster? - why? I leave the cluster)\r
231     iniTower->SetNExMax(-1) ;\r
232     delete[] fitparameters ;\r
233     return kFALSE;\r
234   }\r
235   \r
236   // create unfolded rec points and fill them with new energy lists\r
237   // First calculate energy deposited in each sell in accordance with\r
238   // fit (without fluctuations): efit[]\r
239   // and later correct this number in acordance with actual energy\r
240   // deposition\r
241   \r
242   Int_t nDigits = iniTower->GetMultiplicity() ;\r
243   Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
244   Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r
245   \r
246   AliEMCALDigit * digit = 0 ;\r
247   Int_t * digitsList = iniTower->GetDigitsList() ;\r
248   \r
249   Int_t iSupMod =  0 ;\r
250   Int_t iTower  =  0 ;\r
251   Int_t iIphi   =  0 ;\r
252   Int_t iIeta   =  0 ;\r
253   Int_t iphi    =  0 ;//x direction\r
254   Int_t ieta    =  0 ;//z direstion\r
255   \r
256   Int_t iparam = 0 ;\r
257   Int_t iDigit = 0 ;\r
258   \r
259   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
260     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
261     if(digit){\r
262       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
263       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
264                                          iIphi, iIeta,iphi,ieta);\r
265       EvalParsPhiDependence(digit->GetId(),fGeom);\r
266       \r
267       efit[iDigit] = 0.;\r
268       iparam = 0;\r
269       while(iparam < nPar ){\r
270         xpar = fitparameters[iparam] ;\r
271         zpar = fitparameters[iparam+1] ;\r
272         epar = fitparameters[iparam+2] ;\r
273         iparam += 3 ;\r
274         \r
275         efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
276       }\r
277     } else AliError("Digit NULL!");\r
278     \r
279   }//digit loop\r
280   \r
281   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
282   // so that energy deposited in each cell is distributed between new clusters proportionally\r
283   // to its contribution to efit\r
284   \r
285   Float_t * energiesList = iniTower->GetEnergiesList() ;\r
286   Float_t ratio = 0 ;\r
287   \r
288   iparam = 0 ;\r
289   while(iparam < nPar ){\r
290     xpar = fitparameters[iparam] ;\r
291     zpar = fitparameters[iparam+1] ;\r
292     epar = fitparameters[iparam+2] ;\r
293     iparam += 3 ;\r
294     \r
295     AliEMCALRecPoint * recPoint = 0 ;\r
296     \r
297     if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
298       fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
299     \r
300     //add recpoint\r
301     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
302     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
303     \r
304     if(recPoint){\r
305       \r
306       fNumberOfECAClusters++ ;\r
307       recPoint->SetNExMax((Int_t)nPar/3) ;\r
308       \r
309       Float_t eDigit = 0. ;\r
310       for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
311         digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
312         if(digit){\r
313           fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
314           fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
315                                              iIphi, iIeta,iphi,ieta);\r
316           EvalParsPhiDependence(digit->GetId(),fGeom);\r
317           if(efit[iDigit]==0) continue;//just for sure\r
318           ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
319           eDigit = energiesList[iDigit] * ratio ;\r
320           recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
321         } else AliError("NULL digit");\r
322       }//digit loop \r
323     } else AliError("NULL RecPoint");\r
324   }//while\r
325   \r
326   delete[] fitparameters ;\r
327   delete[] efit ;\r
328   \r
329   return kTRUE;\r
330 }\r
331 \r
332 //____________________________________________________________________________\r
333 Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, \r
334                                         const Float_t* maxAtEnergy,\r
335                                         Int_t nPar, Float_t * fitparameters) const\r
336 {\r
337   // Calls TMinuit to fit the energy distribution of a cluster with several maxima\r
338   // The initial values for fitting procedure are set equal to the\r
339   // positions of local maxima.       \r
340   // Cluster will be fitted as a superposition of nPar/3\r
341   // electromagnetic showers\r
342 \r
343   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r
344         \r
345   if(!gMinuit)\r
346     gMinuit = new TMinuit(100) ;//max 100 parameters\r
347 \r
348   gMinuit->mncler();                     // Reset Minuit's list of paramters\r
349   gMinuit->SetPrintLevel(-1) ;           // No Printout\r
350   gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;\r
351   // To set the address of the minimization function\r
352   TList * toMinuit = new TList();\r
353   toMinuit->AddAt(recPoint,0) ;\r
354   toMinuit->AddAt(fDigitsArr,1) ;\r
355   toMinuit->AddAt(fGeom,2) ;\r
356 \r
357   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare\r
358 \r
359   // filling initial values for fit parameters\r
360   AliEMCALDigit * digit ;\r
361 \r
362   Int_t ierflg  = 0;\r
363   Int_t index   = 0 ;\r
364   Int_t nDigits = (Int_t) nPar / 3 ;\r
365 \r
366   Int_t iDigit ;\r
367 \r
368   Int_t iSupMod =  0 ;\r
369   Int_t iTower  =  0 ;\r
370   Int_t iIphi   =  0 ;\r
371   Int_t iIeta   =  0 ;\r
372   Int_t iphi    =  0 ;//x direction\r
373   Int_t ieta    =  0 ;//z direstion\r
374 \r
375   for(iDigit = 0; iDigit < nDigits; iDigit++){\r
376     digit = maxAt[iDigit];\r
377     fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
378     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
379                                        iIphi, iIeta,iphi,ieta);\r
380 \r
381     Float_t energy = maxAtEnergy[iDigit] ;\r
382 \r
383     //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original\r
384     gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;\r
385     index++ ;\r
386     if(ierflg != 0){\r
387       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;\r
388       toMinuit->Clear();\r
389       delete toMinuit ;\r
390       return kFALSE;\r
391     }\r
392     //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original\r
393     gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;\r
394     index++ ;\r
395     if(ierflg != 0){\r
396       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;\r
397       toMinuit->Clear();\r
398       delete toMinuit ;\r
399       return kFALSE;\r
400     }\r
401     //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original\r
402     gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r
403     index++ ;\r
404     if(ierflg != 0){\r
405       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;\r
406       toMinuit->Clear();\r
407       delete toMinuit ;\r
408       return kFALSE;\r
409     }\r
410   }\r
411 \r
412   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; \r
413                       // The number of function call slightly depends on it.\r
414   //  Double_t p1 = 1.0 ;// par to gradient \r
415   Double_t p2 = 0.0 ;\r
416   //  Double_t p3 = 3.0 ;\r
417   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls\r
418   //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient\r
419   gMinuit->SetMaxIterations(5);//was 5\r
420   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings\r
421   //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts\r
422 \r
423   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize\r
424   //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize\r
425   if(ierflg == 4){  // Minimum not found\r
426     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;\r
427     toMinuit->Clear();\r
428     delete toMinuit ;\r
429     return kFALSE ;\r
430   }\r
431   for(index = 0; index < nPar; index++){\r
432     Double_t err = 0. ;\r
433     Double_t val = 0. ;\r
434     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index\r
435     fitparameters[index] = val ;\r
436   }\r
437 \r
438   toMinuit->Clear();\r
439   delete toMinuit ;\r
440   return kTRUE;\r
441 \r
442 }\r
443 \r
444 //____________________________________________________________________________\r
445 Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)\r
446\r
447   // extended to whole EMCAL \r
448   // Shape of the shower\r
449   // If you change this function, change also the gradient evaluation in ChiSquare()\r
450 \r
451   Double_t r = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r
452   Double_t rp1  = TMath::Power(r, fSSPars[1]) ;\r
453   Double_t rp5  = TMath::Power(r, fSSPars[5]) ;\r
454   Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[6] * rp5) ) ) ;\r
455   return shape ;\r
456 }\r
457 \r
458 //____________________________________________________________________________\r
459 void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r
460                                                  Double_t & fret,\r
461                                                  Double_t * x, Int_t iflag)\r
462 {\r
463   // Calculates the Chi square for the cluster unfolding minimization\r
464   // Number of parameters, Gradient, Chi squared, parameters, what to do\r
465   \r
466   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
467   if(toMinuit){\r
468     AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;\r
469     TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;\r
470     // A bit buggy way to get an access to the geometry\r
471     // To be revised!\r
472     AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
473     \r
474     if(recPoint && digits && geom){\r
475       \r
476       Int_t * digitsList     = recPoint->GetDigitsList() ;\r
477       \r
478       Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
479       \r
480       Float_t * energiesList = recPoint->GetEnergiesList() ;\r
481       \r
482       fret = 0. ;\r
483       Int_t iparam = 0 ;\r
484       \r
485       if(iflag == 2)\r
486         for(iparam = 0 ; iparam < nPar ; iparam++)\r
487           Grad[iparam] = 0 ; // Will evaluate gradient\r
488       \r
489       Double_t efit = 0. ;\r
490       \r
491       AliEMCALDigit * digit ;\r
492       Int_t iDigit ;\r
493       \r
494       Int_t iSupMod =  0 ;\r
495       Int_t iTower  =  0 ;\r
496       Int_t iIphi   =  0 ;\r
497       Int_t iIeta   =  0 ;\r
498       Int_t iphi    =  0 ;//x direction\r
499       Int_t ieta    =  0 ;//z direstion\r
500       \r
501       \r
502       for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
503         if(energiesList[iDigit]==0) continue;\r
504         \r
505         digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
506         \r
507         if(digit){\r
508         geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
509         geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
510                                           iIphi, iIeta,iphi,ieta);\r
511         EvalParsPhiDependence(digit->GetId(),geom);\r
512         \r
513         if(iflag == 2){  // calculate gradient\r
514           Int_t iParam = 0 ;\r
515           efit = 0. ;\r
516           while(iParam < nPar ){\r
517             Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
518             iParam++ ;\r
519             Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
520             iParam++ ;\r
521             efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
522             iParam++ ;\r
523           }\r
524           \r
525           Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
526           iParam = 0 ;\r
527           while(iParam < nPar ){\r
528             Double_t xpar = x[iParam] ;\r
529             Double_t zpar = x[iParam+1] ;\r
530             Double_t epar = x[iParam+2] ;\r
531             \r
532             Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r
533             Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
534             Double_t rp1  = TMath::Power(dr, fSSPars[1]) ;\r
535             Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;\r
536             \r
537             Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
538             (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
539              (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
540               fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\r
541             \r
542             //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )\r
543             //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
544             \r
545             Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x\r
546             iParam++ ;\r
547             Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z\r
548             iParam++ ;\r
549             Grad[iParam] += shape ;                                  // Derivative over energy\r
550             iParam++ ;\r
551           }\r
552         }\r
553         efit = 0;\r
554         iparam = 0 ;\r
555         \r
556         while(iparam < nPar ){\r
557           Double_t xpar = x[iparam] ;\r
558           Double_t zpar = x[iparam+1] ;\r
559           Double_t epar = x[iparam+2] ;\r
560           iparam += 3 ;\r
561           efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
562         }\r
563         \r
564         fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
565         // Here we assume, that sigma = sqrt(E) \r
566         } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n");\r
567       } // digit loop\r
568     } // recpoint, digits and geom not NULL\r
569   }// List is not NULL\r
570   \r
571 }\r
572 \r
573 \r
574 //____________________________________________________________________________\r
575 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r
576   for(UInt_t i=0;i<7;++i)\r
577     fSSPars[i]=pars[i];\r
578   if(pars[2]==0. && pars[3]==0.) fSSPars[2]=1.;//to avoid dividing by 0\r
579 }\r
580 \r
581 //____________________________________________________________________________\r
582 void AliEMCALUnfolding::SetPar5(Double_t *pars){\r
583   for(UInt_t i=0;i<3;++i)\r
584     fPar5[i]=pars[i];\r
585 }\r
586 \r
587 //____________________________________________________________________________\r
588 void AliEMCALUnfolding::SetPar6(Double_t *pars){\r
589   for(UInt_t i=0;i<3;++i)\r
590     fPar6[i]=pars[i];\r
591 }\r
592 \r
593 //____________________________________________________________________________\r
594 void AliEMCALUnfolding::EvalPar5(Double_t phi){\r
595   //\r
596   //Evaluate the 5th parameter of the shower shape function\r
597   //phi in degrees range (-10,10)\r
598   //\r
599   //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r
600   fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[2];\r
601 }\r
602 \r
603 //____________________________________________________________________________\r
604 void AliEMCALUnfolding::EvalPar6(Double_t phi){\r
605   //\r
606   //Evaluate the 6th parameter of the shower shape function\r
607   //phi in degrees range (-10,10)\r
608   //\r
609   //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r
610   fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r
611 }\r
612 \r
613 //____________________________________________________________________________\r
614 void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom){\r
615   //\r
616   // calculate params p5 and p6 depending on the phi angle in global coordinate\r
617   // for the cell with given absId index\r
618   //\r
619   Double_t etaGlob = 0.;//eta in global c.s. - unused\r
620   Double_t phiGlob = 0.;//phi in global c.s. in radians\r
621   geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);\r
622   phiGlob*=180./TMath::Pi();\r
623   phiGlob-=90.;\r
624   phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);\r
625 \r
626   EvalPar5(phiGlob);\r
627   EvalPar6(phiGlob);\r
628 }\r
629 \r