/*
 *  Program: Analyze Lambda Electroproduction
 *
 *  Author : Steve Barrow 
 *
 *  Date :   Summer, 1998. 
 *
*/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <signal.h>
#include <stdlib.h>
#include <errno.h>
#include <ntypes.h>
#include <bostypes.h>
#include <clas_cern.h>
#include <kinematics.h>
#include <map_manager.h>
#include <sc.h>
#include <ec.h>
#include <utility.h>
#include <particleType.h>
#include <pid.h>
#include "accept.h"

/* ----------- Function prototypes ------------------ */

int pe_sbarrow_(float beam_energy);
void bh_sbarrow_(int runno, float beam_energy, float Q_live);
void pl_1520_sbarrow_(float beam_energy);

#define P_kaon_min 0.283
#define P_kaon_max 2.017
#define N_KAON_SPA_P 102
int kaon_in[N_KAON_SPA_P];
int kaon_index;
#include "rho.h"

/* --------------------------------------------------- */

void pl_1520_sbarrow_(float beam_energy) 
{
  int ind_1,ind_2,ind_3,ind_4,ind_5,ind_6;
  float theta_kaon_minus;

  hf1(1,RAD2DEG*get_rel_angle_spb(Kaon,Proton_SPB),1.0);

  if ((mag_spb(Kaon.space)>P_kaon_min)&&(mag_spb(Kaon.space)<P_kaon_max)) {
    kaon_index=get_index(mag_spb(Kaon.space),P_kaon_min,P_kaon_max,N_KAON_SPA_P);
    kaon_in[kaon_index] = kaon_in[kaon_index] +1;
  }

  theta_rel=get_rel_plane2(Electron_SPB.space,Kaon.space,Missing.space); 
  lambda = Missing;
  test_cm = Lab2CMFrame(beam_energy,Electron_SPB,Kaon,Proton_SPB);
  Target_cm = Lab2CMFrame(beam_energy,Electron_SPB,Kaon,Target);
  k_minus_cm = Lab2CMFrame(beam_energy,Electron_SPB,Kaon,v4sub(Missing,Proton_SPB));
  lambda_cm = Lab2CMFrame(beam_energy,Electron_SPB,Kaon,Missing);
  
  t = v4sub(Target,Missing);

  cosines = get_lambda_decay(lambda_cm, test_cm);
  
  kaon_cm = Lab2CMFrame(beam_energy,Electron_SPB,Kaon,Kaon);
  theta_cm_kaon = (float)RAD2DEG*acos(kaon_cm.space.z/\
              sqrt(pow(kaon_cm.space.x,2)+pow(kaon_cm.space.z,2)));
  

  beta_lambda_cm.x=lambda_cm.space.x/lambda_cm.t;
  beta_lambda_cm.y=lambda_cm.space.y/lambda_cm.t;
  beta_lambda_cm.z=lambda_cm.space.z/lambda_cm.t;

  k_minus_t = v4boost(beta_lambda_cm, k_minus_cm);
  Target_t =  v4boost(beta_lambda_cm, Target_cm);
  k_minus_t = rotate_CM_SPB(k_minus_t,Target_t.space);

  if (k_minus_t.space.y > 0.0) {
     phi_kaon_minus=(float)RAD2DEG*acos((double)(k_minus_t.space.x/\
      sqrt(pow(k_minus_t.space.x,2)+pow(k_minus_t.space.y,2))));
  }
  if (k_minus_t.space.y < 0.0) {
     phi_kaon_minus=360-(float)RAD2DEG*acos((double)(k_minus_t.space.x/\
      sqrt(pow(k_minus_t.space.x,2)+pow(k_minus_t.space.y,2))));
  }

  /* end of sept. 16th stuff */

    if ((mag_spb(Electron_SPB.space)>E_CUT)) {
      theta_kaon_minus = get_theta_rel(Missing,Missing_prot);
      ind_1 = get_index(W,W_MIN,W_MAX,N_W);
      ind_2 = get_index(Q2,Q2_MIN,Q2_MAX,N_Q2);
      ind_3 = get_index(theta_rel,PHI_REL_MIN,PHI_REL_MAX,N_PHI_REL);
      ind_4 = get_index(180 - theta_cm_kaon,THETA_CM_MIN,THETA_CM_MAX,N_THETA_CM);
      ind_5 = get_index(theta_kaon_minus,THETA_LAM_MIN,THETA_LAM_MAX,N_THETA_LAM);
      ind_6 = get_index(phi_kaon_minus,PHI_LAM_MIN,PHI_LAM_MAX,N_PHI_LAM);
      if ((current_run ==1)&&(mag_spb(Electron_SPB.space)>E_CUT)&&(Q2<Q2_MAX)&&(Q2>Q2_MIN)&&\
	(W>W_MIN)&&(W<W_MAX)) {
	accept[ind_1][ind_2][ind_3][ind_4][ind_5][ind_6]=accept[ind_1]\
              [ind_2][ind_3][ind_4][ind_5][ind_6]+1; 
      }

    }

} 

