/******************************************************************************
*  $Id: CbmRichRingQa.cxx,v 1.3 2006/06/19 11:19:52 turany Exp $
*
*  Class  : CbmRichRingQa
*  Description : Quality checks or ring finders:
*                efficiency calculation etc.
*
*  Author : Simeon Lebedev
*
*******************************************************************************
*  $Log: CbmRichRingQa.cxx,v $
*  Revision 1.3  2006/06/19 11:19:52  turany
*  map<Int_t, Int_t>iterator it do not need to be set to zero, it will be reseted by calling it.begin(), putting it=0 gives a compiler error with the new compilers
*
*  Revision 1.2  2006/04/17 22:11:04  sgorboun
*  changes in L1 ENN Ring Finder
*
*  Revision 1.1  2006/04/10 08:21:03  hoehne
*  routine for checking the ring finders (quality, efficiency)
*  initial version
*
*
*
*******************************************************************************/

#include <map>
#include <iostream>

#include "CbmRichRingQa.h"
#include "CbmTask.h"
#include "CbmRootManager.h"
#include "CbmMCTrack.h"
#include "CbmRichPoint.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "CbmTrackParam.h"
#include "CbmGlobalTrack.h"
#include "CbmMCTrack.h"
#include "CbmRichRingMatch.h"
#include "CbmStsTrackMatch.h"
// -----   Default constructor   -------------------------------------------
CbmRichRingQa::CbmRichRingQa() :CbmTask("RichRingQa")
{
  fNormType = 0;
}
// -------------------------------------------------------------------------

