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