\r\r
// --- ROOT system ---\r\r
#include "TClonesArray.h"\r\r
-//#include "TTree.h"\r\r
-//#include <TFile.h> \r\r
-//class TFolder;\r\r
#include <TMath.h> \r\r
#include <TMinuit.h>\r\r
-//#include <TTree.h> \r\r
-//class TSystem; \r\r
-//#include <TBenchmark.h>\r\r
-//#include <TBrowser.h>\r\r
-//#include <TROOT.h>\r\r
\r\r
// --- Standard library ---\r\r
#include <cassert>\r\r
#include "AliEMCALRecPoint.h"\r\r
#include "AliEMCALDigit.h"\r\r
#include "AliEMCALReconstructor.h"\r\r
-//#include "AliEMCALClusterizer.h"\r\r
-\r\r
-\r\r
\r\r
#include "AliLog.h"\r\r
-\r\r
#include "AliCDBManager.h"\r\r
-//#include "AliCaloCalibPedestal.h"\r\r
-//#include "AliEMCALCalibData.h"\r\r
class AliCDBStorage;\r\r
#include "AliCDBEntry.h"\r\r
\r\r
-Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r\r
-Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r\r
-Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};\r\r
+Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r\r
+Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};\r\r
+Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};\r\r
\r\r
ClassImp(AliEMCALUnfolding)\r\r
\r\r
fDigitsArr(NULL)\r\r
{\r\r
// ctor with the indication of the file where header Tree and digits Tree are stored\r\r
- \r\r
Init() ;\r\r
}\r\r
\r\r
{\r\r
AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r\r
}\r\r
- \r\r
+\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r\r
}\r\r
Int_t i=0;\r\r
- for (i = 0; i < 8; i++) fSSPars[i] = SSPars[i];\r\r
+ for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];\r\r
for (i = 0; i < 3; i++) {\r\r
- fPar5[i] = Par5[i];\r\r
- fPar6[i] = Par6[i];\r\r
+ fgPar5[i] = Par5[i];\r\r
+ fgPar6[i] = Par6[i];\r\r
}\r\r
\r\r
}\r\r
AliDebug(1,Form("geom %p",fGeom));\r\r
\r\r
if(!gMinuit) \r\r
- gMinuit = new TMinuit(100) ;\r\r
+ // gMinuit = new TMinuit(100) ;//the same is in FindFitV2\r\r
+ gMinuit = new TMinuit(30) ;//the same is in FindFitV2\r\r
\r\r
}\r\r
\r\r
{\r\r
//\r\r
//Set input for unfolding purposes\r\r
+ //\r\r
SetNumberOfECAClusters(numberOfECAClusters);\r\r
SetRecPoints(recPoints);\r\r
SetDigitsArr(digitsArr);\r\r
AliFatal("Did not get geometry from EMCALLoader") ;\r\r
//Int_t nModulesToUnfold = fGeom->GetNCells();\r\r
\r\r
- Int_t numberofNotUnfolded = fNumberOfECAClusters ;\r\r
+ Int_t numberOfClustersToUnfold=fNumberOfECAClusters;\r\r
+ //we unfold only clusters present in the array untill now\r\r
+ //fNumberOfECAClusters may change due to unfilded clusters\r\r
+ //so 0 to numberOfClustersToUnfold-1: clusters before unfolding\r\r
+ //numberOfClustersToUnfold to the end: new clusters from unfolding\r\r
+ //of course numberOfClustersToUnfold also is decreased but we don't loop over clusters added in UF \r\r
Int_t index ;\r\r
- for(index = 0 ; index < numberofNotUnfolded ; index++){\r\r
+ for(index = 0 ; index < numberOfClustersToUnfold ; index++){\r\r
AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r\r
if(recPoint){\r\r
- //do we really need it?\r\r
- // TVector3 gpos;\r\r
- // Int_t absId = -1;\r\r
- // recPoint->GetGlobalPosition(gpos);\r\r
- // fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r\r
- // if(absId > nModulesToUnfold)\r\r
- // break ;\r\r
- \r\r
Int_t nMultipl = recPoint->GetMultiplicity() ;\r\r
AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r\r
Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r\r
Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r\r
- \r\r
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0\r\r
if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r\r
+ //if unfolding correct remove old recPoint\r\r
fRecPoints->Remove(recPoint);\r\r
fRecPoints->Compress() ;//is it really needed\r\r
index-- ;\r\r
fNumberOfECAClusters-- ;\r\r
- numberofNotUnfolded-- ;\r\r
+ numberOfClustersToUnfold--;\r\r
}\r\r
- }\r\r
- else{\r\r
+ } else{\r\r
recPoint->SetNExMax(1) ; //Only one local maximum\r\r
}\r\r
\r\r
delete[] maxAt ;\r\r
delete[] maxAtEnergy ;\r\r
- } else AliError("RecPoint NULL");\r\r
+ } else {\r\r
+ //AliError("RecPoint NULL"); //end of check if recPoint exist\r\r
+ Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ;\r\r
+ }\r\r
} // rec point loop\r\r
- }\r\r
+ }//end of check fNumberOfECAClusters\r\r
// End of Unfolding of clusters\r\r
}\r\r
\r\r
Float_t * maxAtEnergy)\r\r
{\r\r
// Extended to whole EMCAL \r\r
-\r\r
+ \r\r
//**************************** part 1 *******************************************\r\r
// Performs the unfolding of a cluster with nMax overlapping showers \r\r
\r\r
+ //cout<<"unfolding check here part 1"<<endl;\r\r
+ //printf("Original cluster E %f\n",iniTower->GetEnergy());\r\r
+ \r\r
Int_t nPar = 3 * nMax ;\r\r
Float_t * fitparameters = new Float_t[nPar] ;\r\r
\r\r
AliFatal("Did not get geometry from EMCALLoader") ;\r\r
\r\r
Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r\r
- if( !rv ) {\r\r
+ if( !rv ) \r\r
+ {\r\r
// Fit failed, return (and remove cluster? - why? I leave the cluster)\r\r
iniTower->SetNExMax(-1) ;\r\r
delete[] fitparameters ;\r\r
return kFALSE;\r\r
}\r\r
\r\r
+ //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold\r\r
+ if(nMax==2){\r\r
+ if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){\r\r
+ AliDebug(1,"One of fitted energy below threshold");\r\r
+ iniTower->SetNExMax(1) ;\r\r
+ delete[] fitparameters ;\r\r
+ return kFALSE;\r\r
+ }\r\r
+ }\r\r
+\r\r
//**************************** part 2 *******************************************\r\r
// create unfolded rec points and fill them with new energy lists\r\r
// First calculate energy deposited in each sell in accordance with\r\r
// and later correct this number in acordance with actual energy\r\r
// deposition\r\r
\r\r
+ // cout<<"unfolding check here part 2"<<endl;\r\r
Int_t nDigits = iniTower->GetMultiplicity() ;\r\r
Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r\r
Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units\r\r
Int_t iparam = 0 ;\r\r
Int_t iDigit = 0 ;\r\r
\r\r
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r\r
+ {\r\r
digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r\r
- if(digit){\r\r
+ if(digit)\r\r
+ {\r\r
fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
iIphi, iIeta,iphi,ieta);\r\r
\r\r
efit[iDigit] = 0.;\r\r
iparam = 0;\r\r
- while(iparam < nPar ){\r\r
+ while(iparam < nPar )\r\r
+ {\r\r
xpar = fitparameters[iparam] ;\r\r
zpar = fitparameters[iparam+1] ;\r\r
epar = fitparameters[iparam+2] ;\r\r
- iparam += 3 ;\r\r
- \r\r
+\r\r
efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
+ iparam += 3 ;\r\r
}\r\r
- } else AliError("Digit NULL!");\r\r
+\r\r
+ } else AliDebug(1,"Digit NULL part 2!");\r\r
\r\r
}//digit loop\r\r
\r\r
// to its contribution to efit\r\r
\r\r
Float_t * energiesList = iniTower->GetEnergiesList() ;\r\r
- Float_t ratio = 0 ;\r\r
+ Float_t ratio = 0. ;//0 -> 0. changed\r\r
Float_t eDigit = 0. ;\r\r
Int_t nSplittedClusters=(Int_t)nPar/3;\r\r
\r\r
Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters];\r\r
//above - temporary table with energies after unfolding.\r\r
- //the orderis following: \r\r
+ //the order is following: \r\r
//first cluster <first cell - last cell>, \r\r
//second cluster <first cell - last cell>, etc.\r\r
-\r\r
+ \r\r
//**************************** sub-part 3.1 *************************************\r\r
//here we check if energy of the cell in the cluster after unfolding is above threshold. \r\r
//If not the energy from a given cell in the cluster is divided in correct proportions \r\r
//in accordance to the other clusters and added to them and set to 0.\r\r
+ \r\r
+ // cout<<"unfolding check here part 3.1"<<endl;\r\r
\r\r
iparam = 0 ;\r\r
- while(iparam < nPar ){\r\r
+ while(iparam < nPar )\r\r
+ {\r\r
xpar = fitparameters[iparam] ;\r\r
zpar = fitparameters[iparam+1] ;\r\r
epar = fitparameters[iparam+2] ;\r\r
-\r\r
-\r\r
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
+ \r\r
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r\r
+ {\r\r
digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
- if(digit){\r\r
- fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
- fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
- iIphi, iIeta,iphi,ieta);\r\r
- EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
- if(efit[iDigit]==0) {//just for sure\r\r
- correctedEnergyList[iparam/3+iDigit] = 0;\r\r
- continue;\r\r
- }\r\r
- ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r\r
- eDigit = energiesList[iDigit] * ratio ;\r\r
-\r\r
- //add energy to temporary matrix\r\r
- correctedEnergyList[iparam/3+iDigit] = eDigit;\r\r
-\r\r
- } else AliError("NULL digit");\r\r
+ if(digit)\r\r
+ {\r\r
+ fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
+ fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
+ iIphi, iIeta,iphi,ieta);\r\r
+ \r\r
+ EvalParsPhiDependence(digit->GetId(),fGeom);\r\r
+ \r\r
+ if(efit[iDigit]==0) \r\r
+ {//just for sure\r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here\r\r
+ continue;\r\r
+ }\r\r
+ \r\r
+ ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r\r
+ eDigit = energiesList[iDigit] * ratio ;\r\r
+ \r\r
+ //add energy to temporary matrix\r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;\r\r
+ \r\r
+ } else AliDebug(1,"NULL digit part 3");\r\r
}//digit loop \r\r
iparam += 3 ;\r\r
}//while\r\r
-\r\r
+ \r\r
//**************************** sub-part 3.2 *************************************\r\r
//here we correct energy for each cell and cluster\r\r
- Float_t maximumEne=0;\r\r
+ // cout<<"unfolding check here part 3.2"<<endl;\r\r
+\r\r
+ Float_t maximumEne=0.;//0 -> 0. changed\r\r
Int_t maximumIndex=0;\r\r
Bool_t isAnyBelowThreshold=kFALSE;\r\r
// Float_t Threshold=0.01;\r\r
Float_t * energyFraction = new Float_t[nSplittedClusters];\r\r
Int_t iparam2 = 0 ;\r\r
- for(iDigit = 0 ; iDigit < nDigits ; ++iDigit){\r\r
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit++)\r\r
+ {\r\r
isAnyBelowThreshold=kFALSE;\r\r
- maximumEne=0;\r\r
- for(iparam = 0 ; iparam < nPar ; iparam+=3){\r\r
-\r\r
- if(correctedEnergyList[iparam/3+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;\r\r
- if(correctedEnergyList[iparam/3+iDigit] > maximumEne) {\r\r
- maximumEne = correctedEnergyList[iparam/3+iDigit];\r\r
- maximumIndex = iparam;\r\r
+ maximumEne=0.;//0 -> 0. changed\r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
+ {\r\r
+\r\r
+ if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;\r\r
+ if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne) \r\r
+ {\r\r
+ maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
+ maximumIndex = iparam;\r\r
}\r\r
}//end of loop over clusters after unfolding\r\r
-\r\r
+ \r\r
if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold \r\r
- if(maximumEne < fThreshold) {//add all cluster cells and put energy into max index, other set to 0\r\r
+\r\r
+\r\r
+ if(maximumEne < fThreshold) \r\r
+ {//add all cluster cells and put energy into max index, other set to 0\r\r
maximumEne=0.;\r\r
- for(iparam = 0 ; iparam < nPar ; iparam+=3){\r\r
- maximumEne+=correctedEnergyList[iparam/3+iDigit];\r\r
- correctedEnergyList[iparam/3+iDigit]=0;\r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
+ {\r\r
+ maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit]=0;\r\r
}\r\r
- correctedEnergyList[maximumIndex/3+iDigit]=maximumEne;\r\r
+ correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;\r\r
continue;\r\r
}//end if\r\r
-\r\r
+ \r\r
//divide energy of cell below threshold in the correct proportion and add to other cells\r\r
- maximumEne=0;//not used any more so use it for the energy sum\r\r
- for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate energy sum\r\r
- if(correctedEnergyList[iparam/3+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r\r
- else {\r\r
- energyFraction[iparam/3]=1;\r\r
- maximumEne+=correctedEnergyList[iparam/3+iDigit];\r\r
+ maximumEne=0.;//not used any more so use it for the energy sum \r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
+ {//calculate energy sum\r\r
+ if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r\r
+ else \r\r
+ {\r\r
+ energyFraction[iparam/3]=1.;\r\r
+ maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r\r
}\r\r
}//end of loop over clusters after unfolding\r\r
- for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r\r
- energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3+iDigit] / maximumEne;\r\r
- }\r\r
- for(iparam = 0 ; iparam < nPar ; iparam+=3){//add energy from cells below threshold to others\r\r
- if(energyFraction[iparam/3]>0) continue;\r\r
- else{\r\r
- for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3){\r\r
- correctedEnergyList[iparam2/3+iDigit] += (energyFraction[iparam2/3] * \r\r
- correctedEnergyList[iparam/3+iDigit]) ;\r\r
- }//inner loop\r\r
- correctedEnergyList[iparam/3+iDigit] = 0;\r\r
+ if(maximumEne>0.)\r\r
+ {\r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r\r
+ energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;\r\r
+ }\r\r
+ \r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
+ {//add energy from cells below threshold to others\r\r
+ if(energyFraction[iparam/3]>0.) continue;\r\r
+ else\r\r
+ {\r\r
+ for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)\r\r
+ {\r\r
+ correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * \r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit]) ;\r\r
+ }//inner loop\r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;\r\r
+ }\r\r
}\r\r
}\r\r
-\r\r
+ else\r\r
+ {\r\r
+ //digit energy to be set to 0\r\r
+ for(iparam = 0 ; iparam < nPar ; iparam+=3)\r\r
+ {\r\r
+ correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;\r\r
+ }\r\r
+ }//correction for: is energy>0\r\r
+ \r\r
}//end of loop over digits\r\r
delete[] energyFraction;\r\r
-\r\r
+ \r\r
//**************************** sub-part 3.3 *************************************\r\r
//here we add digits to recpoints with corrected energy\r\r
+ // cout<<"unfolding check here part 3.3"<<endl;\r\r
+\r\r
iparam = 0 ;\r\r
- while(iparam < nPar ){\r\r
+ while(iparam < nPar )\r\r
+ {\r\r
AliEMCALRecPoint * recPoint = 0 ;\r\r
\r\r
if(fNumberOfECAClusters >= fRecPoints->GetSize())\r\r
fRecPoints->Expand(2*fNumberOfECAClusters) ;\r\r
\r\r
//add recpoint\r\r
- (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r\r
+ (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;//fNumberOfECAClusters-1 is old cluster before unfolding\r\r
recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r\r
\r\r
- if(recPoint){\r\r
- \r\r
- fNumberOfECAClusters++ ;\r\r
+ if(recPoint){//recPoint present -> good\r\r
recPoint->SetNExMax(nSplittedClusters) ;\r\r
\r\r
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) {\r\r
digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
+ if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. ){\r\r
+ //if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold);\r\r
+ recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case\r\r
+ } else {\r\r
+ AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));\r\r
+ }\r\r
+ }//digit loop\r\r
+\r\r
+ if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list\r\r
+ delete (*fRecPoints)[fNumberOfECAClusters];\r\r
+ fRecPoints->RemoveAt(fNumberOfECAClusters);\r\r
+ //fNumberOfECAClusters--; //it should be like that we don't want to decrease cluster number\r\r
+ nSplittedClusters--;\r\r
+ } else {//recPoint exists and has digits associated -> very good increase number of clusters \r\r
+ fNumberOfECAClusters++ ; \r\r
+ }\r\r
+ \r\r
+ } else {//recPoint empty -> remove from list\r\r
+ AliError("NULL RecPoint");\r\r
+ //protection from recpoint with no digits\r\r
+// if(recPoint->GetMultiplicity()==0)\r\r
+// {\r\r
+ delete (*fRecPoints)[fNumberOfECAClusters];\r\r
+ fRecPoints->RemoveAt(fNumberOfECAClusters);\r\r
+ // fNumberOfECAClusters--;//it should be like that we don't want to decrease cluster number\r\r
+ nSplittedClusters--;\r\r
+// }\r\r
+ \r\r
+ }\r\r
\r\r
- if(digit && correctedEnergyList[iparam/3+iDigit]>0. ){\r\r
- recPoint->AddDigit( *digit, correctedEnergyList[iparam/3+iDigit], kFALSE ) ; //FIXME, need to study the shared case\r\r
- } else AliError("NULL digit");\r\r
- }//digit loop \r\r
- } else AliError("NULL RecPoint");\r\r
iparam += 3 ;\r\r
}//while\r\r
\r\r
delete[] fitparameters ;\r\r
delete[] efit ;\r\r
delete[] correctedEnergyList ;\r\r
-\r\r
+ \r\r
+ // cout<<"end of unfolding check part 3.3"<<endl;\r\r
return kTRUE;\r\r
}\r\r
\r\r
\r\r
if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r\r
\r\r
- if(!gMinuit)\r\r
- gMinuit = new TMinuit(100) ;//max 100 parameters\r\r
+ if(!gMinuit){\r\r
+ // gMinuit = new TMinuit(100) ;//max 100 parameters\r\r
+ if(nPar<30) gMinuit = new TMinuit(30);\r\r
+ else gMinuit = new TMinuit(nPar) ;//max nPar parameters\r\r
+ //\r\r
+ } else {\r\r
+ if(gMinuit->fMaxpar < nPar) {\r\r
+ delete gMinuit;\r\r
+ gMinuit = new TMinuit(nPar);\r\r
+ }\r\r
+ }\r\r
\r\r
gMinuit->mncler(); // Reset Minuit's list of paramters\r\r
gMinuit->SetPrintLevel(-1) ; // No Printout\r\r
\r\r
for(iDigit = 0; iDigit < nDigits; iDigit++){\r\r
digit = maxAt[iDigit];\r\r
+ if(digit==0) AliError("energy of digit = 0!");\r\r
fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r\r
fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r\r
iIphi, iIeta,iphi,ieta);\r\r
gMinuit->mnparm(index, "x", iphi, 0.05, 0, 0, ierflg) ;\r\r
index++ ;\r\r
if(ierflg != 0){\r\r
- Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %d", iphi ) ;\r\r
+ Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: x=%d, param.id=%d, nMaxima=%d",iphi,index-1,nPar/3 ) ;\r\r
toMinuit->Clear();\r\r
delete toMinuit ;\r\r
return kFALSE;\r\r
gMinuit->mnparm(index, "z", ieta, 0.05, 0, 0, ierflg) ;\r\r
index++ ;\r\r
if(ierflg != 0){\r\r
- Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %d", ieta) ;\r\r
+ Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: z=%d, param.id=%d, nMaxima=%d", ieta, index-1,nPar/3) ;\r\r
toMinuit->Clear();\r\r
delete toMinuit ;\r\r
return kFALSE;\r\r
gMinuit->mnparm(index, "Energy", energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r\r
index++ ;\r\r
if(ierflg != 0){\r\r
- Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;\r\r
+ Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: energy = %f, param.id=%d, nMaxima=%d", energy, index-1, nPar/3) ;\r\r
toMinuit->Clear();\r\r
delete toMinuit ;\r\r
return kFALSE;\r\r
gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize\r\r
//gMinuit->mnexcm("MINI", &p0, 0, ierflg) ; // minimize\r\r
if(ierflg == 4){ // Minimum not found\r\r
- Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;\r\r
+ AliDebug(1,"EMCAL Unfolding Fit not converged, cluster abandoned " ) ;\r\r
toMinuit->Clear();\r\r
delete toMinuit ;\r\r
return kFALSE ;\r\r
for(index = 0; index < nPar; index++){\r\r
Double_t err = 0. ;\r\r
Double_t val = 0. ;\r\r
- gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index\r\r
+ gMinuit->GetParameter(index, val, err) ; // Returns value and error ofOA parameter index\r\r
fitparameters[index] = val ;\r\r
}\r\r
\r\r
toMinuit->Clear();\r\r
delete toMinuit ;\r\r
+\r\r
+ if(gMinuit->fMaxpar>30) delete gMinuit;\r\r
+\r\r
return kTRUE;\r\r
\r\r
}\r\r
// Shape of the shower\r\r
// If you change this function, change also the gradient evaluation in ChiSquare()\r\r
\r\r
- Double_t r = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r\r
- Double_t rp1 = TMath::Power(r, fSSPars[1]) ;\r\r
- Double_t rp5 = TMath::Power(r, fSSPars[5]) ;\r\r
- Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[6] * rp5) ) ) ;\r\r
+ Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);\r\r
+ Double_t rp1 = TMath::Power(r, fgSSPars[1]) ;\r\r
+ Double_t rp5 = TMath::Power(r, fgSSPars[5]) ;\r\r
+ Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;\r\r
return shape ;\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r\r
- Double_t & fret,\r\r
- Double_t * x, Int_t iflag)\r\r
+ Double_t & fret,\r\r
+ Double_t * x, Int_t iflag)\r\r
{\r\r
// Calculates the Chi square for the cluster unfolding minimization\r\r
// Number of parameters, Gradient, Chi squared, parameters, what to do\r\r
Double_t zpar = x[iParam+1] ;\r\r
Double_t epar = x[iParam+2] ;\r\r
\r\r
- 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
+ Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r\r
Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r\r
- Double_t rp1 = TMath::Power(dr, fSSPars[1]) ;\r\r
- Double_t rp5 = TMath::Power(dr, fSSPars[5]) ;\r\r
+ Double_t rp1 = TMath::Power(dr, fgSSPars[1]) ;\r\r
+ Double_t rp5 = TMath::Power(dr, fgSSPars[5]) ;\r\r
\r\r
- Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r\r
- (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r\r
- (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r\r
- fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\r\r
+ Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * \r\r
+ (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - \r\r
+ (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + \r\r
+ fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );\r\r
\r\r
//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
// - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r\r
\r\r
fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r\r
// Here we assume, that sigma = sqrt(E) \r\r
- } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n");\r\r
+ } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker\r\r
} // digit loop\r\r
} // recpoint, digits and geom not NULL\r\r
}// List is not NULL\r\r
//____________________________________________________________________________\r\r
void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r\r
for(UInt_t i=0;i<7;++i)\r\r
- fSSPars[i]=pars[i];\r\r
- if(pars[2]==0. && pars[3]==0.) fSSPars[2]=1.;//to avoid dividing by 0\r\r
+ fgSSPars[i]=pars[i];\r\r
+ if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
void AliEMCALUnfolding::SetPar5(Double_t *pars){\r\r
for(UInt_t i=0;i<3;++i)\r\r
- fPar5[i]=pars[i];\r\r
+ fgPar5[i]=pars[i];\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
void AliEMCALUnfolding::SetPar6(Double_t *pars){\r\r
for(UInt_t i=0;i<3;++i)\r\r
- fPar6[i]=pars[i];\r\r
+ fgPar6[i]=pars[i];\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
//phi in degrees range (-10,10)\r\r
//\r\r
//fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r\r
- fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[2];\r\r
+ fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
//phi in degrees range (-10,10)\r\r
//\r\r
//fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r\r
- fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r\r
+ fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];\r\r
}\r\r
\r\r
//____________________________________________________________________________\r\r
-void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom){\r\r
+void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){\r\r
//\r\r
// calculate params p5 and p6 depending on the phi angle in global coordinate\r\r
// for the cell with given absId index\r\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff