//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
///  TCHIGraphCsI  implementation
///  contains methods to calculate position of point respect to a segment 
///  adapted from KaliVeda class KVIDLine, KVIDGrid, KVIDZAGrid and other classes. 
///  e.d.f. 2011   v. 1.1
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
#include "TCHIGridCsI.h"

void TCHIGraphCsI::Init()
{
//initialisation
 fZMax = 0;
 fZMaxLine = 0;
 finfi = finf = fsup = fsups = -1;
 fdinf = fdsup = fdinfi = fdsups = 0.;
 fwinf = fwsup = fwinfi = fwsups = 0.;
 fZinfi = fZinf = fZsup = fZsups = fAinfi = fAinf = fAsup = fAsups = 0;
 fDistanceClosest = -1.;
 fClosest = fLsups = fLsup = fLinf = fLinfi = 0;
 fIdxClosest = -1;
 fIMFlineadded = false;
 fAint = fZint = 0;
 fICode = kICODE0;  //just to initialize
} 

// based on KVIDZAGrid::Initialize() and KVIDGrid::Initialize();
void TCHIGraphCsI::Initialize()
{
 //Sort the list 
 fgrid->Sort();
 // Calculate the line width 
 // Method KVIDZAGrid::CalculateLineWidths() not implemented.  
 // Zmax should be Z of last line in sorted list
 fZMaxLine = (TGraph *)fgrid->Last();
 if(fZMaxLine) 
  fZMax = GetZ(fZMaxLine);
 else 				
  fZMax=0; 
}
//Returns a pointer to line corresponding to particle (Z,A) and its index
TGraph *TCHIGraphCsI::GetZALine(int Z, int A, int &index)
{
 index = 0;
 bool found=false;
 if(fgrid->IsEmpty()) {
  return 0;
 } 

 int nlines = fgrid->GetEntries();
 TGraph *line = dynamic_cast <TGraph *>(fgrid->At(index));
 while(index<nlines) {
  if(GetZ(line) == Z && GetA(line) == A) {
   found = true;
   break;
  }
  line = dynamic_cast <TGraph * >(fgrid->At(++index));
 }

 if(!found) {
  index = -1;
  return 0;
 }
 return line;
}

//Return the grid Z value scanning the name 
int TCHIGraphCsI::GetZ(TGraph *pg)
{
 char s[3];
 const char *name = pg->GetName();
 const char *p=name;
 //special case 8Be [saved as Z=2,A=8, name=8Be]
 string sb(name);
// cout << "sb = " << sb << endl;
 int pos = sb.find("8Be");
 if(pos!=-1) {
  return 4;
 }
 if(*p=='Z') {
  p++;
  int j=0;
  while(*p != 'A') {
   if(*p>='0' && *p<='9') {
    s[j] = *p;
//    cout << j << " " << s[j] << endl; 
    j++;
   }
   p++;
  }
  s[j] = '\0';
//   cout << "name = " << name << " *** " << j << " " << s[0] << s[1] << " ***" ;
  if(atoi(s)>10){
   cout << "dentro TCHIGraphCsI.cc linea 90 GetZ return= " << atoi(s) << endl;
   getchar();
  }
  else {
//  cout << " OK GetZ return= " << atoi(s) << endl;
  }
  return atoi(s);
 }
 else {
  return -1;
 }
}

//Return the grid A value scanning the name
int TCHIGraphCsI::GetA(TGraph *pg)
{
 char s[3];
 const char *name = pg->GetName();
 const char *p=name;
 while(*p && *p!='A')
  p++;
 if(*p!='A')return -1;
 int j=0;
 p++;
 while(*p!='_') {
  if(*p>='0' && *p<='9') {
   s[j] = *p;
   j++;
  }
  p++;
 }
 return atoi(s);
} 


//Compute the closest distance of approach from point px,py to this line.
//The units of px, py are the same as the coordinates of the graph
//
//   ** Adapted from kaliveda:  KVIDLine::DistanceToLine **
//
//The distance from the point to each segment is calculated in turn.
//If the point lies between the endpoints of at least one segment, the returned
//distance will be the closest of these segments. In this case i_near gives the
//index of the first point of the closest segment.
//
//If not then the distance is the smallest distance between the endpoint
//of a segment and the point. In this case i_near gives the index of the
//closest endpoint.
//If the closest segment-endpoint is not one of the endpoints of the whole line,
//the returned distance is positive.
//On the other hand, if the closest part of the line to the point is one of the
//two endpoints of the line, the returned distance is NEGATIVE.
Double_t TCHIGraphCsI::DistanceToLine(TGraph *gr, Double_t px, Double_t py, Int_t &i_near)
{
 Double_t distance, dist2;
 Int_t i;
 Double_t d;
   
 distance = dist2 = 9999.;
 Int_t i_nearest_point = 0, inear1 = 0, inear2 = 0;
 //distance = closest approach to a segment with point lying between endpoints
 //dist2 = closest approach to a segment with point outside endpoints
 //if distance == 9999 at end, then dist2 is used
 Int_t NPoints = gr->GetN();
 Double_t *fx = gr->GetX();
 Double_t *fy = gr->GetY();

 for (i = 0; i < NPoints - 1; i++) {
  d = DistanceToLine(px, py, fx[i], fy[i], fx[i + 1], fy[i + 1], i_nearest_point);
  if (d >= 0.) {
   if (d < distance) {
    distance = d;
    inear1 = i;
   }
  } 
  else {
  //point not between endpoints of segment
   if (-d < dist2) {
    dist2 = -d;
    inear2 = i + i_nearest_point;
   }
  }
 }

 i_near = inear1;
 if (distance < 9999.)
  return distance;
 i_near = inear2;
 //closest point is an 'internal' point of line - POSITIVE distance
 if (inear2 > 0 && inear2 < (NPoints - 1))
  return dist2;
 //closest point is one of endpoints of line - NEGATIVE distance
 return -dist2;
}

