//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of resCalc
//      Uses existing GFTracks and calculates
//      residua in XY.
//      Developed as part of the GEM-TPC project
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Felix Boehmer                     (original author)
//      Physik Department E18, TUM
//
//-----------------------------------------------------------

#include "TpcRefGFTrkTpcClusterResCalc_Alignment.h"
#include "TClonesArray.h"
#include "TpcCluster.h"
#include "TpcResidual.h"
#include "TpcRefResidualCollection.h"
#include "TpcSPHit.h"
#include "TpcAlignmentManager.h"

#include "GFTrack.h"
#include "GFException.h"
#include "GFAbsTrackRep.h"
#include "GFRecoHitFactory.h"
#include "GFRecoHitProducer.h"
#include "RKTrackRep.h"
#include "GFTrackCand.h"
#include <FairRootManager.h>

#include "TpcRefGFTrkResCalc_Alignment.h"
#include <TClass.h>


#include <string>
#include <iostream>


ClassImp(TpcRefGFTrkTpcClusterResCalc_Alignment)

TpcRefGFTrkTpcClusterResCalc_Alignment::TpcRefGFTrkTpcClusterResCalc_Alignment()
    : fGFTrackBranchName("GFTrack"),
      fTpcClusterBranchName("TpcCluster"),
      fTpcSPHitBranchName("TpcCluster"),
      fNRep(-1),
      fDetID(-1),
      fRequireGoodFit(true),
      fMomCutRefTrack(0.5),
      fQAFileName("TpcRefGFTrkTpcClusterResCalc_AlignmentQA.root"),
      fIdealTpcTrackTpcClusterMatching(false),
      fSimulateTrackMatching(false),
      fdPhiCut(10),
      fMatchOnCylinder(false),
      fdebug(false)

{
    fNExpectedBranches = 2;
    TH1D* rHist=new TH1D("rHist","rHist",1000,0,100);
    hists["rHist"]=rHist;
    TH1D* rHistCut=new TH1D("rHistCut","rHistCut",1000,0,100);
    hists["rHistCut"]=rHistCut;
    TH1D* momHist=new TH1D("momHist","momHist",1000,0,5);
    hists["momHist"]=momHist;
    TH1D* momHistCut=new TH1D("momHistCut","momHistCut",1000,0,5);
    hists["momHistCut"]=momHistCut;

    TH2D* histdPhiVsTpcPhi=new TH2D("dPhiVsTpcPhi","dPhiVsTpcPhi",2500,-180,180,2500,-10,10);
    hists["dPhiVsTpcPhi"]=histdPhiVsTpcPhi;
    TH2D* histdPhiVsMisPhi=new TH2D("dPhiVsMisPhi","dPhiVsMisPhi",2500,-180,180,2500,-10,10);
    hists["dPhiVsMisPhi"]=histdPhiVsMisPhi;
    TH2D* histdHitPhiVsTpcHitPhi=new TH2D("dHitPhiVsTpcHitPhi","dHitPhiVsTpcHitPhi",2500,-180,180,2500,-10,10);
    hists["dHitPhiVsTpcHitPhi"]=histdHitPhiVsTpcHitPhi;
    TH2D* histdHitPhiVsMisHitPhi=new TH2D("dHitPhiVsMisHitPhi","dHitPhiVsMisHitPhi",2500,-180,180,2500,-10,10);
    hists["dHitPhiVsMisHitPhi"]=histdHitPhiVsMisHitPhi;
    TH1D* histdPhi=new TH1D("dPhi","dPhi",2500,-180,180);
    hists["dPhi"]=histdPhi;
    TH1D* histdHitPhi=new TH1D("dHitPhi","dHitPhi",2500,-180,180);
    hists["dHitPhi"]=histdHitPhi;
    //Positive************************************************************************************
    TH2D* histPositivedPhiVsTpcPhi=new TH2D("PositivedPhiVsTpcPhi","PositivedPhiVsTpcPhi",2500,-180,180,2500,-10,10);
    hists["PositivedPhiVsTpcPhi"]=histPositivedPhiVsTpcPhi;
    TH2D* histPositivedPhiVsMisPhi=new TH2D("PositivedPhiVsMisPhi","PositivedPhiVsMisPhi",2500,-180,180,2500,-10,10);
    hists["PositivedPhiVsMisPhi"]=histPositivedPhiVsMisPhi;
    TH2D* histPositivedHitPhiVsTpcHitPhi=new TH2D("PositivedHitPhiVsTpcHitPhi","PositivedHitPhiVsTpcHitPhi",2500,-180,180,2500,-10,10);
    hists["PositivedHitPhiVsTpcHitPhi"]=histPositivedHitPhiVsTpcHitPhi;
    TH2D* histPositivedHitPhiVsMisHitPhi=new TH2D("PositivedHitPhiVsMisHitPhi","PositivedHitPhiVsMisHitPhi",2500,-180,180,2500,-10,10);
    hists["PositivedHitPhiVsMisHitPhi"]=histPositivedHitPhiVsMisHitPhi;
    TH1D* histPositivedPhi=new TH1D("PositivedPhi","PositivedPhi",2500,-180,180);
    hists["PositivedPhi"]=histPositivedPhi;
    TH1D* histPositivedHitPhi=new TH1D("PositivedHitPhi","PositivedHitPhi",2500,-180,180);
    hists["PositivedHitPhi"]=histPositivedHitPhi;
    //Negative ************************************************************************************
    TH2D* histNegativedPhiVsTpcPhi=new TH2D("NegativedPhiVsTpcPhi","NegativedPhiVsTpcPhi",2500,-180,180,2500,-10,10);
    hists["NegativedPhiVsTpcPhi"]=histNegativedPhiVsTpcPhi;
    TH2D* histNegativedPhiVsMisPhi=new TH2D("NegativedPhiVsMisPhi","NegativedPhiVsMisPhi",2500,-180,180,2500,-10,10);
    hists["NegativedPhiVsMisPhi"]=histNegativedPhiVsMisPhi;
    TH2D* histNegativedHitPhiVsTpcHitPhi=new TH2D("NegativedHitPhiVsTpcHitPhi","NegativedHitPhiVsTpcHitPhi",2500,-180,180,2500,-10,10);
    hists["NegativedHitPhiVsTpcHitPhi"]=histNegativedHitPhiVsTpcHitPhi;
    TH2D* histNegativedHitPhiVsMisHitPhi=new TH2D("NegativedHitPhiVsMisHitPhi","NegativedHitPhiVsMisHitPhi",2500,-180,180,2500,-10,10);
    hists["NegativedHitPhiVsMisHitPhi"]=histNegativedHitPhiVsMisHitPhi;
    TH1D* histNegativedPhi=new TH1D("NegativedPhi","NegattivedPhi",2500,-180,180);
    hists["NegativedPhi"]=histNegativedPhi;
    TH1D* histNegativedHitPhi=new TH1D("NegativedHitPhi","NegattivedHitPhi",2500,-180,180);
    hists["NegativedHitPhi"]=histNegativedHitPhi;

    TH1D* histMisPhi=new TH1D("MisPhi","MisPhi",2500,-180,180);
    hists["MisPhi"]=histMisPhi;
    TH1D* histTpcPhi=new TH1D("TpcPhi","TpcPhi",2500,-180,180);
    hists["TpcPhi"]=histTpcPhi;
    TH1D* histMisHitPhi=new TH1D("MisHitPhi","MisHitPhi",2500,-180,180);
    hists["MisHitPhi"]=histMisHitPhi;
    TH1D* histTpcHitPhi=new TH1D("TpcHitPhi","TpcHitPhi",2500,-180,180);
    hists["TpcHitPhi"]=histTpcHitPhi;


    TH1D* dTheta=new TH1D("dTheta","dTheta",1000,-180,180);
    hists["dTheta"]=dTheta;
    TH2D* dThetaVTheta=new TH2D("dThetaVTheta","dThetaVTheta",2500,-180,180,2500,-180,180);
    hists["dThetaVTheta"]=dThetaVTheta;
    TH2D* dThetaVRefTheta=new TH2D("dThetaVRefTheta","dThetaVRefTheta",2500,-180,180,2500,-180,180);
    hists["dThetaVRefTheta"]=dThetaVRefTheta;

    TH1D* dThetaCut=new TH1D("dThetaCut","dThetaCut",1000,-180,180);
    hists["dThetaCut"]=dThetaCut;
    TH2D* dThetaVThetaCut=new TH2D("dThetaVThetaCut","dThetaVThetaCut",2500,-180,180,2500,-180,180);
    hists["dThetaVThetaCut"]=dThetaVThetaCut;
    TH2D* dThetaVRefThetaCut=new TH2D("dThetaVRefThetaCut","dThetaVRefThetaCut",2500,-180,180,2500,-180,180);
    hists["dThetaVRefThetaCut"]=dThetaVRefThetaCut;

    TH1D* dThetaPos=new TH1D("dThetaPos","dThetaPos",1000,-180,180);
    hists["dThetaPos"]=dThetaPos;
    TH2D* dThetaPosVThetaPos=new TH2D("dThetaPosVThetaPos","dThetaPosVThetaPos",2500,-180,180,2500,-180,180);
    hists["dThetaPosVThetaPos"]=dThetaPosVThetaPos;
    TH2D* dThetaPosVRefThetaPos=new TH2D("dThetaPosVRefThetaPos","dThetaPosVRefThetaPos",2500,-180,180,2500,-180,180);
    hists["dThetaPosVRefThetaPos"]=dThetaPosVRefThetaPos;
    TH1D* dThetaPosCut=new TH1D("dThetaPosCut","dThetaPosCut",1000,-180,180);
    hists["dThetaPosCut"]=dThetaPosCut;
    TH2D* dThetaPosVThetaPosCut=new TH2D("dThetaPosVThetaPosCut","dThetaPosVThetaPosCut",2500,-180,180,2500,-180,180);
    hists["dThetaPosVThetaPosCut"]=dThetaPosVThetaPosCut;
    TH2D* dThetaPosVRefThetaPosCut=new TH2D("dThetaPosVRefThetaPosCut","dThetaPosVRefThetaPosCut",2500,-180,180,2500,-180,180);
    hists["dThetaPosVRefThetaPosCut"]=dThetaPosVRefThetaPosCut;

    zSlices=10;
    zStart=-30;
    zEnd=40;
    double zRange=zEnd-zStart;
    double zBin=zRange/((double)zSlices);
    double zBinStart=zStart;
    double zBinEnd=zBinStart;
    for(int i =0;i<zSlices;++i){
        std::cout<<"i: "<<i<<" ";
        zBinEnd+=zBin;
        TString name=TString::Format("dPhiPlusPos_vsTpcPhi_ZSlice_%i",i);
        TString name2=TString::Format("dPhiPlusMom_vsTpcPhi_ZSlice_%i",i);
        TString name3=TString::Format("dPhiMinusPos_vsTpcPhi_ZSlice_%i",i);
        TString name4=TString::Format("dPhiMinusMom_vsTpcPhi_ZSlice_%i",i);
        TString prettyName=TString::Format("dPhiPlusPos_vsTpcPhi_ZSlice %f - %f",zBinStart,zBinEnd);
        TString prettyName2=TString::Format("dPhiPlusMom_vsTpcPhi_ZSlice  %f - %f",zBinStart,zBinEnd);
        TString prettyName3=TString::Format("dPhiMinusPos_vsTpcPhi_ZSlice  %f - %f",zBinStart,zBinEnd);
        TString prettyName4=TString::Format("dPhiMinusMom_vsTpcPhi_ZSlice  %f - %f",zBinStart,zBinEnd);
        zBinStart=zBinEnd;
        TH2D* hist=new TH2D(name,prettyName,500,-180,180,500,-10,10);
        TH2D* hist2=new TH2D(name2,prettyName2,500,-180,180,500,-10,10);
        TH2D* hist3=new TH2D(name3,prettyName3,500,-180,180,500,-10,10);
        TH2D* hist4=new TH2D(name4,prettyName4,500,-180,180,500,-10,10);
        hists[name.Data()]=hist;
        hists[name2.Data()]=hist2;
        hists[name3.Data()]=hist3;
        hists[name4.Data()]=hist4;
        sliceNamesPlusPos.push_back(name.Data());
        sliceNamesPlusMom.push_back(name2.Data());
        sliceNamesMinusPos.push_back(name3.Data());
        sliceNamesMinusMom.push_back(name4.Data());
        std::cout<<name<<std::endl;
    }
}

TpcRefGFTrkTpcClusterResCalc_Alignment::~TpcRefGFTrkTpcClusterResCalc_Alignment() {
}

bool TpcRefGFTrkTpcClusterResCalc_Alignment::init() {
    if(fIdealTpcTrackTpcClusterMatching){
        fNExpectedBranches=3;
    }
    if(fBranchMap.size() != fNExpectedBranches) {
        std::cout<<"TpcRefGFTrkTpcClusterResCalc_Alignment::init(): "<<std::endl
                <<"   ERROR: Wrong number of input branches ("
               <<fBranchMap.size()<<")!"<<std::endl;
        return true;
    }
    //check for proper data types:
    bool foundGFTrack=false;
    bool foundTpcCluster=false;
    bool foundTpcSPHit=false;
    std::string GFTrackRefType("GFTrack");
    std::string tpcClusterRefType("TpcCluster");
    std::string tpcSPHitType("TpcSPHit");
    //identify the needed branches in the map
    std::map<TString, TClonesArray*>::const_iterator it;
    for(it=fBranchMap.begin(); it!=fBranchMap.end(); it++) {
        std::string type((it->second)->GetClass()->GetName());
        if(type.find(GFTrackRefType)<type.size()) {
            fGFTrackArray=it->second;
            fGFTrackBranchName = it->first;
            foundGFTrack=true;
            std::cout<<"TpcRefGFTrkTpcClusterResCalc_Alignment::init(): "<<std::endl
                    <<"   INFO: got "<<type<<" as Object type for GFTracks with branch name "<<it->first
                   <<std::endl;
        }
        if(type.find(tpcClusterRefType)<type.size()) {
            fTpcClusterArray=it->second;
            foundTpcCluster=true;
            std::cout<<"TpcRefGFTrkTpcClusterResCalc_Alignment::init(): "<<std::endl
                    <<"   INFO: got "<<type<<" as Object type for TPC Clusters with branch name "<<it->first
                   <<std::endl;
        }
        if(fIdealTpcTrackTpcClusterMatching){
            if(type.find(tpcSPHitType)<type.size()) {
                fTpcSPHitArray=it->second;
                foundTpcSPHit=true;
                std::cout<<"TpcRefGFTrkTpcClusterResCalc_Alignment::init(): "<<std::endl
                        <<"   INFO: got "<<type<<" as Object type for TPC SP Hits with branch name "<<it->first
                       <<std::endl;
            }
        }
    }
    if((!(foundGFTrack&&foundTpcCluster))&&(fIdealTpcTrackTpcClusterMatching||foundTpcSPHit)) {
        std::cout<<"TpcRefGFTrkTpcClusterResCalc_Alignment::init(): "<<std::endl
                <<"   ERROR: Wrong Class Types in input branches!"<<std::endl;
        return true;
    }

    //all ok
    return false;
}


