ftp.nice.ch/pub/next/science/mathematics/HippoDraw.2.0.s.tar.gz#/HippoDraw/Hippo.bproj/Zshape.m

This is Zshape.m in view mode; [Download] [Up]

/* Zshape.m  	by Paul Kunz	January 1993
 * Subclass of PFunction to manage a Z shape function
 *
 * Copyright (C)  1993  The Board of Trustees of
 * The Leland Stanford Junior University. All Rights Reserved.
 */

#import "Zshape.h"
#import <objc/Storage.h>
#import <math.h>

const char Zshape_h_rcsid[] = ZSHAPE_H_ID;
const char Zshape_m_rcsid[] = "$Id: Zshape.m,v 1.2 1993/02/01 22:14:59 pfkeb Exp $";

/* The form of this function was taken from Jordan Nash's thesis, 
 * SLAC-Rep-356, p 22.
 */

@implementation Zshape : PFunction

static double phi(double cos_beta, double nu)
{ 
    return M_PI*nu*sin((1-nu)*acos(cos_beta))/
            (sin(M_PI*nu)*sin(acos(cos_beta)));
}

static double Z_shape(double x, double binW, double *par)
{
    double y, gamma, a, cos_beta, t=0.108, cross_section;
    double xsect = par[0];
    double mass  = par[1];
    double width = par[2];

    gamma = width/mass;
    y = pow( (x/mass), 2. ) - 1.;
    a = sqrt( (y*y + gamma*gamma*pow((1+y),2.))/(1+gamma*gamma) );
    cos_beta = ( y + gamma*gamma*(1+y))/(a*(1+gamma*gamma));
    cos_beta *= -1.;
	
    cross_section = xsect*(1+0.75*t)*gamma*gamma/(1+gamma*gamma);
    cross_section *= (1+y)*pow(a,t-2)*phi(cos_beta,t) 
                     -     pow(a,t-1)*t*phi(cos_beta,1+t)/(1-t);

// add in the correction for hard photon radiation :
// MOD 900723 WSL multiply correction by peak cross section before adding. 

    cross_section -= xsect*t*width*( atan(2/gamma)-atan(2*(mass-x)/gamma) )/x ;
	
    return cross_section;
}

- init
{
    [super init];
    [[[self setTitle:"Z_shape"] setFunctionPtr:Z_shape] setNumberArgs:3];
    [self registerFunc];
    [self setArgName:"Peak_x-sect" at:0];
    [self setArgName:"Z_mass" at:1];
    [self setArgName:"Z_width" at:2];
    return self;
}
- setInitialValues
{
    fitParm             *parm;
    graphtype_t		type;
    
    [super setInitialValues];
    if ( !disp ) return self;
    
    type = h_getDispType(disp);
    
    parm = [variedParms elementAt:0];
    if ( type == HISTOGRAM ) {
	parm->value       = disp->bins.totals[1][0];
	parm->lower_limit = 0.0;
	parm->upper_limit = 2.0 * parm->value; 
	parm->step_size   = 0.01 * parm->upper_limit;
    } else {
        parm->lower_limit = disp->yAxis.low;
        parm->upper_limit = 4.0 * disp->yAxis.high;
        parm->value       = 0.5 * parm->upper_limit;
	parm->step_size   = 0.01 * (disp->yAxis.high - disp->yAxis.low);
    }
    
    parm = [variedParms elementAt:1];
    parm->lower_limit = disp->xAxis.low;
    parm->upper_limit = disp->xAxis.high;
    parm->step_size   = 0.01 * (disp->xAxis.high - disp->xAxis.low);
    parm->value       = parm->lower_limit 
                        + 0.50 * (disp->xAxis.high - disp->xAxis.low);
    
    parm = [variedParms elementAt:2];
    parm->upper_limit = 0.2 * (disp->xAxis.high - disp->xAxis.low);
    parm->lower_limit = 0.02 * parm->upper_limit;
    parm->step_size   = 0.01 * (parm->upper_limit - parm->lower_limit);
    parm->value = 0.5 * (parm->upper_limit - parm->lower_limit );
    return self;
}

@end  
 

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.