remove exec
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALUnfolding.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 //  Base class for the cluster unfolding algorithm 
18 //*-- Author: Adam Matyja (SUBATECH)
19 //  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis
20 //-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1
21 //-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)
22 //
23 //  unfolds the clusters having several local maxima. 
24 //////////////////////////////////////////////////////////////////////////////
25
26 // --- ROOT system ---
27 #include "TClonesArray.h"
28 #include <TMath.h> 
29 #include <TMinuit.h>
30
31 // --- Standard library ---
32 #include <cassert>
33
34 // --- AliRoot header files ---
35 #include "AliEMCALUnfolding.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliRunLoader.h"
38 #include "AliRun.h"
39 #include "AliEMCAL.h"
40 #include "AliEMCALRecParam.h"
41 #include "AliEMCALRecPoint.h"
42 #include "AliEMCALDigit.h"
43 #include "AliEMCALReconstructor.h"
44
45 #include "AliLog.h"
46 #include "AliCDBManager.h"
47 class AliCDBStorage;
48 #include "AliCDBEntry.h"
49
50 Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};
51 Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};
52 Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};
53
54 ClassImp(AliEMCALUnfolding)
55   
56 //____________________________________________________________________________
57 AliEMCALUnfolding::AliEMCALUnfolding():
58   fNumberOfECAClusters(0),
59   fECALocMaxCut(0),
60   fThreshold(0.01),//10 MeV
61   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 //____________________________________________________________________________
71 AliEMCALUnfolding::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 //____________________________________________________________________________
90 AliEMCALUnfolding::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++) fgSSPars[i] = SSPars[i];
107   for (i = 0; i < 3; i++) {
108     fgPar5[i] = Par5[i];
109     fgPar6[i] = Par6[i];
110   }
111
112 }
113
114 //____________________________________________________________________________
115 void 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 //____________________________________________________________________________
143 void 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 //____________________________________________________________________________
153 void 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 //____________________________________________________________________________
196 Bool_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 //____________________________________________________________________________
435 Bool_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 //____________________________________________________________________________
555 Bool_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 //____________________________________________________________________________
668 Double_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 = fgSSPars[7]*TMath::Sqrt(x*x+y*y);
675   Double_t rp1  = TMath::Power(r, fgSSPars[1]) ;
676   Double_t rp5  = TMath::Power(r, fgSSPars[5]) ;
677   Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;
678   return shape ;
679 }
680
681 //____________________________________________________________________________
682 void 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 = fgSSPars[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, fgSSPars[1]) ;
758             Double_t rp5  = TMath::Power(dr, fgSSPars[5]) ;
759             
760             Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * 
761             (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - 
762              (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + 
763               fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[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!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker
790       } // digit loop
791     } // recpoint, digits and geom not NULL
792   }// List is not NULL
793   
794 }
795
796
797 //____________________________________________________________________________
798 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){
799   for(UInt_t i=0;i<7;++i)
800     fgSSPars[i]=pars[i];
801   if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0
802 }
803
804 //____________________________________________________________________________
805 void AliEMCALUnfolding::SetPar5(Double_t *pars){
806   for(UInt_t i=0;i<3;++i)
807     fgPar5[i]=pars[i];
808 }
809
810 //____________________________________________________________________________
811 void AliEMCALUnfolding::SetPar6(Double_t *pars){
812   for(UInt_t i=0;i<3;++i)
813     fgPar6[i]=pars[i];
814 }
815
816 //____________________________________________________________________________
817 void 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   fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];
824 }
825
826 //____________________________________________________________________________
827 void 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   fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];
834 }
835
836 //____________________________________________________________________________
837 void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const 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