]>
Commit | Line | Data |
---|---|---|
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 | |
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 | |
d83af7ae | 79 | \r |
65bec413 | 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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 |