#include "mwdc_defs.h"
#include "std_funcs.h"
#include <iostream>
#include <fstream>

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

  Set of functions which to be accessed by analysis scripts

  11/1/05 B. Craver

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

plane_type paired_plane( plane_type plane )
{
  Int_t pair_counter;
  plane_type paired_plane;

  paired_plane = -1;

  if (plane == 0)   paired_plane = 1;
  if (plane == 1)  paired_plane = 0;

  if (plane == 3)   paired_plane = 4;
  if (plane == 4)  paired_plane = 3;

  if (plane == 5)   paired_plane = 6;
  if (plane == 6)  paired_plane = 5;

  if (plane == 8)   paired_plane = 9;
  if (plane == 9)  paired_plane = 8;

  if (plane == 10)   paired_plane = 11;
  if (plane == 11)  paired_plane = 10;

  if (plane == 13)   paired_plane = 14;
  if (plane == 14)  paired_plane = 13;

  return paired_plane;
}


/* Returns the orientation of the wires of the current plane with the following format:
0 - U plane
1 - X plane
2 - V plane
*/
Int_t get_pair_type(plane_type plane)
{
  Int_t pair_type;
  Int_t chamber_number = chamber[plane];

  if ( (plane == upair[chamber_number][0]) || (plane == upair[chamber_number][1]) )
    pair_type = 0;
  if ( (plane == xpair[chamber_number][0]) || (plane == xpair[chamber_number][1]) )
    pair_type = 1;
  if ( (plane == vpair[chamber_number][0]) || (plane == vpair[chamber_number][1]) )
    pair_type = 2;

  return (pair_type);
}



