//-----------------------------------------------------------
// 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 "TpcShieldPointTpcHitResCalc.h"
#include "TClonesArray.h"
#include "TpcCluster.h"
#include "TpcResidual.h"
#include "TpcRefResidualCollection.h"
#include "TpcSPHit.h"
#include "TpcAlignmentManager.h"
#include "TpcShieldPoint.h"

#include "GFTools.h"
#include "GFTrack.h"
#include "GFException.h"
#include "GFAbsTrackRep.h"
#include "GFRecoHitFactory.h"
#include "GFRecoHitProducer.h"
#include "RKTrackRep.h"

#include <FairRootManager.h>

#include <PndMCTrack.h>

#include <TClass.h>
#include <TString.h>
#include <TProfile.h>
#include <string>
#include <iostream>
#include <fstream>
#include <TRandom3.h>

#include <set>


double TpcShieldPointTpcHitResCalc::phiDifference(double angle1, double angle2){
    double difference=angle1-angle2;
    difference=difference*180/TMath::Pi();
    if(difference>180){
        difference-=360;
    }else if(difference<-180){
        difference+=360;
    }
    return difference;
}

TpcShieldPointTpcHitResCalc::TpcShieldPointTpcHitResCalc()
    : fTpcTrackBranchName("TpcTrackPostFit"),
      fTpcSPHitBranchName("TpcSPHit"),
      fShieldPointBranchName("TpcShieldPoint"),
      fMCTrackBranchName("MCTrack"),
      fNRep(-1),
      fDetID(-1),
      fAlignedMatching(true),
      fShieldTrackPtCut(0.5),
      fTpcTrackPtCut(0.5),
      fTrackdPhiCut(10),
      fOrigTransfName("tpc"),
      fAlMan(NULL),
      fQAFileName("TpcShieldPointTpcHitResCalcQA.root"),
      fDistortShield(false)
{
    
    initHists();
    fNExpectedBranches = 4;
    fPdg=new TDatabasePDG();
    fQA=new CdcTpcMatchingQA("TpcShieldPointTpcHitResCalcQA_matching.root");
    fQA->initHists();
    
    fMomPhiDistortions.push_back(-11.25/2);
    fMomPhiDistortions.push_back(0.4);
    fMomPhiDistortions.push_back(0.0);
    fPosPhiDistortions.push_back(0.);
    fPosPhiDistortions.push_back(0.4);
    fPosPhiDistortions.push_back(0.0);
}

TpcShieldPointTpcHitResCalc::~TpcShieldPointTpcHitResCalc() {
    
}

bool TpcShieldPointTpcHitResCalc::init() {
    if(fBranchMap.size() != fNExpectedBranches) {
        std::cout<<"TpcShieldPointTpcHitResCalc::init(): "<<std::endl
                <<"   ERROR: Wrong number of input branches ("
               <<fBranchMap.size()<<")!"<<std::endl;
        return true;
    }
    //check for proper data types:
    bool foundTpcTrack=false;
    bool foundShieldPoints=false;
    bool foundSPHitArray=false;
    bool foundMCTrackArray=false;
    std::string TpcTrackRefType("GFTrack");
    std::string tpcSPHitType("TpcSPHit");
    std::string tpcShieldType("TpcShieldPoint");
    std::string mcTrackType("PndMCTrack");
    
    if(fBranchMap.count(fShieldPointBranchName)){
        std::cout<<"shield branch found:"<<fShieldPointBranchName<<std::endl;
        std::string type(fBranchMap[fShieldPointBranchName]->GetClass()->GetName());
        if(type.find(tpcShieldType)<type.size()){
            fShieldPointArray=fBranchMap[fShieldPointBranchName];
            foundShieldPoints=true;
        }
    }
    if(fBranchMap.count(fTpcTrackBranchName)){
        std::cout<<"track branch found:"<<fTpcTrackBranchName<<std::endl;
        std::string type(fBranchMap[fTpcTrackBranchName]->GetClass()->GetName());
        if(type.find(TpcTrackRefType)<type.size()){
            fTpcTrackArray=fBranchMap[fTpcTrackBranchName];
            foundTpcTrack=true;
        }
    }
    if(fBranchMap.count(fTpcSPHitBranchName)){
        std::cout<<"SPHit branch found:"<<fTpcSPHitBranchName<<std::endl;
        std::string type(fBranchMap[fTpcSPHitBranchName]->GetClass()->GetName());
        if(type.find(tpcSPHitType)<type.size()){
            fTpcSPHitArray=fBranchMap[fTpcSPHitBranchName];
            foundSPHitArray=true;
        }
    }
    if(fBranchMap.count(fMCTrackBranchName)){
        std::cout<<"MCTrack branch found:"<<fMCTrackBranchName<<std::endl;
        std::string type(fBranchMap[fMCTrackBranchName]->GetClass()->GetName());
        if(type.find(mcTrackType)<type.size()){
            fMCTrackArray=fBranchMap[fMCTrackBranchName];
            foundMCTrackArray=true;
        }
    }
    
    if(!(foundTpcTrack&&foundSPHitArray&&foundShieldPoints&&foundMCTrackArray)) {
        std::cout<<"TpcShieldPointTpcHitResCalc::init(): "<<std::endl
                <<"   ERROR: Wrong Class Types in input branches!"<<std::endl;
        return true;
    }
    fAlMan=TpcAlignmentManager::getInstance();
    
    std::ofstream txtoutfile("TpcShieldPointTpcHitResCalc_settings.txt");
    txtoutfile<<"Settings for TpcShieldPointTpcHitResCalc"<<std::endl;
    txtoutfile<<"Original applied alignment: "<<fOrigTransfName<<std::endl;
    txtoutfile<<"dPhi cuts on the track tangentials used for matching: "<<fTrackdPhiCut<<std::endl;
    txtoutfile<<"dPhi cuts on the track positions used for matching: "<<fTrackdPhiCut/2<<std::endl;
    txtoutfile<<"Momentum cut, minimum momentum for shield tracks: "<<fShieldTrackPtCut<<std::endl;
    txtoutfile<<"Momentum cut, minimum momentum for Tpc tracks: "<<fTpcTrackPtCut<<std::endl;
    //All ok
    return false;
}


