//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class TpcQATask
//      see TpcQATask.h for details
//
// Environment:
//      Software developed for the FOPI GEM-TPC.
//
// Author List:
//      Felix Boehmer          TUM         
//
//
//-----------------------------------------------------------


// This Class' Header ------------------
#include "TpcQATask.h"

#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairRootManager.h"
#include "TClonesArray.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TVector3.h"
#include "TMatrixT.h"
#include "TMath.h"

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

#include "TpcSample.h"
#include "TpcCluster.h"
#include "TpcSPHit.h"
#include "TpcdEdx.h"
#include "MatchingPair.h"
#include "MatchingCriterion.h"
#include "MatchingTuple.h"

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

#include <iostream>


// Class Member definitions -----------

static double fTpcdEdxConversion = 1.e-3;
static double fCdcdEdxConversion = 1./1.5e3;

TpcQATask::TpcQATask()
  : FairTask("QATask")
{
  fTpcTrackBranchName = "TpcTrackPostFit";
  fTpcCdcTrackBranchName = "FopiTrackPostFit";
  fTpcSampleBranchName = "TpcSample";
  fTpcDigiBranchName = "TpcDigi";
  fTpcClusterBranchName = "TpcCluster";
  fTpcSPHitBranchName = "TpcSPHit";
  fCdcTrackBranchName = "CdcTrack";
  fCdcGFTrackBranchName = "CdcTrackPostFit";
  fCdcHitBranchName = "CdcHit";
  fTpcCdcMatchingBranchName = "TpcCdcMatchingPair";
  fTpcdEdxBranchName = "dEdx";
  fFopiTupleBranchName = "FopiTrackTuples";
  fRpcTrackBranchName = "RpcTrack";
  fTpcDec = "";
  fTpcDecNoiseName = "SigmaPadIDGlobalRef";
}


TpcQATask::~TpcQATask()
{
}

