Using decays to ϕ -meson pairs, the inclusive production of charmonium states in b -hadron decays is studied with pp collision data corresponding to an integrated luminosity of 3.0fb-1 , collected by the LHCb experiment at centre-of-mass energies of 7 and 8 TeV. Denoting by BC≡B(b→CX)×B(C→ϕϕ) the inclusive branching fraction of a b hadron to a charmonium state C that decays into a pair of ϕ mesons, ratios RC2C1≡BC1/BC2 are determined as Rηc(1S)χc0=0.147±0.023±0.011 , Rηc(1S)χc1=0.073±0.016±0.006 , Rηc(1S)χc2=0.081±0.013±0.005 , Rχc0χc1=0.50±0.11±0.01 , Rχc0χc2=0.56±0.10±0.01 and Rηc(1S)ηc(2S)=0.040±0.011±0.004 . Here and below the first uncertainties are statistical and the second systematic. Upper limits at 90% confidence level for the inclusive production of X(3872), X(3915) and χc2(2P) states are obtained as Rχc1X(3872)<0.34 , Rχc0X(3915)<0.12 and Rχc2χc2(2P)<0.16 . Differential cross-sections as a function of transverse momentum are measured for the ηc(1S) and χc states. The branching fraction of the decay Bs0→ϕϕϕ is measured for the first time, B(Bs0→ϕϕϕ)=(2.15±0.54±0.28±0.21B)×10-6 . Here the third uncertainty is due to the branching fraction of the decay Bs0→ϕϕ , which is used for normalization. No evidence for intermediate resonances is seen. A preferentially transverse ϕ polarization is observed. The measurements allow the determination of the ratio of the branching fractions for the ηc(1S) decays to ϕϕ and pp¯ as B(ηc(1S)→ϕϕ)/B(ηc(1S)→pp¯)=1.79±0.14±0.32 .