//-----------------------------------------------------------
// 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 "CdcTpcMatchingQA.h"
#include "TClonesArray.h"
#include "TpcCluster.h"
#include "TpcResidual.h"
#include "TpcRefResidualCollection.h"
#include "TpcSPHit.h"
#include "TpcAlignmentManager.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 <TClass.h>
#include <TString.h>
#include <TProfile.h>
#include <string>
#include <iostream>
#include <fstream>




double CdcTpcMatchingQA::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;
}

CdcTpcMatchingQA::CdcTpcMatchingQA()
  : fQAFileName("CdcTpcMatchingQA.root")
{
    
  initHists();

}
CdcTpcMatchingQA::CdcTpcMatchingQA(TString name)
  : fQAFileName(name)
{
    
  initHists();

}

CdcTpcMatchingQA::~CdcTpcMatchingQA() {

}


void CdcTpcMatchingQA::writeHists(){
    
  fHistoFile->cd();
  
  for(std::map<TString,TH1*>::iterator it =histMap.begin();it!=histMap.end();++it){
    fHistoFile->cd();
    fHistoFile->cd(folderNameMap[it->first]);
    it->second->Write();
    if(it->second->InheritsFrom("TH2D")){
      TProfile* temp = ((TH2D*)it->second)->ProfileX();
      temp->Write();
    }
  }
  fHistoFile->Close();
}

void CdcTpcMatchingQA::initHists(){
  fHistoFile = new TFile(fQAFileName,"recreate");
  
  initBasicHists("tpc","",0);
  initBasicHists("tpc","",1);
  initBasicHists("tpc","",-1);
  initBasicHists("tpc","Pt",0);
  initBasicHists("tpc","Pt",-1);
  initBasicHists("tpc","Pt",1);
  initBasicHists("tpc","Match",0);
  initBasicHists("tpc","Match",-1);
  initBasicHists("tpc","Match",1);
  initBasicHists("cdc","",0);
  initBasicHists("cdc","",-1);
  initBasicHists("cdc","",1);
  initBasicHists("cdc","Pt",0);
  initBasicHists("cdc","Pt",-1);
  initBasicHists("cdc","Pt",1);
  initBasicHists("cdc","Match",0);
  initBasicHists("cdc","Match",-1);
  initBasicHists("cdc","Match",1);
  //for(std::map<TString,TH1*>::iterator it =histMap.begin();it!=histMap.end();++it){
  //it->second->Print();
  //}
  initMatchingHists(true,0);
  initMatchingHists(true,1);
  initMatchingHists(true,-1);
  initMatchingHists(false,0);
  initMatchingHists(false,1);
  initMatchingHists(false,-1);
  //################################################################################
  
  TH1D* cdcExtrapZ_before=new TH1D("cdc_poca_z_before","cdc_poca_z z;z (cm)",1000,0,150);
  histMap["cdc_poca_z_before"]=cdcExtrapZ_before;
  folderNameMap["cdc_poca_z_before"]="";

  TH1D* cdcExtrapZ_after=new TH1D("cdc_poca_z_after","cdc_poca_z z;z (cm)",1000,0,150);
  histMap["cdc_poca_z_after"]=cdcExtrapZ_after;
  folderNameMap["cdc_poca_z_after"]="";

  TH1D* cdcExtrapZ_perp_before=new TH1D("cdc_poca_z_perp_before","cdc_poca_z perp;perp (cm)",1000,0,150);
  histMap["cdc_poca_z_perp_before"]=cdcExtrapZ_perp_before;
  folderNameMap["cdc_poca_z_before"]="";

  TH1D* cdcExtrapZ_perp_after=new TH1D("cdc_poca_z_perp_after","cdc_poca_z perp;perp (cm)",1000,0,150);
  histMap["cdc_poca_z_perp_after"]=cdcExtrapZ_perp_after;
  folderNameMap["cdc_poca_z_perp_after"]="";

}


void CdcTpcMatchingQA::fillBasicHists(TString prefix, TVector3 mom, TVector3 pos, int charge){
  if(!(prefix=="cdc" or prefix=="cdc_cutPt" or prefix=="tpc" or prefix=="tpc_cutPt" or
       prefix=="cdc_cutMatch" or  prefix=="tpc_cutMatch")){
    std::cerr<<"void CdcTpcMatchingQA::fillBasicHists() called with invalid prefix: "
	     <<prefix<<std::endl
	     <<"Valid prefixes are: cdc, cdc_cutPt, tpc, tpc_cutPt"<<std::endl;
    throw;
  }
  if(charge>0){
    prefix+="_plus";
  }else if(charge==0){
    
  }else{
    prefix+="_minus";
  }
  prefix+="_";

  histMap[prefix+"R"]->Fill(pos.Perp());
  histMap[prefix+"XY"]->Fill(pos.X(),pos.Y());
  histMap[prefix+"Mom"]->Fill(mom.Mag());
  histMap[prefix+"Pt"]->Fill(mom.Perp());
  histMap[prefix+"PhiPos"]->Fill(pos.Phi()*TMath::RadToDeg());
  histMap[prefix+"PhiMom"]->Fill(mom.Phi()*TMath::RadToDeg());
  histMap[prefix+"ThetaPos"]->Fill(pos.Theta()*TMath::RadToDeg());
  histMap[prefix+"ThetaMom"]->Fill(mom.Theta()*TMath::RadToDeg());
  double dPhi=phiDifference(pos.Phi(),mom.Phi()); //in degree!
  histMap[prefix+"dPhiMomPos"]->Fill(dPhi);
  
}


void CdcTpcMatchingQA::initBasicHists(TString detector, TString cut, int charge){
  TString folderName;
  TString cutstring;
  if(cut==""){
    cutstring="";
    folderName=detector+"_before_cuts";
  }else{
    cutstring="_cut"+cut;
    folderName=detector+"_after_cut"+cut;
  }
  if(charge>0){
    cutstring+="_plus";
    folderName+="/plus";
  }else if(charge==0){
    
  }else{
    cutstring+="_minus";
    folderName+="/minus";
  }
  fHistoFile->cd();
  fHistoFile->mkdir(folderName);
  std::cout<<"detector: "<<detector<<", cut: "<<cutstring<<", "<<cut<<std::endl;
    
  TH1D* R=new TH1D(detector+cutstring+"_"+"R","r_{"+detector+"};r (cm)",1000,0,150);
  histMap[detector+cutstring+"_"+"R"]=R;
  folderNameMap[detector+cutstring+"_"+"R"]=folderName;
  TH2D* XY=new TH2D(detector+cutstring+"_"+"XY","XY_{"+detector+"};x (cm);y (cm)",1000,-150,150,1000,-150,150);
  histMap[detector+cutstring+"_"+"XY"]=XY;
  folderNameMap[detector+cutstring+"_"+"XY"]=folderName;
  TH1D* Mom =new TH1D(detector+cutstring+"_"+"Mom","p_{"+detector+"};p (GeV/c)",1000,0,10);
  histMap[detector+cutstring+"_"+"Mom"]=Mom;
  folderNameMap[detector+cutstring+"_"+"Mom"]=folderName;
  TH1D* Pt =new TH1D(detector+cutstring+"_"+"Pt","p_{t,"+detector+"};p_{t} (GeV/c)",1000,0,10);
  histMap[detector+cutstring+"_"+"Pt"]=Pt;
  folderNameMap[detector+cutstring+"_"+"Pt"]=folderName;
  TH1D* PhiPos=new TH1D(detector+cutstring+"_"+"PhiPos","#phi_{pos,"+detector+"}; #phi (deg)",1000,-180,180);
  histMap[detector+cutstring+"_"+"PhiPos"]=PhiPos;
  folderNameMap[detector+cutstring+"_"+"PhiPos"]=folderName;
  TH1D* PhiMom=new TH1D(detector+cutstring+"_"+"PhiMom","#phi_{mom,"+detector+"}; #phi (deg)",1000,-180,180);
  histMap[detector+cutstring+"_"+"PhiMom"]=PhiMom;
  folderNameMap[detector+cutstring+"_"+"PhiMom"]=folderName;
  TH1D* dPhiMomPos=new TH1D(detector+cutstring+"_"+"dPhiMomPos",detector+" #phi_{pos}-#phi_{mom}; d#phi_{"+detector+"} (deg)",1000,-180,180);
  histMap[detector+cutstring+"_"+"dPhiMomPos"]=dPhiMomPos;
  folderNameMap[detector+cutstring+"_"+"dPhiMomPos"]=folderName;
  TH1D* ThetaPos=new TH1D(detector+cutstring+"_"+"ThetaPos","#theta_{pos,"+detector+"}; #theta (deg)",1000,0,180);
  histMap[detector+cutstring+"_"+"ThetaPos"]=ThetaPos;
  folderNameMap[detector+cutstring+"_"+"ThetaPos"]=folderName;
  TH1D* ThetaMom=new TH1D(detector+cutstring+"_"+"ThetaMom","#theta_{mom,"+detector+"}; #theta (deg)",1000,0,180);
  histMap[detector+cutstring+"_"+"ThetaMom"]=ThetaMom;
  folderNameMap[detector+cutstring+"_"+"ThetaMom"]=folderName;
}

void CdcTpcMatchingQA::initMatchingHists(bool cut, int charge){
  TString folderName;
  TString postfix="";
  if(cut){
    folderName="Matching_after_cuts";
    postfix="_cut";
  }else{
    folderName="Matching_before_cuts";
  }
  if(charge>0){
    postfix+="_plus";
    folderName+="/plus";
  }else if(charge==0){
    
  }else{
    postfix+="_minus";
    folderName+="/minus";
  }
  fHistoFile->mkdir(folderName);
  TH1D* dPhiPos=new TH1D("dPhiPos"+postfix,"pos: #phi_{tpc}-#phi_{cdc}; d#phi_{pos} (deg)",2000,-180,180);
  histMap["dPhiPos"+postfix]=dPhiPos;
  folderNameMap["dPhiPos"+postfix]=folderName;
  TH1D* dPhiMom=new TH1D("dPhiMom"+postfix,"mom: #phi_{tpc}-#phi_{cdc}; d#phi_{pos} (deg)",2000,-180,180);
  histMap["dPhiMom"+postfix]=dPhiMom;
  folderNameMap["dPhiMom"+postfix]=folderName;  
  TH1D* dThetaPos=new TH1D("dThetaPos"+postfix,"pos: #theta_{tpc}-#theta_{cdc}; d#theta_{pos} (deg)",2000,-180,180);
  histMap["dThetaPos"+postfix]=dThetaPos;
  folderNameMap["dThetaPos"+postfix]=folderName;
  TH1D* dThetaMom=new TH1D("dThetaMom"+postfix,"mom: #theta_{tpc}-#theta_{cdc}; d#theta_{pos} (deg)",2000,-180,180);
  histMap["dThetaMom"+postfix]=dThetaMom;
  folderNameMap["dThetaMom"+postfix]=folderName;
  
  TH1D* dZ=new TH1D("dZ"+postfix,"z_{tpc}-z_{cdc}; dz (cm)",2000,-50,50);
  histMap["dZ"+postfix]=dZ;
  folderNameMap["dZ"+postfix]=folderName;

  TH2D* dPhivPhiPos=new TH2D("dPhivPhiPos"+postfix,
			     "pos: #phi_{tpc}-#phi_{cdc} vs #phi_{tpc}; #phi_{tpc}; d#phi_{pos} (deg)",
			     500,-180,180,2000,-180,180);
  histMap["dPhivPhiPos"+postfix]=dPhivPhiPos;
  folderNameMap["dPhivPhiPos"+postfix]=folderName;
  TH2D* dPhivPhiMom=new TH2D("dPhivPhiMom"+postfix,
			     "mom: #phi_{tpc}-#phi_{cdc} vs #phi_{tpc}; #phi_{tpc}; d#phi_{mom} (deg)",
			     500,-180,180,2000,-180,180);
  histMap["dPhivPhiMom"+postfix]=dPhivPhiMom;
  folderNameMap["dPhivPhiMom"+postfix]=folderName;  
  TH2D* dThetavThetaPos=new TH2D("dThetavThetaPos"+postfix,
				 "pos: #theta_{tpc}-#theta_{cdc} vs #theta_{tpc};#theta_{tpc}; d#theta_{pos} (deg)",
				 500,-180,180,2000,-180,180);
  histMap["dThetavThetaPos"+postfix]=dThetavThetaPos;
  folderNameMap["dThetavThetaPos"+postfix]=folderName;
  TH2D* dThetavThetaMom=new TH2D("dThetavThetaMom"+postfix,
				 "mom: #theta_{tpc}-#theta_{cdc} vs #theta_{tpc};#theta_{tpc}; d#theta_{mom} (deg)",
				 500,-180,180,2000,-180,180);
  histMap["dThetavThetaMom"+postfix]=dThetavThetaMom;
  folderNameMap["dThetavThetaMom"+postfix]=folderName;


  TH1D* momRes=new TH1D("momRes"+postfix,"mom_{tpc}-mom_{cdc}; momres (MeV/c)",2000,-4,4);
  histMap["momRes"+postfix]=momRes;
  folderNameMap["momRes"+postfix]=folderName;
  TH1D* ptRes=new TH1D("ptRes"+postfix,"pt_{tpc}-pt_{cdc}; ptres (MeV/c)",2000,-4,4);
  histMap["ptRes"+postfix]=ptRes;
  folderNameMap["ptRes"+postfix]=folderName;

  TH2D* momResvMom=new TH2D("momResvMom"+postfix,"mom_{tpc}-mom_{cdc} vs mom_{cdc};mom_{cdc} (MeV/C);momres (MeV/c)",500,0,4,1000,-4,4);
  histMap["momResvMom"+postfix]=momResvMom;
  folderNameMap["momResvMom"+postfix]=folderName;
  TH2D* ptResvPt=new TH2D("ptResvPt"+postfix,"pt_{tpc}-pt_{cdc} vs pt_{cdc};pt_{cdc} (MeV/C); ptres (MeV/c)",500,0,4,1000,-4,4);
  histMap["ptResvPt"+postfix]=ptResvPt;
  folderNameMap["ptResvPt"+postfix]=folderName;

}

