Detect the correct distribution from a small sample size by using fitdistrplus in R

The simplest version of the issue that I am looking for help is:

I am generating some random numbers with Gamma distribution and fit these random numbers to different distributions (Lognormal, Weibull, Exp, Gamma), but unfortunately the AIC obtained from Gamma is not always the minimum AIC, I appreciate if someone could help me to find an approach to detect Gamma after fitting even for a small sample size. I am using fitdist library.

Longer version: I am generating some random numbers with Gamma distribution. In order to do that I need three parameters:Number of random numbers, shape and scale. I am using pre-calculated shape and scale and for the number of random numbers, this varies from 40 to 120. I know that if shape gets close to 1 it is difficult to distinguish exp, gamma and weibull as stated here, so I am trying to keep shape far away to 1. Sample size is another matter and if I increase my sample size the results will be much better, but I have to keep my sample size small. I am trying to detect Gamma after fitting with a high accuracy but seems it is not possible. I am thinking of changing the method from mle to qme or something else but not sure which one I shall go for. I have tried a few of them but no success as I am not a statistician. Anothe issue that I considered is not only detect the best fitted distribution by the lowest aic, but also by some other parameters such as std error of fitting, but no success. I appreciate any help especially in simple terms:). Please let me know if you need more information.

So this is my code in R:

library(fitdistrplus)
require(distr)

shapegoriginal=0.769230769230769
scalegoriginal=78
numberofrandomnumbers=60
numbeoftrial=100
counter_AIC_fitw=0;
counter_AIC_fitl=0;
counter_AIC_fitgamma=0;
counter_AIC_fite=0;
out <- matrix(NA, nrow=numbeoftrial, ncol=13)
for(i in 1:numbeoftrial)
{

    nn=(rgamma(numberofrandomnumbers, shape = shapegoriginal, scale = scalegoriginal))
    fite=fitdist (nn ,'exp')
    lambda=fite[1]$estimate[1]

    fitl=fitdist (nn ,'lnorm')
    meanl=fitl[1]$estimate[1] 
    sdl=fitl[1]$estimate[2]   

    fitw=fitdist (nn  ,'weibull')
    shape=fitw[1]$estimate[1] 
    scale=fitw[1]$estimate[2]

    fitgamma=fitdist (nn ,'gamma')
    shapeg=fitgamma[1]$estimate[1] 
    scaleg=1/fitgamma[1]$estimate[2]

    AIC_fitw=summary(fitw)$aic
    AIC_fitl=summary(fitl)$aic
    AIC_fitgamma=summary(fitgamma)$aic
    AIC_fite=summary(fite)$aic

    min_AIC=min(AIC_fitw,AIC_fitl,AIC_fitgamma,AIC_fite)

    if(min_AIC==AIC_fitw){counter_AIC_fitw=counter_AIC_fitw+1 }
    if(min_AIC==AIC_fitl){counter_AIC_fitl=counter_AIC_fitl+1}
    if(min_AIC==AIC_fitgamma){counter_AIC_fitgamma=counter_AIC_fitgamma+1}
    if(min_AIC==AIC_fite){counter_AIC_fite=counter_AIC_fite+1}
        out[i,]=c(i,lambda,meanl,sdl,shape,scale,shapeg,scaleg,AIC_fitw,AIC_fitl,AIC_fitgamma,AIC_fite,min_AIC)

}
print('#when Weibull detected')
print(counter_AIC_fitw)
print('#when Lognormal detected')
print(counter_AIC_fitl)
print('#when Gamma detected')
print(counter_AIC_fitgamma)
print('#when Exp detected')
print(counter_AIC_fite)
colnames(out)=c('i','lambda','meanl','sdl','shape','scale','shapeg','scaleg','AIC_fitw','AIC_fitl','AIC_fitgamma','AIC_fite','min_AIC')

out