//main functionality
int TpcShieldPointTpcHitResCalc::calc() {
    
    //clear results
    for(unsigned int i=0; i<fResults.size(); i++){
        delete fResults[i];
    }
    
    fResults.clear();
   
    //create reference tracks
    std::set<int> usedTracks;
   int nShieldPoints  = fShieldPointArray->GetEntriesFast();
   for(int iRef=0;iRef<nShieldPoints;++iRef){
       TpcShieldPoint* shieldP = (TpcShieldPoint*) fShieldPointArray->At(iRef);
       int trackid=shieldP->GetTrackID();
       //Use only first MC-point of one ID
       if(!(usedTracks.insert(trackid).second)){
           continue;
       }
       PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(trackid);
       //remove ions++
       if(fabs(mcTrack->GetPdgCode())>10000){
           continue;
       }
       
       double refCharge=fPdg->GetParticle(mcTrack->GetPdgCode())->Charge()/3;
       
       TVector3 shieldPos(shieldP->GetX(),shieldP->GetY(),shieldP->GetZ());
       TVector3 shieldMom(shieldP->GetPx(),shieldP->GetPy(),shieldP->GetPz());
       TVector3 mcMom;
       mcMom=mcTrack->GetMomentum();
//       std::cout<<"#############################"<<std::endl;
//       std::cout<<"mcMom: "<<mcMom.Perp()<<std::endl;
//       std::cout<<"shieldMom: "<<shieldMom.Perp()<<std::endl;
                  
       fQA->fillBasicHists("cdc",shieldMom,shieldPos,refCharge);
       fQA->fillBasicHists("cdc",shieldMom,shieldPos,0);
       
       if(shieldMom.Perp()<fShieldTrackPtCut){
           continue;
       }
       
       double phi=shieldPos.Phi()*TMath::RadToDeg();
       histMap["rawPhiHist"]->Fill(phi);
             
       
       TVector3 shieldPosNew=shieldPos;
       TVector3 shieldMomNew=shieldMom;
       double momPhi=shieldMom.Phi()*TMath::RadToDeg();
       
       //has to be a function of phiInSector;
       //same direction as phiInSector, move away from center;
       double phiInSector=getPhiInSector(phi);
       double momPhiInSector=getPhiInSector(momPhi);
       double phiCor=0;
       double phiMomCor=0;
       
       if(fDistortShield){
           phiCor=getDistortion(phiInSector,fMomPhiDistortions[0],fMomPhiDistortions[1],fMomPhiDistortions[2]);
           phiMomCor=getDistortion(momPhiInSector,fPosPhiDistortions[0],fPosPhiDistortions[1],fPosPhiDistortions[2]);
       }
       
       histMap["phiCorVsPhi"]->Fill(phi,phiCor);
       histMap["momPhiCorVsPhi"]->Fill(momPhi,phiMomCor);
       histMap["rawMomPhiHist"]->Fill(momPhi);
       
       shieldPosNew.RotateZ(phiCor*TMath::DegToRad());
       shieldMomNew.RotateZ(phiMomCor*TMath::DegToRad());
       
       histMap["momPhiHist"]->Fill(shieldMomNew.Phi()*TMath::RadToDeg());
       histMap["phiHist"]->Fill(shieldPosNew.Phi()*TMath::RadToDeg());
       double phiDiff=phiDifference(shieldPos.Phi(),shieldPosNew.Phi());
       double phiDiff2=phiDifference(shieldMom.Phi(),shieldMomNew.Phi());
       histMap["InDirPhiCorVsPhi"]->Fill(phi,phiDiff);
       histMap["InDirMomPhiCorVsPhi"]->Fill(momPhi,phiDiff2);
       if(fabs(refCharge)>1){
           std::cout<<refCharge<<" "<<mcTrack->GetPdgCode()<<std::endl;
           continue;
       }
       //create ref rep using pion hypothesis!
       RKTrackRep* refRep = new RKTrackRep(shieldPosNew,shieldMomNew,refCharge*211);
       
       try{
           refRep->extrapolateToCylinder(15,shieldPosNew,shieldMomNew);
       }catch(GFException& ex){
           std::cout<<"TpcShieldPointTpcHitResCalc::calc(): Error during CDC extrap to cylinder: \n"
                   <<"     "<<ex.what()<<std::endl;
           continue;
       }
       //second mom cut
       if(shieldMomNew.Perp()<fShieldTrackPtCut){
           continue;
       }
       fQA->fillBasicHists("cdc_cutPt",shieldMomNew,shieldPosNew,refCharge);
       fQA->fillBasicHists("cdc_cutPt",shieldMomNew,shieldPosNew,0);
       
       //find TPC match
       fResults.push_back(new TpcRefResidualCollection(fTpcSPHitBranchName.Data(),iRef));
       unsigned int nTpcTracks = fTpcTrackArray->GetEntriesFast();
       for(unsigned int itr=0;itr<nTpcTracks;++itr){
           GFTrack* iTpcTrack=(GFTrack*)((*fTpcTrackArray)[itr]);
           GFAbsTrackRep* tpcRep= iTpcTrack->getCardinalRep();
           GFTrackCand cand=iTpcTrack->getCand();
           int tpcTrackId=cand.getMcTrackId();
           //MC matching
           if(tpcTrackId!=trackid){
               continue;
           }
           if(tpcRep->getStatusFlag()!=0){
               continue;
           }
           histMap["secID"]->Fill(shieldP->GetSecID());
           if(shieldP->GetSecID()!=0){
               std::cout<<"Matched track secId: "<<shieldP->GetSecID()<<std::endl;
               std::cout<<"mcTrack->GetStartVertex().Perp() "<<mcTrack->GetStartVertex().Perp() <<std::endl;
           }
           
           histMap["pdgHist_noPtCut"]->Fill(mcTrack->GetPdgCode());
           
           TVector3 tpcMom, tpcPos;
           TMatrixDSym tpcCov_;
           int tpcCharge=tpcRep->getCharge(); //not used, not interesting (tpc has no charge resolution)
           tpcRep->getPosMomCov(tpcPos,tpcMom,tpcCov_);
           fQA->fillBasicHists("tpc",tpcMom,tpcPos,refCharge);
           fQA->fillBasicHists("tpc",tpcMom,tpcPos,0);
           //should be set to zero, still there to reproduce old results
           if(tpcMom.Perp()<fTpcTrackPtCut){
               continue;
           }
           if(false){
           //Making a copy of the track rep, and giving it the last state.
//           RKTrackRep* lastTpcRep = (RKTrackRep*) tpcRep->clone(); //      delete lastTpcRep;!!!
//           TMatrixDSym lastTpcCov=lastTpcRep->getLastCov();      
//           GFDetPlane lastTpcPl = lastTpcRep->getLastPlane();
//           TVectorD lastTpcState = lastTpcRep->getLastState();
//           try{
//               lastTpcRep->setData(lastTpcState, lastTpcPl, &lastTpcCov,
//                                   lastTpcRep->getAuxInfo(lastTpcRep->getReferencePlane()));
//           }catch(GFException& ex){
//               std::cout<<"CdcTrackTpcHitResCalc::calc():during getting of AuxInfo: \n"
//                       <<"     "<<ex.what()<<std::endl;
//               delete lastTpcRep;
//               continue;
//           }
//           TVector3 lastTpcMom, lastTpcPos;
//           ((GFAbsTrackRep*)lastTpcRep)->getPosMomCov(lastTpcPos,lastTpcMom,lastTpcCov);
//           //LastMom << firstMom
//           if(lastTpcMom.Mag()<tpcMom.Mag()/4){
//               std::cout<<"################################################################################"<<std::endl;
//               std::cout<<"Tpc last state broken continue shield#: "<<iRef<<", tpc# "<<itr<<std::endl;
//               std::cout<<"First Mom: "<<tpcMom.Mag()<<", lastMom: "<<lastTpcMom.Mag()<<std::endl;
//               std::cout<<"pdg code: "<<mcTrack->GetPdgCode()<<std::endl;
//               std::cout<<"motherId: "<<mcTrack->GetMotherID()<<std::endl;
//               std::cout<<"points: "<<mcTrack->GetNPoints()<<std::endl;
//               std::cout<<"tpcMCid "<<cand.getMcTrackId()<<", shield "<<trackid<<std::endl;
               
//               TVector3 mom, pos;
//               shieldP->Momentum(mom);
//               shieldP->Position(pos);
//               std::cout<<"mom: "<<mom.Mag()<<std::endl;
//               std::cout<<"pt: "<<mom.Perp()<<std::endl;
//               std::cout<<"theta: "<<mom.Theta()*TMath::RadToDeg()<<std::endl;
//               std::cout<<" spos: ";pos.Print();
//               std::cout<<"red-chi2: "<<tpcRep->getChiSqu()/tpcRep->getNDF()<<std::endl;
//               std::cout<<"ndf: "<<tpcRep->getNDF()<<" npoints: "<<cand.getNHits()<<std::endl;
//               std::cout<<"################################################################################"<<std::endl;
//               delete lastTpcRep;
//               histMap["lastStateBroken"]->Fill(1);
//               lastTpcRep=(RKTrackRep*)tpcRep->clone();
//           }else{
//               histMap["lastStateBroken"]->Fill(0);
//           }
//           //rotate and shift tpctrack using alignment 
//           if(fOrigTransfName!="tpc"){
//               lastTpcMom=fAlMan->masterToLocalVect(fOrigTransfName,lastTpcMom);
//               lastTpcMom=fAlMan->localToMasterVect("tpc",lastTpcMom);
//               lastTpcPos=fAlMan->masterToLocal(fOrigTransfName,lastTpcPos);
//               lastTpcPos=fAlMan->localToMaster("tpc",lastTpcPos);
//               lastTpcRep->setPosMomCov(lastTpcPos,lastTpcMom,lastTpcCov);
//           }
           }
           if(fOrigTransfName!="tpc"){
               tpcMom=fAlMan->masterToLocalVect(fOrigTransfName,tpcMom);
               tpcMom=fAlMan->localToMasterVect("tpc",tpcMom);
               tpcPos=fAlMan->masterToLocal(fOrigTransfName,tpcPos);
               tpcPos=fAlMan->localToMaster("tpc",tpcPos);
               tpcRep->setPosMomCov(tpcPos,tpcMom,tpcCov_);
           }
           
           //std::cout<<" tpc pdg "<<((RKTrackRep*)tpcRep)->getPDG()<<std::endl;
           TVector3 tpcPos15=tpcPos;
           TVector3 tpcMom15=tpcMom;
           try{
               //lastTpcRep->extrapolateToCylinder(15,tpcPos,tpcMom);
               tpcRep->extrapolateToCylinder(15,tpcPos15,tpcMom15);
           }catch(GFException& ex){
               std::cout<<"TpcShieldPointTpcHitResCalc::calc(): Error during TPC extrap to cylinder: \n"
                       <<"pt: "<<tpcMom.Perp()<<", theta: "<<tpcMom.Theta()*TMath::RadToDeg()<<" pos: "<<tpcPos.Perp()<<" mcPt: "<<shieldMom.Perp()
                       <<"     \n"<<ex.what()<<std::endl;
               //delete lastTpcRep;
               continue;
           }           
           if(tpcMom.Perp()<fTpcTrackPtCut){
               //delete lastTpcRep;
               continue;
           }
           fQA->fillBasicHists("tpc_cutPt",tpcMom15,tpcPos15,refCharge);
           fQA->fillBasicHists("tpc_cutPt",tpcMom15,tpcPos15,0);
           fQA->fillMatchingHists(tpcMom15,tpcPos15,shieldMomNew,shieldPosNew,false,0);
           fQA->fillMatchingHists(tpcMom15,tpcPos15,shieldMomNew,shieldPosNew,false,refCharge);
           //ONLY MC matching
           
           //std::cout<<mcTrack->GetPdgCode()<<std::endl;
           //           double dPhiMom=phiDifference(tpcMom.Phi(),shieldMomNew.Phi());
           //           double dPhiPos=phiDifference(tpcPos.Phi(),shieldPosNew.Phi());
           //           if((fabs(dPhiPos)>(this->fTrackdPhiCut/2)) or
           //                   (fabs(dPhiMom)>this->fTrackdPhiCut)){
           //               //	std::cout<<"failing cut "<<this->fTrackdPhiCut<<std::endl;
           //               delete lastTpcRep;
           //               continue;
           //           }
           
           fQA->fillMatchingHists(tpcMom,tpcPos,shieldMomNew,shieldPosNew,true,0);
           fQA->fillMatchingHists(tpcMom,tpcPos,shieldMomNew,shieldPosNew,true,refCharge);
           histMap["pdgHist"]->Fill(mcTrack->GetPdgCode());
           
           //match, lets calculate residuals
           unsigned int nHits = cand.getNHits();
           for(int iHit=(nHits-1);iHit>=0;iHit=(iHit-1)){
               int hitId, detId;
               cand.getHit(iHit,detId,hitId);
               TpcSPHit* recoHit = (TpcSPHit*) fTpcSPHitArray->At(hitId);
               if(recoHit==NULL){
                   std::cout<<"nooooooooooooooooooooooooooooooooo "<<iHit<<" "<<hitId<<" "<<detId<<std::endl;
                   std::cout<<fTpcSPHitArray->GetEntries()<<std::endl;
               }
               TVector3 hPos;
               TMatrixT<double> cov(3,3) ;
               TVectorT<double> rc = ((GFAbsRecoHit*)recoHit)->getRawHitCoord();
               cov =((GFAbsRecoHit*)recoHit)->getRawHitCov();
               hPos = TVector3(rc[0], rc[1], rc[2]);
               //recalculate position based on initial alignment and current alignment
               if(fOrigTransfName!="tpc"){
                   hPos=fAlMan->masterToLocal(fOrigTransfName,hPos);
                   // Not needed in sim, as no TPC distortions are used
                   //                   if(hPos.Perp()>13 or hPos.Perp()<7){
                   //                       continue;
                   //                   }
                   hPos=fAlMan->localToMaster("tpc",hPos);
                   cov=fAlMan->masterToLocal(fOrigTransfName,cov);
                   cov=fAlMan->localToMaster("tpc",cov);
               }
               //extrapolations
               
               double tracklength=0;
               bool ok=true;
               TVector3 poca, norm;
               TVector3 pos1,mom1;
               TMatrixDSym cov6;
               try {
                   refRep->extrapolateToCylinder(hPos.Perp(),poca,norm);
                   TVectorT<double> dummy;
                   GFDetPlane plane(poca,norm);
                   tracklength=((GFAbsTrackRep*)refRep)->extrapolate(plane);
                   ((GFAbsTrackRep*)refRep)->getPosMomCov(pos1,mom1,cov6);
                   
               }
               catch(GFException& ex) {
                   std::cout<<"TpcRefGFTrackResCalc::calc(): Error during extrapolation: \n"
                           <<"     "<<ex.what()<<std::endl;
                   ok=false;
               }
               if(ok) {  //add residual
                   TpcResidual* theRes = new TpcResidual();
                   theRes->setResXYZ(poca-hPos);
                   theRes->setRefPos(poca);
                   theRes->setHitPos(hPos);
                   theRes->setTrackDir(mom1.Unit());
                   //TODO, use correct assignment (cluster/SPHit)
                   theRes->setHitIndex(hitId);
                   TMatrixDSym cov3=cov6.GetSub(0,2,0,2);
                   //cov3.Print();
                   theRes->setHitCov(cov);
                   theRes->setTrackCov(cov3);
                   theRes->setTrackLength(tracklength);
                   theRes->setCharge(refCharge);
                   fResults.back()->addResidual(theRes);
               }
               
           }//end looping over hits
           
       }//end looping over tpc tracks
       delete refRep;
   }//end looping over shield-points
   
   
}
double TpcShieldPointTpcHitResCalc::getPhiInSector(double phi){
    double shiftedPhi360=phi-11.25;
    if(shiftedPhi360<0){
       shiftedPhi360+=360;
    }
    return fmod(shiftedPhi360,22.5)-11.25;
}