/* end pl_1520_sbarrow_() */
/* --------------------------------------------------- */



void bh_sbarrow_( int runno, float beam_energy, float Q_live){
  int i, id, bins;
  float min, max, v;
  char title[100];
  char runmes[50];
  int sec;
  
  sprintf(runmes,"Run %d", runno);

  /* booking  histos */
  
  sprintf(title, "rel-angle kplus, proton, all events", runmes);
  hbook1(1, title, 60, 0.0, 80.0, 0);
  sprintf(title, "M versus phi, proton tracks", runmes);
  hbook2(2, title, 20, 0.4, 0.6, 80, 0.0, 360.0, 0); 
  sprintf(title, "Q2, e-k events", runmes);
  hbook1(3, title, 40, 0.5, 3.5, 0);
  sprintf(title, "W, e-k events", runmes);
  hbook1(4, title, 30, 0.0, 3.0, 0);
  sprintf(title, "Q2 versus W, e-k events", runmes);
  hbook2(5, title, 20, 0.5, 3.5, 30, 1.0, 3.0, 0); 
  sprintf(title, "M versus P hadron tracks", runmes);
  hbook2(10, title, 40, 0.0, 1.2, 40, 0.0, 4.0, 0); 
  sprintf(title, "e-k missing mass", runmes);
  hbook1(11, title, 40, 1.0, 2.0, 0);

}


