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

#include "hrich700digiparcreator.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 "hrichhitsim.h"
#include "htool.h"
#include "hhistconverter.h"
#include "hparticlecandfillerpar.h"
#include "hparticleanglecor.h"
#include "hrich700pmttypeenum.h"

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

#include "hgeantkine.h"

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

#include <iostream>
#include <sstream>
#include <fstream>

using namespace std;


ClassImp(HRich700DigiParCreator)

HRich700DigiParCreator::HRich700DigiParCreator():
fEventNum(0),
fOutputHistFilePath("")
{

}

HRich700DigiParCreator::~HRich700DigiParCreator()
{

}

Bool_t HRich700DigiParCreator::init()
{
    fCatKine = gHades->getCurrentEvent()->getCategory(catGeantKine);
    if (NULL == fCatKine) {
    	Error("init", "Initializatin of kine category failed, returning...");
    	return kFALSE;
    }

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

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

    fRingFinderPar = (HRich700RingFinderPar*) gHades->getRuntimeDb()->getContainer("Rich700RingFinderPar");
    if(fRingFinderPar == NULL) {
        Error("init", "Can not retrieve fRingFinderPar");
        return kFALSE;
    }

    fCandFillerPar =(HParticleCandFillerPar*) gHades->getRuntimeDb()->getContainer("ParticleCandFillerPar");
    if(fCandFillerPar == NULL) {
	Error ("init", "Can not retrieve fCandFillerPar");
	return kFALSE;
    }

    initHist();
    initZVertex();
    return kTRUE;
}

void HRich700DigiParCreator::initZVertex()
{
    // IMPORTANT!!!
    // the values should correspond to simulated Z vertices
    // see run_ascii.C macro
    for (int i = 0; i < knZ; i++) {
       kZ[i] = -25. * i;
    }
}

void HRich700DigiParCreator::initHist()
{
    fHM = new HRich700HistManager();

    //XY[mm]
    const Double_t hMinXY = -660.;
    const Double_t hMaxXY = 660.;
    //XY pmtXInd, pmtYInd
    const Double_t hMinInd = -0.5;
    const Double_t hMaxInd = 23.5;
    const Int_t nofBinsInd = 24;
    for (int iZ = 0; iZ < knZ; iZ++) {
	fHM->AddTNamedObject(new TProfile2D(Form("fpPhiXY[%i]", iZ), Form("fpPhiXY[%i];X [mm];Y [mm]", iZ), knXY, hMinXY, hMaxXY, knXY, hMinXY, hMaxXY, 0., 360.));
        fHM->P2(Form("fpPhiXY[%i]", iZ))->SetErrorOption("s");
	fHM->AddTNamedObject( new TProfile2D(Form("fpThetaXY[%i]", iZ), Form("fpThetaXY[%i];X [mm];Y [mm]", iZ), knXY, hMinXY, hMaxXY, knXY, hMinXY, hMaxXY, 0., 90.));
        fHM->P2(Form("fpThetaXY[%i]", iZ))->SetErrorOption("s");

	fHM->AddTNamedObject(new TProfile2D(Form("fpPhiXYInd[%i]", iZ), Form("fpPhiXYInd[%i];PMT X_{ind};PMT Y_{ind};", iZ), nofBinsInd, hMinInd, hMaxInd, nofBinsInd, hMinInd, hMaxInd, 0., 360.));
        fHM->P2(Form("fpPhiXYInd[%i]", iZ))->SetErrorOption("s");
	fHM->AddTNamedObject(new TProfile2D(Form("fpThetaXYInd[%i]", iZ), Form("fpThetaXYInd[%i];PMT X_{ind};PMT Y_{ind};", iZ), nofBinsInd, hMinInd, hMaxInd, nofBinsInd, hMinInd, hMaxInd, 0., 90.));
        fHM->P2(Form("fpThetaXYInd[%i]", iZ))->SetErrorOption("s");
    }

    fHM->Create1<TH1F>("fhNofKineInEvent", "fhNofKineInEvent;Nof Kine in event;Counter", 150, -0.5, 149.5);
    fHM->Create1<TH1F>("fhNofRichRingsInEvent", "fhNofRichRingsInEvent;Nof Kine in event;Counter", 150, -0.5, 149.5);

    fHM->Create2<TH2F>("fhXYVertexKine", "fhXYVertexKine;X [mm];Y[mm];Counter", 200, -1, 1, 200, -1, 1);
    fHM->Create1<TH1F>("fhZVertexKine", "fhZVertexKine;Z [mm];Counter", 90, -80, 10);
    fHM->Create1<TH1F>("fhMomKine", "fhMomKine;Momentum [MeV/c];Counter", 101, 0., 1010.);
    fHM->Create2<TH2F>("fhThetaPhiKine", "fhThetaPhiKine;Theta [deg];Phi [deg];Counter", 90, 0., 90., 360, 0., 360.);

    fHM->Create2<TH2F>("fhZIndVsVz", "fhZIndVsVz;Z_{ind};Z_{vertex} [mm];Counter", knZ, -0.5, knZ - 0.5, 100, -100., 0.);
}

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


Int_t HRich700DigiParCreator::execute()
{
    HRichDrawHist::SetDefaultDrawStyle();
    fEventNum++;
    if (fEventNum % 10000 == 0) cout << "HRich700DigiParCreator::execute eventNum " << fEventNum << endl;
    processEvent();


    return 0;
}

Int_t HRich700DigiParCreator::getZIndex(double z)
{
    for (Int_t i = 0; i < knZ; i++) {
	if (std::abs(kZ[i] - z) < 0.01) return i;
    }
    return 0;
}

string HRich700DigiParCreator::getStringByZIndex(Int_t zInd)
{
    if (zInd < 0 || zInd >= knZ) return "";
    stringstream ss;
    ss << "Z_{vertex} = "<< kZ[zInd] << " mm";
    return ss.str();
}


void HRich700DigiParCreator::processEvent()
{
    Int_t nofKine = fCatKine->getEntries();
    Int_t nofRichHits = fCatRichHit->getEntries();

    fHM->H1("fhNofKineInEvent")->Fill(nofKine);
    fHM->H1("fhNofRichRingsInEvent")->Fill(nofRichHits);

    if (nofKine >= 4 || nofRichHits >= 4) return;

    for (Int_t i = 0; i < nofKine; i++) {
        HGeantKine* kine = (HGeantKine*)fCatKine->getObject(i);
        if (!isPrimaryElectron(kine)) continue;
        Float_t vx, vy, vz;
        kine->getVertex(vx, vy, vz);
        fHM->H2("fhXYVertexKine")->Fill(vx, vy);
        fHM->H1("fhZVertexKine")->Fill(vz);
        fHM->H1("fhMomKine")->Fill(kine->getTotalMomentum());
        fHM->H2("fhThetaPhiKine")->Fill(kine->getThetaDeg(), kine->getPhiDeg());
    }


    
    for (Int_t i = 0; i < nofRichHits; i++) {
    	HRichHitSim* richHit = (HRichHitSim*) fCatRichHit->getObject(i);
	if (richHit == NULL
	    || richHit->fRich700CircleRadius < 12. || richHit->fRich700CircleRadius > 32.
	    || richHit->fRich700NofRichCals < 3 || richHit->fRich700NofRichCals > 35) continue;

        HGeantKine* kine = (HGeantKine*)fCatKine->getObject(richHit->track1 - 1);
        // select only rings from primary electrons
    	if (!isPrimaryElectron(kine)) continue;

	Float_t vx, vy, vz;
	kine->getVertex(vx, vy, vz);
	Int_t iZ = getZIndex(vz);
	fHM->H2("fhZIndVsVz")->Fill(iZ, vz);

         // x,y [mm]
	Double_t xR = richHit->fRich700CircleCenterX;
	Double_t yR = richHit->fRich700CircleCenterY;

	Float_t theta = 0, phi = 0;

        // remove bad rings
        if ((phi < 90 && yR < -20.) || (phi > 270 && yR > 20)){
            continue;
	}

	// correct phi for 0-360 deg problem
	if ((phi < 90 && yR < 0) || (phi > 270 && yR > 0)){
	    yR = -yR;
	}

	//PMT
	Int_t pmtId = fDigiPar->getPMTId((Float_t)xR, yR);
	if(pmtId >= 0){
	    HRich700PmtData* rd = fDigiPar->getPMTData(pmtId);
	    fHM->P2(Form("fpPhiXYInd[%i]", iZ))->Fill(rd->fIndX, rd->fIndY, phi);
	    fHM->P2(Form("fpThetaXYInd[%i]", iZ))->Fill(rd->fIndX, rd->fIndY, theta);
	}

	fHM->P2(Form("fpPhiXY[%i]", iZ))->Fill(xR, yR, phi);
	fHM->P2(Form("fpThetaXY[%i]", iZ))->Fill(xR, yR, theta);
    }
}

void HRich700DigiParCreator::createOutputFile()
{
    string rfStr = fRingFinderPar->getStringForParTxtFile();
    string digiStr = fDigiPar->getStringForParTxtFile();

    string digiPmtIdStr = createStringPhiThetaPmtId();
    string digiXYStr = createStringPhiThetaXY();
    string thetaPolyStr = createStringThetaTransParamsPoly();

    ofstream fout;
    fout.open(fOutputFilePath.c_str());
    fout << rfStr << digiStr << digiPmtIdStr << digiXYStr << thetaPolyStr;
    fout << "##############################################################################" <<endl;
    fout.close();

}

string HRich700DigiParCreator::createStringPhiThetaPmtId()
{
    cout <<"createStringPhiThetaPmtId" << endl;
    ostringstream ss;
    ss <<"fArrayPmt:  Double_t \\" << endl;
    Int_t iZ = 0;

    // WLS pmt for mar19 beamtime
    const Int_t nofWlsIds = 48;
    Int_t wlsIds[nofWlsIds] = {204, 205, 225, 226, 227, 228, 229, 230, 231, 232, 249, 250, 251, 252, 253, 254, 255, 256,
    273, 274, 275, 278, 279, 280, 297, 298, 299, 302, 303, 304, 321, 322, 323, 324, 325, 326, 327, 328,
    345, 346, 347, 348, 349, 350, 351, 352, 372, 373};

    for(Int_t i = 0; i < knPmt; i++){
	HRich700PmtData* pd = fDigiPar->getPMTData(i);
	if(pd == NULL) continue;
	Double_t phiMean = fHM->P2(Form("fpPhiXYInd[%i]", iZ))->GetBinContent(pd->fIndX + 1, pd->fIndY + 1);
	Double_t thetaMean = fHM->P2(Form("fpThetaXYInd[%i]", iZ))->GetBinContent(pd->fIndX + 1, pd->fIndY + 1);
	Int_t pmtType = pd->fPmtType;
	for (int i = 0; i < nofWlsIds; i++) {
	    if (wlsIds[i] == pd->fPmtId){
		pmtType = HRich700PmtTypeCosy17WithWls;
                break;
	    }
	}
	ss << pd->fPmtId << " " << pd->fIndX << " " << pd->fIndY << " " << pd->fX << " "
	    << pd->fY << " " << pd->fZ << " " << pmtType << " " << thetaMean << " "
	    << phiMean << " \\" << endl;
    }
    // remove last " \" and endl symbols
    string str = ss.str().substr(0, ss.str().size() - 3);
    str = string(str + "\n");
    //cout << str << endl;

    return str;
    //fclose(pFile);
}

string HRich700DigiParCreator::createStringPhiThetaXY()
{
    cout << "createStringPhiThetaXY" << endl;
    stringstream ss;
    TArrayD zAr;
    zAr.Set(knZ, kZ);
    HHistConverter::writeArray(ss, "fArrayZVertices", zAr, 10);

    for (Int_t i = 0; i < knZ; i++) {
	//TArrayD linDataPhi;
	TArrayD linDataTheta;
	//HHistConverter::fillArray(fHM->H2(Form("fpPhiXY[%i]",i)), linDataPhi, "Form(fArrayPhiMean[%i],i)", 10, 10, kFALSE);
	HHistConverter::fillArray(fHM->H2(Form("fpThetaXY[%i]",i)), linDataTheta, Form("fArrayThetaMean[%i]", i), 10, 10, kFALSE);

	//HHistConverter::writeArray(ss, Form("fArrayPhiMean[%i]", i), linDataPhi, 10);
	HHistConverter::writeArray(ss, Form("fArrayThetaMean[%i]", i), linDataTheta, 10);
    }
    return ss.str();
}

