{\r
// Unfolds clusters using the shape of an ElectroMagnetic shower\r
// Performs unfolding of all clusters\r
-\r
+ \r
if(fNumberOfECAClusters > 0){\r
if (fGeom==0)\r
AliFatal("Did not get geometry from EMCALLoader") ;\r
//Int_t nModulesToUnfold = fGeom->GetNCells();\r
-\r
+ \r
Int_t numberofNotUnfolded = fNumberOfECAClusters ;\r
Int_t index ;\r
for(index = 0 ; index < numberofNotUnfolded ; index++){\r
AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r
-\r
- //do we really need it?\r
-// TVector3 gpos;\r
-// Int_t absId = -1;\r
-// recPoint->GetGlobalPosition(gpos);\r
-// fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r
-// if(absId > nModulesToUnfold)\r
-// break ;\r
-\r
- Int_t nMultipl = recPoint->GetMultiplicity() ;\r
- AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
- Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
- Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
-\r
- if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0\r
- if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
- fRecPoints->Remove(recPoint);\r
- fRecPoints->Compress() ;//is it really needed\r
- index-- ;\r
- fNumberOfECAClusters-- ;\r
- numberofNotUnfolded-- ;\r
- }\r
- }\r
- else{\r
- recPoint->SetNExMax(1) ; //Only one local maximum\r
- }\r
-\r
- delete[] maxAt ;\r
- delete[] maxAtEnergy ;\r
- }\r
+ if(recPoint){\r
+ //do we really need it?\r
+ // TVector3 gpos;\r
+ // Int_t absId = -1;\r
+ // recPoint->GetGlobalPosition(gpos);\r
+ // fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r
+ // if(absId > nModulesToUnfold)\r
+ // break ;\r
+ \r
+ Int_t nMultipl = recPoint->GetMultiplicity() ;\r
+ AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
+ Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
+ Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
+ \r
+ if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0\r
+ if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
+ fRecPoints->Remove(recPoint);\r
+ fRecPoints->Compress() ;//is it really needed\r
+ index-- ;\r
+ fNumberOfECAClusters-- ;\r
+ numberofNotUnfolded-- ;\r
+ }\r
+ }\r
+ else{\r
+ recPoint->SetNExMax(1) ; //Only one local maximum\r
+ }\r
+ \r
+ delete[] maxAt ;\r
+ delete[] maxAtEnergy ;\r
+ } else AliError("RecPoint NULL");\r
+ } // rec point loop\r
}\r
// End of Unfolding of clusters\r
}\r
\r
Int_t nPar = 3 * nMax ;\r
Float_t * fitparameters = new Float_t[nPar] ;\r
-\r
+ \r
if (fGeom==0)\r
AliFatal("Did not get geometry from EMCALLoader") ;\r
-\r
+ \r
Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
if( !rv ) {\r
// Fit failed, return (and remove cluster? - why? I leave the cluster)\r
delete[] fitparameters ;\r
return kFALSE;\r
}\r
-\r
+ \r
// create unfolded rec points and fill them with new energy lists\r
// First calculate energy deposited in each sell in accordance with\r
// fit (without fluctuations): efit[]\r
// and later correct this number in acordance with actual energy\r
// deposition\r
-\r
+ \r
Int_t nDigits = iniTower->GetMultiplicity() ;\r
Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units\r
-\r
+ \r
AliEMCALDigit * digit = 0 ;\r
Int_t * digitsList = iniTower->GetDigitsList() ;\r
\r
Int_t iIeta = 0 ;\r
Int_t iphi = 0 ;//x direction\r
Int_t ieta = 0 ;//z direstion\r
-\r
+ \r
Int_t iparam = 0 ;\r
Int_t iDigit = 0 ;\r
-\r
+ \r
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
- fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
- fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
- iIphi, iIeta,iphi,ieta);\r
- EvalParsPhiDependence(digit->GetId(),fGeom);\r
-\r
- efit[iDigit] = 0.;\r
- iparam = 0;\r
- while(iparam < nPar ){\r
- xpar = fitparameters[iparam] ;\r
- zpar = fitparameters[iparam+1] ;\r
- epar = fitparameters[iparam+2] ;\r
- iparam += 3 ;\r
-\r
- efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
- }\r
-\r
- }\r
-\r
+ if(digit){\r
+ fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+ fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+ iIphi, iIeta,iphi,ieta);\r
+ EvalParsPhiDependence(digit->GetId(),fGeom);\r
+ \r
+ efit[iDigit] = 0.;\r
+ iparam = 0;\r
+ while(iparam < nPar ){\r
+ xpar = fitparameters[iparam] ;\r
+ zpar = fitparameters[iparam+1] ;\r
+ epar = fitparameters[iparam+2] ;\r
+ iparam += 3 ;\r
+ \r
+ efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+ }\r
+ } else AliError("Digit NULL!");\r
+ \r
+ }//digit loop\r
+ \r
// Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
// so that energy deposited in each cell is distributed between new clusters proportionally\r
// to its contribution to efit\r
-\r
+ \r
Float_t * energiesList = iniTower->GetEnergiesList() ;\r
Float_t ratio = 0 ;\r
-\r
+ \r
iparam = 0 ;\r
while(iparam < nPar ){\r
xpar = fitparameters[iparam] ;\r
zpar = fitparameters[iparam+1] ;\r
epar = fitparameters[iparam+2] ;\r
iparam += 3 ;\r
-\r
+ \r
AliEMCALRecPoint * recPoint = 0 ;\r
-\r
+ \r
if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
-\r
+ \r
//add recpoint\r
(*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
- fNumberOfECAClusters++ ;\r
- recPoint->SetNExMax((Int_t)nPar/3) ;\r
-\r
- Float_t eDigit = 0. ;\r
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
- digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
- fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
- fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
- iIphi, iIeta,iphi,ieta);\r
- EvalParsPhiDependence(digit->GetId(),fGeom);\r
- if(efit[iDigit]==0) continue;//just for sure\r
- ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
- eDigit = energiesList[iDigit] * ratio ;\r
- recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
- }\r
-\r
- }\r
-\r
+ \r
+ if(recPoint){\r
+ \r
+ fNumberOfECAClusters++ ;\r
+ recPoint->SetNExMax((Int_t)nPar/3) ;\r
+ \r
+ Float_t eDigit = 0. ;\r
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
+ digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
+ if(digit){\r
+ fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+ fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+ iIphi, iIeta,iphi,ieta);\r
+ EvalParsPhiDependence(digit->GetId(),fGeom);\r
+ if(efit[iDigit]==0) continue;//just for sure\r
+ ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
+ eDigit = energiesList[iDigit] * ratio ;\r
+ recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
+ } else AliError("NULL digit");\r
+ }//digit loop \r
+ } else AliError("NULL RecPoint");\r
+ }//while\r
+ \r
delete[] fitparameters ;\r
delete[] efit ;\r
-\r
+ \r
return kTRUE;\r
}\r
\r
// Cluster will be fitted as a superposition of nPar/3\r
// electromagnetic showers\r
\r
- cout<<"inside FindFitV2"<<endl;\r
-\r
if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r
\r
if(!gMinuit)\r
{\r
// Calculates the Chi square for the cluster unfolding minimization\r
// Number of parameters, Gradient, Chi squared, parameters, what to do\r
-\r
+ \r
TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
-\r
- AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;\r
- TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;\r
- // A bit buggy way to get an access to the geometry\r
- // To be revised!\r
- AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
-\r
- Int_t * digitsList = recPoint->GetDigitsList() ;\r
-\r
- Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
-\r
- Float_t * energiesList = recPoint->GetEnergiesList() ;\r
-\r
- fret = 0. ;\r
- Int_t iparam = 0 ;\r
-\r
- if(iflag == 2)\r
- for(iparam = 0 ; iparam < nPar ; iparam++)\r
- Grad[iparam] = 0 ; // Will evaluate gradient\r
-\r
- Double_t efit = 0. ;\r
-\r
- AliEMCALDigit * digit ;\r
- Int_t iDigit ;\r
-\r
- Int_t iSupMod = 0 ;\r
- Int_t iTower = 0 ;\r
- Int_t iIphi = 0 ;\r
- Int_t iIeta = 0 ;\r
- Int_t iphi = 0 ;//x direction\r
- Int_t ieta = 0 ;//z direstion\r
-\r
-\r
- for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
- if(energiesList[iDigit]==0) continue;\r
-\r
- digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
-\r
- geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
- geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
- iIphi, iIeta,iphi,ieta);\r
- EvalParsPhiDependence(digit->GetId(),geom);\r
-\r
- if(iflag == 2){ // calculate gradient\r
- Int_t iParam = 0 ;\r
- efit = 0. ;\r
- while(iParam < nPar ){\r
- Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
- iParam++ ;\r
- Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
- iParam++ ;\r
- efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
- iParam++ ;\r
- }\r
-\r
- Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
- iParam = 0 ;\r
- while(iParam < nPar ){\r
- Double_t xpar = x[iParam] ;\r
- Double_t zpar = x[iParam+1] ;\r
- Double_t epar = x[iParam+2] ;\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
- Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
- Double_t rp1 = TMath::Power(dr, fSSPars[1]) ;\r
- Double_t rp5 = TMath::Power(dr, fSSPars[5]) ;\r
-\r
- Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
- (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
- (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
- fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\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
- // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
-\r
- Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x\r
- iParam++ ;\r
- Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z\r
- iParam++ ;\r
- Grad[iParam] += shape ; // Derivative over energy\r
- iParam++ ;\r
- }\r
- }\r
- efit = 0;\r
- iparam = 0 ;\r
-\r
- while(iparam < nPar ){\r
- Double_t xpar = x[iparam] ;\r
- Double_t zpar = x[iparam+1] ;\r
- Double_t epar = x[iparam+2] ;\r
- iparam += 3 ;\r
- efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
- }\r
-\r
- fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
- // Here we assume, that sigma = sqrt(E) \r
- }\r
-\r
+ if(toMinuit){\r
+ AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;\r
+ TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;\r
+ // A bit buggy way to get an access to the geometry\r
+ // To be revised!\r
+ AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
+ \r
+ if(recPoint && digits && geom){\r
+ \r
+ Int_t * digitsList = recPoint->GetDigitsList() ;\r
+ \r
+ Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
+ \r
+ Float_t * energiesList = recPoint->GetEnergiesList() ;\r
+ \r
+ fret = 0. ;\r
+ Int_t iparam = 0 ;\r
+ \r
+ if(iflag == 2)\r
+ for(iparam = 0 ; iparam < nPar ; iparam++)\r
+ Grad[iparam] = 0 ; // Will evaluate gradient\r
+ \r
+ Double_t efit = 0. ;\r
+ \r
+ AliEMCALDigit * digit ;\r
+ Int_t iDigit ;\r
+ \r
+ Int_t iSupMod = 0 ;\r
+ Int_t iTower = 0 ;\r
+ Int_t iIphi = 0 ;\r
+ Int_t iIeta = 0 ;\r
+ Int_t iphi = 0 ;//x direction\r
+ Int_t ieta = 0 ;//z direstion\r
+ \r
+ \r
+ for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
+ if(energiesList[iDigit]==0) continue;\r
+ \r
+ digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
+ \r
+ if(digit){\r
+ geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+ geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+ iIphi, iIeta,iphi,ieta);\r
+ EvalParsPhiDependence(digit->GetId(),geom);\r
+ \r
+ if(iflag == 2){ // calculate gradient\r
+ Int_t iParam = 0 ;\r
+ efit = 0. ;\r
+ while(iParam < nPar ){\r
+ Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
+ iParam++ ;\r
+ Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
+ iParam++ ;\r
+ efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
+ iParam++ ;\r
+ }\r
+ \r
+ Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
+ iParam = 0 ;\r
+ while(iParam < nPar ){\r
+ Double_t xpar = x[iParam] ;\r
+ Double_t zpar = x[iParam+1] ;\r
+ Double_t epar = x[iParam+2] ;\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
+ Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+ Double_t rp1 = TMath::Power(dr, fSSPars[1]) ;\r
+ Double_t rp5 = TMath::Power(dr, fSSPars[5]) ;\r
+ \r
+ Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
+ (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
+ (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
+ fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\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
+ // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
+ \r
+ Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x\r
+ iParam++ ;\r
+ Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z\r
+ iParam++ ;\r
+ Grad[iParam] += shape ; // Derivative over energy\r
+ iParam++ ;\r
+ }\r
+ }\r
+ efit = 0;\r
+ iparam = 0 ;\r
+ \r
+ while(iparam < nPar ){\r
+ Double_t xpar = x[iparam] ;\r
+ Double_t zpar = x[iparam+1] ;\r
+ Double_t epar = x[iparam+2] ;\r
+ iparam += 3 ;\r
+ efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+ }\r
+ \r
+ fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
+ // Here we assume, that sigma = sqrt(E) \r
+ } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n");\r
+ } // digit loop\r
+ } // recpoint, digits and geom not NULL\r
+ }// List is not NULL\r
+ \r
}\r
\r
\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff