]>
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 | |
2037b85d | 202 | |
35abc2bd | 203 | //**************************** part 1 ******************************************* |
204 | // Performs the unfolding of a cluster with nMax overlapping showers | |
205 | ||
2037b85d | 206 | //printf("Original cluster E %f\n",iniTower->GetEnergy()); |
207 | ||
35abc2bd | 208 | Int_t nPar = 3 * nMax ; |
209 | Float_t * fitparameters = new Float_t[nPar] ; | |
210 | ||
211 | if (fGeom==0) | |
212 | AliFatal("Did not get geometry from EMCALLoader") ; | |
213 | ||
214 | Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; | |
2037b85d | 215 | if( !rv ) |
216 | { | |
35abc2bd | 217 | // Fit failed, return (and remove cluster? - why? I leave the cluster) |
218 | iniTower->SetNExMax(-1) ; | |
219 | delete[] fitparameters ; | |
220 | return kFALSE; | |
221 | } | |
222 | ||
223 | //**************************** part 2 ******************************************* | |
224 | // create unfolded rec points and fill them with new energy lists | |
225 | // First calculate energy deposited in each sell in accordance with | |
226 | // fit (without fluctuations): efit[] | |
227 | // and later correct this number in acordance with actual energy | |
228 | // deposition | |
229 | ||
230 | Int_t nDigits = iniTower->GetMultiplicity() ; | |
231 | Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells | |
232 | Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units | |
233 | ||
234 | AliEMCALDigit * digit = 0 ; | |
235 | Int_t * digitsList = iniTower->GetDigitsList() ; | |
236 | ||
237 | Int_t iSupMod = 0 ; | |
238 | Int_t iTower = 0 ; | |
239 | Int_t iIphi = 0 ; | |
240 | Int_t iIeta = 0 ; | |
241 | Int_t iphi = 0 ;//x direction | |
242 | Int_t ieta = 0 ;//z direstion | |
243 | ||
244 | Int_t iparam = 0 ; | |
245 | Int_t iDigit = 0 ; | |
246 | ||
2037b85d | 247 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) |
248 | { | |
35abc2bd | 249 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ; |
2037b85d | 250 | if(digit) |
251 | { | |
35abc2bd | 252 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); |
253 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
254 | iIphi, iIeta,iphi,ieta); | |
255 | EvalParsPhiDependence(digit->GetId(),fGeom); | |
256 | ||
257 | efit[iDigit] = 0.; | |
258 | iparam = 0; | |
2037b85d | 259 | while(iparam < nPar ) |
260 | { | |
35abc2bd | 261 | xpar = fitparameters[iparam] ; |
262 | zpar = fitparameters[iparam+1] ; | |
263 | epar = fitparameters[iparam+2] ; | |
264 | iparam += 3 ; | |
265 | ||
266 | efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; | |
267 | } | |
2037b85d | 268 | } else AliDebug(1,"Digit NULL part 2!"); |
35abc2bd | 269 | |
270 | }//digit loop | |
271 | ||
272 | //**************************** part 3 ******************************************* | |
273 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations | |
274 | // so that energy deposited in each cell is distributed between new clusters proportionally | |
275 | // to its contribution to efit | |
276 | ||
277 | Float_t * energiesList = iniTower->GetEnergiesList() ; | |
278 | Float_t ratio = 0 ; | |
279 | Float_t eDigit = 0. ; | |
280 | Int_t nSplittedClusters=(Int_t)nPar/3; | |
281 | ||
282 | Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters]; | |
283 | //above - temporary table with energies after unfolding. | |
284 | //the orderis following: | |
285 | //first cluster <first cell - last cell>, | |
286 | //second cluster <first cell - last cell>, etc. | |
2037b85d | 287 | |
35abc2bd | 288 | //**************************** sub-part 3.1 ************************************* |
289 | //here we check if energy of the cell in the cluster after unfolding is above threshold. | |
290 | //If not the energy from a given cell in the cluster is divided in correct proportions | |
291 | //in accordance to the other clusters and added to them and set to 0. | |
2037b85d | 292 | |
35abc2bd | 293 | iparam = 0 ; |
2037b85d | 294 | while(iparam < nPar ) |
295 | { | |
35abc2bd | 296 | xpar = fitparameters[iparam] ; |
297 | zpar = fitparameters[iparam+1] ; | |
298 | epar = fitparameters[iparam+2] ; | |
2037b85d | 299 | |
300 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) | |
301 | { | |
35abc2bd | 302 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
2037b85d | 303 | if(digit) |
304 | { | |
305 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); | |
306 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
307 | iIphi, iIeta,iphi,ieta); | |
308 | ||
309 | EvalParsPhiDependence(digit->GetId(),fGeom); | |
310 | ||
311 | if(efit[iDigit]==0) | |
312 | {//just for sure | |
313 | correctedEnergyList[iparam/3+iDigit] = 0; | |
314 | continue; | |
315 | } | |
316 | ||
317 | ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; | |
318 | eDigit = energiesList[iDigit] * ratio ; | |
319 | ||
320 | //add energy to temporary matrix | |
321 | correctedEnergyList[iparam/3+iDigit] = eDigit; | |
322 | ||
323 | } else AliDebug(1,"NULL digit part 3"); | |
35abc2bd | 324 | }//digit loop |
325 | iparam += 3 ; | |
326 | }//while | |
2037b85d | 327 | |
35abc2bd | 328 | //**************************** sub-part 3.2 ************************************* |
329 | //here we correct energy for each cell and cluster | |
330 | Float_t maximumEne=0; | |
331 | Int_t maximumIndex=0; | |
332 | Bool_t isAnyBelowThreshold=kFALSE; | |
333 | // Float_t Threshold=0.01; | |
334 | Float_t * energyFraction = new Float_t[nSplittedClusters]; | |
335 | Int_t iparam2 = 0 ; | |
2037b85d | 336 | for(iDigit = 0 ; iDigit < nDigits ; iDigit++) |
337 | { | |
35abc2bd | 338 | isAnyBelowThreshold=kFALSE; |
339 | maximumEne=0; | |
2037b85d | 340 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
341 | { | |
35abc2bd | 342 | |
343 | if(correctedEnergyList[iparam/3+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE; | |
2037b85d | 344 | if(correctedEnergyList[iparam/3+iDigit] > maximumEne) |
345 | { | |
346 | maximumEne = correctedEnergyList[iparam/3+iDigit]; | |
347 | maximumIndex = iparam; | |
35abc2bd | 348 | } |
349 | }//end of loop over clusters after unfolding | |
2037b85d | 350 | |
35abc2bd | 351 | if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold |
2037b85d | 352 | |
353 | //printf("Correct E cell %f < %f, org Digit index %d, e = %f\n",correctedEnergyList[iparam/3+iDigit],fThreshold,iDigit, energiesList[iDigit]); | |
354 | //if( energiesList[iDigit] < correctedEnergyList[iparam/3+iDigit]) printf("\t What? \n"); | |
355 | ||
356 | if(maximumEne < fThreshold) | |
357 | {//add all cluster cells and put energy into max index, other set to 0 | |
35abc2bd | 358 | maximumEne=0.; |
2037b85d | 359 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
360 | { | |
361 | maximumEne+=correctedEnergyList[iparam/3+iDigit]; | |
362 | correctedEnergyList[iparam/3+iDigit]=0; | |
35abc2bd | 363 | } |
364 | correctedEnergyList[maximumIndex/3+iDigit]=maximumEne; | |
365 | continue; | |
366 | }//end if | |
2037b85d | 367 | |
35abc2bd | 368 | //divide energy of cell below threshold in the correct proportion and add to other cells |
369 | maximumEne=0;//not used any more so use it for the energy sum | |
2037b85d | 370 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
371 | {//calculate energy sum | |
35abc2bd | 372 | if(correctedEnergyList[iparam/3+iDigit] < fThreshold) energyFraction[iparam/3]=0; |
2037b85d | 373 | else |
374 | { | |
375 | energyFraction[iparam/3]=1; | |
376 | maximumEne+=correctedEnergyList[iparam/3+iDigit]; | |
35abc2bd | 377 | } |
378 | }//end of loop over clusters after unfolding | |
2037b85d | 379 | if(maximumEne>0) |
380 | { | |
35abc2bd | 381 | for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction |
2037b85d | 382 | energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3+iDigit] / maximumEne; |
35abc2bd | 383 | } |
384 | ||
2037b85d | 385 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
386 | {//add energy from cells below threshold to others | |
387 | if(energyFraction[iparam/3]>0) continue; | |
388 | else | |
389 | { | |
390 | for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3) | |
391 | { | |
392 | correctedEnergyList[iparam2/3+iDigit] += (energyFraction[iparam2/3] * | |
393 | correctedEnergyList[iparam/3+iDigit]) ; | |
394 | }//inner loop | |
395 | correctedEnergyList[iparam/3+iDigit] = 0; | |
396 | } | |
35abc2bd | 397 | } |
398 | } | |
2037b85d | 399 | else |
400 | { | |
35abc2bd | 401 | //digit energy to be set to 0 |
2037b85d | 402 | for(iparam = 0 ; iparam < nPar ; iparam+=3) |
403 | { | |
404 | correctedEnergyList[iparam/3+iDigit] = 0; | |
35abc2bd | 405 | } |
406 | }//new adam correction for is energy>0 | |
2037b85d | 407 | |
35abc2bd | 408 | }//end of loop over digits |
409 | delete[] energyFraction; | |
2037b85d | 410 | |
35abc2bd | 411 | //**************************** sub-part 3.3 ************************************* |
412 | //here we add digits to recpoints with corrected energy | |
413 | iparam = 0 ; | |
2037b85d | 414 | while(iparam < nPar ) |
415 | { | |
35abc2bd | 416 | AliEMCALRecPoint * recPoint = 0 ; |
417 | ||
418 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) | |
419 | fRecPoints->Expand(2*fNumberOfECAClusters) ; | |
420 | ||
421 | //add recpoint | |
422 | (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; | |
423 | recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ; | |
424 | ||
2037b85d | 425 | if(recPoint) |
426 | { | |
35abc2bd | 427 | fNumberOfECAClusters++ ; |
428 | recPoint->SetNExMax(nSplittedClusters) ; | |
429 | ||
2037b85d | 430 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) |
431 | { | |
35abc2bd | 432 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; |
2037b85d | 433 | |
434 | ||
435 | if(digit && correctedEnergyList[iparam/3+iDigit]>0. ) | |
436 | { | |
437 | //if(correctedEnergyList[iparam/3+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3+iDigit],fThreshold); | |
35abc2bd | 438 | recPoint->AddDigit( *digit, correctedEnergyList[iparam/3+iDigit], kFALSE ) ; //FIXME, need to study the shared case |
2037b85d | 439 | } else |
440 | { | |
441 | AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3+iDigit])); | |
442 | //cout<<"nDigits "<<nDigits<<" iParam/3 "<<iparam/3<< endl; | |
443 | } | |
35abc2bd | 444 | }//digit loop |
445 | } else AliError("NULL RecPoint"); | |
2037b85d | 446 | |
35abc2bd | 447 | //protection from recpoint with no digits |
448 | //cout<<"multi rec "<<recPoint->GetMultiplicity()<<endl; | |
2037b85d | 449 | if(recPoint->GetMultiplicity()==0) |
450 | { | |
35abc2bd | 451 | delete (*fRecPoints)[fNumberOfECAClusters]; |
452 | //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl; | |
453 | fRecPoints->RemoveAt(fNumberOfECAClusters); | |
454 | //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl; | |
455 | fNumberOfECAClusters--; | |
456 | nSplittedClusters--; | |
2037b85d | 457 | |
35abc2bd | 458 | } |
2037b85d | 459 | |
35abc2bd | 460 | iparam += 3 ; |
461 | }//while | |
462 | ||
463 | delete[] fitparameters ; | |
464 | delete[] efit ; | |
465 | delete[] correctedEnergyList ; | |
2037b85d | 466 | |
35abc2bd | 467 | return kTRUE; |
468 | } | |
469 | ||
470 | ||
471 | //____________________________________________________________________________ | |
472 | Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, | |
473 | Int_t nMax, | |
474 | AliEMCALDigit ** maxAt, | |
475 | Float_t * maxAtEnergy) | |
476 | { | |
477 | // Extended to whole EMCAL | |
478 | // Performs the unfolding of a cluster with nMax overlapping showers | |
479 | ||
480 | Int_t nPar = 3 * nMax ; | |
481 | Float_t * fitparameters = new Float_t[nPar] ; | |
482 | ||
483 | if (fGeom==0) | |
484 | AliFatal("Did not get geometry from EMCALLoader") ; | |
485 | ||
486 | Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; | |
487 | if( !rv ) { | |
488 | // Fit failed, return (and remove cluster? - why? I leave the cluster) | |
489 | iniTower->SetNExMax(-1) ; | |
490 | delete[] fitparameters ; | |
491 | return kFALSE; | |
492 | } | |
493 | ||
494 | // create unfolded rec points and fill them with new energy lists | |
495 | // First calculate energy deposited in each sell in accordance with | |
496 | // fit (without fluctuations): efit[] | |
497 | // and later correct this number in acordance with actual energy | |
498 | // deposition | |
499 | ||
500 | Int_t nDigits = iniTower->GetMultiplicity() ; | |
501 | Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells | |
502 | Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units | |
503 | ||
504 | AliEMCALDigit * digit = 0 ; | |
505 | Int_t * digitsList = iniTower->GetDigitsList() ; | |
506 | ||
507 | Int_t iSupMod = 0 ; | |
508 | Int_t iTower = 0 ; | |
509 | Int_t iIphi = 0 ; | |
510 | Int_t iIeta = 0 ; | |
511 | Int_t iphi = 0 ;//x direction | |
512 | Int_t ieta = 0 ;//z direstion | |
513 | ||
514 | Int_t iparam = 0 ; | |
515 | Int_t iDigit = 0 ; | |
516 | ||
517 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
518 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ; | |
519 | if(digit){ | |
520 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); | |
521 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
522 | iIphi, iIeta,iphi,ieta); | |
523 | EvalParsPhiDependence(digit->GetId(),fGeom); | |
524 | ||
525 | efit[iDigit] = 0.; | |
526 | iparam = 0; | |
527 | while(iparam < nPar ){ | |
528 | xpar = fitparameters[iparam] ; | |
529 | zpar = fitparameters[iparam+1] ; | |
530 | epar = fitparameters[iparam+2] ; | |
531 | iparam += 3 ; | |
532 | ||
533 | efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; | |
534 | } | |
535 | } else AliError("Digit NULL!"); | |
536 | ||
537 | }//digit loop | |
538 | ||
539 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations | |
540 | // so that energy deposited in each cell is distributed between new clusters proportionally | |
541 | // to its contribution to efit | |
542 | ||
543 | Float_t * energiesList = iniTower->GetEnergiesList() ; | |
544 | Float_t ratio = 0 ; | |
545 | ||
546 | iparam = 0 ; | |
547 | while(iparam < nPar ){ | |
548 | xpar = fitparameters[iparam] ; | |
549 | zpar = fitparameters[iparam+1] ; | |
550 | epar = fitparameters[iparam+2] ; | |
551 | iparam += 3 ; | |
552 | ||
553 | AliEMCALRecPoint * recPoint = 0 ; | |
554 | ||
555 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) | |
556 | fRecPoints->Expand(2*fNumberOfECAClusters) ; | |
557 | ||
558 | //add recpoint | |
559 | (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; | |
560 | recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ; | |
561 | ||
562 | if(recPoint){ | |
563 | ||
564 | fNumberOfECAClusters++ ; | |
565 | recPoint->SetNExMax((Int_t)nPar/3) ; | |
566 | ||
567 | Float_t eDigit = 0. ; | |
568 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
569 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; | |
570 | if(digit){ | |
571 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); | |
572 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
573 | iIphi, iIeta,iphi,ieta); | |
574 | EvalParsPhiDependence(digit->GetId(),fGeom); | |
575 | if(efit[iDigit]==0) continue;//just for sure | |
576 | ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; | |
577 | eDigit = energiesList[iDigit] * ratio ; | |
578 | recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case | |
579 | } else AliError("NULL digit"); | |
580 | }//digit loop | |
581 | } else AliError("NULL RecPoint"); | |
582 | }//while | |
583 | ||
584 | delete[] fitparameters ; | |
585 | delete[] efit ; | |
586 | ||
587 | return kTRUE; | |
588 | } | |
589 | ||
590 | ||
591 | //____________________________________________________________________________ | |
592 | Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, | |
593 | const Float_t* maxAtEnergy, | |
594 | Int_t nPar, Float_t * fitparameters) const | |
595 | { | |
596 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima | |
597 | // The initial values for fitting procedure are set equal to the | |
598 | // positions of local maxima. | |
599 | // Cluster will be fitted as a superposition of nPar/3 | |
600 | // electromagnetic showers | |
601 | ||
602 | if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); | |
603 | ||
604 | if(!gMinuit) | |
605 | gMinuit = new TMinuit(100) ;//max 100 parameters | |
606 | ||
607 | gMinuit->mncler(); // Reset Minuit's list of paramters | |
608 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
609 | gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ; | |
610 | // To set the address of the minimization function | |
611 | TList * toMinuit = new TList(); | |
612 | toMinuit->AddAt(recPoint,0) ; | |
613 | toMinuit->AddAt(fDigitsArr,1) ; | |
614 | toMinuit->AddAt(fGeom,2) ; | |
615 | ||
616 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare | |
617 | ||
618 | // filling initial values for fit parameters | |
619 | AliEMCALDigit * digit ; | |
620 | ||
621 | Int_t ierflg = 0; | |
622 | Int_t index = 0 ; | |
623 | Int_t nDigits = (Int_t) nPar / 3 ; | |
624 | ||
625 | Int_t iDigit ; | |
626 | ||
627 | Int_t iSupMod = 0 ; | |
628 | Int_t iTower = 0 ; | |
629 | Int_t iIphi = 0 ; | |
630 | Int_t iIeta = 0 ; | |
631 | Int_t iphi = 0 ;//x direction | |
632 | Int_t ieta = 0 ;//z direstion | |
633 | ||
634 | for(iDigit = 0; iDigit < nDigits; iDigit++){ | |
635 | digit = maxAt[iDigit]; | |
636 | if(digit==0) AliError("energy of digit = 0!"); | |
637 | fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); | |
638 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
639 | iIphi, iIeta,iphi,ieta); | |
640 | ||
641 | Float_t energy = maxAtEnergy[iDigit] ; | |
642 | ||
643 | //gMinuit->mnparm(index, "x", iphi, 0.1, 0, 0, ierflg) ;//original | |
644 | gMinuit->mnparm(index, "x", iphi, 0.05, 0, 0, ierflg) ; | |
645 | index++ ; | |
646 | if(ierflg != 0){ | |
647 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %d", iphi ) ; | |
648 | toMinuit->Clear(); | |
649 | delete toMinuit ; | |
650 | return kFALSE; | |
651 | } | |
652 | //gMinuit->mnparm(index, "z", ieta, 0.1, 0, 0, ierflg) ;//original | |
653 | gMinuit->mnparm(index, "z", ieta, 0.05, 0, 0, ierflg) ; | |
654 | index++ ; | |
655 | if(ierflg != 0){ | |
656 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %d", ieta) ; | |
657 | toMinuit->Clear(); | |
658 | delete toMinuit ; | |
659 | return kFALSE; | |
660 | } | |
661 | //gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original | |
662 | gMinuit->mnparm(index, "Energy", energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05 | |
663 | index++ ; | |
664 | if(ierflg != 0){ | |
665 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ; | |
666 | toMinuit->Clear(); | |
667 | delete toMinuit ; | |
668 | return kFALSE; | |
669 | } | |
670 | } | |
671 | ||
672 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; | |
673 | // The number of function call slightly depends on it. | |
674 | // Double_t p1 = 1.0 ;// par to gradient | |
675 | Double_t p2 = 0.0 ; | |
676 | // Double_t p3 = 3.0 ; | |
677 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls | |
678 | // gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient | |
679 | gMinuit->SetMaxIterations(5);//was 5 | |
680 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
681 | //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ; // printouts | |
682 | ||
683 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize | |
684 | //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ; // minimize | |
685 | if(ierflg == 4){ // Minimum not found | |
2037b85d | 686 | AliDebug(1,"EMCAL Unfolding Fit not converged, cluster abandoned " ) ; |
35abc2bd | 687 | toMinuit->Clear(); |
688 | delete toMinuit ; | |
689 | return kFALSE ; | |
690 | } | |
691 | for(index = 0; index < nPar; index++){ | |
692 | Double_t err = 0. ; | |
693 | Double_t val = 0. ; | |
694 | gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index | |
695 | fitparameters[index] = val ; | |
696 | } | |
697 | ||
698 | toMinuit->Clear(); | |
699 | delete toMinuit ; | |
700 | return kTRUE; | |
701 | ||
702 | } | |
703 | ||
704 | //____________________________________________________________________________ | |
705 | Double_t AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y) | |
706 | { | |
707 | // extended to whole EMCAL | |
708 | // Shape of the shower | |
709 | // If you change this function, change also the gradient evaluation in ChiSquare() | |
710 | ||
1186cd2b | 711 | Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y); |
712 | Double_t rp1 = TMath::Power(r, fgSSPars[1]) ; | |
713 | Double_t rp5 = TMath::Power(r, fgSSPars[5]) ; | |
714 | Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ; | |
35abc2bd | 715 | return shape ; |
716 | } | |
717 | ||
718 | //____________________________________________________________________________ | |
719 | void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, | |
1186cd2b | 720 | Double_t & fret, |
721 | Double_t * x, Int_t iflag) | |
35abc2bd | 722 | { |
723 | // Calculates the Chi square for the cluster unfolding minimization | |
724 | // Number of parameters, Gradient, Chi squared, parameters, what to do | |
725 | ||
726 | TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ; | |
727 | if(toMinuit){ | |
728 | AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ; | |
729 | TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ; | |
730 | // A bit buggy way to get an access to the geometry | |
731 | // To be revised! | |
732 | AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2)); | |
733 | ||
734 | if(recPoint && digits && geom){ | |
735 | ||
736 | Int_t * digitsList = recPoint->GetDigitsList() ; | |
737 | ||
738 | Int_t nOdigits = recPoint->GetDigitsMultiplicity() ; | |
739 | ||
740 | Float_t * energiesList = recPoint->GetEnergiesList() ; | |
741 | ||
742 | fret = 0. ; | |
743 | Int_t iparam = 0 ; | |
744 | ||
745 | if(iflag == 2) | |
746 | for(iparam = 0 ; iparam < nPar ; iparam++) | |
747 | Grad[iparam] = 0 ; // Will evaluate gradient | |
748 | ||
749 | Double_t efit = 0. ; | |
750 | ||
751 | AliEMCALDigit * digit ; | |
752 | Int_t iDigit ; | |
753 | ||
754 | Int_t iSupMod = 0 ; | |
755 | Int_t iTower = 0 ; | |
756 | Int_t iIphi = 0 ; | |
757 | Int_t iIeta = 0 ; | |
758 | Int_t iphi = 0 ;//x direction | |
759 | Int_t ieta = 0 ;//z direstion | |
760 | ||
761 | ||
762 | for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { | |
763 | if(energiesList[iDigit]==0) continue; | |
764 | ||
765 | digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) ); | |
766 | ||
767 | if(digit){ | |
768 | geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); | |
769 | geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, | |
770 | iIphi, iIeta,iphi,ieta); | |
771 | EvalParsPhiDependence(digit->GetId(),geom); | |
772 | ||
773 | if(iflag == 2){ // calculate gradient | |
774 | Int_t iParam = 0 ; | |
775 | efit = 0. ; | |
776 | while(iParam < nPar ){ | |
777 | Double_t dx = ((Float_t)iphi - x[iParam]) ; | |
778 | iParam++ ; | |
779 | Double_t dz = ((Float_t)ieta - x[iParam]) ; | |
780 | iParam++ ; | |
781 | efit += x[iParam] * ShowerShapeV2(dx,dz) ; | |
782 | iParam++ ; | |
783 | } | |
784 | ||
785 | Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E) | |
786 | iParam = 0 ; | |
787 | while(iParam < nPar ){ | |
788 | Double_t xpar = x[iParam] ; | |
789 | Double_t zpar = x[iParam+1] ; | |
790 | Double_t epar = x[iParam+2] ; | |
791 | ||
1186cd2b | 792 | Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) ); |
35abc2bd | 793 | Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; |
1186cd2b | 794 | Double_t rp1 = TMath::Power(dr, fgSSPars[1]) ; |
795 | Double_t rp5 = TMath::Power(dr, fgSSPars[5]) ; | |
35abc2bd | 796 | |
1186cd2b | 797 | Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * |
798 | (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - | |
799 | (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + | |
800 | fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) ); | |
35abc2bd | 801 | |
802 | //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) | |
803 | // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; | |
804 | ||
805 | Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x | |
806 | iParam++ ; | |
807 | Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z | |
808 | iParam++ ; | |
809 | Grad[iParam] += shape ; // Derivative over energy | |
810 | iParam++ ; | |
811 | } | |
812 | } | |
813 | efit = 0; | |
814 | iparam = 0 ; | |
815 | ||
816 | while(iparam < nPar ){ | |
817 | Double_t xpar = x[iparam] ; | |
818 | Double_t zpar = x[iparam+1] ; | |
819 | Double_t epar = x[iparam+2] ; | |
820 | iparam += 3 ; | |
821 | efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; | |
822 | } | |
823 | ||
824 | fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ; | |
825 | // Here we assume, that sigma = sqrt(E) | |
f11b2cf8 | 826 | } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker |
35abc2bd | 827 | } // digit loop |
828 | } // recpoint, digits and geom not NULL | |
829 | }// List is not NULL | |
830 | ||
831 | } | |
832 | ||
833 | ||
834 | //____________________________________________________________________________ | |
835 | void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){ | |
836 | for(UInt_t i=0;i<7;++i) | |
1186cd2b | 837 | fgSSPars[i]=pars[i]; |
838 | if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0 | |
35abc2bd | 839 | } |
840 | ||
841 | //____________________________________________________________________________ | |
842 | void AliEMCALUnfolding::SetPar5(Double_t *pars){ | |
843 | for(UInt_t i=0;i<3;++i) | |
1186cd2b | 844 | fgPar5[i]=pars[i]; |
35abc2bd | 845 | } |
846 | ||
847 | //____________________________________________________________________________ | |
848 | void AliEMCALUnfolding::SetPar6(Double_t *pars){ | |
849 | for(UInt_t i=0;i<3;++i) | |
1186cd2b | 850 | fgPar6[i]=pars[i]; |
35abc2bd | 851 | } |
852 | ||
853 | //____________________________________________________________________________ | |
854 | void AliEMCALUnfolding::EvalPar5(Double_t phi){ | |
855 | // | |
856 | //Evaluate the 5th parameter of the shower shape function | |
857 | //phi in degrees range (-10,10) | |
858 | // | |
859 | //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936; | |
1186cd2b | 860 | fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2]; |
35abc2bd | 861 | } |
862 | ||
863 | //____________________________________________________________________________ | |
864 | void AliEMCALUnfolding::EvalPar6(Double_t phi){ | |
865 | // | |
866 | //Evaluate the 6th parameter of the shower shape function | |
867 | //phi in degrees range (-10,10) | |
868 | // | |
869 | //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361; | |
1186cd2b | 870 | fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2]; |
35abc2bd | 871 | } |
872 | ||
873 | //____________________________________________________________________________ | |
1186cd2b | 874 | void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){ |
35abc2bd | 875 | // |
876 | // calculate params p5 and p6 depending on the phi angle in global coordinate | |
877 | // for the cell with given absId index | |
878 | // | |
879 | Double_t etaGlob = 0.;//eta in global c.s. - unused | |
880 | Double_t phiGlob = 0.;//phi in global c.s. in radians | |
881 | geom->EtaPhiFromIndex(absId, etaGlob, phiGlob); | |
882 | phiGlob*=180./TMath::Pi(); | |
883 | phiGlob-=90.; | |
884 | phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20); | |
885 | ||
886 | EvalPar5(phiGlob); | |
887 | EvalPar6(phiGlob); | |
888 | } | |
889 |