New Clusterizer that makes cluster grouping neighbour cells around highest energy...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.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 /* $Id$ */
17
18 //-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
19 //--         Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
20 //////////////////////////////////////////////////////////////////////////////
21 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
22 //  unfolds the clusters having several local maxima.  
23 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
24 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all 
25 //  parameters including input digits branch title, thresholds etc.)
26 //
27
28 // --- ROOT system ---
29
30 #include <TFile.h> 
31 #include <TMath.h> 
32 #include <TMinuit.h>
33 #include <TTree.h> 
34 #include <TBenchmark.h>
35 #include <TBrowser.h>
36 #include <TROOT.h>
37 #include <TList.h>
38 #include <TClonesArray.h>
39
40 // --- Standard library ---
41 #include <cassert>
42
43 // --- AliRoot header files ---
44 #include "AliLog.h"
45 #include "AliEMCALClusterizerv1.h"
46 #include "AliEMCALRecPoint.h"
47 #include "AliEMCALDigit.h"
48 #include "AliEMCALGeometry.h"
49 #include "AliCaloCalibPedestal.h"
50 #include "AliEMCALCalibData.h"
51 #include "AliESDCaloCluster.h"
52
53 ClassImp(AliEMCALClusterizerv1)
54
55 //____________________________________________________________________________
56 AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
57 {
58   // ctor with the indication of the file where header Tree and digits Tree are stored
59   
60   Init() ;
61 }
62
63 //____________________________________________________________________________
64 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
65   : AliEMCALClusterizer(geometry)
66 {
67   // ctor with the indication of the file where header Tree and digits Tree are stored
68   // use this contructor to avoid usage of Init() which uses runloader
69   // change needed by HLT - MP
70
71 }
72
73 //____________________________________________________________________________
74 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
75 : AliEMCALClusterizer(geometry, calib, caloped)
76 {
77         // ctor, geometry and calibration are initialized elsewhere.
78                                 
79 }
80
81
82 //____________________________________________________________________________
83   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
84 {
85   // dtor
86 }
87
88 //____________________________________________________________________________
89 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
90 {
91   // Steering method to perform clusterization for the current event 
92   // in AliEMCALLoader
93
94   if(strstr(option,"tim"))
95     gBenchmark->Start("EMCALClusterizer"); 
96   
97   if(strstr(option,"print"))
98     Print("") ; 
99  
100   //Get calibration parameters from file or digitizer default values.
101   GetCalibrationParameters() ;
102
103   //Get dead channel map from file or digitizer default values.
104   GetCaloCalibPedestal() ;
105         
106   fNumberOfECAClusters = 0;
107
108   MakeClusters() ;  //only the real clusters
109
110   if(fToUnfold)
111     MakeUnfolding() ;
112
113   Int_t index ;
114
115   //Evaluate position, dispersion and other RecPoint properties for EC section                      
116   for(index = 0; index < fRecPoints->GetEntries(); index++) {
117       dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
118           //For each rec.point set the distance to the nearest bad crystal
119           dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
120   }
121
122   fRecPoints->Sort() ;
123
124   for(index = 0; index < fRecPoints->GetEntries(); index++) {
125     (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
126     (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
127   }
128
129   fTreeR->Fill();
130   
131   if(strstr(option,"deb") || strstr(option,"all"))  
132     PrintRecPoints(option) ;
133
134   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
135
136   fRecPoints->Delete();
137
138   if(strstr(option,"tim")){
139     gBenchmark->Stop("EMCALClusterizer");
140     printf("Exec took %f seconds for Clusterizing", 
141            gBenchmark->GetCpuTime("EMCALClusterizer"));
142   }    
143 }
144
145 //____________________________________________________________________________
146 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, 
147                                       const Float_t* maxAtEnergy,
148                                       Int_t nPar, Float_t * fitparameters) const
149 {
150   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
151   // The initial values for fitting procedure are set equal to the
152   // positions of local maxima.       
153   // Cluster will be fitted as a superposition of nPar/3
154   // electromagnetic showers
155
156   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
157         
158   if(!gMinuit)
159      gMinuit = new TMinuit(100) ;
160
161   gMinuit->mncler();                     // Reset Minuit's list of paramters
162   gMinuit->SetPrintLevel(-1) ;           // No Printout
163   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
164   // To set the address of the minimization function
165   TList * toMinuit = new TList();
166   toMinuit->AddAt(recPoint,0) ;
167   toMinuit->AddAt(fDigitsArr,1) ;
168   toMinuit->AddAt(fGeom,2) ;
169
170   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
171
172   // filling initial values for fit parameters
173   AliEMCALDigit * digit ;
174
175   Int_t ierflg  = 0;
176   Int_t index   = 0 ;
177   Int_t nDigits = (Int_t) nPar / 3 ;
178
179   Int_t iDigit ;
180
181   for(iDigit = 0; iDigit < nDigits; iDigit++){
182     digit = maxAt[iDigit];
183     Double_t x = 0.;
184     Double_t y = 0.;
185     Double_t z = 0.;
186
187     fGeom->RelPosCellInSModule(digit->GetId(), y, x, z);
188
189     Float_t energy = maxAtEnergy[iDigit] ;
190
191     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
192     index++ ;
193     if(ierflg != 0){
194       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f", x ) ;
195       return kFALSE;
196     }
197     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
198     index++ ;
199     if(ierflg != 0){
200       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
201       return kFALSE;
202     }
203     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
204     index++ ;
205     if(ierflg != 0){
206       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
207       return kFALSE;
208     }
209   }
210
211   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; 
212                       // The number of function call slightly depends on it.
213   //Double_t p1 = 1.0 ;
214   Double_t p2 = 0.0 ;
215
216   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
217   //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
218   gMinuit->SetMaxIterations(5);
219   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
220   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
221
222   if(ierflg == 4){  // Minimum not found
223     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
224     return kFALSE ;
225   }
226   for(index = 0; index < nPar; index++){
227     Double_t err = 0. ;
228     Double_t val = 0. ;
229     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
230     fitparameters[index] = val ;
231   }
232
233   delete toMinuit ;
234   return kTRUE;
235
236 }
237
238 //____________________________________________________________________________
239 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
240 {
241         // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
242         //                                       = 1 are neighbour
243         //                                       = 2 is in different SM; continue searching 
244         // In case it is in different SM, but same phi rack, check if neigbours at eta=0
245         // neighbours are defined as digits having at least a common side 
246         // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
247         //                                      which is compared to a digit (d2)  not yet in a cluster  
248         
249         static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
250         static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
251
252         shared = kFALSE;
253         
254         fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
255         fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
256         fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
257         fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
258         
259         //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
260         if(nSupMod1 != nSupMod2 ) {
261                 //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
262                 Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
263                 Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
264                 
265                 if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
266                                 
267                 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
268                 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
269                 if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
270                 else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
271                 
272                 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
273                 
274         }//Different SM, same phi
275         
276         Int_t rowdiff = TMath::Abs(iphi1 - iphi2);  
277         Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;  
278
279         // neighbours with at least common side; May 11, 2007
280         if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {  
281         //Diagonal?
282         //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
283         
284         if (gDebug == 2) 
285                 printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
286                                 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);   
287         
288                 return 1;
289         }//Neighbours
290         else {
291                 shared = kFALSE;
292                 return 2 ; 
293         }//Not neighbours
294 }
295
296 //____________________________________________________________________________
297 void AliEMCALClusterizerv1::MakeClusters()
298 {
299   // Steering method to construct the clusters stored in a list of Reconstructed Points
300   // A cluster is defined as a list of neighbour digits
301   // Mar 03, 2007 by PAI
302
303   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
304
305   fRecPoints->Clear();
306
307   // Set up TObjArray with pointers to digits to work on 
308   TObjArray *digitsC = new TObjArray();
309   TIter nextdigit(fDigitsArr);
310   AliEMCALDigit *digit;
311   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
312     digitsC->AddLast(digit);
313   }
314
315   double e = 0.0, ehs = 0.0;
316   TIter nextdigitC(digitsC);
317   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
318     e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR?
319     if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate
320       digitsC->Remove(digit);
321     else    
322       ehs += e;
323   } 
324   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
325                   fDigitsArr->GetEntries(),fMinECut,ehs));
326
327   nextdigitC.Reset();
328
329   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
330     TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
331
332     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold  ) ){
333       // start a new Tower RecPoint
334       if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
335
336       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ; 
337       fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
338       recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ; 
339       fNumberOfECAClusters++ ; 
340
341       recPoint->SetClusterType(AliESDCaloCluster::kEMCALClusterv1);
342
343       recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE) ; //Time or TimeR?
344       TObjArray clusterDigits;
345       clusterDigits.AddLast(digit);     
346       digitsC->Remove(digit) ; 
347
348       AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
349       Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold));  //Time or TimeR?
350           Float_t time = digit->GetTime();//Time or TimeR?
351       // Grow cluster by finding neighbours
352       TIter nextClusterDigit(&clusterDigits);
353       while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster 
354         TIter nextdigitN(digitsC); 
355         AliEMCALDigit *digitN = 0; // digi neighbor
356         while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
357                         
358                         //Do not add digits with too different time 
359                         Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
360                         if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
361                         if (AreNeighbours(digit, digitN, shared)==1) {      // call (digit,digitN) in THAT order !!!!! 
362                                 recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared) ;//Time or TimeR?
363                                 clusterDigits.AddLast(digitN) ; 
364                                 digitsC->Remove(digitN) ; 
365                         } // if(ineb==1)
366                 } // scan over digits
367       } // scan over digits already in cluster
368         
369       if(recPoint)
370         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); 
371     } // If seed found
372   } // while digit 
373
374   delete digitsC ;
375   
376   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
377 }
378
379 //____________________________________________________________________________
380 void AliEMCALClusterizerv1::MakeUnfolding()
381 {
382   // Unfolds clusters using the shape of an ElectroMagnetic shower
383   // Performs unfolding of all clusters
384                 
385   if(fNumberOfECAClusters > 0){
386     if (fGeom==0)
387       AliFatal("Did not get geometry from EMCALLoader") ;
388     Int_t nModulesToUnfold = fGeom->GetNCells();
389
390     Int_t numberofNotUnfolded = fNumberOfECAClusters ;
391     Int_t index ;
392     for(index = 0 ; index < numberofNotUnfolded ; index++){
393
394       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
395
396       TVector3 gpos;
397       Int_t absId = -1;
398       recPoint->GetGlobalPosition(gpos);
399       fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
400       if(absId > nModulesToUnfold)
401         break ;
402
403       Int_t nMultipl = recPoint->GetMultiplicity() ;
404       AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
405       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
406       Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
407
408       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
409         UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
410         fRecPoints->Remove(recPoint);
411         fRecPoints->Compress() ;
412         index-- ;
413         fNumberOfECAClusters-- ;
414         numberofNotUnfolded-- ;
415       }
416       else{
417         recPoint->SetNExMax(1) ; //Only one local maximum
418       }
419
420       delete[] maxAt ;
421       delete[] maxAtEnergy ;
422     }
423   }
424   // End of Unfolding of clusters
425 }
426
427 //____________________________________________________________________________
428 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t x, Double_t y)
429
430   // Shape of the shower
431   // If you change this function, change also the gradient evaluation in ChiSquare()
432
433   Double_t r = sqrt(x*x+y*y);
434   Double_t r133  = TMath::Power(r, 1.33) ;
435   Double_t r669  = TMath::Power(r, 6.69) ;
436   Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
437   return shape ;
438 }
439
440 //____________________________________________________________________________
441 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, 
442                                            Int_t nMax, 
443                                            AliEMCALDigit ** maxAt, 
444                                            Float_t * maxAtEnergy)
445 {
446   // Performs the unfolding of a cluster with nMax overlapping showers 
447   Int_t nPar = 3 * nMax ;
448   Float_t * fitparameters = new Float_t[nPar] ;
449
450   if (fGeom==0)
451     AliFatal("Did not get geometry from EMCALLoader") ;
452
453   Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
454   if( !rv ) {
455     // Fit failed, return and remove cluster
456     iniTower->SetNExMax(-1) ;
457     delete[] fitparameters ;
458     return ;
459   }
460
461   // create unfolded rec points and fill them with new energy lists
462   // First calculate energy deposited in each sell in accordance with
463   // fit (without fluctuations): efit[]
464   // and later correct this number in acordance with actual energy
465   // deposition
466
467   Int_t nDigits = iniTower->GetMultiplicity() ;
468   Float_t * efit = new Float_t[nDigits] ;
469   Double_t xDigit=0.,yDigit=0.,zDigit=0. ;
470   Float_t xpar=0.,zpar=0.,epar=0.  ;
471
472   AliEMCALDigit * digit = 0 ;
473   Int_t * digitsList = iniTower->GetDigitsList() ;
474   
475   Int_t iparam = 0 ;
476   Int_t iDigit ;
477   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
478     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
479     fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
480     efit[iDigit] = 0;
481
482     while(iparam < nPar ){
483       xpar = fitparameters[iparam] ;
484       zpar = fitparameters[iparam+1] ;
485       epar = fitparameters[iparam+2] ;
486       iparam += 3 ;
487       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
488     }
489   }
490
491
492   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
493   // so that energy deposited in each cell is distributed between new clusters proportionally
494   // to its contribution to efit
495
496   Float_t * energiesList = iniTower->GetEnergiesList() ;
497   Float_t ratio = 0 ;
498
499   iparam = 0 ;
500   while(iparam < nPar ){
501     xpar = fitparameters[iparam] ;
502     zpar = fitparameters[iparam+1] ;
503     epar = fitparameters[iparam+2] ;
504     iparam += 3 ;
505
506     AliEMCALRecPoint * recPoint = 0 ;
507
508     if(fNumberOfECAClusters >= fRecPoints->GetSize())
509       fRecPoints->Expand(2*fNumberOfECAClusters) ;
510
511     (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
512     recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
513     fNumberOfECAClusters++ ;
514     recPoint->SetNExMax((Int_t)nPar/3) ;
515
516     Float_t eDigit = 0. ;
517     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
518       digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
519       fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
520
521       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
522       eDigit = energiesList[iDigit] * ratio ;
523       recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case
524     }
525   }
526
527   delete[] fitparameters ;
528   delete[] efit ;
529
530 }
531
532 //_____________________________________________________________________________
533 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
534                                                Double_t & fret,
535                                                Double_t * x, Int_t iflag)
536 {
537   // Calculates the Chi square for the cluster unfolding minimization
538   // Number of parameters, Gradient, Chi squared, parameters, what to do
539
540   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
541
542   AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;
543   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
544   // A bit buggy way to get an access to the geometry
545   // To be revised!
546   AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
547
548   Int_t * digitsList     = recPoint->GetDigitsList() ;
549
550   Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
551
552   Float_t * energiesList = recPoint->GetEnergiesList() ;
553
554   fret = 0. ;
555   Int_t iparam ;
556
557   if(iflag == 2)
558     for(iparam = 0 ; iparam < nPar ; iparam++)
559       Grad[iparam] = 0 ; // Will evaluate gradient
560
561   Double_t efit = 0. ;
562
563   AliEMCALDigit * digit ;
564   Int_t iDigit ;
565
566   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
567
568     digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
569
570     Double_t xDigit=0 ;
571     Double_t zDigit=0 ;
572     Double_t yDigit=0 ;//not used yet, assumed to be 0
573
574     geom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
575
576     if(iflag == 2){  // calculate gradient
577       Int_t iParam = 0 ;
578       efit = 0. ;
579       while(iParam < nPar ){
580         Double_t dx = (xDigit - x[iParam]) ;
581         iParam++ ;
582         Double_t dz = (zDigit - x[iParam]) ;
583         iParam++ ;
584         efit += x[iParam] * ShowerShape(dx,dz) ;
585         iParam++ ;
586       }
587       Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
588       iParam = 0 ;
589       while(iParam < nPar ){
590         Double_t xpar = x[iParam] ;
591         Double_t zpar = x[iParam+1] ;
592         Double_t epar = x[iParam+2] ;
593         Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
594         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
595         Double_t r133 =  TMath::Power(dr, 1.33);
596         Double_t r669 = TMath::Power(dr,6.69);
597         Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
598                                                              - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
599
600         Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
601         iParam++ ;
602         Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
603         iParam++ ;
604         Grad[iParam] += shape ;                                  // Derivative over energy
605         iParam++ ;
606       }
607     }
608     efit = 0;
609     iparam = 0 ;
610
611
612     while(iparam < nPar ){
613       Double_t xpar = x[iparam] ;
614       Double_t zpar = x[iparam+1] ;
615       Double_t epar = x[iparam+2] ;
616       iparam += 3 ;
617       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
618     }
619
620     fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
621     // Here we assume, that sigma = sqrt(E) 
622   }
623 }