]> 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
Getting keyword substitution to work.
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALUnfolding.cxx
CommitLineData
65bec413 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
61class AliCDBStorage;\r
62#include "AliCDBEntry.h"\r
63\r
64Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r
65Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r
66Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};\r
67\r
68ClassImp(AliEMCALUnfolding)\r
69 \r
70//____________________________________________________________________________\r
71AliEMCALUnfolding::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
d83af7ae 79 \r
65bec413 80 Init() ;\r
81}\r
82\r
83//____________________________________________________________________________\r
84AliEMCALUnfolding::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
102AliEMCALUnfolding::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
126void 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
d83af7ae 130\r
65bec413 131 AliRunLoader *rl = AliRunLoader::Instance();\r
d83af7ae 132 if (rl && rl->GetAliRun()){\r
65bec413 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
154void 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
164void AliEMCALUnfolding::MakeUnfolding()\r
165{\r
166 // Unfolds clusters using the shape of an ElectroMagnetic shower\r
167 // Performs unfolding of all clusters\r
622e10be 168 \r
65bec413 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
622e10be 173 \r
65bec413 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
622e10be 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
65bec413 209 }\r
210 // End of Unfolding of clusters\r
211}\r
212\r
213//____________________________________________________________________________\r
214Bool_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
622e10be 224 \r
65bec413 225 if (fGeom==0)\r
226 AliFatal("Did not get geometry from EMCALLoader") ;\r
622e10be 227 \r
65bec413 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
622e10be 235 \r
65bec413 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
622e10be 241 \r
65bec413 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
622e10be 245 \r
65bec413 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
622e10be 255 \r
65bec413 256 Int_t iparam = 0 ;\r
257 Int_t iDigit = 0 ;\r
622e10be 258 \r
65bec413 259 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
260 digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
622e10be 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
65bec413 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
622e10be 284 \r
65bec413 285 Float_t * energiesList = iniTower->GetEnergiesList() ;\r
286 Float_t ratio = 0 ;\r
622e10be 287 \r
65bec413 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
622e10be 294 \r
65bec413 295 AliEMCALRecPoint * recPoint = 0 ;\r
622e10be 296 \r
65bec413 297 if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
298 fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
622e10be 299 \r
65bec413 300 //add recpoint\r
301 (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
302 recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
622e10be 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
65bec413 326 delete[] fitparameters ;\r
327 delete[] efit ;\r
622e10be 328 \r
65bec413 329 return kTRUE;\r
330}\r
331\r
332//____________________________________________________________________________\r
333Bool_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
65bec413 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
445Double_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
459void 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
622e10be 465 \r
65bec413 466 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
622e10be 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
65bec413 571}\r
572\r
573\r
574//____________________________________________________________________________\r
575void 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
582void 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
588void 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
594void 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
604void 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
614void 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