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