//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//     Class to align tpc with respect to external tracks
//
//
// Environment:
//      Software NOT developed for the PANDA Detector at FAIR.
//
// Author List:
//      Alexander Schmah      TUM            (original author)
//      Maxence Vandenbroucke TUM            (author)
//      Francesco Cusanno     TUM            (author)
//      Sverre Doerheim       TUM            (author, rewrite from macro to
//      compiled object)
//
//
//-----------------------------------------------------------
#include <stdlib.h>
#include <omp.h>
#include <sstream>
#include <fstream>
#include <iomanip>


#include <TH1D.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TProfile.h>
#include <TMinuit.h>
#include <TString.h>
#include <TFile.h>
#include <TTree.h>
#include <TStopwatch.h>
#include <TTreeCache.h>

#include "GFTrack.h"
#include "GFFieldManager.h"
#include "RKTrackRep.h"
#include "GFException.h"

#include "FairField.h"
#include "FairRunAna.h"

#include "PndFieldAdaptor.h"

#include "TpcStandaloneResCalc.h"
#include "TpcRefGFTrkResCalc.h"
#include "TpcRefGFTrkTpcClusterResCalc.h"
#include "TpcCdcFit2DResCalc_Alignment.h"
#include "TpcAlignmentHelper.h"
#include "TpcEventIdentifier.h"

#include "CdcHit.h"
#include "CdcTrack.h"

extern double getPhiInSector(double phi);
TpcStandaloneResCalc::TpcStandaloneResCalc()
    : fVerbose(0),
      fDetectorName("tpc"),
      fAlignmentFile("alinput.txt"),
      fResCalcName(""),
      fMaxEvents(0),
      fChargeFiltering(0),
      fOutFileName("testResiduals.root"),
      fSaveHistos(false),
      fHistoFileName("AlignmentHists.root"),
      fzSlicesHists(10),
      fMaxRangeXY(15),
      fMaxRangeZ(40),
      fzStart(-25),
      fzEnd(30),
      fFopiEventBranchName(""),
      fTpcEventIdentifierBranchName(""){
    fAlMan = NULL;
    fInputChain = NULL;
    fFopiTree = NULL;
    fResidualCalculator = NULL;
}

void TpcStandaloneResCalc::setTpcEventBranchName(const TString name)
{
    fTpcEventIdentifierBranchName=name;
    fBranches[name]=new TClonesArray("TpcEventIdentifier");
}

void TpcStandaloneResCalc::setFopiEventBranchName(const TString name)
{
    fFopiEventBranchName=name;
    fFopiBranches[name]=new TClonesArray("FopiEvent");
}

void TpcStandaloneResCalc::addBranch(TString name, TString classType) {
    std::cout<<"Adding branch "<<name<<", class: "<<classType<<std::endl;
    fBranches[name] = new TClonesArray(classType);
}
void TpcStandaloneResCalc::addFopiBranch(TString name, TString classType) {
    std::cout<<"Adding fopi branch "<<name<<", class: "<<classType<<std::endl;
    fFopiBranches[name] = new TClonesArray(classType);
}
bool TpcStandaloneResCalc::init() {
    if (fResidualCalculator == NULL) {
        std::cerr << "TpcStandaloneResCalc::init residual calculator missing"
                  << std::endl;
        return false;
    }
    if (fInputChain == NULL) {
        std::cerr << "TpcStandaloneResCalc::init() \n "
                  << "Input chain missing, please execuite setInputChain(chain) "
                     "before running init" << std::endl;
        return false;
    }
    
    fInputChain->SetCacheSize(100e6);
    TTreeCache::SetLearnEntries(10);
    
    TpcAlignmentManager::init(fAlignmentFile.c_str());
    fAlMan = TpcAlignmentManager::getInstance();
    if (fAlMan == 0) {
        std::cerr << "TpcStandaloneResCalc::init() Init TpcAlignmentManager:"
                  << fAlignmentFile.c_str() << " not found!" << std::endl;
        return false;
    }
    fAlMan->write("alignment_rescalc.txt");
    // TpcAlignmentHelper a;
    
    std::map<TString, TClonesArray*>::iterator it1;
    
    for (it1 = fBranches.begin(); it1 != fBranches.end(); ++it1) {
        std::cout << "setting branch address for branch: "<<it1->first<<std::endl;
        int test = fInputChain->SetBranchAddress(it1->first, &(it1->second));
        if (test < 0) {
            std::cerr << "TpcStandaloneResCalc::init() Branch " << it1->first
                      << " not found " << test << std::endl;
            return 1;
        }
    }
    
    for (it1 = fFopiBranches.begin(); it1 != fFopiBranches.end(); ++it1) {
        int test = fFopiTree->SetBranchAddress(it1->first, &(it1->second));
        if (test < 0) {
            std::cerr << "TpcStandaloneResCalc::init() Branch " << it1->first
                      << " not found " << test << std::endl;
            return 1;
        }
    }
    std::cout << "adding branches to residual calculator branch-map" << std::endl;
    // Adding Branches to the residual calculator
    std::map<TString, TClonesArray*>* branches =
            fResidualCalculator->getBranchMap();
    std::map<TString, TClonesArray*>::iterator it;
    std::cout << "looking for " << branches->size() << " branches" << std::endl;
    
    for (it = branches->begin(); it != branches->end(); ++it) {
        TString branchName = it->first;
        TClonesArray* br=NULL;
        if(fBranches.count(branchName)>0){
            br = fBranches[branchName];
            std::cout << it->first << std::endl;
            std::cout<<"fInputChain->SetBranchStatus("<<branchName<<"+\".*\",1);"<<std::endl;
            fInputChain->SetBranchStatus(branchName+".*",1);
        }  else{
            std::cout<<"Branch not in main tree, trying to find it in fopi tree "<<branchName<<std::endl;
            br=fFopiBranches[branchName];
            if(br==NULL){
                std::string s("Branch ");
                s.append(it->first.Data());
                s.append(" requested by residual calculator, but not found in tree!");
                std::cerr << "TpcStandaloneResCalc::Init() " << s << std::endl;
                return 1; 
            }
            fFopiTree->SetBranchStatus(branchName+".*",1);
        }
        it->second = br;
        
    }
    // Initializing, searching for set branches
    std::cout << "TpcStandaloneResCalc:init() initializing residual calculator: "
              << std::endl;
    std::cout << "Self awareness ... " << std::endl;
    if (fResidualCalculator->init()) {
        throw;
    }
    fTCResiduals=new TClonesArray("TpcResidual");
    unsigned int nEvents = fInputChain->GetEntries();
    std::cout << "Number of events in the input file: " << nEvents << std::endl;
    if(fSaveHistos){
        initHistos();
    }
    return false;
}

bool TpcStandaloneResCalc::execute() {
    std::cout << "TpcStandaloneResCalc::execute()" << std::endl;
    if (fInputChain==NULL){
        std::cerr<<"TpcStandaloneResCalc::execute():"
                <<"Input chain missing. Aborting"<<std::endl;
        return true;
    }
    if(fFopiBranches.size() >0 and fFopiTree==0){
        std::cerr<<"TpcStandaloneResCalc::execute(): "
                <<" fopiBranches selected but no fopi chain"
               <<"Aborting"<<std::endl;
        return true;
    }
    
    if(fFopiBranches.size()>0){
        TString major = fFopiEventBranchName;
        TString minor = major;
        fFopiTree->SetBranchStatus(fFopiEventBranchName+".*",1);
        fInputChain->SetBranchStatus(fTpcEventIdentifierBranchName+".*",1);
        std::cout<<fFopiTree->GetEntries()<<std::endl;
        //please let us never change the accessor names
        major += ".GetSpillNum()";
        minor += ".GetEvtNum()";
        std::cout << "major: "<<major<<", minor: "<<minor<<std::endl;
        int nEv=fFopiTree->BuildIndex(major,minor);
        std::cout<<"   ... Index has "<<nEv<<" entries"<<std::endl;
        fFopiIndex=(TTreeIndex*) fFopiTree->GetTreeIndex();
        
    }
    int nEvents = fInputChain->GetEntries();
    
    if (fMaxEvents != 0 && nEvents > fMaxEvents) {
        nEvents = fMaxEvents;
    }
    
    TFile* outfile;
    TTree* outTree;
    
    outfile = new TFile(fOutFileName, "recreate");
    outfile->cd();
    outTree = new TTree("res", "res");
    outTree->Branch("residual", "TClonesArray", &fTCResiduals, 8000,99);
    
    // clock.Start();
    int nFailedEvents=0;
    std::cout << "TpcStandaloneResCalc::execute() nEvents:" << nEvents
              << std::endl;
    for (unsigned int iEv = 0; iEv < nEvents; ++iEv) {
    if (iEv % 1000 == 0) {
            std::cout << " " << iEv << " " <<std::flush;
        }
        if (iEv % 10000 == 0) {
            std::cout << "###########################################################"
                         "######### " << iEv << " / " << nEvents << " ("
                      << ((double)iEv / (double)nEvents) * 100 << " %)"
                      << " ##########################################################"
                         "###########" << std::endl;
        }
        
        std::map<TString, TClonesArray*>::iterator it1;
        for (it1 = fBranches.begin(); it1 != fBranches.end(); ++it1) {
            //      
            if(it1->first.Contains("track") or it1->first.Contains("Track")){
                   (it1->second)->Delete();    
            }
            
        }
        fTCResiduals->Clear();
        //fInputChain->Print();
        fInputChain->GetEvent(iEv);
        
        if(fFopiBranches.size()>0){
            TpcEventIdentifier* tpcid=(TpcEventIdentifier*) 
                    fBranches[fTpcEventIdentifierBranchName]->At(0);
            unsigned int tpcEvent = tpcid->getEventInSpill();
            unsigned int tpcSpill = tpcid->getSpill();
            //fFopiIndex->Print();
            int treeEv = fFopiIndex->GetEntryNumberWithIndex((int)tpcSpill,
                                                             (int)tpcEvent);
            if (treeEv< 0 ){
                std::cout<<"error looking for entry with"<<"\n"
                        <<"   Spill: "<<tpcSpill<<"\n"
                       <<"   Event: "<<tpcEvent<<"\n"
                      <<"# "<<++nFailedEvents<<"/ "<<iEv<<" "<<(double) nFailedEvents/(double) iEv<<"\n";
                //                Warning("FopiDataReaderTask::Exec()",
                //                        "No matching Entry in Fopi Tree could be found");
                continue;
            }
            fFopiTree->GetEvent(treeEv);
        }
        fResidualCalculator->calc();
        
        std::vector<TpcRefResidualCollection*>* resCols =
                fResidualCalculator->getResiduals();
        int nRes = resCols->size();
        //  std::cout<<nRes<<std::endl;
        for (int ires = 0; ires < nRes; ++ires) {
            TpcRefResidualCollection* resCol = (*resCols)[ires];
            for (int ii = 0; ii < resCol->getSize(); ii++) {
                const TpcResidual* miaou;
                
                miaou = resCol->getResidual(ii);
                int charge = miaou->getCharge();
                if (fChargeFiltering == 0 or charge * fChargeFiltering > 0) {
                    //	  std::cout<<ires<<", "<<ii<<", "<<fTCResiduals->GetEntriesFast()<<std::endl;
                    TpcResidual* miaou2=new ((*fTCResiduals)[fTCResiduals->GetEntriesFast()]) TpcResidual(*miaou);
                    if(fSaveHistos){
                        fillHistos(miaou, miaou->getRefCylPos(),miaou->getRefCylMom());
                    }
                }
                
                
            }  // loop over residuals in this collection
        }    // loop over residual collections for this event
        //std::cout << fTCResiduals->GetEntries() << std::endl;
        outfile->cd();
        
        outTree->Fill();
        
    }  // Event loop end
    
    outTree->Write();
    outfile->Close();
    
    return false;
}
double TpcStandaloneResCalc::getzEnd() const
{
    return fzEnd;
}