void CdcTpcMatchingQA::fillMatchingHists(TVector3 tpcMom, TVector3 tpcPos,TVector3 cdcMom, TVector3 cdcPos, bool cut, int charge){
  TString postfix="";
  if(cut){
    postfix="_cut";
  }
  if(charge>0){
    postfix+="_plus";
  }else if(charge==0){
    
  }else{
    postfix+="_minus";
  }
  if(cut){
    fillBasicHists("tpc_cutMatch",tpcMom,tpcPos,charge);
    fillBasicHists("cdc_cutMatch",cdcMom,cdcPos,charge);
  }
  
  double dMomPhi=phiDifference(tpcMom.Phi(),cdcMom.Phi()); //in degree!
  double dPosPhi=phiDifference(tpcPos.Phi(),cdcPos.Phi()); //in degree!
  double dMomTheta=phiDifference(tpcMom.Theta(),cdcMom.Theta()); //in degree!
  double dPosTheta=phiDifference(tpcPos.Theta(),cdcPos.Theta()); //in degree!

  histMap["dPhiPos"+postfix]->Fill(dPosPhi);
  histMap["dPhiMom"+postfix]->Fill(dMomPhi);
  histMap["dThetaPos"+postfix]->Fill(dPosTheta);
  histMap["dThetaMom"+postfix]->Fill(dMomTheta);
  histMap["dZ"+postfix]->Fill(tpcPos.Z()-cdcPos.Z());


  histMap["dPhivPhiPos"+postfix]->Fill(tpcPos.Phi()*TMath::RadToDeg(),dPosPhi);
  histMap["dPhivPhiMom"+postfix]->Fill(tpcMom.Phi()*TMath::RadToDeg(),dMomPhi);
  histMap["dThetavThetaPos"+postfix]->Fill(tpcPos.Theta()*TMath::RadToDeg(),dPosTheta);
  histMap["dThetavThetaMom"+postfix]->Fill(tpcMom.Theta()*TMath::RadToDeg(),dMomTheta);



  histMap["momRes"+postfix]->Fill(tpcMom.Mag()-cdcMom.Mag());
  histMap["ptRes"+postfix]->Fill(tpcMom.Perp()-cdcMom.Perp());
  histMap["momResvMom"+postfix]->Fill(cdcMom.Mag(),tpcMom.Mag()-cdcMom.Mag()); 
  histMap["ptResvPt"+postfix]->Fill(cdcMom.Perp(),tpcMom.Perp()-cdcMom.Perp());

}



ClassImp(CdcTpcMatchingQA)