//Given a line segment defined by endpoints (xp1,yp1) and (xp2,yp2) find the
//shortest distance between point (px,py) and the line.
//
//                                                     M
//      P1 (xp1,yp1) x---------------------------------------x P2 (xp2,yp2)
//                                                      |
//                                                      |
//                                                      |
//                                                      |
//                                                      |
//                                                      x
//                                                                P (px,py)
//This is simply the magnitude of the component of vector P1->P  (or P2->P)
//perpendicular to P1->P2.
//If the point is indeed between the two endpoints as shown, then
//P1M + P2M = P1P2
//If not, P1M + P2M > P1P2. In this case the closest distance is that to the nearest
//of the two endpoints, but the value returned is negative; i_nearest_point gives the index (0 or 1)
//of the nearer endpoint
//   ** from kaliveda:  KVIDLine::DistanceToLine [no modifications] **
Double_t TCHIGraphCsI::DistanceToLine(Double_t px, Double_t py,
                                      Double_t xp1, Double_t yp1,
                                      Double_t xp2, Double_t yp2,
                                      Int_t &i_nearest_point)
{
 TVector2 P1(xp1, yp1), P2(xp2, yp2), P(px, py);
 TVector2 P1P2 = P2 - P1;
 TVector2 P1P = P - P1;
 TVector2 P2P = P - P2;
 TVector2 MP = P1P.Norm((P1P2));
 i_nearest_point = 0;

 //check point lies between endpoints
 Double_t sum = (P1P - MP).Mod() + (P2P - MP).Mod();
 //cout << "sum = "<<sum<<"  P1P2 = " << P1P2.Mod() << " difference = " << sum-P1P2.Mod() << endl;
 if (TMath::Abs(sum - P1P2.Mod()) > 1.e-06) {
 //cout << "take endpoint" << endl;
  if (P1P.Mod() < P2P.Mod()) {
   i_nearest_point = 0;
   return -(P1P.Mod());
  }
  i_nearest_point = 1;
  return -(P2P.Mod());
 }
 //cout << "point is close to segment" << endl;
 return MP.Mod();
}

//Given a line segment defined by endpoints (xp1,yp1) and (xp2,yp2) test the
//  direction of (px,py) with respect to the line in the (x,y) plane.
//
//                                                        M
//      P1 (xp1,yp1) x---------------------------------------x P2 (xp2,yp2)
//                                                         |
//                                                         |
//                                                         |
//                                                         |
//                                                         |
//                                                         x
//                                                                P (px,py)
//This is given by the angle phi between vector MP (component of vector P1->P  (or P2->P)
//perpendicular to P1->P2) and the +ve x-axis.
//Point is left of line if 90<phi<270
//Point is right of line if 270<phi<90
//Point is above line if 0<phi<180
//Point is below line if 180<phi<360
//The option string is "left", "right", "above" or "below", and the result is either kTRUE or kFALSE
//depending on the relative position of the point and the line.
//E.g. in the diagram shown above, PosRelToLine("below", ...) would give kTRUE - the point is
//below the line.
//   ** from kaliveda:  KVIDLine::PosRelToLine [no modifications] **
Bool_t TCHIGraphCsI::PosRelToLine(Option_t * opt, Double_t px,
                              Double_t py, Double_t xp1,
                              Double_t yp1, Double_t xp2,
                              Double_t yp2)
{
 TVector2 P1(xp1, yp1), P2(xp2, yp2), P(px, py);
 TVector2 P1P2 = P2 - P1;
 TVector2 P1P = P - P1;
 TVector2 MP = P1P.Norm((P1P2));
 Double_t phi = MP.Phi() * TMath::RadToDeg();
   //Point is left of line if 90<phi<270
   //Point is right of line if 270<phi<90
   //Point is above line if 0<phi<180
   //Point is below line if 180<phi<360
 Bool_t result = kFALSE;
 if (!strcmp(opt, "left")) {
  result = (phi > 90. && phi < 270.);
 } else if (!strcmp(opt, "right")) {
  result = (phi < 90. || phi > 270.);
 } else if (!strcmp(opt, "above")) {
  result = (phi > 0. && phi < 180.);
 } else if (!strcmp(opt, "below")) {
  result = (phi > 180. && phi < 360.);
 }
 return result;
}




//The relative position of point (px,py) with respect to the line is tested.
//The option string can be "left", "right", "above" or "below".
//WhereAmI( px, py, "above") returns kTRUE if the point is above the line etc. etc.
//
//First of all, the closest segment/point of the line to the point is found.
//Then the relative position of the point and this segment/point is tested.
// ** Adapted from kaliveda:  KVIDLine::WhereAmI **
Bool_t TCHIGraphCsI::WhereAmI(TGraph *gr, Double_t px, Double_t py, Option_t *opt)
{
//Find closest point/segment
 Int_t i_close = 0;
 Double_t d = DistanceToLine(gr, px, py, i_close);
 Double_t *fx = gr->GetX();
 Double_t *fy = gr->GetY();
 
 if (d < 0) {
 //negative distance => closest point
  if (!strcmp(opt, "left")) {
   return (px < fx[i_close]);
   } else if (!strcmp(opt, "right")) {
   return (px > fx[i_close]);
   } else if (!strcmp(opt, "above")) {
   return (py > fy[i_close]);
   } else if (!strcmp(opt, "below")) {
   return (py < fy[i_close]);
   }
   return kFALSE;
  } 
 else {
  //positive distance => closest segment
  return PosRelToLine(opt, px, py, fx[i_close], fy[i_close],
                          fx[i_close + 1], fy[i_close + 1]);
 }
 return kFALSE;
}

//Returns kTRUE for point (x,y) if :
// axis = "" (default) and both x and y lie inside the endpoints of the line x1 < x < x2, y1 < y < y2
// axis = "x" and x lies inside the endpoints x1 < x < x2
// axis = "y" and y lies inside the endpoints y1 < y < y2
// * Adapted from kaliveda:  KVIDLine::IsBetweenEndPoints **
Bool_t TCHIGraphCsI::IsBetweenEndPoints(TGraph *gr, Double_t x, Double_t y, const Char_t *axis) 
{
 //Note: here we do not transform to uppercase as in the original routine
 //      but we accept also lowercase string. 
 string ax(axis);
 Double_t x1, y1, x2, y2;
 GetStartPoint(gr, x1, y1);
 GetEndPoint(gr, x2, y2);
 Double_t xmin = TMath::Min(x1, x2);
 Double_t xmax = TMath::Max(x1, x2);
 Double_t ymin = TMath::Min(y1, y2);
 Double_t ymax = TMath::Max(y1, y2);
 Bool_t in_range_x = (x <= xmax && x >= xmin);
 Bool_t in_range_y = (y <= ymax && y >= ymin);
 if (ax == "X" || ax=="x") {
  return in_range_x;
 } else if (ax == "Y" || ax=="y") {
  return in_range_y;
 }

 return (in_range_x && in_range_y);
}

// Find the coordinates of first point in line
// Adapted from kaliveda  KVIDLine::GetStartPoint
void TCHIGraphCsI::GetStartPoint(TGraph *gr, Double_t &x, Double_t &y)
{
 gr->GetPoint(0,x,y);
}

// Find the coordinates of last point in line
// Adapted from kaliveda  KVIDLine::GetEndPoint
void TCHIGraphCsI::GetEndPoint(TGraph *gr, Double_t &x, Double_t &y)
{
 Int_t n = gr->GetN() - 1; //last point
 gr->GetPoint(n,x,y);
}

//Performs rejection of gammas - returns kTRUE if point above gamma line
//(note - position w.r.t. IMF line is not considered)
//Returns kFALSE for points below gamma line i.e. for gammas ;-)
//If no gamma line defined, returns kTRUE for all points
Bool_t TCHIGraphCsI::IsIdentifiable(Double_t x, Double_t y)
{
 if(fgamma->GetN()==0)return kTRUE;
 if(WhereAmI(fgamma, x, y, "above")==kTRUE) {
  return kTRUE;
 }
 else {
  return kFALSE;
 }
}

// This method will locate (at most) four lines close to the point (x,y), the point must
// lie within the endpoints (in X) of each line (the lines "embrace" the point).
// Returns kTRUE if at least one line is found. Identification can then be carried out with
// either IdentZA or IdentZ (see Identify).
// Returns kFALSE if no lines are found (not even a closest embracing line).
//
// We look for two lines above the point and two lines below the point, as in one of 
// the following two cases:
//
//                      ------------------------     ksups       ---------------------
//
//  closest --->  ------------------------      ksup        ---------------------
//                                     X
//
//                                                                                             X
//                      ------------------------      kinf        ---------------------     <--- closest
//
//                      ------------------------     kinfi        ---------------------
//
// First we find the closest embracing line to the point, using FindNearestEmbracingIDLine.
// Then we search above and below for the other 'embracing' lines. Note that no condition is
// applied regarding the distances to these lines: the lines must have been sorted in order of increasing
// ordinate before hand in Initialize(), we simply use the order of lines in the list of identifiers.
// The Z, A, width and distance to each of these lines are stored in the variables
// Zsups, Asups, wsups, dsups
// etc. etc. to be used by IdentZA or IdentZ.
// Kaliveda * adapted from KVIDZAGrid class method * 
Bool_t TCHIGraphCsI::FindFourEmbracingLines(Double_t x, Double_t y, const Char_t* position)
{
 finfi = finf = fsup = fsups = -1;
 fdinf = fdsup = fdinfi = fdsups = 0.;
 fwinf = fwsup = fwinfi = fwsups = 16000.;
 fZinfi = fZinf = fZsup = fZsups = fAinfi = fAinf = fAsup = fAsups = -1;
 fClosest = fLsups = fLsup = fLinf = fLinfi = 0;
 fDistanceClosest = -1.;
 fIdxClosest = -1;


 fClosest = FindNearestEmbracingIDLine(x,y,position,"x", fIdxClosest, finf, fsup, fDistanceClosest, 
            fdinf, fdsup);

 if (!fClosest) return kFALSE;   //no lines found

 Int_t dummy=0;
 if (finf>-1 && finf==fIdxClosest) {
  //point is above closest line, closest line is "kinf"
  //need to look for 2 lines above (ksup, ksups) and 1 line below (kinfi)
  fLinf = fClosest;
//  fwinf = ((KVIDZALine*)fLinf)->GetWidth();   /* skipped: not implemented */ 
   fZinf = GetZ(fLinf);
   if(fZinf>10)cout << "dentro FindFourEmbracingLines fZinf=" << fZinf<< endl;
   fAinf = GetA(fLinf);
  
  if (fsup>-1) {
   fLsup = (TGraph *)fgrid->At(fsup);
// fLsup = (KVIDLine*)GetIdentifierAt(ksup);    /* original code */
//   fwsup = ((KVIDZALine*)fLsup)->GetWidth();  /* skipped: not implemented */
   fZsup = GetZ(fLsup);
   fAsup = GetA(fLsup);
  }
 }
 else if (fsup>-1 && fsup==fIdxClosest) {
 //point is below closest line, closest line is "ksup"
 //need to look for 1 line above (ksups) and 2 lines below (kinf, kinfi)
  fLsup = fClosest;        
//  fwsup = ((KVIDZALine*)fLsup)->GetWidth();  /* skipped: not implemented */ 
  fZsup = GetZ(fLsup);
  fAsup = GetA(fLsup);
  if (finf>-1) {
   fLinf = (TGraph *)fgrid->At(finf);
   //fLinf = (KVIDLine*)GetIdentifierAt(kinf);  /* original code */
   //fwinf = ((KVIDZALine*)fLinf)->GetWidth();  /* skipped: not implemented */ 
   fZinf = GetZ(fLinf);
   if(fZinf>10)cout << "dentro FindFourEmbracingLines 2nd fZinf=" << fZinf<< endl;
   fAinf = GetA(fLinf);                   
  }
 }
 else {
  cout<<"TCHIGraphCsI>> FindFourEmbracingLines: "<<
          "I do not understand the result of FindNearestEmbracingIDLine!!!"<<endl;
  return kFALSE;
 }

 if ( finf > -1 ) {
  // look for finfi line -> next line below 'inf' line
  finfi = finf;
  fLinfi = FindNextEmbracingLine(finfi, -1, x, y, "x");
  if ( !fLinfi ) finfi = -1;  //no 'infi' line found
  else {
   fdinfi = TMath::Abs(DistanceToLine(fLinfi,x,y,dummy) );
   // winfi = ((KVIDZALine*)fLinfi)->GetWidth();  /* skipped: not well understood what is fwidth */ 
   fZinfi = GetZ(fLinfi);
   fAinfi = GetA(fLinfi);         
  }
 }

 if ( fsup > -1 ) {
 // look for ksups line -> next line above 'sup' line
  fsups = fsup;
  fLsups = FindNextEmbracingLine(fsups, 1, x, y, "x");
  if ( !fLsups )fsups = -1; // no 'sups' line found
  else {
   fdsups = TMath::Abs(DistanceToLine(fLsups,x,y,dummy) );
   //wsups = ((KVIDZALine*)fLsups)->GetWidth();  /* skipped: not well understood what is fwidth */ 
   fZsups = GetZ(fLsups);
   fAsups = GetA(fLsups);         
  }
 }
 return kTRUE;
}

//
// * KALIVEDA adapted from KVIDLine::FindNearestEmbracingIDLine *
//   slight modifications
//
TGraph* TCHIGraphCsI::FindNearestEmbracingIDLine(Double_t x, Double_t y, const Char_t *position, 
                      const Char_t* axis, Int_t &idx, Int_t & idx_min, Int_t & idx_max, 
                      Double_t &dist, Double_t &dist_min, Double_t &dist_max)
{
 static TList emL;

 Int_t nlines = GetIDLinesEmbracingPoint(axis, x, y, emL);
 if ( !nlines ) return 0;     // no lines

 idx_min = 0;                 //minimum index
 idx_max = nlines - 1;        // maximum index
 idx = idx_max / 2;           //current index i.e. we begin in the middle
 dist = dist_min = dist_max = -1.;

 while (idx_max > idx_min + 1) {
  TGraph *line = (TGraph *)emL.At(idx);
  Bool_t point_above_line = WhereAmI(line, x, y, position);
/*  original code
  KVIDLine *line = (KVIDLine*)emL.At(idx);
  Bool_t point_above_line = line->WhereAmI(x, y, position);
*/
  if (point_above_line) {
  //increase index
   idx_min = idx;
   idx += (Int_t) ((idx_max - idx) / 2 + 0.5);
  }
  else {
  //decrease index
   idx_max = idx;
   idx -= (Int_t) ((idx - idx_min) / 2 + 0.5);
  }
 }

 //calculate distance of point to the two lines above and below
 TGraph *upper = (TGraph *)emL.At(idx_max);
 TGraph *lower = (TGraph *)emL.At(idx_min);
 Int_t dummy = 0;
 //if idx_max = nlines-1, the point may be above the last line
 //in which case we put idx_max = -1 (no line above point)
 if(idx_max == nlines-1) {
  if(WhereAmI(upper, x, y, position)) {
  // above last line
/*
  original code: get the index of "upper" in list  
  idx = idx_min = fIdentifiers->IndexOf(upper);   //index of last line
  critical point: it is supposed here and in the following that in the original 
  code we are looking at line index in the main grid list "fgrid" for a given telescope
*/
   idx = idx_min = fgrid->IndexOf(upper);
   idx_max = -1;
   dist = dist_min = TMath::Abs(DistanceToLine(upper, x, y, dummy));
   return upper;
  }
 }


//if idx_min = 0, the point may be below the first line
//in which case we put idx_min = -1 (no line below point)
 if(idx_min == 0) {
  if(!WhereAmI(lower, x, y, position)) {
   // below first line
   idx_min = -1;
   idx = idx_max = fgrid->IndexOf(lower); // index of first line
   dist = dist_max = TMath::Abs(DistanceToLine(lower, x, y, dummy));
   return lower;
  }
 }

 dist_max = TMath::Abs(DistanceToLine(upper, x, y, dummy));
 dist_min = TMath::Abs(DistanceToLine(lower, x, y, dummy));
 // convert indices back to index in main list
 idx_min = fgrid->IndexOf(lower);
 idx_max = fgrid->IndexOf(upper);
 if ( dist_max < dist_min ) {
  dist = dist_max;
  idx = idx_max;
  return upper;
 }
 dist = dist_min;
 idx = idx_min;
 return lower;
}

// Replaces contents of TList 'tmp' with subset of ID lines 
// for which IsBetweenEndPoints(x,y,direction) == kTRUE.
// Kaliveda * adapted from KVIDGrid class method [warning: heavily modified] * 
Int_t TCHIGraphCsI::GetIDLinesEmbracingPoint(const Char_t *direction,
        Double_t x, Double_t y, TList& tmp)
{
 //nlines = new number of lines in list

 TIter next(fgrid);
 Int_t nlines = 0;
 TGraph *line;
 tmp.Clear();
 while ( (line = (TGraph *)next()) ) {
  if (IsBetweenEndPoints(line, x, y, direction)) {
   tmp.Add(line);
   nlines++;
  }  
 }
//the original loop: KVIDGrid::GetIDLinesEmbracingPoint
/*  
    TIter next(GetIdentifiers());
    Int_t nlines = 0;
    KVIDLine *line;
    tmp.Clear();
    while ((line = (KVIDLine *) next()))
    {
        if (line->IsBetweenEndPoints(x, y, direction))
        {
            tmp.Add(line);
            nlines++;
        }
    }
*/
 return nlines;
}


// Starting from the line with given 'index', we search for the next line in the list of identifiers 
// for which KVIDLine::IsBetweenEndPoints(x,y,axis) returns kTRUE.
// 'inc_index' is the step used to scan the list, i.e. inc_index=1 will scan index+1, index+2, ...
// inc_index=-1 will scan index-1, index-2, ...
// Returns pointer to line (or 0x0 if not found) and 'index' contains index of this line 
// (or -1 if no line found)
// Kaliveda * adapted from KVIDGrid class method [slightly modified] * 
TGraph* TCHIGraphCsI::FindNextEmbracingLine(Int_t &index, Int_t inc_index, Double_t x, Double_t y, 
                              const Char_t* axis)
{
 int i = index + inc_index;
 Int_t nlines = fgrid->GetSize();
 TGraph *l = 0;
 while ( (i > -1 && i < nlines) ) {
  l = (TGraph *)fgrid->At(i);
  if (IsBetweenEndPoints(l, x,y,axis) ) break;
  i+=inc_index;
 }
 if (i<0 || i>=nlines) {
 // no line found
  index=-1;
  return 0;
 }
 index = i;
 return l;
/*   original code KVIDGrid::FindNextEmbracingLine
     register int i=index+inc_index;
     Int_t nlines = GetNumberOfIdentifiers();
     KVIDLine* l=0;
     while ( (i > -1 && i < nlines) )
     {
       l = (KVIDLine*)GetIdentifierAt(i);
       if ( l->IsBetweenEndPoints(x,y,axis) ) break;
       i+=inc_index;
     }
*/
}

// Find Z, A and 'real A' for point (x,y) once closest lines to point have been found.
// Double_t A = mass calculated by interpolation
// Adapted from  KVIDGCsI::IdentZA. This large routine must be absolutely unchanged and
// is adapted "as the original" one just changing the variables names to follow 
// TCHIEvent data member conventions. It is 
// adapted by Kaliveda from an original Tassan-Got code.  
// **  Note in particular:
//                TGraph *nextline = (TGraph *)GetIdentifierAt(idx);
//                if (GetZ(nextline) == Z
//                        && !IsBetweenEndPoints(nextline, x, y, "x"))
//                {
//                    fICode = ECode(fICode+1);    //dirty casting but ok
//  respect to the original code 
//                KVIDCsIRLLine *nextline = (KVIDCsIRLLine *) GetIdentifierAt(idx);
//                if (nextline->GetZ() == Z
//                        && !nextline->IsBetweenEndPoints(x, y, "x"))
//                {
//                    fICode++;
////////////////////////////////////////////////////////////////////////////////////////
void TCHIGraphCsI::IdentZA(Double_t x, Double_t y, Int_t& Z, Double_t& A)
{
    fICode = kICODE0;
    A = -1.;
    Z = -1;
    fAint = 0;

//if(fIdxClosest==ksups) cout << "*** ";
//cout<<"ksups="<<ksups<<" Zsups = "<<Zsups<<"  Asups = "<<Asups<<"  wsups = "<<wsups<< "  dsups = " 
//    << dsups <<endl;
//if(fIdxClosest==ksup) cout << "*** ";
//cout << "ksup = " << ksup << " Zsup = " << Zsup << "  Asup = " << Asup << "  wsup = " 
//     << wsup << "  dsup = " << dsup << endl;
//if(fIdxClosest==kinf) cout << "*** ";
//   cout << "kinf = " << kinf << " Zinf = " << Zinf << "  Ainf = " << Ainf << "  winf = " 
//         << winf << "  dinf = " << dinf << endl;
//if(fIdxClosest==kinfi) cout << "*** ";
//cout << "kinfi = " << kinfi << " Zinfi = " << Zinfi << "  Ainfi = " << Ainfi << "  winfi = " 
//      << winfi << "  dinfi = " << dinfi << endl;

 Int_t ibif = 0;
 Int_t k = -1;
 Double_t yy, y1, y2;
 Int_t ix1, ix2;
 yy = y1 = y2 = 0;
 ix1 = ix2 = 0;
 if (fsup > -1)
 {
  if (finf > -1)
  {
//cout << " /******************* 2 lignes encadrant le point ont ete trouvees ***************/" << endl;
   Double_t dt = fdinf + fdsup;     //distance between the 2 lines
   if (fZinf == fZsup)
   {
    //   cout << "      /****************meme Z**************/" << endl;
    Z = fZinf;
    if(Z>10)cout << "linea 700 Z= " << Z<< endl; 
    Int_t dA = fAsup - fAinf;
    Double_t dist = dt / dA;    //distance between the 2 lines normalised to difference in A of lines
    /*** A = Asup ***/
    if (fdinf > fdsup)    //point is closest to upper line, 'sup' line
    {
     ibif = 1;
     k = fsup;
     yy = -fdsup;
     A = fAsup;
     fAint = fAsup;
     if (fsups > -1)          // there is a 'sups' line above the 2 which encadrent le point
     {
      y2 = fdsups - fdsup;
      if (fZsups == fZsup)
      {
       ibif = 0;
       y2 /= 2.;
       ix2 = fAsups - fAsup;
      }
      else
      {
       if (fZsups > 0)
        y2 /= 2.;       // 'sups' line is not IMF line
       Double_t x2 = fwsup;
       x2 = 0.5 * TMath::Max(x2, dist);
       y2 = TMath::Min(y2, x2);
       ix2 = 1;
      }
     }
     else           // ksups == -1 i.e. no 'sups' line
     {
      y2 = fwsup;
      y2 = 0.5 * TMath::Max(y2, dist);
      ix2 = 1;
     }
     y1 = -dt / 2.;
     ix1 = -dA;
    }
    /*** A = Ainf ***/
    else                //point is closest to lower line, 'inf' line
    {
     ibif = 2;
     k = finf;
     yy = fdinf;
     A = fAinf;
     fAint = fAinf;
     if (finfi > -1)          // there is a 'infi' line below the 2 which encadrent le point
     {
      y1 = 0.5 * (fdinfi - fdinf);
      if (fZinfi == fZinf)
      {
       ibif = 0;
       ix1 = fAinfi - fAinf;
       y1 = -y1;
      }
      else
      {
       Double_t x1 = fwinf;
       x1 = 0.5 * TMath::Max(x1, dist);
       y1 = -TMath::Min(y1, x1);
       ix1 = -1;
      }
     }
     else           // kinfi = -1 i.e. no 'infi' line
     {
      y1 = fwinf;
      y1 = -0.5 * TMath::Max(y1, dist);
      ix1 = -1;
     }
     y2 = dt / 2.;
     ix2 = dA;
    }
   }
   else
   {
    //cout << "         /****************Z differents**************/ " << endl;
    if (fZsup == -1)     //'sup' is the IMF line
    {
     dt *= 2.;
     fdsup = dt - fdinf;
    }
    /*** Z = Zsup ***/
    ibif = 3;
    if (fdinf > fdsup)    // closest to upper 'sup' line
    {
     k = fsup;
     yy = -fdsup;
     Z = fZsup;
     if(Z>10)cout << "linea 789 Z= " << Z<< endl; 
     A = fAsup;
     fAint = fAsup;
     y1 = 0.5 * fwsup;
     if (fsups > -1)          // there is a 'sups' line above the 2 which encadrent the point
     {
      y2 = fdsups - fdsup;
      if (fZsups == fZsup)
      {
       ibif = 2;
       ix2 = fAsups - fAsup;
       Double_t x1 = y2 / ix2 / 2.;
       y1 = TMath::Max(y1, x1);
       y1 = -TMath::Min(y1, dt / 2.);
       ix1 = -1;
       y2 /= 2.;
      }
      else
      {
       if (fZsups > 0)
        y2 /= 2.;       // 'sups" is not IMF line
       y2 = TMath::Min(y1, y2);
       ix2 = 1;
       y1 = -TMath::Min(y1, dt / 2.);
       ix1 = -1;
      }
     }
     else           // ksups == -1, i.e. no 'sups' line
     {
      fICode = kICODE7;   //a gauche de la ligne fragment, Z est alors un Zmin et le plus probable
      y2 = y1;
      ix2 = 1;
      y1 = -TMath::Min(y1, dt / 2.);
      ix1 = -1;
     }
    }
    /*** Z = Zinf ***/
    else                // closest to lower 'inf' line
    {
     k = finf;
     yy = fdinf;
     Z = fZinf;
     if(Z>10)cout << "linea 830 Z= " << Z<< endl; 
     A = fAinf;
     fAint = fAinf;
     y2 = 0.5 * fwinf;
     if (finfi > -1)        // there is a 'infi' line below the 2 which encadrent the point
     {
      y1 = fdinfi - fdinf;
      if (fZinfi == fZinf)
      {
       ibif = 1;
       ix1 = fAinfi - fAinf;
       Double_t x2 = -y1 / ix1 / 2.;
       y2 = TMath::Max(y2, x2);
       y2 = TMath::Min(y2, dt / 2.);
       ix2 = 1;
       y1 /= -2.;
      }
      else
      {
       y1 = -TMath::Min(y2, y1 / 2.);
       ix1 = -1;
       y2 = TMath::Min(y2, dt / 2.);
       ix2 = 1;
      }
     }
     else           // no kinfi line, kinfi==-1
     {
      y1 = -y2;
      ix1 = -1;
      y2 = TMath::Min(y2, dt / 2.);
      ix1 = 1;
     }
    }
   }
  }   // if(kinf>-1)...
  else if (fZsup > 0)        // 'sup' is not IMF line
  {
  //cout<<" /****************** Seule une ligne superieure a ete trouvee *********************/" << endl;
   ibif = 3;
   k = fsup;
   yy = -fdsup;
   Z = fZsup;
   if(Z>10)cout << "linea 873 Z= " << Z<< endl; 
   A = fAsup;
   fAint = fAsup;
   y1 = 0.5 * fwsup;
   if (fsups > -1)        // there is a 'sups' line above the closest line to the point
   {
    y2 = fdsups - fdsup;
    if (fZsups == fZsup)
    {
     ibif = 2;
     ix2 = fAsups - fAsup;
     Double_t x1 = y2 / ix2 / 2.;
     y1 = -TMath::Max(y1, x1);
     ix1 = -1;
     y2 /= 2.;
    }
    else
    {
     if (fZsups > 0)
      y2 /= 2.;     // 'sups' is not IMF line
     y2 = TMath::Min(y1, y2);
     ix2 = 1;
     y1 = -y1;
     ix1 = -1;
    }
   }
   else                 // no 'sups' line above closest line
   {
    fICode = kICODE7;   //a gauche de la ligne fragment, Z est alors un Zmin et le plus probable
    y2 = y1;
    ix2 = 1;
    y1 = -y1;
    ix1 = -1;
   }
  }
  else
  {
   fICode = kICODE8;      //  Z indetermine ou (x,y) hors limites
  }
 }  //if (fsup > -1)
 else if (finf > -1)
 {
  //cout <<"/****************** Seule une ligne inferieure a ete trouvee ***********************/" << endl;
  /*** Sep. fragment ***/
  if (fZinf == -1)           // 'inf' is IMF line
  {
   //point is above IMF line. Z = Z of last line in grid, A = -1
   k = -1;
   Z = GetZmax();
   if(Z>10)cout << "linea 922 Z= " << Z<< endl; 
   A = -1;
   fAint = 0;
   fICode = kICODE6;      // au-dessus de la ligne fragment, Z est alors un Zmin
  }
  /*** Ligne de crete (Z,A line)***/
  else
  {
   ibif = 3;
   k = finf;
   Z = fZinf;
   if(Z>10)cout << "linea 933 Z= " << Z<< endl; 
   A = fAinf;
   fAint = fAinf;
   yy = fdinf;
   y2 = 0.5 * fwinf;
   if (finfi > -1)
   {
    y1 = fdinfi - fdinf;
    if (fZinfi == fZinf)
    {
     ibif = 1;
     ix1 = fAinfi - fAinf;
     Double_t x2 = -y1 / ix1 / 2.;
     y2 = TMath::Max(y2, x2);
     ix2 = 1;
     y1 /= -2.;
    }
    else
    {
     y1 = -TMath::Min(y2, y1 / 2.);
     ix1 = -1;
     ix2 = 1;
    }
   }
   else
   {
    y1 = -y2;
    ix1 = -1;
    ix2 = 1;
   }
   fICode = kICODE7;      // a gauche de la ligne fragment, Z est alors un Zmin et le plus probable
  }
 }
 /*****************Aucune ligne n'a ete trouvee*********************************/
 else
 {
  fICode = kICODE8;         // Z indetermine ou (x,y) hors limites
 }
 /****************Test des bornes********************************************/
 if (k > -1 && fICode == kICODE0)
 {
  if (yy > y2)
   fICode = kICODE4;      // Z ok, masse hors limite superieure ou egale a A
 }
 if (k > -1 && (fICode == kICODE0 || fICode == kICODE7))
 {
  if (yy < y1)
   fICode = kICODE5;      // Z ok, masse hors limite inferieure ou egale a A
 }
 if (fICode == kICODE4 || fICode == kICODE5)
  { A = -1; fAint = 0; }
            
 /****************Interpolation de la masse: da = f*log(1+b*dy)********************/
 if (fICode == kICODE0 || (fICode == kICODE7 && yy <= y2))
 {
  Double_t deltaA = 0.;
  Bool_t i = kFALSE;
  Double_t dt, dist = y1 * y2;
  dt = 0.;
  if (ix2 == -ix1)          //dA1 = dA2
  {
   if(dist!=0) {
    dt = -(y1 + y2) / dist;
    i = kTRUE;
   }
            /*else
            	Warning("IdentZA","%s : cannot calculate interpolated mass (dist=%f), Areal will equal Aint (Z=%d Aint=%d fICode=%d)", GetName(), dist, Z, Aint, fICode);*/
  }
  else if (ix2 == -ix1 * 2)       // dA2 = 2*dA1
  {
   Double_t tmp = y1 * y1 - 4. * dist;
   if(tmp>0. && dist!=0) {
    dt = -(y1 + 2. * y2 -
    TMath::Sqrt(tmp)) / dist / 2.;
    i = kTRUE;
   }
     /*else
     Warning("IdentZA","%s : cannot calculate interpolated mass (y1*y1-4*dist=%f), Areal will equal Aint (Z=%d Aint=%d fICode=%d)", GetName(), tmp, Z, Aint, fICode);*/
  }
  else if (ix1 == -ix2 * 2)       // dA1 = 2*dA2
  {
   Double_t tmp = y2 * y2 - 4. * dist;
   if(tmp>0. && dist!=0) {
    dt = -(y2 + 2. * y1 +
    TMath::Sqrt(tmp)) / dist / 2.;
    i = kTRUE;
   }
   /*else
   Warning("IdentZA","%s : cannot calculate interpolated mass (y2*y2-4*dist=%f), Areal will equal Aint (Z=%d Aint=%d fICode=%d)", GetName(), tmp, Z, Aint, fICode);*/
  }
  if (i)
  {
   dist = dt * y2;
   if (TMath::Abs(dist) < 0.001)
   {
    if(y2!=0)
     deltaA = yy * ix2 / y2 / 2.;
     /*else
     Warning("IdentZA","%s : cannot calculate interpolated mass (y2=%f), Areal will equal Aint (Z=%d Aint=%d fICode=%d)", GetName(), y2, Z, Aint, fICode);*/
   }
   else
   {
    if(dist>-1. && dt*yy>-1.)
     deltaA = ix2 / 2. / TMath::Log(1. + dist) * TMath::Log(1. + dt * yy);
    /*else
    Warning("IdentZA","%s : cannot calculate interpolated mass (dist=%f dt*yy=%f), Areal will equal Aint (Z=%d Aint=%d fICode=%d)", GetName(), dist, dt*yy, Z, Aint, fICode);*/
   }
   A += deltaA;
  }
 }
 /***************D'autres masses sont-elles possibles ?*************************/
 if (fICode == kICODE0)
 {
  //cout << "icode = 0, ibif = " << ibif << endl;
  /***Masse superieure***/
  if (ibif == 1 || ibif == 3)
  {
   //We look at next line in the complete list of lines, after the closest line.
   //If it has the same Z as the closest line, but was excluded from research for closest line
   //because the point lies outside the endpoints, there remains a doubt about the mass:
   //on rajoute 1 a fICode, effectivement on le met = kICODE1
   Int_t idx = fIdxClosest;
// Warning: slight modified respect to original code [e.d.f.]
   if (idx > -1 && ++idx < GetNumberOfIdentifiers())
   {
    TGraph *nextline = (TGraph *)GetIdentifierAt(idx);
    if (GetZ(nextline) == Z
                       && !IsBetweenEndPoints(nextline, x, y, "x"))
    {
     fICode = ECode(fICode+1);   
     // Z ok, mais les masses superieures a A sont possibles
     //cout <<"//on rajoute 1 a fICode, effectivement on le met = kICODE1" << endl;
    }
   }
  }

  /***Masse inferieure***/
  if (ibif == 2 || ibif == 3)
  {
   //We look at line below the closest line in the complete list of lines.
   //If it has the same Z as the closest line, but was excluded from research for closest line
   //because the point lies outside the endpoints, there remains a doubt about the mass:
   //on rajoute 2 a fICode, so it can be = kICODE2 or kICODE3
   Int_t idx = fIdxClosest;
   if (idx > -1 && --idx >= 0)
   {
    TGraph *nextline = (TGraph *) GetIdentifierAt(idx);
    if (GetZ(nextline) == Z
                       && IsBetweenEndPoints(nextline, x, y, "x"))
    {
     fICode = ECode(fICode + 2);
     //cout << "//on rajoute 2 a fICode, so it can be = kICODE2 or kICODE3" << endl;
    }
   }
  }
// Warning: end modification
 }
 if(Z>10){
  cout << "*******************************"<< endl;
  cout << "dentro TCHIGraphCsI linea 1084  Z= "<< Z << endl;
  cout << "Z = " << Z << " A = " << A << " icode = " << fICode << endl;
 }
}