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