/// merge N events into 1.
 /// to run:  convert NMergeEvents NGeneratedEvents. If you want merge 20 events into 10, then print: convert 2 10

// #define DEBUG

 
#include <iostream>
#include <fstream>
#include <vector>
#include <map>

#include <cstring>
#include <stdio.h> // for atoi
#include <stdlib.h>

using namespace std;


typedef  float fscal;

// #include "code/Algo/L1Algo/L1Types.h"
#include "code/Performance/CbmL1MC.h"
#include "code/Algo/L1Algo/L1StsHit.h"
#include "code/Algo/L1Algo/L1Strip.h"

const int MCTrStep = 10000;
char work_dir[100] = "/u/ikulakov/STAP/data/mbias_1000ev_25022010/";  // path to input and output files
// char work_dir[100] = "/u/ikulakov/STAP/data/centr_100ev_07042009/";  // path to input and output files
char out_dir[100] = "/u/ikulakov/STAP/data/mbias_1000ev_25022010_merged/";  // path to input and output files
const int NStations = 10;


class CbmL1HitStore{
  public:
    int ExtIndex;
    int iStation;
    double x, y;
//     bool operator>(CbmL1HitStore& a1) { return a1.y > y;}
    int IsLessThen(CbmL1HitStore& a1){ if (a1.y > y) return 1; if (a1.y == y) return 0; return -1; }
};

struct TData{

  vector< L1Strip > vStsStrips, vStsStripsB;
  vector< L1StsHit > vStsHits;
  vector< unsigned char > vSFlag;  // = iStation*4 + used*2 + used_by_duplets;
  vector< unsigned char > vSFlagB;
  int StsHitsStartIndex[20], StsHitsStopIndex[20];

  vector<CbmL1MCPoint> vMCPoints;  // list of MonteCarlo(=simulated) point
  vector<CbmL1MCTrack> vMCTracks; // list of MonteCarlo(=simulated) track
  vector<int> vHitMCRef;            // index the registred(effective(reference??)) MonteCarlo points(=hits) in vMCPoints
  vector<CbmL1HitStore> vHitStore;  // array of all hits. all hits must have corespondend point, which marked in vHitMCRef

  void Add(TData& data_new, int iEv);
  int Correct();
//   TData& operator= (TData& a){
//     vStsStrips = a.vStsStrips;
//     vStsStripsB = a.vStsStripsB;
//     vStsHits = a.vStsHits;
//     vSFlag = a.vSFlag;
//     vSFlagB = a.vSFlagB;
// 
//     for (int i = 0; i < 20; i++){
//       StsHitsStartIndex[i] = a.StsHitsStartIndex[i];
//       StsHitsStopIndex[i] = a.StsHitsStopIndex[i];
//     };
//     
//     vMCPoints = a.vMCPoints;
//     vMCTracks = a.vMCTracks;
//     vHitMCRef = a.vHitMCRef;
//     vHitStore = a.vHitStore;
//   };
};

int TData::Correct()
{
  int flag = 1;
  int curStart = 0, curStop = 0;
  for (int i = 0; i < NStations; i++){
#ifdef DEBUG
    cout << i << " : " << " " << StsHitsStartIndex[i] << endl;
    cout << i << " : " << " " << StsHitsStopIndex[i] << endl;
#endif
    if (curStart > StsHitsStartIndex[i]){
#ifdef DEBUG
      cout << " : StsHitsStartIndex " << StsHitsStartIndex[i] << " -> ";
      StsHitsStartIndex[i] = curStart;
      cout << StsHitsStartIndex[i] << endl;
#else
      StsHitsStartIndex[i] = curStart;
#endif
      flag = 0;
    }
    else curStart = StsHitsStartIndex[i];
    if (curStart > StsHitsStopIndex[i]){
#ifdef DEBUG
      cout << " : StsHitsStopIndex " << StsHitsStopIndex[i] << " -> ";
      StsHitsStopIndex[i] = curStart-1;
      cout << StsHitsStopIndex[i] << endl;
#else
      StsHitsStopIndex[i] = curStart-1;
#endif
      flag = 0;
    }
    else curStart = StsHitsStopIndex[i]+1;
  };
  return flag;
};

void TData::Add(TData& data_new, int iEv)
{
  const int stSize =  vStsStrips.size();
  const int stBSize =  vStsStripsB.size();
  const int mcPointsSize =  vMCPoints.size();

  cout << data_new.vStsStrips.size()+data_new.vStsStripsB.size() << " strips and ";
  cout << data_new.vStsHits.size() << " hits for ";
  cout << data_new.vMCTracks.size() << " tracks ";
  cout << " have been added." << endl;
  for (int i = 0; i < data_new.vStsStrips.size(); i++){
     vStsStrips.push_back(L1Strip(data_new.vStsStrips[i], iEv));
  }
  for (int i = 0; i < data_new.vStsStripsB.size(); i++)
     vStsStripsB.push_back(L1Strip(data_new.vStsStripsB[i], iEv));
  for (int i = 0; i < data_new.vSFlag.size(); i++)
    vSFlag.push_back(data_new.vSFlag[i]);
  for (int i = 0; i < data_new.vSFlagB.size(); i++)
    vSFlagB.push_back(data_new.vSFlagB[i]);
  for (int i = 0; i < data_new.vMCPoints.size(); i++){
    CbmL1MCPoint p = data_new.vMCPoints[i];
    p.ID += MCTrStep*iEv;
    p.mother_ID += MCTrStep*iEv;
    vMCPoints.push_back(p);
  }

  for (int i = 0; i < data_new.vHitMCRef.size(); i++){
    data_new.vHitMCRef[i] += mcPointsSize;
  }
  for (int i = 0; i < data_new.vStsHits.size(); i++)
    data_new.vStsHits[i] = L1StsHit(data_new.vStsHits[i], iEv, stSize, stBSize);

  map<int, int> hitMoveMap, hitMoveMap_new;
  for (int i = 0; i < vHitStore.size(); i++){
    hitMoveMap[i] = i;
  }
  
//   for (int iS = 0; iS < NStations; iS++)
//     cout << "sta: " << iS << " " << StsHitsStopIndex[iS] << " " <<  StsHitsStartIndex[iS]  << endl;
  
  for (int iS = 0; iS < NStations; iS++){
    int beginH = data_new.StsHitsStartIndex[iS];
    int endH = data_new.StsHitsStopIndex[iS];
    int iH2 = StsHitsStartIndex[iS];
    int endH2 = StsHitsStopIndex[iS];

    // sort in an increasing order
    for (int iH = beginH; iH <= endH; iH++){
      CbmL1HitStore hs = data_new.vHitStore[iH];
        // check
//       int i = iH;
//       if (!( (data_new.StsHitsStartIndex[data_new.vHitStore[i].iStation] <= i) && (data_new.StsHitsStopIndex[data_new.vHitStore[i].iStation] >= i) )){
//         cout << " Hit fault. " << endl;
//         cout << data_new.StsHitsStartIndex[data_new.vHitStore[i].iStation] << " " << i << " " << data_new.StsHitsStopIndex[data_new.vHitStore[i].iStation] << " " << endl;
//       }


        // find right place in our array
      CbmL1HitStore hs2 = vHitStore[iH2];
      while ( (hs2.IsLessThen(hs) > 0) && (iH2 <= endH2) ){
        hitMoveMap[iH2 - iH] = iH2;

        // check

        if ( vHitStore[iH2].iStation != iS || data_new.vHitStore[iH].iStation != iS  )
          cout << " Hit fault2. " << endl;
//         int i = iH2;
//       if (!( (StsHitsStartIndex[iS] <= i) && (endH2 >= i) )){
//         cout << " Hit fault2. " << endl;
//         cout << StsHitsStartIndex[iS] << " " << i << " " << endH2 << " " << endl;
//         static int k = 0;
//         k++;
//         cout << i << " " << k << endl;
//       }
//       if (!( (StsHitsStartIndex[vHitStore[i].iStation] <= i) && (endH2 >= i) )){
//         cout << " Hit fault2. " << endl;
//         cout << StsHitsStartIndex[vHitStore[i].iStation] << " " << i << " " << endH2 << " " << endl;
//         static int k = 0;
//         k++;
//         cout << i << " " << k << endl;
//       }
        
        iH2++;
        hs2 = vHitStore[iH2];
      }

//       cout << iH2 << endl; // gdb 0
      vHitStore.insert(vHitStore.begin() + iH2, hs);
      vHitMCRef.insert(vHitMCRef.begin() + iH2, data_new.vHitMCRef[iH]);
      vStsHits.insert (vStsHits.begin()  + iH2, data_new.vStsHits[iH]);
      hitMoveMap_new[iH] = iH2;
//       if ( vHitStore[iH2].iStation != iS )
//         cout << " Hit fault2. iH2 " << endl;
//       if ( data_new.vHitStore[iH].iStation != iS )
//         cout << " Hit fault2. iH " << endl;
      endH2++;
    }; // iH

    int nH = endH - beginH + 1;
    StsHitsStopIndex[iS] += nH;
//     cout << endH2 << " " << StsHitsStopIndex[iS] << endl;

    
    for (int iS2 = iS + 1; iS2 < NStations; iS2++){
      StsHitsStartIndex[iS2] += nH;
      StsHitsStopIndex[iS2] += nH;
    }

//     cout << "sta: " << iS << " " << data_new.StsHitsStopIndex[iS] << " " <<  data_new.StsHitsStartIndex[iS] /*<< nH data_new.StsHitsStopIndex[iS] - data_new.StsHitsStartIndex[iS] + 1*/ << " hits have been added" << endl;
//     cout << "sta: " << iS << " " << StsHitsStopIndex[iS] << " " <<  StsHitsStartIndex[iS]  << endl;
  } // iS

    // correct all tracks
  for (int i = 0; i < vMCTracks.size(); i++){
    CbmL1MCTrack &t = vMCTracks[i];
    for (int j = 0; j < t.StsHits.size(); j++){
      t.StsHits[j] = hitMoveMap[t.StsHits[j]];
    }
  }
    // correct and add new tracks
  for (int i = 0; i < data_new.vMCTracks.size(); i++){
    CbmL1MCTrack t = data_new.vMCTracks[i];
    t.ID += MCTrStep*iEv;
    if (t.mother_ID >= 0) t.mother_ID += MCTrStep*iEv;
    for (int j = 0; j < t.StsHits.size(); j++){
//       cout << t.StsHits[j] << " ";
      t.StsHits[j] = hitMoveMap_new[t.StsHits[j]];
//       cout << t.StsHits[j] << endl;
    }
    vMCTracks.push_back(t);
  }
}



void ReadStAPAlgoData(TData &data)
{
  static int iEvent = 1;
  static ifstream fadata;
  static char fname[100];
  
  {

    if ( iEvent == 1 ){
      strcpy(fname, work_dir);
      strcat(fname, "data_algo.txt");
      fadata.open(fname);
    };

    cout << " -I- begin read " <<  iEvent << " event algo data from file " << fname << endl;
        
    data.vStsHits.clear();
    data.vStsStrips.clear();
    data.vStsStripsB.clear();
    data.vSFlag.clear();
    data.vSFlagB.clear();
      // check correct position in file
    char s[] = "EVENT:  ";
    int nEv;
    fadata >> s;
    fadata >> nEv;
    if (nEv != iEvent)  cout << " -E- Can't read event number " << iEvent << " from file " << fname << endl;
//  cout << s << "|" << nEv << endl; int kk; cin >> kk;
    int n;  // number of elements
      // read data.vStsStrips
    fadata >> n;
    for (int i = 0; i < n; i++){
      fscal element;
      fadata >> element;
      data.vStsStrips.push_back(element);
    };
      // read data.vStsStripsB
    fadata >> n;
    for (int i = 0; i < n; i++){
      fscal element;
      fadata >> element;
      data.vStsStripsB.push_back(element);
    };
      // read data.vSFlag
    fadata >> n;
    for (int i = 0; i < n; i++){
      int element;
      fadata >> element;
      data.vSFlag.push_back((unsigned char)element);
    };
      // read data.vSFlagB
    fadata >> n;
    for (int i = 0; i < n; i++){
      int element;
      fadata >> element;
      data.vSFlagB.push_back((unsigned char)element);
    };
      // read data.vStsHits
    fadata >> n;
    int element_f;  // for convert
    int element_b;
    for (int i = 0; i < n; i++){
      L1StsHit element;
      fadata >> element_f >> element_b;
      element.f = (unsigned short int) element_f;
      element.b = (unsigned short int) element_b;
      data.vStsHits.push_back(element);
    };
      // read StsHitsStartIndex and StsHitsStopIndex
    n = 20;
    for (int i = 0; i < n; i++){
      fadata >> data.StsHitsStartIndex[i];
    };
    for (int i = 0; i < n; i++){
      fadata >> data.StsHitsStopIndex[i];
    };

    cout << iEvent << " event data read from file " << fname << " successfully" << endl;
//    if (iEvent == maxNEvent) fadata.close();  
  }
  iEvent++;
};

    /// read data for performance(vMCPoints, vMCTracks, vHitMCRef, vHitStore) of one event from file
void ReadStAPPerfoData(TData &data)
{

  static int iEvent = 1;
  static ifstream fadata;
  static char fname[100];

  {
  
    if ( iEvent == 1 ){       // file open on begin of all work class and close at end
      strcpy(fname, work_dir);
      strcat(fname, "data_perfo.txt");
      fadata.open(fname);
    };
    data.vMCPoints.clear();
    data.vMCTracks.clear();
    data.vHitStore.clear();
    data.vHitMCRef.clear();
      // check if it is right position in file
    char s[] = "Event:  ";  // buffer
    int nEv;                // event number
    fadata >> s;
    fadata >> nEv;
    if (nEv != iEvent)  cout << " -E- Performance: can't read event number " << iEvent << " from file " << fname << endl;
      // vMCPoints
    int n;  // number of elements
    fadata >> n;
    for (int i = 0; i < n; i++){
      CbmL1MCPoint element;
  
      fadata >> element.x;
      fadata >> element.y;
      fadata >> element.z;
      fadata >> element.px;
      fadata >> element.py;
      fadata >> element.pz;
      fadata >> element.p;
      fadata >> element.q;
      fadata >> element.mass;
      
      fadata >> element.pdg;
      fadata >> element.ID;
      fadata >> element.mother_ID;
      fadata >> element.iStation;
      data.vMCPoints.push_back(element);
    };
    
      // vMCTracks . without Points
    fadata >> n;
    for (int i = 0; i < n; i++){
      CbmL1MCTrack element;
  
      fadata >> element.x;
      fadata >> element.y;
      fadata >> element.z;
      fadata >> element.px;
      fadata >> element.py;
      fadata >> element.pz;
      fadata >> element.p;
      fadata >> element.q;
      fadata >> element.mass;
      
      fadata >> element.pdg;
      fadata >> element.ID;
      fadata >> element.mother_ID;
      
      int nhits;
      fadata >> nhits;
      for (int k = 0; k < nhits; k++){
        int helement;
        fadata >> helement;
        element.StsHits.push_back(helement);
//         cout << helement << endl; // gdb 1
      };
      
      fadata >> element.nMCContStations;
      fadata >> element.nHitContStations;
      data.vMCTracks.push_back(element);
    };

//         // find the corrspondece between MCPoint MCTrack
//     
//     for (unsigned k = 0; k < data.vMCPoints.size(); k++){
//           // find track with ID
//       int itr = 0;
//       for (;data.vMCTracks[itr].ID != data.vMCPoints[k].ID;itr++);
//       data.vMCTracks[itr].Points.push_back(&data.vMCPoints[k]);
//     };
    
        // vHitMCRef
    fadata >> n;
    for (int i = 0; i < n; i++){
      int element;
      fadata >> element;
      data.vHitMCRef.push_back(element);
    };

      // vHitStore
    fadata >> n;
    for (int i = 0; i < n; i++){
      CbmL1HitStore element;
      fadata >> element.ExtIndex;
      fadata >> element.iStation;
      
      fadata >> element.x;
      fadata >> element.y;
      data.vHitStore.push_back(element);
    };
      
    cout << " -I- Performance: data read from file " << fname << " successfully"<< endl;
  }
  iEvent++;
};  // Performance::ReadOneEvData()

void WriteStAPAlgoData(TData &data_out)  // must be called after ReadEvent
{
    // write algo data in file
  static int vNEvent = 1;
  fstream fadata;
  static char fadata_name[100];

  strcpy(fadata_name, out_dir);
  strcat(fadata_name, "data_algo.txt");

//    if ( vNEvent <= maxNEvent ) {
  if ( 1 ) {
 
    if (vNEvent == 1)
      fadata.open(fadata_name,fstream::out);  // begin new file
    else
      fadata.open(fadata_name,fstream::out | fstream::app);
        
    fadata << "Event:" << " ";
    fadata << vNEvent << endl;
      // write vStsStrips
    int n = data_out.vStsStrips.size();  // number of elements
    fadata << n << endl;
    for (int i = 0; i < n; i++){
//       fadata << data_out.vStsStrips[i] << endl;
      fadata << data_out.vStsStrips[i].f  << " " << data_out.vStsStrips[i].n << endl;
    };
      // write vStsStripsB
    n = data_out.vStsStripsB.size();
    fadata << n << endl;
    for (int i = 0; i < n; i++){
//       fadata << data_out.vStsStripsB[i] << endl;
      fadata << data_out.vStsStripsB[i].f  << " " << data_out.vStsStripsB[i].n << endl;
    };
//       // write vStsZPos
//     n = data_out.vStsZPos.size();
//     fadata << n << endl;
//     for (int i = 0; i < n; i++){
//       fadata << data_out.vStsZPos[i] << endl;
//     };
      // write vSFlag
    n = data_out.vSFlag.size();
    fadata << n << endl;
    unsigned char element;
    for (int i = 0; i < n; i++){
      element = data_out.vSFlag[i];
      fadata << (int)element << endl;
    };
      // write vSFlagB
    n = data_out.vSFlagB.size();
    fadata << n << endl;
    for (int i = 0; i < n; i++){
      element = data_out.vSFlagB[i];
      fadata << (int)element << endl;
    };
      // write vStsHits
    n = data_out.vStsHits.size();
    fadata << n << endl;
    unsigned /*short*/ int element2;
//     char element3;
    for (int i = 0; i < n; i++){
      element2 = data_out.vStsHits[i].f;
      fadata  << (int)element2 << " ";
      element2 = data_out.vStsHits[i].b;
      fadata  << (int)element2 << " ";
      element2 = data_out.vStsHits[i].n;
      fadata  << (int)element2 << " ";
//       element3 = data_out.vStsHits[i].iz;
//       fadata  << (int)element3 << endl;
    };
      // write StsHitsStartIndex and StsHitsStopIndex
    n = 20;
    for (int i = 0; i < n; i++){
      fadata  << data_out.StsHitsStartIndex[i] << endl;
    };
    for (int i = 0; i < n; i++){
      fadata  << data_out.StsHitsStopIndex[i] << endl;
    };
    
    fadata.close();
  }
  cout << "-I- CbmL1: CATrackFinder data for event number " << vNEvent << " have been written in file " << fadata_name << endl;
  vNEvent++;
}

void WriteStAPPerfData(TData &data_out)
{
  fstream fpdata;
  static int vNEvent = 1;
  
  static char fpdata_name[100];

  strcpy(fpdata_name, out_dir);
  strcat(fpdata_name, "data_perfo.txt");


        // write data for performance in file
  //   if ( vNEvent <= maxNEvent )  {
  if ( 1 ) {

    if (vNEvent == 1)
      fpdata.open(fpdata_name,fstream::out);  // begin new file
    else
      fpdata.open(fpdata_name,fstream::out | fstream::app);
  
    fpdata << "Event: " ;
    fpdata << vNEvent << endl;
      // write vMCPoints
    int n = data_out.vMCPoints.size();  // number of elements
    fpdata << n << endl;
    for (int i = 0; i < n; i++){
      fpdata << data_out.vMCPoints[i].x << " ";
      fpdata << data_out.vMCPoints[i].y << " ";
      fpdata << data_out.vMCPoints[i].z << "  ";
      fpdata << data_out.vMCPoints[i].px << " ";
      fpdata << data_out.vMCPoints[i].py << " ";
      fpdata << data_out.vMCPoints[i].pz << " ";
      fpdata << data_out.vMCPoints[i].p << "  ";
      fpdata << data_out.vMCPoints[i].q << " ";
      fpdata << data_out.vMCPoints[i].mass << "   ";
      
      fpdata << data_out.vMCPoints[i].pdg << " ";
      fpdata << data_out.vMCPoints[i].ID << " ";
      fpdata << data_out.vMCPoints[i].mother_ID << " ";
      fpdata << data_out.vMCPoints[i].iStation << endl;
    };
    
          // write vMCTracks  . without Points
    n = data_out.vMCTracks.size();  // number of elements
    fpdata << n << endl;
    for (int i = 0; i < n; i++){
      fpdata << data_out.vMCTracks[i].x << " ";
      fpdata << data_out.vMCTracks[i].y << " ";
      fpdata << data_out.vMCTracks[i].z << "  ";
      fpdata << data_out.vMCTracks[i].px << " ";
      fpdata << data_out.vMCTracks[i].py << " ";
      fpdata << data_out.vMCTracks[i].pz << " ";
      fpdata << data_out.vMCTracks[i].p << "  ";
      fpdata << data_out.vMCTracks[i].q << " ";
      fpdata << data_out.vMCTracks[i].mass << "   ";
      
      fpdata << data_out.vMCTracks[i].pdg << " ";
      fpdata << data_out.vMCTracks[i].ID << " ";
      fpdata << data_out.vMCTracks[i].mother_ID << endl;
      
      int nhits = data_out.vMCTracks[i].StsHits.size();
      fpdata  << "   " << nhits << endl << "   ";
      for (int k = 0; k < nhits; k++){
        fpdata << data_out.vMCTracks[i].StsHits[k] << " ";
      };
      fpdata << endl;

      fpdata << data_out.vMCTracks[i].nMCContStations << " ";
      fpdata << data_out.vMCTracks[i].nHitContStations << endl;
    };
    
        // write vHitMCRef
    n = data_out.vHitMCRef.size();  // number of elements
    fpdata << n << endl;
    for (int i = 0; i < n; i++){
      fpdata << data_out.vHitMCRef[i] << endl;
    };
        
      // write vHitStore
    n = data_out.vHitStore.size();  // number of elements
    fpdata << n << endl;
    int stat_sf = 0; // check
    for (int i = 0; i < n; i++){
      fpdata << data_out.vHitStore[i].ExtIndex << "  ";
      fpdata << data_out.vHitStore[i].iStation << "  ";
      
      fpdata << data_out.vHitStore[i].x << " ";
      fpdata << data_out.vHitStore[i].y << endl;

        // check
//       cout << data_out.vHitStore[i].iStation << " ";
      if (!( (data_out.StsHitsStartIndex[data_out.vHitStore[i].iStation] <= i) && (data_out.StsHitsStopIndex[data_out.vHitStore[i].iStation] >= i) )){
//         cout << " Hit sort fault. " << " ";
//         cout << data_out.StsHitsStartIndex[data_out.vHitStore[i].iStation] << " " << i << " " << data_out.StsHitsStopIndex[data_out.vHitStore[i].iStation] << " " << endl;
        
        stat_sf++;
      }
    };
    cout << " NSort faults: " << stat_sf << " from " << n << endl;// check

    fpdata.close();
  }
  cout << "-I- CbmL1: Data for performance of event number " << vNEvent << " have been written in file " << fpdata_name << endl;
  vNEvent++;
    
}

  
int main(int argc, char *argv[])
{
  int NMergeEvents = 2;
  int NEvents = 1;
  
  if (argc >= 2){
    NMergeEvents = atoi(argv[1]);
  }

  if (argc >= 3){
    NEvents = atoi(argv[2]);
  }
    // skip
//   for (int igEv = 0; igEv < 15; igEv++ ){
//     TData data;
//     ReadStAPAlgoData(data);
//     ReadStAPPerfoData(data);
//   }
  
  for (int igEv = 0; igEv < NEvents; igEv++ ){
    TData data;
    TData data_out;

    ReadStAPAlgoData(data);
    ReadStAPPerfoData(data);
    data.Correct();
//     while (! data.Correct()){
//       ReadStAPAlgoData(data);
//       ReadStAPPerfoData(data);
//     }
    data_out = data;
    for (int iEv = 1; iEv < NMergeEvents; iEv++ ){
      ReadStAPAlgoData(data);
      ReadStAPPerfoData(data);
      data.Correct();
//       while (! data.Correct()){
//         ReadStAPAlgoData(data);
//         ReadStAPPerfoData(data);
//       }
      data_out.Add(data,iEv);
    }

    cout << endl; cout << endl; cout << endl;
    WriteStAPAlgoData(data_out);
    WriteStAPPerfData(data_out);
    
  } // igEv
  
  return 1;
}