-
- delete digitsC ;
-
- AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
-}
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::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<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
-
- TVector3 gpos;
- Int_t absId;
- 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
- UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
- fRecPoints->Remove(recPoint);
- fRecPoints->Compress() ;
- index-- ;
- fNumberOfECAClusters-- ;
- numberofNotUnfolded-- ;
- }
- else{
- recPoint->SetNExMax(1) ; //Only one local maximum
- }
-
- delete[] maxAt ;
- delete[] maxAtEnergy ;
- }
- }
- // End of Unfolding of clusters
-}
-
-//____________________________________________________________________________
-Double_t AliEMCALClusterizerv1::ShowerShape(Double_t x, Double_t y)
-{
- // Shape of the shower
- // If you change this function, change also the gradient evaluation in ChiSquare()
-
- Double_t r = sqrt(x*x+y*y);
- Double_t r133 = TMath::Power(r, 1.33) ;
- Double_t r669 = TMath::Power(r, 6.69) ;
- Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
- return shape ;
-}
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
- Int_t nMax,
- AliEMCALDigit ** 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] ;
-
- if (fGeom==0)
- AliFatal("Did not get geometry from EMCALLoader") ;
-
- Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
- if( !rv ) {
- // Fit failed, return and remove cluster
- iniTower->SetNExMax(-1) ;
- delete[] fitparameters ;
- return ;
- }
-
- // 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] ;
- Double_t xDigit=0.,yDigit=0.,zDigit=0. ;
- Float_t xpar=0.,zpar=0.,epar=0. ;
-
- AliEMCALDigit * digit = 0 ;
- Int_t * digitsList = iniTower->GetDigitsList() ;
-
- Int_t iparam ;
- Int_t iDigit ;
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
- fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
- efit[iDigit] = 0;
-
- iparam = 0 ;
- while(iparam < nPar ){
- xpar = fitparameters[iparam] ;
- zpar = fitparameters[iparam+1] ;
- epar = fitparameters[iparam+2] ;
- iparam += 3 ;
- 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 between new clusters proportionally
- // to its contribution to efit
-
- Float_t * energiesList = iniTower->GetEnergiesList() ;
- Float_t ratio ;
-
- 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) ;
-
- (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
- recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
- fNumberOfECAClusters++ ;
- recPoint->SetNExMax((Int_t)nPar/3) ;
-
- Float_t eDigit ;
- for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
- fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
-
- ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
- eDigit = energiesList[iDigit] * ratio ;
- recPoint->AddDigit( *digit, eDigit ) ;
- }
- }
-
- delete[] fitparameters ;
- delete[] efit ;
-
-}
-
-//_____________________________________________________________________________
-void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
- Double_t & fret,
- Double_t * x, Int_t iflag)
-{
- // Calculates the Chi square for the cluster unfolding minimization
- // Number of parameters, Gradient, Chi squared, parameters, what to do
-
- TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
-
- AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ;
- TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
- // A bit buggy way to get an access to the geometry
- // To be revised!
- AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
-
- Int_t * digitsList = recPoint->GetDigitsList() ;
-
- Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
-
- Float_t * energiesList = recPoint->GetEnergiesList() ;
-
- fret = 0. ;
- Int_t iparam ;
-
- if(iflag == 2)
- for(iparam = 0 ; iparam < nPar ; iparam++)
- Grad[iparam] = 0 ; // Will evaluate gradient
-
- Double_t efit ;
-
- AliEMCALDigit * digit ;
- Int_t iDigit ;
-
- for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
-
- digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
-
- Double_t xDigit=0 ;
- Double_t zDigit=0 ;
- Double_t yDigit=0 ;//not used yet, assumed to be 0
-
- geom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
-
- if(iflag == 2){ // calculate gradient
- Int_t iParam = 0 ;
- efit = 0 ;
- while(iParam < nPar ){
- Double_t dx = (xDigit - x[iParam]) ;
- iParam++ ;
- Double_t dz = (zDigit - x[iParam]) ;
- iParam++ ;
- efit += x[iParam] * ShowerShape(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 = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
- Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
- Double_t r133 = TMath::Power(dr, 1.33);
- Double_t r669 = TMath::Power(dr,6.69);
- 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 * (xpar - xDigit) ; // Derivative over x
- iParam++ ;
- Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // 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 * ShowerShape(xDigit - xpar,zDigit - zpar) ;
- }
-
- fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
- // Here we assume, that sigma = sqrt(E)
- }
-}
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
-{
- // Print clusterizer parameters
-
- TString message("\n") ;