///////////////////////////////////////////////////////////
//  Chimera Event UnPacker
//  CsI Identification 
//  partially adapted upon KaliVeda (C) 2005 by J.D. Frankland
//  identification methods for CsI(Tl)
//
//  E. De Filippo  v 1.1
//  TCHILNS / Exochim - AsyEos collaboration 
//  TCHIGridCsI implementation 
///////////////////////////////////////////////////////////
#include "TCHIGridCsI.h"

// Constructor
TCHIGridCsI::TCHIGridCsI(string filename) : 
 fgridfilename(filename),
 fgridOK(false)
{
 for(int i=0; i<MAXTEL; i++) {
  fgraph[i] = 0;
  frmin[i] = frmax[i] = 0;
 }

 // A simple light particle's map based on key Z*100 + A 
 felem[101] = "1H";
 felem[102] = "2H";
 felem[103] = "3H";
 felem[203] = "3He";
 felem[204] = "4He";
 felem[206] = "6He";
 felem[306] = "6Li";
 felem[307] = "7Li";
 felem[308] = "8Li";
 felem[309] = "9Li";
 felem[407] = "7Be";
 felem[408] = "8Be";
 felem[208] = "8Be";   //Z2A8 used for correct 8Be sorting (8Be is below Li and over 6He)
 felem[409] = "9Be"; 
 felem[410] = "10Be"; 

 InitChimera();
}

//Destructor
TCHIGridCsI::~TCHIGridCsI()
{
 for(int i=0; i<MAXTEL; i++) {
  if(fgraph[i] != 0) {
   delete fgraph[i];
   fgraph[i] = 0;
  }
 }
}  

// calculate ring-module value for each telescope number in Chimera
// The ring value is here the same for internal and external rings
// as required by the isospin program for CsI identification
void TCHIGridCsI::InitChimera()
{
 Int_t modules[26] = {32,48,64,80,80,96,96,96,96,
                      32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,16,8};
 
 Int_t ntel=0;
 for(Int_t i=0; i<26; i++) {
  Int_t k=modules[i]/2;
  for(Int_t j=0; j<modules[i]; j++) {
   fring[ntel] = i+1;
   if(i<9 && j>=k) {
    fmod[ntel] = j-k;
   }
   else {
    fmod[ntel] = j;
   }  
   ntel++;
  }
 }
}

//read grids file
bool TCHIGridCsI::ReadAsciiFile()
{
 int row = 0;
 int ntel = 0;
 int ncut = 0;
 int ngri = 0;

 fgridOK = false;
 ifstream gridf(fgridfilename.c_str()); 
 if (!gridf.good()) {
  cout<<"TCHIGridCsI>> "<< "File "<<fgridfilename<<
        " cannot be opened"<<endl; 
  return fgridOK;
 }
 
 string s;
 int telnum = -1;
 int rmin=0;
 int rmax=0;
 while (!gridf.eof())
 {
  //read a line 
  getline(gridf, s);
  row++;
  if(s.find("++")==0) {   //find '++' at beginning 
   s.erase(0,2);
   if(s != "KVIDGCsI") {
    cout<<"TCHIGridCsI>> The ascii file for Csi Chimera Identification should\n"<<
          "follow the KaliVeda grid format standard: missing KVIDGCsi label"
          <<endl;
    return fgridOK;
   }
  }

  // run number pseudo-command (out of kaliveda standard)
  if(s.find("#@")==0) { 
   s.erase(0,2);
   istringstream os(s);
   os>>rmin>>rmax;
  }  

  // End of telescope grid, next telescope starts 
  if(s.find("!")==0) {
   ntel++;
   telnum = -1;
   continue; 
  }

  //read the telescope number and create the relative object
  if(s.find("<PARAMETER> IDTelescopes=CSI_R_L_")==0) {
   s.erase(0, 33);
   istringstream os(s);
   os>>telnum;
   if(telnum>=0 && telnum<= MAXTEL) {
    if(rmin>0 && rmax>0) {
     frmin[telnum] = rmin;
     frmax[telnum] = rmax;
     rmin = 0;
     rmax = 0;
    }
    if(fgraph[telnum] != 0) {
     cout<<"TCHIGridCsI>> Warning telescope "<<telnum<<
     " already existing...replaced"<<endl; 
     delete fgraph[telnum];
     fgraph[telnum] = new TCHIGraphCsI(telnum);
    }
    else {
     fgraph[telnum] = new TCHIGraphCsI(telnum);
    }
   }
  }

  //Get name of class by stripping off the '+' at the start of the line
  if(s.find("+")==0) {  
   s.erase(0, 1);
   if(s == "KVIDCsIRLLine") {  
   //next line is type ("ID" or "OK") followed by ':', followed by name of line
    getline(gridf,s);
    row++;    
    if(s.find("ID:")==0) {
     ngri++;
     AddLine(gridf, &row, telnum);  
    }
    else {
     cout<<"TCHIGridCsI>> Severe error ID not found at row "<<row<<
       " grid line skipped"<<endl;
    }
   }

   if(s == "KVIDCutLine") {  
   //next line is type ("ID" or "OK") followed by ':', followed by name of line
    getline(gridf,s);
    row++;  
    ncut++;  
    if(s.find("OK:gamma_line")==0) {
     AddCut(gridf, &row, telnum, kgamma);  
    }
    else if(s.find("OK:IMF_line")==0) {
     AddCut(gridf, &row, telnum, kIMF);  
    }
    else {
     cout<<"TCHIGridCsI>> Severe error OK:cut not found at row "<<row<<
       " cut line skipped"<<endl;
    }
   }
  }
 }
 fgridOK = true;
 cout<<"TCHIGridCsI>> Read "<<row-1<<" lines from file "<<fgridfilename<<endl;
 cout<<"TCHIGridCsI>> Processed "<<ngri<<" grids for "<<ntel<<" telescopes"<<endl;
 cout<<"TCHIGridCsI>> Stored "<<ncut<<" delimiters lines"<<endl;

 gridf.close();
 return fgridOK; 
}

// Add grids lines
int TCHIGridCsI::AddLine(ifstream &gridf, int *row, int tel)
{
 string s;
 stringstream name;
 int Z, A, npar;
 float x,y;
 

 if(tel==-1)return -1; 
 getline(gridf,s);
 (*row)++;
 istringstream os(s);
 os>>Z>>A;
 // special case for 8Be. We write Z=2 A=8 in order to 
 // maintain 8Be between alpha and Z=3 during "sorting"
 // GetZ method know this and will give the 
 // correct Z=4 for this line. 
 if(Z==4 && A==8) {
  Z=2;
 }
 name<<"Z"<<Z<<"A"<<A<<"_"<<"T"<<tel<<"_"<<felem[Z*100+A];
 getline(gridf,s);
 (*row)++;
 istringstream osn(s);
 osn>>npar;
 TGraph *pg = new TGraph(npar);   
 pg->SetName(name.str().c_str());    //set name to grid
 pg->SetTitle(name.str().c_str());   //set title to grid 
 for(int i=0; i<npar; i++) {
  getline(gridf,s);
  (*row)++;
  istringstream ns(s);
  ns>>x>>y;
  pg->SetPoint(i, x, y);
 }
 if(fgraph[tel] != 0) {
  fgraph[tel]->GetGrid()->Add(pg);   // Add the grid to the list
 }
 fgraph[tel]->SetValidGrid(true);
 return npar;
}

// Add cuts lines
int  TCHIGridCsI::AddCut(ifstream &gridf, int *row, int tel, ECut cuttype) 
{
 string s;
 stringstream name;
 int Z, A, npar;
 float x,y;
 TGraph *pg;
 TGraph *gamma = fgraph[tel]->GetGammaL();
 TGraph *imf = fgraph[tel]->GetIMFL();

 if(tel==-1)return -1; 
 getline(gridf,s);
 (*row)++;
 if(cuttype == kgamma) {
  name<<"gamma_line_"<<s<<"_"<<"T"<<tel;
 }
 if(cuttype == kIMF) {
  name<<"IMF_line_"<<s<<"_"<<"T"<<tel;
 }
 const char *ntit = name.str().c_str();
 getline(gridf,s);
 (*row)++;
 istringstream osn(s);
 osn>>npar;

 if(cuttype == kgamma) {
  pg = gamma;
  if(pg->GetN() != 0) {
   cout<<"TCHIGridCsI>> Warning gamma line already existing for tel "<<tel<<
  "...replaced"<<endl;
   delete pg;
   pg = new TGraph;
   gamma = pg;
  }
 }

 if(cuttype == kIMF) {
  pg = imf;
  if(pg->GetN() != 0) {
   cout<<"TCHIGridCsI>> Warning IMF line already existing for tel "<<tel<<
  "...replaced"<<endl;
   delete pg;
   pg = new TGraph;
   imf = pg;
  }
 }

 pg->SetName(ntit);    //set name to grid
 pg->SetTitle(ntit);   //set title to grid 
 for(int i=0; i<npar; i++) {
  getline(gridf,s);
  (*row)++;
  istringstream ns(s);
  ns>>x>>y;
  pg->SetPoint(i, x, y);
 }
 fgraph[tel]->SetValidCuts(true);
 return npar;
}

// Print grid contents on stdout
void TCHIGridCsI::PrintLine(int ntel, int opt)
{
 TList *lg;

 //Security pointers check  
 if(fgraph[ntel] == 0 ) {
  cout<<"TCHIGridCsI>> Grid line not existing for telescope "<<ntel<<endl;
  return;
 }
 // now iteration over grids 
 if(fgraph[ntel]->HasValidGrid()) {
  lg = fgraph[ntel]->GetGrid();
  TIter next(lg);
  TGraph *gr = 0;
  int i=0;
  while ( (gr = (TGraph *)next()) ) {
   int np = gr->GetN();
   cout<<"Name:  "<<gr->GetName()<<"  Points "<<np<<
   " Z="<<fgraph[ntel]->GetZ(gr)<<" A="<<fgraph[ntel]->GetA(gr)<<endl;
   if(opt==1) 
    gr->Print();
  } 
 }
}

// Export all valid grids in an ascii file 
// compatible with the Chimera-Isospin program 
// fortran CsI routines and Chimera graphical 
// programs for CsI. Utility method.  
void TCHIGridCsI::Export2Isospin(string chname)
{
 TList *lg;

 ofstream chiout(chname.c_str(), ios::out); 
 if (!chiout.good()) {
  cout<<"TCHIGridCsI>> "<< "File "<<chname<<
        " cannot be opened for output"<<endl; 
  return;
 }
 
 for(Int_t ntel=0; ntel<MAXTEL; ntel++) { 
  //Security pointers check  
  if(fgraph[ntel] == 0 ) {
   continue;
  }

  TCHIGraphCsI *currg = fgraph[ntel];
  // now iteration over grids 
  if(currg->HasValidGrid()) {
   chiout<<setw(5)<<fring[ntel]<<setw(5)<<ntel<<
   setw(6)<<frmin[ntel]<<setw(6)<<frmax[ntel]<<endl;
   lg = currg->GetGrid();
   int nlines = lg->GetEntries();
   if(lg>0) { 
    if(currg->IsIMFlineadded())nlines--;
    chiout<<setw(5)<<nlines+2<<endl;    
   }
   // Add gamma
   TGraph *gamma = currg->GetGammaL();
   int ng = gamma->GetN();
   chiout<<setw(5)<<"0"<<setw(5)<<"0"<<setw(5)<<ng<<endl;
   // Add Isotopes 
   TIter next(lg);
   TGraph *gr = 0;
   int k=nlines;
   while ( (gr = (TGraph *)next()) && k--) {
    int np = gr->GetN();
    chiout<<setw(5)<<currg->GetZ(gr)<<setw(5)<<currg->GetA(gr)<<setw(5)<<np<<endl;
   }
   // Add IMF 
   TGraph *IMFLine = currg->GetIMFL(); 
   int nm = IMFLine->GetN();
   chiout<<setw(5)<<"-1"<<setw(5)<<"-1"<<setw(5)<<nm<<endl;

   // Chimera-isospin grids ends with a complete set of floating points
   // in rows of five couples (x,y) following the order 
   // 1) gamma opening/ 2) particles / 3) IMF closing

   //gamma line
   Double_t *fx = gamma->GetX();
   Double_t *fy = gamma->GetY(); 
   int j=5;
   chiout.setf(ios::fixed);
   chiout<<setprecision(1); 
   for(Int_t i=0; i<ng; i++) {
    j--;
    chiout<<right<<setw(8)<<fx[i]<<setw(8)<<fy[i];
    if(j==0) { 
     chiout<<endl;
     j=5;
    }
   }

   // particles grids lines 
   next=lg;
   k=nlines;
   while ( (gr = (TGraph *)next()) && k--) {
    int np = gr->GetN();
    fx = gr->GetX();
    fy = gr->GetY();
    for(Int_t i=0; i<np; i++) {
    j--;
    chiout<<right<<setw(8)<<fx[i]<<setw(8)<<fy[i];
     if(j==0) { 
      chiout<<endl;
      j=5;
     }
    }    
   } 

   //IMF line
   fx = IMFLine->GetX();
   fy = IMFLine->GetY(); 
   for(Int_t i=0; i<nm; i++) {
    j--;
    chiout<<right<<setw(8)<<fx[i]<<setw(8)<<fy[i];
    if(j==0) { 
     chiout<<endl;
     j=5;
    }
   }
   // closing telescope
   if(j!=5)chiout<<endl;
  }
 }
 chiout.close();
}

