+ } // energy theshold
+ }
+
+}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::MakeUnfolding()
+{
+ // Unfolds clusters using the shape of an ElectroMagnetic shower
+ // Performs unfolding of all EMC/CPV clusters
+
+ // Unfold first EMC clusters
+ if(fNumberOfEmcClusters > 0){
+
+ Int_t nModulesToUnfold = fGeom->GetNModules() ;
+
+ Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
+ Int_t index ;
+ for(index = 0 ; index < numberofNotUnfolded ; index++){
+
+ AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
+ if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
+ break ;
+
+ Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
+ AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
+ Float_t * maxAtEnergy = new Float_t[nMultipl] ;
+ Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
+
+ if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
+ UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+
+ fEMCRecPoints->Remove(emcRecPoint);
+ fEMCRecPoints->Compress() ;
+ index-- ;
+ fNumberOfEmcClusters -- ;
+ numberofNotUnfolded-- ;
+ }
+ else{
+ emcRecPoint->SetNExMax(1) ; //Only one local maximum
+ }
+
+ delete[] maxAt ;
+ delete[] maxAtEnergy ;
+ }
+ }
+ // Unfolding of EMC clusters finished
+
+
+ // Unfold now CPV clusters
+ if(fNumberOfCpvClusters > 0){
+
+ Int_t nModulesToUnfold = fGeom->GetNModules() ;
+
+ Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
+ Int_t index ;
+ for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
+
+ AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
+
+ if(recPoint->GetPHOSMod()> nModulesToUnfold)
+ break ;
+
+ AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
+
+ Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
+ AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
+ Float_t * maxAtEnergy = new Float_t[nMultipl] ;
+ Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
+
+ if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
+ UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+ fCPVRecPoints->Remove(emcRecPoint);
+ fCPVRecPoints->Compress() ;
+ index-- ;
+ numberofCpvNotUnfolded-- ;
+ fNumberOfCpvClusters-- ;
+ }
+
+ delete[] maxAt ;
+ delete[] maxAtEnergy ;
+ }
+ }
+ //Unfolding of Cpv clusters finished
+
+}
+
+//____________________________________________________________________________
+Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
+{
+ // Shape of the shower (see PHOS TDR)
+ // If you change this function, change also the gradient evaluation in ChiSquare()
+
+ //for the moment we neglect dependence on the incident angle.
+
+ Double_t r2 = x*x + z*z ;
+ Double_t r4 = r2*r2 ;
+ Double_t r295 = TMath::Power(r2, 2.95/2.) ;
+ Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
+ return shape ;
+}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
+ Int_t nMax,
+ AliPHOSDigit ** maxAt,
+ Float_t * maxAtEnergy)
+{
+ // Performs the unfolding of a cluster with nMax overlapping showers
+
+ Int_t nPar = 3 * nMax ;
+ Float_t * fitparameters = new Float_t[nPar] ;
+
+ Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
+
+ if( !rv ) {
+ // Fit failed, return and remove cluster
+ iniEmc->SetNExMax(-1) ;
+ delete[] fitparameters ;
+ return ;
+ }
+
+ // create ufolded 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 = iniEmc->GetMultiplicity() ;
+ Float_t * efit = new Float_t[nDigits] ;
+ Float_t xDigit=0.,zDigit=0. ;
+ Float_t xpar=0.,zpar=0.,epar=0. ;
+ Int_t relid[4] ;
+ AliPHOSDigit * digit = 0 ;
+ Int_t * emcDigits = iniEmc->GetDigitsList() ;
+
+ TVector3 vIncid ;
+
+ Int_t iparam ;
+ Int_t iDigit ;
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+ digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+ efit[iDigit] = 0;
+
+ iparam = 0 ;
+ while(iparam < nPar ){
+ xpar = fitparameters[iparam] ;
+ zpar = fitparameters[iparam+1] ;
+ epar = fitparameters[iparam+2] ;
+ iparam += 3 ;
+// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
+// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
+ efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
+ }
+ }
+
+ // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
+ // so that energy deposited in each cell is distributed betwin new clusters proportionally
+ // to its contribution to efit
+
+ Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
+ Float_t ratio ;
+
+ iparam = 0 ;
+ while(iparam < nPar ){
+ xpar = fitparameters[iparam] ;
+ zpar = fitparameters[iparam+1] ;
+ epar = fitparameters[iparam+2] ;
+ iparam += 3 ;
+// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
+
+ AliPHOSEmcRecPoint * emcRP = 0 ;
+
+ if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
+
+ if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
+ fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
+
+ (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
+ fNumberOfEmcClusters++ ;
+ emcRP->SetNExMax((Int_t)nPar/3) ;
+ }
+ else{//create new entries in fCPVRecPoints
+ if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
+ fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
+
+ (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
+ fNumberOfCpvClusters++ ;
+ }
+
+ Float_t eDigit ;
+ for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+ digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
+ ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
+ eDigit = emcEnergies[iDigit] * ratio ;
+ emcRP->AddDigit( *digit, eDigit ) ;
+ }
+ }