/// merge N events into 1.
  /// to compile: g++ convert.cpp -O3 -Icode/Common -o convert
  /// to run:  convert NMergeEvents NGeneratedEvents OutputDir InputDir. If you want merge 20 events into 10, then print: convert 2 10

#define DEBUG
#define MUTE1
#define MUTE2
 
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cstring>
#include <stdio.h> // for atoi
#include <stdlib.h>

#include "TRandom.h"

using namespace std;


typedef  float fscal;

  // #include "code/Algo/L1Algo/L1Types.h"
#include "code/Performance/CbmL1MCPoint.h"
#include "code/L1Algo/L1StsHit.h"
#include "code/L1Algo/L1Strip.h"
#include "code/Performance/CbmL1StsHit.h"
#include "code/Performance/CbmL1MCTrack.h"

  char work_dir[100] = "/d/cbm05/akishina/STAP/centr_25_1e/";  // path to input and output files

//char work_dir[100] = "/d/cbm05/akishina/STAP/mbias_1k/";
  //char work_dir[100] = "/d/cbm02/ikulakov/STAP/centr_0maps_015deg_100ev_MCInOut_13032012/";  // path to input and output files

  //char work_dir[100] = "/d/cbm02/ikulakov/STAP/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] = "/d/cbm05/akishina/STAP/mbias_0maps_015deg_1000ev_13032012_30e_30g_new_track/";
char out_dir[100] = "/d/cbm05/akishina/STAP/centr_25_convert_1event/";
const int NStations = 8;

float EventRate = 1e7f;
float tau = 1.e9f/EventRate;  // average time between events in [ns]
float T0 = 0.f;

float T0_shift = 0.f;

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

struct TData{

  TData() { for(int i=0; i<20; i++) {StsHitsStartIndex[i]=-1; StsHitsStopIndex[i]=0;}}
  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
  
  vector<CbmL1StsHit> vStsHitsPerfo;   
  vector< fscal >   vStsZPos;  
  const THitI* StsHitsStartIndex_;
  const THitI* StsHitsStopIndex_;
  
  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;
    //   };

  void Reserve(int nEv);
};

void TData::Reserve(int nEv)
{
  const int k = 1.2*nEv;
  
#define res(AAAA) AAAA.reserve(AAAA.size()*k)
  
  res(vStsStrips);
  res(vStsStripsB);
  res(vStsHits);
  res(vSFlag);
  res(vSFlagB);
  res(vMCPoints);
  res(vMCTracks);
  res(vHitMCRef);
  res(vHitStore);
  res(vStsHitsPerfo);
#undef res
}

int TData::Correct()
{
  int flag = 1;
    //  int curStart = -1, 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];
    //  };

  int curStart = 0;
  for(int i = 0; i<NStations; i++ )
  {
    if( StsHitsStartIndex[i] < curStart )
    {
      StsHitsStartIndex[i] = curStart;
    }
    else curStart = StsHitsStartIndex[i];

    if( StsHitsStopIndex[i] < curStart )
    {
        
      StsHitsStopIndex[i] = curStart;
    }
    else curStart = StsHitsStopIndex[i];
  }

  return flag;
}

void TData::Add(TData& data_new, int iEv)
{
#ifndef MUTE0
  cout <<  " Add"<<endl;
#endif
  const int stSize =  vStsStrips.size();
  const int stBSize =  vStsStripsB.size();
  const int mcPointsSize =  vMCPoints.size();
  vector< fscal >   vStsZPos_old;  

#ifndef MUTE0
  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 << " will be added." << endl;
#endif  
  
  for (unsigned int i = 0; i < data_new.vStsStrips.size(); i++)
    vStsStrips.push_back(L1Strip(data_new.vStsStrips[i], iEv));
  for (unsigned int i = 0; i < data_new.vStsStripsB.size(); i++)
    vStsStripsB.push_back(L1Strip(data_new.vStsStripsB[i], iEv));
  for (unsigned int i = 0; i < data_new.vSFlag.size(); i++)
    vSFlag.push_back(data_new.vSFlag[i]);
  for (unsigned int i = 0; i < data_new.vSFlagB.size(); i++)
    vSFlagB.push_back(data_new.vSFlagB[i]);
  vStsZPos_old.resize(vStsZPos.size());
  for (unsigned int i = 0; i < vStsZPos.size(); i++)
    vStsZPos_old[i] = vStsZPos[i];

  for (unsigned int i = 0; i < data_new.vStsZPos.size(); i++) {
      //   cout<<data_new.vStsZPos[i]<<" data_new.vStsZPos[i]"<<endl;
    bool flag=0;
      // find between old z-poz
    for (unsigned int m = 0; m < vStsZPos_old.size(); m++) {
      if (vStsZPos_old[m]==data_new.vStsZPos[i]) {
	flag =1;
	break;
      }
    }
    if (!flag)
      vStsZPos.push_back(data_new.vStsZPos[i]);
  }
   
  std::sort(vStsZPos.begin(), vStsZPos.end());
  
    //   for (unsigned int i = 0; i < vStsZPos.size(); i++) {
    //    cout<<vStsZPos[i]<<" vStsZPos[i]"<<endl;}
  
  int mcTrackIndexOffset = vMCTracks.size();
  
  for (unsigned int i = 0; i < data_new.vMCPoints.size(); i++){
    CbmL1MCPoint &p = data_new.vMCPoints[i];
    p.ID += mcTrackIndexOffset;
    p.mother_ID += mcTrackIndexOffset;
    p.time += T0;
    for (unsigned int j = 0; j < p.hitIds.size(); j++){
      p.hitIds[j] += vStsHitsPerfo.size();//PointMoveMap[p.hitIds[j]+1000000];
    }
    vMCPoints.push_back(p);
  }
    //   for (unsigned int i = 0; i < vMCPoints.size(); i++){
    //     if vMCPoints.hitIds[i]< 1000000
    //       
    //   
    //   }

  
  for (unsigned int i = 0; i < data_new.vStsHits.size(); i++)
  {
    int iz_cur=data_new.vStsHits[i].iz;
    float time = data_new.vStsHits[i].t_reco;// + T0;
    
    data_new.vStsHits[i] = L1StsHit(data_new.vStsHits[i], stSize, stBSize);
    data_new.vStsHits[i].n = iEv;
    data_new.vStsHits[i].iz=iz_cur; // ???
    data_new.vStsHits[i].t_reco = time;
  }

  for (unsigned int i = 0; i < data_new.vHitMCRef.size(); i++){
    data_new.vHitMCRef[i] += mcPointsSize;
  }
 

  
 
 
  
 
  for(unsigned int iH=0; iH<vStsHits.size(); iH++){
    fscal zCur = vStsZPos_old[ vStsHits[iH].iz ];
    int zIndNew = -1;
    for(unsigned int iZ=0; iZ<vStsZPos.size(); iZ++){
        //        std::cout << "tatata " << zCur << " " << vStsZPos[iZ] << std::endl;
      if(vStsZPos[iZ]==zCur) {
        zIndNew=iZ;
        break;
      }
    }
    vStsHits[iH].iz=zIndNew;
  }
  
  for(unsigned int iH=0; iH<data_new.vStsHits.size(); iH++){
    fscal zCur = data_new.vStsZPos[ data_new.vStsHits[iH].iz ];
    int zIndNew = -1;
    for(unsigned int iZ=0; iZ<vStsZPos.size(); iZ++){
        //        std::cout << "tatatatatatatata " << zCur << " " << vStsZPos[iZ] << std::endl;

      if(vStsZPos[iZ]==zCur) {
        zIndNew=iZ;
        break;
      }
    }
    data_new.vStsHits[iH].iz=zIndNew;
    assert (zIndNew>=0);
  }



  for (unsigned int i = 0; i < data_new.vStsHitsPerfo.size(); i++){
    CbmL1StsHit &p = data_new.vStsHitsPerfo[i];
      //p.hitId = hitMoveMap_new[p.hitId];

    for (unsigned int j = 0; j < p.mcPointIds.size(); j++){
      p.mcPointIds[j] += mcPointsSize;//StsPerfoMoveMap[p.mcPointIds[j]];
    }
      // vStsHitsPerfo.push_back(p);
  }

  vector<int> hitMoveMap(vHitStore.size()), hitMoveMap_new(data_new.vHitStore.size());
  for (unsigned int i = 0; i < vHitStore.size(); i++){
    hitMoveMap[i] = i;
  }

  vector<CbmL1HitStore> vHitStore_cmb(vHitStore.size()+data_new.vHitStore.size());
  vector<int> vHitMCRef_cmb(vHitMCRef.size()+data_new.vHitMCRef.size());
  vector<L1StsHit> vStsHits_cmb(vStsHits.size()+data_new.vStsHits.size());
  vector<CbmL1StsHit> vStsHitsPerfo_cmb(vStsHitsPerfo.size()+data_new.vStsHitsPerfo.size());
  for (int iS = 0; iS < NStations; iS++){
    int beginH = data_new.StsHitsStartIndex[iS];
    int endH = data_new.StsHitsStopIndex[iS];
    int beginH2 = StsHitsStartIndex[iS];
    int endH2 = StsHitsStopIndex[iS];

      // sort in an increasing order
    int iH = beginH, iH2 = beginH2;
    for (; (iH < endH) && (iH2 < endH2); ){
      L1_ASSERT(data_new.vHitStore[iH].y <= data_new.vHitStore[iH+1].y || iH + 1 >= endH, data_new.vHitStore[iH].y << " " << data_new.vHitStore[iH+1].y << " " << iH << " " << endH);
      int iHc = iH+iH2;
      
      CbmL1HitStore hs = data_new.vHitStore[iH];
      CbmL1HitStore hs2 = vHitStore[iH2];

      if (hs2.IsLessThan(hs) > 0){
          // add hs2
        hitMoveMap[iH2] = iHc;
        vHitStore_cmb[iHc] = hs2;
        vHitMCRef_cmb[iHc] = vHitMCRef[iH2];
        vStsHits_cmb[iHc] = vStsHits[iH2];
        vStsHitsPerfo_cmb[iHc] = vStsHitsPerfo[iH2];
                
        iH2++;
      }
      else {
          // add hs
        hitMoveMap_new[iH] = iHc;
        vHitStore_cmb[iHc] = hs;
        vHitMCRef_cmb[iHc] = data_new.vHitMCRef[iH];
        vStsHits_cmb[iHc] = data_new.vStsHits[iH];
        vStsHitsPerfo_cmb[iHc] = data_new.vStsHitsPerfo[iH];
        
        iH++;
      }

    }; // iH
      // add the rest
    for (; (iH < endH) || (iH2 < endH2); ){
      int iHc = iH+iH2;

      if (iH2 < endH2){
          // add hs2
        hitMoveMap[iH2] = iHc;
        vHitStore_cmb[iHc] = vHitStore[iH2];
        vHitMCRef_cmb[iHc] = vHitMCRef[iH2];
        vStsHits_cmb[iHc] = vStsHits[iH2];
        vStsHitsPerfo_cmb[iHc] = vStsHitsPerfo[iH2];
                
        iH2++;
      }
      else {
          // add hs
        hitMoveMap_new[iH] = iHc;
        vHitStore_cmb[iHc] = data_new.vHitStore[iH];
        vHitMCRef_cmb[iHc] = data_new.vHitMCRef[iH];
        vStsHits_cmb[iHc] = data_new.vStsHits[iH];
        vStsHitsPerfo_cmb[iHc] = data_new.vStsHitsPerfo[iH];
        
        iH++;
      }

    }; // iH
  } // iS

  for (int iS = 0; iS < NStations; iS++){
    int beginH = data_new.StsHitsStartIndex[iS];
    int endH = data_new.StsHitsStopIndex[iS];

      //    cout <<(endH2-StsHitsStopIndex[iS])<<" (endH2-StsHitsStopIndex[iS])" <<endl;

    int nH = endH - beginH;
      //    cout << "!!!!!!!   " << iS << " " << StsHitsStartIndex[iS] << " " << StsHitsStopIndex[iS] << " " << nH << endl;

      //    if( (StsHitsStartIndex[iS] == -1) && (nH >= 0) )
      //    {
      //      StsHitsStartIndex[iS] = 0;
      //      for(int iS2 = iS + 1; iS2 < NStations; iS2++)
      //        if(StsHitsStartIndex[iS2] == -1) StsHitsStartIndex[iS2] = 0;
      //    }
      //   if(nH>0)
    {
      StsHitsStopIndex[iS] += nH;
        //    cout <<nH<<" nH" <<endl;
        //     cout << endH2 << " " << StsHitsStopIndex[iS] << endl;

    
      for (int iS2 = iS + 1; iS2 < NStations; iS2++){
        StsHitsStartIndex[iS2] += nH;
        StsHitsStopIndex[iS2] += nH;
          //        cout << "???   " << StsHitsStartIndex[iS2] << " " << StsHitsStopIndex[iS2] << " " << nH << endl;
      }
    }
      // cout<< vHitStore[StsHitsStartIndex[iS] ].iStation<<" vHitStore[StsHitsStartIndex[iS] ].iStation "<<iS<<" iS "<<endl;
      // cout<< vHitStore[StsHitsStopIndex[iS] ].iStation<<" vHitStore[StsHitsStopIndex[iS] ].iStation "<<iS<<" iS "<<endl;
      // cout<< vHitStore[StsHitsStopIndex[iS]-1 ].iStation<<" vHitStore[StsHitsStopIndex[iS] -1].iStation "<<iS<<" iS "<<endl;
      //     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

  vHitStore = vHitStore_cmb;
  vHitMCRef = vHitMCRef_cmb;
  vStsHits = vStsHits_cmb;
  vStsHitsPerfo = vStsHitsPerfo_cmb;


    //  cout <<" GEST"<<endl; 
    // correct all tracks

  for (unsigned int i = 0; i < vStsHitsPerfo.size(); i++){
    CbmL1StsHit &p = vStsHitsPerfo[i];
    p.hitId = i;//hitMoveMap[p.hitId];
  }



  for (unsigned int i = 0; i < vMCTracks.size(); i++){
    CbmL1MCTrack &t = vMCTracks[i];
    for (unsigned int j = 0; j < t.StsHits.size(); j++){
      t.StsHits[j] = hitMoveMap[t.StsHits[j]];
    }
      //     for (int l = 0; l < t.Points.size(); l++){
      //       t.Points[l] = pointMoveMap_new[t.Points[l]];
      //     }
  }

    // correct and add new tracks
  for (unsigned int i = 0; i < data_new.vMCTracks.size(); i++){
    CbmL1MCTrack &t = data_new.vMCTracks[i];
    t.time += T0;
    t.ID += mcTrackIndexOffset;
    if (t.mother_ID >= 0) t.mother_ID += mcTrackIndexOffset;
    else t.mother_ID = -iEv-1;
    for (unsigned 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;
    }
    
    for (unsigned int l = 0; l < t.Points.size(); l++){
      t.Points[l] += mcPointsSize;
    }
    
    vMCTracks.push_back(t);
  }
}



void ReadStAPAlgoData(TData &data)
{
#ifndef MUTE0
  cout <<  " ReadStAPAlgoData"<<endl;
#endif
  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);
    };

#ifndef MUTE1
    cout << " -I- begin read " <<  iEvent << " event algo data from file " << fname << endl;
#endif
        
    data.vStsHits.clear();
    data.vStsStrips.clear();
    data.vStsStripsB.clear();
    data.vSFlag.clear();
    data.vSFlagB.clear();
    data.vStsZPos.clear();

    for(int i=0; i<20; i++) {data.StsHitsStartIndex[i]=-2; data.StsHitsStopIndex[i]=0;}

      // check correct position in file
    char s[] = "EVENT:  ";
    int nEv;
    fadata >> s;
      // cout << s<<  " s"<<endl;
    fadata >> nEv;
      //  cout << nEv<<  " nEv"<<endl;
    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;
      //       cout << n<<  " vStsStrips"<<endl;
    for (int i = 0; i < n; i++){
      fscal element;
    
        //  unsigned short int element2;
        //  fadata >> element>> element2;
        // vStsStrips.push_back(L1Strip(element,element2));
      fadata >> element;
      data.vStsStrips.push_back(L1Strip(element));
    };
      // read data.vStsStripsB
    fadata >> n;
      //  cout << n<<  " vStsStripsB"<<endl;
    
    for (int i = 0; i < n; i++){
      fscal element;
        //   unsigned short int element2;
        //  fadata >> element >> element2;
        // vStsStripsB.push_back(L1Strip(element,element2));
      fadata >> element;
      data.vStsStripsB.push_back(L1Strip(element));
    };
      // read data.vSFlag

    fadata >> n;
      //    cout << n<<  " vStsZPos"<<endl;
       
    for (int i = 0; i < n; i++){
      fscal element;
      fadata >> element;
        //   cout<<element<<" element Z"<<endl;
      data.vStsZPos.push_back(element);
    };

#ifndef MUTE2
    printf("In main after ReadStAPAlgoData: data.vStsZPos.size() = %d\n", data.vStsZPos.size() );
#endif
    fadata >> n;
      //   cout << n<<  " vSFlag"<<endl;
    for (int i = 0; i < n; i++){
      int element;
      fadata >> element;
        //   cout << element<< " vSFlag"<<endl;
      data.vSFlag.push_back((unsigned char)element);
    };
   
      // read data.vSFlagB
    fadata >> n;
      //    cout << n<<  " vSFlagB"<<endl;
    for (int i = 0; i < n; i++){
      int element;
      fadata >> element;
        //     cout << element<< " vSFlagB"<<endl;
      data.vSFlagB.push_back((unsigned char)element);
    };
    
    fadata >> n;
      //    cout << n<<  " vStsHits"<<endl;
    int element_f;  // for convert
    int element_b;
    int element_iz;
    float time;
    for (int i = 0; i < n; i++){
      L1StsHit element;
      fadata >> element_f >> element_b >> element_iz >> time;
        //     cout << element_f<< " "<<element_b<<" "<<element_iz<< " element_iz"<<endl;
      element.f = static_cast<THitI>(element_f);
      element.b = static_cast<THitI>(element_b);
      element.iz = static_cast<TZPosI>(element_iz);
      element.t_reco = time;
      data.vStsHits.push_back(element);
    }
     
    n = 20;
    int MaxNStations = NStations;
    for (int i = 0; i < n; i++){
      int tmp;
      fadata >> tmp;
      if (MaxNStations+1 > i) data.StsHitsStartIndex[i] = tmp;
#ifndef MUTE2
      cout << tmp<<  " tmp"<<endl;
#endif
    }
    for (int i = 0; i < n; i++){
      int tmp;
      fadata >> tmp;
      if (MaxNStations+1 > i) data.StsHitsStopIndex[i] = tmp;
#ifndef MUTE2
      cout << tmp<<  " tmp"<<endl;
#endif
    }
   
    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)
{
#ifndef MUTE0
  cout <<  "  ReadStAPPerfoDataa"<<endl;
#endif
  static int iEvent = 1;
  static ifstream fpdata;
  static char fname[100];
  const int iVerbose = 0;

  {
  
    if ( iEvent == 1 ){       // file open on begin of all work class and close at end
      strcpy(fname, work_dir);
      strcat(fname, "data_perfo.txt");
      fpdata.open(fname);
    };        
    data.vMCPoints.clear();
    data.vMCTracks.clear();
    data.vHitStore.clear();
    data.vHitMCRef.clear();
    data.vStsHitsPerfo.clear();
      // check if it is right position in file
    char s[] = "Event:  ";  // buffer
    int nEv=0;                // event number
    fpdata >> s;
    
    fpdata >> nEv;
      //   cout << nEv<< " nEv"<<endl;
    if (nEv != iEvent)  cout << " -E- Performance: can't read event number " << iEvent << " from file " << fname << endl;
      // vMCPoints
    int n;  // number of elements
    fpdata >> n;
      //    cout <<n<<  " vMCPoints"<<endl; 
    for (int i = 0; i < n; i++){
      CbmL1MCPoint element;
  
   

      fpdata >> element.xIn;
      fpdata >> element.yIn;
      fpdata >> element.zIn;
      fpdata >> element.pxIn;
      fpdata >> element.pyIn;
      fpdata >> element.pzIn;
            
      fpdata >> element.xOut;
      fpdata >> element.yOut;
      fpdata >> element.zOut;
      fpdata >> element.pxOut;
      fpdata >> element.pyOut;
      fpdata >> element.pzOut;

      element.x = (element.xIn + element.xOut)/2;
      element.y = (element.yIn + element.yOut)/2;
      element.z = (element.zIn + element.zOut)/2;
      element.px = (element.pxIn + element.pxOut)/2;
      element.py = (element.pyIn + element.pyOut)/2;
      element.pz = (element.pzIn + element.pzOut)/2;
      
      fpdata >> element.p;
      fpdata >> element.q;
      fpdata >> element.mass;
      fpdata >> element.time;

      fpdata >> element.pdg;
      fpdata >> element.ID;
      fpdata >> element.mother_ID;
      fpdata >> element.iStation;

   
      int nPoints;
        
      fpdata >> nPoints;
    

      for (int k = 0; k < nPoints; k++){
        int helement;
	
        fpdata >> helement;
	 
        element.hitIds.push_back(helement);
	 
      };
      
      data.vMCPoints.push_back(element);
    };
    
    
    
    
    
      // vMCTracks . without Points
    fpdata >> n;
 
      //    cout << n<<" vMCTracks"<< endl;
    for (int i = 0; i < n; i++){
      CbmL1MCTrack element;
  
      fpdata >> element.x;
      fpdata >> element.y;
      fpdata >> element.z;
      fpdata >> element.px;
      fpdata >> element.py;
      fpdata >> element.pz;
      fpdata >> element.p;
      fpdata >> element.q;
      fpdata >> element.mass;
      fpdata >> element.time;

      fpdata >> element.pdg;
      fpdata >> element.ID;
      fpdata >> element.mother_ID;
      
      int nhits;
      fpdata >> nhits;

      for (int k = 0; k < nhits; k++){
        int helement;
        fpdata >> helement;
        element.StsHits.push_back(helement);
          //         cout << helement << endl; // gdb 1
      };
     
    
      fpdata >> nhits;

      for (int k = 0; k < nhits; k++){
        int helement;
        fpdata >> helement;
	
        element.Points.push_back(helement);
      };
        //         // 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]);
        //     };
      fpdata >> element.nMCContStations;
      fpdata >> element.nHitContStations;
      fpdata >> element.maxNStaMC;
      fpdata >> element.maxNSensorMC;
      fpdata >> element.maxNStaHits;
      fpdata >> element.nStations;
      

      
      data.vMCTracks.push_back(element);
    };
    if (iVerbose >= 4) cout << "vMCTracks[" << n << "]" << " have been read." << endl;
      // vHitMCRef
    fpdata >> n;
#ifndef MUTE2
    cout << n<<" vHitMCRef"<< endl;
#endif
    for (int i = 0; i < n; i++){
      int element;
      fpdata >> element;
      data.vHitMCRef.push_back(element);
    };
    if (iVerbose >= 4) cout << "vHitMCRef[" << n << "]" << " have been read." << endl;
      // vHitStore
    fpdata >> n;
#ifndef MUTE2
    cout << n<<" vHitStore"<< endl;
#endif
    for (int i = 0; i < n; i++){
      CbmL1HitStore element;
      
      fpdata >> element.ExtIndex;
      fpdata >> element.iStation;
      
      fpdata >> element.x;
      fpdata >> element.y;
      data.vHitStore.push_back(element);
    };
    
    if (iVerbose >= 4) cout << "vHitStore[" << n << "]" << " have been read." << endl;
    
      // vStsHits
    fpdata >> n;
#ifndef MUTE2
    cout << n<<" vStsHits"<< endl;
#endif
    for (int i = 0; i < n; i++){
      CbmL1StsHit element;
      fpdata >> element.hitId;
      fpdata >> element.extIndex;

      int nPoints;
      fpdata >> nPoints;
      element.mcPointIds.clear();
      
      for (int k = 0; k < nPoints; k++){
        int id;
        fpdata >> id;

        element.mcPointIds.push_back(id);
      };
      data.vStsHitsPerfo.push_back(element);
    };
    if (iVerbose >= 4) cout << "vStsHits[" << n << "]" << " have been read." << endl;

#ifndef MUTE1
    cout << " -I- Performance: data read from file " << fname << " successfully"<< endl;
#endif
  }
  iEvent++;
}  // Performance::ReadOneEvData()
     

void WriteStAPAlgoData(TData &data_out)  // must be called after ReadEvent
    
{
#ifndef MUTE0
  cout <<  "  WriteStAPAlgoData"<<endl;
#endif
    // 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;
#ifndef MUTE2
    cout << vNEvent <<  " vNEvents"<<endl;
#endif
      // write vStsStrips
      
      
    
      
    int n = data_out.vStsStrips.size();  // number of elements
   
    fadata << n << endl;
      //    cout << n <<  " Strip size"<<endl;  
    
    for (int i = 0; i < n; i++){
      fadata << data_out.vStsStrips[i].f << endl;
//       fadata << data_out.vStsStrips[i].f  << " " << data_out.vStsStrips[i].n << endl;
        //  fadata << data_out.vStsStrips[i].f<<endl;
    };
      // write vStsStripsB
    n = data_out.vStsStripsB.size();
    fadata << n << endl;
    
      //  cout << n <<  " StripB size"<<endl;    
    
    for (int i = 0; i < n; i++){
      fadata << data_out.vStsStripsB[i].f << endl;
//       fadata << data_out.vStsStripsB[i].f  << " " << data_out.vStsStripsB[i].n << endl;
        // fadata << data_out.vStsStripsB[i].f<< endl;
    };
    
      //       // write vStsZPos
    n = data_out.vStsZPos.size();
    
    fadata << n << endl;
      //     cout << n <<  "vStsZPos.size"<<endl;  
    for (int i = 0; i < n; i++){
      fadata << data_out.vStsZPos[i] << endl;
        //     cout << data_out.vStsZPos[i] <<  " data_out.vStsZPos[i]"<<endl;  
    };
      // write vSFlag
    n = data_out.vSFlag.size();
    fadata << n << endl;
      //  cout << n <<  "vSFlag.size"<<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;
      //  cout << n <<  "vSFlagB.size"<<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;
      //   cout << n <<  "vStsHits.size"<<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 << " ";
      element2 = data_out.vStsHits[i].iz;
      fadata  << (int)element2 << " ";
      data_out.vStsHits[i].t_reco += T0_shift;
      fadata  << data_out.vStsHits[i].t_reco << endl;
    };
      // write StsHitsStartIndex and StsHitsStopIndex
    n = 20;
    for (int i = 0; i < n; i++){

      fadata  << data_out.StsHitsStartIndex[i] << endl;
#ifndef MUTE2
      cout<< data_out.StsHitsStartIndex[i] << endl;
#endif
    };
    for (int i = 0; i < n; i++){

      fadata  << data_out.StsHitsStopIndex[i] << endl;
#ifndef MUTE2
      cout<< data_out.StsHitsStopIndex[i] << endl;
#endif
    };
    
    fadata.close();
  }
#ifndef MUTE1
  cout << "-I- CbmL1: CATrackFinder data for event number " << vNEvent << " have been written in file " << fadata_name << endl;
#endif
  vNEvent++;
}

void WriteStAPPerfData(TData &data_out)
{
#ifndef MUTE0
cout <<  "  WriteStAPPerfData"<<endl;
#endif
  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;
#ifndef MUTE2
    cout << vNEvent <<"Event: "<< endl;
#endif
    
      // write vMCPoints
    int n = data_out.vMCPoints.size();  // number of elements
    fpdata << n << endl;
#ifndef MUTE2
    cout << n <<"vMCPoints: "<< endl;
#endif
    for (int i = 0; i < n; i++){
    
      fpdata << data_out.vMCPoints[i].xIn << " ";
      fpdata << data_out.vMCPoints[i].yIn << " ";
      fpdata << data_out.vMCPoints[i].zIn << "  ";
      fpdata << data_out.vMCPoints[i].pxIn << " ";
      fpdata << data_out.vMCPoints[i].pyIn << " ";
      fpdata << data_out.vMCPoints[i].pzIn << " ";
      fpdata << data_out.vMCPoints[i].xOut << " ";
      fpdata << data_out.vMCPoints[i].yOut << " ";
      fpdata << data_out.vMCPoints[i].zOut << "  ";
      fpdata << data_out.vMCPoints[i].pxOut << " ";
      fpdata << data_out.vMCPoints[i].pyOut << " ";
      fpdata << data_out.vMCPoints[i].pzOut << " ";
      fpdata << data_out.vMCPoints[i].p << "  ";
      fpdata << data_out.vMCPoints[i].q << " ";
      fpdata << data_out.vMCPoints[i].mass << " ";  
      data_out.vMCPoints[i].time += T0_shift;
      fpdata << data_out.vMCPoints[i].time << "   ";
      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;
      
      int nhits = data_out.vMCPoints[i].hitIds.size();
      
      fpdata << nhits<< endl<<" ";
      for (int k = 0; k < nhits; k++){
        
        fpdata << data_out.vMCPoints[i].hitIds[k]<<" ";
       
      };
      fpdata<< endl;
    };
    
      //   cout <<  " write vMCPoints"<<endl;
    
      // write vMCTracks  . without Points
      // number of elements
    n = data_out.vMCTracks.size(); 
      //  cout << n <<"vMCTracks: "<< endl;
    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 << " ";
      data_out.vMCTracks[i].time += T0_shift;
      fpdata << data_out.vMCTracks[i].time << "   ";
      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;
      
      int npoints = data_out.vMCTracks[i].Points.size();
      fpdata  << "   " << npoints << endl << "   ";
      for (int k = 0; k < npoints; k++){
        fpdata << data_out.vMCTracks[i].Points[k] << " ";
      };
      
      
 
      fpdata << data_out.vMCTracks[i].nMCContStations << " ";
      fpdata << data_out.vMCTracks[i].nHitContStations<< " ";
      fpdata << data_out.vMCTracks[i].maxNStaMC<< " ";
      fpdata << data_out.vMCTracks[i].maxNSensorMC<< " ";
      fpdata << data_out.vMCTracks[i].maxNStaHits<< endl;
      fpdata << data_out.vMCTracks[i].nStations<< endl;
  
    };
    
      //    cout <<  " write vMCTracks"<<endl;  
    
      // write vHitMCRef
    n = data_out.vHitMCRef.size();  // number of elements
      //   cout << n <<"vHitMCRef: "<< endl;
    fpdata << n << endl;
    for (int i = 0; i < n; i++){
      fpdata << data_out.vHitMCRef[i] << endl;
    };
      //      cout <<  " write vHitMCRef"<<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) ))
	stat_sf++;
        //         cout << " Hit sort fault. " << " ";
        //         cout << data_out.StsHitsStartIndex[data_out.vHitStore[i].iStation] << " " << i << " " << data_out.StsHitsStopIndex[data_out.vHitStore[i].iStation] << " " << endl;
        
       
      
    };
#ifndef MUTE1
    cout << data_out.vHitStore.size()<< " vHitStore.size("<<endl; 
    cout << " NSort faults: " << stat_sf << " from " << n << endl;// check
    cout <<  " write vHitStore"<<endl; 
#endif
            
    int s = data_out.vStsHitsPerfo.size();  // number of elements
    fpdata << s << endl;

    
      // write  vStsHits
    for (int i = 0; i < s; i++){

      fpdata << data_out.vStsHitsPerfo[i].hitId << " ";  
      fpdata << data_out.vStsHitsPerfo[i].extIndex << endl; 
       
      int nPointIds = data_out.vStsHitsPerfo[i].mcPointIds.size();
      fpdata  << "   " << nPointIds << endl << "   ";
      
      for (int j = 0; j < nPointIds;j++)
      {
        fpdata << data_out.vStsHitsPerfo[i].mcPointIds[j]<< " ";
      };
      fpdata << endl;       
    };
      //    cout << s <<"vStsHits: "<< endl;
      //    cout <<  " write vStsHits"<<endl; 
   
    
   

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

  
int main(int argc, char *argv[])
{
  gRandom->SetSeed(1);

  int NMergeEvents = 2;
  int NEvents = 1;

  if (argc >= 2){
    NMergeEvents = atoi(argv[1]);
  }

  if (argc >= 3){
    NEvents = atoi(argv[2]);
  }

  if (argc >= 4){
    strcpy(out_dir,argv[3]);
  }

  if (argc >= 5){
    strcpy(work_dir,argv[4]);
  }
  std::cout << " Create " << NEvents << " groups of events. " << NMergeEvents << " events in group " << std::endl;
  std::cout << " Output Dir = " << out_dir << std::endl;
  std::cout << " Input  Dir = " << work_dir << std::endl;
    // skip
    //   for (int igEv = 0; igEv < 15; igEv++ ){
    //     TData data;
    //     ReadStAPAlgoData(data);
    //     ReadStAPPerfoData(data);
    //   }
  
  for (int igEv = 0; igEv < NEvents; igEv++ ){
    
    T0 = 0;
    
    TData data;
    TData data_out;

    ReadStAPAlgoData(data);
    
    //find shift for timecalculations to avoid negative time values
    
    T0_shift = 1.e8f;
    if(data.vStsHits.size() == 0)
      T0_shift = 0.f;
    for(int iHit=0; iHit<data.vStsHits.size(); iHit++)
    {
      if(data.vStsHits[iHit].t_reco < T0_shift)
        T0_shift = data.vStsHits[iHit].t_reco;
    }
    T0_shift = -T0_shift;
    
#ifndef MUTE1
    cout<<data.vHitStore.size()<< " data.vHitStore.size"<<endl;
#endif
    ReadStAPPerfoData(data);
#ifndef MUTE1
    cout<<data.vHitStore.size()<< " data.vHitStore.size"<<endl;
#endif

    data.Correct();
      //     while (! data.Correct()){
      //       ReadStAPAlgoData(data);
      //       ReadStAPPerfoData(data);
      //     }

   
    data_out = data;
      //     for (int iS = 0; iS < NStations; iS++ )
      //         cout<< data_out.StsHitsStartIndex[iS]<<" data_out.StsHitsStartIndex[iS] "<< data_out.StsHitsStopIndex[iS]<< " data_out.StsHitsStopIndex[iS]"<<endl;

#ifndef MUTE0
  cout << data_out.vStsStrips.size()+data_out.vStsStripsB.size() << " strips and ";
  cout << data_out.vStsHits.size() << " hits for ";
  cout << data_out.vMCTracks.size() << " tracks ";
  cout << " have been read in 1st event in " << igEv << "-th group." << endl;
#endif  
    
    data_out.Reserve(NMergeEvents);
    for (int iEv = 1; iEv < NMergeEvents; iEv++ ){
      
      float deltaT = gRandom->Exp(tau);
      T0 += deltaT;
      
      ReadStAPAlgoData(data);
        //           for (int iS = 0; iS < NStations; iS++ )
        //         cout<< data_out.StsHitsStartIndex[iS]<<" Algo data_out.StsHitsStartIndex[iS] "<< data_out.StsHitsStopIndex[iS]<< " data_out.StsHitsStopIndex[iS]"<<endl;
        //
      ReadStAPPerfoData(data);
        //           for (int iS = 0; iS < NStations; iS++ )
        //         cout<< data_out.StsHitsStartIndex[iS]<<" Perfo data_out.StsHitsStartIndex[iS] "<< data_out.StsHitsStopIndex[iS]<< " data_out.StsHitsStopIndex[iS]"<<endl;
   
      
      data.Correct();
        //     for (int iS = 0; iS < NStations; iS++ )
        //     cout<< data_out.StsHitsStartIndex[iS]<<" Correct data_out.StsHitsStartIndex[iS] "<< data_out.StsHitsStopIndex[iS]<< " data_out.StsHitsStopIndex[iS]"<<endl;
        //       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;
}