// This method initializes the identifications grids for all telescopes. 
// This procedure is drastically different from the Kaliveda one. It is 
// done globally (before identification for all telescopes for which 
// a valid grid is found. As in Kaliveda it makes: 1) sorting, 2) not implemented : 
// calculate line width    
// 3) find Zmax of sorted lines. Several KaliVeda "Init" methods are merged here.   
// It has to be mandatory called once before using the grids. This means that, 
// following the guide-lines of the Chimera-Isospin program, all the grids have to be ready 
// and initialized for all telescopes before identification.  
void TCHIGridCsI::Initialize()
{
 for(int ntel=0; ntel<MAXTEL; ntel++) {
  if(fgraph[ntel] == 0 )continue;
  if(fgraph[ntel]->HasValidGrid()) {
   fgraph[ntel]->Initialize();
  } 
 }
}


// Set Z and A of nucleus based on position in R-L grid
// The identification of gammas (kICODE10) and charged particles is performed
// Main Identification method [heavily modified respect to kaliveda but it is 
// supposed doing the same things]
// * adapted from  KVIDGCsI::Identify
void TCHIGridCsI::Identify(int ntel, Double_t x, Double_t y, TCHIResult *idr)
{
 TCHIGraphCsI *currg = fgraph[ntel];  // the current lines+cuts object for tel ntel

 idr->Init();
 if(!currg) {            //no pointer, no grids
  idr->fcod = kICODE_nogrid;
  return;
 }
 if(!currg->HasValidGrid()) {
  idr->fcod = kICODE9;  //incomplete grids
  return;
 }

 TList *Identifiers = currg->GetGrid();   //the current list of identifiers (grids)
 TGraph *IMFLine = currg->GetIMFL();      //the IMF closing line 
// before first use of this method, we add the IMF line to the list of identifiers
// which is necessary for IdentZA to work correctly. IMF line is cloned in order 
// to avoid problems with class destructors
 if(!currg->IsIMFlineadded()) {
  TGraph *newIMF = (TGraph *)IMFLine->Clone();
  if(IMFLine)Identifiers->AddLast(newIMF);
  currg->SetIMFLineAdded(kTRUE);
 }

 if(!currg->IsIdentifiable(x, y)) {
  //point below gamma line
  idr->fcod = kICODE10;
  idr->fZ = 0;
  idr->fA = 0;
  idr->fIDOK = kTRUE;
  return;
 }
 
 if (!currg->FindFourEmbracingLines(x,y,"above")) {
  //no lines corresponding to point were found
  idr->fcod = kICODE8;        // Z indetermine ou (x,y) hors limites
  return;
 }

 Int_t Z;
 Double_t A;
 currg->IdentZA(x, y, Z, A);    // Main identification Tassan-Got-like routine 
 if(Z>10) cout << "dentro TCHIGridCsI linea 496 ntel= " << ntel << " Z= "<< Z << endl;
 idr->fZ = Z;
 idr->fA = currg->GetAint();
 idr->fPID = A;
 idr->fcod = currg->GetICode();
 if(idr->fcod < kICODE4){
  idr->fIDOK = kTRUE;
  idr->fZident = kTRUE;
  idr->fAident = kTRUE;
 }
}