void TpcStandaloneResCalc::setzEnd(double value)
{
    fzEnd = value;
}
bool TpcStandaloneResCalc::getSaveHistos() const
{
    return fSaveHistos;
}

void TpcStandaloneResCalc::setSaveHistos(bool value)
{
    fSaveHistos = value;
}

TString TpcStandaloneResCalc::getHistoFileName() const
{
    return fHistoFileName;
}

void TpcStandaloneResCalc::setHistoFileName(const TString &value)
{
    fHistoFileName = value;
}

double TpcStandaloneResCalc::getzStart() const
{
    return fzStart;
}
void TpcStandaloneResCalc::setzStart(double value)
{
    fzStart = value;
}
double TpcStandaloneResCalc::getMaxRangeZ() const
{
    return fMaxRangeZ;
}
void TpcStandaloneResCalc::setMaxRangeZ(double value)
{
    fMaxRangeZ = value;
}
double TpcStandaloneResCalc::getMaxRangeXY() const
{
    return fMaxRangeXY;
}
void TpcStandaloneResCalc::setMaxRangeXY(double value)
{
    fMaxRangeXY = value;
}
int TpcStandaloneResCalc::getzSlicesHists() const
{
    return fzSlicesHists;
}
void TpcStandaloneResCalc::setzSlicesHists(int value)
{
    fzSlicesHists = value;
}


