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