InitStatus
TpcQATask::Init()
{
  //Get ROOT Manager
  FairRootManager* ioman= FairRootManager::Instance();
  
  if(ioman==0)
    {
      Error("TpcQATask::Init","RootManager not instantiated!");
      return kERROR;
    }
  
  // Get input collection
  fTpcTrackArray=(TClonesArray*) ioman->GetObject(fTpcTrackBranchName);
  if(fTpcTrackArray==0)
    {
      Error("TpcQATask::Init","TPC Track-array not found!");
      return kERROR;
    }
  
  fCdcTrackArray=(TClonesArray*) ioman->GetObject(fCdcTrackBranchName);
  if(fCdcTrackArray==0)
    {
      Error("TpcQATask::Init","CDC Track-array not found!");
      return kERROR;
    }

  fCdcGFTrackArray=(TClonesArray*) ioman->GetObject(fCdcGFTrackBranchName);
  if(fCdcGFTrackArray==0)
    {
      Error("TpcQATask::Init","CDC GFTrack-array not found!");
      return kERROR;
    }

  fTpcCdcTrackArray=(TClonesArray*) ioman->GetObject(fTpcCdcTrackBranchName);
  if(fTpcCdcTrackArray==0)
    {
      Error("TpcQATask::Init","TPC-CDC Track-array not found!");
      return kERROR;
    }

  fCdcHitArray=(TClonesArray*) ioman->GetObject(fCdcHitBranchName);
  if(fCdcHitArray==0)
    {
      Error("TpcQATask::Init","CDC Hit-array not found!");
      return kERROR;
    }

  fTpcCdcMatchingArray=(TClonesArray*) ioman->GetObject(fTpcCdcMatchingBranchName);
  if(fTpcCdcMatchingArray==0)
    {
      Error("TpcQATask::Init","TPC-CDC Matching-array not found!");
      return kERROR;
    }
  
  fTpcClusterArray=(TClonesArray*) ioman->GetObject(fTpcClusterBranchName);
  if(fTpcClusterArray==0)
    {
      Error("TpcQATask::Init","Cluster-array not found!");
      return kERROR;
    }
  
  fTpcDigiArray=(TClonesArray*) ioman->GetObject(fTpcDigiBranchName);
  if(fTpcDigiArray==0)
    {
      Error("TpcQATask::Init","Digi-array not found!");
      return kERROR;
    }

  fTpcSampleArray=(TClonesArray*) ioman->GetObject(fTpcSampleBranchName);
  if(fTpcSampleArray==0)
    {
      Error("TpcQATask::Init","TpcSample-array not found!");
      return kERROR;
    }

  fTpcSPHitArray=(TClonesArray*) ioman->GetObject(fTpcSPHitBranchName);
  if(fTpcSPHitArray==0)
    {
      Error("TpcQATask::Init","SPHit-array not found!");
      return kERROR;
    }

  fTpcdEdxArray=(TClonesArray*) ioman->GetObject(fTpcdEdxBranchName);
  if(fTpcdEdxArray==0)
    {
      Error("TpcQATask::Init","dEdx-array not found!");
      return kERROR;
    }
  
  fFopiTupleArray=(TClonesArray*) ioman->GetObject(fFopiTupleBranchName);
  if(fFopiTupleArray==0)
    {
      Error("TpcQATask::Init","FopiTuple-array not found!");
      return kERROR;
    }

  fRpcTrackArray=(TClonesArray*) ioman->GetObject(fRpcTrackBranchName);
  if(fRpcTrackArray==0)
    {
      Error("TpcQATask::Init","RpcTrack-array not found!");
      return kERROR;
    }
  
  // setup histograms -------------------------------------------------------------------
  
  TpcTotAmpsPerPadHist = new TH1D("TpcTotAmpsPerPadHist","TPC total Amplitudes per pad",10254,0,10254);
  TpcTotAmpsPerPadHist->GetXaxis()->SetTitle("Pad Number");
  fHistos.push_back(TpcTotAmpsPerPadHist);

  if(fTpcDec.Length() > 1) {
    TpcSampleSigPerPadHist = new TH1D("TpcSampleSigPerPadHist","TPC Noise per pad",10254,0,10254);
    TpcSampleSigPerPadHist->GetXaxis()->SetTitle("Pad Number");
    fHistos.push_back(TpcSampleSigPerPadHist);
  }

  TpcSampleAmpHist = new TH1D("TpcSampleAmps","TPC Sample Amplitude distribution",2000,0,2000);
  TpcSampleAmpHist->GetXaxis()->SetTitle("Amplitude (ADC channels)");
  fHistos.push_back(TpcSampleAmpHist);

  TpcNSamplesPerDigiHist = new TH1D("TpcNSamplesPerDigi","TPC Digi Length distribution",50,0,50);
  TpcNSamplesPerDigiHist->GetXaxis()->SetTitle("# samples");
  fHistos.push_back(TpcNSamplesPerDigiHist);
  
  TpcXYSPHitHist = new TH2D("TpcXYSPHits","TPC XY SPHit distribution",100,-15,15,100,-15.5,15.5);
  TpcXYSPHitHist->GetXaxis()->SetTitle("X (cm)");
  TpcXYSPHitHist->GetYaxis()->SetTitle("Y (cm)");
  fHistos.push_back(TpcXYSPHitHist);
  
  TpcZSPHitHist = new TH1D("TpcZSPHits","TPC Z SPHit distribution",200,-70,70);
  TpcZSPHitHist->GetXaxis()->SetTitle("Z (cm)");
  fHistos.push_back(TpcZSPHitHist);

  TpcLastSPHitHist = new TH2D("TpcLastSPHits",
			      "TPC Last SPHit Position Distribution (TPC tracks)",
			      300,-70,70,100,5.,15.5);
  TpcLastSPHitHist->GetXaxis()->SetTitle("X (cm)");
  TpcLastSPHitHist->GetYaxis()->SetTitle("Radius (cm)");
  fHistos.push_back(TpcLastSPHitHist);

  
  TpcClusterSize2DHist = new TH1D("TpcClusterSize2D","TPC ClusterSize2D",50,0,50);
  TpcClusterSize2DHist->GetXaxis()->SetTitle("# Digis");
  fHistos.push_back(TpcClusterSize2DHist);

  TpcClusterSize3DHist = new TH1D("TpcClusterSize3D","TPC ClusterSize3D",50,0,50);
  TpcClusterSize3DHist->GetXaxis()->SetTitle("# Digis");
  fHistos.push_back(TpcClusterSize3DHist);
  
  TpcClusterMultHist = new TH1D("TpcClusterMults","TPC Cluster Multiplicities",300,0,300);
  TpcClusterMultHist->GetXaxis()->SetTitle("# TPC Clusters per event");
  fHistos.push_back(TpcClusterMultHist);
    
  TpcTrackMultHist = new TH1D("TpcTrackMults","TPC Track Multiplicities",30,0,30);
  TpcTrackMultHist->GetXaxis()->SetTitle("# TPC tracks");
  fHistos.push_back(TpcTrackMultHist);
  
  TpcHitMultHist = new TH1D("TpcHitMults","TPC Hit Multiplicities",100,0,100);
  TpcHitMultHist->GetXaxis()->SetTitle("# TPC hits per track");
  fHistos.push_back(TpcHitMultHist);
 
  CdcTrackMultHist = new TH1D("CdcTrackMults","CDC Track Multiplicities",30,0,30);
  CdcTrackMultHist->GetXaxis()->SetTitle("# CDC tracks");
  fHistos.push_back(CdcTrackMultHist);
  
  CdcdEdxVsMomHist = new TH2D("CdcdEdxVsMomHist","CDC dE/dx Vs. Momentum",300,0,2,300,0,10);
  CdcdEdxVsMomHist->GetXaxis()->SetTitle("CDC Momentum (GeV/c)");
  CdcdEdxVsMomHist->GetYaxis()->SetTitle("dE/dx (a.u.)");
  fHistos.push_back(CdcdEdxVsMomHist);
  
  TpcdEdxVsMomHist = new TH2D("TpcdEdxVsMomHist","TPC dE/dx [5,25] Vs. Momentum",300,0,2,300,0,10);
  TpcdEdxVsMomHist->GetXaxis()->SetTitle("Momentum (GeV/c)");
  TpcdEdxVsMomHist->GetYaxis()->SetTitle("dE/dx_{trunc} (a.u.)");
  fHistos.push_back(TpcdEdxVsMomHist);

  CdcHitMap = new TH2D("CdcHitMap","CDC Hit Map (from CDC tracks)",300,-100,100,300,-100,100);
  CdcHitMap->GetXaxis()->SetTitle("X (cm)");
  CdcHitMap->GetYaxis()->SetTitle("Y (cm)");
  fHistos.push_back(CdcHitMap);
    
  CdcHitMultHist = new TH1D("CdcHitMults","CDC Hit Multiplicities",100,0,100);
  CdcHitMultHist->GetXaxis()->SetTitle("# CDC hits per track");
  fHistos.push_back(CdcHitMultHist);

  CdcMassHist = new TH1D("CdcMasss","CDC Raw Mass*Charge Spectrum",300,-3,3);
  CdcMassHist->GetXaxis()->SetTitle("Mass*Charge (GeV/c^{2})");
  fHistos.push_back(CdcMassHist);
  
  TpcThetaHist = new TH1D("TpcThetas","TpcTrack #Theta Distribution",300,-1,1);
  TpcThetaHist->GetXaxis()->SetTitle("Tpc Track Cos(#Theta)");
  fHistos.push_back(TpcThetaHist);

  CdcThetaHist = new TH1D("CdcThetas","CdcTrack #Theta Distribution",300,-1,1);
  CdcThetaHist->GetXaxis()->SetTitle("Cdc Track Cos(#Theta)");
  fHistos.push_back(CdcThetaHist);
  
  //CdcRadiusHist = new TH1D("CdcRadiuss",
  //                         "CDC Track Radius Distribution",300,0,1000);
  //CdcRadiusHist->GetXaxis()->SetTitle("Radius (cm)");
  //fHistos.push_back(CdcRadiusHist);
  
  TpcStatusFlag = new TH1D("TpcStatusFlag",
			   "TPC CardinalRep StatusFlag",2,0,2);
  TpcStatusFlag->GetXaxis()->SetTitle("Status flag");
  fHistos.push_back(TpcStatusFlag);
  
  CdcStatusFlag = new TH1D("CdcStatusFlag",
			   "CDC CardinalRep StatusFlag",2,0,2);
  CdcStatusFlag->GetXaxis()->SetTitle("Status flag");
  fHistos.push_back(CdcStatusFlag);
  
  TpcCdcStatusFlag = new TH1D("TpcCdcStatusFlag","TPC-CDC CardinalRep StatusFlag",2,0,2);
  TpcCdcStatusFlag->GetXaxis()->SetTitle("Status flag");
  fHistos.push_back(TpcCdcStatusFlag);

  TpcChi2NDF = new TH1D("TpcChi2NDF","TPC Genfit Chi2/NDF",300,0,10);
  TpcChi2NDF->GetXaxis()->SetTitle("#chi^{2} / NDF");
  fHistos.push_back(TpcChi2NDF);

  CdcChi2NDF = new TH1D("CdcChi2NDF","CDC Genfit Chi2/NDF",300,0,10);
  CdcChi2NDF->GetXaxis()->SetTitle("#chi^{2} / NDF");
  fHistos.push_back(CdcChi2NDF);
  
  TpcCdcChi2NDF = new TH1D("TpcCdcChi2NDF","TPC-CDC Genfit Chi2/NDF",300,0,10);
  TpcCdcChi2NDF->GetXaxis()->SetTitle("#chi^{2} / NDF");
  fHistos.push_back(TpcCdcChi2NDF);

  TpcPval = new TH1D("TpcPval","TPC Genfit p-value",100,0,1);
  TpcPval->GetXaxis()->SetTitle("P-value");
  fHistos.push_back(TpcPval);

  CdcPval = new TH1D("CdcPval","CDC Genfit p-value",100,0,1);
  CdcPval->GetXaxis()->SetTitle("P-value");
  fHistos.push_back(CdcPval);
  
  TpcCdcPval = new TH1D("TpcCdcPval","TPC-CDC Genfit p-value",100,0,1);
  TpcCdcPval->GetXaxis()->SetTitle("P-value");
  fHistos.push_back(TpcCdcPval);

  TpcCdcMatchingEff = new TH1D("TpcCdcMatchingEff","TPC-CDC Track Matching Efficiency",2,0,2);
  TpcCdcMatchingEff->GetXaxis()->SetTitle("Was Matched?");
  fHistos.push_back(TpcCdcMatchingEff);

  CdcMomDiffVsTheta = new TH2D("CdcMomDiffVsTheta",
			       "CDC GF - CDC FOPI Momentum Vs. Track #Theta",
			       300,-1,1,300,-0.3,0.3);
  CdcMomDiffVsTheta->GetXaxis()->SetTitle("FOPI CDC cos(#Theta)");
  CdcMomDiffVsTheta->GetYaxis()->SetTitle("#Delta p (GeV/c)");
  fHistos.push_back(CdcMomDiffVsTheta);  

  TpcCdcDeltaZVsZ = new TH2D("TpcCdcDeltaZVsZ",
			     "CDC-TPC TrackPos #Delta z (first CDC DetPlane) Vs. Tpc z",
			     300,-40,60,300,-30,30);
  TpcCdcDeltaZVsZ->GetXaxis()->SetTitle("Tpc track z (cm)");
  TpcCdcDeltaZVsZ->GetYaxis()->SetTitle("#Delta z (cm)");
  fHistos.push_back(TpcCdcDeltaZVsZ);  

  TpcCdcDeltaThetaVsZ = new TH2D("TpcCdcDeltaThetaVsZ",
				 "CDC-TPC Track #Delta #theta Vs. Tpc z",
				 300,-40,60,300,-20,20);
  TpcCdcDeltaThetaVsZ->GetXaxis()->SetTitle("Tpc track z (cm)");
  TpcCdcDeltaThetaVsZ->GetYaxis()->SetTitle("#Delta #theta (deg.)");
  fHistos.push_back(TpcCdcDeltaThetaVsZ);  

  TpcCdcDeltaPhiVsZ = new TH2D("TpcCdcDeltaPhiVsZ",
			       "CDC-TPC Track #Delta#phi Vs. Tpc z",
			       300,-40,60,300,-20,20);
  TpcCdcDeltaPhiVsZ->GetXaxis()->SetTitle("CDC track z (cm)");
  TpcCdcDeltaPhiVsZ->GetYaxis()->SetTitle("#Delta #phi (deg.)");
  fHistos.push_back(TpcCdcDeltaPhiVsZ);    
  
  TpcCdcDeltaPhiVsPhi = new TH2D("TpcCdcDeltaPhiVsPhi","CDC-TPC Track #Delta#phi Vs. Tpc #phi",
                                 300,-180,180,300,-20,20);
  TpcCdcDeltaPhiVsPhi->GetXaxis()->SetTitle("TPC track #phi (deg)");
  TpcCdcDeltaPhiVsPhi->GetYaxis()->SetTitle("#Delta#phi (deg.)");
  fHistos.push_back(TpcCdcDeltaPhiVsPhi);
  
  TpcCdcDeltaPhiPosVsPhiPos = new TH2D("TpcCdcDeltaPhiPosVsPhiPos","CDC-TPC Track #Delta#phi_{Pos} Vs. Tpc #phi_{Pos}",
                                 300,-180,180,300,-20,20);
  TpcCdcDeltaPhiPosVsPhiPos->GetXaxis()->SetTitle("TPC track #phi_{Pos} (deg)");
  TpcCdcDeltaPhiPosVsPhiPos->GetYaxis()->SetTitle("#Delta#phi_{Pos} (deg.)");
  fHistos.push_back(TpcCdcDeltaPhiPosVsPhiPos);
  
  CdcPhiPosHist = new TH1D("CdcPhiPos", "Cdc track #phi_{pos}",300,-180,180);
  CdcPhiPosHist->GetXaxis()->SetTitle("CDC track #phi_{pos} (deg)");
  fHistos.push_back(CdcPhiPosHist);
  
  TpcPhiPosHist = new TH1D("TpcPhiPos", "TPC track #phi_{pos}",300,-180,180);
  TpcPhiPosHist->GetXaxis()->SetTitle("TPC track #phi_{pos} (deg)");
  fHistos.push_back(TpcPhiPosHist);
  
  TpcCdcPhaseSpace = new TH2D("TpcCdcPhaseSpace","TPC-CDC PhaseSpace",300,-1,1,300,0,3);
  TpcCdcPhaseSpace->GetXaxis()->SetTitle("Cos #Theta");
  TpcCdcPhaseSpace->GetYaxis()->SetTitle("Comb. momentum (GeV/c)");
  fHistos.push_back(TpcCdcPhaseSpace);
  
  RpcMomVsVelHist = new TH2D("RpcMomVsVelHist","CDC Momentum Vs RPC Velocity",300,0,2,300,0,35);
  RpcMomVsVelHist->GetYaxis()->SetTitle("RPC Velocity (cm/ns)");
  RpcMomVsVelHist->GetXaxis()->SetTitle("CDC Momentum (GeV/c)");
  fHistos.push_back(RpcMomVsVelHist);

  BarMomVsVelHist = new TH2D("BarMomVsVelHist","CDC Momentum Vs BAR Velocity",300,0,2,300,0,35);
  BarMomVsVelHist->GetYaxis()->SetTitle("BAR Velocity (cm/ns)");
  BarMomVsVelHist->GetXaxis()->SetTitle("CDC Momentum (GeV/c)");
  fHistos.push_back(BarMomVsVelHist);

  

  //set up MatchingPair info:
  fMatchVals[TString("meanResidualXY")] = std::vector<double>();
  fMatchVals[TString("dPhiHitTransition")] = std::vector<double>();
  fMatchVals[TString("dZHits")] = std::vector<double>();
  fMatchVals[TString("dPhiTrackTransition")] = std::vector<double>();
  fMatchVals[TString("dThetaTrack")] = std::vector<double>();
  fMatchVals[TString("dCharge")] = std::vector<double>();
  
  //only from matched tracks -- HAS TO HAVE EQUAL ENTRIES AS fMatchVals !!!!!!!!!
  fMatchValsM[TString("meanResidualXY")] = std::vector<double>();
  fMatchValsM[TString("dPhiHitTransition")] = std::vector<double>();
  fMatchValsM[TString("dZHits")] = std::vector<double>();
  fMatchValsM[TString("dPhiTrackTransition")] = std::vector<double>();
  fMatchValsM[TString("dThetaTrack")] = std::vector<double>();
  fMatchValsM[TString("dCharge")] = std::vector<double>();


  fMatchBins[TString("meanResidualXY")] =  300;
  fMatchBins[TString("dPhiHitTransition")] =  300;
  fMatchBins[TString("dZHits")] =  300;
  fMatchBins[TString("dPhiTrackTransition")] =  300;
  fMatchBins[TString("dThetaTrack")] =  300;
  fMatchBins[TString("dCharge")] =  5;

  std::pair<double,double> pair;
  pair.first = -18.;
  pair.second = 18.;
  fMatchRange[TString("meanResidualXY")] =  pair;
  fMatchRange[TString("dPhiHitTransition")] =  pair;
  fMatchRange[TString("dPhiTrackTransition")] =  pair;
  fMatchRange[TString("dThetaTrack")] =  pair;

  pair.first = -30.;
  pair.second = 30.;
  fMatchRange[TString("dZHits")] =  pair;

  pair.first = -2.1;
  pair.second = 2.1;
  fMatchRange[TString("dCharge")] =  pair;
  
  
 // _hitBranchID = FairRootManager::Instance()->GetBranchId(_clusterBranchName);

 //  std::cout<<"Initializing PadPlane ... "<<std::endl;
 //  //TpcPadShapePool* pool = fPar->getPadShapes();
 //  fPlane = fPar->getPadPlane();
 //  fRMin = fPar->getRMin();
 //  fRMax = fPar->getRMax();

 //  if(fIgnoreEdgeDigis) {
 //    fRMin = fRCutoff1;
 //    fRMax = fRCutoff2;
 //  }
 
  return kSUCCESS;
}


