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