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