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