In order to decrease the tails of heavy ion peak I processed 20 day data sample using the same cuts as for proton selection:
1) multiplicity (number of crystals with energy above the threshold) is calculated for whole CAL layer (not within one tower)
2) multiplicity == 1 in the layer of crystal being calibrated
3) multiplicity <= 2 in all other layers parallel to the one being calibrated
4) multiplicity <=4 in all other layers perpendicular to the one being calibrated (to allow big angles)
5) extention of tracker track passes through top and bottom of the calibrated crystal
6) tkr1NumHits>5
7) 4-range readout
The 2-d plot of GCR spectrum vs theta angle is show on the following plot: 
I've made four slices of this for different theta angle bands and fit B,C,N and O peaks on each slice:
0< theta < 10 deg
25 deg < theta < 35 deg
40 deg < theta < 50 deg
55 deg < theta < 65 deg
Using peak position for theta < 10 deg as normalization I found the deviation for linearity as a function of crystal energy (not corrected for angle):
While there is some systematics ~ 0.5% in curves for all peaks, we can conclude that deviation from linearity between 300 MeV and 1700 MeV is better than 1%.
Systematics is probably due low energy tail seen in simulation and on data histograms above - could be corrected using simulated shape.
The method could be extended down to Be peak (200 MeV), but not to Li peak because ratio peakBe/peakLi ~2 and Li is absent at small angles.
After adding statistics, we could try to use Ne and Mg (2.0 GeV)
I'm developping alternative method, which could be applied to all peaks, including p, He, and Li.
Method is based on energy split between adjacent crystals.