//main functionality
int TpcRefGFTrkTpcClusterResCalc_Alignment::calc() {
    double rad2deg=180/TMath::Pi();
    //clear results
    for(unsigned int i=0; i<fResults.size(); i++)
        delete fResults[i];

    fResults.clear();


    unsigned int nTracks = fGFTrackArray->GetEntriesFast();
    // GF trackfit loop
    for(unsigned int itr=0; itr<nTracks; itr++) {
        GFTrack* itrack = (GFTrack*)((*fGFTrackArray)[itr]);
        const GFAbsTrackRep* rep;
        if(fNRep>=0){
            rep = itrack->getTrackRep(fNRep);
        }else{
            rep = itrack->getCardinalRep();
        }
        if(abs(rep->getPDG())>10000){
            continue;
        }
        GFTrackCand cand = itrack->getCand();
        if(rep->getStatusFlag()!=0 && fRequireGoodFit){
            continue;
        }


        GFDetPlane tpcLastPl = rep->getLastPlane();

        RKTrackRep* lastRep = (RKTrackRep*) rep->clone();
        const TMatrixDSym lastCov=lastRep->getLastCov();
        TVectorD lastState = lastRep->getLastState();
        TVector3 refMom=((GFAbsTrackRep*)lastRep)->getMom();
        TVector3 refPos=((GFAbsTrackRep*)lastRep)->getPos();
        hists["rHist"]->Fill(refPos.Perp());
        hists["momHist"]->Fill(refMom.Mag());
        double p_firstPoint=refMom.Mag();
        double pt_firstPoint=refMom.Perp();
        double p2=0;
        double pt2=0;
        if(refMom.Mag()<fMomCutRefTrack){
            delete lastRep;
            continue;
        }
        double len=0;
        double q_over_p=lastState[0];
        if(!fdebug){
            if((1/q_over_p)<0.05){
                delete lastRep;
                continue;
            }
        }
        if(fMatchOnCylinder){
            try{
                len=((GFAbsTrackRep*)lastRep)->extrapolate(tpcLastPl);
            }catch(GFException& ex){
                std::cout<<"TpcRefGFTrackResCalc::calc(): Error during 1st Tpc extrap last point: \n"
                        <<"     "<<ex.what()<<std::endl;
                delete lastRep;
                continue;
            }
            p2=((GFAbsTrackRep*)lastRep)->getMom().Mag();
            pt2=((GFAbsTrackRep*)lastRep)->getMom().Perp();
            lastRep->setData(lastState,tpcLastPl,&lastCov,lastRep->getAuxInfo(tpcLastPl));
        }else{
            lastRep->setData(lastState,tpcLastPl,&lastCov,lastRep->getAuxInfo(lastRep->getFirstPlane()));
        }

        refMom=((GFAbsTrackRep*)lastRep)->getMom();
        refPos=((GFAbsTrackRep*)lastRep)->getPos();
        double pt=refMom.Perp();
        double p=refMom.Mag();
        if(fMatchOnCylinder){
            try{
                lastRep->extrapolateToCylinder(15,refPos,refMom);
            }catch(GFException& ex){

                std::cout<<"TpcRefGFTrackResCalc::calc(): Error during 1st Tpc extrap to cylinder: \n"
                        <<"pt lastState: "<<pt<<" p last state: "<<p<<"\n"
                       <<"pt firstPoint: "<<pt_firstPoint<<" p1: "<<p_firstPoint<<" length "<<len<<"\n"
                      <<"pt lastRep before set data: "<<pt2<<" p2: "<<p2<<" p over q: "<<(1/q_over_p)<<"\n"
                     <<tpcLastPl.getO().X()<<", "<<tpcLastPl.getO().Y()<<", "<<tpcLastPl.getO().Z()<<"\n"
                    <<"pdgId: "<<rep->getPDG()<<"\n"
                   <<"mcId: "<<cand.getMcTrackId()<<"\n"
                  <<"     "<<ex.what()<<"\n";

                if(fdebug){
                    double ncl=cand.getNHits();
                    double lastPhi=0;
                    for(unsigned int icl=0; icl<(ncl); icl++) {
                        TMatrixT<double> cov(3,3);
                        TpcCluster* theCl;
                        TVector3 posSP;
                        if(fIdealTpcTrackTpcClusterMatching){
                            int detID,hitID;
                            cand.getHit(icl,detID,hitID);
                            TpcSPHit* hit= (TpcSPHit*) fTpcSPHitArray->At(hitID);
                            posSP=hit->getPos();
                            theCl=(TpcCluster*) fTpcClusterArray->At(hit->getClusterID());
                            if(theCl==NULL){
                                std::cout<<"########################## "
                                        <<"detID: "<<detID<<" hitID: "<<hitID<<" totCl: "<<fTpcClusterArray->GetEntriesFast()
                                       <<" ###########################"<<std::endl;
                                continue;
                            }
                        }
                        TVector3 pos = theCl->pos();

                        posSP.Print();
                        //pos=TpcAlignmentManager::getInstance()->localToMaster("tpc",pos);
                        if(icl==0){
                            lastPhi=pos.Phi();
                        }
                        std::cout<<"McId: "<<theCl->mcId().DominantID().mctrackID()<<", dPhi = "<<pos.Phi()-lastPhi<<"\n ";pos.Print();
                        lastPhi=pos.Phi();
                    }//end loop over clusters
                }
                delete lastRep;
                continue;
            }

        }



        hists["rHistCut"]->Fill(refPos.Perp());
        hists["momHistCut"]->Fill(refMom.Mag());

        //create residualCollection for this track
        fResults.push_back(new TpcRefResidualCollection(fTpcClusterBranchName.Data(),itr));

        //calculate residuals in XYZ --------------------------
        unsigned int ncl =0;
        if(fIdealTpcTrackTpcClusterMatching){
            ncl=cand.getNHits();

        }else{
            ncl = fTpcClusterArray->GetEntriesFast();
        }
        double refPhi=refMom.Phi();
        double refHitPhi=refPos.Phi();
        double refTheta=refMom.Theta();
        double refHitTheta=refPos.Theta();

        hists["TpcPhi"]->Fill(refPhi*rad2deg);
        hists["TpcHitPhi"]->Fill(refHitPhi*rad2deg);

        TVector3 poca, norm;
        int nTracks2=nTracks;
        if(!fSimulateTrackMatching){
            nTracks2=1;
        }
        for(unsigned int itr2=0;itr2<nTracks2;++itr2){
            //misaligned tracks
            if(fSimulateTrackMatching){
                GFTrack* itrack2 = (GFTrack*)((*fGFTrackArray)[itr2]);
                const GFAbsTrackRep* rep2;
                if(fNRep>=0){
                    rep2 = itrack2->getTrackRep(fNRep);
                }else{
                    rep2 = itrack2->getCardinalRep();
                }
                if(abs(rep2->getPDG())>10000){
                    continue;
                }
                if(((GFAbsTrackRep*)rep2)->getMom().Mag()<fMomCutRefTrack){
                    continue;
                }

                GFDetPlane tpcLastPl2 = rep2->getLastPlane();
                GFDetPlane tpcFirstPl2 = rep2->getFirstPlane();
                RKTrackRep* lastRep2 = (RKTrackRep*) rep2->clone();
                const TMatrixDSym lastCov2=lastRep2->getLastCov();
                TVectorD lastState2 = lastRep2->getLastState();
                if(fMatchOnCylinder){
                    try{
                        ((GFAbsTrackRep*)lastRep2)->extrapolate(tpcLastPl2);
                    }catch(GFException& ex){
                        std::cout<<"TpcRefGFTrackResCalc::calc(): Error during 2st Tpc extrap last point: \n"
                                <<"     "<<ex.what()<<std::endl;
                        delete lastRep2;
                        continue;
                    }
                    lastRep2->setData(lastState2,tpcLastPl2,&lastCov2,lastRep2->getAuxInfo(tpcLastPl2));
                }else{
                    lastRep2->setData(lastState2,tpcLastPl2,&lastCov2,lastRep2->getAuxInfo(lastRep2->getFirstPlane()));
                }

                TVector3 lastTpcMom2=((GFAbsTrackRep*)lastRep2)->getMom();
                TVector3 lastTpcPos2=((GFAbsTrackRep*)lastRep2)->getPos();

                if(rep2->getStatusFlag()!=0 && fRequireGoodFit){
                    continue;
                }
                if(lastTpcMom2.Mag()<fMomCutRefTrack){
                    continue;
                }
                if(fMatchOnCylinder){
                    try{

                        lastRep2->extrapolateToCylinder(15,lastTpcPos2,lastTpcMom2);
                    }catch(GFException& ex){
                        lastTpcMom2.Print();
                        std::cout<<"TpcRefGFTrackResCalc::calc(): Error during 2nd Tpc extrap to cylinder: \n"
                                <<"     "<<ex.what()<<std::endl;
                        continue;
                    }
                }
                lastTpcPos2 = TpcAlignmentManager::getInstance()->localToMaster("tpc",lastTpcPos2);
                lastTpcMom2 = TpcAlignmentManager::getInstance()->localToMasterVect("tpc",lastTpcMom2);


                GFTrackCand cand2 = itrack2->getCand();


                double phi=lastTpcMom2.Phi();
                double hitPhi=lastTpcPos2.Phi();
                double theta=lastTpcMom2.Theta();
                double hitTheta=lastTpcPos2.Theta();
                double dPhi=TpcRefGFTrkResCalc_Alignment::phiDifference(refPhi,phi);
                double dHitPhi=TpcRefGFTrkResCalc_Alignment::phiDifference(refHitPhi,hitPhi);

                hists["dPhiVsTpcPhi"]->Fill(refPhi*rad2deg,dPhi);
                hists["dPhiVsMisPhi"]->Fill(phi*rad2deg,dPhi);
                hists["dHitPhiVsTpcHitPhi"]->Fill(refHitPhi*rad2deg,dHitPhi);
                hists["dHitPhiVsMisHitPhi"]->Fill(hitPhi*rad2deg,dHitPhi);
                hists["dPhi"]->Fill(dPhi);
                hists["dHitPhi"]->Fill(dHitPhi);
                hists["MisPhi"]->Fill(phi*rad2deg);
                hists["MisHitPhi"]->Fill(hitPhi*rad2deg);
                if(rep->getCharge()>0){
                    hists["PositivedPhiVsTpcPhi"]->Fill(refPhi*rad2deg,dPhi);
                    hists["PositivedPhiVsMisPhi"]->Fill(phi*rad2deg,dPhi);
                    hists["PositivedHitPhiVsTpcHitPhi"]->Fill(refHitPhi*rad2deg,dHitPhi);
                    hists["PositivedHitPhiVsMisHitPhi"]->Fill(hitPhi*rad2deg,dHitPhi);
                    hists["PositivedPhi"]->Fill(dPhi);
                    hists["PositivedHitPhi"]->Fill(dHitPhi);
                    double zBinEnd=zStart;
                    double zBin=(zEnd-zStart)/((double)zSlices);
                    double zPos=lastTpcPos2.Z();
                    for(int i=0;i<zSlices;++i){
                        zBinEnd+=zBin;
                        //std::cout<<zBinEnd<<std::endl;
                        if(zPos<zBinEnd){
                            hists[sliceNamesPlusPos[i]]->Fill(hitPhi*rad2deg,dHitPhi);
                            hists[sliceNamesPlusMom[i]]->Fill(phi*rad2deg,dPhi);
                            break;
                        }
                    }

                }else if (rep->getCharge()<0){
                    hists["NegativedPhiVsTpcPhi"]->Fill(refPhi*rad2deg,dPhi);
                    hists["NegativedPhiVsMisPhi"]->Fill(phi*rad2deg,dPhi);
                    hists["NegativedHitPhiVsTpcHitPhi"]->Fill(refHitPhi*rad2deg,dHitPhi);
                    hists["NegativedHitPhiVsMisHitPhi"]->Fill(hitPhi*rad2deg,dHitPhi);
                    hists["NegativedPhi"]->Fill(dPhi);
                    hists["NegativedHitPhi"]->Fill(dHitPhi);

                    double zBinEnd=zStart;
                    double zBin=(zEnd-zStart)/((double)zSlices);
                    double zPos=lastTpcPos2.Z();
                    for(int i=0;i<zSlices;++i){
                        zBinEnd+=zBin;
                        //std::cout<<zBinEnd<<std::endl;
                        if(zPos<zBinEnd){
                            hists[sliceNamesMinusPos[i]]->Fill(hitPhi*rad2deg,dHitPhi);
                            hists[sliceNamesMinusMom[i]]->Fill(phi*rad2deg,dPhi);
                            break;
                        }
                    }

                }
                double dTheta=refTheta-theta;
                double dHitTheta=refHitTheta-hitTheta;

                hists["dTheta"]->Fill(dTheta*rad2deg);
                hists["dThetaVTheta"]->Fill(theta*rad2deg,dTheta*rad2deg);
                hists["dThetaVRefTheta"]->Fill(refTheta*rad2deg,dTheta*rad2deg);
                hists["dThetaPos"]->Fill(dHitTheta*rad2deg);
                hists["dThetaPosVThetaPos"]->Fill(hitTheta*rad2deg,dHitTheta*rad2deg);
                hists["dThetaPosVRefThetaPos"]->Fill(refHitTheta*rad2deg,dHitTheta*rad2deg);

                if(fabs(dPhi)<fdPhiCut&&fabs(dHitPhi)<(fdPhiCut/2)){
                    if(itr!=itr2){
                        std::cout<<"dPhi: "<<dPhi<<" dHitPhi: "<<dHitPhi<<" dTheta: :"<<dTheta<<" dHitTheta "<<dHitTheta<<" itr: "<<itr<<" itr2: " <<itr2<<std::endl;
                    }
                    cand=itrack2->getCand();
                    ncl=cand.getNHits();
                    if(abs(rep2->getPDG())>10000){
                        continue;
                    }
                    hists["dThetaCut"]->Fill(dTheta*rad2deg);
                    hists["dThetaVThetaCut"]->Fill(theta*rad2deg,dTheta*rad2deg);
                    hists["dThetaVRefThetaCut"]->Fill(refTheta*rad2deg,dTheta*rad2deg);
                    hists["dThetaPosCut"]->Fill(dHitTheta*rad2deg);
                    hists["dThetaPosVThetaPosCut"]->Fill(hitTheta*rad2deg,dHitTheta*rad2deg);
                    hists["dThetaPosVRefThetaPosCut"]->Fill(refHitTheta*rad2deg,dHitTheta*rad2deg);
                }else{
                    continue;
                }
                delete lastRep2;
            }
            for(unsigned int icl=0; icl<(ncl); icl++) {
                TMatrixT<double> cov(3,3);
                TpcCluster* theCl;
                if(fIdealTpcTrackTpcClusterMatching){
                    int detID,hitID;
                    cand.getHit(icl,detID,hitID);
                    TpcSPHit* hit= (TpcSPHit*) fTpcSPHitArray->At(hitID);
                    theCl=(TpcCluster*) fTpcClusterArray->At(hit->getClusterID());
                    if(theCl==NULL){
                        std::cout<<"########################## "
                                <<"detID: "<<detID<<" hitID: "<<hitID<<" totCl: "<<fTpcClusterArray->GetEntriesFast()
                               <<" ###########################"<<std::endl;
                        continue;
                    }
                }else{
                    theCl= (TpcCluster*) fTpcClusterArray->At(icl);
                }
                TVector3 pos =theCl->pos();
                TVector3 sig = theCl->sig();
                cov[0][0] = sig.X()*sig.X();
                cov[1][1] = sig.Y()*sig.Y();
                cov[2][2] = sig.Z()*sig.Z();
                pos = TpcAlignmentManager::getInstance()->localToMaster("tpc",pos);
                cov = TpcAlignmentManager::getInstance()->localToMaster("tpc",cov);

                double tracklength=0;
                bool ok=true;
                TVector3 pos1, mom1;
                TMatrixDSym cov6;
                try {
                    lastRep->extrapolateToPoint(pos, poca, norm);
                    TVectorT<double> dummy;
                    GFDetPlane plane(poca,norm);
                    tracklength=((GFAbsTrackRep*)lastRep)->extrapolate(plane);
                    ((GFAbsTrackRep*)lastRep)->getPosMomCov(pos1,mom1,cov6);
                }
                catch(GFException& ex) {
                    std::cout<<" TpcRefGFTrkTpcClusterResCalc_Alignment::calc(): Error during extrapolation: to point #"
                            <<icl<<", r="<<theCl->pos().Perp()<<"\n"
                           <<"     "<<ex.what()<<"\n";
                    ok=false;
                }
                if(ok) {  //add residual
                    TpcResidual* theRes = new TpcResidual();
                    theRes->setResXYZ(poca-pos);
                    theRes->setRefPos(poca);
                    theRes->setHitPos(pos);
                    theRes->setHitIndex(icl);
                    theRes->setTrackDir(mom1.Unit());
                    theRes->setHitCov(cov);
                    TMatrixD cov3=cov6.GetSub(0,2,0,2);
                    theRes->setTrackCov(cov3);
                    theRes->setTrackLength(tracklength);
                    fResults.back()->addResidual(theRes);
                }
            } //end Cluster loop
        }//end second trackfit loop
        delete lastRep;
    }//end GF trackfit loop

}
void TpcRefGFTrkTpcClusterResCalc_Alignment::writeHists(){
    TFile* file = new TFile(fQAFileName,"recreate");
    file->cd();
    for(std::map<std::string,TH1*>::iterator it=hists.begin();it!=hists.end();++it){
        (*it).second->Write();
    }
    file->Close();
}