/** CbmMuchFindVectorsQA.cxx
 *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
 *@since 2014
 **/
#include "CbmMuchFindVectorsQA.h"
#include "CbmMuchCluster.h"
#include "CbmMuchFindVectors.h"
#include "CbmMuchDigiMatch.h"
#include "CbmMuchFindHitsStraws.h"
#include "CbmMuchModule.h"
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
#include "CbmMuchStation.h"
#include "CbmMuchStrawHit.h"
#include "CbmMuchTrack.h"

#include "FairEventHeader.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include <TClonesArray.h>
#include <TF1.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TMath.h>
#include <TROOT.h>

#include <iostream>

using std::map;
using std::pair;
using std::set;
using std::multiset;
using std::cout;
using std::endl;

// -----   Default (stabdard) constructor   --------------------------------
CbmMuchFindVectorsQA::CbmMuchFindVectorsQA()
  : FairTask("MuchFindVectorsQA"),
    fGeoScheme(CbmMuchGeoScheme::Instance()),
    fStatFirst(-1),
    fNstat(0)
{
}
// -------------------------------------------------------------------------

// -----   Destructor   ----------------------------------------------------
CbmMuchFindVectorsQA::~CbmMuchFindVectorsQA()
{
}
// -------------------------------------------------------------------------

// -----   Public method Init (abstract in base class)  --------------------
InitStatus CbmMuchFindVectorsQA::Init()
{

  // Get and check FairRootManager
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman == NULL) Fatal("CbmMuchFindTracksQA::Init", "RootManager not instantiated!");

  //CbmMuchFindHitsStraws *hitFinder = (CbmMuchFindHitsStraws*) 
  //FairRun::Instance()->GetTask("CbmMuchFindHitsStraws");
  //if (hitFinder == NULL) Fatal("CbmMuchFindTracks::Init", "CbmMuchFindHitsStraws not run!");

  fVectors     = static_cast<TClonesArray*> (ioman->GetObject("MuchVector"));
  fMCTracks    = static_cast<TClonesArray*> (ioman->GetObject("MCTrack"));
  fPoints      = static_cast<TClonesArray*> (ioman->GetObject("MuchPoint"));
  fHits        = static_cast<TClonesArray*> (ioman->GetObject("MuchStrawHit"));
  fHitsGem     = static_cast<TClonesArray*> (ioman->GetObject("MuchPixelHit"));
  fDigis       = static_cast<TClonesArray*> (ioman->GetObject("MuchStrawDigi"));
  fDigisGem    = static_cast<TClonesArray*> (ioman->GetObject("MuchDigi"));
  fDigiMatches = static_cast<TClonesArray*> (ioman->GetObject("MuchStrawDigiMatch"));
  fDigiMatchesGem = static_cast<TClonesArray*> (ioman->GetObject("MuchDigiMatch"));
  fClusters    = static_cast<TClonesArray*> (ioman->GetObject("MuchCluster"));

  // Get straw system configuration
  Int_t nSt = fGeoScheme->GetNStations();

  for (Int_t i = 0; i < nSt; ++i) {
    CbmMuchStation* st = fGeoScheme->GetStation(i);
    CbmMuchModule* mod = fGeoScheme->GetModule(i, 0, 0, 0);
    Int_t stat = CbmMuchAddress::GetStationIndex(st->GetDetectorId()); 
    if (mod->GetDetectorType() == 2 && fStatFirst < 0) fStatFirst = stat;
    ++fNstat;
    Int_t nLays = st->GetNLayers();
    fNdoubl[stat] = nLays;

    for (Int_t lay = 0; lay < nLays; ++lay) {
      CbmMuchLayer* layer = st->GetLayer(lay);
      
      for (Int_t iside = 0; iside < 2; ++iside) {
	CbmMuchLayerSide* side = layer->GetSide(iside);
	Int_t plane = lay * 2 + iside;
	fZpos[stat][plane] = side->GetZ();
      }
    }
  }
  if (fStatFirst < 0) fStatFirst = 99; // all GEM stations
 
  // Create directory for histograms
  TDirectory *dir = new TDirectory("muchQA","","");
  dir->cd();

  BookHistos();

  return kSUCCESS;
}
// -------------------------------------------------------------------------

// -----  SetParContainers -------------------------------------------------
void CbmMuchFindVectorsQA::SetParContainers()
{
    FairRunAna* ana = FairRunAna::Instance();
    FairRuntimeDb* rtdb = ana->GetRuntimeDb();

    rtdb->getContainer("FairBaseParSet");
    rtdb->getContainer("CbmGeoPassivePar");
    //rtdb->getContainer("CbmGeoStsPar");
    //rtdb->getContainer("CbmGeoRichPar");
    rtdb->getContainer("CbmGeoMuchPar");
    //rtdb->getContainer("CbmGeoTrdPar");
    rtdb->getContainer("CbmFieldPar");
}
// -------------------------------------------------------------------------

// -----   Private method BookHistos   -------------------------------------
void CbmMuchFindVectorsQA::BookHistos()
{
  // Book histograms

  fhNvec      = new TH1D*[fNstat];
  fhNdoubl    = new TH1D*[fNstat];
  fhNhits     = new TH1D*[fNstat];
  fhNhitsOk   = new TH1D*[fNstat];
  fhChi2      = new TH1D*[fNstat];
  fhNgood     = new TH1D*[fNstat];
  fhNghost    = new TH1D*[fNstat];
  fhChi2ok    = new TH1D*[fNstat];
  fhChi2bad   = new TH1D*[fNstat];
  fhDx        = new TH1D*[fNstat];
  fhDy        = new TH1D*[fNstat];
  fhDtx       = new TH1D*[fNstat];
  fhDty       = new TH1D*[fNstat];
  fhIds       = new TH1D*[fNstat];
  fhIdVsEv    = new TH2D*[fNstat];
  fhDtxAll    = new TH1D*[fNstat];
  fhDtyAll    = new TH1D*[fNstat];
  fhDtxOk     = new TH1D*[fNstat];
  fhDtyOk     = new TH1D*[fNstat];
  fhShort     = new TH1D*[fNstat];
  fhOverlap   = new TH2D*[fNstat];
  fhChi2mat   = new TH1D*[fNstat];
  fhMatchMult = new TH1D*[fNstat];
  fhOccup     = new TH1D*[fNstat];
  //fhSim     = new TH1D*[fNstat];
  //fhRec     = new TH1D*[fNstat];
  fhDx12    = new TH2D*[fNstat];
  fhDx23    = new TH2D*[fNstat];
  fhDy12    = new TH2D*[fNstat];
  fhDy23    = new TH2D*[fNstat];

  for (Int_t ist = 0; ist < fNstat; ++ist) {
    Int_t stat = ist;
    fhNvec[ist] = new TH1D(Form("hNvec%i",stat),Form("Number of vectors in station %i",stat),200,0,1000);
    fhNdoubl[ist] = new TH1D(Form("hNdoubl%i",stat),Form("Number of tubes in lay. 1 in station %i",stat),200,0,200);
    fhNhits[ist] = new TH1D(Form("hNhits%i",stat),Form("Number of hits/vector in station %i",stat),20,0,20);
    fhNhitsOk[ist] = new TH1D(Form("hNhitsOk%i",stat),Form("Number of hits/good vector in station %i",stat),20,0,20);
    fhChi2[ist] = new TH1D(Form("hChi2all%i",stat),Form("Chi2 of all vectors in station %i",stat),200,0,200);
    fhNgood[ist] = new TH1D(Form("hNgood%i",stat),Form("Number of good vectors in station %i",stat),100,0,100);
    fhNghost[ist] = new TH1D(Form("hNghost%i",stat),Form("Number of ghost vectors in station %i",stat),200,0,200);
    fhChi2ok[ist] = new TH1D(Form("hChi2ok%i",stat),Form("Chi2 of good vectors in station %i",stat),200,0,200);
    fhChi2bad[ist] = new TH1D(Form("hChi2bad%i",stat),Form("Chi2 of ghost vectors in station %i",stat),200,0,200);
    fhDx[ist] = new TH1D(Form("hDx%i",stat),Form("Xrec-Xmc in station %i",stat),100,-2,2);
    fhDy[ist] = new TH1D(Form("hDy%i",stat),Form("Yrec-Ymc in station %i",stat),100,-15,15);
    fhDtx[ist] = new TH1D(Form("hDtx%i",stat),Form("TXrec-TXmc in station %i",stat),100,-0.2,0.2);
    fhDty[ist] = new TH1D(Form("hDty%i",stat),Form("TYrec-TYmc in station %i",stat),100,-1,1);
    fhIds[ist] = new TH1D(Form("hIds%i",stat),Form("Good track ID in station %i",stat),100,0,100);
    fhIdVsEv[ist] = new TH2D(Form("hIdVsEv%i",stat),Form("Good track ID vs ev. No. in station %i",stat),
			     100,-1,99,20,-1,19);
    fhIdVsEv[ist]->SetOption("box");
    fhDtxAll[ist] = new TH1D(Form("hDtxAll%i",stat),Form("TXvec-TXhit in station %i",stat),100,-2,2);
    fhDtyAll[ist] = new TH1D(Form("hDtyAll%i",stat),Form("TYvec-TYhit in station %i",stat),100,-2,2);
    fhDtxOk[ist] = new TH1D(Form("hDtxOk%i",stat),Form("TXvec-TXhit in station %i (ok)",stat),100,-2,2);
    fhDtyOk[ist] = new TH1D(Form("hDtyOk%i",stat),Form("TYvec-TYhit in station %i (ok)",stat),100,-2,2);
    fhDtube[ist][0] = new TH1D(Form("hDtube21_%i",stat),Form("Dtube 21 vs tube in station %i",stat),100,-50,50);
    fhDtube[ist][1] = new TH1D(Form("hDtube32_%i",stat),Form("Dtube 32 vs tube in station %i",stat),100,-100,100);
    fhDtube[ist][2] = new TH1D(Form("hDtube43_%i",stat),Form("Dtube 43 vs tube in station %i",stat),100,-200,200);
    fhShort[ist] = new TH1D(Form("hShort%i",stat),Form("Number of shared hits in station %i",stat),20,0,20);
    fhOverlap[ist] = new TH2D(Form("hOverlap%i",stat),Form("Overlap multipl. in station %i",stat),20,0,20,20,0,20);
    fhChi2mat[ist] = new TH1D(Form("hChi2mat%i",stat),Form("Chi2 of matching in station %i",stat),100,0,50);
    fhMatchMult[ist] = new TH1D(Form("hMatchMult%i",stat),Form("Multiplicity of matching in station %i",stat),10,0,10);
    fhOccup[ist] = new TH1D(Form("hOccup%i",stat),Form("Occupancy in station %i",stat),600,-300,300);
    fhDtube2[ist][0] = new TH2D(Form("hDtube31_%i",stat),Form("Dtube 31 vs tube in station %i",stat),100,-300,300,100,-200,200);
    fhDtube2[ist][1] = new TH2D(Form("hDtube42_%i",stat),Form("Dtube 42 vs tube in station %i",stat),100,-300,300,100,-100,100);
    fhMCFit[ist][0] = new TH1D(Form("hMCFitX_%i",stat),Form("Chi2 of fit in X in station %i",stat),100,0,50);
    fhMCFit[ist][1] = new TH1D(Form("hMCFitY_%i",stat),Form("Chi2 of fit in Y in station %i",stat),100,0,50);
    fhDx12[ist] = new TH2D(Form("hDx12_%i",stat),Form("Dx2 vs Dx1 in station %i",stat),50,-10,10,50,-10,10);
    fhDx23[ist] = new TH2D(Form("hDx23_%i",stat),Form("Dx3 vs Dx2 in station %i",stat),50,-10,10,50,-10,10);
    fhDy12[ist] = new TH2D(Form("hDy12_%i",stat),Form("Dy2 vs Dy1 in station %i",stat),50,-10,10,50,-10,10);
    fhDy23[ist] = new TH2D(Form("hDy23_%i",stat),Form("Dy3 vs Dy2 in station %i",stat),50,-10,10,50,-10,10);
  }
  fhSim = new TH1D("hSim","Number of reconstructable muons",10,0,10);
  fhRec = new TH1D("hRec","Number of reconstructed muons",10,0,10);
  fhZXY[0] = new TH1D("hZX","Z - X",180,0,45);
  fhZXY[1] = new TH1D("hZY","Z - Y",180,0,45);
  fhEvents = new TH1D("hEvents","Number of processed events",5,0,5);
}
// -------------------------------------------------------------------------