void
TpcQATask::SetParContainers() {

  std::cout<<"TpcQATask::SetParContainers"<<std::endl;
  std::cout.flush();
  
  // Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get Tpc digitisation parameter container
  fPar= (TpcDigiPar*) db->getContainer("TpcDigiPar");
  if (! fPar ) Fatal("SetParContainers", "TpcDigiPar not found");

}



void
TpcQATask::Exec(Option_t* opt)
{
  std::cout<<"TpcQATask::Exec() ... "<<std::endl;
  

  //loop over TpcSamples
  for(unsigned int is=0; is<fTpcSampleArray->GetEntriesFast(); is++) {
    TpcSample* isam = (TpcSample*) fTpcSampleArray->At(is);
    TpcSampleAmpHist->Fill(isam->amp());
    if(isam->padId()>=1.e-9) {
      TpcTotAmpsPerPadHist->Fill(isam->padId(), isam->amp());
    }
  } //end loop over TpcSamples

  
  //loop over TpcDigis
  for(unsigned int id=0; id<fTpcDigiArray->GetEntriesFast(); id++) {
    TpcDigi* idigi = (TpcDigi*) fTpcDigiArray->At(id);
    TpcNSamplesPerDigiHist->Fill(idigi->nSample());
  } //end loop over TpcDigis
  
 

  //loop over TpcClusters
  TpcClusterMultHist->Fill(fTpcClusterArray->GetEntriesFast());
  for(unsigned int ic=0; ic<fTpcClusterArray->GetEntriesFast(); ic++) {
    TpcCluster* icl = (TpcCluster*) fTpcClusterArray->At(ic);
    TpcClusterSize2DHist->Fill(icl->get2DSize());
    TpcClusterSize3DHist->Fill(icl->size());
  }
  

  //loop over TPC SPhits
  for(unsigned int isp=0; isp<fTpcSPHitArray->GetEntriesFast(); isp++) {
    TpcSPHit* isph = (TpcSPHit*)fTpcSPHitArray->At(isp);
    TVectorD iraw(3);
    iraw = isph->getRawHitCoord();
    TVector3 ipos(iraw[0],iraw[1],iraw[2]);
    TpcXYSPHitHist->Fill(ipos.X(),ipos.Y());
    TpcZSPHitHist->Fill(ipos.Z());
  } //end loop over TpcSPHits
  
  //loop over Tpc tracks
  for(unsigned int itr=0; itr<fTpcTrackArray->GetEntriesFast(); itr++) {
    TpcTrackMultHist->Fill(fTpcTrackArray->GetEntriesFast());
    GFTrack* itpctr = (GFTrack*) fTpcTrackArray->At(itr);
    TpcHitMultHist->Fill(itpctr->getHits().size());
    std::vector<GFAbsRecoHit*> hits = itpctr->getHits();
    TpcSPHit* lhit = (TpcSPHit*) hits.back();
    TVectorD lraw(3);
    lraw = lhit->getRawHitCoord();
    TVector3 lpos(lraw[0],lraw[1],lraw[2]);
    TpcLastSPHitHist->Fill(lpos.Z(),lpos.Perp());    
    TpcThetaHist->Fill(TMath::Cos(itpctr->getCardinalRep()->getMom().Theta()));
    TpcChi2NDF->Fill(itpctr->getCardinalRep()->getRedChiSqu());
    TpcPval->Fill(itpctr->getCardinalRep()->getPVal());
    TpcStatusFlag->Fill(itpctr->getCardinalRep()->getStatusFlag());
  } //end loop over Tpc tracks

  //loop over TPC dE/dx
  for(unsigned int ide=0; ide<fTpcdEdxArray->GetEntriesFast(); ide++) {
    TpcdEdx* idEdx = (TpcdEdx*) fTpcdEdxArray->At(ide);
    double mom = idEdx->getMom().Mag();
    double val = idEdx->truncMean(0.05,0.25);
    TpcdEdxVsMomHist->Fill(mom,val*fTpcdEdxConversion);
  }
  //end loop over TPC dE/dx
  

  //loop over Cdc GF tracks
  for(unsigned int itr=0; itr<fCdcGFTrackArray->GetEntriesFast(); itr++) {
    GFTrack* icdcgf = (GFTrack*) fCdcGFTrackArray->At(itr);
    CdcChi2NDF->Fill(icdcgf->getCardinalRep()->getRedChiSqu());
    CdcPval->Fill(icdcgf->getCardinalRep()->getPVal());
    CdcStatusFlag->Fill(icdcgf->getCardinalRep()->getStatusFlag());
  }
  //end loop over Cdc GF tracks

  //loop over Cdc tracks (FOPI reco output)
  for(unsigned int itr=0; itr<fCdcTrackArray->GetEntriesFast(); itr++) {
    CdcTrackMultHist->Fill(fCdcTrackArray->GetEntriesFast());
    CdcTrack* icdctr = (CdcTrack*) fCdcTrackArray->At(itr);
    CdcHitMultHist->Fill(icdctr->GetCdcIdHits().size());
    double charge = icdctr->GetCharge();
    double mass = icdctr->GetMass();
    CdcMassHist->Fill(charge*mass);
    //CdcRadiusHist->Fill(icdctr->GetRadius());
    double convdEdx = icdctr->GetEloss()*fCdcdEdxConversion;
    CdcdEdxVsMomHist->Fill(icdctr->GetMom(),convdEdx);
    CdcThetaHist->Fill(TMath::Cos(icdctr->GetTheta()));
    if(icdctr->GetVel() > 1.e-3) {
      int rpcid = icdctr->GetRpcHitID();
      int barid = icdctr->GetBarHitID();
      if(rpcid >= 0 && barid == -1)
	RpcMomVsVelHist->Fill(icdctr->GetMom(),icdctr->GetVel());
      if(barid >= 0 && rpcid == -1)
	BarMomVsVelHist->Fill(icdctr->GetMom(),icdctr->GetVel());
    }
	
    
    
    //Get CDC hits from tracks:
    std::vector<unsigned int> hits = icdctr->GetCdcIdHits();
    for(unsigned int ih=0; ih<hits.size(); ih++) {
      CdcHit* ihit = (CdcHit*) fCdcHitArray->At(hits[ih]);
      CdcHitMap->Fill(ihit->GetHitPos().X(),ihit->GetHitPos().Y());
    }
    
    
    
  } //end loop over Cdc tracks
  
  //loop over combined trackfits :
  for(unsigned int itr=0; itr<fTpcCdcTrackArray->GetEntriesFast(); itr++) {
    GFTrack* itpccdctr = (GFTrack*) fTpcCdcTrackArray->At(itr);
    TpcCdcStatusFlag->Fill(itpccdctr->getCardinalRep()->getStatusFlag());
    if(itpccdctr->getCardinalRep()->getStatusFlag())
      continue;

    TVector3 mom3 = itpccdctr->getCardinalRep()->getMom();
    //TpcCdcMomHist->Fill(mom3.Mag());
    TpcCdcPhaseSpace->Fill(TMath::Cos(mom3.Theta()),mom3.Mag());
    TpcCdcChi2NDF->Fill(itpccdctr->getCardinalRep()->getRedChiSqu());    
    TpcCdcPval->Fill(itpccdctr->getCardinalRep()->getPVal());
  } //end loop over combined TPC CDC tracks

  //loop over MatchingPairs
  for(unsigned int imp=0; imp<fTpcCdcMatchingArray->GetEntriesFast(); imp++) {
    MatchingPair* mp = (MatchingPair*) fTpcCdcMatchingArray->At(imp);
    std::vector<TString> cnames = mp->getListOfCritNames();
    for(unsigned int ic=0; ic<cnames.size(); ic++) {
      if(fMatchVals.count(cnames[ic]))  {
	fMatchVals[cnames[ic]].push_back(mp->getCriterion(cnames[ic]).getValue());
	//if(mp->getCriterion(cnames[ic]).isConstrained() && fMatchCuts.count(cnames[ic])==0) 
	//  fMatchCuts[cnames[ic]] = mp->getCriterion(cnames[ic]).getCuts();
	if(mp->isMatched()) 
	  fMatchValsM[cnames[ic]].push_back(mp->getCriterion(cnames[ic]).getValue());
	if(mp->getCriterion(cnames[ic]).isConstrained() && !fMatchCuts.count(cnames[ic])) {
	  fMatchCuts[cnames[ic]] = mp->getCriterion(cnames[ic]).getCuts();
	}
      }
    }
  } //end loop over MatchingPairs

  //loop over track tuples ... associated information
  for(unsigned int it=0; it<fFopiTupleArray->GetEntriesFast(); it++) {
    MatchingTuple* itup = (MatchingTuple*)fFopiTupleArray->At(it);
    TpcCdcMatchingEff->Fill(itup->isMatched());
    if(itup->isMatched()) {
      unsigned int tpcGFId = itup->getIdFromBranch(fTpcTrackBranchName);
      unsigned int cdcFopiId = itup->getIdFromBranch(fCdcTrackBranchName);
      GFTrack* tpcTrack = (GFTrack*) fTpcTrackArray->At(tpcGFId);
      CdcTrack* cdcTrack = (CdcTrack*) fCdcTrackArray->At(cdcFopiId);
      //try to get CDC GF Track: relies on same sorting of FOPI reco and CDCGF TCA
      GFTrack* cdcGFTrack = (GFTrack*) fCdcGFTrackArray->At(cdcFopiId);

      //TODO: investigate: How can this be empty? Failed fit somewhere?
      if(cdcGFTrack==NULL || tpcTrack==NULL) 
	continue;

      GFDetPlane tpcLastPl = tpcTrack->getCardinalRep()->getLastPlane();
      GFDetPlane cdcFirstPl = cdcGFTrack->getCardinalRep()->getFirstPlane();
      
      
      TVector3 cdcTrackPos; //track position of CDC track
      TVector3 tpcTrackPos; //track position of TPC track
      TVector3 cdcTrackMom; //track momentum of CDC track
      TVector3 tpcTrackMom; //track momentum of TPC track
      double radius=15;
      try { 
	RKTrackRep* lastRep = (RKTrackRep*) tpcTrack->getCardinalRep()->clone();
	if(!tpcTrack->getSmoothing()) {
	  lastRep->setData(lastRep->getLastState(),tpcLastPl,
                           &lastRep->getLastCov());
	  //tpcTrack->getCardinalRep()->getPosMom(tpcLastPl,
	  //				      tpcTrackPos,tpcTrackMom);
	}
	else {
	  TVectorD smoothed_state;
	  TMatrixDSym smoothed_cov;
	  TMatrixD auxInfo;
	  GFDetPlane smoothing_plane;
	  unsigned int nHits = tpcTrack->getNumHits();
	  if(nHits==0) 
	    nHits=1;
	  
	  GFTools::getBiasedSmoothedData(tpcTrack, 0, nHits-1, smoothed_state, smoothed_cov,
					 smoothing_plane, auxInfo);
	  
	  if(lastRep->hasAuxInfo()) {
	    lastRep->setData(smoothed_state, smoothing_plane, &smoothed_cov, &auxInfo);
	  } else {
	    lastRep->setData(smoothed_state, smoothing_plane, &smoothed_cov);
	  }
	}
	lastRep->extrapolateToCylinder(radius,tpcTrackPos,tpcTrackMom);
	delete lastRep;
      }
      catch (GFException ex) {
	if (ex.isFatal()) 
	  continue;
      }
      try { 
          ((RKTrackRep*)cdcGFTrack->getCardinalRep())->extrapolateToCylinder(radius,cdcTrackPos,cdcTrackMom);
      }
      catch (GFException ex) {
          if (ex.isFatal()) 
              continue;
      }
      //this shouldn't actually extrapolate anything, but let's be safe

      TVector3 delta = cdcTrackPos - tpcTrackPos;
      double cdcTheta = cdcTrackMom.Theta();
      double tpcTheta = tpcTrackMom.Theta();
      double cdcPhi = cdcTrackMom.Phi();
      double tpcPhi = tpcTrackMom.Phi();
      double cdcPosPhi = cdcTrackPos.Phi();
      double tpcPosPhi = tpcTrackPos.Phi();
      
      double deltaTheta = cdcTheta-tpcTheta;
      double deltaPhi = cdcPhi-tpcPhi;
      double deltaPhiPos=cdcPosPhi-tpcPosPhi;
      if(deltaPhi>180){
          deltaPhi-=360;
      }else if(deltaPhi<-180){
          deltaPhi+=360;
      }
      if(deltaPhiPos>180){
          deltaPhiPos-=360;
      }else if(deltaPhiPos<-180){
          deltaPhiPos+=360;
      }
      double momDiff = cdcTrackMom.Mag()-cdcTrack->GetMom();
      TpcCdcDeltaZVsZ->Fill(tpcTrackPos.Z(),delta.Z());
      TpcCdcDeltaThetaVsZ->Fill(tpcTrackPos.Z(),deltaTheta*180./TMath::Pi());
      TpcCdcDeltaPhiVsZ->Fill(tpcTrackPos.Z(),deltaPhi*180./TMath::Pi());
      CdcMomDiffVsTheta->Fill(TMath::Cos(cdcTrack->GetTheta()),momDiff);
      TpcCdcDeltaPhiVsPhi->Fill(tpcPhi*180,deltaPhi*180);
      TpcCdcDeltaPhiPosVsPhiPos->Fill(tpcPosPhi*180,deltaPhiPos*180);
      TpcPhiPosHist->Fill(tpcPosPhi*180);
      CdcPhiPosHist->Fill(cdcPosPhi*180);
    }
  }
  //end loop over track tuples

  

  //loop over RPC tracks
  
  //COMMENTED out for now since RPC vel is stored in CdcTrack since rev155

  // for(unsigned int it=0; it<fRpcTrackArray->GetEntriesFast(); it++) {
  //   RpcTrack* irpc = (RpcTrack*)fRpcTrackArray->At(it);
  //   if(irpc->GetVel() > 1.e-2) {
  //     CdcTrack* tt = (CdcTrack*) fCdcTrackArray->At(it);
  //     if(tt!=NULL) {
  // 	RpcMomVsVelHist->Fill(irpc->GetVel(),tt->GetMom());
  //     }
  //   }
  // }
  //end loop over RPC tracks
  								       
}