void TpcStandaloneResCalc::fillHistos(const TpcResidual *theRes, 
                                      const TVector3 pos20, 
                                      const TVector3 mom20) {
    TVector3 dir = theRes->getTrackDir();
    TVector3 ref = theRes->getRefPos();
    
    TVector3 pos = theRes->getHitPos();
    
    TVector3 res = theRes->getResXYZ();
    
    TVector3 dirLocal = theRes->getLocalTrackDir();
    TVector3 uvw = fAlMan->masterToLocal("tpc", theRes->getHitPos());
    //TVector3 uvw = alMan->masterToLocal("tpc", theRes->getHitPos());
    
    //double refPhi=ref.Phi()* TMath::RadToDeg();
    double phiPos = uvw.Phi() * TMath::RadToDeg();
    double phiPosInSector = getPhiInSector(phiPos);
    double phiDir = dir.Phi() * TMath::RadToDeg();
    double phiLocalDir = dirLocal.Phi() * TMath::RadToDeg();
    
    double phiOuter=pos20.Phi()*TMath::RadToDeg();
    
    double charge = theRes->getCharge();
    
    double resMag = res.Mag();
    double resPerp = res.Perp();
    double posPerp = uvw.Perp();
    
    double dPhiHist = pos.Phi() - ref.Phi();
    dPhiHist = dPhiHist * 180 / TMath::Pi();
    if (dPhiHist > 180) {
        dPhiHist -= 360;
    } else if (dPhiHist < -180) {
        dPhiHist += 360;
    }
    
    double dPhiTrackDir = phiLocalDir - phiDir;
    if (dPhiTrackDir > 180) {
        dPhiTrackDir -= 360;
    } else if (dPhiTrackDir < -180) {
        dPhiTrackDir += 360;
    }
    double perpDir = dPhiHist / fabs(dPhiHist);
    hists["hresx"]->Fill(res.X());
    hists["hresXvX"]->Fill(pos.X(), res.X());
    hists["hresXvY"]->Fill(pos.Y(), res.X());
    hists["hresXvZ"]->Fill(pos.Z(), res.X());
    hists["hresXvPerp"]->Fill(posPerp, res.X());
    hists["hresXvPhi"]->Fill(phiDir, res.X());
    hists["hresXvPhiPos"]->Fill(phiOuter, res.X());
    
    hists["hresy"]->Fill(res.Y());
    hists["hresYvX"]->Fill(pos.X(), res.Y());
    hists["hresYvY"]->Fill(pos.Y(), res.Y());
    hists["hresYvZ"]->Fill(pos.Z(), res.Y());
    hists["hresYvPerp"]->Fill(posPerp, res.Y());
    hists["hresYvPhi"]->Fill(phiPos, res.Y());
    hists["hresYvPhiPos"]->Fill(phiOuter, res.Y());
    
    hists["hresz"]->Fill(res.Z());
    hists["hresZvX"]->Fill(pos.X(), res.Z());
    hists["hresZvY"]->Fill(pos.Y(), res.Z());
    hists["hresZvZ"]->Fill(pos.Z(), res.Z());
    hists["hresZvPerp"]->Fill(posPerp, res.Z());
    hists["hresZvPhi"]->Fill(phiDir, res.Z());
    hists["hresZvPhiPos"]->Fill(phiOuter, res.Z());
    hists["hresZvPhiCdc"]->Fill(phiOuter, res.Z());
    if (charge > 0) {
        hists["hresZvPhiCdcPlus"]->Fill(phiOuter, res.Z());
    }else{
        hists["hresZvPhiCdcMinus"]->Fill(phiOuter, res.Z());
    }
    hists["hresPerp"]->Fill(resPerp);
    hists["hresPerpvX"]->Fill(pos.X(), resPerp);
    hists["hresPerpvY"]->Fill(pos.Y(), resPerp);
    hists["hresPerpvZ"]->Fill(pos.Z(), resPerp);
    hists["hresPerpvPerp"]->Fill(posPerp, resPerp);
    hists["hresPerpvPhi"]->Fill(phiPos, resPerp);
    hists["hresPerpvPhiPos"]->Fill(phiOuter, resPerp);
    
    
    hists["hresSignedPerp"]->Fill(resPerp * perpDir);
    hists["hresSignedPerpvX"]->Fill(pos.X(), resPerp * perpDir);
    hists["hresSignedPerpvY"]->Fill(pos.Y(), resPerp * perpDir);
    hists["hresSignedPerpvZ"]->Fill(pos.Z(), resPerp * perpDir);
    hists["hresSignedPerpvPerp"]->Fill(posPerp, resPerp * perpDir);
    hists["hresSignedPerpvPhi"]->Fill(phiPos,
                                      resPerp * perpDir);
    hists["hresSignedPerpvPhiPos"]->Fill(phiPos,
                                         resPerp * perpDir);
    hists["hresSignedPerpvPhiCdc"]->Fill(phiOuter,
                                         resPerp * perpDir);
    hists["hresSignedPerpvDafWeight"]->Fill(theRes->getTrackLength(),
                                            resPerp * perpDir);
    if (charge > 0) {
        hists["hresSignedPerpPlus"]->Fill(resPerp * perpDir);
        hists["hresSignedPerpvPhiPosPlus"]->Fill(phiPos,
                                                 resPerp * perpDir);
        hists["hresSignedPerpvPhiCdcPlus"]->Fill(phiOuter,
                                                 resPerp * perpDir);
        hists["hresSignedPerpvPerpPlus"]->Fill(posPerp, resPerp * perpDir);
    } else {
        hists["hresSignedPerpMinus"]->Fill(resPerp * perpDir);
        hists["hresSignedPerpvPhiPosMinus"]->Fill(phiPos,
                                                  resPerp * perpDir);
        hists["hresSignedPerpvPhiCdcMinus"]->Fill(phiOuter,
                                                  resPerp * perpDir);
        hists["hresSignedPerpvPerpMinus"]->Fill(posPerp, resPerp * perpDir);
    }
    
    hists["hresMag"]->Fill(resMag);
    hists["hresMagvX"]->Fill(pos.X(), resMag);
    hists["hresMagvY"]->Fill(pos.Y(), resMag);
    hists["hresMagvZ"]->Fill(pos.Z(), resMag);
    hists["hresMagvPerp"]->Fill(posPerp, resMag);
    hists["hresMagvPhi"]->Fill(phiPos, resMag);
    hists["hresMagvPhiPos"]->Fill(phiOuter, resMag);
    
    hists["hresPhi"]->Fill(dPhiHist);
    hists["hresPhivX"]->Fill(pos.X(), dPhiHist);
    hists["hresPhivY"]->Fill(pos.Y(), dPhiHist);
    hists["hresPhivZ"]->Fill(pos.Z(), dPhiHist);
    hists["hresPhivPerp"]->Fill(posPerp, dPhiHist);
    hists["hresPhivPhi"]->Fill(phiPos, dPhiHist);
    hists["hresPhivPhiPos"]->Fill(phiOuter, dPhiHist);
    if (charge > 0) {
        hists["hresPhiPlus"]->Fill(dPhiHist);
        hists["hresPhivPhiPosPlus"]->Fill(phiOuter, dPhiHist);
        hists["hresPhivPerpPlus"]->Fill(posPerp, dPhiHist);
    } else {
        hists["hresPhiMinus"]->Fill(dPhiHist);
        hists["hresPhivPhiPosMinus"]->Fill(phiOuter, dPhiHist);
        hists["hresPhivPerpMinus"]->Fill(posPerp, dPhiHist);
    }
    
    hists["hres_trDist"]->Fill(theRes->getTrackLength());
    hists["hresPerpvDist"]->Fill(theRes->getTrackLength(), resPerp);
    hists["hitsX"]->Fill(pos.X());
    hists["hitsY"]->Fill(pos.Y());
    hists["hitsZ"]->Fill(pos.Z());
    hists["hitsPhi"]->Fill(phiPos);
    
    hists["hitsW"]->Fill(uvw.Z());
    hists["hitsRl"]->Fill(uvw.Perp());
    hists["hitsPerpl"]->Fill(uvw.x(), uvw.y());
    hists["hitsPerp"]->Fill(pos.X(), pos.Y());
    hists["hitsZR"]->Fill(pos.Z(), posPerp);
    
    hists["hresPhivPhiPosModulo"]->Fill(phiPosInSector, dPhiHist);
    if (charge > 0) {
        hists["hresPhivPhiPosPlusModulo"]->Fill(phiPosInSector, dPhiHist);
        hists["hresSignedPerpvPhiPosPlusModulo"]->Fill(phiPosInSector,
                                                       resPerp * perpDir);
    } else {
        hists["hresPhivPhiPosMinusModulo"]->Fill(phiPosInSector, dPhiHist);
        hists["hresSignedPerpvPhiPosMinusModulo"]->Fill(phiPosInSector,
                                                        resPerp * perpDir);
    }
    hists["hresSignedPerpvPhiPosModulo"]->Fill(phiPosInSector, resPerp * perpDir);
    
    double zRange = fzEnd - fzStart;
    double zBin = zRange / ((double)fzSlicesHists);
    double zBinStart = fzStart;
    double zBinEnd = zBinStart;
    for (int i = 0; i < fzSlicesHists; ++i) {
        zBinEnd += zBin;
        if (pos.Z() < zBinEnd) {
            TString histPostFix = TString::Format("_zSlice_%i", i);
            hists["hresPhivPerp" + histPostFix]->Fill(posPerp, dPhiHist);
            hists["hresPhivPhiTemp" + histPostFix]->Fill(
                        phiOuter, dPhiHist);
            hists["hresPhivPhiPos" + histPostFix]->Fill(phiOuter,
                                                        dPhiHist);
            hists["hresSignedPerpvPhiPos" + histPostFix]->Fill(
                        phiOuter , resPerp * perpDir);
            
            if (charge > 0) {
                hists["hresPhivPhiPosPlus" + histPostFix]->Fill(
                            phiOuter, dPhiHist);
                hists["hresSignedPerpvPhiPosPlus" + histPostFix]->Fill(
                            phiOuter, resPerp * perpDir);
            } else {
                hists["hresPhivPhiPosMinus" + histPostFix]->Fill(
                            phiOuter, dPhiHist);
                hists["hresSignedPerpvPhiPosMinus" + histPostFix]->Fill(
                            phiOuter , resPerp * perpDir);
            }
            break;
        }
    }
    double phiRange = 22.5;
    int phiSlices = 20;
    double phiBinStart = -11.25;
    double phiBin = phiRange / phiSlices;
    double phiBinEnd = phiBinStart;
    
    for (int i = 0; i < phiSlices; ++i) {
        phiBinEnd += phiBin;
        
        if (phiPosInSector < phiBinEnd) {
            TString histPostFix = TString::Format("_phiSlice_%i", i);
            hists["hresSignedPerpvPerp" + histPostFix]->Fill(posPerp,
                                                             resPerp * perpDir);
            if (charge > 0) {
                hists["hresSignedPerpvPerpPlus" + histPostFix]->Fill(posPerp,
                                                                     resPerp * perpDir);
            } else {
                hists["hresSignedPerpvPerpMinus" + histPostFix]->Fill(
                            posPerp, resPerp * perpDir);
            }
            break;
        }
        phiBinStart = phiBinEnd;
    }
    
    TVectorD normResD(3);
    TVectorD normRes2D(3);
    TVector3 normRes = res.Unit();
    if(theRes->getHitCov()[0][0]!=0){
        TMatrixDSym covHitSym(3, theRes->getHitCov().GetMatrixArray());
        TMatrixDSym covHitSymInv = covHitSym.Invert();
        TMatrixDSym covTrackSym(5, theRes->getTrackCov().GetMatrixArray());
        
        normResD[0] = normRes.X();
        normResD[1] = normRes.Y();
        normResD[2] = normRes.Z();
        normRes2D[0] = res.XYvector().Unit().X();
        normRes2D[1] = res.XYvector().Unit().Y();
        normRes2D[2] = 0;
        TVectorD xD(3);
        xD[0] = 1;    xD[1] = 0;    xD[2] = 0;
        TVectorD yD(3);
        yD[0] = 0;    yD[1] = 1;    yD[2] = 0;
        TVectorD zD(3);
        zD[0] = 0;    zD[1] = 0;    zD[2] = 1;
        
        // lambda2= res*cov*res // projection of cov
        //##########################################################
        // # Additional separate histograms for chi2 method:
        //##########################################################matrix onto
        //residual vector
        // 1/lambda1=resT * inverse(cov) * residual //cut of cov-matrix onto residual
        // vector
        
        double covRes = covHitSym.Similarity(normResD);
        double covX = covHitSym.Similarity(xD);
        double covY = covHitSym.Similarity(yD);
        double covZ = covHitSym.Similarity(zD);
        double cov2D = covHitSym.Similarity(normRes2D);
        // Trackerr
        //  double covTrackRes = covTrackSym.Similarity(normResD);
        //  double covTrackX = covTrackSym.Similarity(xD);
        //  double covTrackY = covTrackSym.Similarity(yD);
        //  double covTrackZ = covTrackSym.Similarity(zD);
        //  double covTrack2D = covTrackSym.Similarity(normRes2D);
        
        double covTrackRes = 0;
        double covTrackX = 0;
        double covTrackY = 0;
        double covTrackZ = sqrt(covTrackSym[3][3]);
        ;
        double covTrack2D = sqrt(covTrackSym[4][4]);
        
        //    double covRes = 1/covHitSymInv.Similarity(normResD);
        //    double covX =  1/covHitSymInv.Similarity(xD);
        //    double covY =  1/covHitSymInv.Similarity(yD);
        //    double covZ =  1/covHitSymInv.Similarity(zD);
        //    double cov2D = 1/covHitSymInv.Similarity(normRes2D);
        //    // Trackerr
        //    double covTrackRes = 1/covTrackSymInv.Similarity(normResD);
        //    double covTrackX = 1/covTrackSymInv.Similarity(xD);
        //    double covTrackY = 1/covTrackSymInv.Similarity(yD);
        //    double covTrackZ = 1/covTrackSymInv.Similarity(zD);
        //    double covTrack2D = 1/covTrackSymInv.Similarity(normRes2D);
        
        // std::cout<<covTrackRes<<", "<<covTrack2D<<"\n";
        
        hists["hresCovX"]->Fill(covX);
        hists["hresCovY"]->Fill(covY);
        hists["hresCovZ"]->Fill(covZ);
        hists["hresCov"]->Fill(covRes);
        hists["hresCovXY"]->Fill(cov2D);
        hists["hresMagdCov"]->Fill(resMag / covRes);
        hists["hresSignedPerpdCov"]->Fill(resPerp * perpDir / cov2D);
        hists["hresMagVCov"]->Fill(covRes, resMag);
        hists["hresXVCov"]->Fill(covX, res.x());
        hists["hresYVCov"]->Fill(covY, res.y());
        hists["hresZVCov"]->Fill(covZ, res.z());
        hists["hresSignedPerpVCov"]->Fill(cov2D, resPerp * perpDir);
        hists["hresCovVPerp"]->Fill(posPerp, covRes);
        hists["hresCovPerpVPerp"]->Fill(posPerp, cov2D);
        
        hists["hresCovTrackX"]->Fill(covTrackX);
        hists["hresCovTrackY"]->Fill(covTrackY);
        hists["hresCovTrackZ"]->Fill(covTrackZ);
        hists["hresCovTrack"]->Fill(covTrackRes);
        hists["hresCovTrackXY"]->Fill(covTrack2D);
        hists["hresMagdCovTrack"]->Fill(resMag / covTrackRes);
        hists["hresSignedPerpdCovTrack"]->Fill(resPerp * perpDir / covTrack2D);
        hists["hresMagVCovTrack"]->Fill(covTrackRes, resMag);
        hists["hresXVCovTrack"]->Fill(covTrackX, res.x());
        hists["hresYVCovTrack"]->Fill(covTrackY, res.y());
        hists["hresZVCovTrack"]->Fill(covTrackZ, res.z());
        hists["hresSignedPerpVCovTrack"]->Fill(covTrack2D, resPerp * perpDir);
        hists["hresCovTrackVPerp"]->Fill(posPerp, covTrackRes);
        hists["hresCovTrackPerpVPerp"]->Fill(posPerp, covTrack2D);
    }
    hists["hresdTrackPhi"]->Fill(dPhiTrackDir);
    hists["hresdTrackPhiVPhiPos"]->Fill(phiOuter, dPhiTrackDir);
    hists["hresdTrackPhiVPerp"]->Fill(posPerp, dPhiTrackDir);
    if (charge > 0) {
        hists["hresdTrackPhiPlus"]->Fill(dPhiTrackDir);
        hists["hresdTrackPhiVPhiPosPlus"]->Fill(phiOuter, dPhiTrackDir);
        hists["hresdTrackPhiVPerpPlus"]->Fill(posPerp, dPhiTrackDir);
    } else {
        hists["hresdTrackPhiMinus"]->Fill(dPhiTrackDir);
        hists["hresdTrackPhiVPhiPosMinus"]->Fill(phiOuter, dPhiTrackDir);
        hists["hresdTrackPhiVPerpMinus"]->Fill(posPerp, dPhiTrackDir);
    }
}
void TpcStandaloneResCalc::initHistos() {
    fHistoFile=new TFile(fHistoFileName,"recreate");
    fHistoFile->mkdir("1d-histograms");
    fHistoFile->cd("1d-histograms");
    TH1D *hresx =
            new TH1D("hresX", "residual X; dx (cm)", 8000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresx"] = hresx;
    histDirName["hresx"] = "1d-histograms";
    TH1D *hresy =
            new TH1D("hresY", "residual Y; dy (cm)", 8000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresy"] = hresy;
    histDirName["hresy"] = "1d-histograms";
    TH1D *hresz =
            new TH1D("hresZ", "residual Z; dz (cm)", 8000, -fMaxRangeZ, fMaxRangeZ);
    hists["hresz"] = hresz;
    histDirName["hresz"] = "1d-histograms";
    TH1D *hresPerp = new TH1D("hresPerp", "Resolution on Perp plane;; dperp (cm)",
                              8000, 0, fMaxRangeXY);
    hists["hresPerp"] = hresPerp;
    histDirName["hresPerp"] = "1d-histograms";
    TH1D *hresMag = new TH1D("hresMag", "residual; dMag (cm)", 8000, 0, 10);
    hists["hresMag"] = hresMag;
    histDirName["hresMag"] = "1d-histograms";
    TH1D *hresPhi =
            new TH1D("hresPhi", "#phi residual; dPhi (deg)", 8000, -10, 10);
    hists["hresPhi"] = hresPhi;
    histDirName["hresPhi"] = "1d-histograms";
    TH1D *hresPhiPlus =
            new TH1D("hresPhiPlus", "#phi residual, plus; dPhi (deg)", 8000, -10, 10);
    hists["hresPhiPlus"] = hresPhiPlus;
    histDirName["hresPhiPlus"] = "1d-histograms";
    TH1D *hresPhiMinus = new TH1D(
                "hresPhiMinus", "#phi residual, minus; dPhi (deg)", 8000, -10, 10);
    hists["hresPhiMinus"] = hresPhiMinus;
    histDirName["hresPhiMinus"] = "1d-histograms";
    TH1D *hres_trDist =
            new TH1D("hres_trDist", "extrapolation length; cm", 8000, -100, 100);
    hists["hres_trDist"] = hres_trDist;
    histDirName["hres_trDist"] = "1d-histograms";
    TH1D *hresSignedPerp =
            new TH1D("hresSignedPerp", "Resolution on Perp plane; dperp (cm)", 8000,
                     -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerp"] = hresSignedPerp;
    histDirName["hresSignedPerp"] = "1d-histograms";
    TH1D *hresSignedPerpPlus = new TH1D(
                "hresSignedPerpPlus", "Resolution on Perp plane, plus; dperp (cm)", 8000,
                -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpPlus"] = hresSignedPerpPlus;
    histDirName["hresSignedPerpPlus"] = "1d-histograms";
    TH1D *hresSignedPerpMinus = new TH1D(
                "hresSignedPerpMinus", "Resolution on Perp plane, minus; dperp (cm)",
                8000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpMinus"] = hresSignedPerpMinus;
    histDirName["hresSignedPerpMinus"] = "1d-histograms";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("x_residual");
    fHistoFile->cd("x_residual");
    // X-residuals versus X,Y,Z, perp
    TH2D *hresXvX = new TH2D("hresXvX", "residual X v X;x (cm);dx (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvX"] = hresXvX;
    histDirName["hresXvX"] = "x_residual";
    TH2D *hresXvY = new TH2D("hresXvY", "residual X v Y;y (cm);dx (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvY"] = hresXvY;
    histDirName["hresXvY"] = "x_residual";
    TH2D *hresXvZ = new TH2D("hresXvZ", "residual X v Z;z (cm);dx (cm)", 500, -40,
                             40, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvZ"] = hresXvZ;
    histDirName["hresXvZ"] = "x_residual";
    TH2D *hresXvPerp =
            new TH2D("hresXvPerp", "residual X v Perp;perp (cm);dx(cm)", 500, -0,
                     fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvPerp"] = hresXvPerp;
    histDirName["hresXvPerp"] = "x_residual";
    TH2D *hresXvPhi = new TH2D(
                "hresXvPhi", "residual X v #phi_{track};#phi_{track} (deg);dx(cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvPhi"] = hresXvPhi;
    histDirName["hresXvPhi"] = "x_residual";
    TH2D *hresXvPhiPos = new TH2D(
                "hresXvPhiPos", "residual X v #phi_{pos};#phi_{pos} (deg);dx(cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresXvPhiPos"] = hresXvPhiPos;
    histDirName["hresXvPhiPos"] = "x_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("y_residual");
    fHistoFile->cd("y_residual");
    // Y-residuals versus X, Y, Z, perp
    TH2D *hresYvX = new TH2D("hresYvX", "residual Y v X;x (cm);dy (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvX"] = hresYvX;
    histDirName["hresYvX"] = "y_residual";
    TH2D *hresYvY = new TH2D("hresYvY", "residual Y v Y;y (cm);dy (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvY"] = hresYvY;
    histDirName["hresYvY"] = "y_residual";
    TH2D *hresYvZ = new TH2D("hresYvZ", "residual Y v Z;z (cm);dy (cm)", 500, -40,
                             40, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvZ"] = hresYvZ;
    histDirName["hresYvZ"] = "y_residual";
    TH2D *hresYvPerp =
            new TH2D("hresYvPerp", "residual Y v Perp;perp (cm);dy (cm)", 500, 0,
                     fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvPerp"] = hresYvPerp;
    histDirName["hresYvPerp"] = "y_residual";
    TH2D *hresYvPhi = new TH2D(
                "hresYvPhi", "residual Y v #phi_{track};#phi_{track} (deg);dy (cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvPhi"] = hresYvPhi;
    histDirName["hresYvPhi"] = "y_residual";
    TH2D *hresYvPhiPos = new TH2D(
                "hresYvPhiPos", "residual Y v #phi_{pos};#phi_{pos} (deg);dy (cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresYvPhiPos"] = hresYvPhiPos;
    histDirName["hresYvPhiPos"] = "y_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("z_residual");
    fHistoFile->cd("z_residual");
    // Z-residual versus X, Y, Z, perp
    TH2D *hresZvX = new TH2D("hresZvX", "residual Z v X;x (cm);dz (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvX"] = hresZvX;
    histDirName["hresZvX"] = "z_residual";
    TH2D *hresZvY = new TH2D("hresZvY", "residual Z v Y;y (cm);dz (cm)", 500,
                             -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvY"] = hresZvY;
    histDirName["hresZvY"] = "z_residual";
    TH2D *hresZvZ = new TH2D("hresZvZ", "residual Z v Z;z(cm);dz (cm)", 500, -40,
                             40, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvZ"] = hresZvZ;
    histDirName["hresZvZ"] = "z_residual";
    TH2D *hresZvPerp =
            new TH2D("hresZvPerp", "residual Z v Perp; perp(cm);dz (cm)", 500, 0,
                     fMaxRangeXY, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPerp"] = hresZvPerp;
    histDirName["hresZvPerp"] = "z_residual";
    TH2D *hresZvPhi = new TH2D(
                "hresZvPhi", "residual Z v #phi_{track}; #phi_{track} (deg);dz (cm)", 500,
                -180, 180, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPhi"] = hresZvPhi;
    histDirName["hresZvPhi"] = "z_residual";
    TH2D *hresZvPhiPos = new TH2D(
                "hresZvPhiPos", "residual Z v #phi_{pos}; #phi_{pos} (deg);dz (cm)", 500,
                -180, 180, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPhiPos"] = hresZvPhiPos;
    histDirName["hresZvPhiPos"] = "z_residual";
    TH2D *hresZvPhiCdc = new TH2D(
                "hresZvPhiCdc", "residual Z v #phi_{cdc}; #phi_{cdc} (deg);dz (cm)", 500,
                -180, 180, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPhiCdc"] = hresZvPhiCdc;
    histDirName["hresZvPhiCdc"] = "z_residual";
    TH2D *hresZvPhiCdcMinus = new TH2D(
                "hresZvPhiCdcMinus", "residual Z v #phi_{cdc} minus; #phi_{cdc} (deg);dz (cm)", 500,
                -180, 180, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPhiCdcMinus"] = hresZvPhiCdcMinus;
    histDirName["hresZvPhiCdcMinus"] = "z_residual";
    TH2D *hresZvPhiCdcPlus = new TH2D(
                "hresZvPhiCdcPlus", "residual Z v #phi_{cdc} Plus; #phi_{cdc} (deg);dz (cm)", 500,
                -180, 180, 500, -fMaxRangeZ, fMaxRangeZ);
    hists["hresZvPhiCdcPlus"] = hresZvPhiCdcPlus;
    histDirName["hresZvPhiCdcPlus"] = "z_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("perp_residual");
    fHistoFile->cd("perp_residual");
    // Perp-residual versus X, Y, Z, perp
    TH2D *hresPerpvX =
            new TH2D("hresPerpvX", "residual Perp v X;x (cm);dPerp (cm)", 500,
                     -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvX"] = hresPerpvX;
    histDirName["hresPerpvX"] = "perp_residual";
    TH2D *hresPerpvY =
            new TH2D("hresPerpvY", "residual Perp v Y;y (cm);dPerp (cm)", 500,
                     -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvY"] = hresPerpvY;
    histDirName["hresPerpvY"] = "perp_residual";
    TH2D *hresPerpvZ =
            new TH2D("hresPerpvZ", "residual Perp v Z;z (cm);dPerp (cm)", 500, -40,
                     40, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvZ"] = hresPerpvZ;
    histDirName["hresPerpvZ"] = "perp_residual";
    TH2D *hresPerpvPerp =
            new TH2D("hresPerpvPerp", "residual PerpvPerp;perp (cm);dPerp (cm)", 500,
                     0, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvPerp"] = hresPerpvPerp;
    histDirName["hresPerpvPerp"] = "perp_residual";
    TH2D *hresPerpvPhi =
            new TH2D("hresPerpvPhi",
                     "residual Perpv #phi_{track};#phi_{track} (deg);dPerp (cm)", 500,
                     -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvPhi"] = hresPerpvPhi;
    histDirName["hresPerpvPhi"] = "perp_residual";
    TH2D *hresPerpvPhiPos =
            new TH2D("hresPerpvPhiPos",
                     "residual Perpv #phi_{pos};#phi_{pos} (deg);dPerp (cm)", 500,
                     -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPerpvPhiPos"] = hresPerpvPhiPos;
    histDirName["hresPerpvPhiPos"] = "perp_residual";
    TH2D *hresPerpvDist =
            new TH2D("hresPerpvDist", "Residual_{xy} vs extrap length;cm;dperp (cm)",
                     500, -100, 100, 500, 0, 100);
    hists["hresPerpvDist"] = hresPerpvDist;
    histDirName["hresPerpvDist"] = "perp_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("signedPerp_residual");
    fHistoFile->cd("signedPerp_residual");
    
    TH2D *hresSignedPerpvX = new TH2D(
                "hresSignedPerpvX", "residual SignedPerp v X;x (cm);dSignedPerp (cm)",
                500, -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvX"] = hresSignedPerpvX;
    histDirName["hresSignedPerpvX"] = "signedPerp_residual";
    TH2D *hresSignedPerpvY = new TH2D(
                "hresSignedPerpvY", "residual SignedPerp v Y;y (cm);dSignedPerp (cm)",
                500, -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvY"] = hresSignedPerpvY;
    histDirName["hresSignedPerpvY"] = "signedPerp_residual";
    TH2D *hresSignedPerpvZ = new TH2D(
                "hresSignedPerpvZ", "residual SignedPerp v Z;z (cm);dSignedPerp (cm)",
                500, -fMaxRangeZ, fMaxRangeZ, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvZ"] = hresSignedPerpvZ;
    histDirName["hresSignedPerpvZ"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPerp =
            new TH2D("hresSignedPerpvPerp",
                     "residual SignedPerpvPerp;Perp (cm);dSignedPerp (cm)", 500, 0,
                     15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPerp"] = hresSignedPerpvPerp;
    histDirName["hresSignedPerpvPerp"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPerpPlus =
            new TH2D("hresSignedPerpvPerpPlus",
                     "residual SignedPerpvPerp Plus;Perp (cm);dSignedPerp (cm)", 500,
                     0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPerpPlus"] = hresSignedPerpvPerpPlus;
    histDirName["hresSignedPerpvPerpPlus"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPerpMinus =
            new TH2D("hresSignedPerpvPerpMinus",
                     "residual SignedPerpvPerp Minus;Perp (cm);dSignedPerp (cm)", 500,
                     0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPerpMinus"] = hresSignedPerpvPerpMinus;
    histDirName["hresSignedPerpvPerpMinus"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPhi = new TH2D(
                "hresSignedPerpvPhi",
                "residual SignedPerpv #phi_{track};#phi_{track} (deg);dSignedPerp (cm)",
                500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhi"] = hresSignedPerpvPhi;
    histDirName["hresSignedPerpvPhi"] = "signedPerp_residual";
    
    TH2D *hresSignedPerpvPhiPos = new TH2D(
                "hresSignedPerpvPhiPos",
                "residual SignedPerpv #phi_{pos};#phi_{pos} (deg);dSignedPerp (cm)", 1000,
                -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPos"] = hresSignedPerpvPhiPos;
    histDirName["hresSignedPerpvPhiPos"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPhiPosPlus = new TH2D(
                "hresSignedPerpvPhiPosPlus",
                "residual SignedPerpv #phi_{pos} Plus;#phi_{pos} (deg);dSignedPerp (cm)",
                1000, -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPosPlus"] = hresSignedPerpvPhiPosPlus;
    histDirName["hresSignedPerpvPhiPosPlus"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPhiPosMinus = new TH2D(
                "hresSignedPerpvPhiPosMinus",
                "residual SignedPerpv #phi_{pos} Minus;#phi_{pos} (deg);dSignedPerp (cm)",
                1000, -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPosMinus"] = hresSignedPerpvPhiPosMinus;
    histDirName["hresSignedPerpvPhiPosMinus"] = "signedPerp_residual";
    
    TH2D *hresSignedPerpvPhiCdc = new TH2D(
                "hresSignedPerpvPhiCdc",
                "residual SignedPerpv #phi_{cdc};#phi_{cdc} (deg);dSignedPerp (cm)", 360,
                -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiCdc"] = hresSignedPerpvPhiCdc;
    histDirName["hresSignedPerpvPhiCdc"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPhiCdcPlus = new TH2D(
                "hresSignedPerpvPhiCdcPlus",
                "residual SignedPerpv #phi_{cdc} Plus;#phi_{cdc} (deg);dSignedPerp (cm)",
                360, -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiCdcPlus"] = hresSignedPerpvPhiCdcPlus;
    histDirName["hresSignedPerpvPhiCdcPlus"] = "signedPerp_residual";
    TH2D *hresSignedPerpvPhiCdcMinus = new TH2D(
                "hresSignedPerpvPhiCdcMinus",
                "residual SignedPerpv #phi_{cdc} Minus;#phi_{cdc} (deg);dSignedPerp (cm)",
                360, -180, 180, 1000, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiCdcMinus"] = hresSignedPerpvPhiCdcMinus;
    histDirName["hresSignedPerpvPhiCdcMinus"] = "signedPerp_residual";
    
    
    TH2D *hresSignedPerpvDafWeight =
            new TH2D("hresSignedPerpvDafWeight",
                     "residual SignedPerp v daf weight;daf weightdSignedPerp (cm)",
                     500, 0, 1, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvDafWeight"] = hresSignedPerpvDafWeight;
    histDirName["hresSignedPerpvDafWeight"] = "signedPerp_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("3D_residual");
    fHistoFile->cd("3D_residual");
    // Residual versus X, Y, Z, perp
    TH2D *hresMagvX = new TH2D("hresMagvX", "residual v X;x (cm);dMag (cm)", 500,
                               -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvX"] = hresMagvX;
    histDirName["hresMagvX"] = "3D_residual";
    TH2D *hresMagvY = new TH2D("hresMagvY", "residual v Y;y (cm);dMag (cm)", 500,
                               -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvY"] = hresMagvY;
    histDirName["hresMagvY"] = "3D_residual";
    TH2D *hresMagvZ = new TH2D("hresMagvZ", "residual v Z;z (cm);dMag (cm)", 500,
                               -40, 40, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvZ"] = hresMagvZ;
    histDirName["hresMagvZ"] = "3D_residual";
    TH2D *hresMagvPerp =
            new TH2D("hresMagvPerp", "residual vPerp;perp (cm);dMag (cm)", 500, 0,
                     fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvPerp"] = hresMagvPerp;
    histDirName["hresMagvPerp"] = "3D_residual";
    TH2D *hresMagvPhi = new TH2D(
                "hresMagvPhi", "residual v #phi_{track};#phi_{track} (cm);dMag (cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvPhi"] = hresMagvPhi;
    histDirName["hresMagvPhi"] = "3D_residual";
    TH2D *hresMagvPhiPos = new TH2D(
                "hresMagvPhiPos", "residual v #phi_{pos};#phi_{pos} (cm);dMag (cm)", 500,
                -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresMagvPhiPos"] = hresMagvPhiPos;
    histDirName["hresMagvPhiPos"] = "3D_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("phi_residual");
    fHistoFile->cd("phi_residual");
    // Phi-residual versus X, Y, Z, perp
    TH2D *hresPhivX = new TH2D("hresPhivX", "#phi residual v X;x (cm);dPhi (deg)",
                               500, -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivX"] = hresPhivX;
    histDirName["hresPhivX"] = "phi_residual";
    TH2D *hresPhivY = new TH2D("hresPhivY", "#phi residual v Y;y (cm);dPhi (deg)",
                               500, -fMaxRangeXY, fMaxRangeXY, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivY"] = hresPhivY;
    histDirName["hresPhivY"] = "phi_residual";
    TH2D *hresPhivZ = new TH2D("hresPhivZ", "#phi residual v Z;z (cm);dPhi (deg)",
                               500, -40, 40, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivZ"] = hresPhivZ;
    histDirName["hresPhivZ"] = "phi_residual";
    TH2D *hresPhivPerp =
            new TH2D("hresPhivPerp", "#phi residual vPerp;perp (cm);dPhi (deg)", 500,
                     0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPerp"] = hresPhivPerp;
    histDirName["hresPhivPerp"] = "phi_residual";
    
    TH2D *hresPhivPerpPlus = new TH2D(
                "hresPhivPerpPlus", "#phi residual vPerpPlus;perp (cm);dPhi (deg)", 500,
                0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPerpPlus"] = hresPhivPerpPlus;
    histDirName["hresPhivPerpPlus"] = "phi_residual";
    TH2D *hresPhivPerpMinus = new TH2D(
                "hresPhivPerpMinus", "#phi residual vPerpMinus;perp (cm);dPhi (deg)", 500,
                0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPerpMinus"] = hresPhivPerpMinus;
    histDirName["hresPhivPerpMinus"] = "phi_residual";
    
    TH2D *hresPhivPhi =
            new TH2D("hresPhivPhi",
                     "#phi residual v #phi_{track};#phi_{track} (cm);dPhi (deg)", 500,
                     -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhi"] = hresPhivPhi;
    histDirName["hresPhivPhi"] = "phi_residual";
    TH2D *hresPhivPhiPos = new TH2D(
                "hresPhivPhiPos", "#phi residual v #phi_{pos};#phi_{pos} (cm);dPhi (deg)",
                500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPos"] = hresPhivPhiPos;
    histDirName["hresPhivPhiPos"] = "phi_residual";
    TH2D *hresPhivPhiPosPlus =
            new TH2D("hresPhivPhiPosPlus",
                     "#phi residual v #phi_{pos} Plus;#phi_{pos} (cm);dPhi (deg)",
                     500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPosPlus"] = hresPhivPhiPosPlus;
    histDirName["hresPhivPhiPosPlus"] = "phi_residual";
    TH2D *hresPhivPhiPosMinus =
            new TH2D("hresPhivPhiPosMinus",
                     "#phi residual v #phi_{pos} Minus;#phi_{pos} (cm);dPhi (deg)",
                     500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPosMinus"] = hresPhivPhiPosMinus;
    histDirName["hresPhivPhiPosMinus"] = "phi_residual";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("occupancy");
    fHistoFile->cd("occupancy");
    TH1D *hitsX = new TH1D("hitsX", "hitsX;x (cm)", 800, -20, 20);
    hists["hitsX"] = hitsX;
    histDirName["hitsX"] = "occupancy";
    TH1D *hitsY = new TH1D("hitsY", "hitsY;y (cm)", 800, -20, 20);
    hists["hitsY"] = hitsY;
    histDirName["hitsY"] = "occupancy";
    TH1D *hitsZ = new TH1D("hitsZ", "hitsZ;z (cm)", 1600, -200, 200);
    hists["hitsZ"] = hitsZ;
    histDirName["hitsZ"] = "occupancy";
    
    TH1D *hitsW = new TH1D("hitsW", "hitsW;W (cm)", 1600, -200, 200);
    hists["hitsW"] = hitsW;
    histDirName["hitsW"] = "occupancy";
    TH1D *hitsRl =
            new TH1D("hitsRl", "hitsR_{local};R_{local} (cm)", 1600, 0, 15);
    hists["hitsRl"] = hitsRl;
    histDirName["hitsRl"] = "occupancy";
    TH2D *hitsPerpl = new TH2D("hitsPerpl", "hitsPerp_{local};u (cm);v (cm)", 800,
                               -20, 20, 800, -20, 20);
    hists["hitsPerpl"] = hitsPerpl;
    histDirName["hitsPerpl"] = "occupancy";
    
    TH1D *hitsPhi = new TH1D("hitsPhi", "hits #phi; #phi (deg)", 1600, -180, 180);
    hists["hitsPhi"] = hitsPhi;
    histDirName["hitsPhi"] = "occupancy";
    TH2D *hitsPerp = new TH2D("hitsPerp", "hitsPerp;x (cm);y (cm)", 800, -20,
                              20, 800, -20, 20);
    hists["hitsPerp"] = hitsPerp;
    histDirName["hitsPerp"] = "occupancy";
    TH2D *hitsZR = new TH2D("hitsZR", "hitsZR;z (cm);r (cm)", 1600, -200, 200,
                            400, 0, 15);
    hists["hitsZR"] = hitsZR;
    histDirName["hitsZR"] = "occupancy";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("covariance");
    fHistoFile->cd("covariance");
    TH1D *hresCovX =
            new TH1D("hresCovX", "hit covariance x projection", 1000, 0, 1.05);
    hists["hresCovX"] = hresCovX;
    histDirName["hresCovX"] = "covariance";
    TH1D *hresCovY =
            new TH1D("hresCovY", "hit covariance Y projection", 1000, 0, 1.05);
    hists["hresCovY"] = hresCovY;
    histDirName["hresCovY"] = "covariance";
    TH1D *hresCovZ =
            new TH1D("hresCovZ", "hit covariance Z projection", 1000, 0, 1.05);
    hists["hresCovZ"] = hresCovZ;
    histDirName["hresCovZ"] = "covariance";
    TH1D *hresCov =
            new TH1D("hresCov", "Covariance residual projection", 1000, 0, 1.05);
    hists["hresCov"] = hresCov;
    histDirName["hresCov"] = "covariance";
    TH1D *hresCovXY =
            new TH1D("hresCovXY", "Covariance residualXY projection", 1000, 0, 1.05);
    hists["hresCovXY"] = hresCovXY;
    histDirName["hresCovXY"] = "covariance";
    TH1D *hresMagdCov = new TH1D(
                "hresMagdCov", "res.Mag weighted with hit covariance", 5000, 0, 50000);
    hists["hresMagdCov"] = hresMagdCov;
    histDirName["hresMagdCov"] = "covariance";
    TH1D *hresSignedPerpdCov = new TH1D(
                "hresSignedPerpdCov", "signed res.Perp weighted with hit covariance",
                5000, -50000, 50000);
    hists["hresSignedPerpdCov"] = hresSignedPerpdCov;
    histDirName["hresSignedPerpdCov"] = "covariance";
    
    TH2D *hresMagVCov = new TH2D("hresMagVCov", "residual v covariance", 500, 0,
                                 1.05, 500, 0, 10);
    hists["hresMagVCov"] = hresMagVCov;
    histDirName["hresMagVCov"] = "covariance";
    TH2D *hresXVCov = new TH2D("hresXVCov", "residual X v covarianceX", 500, 0,
                               1.05, 500, -3, 3);
    hists["hresXVCov"] = hresXVCov;
    histDirName["hresXVCov"] = "covariance";
    TH2D *hresYVCov = new TH2D("hresYVCov", "residual y v covarianceY", 500, 0,
                               1.05, 500, -3, 3);
    hists["hresYVCov"] = hresYVCov;
    histDirName["hresYVCov"] = "covariance";
    TH2D *hresZVCov = new TH2D("hresZVCov", "residual z v covarianceZ", 500, 0,
                               1.05, 500, -10, 10);
    hists["hresZVCov"] = hresZVCov;
    histDirName["hresZVCov"] = "covariance";
    TH2D *hresSignedPerpVCov =
            new TH2D("hresSignedPerpVCov", "signed residual xy v covarianceXY", 500,
                     0, 1.05, 500, -10, 10);
    hists["hresSignedPerpVCov"] = hresSignedPerpVCov;
    histDirName["hresSignedPerpVCov"] = "covariance";
    
    TH2D *hresCovVPerp =
            new TH2D("hresCovVPerp", "hit cov vs hit_xy", 500, 0, 15, 500, 0, 0.4);
    hists["hresCovVPerp"] = hresCovVPerp;
    histDirName["hresCovVPerp"] = "covariance";
    TH2D *hresCovPerpVPerp = new TH2D("hresCovPerpVPerp", "hit cov_xy vs hit_xy",
                                      500, 0, 15, 500, 0, 0.4);
    hists["hresCovPerpVPerp"] = hresCovPerpVPerp;
    histDirName["hresCovPerpVPerp"] = "covariance";
    
    fHistoFile->cd("");
    fHistoFile->mkdir("covarianceTrack");
    fHistoFile->cd("covarianceTrack");
    TH1D *hresCovTrackX =
            new TH1D("hresCovTrackX", "track covariance x projection", 1000, 0, 1.);
    hists["hresCovTrackX"] = hresCovTrackX;
    histDirName["hresCovTrackX"] = "covarianceTrack";
    TH1D *hresCovTrackY =
            new TH1D("hresCovTrackY", "track covariance Y projection", 1000, 0, 1.);
    hists["hresCovTrackY"] = hresCovTrackY;
    histDirName["hresCovTrackY"] = "covarianceTrack";
    TH1D *hresCovTrackZ =
            new TH1D("hresCovTrackZ", "track covariance Z projection", 1000, 0, 20.);
    hists["hresCovTrackZ"] = hresCovTrackZ;
    histDirName["hresCovTrackZ"] = "covarianceTrack";
    TH1D *hresCovTrack = new TH1D(
                "hresCovTrack", "Track covariance residual projection", 1000, 0, 1.);
    hists["hresCovTrack"] = hresCovTrack;
    histDirName["hresCovTrack"] = "covarianceTrack";
    TH1D *hresCovTrackXY = new TH1D(
                "hresCovTrackXY", "Track Covariance residualXY projection", 1000, 0, .05);
    hists["hresCovTrackXY"] = hresCovTrackXY;
    histDirName["hresCovTrackXY"] = "covarianceTrack";
    TH1D *hresMagdCovTrack =
            new TH1D("hresMagdCovTrack", "res.Mag weighted with track covariance",
                     5000, 0, 50000);
    hists["hresMagdCovTrack"] = hresMagdCovTrack;
    histDirName["hresMagdCovTrack"] = "covarianceTrack";
    TH1D *hresSignedPerpdCovTrack = new TH1D(
                "hresSignedPerpdCovTrack",
                "signed res.Perp weighted with track covariance", 5000, -5000, 5000);
    hists["hresSignedPerpdCovTrack"] = hresSignedPerpdCovTrack;
    histDirName["hresSignedPerpdCovTrack"] = "covarianceTrack";
    
    TH2D *hresMagVCovTrack =
            new TH2D("hresMagVCovTrack", "residual v track covariance", 500, 0, 20.0,
                     500, 0, 10);
    hists["hresMagVCovTrack"] = hresMagVCovTrack;
    histDirName["hresMagVCovTrack"] = "covarianceTrack";
    TH2D *hresXVCovTrack =
            new TH2D("hresXVCovTrack", "residual X v track covarianceX", 500, 0, 0.2,
                     500, -3, 3);
    hists["hresXVCovTrack"] = hresXVCovTrack;
    histDirName["hresXVCovTrack"] = "covarianceTrack";
    TH2D *hresYVCovTrack =
            new TH2D("hresYVCovTrack", "residual y v track covarianceY", 500, 0, 0.2,
                     500, -3, 3);
    hists["hresYVCovTrack"] = hresYVCovTrack;
    histDirName["hresYVCovTrack"] = "covarianceTrack";
    TH2D *hresZVCovTrack =
            new TH2D("hresZVCovTrack", "residual z v track covarianceZ", 500, 0, 30.,
                     500, -30, 30);
    hists["hresZVCovTrack"] = hresZVCovTrack;
    histDirName["hresZVCovTrack"] = "covarianceTrack";
    TH2D *hresSignedPerpVCovTrack = new TH2D(
                "hresSignedPerpVCovTrack", "signed residual xy v track covarianceXY", 500,
                0, 0.2, 500, -3, 3);
    hists["hresSignedPerpVCovTrack"] = hresSignedPerpVCovTrack;
    histDirName["hresSignedPerpVCovTrack"] = "covarianceTrack";
    
    TH2D *hresCovTrackVPerp = new TH2D("hresCovTrackVPerp", "track cov vs hit_xy",
                                       500, 0, 15, 500, 0, 0);
    hists["hresCovTrackVPerp"] = hresCovTrackVPerp;
    histDirName["hresCovTrackVPerp"] = "covarianceTrack";
    TH2D *hresCovTrackPerpVPerp =
            new TH2D("hresCovTrackPerpVPerp", "track cov_xy vs hit_xy", 500, 0, 15,
                     500, 0, 0.05);
    hists["hresCovTrackPerpVPerp"] = hresCovTrackPerpVPerp;
    histDirName["hresCovTrackPerpVPerp"] = "covarianceTrack";
    // zSliced hists
    double zRange = fzEnd - fzStart;
    double zBin = zRange / ((double)fzSlicesHists);
    double zBinStart = fzStart;
    double zBinEnd = zBinStart;
    fHistoFile->cd("");
    fHistoFile->mkdir("hresPhivPerp_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresPhivPhi_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresPhivPhiPos_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresPhivPhiPosPlus_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresPhivPhiPosMinus_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresSignedPerpvPhiPos_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresSignedPerpvPhiPosPlus_zslice");
    fHistoFile->cd("");
    fHistoFile->mkdir("hresSignedPerpvPhiPosMinus_zslice");
    for (int i = 0; i < fzSlicesHists; ++i) {
        zBinEnd += zBin;
        TString histPostFix = TString::Format("_zSlice_%i", i);
        TString histPrettyPostFix =
                TString::Format(" zSlice %f - %f", zBinStart, zBinEnd);
        
        TH2D *hresPhivPerpTemp = new TH2D(
                    "hresPhivPerp" + histPostFix,
                    "#phi residual vPerp" + histPrettyPostFix + ";perp (cm);dPhi (deg)",
                    500, 0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresPhivPerp" + histPostFix] = hresPhivPerpTemp;
        histDirName["hresPhivPerp" + histPostFix] = "hresPhivPerp_zslice";
        
        TH2D *hresPhivPhiTemp =
                new TH2D("hresPhivPhi" + histPostFix,
                         "#phi residual v #phi_{track}" + histPrettyPostFix +
                         ";#phi_{track} (cm);dPhi (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresPhivPhiTemp" + histPostFix] = hresPhivPhiTemp;
        histDirName["hresPhivPhiTemp" + histPostFix] = "hresPhivPhi_zslice";
        
        TH2D *hresPhivPhiPosTemp =
                new TH2D("hresPhivPhiPos" + histPostFix,
                         "#phi residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dPhi (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresPhivPhiPos" + histPostFix] = hresPhivPhiPosTemp;
        histDirName["hresPhivPhiPos" + histPostFix] = "hresPhivPhiPos_zslice";
        
        TH2D *hresPhivPhiPosPlusTemp =
                new TH2D("hresPhivPhiPosPlus" + histPostFix,
                         "plus #phi residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dPhi (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresPhivPhiPosPlus" + histPostFix] = hresPhivPhiPosPlusTemp;
        histDirName["hresPhivPhiPosPlus" + histPostFix] =
                "hresPhivPhiPosPlus_zslice";
        
        TH2D *hresPhivPhiPosMinusTemp =
                new TH2D("hresPhivPhiPosMinus" + histPostFix,
                         "Minus #phi residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dPhi (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresPhivPhiPosMinus" + histPostFix] = hresPhivPhiPosMinusTemp;
        histDirName["hresPhivPhiPosMinus" + histPostFix] =
                "hresPhivPhiPosMinus_zslice";
        
        // PERP!
        TH2D *hresSignedPerpvPhiPosTemp =
                new TH2D("hresSignedPerpvPhiPos" + histPostFix,
                         "SignedPerp residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dSignedPerp (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresSignedPerpvPhiPos" + histPostFix] = hresSignedPerpvPhiPosTemp;
        histDirName["hresSignedPerpvPhiPos" + histPostFix] =
                "hresSignedPerpvPhiPos_zslice";
        
        TH2D *hresSignedPerpvPhiPosPlusTemp =
                new TH2D("hresSignedPerpvPhiPosPlus" + histPostFix,
                         "plus SignedPerp residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dSignedPerp (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresSignedPerpvPhiPosPlus" + histPostFix] =
                hresSignedPerpvPhiPosPlusTemp;
        histDirName["hresSignedPerpvPhiPosPlus" + histPostFix] =
                "hresSignedPerpvPhiPosPlus_zslice";
        
        TH2D *hresSignedPerpvPhiPosMinusTemp =
                new TH2D("hresSignedPerpvPhiPosMinus" + histPostFix,
                         "Minus SignedPerp residual v #phi_{pos}" + histPrettyPostFix +
                         ";#phi_{pos} (cm);dSignedPerp (deg)",
                         500, -180, 180, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresSignedPerpvPhiPosMinus" + histPostFix] =
                hresSignedPerpvPhiPosMinusTemp;
        histDirName["hresSignedPerpvPhiPosMinus" + histPostFix] =
                "hresSignedPerpvPhiPosMinus_zslice";
        
        zBinStart = zBinEnd;
    }
    
    double phiRange = 22.5;
    int phiSlices = 20;
    double phiBinStart = -11.25;
    double phiBin = phiRange / phiSlices;
    double phiBinEnd = phiBinStart;
    fHistoFile->cd("");
    fHistoFile->mkdir("hresSignedPerpVSignedPerp_phiSlice");
    fHistoFile->cd("hresSignedPerpVSignedPerp_phiSlice");
    for (int i = 0; i < phiSlices; ++i) {
        phiBinEnd += phiBin;
        TString histPostFix = TString::Format("_phiSlice_%i", i);
        TString histPrettyPostFix =
                TString::Format("#phi-slice %f - %f", phiBinStart, phiBinEnd);
        TH2D *hresSignedPerpvPerpTmp =
                new TH2D("hresSignedPerpvSignedPerp" + histPostFix,
                         "residual SignedPerpvPerp" + histPrettyPostFix +
                         ";Perp (cm);dSignedPerp (cm)",
                         500, 0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
        TH2D *hresSignedPerpvPerpPlusTmp =
                new TH2D("hresSignedPerpvSignedPerpPlus" + histPostFix,
                         "residual SignedPerpvPerp" + histPrettyPostFix +
                         ", plus;Perp (cm);dSignedPerp (cm)",
                         500, 0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
        TH2D *hresSignedPerpvPerpMinusTmp =
                new TH2D("hresSignedPerpvSignedPerpMinus" + histPostFix,
                         "residual SignedPerpvPerp" + histPrettyPostFix +
                         ", minus;Perp (cm);dSignedPerp (cm)",
                         500, 0, 15, 500, -fMaxRangeXY, fMaxRangeXY);
        hists["hresSignedPerpvPerp" + histPostFix] = hresSignedPerpvPerpTmp;
        hists["hresSignedPerpvPerpPlus" + histPostFix] = hresSignedPerpvPerpPlusTmp;
        hists["hresSignedPerpvPerpMinus" + histPostFix] =
                hresSignedPerpvPerpMinusTmp;
        histDirName["hresSignedPerpvPerp" + histPostFix] =
                "hresSignedPerpVSignedPerp_phiSlice";
        histDirName["hresSignedPerpvPerpPlus" + histPostFix] =
                "hresSignedPerpVSignedPerp_phiSlice";
        histDirName["hresSignedPerpvPerpMinus" + histPostFix] =
                "hresSignedPerpVSignedPerp_phiSlice";
        phiBinStart = phiBinEnd;
    }
    
    // Hists modulo 22.5
    fHistoFile->cd();
    fHistoFile->mkdir("modulo");
    
    TH2D *hresPhivPhiPosModulo = new TH2D(
                "hresPhivPhiPos", "#phi residual v $phi modulo 22.5;#phi(cm);d#Phi (deg)",
                500, -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPosModulo"] = hresPhivPhiPosModulo;
    histDirName["hresPhivPhiPosModulo"] = "modulo";
    TH2D *hresPhivPhiPosPlusModulo =
            new TH2D("hresPhivPhiPosPlus",
                     "#phi residual v $phi modulo 22.5 plus;#phi(cm);d#Phi (deg)",
                     500, -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPosPlusModulo"] = hresPhivPhiPosPlusModulo;
    histDirName["hresPhivPhiPosPlusModulo"] = "modulo";
    TH2D *hresPhivPhiPosMinusModulo =
            new TH2D("hresPhivPhiPosMinus",
                     "#phi residual v $phi modulo 22.5 minus;#phi(cm);d#Phi (deg)",
                     500, -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresPhivPhiPosMinusModulo"] = hresPhivPhiPosMinusModulo;
    histDirName["hresPhivPhiPosMinusModulo"] = "modulo";
    
    TH2D *hresSignedPerpvPhiPosModulo =
            new TH2D("hresSignedPerpvPhiPos",
                     "signed perp residual v $phi modulo 22.5;#phi(cm);d#Phi (deg)",
                     500, -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPosModulo"] = hresSignedPerpvPhiPosModulo;
    histDirName["hresSignedPerpvPhiPosModulo"] = "modulo";
    TH2D *hresSignedPerpvPhiPosPlusModulo = new TH2D(
                "hresSignedPerpvPhiPosPlus",
                "signed perp residual v $phi modulo 22.5 plus;#phi(cm);d#Phi (deg)", 500,
                -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPosPlusModulo"] = hresSignedPerpvPhiPosPlusModulo;
    histDirName["hresSignedPerpvPhiPosPlusModulo"] = "modulo";
    TH2D *hresSignedPerpvPhiPosMinusModulo =
            new TH2D("hresSignedPerpvPhiPosMinus",
                     "signed residual v $phi modulo 22.5 minus;#phi(cm);d#Phi (deg)",
                     500, -11.25, 11.25, 500, -fMaxRangeXY, fMaxRangeXY);
    hists["hresSignedPerpvPhiPosMinusModulo"] = hresSignedPerpvPhiPosMinusModulo;
    histDirName["hresSignedPerpvPhiPosMinusModulo"] = "modulo";
    
    fHistoFile->cd();
    fHistoFile->mkdir("TrackDirs");
    
    TH1D *hresdTrackPhi =
            new TH1D("hresdTrackPhi", "d#phi_{track};#dphi", 1000, -5, 5);
    hists["hresdTrackPhi"] = hresdTrackPhi;
    histDirName["hresdTrackPhi"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPhiPos =
            new TH2D("hresdTrackPhiVPhisPos",
                     "d#phi_{track} vs #phi_{pos} (deg);#phi_{pos} (deg);#dphi", 1000,
                     -180, 180, 1000, -5, 5);
    hists["hresdTrackPhiVPhiPos"] = hresdTrackPhiVPhiPos;
    histDirName["hresdTrackPhiVPhiPos"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPerp = new TH2D(
                "hresdTrackPhiVPerp", "d#phi_{track} vs perp;perp (cm);#dphi (deg)", 1000,
                0, 15, 1000, -5, 5);
    hists["hresdTrackPhiVPerp"] = hresdTrackPhiVPerp;
    histDirName["hresdTrackPhiVPerp"] = "TrackDirs";
    
    // plus
    TH1D *hresdTrackPhiPlus =
            new TH1D("hresdTrackPhiPlus", "d#phi_{track} Plus;#dphi", 1000, -5, 5);
    hists["hresdTrackPhiPlus"] = hresdTrackPhiPlus;
    histDirName["hresdTrackPhiPlus"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPhiPosPlus =
            new TH2D("hresdTrackPhiVPhisPosPlus",
                     "d#phi_{track} vs #phi_{pos} (deg) Plus;#phi_{pos} (deg);#dphi",
                     1000, -180, 180, 1000, -5, 5);
    hists["hresdTrackPhiVPhiPosPlus"] = hresdTrackPhiVPhiPosPlus;
    histDirName["hresdTrackPhiVPhiPosPlus"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPerpPlus =
            new TH2D("hresdTrackPhiVPerpPlus",
                     "d#phi_{track} vs perp Plus;perp (cm);#dphi (deg)", 1000, 0, 15,
                     1000, -5, 5);
    hists["hresdTrackPhiVPerpPlus"] = hresdTrackPhiVPerpPlus;
    histDirName["hresdTrackPhiVPerpPlus"] = "TrackDirs";
    // minus
    TH1D *hresdTrackPhiMinus =
            new TH1D("hresdTrackPhiMinus", "d#phi_{track} Minus;#dphi", 1000, -5, 5);
    hists["hresdTrackPhiMinus"] = hresdTrackPhiMinus;
    histDirName["hresdTrackPhiMinus"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPhiPosMinus =
            new TH2D("hresdTrackPhiVPhisPosMinus",
                     "d#phi_{track} vs #phi_{pos} (deg) Minus;#phi_{pos} (deg);#dphi",
                     1000, -180, 180, 1000, -5, 5);
    hists["hresdTrackPhiVPhiPosMinus"] = hresdTrackPhiVPhiPosMinus;
    histDirName["hresdTrackPhiVPhiPosMinus"] = "TrackDirs";
    
    TH2D *hresdTrackPhiVPerpMinus =
            new TH2D("hresdTrackPhiVPerpMinus",
                     "d#phi_{track} vs perp Minus;perp (cm);#dphi (deg)", 1000, 0, 15,
                     1000, -5, 5);
    hists["hresdTrackPhiVPerpMinus"] = hresdTrackPhiVPerpMinus;
    histDirName["hresdTrackPhiVPerpMinus"] = "TrackDirs";
    
}
void TpcStandaloneResCalc::writeHistos() {
    std::cout << "writing hists to file: " << fHistoFileName << std::endl;
    if(!fSaveHistos){
     return;   
    }
    fHistoFile->cd();
    std::map<TString, TH1 *>::iterator it;
    for (it = hists.begin(); it != hists.end(); ++it) {
        std::cout << "." << std::flush;
        fHistoFile->cd();
        fHistoFile->cd(histDirName[(*it).first]);
        (*it).second->Write();
        if (it->second->InheritsFrom("TH2D")) {
            TString name = it->second->GetName();
            TProfile *temp = ((TH2D *)it->second)->ProfileX();
            temp->Write();
        }
    }
}