// -----   Public method Exec   --------------------------------------------
void CbmMuchFindVectorsQA::Exec(Option_t* opt)
{
  // Do all processing

  Int_t *mult = new Int_t [fNstat];
  Int_t *multOk = new Int_t [fNstat];
  Int_t *multBad = new Int_t [fNstat];
  for (Int_t i = 0; i < fNstat; ++i) mult[i] = multOk[i] = multBad[i] = 0;
  fhEvents->Fill(1);

  Int_t nvec = fVectors->GetEntriesFast();
  cout << " nvec " << nvec << endl;
  for (Int_t iv = 0; iv < nvec; ++iv) {
    CbmMuchTrack *vec = (CbmMuchTrack*) fVectors->UncheckedAt(iv);
    Int_t ista = vec->GetUniqueID();
    const FairTrackParam *params = vec->GetParamFirst();
    ++mult[ista];
    fhChi2[ista]->Fill(vec->GetChiSq());
    fhNhits[ista]->Fill(vec->GetNofHits());
    fhDtxAll[ista]->Fill(params->GetTx()-params->GetX()/params->GetZ());
    fhDtyAll[ista]->Fill(params->GetTy()-params->GetY()/params->GetZ());
    if (CheckMatch(vec)) {
      // "Good" vector
      ++multOk[ista];
      fhChi2ok[ista]->Fill(vec->GetChiSq());
    } else {
      // Ghost
      ++multBad[ista];
      fhChi2bad[ista]->Fill(vec->GetChiSq());
    }
  }
  cout << mult[0] << " " << mult[1] << endl;
  
  for (Int_t i = 0; i < fNstat; ++i) {
    fhNvec[i]->Fill(mult[i]);
    fhNgood[i]->Fill(multOk[i]);
    fhNghost[i]->Fill(multBad[i]);
  }

  // Check short tracks
  CheckShorts(fHits);

  // Check efficiency
  CheckEffic();

  // For straws: fill occupancy histos as well (for completeness)
  // and number of hits in the first layer of the first doublet
  Int_t nhits = fHits->GetEntriesFast(), nDouble[10] = {0};
  for (Int_t i = 0; i < nhits; ++i) {
    CbmMuchStrawHit *hit = (CbmMuchStrawHit*) fHits->UncheckedAt(i);
    Int_t ista = fGeoScheme->GetStationIndex(hit->GetAddress());
    fhOccup[ista]->Fill(hit->GetTube());
    Int_t lay = fGeoScheme->GetLayerIndex(hit->GetAddress());
    Int_t side = fGeoScheme->GetLayerSideIndex(hit->GetAddress());
    if (lay == 0 && side == 0) ++nDouble[ista];
  }
  //for (Int_t ista = 0; ista < 2; ++ista) fhNdoubl[ista]->Fill(nDouble[ista]);
  for (Int_t ista = 0; ista < fNstat; ++ista) fhNdoubl[ista]->Fill(nDouble[ista]);

  delete [] mult;
  delete [] multOk;
  delete [] multBad;
}
// -------------------------------------------------------------------------