double TpcShieldPointTpcHitResCalc::getDistortion(double phiInSector,double phase, double magnitude,double offset){
    //double phiCor=0.5*TMath::Cos(TMath::Pi()*2/22.5*(phiInSector));
    double phiCor=(magnitude*TMath::Cos(TMath::Pi()*2/22.5*(phiInSector-phase))+offset);
    return phiCor;//+gRandom->Gaus(phiCor,0.2)+gRandom->Uniform(2)-1;
    
}

void TpcShieldPointTpcHitResCalc::initHists(){
    TString histname="phiCorVsPhi";
    TH2D* phiCorVsPhi=new TH2D(histname,"#phi_{correction} Vs #phi;#phi (deg);d#phi (deg)",500,-180,180,500,-5,5);
    histMap[histname]=phiCorVsPhi;
    
    histname="InDirPhiCorVsPhi";
    TH2D* phiCorVsPhi2=new TH2D(histname,"#phi_{Raw}-#phi Vs #phi;#phi (deg);d#phi (deg)",500,-180,180,500,-5,5);
    histMap[histname]=phiCorVsPhi2;
    
    histname="momPhiCorVsPhi";
    TH2D* momPhiCorVsPhi=new TH2D(histname,"mom#phi_{correction} Vs mom #phi;mom d#phi (deg);mom #phi (deg)",500,-180,180,500,-5,5);
    histMap[histname]=momPhiCorVsPhi;
    
    histname="InDirMomPhiCorVsPhi";
    TH2D* momPhiCorVsPhi2=new TH2D(histname,"mom#phi_{Raw}-#phi Vs mom#phi;mom#phi (deg);momd#phi (deg)",500,-180,180,500,-5,5);
    histMap[histname]=momPhiCorVsPhi2;
    
    histname="rawPhiHist";
    TH1D* rawPhiHist = new TH1D(histname,"raw #phi; raw #phi (deg)",500,-180,180);
    histMap[histname]=rawPhiHist;
    
    histname="phiHist";
    TH1D* phiHist = new TH1D(histname,"#phi;#phi (deg)",500,-180,180);
    histMap[histname]=phiHist;
    
    histname="rawMomPhiHist";
    TH1D* rawMomPhiHist = new TH1D(histname,"raw Mom#phi; raw Mom#phi (deg)",500,-180,180);
    histMap[histname]=rawMomPhiHist;
    
    histname="momPhiHist";
    TH1D* momPhiHist = new TH1D(histname,"Mom#phi;Mom#phi (deg)",500,-180,180);
    histMap[histname]=momPhiHist;    
    
    histname="pdgHist";
    TH1D* pdgHist= new TH1D(histname,"pdg code",10000,-5000,5000);
    histMap[histname]=pdgHist;
    
    histname="pdgHist_noPtCut";
    TH1D* pdgHist2= new TH1D(histname,"pdg code no pt cut",10000,-5000,5000);
    histMap[histname]=pdgHist2;
    
    histname="secID";
    TH1D* secID= new TH1D(histname,"secId",10,0,10);
    histMap[histname]=secID;
    
    histname="lastStateBroken";
    TH1D* lastStateBroken=new TH1D(histname,histname,11,0,10);
    histMap[histname]=lastStateBroken;
}

void TpcShieldPointTpcHitResCalc::writeHists(){
    TFile* fHistoFile= new TFile(fQAFileName,"recreate");
    for(std::map<TString,TH1*>::iterator it =histMap.begin();it!=histMap.end();++it){
        it->second->Write();
        if(it->second->InheritsFrom("TH2D")){
            TProfile* temp = ((TH2D*)it->second)->ProfileX();
            temp->Write();
        }
    }
    fHistoFile->Close();
    fQA->writeHists();
}
  




ClassImp(TpcShieldPointTpcHitResCalc)