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