1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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)
23 // unfolds the clusters having several local maxima.
24 //////////////////////////////////////////////////////////////////////////////
26 // --- ROOT system ---
27 #include "TClonesArray.h"
31 // --- Standard library ---
34 // --- AliRoot header files ---
35 #include "AliEMCALUnfolding.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliRunLoader.h"
40 #include "AliEMCALRecParam.h"
41 #include "AliEMCALRecPoint.h"
42 #include "AliEMCALDigit.h"
43 #include "AliEMCALReconstructor.h"
46 #include "AliCDBManager.h"
48 #include "AliCDBEntry.h"
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};
54 ClassImp(AliEMCALUnfolding)
56 //____________________________________________________________________________
57 AliEMCALUnfolding::AliEMCALUnfolding():
58 fNumberOfECAClusters(0),
60 fThreshold(0.01),//10 MeV
61 fRejectBelowThreshold(0),//split
66 // ctor with the indication of the file where header Tree and digits Tree are stored
70 //____________________________________________________________________________
71 AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):
72 fNumberOfECAClusters(0),
74 fThreshold(0.01),//10 MeV
75 fRejectBelowThreshold(0),//split
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
85 AliFatal("AliEMCALUnfolding: Geometry not initialized.");
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
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
105 AliFatal("AliEMCALUnfolding: Geometry not initialized.");
108 for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];
109 for (i = 0; i < 3; i++) {
116 //____________________________________________________________________________
117 void AliEMCALUnfolding::Init()
119 // Make all memory allocations which can not be done in default constructor.
120 // Attach the Clusterizer task to the list of EMCAL tasks
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();
129 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
131 AliDebug(1,Form("geom %p",fGeom));
134 // gMinuit = new TMinuit(100) ;//the same is in FindFitV2
135 gMinuit = new TMinuit(30) ;//the same is in FindFitV2
139 //____________________________________________________________________________
140 AliEMCALUnfolding::~AliEMCALUnfolding()
145 //____________________________________________________________________________
146 void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)
149 //Set input for unfolding purposes
151 SetNumberOfECAClusters(numberOfECAClusters);
152 SetRecPoints(recPoints);
153 SetDigitsArr(digitsArr);
156 //____________________________________________________________________________
157 void AliEMCALUnfolding::MakeUnfolding()
159 // Unfolds clusters using the shape of an ElectroMagnetic shower
160 // Performs unfolding of all clusters
162 AliDebug(4,Form(" V1: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
163 if(fNumberOfECAClusters > 0){
165 AliFatal("Did not get geometry from EMCALLoader") ;
166 //Int_t nModulesToUnfold = fGeom->GetNCells();
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
175 for(index = 0 ; index < numberOfClustersToUnfold ; index++){
176 AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
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
189 fNumberOfECAClusters-- ;
190 numberOfClustersToUnfold--;
192 AliDebug(4,Form(" Cluster index after UF %d",fNumberOfECAClusters));
194 recPoint->SetNExMax(1) ; //Only one local maximum
198 delete[] maxAtEnergy ;
200 //AliError("RecPoint NULL"); //end of check if recPoint exist
201 Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ;
204 }//end of check fNumberOfECAClusters
205 // End of Unfolding of clusters
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);
222 //____________________________________________________________________________
223 Int_t AliEMCALUnfolding::UnfoldOneCluster(AliEMCALRecPoint * iniTower,
225 AliEMCALDigit ** maxAt,
226 Float_t * maxAtEnergy,
230 // Output list of clusters
231 // returns number of clusters
232 // if fit failed or unfolding is not applicable returns 0 and empty list
234 //**************************** part 1 *******************************************
235 // Performs the unfolding of a cluster with nMax overlapping showers
237 //cout<<"unfolding check here part 1"<<endl;
238 AliDebug(5,Form(" Original cluster E %f, nMax = %d",iniTower->GetEnergy(),nMax ));
240 Int_t nPar = 3 * nMax ;
241 Float_t * fitparameters = new Float_t[nPar] ;
244 AliFatal("Did not get geometry from EMCALLoader") ;
246 Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
249 // Fit failed, return (and remove cluster? - why? I leave the cluster)
250 iniTower->SetNExMax(-1) ;
251 delete[] fitparameters ;
252 return 0;//changed here
255 //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold
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
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
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
277 AliEMCALDigit * digit = 0 ;
278 Int_t * digitsList = iniTower->GetDigitsList() ;
284 Int_t iphi = 0 ;//x direction
285 Int_t ieta = 0 ;//z direstion
290 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
292 digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
295 fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
296 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
297 iIphi, iIeta,iphi,ieta);
298 EvalParsPhiDependence(digit->GetId(),fGeom);
302 while(iparam < nPar )
304 xpar = fitparameters[iparam] ;
305 zpar = fitparameters[iparam+1] ;
306 epar = fitparameters[iparam+2] ;
308 efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
312 } else AliDebug(1,"Digit NULL part 2!");
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
321 Float_t * energiesList = iniTower->GetEnergiesList() ;
323 Float_t eDigit = 0. ;
324 Int_t nSplittedClusters=(Int_t)nPar/3;
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.
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.
336 // cout<<"unfolding check here part 3.1"<<endl;
339 while(iparam < nPar )
341 xpar = fitparameters[iparam] ;
342 zpar = fitparameters[iparam+1] ;
343 epar = fitparameters[iparam+2] ;
345 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
347 digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
350 fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
351 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
352 iIphi, iIeta,iphi,ieta);
354 EvalParsPhiDependence(digit->GetId(),fGeom);
358 correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here
362 ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
363 eDigit = energiesList[iDigit] * ratio ;
365 //add energy to temporary matrix
366 correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;
368 } else AliDebug(1,"NULL digit part 3");
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;
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
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.;
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];
400 for(iDigit = 0 ; iDigit < nDigits ; iDigit++){
401 isAnyBelowThreshold=kFALSE;
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)
407 maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];
408 maximumIndex = iparam;
410 }//end of loop over clusters after unfolding
412 if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold
414 if(maximumEne < fThreshold)
415 {//add all cluster cells and put energy into max index, other set to 0
417 for(iparam = 0 ; iparam < nPar ; iparam+=3)
419 maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
420 correctedEnergyList[iparam/3*nDigits+iDigit]=0;
422 correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;
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;
433 energyFraction[iparam/3]=1.;
434 maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
436 }//end of loop over clusters after unfolding
438 for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction
439 energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;
442 for(iparam = 0 ; iparam < nPar ; iparam+=3)
443 {//add energy from cells below threshold to others
444 if(energyFraction[iparam/3]>0.) continue;
447 for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)
449 correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] *
450 correctedEnergyList[iparam/3*nDigits+iDigit]) ;
452 correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
456 //digit energy to be set to 0
457 for(iparam = 0 ; iparam < nPar ; iparam+=3)
459 correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
461 }//correction for: is energy>0
463 }//end of loop over digits
464 delete[] energyFraction;
466 }//end of option 2 or 3
472 //**************************** sub-part 3.3 *************************************
473 //here we add digits to recpoints with corrected energy
474 // cout<<"unfolding check here part 3.3"<<endl;
476 Int_t newClusterIndex=0;
478 while(iparam < nPar )
480 AliEMCALRecPoint * recPoint = 0 ;
482 if(nSplittedClusters >= list->GetSize())
483 list->Expand(nSplittedClusters);
486 (*list)[newClusterIndex] = new AliEMCALRecPoint("") ;
487 recPoint = dynamic_cast<AliEMCALRecPoint *>( list->At(newClusterIndex) ) ;
489 if(recPoint){//recPoint present -> good
490 recPoint->SetNExMax(nSplittedClusters) ;//can be wrong number, to be corrected in outer method
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
498 AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));
502 if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list
503 delete (*list)[newClusterIndex];
504 list->RemoveAt(newClusterIndex);
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]));
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);
517 newClusterIndex--;//decrease cluster number
524 delete[] fitparameters ;
526 delete[] correctedEnergyList ;
529 AliDebug(5,Form(" nSplittedClusters %d, fNumberOfECAClusters %d, newClusterIndex %d,list->Entries() %d\n",nSplittedClusters,fNumberOfECAClusters,newClusterIndex,list->GetEntriesFast() ));
531 // cout<<"end of unfolding check part 3.3"<<endl;
532 return nSplittedClusters;
535 //____________________________________________________________________________
536 Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
538 AliEMCALDigit ** maxAt,
539 Float_t * maxAtEnergy)
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)
545 TObjArray *list =new TObjArray(2);//temporary object
546 Int_t nUnfoldedClusters=UnfoldOneCluster(iniTower,nMax,maxAt,maxAtEnergy,list);
548 // here we write new clusters from list to fRecPoints
549 AliDebug(5,Form("Number of clusters after unfolding %d",list->GetEntriesFast()));
551 AliEMCALDigit * digit = 0 ;
552 for(Int_t i=0;i<list->GetEntriesFast();i++) {
553 AliEMCALRecPoint * recPoint = 0 ;
555 if(fNumberOfECAClusters >= fRecPoints->GetSize())
556 fRecPoints->Expand(2*fNumberOfECAClusters) ;
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)) ;
563 if( recPoint && rpUFOne ){//recPoint present -> good
565 recPoint->SetNExMax(list->GetEntriesFast()) ;
567 Int_t *digitsList = rpUFOne->GetDigitsList();
568 Float_t *energyList = rpUFOne->GetEnergiesList();
570 if(!digitsList || ! energyList)
572 AliDebug(-1,"No digits index or energy available");
573 delete (*fRecPoints)[fNumberOfECAClusters];
574 fRecPoints->RemoveAt(fNumberOfECAClusters);
578 AliDebug(5,Form("cluster %d, digit no %d, energy %f\n",i,digitsList[0],energyList[0]));
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
584 fNumberOfECAClusters++ ;
585 } else {//recPoint empty -> remove from list
586 AliError("NULL RecPoint");
587 delete (*fRecPoints)[fNumberOfECAClusters];
588 fRecPoints->RemoveAt(fNumberOfECAClusters);
591 }//loop over unfolded clusters
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() ));
601 list->SetOwner(kTRUE);
604 if(nUnfoldedClusters>1) return kTRUE;
610 //____________________________________________________________________________
611 Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower,
613 AliEMCALDigit ** maxAt,
614 Float_t * maxAtEnergy)
616 // Extended to whole EMCAL
617 // Performs the unfolding of a cluster with nMax overlapping showers
619 Int_t nPar = 3 * nMax ;
620 Float_t * fitparameters = new Float_t[nPar] ;
623 AliFatal("Did not get geometry from EMCALLoader") ;
625 Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
627 // Fit failed, return (and remove cluster? - why? I leave the cluster)
628 iniTower->SetNExMax(-1) ;
629 delete[] fitparameters ;
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
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
643 AliEMCALDigit * digit = 0 ;
644 Int_t * digitsList = iniTower->GetDigitsList() ;
650 Int_t iphi = 0 ;//x direction
651 Int_t ieta = 0 ;//z direstion
656 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
657 digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
659 fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
660 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
661 iIphi, iIeta,iphi,ieta);
662 EvalParsPhiDependence(digit->GetId(),fGeom);
666 while(iparam < nPar ){
667 xpar = fitparameters[iparam] ;
668 zpar = fitparameters[iparam+1] ;
669 epar = fitparameters[iparam+2] ;
672 efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
674 } else AliError("Digit NULL!");
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
682 Float_t * energiesList = iniTower->GetEnergiesList() ;
686 while(iparam < nPar ){
687 xpar = fitparameters[iparam] ;
688 zpar = fitparameters[iparam+1] ;
689 epar = fitparameters[iparam+2] ;
692 AliEMCALRecPoint * recPoint = 0 ;
694 if(fNumberOfECAClusters >= fRecPoints->GetSize())
695 fRecPoints->Expand(2*fNumberOfECAClusters) ;
698 (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
699 recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
703 fNumberOfECAClusters++ ;
704 recPoint->SetNExMax((Int_t)nPar/3) ;
706 Float_t eDigit = 0. ;
707 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
708 digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
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");
720 } else AliError("NULL RecPoint");
723 delete[] fitparameters ;
730 //____________________________________________________________________________
731 Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt,
732 const Float_t* maxAtEnergy,
733 Int_t nPar, Float_t * fitparameters) const
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
741 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
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
749 if(gMinuit->fMaxpar < nPar) {
751 gMinuit = new TMinuit(nPar);
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) ;
764 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
766 // filling initial values for fit parameters
767 AliEMCALDigit * digit ;
771 Int_t nDigits = (Int_t) nPar / 3 ;
779 Int_t iphi = 0 ;//x direction
780 Int_t ieta = 0 ;//z direstion
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);
789 Float_t energy = maxAtEnergy[iDigit] ;
791 //gMinuit->mnparm(index, "x", iphi, 0.1, 0, 0, ierflg) ;//original
792 gMinuit->mnparm(index, "x", iphi, 0.05, 0, 0, ierflg) ;
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 ) ;
800 //gMinuit->mnparm(index, "z", ieta, 0.1, 0, 0, ierflg) ;//original
801 gMinuit->mnparm(index, "z", ieta, 0.05, 0, 0, ierflg) ;
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) ;
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
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) ;
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
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
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 " ) ;
839 for(index = 0; index < nPar; index++){
842 gMinuit->GetParameter(index, val, err) ; // Returns value and error ofOA parameter index
843 fitparameters[index] = val ;
849 if(gMinuit->fMaxpar>30) delete gMinuit;
855 //____________________________________________________________________________
856 Double_t AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)
858 // extended to whole EMCAL
859 // Shape of the shower
860 // If you change this function, change also the gradient evaluation in ChiSquare()
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) ) ) ;
869 //____________________________________________________________________________
870 void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
872 Double_t * x, Int_t iflag)
874 // Calculates the Chi square for the cluster unfolding minimization
875 // Number of parameters, Gradient, Chi squared, parameters, what to do
877 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
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
883 AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
885 if(recPoint && digits && geom){
887 Int_t * digitsList = recPoint->GetDigitsList() ;
889 Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
891 Float_t * energiesList = recPoint->GetEnergiesList() ;
897 for(iparam = 0 ; iparam < nPar ; iparam++)
898 Grad[iparam] = 0 ; // Will evaluate gradient
902 AliEMCALDigit * digit ;
909 Int_t iphi = 0 ;//x direction
910 Int_t ieta = 0 ;//z direstion
913 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
914 if(energiesList[iDigit]==0) continue;
916 digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
919 geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
920 geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
921 iIphi, iIeta,iphi,ieta);
922 EvalParsPhiDependence(digit->GetId(),geom);
924 if(iflag == 2){ // calculate gradient
927 while(iParam < nPar ){
928 Double_t dx = ((Float_t)iphi - x[iParam]) ;
930 Double_t dz = ((Float_t)ieta - x[iParam]) ;
932 efit += x[iParam] * ShowerShapeV2(dx,dz) ;
936 Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
938 while(iParam < nPar ){
939 Double_t xpar = x[iParam] ;
940 Double_t zpar = x[iParam+1] ;
941 Double_t epar = x[iParam+2] ;
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]) ;
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) ) ) );
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) ) ) ;
956 Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x
958 Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z
960 Grad[iParam] += shape ; // Derivative over energy
967 while(iparam < nPar ){
968 Double_t xpar = x[iparam] ;
969 Double_t zpar = x[iparam+1] ;
970 Double_t epar = x[iparam+2] ;
972 efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
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
979 } // recpoint, digits and geom not NULL
985 //____________________________________________________________________________
986 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){
987 for(UInt_t i=0;i<7;++i)
989 if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0
992 //____________________________________________________________________________
993 void AliEMCALUnfolding::SetPar5(Double_t *pars){
994 for(UInt_t i=0;i<3;++i)
998 //____________________________________________________________________________
999 void AliEMCALUnfolding::SetPar6(Double_t *pars){
1000 for(UInt_t i=0;i<3;++i)
1004 //____________________________________________________________________________
1005 void AliEMCALUnfolding::EvalPar5(Double_t phi){
1007 //Evaluate the 5th parameter of the shower shape function
1008 //phi in degrees range (-10,10)
1010 //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;
1011 fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];
1014 //____________________________________________________________________________
1015 void AliEMCALUnfolding::EvalPar6(Double_t phi){
1017 //Evaluate the 6th parameter of the shower shape function
1018 //phi in degrees range (-10,10)
1020 //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;
1021 fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];
1024 //____________________________________________________________________________
1025 void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){
1027 // calculate params p5 and p6 depending on the phi angle in global coordinate
1028 // for the cell with given absId index
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();
1035 phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);