//------------  standard constructor (with verbosity level)  ---------------------------------
CbmRichRingQa::CbmRichRingQa(const char *name, const char *title, Int_t verbose, Int_t normType)
  :CbmTask(name)
{

  // verbosity level
  fVerbose = verbose;
  fNormType = normType;


  // count events
  nevents = 0;
  fh_nAllRings = new TH1D("fh_nAllRings","Number of All RICH rings",50,0,150);
  fh_n5HitsAllRings = new TH1D("fh_n5HitsAllRings","Number of RICH rings. >= 5 hits",50,0,150);
  fh_n10HitsAllRings = new TH1D("fh_n10HitsAllRings","Number of RICH rings. >= 10 hits",50,0,150);
  fh_n15HitsAllRings = new TH1D("fh_n15HitsAllRings","Number of RICH rings. >= 15hits",50,0,150);

  fh_nERings = new TH1D("fh_nERings","Number of all electron RICH rings",50,0,150);
  fh_n5HitsERings = new TH1D("fh_n5HitsERings","Number of electron RICH rings. >= 5 hits",50,0,150);
  fh_n10HitsERings = new TH1D("fh_n10HitsERings","Number of electron RICH rings. >= 10 hits",50,0,150);
  fh_n15HitsERings = new TH1D("fh_n15HitsERings","Number of electron RICH rings. >= 15hits",50,0,150);

  fh_nPiRings = new TH1D("fh_nPiRings","Number of all pion RICH rings",50,0,150);
  fh_n5HitsPiRings = new TH1D("fh_n5HitsPiRings","Number of pion RICH rings. >= 5 hits",50,0,150);
  fh_n10HitsPiRings = new TH1D("fh_n10HitsPiRings","Number of pion RICH rings. >= 10 hits",50,0,150);
  fh_n15HitsPiRings = new TH1D("fh_n15HitsPiRings","Number of pion RICH rings. >= 15hits",50,0,150);


  fh_n6StsAllRings = new TH1D("fh_n6StsAllRings","Number of All RICH rings. >=6 STS Points",50,0,150);
  fh_n5Hits6StsAllRings = new TH1D("fh_n5Hits6StsAllRings","Number of RICH rings. >=6 STS Points. >= 5 hits",50,0,150);
  fh_n10Hits6StsAllRings = new TH1D("fh_n10Hits6StsAllRings","Number of RICH rings. >=6 STS Points. >= 10 hits",50,0,150);
  fh_n15Hits6StsAllRings = new TH1D("fh_n15Hits6StsAllRings","Number of RICH rings. >=6 STS Points. >= 15hits",50,0,150);

  fh_n6StsERings = new TH1D("fh_n6StsERings","Number of all electron RICH rings. >=6 STS Points",50,0,150);
  fh_n5Hits6StsERings = new TH1D("fh_n5Hits6StsERings","Number of electron RICH rings. >=6 STS Points. >= 5 hits",50,0,150);
  fh_n10Hits6StsERings = new TH1D("fh_n10Hits6StsERings","Number of electron RICH rings. >=6 STS Points. >= 10 hits",50,0,150);
  fh_n15Hits6StsERings = new TH1D("fh_n15Hits6StsERings","Number of electron RICH rings. >=6 STS Points. >= 15hits",50,0,150);

  fh_n6StsPiRings = new TH1D("fh_n6StsPiRings","Number of all pion RICH rings. >=6 STS Points",50,0,150);
  fh_n5Hits6StsPiRings = new TH1D("fh_n5Hits6StsPiRings","Number of pion RICH rings. >=6 STS Points. >= 5 hits",50,0,150);
  fh_n10Hits6StsPiRings = new TH1D("fh_n10Hits6StsPiRings","Number of pion RICH rings. >=6 STS Points. >= 10 hits",50,0,150);
  fh_n15Hits6StsPiRings = new TH1D("fh_n15Hits6StsPiRings","Number of pion RICH rings. >=6 STS Points. >= 15hits",50,0,150);


  fh_n6StsProjAllRings = new TH1D("fh_n6StsProjAllRings","Number of All RICH rings. >=6 STS Points. Proj.",50,0,150);
  fh_n5Hits6StsProjAllRings = new TH1D("fh_n5Hits6StsProjAllRings","Number of RICH rings. >=6 STS Points. >= 5 hits. Proj.",50,0,150);
  fh_n10Hits6StsProjAllRings = new TH1D("fh_n10Hits6StsProjAllRings","Number of RICH rings. >=6 STS Points. >= 10 hits. Proj.",50,0,150);
  fh_n15Hits6StsProjAllRings = new TH1D("fh_n15Hits6StsProjAllRings","Number of RICH rings. >=6 STS Points. >= 15hits. Proj.",50,0,150);

  fh_n6StsProjERings = new TH1D("fh_n6StsProjERings","Number of all electron RICH rings. >=6 STS Points. Proj.",50,0,150);
  fh_n5Hits6StsProjERings = new TH1D("fh_n5Hits6StsProjERings","Number of electron RICH rings. >=6 STS Points. >= 5 hits. Proj.",50,0,150);
  fh_n10Hits6StsProjERings = new TH1D("fh_n10Hits6StsProjERings","Number of electron RICH rings. >=6 STS Points. >= 10 hits. Proj.",50,0,150);
  fh_n15Hits6StsProjERings = new TH1D("fh_n15Hits6StsProjERings","Number of electron RICH rings. >=6 STS Points. >= 15hits. Proj.",50,0,150);

  fh_n6StsProjPiRings = new TH1D("fh_n6StsProjPiRings","Number of all pion RICH rings. >=6 STS Points. Proj.",50,0,150);
  fh_n5Hits6StsProjPiRings = new TH1D("fh_n5Hits6StsProjPiRings","Number of pion RICH rings. >=6 STS Points. >= 5 hits. Proj.",50,0,150);
  fh_n10Hits6StsProjPiRings = new TH1D("fh_n10Hits6StsProjPiRings","Number of pion RICH rings. >=6 STS Points. >= 10 hits. Proj.",50,0,150);
  fh_n15Hits6StsProjPiRings = new TH1D("fh_n15Hits6StsProjPiRings","Number of pion RICH rings. >=6 STS Points. >= 15hits. Proj.",50,0,150);

  fh_percOfTrueHitsAll = new TH1D("fh_percOfTrueHitsAll","Percent of true hits.",50,0,1.01);
  fh_percOfTrueMCHitsAll = new TH1D("fh_percOfTrueMCHitsAll","Percent of true hits.",50,0,1.01);

  fh_percOfTrueHitsE = new TH1D("fh_percOfTrueHitsE","Percent of true hits. Electrons. Proj",50,0,1.01);
  fh_percOfTrueMCHitsE = new TH1D("fh_percOfTrueMCHitsE","Percent of true hits. Electrons. Proj",50,0,1.01);

  fh_percOfTrueHitsPi = new TH1D("fh_percOfTrueHitsPi","Percent of true hits. Pions. Proj",50,0,1.01);
  fh_percOfTrueMCHitsPi = new TH1D("fh_percOfTrueMCHitsPi","Percent of true hits. Pions. Proj",50,0,1.01);


  fh_True10HitsFoundRingsMomAll = new TH1D("fh_True10HitsFoundRingsMomAll","Number of rings vs momentum.>=10 hits.Proj",40,0,15);
  fh_10HitsMCRingsMomAll = new TH1D("fh_10HitsMCRingsMomAll","Number of rings vs momentum..>=10 hits.Proj",40,0,15);

  fh_True10HitsFoundRingsMomE = new TH1D("fh_True10HitsFoundRingsMomE","Number of rings vs momentum. Electrons. Proj.>=10 hits",40,0,10);
  fh_10HitsMCRingsMomE = new TH1D("fh_10HitsMCRingsMomE","Number of rings vs momentum. Electrons. Proj.>=10 hits",40,0,10);

  fh_True10HitsFoundRingsMomPi = new TH1D("fh_True10HitsFoundRingsMomPi","Number of rings vs momentum. Pions. Proj.>=10 hits",40,0,15);
  fh_10HitsMCRingsMomPi = new TH1D("fh_10HitsMCRingsMomPi","Number of rings vs momentum. Pions. Proj.>=10 hits",40,0,15);

  fh_CloneMom = new TH1D("fh_CloneMom","Number of clone rings vs momentum",40,0,15);

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

  fh_Fake7RingsNhits = new TH1D("fh_Fake7RingsNhits","Number of hits in found rings",50,0,50);
  fh_Fake8RingsNhits = new TH1D("fh_Fake8RingsNhits","Number of hits in found rings",50,0,50);
  fh_TrueRingsNhitsE = new TH1D("fh_TrueRingsNhitsE","Number of hits in found rings",50,0,50);
  fh_TrueRingsNhitsPi = new TH1D("fh_TrueRingsNhitsPi","Number of hits in found rings",50,0,50);


  fh_Fake7RingTrackDist = new TH1D("fh_Fake7RingTrackDist","Distance between track and ring center",50,0,5);
  fh_Fake8RingTrackDist = new TH1D("fh_Fake8RingTrackDist","Distance between track and ring center",50,0,5);
  fh_TrueRingTrackDistE = new TH1D("fh_TrueRingTrackDistE","Distance between track and ring center",50,0,5);
  fh_TrueRingTrackDistPi = new TH1D("fh_TrueRingTrackDistPi","Distance between track and ring center",50,0,5);

  fh_Fake7Angle = new TH1D("fh_Fake7Angle","Biggest angle in found ring",50,0,6.5);
  fh_Fake8Angle = new TH1D("fh_Fake8Angle","Biggest angle in found ring",50,0,6.5);
  fh_TrueAngleE = new TH1D("fh_TrueAngleE","Biggest angle in found ring",50,0,6.5);
  fh_TrueAnglePi = new TH1D("fh_TrueAnglePi","Biggest angle in found ring",50,0,6.5);

  fh_Fake7RingsTBsum = new TH1D("fh_Fake7RingsTBsum","TBsum in found rings",50,0,50);
  fh_TrueRingsTBsumPi = new TH1D("fh_TrueRingsTBsumPi","TBsum in found rings",50,0,50);
  fh_TrueRingsTBsumE = new TH1D("fh_TrueRingsTBsumE","TBsum in found rings",50,0,50);
  fh_Fake8RingsTBsum = new TH1D("fh_Fake8RingsTBsum","TBsum in found rings",50,0,50);

  fh_Fake7RingsTBsumNhits = new TH1D("fh_Fake7RingsTBsumNhits","TBsum/Nhits in found rings",50,0,1);
  fh_TrueRingsTBsumNhitsPi = new TH1D("fh_TrueRingsTBsumNhitsPi","TBsum/Nhits in found rings",50,0,1);
  fh_TrueRingsTBsumNhitsE = new TH1D("fh_TrueRingsTBsumNhitsE","TBsum/Nhits in found rings",50,0,1);
  fh_Fake8RingsTBsumNhits = new TH1D("fh_Fake8RingsTBsumNhits","TBsum/Nhits in found rings",50,0,1);

  fh_Fake7RingsChi2 = new TH1D("fh_Fake7RingsChi2","Chi2 in found rings",50,0,1.5);
  fh_TrueRingsChi2Pi = new TH1D("fh_TrueRingsChi2Pi","Chi2 in found rings",50,0,1.5);
  fh_TrueRingsChi2E = new TH1D("fh_TrueRingsChi2E","Chi2 in found rings",50,0,1.5);
  fh_Fake8RingsChi2 = new TH1D("fh_Fake8RingsChi2","Chi2 in found rings",50,0,1.5);

  fh_RvsMomTrue =new TH2D("fh_RvsMomTrue","momentum vs radius true ring",40,0,10,50,0,7);
  fh_RvsMomFake = new TH2D("fh_RvsMomFake","momentum vs radius fake ring",40,0,10,50,0,7);

  //hits distribution (x,y)
  fh_HitsXY = new TH2D("fh_HitsXY","Hits distribution (x,y)",100,-200,200,125,-250,250);
  // (x,y) of true primary electron rings
  fh_TrueFoundRingsXYE = new TH2D("fh_TrueFoundRingsXYE","(x,y) true electrons rings. >=10 Hits. Proj",100,-200,200,125,-250,250);
  // (x,y) of true primary pions rings
  fh_TrueFoundRingsXYPi = new TH2D("fh_TrueFoundRingsXYPi","(x,y) true pions rings. >=10 Hits. Proj",100,-200,200,125,-250,250);
  //number of hits per event
  fh_NhitsPerEvent = new TH1D("fh_NhitsPerEvent","Number of hits per event",100,1000,3000);
  //number of projections per event
  fh_NprojPerEvent = new TH1D("fh_NprojPerEvent","Number of projections per event",50,200,450);

  fh_NClonesPerEvent = new TH1D("fh_NClonesPerEvent","Number of clone rings per event",30,0,30);
  fh_NFakesPerEvent = new TH1D("fh_NFakesPerEvent","Number of fake rings per event",50,0,50);

  fh_MomSlice1GeVE = new TH1D("fh_MomSlice1GeVE","Slice on p=1+-0.2 GeV",50,0,7);
  fh_MomSlice7GeVE = new TH1D("fh_MomSlice7GeVE","Slice on p=7+-0.2 GeV",50,0,7);
  fh_MomSlice9GeVE = new TH1D("fh_MomSlice9GeVE","Slice on p=9+-0.2 GeV",50,0,7);

  fh_MomSlice1GeVPi = new TH1D("fh_MomSlice1GeVPi","Slice on p=1+-0.2 GeV",50,0,7);
  fh_MomSlice7GeVPi = new TH1D("fh_MomSlice7GeVPi","Slice on p=7+-0.2 GeV",50,0,7);
  fh_MomSlice9GeVPi = new TH1D("fh_MomSlice9GeVPi","Slice on p=9+-0.2 GeV",50,0,7);

}

  // -----   Destructor   ----------------------------------------------------
