my_star/crab/binnings/crab_bindefs_qinv.cpp

/* Set ABSOLUTE_MAXMOM so that program won't bother calculating correlations
   when rel. mom. is greater than ABSOLUTE_MAXMOM. This is ignored if
   "MIXED_PAIRS_FOR_DENOM" is defined in crab.cpp below */
#define ABSOLUTE_MAXMOM 800

/* Four Examples are presented here (A-D). For each binning, modifications
   must be made in four locations in this file as illustrated. One can
   add as many binnings as desired. Binnings are of two types:
   1.) one-dimensional, where the answer will be given as C(x).
   2.) three-dimensional, where the answer will be given as C(x,y,z). */

/* One of the included routines is 'decompose(p1,p2,vs,
   long_comoving_boost,x_along_p,outwards_boost,
   &kmag,&kx,&ky,&kz)'. This routine returns the 3-dimensions of the
   relative momentum. ( If one only uses Qinv, one can skip calling this
   function since 'mom' equals 'Qinv' and is already fed into 'binner'. )
   When one calls decompose, three additional arguments,
   long_comoving_boost,x_along_p,outwards_boost, are used to determine the
   reference frame, (both relativistically and rotationally).
   1.) If one chooses 'long_comoving_boost' equal to one,
   one boosts along the beam axis to the frame where pz1+pz2=0.
   If 'long_comoving_boost' is not equal to one, the longitudinal frame
   will be set by the value of 'vs'.
   2.) If one chooses 'x_along_p' equal to one,
   one performs a rotation, in the reaction plane defined by the pairs
   momentum, such that the x-axis is parallel to the total momentum.
   3.) If one chooses outwards_boost equal to one,
   one also boost along the direction of p1+p2 to the pair's c.o.m. frame. */

/* Example A: Bin the correlation function according to Qinv */
/* Example B: Bin the correlation function according to Qinv
   in 10 different pt bins */
/* Example C: Bin the correlation function according to |Q|, where Q is
   determined in the longitudinally comoving frame. Make three binnings for
   three different directional cuts. */
/* Example D: Save the correlation function in a 3-d binning, according
   to Qx, Qy, and Qz. */

/* ********************************************************* */

/* First define the binnings you will do. */

/* A. */
onedbin_class *qinv;

/* B. */
/*
#define NPTCUTS 10
#define DELPT 40
onedbin_class *qinv_ptcut[NPTCUTS];
*/

/* C */
/*
onedbin_class *q_outwards;
onedbin_class *q_sidewards;
onedbin_class *q_beam;
*/

/* D. */
threedbin_class *qinv3d;

/* E */
#define NMTCUTS 5
threedbin_class *qinv3d_mt[NMTCUTS];
const float mTarr[NMTCUTS] = {0.2, 0.4, 0.6, 0.7, 1.};
const float partMass = 0.13957018; // pi 


/* ********************************************************* */

/* Initialize binning objects with desired mesh sizes and momentum ranges */

void bininit(void){
  int ipt,nkmax,nkxmax,nkymax,nkzmax;
  int nptmax;
  double maxmom,kxmax,kymax,kzmax;
  double pt;
  char title[80];


  /* A */
  nkmax=80;
  maxmom=800.0;
  qinv=new onedbin_class(nkmax,maxmom,"qinv.dat");

  /* B */
  /*
  nkmax=40;
  maxmom=400.0;
  for(ipt=0;ipt<NPTCUTS;ipt++){
    pt=(0.5+ipt)*maxmom/NPTCUTS;
    pt=floor(pt+0.5);
    sprintf(title,"qinv_pt%03g.dat",pt);
    qinv_ptcut[ipt]=new onedbin_class(nkmax,maxmom,title);
  }
  */
  
  /* C */
  /*
  nkmax=40;
  maxmom=400.0;
  q_outwards=new onedbin_class(nkmax,maxmom,"q_outwards.dat");
  q_sidewards=new onedbin_class(nkmax,maxmom,"q_sidewards.dat");
  q_beam=new onedbin_class(nkmax,maxmom,"q_beam.dat");
  */
  
  /* D */
  nkxmax = nkymax = nkzmax = 80;
  kxmax = kymax = kzmax = 800.0;
  qinv3d=new threedbin_class(nkxmax,nkymax,nkzmax,
			     kxmax,kymax,kzmax,"qinv3d.dat");

  /* E */
  nkxmax = nkymax = nkzmax = 80;
  kxmax = kymax = kzmax = 800.0;
  for (int iMt = 0; iMt < NMTCUTS; iMt++) {
    sprintf(title, "qinv3d_%i.dat", iMt);
    qinv3d_mt[iMt]=new threedbin_class(nkxmax, nkymax, nkzmax,
				       kxmax, kymax, kzmax, title);
  }
  

}

