Line data Source code
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 42 : ClassImp(AliEMCALUnfolding)
55 :
56 : //____________________________________________________________________________
57 0 : AliEMCALUnfolding::AliEMCALUnfolding():
58 0 : fNumberOfECAClusters(0),
59 0 : fECALocMaxCut(0),
60 0 : fThreshold(0.01),//10 MeV
61 0 : fRejectBelowThreshold(0),//split
62 0 : fGeom(NULL),
63 0 : fRecPoints(NULL),
64 0 : fDigitsArr(NULL)
65 0 : {
66 : // ctor with the indication of the file where header Tree and digits Tree are stored
67 0 : Init() ;
68 0 : }
69 :
70 : //____________________________________________________________________________
71 0 : AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):
72 0 : fNumberOfECAClusters(0),
73 0 : fECALocMaxCut(0),
74 0 : fThreshold(0.01),//10 MeV
75 0 : fRejectBelowThreshold(0),//split
76 0 : fGeom(geometry),
77 0 : fRecPoints(NULL),
78 0 : fDigitsArr(NULL)
79 0 : {
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 0 : if (!fGeom)
84 : {
85 0 : AliFatal("AliEMCALUnfolding: Geometry not initialized.");
86 : }
87 :
88 0 : }
89 :
90 : //____________________________________________________________________________
91 0 : AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):
92 0 : fNumberOfECAClusters(0),
93 0 : fECALocMaxCut(ECALocMaxCut),
94 0 : fThreshold(0.01),//10 MeV
95 0 : fRejectBelowThreshold(0),//split
96 0 : fGeom(geometry),
97 0 : fRecPoints(NULL),
98 0 : fDigitsArr(NULL)
99 0 : {
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 0 : if (!fGeom)
104 : {
105 0 : AliFatal("AliEMCALUnfolding: Geometry not initialized.");
106 : }
107 : Int_t i=0;
108 0 : for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];
109 0 : for (i = 0; i < 3; i++) {
110 0 : fgPar5[i] = Par5[i];
111 0 : fgPar6[i] = Par6[i];
112 : }
113 :
114 0 : }
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 0 : AliRunLoader *rl = AliRunLoader::Instance();
123 0 : if (rl && rl->GetAliRun()){
124 0 : AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
125 0 : if(emcal)fGeom = emcal->GetGeometry();
126 0 : }
127 :
128 0 : if(!fGeom)
129 0 : fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
130 :
131 0 : AliDebug(1,Form("geom %p",fGeom));
132 :
133 0 : if(!gMinuit)
134 : // gMinuit = new TMinuit(100) ;//the same is in FindFitV2
135 0 : gMinuit = new TMinuit(30) ;//the same is in FindFitV2
136 :
137 0 : }
138 :
139 : //____________________________________________________________________________
140 : AliEMCALUnfolding::~AliEMCALUnfolding()
141 0 : {
142 : // dtor
143 0 : }
144 :
145 : //____________________________________________________________________________
146 : void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)
147 : {
148 : //
149 : //Set input for unfolding purposes
150 : //
151 0 : SetNumberOfECAClusters(numberOfECAClusters);
152 0 : SetRecPoints(recPoints);
153 0 : SetDigitsArr(digitsArr);
154 0 : }
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 0 : AliDebug(4,Form(" V1: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
163 0 : if(fNumberOfECAClusters > 0){
164 0 : if (fGeom==0)
165 0 : AliFatal("Did not get geometry from EMCALLoader") ;
166 : //Int_t nModulesToUnfold = fGeom->GetNCells();
167 :
168 0 : 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 0 : for(index = 0 ; index < numberOfClustersToUnfold ; index++){
176 0 : AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
177 0 : if(recPoint){
178 0 : Int_t nMultipl = recPoint->GetMultiplicity() ;
179 0 : AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
180 0 : Float_t * maxAtEnergy = new Float_t[nMultipl] ;
181 0 : Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
182 0 : if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
183 0 : AliDebug(4,Form(" *** V1+UNFOLD *** Cluster index before UF %d",fNumberOfECAClusters));
184 0 : if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){
185 : //if unfolding correct remove old recPoint
186 0 : fRecPoints->Remove(recPoint);
187 0 : fRecPoints->Compress() ;//is it really needed
188 0 : index-- ;
189 0 : fNumberOfECAClusters-- ;
190 0 : numberOfClustersToUnfold--;
191 0 : }
192 0 : AliDebug(4,Form(" Cluster index after UF %d",fNumberOfECAClusters));
193 : } else{
194 0 : recPoint->SetNExMax(1) ; //Only one local maximum
195 : }
196 :
197 0 : delete[] maxAt ;
198 0 : delete[] maxAtEnergy ;
199 0 : } else {
200 : //AliError("RecPoint NULL"); //end of check if recPoint exist
201 0 : Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ;
202 : }
203 : } // rec point loop
204 0 : }//end of check fNumberOfECAClusters
205 : // End of Unfolding of clusters
206 :
207 0 : 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 0 : }
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 0 : AliDebug(5,Form(" Original cluster E %f, nMax = %d",iniTower->GetEnergy(),nMax ));
239 :
240 0 : Int_t nPar = 3 * nMax ;
241 0 : Float_t * fitparameters = new Float_t[nPar] ;
242 :
243 0 : if (fGeom==0)
244 0 : AliFatal("Did not get geometry from EMCALLoader") ;
245 :
246 0 : Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
247 0 : if( !rv )
248 : {
249 : // Fit failed, return (and remove cluster? - why? I leave the cluster)
250 0 : iniTower->SetNExMax(-1) ;
251 0 : delete[] fitparameters ;
252 0 : return 0;//changed here
253 : }
254 :
255 : //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold
256 0 : if(nMax==2){
257 0 : if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){
258 0 : AliDebug(1,"One of fitted energy below threshold");
259 0 : iniTower->SetNExMax(1) ;
260 0 : delete[] fitparameters ;
261 0 : 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 0 : Int_t nDigits = iniTower->GetMultiplicity() ;
274 0 : 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 0 : Int_t * digitsList = iniTower->GetDigitsList() ;
279 :
280 0 : Int_t iSupMod = 0 ;
281 0 : Int_t iTower = 0 ;
282 0 : Int_t iIphi = 0 ;
283 0 : Int_t iIeta = 0 ;
284 0 : Int_t iphi = 0 ;//x direction
285 0 : Int_t ieta = 0 ;//z direstion
286 :
287 : Int_t iparam = 0 ;
288 : Int_t iDigit = 0 ;
289 :
290 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
291 : {
292 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
293 0 : if(digit)
294 : {
295 0 : fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
296 0 : fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
297 0 : iIphi, iIeta,iphi,ieta);
298 0 : EvalParsPhiDependence(digit->GetId(),fGeom);
299 :
300 0 : efit[iDigit] = 0.;
301 : iparam = 0;
302 0 : while(iparam < nPar )
303 : {
304 0 : xpar = fitparameters[iparam] ;
305 0 : zpar = fitparameters[iparam+1] ;
306 0 : epar = fitparameters[iparam+2] ;
307 :
308 0 : efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
309 0 : iparam += 3 ;
310 : }
311 :
312 0 : } 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 0 : 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 0 : 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 0 : while(iparam < nPar )
340 : {
341 0 : xpar = fitparameters[iparam] ;
342 0 : zpar = fitparameters[iparam+1] ;
343 0 : epar = fitparameters[iparam+2] ;
344 :
345 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
346 : {
347 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
348 0 : if(digit)
349 : {
350 0 : fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
351 0 : fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
352 0 : iIphi, iIeta,iphi,ieta);
353 :
354 0 : EvalParsPhiDependence(digit->GetId(),fGeom);
355 :
356 0 : if(efit[iDigit]==0)
357 : {//just for sure
358 0 : correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here
359 0 : continue;
360 : }
361 :
362 0 : ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
363 0 : eDigit = energiesList[iDigit] * ratio ;
364 :
365 : //add energy to temporary matrix
366 0 : correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;
367 :
368 0 : } else AliDebug(1,"NULL digit part 3");
369 : }//digit loop
370 0 : 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 0 : if(fThreshold > 0){//option 2 or 3
387 0 : if(fRejectBelowThreshold){//option 3
388 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit++){//digit loop
389 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3){//param0 loop = energy loop
390 0 : 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 0 : Float_t * energyFraction = new Float_t[nSplittedClusters];
399 : Int_t iparam2 = 0 ;
400 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit++){
401 : isAnyBelowThreshold=kFALSE;
402 : maximumEne=0.;
403 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3){
404 0 : if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;
405 0 : if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne)
406 : {
407 : maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];
408 : maximumIndex = iparam;
409 0 : }
410 : }//end of loop over clusters after unfolding
411 :
412 0 : if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold
413 :
414 0 : if(maximumEne < fThreshold)
415 : {//add all cluster cells and put energy into max index, other set to 0
416 : maximumEne=0.;
417 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3)
418 : {
419 0 : maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
420 0 : correctedEnergyList[iparam/3*nDigits+iDigit]=0;
421 : }
422 0 : correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;
423 0 : 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 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3)
429 : {//calculate energy sum
430 0 : if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;
431 : else
432 : {
433 0 : energyFraction[iparam/3]=1.;
434 0 : maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
435 : }
436 : }//end of loop over clusters after unfolding
437 0 : if(maximumEne>0.) {
438 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction
439 0 : energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;
440 : }
441 :
442 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3)
443 : {//add energy from cells below threshold to others
444 0 : if(energyFraction[iparam/3]>0.) continue;
445 : else
446 : {
447 0 : for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)
448 : {
449 0 : correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] *
450 0 : correctedEnergyList[iparam/3*nDigits+iDigit]) ;
451 : }//inner loop
452 0 : correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
453 : }
454 0 : }
455 : } else {
456 : //digit energy to be set to 0
457 0 : for(iparam = 0 ; iparam < nPar ; iparam+=3)
458 : {
459 0 : correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
460 : }
461 : }//correction for: is energy>0
462 :
463 : }//end of loop over digits
464 0 : 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 0 : while(iparam < nPar )
479 : {
480 : AliEMCALRecPoint * recPoint = 0 ;
481 :
482 0 : if(nSplittedClusters >= list->GetSize())
483 0 : list->Expand(nSplittedClusters);
484 :
485 : //add recpoint
486 0 : (*list)[newClusterIndex] = new AliEMCALRecPoint("") ;
487 0 : recPoint = dynamic_cast<AliEMCALRecPoint *>( list->At(newClusterIndex) ) ;
488 :
489 0 : if(recPoint){//recPoint present -> good
490 0 : recPoint->SetNExMax(nSplittedClusters) ;//can be wrong number, to be corrected in outer method
491 :
492 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) {
493 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
494 0 : 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 0 : recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case
497 0 : } else {
498 0 : AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));
499 : }
500 : }//digit loop
501 :
502 0 : if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list
503 0 : delete (*list)[newClusterIndex];
504 0 : list->RemoveAt(newClusterIndex);
505 0 : nSplittedClusters--;
506 0 : newClusterIndex--;//decrease cluster number
507 0 : }else {//recPoint exists and has digits associated -> very good increase number of clusters
508 0 : 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 0 : AliError("NULL RecPoint");
513 : //protection from recpoint with no digits
514 0 : delete (*list)[newClusterIndex];
515 0 : list->RemoveAt(newClusterIndex);
516 0 : nSplittedClusters--;
517 0 : newClusterIndex--;//decrease cluster number
518 : }
519 :
520 0 : iparam += 3 ;
521 0 : newClusterIndex++;
522 : }//while
523 :
524 0 : delete[] fitparameters ;
525 0 : delete[] efit ;
526 0 : delete[] correctedEnergyList ;
527 :
528 : // print
529 0 : 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 0 : }
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 0 : TObjArray *list =new TObjArray(2);//temporary object
546 0 : Int_t nUnfoldedClusters=UnfoldOneCluster(iniTower,nMax,maxAt,maxAtEnergy,list);
547 :
548 : // here we write new clusters from list to fRecPoints
549 0 : AliDebug(5,Form("Number of clusters after unfolding %d",list->GetEntriesFast()));
550 : Int_t iDigit=0;
551 : AliEMCALDigit * digit = 0 ;
552 0 : for(Int_t i=0;i<list->GetEntriesFast();i++) {
553 : AliEMCALRecPoint * recPoint = 0 ;
554 :
555 0 : if(fNumberOfECAClusters >= fRecPoints->GetSize())
556 0 : fRecPoints->Expand(2*fNumberOfECAClusters) ;
557 :
558 : //add recpoint
559 0 : (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;//fNumberOfECAClusters-1 is old cluster before unfolding
560 0 : recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
561 0 : AliEMCALRecPoint * rpUFOne = dynamic_cast<AliEMCALRecPoint *>(list->At(i)) ;
562 :
563 0 : if( recPoint && rpUFOne ){//recPoint present -> good
564 :
565 0 : recPoint->SetNExMax(list->GetEntriesFast()) ;
566 :
567 0 : Int_t *digitsList = rpUFOne->GetDigitsList();
568 0 : Float_t *energyList = rpUFOne->GetEnergiesList();
569 :
570 0 : if(!digitsList || ! energyList)
571 : {
572 0 : AliDebug(-1,"No digits index or energy available");
573 0 : delete (*fRecPoints)[fNumberOfECAClusters];
574 0 : fRecPoints->RemoveAt(fNumberOfECAClusters);
575 0 : continue;
576 : }
577 :
578 0 : AliDebug(5,Form("cluster %d, digit no %d, energy %f\n",i,digitsList[0],energyList[0]));
579 :
580 0 : for(iDigit = 0 ; iDigit < rpUFOne->GetMultiplicity(); iDigit ++) {
581 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
582 0 : if(digit) recPoint->AddDigit( *digit, energyList[iDigit], kFALSE ) ; //FIXME, need to study the shared case
583 : }//digit loop
584 0 : fNumberOfECAClusters++ ;
585 0 : } else {//recPoint empty -> remove from list
586 0 : AliError("NULL RecPoint");
587 0 : delete (*fRecPoints)[fNumberOfECAClusters];
588 0 : fRecPoints->RemoveAt(fNumberOfECAClusters);
589 : }
590 :
591 0 : }//loop over unfolded clusters
592 :
593 : //print energy of new unfolded clusters
594 0 : AliDebug(5,Form(" nUnfoldedClusters %d, fNumberOfECAClusters %d",nUnfoldedClusters,fNumberOfECAClusters ));
595 0 : for(Int_t inewclus=0; inewclus<nUnfoldedClusters;inewclus++){
596 0 : AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters-1-inewclus));
597 0 : if(rp) AliDebug(5,Form(" Unfolded cluster %d E %f",inewclus, rp->GetEnergy() ));
598 : }
599 :
600 : //clear tables
601 0 : list->SetOwner(kTRUE);
602 0 : list->Delete();
603 0 : delete list;
604 0 : if(nUnfoldedClusters>1) return kTRUE;
605 0 : return kFALSE;
606 0 : }
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 0 : Int_t nPar = 3 * nMax ;
620 0 : Float_t * fitparameters = new Float_t[nPar] ;
621 :
622 0 : if (fGeom==0)
623 0 : AliFatal("Did not get geometry from EMCALLoader") ;
624 :
625 0 : Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
626 0 : if( !rv ) {
627 : // Fit failed, return (and remove cluster? - why? I leave the cluster)
628 0 : iniTower->SetNExMax(-1) ;
629 0 : delete[] fitparameters ;
630 0 : 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 0 : Int_t nDigits = iniTower->GetMultiplicity() ;
640 0 : 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 0 : Int_t * digitsList = iniTower->GetDigitsList() ;
645 :
646 0 : Int_t iSupMod = 0 ;
647 0 : Int_t iTower = 0 ;
648 0 : Int_t iIphi = 0 ;
649 0 : Int_t iIeta = 0 ;
650 0 : Int_t iphi = 0 ;//x direction
651 0 : Int_t ieta = 0 ;//z direstion
652 :
653 : Int_t iparam = 0 ;
654 : Int_t iDigit = 0 ;
655 :
656 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
657 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
658 0 : if(digit){
659 0 : fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
660 0 : fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
661 0 : iIphi, iIeta,iphi,ieta);
662 0 : EvalParsPhiDependence(digit->GetId(),fGeom);
663 :
664 0 : efit[iDigit] = 0.;
665 : iparam = 0;
666 0 : while(iparam < nPar ){
667 0 : xpar = fitparameters[iparam] ;
668 0 : zpar = fitparameters[iparam+1] ;
669 0 : epar = fitparameters[iparam+2] ;
670 0 : iparam += 3 ;
671 :
672 0 : efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
673 : }
674 0 : } 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 0 : Float_t * energiesList = iniTower->GetEnergiesList() ;
683 : Float_t ratio = 0 ;
684 :
685 : iparam = 0 ;
686 0 : while(iparam < nPar ){
687 0 : xpar = fitparameters[iparam] ;
688 0 : zpar = fitparameters[iparam+1] ;
689 0 : epar = fitparameters[iparam+2] ;
690 0 : iparam += 3 ;
691 :
692 : AliEMCALRecPoint * recPoint = 0 ;
693 :
694 0 : if(fNumberOfECAClusters >= fRecPoints->GetSize())
695 0 : fRecPoints->Expand(2*fNumberOfECAClusters) ;
696 :
697 : //add recpoint
698 0 : (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
699 0 : recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
700 :
701 0 : if(recPoint){
702 :
703 0 : fNumberOfECAClusters++ ;
704 0 : recPoint->SetNExMax((Int_t)nPar/3) ;
705 :
706 : Float_t eDigit = 0. ;
707 0 : for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
708 0 : digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
709 0 : if(digit){
710 0 : fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
711 0 : fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
712 0 : iIphi, iIeta,iphi,ieta);
713 0 : EvalParsPhiDependence(digit->GetId(),fGeom);
714 0 : if(efit[iDigit]==0) continue;//just for sure
715 0 : ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
716 0 : eDigit = energiesList[iDigit] * ratio ;
717 0 : recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case
718 0 : } else AliError("NULL digit");
719 : }//digit loop
720 0 : } else AliError("NULL RecPoint");
721 : }//while
722 :
723 0 : delete[] fitparameters ;
724 0 : delete[] efit ;
725 :
726 : return kTRUE;
727 0 : }
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 0 : if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
742 :
743 0 : if(!gMinuit){
744 : // gMinuit = new TMinuit(100) ;//max 100 parameters
745 0 : if(nPar<30) gMinuit = new TMinuit(30);
746 0 : else gMinuit = new TMinuit(nPar) ;//max nPar parameters
747 : //
748 : } else {
749 0 : if(gMinuit->fMaxpar < nPar) {
750 0 : delete gMinuit;
751 0 : gMinuit = new TMinuit(nPar);
752 0 : }
753 : }
754 :
755 0 : gMinuit->mncler(); // Reset Minuit's list of paramters
756 0 : gMinuit->SetPrintLevel(-1) ; // No Printout
757 0 : gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;
758 : // To set the address of the minimization function
759 0 : TList * toMinuit = new TList();
760 0 : toMinuit->AddAt(recPoint,0) ;
761 0 : toMinuit->AddAt(fDigitsArr,1) ;
762 0 : toMinuit->AddAt(fGeom,2) ;
763 :
764 0 : gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
765 :
766 : // filling initial values for fit parameters
767 : AliEMCALDigit * digit ;
768 :
769 0 : Int_t ierflg = 0;
770 : Int_t index = 0 ;
771 0 : Int_t nDigits = (Int_t) nPar / 3 ;
772 :
773 : Int_t iDigit ;
774 :
775 0 : Int_t iSupMod = 0 ;
776 0 : Int_t iTower = 0 ;
777 0 : Int_t iIphi = 0 ;
778 0 : Int_t iIeta = 0 ;
779 0 : Int_t iphi = 0 ;//x direction
780 0 : Int_t ieta = 0 ;//z direstion
781 :
782 0 : for(iDigit = 0; iDigit < nDigits; iDigit++)
783 : {
784 0 : digit = maxAt[iDigit];
785 0 : if(!digit)
786 : {
787 0 : AliError("Null digit pointer");
788 0 : continue;
789 : }
790 :
791 0 : fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
792 0 : fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
793 0 : iIphi, iIeta,iphi,ieta);
794 :
795 0 : Float_t energy = maxAtEnergy[iDigit] ;
796 :
797 : //gMinuit->mnparm(index, "x", iphi, 0.1, 0, 0, ierflg) ;//original
798 0 : gMinuit->mnparm(index, "x", iphi, 0.05, 0, 0, ierflg) ;
799 0 : index++ ;
800 0 : if(ierflg != 0){
801 0 : Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: x=%d, param.id=%d, nMaxima=%d",iphi,index-1,nPar/3 ) ;
802 0 : toMinuit->Clear();
803 0 : delete toMinuit ;
804 0 : return kFALSE;
805 : }
806 : //gMinuit->mnparm(index, "z", ieta, 0.1, 0, 0, ierflg) ;//original
807 0 : gMinuit->mnparm(index, "z", ieta, 0.05, 0, 0, ierflg) ;
808 0 : index++ ;
809 0 : if(ierflg != 0){
810 0 : Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: z=%d, param.id=%d, nMaxima=%d", ieta, index-1,nPar/3) ;
811 0 : toMinuit->Clear();
812 0 : delete toMinuit ;
813 0 : return kFALSE;
814 : }
815 : //gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original
816 0 : gMinuit->mnparm(index, "Energy", energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05
817 0 : index++ ;
818 0 : if(ierflg != 0){
819 0 : Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: energy = %f, param.id=%d, nMaxima=%d", energy, index-1, nPar/3) ;
820 0 : toMinuit->Clear();
821 0 : delete toMinuit ;
822 0 : return kFALSE;
823 : }
824 0 : }
825 :
826 0 : Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ;
827 : // The number of function call slightly depends on it.
828 : // Double_t p1 = 1.0 ;// par to gradient
829 0 : Double_t p2 = 0.0 ;
830 : // Double_t p3 = 3.0 ;
831 0 : gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
832 : // gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
833 0 : gMinuit->SetMaxIterations(5);//was 5
834 0 : gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
835 : //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ; // printouts
836 :
837 0 : gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
838 : //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ; // minimize
839 0 : if(ierflg == 4){ // Minimum not found
840 0 : AliDebug(1,"EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
841 0 : toMinuit->Clear();
842 0 : delete toMinuit ;
843 0 : return kFALSE ;
844 : }
845 0 : for(index = 0; index < nPar; index++){
846 0 : Double_t err = 0. ;
847 0 : Double_t val = 0. ;
848 0 : gMinuit->GetParameter(index, val, err) ; // Returns value and error ofOA parameter index
849 0 : fitparameters[index] = val ;
850 0 : }
851 :
852 0 : toMinuit->Clear();
853 0 : delete toMinuit ;
854 :
855 0 : if(gMinuit->fMaxpar>30) delete gMinuit;
856 :
857 0 : return kTRUE;
858 :
859 0 : }
860 :
861 : //____________________________________________________________________________
862 : Double_t AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)
863 : {
864 : // extended to whole EMCAL
865 : // Shape of the shower
866 : // If you change this function, change also the gradient evaluation in ChiSquare()
867 :
868 0 : Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);
869 0 : Double_t rp1 = TMath::Power(r, fgSSPars[1]) ;
870 0 : Double_t rp5 = TMath::Power(r, fgSSPars[5]) ;
871 0 : Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;
872 0 : return shape ;
873 : }
874 :
875 : //____________________________________________________________________________
876 : void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
877 : Double_t & fret,
878 : Double_t * x, Int_t iflag)
879 : {
880 : // Calculates the Chi square for the cluster unfolding minimization
881 : // Number of parameters, Gradient, Chi squared, parameters, what to do
882 :
883 0 : TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
884 0 : if(toMinuit){
885 0 : AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;
886 0 : TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
887 : // A bit buggy way to get an access to the geometry
888 : // To be revised!
889 0 : AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
890 :
891 0 : if(recPoint && digits && geom){
892 :
893 0 : Int_t * digitsList = recPoint->GetDigitsList() ;
894 :
895 0 : Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
896 :
897 0 : Float_t * energiesList = recPoint->GetEnergiesList() ;
898 :
899 0 : fret = 0. ;
900 : Int_t iparam = 0 ;
901 :
902 0 : if(iflag == 2)
903 0 : for(iparam = 0 ; iparam < nPar ; iparam++)
904 0 : Grad[iparam] = 0 ; // Will evaluate gradient
905 :
906 : Double_t efit = 0. ;
907 :
908 : AliEMCALDigit * digit ;
909 : Int_t iDigit ;
910 :
911 0 : Int_t iSupMod = 0 ;
912 0 : Int_t iTower = 0 ;
913 0 : Int_t iIphi = 0 ;
914 0 : Int_t iIeta = 0 ;
915 0 : Int_t iphi = 0 ;//x direction
916 0 : Int_t ieta = 0 ;//z direstion
917 :
918 :
919 0 : for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
920 0 : if(energiesList[iDigit]==0) continue;
921 :
922 0 : digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
923 :
924 0 : if(digit){
925 0 : geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
926 0 : geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
927 0 : iIphi, iIeta,iphi,ieta);
928 0 : EvalParsPhiDependence(digit->GetId(),geom);
929 :
930 0 : if(iflag == 2){ // calculate gradient
931 : Int_t iParam = 0 ;
932 : efit = 0. ;
933 0 : while(iParam < nPar ){
934 0 : Double_t dx = ((Float_t)iphi - x[iParam]) ;
935 0 : iParam++ ;
936 0 : Double_t dz = ((Float_t)ieta - x[iParam]) ;
937 0 : iParam++ ;
938 0 : efit += x[iParam] * ShowerShapeV2(dx,dz) ;
939 0 : iParam++ ;
940 : }
941 :
942 0 : Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
943 : iParam = 0 ;
944 0 : while(iParam < nPar ){
945 0 : Double_t xpar = x[iParam] ;
946 0 : Double_t zpar = x[iParam+1] ;
947 0 : Double_t epar = x[iParam+2] ;
948 :
949 0 : Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );
950 0 : Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
951 0 : Double_t rp1 = TMath::Power(dr, fgSSPars[1]) ;
952 0 : Double_t rp5 = TMath::Power(dr, fgSSPars[5]) ;
953 :
954 0 : Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] *
955 0 : (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) -
956 0 : (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) +
957 0 : fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );
958 :
959 : //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
960 : // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
961 :
962 0 : Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x
963 : iParam++ ;
964 0 : Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z
965 0 : iParam++ ;
966 0 : Grad[iParam] += shape ; // Derivative over energy
967 0 : iParam++ ;
968 : }
969 0 : }
970 : efit = 0;
971 : iparam = 0 ;
972 :
973 0 : while(iparam < nPar ){
974 0 : Double_t xpar = x[iparam] ;
975 0 : Double_t zpar = x[iparam+1] ;
976 0 : Double_t epar = x[iparam+2] ;
977 0 : iparam += 3 ;
978 0 : efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
979 : }
980 :
981 0 : fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
982 : // Here we assume, that sigma = sqrt(E)
983 0 : } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker
984 : } // digit loop
985 0 : } // recpoint, digits and geom not NULL
986 0 : }// List is not NULL
987 :
988 0 : }
989 :
990 :
991 : //____________________________________________________________________________
992 : void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){
993 0 : for(UInt_t i=0;i<7;++i)
994 0 : fgSSPars[i]=pars[i];
995 0 : if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0
996 0 : }
997 :
998 : //____________________________________________________________________________
999 : void AliEMCALUnfolding::SetPar5(Double_t *pars){
1000 0 : for(UInt_t i=0;i<3;++i)
1001 0 : fgPar5[i]=pars[i];
1002 0 : }
1003 :
1004 : //____________________________________________________________________________
1005 : void AliEMCALUnfolding::SetPar6(Double_t *pars){
1006 0 : for(UInt_t i=0;i<3;++i)
1007 0 : fgPar6[i]=pars[i];
1008 0 : }
1009 :
1010 : //____________________________________________________________________________
1011 : void AliEMCALUnfolding::EvalPar5(Double_t phi){
1012 : //
1013 : //Evaluate the 5th parameter of the shower shape function
1014 : //phi in degrees range (-10,10)
1015 : //
1016 : //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;
1017 0 : fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];
1018 0 : }
1019 :
1020 : //____________________________________________________________________________
1021 : void AliEMCALUnfolding::EvalPar6(Double_t phi){
1022 : //
1023 : //Evaluate the 6th parameter of the shower shape function
1024 : //phi in degrees range (-10,10)
1025 : //
1026 : //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;
1027 0 : fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];
1028 0 : }
1029 :
1030 : //____________________________________________________________________________
1031 : void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){
1032 : //
1033 : // calculate params p5 and p6 depending on the phi angle in global coordinate
1034 : // for the cell with given absId index
1035 : // output phiglob should be in range -10 to 10 degree
1036 : //
1037 0 : Double_t etaGlob = 0.;//eta in global c.s. - unused
1038 0 : Double_t phiGlob = 0.;//phi in global c.s. in radians
1039 0 : geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);
1040 0 : if(phiGlob<0)phiGlob+=TMath::TwoPi();
1041 0 : phiGlob*=180./TMath::Pi();
1042 0 : phiGlob-=20.*(Int_t)(phiGlob/20.);
1043 0 : Int_t superModule=geom->GetSuperModuleNumber(absId);
1044 0 : if(superModule==10 || superModule==11 || superModule==18 || superModule==19) //sm 10,11,18,19 shift by 3.5 deg other 10 deg
1045 0 : phiGlob-=3.5;
1046 : else
1047 0 : phiGlob-=10.;
1048 :
1049 0 : EvalPar5(phiGlob);
1050 0 : EvalPar6(phiGlob);
1051 0 : }
1052 :
|
