From 622e10be4fc3011cab9894a487f5b08593658524 Mon Sep 17 00:00:00 2001 From: gconesab Date: Tue, 19 Oct 2010 08:30:17 +0000 Subject: [PATCH] Correct coverity FORWARD_NULL Remove prints in case of unfolding --- EMCAL/AliEMCALClusterizer.cxx | 48 ++-- EMCAL/AliEMCALRecParam.h | 2 +- EMCAL/AliEMCALSensorTempArray.cxx | 7 +- EMCAL/AliEMCALTrigger.cxx | 22 +- EMCAL/AliEMCALTriggerDCSConfigDB.cxx | 32 ++- EMCAL/AliEMCALUnfolding.cxx | 374 ++++++++++++++------------- 6 files changed, 256 insertions(+), 229 deletions(-) diff --git a/EMCAL/AliEMCALClusterizer.cxx b/EMCAL/AliEMCALClusterizer.cxx index c3a6d93d7e2..ad7a237941e 100644 --- a/EMCAL/AliEMCALClusterizer.cxx +++ b/EMCAL/AliEMCALClusterizer.cxx @@ -319,7 +319,6 @@ void AliEMCALClusterizer::InitParameters() fECAW0 = recParam->GetW0(); fMinECut = recParam->GetMinECut(); fToUnfold = recParam->GetUnfold(); - if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); fECALocMaxCut = recParam->GetLocMaxCut(); fTimeCut = recParam->GetTimeCut(); fTimeMin = recParam->GetTimeMin(); @@ -327,30 +326,31 @@ void AliEMCALClusterizer::InitParameters() AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s", fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax)); - } - - if(fToUnfold){ - Int_t i=0; - for (i = 0; i < 8; i++) { - fSSPars[i] = recParam->GetSSPars(i); - }//end of loop over parameters - for (i = 0; i < 3; i++) { - fPar5[i] = recParam->GetPar5(i); - fPar6[i] = recParam->GetPar6(i); - }//end of loop over parameters - - fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6); - for (i = 0; i < 8; i++) { - AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i])); - } - for (i = 0; i < 3; i++) { - AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i])); - AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i])); - } - - } - + if(fToUnfold){ + + Int_t i=0; + for (i = 0; i < 8; i++) { + fSSPars[i] = recParam->GetSSPars(i); + }//end of loop over parameters + for (i = 0; i < 3; i++) { + fPar5[i] = recParam->GetPar5(i); + fPar6[i] = recParam->GetPar6(i); + }//end of loop over parameters + + fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6); + + for (i = 0; i < 8; i++) { + AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i])); + } + for (i = 0; i < 3; i++) { + AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i])); + AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i])); + } + + }// to unfold + }// recparam not null + } diff --git a/EMCAL/AliEMCALRecParam.h b/EMCAL/AliEMCALRecParam.h index 3e2be8bfea0..b60fe412bd8 100644 --- a/EMCAL/AliEMCALRecParam.h +++ b/EMCAL/AliEMCALRecParam.h @@ -52,7 +52,7 @@ class AliEMCALRecParam : public AliDetectorRecoParam void SetTimeCut (Float_t t) {fTimeCut = t ;} void SetTimeMin (Float_t t) {fTimeMin = t ;} void SetTimeMax (Float_t t) {fTimeMax = t ;} - void SetUnfold (Bool_t unfold) {fUnfold = unfold ; if(fUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");} + void SetUnfold (Bool_t unfold) {fUnfold = unfold ;} //PID (Guenole) Double_t GetGamma(Int_t i, Int_t j) const {return fGamma[i][j];} diff --git a/EMCAL/AliEMCALSensorTempArray.cxx b/EMCAL/AliEMCALSensorTempArray.cxx index 7976bb3c725..ae47d5cd997 100644 --- a/EMCAL/AliEMCALSensorTempArray.cxx +++ b/EMCAL/AliEMCALSensorTempArray.cxx @@ -116,8 +116,11 @@ void AliEMCALSensorTempArray::ReadSensors(const char *dbEntry) // Read list of temperature sensors from text file // AliCDBEntry *entry = AliCDBManager::Instance()->Get(dbEntry); - TTree *tree = (TTree*) entry->GetObject(); - fSensors = AliEMCALSensorTemp::ReadTree(tree); + if(entry){ + TTree *tree = (TTree*) entry->GetObject(); + fSensors = AliEMCALSensorTemp::ReadTree(tree); + } + else AliFatal("NULL CDB entry!"); } //_____________________________________________________________________________ diff --git a/EMCAL/AliEMCALTrigger.cxx b/EMCAL/AliEMCALTrigger.cxx index 933d7fb0d3b..6e1f8bd8797 100644 --- a/EMCAL/AliEMCALTrigger.cxx +++ b/EMCAL/AliEMCALTrigger.cxx @@ -244,7 +244,8 @@ Bool_t AliEMCALTrigger::IsPatchIsolated(Int_t iPatchType, const TClonesArray * a // 1 4x4 8x8 Bool_t b = kFALSE; - + if(!ampmatrixes) return kFALSE; + // Get matrix of TRU or Module with maximum amplitude patch. Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36. TMatrixD * ampmatrix = 0x0; @@ -858,17 +859,18 @@ void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmat if(!amptrus || timeRtrus){ AliError("Could not recover the TRU matrix with amplitudes or times"); - continue; } - //Calculate row and column of the module inside the TRU with number itru - Int_t irow = iphim - row * nModulesPhi; - if(iSupMod > 9) - irow = iphim - row * nModulesPhi2; // size of matrix the same - Int_t icol = ietam - col * nModulesEta; + else{ + //Calculate row and column of the module inside the TRU with number itru + Int_t irow = iphim - row * nModulesPhi; + if(iSupMod > 9) + irow = iphim - row * nModulesPhi2; // size of matrix the same + Int_t icol = ietam - col * nModulesEta; - (*amptrus)(irow,icol) += amp ; - if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ?? - (*timeRtrus)(irow,icol) = timeR ; + (*amptrus)(irow,icol) += amp ; + if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ?? + (*timeRtrus)(irow,icol) = timeR ; + } } //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n", // ieta, iphi, iSupMod, col, row, itru, amp); diff --git a/EMCAL/AliEMCALTriggerDCSConfigDB.cxx b/EMCAL/AliEMCALTriggerDCSConfigDB.cxx index 9d243014205..666255d0f97 100644 --- a/EMCAL/AliEMCALTriggerDCSConfigDB.cxx +++ b/EMCAL/AliEMCALTriggerDCSConfigDB.cxx @@ -291,15 +291,19 @@ Int_t AliEMCALTriggerDCSConfigDB::GetTRUSegmentation(Int_t iTRU) { // const AliEMCALTriggerDCSConfig* dcsConf = dynamic_cast(GetCachedCDBObject(kIDTriggerConfig)); + if(dcsConf){ + AliEMCALTriggerTRUDCSConfig* truConf = dcsConf->GetTRUDCSConfig(iTRU); + if(truConf){ + Int_t sel = truConf->GetL0SEL(); - AliEMCALTriggerTRUDCSConfig* truConf = dcsConf->GetTRUDCSConfig(iTRU); - - Int_t sel = truConf->GetL0SEL(); - - if (sel & 0x0001) - return 2; - else - return 1; + if (sel & 0x0001) + return 2; + else + return 1; + } else AliFatal("TRUDCSConf Null!") ; + }else AliFatal("TriggerDCSConf Null!") ; + + return -1; } //_____________________________________________________________________________ @@ -309,8 +313,12 @@ Int_t AliEMCALTriggerDCSConfigDB::GetTRUGTHRL0(Int_t iTRU) // // const AliEMCALTriggerDCSConfig* dcsConf = dynamic_cast(GetCachedCDBObject(kIDTriggerConfig)); - - AliEMCALTriggerTRUDCSConfig* truConf = dcsConf->GetTRUDCSConfig(iTRU); - - return truConf->GetGTHRL0(); + if(dcsConf){ + AliEMCALTriggerTRUDCSConfig* truConf = dcsConf->GetTRUDCSConfig(iTRU); + if(truConf){ + return truConf->GetGTHRL0(); + } else AliFatal("TRUDCSConf Null!") ; + }else AliFatal("TriggerDCSConf Null!") ; + + return -1; } diff --git a/EMCAL/AliEMCALUnfolding.cxx b/EMCAL/AliEMCALUnfolding.cxx index f9cd63b8127..f0c7d3172c0 100644 --- a/EMCAL/AliEMCALUnfolding.cxx +++ b/EMCAL/AliEMCALUnfolding.cxx @@ -165,46 +165,47 @@ void AliEMCALUnfolding::MakeUnfolding() { // Unfolds clusters using the shape of an ElectroMagnetic shower // Performs unfolding of all clusters - + if(fNumberOfECAClusters > 0){ if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader") ; //Int_t nModulesToUnfold = fGeom->GetNCells(); - + Int_t numberofNotUnfolded = fNumberOfECAClusters ; Int_t index ; for(index = 0 ; index < numberofNotUnfolded ; index++){ AliEMCALRecPoint * recPoint = dynamic_cast( fRecPoints->At(index) ) ; - - //do we really need it? -// TVector3 gpos; -// Int_t absId = -1; -// recPoint->GetGlobalPosition(gpos); -// fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId); -// if(absId > nModulesToUnfold) -// break ; - - Int_t nMultipl = recPoint->GetMultiplicity() ; - AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ; - Float_t * maxAtEnergy = new Float_t[nMultipl] ; - Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ; - - if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 - if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){ - fRecPoints->Remove(recPoint); - fRecPoints->Compress() ;//is it really needed - index-- ; - fNumberOfECAClusters-- ; - numberofNotUnfolded-- ; - } - } - else{ - recPoint->SetNExMax(1) ; //Only one local maximum - } - - delete[] maxAt ; - delete[] maxAtEnergy ; - } + if(recPoint){ + //do we really need it? + // TVector3 gpos; + // Int_t absId = -1; + // recPoint->GetGlobalPosition(gpos); + // fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId); + // if(absId > nModulesToUnfold) + // break ; + + Int_t nMultipl = recPoint->GetMultiplicity() ; + AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ; + Float_t * maxAtEnergy = new Float_t[nMultipl] ; + Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ; + + if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 + if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){ + fRecPoints->Remove(recPoint); + fRecPoints->Compress() ;//is it really needed + index-- ; + fNumberOfECAClusters-- ; + numberofNotUnfolded-- ; + } + } + else{ + recPoint->SetNExMax(1) ; //Only one local maximum + } + + delete[] maxAt ; + delete[] maxAtEnergy ; + } else AliError("RecPoint NULL"); + } // rec point loop } // End of Unfolding of clusters } @@ -220,10 +221,10 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, Int_t nPar = 3 * nMax ; Float_t * fitparameters = new Float_t[nPar] ; - + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader") ; - + Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; if( !rv ) { // Fit failed, return (and remove cluster? - why? I leave the cluster) @@ -231,17 +232,17 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, delete[] fitparameters ; return kFALSE; } - + // create unfolded rec points and fill them with new energy lists // First calculate energy deposited in each sell in accordance with // fit (without fluctuations): efit[] // and later correct this number in acordance with actual energy // deposition - + Int_t nDigits = iniTower->GetMultiplicity() ; Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells Float_t xpar=0.,zpar=0.,epar=0. ;//center of gravity in cell units - + AliEMCALDigit * digit = 0 ; Int_t * digitsList = iniTower->GetDigitsList() ; @@ -251,73 +252,80 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, Int_t iIeta = 0 ; Int_t iphi = 0 ;//x direction Int_t ieta = 0 ;//z direstion - + Int_t iparam = 0 ; Int_t iDigit = 0 ; - + for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ digit = dynamic_cast( fDigitsArr->At(digitsList[iDigit] ) ) ; - fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); - fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, - iIphi, iIeta,iphi,ieta); - EvalParsPhiDependence(digit->GetId(),fGeom); - - efit[iDigit] = 0.; - iparam = 0; - while(iparam < nPar ){ - xpar = fitparameters[iparam] ; - zpar = fitparameters[iparam+1] ; - epar = fitparameters[iparam+2] ; - iparam += 3 ; - - efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; - } - - } - + if(digit){ + fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, + iIphi, iIeta,iphi,ieta); + EvalParsPhiDependence(digit->GetId(),fGeom); + + efit[iDigit] = 0.; + iparam = 0; + while(iparam < nPar ){ + xpar = fitparameters[iparam] ; + zpar = fitparameters[iparam+1] ; + epar = fitparameters[iparam+2] ; + iparam += 3 ; + + efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; + } + } else AliError("Digit NULL!"); + + }//digit loop + // Now create new RecPoints and fill energy lists with efit corrected to fluctuations // so that energy deposited in each cell is distributed between new clusters proportionally // to its contribution to efit - + Float_t * energiesList = iniTower->GetEnergiesList() ; Float_t ratio = 0 ; - + iparam = 0 ; while(iparam < nPar ){ xpar = fitparameters[iparam] ; zpar = fitparameters[iparam+1] ; epar = fitparameters[iparam+2] ; iparam += 3 ; - + AliEMCALRecPoint * recPoint = 0 ; - + if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters) ; - + //add recpoint (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; recPoint = dynamic_cast( fRecPoints->At(fNumberOfECAClusters) ) ; - fNumberOfECAClusters++ ; - recPoint->SetNExMax((Int_t)nPar/3) ; - - Float_t eDigit = 0. ; - for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = dynamic_cast( fDigitsArr->At( digitsList[iDigit] ) ) ; - fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); - fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, - iIphi, iIeta,iphi,ieta); - EvalParsPhiDependence(digit->GetId(),fGeom); - if(efit[iDigit]==0) continue;//just for sure - ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; - eDigit = energiesList[iDigit] * ratio ; - recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case - } - - } - + + if(recPoint){ + + fNumberOfECAClusters++ ; + recPoint->SetNExMax((Int_t)nPar/3) ; + + Float_t eDigit = 0. ; + for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ + digit = dynamic_cast( fDigitsArr->At( digitsList[iDigit] ) ) ; + if(digit){ + fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, + iIphi, iIeta,iphi,ieta); + EvalParsPhiDependence(digit->GetId(),fGeom); + if(efit[iDigit]==0) continue;//just for sure + ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ; + eDigit = energiesList[iDigit] * ratio ; + recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case + } else AliError("NULL digit"); + }//digit loop + } else AliError("NULL RecPoint"); + }//while + delete[] fitparameters ; delete[] efit ; - + return kTRUE; } @@ -456,106 +464,112 @@ void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, { // Calculates the Chi square for the cluster unfolding minimization // Number of parameters, Gradient, Chi squared, parameters, what to do - + TList * toMinuit = dynamic_cast( gMinuit->GetObjectFit() ) ; - - AliEMCALRecPoint * recPoint = dynamic_cast( toMinuit->At(0) ) ; - TClonesArray * digits = dynamic_cast( toMinuit->At(1) ) ; - // A bit buggy way to get an access to the geometry - // To be revised! - AliEMCALGeometry *geom = dynamic_cast(toMinuit->At(2)); - - Int_t * digitsList = recPoint->GetDigitsList() ; - - Int_t nOdigits = recPoint->GetDigitsMultiplicity() ; - - Float_t * energiesList = recPoint->GetEnergiesList() ; - - fret = 0. ; - Int_t iparam = 0 ; - - if(iflag == 2) - for(iparam = 0 ; iparam < nPar ; iparam++) - Grad[iparam] = 0 ; // Will evaluate gradient - - Double_t efit = 0. ; - - AliEMCALDigit * digit ; - Int_t iDigit ; - - Int_t iSupMod = 0 ; - Int_t iTower = 0 ; - Int_t iIphi = 0 ; - Int_t iIeta = 0 ; - Int_t iphi = 0 ;//x direction - Int_t ieta = 0 ;//z direstion - - - for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { - if(energiesList[iDigit]==0) continue; - - digit = dynamic_cast( digits->At( digitsList[iDigit] ) ); - - geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); - geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, - iIphi, iIeta,iphi,ieta); - EvalParsPhiDependence(digit->GetId(),geom); - - if(iflag == 2){ // calculate gradient - Int_t iParam = 0 ; - efit = 0. ; - while(iParam < nPar ){ - Double_t dx = ((Float_t)iphi - x[iParam]) ; - iParam++ ; - Double_t dz = ((Float_t)ieta - x[iParam]) ; - iParam++ ; - efit += x[iParam] * ShowerShapeV2(dx,dz) ; - iParam++ ; - } - - Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E) - iParam = 0 ; - while(iParam < nPar ){ - Double_t xpar = x[iParam] ; - Double_t zpar = x[iParam+1] ; - Double_t epar = x[iParam+2] ; - - Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) ); - Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; - Double_t rp1 = TMath::Power(dr, fSSPars[1]) ; - Double_t rp5 = TMath::Power(dr, fSSPars[5]) ; - - Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * - (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - - (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + - fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) ); - - //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) - // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; - - Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x - iParam++ ; - Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z - iParam++ ; - Grad[iParam] += shape ; // Derivative over energy - iParam++ ; - } - } - efit = 0; - iparam = 0 ; - - while(iparam < nPar ){ - Double_t xpar = x[iparam] ; - Double_t zpar = x[iparam+1] ; - Double_t epar = x[iparam+2] ; - iparam += 3 ; - efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; - } - - fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ; - // Here we assume, that sigma = sqrt(E) - } - + if(toMinuit){ + AliEMCALRecPoint * recPoint = dynamic_cast( toMinuit->At(0) ) ; + TClonesArray * digits = dynamic_cast( toMinuit->At(1) ) ; + // A bit buggy way to get an access to the geometry + // To be revised! + AliEMCALGeometry *geom = dynamic_cast(toMinuit->At(2)); + + if(recPoint && digits && geom){ + + Int_t * digitsList = recPoint->GetDigitsList() ; + + Int_t nOdigits = recPoint->GetDigitsMultiplicity() ; + + Float_t * energiesList = recPoint->GetEnergiesList() ; + + fret = 0. ; + Int_t iparam = 0 ; + + if(iflag == 2) + for(iparam = 0 ; iparam < nPar ; iparam++) + Grad[iparam] = 0 ; // Will evaluate gradient + + Double_t efit = 0. ; + + AliEMCALDigit * digit ; + Int_t iDigit ; + + Int_t iSupMod = 0 ; + Int_t iTower = 0 ; + Int_t iIphi = 0 ; + Int_t iIeta = 0 ; + Int_t iphi = 0 ;//x direction + Int_t ieta = 0 ;//z direstion + + + for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { + if(energiesList[iDigit]==0) continue; + + digit = dynamic_cast( digits->At( digitsList[iDigit] ) ); + + if(digit){ + geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); + geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower, + iIphi, iIeta,iphi,ieta); + EvalParsPhiDependence(digit->GetId(),geom); + + if(iflag == 2){ // calculate gradient + Int_t iParam = 0 ; + efit = 0. ; + while(iParam < nPar ){ + Double_t dx = ((Float_t)iphi - x[iParam]) ; + iParam++ ; + Double_t dz = ((Float_t)ieta - x[iParam]) ; + iParam++ ; + efit += x[iParam] * ShowerShapeV2(dx,dz) ; + iParam++ ; + } + + Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E) + iParam = 0 ; + while(iParam < nPar ){ + Double_t xpar = x[iParam] ; + Double_t zpar = x[iParam+1] ; + Double_t epar = x[iParam+2] ; + + Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) ); + Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; + Double_t rp1 = TMath::Power(dr, fSSPars[1]) ; + Double_t rp5 = TMath::Power(dr, fSSPars[5]) ; + + Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * + (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - + (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + + fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) ); + + //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) + // - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; + + Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ; // Derivative over x + iParam++ ; + Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ; // Derivative over z + iParam++ ; + Grad[iParam] += shape ; // Derivative over energy + iParam++ ; + } + } + efit = 0; + iparam = 0 ; + + while(iparam < nPar ){ + Double_t xpar = x[iparam] ; + Double_t zpar = x[iparam+1] ; + Double_t epar = x[iparam+2] ; + iparam += 3 ; + efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ; + } + + fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ; + // Here we assume, that sigma = sqrt(E) + } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n"); + } // digit loop + } // recpoint, digits and geom not NULL + }// List is not NULL + } -- 2.43.0