int pe_sbarrow_(float beam_energy){

  FILE *fp_part = NULL;
  
  clasHEAD_t *HEAD = getBank(&bcs_, "HEAD");
  clasTBID_t *TBID = NULL;
  clasPART_t *PART = NULL;
  clasTRGS_t *TRGS = NULL;
  clasECHB_t *ECHB = NULL;
  clasSCPB_t *SCPB = NULL;
  clasSCR_t  *SCR = NULL;
  clasSCRC_t  *SCRC = NULL;
  clasVERT_t *VERT = NULL;
  clasEC_t *EC = NULL;
  clasEC01_t *EC01 = NULL;  /* added oct. 3rd */
  clasCC0_t *CC = NULL;  /*added oct 6th */
  vert_t *testvert;
  float testvert_x,testvert_y,testvert_z;
  int jj_ec, jj_cc, sect_high;   /* added oct. 3rd, 6th */
  float E_trig_sector; /* added oct. 3rd */
  float E_trig_sector2, E_trig_highest;
  int EC_id; 
  float EC_energy; /* added oct. 3rd */  
  int head_test1;
  int head_test2;
  int head_test3;

  dec97 = 0;
  feb98 = 0;
  bad_sc = 0;
  npart = 0;
  nelec = 0;
  nprot = 0;
  nkaon = 0;
  npion = 0;
  nkplus = 0;
  nkminus = 0;
  npiplus = 0;
  npiminus = 0;
  
  Nep = 0;
  Q2 = 0.0;
  W = 0.0;

  /* dropAllBanks(&bcs_,"ECHB"); 
  ec_evnt_();  */     

  /* dropAllBanks(&bcs_,"CC01");
  dropAllBanks(&bcs_,"SC1 SCR SCRCEC01");   
  for(sec=1; sec <= 6; sec++) {   
    make_CC01_bank(sec);     
    make_SC1_bank(sec);   
    make_SCR_bank(sec,"TDPL");
    make_SCRC_bank(sec);
    make_EC01_bank(sec);
  }


  dec97 = 0;
  feb98 = 0;
  bad_sc = 0;
  npart = 0;
  nelec = 0;
  nprot = 0;
  nkaon = 0;
  npion = 0;
  nkplus = 0;
  nkminus = 0;
  npiplus = 0;
  npiminus = 0;
  
  Nep = 0;
  
  make_CL01_bank();
  dropAllBanks(&bcs_,"ECHB");  
  ec_evnt_(); 
  
  dropAllBanks(&bcs_,"TBIDPART");  */
  /* make_TBID_group(1); 
  make_PART_group(1); */
  
  TBID = getGroup(&bcs_, "TBID",1);  
  PART = getGroup(&bcs_, "PART",1); 
  /* printf("HEAD:%i, TBID:%i, PART:%i \n.", HEAD, TBID, PART);
  printf("HEAD:%i, TBID:%i, PART:%i \n.", HEAD, TBID, PART); */

  /* head_test1 = 0;
  head_test2 = 1;
  head_test3 = -1;
    if (head_test1) { 
      printf("HEAD1:%i\n .",head_test1 );
    }
    if (head_test2) { 
      printf("HEAD2:%i\n .",head_test2 );
    }
    if (head_test3) { 
      printf("HEAD3:%i\n .",head_test3 );
    } */

  
  /* ------------------ RUN FLAGS ------------------------------ */
  
      if(HEAD){
	class_evt = HEAD->head[0].evtclass;
      }	
  /* ------------------------------------------------------------ */
  

  if(HEAD && TBID && PART ){

   tracks = 0; 
    Beam.space.x = Beam.space.y = 0.;
    Beam.space.z = beam_energy;
    Beam.t = beam_energy;
    Target.space.x = Target.space.y = Target.space.z = 0.0;
    Target.t = PROTON_MASS;
    beta_cm.x =0;
    beta_cm.y =0;
    beta_cm.z = (Beam.space.z)/(Beam.t+Target.t);
    
    Total = v4add(Beam,Target);
    
    
    for (j = 0; j < PART->bank.nrow; j++){
      npart++;
      part_id = PART->part[j].pid;
      index_tbid = PART->part[j].trkid - 1;
      stat_sc = TBID->tbid[index_tbid].sc_stat;
      stat_ec = TBID->tbid[index_tbid].ec_stat;
      stat_cc = TBID->tbid[index_tbid].cc_stat;
      stat_trk = PART->part[j].qtrk;
      sector = TBID->tbid[index_tbid].sec;
      sc_nn = TBID->tbid[index_tbid].sc_id-1;
     if (PART->part[j].q != 0) tracks++; 
      paddle = 0;
      if (SCRC = getGroup(&bcs_, "SCRC",sector)) {
	if (TBID->tbid[j].sc_stat)
	  paddle = SCRC->scrc[sc_nn].id; 
      }
     
      /* bad_sc_spb = 0;
      if ((paddle < 51) && (0< paddle)){	
	bad_sc_spb = bad_sc_arr[sector][paddle];
      } */

      /* if((part_id==3) && (bad_sc_spb != -13) && (paddle < 40) && stat_ec == (GOOD_MATCH)){ */
      if((part_id==3) && (paddle < 40) && stat_ec == (GOOD_MATCH)){
	/* ELECTRON STUFF */

	Elec_p = mag_spb(PART->part[j].p.space);
	nelec++;
	Electron_SPB = PART->part[j].p;
	/* Vertex_elec = PART->part[j].vert; */
	jj = PART->part[j].trkid;

	make_vert();
	VERT = getBank(&bcs_,"VERT");
	testvert = VERT->vert;
	testvert_x = (testvert + index_tbid)->vert.x;
	testvert_y = (testvert + index_tbid)->vert.y;
	testvert_z = (testvert + index_tbid)->vert.z;

	  
	hf1(2,VERT->vert[PART->part[j].trkid].vert.z,1.0); 
	  
	if ((fabs(testvert_z-VERT_OFF)<VERT_ELEC)&&(testvert_x*testvert_x+\
                                         testvert_y*testvert_y < 0.85)) {
	  vert_good = 100; 
	}  /* end of vertex cut for electron track */
	

	Photon = v4sub(Beam, Electron_SPB);
	Resonance = v4add(Photon, Target);
	Particle = PART->part[j].p;
	
	Nu = beam_energy - Electron_SPB.t;    
	Q2 = -v4magsq(Photon);
	W  =  v4mass(Resonance);  
	
	Xbrk = Q2/(2.0*PROTON_MASS*Nu); 
	Ybrk = Nu/beam_energy;
	
	Particle_x = Particle.space.x;
	Particle_y = Particle.space.y;

	if (Particle_y > 0.0) {
	  phi_elec=(float)RAD2DEG*acos((double)(Particle_x/\
                          sqrt(pow(Particle_x,2)+pow(Particle_y,2))));
	}
	if (Particle_y < 0.0) {
	  phi_elec=360-(float)RAD2DEG*acos((double)(Particle_x/\
                          sqrt(pow(Particle_x,2)+pow(Particle_y,2))));
	}

	sect_elec = sector;
	sc_nn = TBID->tbid[index_tbid].sc_id; 
      }
    }
    
    
    if(W > 1.5 && Q2 > 0.5 &&  stat_ec == (GOOD_MATCH) ){


      printf("**************************************\n");

      nkaon = 0;
      for (jj = 0; jj < PART->bank.nrow; jj++){
	
	part_id = PART->part[jj].pid;
	index_tbid = PART->part[jj].trkid - 1;
	stat_sc = TBID->tbid[index_tbid].sc_stat;
	stat_ec = TBID->tbid[index_tbid].ec_stat;
	stat_cc = TBID->tbid[index_tbid].cc_stat;
	stat_trk = PART->part[jj].qtrk;
	sector = TBID->tbid[index_tbid].sec;
	sc_nn = TBID->tbid[index_tbid].sc_id -1;
	bad_sc_spb = 0;
	paddle = 0;
	if (SCRC = getGroup(&bcs_, "SCRC",sector)) {
	  if (TBID->tbid[index_tbid].sc_stat)
	    paddle = SCRC->scrc[sc_nn].id; 
	}

	/* if ((paddle < 51) && (0< paddle)){
	  bad_sc_spb = bad_sc_arr[sector][paddle];
	} */
	
	/* ---------------- VERTEX ----------------------------------*/
	
	/* Vertex_part = PART->part[jj].vert; */
	
	/*-------------------------------------------------------------*/
	
	
	Particle = PART->part[jj].p;
	charge = PART->part[jj].q;
	index_tbid = PART->part[jj].trkid - 1;
	Part_beta = TBID->tbid[index_tbid].beta;
	Particle_x = Particle.space.x;
	Particle_y = Particle.space.y;
	Particle_z = Particle.space.z;
	Particle_e = Particle.t;
	Part_mom = sqrt(Particle_x*Particle_x+Particle_y*Particle_y+Particle_z*Particle_z);
	
       	if(Part_beta > 0.0 && Part_beta < 1.1 && charge != 0){
	  Part_mass = Part_mom * sqrt(1.0 - Part_beta*Part_beta)/Part_beta;
	
	  if (charge == 1) hf2(10,Part_mass, Part_mom,1.0);


	  /* if((bad_sc_spb != -13) && (paddle < 40) && (vert_good = 100)){ */
	  if((paddle < 40) && (vert_good = 100)){ 
	  
	    /* KAONS */

	    /*              0.42                  0.58     */
	    if((Part_mass > 0.42) && (Part_mass < 0.58) && (charge == 1) && stat_sc == (GOOD_MATCH)){
	      p_kaon = Part_mom;
              beta_kaon2 = Part_beta;
	      kaon_mass = Part_mass; 
	      sect_kaon = sector;
	
	      hf1(3,Q2,1.0);
	      hf1(4,W,1.0);
	      hf2(5,Q2,W,1.0);


	      testvert_x = (testvert + index_tbid)->vert.x;
	      testvert_y = (testvert + index_tbid)->vert.y;
	      testvert_z = (testvert + index_tbid)->vert.z;
	      
              if((fabs(testvert_x)<2.5)&&(fabs(testvert_y)<2.5)&&\
                 (fabs(testvert_z-VERT_OFF)<VERT_KAON)&&(vert_good==100)) {
		
		nkaon=100;
		
		Kaon = PART->part[jj].p;
		temp = v4add(Kaon,Electron_SPB);
		Missing = v4sub(Total,temp);

		hf1(11,v4mass(Missing),1.0);
		
		t = v4sub(Target,Missing);

		if (Particle_y > 0.0) {
		  phi_kaon=(float)RAD2DEG*acos((double)(Particle_x/\
                           sqrt(pow(Particle_x,2)+pow(Particle_y,2))));
		}
		if (Particle_y < 0.0) 
		  phi_kaon=360-(float)RAD2DEG*acos((double)(Particle_x/\
                           sqrt(pow(Particle_x,2)+pow(Particle_y,2))));
		}

		hf2(2,Part_mass,phi_kaon,1.0); 
	      } /* end of cut on kaon vertex */
	    } /* end of cut on kaon mass */
	  } /* end of cut on good TOF paddles */   
	} /* cut on part_beta */
      } /* end of loop over part bank */
    } /* end of W, Q2 cut */
    
     if (nkaon == 100) {
       
       for (j = 0; j < PART->bank.nrow; j++){
	
	part_id = PART->part[j].pid;
	index_tbid = PART->part[j].trkid - 1;
	stat_sc = TBID->tbid[index_tbid].sc_stat;
	stat_ec = TBID->tbid[index_tbid].ec_stat;
	stat_cc = TBID->tbid[index_tbid].cc_stat;
	stat_trk = PART->part[j].qtrk;
	sector = TBID->tbid[index_tbid].sec;
	sc_nn = TBID->tbid[index_tbid].sc_id - 1;
	bad_sc_spb = 0;
	paddle = 0;
	if (SCRC = getGroup(&bcs_, "SCRC",sector)) {
	  if (TBID->tbid[index_tbid].sc_stat)
	    paddle = SCRC->scrc[sc_nn].id; 
	}
	
	
	/* ------------------ VERTEX ----------------------------------*/
	
	/* Vertex_part = PART->part[j].vert; */

	/*-------------------------------------------------------------*/
	
	/* if(bad_sc_spb != -13){ */
	  
	  Particle = PART->part[j].p;
	  charge = PART->part[j].q;
	  Part_beta = TBID->tbid[index_tbid].beta;
	  Particle_x = Particle.space.x;
	  Particle_y = Particle.space.y;
	  Particle_z = Particle.space.z;
	  Particle_e = Particle.t;
	  Part_mom = sqrt(Particle_x*Particle_x+Particle_y*Particle_y+Particle_z*Particle_z);
	  
	  if(Part_beta > 0.0 && Part_beta < 1.1 && charge != 0){
	    Part_mass = Part_mom * sqrt(1.0 - Part_beta*Part_beta)/Part_beta;
	    
	    /* PROTON */
	    /*             0.905                1.050  */
	    if((Part_mass > 0.905) && (Part_mass < 1.050) && (charge == 1) &&\
               (stat_sc == (GOOD_MATCH))&&(tracks==3)){ 
	      beta_proton2 = Part_beta;
	      proton_mass = Part_mass;
	      sector = TBID->tbid[index_tbid].sec;
	      sect_prot = sector;

	      if ((nkaon == 100)&&(vert_good==100)) { 
		nprot++;
		Proton_SPB = PART->part[j].p;
		Missing_prot = v4sub(Missing,Proton_SPB);
		m_miss = v4magsq(Missing_prot);

		if((m_miss<0.75) && (m_miss> -0.1) ) {

		  testvert_x = (testvert + index_tbid)->vert.x;
		  testvert_y = (testvert + index_tbid)->vert.y;
		  testvert_z = (testvert + index_tbid)->vert.z;
		  
		  if((fabs(testvert_x)<2.5)&&(fabs(testvert_y)<2.5)&&\
                     (fabs(testvert_z-VERT_OFF)<VERT_PROT)&&(vert_good==100)) {
		    

		    if((1.46<v4mass(Missing))&&(1.58>v4mass(Missing))&&(W>1.95)\
                       &&(sqrt(m_miss)>0.42)&&(sqrt(m_miss)<0.58)\
		       &&mag_spb(Electron_SPB.space)>E_CUT&&(side_low !=100)&&(side_high !=100)) { 
		      p_proton = Part_mom;
		      S = v4add(Missing,Kaon);
		      pl_1520_sbarrow_(beam_energy);
     		    }
		  } /* end of cut on proton vertex */
		}
	      } /* end of nkaon==100, vert_good==100 */
	    } /*end of cut on proton mass */
	  } /*end of beta cut */
      } /*end of loop over part bank */
     p_kaon = 0;
     p_proton = 0;
     nkaon=0;
     npion=0;
     vert_good = 0;
     side_low = 0;
     side_high = 0;

     }
  /* } */
  /* fclose(fp_part); */
  return(1);  
}

/* end file */


<Hit "Back" button on your browser>