18a19,60
> 
> 
> static double IMP_PROTEIN_DENSITY=0.826446; // 0.826446=1/1.21   don't know exactly where that constant comes from
> static double protein_density_ = IMP_PROTEIN_DENSITY;
> double get_used_protein_density(){
>   return protein_density_;
>  };
> 
> double set_used_protein_density(double densityValue){
>   IMP_USAGE_CHECK(densityValue > 0,
>         "A protein density value is most probably a positive value");
>   protein_density_=densityValue;
>   return protein_density_;
>  }
> 
> double set_used_protein_density(ProteinDensityReference densityReference){
>   double density=0.0;
>   switch (densityReference) {
>     case HARPAZ:
>       density = IMP_PROTEIN_DENSITY;
>       break;
>     case ANDERSSON :
>       density = 0.7347;  // Andersson and Hovmšller (1998) Theoretical 1.22 g/cm3
>       break;
>     case TSAI:
>       density = 0.84309; // Tsai et al. (1999) Theoretical 1.40 g/cm3
>       break;
>     case QUILLIN:
>       density = 0.86116; // Quillin and Matthews (2000) Theoretical 1.43 g/cm3
>       break;
>     case SQUIRE:
>       density = 0.82503; // Squire and Himmel (1979) and Gekko and Noguchi (1979) Experimental 1.37
>       break;
>     default :
>       // should warn IMP_WARNING("unknown density reference... density set to default");
>       density = IMP_PROTEIN_DENSITY;
>    }
>   density=set_used_protein_density(density);
>   return density;
>  }
> 
> 
725a768,785
> Float compute_threshold_for_approximate_volume(DensityMap* d, Float volumeDesired) {
>   Float voxelVolume     = d->get_spacing()*d->get_spacing()*d->get_spacing();
>   long numVoxelsNeeded  = volumeDesired / voxelVolume;
>   long mapSizeInVoxels  = d->get_number_of_voxels();
>   std::vector<emreal> data(mapSizeInVoxels);
>   for (long l=0;l<mapSizeInVoxels;++l) {
>     // This ugly loop should be replaced by a memcopy if only I had access to d->data_
>     data[l]=(d->get_value(l));
>    }
>   std::sort(data.begin(),data.end());
>   emreal threshold = data[mapSizeInVoxels-numVoxelsNeeded];
>   return static_cast<Float>(threshold);
>  }
> 
> Float compute_threshold_for_approximate_mass(DensityMap* d, Float desiredMass){
>   Float desiredVolume = desiredMass / get_used_protein_density();
>   return compute_threshold_for_approximate_volume(d,desiredVolume);
> }
728c788
<   long counter=0;//number of voxles above the threshold
---
>   long counter=0;//number of voxels above the threshold
734c794,795
<   return d->get_spacing()*counter/1.21;
---
>   Float s=d->get_spacing();
>   return  s*s*s*counter*get_used_protein_density();