//////////////////////////////////////////////////////////////////////////////
//
// $Id: $
//
//*-- Author  : S. Lebedev
//
//_HADES_CLASS_DESCRIPTION
//////////////////////////////////////////////////////////////////////////////
//
//  HRich700GeoTestRealQa
//
//
//////////////////////////////////////////////////////////////////////////////

#include "hrich700geotestrealqa.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hcategory.h"
#include "hevent.h"
#include "hgeantrich.h"
#include "hlinearcatiter.h"
#include "hmatrixcatiter.h"
#include "hparset.h"
#include "hspectrometer.h"
#include "richdef.h"
#include "hrich700drawhist.h"
#include "hrich700pmt.h"
#include "hrichhit.h"
#include "hrichcal.h"
#include "hrich700raw.h"
#include "hparticletool.h"
#include "hparticledef.h"
#include "hparticlecand.h"
#include "hparticleevtinfo.h"

#include "TCanvas.h"
#include "TH2D.h"
#include "TPad.h"
#include "TEllipse.h"
#include "TRandom.h"
#include "TLatex.h"

#include "hgeantkine.h"

#include "hrich700histmanager.h"
#include "hrich700digipar.h"
#include "hrich700utils.h"

#include <iostream>
#include <sstream>

using namespace std;


ClassImp(HRich700GeoTestRealQa)

    HRich700GeoTestRealQa::HRich700GeoTestRealQa():
fEventNum(0),
fOutputDir(""),
fOutputRootFile("")
{

}

HRich700GeoTestRealQa::~HRich700GeoTestRealQa()
{

}

Bool_t HRich700GeoTestRealQa::init()
{
    fCatRichCal = gHades->getCurrentEvent()->getCategory(catRichCal);
    if (NULL == fCatRichCal) {
       Error("init()", "Initialization of catRicCal category failed, returning...");
       return kFALSE;
    }

    fCatRichHit = gHades->getCurrentEvent()->getCategory(catRichHit);
    if (NULL == fCatRichHit) {
	   Error("init", "Initializatin of RICH hit category failed, returning...");
	   return kFALSE;
    }

    fCatRichRaw = gHades->getCurrentEvent()->getCategory(catRich700Raw);
    // if (NULL == fCatRichRaw) {
    //     Error("init", "Initializatin of RICH raw category failed, returning...");
    //     return kFALSE;
    // }

    fCatParticleCand = gHades->getCurrentEvent()->getCategory(catParticleCand);
    if (NULL == fCatParticleCand) {
	Error("init", "Initializatin of catParticleCand category failed, returning...");
	return kFALSE;
    }

    fDigiPar = (HRich700DigiPar*) gHades->getRuntimeDb()->getContainer("Rich700DigiPar");
    if(!fDigiPar) {
	Error("init", "Can not retrieve HRich700DigiPar");
        return kFALSE;
    }

    fCatEvtInfo = gHades->getCurrentEvent()->getCategory(catParticleEvtInfo);
    if (NULL == fCatEvtInfo) {
	Error("init", "Initializatin of catParticleEvtInfo category failed, returning...");
	return kFALSE;
    }

    initHist();
    return kTRUE;

}