/* ********************************************************* */

/* Bin |phi|^2 according to p1 and p2 */
/* Use the routine "decompose" to get projections of rel. mom. if desired */

void binner(int inumden,double mom,double corr,double *p1,double *p2){
  double pt,kx,ky,kz,kmag;
  int ipt;
  /* These 3 variables define conv.s used by decompose, 0=no, 1=yes */
  int long_comoving_boost,x_along_p,outwards_boost;
  double vs=0.0;

  /* A */
  qinv->bincorr(inumden,mom,corr);

  /* B */
  pt=sqrt(pow(p1[1]+p2[1], 2.)+pow(p1[2]+p2[2], 2.));
  /*
  ipt=(int)floor(pt/DELPT);
  if(ipt<NPTCUTS){
    qinv_ptcut[ipt]->bincorr(inumden,mom,corr);
  }
  */
  
  /* C */
  /*
  long_comoving_boost=1;
  x_along_p=0;
  outwards_boost=0;
  decompose(p1,p2,vs,
	    long_comoving_boost,x_along_p,outwards_boost,
	    &kmag,&kx,&ky,&kz);
  if(fabs(ky)<50.0 && fabs(kz)<50.0) q_outwards->bincorr(inumden,kmag,corr);
  if(fabs(kx)<50.0 && fabs(kz)<50.0) q_sidewards->bincorr(inumden,kmag,corr);
  if(fabs(kx)<50.0 && fabs(ky)<50.0) q_beam->bincorr(inumden,kmag,corr);
  */
  
  /* D */
  long_comoving_boost=1;
  x_along_p=0;
  outwards_boost=0;
  decompose(p1,p2,vs,
	    long_comoving_boost,x_along_p,outwards_boost,
	    &kmag,&kx,&ky,&kz);
  qinv3d->bincorr(inumden,kx,ky,kz,corr);
  
  /* E */
  float kT = pt*1.e-3/2.;
  float mT = sqrt(kT*kT + partMass*partMass);
  for (int iMt = 0; iMt < NMTCUTS; iMt++) {
    if (mT <= mTarr[iMt]) {
      long_comoving_boost=1;
      x_along_p=0;
      outwards_boost=0;
      decompose(p1,p2,vs,
		long_comoving_boost,x_along_p,outwards_boost,
		&kmag,&kx,&ky,&kz);
      qinv3d_mt[iMt]->bincorr(inumden,kx,ky,kz,corr);
      break;
    }
  }
  
}

void printstatements(int onlyzerob){
  int ipt;

  /* A */
  qinv->printcorr(onlyzerob);

  /* B */
  /*
  for(ipt=0;ipt<10;ipt++){
    qinv_ptcut[ipt]->printcorr(onlyzerob);
  }
  */
  
  /* C */
  /*
  q_outwards->printcorr(onlyzerob);
  q_sidewards->printcorr(onlyzerob);
  q_beam->printcorr(onlyzerob);
  */

  /* D */
  qinv3d->printcorr(onlyzerob);

  /* E */
  for (int iMt = 0; iMt < NMTCUTS; iMt++) {
    qinv3d_mt[iMt]->printcorr(onlyzerob);
  }
}