string HRich700DigiParCreator::createStringThetaTransParamsPoly()
{
    cout << "createStringThetaTransParamsPoly" << endl;
    stringstream ss;

    // The parameters for analytical X,Y -> Theta transformation developed by Joerg Foertsch (Wuppertal Uni)
    Double_t polyParams[16][5] = {
	{+1.1369326537813567e+00, -1.0632571708885562e-02, +3.7328828043546198e-05, -5.8129388919359630e-08, +3.3736510274371330e-11},
	{+1.6921292491301611e-02, +1.1367376718416013e-03, -4.0964484469121487e-06, +6.4467554002406032e-09, -3.7779327749509494e-12},
	{+3.0432354338100452e-03, -2.9108617521613769e-05, +1.0459414484497346e-07, -1.6669710778113849e-10, +9.8925641913421989e-14},
	{-2.9670946160677261e-05, +2.8765294196925021e-07, -1.0499767058846625e-09, +1.7028792672246029e-12, -1.0285451414655484e-15},
	{+1.2879373308389344e-07, -1.2951920086928232e-09, +4.8875634831362772e-12, -8.1610531442154027e-15, +5.0545063006800175e-18},
	{-2.5055230014004777e-10, +2.6791881754382908e-12, -1.0687796814869330e-14, +1.8742115249294953e-17, -1.2103612918118934e-20},
	{+1.7074396225102669e-13, -2.1651432517492696e-15, +9.7827195539292294e-18, -1.8858059984766831e-20, +1.3098563438069365e-23},
	{+8.8070344305834701e-18, +2.7528980557118093e-19, -2.3366811003398705e-21, +5.9307162038263844e-24, -4.8215115873206042e-27},
	{+1.0255107536596604e+04, -8.6480386267682249e+01, +2.7013516608647020e-01, -3.6897803501428233e-04, +1.8506083880516984e-07},
	{-2.3507656748940855e+02, +1.9916941257550129e+00, -6.2518906742161460e-03, +8.5904274471498453e-06, -4.3406788029937180e-09},
	{+2.2254320102725140e+00, -1.8918298862692612e-02, +5.9674149619642513e-05, -8.2491319710548530e-08, +4.1994575592520704e-11},
	{-1.1234339502120968e-02, +9.5897773382736873e-05, -3.0400501529305233e-07, +4.2280052379232953e-10, -2.1683822369569551e-13},
	{+3.2649936539464306e-05, -2.7985172625980123e-07, +8.9157752086558108e-10, -1.2474411701236813e-12, +6.4441703761595720e-16},
	{-5.4624674421154949e-08, +4.7002418297894819e-10, -1.5045582409448433e-12, +2.1172251767103056e-15, -1.1013422467825424e-18},
	{+4.8764885686986162e-11, -4.2109646679116353e-13, +1.3538273251623813e-15, -1.9152849190112116e-18, +1.0027243308189633e-21},
	{-1.7958902985195771e-14, +1.5557460663233240e-16, -5.0213593687252007e-19, +7.1380360928510676e-22, -3.7588775710427317e-25}
    };
    TArrayD polyAr;
    polyAr.Set(16 * 5);
    Int_t i = 0;
    ss << std::scientific;
    ss.precision(16);
    for (Int_t i1 = 0; i1 < 16; i1++) {
	for (Int_t i2 = 0; i2 < 5; i2++) {
	    polyAr.SetAt(polyParams[i1][i2], i);
            i++;
	}
    }
    HHistConverter::writeArray(ss, "fArrayThetaTransParamsPoly", polyAr, 5);

    // These geometry parametrs are used for the transformation algorithm [mm]
    const Double_t innerPlaneWidth = 371.;
    const Double_t meanTargetWRTMirror = 533.3; //This is actually more like the hydra z=0 in the mirror coordinate system
    const Double_t TargetCenter = -40.; //This is the target center in the hydra coordinate system
    const Double_t mirrorRadius = 872.;
    const Double_t PMTPlane1Z = 533.3 - 122.2;
    const Double_t PMTPlane2Z = 533.3 - 122.2 - 128.;
    Double_t geoParams[6];
    geoParams[0] = innerPlaneWidth;
    geoParams[1] = meanTargetWRTMirror;
    geoParams[2] = TargetCenter;
    geoParams[3] = mirrorRadius;
    geoParams[4] = PMTPlane1Z;
    geoParams[5] = PMTPlane2Z;

    ss << std::fixed;
    ss.precision(5);
    TArrayD geoAr;
    geoAr.Set(6, geoParams);
    HHistConverter::writeArray(ss, "fArrayThetaTransParamsGeo", geoAr, 10);

    return ss.str();
}

Bool_t HRich700DigiParCreator::isPrimaryElectron(
					    HGeantKine* kine)
{
    if (kine == NULL) return kFALSE;
    return (kine->getID() == 2 || kine->getID() == 3) && kine->isPrimary();
}

void  HRich700DigiParCreator::convertProfile(const string& profileName, const string& zTitle)
{
    TProfile2D* profile = fHM->P2(profileName);
    TH2D* rms = profile->ProjectionXY((string(profile->GetName()) + "_RMS").c_str(), "C=E");
    rms->GetXaxis()->SetTitle(profile->GetXaxis()->GetTitle());
    rms->GetYaxis()->SetTitle(profile->GetYaxis()->GetTitle());
    rms->GetZaxis()->SetTitle(("RMS. " + zTitle).c_str());
    fHM->AddTNamedObject(rms);

    TH2D* counter = profile->ProjectionXY((string(profile->GetName()) + "_Counter").c_str(), "B");
    counter->GetXaxis()->SetTitle(profile->GetXaxis()->GetTitle());
    counter->GetYaxis()->SetTitle(profile->GetYaxis()->GetTitle());
    counter->GetZaxis()->SetTitle(("Counter. " + zTitle).c_str());
    fHM->AddTNamedObject(counter);

    TH2D* mean = profile->ProjectionXY((string(profile->GetName()) + "_Mean").c_str());
    mean->GetXaxis()->SetTitle(profile->GetXaxis()->GetTitle());
    mean->GetYaxis()->SetTitle(profile->GetYaxis()->GetTitle());
    mean->GetZaxis()->SetTitle(("Mean. " + zTitle).c_str());
    fHM->AddTNamedObject(mean);
}

void  HRich700DigiParCreator::convertProfileAll()
{
    for (int iZ = 0; iZ < knZ; iZ++) {
	convertProfile(Form("fpPhiXY[%i]", iZ), "Phi [deg]");
	convertProfile(Form("fpThetaXY[%i]", iZ), "Theta [deg]");
	convertProfile(Form("fpPhiXYInd[%i]", iZ), "Phi [deg]");
	convertProfile(Form("fpThetaXYInd[%i]", iZ), "Theta [deg]");
    }
}

Bool_t HRich700DigiParCreator::finalize()
{
    TFile* out = new TFile(fOutputHistFilePath.c_str(), "RECREATE");
    out->cd();
    fHM->WriteToFile();
    out->Close();
    return kTRUE;
}

void HRich700DigiParCreator::drawZDiffFromRef( TCanvas* c, const string& histNameTemplate, Int_t indMin, Int_t indMax, Int_t refInd,
					      Double_t zMin, Double_t zMax )
{
    TH2D* diff0 = (TH2D*)fHM->H2(Form(histNameTemplate.c_str(), refInd))->Clone();
    for (int iZ = indMin; iZ < indMax; iZ++){
	c->cd(iZ+1);
	TH2D* diff1 = (TH2D*)fHM->H2(Form(histNameTemplate.c_str(), iZ))->Clone();
	diff1->Add(diff0, -1.);
	diff1->GetZaxis()->SetRangeUser(zMin, zMax);
	HRichDrawHist::DrawH2(diff1);
	HRichDrawHist::DrawTextOnPad(getStringByZIndex(iZ), 0.1, 0.9, 0.7, 0.99);
    }
}