void HRich700GeoTestRealQa::initHist()
{
    fHM = new HRich700HistManager();
    fHM->Create1<TH1D>("fhNofEvents", "fhNofEvents;# Events;Entries", 1, 0, 1);

    fHM->Create2<TH2D>("fhCalsColRow", "fhCalsColRow;Col index;Row index;# Cals/event", 194,-0.5, 193.5, 194, -0.5, 193.5);
    fHM->Create2<TH2D>("fhCalsXY", "fhCalsXY;X [mm];Y [mm];# Cals/event", 175,-700, 700, 175, -700, 700);
    fHM->Create1<TH1D>("fhCalsPmtId", "fhCalsPmtId;PMT ID;# Cals/event", 600,-0.5, 599.5);
    fHM->Create2<TH2D>("fhCalsPmtXPmtY", "fhCalsPmtXPmtY;PMT X;PMT Y;# Cals/event", 25,-0.5, 24.5, 25, -0.5, 24.5);
    fHM->Create1<TH1D>("fhCalsSectors", "fhCalsSectors;Sector;# Cals/event", 6,-0.5, 5.5);
    fHM->Create1<TH1D>("fhNofCalsPerPmt", "fhNofCalsPerPmt;# Cals/PMT/event;# PMT/event", 67,-0.5, 66.5);
    fHM->Create3<TH3D>("fhNofCalsPerPmtPmtXPmtY", "fhNofCalsPerPmtPmtXPmtY;PMT X;PMT Y;# Cals/PMT/event", 25,-0.5, 24.5, 25, -0.5, 24.5, 67,-0.5, 66.5);

    //blobs
    fHM->Create1<TH1D>("fhNofBlobsPerEvent", "fhNofBlobsPerEvent;Nof blobs per event;% of total events", 10,-0.5, 9.5);
    fHM->Create2<TH2D>("fhBlobsPmtXPmtY", "fhBlobsPmtXPmtY;PMT X;PMT Y;Nof blobs/PMT/event", 25,-0.5, 24.5, 25, -0.5, 24.5);
    fHM->Create1<TH1D>("fhBlobsMatchQa", "fhBlobsMatchQa;RICH Matching Qa;Entries/event", 60, 0, 60.);
    fHM->Create2<TH2D>("fhBlobsDThetaVsDPhi", "fhBlobsDThetaVsDPhi;dTheta [deg];dPhi [deg];Entries/event", 60, -30, 30, 60, -30, 30);
    fHM->Create2<TH2D>("fhNofCandsVsNofBlobs", "fhNofCandsVsNofBlobs;Nof cands (used) per event;Nof blobs per event;Entries/event", 20, -0.5, 39.5, 10,-0.5, 9.5);

    fHM->Create1<TH1D>("fhCalsTheta", "fhCalsTheta;Theta [deg];# Cals/event", 100, 0., 100.);
    fHM->Create1<TH1D>("fhCalsPhi", "fhCalsPhi;Phi [deg];# Cals/event", 360, 0., 360.);
    fHM->Create1<TH1D>("fhCalsSectorsInterp", "fhCalsSectorsInterp;Sector;# Cals/event", 6,-0.5, 5.5);

    //ring center XY
    fHM->Create2<TH2D>("fhRingXYMore18mm", "fhRingXYMore18mm;X [mm];Y [mm];Entries/event", 140,-700, 700, 140, -700, 700);
    fHM->Create2<TH2D>("fhRingXYLess18mm", "fhRingXYLess18mm;X [mm];Y [mm];Entries/event", 140,-700, 700, 140, -700, 700);
    fHM->Create2<TH2D>("fhRingXY", "fhRingXY;X [mm];Y [mm];Entries/event", 140,-700, 700, 140, -700, 700);
    fHM->Create2<TH2D>("fhRingXYCand", "fhRingXYCand;X [mm];Y [mm];Entries/event", 140,-700, 700, 140, -700, 700);
    fHM->Create2<TH2D>("fhRingXYTofId", "fhRingXYTofId;X [mm];Y [mm];Entries/event", 140,-700, 700, 140, -700, 700);


    fHM->Create2<TH2D>("fhMomVsBetaMore18mm", "fhMomVsBetaMore18mm;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaLess18mm", "fhMomVsBetaLess18mm;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaCand", "fhMomVsBetaCand;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaRichMatch<1", "fhMomVsBetaRichMatch<1;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaRichMatch<2", "fhMomVsBetaRichMatch<2;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaRichMatch<3", "fhMomVsBetaRichMatch<3;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);
    fHM->Create2<TH2D>("fhMomVsBetaRichMatch<4", "fhMomVsBetaRichMatch<4;p x q [MeV/c];#beta;Entries/event", 50, -1000, 1000, 50, 0., 1.2);

    // Radius
    fHM->Create2<TH2D>("fhThetaVsRadius", "fhThetaVsRadius;Theta [deg];Radius [mm];Entries/event", 15, 0, 90, 60, 0, 30);
    fHM->Create2<TH2D>("fhThetaVsRadiusCand", "fhThetaVsRadiusCand;Theta [deg];Radius [mm];Entries/event", 15, 0, 90, 60, 0, 30);
    fHM->Create2<TH2D>("fhThetaVsRadiusTofId", "fhThetaVsRadiusTofId;Theta [deg];Radius [mm];Entries/event", 15, 0, 90, 60, 0, 30);

    fHM->Create2<TH2D>("fhNofCalsInRingVsRadius", "fhNofCalsInRingVsRadius;Nof cals per ring;Radius [mm];Entries/event", 30, -.5, 29.5, 60, 0, 30);
    fHM->Create2<TH2D>("fhNofCalsInRingVsRadiusCand", "fhNofCalsInRingVsRadiusCand;Nof cals per ring;Radius [mm];Entries/event", 30, -.5, 29.5, 60, 0, 30);
    fHM->Create2<TH2D>("fhNofCalsInRingVsRadiusTofId", "fhNofCalsInRingVsRadiusTofId;Nof cals per ring;Radius [mm];Entries/event", 30, -.5, 29.5, 60, 0, 30);

    fHM->Create2<TH2D>("fhBetaVsRadiusSys0", "fhBetaVsRadiusSys0;#beta;Radius [mm];Entries/event", 50, 0., 1.2, 60, 0, 30);
    fHM->Create2<TH2D>("fhBetaVsRadiusSys1", "fhBetaVsRadiusSys1;#beta;Radius [mm];Entries/event", 50, 0., 1.2, 60, 0, 30);
    fHM->Create2<TH2D>("fhMomVsRadius", "fhMomVsRadius;p x q [MeV/c];Radius [mm];Entries/event", 50, -1000, 1000, 60, 0, 30);
    fHM->Create2<TH2D>("fhNofUsedFlagCandVsRadius", "fhNofUsedFlagCandVsRadius;Nof cands (used) per event;Radius [mm];Entries/event", 15, -0.5, 29.5, 60, 0, 30);
    fHM->Create2<TH2D>("fhRichMatchQaVsRadiusCand", "fhRichMatchQaVsRadiusCand;RICH Matching Qa;Radius [mm];Entries/event", 60, 0., 12., 60, 0, 30);
    fHM->Create2<TH2D>("fhRichMatchQaVsRadiusTofId", "fhRichMatchQaVsRadiusTofId;RICH Matching Qa;Radius [mm];Entries/event", 60, 0., 12., 60, 0, 30);

    // Nof cals per ring
    fHM->Create2<TH2D>("fhThetaVsNofCalsInRing", "fhThetaVsNofCalsInRing;Theta [deg];Nof cals per ring;Entries/event", 15, 0, 90, 30, -0.5, 29.5);
    fHM->Create2<TH2D>("fhThetaVsNofCalsInRingCand", "fhThetaVsNofCalsInRingCand;Theta [deg];Nof cals per ring;Entries/event", 15, 0, 90, 30, -0.5, 29.5);
    fHM->Create2<TH2D>("fhThetaVsNofCalsInRingTofId", "fhThetaVsNofCalsInRingTofId;Theta [deg];Nof cals per ring;Entries/event", 15, 0, 90, 30, -0.5, 29.5);

    // nof Cals in ring vs XY
    fHM->Create3<TH3D>("fhNofCalsInRingVsXY", "fhNofCalsInRingVsXY;PMT X;PMT Y;Nof cals per ring", 25,-0.5, 24.5, 25, -0.5, 24.5, 30, -0.5, 29.5);
    fHM->Create3<TH3D>("fhNofCalsInRingVsXYCand", "fhNofCalsInRingVsXYCand;PMT X;PMT Y;Nof cals per ring", 25,-0.5, 24.5, 25, -0.5, 24.5, 30, -0.5, 29.5);
    fHM->Create3<TH3D>("fhNofCalsInRingVsXYTofId", "fhNofCalsInRingVsXYTofId;PMT X;PMT Y;Nof cals per ring", 25,-0.5, 24.5, 25, -0.5, 24.5, 30, -0.5, 29.5);

    // nof Cals in ring vs Theta Phi
    fHM->Create3<TH3D>("fhNofCalsInRingVsThetaPhi", "fhNofCalsInRingVsThetaPhi;Theta [deg];Phi [deg];Nof cals per ring",  30, 0, 90, 30, 0, 360, 30, -0.5, 29.5);
    fHM->Create3<TH3D>("fhNofCalsInRingVsThetaPhiCand", "fhNofCalsInRingVsThetaPhiCand;Theta [deg];Phi [deg];Nof cals per ring",  30, 0, 90, 30, 0, 360, 30, -0.5, 29.5);
    fHM->Create3<TH3D>("fhNofCalsInRingVsThetaPhiTofId", "fhNofCalsInRingVsThetaPhiTofId;Theta [deg];Phi [deg];Nof cals per ring",  30, 0, 90, 30, 0, 360, 30, -0.5, 29.5);

    fHM->Create1<TH1D>("fhNofCals", "fhNofCals;Nof cals per event;Entries/event", 50, 0, 1500);
    fHM->Create1<TH1D>("fhNofRings", "fhNofRings;Nof rings per event;Entries/event", 20, -0.5, 19.5);
    fHM->Create1<TH1D>("fhNofCands", "fhNofCands;Nof particle cands per event;Entries/event", 50, -0.5, 149.5);
    fHM->Create2<TH2D>("fhNofCalsVsNofRings", "fhNofCalsVsNofRings;Nof cals per event;Nof rings per event;Entries/event", 30, 0, 600, 20, -0.5, 19.5);
    fHM->Create2<TH2D>("fhNofCandsVsNofRings", "fhNofCandsVsNofRings;Nof cands (used) per event;Nof rings per event;Entries/event", 20, -0.5, 39.5, 20, -0.5, 19.5);

    fHM->Create1<TH1D>("fhNofCandsUsedFlag", "fhNofCandsUsedFlag;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);
    fHM->Create1<TH1D>("fhNofCandsTof", "fhNofCandsTofId;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);
    fHM->Create1<TH1D>("fhNofCandsRich", "fhNofCandsRich;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);
    fHM->Create1<TH1D>("fhNofCandsRichTof", "fhNofCandsRichTof;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);
    fHM->Create1<TH1D>("fhNofCandsTofId", "fhNofCandsTofId;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);
    fHM->Create1<TH1D>("fhNofCandsRichTofId", "fhNofCandsRichTofId;Nof particle cands per event;Entries/event", 30, -0.5, 29.5);

    fHM->Create1<TH1D>("fhMomUsedFlag", "fhMomUsedFlag;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomTof", "fhMomTof;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomRich", "fhMomRich;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomRichTof", "fhMomRichTof;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomTofId", "fhMomTofId;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomRichTofId", "fhMomRichTofId;Momentum [MeV/c];Entries/event", 40, 0, 1000);

    fHM->Create2<TH2D>("fhThetaVsPhiUsedFlag", "fhThetaVsPhiUsedFlag;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiTof", "fhThetaVsPhiTof;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiRich", "fhThetaVsPhiRich;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiRichTof", "fhThetaVsPhiRichTof;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiTofId", "fhThetaVsPhiTofId;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiRichTofId", "fhThetaVsPhiRichTofId;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);

    //Plus and Minus efficiency
    fHM->Create1<TH1D>("fhMomTofIdMinus", "fhMomTofIdMinus;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomRichTofIdMinus", "fhMomRichTofIdMinus;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomTofIdPlus", "fhMomTofIdPlus;Momentum [MeV/c];Entries/event", 40, 0, 1000);
    fHM->Create1<TH1D>("fhMomRichTofIdPlus", "fhMomRichTofIdPlus;Momentum [MeV/c];Entries/event", 40, 0, 1000);

    fHM->Create2<TH2D>("fhThetaVsPhiTofIdMinus", "fhThetaVsPhiTofIdMinus;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiRichTofIdMinus", "fhThetaVsPhiRichTofIdMinus;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiTofIdPlus", "fhThetaVsPhiTofIdPlus;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);
    fHM->Create2<TH2D>("fhThetaVsPhiRichTofIdPlus", "fhThetaVsPhiRichTofIdPlus;Theta [deg];Phi [deg];Entries/event", 30, 0, 90, 30, 0, 360);


    fHM->Create1<TH1D>("fhRichMatchQaCand", "fhRichMatchQaCand;RICH Matching Qa;Entries/event", 40, -2, 12.);
    fHM->Create2<TH2D>("fhDThetaVsDPhiCand", "fhDThetaVsDPhiCand;dTheta [deg];dPhi [deg];Entries/event", 60, -12, 12, 60, -12, 12);
    fHM->Create1<TH1D>("fhRichMatchQaTofId", "fhRichMatchQaTofId;RICH Matching Qa;Entries/event", 40, -2, 12.);
    fHM->Create2<TH2D>("fhDThetaVsDPhiTofId", "fhDThetaVsDPhiTofId;dTheta [deg];dPhi [deg];Entries/event", 60, -12, 12, 60, -12, 12);


    for (Int_t iS = 0; iS < 6; iS++) {
        stringstream ss;
        ss << "fhCalsThetaForSector" << iS;
        fHM->Create1<TH1D>(ss.str().c_str() , (ss.str() + ";Theta [deg];# Cals/event").c_str(), 100, 0., 100.);
    }

    fHM->Create2<TH2D>("fhRawPmtPixel", "fhRawPmtPixel;PMT ID;Pixel ID;Entries", 600,-0.5, 599.5, 65, -0.5, 64.5);
    fHM->Create1<TH1D>("fhRawMult", "fhRawMult;Raw hit multiplicity;Entries", 10,-0.5, 9.5);
    fHM->Create1<TH1D>("fhRawToT", "fhRawToT;ToT [ns];Entries", 100, 0., 10);
}

Bool_t HRich700GeoTestRealQa::reinit()
{
    return kTRUE;
}


Int_t HRich700GeoTestRealQa::execute()
{
    if (!isCurrentEventGood()) return 0;
    fEventNum++;
    fHM->H1("fhNofEvents")->Fill(0.5);
    //if (fEventNum%1000 == 0) cout << "HRich700GeoTestRealQa::execute nofGoodEvents " << fEventNum << endl;

    fillHistCals();
    fillHistRaw();
    fillHistRings();
    fillHistParticleCands();
    return 0;
}

Bool_t HRich700GeoTestRealQa::isCurrentEventGood()
{
    HParticleEvtInfo* evtInfo = static_cast<HParticleEvtInfo*>(fCatEvtInfo->getObject(0));
    if(evtInfo&&!evtInfo->isGoodEvent(
        Particle::kGoodTRIGGER
    	// | Particle::kGoodVertexClust
        // | Particle::kGoodVertexCand
        // | Particle::kGoodSTART
        // | Particle::kNoPileUpSTART
        // | Particle::kNoVETO
        // | Particle::kGoodSTARTVETO
        // | Particle::kGoodSTARTMETA
    )) return false;

    // HEventHeader* eventHeader = gHades->getCurrentEvent()->getEventHeader();
    // eventHeader->isTBit(13); //PT3
    return true;
}

void HRich700GeoTestRealQa::fillHistCals()
{
    Int_t blobsCut = 28;
    Int_t nofRichCals = fCatRichCal->getEntries();
    Int_t nofRichRings = fCatRichHit->getEntries();
    fHM->H1("fhNofCals")->Fill(nofRichCals);
    fHM->H1("fhNofRings")->Fill(nofRichRings);
    fHM->H1("fhNofCalsVsNofRings")->Fill(nofRichCals, nofRichRings);
    map<Int_t, Int_t> nofCalsPerPmtMap;
    nofCalsPerPmtMap.clear();
    for (Int_t iPmt = 0; iPmt < 600; iPmt++) {
        nofCalsPerPmtMap[iPmt] = 0;
    }
    for (Int_t iC = 0; iC < nofRichCals; iC++) {
       HRichCal* cal = static_cast<HRichCal*>(fCatRichCal->getObject(iC));
       if (NULL == cal) continue;
       Int_t col = cal->getCol();
       Int_t row = cal->getRow();
       Int_t sector = cal->getSector();
       fHM->H2("fhCalsColRow")->Fill(col, row);
       fHM->H1("fhCalsSectors")->Fill(sector);

       Int_t loc[3];
       loc[0] = cal->getSector();
       loc[1] = col;
       loc[2] = row;
       pair<Double_t, Double_t> xy = fDigiPar->getXY(loc);
       fHM->H2("fhCalsXY")->Fill(xy.first, xy.second);

       Int_t pmtId = fDigiPar->getPMTId(col, row);
       fHM->H1("fhCalsPmtId")->Fill(pmtId);
       nofCalsPerPmtMap[pmtId]++;

       HRich700PmtData* pmtData = fDigiPar->getPMTData(pmtId);
       if (pmtData != NULL) {
           fHM->H2("fhCalsPmtXPmtY")->Fill(pmtData->fIndX, pmtData->fIndY);
       }

       Float_t phiDeg, thetaDeg;
       Int_t sectorInterp = fDigiPar->getSectorPhiThetaDegPixels(col, row, phiDeg, thetaDeg);
       fHM->H1("fhCalsPhi")->Fill(phiDeg);
       fHM->H1("fhCalsTheta")->Fill(thetaDeg);
       fHM->H1("fhCalsSectorsInterp")->Fill(sectorInterp);

       stringstream ss;
       ss << "fhCalsThetaForSector" << sectorInterp;
       fHM->H1(ss.str())->Fill(thetaDeg);
    }

    Int_t blobsCounter = 0;
    map<Int_t, Int_t>::iterator it;
    for ( it = nofCalsPerPmtMap.begin(); it != nofCalsPerPmtMap.end(); it++ ) {
        HRich700PmtData* pmtData = fDigiPar->getPMTData(it->first);
        if (pmtData != NULL) {
            fHM->H1("fhNofCalsPerPmt")->Fill(it->second);
            fHM->H3("fhNofCalsPerPmtPmtXPmtY")->Fill(pmtData->fIndX, pmtData->fIndY, it->second);
            // select only blob PMT
            if (it->second >= blobsCut){
                blobsCounter++;
                fHM->H2("fhBlobsPmtXPmtY")->Fill(pmtData->fIndX, pmtData->fIndY);

                // find the closest track to blob PMT
                Int_t candId = getClosestCand(pmtData->fTheta, pmtData->fPhi );
                HParticleCand* cand = (candId >=0)?static_cast<HParticleCand*>(fCatParticleCand->getObject(candId)):NULL;
                if (NULL != cand) {
                    Double_t candTheta = cand->getTheta();
                    Double_t candPhi = cand->getPhi();
                    Double_t dPhi = ( pmtData->fPhi - candPhi ) * TMath::Sin(TMath::DegToRad() * candTheta);
                    Float_t dTheta = pmtData->fTheta - candTheta;
                    Double_t dQa = sqrt(dPhi*dPhi + dTheta*dTheta);
                    fHM->H1("fhBlobsMatchQa")->Fill(dQa);
                    fHM->H2("fhBlobsDThetaVsDPhi")->Fill(dTheta, dPhi);
                }
            }
        }
    }
    Int_t nofCandsUsedFlag = 0;
    for (Int_t iP = 0; iP < fCatParticleCand->getEntries(); iP++) {
        HParticleCand* cand = static_cast<HParticleCand*>(fCatParticleCand->getObject(iP));
        if (NULL != cand && cand->isFlagBit(kIsUsed)) nofCandsUsedFlag++;
    }
    fHM->H1("fhNofBlobsPerEvent")->Fill(blobsCounter);
    fHM->H1("fhNofCandsVsNofBlobs")->Fill(nofCandsUsedFlag, blobsCounter);
}

Int_t HRich700GeoTestRealQa::getClosestCand(Double_t theta, Double_t phi)
{
    Double_t closestQa = 999.;
    Double_t closestCandId = -1;
    Int_t nofPartCands = fCatParticleCand->getEntries();
    for (Int_t iP = 0; iP < nofPartCands; iP++) {
        HParticleCand* cand = static_cast<HParticleCand*>(fCatParticleCand->getObject(iP));
        if (NULL == cand) continue;
        if (!cand->isFlagBit(kIsUsed)) continue;
        Double_t candTheta = cand->getTheta();
        Double_t candPhi = cand->getPhi();
        Double_t dPhi = ( phi - candPhi ) * TMath::Sin(TMath::DegToRad() * candTheta);
        Float_t dTheta = theta - candTheta;
        Double_t dQa = sqrt(dPhi*dPhi + dTheta*dTheta);

        if (closestQa >= dQa) {
            closestQa = dQa;
            closestCandId = iP;
        }
    }
    return closestCandId;
}

void HRich700GeoTestRealQa::fillHistRings()
{
    Int_t nofRichRings = fCatRichHit->getEntries();
    fHM->H1("fhNofRings")->Fill(nofRichRings);
    for (Int_t iR = 0; iR < nofRichRings; iR++) {
       HRichHit* ring = static_cast<HRichHit*>(fCatRichHit->getObject(iR));
       if (NULL == ring) continue;

       Int_t nofCalsInRing = ring->fRich700NofRichCals;
       Double_t radius = ring->fRich700CircleRadius;

       Double_t xc = ring->fRich700CircleCenterX;
       Double_t yc = ring->fRich700CircleCenterY;

       HRich700PmtData* pmtData = fDigiPar->getPMTData( fDigiPar->getPMTId((Float_t)xc, (Float_t)yc) );

       if (radius >= 18){
           fHM->H2("fhRingXYMore18mm")->Fill(xc, yc);
       } else {
           fHM->H2("fhRingXYLess18mm")->Fill(xc, yc);
       }
       fHM->H2("fhRingXY")->Fill(xc, yc);

       fHM->H2("fhThetaVsNofCalsInRing")->Fill(ring->getTheta(), nofCalsInRing);
       if (pmtData != NULL) fHM->H3("fhNofCalsInRingVsXY")->Fill(pmtData->fIndX, pmtData->fIndY, nofCalsInRing);
       fHM->H3("fhNofCalsInRingVsThetaPhi")->Fill(ring->getTheta(), ring->getPhi(), nofCalsInRing);
       fHM->H2("fhThetaVsRadius")->Fill(ring->getTheta(), radius);
       fHM->H2("fhNofCalsInRingVsRadius")->Fill(nofCalsInRing, radius);
    }
}

void HRich700GeoTestRealQa::fillHistParticleCands()
{
    Int_t nofPartCands = fCatParticleCand->getEntries();
    Int_t nofRichRings = fCatRichHit->getEntries();
    fHM->H1("fhNofCands")->Fill(nofPartCands);

    Int_t nofCandsUsedFlag = 0;
    Int_t nofCandsTof = 0;
    Int_t nofCandsRich = 0;
    Int_t nofCandsRichTof = 0;
    Int_t nofCandsTofId = 0;
    Int_t nofCandsRichTofId = 0;

    for (Int_t iP = 0; iP < nofPartCands; iP++) {
        HParticleCand* cand = static_cast<HParticleCand*>(fCatParticleCand->getObject(iP));
        if (NULL == cand) continue;
        if (!cand->isFlagBit(kIsUsed)) continue;

        nofCandsUsedFlag++;
        Double_t candTheta = cand->getTheta();
        Double_t candPhi = cand->getPhi();
        Double_t candBeta = cand->getBeta();
        Double_t candMom = cand->getMomentum() * cand->getCharge();
        Double_t candMomAbs = std::abs(candMom);
        Double_t candDeltaPhi = cand->getDeltaPhi();;
        Double_t candDeltaTheta = cand->getDeltaTheta();
        Double_t candRichMatchQa = cand->getRichMatchingQuality();;

        Bool_t isRich = cand->getRichInd() >= 0;
        Bool_t isTof = cand->getSystemUsed() >= 0;
        Bool_t isTofIdOk =  (candBeta >= 0.95 && candMomAbs < 250.);
        //Bool_t isTofIdOk =  (candBeta >= 0.9 && candMom < 0. && candMom > -250);

        fHM->H1("fhDThetaVsDPhiCand")->Fill(candDeltaTheta, candDeltaPhi);
        fHM->H1("fhRichMatchQaCand")->Fill(candRichMatchQa);
        fHM->H1("fhMomVsBetaCand")->Fill(candMom, candBeta);
        if (candRichMatchQa > 0){
            if (candRichMatchQa < 1) fHM->H1("fhMomVsBetaRichMatch<1")->Fill(candMom, candBeta);
            if (candRichMatchQa < 2) fHM->H1("fhMomVsBetaRichMatch<2")->Fill(candMom, candBeta);
            if (candRichMatchQa < 3) fHM->H1("fhMomVsBetaRichMatch<3")->Fill(candMom, candBeta);
            if (candRichMatchQa < 4) fHM->H1("fhMomVsBetaRichMatch<4")->Fill(candMom, candBeta);
        }

        fHM->H1("fhMomUsedFlag")->Fill(candMomAbs);
        fHM->H2("fhThetaVsPhiUsedFlag")->Fill(candTheta, candPhi);
        if (isTof) {
            nofCandsTof++;
            fHM->H1("fhMomTof")->Fill(candMomAbs);
            fHM->H2("fhThetaVsPhiTof")->Fill(candTheta, candPhi);
        }
        if (isTofIdOk) {
            nofCandsTofId++;
            fHM->H1("fhMomTofId")->Fill(candMomAbs);
            fHM->H2("fhThetaVsPhiTofId")->Fill(candTheta, candPhi);
            if (candMom > 0) {
                fHM->H1("fhMomTofIdPlus")->Fill(candMomAbs);
                fHM->H2("fhThetaVsPhiTofIdPlus")->Fill(candTheta, candPhi);
            } else {
                fHM->H1("fhMomTofIdMinus")->Fill(candMomAbs);
                fHM->H2("fhThetaVsPhiTofIdMinus")->Fill(candTheta, candPhi);
            }
            fHM->H1("fhDThetaVsDPhiTofId")->Fill(candDeltaTheta, candDeltaPhi);
            fHM->H1("fhRichMatchQaTofId")->Fill(candRichMatchQa);
        }

        if (isRich) {
            nofCandsRich++;
            fHM->H1("fhMomRich")->Fill(candMomAbs);
            fHM->H2("fhThetaVsPhiRich")->Fill(candTheta, candPhi);
            HRichHit* ring = static_cast<HRichHit*>(fCatRichHit->getObject(cand->getRichInd()));
            if (NULL != ring) {
                fHM->H2("fhThetaVsRadiusCand")->Fill(ring->getTheta(), ring->fRich700CircleRadius);
                fHM->H2("fhNofCalsInRingVsRadiusCand")->Fill(ring->fRich700NofRichCals, ring->fRich700CircleRadius);
                fHM->H2("fhThetaVsNofCalsInRingCand")->Fill(ring->getTheta(), ring->fRich700NofRichCals);
                HRich700PmtData* pmtData = fDigiPar->getPMTData( fDigiPar->getPMTId((Float_t)ring->fRich700CircleCenterX, (Float_t)ring->fRich700CircleCenterY) );
                if (pmtData!=NULL) fHM->H3("fhNofCalsInRingVsXYCand")->Fill(pmtData->fIndX, pmtData->fIndY, ring->fRich700NofRichCals);
                fHM->H3("fhNofCalsInRingVsThetaPhiCand")->Fill(ring->getTheta(), ring->getPhi(), ring->fRich700NofRichCals);
                if (ring->fRich700CircleRadius >= 18){
                    fHM->H1("fhMomVsBetaMore18mm")->Fill(candMom, candBeta);
                } else {
                    fHM->H1("fhMomVsBetaLess18mm")->Fill(candMom, candBeta);
                }
                fHM->H2("fhRingXYCand")->Fill(ring->fRich700CircleCenterX, ring->fRich700CircleCenterY);

                if (isTofIdOk) {
                    fHM->H2("fhThetaVsRadiusTofId")->Fill(ring->getTheta(), ring->fRich700CircleRadius);
                    fHM->H2("fhNofCalsInRingVsRadiusTofId")->Fill(ring->fRich700NofRichCals, ring->fRich700CircleRadius);
                    fHM->H2("fhThetaVsNofCalsInRingTofId")->Fill(ring->getTheta(), ring->fRich700NofRichCals);
                    HRich700PmtData* pmtData = fDigiPar->getPMTData( fDigiPar->getPMTId((Float_t)ring->fRich700CircleCenterX, (Float_t)ring->fRich700CircleCenterY) );
                    if (pmtData!=NULL) fHM->H3("fhNofCalsInRingVsXYTofId")->Fill(pmtData->fIndX, pmtData->fIndY, ring->fRich700NofRichCals);
                    fHM->H3("fhNofCalsInRingVsThetaPhiTofId")->Fill(ring->getTheta(), ring->getPhi(), ring->fRich700NofRichCals);
                    fHM->H2("fhRingXYTofId")->Fill(ring->fRich700CircleCenterX, ring->fRich700CircleCenterY);
                    fHM->H2("fhRichMatchQaVsRadiusTofId")->Fill(candRichMatchQa, ring->fRich700CircleRadius);
                }
                if (cand->getSystemUsed() == 0){
                    fHM->H2("fhBetaVsRadiusSys0")->Fill(candBeta, ring->fRich700CircleRadius);
                } else if (cand->getSystemUsed() == 1){
                    fHM->H2("fhBetaVsRadiusSys1")->Fill(candBeta, ring->fRich700CircleRadius);
                }
                fHM->H2("fhMomVsRadius")->Fill(candMom, ring->fRich700CircleRadius);
                fHM->H2("fhRichMatchQaVsRadiusCand")->Fill(candRichMatchQa, ring->fRich700CircleRadius);
            }
        }

        if (isRich && isTof) {
            fHM->H1("fhMomRichTof")->Fill(candMomAbs);
            fHM->H2("fhThetaVsPhiRichTof")->Fill(candTheta, candPhi);

            nofCandsRichTof++;
        }
        if (isRich && isTofIdOk) {
            fHM->H1("fhMomRichTofId")->Fill(candMomAbs);
            fHM->H2("fhThetaVsPhiRichTofId")->Fill(candTheta, candPhi);
            if (candMom > 0) {
                fHM->H1("fhMomRichTofIdPlus")->Fill(candMomAbs);
                fHM->H2("fhThetaVsPhiRichTofIdPlus")->Fill(candTheta, candPhi);
            } else {
                fHM->H1("fhMomRichTofIdMinus")->Fill(candMomAbs);
                fHM->H2("fhThetaVsPhiRichTofIdMinus")->Fill(candTheta, candPhi);
            }
            nofCandsRichTofId++;
        }
    }

    fHM->H1("fhNofCandsUsedFlag")->Fill(nofCandsUsedFlag);
    fHM->H1("fhNofCandsRich")->Fill(nofCandsRich);
    fHM->H1("fhNofCandsTof")->Fill(nofCandsTof);
    fHM->H1("fhNofCandsRichTof")->Fill(nofCandsRichTof);
    fHM->H1("fhNofCandsTofId")->Fill(nofCandsTofId);
    fHM->H1("fhNofCandsRichTofId")->Fill(nofCandsRichTofId);

    fHM->H2("fhNofCandsVsNofRings")->Fill(nofCandsUsedFlag, nofRichRings);

    for (Int_t iP = 0; iP < nofPartCands; iP++) {
        HParticleCand* cand = static_cast<HParticleCand*>(fCatParticleCand->getObject(iP));
        if (NULL == cand) continue;
        if (!cand->isFlagBit(kIsUsed)) continue;
        if (cand->getRichInd() >= 0) {
            HRichHit* ring = static_cast<HRichHit*>(fCatRichHit->getObject(cand->getRichInd()));
            if (NULL != ring) {
                fHM->H1("fhNofUsedFlagCandVsRadius")->Fill(nofCandsUsedFlag, ring->fRich700CircleRadius);
            }
        }
    }
}

void HRich700GeoTestRealQa::fillHistRaw()
{
    if (fCatRichRaw == NULL) return;
    Int_t nofRichRaws = fCatRichRaw->getEntries();
    for (Int_t iR = 0; iR < nofRichRaws; iR++){
       HRich700Raw* raw = static_cast<HRich700Raw*>(fCatRichRaw->getObject(iR));
       if (NULL == raw) continue;
       Int_t pmtId = raw->getPMT();
       Int_t pixel = raw->getPixel();
       Int_t mult = raw->getMultiplicity();
       fHM->H2("fhRawPmtPixel")->Fill(pmtId, pixel);
       fHM->H1("fhRawMult")->Fill(mult);
       for (Int_t iH = 0; iH < mult; iH++) {
           const HRich700hit_t* hit = raw->getHit(iH);
           if (hit == nullptr) continue;
           fHM->H1("fhRawToT")->Fill(hit->fToT);
       }
    }
}


void HRich700GeoTestRealQa::drawHist(Bool_t doEventScale)
{
    HRichDrawHist::SetDefaultDrawStyle();
    Double_t nofEvents = fHM->H1("fhNofEvents")->GetEntries();
    Double_t eventScale = (doEventScale)?1./nofEvents:1.;
    cout << "Number of events:" << nofEvents << endl;
    {
	TCanvas* c = fHM->CreateCanvas("hrich_fhCalsXY", "hrich_fhCalsXY", 1800, 600);
    c->Divide(3, 1);
    c->cd(1);
    fHM->Scale("fhCalsColRow", eventScale);
	HRichDrawHist::DrawH2(fHM->H2("fhCalsColRow"));
    c->cd(2);
    fHM->Scale("fhCalsXY", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhCalsXY"));
    c->cd(3);
    fHM->Scale("fhCalsPmtXPmtY", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhCalsPmtXPmtY"));
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhCals", "hrich_fhCals", 1200, 1200);
    c->Divide(2,2);
    c->cd(1);
    fHM->Scale("fhCalsPmtId", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhCalsPmtId"));
    c->cd(2);
    fHM->Scale("fhCalsSectors", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhCalsSectors"));
    c->cd(3);
    fHM->Scale("fhNofCalsPerPmt", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhNofCalsPerPmt"), kLinear, kLog);
    c->cd(4);
    TH2D* h = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsPerPmtPmtXPmtY"));
    HRichDrawHist::DrawH2(h);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhBlobs", "hrich_fhBlobs", 1800, 600);
    c->Divide(3,1);
    c->cd(1);
    fHM->NormalizeToIntegral("fhNofBlobsPerEvent");
    fHM->Scale("fhNofBlobsPerEvent", 100.);
    HRichDrawHist::DrawH1(fHM->H1("fhNofBlobsPerEvent"));
    c->cd(2);
    fHM->Scale("fhBlobsPmtXPmtY", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhBlobsPmtXPmtY"));
    c->cd(3);
    fHM->Scale("fhNofCandsVsNofBlobs", eventScale);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhNofCandsVsNofBlobs"));
    fHM->H2("fhNofCandsVsNofBlobs")->GetYaxis()->SetRangeUser(-.5, 3.5);
    gPad->SetLogz();
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhBlobsClosestCand", "hrich_fhBlobsClosestCand", 1200, 1200);
    c->Divide(2, 2);
    c->cd(1);
    fHM->Scale("fhBlobsDThetaVsDPhi", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhBlobsDThetaVsDPhi"));
    c->cd(2);
    fHM->Scale("fhBlobsMatchQa", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhBlobsMatchQa"), kLinear, kLog);
    c->cd(3);
    TH1D* prY = fHM->H2("fhBlobsDThetaVsDPhi")->ProjectionY();
    prY->GetYaxis()->SetTitle("Entries/event");
    HRichDrawHist::DrawH1(prY);
    c->cd(4);
    TH1D* prX = fHM->H2("fhBlobsDThetaVsDPhi")->ProjectionX();
    prX->GetYaxis()->SetTitle("Entries/event");
    HRichDrawHist::DrawH1(prX);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhNofPerEvent", "hrich_fhNofPerEvent", 1800, 600);
    c->Divide(3, 1);
    c->cd(1);
    fHM->Scale("fhNofCals", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhNofCals"));
    string txtCals = RichUtils::NumberToString(fHM->H1("fhNofCals")->GetMean()) + "/" + RichUtils::NumberToString(fHM->H1("fhNofCals")->GetRMS());
    HRichDrawHist::DrawTextOnPad(txtCals, 0.1, 0.90, 0.9, 0.99);
    c->cd(2);
    fHM->Scale("fhNofRings", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhNofRings"));
    string txtRings = RichUtils::NumberToString(fHM->H1("fhNofRings")->GetMean()) + "/" + RichUtils::NumberToString(fHM->H1("fhNofRings")->GetRMS());
    HRichDrawHist::DrawTextOnPad(txtRings, 0.1, 0.90, 0.9, 0.99);
    c->cd(3);
    fHM->Scale("fhNofCands", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhNofCands"));
    string txtCands = RichUtils::NumberToString(fHM->H1("fhNofCands")->GetMean()) + "/" + RichUtils::NumberToString(fHM->H1("fhNofCands")->GetRMS());
    HRichDrawHist::DrawTextOnPad(txtCands, 0.1, 0.90, 0.9, 0.99);
    }

    {
    fHM->CreateCanvas("hrich_fhNofCands", "hrich_fhNofCands", 1000, 800);
    vector<TH1*> vH;
    vH.push_back(fHM->H1("fhNofCandsUsedFlag"));
    vH.push_back(fHM->H1("fhNofCandsRich"));
    vH.push_back(fHM->H1("fhNofCandsTof"));
    vH.push_back(fHM->H1("fhNofCandsRichTof"));
    vH.push_back(fHM->H1("fhNofCandsTofId"));
    vH.push_back(fHM->H1("fhNofCandsRichTofId"));
    for (UInt_t i = 0; i < vH.size();i++){vH[i]->Scale(eventScale);}
    fHM->H1("fhNofCandsUsedFlag")->SetMinimum(0.95 * eventScale);
    vector<string> vLabel;
    vLabel.push_back("UsedFlag (" + RichUtils::NumberToString(fHM->H1("fhNofCandsUsedFlag")->GetMean()) + ")");
    vLabel.push_back("RICH (" + RichUtils::NumberToString(fHM->H1("fhNofCandsRich")->GetMean()) + ")");
    vLabel.push_back("TOF (" + RichUtils::NumberToString(fHM->H1("fhNofCandsTof")->GetMean()) + ")");
    vLabel.push_back("RICH+TOF (" + RichUtils::NumberToString(fHM->H1("fhNofCandsRichTof")->GetMean()) + ")");
    vLabel.push_back("TOF ID (" + RichUtils::NumberToString(fHM->H1("fhNofCandsTofId")->GetMean()) + ")");
    vLabel.push_back("RICH+TOF ID(" + RichUtils::NumberToString(fHM->H1("fhNofCandsRichTofId")->GetMean()) + ")");
    HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog);
    }

    {
        // calculate efficiency before scaling histograms
        TCanvas* c = fHM->CreateCanvas("hrich_fhRichEff", "hrich_fhRichEff", 1000, 1000);
        c->Divide(2,2);
        TH1D* pxEff1 = RichUtils::DivideH1((TH1D*)fHM->H1("fhMomRichTofId")->Clone(), (TH1D*)fHM->H1("fhMomTofId")->Clone(), "", 100., "Relative efficiency [%]");
        pxEff1->SetMinimum(0.);

        TH1D* pxEff2 = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofId")->ProjectionY()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofId")->ProjectionY()->Clone(), "", 100., "Relative efficiency [%]");
        pxEff2->SetMinimum(0.);

        TH1D* pxEff3 = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofId")->ProjectionX()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofId")->ProjectionX()->Clone(), "", 100., "Relative efficiency [%]");
        pxEff3->SetMinimum(0.);

        TH2D* pxEff4 = RichUtils::DivideH2((TH2D*)fHM->H2("fhThetaVsPhiRichTofId")->Clone(), (TH2D*)fHM->H2("fhThetaVsPhiTofId")->Clone(), "", 100., "Relative efficiency [%]");
        pxEff4->SetMinimum(0.);

        c->cd(1);
        vector<TH1*> hist1;
        hist1.push_back(pxEff1);
        vector<string> labels1;
        stringstream ss1;
        ss1.precision(1);
        ss1 << "Efficiency (" << fixed << 100.*fHM->H1("fhMomRichTofId")->GetEntries() / fHM->H1("fhMomTofId")->GetEntries() << "%)";
        labels1.push_back(ss1.str());
        HRichDrawHist::DrawH1(hist1, labels1);
        hist1[0]->GetXaxis()->SetRangeUser(0, 500);

        c->cd(2);
        vector<TH1*> hist2;
        hist2.push_back(pxEff2);
        vector<string> labels2;
        stringstream ss2;
        ss2.precision(1);
        ss2 << "Efficiency (" << fixed << 100.*fHM->H1("fhThetaVsPhiRichTofId")->GetEntries() / fHM->H1("fhThetaVsPhiTofId")->GetEntries() << "%)";
        labels2.push_back(ss2.str());
        HRichDrawHist::DrawH1(hist2, labels2);

        c->cd(3);
        vector<TH1*> hist3;
        hist3.push_back(pxEff3);
        vector<string> labels3;
        stringstream ss3;
        ss3.precision(1);
        ss3 << "Efficiency (" << fixed << 100.*fHM->H2("fhThetaVsPhiRichTofId")->GetEntries() / fHM->H2("fhThetaVsPhiTofId")->GetEntries() << "%)";
        labels3.push_back(ss3.str());
        HRichDrawHist::DrawH1(hist3, labels3);

        c->cd(4);
        HRichDrawHist::DrawH2(pxEff4);
    }

    {
        TCanvas* c = fHM->CreateCanvas("hrich_fhRichEffPlusMinus", "hrich_fhRichEffPlusMinus", 1800, 600);
        c->Divide(3,1);
        TH1D* pxEff1Plus = RichUtils::DivideH1((TH1D*)fHM->H1("fhMomRichTofIdPlus")->Clone(), (TH1D*)fHM->H1("fhMomTofIdPlus")->Clone(), "", 100., "Relative efficiency [%]");
        TH1D* pxEff1Minus = RichUtils::DivideH1((TH1D*)fHM->H1("fhMomRichTofIdMinus")->Clone(), (TH1D*)fHM->H1("fhMomTofIdMinus")->Clone(), "", 100., "Relative efficiency [%]");
        pxEff1Plus->SetMinimum(0.);

        TH1D* pxEff2Plus = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofIdPlus")->ProjectionY()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofIdPlus")->ProjectionY()->Clone(), "", 100., "Relative efficiency [%]");
        TH1D* pxEff2Minus = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofIdMinus")->ProjectionY()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofIdMinus")->ProjectionY()->Clone(), "", 100., "Relative efficiency [%]");
        pxEff2Plus->SetMinimum(0.);

        TH1D* pxEff3Plus = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofIdPlus")->ProjectionX()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofIdPlus")->ProjectionX()->Clone(), "", 100., "Relative efficiency [%]");
        TH1D* pxEff3Minus = RichUtils::DivideH1((TH1D*)fHM->H2("fhThetaVsPhiRichTofIdMinus")->ProjectionX()->Clone(), (TH1D*)fHM->H2("fhThetaVsPhiTofIdMinus")->ProjectionX()->Clone(), "", 100., "Relative efficiency [%]");
        pxEff3Plus->SetMinimum(0.);

        c->cd(1);
        vector<TH1*> hist1;
        hist1.push_back(pxEff1Plus);
        hist1.push_back(pxEff1Minus);
        vector<string> labels1;
        stringstream ss1Plus;
        ss1Plus.precision(1);
        ss1Plus << "Plus (" << fixed << 100.*fHM->H1("fhMomRichTofIdPlus")->GetEntries() / fHM->H1("fhMomTofIdPlus")->GetEntries() << "%)";
        stringstream ss1Minus;
        ss1Minus.precision(1);
        ss1Minus << "Minus (" << fixed << 100.*fHM->H1("fhMomRichTofIdMinus")->GetEntries() / fHM->H1("fhMomTofIdMinus")->GetEntries() << "%)";
        labels1.push_back(ss1Plus.str());
        labels1.push_back(ss1Minus.str());
        HRichDrawHist::DrawH1(hist1, labels1);
        hist1[0]->GetXaxis()->SetRangeUser(0, 500);

        c->cd(2);
        vector<TH1*> hist2;
        hist2.push_back(pxEff2Plus);
        hist2.push_back(pxEff2Minus);
        vector<string> labels2;
        stringstream ss2Plus;
        ss2Plus.precision(1);
        ss2Plus << "Plus (" << fixed << 100.*fHM->H1("fhThetaVsPhiRichTofIdPlus")->GetEntries() / fHM->H1("fhThetaVsPhiTofIdPlus")->GetEntries() << "%)";
        stringstream ss2Minus;
        ss2Minus.precision(1);
        ss2Minus << "Minus (" << fixed << 100.*fHM->H1("fhThetaVsPhiRichTofIdMinus")->GetEntries() / fHM->H1("fhThetaVsPhiTofIdMinus")->GetEntries() << "%)";
        labels2.push_back(ss2Plus.str());
        labels2.push_back(ss2Minus.str());
        HRichDrawHist::DrawH1(hist2, labels2);

        c->cd(3);
        vector<TH1*> hist3;
        hist3.push_back(pxEff3Plus);
        hist3.push_back(pxEff3Minus);
        vector<string> labels3;
        stringstream ss3Plus;
        ss3Plus.precision(1);
        ss3Plus << "Plus (" << fixed << 100.*fHM->H2("fhThetaVsPhiRichTofIdPlus")->GetEntries() / fHM->H2("fhThetaVsPhiTofIdPlus")->GetEntries() << "%)";
        stringstream ss3Minus;
        ss3Minus.precision(1);
        ss3Minus << "Minus (" << fixed << 100.*fHM->H2("fhThetaVsPhiRichTofIdMinus")->GetEntries() / fHM->H2("fhThetaVsPhiTofIdMinus")->GetEntries() << "%)";
        labels3.push_back(ss3Plus.str());
        labels3.push_back(ss3Minus.str());
        HRichDrawHist::DrawH1(hist3, labels3);
    }

    {
    fHM->CreateCanvas("hrich_fhMomCand", "hrich_fhMomCand", 1000, 800);
    vector<TH1*> vH;
    vH.push_back(fHM->H1("fhMomUsedFlag"));
    vH.push_back(fHM->H1("fhMomRich"));
    vH.push_back(fHM->H1("fhMomTof"));
    vH.push_back(fHM->H1("fhMomRichTof"));
    vH.push_back(fHM->H1("fhMomTofId"));
    vH.push_back(fHM->H1("fhMomRichTofId"));
    for (UInt_t i = 0; i < vH.size();i++){vH[i]->Scale(eventScale);}
    fHM->H1("fhMomUsedFlag")->SetMinimum(0.95 * eventScale);
    vector<string> vLabel;
    vLabel.push_back("UsedFlag (" + RichUtils::NumberToString(fHM->H1("fhMomUsedFlag")->GetEntries()*eventScale) + ")");
    vLabel.push_back("RICH (" + RichUtils::NumberToString(fHM->H1("fhMomRich")->GetEntries()*eventScale) + ")");
    vLabel.push_back("TOF (" + RichUtils::NumberToString(fHM->H1("fhMomTof")->GetEntries()*eventScale) + ")");
    vLabel.push_back("RICH+TOF (" + RichUtils::NumberToString(fHM->H1("fhMomRichTof")->GetEntries()*eventScale) + ")");
    vLabel.push_back("TOF ID (" + RichUtils::NumberToString(fHM->H1("fhMomTofId")->GetEntries()*eventScale) + ")");
    vLabel.push_back("RICH+TOF ID(" + RichUtils::NumberToString(fHM->H1("fhMomRichTofId")->GetEntries()*eventScale) + ")");
    HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog, true, 0.7, 0.8, 0.99, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhThetaVsPhiCand1D", "hrich_fhThetaVsPhiCand1D", 1200, 600);
    c->Divide(2,1);
    c->cd(1);
    {
        vector<TH1*> vH;
        vH.push_back(fHM->H2("fhThetaVsPhiUsedFlag")->ProjectionY());
        vH.push_back(fHM->H2("fhThetaVsPhiRich")->ProjectionY());
        vH.push_back(fHM->H2("fhThetaVsPhiTof")->ProjectionY());
        vH.push_back(fHM->H2("fhThetaVsPhiRichTof")->ProjectionY());
        vH.push_back(fHM->H2("fhThetaVsPhiTofId")->ProjectionY());
        vH.push_back(fHM->H2("fhThetaVsPhiRichTofId")->ProjectionY());
        for (UInt_t i = 0; i < vH.size();i++){vH[i]->Scale(eventScale);}
        vH[0]->SetMinimum(0.95 * eventScale);
        vector<string> vLabel;
        vLabel.push_back("UsedFlag (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiUsedFlag")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRich")->GetEntries()*eventScale) + ")");
        vLabel.push_back("TOF (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTof")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH+TOF (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTof")->GetEntries()*eventScale) + ")");
        vLabel.push_back("TOF ID (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTofId")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH+TOF ID(" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTofId")->GetEntries()*eventScale) + ")");
        HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog, true, 0.7, 0.8, 0.99, 0.99);
    }
    c->cd(2);
    {
        vector<TH1*> vH;
        vH.push_back(fHM->H2("fhThetaVsPhiUsedFlag")->ProjectionX());
        vH.push_back(fHM->H2("fhThetaVsPhiRich")->ProjectionX());
        vH.push_back(fHM->H2("fhThetaVsPhiTof")->ProjectionX());
        vH.push_back(fHM->H2("fhThetaVsPhiRichTof")->ProjectionX());
        vH.push_back(fHM->H2("fhThetaVsPhiTofId")->ProjectionX());
        vH.push_back(fHM->H2("fhThetaVsPhiRichTofId")->ProjectionX());
        for (UInt_t i = 0; i < vH.size();i++){vH[i]->Scale(eventScale);}
        vH[0]->SetMinimum(0.95 * eventScale);
        vector<string> vLabel;
        vLabel.push_back("UsedFlag (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiUsedFlag")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRich")->GetEntries()*eventScale) + ")");
        vLabel.push_back("TOF (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTof")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH+TOF (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTof")->GetEntries()*eventScale) + ")");
        vLabel.push_back("TOF ID (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTofId")->GetEntries()*eventScale) + ")");
        vLabel.push_back("RICH+TOF ID(" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTofId")->GetEntries()*eventScale) + ")");
        HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog, true, 0.7, 0.8, 0.99, 0.99);
    }
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhThetaVsPhiCand2D", "hrich_fhThetaVsPhiCand2D", 1800, 1200);
    c->Divide(3,2);
    c->cd(1);
    fHM->Scale("fhThetaVsPhiUsedFlag", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiUsedFlag"));
    HRichDrawHist::DrawTextOnPad("UsedFlag (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiUsedFlag")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    c->cd(2);
    fHM->Scale("fhThetaVsPhiTof", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiTof"));
    HRichDrawHist::DrawTextOnPad("Tof (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTof")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    c->cd(3);
    fHM->Scale("fhThetaVsPhiRich", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiRich"));
    HRichDrawHist::DrawTextOnPad("Rich (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRich")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    c->cd(4);
    fHM->Scale("fhThetaVsPhiRichTof", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiRichTof"));
    HRichDrawHist::DrawTextOnPad("RichTof (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTof")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    c->cd(5);
    fHM->Scale("fhThetaVsPhiTofId", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiTofId"));
    HRichDrawHist::DrawTextOnPad("TofId (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiTofId")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    c->cd(6);
    fHM->Scale("fhThetaVsPhiRichTofId", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhThetaVsPhiRichTofId"));
    HRichDrawHist::DrawTextOnPad("RichTofId (" + RichUtils::NumberToString(fHM->H2("fhThetaVsPhiRichTofId")->GetEntries()*eventScale) + ")", 0.2, 0.90, 0.7, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhNofPerEvent2D", "hrich_fhNofPerEvent2D", 1600, 700);
    c->Divide(2, 1);
    c->cd(1);
    fHM->Scale("fhNofCandsVsNofRings", eventScale);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhNofCandsVsNofRings"));
    gPad->SetLogz();
    c->cd(2);
    fHM->Scale("fhNofCalsVsNofRings", eventScale);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhNofCalsVsNofRings"));
    gPad->SetLogz();
    }

    {
	TCanvas* c = fHM->CreateCanvas("hrich_fhMomVsBeta", "hrich_fhMomVsBeta", 1800, 1200);
    c->Divide(3, 2);
    c->cd(1);
    fHM->Scale("fhMomVsBetaCand", eventScale);
	HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaCand"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Candidates (isUsed)", 0.2, 0.90, 0.8, 0.99);
    c->cd(2);
    fHM->Scale("fhMomVsBetaMore18mm", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaMore18mm"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Radius > 18 mm", 0.2, 0.90, 0.8, 0.99);
    c->cd(3);
    fHM->Scale("fhMomVsBetaLess18mm", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaLess18mm"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Radius < 18 mm", 0.2, 0.90, 0.8, 0.99);

    Double_t max = std::max(fHM->H2("fhMomVsBetaMore18mm")->GetMaximum(), fHM->H2("fhMomVsBetaLess18mm")->GetMaximum());
    fHM->H2("fhMomVsBetaMore18mm")->SetMaximum(1.2 * max);
    fHM->H2("fhMomVsBetaLess18mm")->SetMaximum(1.2 * max);

    c->cd(4);
    TH2D* eff1 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaMore18mm")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff1);
    eff1->SetMinimum(0.); eff1->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("Radius > 18 mm", 0.2, 0.90, 0.8, 0.99);
    c->cd(5);
    TH2D* eff2 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaLess18mm")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff2);
    eff2->SetMinimum(0.); eff2->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("Radius < 18 mm", 0.2, 0.90, 0.8, 0.99);

    }

    {
	TCanvas* c = fHM->CreateCanvas("hrich_fhMomVsBeta2", "hrich_fhMomVsBeta2", 1800, 900);
    c->Divide(4, 2);
    c->cd(1);
    fHM->Scale("fhMomVsBetaRichMatch<1", eventScale);
	HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaRichMatch<1"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 1", 0.2, 0.90, 0.8, 0.99);
    c->cd(2);
    fHM->Scale("fhMomVsBetaRichMatch<2", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaRichMatch<2"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 2", 0.2, 0.90, 0.8, 0.99);
    c->cd(3);
    fHM->Scale("fhMomVsBetaRichMatch<3", eventScale);
	HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaRichMatch<3"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 3", 0.2, 0.90, 0.8, 0.99);
    c->cd(4);
    fHM->Scale("fhMomVsBetaRichMatch<4", eventScale);
	HRichDrawHist::DrawH2(fHM->H2("fhMomVsBetaRichMatch<4"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 4", 0.2, 0.90, 0.8, 0.99);
    c->cd(5);
    TH2D* eff1 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaRichMatch<1")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff1);
    eff1->SetMinimum(0.); eff1->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 1", 0.2, 0.90, 0.8, 0.99);
    c->cd(6);
    TH2D* eff2 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaRichMatch<2")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff2);
    eff2->SetMinimum(0.); eff2->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 2", 0.2, 0.90, 0.8, 0.99);
    c->cd(7);
    TH2D* eff3 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaRichMatch<3")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff3);
    eff3->SetMinimum(0.); eff3->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 3", 0.2, 0.90, 0.8, 0.99);
    c->cd(8);
    TH2D* eff4 = RichUtils::DivideH2((TH2D*)fHM->H2("fhMomVsBetaRichMatch<4")->Clone(), (TH2D*)fHM->H2("fhMomVsBetaCand")->Clone(), "", 100., "% to all candidates");
    HRichDrawHist::DrawH2(eff4);
    eff4->SetMinimum(0.); eff4->SetMaximum(100.);
    HRichDrawHist::DrawTextOnPad("RICH matching Qa < 4", 0.2, 0.90, 0.8, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRingXY", "hrich_fhRingXY", 1800, 1200);
    c->Divide(3, 2);
    c->cd(1);
    fHM->Scale("fhRingXYMore18mm", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRingXYMore18mm"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Radius > 18 mm", 0.3, 0.90, 0.7, 0.99);
    c->cd(2);
    fHM->Scale("fhRingXYLess18mm", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRingXYLess18mm"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Radius < 18 mm", 0.3, 0.90, 0.7, 0.99);
    c->cd(3);
    fHM->Scale("fhRingXY", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRingXY"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("All rings", 0.3, 0.90, 0.7, 0.99);
    c->cd(4);
    fHM->Scale("fhRingXYCand", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRingXYCand"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Candidates", 0.3, 0.90, 0.7, 0.99);
    c->cd(5);
    fHM->Scale("fhRingXYTofId", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRingXYTofId"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.3, 0.90, 0.7, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRichMatchQa", "hrich_fhRichMatchQa", 1800, 1200);
    c->Divide(3, 2);
    {
        c->cd(1);
        fHM->Scale("fhDThetaVsDPhiCand", eventScale);
        HRichDrawHist::DrawH2(fHM->H2("fhDThetaVsDPhiCand"));
        HRichDrawHist::DrawTextOnPad("Candidates (isUsed)", 0.2, 0.90, 0.8, 0.99);
        c->cd(2);
        TH1D* prY = fHM->H2("fhDThetaVsDPhiCand")->ProjectionY();
        TH1D* prX = fHM->H2("fhDThetaVsDPhiCand")->ProjectionX();
        prY->GetYaxis()->SetTitle("Entries/event");
        prY->GetXaxis()->SetTitle("dAngle [deg]");
        vector<TH1*> vH;
        vH.push_back(prY);
        vH.push_back(prX);
        vector<string> vLabel;
        vLabel.push_back("dPhi");
        vLabel.push_back("dTheta");
        HRichDrawHist::DrawH1(vH, vLabel);
        c->cd(3);
        fHM->Scale("fhRichMatchQaCand", eventScale);
        HRichDrawHist::DrawH1(fHM->H1("fhRichMatchQaCand"), kLinear, kLog);
    }
    {
        c->cd(4);
        fHM->Scale("fhDThetaVsDPhiTofId", eventScale);
        HRichDrawHist::DrawH2(fHM->H2("fhDThetaVsDPhiTofId"));
        HRichDrawHist::DrawTextOnPad("TOF ID", 0.3, 0.90, 0.7, 0.99);
        c->cd(5);
        TH1D* prY = fHM->H2("fhDThetaVsDPhiTofId")->ProjectionY();
        TH1D* prX = fHM->H2("fhDThetaVsDPhiTofId")->ProjectionX();
        prY->GetYaxis()->SetTitle("Entries/event");
        prY->GetXaxis()->SetTitle("dAngle [deg]");
        vector<TH1*> vH;
        vH.push_back(prY);
        vH.push_back(prX);
        vector<string> vLabel;
        vLabel.push_back("dPhi");
        vLabel.push_back("dTheta");
        HRichDrawHist::DrawH1(vH, vLabel);
        c->cd(6);
        fHM->Scale("fhRichMatchQaTofId", eventScale);
        HRichDrawHist::DrawH1(fHM->H1("fhRichMatchQaTofId"), kLinear, kLog);
    }
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRadius", "hrich_fhRadius", 1800, 1200);
    c->Divide(3, 2);
    c->cd(1);
    fHM->Scale("fhBetaVsRadiusSys0", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhBetaVsRadiusSys0"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Sys0", 0.4, 0.90, 0.6, 0.99);
    c->cd(2);
    fHM->Scale("fhBetaVsRadiusSys1", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhBetaVsRadiusSys1"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Sys1", 0.4, 0.90, 0.6, 0.99);
    c->cd(3);
    fHM->Scale("fhMomVsRadius", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhMomVsRadius"), kLinear, kLinear, kLog);
    c->cd(4);
    fHM->Scale("fhRichMatchQaVsRadiusCand", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRichMatchQaVsRadiusCand"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Candidates (isUsed)", 0.1, 0.90, 0.8, 0.99);
    c->cd(5);
    fHM->Scale("fhRichMatchQaVsRadiusTofId", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhRichMatchQaVsRadiusTofId"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.2, 0.90, 0.8, 0.99);
    c->cd(6);
    fHM->Scale("fhNofUsedFlagCandVsRadius", eventScale);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhNofUsedFlagCandVsRadius"));
    gPad->SetLogz();
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRadius2", "hrich_fhRadius2", 1200, 1200);
    fHM->Scale("fhThetaVsRadius", eventScale);
    fHM->Scale("fhThetaVsRadiusCand", eventScale);
    fHM->Scale("fhThetaVsRadiusTofId", eventScale);
    c->Divide(2,2);
    c->cd(1);
    vector<TH1*> vH;
    TH1* hPrY1 = (TH1*)fHM->H2("fhThetaVsRadius")->ProjectionY();
    TH1* hPrY2 = (TH1*)fHM->H2("fhThetaVsRadiusCand")->ProjectionY();
    TH1* hPrY3 = (TH1*)fHM->H2("fhThetaVsRadiusTofId")->ProjectionY();
    hPrY1->ResetStats(); // workaround, some problems with statistics
    hPrY2->ResetStats();
    hPrY3->ResetStats();
    hPrY1->GetYaxis()->SetTitle("Entries/event");
    vH.push_back(hPrY1);
    hPrY1->SetMinimum(0.9 * eventScale);
    vH.push_back(hPrY2);
    vH.push_back(hPrY3);
    vector<string> vLabel;
    vLabel.push_back("All (" + RichUtils::NumberToString(hPrY1->GetMean()) + ")");
    vLabel.push_back("Candidates (" + RichUtils::NumberToString(hPrY2->GetMean()) + ")");
    vLabel.push_back("TOF ID (" + RichUtils::NumberToString(hPrY3->GetMean()) + ")");
    HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog, true, 0.6, 0.89, 0.99, 0.99);
    c->cd(2);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsRadius"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("All rings", 0.4, 0.90, 0.6, 0.99);
    c->cd(3);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsRadiusCand"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("Candidates", 0.4, 0.90, 0.6, 0.99);
    c->cd(4);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsRadiusTofId"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.4, 0.90, 0.6, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRadius3", "hrich_fhRadius3", 1800, 600);
    c->Divide(3,1);
    c->cd(1);
    fHM->Scale("fhNofCalsInRingVsRadius", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhNofCalsInRingVsRadius"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("All rings", 0.4, 0.90, 0.6, 0.99);
    c->cd(2);
    fHM->Scale("fhNofCalsInRingVsRadiusCand", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhNofCalsInRingVsRadiusCand"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("Candidates", 0.4, 0.90, 0.6, 0.99);
    c->cd(3);
    fHM->Scale("fhNofCalsInRingVsRadiusTofId", eventScale);
    HRichDrawHist::DrawH2(fHM->H2("fhNofCalsInRingVsRadiusTofId"), kLinear, kLinear, kLog);
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.4, 0.90, 0.6, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhNofCalsInRing", "hrich_fhNofCalsInRing", 1200, 1200);
    fHM->Scale("fhThetaVsNofCalsInRing", eventScale);
    fHM->Scale("fhThetaVsNofCalsInRingCand", eventScale);
    fHM->Scale("fhThetaVsNofCalsInRingTofId", eventScale);
    c->Divide(2,2);
    c->cd(1);
    vector<TH1*> vH;
    TH1* hPrY1 = (TH1*)fHM->H2("fhThetaVsNofCalsInRing")->ProjectionY()->Clone();
    TH1* hPrY2 = (TH1*)fHM->H2("fhThetaVsNofCalsInRingCand")->ProjectionY()->Clone();
    TH1* hPrY3 = (TH1*)fHM->H2("fhThetaVsNofCalsInRingTofId")->ProjectionY()->Clone();
    hPrY1->ResetStats(); // workaround, some problems with statistics
    hPrY2->ResetStats();
    hPrY3->ResetStats();
    hPrY1->GetYaxis()->SetTitle("Entries/event");
    vH.push_back(hPrY1);
    hPrY1->SetMinimum(0.9 * eventScale);
    vH.push_back(hPrY2);
    vH.push_back(hPrY3);
    vector<string> vLabel;
    vLabel.push_back("All (" + RichUtils::NumberToString(hPrY1->GetMean()) + ")");
    vLabel.push_back("Candidates (" + RichUtils::NumberToString(hPrY2->GetMean()) + ")");
    vLabel.push_back("TOF ID (" + RichUtils::NumberToString(hPrY3->GetMean()) + ")");
    HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog);
    c->cd(2);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsNofCalsInRing"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("All rings", 0.4, 0.90, 0.6, 0.99);
    c->cd(3);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsNofCalsInRingCand"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("Candidates", 0.3, 0.90, 0.7, 0.99);
    c->cd(4);
    HRichDrawHist::DrawH2WithProfile(fHM->H2("fhThetaVsNofCalsInRingTofId"));
    gPad->SetLogz();
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.4, 0.90, 0.6, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhNofCalsInRingXY", "hrich_fhNofCalsInRingXY", 1800, 1200);
    c->Divide(3,2);
    c->cd(1);
    TH2D* pr1 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsXY"));
    pr1->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr1);
    HRichDrawHist::DrawTextOnPad("All rings", 0.4, 0.90, 0.6, 0.99);
    c->cd(2);
    TH2D* pr2 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsXYCand"));
    pr2->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr2);
    HRichDrawHist::DrawTextOnPad("Candidates", 0.3, 0.90, 0.7, 0.99);
    c->cd(3);
    TH2D* pr3 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsXYTofId"));
    pr3->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr3);
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.4, 0.90, 0.6, 0.99);
    c->cd(4);
    TH2D* pr4 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsThetaPhi"));
    pr4->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr4);
    HRichDrawHist::DrawTextOnPad("All rings", 0.4, 0.90, 0.6, 0.99);
    c->cd(5);
    TH2D* pr5 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsThetaPhiCand"));
    pr5->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr5);
    HRichDrawHist::DrawTextOnPad("Candidates", 0.3, 0.90, 0.7, 0.99);
    c->cd(6);
    TH2D* pr6 = HRichDrawHist::DrawH3Profile(fHM->H3("fhNofCalsInRingVsThetaPhiTofId"));
    pr6->GetZaxis()->SetRangeUser(8, 22);
    HRichDrawHist::DrawH2(pr6);
    HRichDrawHist::DrawTextOnPad("TOF ID", 0.4, 0.90, 0.6, 0.99);
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhCalsPhiThetaInterp", "hrich_fhCalsPhiThetaInterp", 1500, 500);
    c->Divide(3, 1);
    c->cd(1);
    fHM->Scale("fhCalsSectorsInterp", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhCalsSectorsInterp"));
    c->cd(2);
    fHM->Scale("fhCalsPhi", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhCalsPhi"));
    c->cd(3);
    fHM->Scale("fhCalsTheta", eventScale);
    HRichDrawHist::DrawH1(fHM->H1("fhCalsTheta"));
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhCalsThetaForSector", "hrich_fhCalsThetaForSector", 1500, 1000);
    c->Divide(3, 2);
    for (Int_t iS = 0; iS < 6; iS++) {
        c->cd(iS+1);
        stringstream ss;
        ss << "fhCalsThetaForSector" << iS;
        fHM->Scale(ss.str(), eventScale);
        HRichDrawHist::DrawH1(fHM->H1(ss.str()));
        stringstream ss2;
        ss2 << "Sector " << iS;
        HRichDrawHist::DrawTextOnPad(ss2.str(), 0.25, 0.9, 0.75, 0.99);
    }
    }

    {
    TCanvas* c = fHM->CreateCanvas("hrich_fhRaw", "hrich_fhRaw", 1800, 600);
    c->Divide(3,1);
    c->cd(1);
    HRichDrawHist::DrawH2(fHM->H2("fhRawPmtPixel"));
    c->cd(2);
    HRichDrawHist::DrawH1(fHM->H1("fhRawMult"));
    c->cd(3);
    HRichDrawHist::DrawH1(fHM->H1("fhRawToT"), kLinear, kLog);
    }

}

//---------------------------------------------------------------------------
Bool_t HRich700GeoTestRealQa::finalize()
{
    //drawHist(false);
    if (fOutputDir != "") {
        fHM->SaveCanvasToImage(string(fOutputDir + "/"));
    }

    if (fOutputRootFile != "") {
        TFile* file = new TFile(fOutputRootFile.c_str(), "RECREATE");
        fHM->WriteToFile();
        file->Close();
    }

    return kTRUE;
}

void HRich700GeoTestRealQa::drawFromFile(const string& fileName, const string& outputDir)
{
    TFile* file = new TFile(fileName.c_str(), "READ");

    fHM = new HRich700HistManager();
    fHM->ReadFromFile(file);

    drawHist(true);

    if (outputDir != "") {
        fHM->SaveCanvasToImage(string(outputDir + "/"));
    }

    file->Close();
}

void HRich700GeoTestRealQa::drawManyFromFile(const vector<string>& fileNames, const vector<string>& legNames, const string& outputDir)
{
    HRichDrawHist::SetDefaultDrawStyle();
    fHM = new HRich700HistManager(); // will be used to store canvases

    if (fileNames.size() != legNames.size()) {
        cout << "ERROR fileNames.size() != legNames.size()" << endl;
        return;
    }
    vector<HRich700HistManager*> hm;
    vector<TFile*> rootFiles;
    hm.resize(fileNames.size());
    rootFiles.resize(fileNames.size());
    for (UInt_t i = 0; i < hm.size(); i++) {
        rootFiles[i] = new TFile(fileNames[i].c_str(), "READ");
        hm[i] = new HRich700HistManager();
        hm[i]->ReadFromFile(rootFiles[i]);
    }

    {
        fHM->CreateCanvas("hrich_m_fhNofCalsInRingTofId", "hrich_m_fhNofCalsInRingTofId", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;
        for (UInt_t i = 0; i < hm.size(); i++){
            TH1* hPrY = (TH1*)hm[i]->H2("fhThetaVsNofCalsInRingTofId")->ProjectionY()->Clone();
            hPrY->Scale(1./hm[i]->H1("fhNofEvents")->GetEntries());
            hPrY->GetYaxis()->SetTitle("Entries/event");
            vH.push_back(hPrY);
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(hPrY->GetMean()) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }

    {
        fHM->CreateCanvas("hrich_m_fhNofCalsPerEvent", "hrich_m_fhNofCalsPerEvent", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;
        for (UInt_t i = 0; i < hm.size(); i++){
            hm[i]->H1("fhNofCals")->Scale(1./hm[i]->H1("fhNofEvents")->GetEntries());
            vH.push_back(hm[i]->H1("fhNofCals"));
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(hm[i]->H1("fhNofCals")->GetMean()) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }

    {
        fHM->CreateCanvas("hrich_m_fhRichMatchQaCand", "hrich_m_fhRichMatchQaCand", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;
        for (UInt_t i = 0; i < hm.size(); i++){
            hm[i]->H1("fhRichMatchQaCand")->Scale(1./hm[i]->H1("fhNofEvents")->GetEntries());
            vH.push_back(hm[i]->H1("fhRichMatchQaCand"));
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(hm[i]->H1("fhRichMatchQaCand")->GetMean()) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel, kLinear, kLog);
    }




    {
        fHM->CreateCanvas("hrich_m_fhRadiusTofId", "hrich_m_fhRadiusTofId", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;
        for (UInt_t i = 0; i < hm.size(); i++){
            TH1* hPrY = (TH1*)hm[i]->H2("fhThetaVsRadiusTofId")->ProjectionY()->Clone();
            hPrY->Scale(1./hm[i]->H1("fhNofEvents")->GetEntries());
            hPrY->GetYaxis()->SetTitle("Entries/event");
            vH.push_back(hPrY);
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(hPrY->GetMean()) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }

    {
        fHM->CreateCanvas("hrich_m_fhRichEffMom", "hrich_m_fhRichEffMom", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;

        for (UInt_t i = 0; i < hm.size(); i++){
            TH1D* pxEff = RichUtils::DivideH1((TH1D*)hm[i]->H1("fhMomRichTofId")->Clone(), (TH1D*)hm[i]->H1("fhMomTofId")->Clone(), "", 100., "Relative efficiency [%]");
            pxEff->SetMinimum(0.);
            pxEff->GetXaxis()->SetRangeUser(0, 500);
            vH.push_back(pxEff);
            stringstream ss1;
            ss1.precision(1);
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(100.*hm[i]->H1("fhMomRichTofId")->GetEntries() / hm[i]->H1("fhMomTofId")->GetEntries() ) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }

    {
        fHM->CreateCanvas("hrich_m_fhRichEffPhi", "hrich_m_fhRichEffPhi", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;

        for (UInt_t i = 0; i < hm.size(); i++){
            TH1D* pxEff = RichUtils::DivideH1((TH1D*)hm[i]->H2("fhThetaVsPhiRichTofId")->ProjectionY()->Clone(), (TH1D*)hm[i]->H2("fhThetaVsPhiTofId")->ProjectionY()->Clone(), "", 100., "Relative efficiency [%]");
            pxEff->SetMinimum(0.);
            vH.push_back(pxEff);
            stringstream ss1;
            ss1.precision(1);
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(100.*hm[i]->H1("fhThetaVsPhiRichTofId")->GetEntries() / hm[i]->H1("fhThetaVsPhiTofId")->GetEntries() ) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }

    {
        fHM->CreateCanvas("hrich_m_fhRichEffTheta", "hrich_m_fhRichEffTheta", 900, 900);
        vector<TH1*> vH;
        vector<string> vLabel;

        for (UInt_t i = 0; i < hm.size(); i++){
            TH1D* pxEff = RichUtils::DivideH1((TH1D*)hm[i]->H2("fhThetaVsPhiRichTofId")->ProjectionX()->Clone(), (TH1D*)hm[i]->H2("fhThetaVsPhiTofId")->ProjectionX()->Clone(), "", 100., "Relative efficiency [%]");
            pxEff->SetMinimum(0.);
            vH.push_back(pxEff);
            stringstream ss1;
            ss1.precision(1);
            vLabel.push_back(legNames[i] + " (" + RichUtils::NumberToString(100.*hm[i]->H1("fhThetaVsPhiRichTofId")->GetEntries() / hm[i]->H1("fhThetaVsPhiTofId")->GetEntries() ) + ")");
        }
        HRichDrawHist::DrawH1(vH, vLabel);
    }



    if (outputDir != "") {
        fHM->SaveCanvasToImage(string(outputDir + "/"));
    }
    // for (int i = 0; i < rootFiles.size(); i++) {
    //     rootFiles[i]->Close();
    // }
}