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