void HRich700DigiParCreator::drawHist()
{
    HRichDrawHist::SetDefaultDrawStyle();

    for (Int_t iZ = 0; iZ < knZ; iZ++) {
	TCanvas* c = fHM->CreateCanvas(Form("hrichdigipar_XY_%imm", (Int_t)kZ[iZ]), Form("hrichdigipar_XY_%imm", (Int_t)kZ[iZ]), 2000, 1200);
	c->Divide(3,2);
	c->cd(1);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXY[%i]_Mean",iZ)));
        drawPmts();
	c->cd(2);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXY[%i]_Mean",iZ)));
        drawPmts();
	c->cd(3);
        fHM->H2(Form("fpPhiXY[%i]_RMS",iZ))->GetZaxis()->SetRangeUser(0., 2.);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXY[%i]_RMS",iZ)));
	c->cd(4);
        fHM->H2(Form("fpThetaXY[%i]_RMS",iZ))->GetZaxis()->SetRangeUser(0., 2.);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXY[%i]_RMS",iZ)));
	c->cd(5);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXY[%i]_Counter",iZ)));
	c->cd(6);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXY[%i]_Counter",iZ)));
    }

    for (Int_t iZ = 0; iZ < knZ; iZ++) {
	TCanvas* c = fHM->CreateCanvas(Form("hrichdigipar_XYInd_%imm",(Int_t)kZ[iZ]), Form("hrichdigipar_XYInd_%imm",(Int_t)kZ[iZ]), 2000, 1200);
	c->Divide(3,2);
	c->cd(1);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXYInd[%i]_Mean",iZ)));
	c->cd(2);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXYInd[%i]_Mean",iZ)));
	c->cd(3);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXYInd[%i]_RMS",iZ)));
	c->cd(4);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXYInd[%i]_RMS",iZ)));
	c->cd(5);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpPhiXYInd[%i]_Counter",iZ)));
	c->cd(6);
	HRichDrawHist::DrawH2(fHM->H2(Form("fpThetaXYInd[%i]_Counter",iZ)));
    }


    {
	TCanvas* c = fHM->CreateCanvas("hrichdigipar_XY_Phi_Diff", "hrichdigipar_XY_Phi_Diff", 1500, 1200);
	c->Divide(2,2);
        drawZDiffFromRef(c,"fpPhiXY[%i]_Mean", 0, knZ, 0, -1., 1.);
    }

    {
	TCanvas* c = fHM->CreateCanvas("hrichdigipar_XY_Theta_Diff", "hrichdigipar_XY_Theta_Diff", 1500, 1200);
	c->Divide(2,2);
        drawZDiffFromRef(c,"fpThetaXY[%i]_Mean", 0, knZ, 0, -3., 3.);
    }

    {
	TCanvas* c = fHM->CreateCanvas("hrichdigipar_KineParams", "hrichdigipar_KineParams", 1200, 1200);
	c->Divide(2,2);
	c->cd(1);
	HRichDrawHist::DrawH2(fHM->H2("fhXYVertexKine"));
	c->cd(2);
	HRichDrawHist::DrawH1(fHM->H1("fhZVertexKine"));
	c->cd(3);
	fHM->H1("fhMomKine")->SetMinimum(0.);
	HRichDrawHist::DrawH1(fHM->H1("fhMomKine"));
	c->cd(4);
	HRichDrawHist::DrawH2(fHM->H2("fhThetaPhiKine"));
    }

    {
	fHM->CreateCanvas("hrichdigipar_ZIndVsVz", "hrichdigipar_ZIndVsVz", 900, 900);
	HRichDrawHist::DrawH2(fHM->H2("fhZIndVsVz"));
    }

    {
	TCanvas* c = fHM->CreateCanvas("hrichdigipar_NofObjectsInEvent", "hrichdigipar_NofObjectsInEvent", 1500, 800);
	c->Divide(2,1);
	c->cd(1);
	HRichDrawHist::DrawH1(fHM->H1("fhNofKineInEvent"), kLinear, kLog);
	c->cd(2);
	HRichDrawHist::DrawH1(fHM->H1("fhNofRichRingsInEvent"), kLinear, kLog);
    }
}

void HRich700DigiParCreator::testZInterpolationAll()
{
    HRichDrawHist::SetDefaultDrawStyle();
    {
	vector<int> z; z.push_back(0); z.push_back(7); z.push_back(14);
	testZInterpolation("fpPhiXY[%i]_Mean", "fhIntTestXYdPhiZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdPhi_Z_0_7_14", z);
        testZInterpolation("fpThetaXY[%i]_Mean", "fhIntTestXYdThetaZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdTheta_Z_0_7_14", z);
    }

    {
	vector<int> z; z.push_back(0); z.push_back(3); z.push_back(6); z.push_back(9); z.push_back(12); z.push_back(14);
	testZInterpolation("fpPhiXY[%i]_Mean", "fhIntTestXYdPhiZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdPhi_Z_0_3_6_9_12_14", z);
        testZInterpolation("fpThetaXY[%i]_Mean", "fhIntTestXYdThetaZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdTheta_Z_0_3_6_9_12_14", z);
    }

    {
	vector<int> z; z.push_back(0); z.push_back(2); z.push_back(4); z.push_back(6); z.push_back(8); z.push_back(10); z.push_back(12); z.push_back(14);
	testZInterpolation("fpPhiXY[%i]_Mean", "fhIntTestXYdPhiZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdPhi_Z_0_2_4_6_8_10_12_14", z);
        testZInterpolation("fpThetaXY[%i]_Mean", "fhIntTestXYdThetaZ0[%i]Z1[%i]_Z[%i]", "hrichdigipar_IntTestXYdTheta_Z_0_2_4_6_8_10_12_14", z);
    }
}

