I finally realized my idea to compare the layer energy depositions for trajectories at the same angles ( with 1/cos(theta) ~1.1 - 1.6) for 2 cases:
1) when trajectory goes completely through one crystal
2) when trajectory goes through two adjacent crystals and signal per crystal is twice smaller.
To avoid zero suppression problem I excluded the events with pathlength in any crystal smaller than 7 mm (corresponding to 3 MeV ).
As in case 2 the trajectory crosses the gap between crystals, the signal is smaller than in case 1, but this decrease has simple connection to the trajectory angle along the crystal ( dx = exitPosX-entryPosX), in fact dE/E = -gap/dx, so the linear fit of dE with linear function p0+p1*(1/dx) has p1 = E*gap and p0 = true proton peak position (corrected for existance of the gap).
The results of the procedure is the following:
for case 1 proton peak position doesn't depend on dx and average value is 10.40 +- 0.02 MeV
for case 2 proton peak indeed vary linearly with 1/dx and linear fit has p1 = -10 which corresponds to crystal to crystal gap = 1 mm (which is rather close to reality) and p0 = 10.69 +- 0.03 MeV, so the nonlinearity between 7 MeV and 14 MeV is ~3 %, but going to opposite direction: relative gain at 7 MeV is 3 % BIGGER than at 14 MeV.
The proton peak in single crystal energy depositions (standard calibration procedure, when crystal is inclined along the crystal direction) is 10.457+-0.003 MeV and consistent with case 1 (10.40 MeV).
I'll follow this pass and will find the same way the nonlinearity near He peak (comparing gains at 30 MeV and at 60 MeV).
To compare 15 MeV with 30 MeV I'll try to use the slow protons (you remember they have a peak at 20-30 MeV in pathlength corrected energy), as in my current procedure I compare strictly the particles with the same angle and energy, just having different position accross the crystals. 
=====================================================================================================
...