#include "CbmRichID.h"
#include "CbmRootManager.h"
#include "CbmMCTrack.h"
#include "CbmRichPoint.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "CbmRichLightSpot.h"
#include "iostream.h"
#include "TParticle.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TArrayD.h"
#include "TGraph.h"
#include "TVector3.h"
#include "TLorentzVector.h"

//------------------------------------------------------------------------------
CbmRichID::CbmRichID() :CbmTask()
{
}

//------------------------------------------------------------------------------
CbmRichID::CbmRichID(const char *name, const char *title)
  :CbmTask(name)
{

  fEvent = 0;
  fDebug = "";
  const Int_t maxHit = 5000;
  fNHits = 0;
  fXHits = new TArrayD(maxHit);
  fYHits = new TArrayD(maxHit);
  fNrgs = 0;
  fNminSTS = 6;        // minimum STS hits (for consistency take same number as for TrackExtrapolation!)
  fNminHits = 15;     // minimum number of RICH hits required for acceptance (e)
  fDistanceCut = 1.;  // maximum distance between ring center and matched track in cm
  fDiffRadiusMin = 2.;  // maximum difference between radii of LightSpot/Found Ring
  fPotentialHitsMin = 0.25;  // minimum fraction of potential hits
  fPotentialHitsMax = 3.5;  // maximum fraction of potential hits
  fDistRingID = 1.2; // maximum distance between light spot and ring

  fRmin = 5.;  // min radius for electrons -> fake ring ID
  fRmax = 6.3; // max radius for electrons

  fvx = 0.5;   // (x,y,z) cut for main-vertex tracks in cm
  fvy = 0.5;
  fvz = 0.1;


  // define different sets:
//----------------------------------------------------------
// (a) relaxed matching cuts + high ring quality:
//  fNminSTS = 6;
//  fNminHits = 20;
//  fDistanceCut = 3.;
//  fDiffRadiusMin = 3.;
//  fPotentialHitsMin = 0.5;
//  fPotentialHitsMax = 1.5;
//  fDistRingID = 2.;
//----------------------------------------------------------
// (b) stronger analysis cuts (-> high purity) + high ring quality:
//  fNminSTS = 6;
//  fNminHits = 20;
//  fDistanceCut = 1.;
//  fDiffRadiusMin = 2.;
//  fPotentialHitsMin = 0.7;
//  fPotentialHitsMax = 1.3;
//  fDistRingID = 1.2;


//  Double_t pi = TMath::Pi();

  mother_testID=-111;

  // generated particles

  fh_gen_e = new TH2D("fh_gen_e","generated e (y,pt)",50,0,5,40,0,2);
  fh_gen_pi = new TH2D("fh_gen_pi","generated pi (y,pt)",50,0,5,40,0,2);
  fh_gen_mu = new TH2D("fh_gen_mu","generated mu (y,pt)",50,0,5,40,0,2);
  fh_gen_k = new TH2D("fh_gen_k","generated K (y,pt)",50,0,5,40,0,2);

  // measured particles: minimum number of STS points required!

  fh_meas_e = new TH2D("fh_meas_e","measured e (y,pt)",50,0,5,40,0,2);
  fh_meas_pi = new TH2D("fh_meas_pi","measured pi (y,pt)",50,0,5,40,0,2);
  fh_meas_mu = new TH2D("fh_meas_mu","measured mu (y,pt)",50,0,5,40,0,2);
  fh_meas_k = new TH2D("fh_meas_k","measured K (y,pt)",50,0,5,40,0,2);

  // accepted particles (N_sts >= fNminSTS)

  fh_acc_e = new TH2D("fh_acc_e","accepted e (y,pt)",50,0,5,40,0,2);
  fh_acc_pi = new TH2D("fh_acc_pi","accepted pi (y,pt)",50,0,5,40,0,2);
  fh_acc_mu = new TH2D("fh_acc_mu","accepted mu (y,pt)",50,0,5,40,0,2);
  fh_acc_k = new TH2D("fh_acc_k","accepted K (y,pt)",50,0,5,40,0,2);

  //----------------------------------------------------------------------------------------
  // primary particles with reconstructed ring (no N_STS cut!)

  fh_ring_e = new TH2D("fh_ring_e","prim reconstructed e (y,pt)",50,0,5,40,0,2);
  fh_ring_pi = new TH2D("fh_ring_pi","prim reconstructed pi (y,pt)",50,0,5,40,0,2);
  fh_ring_mu = new TH2D("fh_ring_mu","prim reconstructed mu (y,pt)",50,0,5,40,0,2);
  fh_ring_k = new TH2D("fh_ring_k","prim reconstructed K (y,pt)",50,0,5,40,0,2);

  // ... with N_STS cut
  fh_ring_nsts_e = new TH2D("fh_ring_nsts_e","prim+STS reconstructed e (y,pt)",50,0,5,40,0,2);
  fh_ring_nsts_pi = new TH2D("fh_ring_nsts_pi","prim+STS reconstructed pi (y,pt)",50,0,5,40,0,2);
  fh_ring_nsts_mu = new TH2D("fh_ring_nsts_mu","prim+STS reconstructed mu (y,pt)",50,0,5,40,0,2);

  // main vertex particles with good rings
  fh_ring_mv_e = new TH2D("fh_ring_mv_e","prim+STS reconstructed e (y,pt)",50,0,5,40,0,2);
  fh_ring_mv_pi = new TH2D("fh_ring_mv_pi","prim+STS reconstructed pi (y,pt)",50,0,5,40,0,2);

  // all particles with reconstructed ring

  fh_ringall_e = new TH2D("fh_ringall_e","all reconstructed e (y,pt)",50,0,5,40,0,2);
  fh_ringall_pi = new TH2D("fh_ringall_pi","all reconstructed pi (y,pt)",50,0,5,40,0,2);
  fh_ringall_mu = new TH2D("fh_ringall_mu","all reconstructed mu (y,pt)",50,0,5,40,0,2);

  // all particles with light spot

  fh_spotall_e = new TH2D("fh_spotall_e","all spots e (y,pt)",50,0,5,40,0,2);
  fh_spotall_pi = new TH2D("fh_spotall_pi","all spots pi (y,pt)",50,0,5,40,0,2);
  fh_spotall_mu = new TH2D("fh_spotall_mu","all spots mu (y,pt)",50,0,5,40,0,2);

  // primary particles with light spot

  fh_spot_e = new TH2D("fh_spot_e","prim spot e (y,pt)",50,0,5,40,0,2);
  fh_spot_pi = new TH2D("fh_spot_pi","prim spot pi (y,pt)",50,0,5,40,0,2);
  fh_spot_mu = new TH2D("fh_spot_mu","prim spot mu (y,pt)",50,0,5,40,0,2);

  // (x,y) of not reconstructed spots

  fh_spot_notrec_e = new TH2D("fh_spot_notrec_e","(x,y) not reconstructed electrons",100,-200,200,125,-250,250);
  fh_spot_notrec_pi = new TH2D("fh_spot_notrec_pi","(x,y) not reconstructed pions",100,-200,200,125,-250,250);
  fh_spot_prim_notrec_e = new TH2D("fh_spot_prim_notrec_e","(x,y) not rec prim electrons",100,-200,200,125,-250,250);
  fh_spot_prim_notrec_pi = new TH2D("fh_spot_prim_notrec_pi","(x,y) not rec prim pions",100,-200,200,125,-250,250);

  fh_spot_prim_notrec_e_r = new TH1D("fh_spot_prim_notrec_e_r","radius not rec prim electrons",100,0,10);
  fh_spot_prim_notrec_pi_r = new TH1D("fh_spot_prim_notrec_pi_r","radius not rec prim pions",100,0,10);
  fh_spot_prim_notrec_e_p = new TH1D("fh_spot_prim_notrec_e_p","p not rec prim electrons",100,0,30);
  fh_spot_prim_notrec_pi_p = new TH1D("fh_spot_prim_notrec_pi_p","p not rec prim pions",100,0,30);

//  fh_p_notrec_pi = new TH2D("fh_p_notrec_pi","momentum p, not reconstructed pions",100,-200,200,125,-250,250);

  // (x,y) of fake rings

  fh_fake_ring = new TH2D("fh_fake_ring","(x,y) fake rings",100,-200,200,125,-250,250);

  fh_r_points_true = new TH2D("fh_r_points_true","R vs Nhits, true",50,0,50,80,0,8);
  fh_r_points_fake = new TH2D("fh_r_points_fake","R vs Nhits, fake",50,0,50,80,0,8);

  fh_fake_p_radius = new TH2D("fh_fake_p_radius","radius in dependence on p, fake rings",80,0.,8.,150,0.,15.);

  //---------------------------------------------------------------------------------------
  // tracks HAVE N_STS cut (estimate "reconstructable" tracks)
  // particles with MC ring (light spot) and correctly matched track

  fh_spotmatchtrue_e = new TH2D("fh_spotmatchtrue_e","MC ring + matched e (y,pt)",50,0,5,40,0,2);
  fh_spotmatchtrue_pi = new TH2D("fh_spotmatchtrue_pi","MC ring + matched pi (y,pt)",50,0,5,40,0,2);
  fh_spotmatchtrue_mu = new TH2D("fh_spotmatchtrue_mu","MC ring + matched mu (y,pt)",50,0,5,40,0,2);

  // particles with MC ring (light spot) and wrongly matched track: (y,pt) of MC ring

  fh_spotmatchfalse_e = new TH2D("fh_spotmatchfalse_e","MC ring + wrongly matched e (y,pt)-ring",50,0,5,40,0,2);
  fh_spotmatchfalse_pi = new TH2D("fh_spotmatchfalse_pi","MC ring + wrongly matched pi (y,pt)-ring",50,0,5,40,0,2);
  fh_spotmatchfalse_mu = new TH2D("fh_spotmatchfalse_mu","MC ring + wrongly matched mu (y,pt)-ring",50,0,5,40,0,2);

  // particles with MC ring (light spot) and wrongly matched track: (y,pt) of track

  fh_spotmatchfalse_track_e = new TH2D("fh_spotmatchfalse_track_e","MC ring + wrongly matched e (y,pt)-track",50,0,5,40,0,2);
  fh_spotmatchfalse_track_pi = new TH2D("fh_spotmatchfalse_track_pi","MC ring + wrongly matched pi (y,pt)-track",50,0,5,40,0,2);
  fh_spotmatchfalse_track_mu = new TH2D("fh_spotmatchfalse_track_mu","MC ring + wrongly matched mu (y,pt)-track",50,0,5,40,0,2);


  // particles with reconstructed ring and correctly matched track

  fh_ringmatchtrue_e = new TH2D("fh_ringmatchtrue_e","reconstructed + matched e (y,pt)",50,0,5,40,0,2);
  fh_ringmatchtrue_pi = new TH2D("fh_ringmatchtrue_pi","reconstructed + matched pi (y,pt)",50,0,5,40,0,2);
  fh_ringmatchtrue_mu = new TH2D("fh_ringmatchtrue_mu","reconstructed + matched mu (y,pt)",50,0,5,40,0,2);

  fh_ringmatchtrue_mv_e = new TH2D("fh_ringmatchtrue_mv_e","reconstructed + matched e (y,pt) mv",50,0,5,40,0,2);
  fh_ringmatchtrue_mv_pi = new TH2D("fh_ringmatchtrue_mv_pi","reconstructed + matched pi (y,pt) mv",50,0,5,40,0,2);

  // particles with reconstructed ring and wrongly matched track: (y,pt) of ring

  fh_ringmatchfalse_e = new TH2D("fh_ringmatchfalse_e","rec + wrongly matched e (y,pt)-ring",50,0,5,40,0,2);
  fh_ringmatchfalse_pi = new TH2D("fh_ringmatchfalse_pi","rec + wrongly matched pi (y,pt)-ring",50,0,5,40,0,2);
  fh_ringmatchfalse_mu = new TH2D("fh_ringmatchfalse_mu","rec + wrongly matched mu (y,pt)-ring",50,0,5,40,0,2);

  // particles with reconstructed ring and wrongly matched track: (y,pt) of track

  fh_ringmatchfalse_track_e = new TH2D("fh_ringmatchfalse_track_e","rec + wrongly matched e (y,pt)-track",50,0,5,40,0,2);
  fh_ringmatchfalse_track_pi = new TH2D("fh_ringmatchfalse_track_pi","rec + wrongly matched pi (y,pt)-track",50,0,5,40,0,2);
  fh_ringmatchfalse_track_mu = new TH2D("fh_ringmatchfalse_track_mu","rec + wrongly matched mu (y,pt)-track",50,0,5,40,0,2);

  fh_ringmatchfalse_track_mv_e = new TH2D("fh_ringmatchfalse_track_mv_e","rec + wrongly matched e (y,pt)-track mv",50,0,5,40,0,2);
  fh_ringmatchfalse_track_mv_fake_e = new TH2D("fh_ringmatchfalse_track_mv_fake_e","rec + wrongly matched fake e (y,pt)-track mv",50,0,5,40,0,2);
  fh_ringmatchfalse_track_mv_pi = new TH2D("fh_ringmatchfalse_track_mv_pi","rec + wrongly matched pi (y,pt)-track mv",50,0,5,40,0,2);


//----------------------------------------------------------------------------------------------------
  // identified particles

  fh_id_e = new TH2D("fh_id_e","identified e (y,pt)",50,0,5,40,0,2);
  fh_id_pi = new TH2D("fh_id_pi","identified pi (y,pt)",50,0,5,40,0,2);
  fh_id_mu = new TH2D("fh_id_mu","identified mu (y,pt)",50,0,5,40,0,2);
  fh_id_k = new TH2D("fh_id_k","identified K (y,pt)",50,0,5,40,0,2);

  // misidentified particles

  fh_misid_e = new TH2D("fh_misid_e","misidentified e (y,pt)",50,0,5,40,0,2);
  fh_misid_pi = new TH2D("fh_misid_pi","misidentified pi (y,pt)",50,0,5,40,0,2);
  fh_misid_mu = new TH2D("fh_misid_mu","misidentified mu (y,pt)",50,0,5,40,0,2);
  fh_misid_k = new TH2D("fh_misid_k","misidentified K (y,pt)",50,0,5,40,0,2);

  // additional (fake)  rings

  fh_ringfake_e = new TH2D("fh_ringfake_e","fake ring e (y,pt)",50,0,5,40,0,2);
  fh_ringfake_pi = new TH2D("fh_ringfake_pi","fake ring pi (y,pt)",50,0,5,40,0,2);
  fh_ringfake_mu = new TH2D("fh_ringfake_mu","fake ring mu (y,pt)",50,0,5,40,0,2);
  fh_ringfake_k = new TH2D("fh_ringfake_k","fake ring K (y,pt)",50,0,5,40,0,2);


  // Distance between ring center and track extrapolation for closest distance match

  fh_disttrue_e = new TH3D("fh_disttrue_e","ring-track distance (true match) in dependence on (p,pt) e ",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_disttrue_pi = new TH3D("fh_disttrue_pi","ring-track distance (true match) in dependence on (p,pt) pi",150,0.,15.,20,0.,2.,60,0.,6.);

  fh_distfalse_e = new TH3D("fh_distfalse_e","ring-track distance (false match) in dependence on (p,pt) e",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_distfalse_fake_e = new TH3D("fh_distfalse_fake_e","ring-track distance (false match, fake rings) in dependence on (p,pt) e",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_distfalse_pi = new TH3D("fh_distfalse_pi","ring-track distance (false match) in dependence on (p,pt) pi",150,0.,15.,20,0.,2.,60,0.,6.);

  fh_disttrue_mv_e = new TH3D("fh_disttrue_mv_e","ring-track distance (true match) in dependence on (p,pt) e mv",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_disttrue_mv_pi = new TH3D("fh_disttrue_mv_pi","ring-track distance (true match) in dependence on (p,pt) pi mv",150,0.,15.,20,0.,2.,60,0.,6.);

  fh_distfalse_mv_e = new TH3D("fh_distfalse_mv_e","ring-track distance (false match) in dependence on (p,pt) e mv",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_distfalse_mv_fake_e = new TH3D("fh_distfalse_mv_fake_e","ring-track distance (false match, fake rings) in dependence on (p,pt) e mv",150,0.,15.,20,0.,2.,60,0.,6.);
  fh_distfalse_mv_pi = new TH3D("fh_distfalse_mv_pi","ring-track distance (false match) in dependence on (p,pt) pi mv",150,0.,15.,20,0.,2.,60,0.,6.);


  // Radius versus p in bins of pt

  fh_radius = new TH3D("fh_radius","radius in dependence on (p,pt)",150,0.,15.,20,0.,2.,80,0.,8.);

}

//------------------------------------------------------------------------------
CbmRichID::~CbmRichID()
{
 // generated particles

  fh_gen_e  -> Write();
  fh_gen_pi -> Write();
  fh_gen_mu -> Write();
  fh_gen_k  -> Write();

  // measured particles

  fh_meas_e  -> Write();
  fh_meas_pi -> Write();
  fh_meas_mu -> Write();
  fh_meas_k  -> Write();

  // accepted particles (N_sts >= fNminSTS)

  fh_acc_e  ->Divide(fh_meas_e,fh_gen_e,1,1);
  fh_acc_pi ->Divide(fh_meas_pi,fh_gen_pi,1,1);
  fh_acc_mu ->Divide(fh_meas_mu,fh_gen_mu,1,1);
  fh_acc_k  ->Divide(fh_meas_k,fh_gen_k,1,1);

  fh_acc_e  -> Write();
  fh_acc_pi -> Write();
  fh_acc_mu -> Write();
  fh_acc_k  -> Write();

  // particles with reconstructed ring

  fh_ring_e  -> Write();
  fh_ring_pi -> Write();
  fh_ring_mu -> Write();
  fh_ring_k  -> Write();

  fh_ring_nsts_e  -> Write();
  fh_ring_nsts_pi -> Write();
  fh_ring_nsts_mu -> Write();

  fh_ring_mv_e  -> Write();
  fh_ring_mv_pi -> Write();

  fh_ringall_e  -> Write();
  fh_ringall_pi -> Write();
  fh_ringall_mu -> Write();

  fh_spot_e  -> Write();
  fh_spot_pi -> Write();
  fh_spot_mu -> Write();

  fh_spotall_e  -> Write();
  fh_spotall_pi -> Write();
  fh_spotall_mu -> Write();

  fh_spot_notrec_e -> Write();
  fh_spot_notrec_pi -> Write();
  fh_spot_prim_notrec_e -> Write();
  fh_spot_prim_notrec_pi -> Write();
  fh_spot_prim_notrec_e_r -> Write();
  fh_spot_prim_notrec_pi_r -> Write();
  fh_spot_prim_notrec_e_p -> Write();
  fh_spot_prim_notrec_pi_p -> Write();

  fh_fake_ring -> Write();

  fh_r_points_true -> Write();
  fh_r_points_fake -> Write();
  fh_fake_p_radius -> Write();

  // particles with MC ring (light spot) and correctly matched track

  fh_spotmatchtrue_e -> Write();
  fh_spotmatchtrue_pi -> Write();
  fh_spotmatchtrue_mu -> Write();

  // particles with MC ring (light spot) and wrongly matched track

  fh_spotmatchfalse_e -> Write();
  fh_spotmatchfalse_pi -> Write();
  fh_spotmatchfalse_mu -> Write();
  fh_spotmatchfalse_track_e -> Write();
  fh_spotmatchfalse_track_pi -> Write();
  fh_spotmatchfalse_track_mu -> Write();


  // particles with reconstructed ring and correctly matched track

  fh_ringmatchtrue_e -> Write();
  fh_ringmatchtrue_pi -> Write();
  fh_ringmatchtrue_mu -> Write();

  fh_ringmatchtrue_mv_e -> Write();
  fh_ringmatchtrue_mv_pi -> Write();

  // particles with reconstructed ring and wrongly matched track

  fh_ringmatchfalse_e -> Write();
  fh_ringmatchfalse_pi -> Write();
  fh_ringmatchfalse_mu -> Write();
  fh_ringmatchfalse_track_e -> Write();
  fh_ringmatchfalse_track_pi -> Write();
  fh_ringmatchfalse_track_mu -> Write();

  fh_ringmatchfalse_track_mv_e -> Write();
  fh_ringmatchfalse_track_mv_fake_e -> Write();
  fh_ringmatchfalse_track_mv_pi -> Write();

  // identified particles

  fh_id_e -> Write();
  fh_id_pi -> Write();
  fh_id_mu -> Write();
  fh_id_k -> Write();

  // misidentified particles

  fh_misid_e -> Write();
  fh_misid_pi -> Write();
  fh_misid_mu -> Write();
  fh_misid_k -> Write();

  // additional (fake)  rings

  fh_ringfake_e -> Write();
  fh_ringfake_pi -> Write();
  fh_ringfake_mu -> Write();
  fh_ringfake_k -> Write();


  // Distance between ring center and track extrapolation for closest distance match

  fh_disttrue_e -> Write();
  fh_disttrue_pi -> Write();

  fh_distfalse_e -> Write();
  fh_distfalse_fake_e -> Write();
  fh_distfalse_pi -> Write();

  fh_disttrue_mv_e -> Write();
  fh_disttrue_mv_pi -> Write();

  fh_distfalse_mv_e -> Write();
  fh_distfalse_mv_fake_e -> Write();
  fh_distfalse_mv_pi -> Write();


  // Radius versus p in bins of pt

  fh_radius -> Write();


}

//------------------------------------------------------------------------------
InitStatus CbmRichID::Init()
{
  cout <<"-------------------------------------------------------------------------------------------"<<endl;
  CbmRootManager *manager= CbmRootManager::Instance();

  // all tracks
  fListStack      = (TClonesArray*)manager->ActivateBranch("MCTrack");

  //RICH MC points
  fListRICHpts    = (TClonesArray*)manager->ActivateBranch("RichPoint");

  //RICH hits
  fListRICHhits   = (TClonesArray*)manager->ActivateBranch("RichHit");

  //RICH ring guidances/ track extrapolation
  fListRICHtracks  = (TClonesArray*)manager->ActivateBranch("RichProjection");

  // found rings
  fListRICHrings  = (TClonesArray*)manager->ActivateBranch("RichRing");

  // light spots
  fListLightSpots = new TClonesArray("CbmRichLightSpot",500);

  return kSUCCESS;
}

//------------------------------------------------------------------------------
void CbmRichID::Exec(Option_t* option)
{
  // make RICH ID analysis

  printf("RICH ID analysis: event %d\n",fEvent);
   
  AnalyzeAcceptance();
  AnalyzeMatching_Spots();
  DefineRingID();
  AnalyzeMatching_Rings();

  fEvent++;
}


//------------------------------------------------------------------------------
void CbmRichID::AnalyzeAcceptance()
{
  // Primary track analysis

  CbmMCTrack *primary=NULL;
  fNprimary = 0;
  TLorentzVector momentum, mom_rich;

  //------------------------------------------------------------------------------------------------
  // Fill TClonesArray with LightSpots

  if (strstr(fDebug,"hit"))
    printf("\tNumber of RICH hits: %d\n",fListRICHhits->GetEntriesFast());

  CbmRichHit* hit=NULL;

  Int_t nLightSpots = 0;
//  fNHits=0;
  Int_t iLightSpot =0;
  fListLightSpots->Clear();

  for(Int_t iHit=0; iHit<fListRICHhits->GetEntriesFast(); iHit++) {
    hit = (CbmRichHit*)fListRICHhits->At(iHit);                                    // RICH hit
    if (hit->GetRefIndex() < 0) continue; // fake hits
    CbmRichPoint* point = (CbmRichPoint*)fListRICHpts->At(hit->GetRefIndex());
    CbmMCTrack* hitTrack = (CbmMCTrack*)fListStack->At(point->GetTrackID());
    if(hitTrack->GetPdgCode() != 50000050) continue;                               // select only Cherenkov photons

    // Fill array with cherenkov photon parents: ALL tracks (no NSTS cut!)
    Int_t  motherTrack = hitTrack->GetMotherID();
    Bool_t dejaVue = kFALSE;
    CbmRichLightSpot *newspot;
    for (iLightSpot=0; iLightSpot<nLightSpots; iLightSpot++) {
      newspot = (CbmRichLightSpot*)fListLightSpots->At(iLightSpot);
      if (newspot->GetTrackID() == motherTrack) {
	dejaVue = kTRUE;
	newspot->AddHitIndex(iHit);                               // amplitudes not considered (e.g. if two photons per hit)
	break;
      }
    }
    if (!dejaVue) {
      newspot = new((*fListLightSpots)[nLightSpots]) CbmRichLightSpot();
      newspot->SetTrackID(motherTrack);
      newspot->AddHitIndex(iHit);
      nLightSpots++;
    }
  }

//  for (iLightSpot=0; iLightSpot<nLightSpots; iLightSpot++)
//    ((CbmRichLightSpot*)fListLightSpots->At(iLightSpot))->RingFitCOP(); //??? noetig?

  cout  <<"================ found "
      <<fListLightSpots->GetEntriesFast()
      <<" LightSpots. "<<endl;

  //-----------TClonesArray filled -------------------------------------------------------------------------------------

  for (Int_t iPrimary=0; iPrimary<fListStack->GetEntriesFast(); iPrimary++) {
    primary = (CbmMCTrack*)fListStack->At(iPrimary);
//    TVector3 vertex = primary->GetStartVertex();
//    if (vertex.X() > fvx || vertex.Y() > fvy || vertex.Z() > fvz || primary->GetMotherID() == -1) continue;
    if (primary->GetMotherID() != -1) continue;
//    if (TMath::Abs(primary->GetPdgCode()) == 11 && primary->GetStsPoints() >= 6) fNprimary++;
    fNprimary++;
    momentum = primary->Get4Momentum();

    if (strstr(fDebug,"prim"))
      printf("track %d: id=%d, p=(%f,%f,%f,%f) GeV/c, pt,y,theta=%f,%f,%f\n",
	     iPrimary,
	     primary->GetPdgCode(),
	     momentum.Px(),
	     momentum.Py(),
	     momentum.Pz(),
	     momentum.Energy(),
	     momentum.Pt(),
	     momentum.Eta(),
	     momentum.Theta());

    if (TMath::Abs(primary->GetPdgCode()) == 11) fh_gen_e->Fill(momentum.Rapidity(),momentum.Pt());
    if (TMath::Abs(primary->GetPdgCode()) == 211) fh_gen_pi->Fill(momentum.Rapidity(),momentum.Pt());
    if (TMath::Abs(primary->GetPdgCode()) == 13) fh_gen_mu->Fill(momentum.Rapidity(),momentum.Pt());
    if (TMath::Abs(primary->GetPdgCode()) == 321) fh_gen_k->Fill(momentum.Rapidity(),momentum.Pt());


    CbmRichLightSpot *spot;
    for ( iLightSpot=0; iLightSpot<nLightSpots; iLightSpot++) {
      spot = (CbmRichLightSpot*)fListLightSpots->At(iLightSpot);
      Int_t fSpotID = spot->GetTrackID();
      if (fSpotID == iPrimary) {
        CbmMCTrack *pm= (CbmMCTrack*) fListStack->At(fSpotID);
         Int_t fnsts=pm->GetStsPoints();
	 Int_t fnhits = spot->GetNhits();
//	 Int_t fnrich=pm->GetRichPoints();  // 10 if RichMirrorPoint, 11 if RichPoint (Cherenkovs only)
//	 mom_rich = pm->Get4Momentum();

         if(fnsts >= fNminSTS && fnhits >= fNminHits){
          if (TMath::Abs(primary->GetPdgCode()) == 11) fh_meas_e->Fill(momentum.Rapidity(),momentum.Pt());
          if (TMath::Abs(primary->GetPdgCode()) == 211) fh_meas_pi->Fill(momentum.Rapidity(),momentum.Pt());
          if (TMath::Abs(primary->GetPdgCode()) == 13) fh_meas_mu->Fill(momentum.Rapidity(),momentum.Pt());
          if (TMath::Abs(primary->GetPdgCode()) == 321) fh_meas_k->Fill(momentum.Rapidity(),momentum.Pt());
	  }  // more than fNminSTS points required  (add here requirement on RICH points?!
         } // primary in RICH?
      }  // loop RICH LightSpots

  } // loop over primaries

  if (strstr(fDebug,"prim"))
    printf("\tNumber of primaries: %d\n",fNprimary);


}


//------------------------------------------------------------------------------
void CbmRichID::AnalyzeMatching_Spots()
{
  // Ring-track matching analysis

  if (strstr(fDebug,"match"))
    printf("\tNumber of extrapolated RICH tracks: %d\n",fListRICHtracks->GetEntriesFast());
  CbmRichRing* track;

  cout  <<"================ found "
      <<fListRICHtracks->GetEntriesFast()
      <<" Track Extrapolations "<<endl;

  // Double loop over RICH hits and tracks
  CbmRichLightSpot *spot;
//  nLightSpots = fListLightSpots->GetEntriesFast();
  for (Int_t iLightSpot=0; iLightSpot<fListLightSpots->GetEntriesFast(); iLightSpot++) {
    spot = (CbmRichLightSpot*)fListLightSpots->At(iLightSpot);
    spot->RingFit();
//    if (spot->GetNhits() < 20) continue;
    Double_t xSpot = spot->GetCenterX();
    Double_t ySpot = spot->GetCenterY();

    TLorentzVector momentum = ((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->Get4Momentum();
    Int_t pdg = ((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetPdgCode();
//    Int_t Nsts =((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetStSPoints();
//    if (Nsts >= fNminSTS) {
    if (pdg == 11 || pdg == -11) fh_spotall_e->Fill(momentum.Rapidity(),momentum.Pt());
    if (pdg == 211 || pdg == -211) fh_spotall_pi->Fill(momentum.Rapidity(),momentum.Pt());
    if (pdg == 13 || pdg == -13) fh_spotall_mu->Fill(momentum.Rapidity(),momentum.Pt());
    if (((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetMotherID() == -1){      // originally generated e,pi,mu
       if (pdg == 11 || pdg == -11) fh_spot_e->Fill(momentum.Rapidity(),momentum.Pt());
       if (pdg == 211 || pdg == -211) fh_spot_pi->Fill(momentum.Rapidity(),momentum.Pt());
       if (pdg == 13 || pdg == -13) fh_spot_mu->Fill(momentum.Rapidity(),momentum.Pt());
     }
//     }

    Double_t rMin = 99.;
    Int_t idRingMin = 0;
    Int_t idTrackMin = 0;
    for (Int_t itrack=0; itrack<fListRICHtracks->GetEntriesFast(); itrack++) {
      track = (CbmRichRing*)fListRICHtracks->At(itrack);
      Double_t xTrack = track->GetCenterX();
      Double_t yTrack = track->GetCenterY();
      Float_t dist_SpotTrack = TMath::Sqrt( (xSpot-xTrack)*(xSpot-xTrack) +
				            (ySpot-yTrack)*(ySpot-yTrack) );

      // search for ring-track combination with shortest distance
      if (dist_SpotTrack < rMin){
	rMin = dist_SpotTrack;
        idRingMin = spot->GetTrackID();
        idTrackMin = track->GetTrackID();
      }
    }  // loop over tracks (ring center guidances)

    if (rMin != 99.){
      TLorentzVector momTrack = ((CbmMCTrack*)fListStack->At(idTrackMin))->Get4Momentum();
      TLorentzVector momRing = ((CbmMCTrack*)fListStack->At(idTrackMin))->Get4Momentum();
      if (idRingMin == idTrackMin){  // correct match
	Int_t pdg = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode();
	if (rMin < fDistanceCut){
	 if (((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetMotherID() == -1){
	  if (pdg == 11 || pdg == -11) fh_spotmatchtrue_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdg == 211 || pdg == -211) fh_spotmatchtrue_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdg == 13 || pdg == -13) fh_spotmatchtrue_mu->Fill(momTrack.Rapidity(),momTrack.Pt());
	  } // primaries!
	}
      }
      else {
	Int_t pdgRing = ((CbmMCTrack*)fListStack->At(idRingMin))->GetPdgCode();
//	Int_t pdgTrack = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode();
        if (rMin < fDistanceCut){
         if (((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetMotherID() == -1){
	  if (pdgRing == 11 || pdgRing == -11) fh_spotmatchfalse_e->Fill(momRing.Rapidity(),momRing.Pt());
	  if (pdgRing == 211 || pdgRing == -211) fh_spotmatchfalse_pi->Fill(momRing.Rapidity(),momRing.Pt());
	  if (pdgRing == 13 || pdgRing == -13) fh_spotmatchfalse_mu->Fill(momRing.Rapidity(),momRing.Pt());

	  if (pdgRing == 11 || pdgRing == -11) fh_spotmatchfalse_track_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdgRing == 211 || pdgRing == -211) fh_spotmatchfalse_track_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdgRing == 13 || pdgRing == -13) fh_spotmatchfalse_track_mu->Fill(momTrack.Rapidity(),momTrack.Pt());
	  }   // primaries!
	}
      }
    }  // plot values for closest distance match
  } // loop over rings (LightSpots)


  }


  //------------------------------------------------------------------------------
void CbmRichID::DefineRingID()
{

// Double loop over RICH hits and tracks
  CbmRichLightSpot *ring;
  CbmRichLightSpot *spot;

  Int_t nIdRing = 0;
  Int_t nring_vertex=0;

  Int_t nRing = fListRICHrings->GetEntriesFast();
  cout  <<"================ found "
      << nRing
      <<" Rings "<<endl;

  for (Int_t iring=0; iring<nRing; iring++) {
    ring = (CbmRichLightSpot*)fListRICHrings->At(iring);
    ring->SetRecFlag(-1);
    Double_t xRing = ring->GetCenterX();
    Double_t yRing = ring->GetCenterY();
    Double_t rRing = ring->GetRadius();
    Int_t hitsRing = ring->GetNhits();
    Int_t nLightSpots = fListLightSpots->GetEntriesFast();

//    cout << xRing << " " << yRing << " " << rRing << endl;

    Double_t rMin = 99.;
    Int_t idSpotMin = 0;
    Double_t potential_hits = 0.;
    Double_t spot_radius = 0.;

    for (Int_t ispot=0; ispot<nLightSpots; ispot++){
     spot = (CbmRichLightSpot*)fListLightSpots->At(ispot);
     spot->SetRecFlag(-1);
     Double_t xSpot = spot->GetCenterX();
     Double_t ySpot = spot->GetCenterY();
     Double_t rSpot = spot->GetRadius();


     if (iring == 0) {
       Int_t dum =((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetMotherID();
       if (dum == -1) {
         nring_vertex++;
//         cout << spot->GetTrackID() << endl;
//	 cout << xSpot << " " << ySpot << " " << rSpot << endl;
	 }
      }


     Float_t dist_SpotRing = TMath::Sqrt( (xSpot-xRing)*(xSpot-xRing) +
				          (ySpot-yRing)*(ySpot-yRing) );
     Float_t diff_radius = TMath::Abs(rSpot - rRing);

     if ((dist_SpotRing < rMin) && (diff_radius < fDiffRadiusMin)){
      rMin = dist_SpotRing;
      idSpotMin = spot->GetTrackID();
      potential_hits = spot->GetNhits();
      spot_radius = rSpot;
      }
     }      // loop spots

//     cout << "ring ID " << rMin << " " << rRing-spot_radius << " " << hitsRing << " " << potential_hits << " " << idSpotMin << endl;

     if (rMin < fDistRingID && (Double_t)(hitsRing)/(Double_t)(potential_hits) > fPotentialHitsMin &&
                               (Double_t)(hitsRing)/(Double_t)(potential_hits) < fPotentialHitsMax){
      ring->SetTrackID(idSpotMin);          // ring can be counted as correctly reconstructed
      ring->SetRecFlag(1);

      // mark corresponding spot as found:
      for (Int_t ispot=0; ispot<nLightSpots; ispot++){
       spot = (CbmRichLightSpot*)fListLightSpots->At(ispot);
       if (spot->GetTrackID() == idSpotMin) {
         spot->SetRecFlag(1);
	 fh_r_points_true->Fill((Double_t)(hitsRing),rRing);
	 if (ring->GetNhits() >= fNminHits && ((CbmMCTrack*)fListStack->At(idSpotMin))->GetStsPoints() >= fNminSTS){
          spot->SetRecFlag(2);
	  ring->SetRecFlag(2);
	  }
	 }
      }

      TLorentzVector momentum = ((CbmMCTrack*)fListStack->At(idSpotMin))->Get4Momentum();
      Int_t pdg = ((CbmMCTrack*)fListStack->At(idSpotMin))->GetPdgCode();
      if (((CbmMCTrack*)fListStack->At(idSpotMin))->GetMotherID() == -1){      // originally generated e,pi,mu
       if (pdg == 11 || pdg == -11) fh_ring_e->Fill(momentum.Rapidity(),momentum.Pt());
       if (pdg == 211 || pdg == -211) fh_ring_pi->Fill(momentum.Rapidity(),momentum.Pt());
       if (pdg == 13 || pdg == -13) fh_ring_mu->Fill(momentum.Rapidity(),momentum.Pt());
       if (ring->GetRecFlag() == 2){
         if (pdg == 11 || pdg == -11) fh_ring_nsts_e->Fill(momentum.Rapidity(),momentum.Pt());
         if (pdg == 211 || pdg == -211) fh_ring_nsts_pi->Fill(momentum.Rapidity(),momentum.Pt());
         if (pdg == 13 || pdg == -13) fh_ring_nsts_mu->Fill(momentum.Rapidity(),momentum.Pt());
	 }
        }
      TVector3 vertex = ((CbmMCTrack*)fListStack->At(idSpotMin))->GetStartVertex();
      if (vertex.X() < fvx && vertex.Y() < fvy && vertex.Z() < fvz && ring->GetRecFlag() == 2){ // select only good main vertex tracks
        if (pdg == 11 || pdg == -11) fh_ring_mv_e->Fill(momentum.Rapidity(),momentum.Pt());
	if (pdg == 211 || pdg == -211) fh_ring_mv_pi->Fill(momentum.Rapidity(),momentum.Pt());
	}

      if (pdg == 11 || pdg == -11) fh_ringall_e->Fill(momentum.Rapidity(),momentum.Pt());
      if (pdg == 211 || pdg == -211) fh_ringall_pi->Fill(momentum.Rapidity(),momentum.Pt());
      if (pdg == 13 || pdg == -13) fh_ringall_mu->Fill(momentum.Rapidity(),momentum.Pt());
//      cout << "ring ID " << rMin << " " << rRing-spot_radius << " " << hitsRing << " " << potential_hits << " " << idSpotMin << endl;
      nIdRing++;
     }
     else
      ring->SetTrackID(-1);    // fake rings - not connected to a light-spot, they also have RecFlag = -1

// todo: include comparison of hitarrays into definition of "correct" Ring ID

    // mark fake rings which look like good e-rings
     if (ring->GetRecFlag() == -1){
       fh_fake_ring->Fill(xRing,yRing);
       fh_r_points_fake->Fill((Double_t)(hitsRing),rRing);
       if (rRing > fRmin && rRing < fRmax && hitsRing >= fNminHits) ring->SetRecFlag(5);  // fakes looking like good e-rings
     }

    }  // loop over rings

    // loop spots and plot not reconstructed ones
     for (Int_t ispot=0; ispot<fListLightSpots->GetEntriesFast(); ispot++){
       spot = (CbmRichLightSpot*)fListLightSpots->At(ispot);
       TLorentzVector momentum = ((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->Get4Momentum();
       Int_t pdg = ((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetPdgCode();
       if (spot->GetRecFlag() == -1) {
        if(pdg ==11 || pdg == -11) fh_spot_notrec_e->Fill(spot->GetCenterX(),spot->GetCenterY());
        if(pdg ==211 || pdg == -211) fh_spot_notrec_pi->Fill(spot->GetCenterX(),spot->GetCenterY());
	if (((CbmMCTrack*)fListStack->At(spot->GetTrackID()))->GetMotherID() == -1){
          if(pdg ==11 || pdg == -11) {
	    fh_spot_prim_notrec_e->Fill(spot->GetCenterX(),spot->GetCenterY());
	    fh_spot_prim_notrec_e_r->Fill(spot->GetRadius());
	    fh_spot_prim_notrec_e_p->Fill(momentum.P());
	    }
          if(pdg ==211 || pdg == -211) {
	    fh_spot_prim_notrec_pi->Fill(spot->GetCenterX(),spot->GetCenterY());
	    fh_spot_prim_notrec_pi_r->Fill(spot->GetRadius());
	    fh_spot_prim_notrec_pi_p->Fill(momentum.P());
	    }
	}
	}
      }


  cout  <<"================ found "
      << nIdRing
      <<" identified rings "
      << nring_vertex
      << " spots from vertex "<<endl;

}
//-------------------------------------------------------------------------------

  //------------------------------------------------------------------------------
void CbmRichID::AnalyzeMatching_Rings()
{
  // Ring-track matching analysis

  if (strstr(fDebug,"match"))
    printf("\tNumber of extrapolated RICH tracks: %d\n",fListRICHtracks->GetEntriesFast());
  CbmRichRing* track;

  // Double loop over RICH hits and tracks
  CbmRichLightSpot *ring;

  Int_t nRing = fListRICHrings->GetEntriesFast();
  for (Int_t iring=0; iring<nRing; iring++) {
    ring = (CbmRichLightSpot*)fListRICHrings->At(iring);
//    if (ring->GetNhits() < fNminHits) continue;
    Double_t xRing = ring->GetCenterX();
    Double_t yRing = ring->GetCenterY();
    Double_t rRing = ring->GetRadius();
    Double_t rMin = 99.;
    Int_t idRingMin = 0;
    Int_t idTrackMin = 0;
    for (Int_t itrack=0; itrack<fListRICHtracks->GetEntriesFast(); itrack++) {
      track = (CbmRichRing*)fListRICHtracks->At(itrack);
      Double_t xTrack = track->GetCenterX();
      Double_t yTrack = track->GetCenterY();
      Float_t dist_RingTrack = TMath::Sqrt( (xRing-xTrack)*(xRing-xTrack) +
				            (yRing-yTrack)*(yRing-yTrack) );

      // search for ring-track combination with shortest distance
      if (dist_RingTrack < rMin){
	rMin = dist_RingTrack;
        idRingMin = ring->GetTrackID();
        idTrackMin = track->GetTrackID();
      }
    }  // loop over tracks (ring center guidances)

    TLorentzVector momTrack = ((CbmMCTrack*)fListStack->At(idTrackMin))->Get4Momentum();
    TVector3 vertex = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetStartVertex();
    if (rMin != 99. && idRingMin > -1){
      TLorentzVector momRing = ((CbmMCTrack*)fListStack->At(idRingMin))->Get4Momentum();
      fh_radius->Fill(momRing.P(),momRing.Pt(),rRing);
      if (idRingMin == idTrackMin){  // correct match!
	Int_t pdg = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode();
	if (pdg == 11 || pdg == -11) fh_disttrue_e ->Fill(momTrack.P(),momTrack.Pt(),rMin);
	if (pdg == 211 || pdg == -211) fh_disttrue_pi ->Fill(momTrack.P(),momTrack.Pt(),rMin);
	if (vertex.X() < fvx && vertex.Y() < fvy && vertex.Z() < fvz && ring->GetRecFlag() == 2){ // select only good main vertex tracks
          if (pdg == 11 || pdg == -11) fh_disttrue_mv_e ->Fill(momTrack.P(),momTrack.Pt(),rMin);
	  if (pdg == 211 || pdg == -211) fh_disttrue_mv_pi ->Fill(momTrack.P(),momTrack.Pt(),rMin);
          if (rMin < fDistanceCut){
            if (pdg == 11 || pdg == -11) fh_ringmatchtrue_mv_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	    if (pdg == 211 || pdg == -211) fh_ringmatchtrue_mv_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  }
	}
	if (((CbmMCTrack*)fListStack->At(idTrackMin))->GetMotherID() == -1){                         // select only primaries
	if (rMin < fDistanceCut && ring->GetRecFlag() == 2){
	  if (pdg == 11 || pdg == -11) fh_ringmatchtrue_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdg == 211 || pdg == -211) fh_ringmatchtrue_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdg == 13 || pdg == -13) fh_ringmatchtrue_mu->Fill(momTrack.Rapidity(),momTrack.Pt());
	}
	}  // select primaries
      }
      else {
	Int_t pdgRing = ((CbmMCTrack*)fListStack->At(idRingMin))->GetPdgCode();
//	Int_t pdgTrack = ((CbmMCTrack*)fListStack->At(idTrackMin))->GetPdgCode();
	if (pdgRing == 11 || pdgRing == -11) fh_distfalse_e ->Fill(momRing.P(),momRing.Pt(),rMin);
	if (pdgRing == 211 || pdgRing == -211) fh_distfalse_pi ->Fill(momRing.P(),momRing.Pt(),rMin);
        if (vertex.X() < fvx && vertex.Y() < fvy && vertex.Z() < fvz && ring->GetRecFlag() == 2){ // select only good main vertex tracks
          if (pdgRing == 11 || pdgRing == -11) fh_distfalse_mv_e ->Fill(momRing.P(),momRing.Pt(),rMin);
	  if (pdgRing == 211 || pdgRing == -211) fh_distfalse_mv_pi ->Fill(momRing.P(),momRing.Pt(),rMin);
          if (rMin < fDistanceCut){
            if (pdgRing == 11 || pdgRing == -11) fh_ringmatchfalse_track_mv_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	    if (pdgRing == 211 || pdgRing == -211) fh_ringmatchfalse_track_mv_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  }
	}
	if (((CbmMCTrack*)fListStack->At(idRingMin))->GetMotherID() == -1){
        if (rMin < fDistanceCut && ring->GetRecFlag() == 2){
	  if (pdgRing == 11 || pdgRing == -11) fh_ringmatchfalse_e->Fill(momRing.Rapidity(),momRing.Pt());
	  if (pdgRing == 211 || pdgRing == -211) fh_ringmatchfalse_pi->Fill(momRing.Rapidity(),momRing.Pt());
	  if (pdgRing == 13 || pdgRing == -13) fh_ringmatchfalse_mu->Fill(momRing.Rapidity(),momRing.Pt());

	  if (pdgRing == 11 || pdgRing == -11) fh_ringmatchfalse_track_e->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdgRing == 211 || pdgRing == -211) fh_ringmatchfalse_track_pi->Fill(momTrack.Rapidity(),momTrack.Pt());
	  if (pdgRing == 13 || pdgRing == -13) fh_ringmatchfalse_track_mu->Fill(momTrack.Rapidity(),momTrack.Pt());
	}
	}  // select primaries
      }
    }  // plot values for closest distance match for "non-fake" rings only!

    if (rMin != 99. && ring->GetRecFlag() == -1) fh_fake_p_radius->Fill(rRing,momTrack.P());   // all other fake rings
    if (rMin != 99. && ring->GetRecFlag() == 5){                              // fake e-rings
     fh_fake_p_radius->Fill(rRing,momTrack.P());
     fh_distfalse_fake_e ->Fill(momTrack.P(),momTrack.Pt(),rMin);
     if (vertex.X() < fvx && vertex.Y() < fvy && vertex.Z() < fvz &&
         ((CbmMCTrack*)fListStack->At(idTrackMin))->GetStsPoints() >= fNminSTS) {
       fh_distfalse_mv_fake_e ->Fill(momTrack.P(),momTrack.Pt(),rMin);
       if (rMin < fDistanceCut)  fh_ringmatchfalse_track_mv_fake_e->Fill(momTrack.Rapidity(),momTrack.Pt());
      }   // select only good main vertex tracks

    }    // plot values for closest distance match for "fake" rings only! here: ring ID is missing!

  } // loop over rings


  }


//------------------------------------------------------------------------------
void CbmRichID::Finish()
{
//  fListRICHpts->Clear();
  fListRICHhits->Clear();
  fListRICHtracks->Clear();
  fListRICHrings->Clear();
  fListLightSpots->Clear();
  fListStack->Clear();
}


//------------------------------------------------------------------------------
ClassImp(CbmRichID)