// -----   Private method CheckMatch   -------------------------------------
Bool_t CbmMuchFindVectorsQA::CheckMatch(CbmMuchTrack *vec)
{
  // Check matching quality of the vector

  Int_t ista = vec->GetUniqueID();
  if (ista < fStatFirst) return CheckMatchGem(vec);

  Int_t nhits = vec->GetNofHits(), nthr = nhits / 2, id = 0;
  //Int_t nhits = vec->GetNofHits(), nthr = nhits - 1, id = 0;
  map<Int_t,Int_t> ids;

  for (Int_t ih = 0; ih < nhits; ++ih) {
    CbmMuchStrawHit *hit = (CbmMuchStrawHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
    //if (hit->GetFlag() > 1) { cout << " hit " << hit->GetFlag() << endl; /*exit(0);*/ }
    if (hit->GetFlag() % 2) {
    //if (0) {
      // Mirror hit
      id = -1; 
      if (ids.find(id) == ids.end()) ids.insert(pair<Int_t,Int_t>(id,1));
      else ++ids[id];
    } else {
      CbmMuchDigiMatch* digiM = (CbmMuchDigiMatch*) fDigiMatches->UncheckedAt(hit->GetRefId());
      Int_t np = digiM->GetNofLinks();
      for (Int_t ip = 0; ip < np; ++ip) {
	//CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(0).GetIndex());
	CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
	id = point->GetTrackID();
	//if (np > 1) cout << ip << " " << id << endl;
	if (ids.find(id) == ids.end()) ids.insert(pair<Int_t,Int_t>(id,1));
	else ++ids[id];
      }
      //if (np > 1) { cout << " digi " << digiM->GetNofLinks() << " " << ista << " " << hit->GetX() << " " << hit->GetY() << endl; exit(0); }
    }
  }
  Int_t maxim = -1, idmax = -9;
  for (map<Int_t,Int_t>::iterator it = ids.begin(); it != ids.end(); ++it) {
    if (it->second > maxim) {
      maxim = it->second;
      idmax = it->first;
    }
  }
  // Set vector ID as its flag
  //vec->SetFlag(idmax);
  if (maxim <= nthr || idmax < 0) return kFALSE;

  if (idmax > 1) return kFALSE; // !!! look only at muons from Pluto

  fhNhitsOk[ista]->Fill(vec->GetNofHits());
  fhIds[ista]->Fill(idmax);
  Int_t evID = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber() + 1;
  fhIdVsEv[ista]->Fill(evID,idmax);

  // Fill true track parameters
  Double_t zp = 0.0, parsTr[4] = {0.0};
  Int_t ok = 0;
  for (Int_t ih = 0; ih < nhits; ++ih) {
    CbmMuchStrawHit *hit = (CbmMuchStrawHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
    CbmMuchDigiMatch* digiM = (CbmMuchDigiMatch*) fDigiMatches->UncheckedAt(hit->GetRefId());
    Int_t np = digiM->GetNofLinks();
    for (Int_t ip = 0; ip < np; ++ip) {
      CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
      id = point->GetTrackID();
      if (id != idmax) continue;
      parsTr[0] = (point->GetXIn() + point->GetXOut()) / 2;
      parsTr[1] = (point->GetYIn() + point->GetYOut()) / 2;
      zp = (point->GetZIn() + point->GetZOut()) / 2;
      parsTr[2] = point->GetPx() / point->GetPz();
      parsTr[3] = point->GetPy() / point->GetPz();
      ok = 1;
      break;
    }
    if (ok) break;
  }
  // Adjust track coordinates
  const FairTrackParam *params = vec->GetParamFirst();
  parsTr[0] += parsTr[2] * (params->GetZ() - zp);
  parsTr[1] += parsTr[3] * (params->GetZ() - zp);
  fhDx[ista]->Fill(params->GetX()-parsTr[0]);
  fhDy[ista]->Fill(params->GetY()-parsTr[1]);
  fhDtx[ista]->Fill(params->GetTx()-parsTr[2]);
  fhDty[ista]->Fill(params->GetTy()-parsTr[3]);
  fhDtxOk[ista]->Fill(params->GetTx()-params->GetX()/params->GetZ());
  fhDtyOk[ista]->Fill(params->GetTy()-params->GetY()/params->GetZ());
  cout << " Good: " << idmax << " " << zp << " " << params->GetZ() << " " << parsTr[0] << " " << params->GetX() << " " << parsTr[1] << " " << params->GetY() << " " << parsTr[2] << " " << params->GetTx() << " " << parsTr[3] << " " << params->GetTy() << " " << vec->GetChiSq() << endl;

  // Fill tube difference between doublets
  map<Int_t,Int_t> tubes[2];
  Int_t id0 = -1;
  // Check configuration
  Int_t combi = 0;
  if (fNdoubl[fNstat-1] == 3) combi = 3; // combined stations for 3-layer config.

  for (Int_t ih = 0; ih < nhits; ++ih) {
    CbmMuchStrawHit *hit = (CbmMuchStrawHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
    CbmMuchDigiMatch* digiM = (CbmMuchDigiMatch*) fDigiMatches->UncheckedAt(hit->GetRefId());
    Int_t np = digiM->GetNofLinks();
    for (Int_t ip = 0; ip < np; ++ip) {
      CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
      id = point->GetTrackID();
      if (id > 1) continue; // look only at muons from vector mesons
      if (id0 < 0) id0 = id;
      break;
    }
    if (id0 < 0 || id != id0) continue;
    Int_t stat = fGeoScheme->GetStationIndex(hit->GetAddress()) - fStatFirst;
    Int_t lay = fGeoScheme->GetLayerIndex(hit->GetAddress());
    tubes[id].insert(pair<Int_t,Int_t>(lay+stat*combi,hit->GetTube()));
  }
  Double_t slope[2] = {0.08, 0.06};
  for (Int_t j = 0; j < 2; ++j) {
    // Loop over trackID
    Int_t nd = tubes[j].size();
    if (nd < 2) continue;
    Int_t tubeNos[6] = {1000, 1000, 1000, 1000, 1000, 1000};
    map<Int_t,Int_t>::iterator it = tubes[j].begin();
    tubeNos[it->first] = it->second;
    ++it;
    for ( ; it != tubes[j].end(); ++it) tubeNos[it->first] = it->second;
    for (Int_t lay = 1; lay < 4; ++lay) {
      if (tubeNos[lay] < 1000 & tubeNos[lay-1] < 1000) 
	fhDtube[ista][lay-1]->Fill(tubeNos[lay]-tubeNos[lay-1]);
      if (combi == 0) {
	// 4-layer config.
	if (lay == 2 && tubeNos[lay] < 1000 & tubeNos[lay-2] < 1000) 
	  fhDtube2[ista][lay-2]->Fill(tubeNos[lay-2],tubeNos[lay]-tubeNos[lay-2]);
	if (lay == 3 && tubeNos[lay] < 1000 & tubeNos[lay-2] < 1000) {
	  //fhDtube2[ista][lay-2]->Fill(tubeNos[lay-2],tubeNos[lay]-tubeNos[lay-2]-slope[ista]*tubeNos[lay-2]);
	  Double_t dzz = (fZpos[ista][lay*2]-fZpos[ista][lay*2-4]) / fZpos[ista][lay*2-4];
	  fhDtube2[ista][lay-2]->Fill(tubeNos[lay-2],tubeNos[lay]-tubeNos[lay-2]-dzz*tubeNos[lay-2]);
	}
      } else {
	// 3-layer config.
	if (lay == 3 && tubeNos[lay] < 1000 & tubeNos[lay-3] < 1000) 
	  fhDtube2[ista][lay-3]->Fill(tubeNos[lay-3],tubeNos[lay]-tubeNos[lay-3]);
	if (lay == 1 && tubeNos[lay] < 1000 & tubeNos[lay+3] < 1000) {
	  //fhDtube2[ista][lay-2]->Fill(tubeNos[lay-2],tubeNos[lay]-tubeNos[lay-2]-slope[ista]*tubeNos[lay-2]);
	  Double_t dzz = (fZpos[ista+1][lay*2]-fZpos[ista][lay*2]) / fZpos[ista][lay*2];
	  fhDtube2[ista][lay]->Fill(tubeNos[lay],tubeNos[lay+3]-tubeNos[lay]-dzz*tubeNos[lay]);
	}
      }
    }
    if (ok) break;
  }
  
  return kTRUE;
}
// -------------------------------------------------------------------------

// -----   Private method CheckMatchGem   ----------------------------------
Bool_t CbmMuchFindVectorsQA::CheckMatchGem(CbmMuchTrack *vec)
{
  // Check matching quality of the vector for GEM stations

  Int_t ista = vec->GetUniqueID();
  Int_t nhits = vec->GetNofHits(), nthr = nhits / 2, id = 0;
  //Int_t nhits = vec->GetNofHits(), nthr = nhits - 1, id = 0;
  map<Int_t,Int_t> ids;

  for (Int_t ih = 0; ih < nhits; ++ih) {
    CbmMuchPixelHit *hit = (CbmMuchPixelHit*) fHitsGem->UncheckedAt(vec->GetHitIndex(ih));
    CbmMuchCluster *clus = (CbmMuchCluster*) fClusters->UncheckedAt(hit->GetRefId());
    Int_t nDigis = clus->GetNofDigis();
    for (Int_t j = 0; j < nDigis; ++j) {
      CbmMuchDigiMatch* digiM = (CbmMuchDigiMatch*) fDigiMatchesGem->UncheckedAt(clus->GetDigi(j));
      Int_t np = digiM->GetNofLinks();
      for (Int_t ip = 0; ip < np; ++ip) {
	//CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(0).GetIndex());
	CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
	id = point->GetTrackID();
	//if (np > 1) cout << ip << " " << id << endl;
	ids[id]++;
      }
      //if (np > 1) { cout << " digi " << digiM->GetNofLinks() << " " << ista << " " << hit->GetX() << " " << hit->GetY() << endl; exit(0); }
    }
  }

  Int_t maxim = -1, idmax = -9;
  for (map<Int_t,Int_t>::iterator it = ids.begin(); it != ids.end(); ++it) {
    if (it->second > maxim) {
      maxim = it->second;
      idmax = it->first;
    }
  }
  //cout << " GEMS: " << ista << " " << nhits << " " << maxim << " " << idmax << endl;
  // Set vector ID as its flag
  //vec->SetFlag(idmax);
  if (maxim <= nthr || idmax < 0) return kFALSE;

  if (idmax > 1) return kFALSE; // !!! look only at muons from Pluto

  fhNhitsOk[ista]->Fill(vec->GetNofHits());
  fhIds[ista]->Fill(idmax);
  Int_t evID = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber() + 1;
  fhIdVsEv[ista]->Fill(evID,idmax);

  // Fill true track parameters
  Double_t zp = 0.0, parsTr[4] = {0.0};
  Int_t ok = 0;

  for (Int_t ih = 0; ih < nhits; ++ih) {
    CbmMuchPixelHit *hit = (CbmMuchPixelHit*) fHitsGem->UncheckedAt(vec->GetHitIndex(ih));
    CbmMuchCluster *clus = (CbmMuchCluster*) fClusters->UncheckedAt(hit->GetRefId());
    Int_t nDigis = clus->GetNofDigis();
    for (Int_t j = 0; j < nDigis; ++j) {
      CbmMuchDigiMatch* digiM = (CbmMuchDigiMatch*) fDigiMatchesGem->UncheckedAt(clus->GetDigi(j));
      Int_t np = digiM->GetNofLinks();
      for (Int_t ip = 0; ip < np; ++ip) {
	CbmMuchPoint* point = (CbmMuchPoint*) fPoints->UncheckedAt(digiM->GetLink(ip).GetIndex());
	id = point->GetTrackID();
	if (id != idmax) continue;
	parsTr[0] = (point->GetXIn() + point->GetXOut()) / 2;
	parsTr[1] = (point->GetYIn() + point->GetYOut()) / 2;
	zp = (point->GetZIn() + point->GetZOut()) / 2;
	parsTr[2] = point->GetPx() / point->GetPz();
	parsTr[3] = point->GetPy() / point->GetPz();
	ok = 1;
	break;
      }
      if (ok) break;
    }
    if (ok) break;
  }
  // Adjust track coordinates
  const FairTrackParam *params = vec->GetParamFirst();
  parsTr[0] += parsTr[2] * (params->GetZ() - zp);
  parsTr[1] += parsTr[3] * (params->GetZ() - zp);
  fhDx[ista]->Fill(params->GetX()-parsTr[0]);
  fhDy[ista]->Fill(params->GetY()-parsTr[1]);
  fhDtx[ista]->Fill(params->GetTx()-parsTr[2]);
  fhDty[ista]->Fill(params->GetTy()-parsTr[3]);
  fhDtxOk[ista]->Fill(params->GetTx()-params->GetX()/params->GetZ());
  fhDtyOk[ista]->Fill(params->GetTy()-params->GetY()/params->GetZ());
  cout << " Good: " << idmax << " " << zp << " " << params->GetZ() << " " << parsTr[0] << " " << params->GetX() << " " << parsTr[1] << " " << params->GetY() << " " << parsTr[2] << " " << params->GetTx() << " " << parsTr[3] << " " << params->GetTy() << " " << vec->GetChiSq() << endl;

  return kTRUE;
}
// -------------------------------------------------------------------------

// -----   Private method CheckShorts   ------------------------------------
void CbmMuchFindVectorsQA::CheckShorts(TClonesArray *hitArray)
{
  // Check short tracks

  Int_t nvec = fVectors->GetEntriesFast();
  Int_t nhits0 = 4;
  for (Int_t iv = 0; iv < nvec; ++iv) {
    CbmMuchTrack *vec = (CbmMuchTrack*) fVectors->UncheckedAt(iv);
    Int_t nhits = vec->GetNofHits();
    if (nhits != nhits0) continue;
    Int_t ista = vec->GetUniqueID();
    set<Int_t> overlap;
    multiset<Int_t> overlap1;

    for (Int_t iv1 = iv + 1; iv1 < nvec; ++iv1) {
      CbmMuchTrack *vec1 = (CbmMuchTrack*) fVectors->UncheckedAt(iv1);
      if (vec1->GetUniqueID() != ista) continue;
      Int_t nhits1 = vec1->GetNofHits();

      // Compare hits
      Int_t ih1 = 0;
      for (Int_t ih = 0; ih < nhits; ++ih) {
	//CbmMuchStrawHit *hit = (CbmMuchStrawHit*) hitArray->UncheckedAt(vec->GetHitIndex(ih));
	CbmHit *hit = (CbmHit*) hitArray->UncheckedAt(vec->GetHitIndex(ih));
	Int_t lay = fGeoScheme->GetLayerIndex(hit->GetAddress());
	Int_t side = fGeoScheme->GetLayerSideIndex(hit->GetAddress());
	Int_t plane = lay * 2 + side;

	for ( ; ih1 < nhits1; ++ih1) { 
	  //CbmMuchStrawHit *hit1 = (CbmMuchStrawHit*) hitArray->UncheckedAt(vec1->GetHitIndex(ih1));
	  CbmHit *hit1 = (CbmHit*) hitArray->UncheckedAt(vec1->GetHitIndex(ih1));
	  Int_t lay1 = fGeoScheme->GetLayerIndex(hit1->GetAddress());
	  Int_t side1 = fGeoScheme->GetLayerSideIndex(hit1->GetAddress());
	  Int_t plane1 = lay1 * 2 + side1;
	  if (plane1 < plane) continue;
	  if (plane1 > plane) break;
	  if (hit == hit1) { overlap.insert(ih); overlap1.insert(ih); }
	}
      }
      if (overlap.size() == nhits0) {
	// All hits are shared with other tracks
	//break;
      }
    } // for (Int_t iv1 = iv + 1;
    fhShort[ista]->Fill(overlap.size());
    if (overlap.size()) fhOverlap[ista]->Fill(overlap.size(),overlap1.size());
  } // for (Int_t iv = 0; iv < nvec;
}
// -------------------------------------------------------------------------

// -----   Private method CheckEffic   -------------------------------------
void CbmMuchFindVectorsQA::CheckEffic()
{
  // Check efficiency of the vector reco

  // Find "reconstructable" vectors - having points in all doublets (GEM or straw)
  // !!! Take only trackID = 0, 1 (from signal mesons)
  const Int_t nMu = 2;
  set<Int_t> doublets[20][nMu], singlets[20][nMu];
  Int_t nPoints = fPoints->GetEntriesFast();
  Double_t xp[7][20][nMu], yp[7][20][nMu];
  Bool_t lstraws = kFALSE;

  for (Int_t ip = 0; ip < nPoints; ++ip) {
    CbmMuchPoint *p = (CbmMuchPoint*) fPoints->UncheckedAt(ip);
    Int_t id = p->GetTrackID();
    if (id > 1) continue;
    if (p->GetZ() < fZpos[0][0] - 2.0) continue;
    Int_t ista = fGeoScheme->GetStationIndex(p->GetDetectorId());
    Int_t lay = fGeoScheme->GetLayerIndex(p->GetDetectorId());
    Int_t side = fGeoScheme->GetLayerSideIndex(p->GetDetectorId());

    doublets[ista][id].insert(lay);
    singlets[ista][id].insert(lay*2+side);
    xp[ista][lay*2+side][id] = (p->GetXIn() + p->GetXOut()) / 2;
    yp[ista][lay*2+side][id] = (p->GetYIn() + p->GetYOut()) / 2;

    if (!lstraws) {
      if (fGeoScheme->GetModuleByDetId(p->GetDetectorId())->GetDetectorType() == 2) lstraws = kTRUE;
    }
  }

  Int_t muons[7][nMu] = {{0},{0}}, nMuVec = 0;
  for (Int_t ista = 0; ista < fNstat; ++ista) {
    if (doublets[ista][0].size() == fNdoubl[ista]) { fhSim->Fill(ista); muons[ista][0] = 1; ++nMuVec; }
    if (doublets[ista][1].size() == fNdoubl[ista]) { fhSim->Fill(ista); muons[ista][1] = 1; ++nMuVec; }
  }
  //if (muons[0][0] == 0 && muons[0][1] == 0 && muons[1][0] == 0 && muons[1][1] == 0) return; 
  if (nMuVec == 0) return;

  // Fit to straight lines in 2 projections
  /*
  for (Int_t ista = 0; ista < fNstat; ++ista) {

    for (Int_t mu = 0; mu < 2; ++ mu) {
      if (muons[ista][mu] == 0) continue;

      fhZXY[0]->Reset();
      fhZXY[1]->Reset();
      for (Int_t plane = 0; plane < nLays2; ++plane) {
	if (singlets[ista][mu].find(plane) == singlets[ista][mu].end()) continue;
	fhZXY[0]->Fill(zPos[ista][plane]-zPos[ista][0]+5, xp[ista][plane][mu]);
	Int_t ib = fhZXY[1]->Fill(zPos[ista][plane]-zPos[ista][0]+5, yp[ista][plane][mu]);
	fhZXY[0]->SetBinError(ib, 0.1);
	fhZXY[1]->SetBinError(ib, 0.1);
      }
      for (Int_t jxy = 0; jxy < 2; ++jxy) {
	fhZXY[jxy]->Fit("pol1","Q0");
	TF1 *f = fhZXY[jxy]->GetFunction("pol1");
	fhMCFit[ista][jxy]->Fill(f->GetChisquare() / f->GetNDF());
      }
    }
  }
  */
  
  // Match reconstructed vectors
  Int_t nvec = fVectors->GetEntriesFast();
  Int_t pluto = 1; // !!! set to 1 to check exact matching for Pluto 
  //Double_t errx = 0.22, erry = 1.41, chi2cut = 10.0; // 10 deg
  //Double_t errx = 0.2, erry = 1.9, chi2cut = 10.0; // 5 doublets
  Double_t errxS[6] = {0.33, 0.50, 0.2, 0.2, 1, 1};
  Double_t erryS[6] = {0.33, 0.50, 1.9, 1.9, 1, 1}, chi2cut = 10.0; // 2 GEMS, 2 straws
  //Double_t errxS[6] = {0.33, 0.50, 0.03, 0.03, 1, 1};
  //Double_t erryS[6] = {0.33, 0.50, 0.26, 0.26, 1, 1}, chi2cut = 10.0; // 2 GEMS, 2 straws
  Double_t errxG[6] = {0.33, 0.54, 1.0, 1.7, 1, 1};
  Double_t erryG[6] = {0.33, 0.54, 1.0, 1.7, 1, 1}; // 4 GEMS
  Int_t nMatch[7][nMu] = {{0},{0}}, combi = 0;
  Double_t *errx = errxG, *erry = erryG;
  if (lstraws) { 
    errx = errxS; 
    erry = erryS; 
    if (fNdoubl[fNstat-1] == 3) combi = 6; // combined stations for 3-layer config.
  }

  for (Int_t iv = 0; iv < nvec; ++iv) {
    CbmMuchTrack *vec = (CbmMuchTrack*) fVectors->UncheckedAt(iv);
    Int_t nhits = vec->GetNofHits();
    Int_t ista = vec->GetUniqueID();
    if (muons[ista][0] == 0 && muons[ista][1] == 0) continue;
    Int_t id = vec->GetFlag();
    if (pluto && id > 1) continue; // !!! this is only for Pluto sample - exact ID match!!!
    const FairTrackParam *params = vec->GetParamFirst();
    TClonesArray *hitArray = (ista < fStatFirst) ? fHitsGem : fHits;

    Double_t chi2[nMu] = {0.0};
    for (Int_t mu = 0; mu < nMu; ++ mu) {
      if (!muons[ista][mu]) { chi2[mu] = 999999.; continue; }
      Int_t nm = 0;
      Double_t dx0 = 0, dy0 = 0;

      for (Int_t ih = 0; ih < nhits; ++ih) {
	//CbmMuchStrawHit *hit = (CbmMuchStrawHit*) fHits->UncheckedAt(vec->GetHitIndex(ih));
	CbmHit *hit = (CbmHit*) hitArray->UncheckedAt(vec->GetHitIndex(ih));
	Int_t stat = fGeoScheme->GetStationIndex(hit->GetAddress()) - fStatFirst;
	Int_t lay = fGeoScheme->GetLayerIndex(hit->GetAddress());
	Int_t side = fGeoScheme->GetLayerSideIndex(hit->GetAddress());
	Int_t plane = lay * 2 + side;
	plane += stat * combi;
	if (singlets[ista][mu].find(plane) == singlets[ista][mu].end()) continue;
	Double_t x = params->GetX() + params->GetTx() * (fZpos[ista][plane] - params->GetZ());
	Double_t y = params->GetY() + params->GetTy() * (fZpos[ista][plane] - params->GetZ());
	Double_t dx = x - xp[ista][plane][mu];
	Double_t dy = y - yp[ista][plane][mu];
	chi2[mu] += dx * dx / errx[ista] / errx[ista];
	chi2[mu] += dy * dy / erry[ista] / erry[ista];
	++nm;
	if (lay) {
	  if (lay == 1) {
	    fhDx12[ista]->Fill(dx0,dx);
	    fhDy12[ista]->Fill(dy0,dy);
	  } else if (lay == 2) {
	    fhDx23[ista]->Fill(dx0,dx);
	    fhDy23[ista]->Fill(dy0,dy);
	  }
	}
	dx0 = dx;
	dy0 = dy;
	//cout << nm << " " << dx << " " << dy << endl;
      }
      nm *= 2;
      if (nm > 4) chi2[mu] /= (nm - 4); // per NDF
    }
    if (pluto) fhChi2mat[ista]->Fill(chi2[id]); // !!! this is only for Pluto sample - exact ID match!!!
    else fhChi2mat[ista]->Fill(TMath::Min(chi2[0],chi2[1])); 
    if (chi2[0] < chi2cut) ++nMatch[ista][0];
    if (chi2[1] < chi2cut) ++nMatch[ista][1];
    cout << chi2[0] << " " << chi2[1] << endl;
  } // for (Int_t iv = 0; iv < nvec;
  for (Int_t ista = 0; ista < fNstat; ++ista) {
    if (combi && ista == fNstat - 1) continue; // 2 combined stations
    for (Int_t id = 0; id < nMu; ++id) {
      if (muons[ista][id]) {
	if (combi && ista == fNstat - 2 && muons[ista+1][id] == 0) continue; // 2 combined stations
	fhMatchMult[ista]->Fill(nMatch[ista][id]);
	if (nMatch[ista][id] == 0) cout << " !!! No match found !!! " << ista << " " << id << endl;
      }
    }
  }
}
// -------------------------------------------------------------------------

// -----   Public method Finish   ------------------------------------------
void CbmMuchFindVectorsQA::Finish()
{
  for (Int_t ista = 0; ista < fNstat; ++ista) {
    if (fNdoubl[ista]) fhOccup[ista]->Scale(1./2/fNdoubl[ista]/2/fhEvents->GetEntries());
  }

  TDirectory *dir = (TDirectory*) gROOT->FindObjectAny("muchQA");
  gDirectory->mkdir("muchQA");
  gDirectory->cd("muchQA");
  dir->GetList()->Write();
}
// -------------------------------------------------------------------------

// -------------------------------------------------------------------------

ClassImp(CbmMuchFindVectorsQA);