CbmRichRingQa::~CbmRichRingQa()
{


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


// -----   Initialization   -----------------------------------------------

InitStatus CbmRichRingQa::Init()
{
 // Get and check CbmRootManager
  CbmRootManager* ioman = CbmRootManager::Instance();
  if (! ioman) {
    cout << "-E- CbmRichRingQa::Init: "
	 << "RootManager not instantised!" << endl;
    return kERROR;
  }

  // Get hit Array
  fHits = (TClonesArray*) ioman->GetObject("RichHit");
  if ( ! fHits) {
    cout << "-W- CbmRichRingQa::Init: No RichHit array!"
	 << endl;
  }

  // Get RichRing Array
  fRings = (TClonesArray*) ioman->GetObject("RichRing");
  if ( ! fRings ) {
    cout << "-E- CbmRichRingQa::Init: No RichRing array!" << endl;
    return kERROR;
  }

  // Get MC Point array
  fPoints = (TClonesArray*) ioman->GetObject("RichPoint");
  if ( ! fPoints ) {
    cout << "-E- CbmRichRingQa::Init: No RichPoint array!" << endl;
    return kERROR;
  }

  // Get MC Point array
  fTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  if ( ! fTracks ) {
    cout << "-E- CbmRichRingQa::Init: No MCTrack array!" << endl;
    return kERROR;
  }

  // Get RichRingMatch array
  fMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
  if ( ! fMatches ) {
    cout << "-E- CbmRichRingQa::Init: No RichRingMatch array!" << endl;
    return kERROR;
  }

  // Get RichProjection array
  fProj = (TClonesArray*) ioman->GetObject("RichProjection");
  if ( ! fProj ) {
    cout << "-E- CbmRichRingQa::Init: No RichProjection array!" << endl;
    return kERROR;
  }


  // get TrackMatch array
  fTrackMatch = (TClonesArray*) ioman->GetObject("STSTrackMatch");
  if ( ! fTrackMatch) {
    cout << "-E- CbmRichRingQa::Init: No track match array!"<< endl;
    return kERROR;
  }
  
  // Get global track array
    gTrackArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
    if ( ! gTrackArray) {
        cout << "-W- CbmRichRingQa::Init: No global track array!"
	    << endl;
	return kERROR;
     }
     
  return kSUCCESS;

}

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

// -----   Execution of Task   ---------------------------------------------

// -------------------------------------------------------------------------
void CbmRichRingQa::Exec(Option_t* option)
{
  cout<<"CbmRichRingQa Event No. "<< nevents++ << endl;
  // Create some pointers and variables
  CbmTrackParam* proj  = NULL;
  CbmRichHit* hit    = NULL;
  CbmMCPoint* point  = NULL;
  CbmMCTrack* track  = NULL;
  CbmStsTrackMatch* trackMatch = NULL;
  CbmGlobalTrack* gtrack = NULL;

  Int_t iPoint       = 0;
  Int_t iMCTrack     = 0;
  Int_t iMother      = 0;

  map<Int_t, Int_t>::iterator it;//iterator for acass to fRingMap

  //calculate number of hits in MC ring to
  //compare its with number of hits in found ring
  fRingMap.clear();
  Int_t nRichHits = fHits->GetEntriesFast();
  fh_NhitsPerEvent->Fill(nRichHits);
  fh_NprojPerEvent->Fill(fProj->GetEntriesFast());
  // Loop over Rich hits
  for (Int_t iHit=0; iHit < nRichHits; iHit++) {
      hit = (CbmRichHit*) fHits->At(iHit);
      if ( ! hit ) {
	cout << "-E- CbmRichMatchRings::Exec: "
	     << "No Hit " << iHit << endl;
	continue;
      }
      fh_HitsXY->Fill(hit->X(), hit->Y());
      iPoint = hit->GetRefIndex();

      if ( iPoint < 0 )continue; // Fake or background hit



      //Get the MC Point corresponding to the hit

      point = (CbmMCPoint*) fPoints->At(iPoint);
      if ( ! point ) {
	cout << "-E- CbmRichRingQa::Exec() "
	     << "Empty MCPoint " << iPoint << " from Hit " << iHit
	     << endl;
	continue;
      }
      //Get the MC Track ID corresponding to the MC Point
      iMCTrack = point->GetTrackID();
     // Get the MC Track corresponding to the ID
     track   = (CbmMCTrack*)fTracks->At(iMCTrack);
     if ( ! track ) {
	cout << "-E- CbmRichRingQa::Exec "
	     << "Empty MCtrack " << iMCTrack
	     << endl;
	continue;
      }
     iMother = track->GetMotherID();
     if (iMother == -1) {
       cout << "-E- MotherID of Cherenkov photon = -1" << endl;
       continue;
     }
     fRingMap[iMother]++;
  }  // Loop over Rich hits

  Int_t nAllRings = 0;      //number of all MC rings
  Int_t n5HitsAllRings = 0; //number of MC rings with 5 or more hits
  Int_t n10HitsAllRings = 0;//number of MC rings with 10 or more hits
  Int_t n15HitsAllRings = 0;//number of MC rings with 15 or more hits

  Int_t nERings = 0;  //number of all electron MC rings
  Int_t n5HitsERings = 0;//number of electron MC rings with 5 or more hits
  Int_t n10HitsERings = 0;//number of electron MC rings with 10 or more hits
  Int_t n15HitsERings = 0;//number of electron MC rings with 15 or more hits

  Int_t nPiRings = 0;  //number of all pion MC rings
  Int_t n5HitsPiRings = 0;//number of pion MC rings with 5 or more hits
  Int_t n10HitsPiRings = 0;//number of pion MC rings with 10 or more hits
  Int_t n15HitsPiRings = 0;//number of pion MC rings with 15 or more hits
//------------------------------------------------------------------------------------------------------------
  Int_t n6StsAllRings = 0;      //number of all MC rings
  Int_t n5Hits6StsAllRings = 0; //number of MC rings with 5 or more hits
  Int_t n10Hits6StsAllRings = 0;//number of MC rings with 10 or more hits
  Int_t n15Hits6StsAllRings = 0;//number of MC rings with 15 or more hits

  Int_t n6StsERings = 0;  //number of all electron MC rings
  Int_t n5Hits6StsERings = 0;//number of electron MC rings with 5 or more hits
  Int_t n10Hits6StsERings = 0;//number of electron MC rings with 10 or more hits
  Int_t n15Hits6StsERings = 0;//number of electron MC rings with 15 or more hits

  Int_t n6StsPiRings = 0;  //number of all pion MC rings
  Int_t n5Hits6StsPiRings = 0;//number of pion MC rings with 5 or more hits
  Int_t n10Hits6StsPiRings = 0;//number of pion MC rings with 10 or more hits
  Int_t n15Hits6StsPiRings = 0;//number of pion MC rings with 15 or more hits
//------------------------------------------------------------------------------------------------------------
  Int_t n6StsProjAllRings = 0;      //number of all MC ringsand projection
  Int_t n5Hits6StsProjAllRings = 0; //number of MC rings with 5 or more hitsand projection
  Int_t n10Hits6StsProjAllRings = 0;//number of MC rings with 10 or more hitsand projection
  Int_t n15Hits6StsProjAllRings = 0;//number of MC rings with 15 or more hitsand projection

  Int_t n6StsProjERings = 0;  //number of all electron MC ringsand projection
  Int_t n5Hits6StsProjERings = 0;//number of electron MC rings with 5 or more hitsand projection
  Int_t n10Hits6StsProjERings = 0;//number of electron MC rings with 10 or more hitsand projection
  Int_t n15Hits6StsProjERings = 0;//number of electron MC rings with 15 or more hitsand projection

  Int_t n6StsProjPiRings = 0;  //number of all pion MC ringsand projection
  Int_t n5Hits6StsProjPiRings = 0;//number of pion MC rings with 5 or more hitsand projection
  Int_t n10Hits6StsProjPiRings = 0;//number of pion MC rings with 10 or more hitsand projection
  Int_t n15Hits6StsProjPiRings = 0;//number of pion MC rings with 15 or more hitsand projection

  Int_t nProj = fProj->GetEntriesFast();
  //Loop for all MC rings
  for (it=fRingMap.begin(); it!=fRingMap.end(); it++) {
    track = (CbmMCTrack*) fTracks->At(it->first);
    if (!track){
       cout << "-E- CbmRichRingQa::Exec() no track"<<
       it->first<<endl;
       continue;

    }
    Int_t nStsPoints = track->GetStsPoints();
    double momentum = track->GetMomentum().Mag();
    Int_t trackID = it->first;//track->GetTrackID();
    bool isProj = false;
    //search for projection with such TrackID
    for (Int_t iProj = 0; iProj < nProj; iProj++){
      proj   = (CbmTrackParam*)fProj->At(iProj);
      if (!proj){
        cout << "-E- CbmRichRingQa::Exec() no projection"<<
        iProj<<endl;
        continue;

      }
      gtrack = (CbmGlobalTrack*)gTrackArray->At(iProj);
      trackMatch = (CbmStsTrackMatch*)fTrackMatch->At(gtrack->GetStsTrackIndex());
      if (!trackMatch) cout << "no matched track!: fake?"<< endl;
      if (!trackMatch) continue;;
      if (trackMatch->GetMCTrackId() == trackID){
        isProj = true;
	break;
      }
    }

    nAllRings++;
    if (it->second >=5) n5HitsAllRings++;
    if (it->second >=10) n10HitsAllRings++;
    if (it->second >=15) n15HitsAllRings++;

    if (nStsPoints >=6){ ///number of points in STS more than 6
      n6StsAllRings++;
      if (it->second >=5) n5Hits6StsAllRings++;
      if (it->second >=10) n10Hits6StsAllRings++;
      if (it->second >=15) n15Hits6StsAllRings++;
      if (isProj){//if we have projection on photodetector plane
        n6StsProjAllRings++;
        if (it->second >=5) n5Hits6StsProjAllRings++;
        if (it->second >=10) n10Hits6StsProjAllRings++;
        if (it->second >=15) n15Hits6StsProjAllRings++;
        if (it->second >=10) fh_10HitsMCRingsMomAll->Fill(momentum);
      }
    }


    Int_t gcode = track->GetPdgCode();

    if ((TMath::Abs(gcode) == 11)) {//electron rings
      nERings++;
      if (it->second >=5) n5HitsERings++;
      if (it->second >=10) n10HitsERings++;
      if (it->second >=15) n15HitsERings++;
      if (nStsPoints >=6){ ///number of points in STS more than 6
        n6StsERings++;
        if (it->second >=5) n5Hits6StsERings++;
        if (it->second >=10) n10Hits6StsERings++;
        if (it->second >=15) n15Hits6StsERings++;
        if (isProj){//if we have projection on photodetector plane
	  n6StsProjERings++;
          if (it->second >=5) n5Hits6StsProjERings++;
          if (it->second >=10) n10Hits6StsProjERings++;
          if (it->second >=15) n15Hits6StsProjERings++;
	  if (it->second >=10) fh_10HitsMCRingsMomE->Fill(momentum);

	}
      }
    }

    if ((TMath::Abs(gcode) == 211)) {//pion rings
      nPiRings++;
      if (it->second >=5) n5HitsPiRings++;
      if (it->second >=10) n10HitsPiRings++;
      if (it->second >=15) n15HitsPiRings++;
      if (nStsPoints >=6){ ///number of points in STS more than 6
        n6StsPiRings++;
        if (it->second >=5) n5Hits6StsPiRings++;
        if (it->second >=10) n10Hits6StsPiRings++;
        if (it->second >=15) n15Hits6StsPiRings++;
        if (isProj){//if we have projection on photodetector plane
	  n6StsProjPiRings++;
          if (it->second >=5) n5Hits6StsProjPiRings++;
          if (it->second >=10) n10Hits6StsProjPiRings++;
          if (it->second >=15) n15Hits6StsProjPiRings++;
	  if (it->second >=10) fh_10HitsMCRingsMomPi->Fill(momentum);
	}
      }
    }


  }  //Loop for all MC rings

  // Fill the histograms
  fh_nAllRings->Fill(nAllRings);
  fh_n5HitsAllRings->Fill(n5HitsAllRings);
  fh_n10HitsAllRings->Fill(n10HitsAllRings);
  fh_n15HitsAllRings->Fill(n15HitsAllRings);

  fh_nERings->Fill(nERings);
  fh_n5HitsERings->Fill(n5HitsERings);
  fh_n10HitsERings->Fill(n10HitsERings);
  fh_n15HitsERings->Fill(n15HitsERings);

  fh_nPiRings->Fill(nPiRings);
  fh_n5HitsPiRings->Fill(n5HitsPiRings);
  fh_n10HitsPiRings->Fill(n10HitsPiRings);
  fh_n15HitsPiRings->Fill(n15HitsPiRings);
  //------------------------------------------------------------------------------------------------------------
  fh_n6StsAllRings->Fill(n6StsAllRings);
  fh_n5Hits6StsAllRings->Fill(n5Hits6StsAllRings);
  fh_n10Hits6StsAllRings->Fill(n10Hits6StsAllRings);
  fh_n15Hits6StsAllRings->Fill(n15Hits6StsAllRings);

  fh_n6StsERings->Fill(n6StsERings);
  fh_n5Hits6StsERings->Fill(n5Hits6StsERings);
  fh_n10Hits6StsERings->Fill(n10Hits6StsERings);
  fh_n15Hits6StsERings->Fill(n15Hits6StsERings);

  fh_n6StsPiRings->Fill(n6StsPiRings);
  fh_n5Hits6StsPiRings->Fill(n5Hits6StsPiRings);
  fh_n10Hits6StsPiRings->Fill(n10Hits6StsPiRings);
  fh_n15Hits6StsPiRings->Fill(n15Hits6StsPiRings);
//------------------------------------------------------------------------------------------------------------
  fh_n6StsProjAllRings->Fill(n6StsProjAllRings);
  fh_n5Hits6StsProjAllRings->Fill(n5Hits6StsProjAllRings);
  fh_n10Hits6StsProjAllRings->Fill(n10Hits6StsProjAllRings);
  fh_n15Hits6StsProjAllRings->Fill(n15Hits6StsProjAllRings);

  fh_n6StsProjERings->Fill(n6StsProjERings);
  fh_n5Hits6StsProjERings->Fill(n5Hits6StsProjERings);
  fh_n10Hits6StsProjERings->Fill(n10Hits6StsProjERings);
  fh_n15Hits6StsProjERings->Fill(n15Hits6StsProjERings);

  fh_n6StsProjPiRings->Fill(n6StsProjPiRings);
  fh_n5Hits6StsProjPiRings->Fill(n5Hits6StsProjPiRings);
  fh_n10Hits6StsProjPiRings->Fill(n10Hits6StsProjPiRings);
  fh_n15Hits6StsProjPiRings->Fill(n15Hits6StsProjPiRings);


  EfficiencyCalc();
  DiffFakeTrue();

}




void CbmRichRingQa::EfficiencyCalc()
{

  CbmRichRingMatch *match;
  CbmMCTrack* track;
  CbmTrackParam* proj;
  CbmRichRing* ring;
  CbmStsTrackMatch* trackMatch;
  CbmGlobalTrack* gtrack;
  Int_t NFakesCount = 0;
  Int_t NClonesCount = 0;

  Int_t nMatches = fMatches->GetEntriesFast();
  cout<<"Number of found rings = "<<nMatches<<endl;

  vector <int> clone;
  clone.clear();

  for (Int_t iMatches = 0; iMatches < nMatches; iMatches++){
    match   = (CbmRichRingMatch*)fMatches->At(iMatches);
    if (!match){
       cout << "-E- CbmRichRingQa::EfficiencyCalc() no match"<<
       iMatches<<endl;
       continue;
    }
    ring = (CbmRichRing*)fRings->At(iMatches);
    if (!ring){
       cout << "-E- CbmRichRingQa::EfficiencyCalc() no ring"<<
       iMatches<<endl;
       continue;
    }
    double radius = ring->GetRadius();
    int lTrueHits = match->GetNofTrueHits();
    int lWrongHits = match->GetNofWrongHits();
    int lFakeHits = match->GetNofFakeHits();
    int lMCHits = match->GetNofMCHits();
    int lFoundHits = lTrueHits + lWrongHits + lFakeHits;
    double lPercTrue = (Double_t)lTrueHits / (Double_t)lFoundHits;
    double lPercMCTrue = (Double_t)lTrueHits / (Double_t)lMCHits;
    //cout <<"lPercTrue = "<< lPercTrue<<"  lPercMCTrue = "<<lPercMCTrue <<endl;
    fh_percOfTrueHitsAll->Fill(lPercTrue);
    fh_percOfTrueMCHitsAll->Fill(lPercMCTrue);

    Int_t trackID = match->GetMCTrackID();
    track = (CbmMCTrack*) fTracks->At(trackID);
    if (!track) {
      cout << "-E- CbmRichRingQa::EfficiencyCalc(). No track "
      << trackID <<endl;
      continue;
    }
    Int_t gcode = track->GetPdgCode();
    double momentum = track->GetMomentum().Mag();
    Int_t nProj = fProj->GetEntriesFast();
    bool isProj = false;
    for (Int_t iProj = 0; iProj < nProj; iProj++){
      proj   = (CbmTrackParam*)fProj->At(iProj);
      if (!proj) {
        cout << "-E- CbmRichRingQa::EfficiencyCalc(). No projection "
	<< iProj <<endl;
        continue;
      }
      gtrack = (CbmGlobalTrack*)gTrackArray->At(iProj);
      trackMatch = (CbmStsTrackMatch*)fTrackMatch->At(gtrack->GetStsTrackIndex());
      if (!trackMatch){
        cout << "-E- CbmRichRingQa::EfficiencyCalc(). no matched track!: fake?"
	<< endl;
        continue;
      }

      if (trackMatch->GetMCTrackId() == trackID){
        isProj = true;
	break;
      }
    }

    //choose the normalisation type!!!!!!!
    double PercOfTrueHits = -1;
    if (fNormType == 0) //normalize by Number of MC hits
      PercOfTrueHits = lPercMCTrue;
    if (fNormType == 1) //normalize by Number of found hits
      PercOfTrueHits = lPercTrue;

    if (PercOfTrueHits == -1){
      cout << "-E- CbmRichRingQa::EfficiencyCalc() "<<
        "PercOfTrueHits == -1"<< endl;
      continue;
    }


    ring->SetRecFlag(-1);
    //fake rings
    if (PercOfTrueHits < 0.7){
      ring->SetRecFlag(7);
      if (radius >=5 && radius <=7) ring->SetRecFlag(8);
    }else{//true rings

      if ((TMath::Abs(gcode) == 11)){ //electron rings
        ring->SetRecFlag(1);
        if (lMCHits >= 10){//more than 10 hits
	  ring->SetRecFlag(2);
	  if (isProj) ring->SetRecFlag(3);
	}
      }
      if ((TMath::Abs(gcode) == 211)){ //pion rings
        ring->SetRecFlag(4);
        if (lMCHits >= 10){//more than 10 hits
	  ring->SetRecFlag(5);
	  if (isProj) ring->SetRecFlag(6);
	}
      }
    }//true rings
    for (unsigned int ci = 0; ci < clone.size(); ci++)
      if (trackID == clone[ci]) ring->SetRecFlag(9);
    clone.push_back(trackID);

    int recFlag = ring->GetRecFlag();


    //clone ring
    if (recFlag ==9){
      NClonesCount++;
      fh_CloneMom->Fill(momentum);
    }
    //fake ring
    if (recFlag == 7){
      fh_FakeFoundRingsXYAll->Fill(ring->GetCenterX(), ring->GetCenterY());
      fh_Fake7RingsNhits->Fill(lFoundHits);
      NFakesCount++;
    }
    //fake E Ring
    if (recFlag == 8){
      fh_Fake8RingsNhits->Fill(lFoundHits);
      fh_FakeFoundRingsXYAll->Fill(ring->GetCenterX(), ring->GetCenterY());
      NFakesCount++;
    }

    if (recFlag == 3){//electrons. >=10 hits. Projection
      fh_True10HitsFoundRingsMomE->Fill(momentum);
      fh_TrueRingsNhitsE->Fill(lFoundHits);
      fh_TrueFoundRingsXYE->Fill(ring->GetCenterX(),ring->GetCenterY());
    }

    if (recFlag == 6){//pions. >=10 hits. Projection
      fh_True10HitsFoundRingsMomPi->Fill(momentum);
      fh_TrueFoundRingsXYPi->Fill(ring->GetCenterX(),ring->GetCenterY());
      fh_TrueRingsNhitsPi->Fill(lFoundHits);
    }

    if (recFlag >=1 && recFlag <=6)// all true El and Pi rings
      fh_RvsMomTrue->Fill(momentum,ring->GetRadius());
    if (recFlag == 7 || recFlag == 8)
      fh_RvsMomFake->Fill(momentum,ring->GetRadius());

      //momentum slice on p=1; 7; 9 GeV; +-0.2GeV
    if (recFlag ==3){//  true El
      if (momentum>=0.8 && momentum <=1.2)
        fh_MomSlice1GeVE->Fill(ring->GetRadius());
      if (momentum>=6.8 && momentum <=7.2)
        fh_MomSlice7GeVE->Fill(ring->GetRadius());
      if (momentum>=8.8 && momentum <=9.2)
       fh_MomSlice9GeVE->Fill(ring->GetRadius());
    }
    if (recFlag ==6){//  true Pi rings
      if (momentum>=0.8 && momentum <=1.2)
        fh_MomSlice1GeVPi->Fill(ring->GetRadius());
      if (momentum>=6.8 && momentum <=7.2)
        fh_MomSlice7GeVPi->Fill(ring->GetRadius());
      if (momentum>=8.8 && momentum <=9.2)
       fh_MomSlice9GeVPi->Fill(ring->GetRadius());
    }


    if (isProj){
      if (PercOfTrueHits > 0.7 && lMCHits >= 10)
        fh_True10HitsFoundRingsMomAll->Fill(momentum);
      if ((TMath::Abs(gcode) == 11)){
        fh_percOfTrueHitsE->Fill(lPercTrue);
        fh_percOfTrueMCHitsE->Fill(lPercMCTrue);
      }
      if ((TMath::Abs(gcode) == 211)){
        fh_percOfTrueHitsPi->Fill(lPercTrue);
        fh_percOfTrueMCHitsPi->Fill(lPercMCTrue);

      }
    }
  }
  fh_NClonesPerEvent->Fill(NClonesCount);
  fh_NFakesPerEvent->Fill(NFakesCount);

}

void CbmRichRingQa::DiffFakeTrue()
{
  CbmRichRingMatch *match;
  CbmTrackParam* proj;
  CbmRichRing* ring;
  CbmRichHit* hitRing;
  double xRing, yRing, rMin;
  Int_t nMatches = fMatches->GetEntriesFast();
  Int_t nProj = fProj->GetEntriesFast();
  for (Int_t iMatches = 0; iMatches < nMatches; iMatches++){
    match   = (CbmRichRingMatch*)fMatches->At(iMatches);
    if (!match){
       cout << "-E- CbmRichRingQa::DiffFakeTrue() no match"<<
       iMatches<<endl;
       continue;
    }
    ring = (CbmRichRing*)fRings->At(iMatches);
    if (!ring){
       cout << "-E- CbmRichRingQa::DiffFakeTrue() no ring"<<
       iMatches<<endl;
       continue;
    }
    int recFlag = ring->GetRecFlag();
    xRing = ring->GetCenterX();
    yRing = ring->GetCenterY();

    rMin = 999.;
    for (Int_t iProj=0; iProj < nProj; iProj++){
      proj = (CbmTrackParam*)fProj->At(iProj);
      if (!proj) {
        cout << "-E- CbmRichRingQa::DiffFakeTrue(). No projection "
	<< iProj <<endl;
        continue;
      }
      double xProj = proj->GetX();
      double yProj = proj->GetY();
      double distRingTrack = TMath::Sqrt( (xRing-xProj)*(xRing-xProj) +
				            (yRing-yProj)*(yRing-yProj) );
      if (distRingTrack < rMin)
	rMin = distRingTrack;
    } // loop projections

    ///biggest angle BEGIN
    vector<double> alpha;
    vector<double> phi;

    double Pi = 3.1415;
    Int_t nHits = ring->GetNoOfHits();
    for(int iHit = 0; iHit < nHits; iHit++){
      hitRing = (CbmRichHit*)ring->GetHit(iHit);
      if (!hitRing) {
        cout << "-E- CbmRichRingQa::DiffFakeTrue(). No Hit "<< iHit <<endl;
        continue;
      }
      if( hitRing->X() > xRing && hitRing->Y() > yRing ){
	alpha.push_back(atan(fabs((hitRing->X()-xRing)/(hitRing->Y()-yRing))));
      }
      if( hitRing->X() < xRing && hitRing->Y() > yRing ){
        alpha.push_back(Pi - atan(fabs((hitRing->X()-xRing)/(hitRing->Y()-yRing))));
      }
      if( hitRing->X() < xRing && hitRing->Y() < yRing ){
        alpha.push_back(Pi + atan(fabs((hitRing->X()-xRing)/(hitRing->Y()-yRing))));
      }
      if( hitRing->X() > xRing && hitRing->Y() < yRing ){
        alpha.push_back(2*Pi - atan(fabs((hitRing->X()-xRing)/(hitRing->Y()-yRing))));
      }
    }
    sort(alpha.begin(),alpha.end());
    for(int i = 0; i < nHits-1; i++) phi.push_back(alpha[i+1] - alpha[i]);
    phi.push_back(2*Pi - alpha[nHits-1] + alpha[0]);
    sort(phi.begin(),phi.end());
    double bigestAngle = phi[nHits-1]+phi[nHits-2]+phi[nHits-3];
    ///biggest angle END

    //TBsum BEGIN
    int histTB[20];
    for (int iH = 0; iH < 20; iH++) histTB[iH] = 0;
    for(int iH = 0; iH < nHits; iH++){
      hitRing = (CbmRichHit*)ring->GetHit(iH);
      double r = pow((hitRing->X() - ring->GetCenterX())*
                     (hitRing->X() - ring->GetCenterX()) +
	             (hitRing->Y() - ring->GetCenterY())*
		     (hitRing->Y() - ring->GetCenterY()), 0.5);

      if (r > 3.5 && r < 6.5) {
	double ir ;
	modf(r / (6.5 / 20.0), &ir);
	if (ir < 20 && ir > 0 ) histTB[ (int) ir ]++;
      }
    }
    int ii;
    int max = 0;
    for (int i = 0; i < 20; i++)
      if (histTB[i] >= max) {
	max = histTB[i];
	ii = i;
      }
    unsigned int maxsum = histTB[ii];
    if (ii+1 < 20) maxsum += histTB[ii+1];
    if (ii-1 > 0) maxsum += histTB[ii-1];
    //TBSum END

    //CHI 2 BEGIN
    double chi2 = 0;
    for(int iH = 0; iH < nHits; iH++){
      hitRing = (CbmRichHit*)ring->GetHit(iH);
      chi2 += pow( sqrt((ring->GetCenterX() - hitRing->X())*
                         (ring->GetCenterX() - hitRing->X()) +
		         (ring->GetCenterY() - hitRing->Y())*
			 (ring->GetCenterY() - hitRing->Y()) )
			 - ring->GetRadius(), 2 );
    }
    chi2 = chi2 / (nHits - 3);
   //CHI 2 END


    if (recFlag == 7) fh_Fake7RingTrackDist->Fill(rMin);
    if (recFlag == 8) fh_Fake8RingTrackDist->Fill(rMin);
    if (recFlag == 3)fh_TrueRingTrackDistE->Fill(rMin);
    if (recFlag == 6)fh_TrueRingTrackDistPi->Fill(rMin);

    if (recFlag == 7) fh_Fake7Angle->Fill(bigestAngle);
    if (recFlag == 8)  fh_Fake8Angle->Fill(bigestAngle);
    if (recFlag == 3)fh_TrueAngleE->Fill(bigestAngle);
    if (recFlag == 6)fh_TrueAnglePi->Fill(bigestAngle);

    if (recFlag == 7) fh_Fake7RingsTBsum->Fill(maxsum);
    if (recFlag == 8) fh_Fake8RingsTBsum->Fill(maxsum);
    if (recFlag == 3) fh_TrueRingsTBsumE->Fill(maxsum);
    if (recFlag == 6) fh_TrueRingsTBsumPi->Fill(maxsum);

    if (recFlag == 7) fh_Fake7RingsTBsumNhits->Fill((double)maxsum/(double)nHits);
    if (recFlag == 8) fh_Fake8RingsTBsumNhits->Fill((double)maxsum/(double)nHits);
    if (recFlag == 3) fh_TrueRingsTBsumNhitsE->Fill((double)maxsum/(double)nHits);
    if (recFlag == 6) fh_TrueRingsTBsumNhitsPi->Fill((double)maxsum/(double)nHits);

    if (recFlag == 7) fh_Fake7RingsChi2->Fill(chi2);
    if (recFlag == 8) fh_Fake8RingsChi2->Fill(chi2);
    if (recFlag == 3) fh_TrueRingsChi2E->Fill(chi2);
    if (recFlag == 6) fh_TrueRingsChi2Pi->Fill(chi2);


  }//loop over matches

}
// -----   Finish Task   ---------------------------------------------------
void CbmRichRingQa::Finish()
{

  fRings->Clear();
  fPoints->Clear();
  fTracks->Clear();
  fHits->Clear();
  fMatches->Clear();
  gTrackArray->Clear();
  
  TDirectory *current = gDirectory;
  TDirectory *rich = current->mkdir("RichRingQaHist");
  rich->cd();

  TDirectory *sub = rich->mkdir("NRings");
  sub->cd();
  fh_nAllRings->Write();
  fh_n5HitsAllRings->Write();
  fh_n10HitsAllRings->Write();
  fh_n15HitsAllRings->Write();

  fh_nERings->Write();
  fh_n5HitsERings->Write();
  fh_n10HitsERings->Write();
  fh_n15HitsERings->Write();

  fh_nPiRings->Write();
  fh_n5HitsPiRings->Write();
  fh_n10HitsPiRings->Write();
  fh_n15HitsPiRings->Write();

  sub = rich->mkdir("N6StsRings");
  sub->cd();

  fh_n6StsAllRings->Write();
  fh_n5Hits6StsAllRings->Write();
  fh_n10Hits6StsAllRings->Write();
  fh_n15Hits6StsAllRings->Write();

  fh_n6StsERings->Write();
  fh_n5Hits6StsERings->Write();
  fh_n10Hits6StsERings->Write();
  fh_n15Hits6StsERings->Write();

  fh_n6StsPiRings->Write();
  fh_n5Hits6StsPiRings->Write();
  fh_n10Hits6StsPiRings->Write();
  fh_n15Hits6StsPiRings->Write();
  
  sub = rich->mkdir("N6StsProjRings");
  sub->cd();

  fh_n6StsProjAllRings->Write();
  fh_n5Hits6StsProjAllRings->Write();
  fh_n10Hits6StsProjAllRings->Write();
  fh_n15Hits6StsProjAllRings->Write();

  fh_n6StsProjERings->Write();
  fh_n5Hits6StsProjERings->Write();
  fh_n10Hits6StsProjERings->Write();
  fh_n15Hits6StsProjERings->Write();

  fh_n6StsProjPiRings->Write();
  fh_n5Hits6StsProjPiRings->Write();
  fh_n10Hits6StsProjPiRings->Write();
  fh_n15Hits6StsProjPiRings->Write();

  sub = rich->mkdir("percOfTrueMCHits");
  sub->cd();

  fh_percOfTrueMCHitsAll->Write();
  fh_percOfTrueHitsAll->Write();

  fh_percOfTrueMCHitsE->Write();
  fh_percOfTrueHitsE->Write();

  fh_percOfTrueMCHitsPi->Write();
  fh_percOfTrueHitsPi->Write();
 
  sub = rich->mkdir("True10HitsFoundRingsMom");
  sub->cd();

  fh_True10HitsFoundRingsMomAll->Write();
  fh_10HitsMCRingsMomAll->Write();

  fh_True10HitsFoundRingsMomE->Write();
  fh_10HitsMCRingsMomE->Write();

  fh_True10HitsFoundRingsMomPi->Write();
  fh_10HitsMCRingsMomPi->Write();

  sub = rich->mkdir("RingsNhits");
  sub->cd();
  fh_Fake7RingsNhits->Write();
  fh_Fake8RingsNhits->Write();
  fh_TrueRingsNhitsE->Write();
  fh_TrueRingsNhitsPi->Write();
 
  sub = rich->mkdir("RingTrackDist");
  sub->cd();
  fh_Fake7RingTrackDist->Write();
  fh_Fake8RingTrackDist->Write();
  fh_TrueRingTrackDistE->Write();
  fh_TrueRingTrackDistPi->Write();
 
  sub = rich->mkdir("Angle");
  sub->cd();
  fh_Fake7Angle->Write();
  fh_Fake8Angle->Write();
  fh_TrueAngleE->Write();
  fh_TrueAnglePi->Write();

  sub = rich->mkdir("RvsMom");
  sub->cd();
  fh_RvsMomTrue->Write();
  fh_RvsMomFake->Write();

  sub = rich->mkdir("XY");
  sub->cd();
  fh_HitsXY->Write();
  fh_TrueFoundRingsXYE->Write();
  fh_TrueFoundRingsXYPi->Write();
  fh_FakeFoundRingsXYAll->Write();

  rich->cd();

  fh_NhitsPerEvent->Write();
  fh_NprojPerEvent->Write();

  sub = rich->mkdir("RingsTBsum");
  sub->cd();

  fh_Fake7RingsTBsum->Write();
  fh_TrueRingsTBsumPi->Write();
  fh_TrueRingsTBsumE->Write();
  fh_Fake8RingsTBsum->Write();

  fh_Fake7RingsTBsumNhits->Write();
  fh_TrueRingsTBsumNhitsPi->Write();
  fh_TrueRingsTBsumNhitsE->Write();
  fh_Fake8RingsTBsumNhits->Write();

  fh_Fake7RingsChi2->Write();
  fh_TrueRingsChi2Pi->Write();
  fh_TrueRingsChi2E->Write();
  fh_Fake8RingsChi2->Write();

  rich->cd();

  fh_CloneMom->Write();

  fh_NClonesPerEvent->Write();
  fh_NFakesPerEvent->Write();

  sub = rich->mkdir("MomSlices");
  sub->cd();

  fh_MomSlice1GeVE->Write();
  fh_MomSlice7GeVE->Write();
  fh_MomSlice9GeVE->Write();

  fh_MomSlice1GeVPi->Write();
  fh_MomSlice7GeVPi->Write();
  fh_MomSlice9GeVPi->Write();

  current->cd();
}

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

ClassImp(CbmRichRingQa)