void HRich700DigiParCreator::testZInterpolation(const string& hTemplate,
                                                const string& hResTemplate,
						const string& cName,
						const vector<int>& indZVec)
{
    TH2* hPar = fHM->H2(Form(hTemplate.c_str(), 0));
    int nBinsX = hPar->GetNbinsX();
    int nBinsY = hPar->GetNbinsY();
    double xMin = hPar->GetXaxis()->GetXmin();
    double xMax = hPar->GetXaxis()->GetXmax();
    double yMin = hPar->GetXaxis()->GetXmin();
    double yMax = hPar->GetXaxis()->GetXmax();

    TCanvas* c = fHM->CreateCanvas(cName.c_str(), cName.c_str(), 2400, 1000);
    c->Divide(5, 3);
    for (UInt_t iV = 0; iV < indZVec.size() - 1; iV++) {
	int indZ0 = indZVec[iV];
        int indZ1 = indZVec[iV+1];
	double z0 = kZ[indZ0];
	double z1 = kZ[indZ1];
	for (int iZ = indZ0 + 1; iZ < indZ1; iZ++) {
	    double z = kZ[iZ];
	    fHM->Create2<TH2D>(Form(hResTemplate.c_str(), indZ0, indZ1, iZ), Form((hResTemplate + ";X [mm];Y [mm];Residual [deg]").c_str(), indZ0, indZ1, iZ),
			       nBinsX, xMin, xMax, nBinsY, yMin, yMax);
	    TH2* hRes = fHM->H2(Form(hResTemplate.c_str(), indZ0, indZ1, iZ));
	    TH2* hZ0 = fHM->H2(Form(hTemplate.c_str(),indZ0));
	    TH2* hZ1 = fHM->H2(Form(hTemplate.c_str(),indZ1));
            TH2* hZ = fHM->H2(Form(hTemplate.c_str(),iZ));

	    for (int xB = 1; xB <= nBinsX; xB++) {
		for (int yB = 1; yB <= nBinsY; yB++) {
		    double v0 = hZ0->GetBinContent(xB, yB);
		    double v1 = hZ1->GetBinContent(xB, yB);
		    double v = v0 + (z - z0) * (v1 - v0) / (z1 - z0);
		    double dV = v - hZ->GetBinContent(xB, yB);
		    hRes->SetBinContent(xB, yB, dV);
		}
	    }
	    c->cd(iZ+1);
	    hRes->GetZaxis()->SetRangeUser(-0.2, 0.2);
	    HRichDrawHist::DrawH2(hRes);
	}
    }
    for (int iZ = 0; iZ < knZ; iZ++) {
        c->cd(iZ+1);
	HRichDrawHist::DrawTextOnPad(getStringByZIndex(iZ), 0.1, 0.9, 0.7, 0.99);
    }


}

void HRich700DigiParCreator::drawPmts()
{
    Double_t pmtHalfSize = fDigiPar->getPmtSize() / 2.;
    map<Int_t,HRich700PmtData> pmtDataMap = fDigiPar->getPmtDataMapPmtId();
    for(map<Int_t,HRich700PmtData>::iterator it = pmtDataMap.begin(); it != pmtDataMap.end(); it++) {
        Double_t pmtX = it->second.fX;
	Double_t pmtY = it->second.fY;
        TBox* box2 = new TBox(pmtX - pmtHalfSize, pmtY - pmtHalfSize, pmtX + pmtHalfSize, pmtY + pmtHalfSize);
        box2->SetLineColor(kBlack);
        box2->SetFillStyle(0);
        box2->SetLineWidth(1);
        box2->DrawClone();
    }
}

void HRich700DigiParCreator::drawFromFile( const string& fileName, const string& outputDir)
{
    cout << "drawFromFile histFile:" << fileName << " outputDir:" << outputDir << endl;
    // this we need to draw PMTs
    fDigiPar = (HRich700DigiPar*) gHades->getRuntimeDb()->getContainer("Rich700DigiPar");

    initZVertex();
    fHM = new HRich700HistManager();
    TFile* file = new TFile(fileName.c_str());
    fHM->ReadFromFile(file);

    convertProfileAll();
    createOutputFile();

    gStyle->SetOptTitle(1);

   // testZInterpolationAll();
    drawHist();

    fHM->SaveCanvasToImage(string(outputDir + "/"));
}