// Calculates "Sigma" = (Avg. U + Avg. V) - Sqrt(3)(Avg. X)
Double_t sigma_calculator(Int_t chamber_number, Int_t number_of_hits[], Double_t hit_wires[][MAX_NUMBER_OF_HITS])
{
  Double_t average_u[NUMBER_OF_CHAMBERS], average_x[NUMBER_OF_CHAMBERS], average_v[NUMBER_OF_CHAMBERS]; 
  Double_t average_uprime[NUMBER_OF_CHAMBERS], average_xprime[NUMBER_OF_CHAMBERS], average_vprime[NUMBER_OF_CHAMBERS]; 
  Double_t sigma = -100.;
  Double_t sigma_offset = 0.;
  plane_type plane;

  if (chamber_number == 0)  sigma_offset = 36.; 
  if ( (chamber_number == 1) || (chamber_number == 2) )  sigma_offset = 50.; 

  plane = upair[chamber_number][0];
  average_u[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  plane = upair[chamber_number][1];
  average_uprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  
  plane = xpair[chamber_number][0];
  average_x[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  plane = xpair[chamber_number][1];
  average_xprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  
  plane = vpair[chamber_number][0];
  average_v[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  plane = vpair[chamber_number][1];
  average_vprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
  
  sigma = (average_u[chamber_number]+ average_uprime[chamber_number] + average_v[chamber_number] + average_vprime[chamber_number]) - (1.7320508) * (average_x[chamber_number] + average_xprime[chamber_number]);

  sigma -= sigma_offset; 
  if (chamber_number == 1) sigma /= 2.; // Value was artificially doubled by adding primed values
  return (sigma);
}


/* Calculates "Sigma" = (Avg. U + Avg. V) - Sqrt(3)(Avg. X) using less than six planes
   The planes to be ignored are specified and are handled in the following way:
   If only one of a pair of planes contains valid hits those hits are duplicated for the other plane
*/
Double_t sigma_calculator_by_plane(Int_t chamber_number, Int_t number_of_hits[], Double_t hit_wires[][MAX_NUMBER_OF_HITS], plane_type ignored_planes[])
{
  Double_t average_u[NUMBER_OF_CHAMBERS], average_x[NUMBER_OF_CHAMBERS], average_v[NUMBER_OF_CHAMBERS]; 
  Double_t average_uprime[NUMBER_OF_CHAMBERS], average_xprime[NUMBER_OF_CHAMBERS], average_vprime[NUMBER_OF_CHAMBERS]; 
  Double_t sigma = -1.;
  Double_t sigma_offset;
  plane_type plane;

  if (chamber_number == 0)  sigma_offset = 36.; 
  if ( (chamber_number == 1) || (chamber_number == 2) )  sigma_offset = 50.; 

  //cout << "chamber number : " <<chamber_number<<endl;
  plane = upair[chamber_number][0];
  //cout << "ignored_planes[" <<plane<< "] : " <<ignored_planes[plane]<<endl;
  if (ignored_planes[plane] != 1) // First plane is not ignored
    {
      average_u[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
      plane = upair[chamber_number][1];
      if (ignored_planes[plane] != 1) // Second plane is not ignored
	{
	  average_uprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	}
      else { // Second plane on pair is ignored: duplicate hits from first pair
	average_uprime[chamber_number] = average_u[chamber_number];
      }
    }
  else { // First plane on pair is ignored: duplicate hits from second pair
    plane = upair[chamber_number][1];
    if (ignored_planes[plane] != 1) // Second plane is not ignored
      {
	average_uprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	average_u[chamber_number] = average_uprime[chamber_number];
      }
    else { // Second plane on pair is ignored: declare averages to be 0
      average_u[chamber_number] = 0.; average_uprime[chamber_number] = 0.;
    }
  }
  //cout << "U : " <<average_u[chamber_number]<< " : U' : " <<average_uprime[chamber_number]<<endl;
  
  plane = xpair[chamber_number][0];
  //cout << "ignored_planes[" <<plane<< "] : " <<ignored_planes[plane]<<endl;
  if (ignored_planes[plane] != 1) // First plane is not ignored
    {
      average_x[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
      //cout << "1 - x = " <<average_x[chamber_number]<<endl;
      plane = xpair[chamber_number][1];
      if (ignored_planes[plane] != 1) // Second plane is not ignored
	{
	  average_xprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	  //cout << "1 - x' = " <<average_xprime[chamber_number]<<endl;
	}
      else { // Second plane on pair is ignored: duplicate hits from first pair (whether or not they exist)
	average_xprime[chamber_number] = average_x[chamber_number];
	//cout << "2 - x' = " <<average_xprime[chamber_number]<<endl;
     }
    }
  else { // First plane on pair is ignored: duplicate hits from second pair (whether or not they exist)
    plane = xpair[chamber_number][1];
    if (ignored_planes[plane] != 1) // Second plane is not ignored
      {
	average_xprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	//cout << "3 - x' = " <<average_xprime[chamber_number]<<endl;
	average_x[chamber_number] = average_xprime[chamber_number];
        //cout << "3 - x = " <<average_x[chamber_number]<<endl;
    }
    else { // Second plane on pair is ignored: declare averages to be 0
      average_x[chamber_number] = 0.; average_xprime[chamber_number] = 0.;
      //cout << "4 - x = " <<average_x[chamber_number]<<endl;
      //cout << "4 - x' = " <<average_xprime[chamber_number]<<endl;
    }
  }
  //cout << "X : " <<average_x[chamber_number]<< " : X' : " <<average_xprime[chamber_number]<<endl;
  
  plane = vpair[chamber_number][0];
  //cout << "ignored_planes[" <<plane<< "] : " <<ignored_planes[plane]<<endl;
  if (ignored_planes[plane] != 1) // First plane is not ignored
    {
      average_v[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
      plane = vpair[chamber_number][1];
      if (ignored_planes[plane] != 1) // Second plane is not ignored
	{
	  average_vprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	}
      else { // Second plane on pair is ignored: duplicate hits from first pair (whether or not they exist)
	average_vprime[chamber_number] = average_v[chamber_number];
      }
    }
  else { // First plane on pair is ignored: duplicate hits from second pair (whether or not they exist)
    plane = vpair[chamber_number][1];
    if (ignored_planes[plane] != 1) // Second plane is not ignored
      {
	average_vprime[chamber_number] = average(number_of_hits[plane], hit_wires[plane]);
	average_v[chamber_number] = average_vprime[chamber_number];
      }
    else { // Second plane on pair is ignored: declare averages to be 0
      average_v[chamber_number] = 0.; average_vprime[chamber_number] = 0.;
    }
  }
  //cout << "V : " <<average_v[chamber_number]<< " : V' : " <<average_vprime[chamber_number]<<endl;
  
  //cout << "average_u[" <<chamber_number<< "] + average_uprime[" <<chamber_number<< "] = " <<average_u[chamber_number] + average_uprime[chamber_number]<<endl;
  //cout << "average_v[" <<chamber_number<< "] + average_vprime[" <<chamber_number<< "] = " <<average_v[chamber_number] + average_vprime[chamber_number]<<endl;
  //cout << "average_x[" <<chamber_number<< "] + average_xprime[" <<chamber_number<< "] = " <<average_x[chamber_number] + average_xprime[chamber_number]<<endl;
  sigma = (average_u[chamber_number]+ average_uprime[chamber_number] + average_v[chamber_number] + average_vprime[chamber_number]) - ((1.7320508) * (average_x[chamber_number] + average_xprime[chamber_number]) );
  
  sigma -= sigma_offset; 
  if (chamber_number == 1) sigma /= 2.; // Value was artificially doubled by adding primed values
  //cout << "sigma = " <<sigma<<endl;
  return (sigma);
}


// Defines "single cluster" as a single wire or 2 adjacent wires firing on the plane
is_single_cluster(plane_type plane, Int_t number_of_hits, Double_t hit_wires[][MAX_NUMBER_OF_HITS])
{
  Int_t hit_counter = 0;  Int_t wire_counter = 0; Int_t number_unique_wires = 0; Int_t last_wire = -1;
  Int_t single_cluster = 0;
  
  for (hit_counter = 0; hit_counter < number_of_hits; hit_counter++)
    {
      wire_counter = hit_wires[plane][hit_counter];
      if (wire_counter != last_wire)
	{
	  number_unique_wires++;
	  if (number_unique_wires == 2)
	    {
	      if ( (wire_counter > last_wire + 1) || (wire_counter < last_wire - 1) )
		{
		  return (0);  // Hit on non-adjacent wire
		}
	    }
	}
      last_wire = wire_counter;
    }

  if ( (number_unique_wires == 1) || (number_unique_wires == 2) )
    {
      return (1);
    }
  else {
    return (0);
  }
}

// One or both planes of similar wire orientation contain a single firing wire
is_single_hit_pair(plane_type plane, Int_t number_of_hits[], Double_t hit_wires[][MAX_NUMBER_OF_HITS])
{
  plane_type paired_plane = paired_plane(plane); plane_type plane_counter;
  Int_t wire_counter = 0;
  Double_t new_wire = 0.;  Double_t old_wire = 0.;

  // Verify that neither of the two planes contain hits on multiple wires
  plane_counter = plane;
  for (wire_counter = 0; wire_counter < number_of_hits[plane_counter]; wire_counter++)
    {
      new_wire = hit_wires[plane_counter][wire_counter];
      if ( (new_wire != old_wire)  && (wire_counter != 0) )
	{
	  //cout << "1 - plane_counter : " <<plane_counter<< " multiple wires" <<endl;
	  return (0);
	}
      old_wire = new_wire;
    }
  plane_counter = paired_plane;
  for (wire_counter = 0; wire_counter < number_of_hits[plane_counter]; wire_counter++)
    {
      new_wire = hit_wires[plane_counter][wire_counter];
      if ( (new_wire != old_wire)  && (wire_counter != 0) )
	{
	  //cout << "2 - plane_counter : " <<plane_counter<< " multiple wires" <<endl;
	  return (0);
	}
      old_wire = new_wire;
    }

  // Allow at most one of the two planes to contain no hits
  if ( (number_of_hits[plane] >= 1) && (number_of_hits[paired_plane] >= 1) )
    return (1);
  if ( (number_of_hits[plane] >= 1) && (number_of_hits[paired_plane] == 0) )
    return (1);
  if ( (number_of_hits[plane] == 0) && (number_of_hits[paired_plane] >= 1) )
    return (1);

  return (0);
}


// Allows at most one plane on the chamber to have zero hits
// and requires that all other planes have exactly one wire firing
is_single_hit_chamber(plane_type plane, Int_t number_of_hits[], Double_t hit_wires[][MAX_NUMBER_OF_HITS])
{
  Int_t chamber_number = chamber[plane]; Int_t zero_counter = 0;
  plane_type plane_counter = 0;
  Int_t wire_counter = 0;
  Double_t new_wire = 0.;  Double_t old_wire = 0.;

  // Verify that no planes contain hits on multiple wires
  for (plane_counter = 0; plane_counter < NUMBER_OF_PLANES; plane_counter++)
    {
      if (chamber[plane_counter] == chamber[plane])
	{
	  for (wire_counter = 0; wire_counter < number_of_hits[plane_counter]; wire_counter++)
	    {
	      new_wire = hit_wires[plane_counter][wire_counter];
	      if ( (new_wire != old_wire) && (wire_counter != 0) )
		{
		  //cout << "plane_counter : " <<plane_counter<< " multiple wires"<<endl;
		  return (0);
		}
	      old_wire = new_wire;
	    }
	}
    }
  
  // Count the planes containing no hits
  for (plane_counter = 0; plane_counter < NUMBER_OF_PLANES; plane_counter++)
    {
      if (chamber[plane_counter] == chamber[plane])
	{
	  if (number_of_hits[plane_counter] == 0)
	    zero_counter++;
	  //cout << "zero_counter = " <<zero_counter<<endl;
	}
    }
  
  if (zero_counter > 1)
    {
      //cout << "returning 0 due to zero's = " <<zero_counter<<endl;
      return (0);
    }
  else {
    return (1);
  }
}

// Requires that all planes on chamber have exactly one wire firing
is_single_hit_chamber_all(plane_type plane, Int_t number_of_hits[], Double_t hit_wires[][MAX_NUMBER_OF_HITS])
{
  Int_t chamber_number = chamber[plane];
  plane_type plane_counter = 0;
  Int_t wire_counter = 0;
  Double_t new_wire = 0.;  Double_t old_wire = 0.;

  // Verify that no planes contain hits on multiple wires
  for (plane_counter = 0; plane_counter < NUMBER_OF_PLANES; plane_counter++)
    {
      if (chamber[plane_counter] == chamber[plane])
	{
	  for (wire_counter = 0; wire_counter < number_of_hits[plane_counter]; wire_counter++)
	    {
	      new_wire = hit_wires[plane_counter][wire_counter];
	      if ( (new_wire != old_wire) && (wire_counter != 0) )
		{
		  //cout << "plane_counter : " <<plane_counter<< " multiple wires"<<endl;
		  return (0);
		}
	      old_wire = new_wire;
	    }
	}
    }
  
  // Count the planes containing no hits
  for (plane_counter = 0; plane_counter < NUMBER_OF_PLANES; plane_counter++)
    {
      if (chamber[plane_counter] == chamber[plane])
	{
	  if (number_of_hits[plane_counter] == 0)
	    return (0);
	}
    }
  
  return (1);
}