void 
TpcQATask::FinishTask() {
  std::cout<<"TpcQATask::FinishTask()  ... writing histograms to subfolder 'TpcQA'"<<std::endl;
  
  //set proper scale for some plots
  std::vector<TH1*> corrs;
  corrs.push_back(TpcStatusFlag);
  corrs.push_back(TpcCdcStatusFlag);
  corrs.push_back(TpcCdcMatchingEff);
  corrs.push_back(CdcStatusFlag);
  for(unsigned int i=0; i<corrs.size(); i++) {
    TH1* icorr = corrs[i];
    double max = icorr->GetMaximum();
    icorr->GetYaxis()->SetRangeUser(0.,1.1*max);     
  }
  
  
  //if decoded file was specified, get noise data from there:
  if(fTpcDec.Length() > 1) {
    TFile* tfile = new TFile(fTpcDec,"read");
    TH1D* noise = (TH1D*) tfile->Get(fTpcDecNoiseName);
    if(noise==NULL) {
      Error("TpcQATask::FinishTask","Could not get noise histogram from TPC decoded file!");
    }
    else 
      TpcSampleSigPerPadHist->Add(noise);
    tfile->Close();
  }
  
  //Create a folder in the output file and store all histograms
  FairRootManager* ioman = FairRootManager::Instance();
  TFile* outf = ioman->GetOutFile();
  TDirectory* qadir = outf->mkdir("TpcQA");
  qadir->cd();
  
  for(unsigned int ih=0; ih<fHistos.size(); ih++) {
    if(fHistos[ih] == NULL) 
      continue;
    (fHistos[ih])->Write();
  }

  //generate histograms for matching criteria:
  std::map<TString, std::vector<double> >::const_iterator it;
  for(it=fMatchVals.begin(); it!=fMatchVals.end(); it++) {
    TString tempname("TpcCdcMatching_");
    TString name = tempname+it->first;
    std::pair<double,double> range = fMatchRange[it->first];
            
    if(fMatchCuts.count(it->first) ) {
      char add[] = "_cutl_%.2f_cuth_%.2f";
      char buff[40];
      std::pair<double,double> cuts = fMatchCuts[it->first];
      sprintf(buff,add,cuts.first,cuts.second);
      name += TString(buff);
    }
    TH1D* hist = new TH1D(name,name,fMatchBins[it->first],range.first,range.second);
    TString nameM = name+"_matched";
    
    TH1D* histm = new TH1D(nameM,nameM,fMatchBins[it->first],range.first,range.second);
    std::vector<double> valsM = fMatchValsM[it->first];
    for(unsigned int ival=0; ival<it->second.size();ival++) 
      hist->Fill((it->second)[ival]);
    
    for(unsigned int ivalm=0; ivalm<valsM.size(); ivalm++ )
      histm->Fill(valsM[ivalm]);
    
    hist->Write();
    histm->Write();
  }
    
}



ClassImp(TpcQATask)