// -------------------------------------------------------------------------
// -----                  CbmLitMuchRecQa source file               -----
// -----                  Created 15/10/07  by A. Lebedev               -----
// -------------------------------------------------------------------------

#include "CbmLitMuchRecQa.h"

#include "CbmMuchTrack.h"
#include "CbmMuchTrackMatch.h"

#include "CbmGeoNode.h"
#include "CbmGeoPassivePar.h"
#include "CbmGeoVector.h"
#include "CbmMCTrack.h"
#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmRuntimeDb.h"
#include "CbmStsTrack.h"
#include "CbmStsTrackMatch.h"

#include "TClonesArray.h"
#include "TDirectory.h"
#include "TH1F.h"
#include "TList.h"
#include "TVector3.h"

#include <iostream>
#include <map>
using std::cout;
using std::endl;
using std::map;

// -----   Default constructor   -------------------------------------------
CbmLitMuchRecQa::CbmLitMuchRecQa()
  : CbmTask("Much Track Finder QA", 1) { 
  fMinPoints = 15;
  fQuota     = 0.7;
  fNormType = 1;
  fVerbose = 1;
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmLitMuchRecQa::CbmLitMuchRecQa(Int_t minPoints, Double_t quota,
                                          Int_t iVerbose) 
  : CbmTask("Much Track Finder QA", iVerbose) {
  fMinPoints = minPoints;
  fQuota     = quota;
  fNormType = 1;
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmLitMuchRecQa::~CbmLitMuchRecQa() { 

  fHistoList->Delete();
  delete fHistoList;
}
// -------------------------------------------------------------------------



// -----   Public method SetParContainers   --------------------------------
void CbmLitMuchRecQa::SetParContainers() {

  // Get Run
  CbmRunAna* run = CbmRunAna::Instance();
  if ( ! run ) {
    cout << "-E- " << GetName() << "::SetParContainers: No CbmRunAna!"
         << endl;
    return;
  }

  // Get Runtime Database
  CbmRuntimeDb* runDb = run->GetRuntimeDb();
  if ( ! run ) {
    cout << "-E- " << GetName() << "::SetParContainers: No runtime database!"
         << endl;
    return;
  }

  // Get passive geometry parameters
  fPassGeo = (CbmGeoPassivePar*) runDb->getContainer("CbmGeoPassivePar");
  if ( ! fPassGeo ) {
    cout << "-E- " << GetName() << "::SetParContainers: "
         << "No passive geometry parameters!" << endl;
    return;
  }
}
// -------------------------------------------------------------------------



// -----   Public method Init   --------------------------------------------
InitStatus CbmLitMuchRecQa::Init() {

  cout << "==========================================================="
       << endl;;
  cout << GetName() << ": Initialising..." << endl;

  // Get RootManager
  CbmRootManager* ioman = CbmRootManager::Instance();
  if (! ioman) {
    cout << "-E- " << GetName() << "::Init: "
         << "RootManager not instantised!" << endl;
    return kFATAL;
  }

  // Get MCTrack array
  fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  if ( ! fMCTracks ) {
    cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl;
    return kFATAL;
  }

  if (fNormType == 2) {
    // Get StsTrack array
    fStsTracks = (TClonesArray*) ioman->GetObject("STSTrack");
    if ( ! fStsTracks ) {
      cout << "-E- " << GetName() << "::Init: No StsTrack array!" << endl;
      return kERROR;
    }

    // Get StsTrackMatch array
    fStsMatches = (TClonesArray*) ioman->GetObject("STSTrackMatch");
    if ( ! fStsMatches ) {
      cout << "-E- " << GetName() << "::Init: No StsTrackMatch array!"
           << endl;
      return kERROR;
    }
  }

  // Get MuchTrack array
  fMuchTracks = (TClonesArray*) ioman->GetObject("MuchTrack");
  if ( ! fMuchTracks ) {
    cout << "-E- " << GetName() << "::Init: No MuchTrack array!" << endl;
    return kERROR;
  }

  // Get MuchTrackMatch array
  fMuchMatches = (TClonesArray*) ioman->GetObject("MuchTrackMatch");
  if ( ! fMuchMatches ) {
    cout << "-E- " << GetName() << "::Init: No MuchTrackMatch array!"
         << endl;
    return kERROR;
  }


  // Get the geometry of target
  InitStatus geoStatus = GetGeometry();
  if ( geoStatus != kSUCCESS ) {
    cout << "-E- " << GetName() << "::Init: Error in reading geometry!"
         << endl;
    return geoStatus;
  }

  // Create histograms
  CreateHistos();
  Reset();

  // Output
  cout << "   Minimum number of Much points : " << fMinPoints << endl;
  cout << "   Matching quota               : " << fQuota << endl;
  if (fActive) cout << "   *****   Task is ACTIVE   *****" << endl;
  cout << "==========================================================="
       << endl << endl;

  return geoStatus;

}
// -------------------------------------------------------------------------



// -----   Public method ReInit   ------------------------------------------
InitStatus CbmLitMuchRecQa::ReInit() {

  cout << "==========================================================="
       << endl;;
  cout << GetName() << ": Reinitialising..." << endl;

  // Get the geometry of target
  InitStatus geoStatus = GetGeometry();
  if ( geoStatus != kSUCCESS ) {
    cout << "-E- " << GetName() << "::ReInit: Error in reading geometry!"
         << endl;
    return geoStatus;
  }

  // Output
  cout << "   Target position ( " << fTargetPos.X() << ", "
       << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl;
  cout << "==========================================================="
       << endl << endl;

  return geoStatus;
}
// -------------------------------------------------------------------------



// -----   Public method Exec   --------------------------------------------
void CbmLitMuchRecQa::Exec(Option_t* opt) {

  fMatchSetSts.clear();

  if (fNormType == 2) { // Normalization to number of Sts AND Much tracks
    Int_t nStsTrack = fStsTracks->GetEntriesFast();
    for(Int_t iStsTrack = 0; iStsTrack < nStsTrack; iStsTrack++) {

      CbmStsTrackMatch* stsTrackM = (CbmStsTrackMatch*) fStsMatches->
                        At(iStsTrack);
      if ( ! stsTrackM ) {
        cout << "-E- " << GetName() << "::Exec: "
	     << "No StsTrackMatch at index " << iStsTrack << endl;
        Fatal("Exec", "No StsTrackMatch in array");
      }

      Int_t iMC = stsTrackM->GetMCTrackId();
      
      if(iMC == -1) continue;

      fMatchSetSts.insert(iMC);
    }
  }

  // Clear matching map and reset eventwise counters
  fMatchMap.clear();
  fQualiMap.clear();

  Int_t nNoSts = 0;
  
  // Loop over MuchTracks. Check matched MCtrack and fill maps.
  Int_t nGhosts = 0;
  Int_t nClones = 0;
  Int_t nRec = fMuchTracks->GetEntriesFast();
  Int_t nMtc = fMuchMatches->GetEntriesFast();
  if ( nMtc != nRec ) {
     Fatal("Exec", "Inequal number of MuchTrack and MuchTrackMatch");
  }
  
  
  Int_t nofRecMuonsSts = 0;
  Int_t nofRecMuonPairsSts = 0;
  for (Int_t iRec = 0; iRec < nRec; iRec++) {

    CbmMuchTrack* MuchTrack = (CbmMuchTrack*) fMuchTracks->At(iRec);
    if ( ! MuchTrack ) {
      Fatal("Exec", "No MuchTrack in array");
    }
    
    
    Int_t stsTrackId = MuchTrack->GetStsTrackID();
    CbmStsTrackMatch* stsTrackM = (CbmStsTrackMatch*) fStsMatches->At(stsTrackId);
    Int_t mcId = stsTrackM->GetMCTrackId();
    CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcId);

    if ( TMath::Abs(mcTrack->GetPdgCode()) == 13 &&
         mcTrack->GetMotherId() == -1 ) {
         
         nofRecMuonsSts++;
         if (nofRecMuonsSts > 1) nofRecMuonPairsSts++;
    } 
    
    
    Int_t nHits = MuchTrack->GetNHits();

    CbmMuchTrackMatch* match = (CbmMuchTrackMatch*) fMuchMatches->At(iRec);
    if ( ! match ) {
      Fatal("Exec", "No MuchTrackMatch in array");
    }
    Int_t nTrue = match->GetNofTrueHits();

    Int_t iMC = match->GetMCTrackId();
    
    if (iMC == -1 ) {       // no common point with MC, really ghastly!
      if ( fVerbose > 1 )
      cout << "-I- " << GetName() << ":"
           << "No MC match for MuchTrack " << iRec << endl;
      fhNhGhosts->Fill(nHits);
      nGhosts++;
      continue;
    }
      

    // Check matching criterion (quota)
    Double_t quali = Double_t(nTrue) / Double_t(nHits);
    if ( quali >= fQuota ) { 
      
      // Good MuchTrack w/o correctly matched STS track
      if ( fNormType == 2 &&
        (fMatchSetSts.find(iMC) == fMatchSetSts.end()) ) nNoSts++;
      
      // No previous match for this MCTrack
      if ( fMatchMap.find(iMC) == fMatchMap.end() ) {
         fMatchMap[iMC] = iRec;
         fQualiMap[iMC] = quali;
      } else {  // Previous match; take the better one
         
         if ( fVerbose > 1 )
          cout << "-I- " << GetName() << ": "
               << "MCTrack " << iMC << " doubly matched."
               << "Current match " << iRec
               << ", previous match " << fMatchMap[iMC]
               << endl;
               
         if ( fQualiMap[iMC] < quali ) {
            CbmMuchTrack* oldTrack
                  = (CbmMuchTrack*) fMuchTracks->At(fMatchMap[iMC]);
            fhNhClones->Fill(Double_t(oldTrack->GetNHits()));
            fMatchMap[iMC] = iRec;
            fQualiMap[iMC] = quali;
         } else fhNhClones->Fill(nHits);
         
         nClones++;
      }
    } else { // If not matched, it's a ghost
      
      if ( fVerbose > 1 )
         cout << "-I- " << GetName() << ":"
              << "MuchTrack " << iRec << " below matching criterion "
              << "(" << quali << ")" << endl;
      fhNhGhosts->Fill(nHits);
      nGhosts++;
      
    }
  }   // Loop over MuchTracks


  // Loop over MCTracks
  Int_t nAll     = 0;
  Int_t nAcc     = 0;
  Int_t nRecAll  = 0;
  Int_t nPrim    = 0;
  Int_t nRecPrim = 0;
  Int_t nRef     = 0;
  Int_t nRecRef  = 0;
  Int_t nSec     = 0;
  Int_t nRecSec  = 0;
  Int_t nRecSignalPair = 0;
  Int_t nAccSignalPair = 0;

  Int_t accSignalCnt = 0;
  Int_t recSignalCnt = 0;
  TVector3 vertex;
  Int_t nMC = fMCTracks->GetEntriesFast();
  for (Int_t iMC = 0; iMC < nMC; iMC++) {
  
    CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(iMC);
    if ( ! mcTrack ) {
       Fatal("Exec", "No MCTrack in array");
    }

    // Check geometrical acceptance; continue only for accepted tracks
    nAll++;
    Int_t nPoints = mcTrack->GetNPoints(kMUCH);

    if ( nPoints < fMinPoints ) continue; /////!!!!!! <

    if ( fNormType == 2 &&
        (fMatchSetSts.find(iMC) == fMatchSetSts.end()) ) continue;

    if ( TMath::Abs(mcTrack->GetPdgCode()) != 13 ) continue;
    if ( mcTrack->GetMotherId() != -1) continue;
    
    accSignalCnt++;
    
    if (accSignalCnt == 2) {
      nAccSignalPair++;
    }
    
    nAcc++;

    // Check origin of MCTrack
    mcTrack->GetStartVertex(vertex);
    Bool_t isPrim = kFALSE;
    if ( (vertex-fTargetPos).Mag() < 1. ) {
      isPrim = kTRUE;
      nPrim++;
    }
    else nSec++;

    // Get momentum
    Double_t mom = mcTrack->GetP();
    Bool_t isRef = kFALSE;
    if ( mom > 1. && isPrim) {
      isRef = kTRUE;
      nRef++;
    }

    // Fill histograms for accepted tracks
    fhMomAccAll->Fill(mom);
    fhNpAccAll->Fill(Double_t(nPoints));
    if ( isPrim) {
      fhMomAccPrim->Fill(mom);
      fhNpAccPrim->Fill(Double_t(nPoints));
    }
    else {
      fhMomAccSec->Fill(mom);
      fhNpAccSec->Fill(Double_t(nPoints));
      fhZAccSec->Fill(vertex.Z());
    }

    // Get matched MuchTrack
    Int_t    iRec  = -1;
    Double_t quali =  0.;
    Bool_t   isRec = kFALSE;
    if (fMatchMap.find(iMC) != fMatchMap.end() ) {
      iRec  = fMatchMap[iMC];
      isRec = kTRUE;
      CbmMuchTrack* MuchTrack = (CbmMuchTrack*) fMuchTracks->At(iRec);
      if ( ! MuchTrack ) {
         cout << "-E- " << GetName() << "::Exec: "
              << "No MuchTrack for matched MCTrack " << iMC << endl;
         Fatal("Exec", "No MuchTrack for matched MCTrack");
      }
      quali = fQualiMap[iMC];
      if ( quali < fQuota ) {
         cout << "-E- " << GetName() << "::Exec: "
              << "Matched MuchTrack " << iRec << " is below matching "
              << "criterion ( " << quali << ")" << endl;
         Fatal("Exec", "Match below matching quota");
      }
      CbmMuchTrackMatch* match = (CbmMuchTrackMatch*) fMuchMatches->At(iRec);
      if ( ! match ) {
         cout << "-E- " << GetName() << "::Exec: "
              << "No MuchTrackMatch for matched MCTrack " << iMC << endl;
         Fatal("Exec", "No MuchTrackMatch for matched MCTrack");
      }
      Int_t nTrue  = match->GetNofTrueHits();
      Int_t nWrong = match->GetNofWrongHits();
      Int_t nFake  = match->GetNofFakeHits();
      Int_t nHits  = MuchTrack->GetNHits();
   //   if ( nTrue + nWrong + nFake != nHits ) {
   //      cout << "True " << nTrue << " wrong " << nWrong << " Fake "
   //           << nFake << " Hits " << nHits << endl;
   //      Fatal("Exec", "Wrong number of hits");
   //   }

      // Verbose output
      if ( fVerbose > 1 )
         cout << "-I- " << GetName() << ": "
              << "MCTrack " << iMC << ", points "
              << nPoints << ", MuchTrack " << iRec << ", hits " << nHits
              << ", true hits " << nTrue << endl;
      
      recSignalCnt++;
    
      if (recSignalCnt == 2) {
         nRecSignalPair++;
      }
      
      
      // Fill histograms for reconstructed tracks
      nRecAll++;
      fhMomRecAll->Fill(mom);
      fhNpRecAll->Fill(Double_t(nPoints));
      if ( isPrim ) {
         nRecPrim++;
         fhMomRecPrim->Fill(mom);
         fhNpRecPrim->Fill(Double_t(nPoints));
         if ( isRef ) nRecRef++;
      }
      else {
         nRecSec++;
         fhMomRecSec->Fill(mom);
         fhNpRecSec->Fill(Double_t(nPoints));
         fhZRecSec->Fill(vertex.Z());
      }

    }  // Match found in map?

  } // Loop over MCTracks

  
  // Calculate efficiencies
  Double_t effAll;
  Double_t effPrim;
  Double_t effRef;
  Double_t effSec;
  Double_t effGhosts;
  Double_t effClones;
  
  if (nAcc > 0) effAll  = Double_t(nRecAll)  / Double_t(nAcc); else effAll = 0.0;
  if (nPrim > 0) effPrim = Double_t(nRecPrim) / Double_t(nPrim); else effPrim = 0.0;
  if (nRef > 0) effRef  = Double_t(nRecRef)  / Double_t(nRef); else effRef = 0.0;
  if (nSec > 0) effSec  = Double_t(nRecSec)  / Double_t(nSec); else effSec = 0.0;
  if (nAcc > 0) effGhosts = Double_t(nGhosts) / Double_t(nAcc); else effGhosts = 0.0;
  if (nAcc > 0) effClones = Double_t(nClones) / Double_t(nAcc); else effClones = 0.0;
  
  


  // Event summary
  if ( fVerbose > 0 ) {
    cout << "----------   MuchFindTracksQaA : Event summary   ------------"
         << endl;
    cout << "MCTracks   : " << nAll << ", accepted: " << nAcc
         << ", reconstructed: " << nRecAll << endl;
    cout << "All        : accepted: " << nAcc << ", reconstructed: "
         << nRecAll << ", efficiency " << effAll*100. << "%" << endl;
    cout << "Vertex     : accepted: " << nPrim << ", reconstructed: "
         << nRecPrim << ", efficiency " << effPrim*100. << "%" << endl;
    cout << "Reference  : accepted: " << nRef  << ", reconstructed: "
         << nRecRef  << ", efficiency " << effRef*100. << "%" << endl;
    cout << "Non-vertex : accepted: " << nSec << ", reconstructed: "
         << nRecSec << ", efficiency " << effSec*100. << "%" << endl;
    cout << "Ghosts : " << nGhosts << ", ghosts/accepted MC tracks: "
         << effGhosts*100. << "%" << endl;
    cout << "Clones : " << nClones << ", clones/accepted MC tracks: "
         << effClones*100. << "%" << endl;
    cout << "Number of MuchTracks " << nRec << ", ghosts " << nGhosts
         << ", clones " << nClones << endl;
    cout << "Number of MuchTracks - " << nRec << endl;
    cout << "Number of good Much Tracks - " << (nRec - nGhosts - nClones) << endl;
    if (fNormType == 2)
       cout << "Number of good Much Tracks w/o matched STS track - " << nNoSts << endl;
    
    cout << "Number of Acc Signal Pairs: " << nAccSignalPair << endl
         << "Number of Rec Signal Pairs: " << nRecSignalPair << endl;
    
    
    cout << "Number of rec signal muons PID from STS: " << nofRecMuonsSts << endl
         << "Number of Rec Signal Pairs PID from STS: " << nofRecMuonPairsSts << endl;    

    
   // cout << "Number of Muon tracks - " << nofMuons << endl;
   // cout << "Number of not Muon tracks - " << nofNoMuons << endl;
    
    cout << "-----------------------------------------------------------"
         << endl;
    cout << endl;
  }
  else cout << "-I- " << GetName() << ": all " << effAll*100. << " %, prim. "
      << effPrim*100. << " %, ref. " << effRef*100. << " %" << endl;

  // Increase counters
  fNNoSts += nNoSts;
  fNAll += nRec;
  fNAllGood += (nRec - nGhosts - nClones);
  fNAccAll  += nAcc;
  fNAccPrim += nPrim;
  fNAccRef  += nRef;
  fNAccSec  += nSec;
  fNRecAll  += nRecAll;
  fNRecPrim += nRecPrim;
  fNRecRef  += nRecRef;
  fNRecSec  += nRecSec;
  fNGhosts  += nGhosts;
  fNClones  += nClones;
  fNAccSignalPairs += nAccSignalPair;
  fNRecSignalPairs += nRecSignalPair;
  fNRecMuonsSts += nofRecMuonsSts;
  fNRecMuonPairsSts += nofRecMuonPairsSts;
  fNEvents++;

}
// -------------------------------------------------------------------------



// -----   Private method Finish   -----------------------------------------
void CbmLitMuchRecQa::Finish() {

  // Divide histograms for efficiency calculation
  DivideHistos(fhMomRecAll,  fhMomAccAll,  fhMomEffAll);
  DivideHistos(fhMomRecPrim, fhMomAccPrim, fhMomEffPrim);
  DivideHistos(fhMomRecSec,  fhMomAccSec,  fhMomEffSec);
  DivideHistos(fhNpRecAll,   fhNpAccAll,   fhNpEffAll);
  DivideHistos(fhNpRecPrim,  fhNpAccPrim,  fhNpEffPrim);
  DivideHistos(fhNpRecSec,   fhNpAccSec,   fhNpEffSec);
  DivideHistos(fhZRecSec,    fhZAccSec,    fhZEffSec);

  // Normalise histos for clones and ghosts to one event
  if ( fNEvents ) {
    fhNhClones->Scale(1./Double_t(fNEvents));
    fhNhGhosts->Scale(1./Double_t(fNEvents));
  }

  // Calculate integrated efficiencies and rates
  Double_t effAll = 0; 
  if (fNAccAll != 0) effAll = Double_t(fNRecAll)  / Double_t(fNAccAll);
  
  Double_t effPrim = 0;
  if (fNAccPrim != 0) effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim);
  
  Double_t effRef = 0;
  if (fNAccRef != 0) effRef = Double_t(fNRecRef)  / Double_t(fNAccRef);
  
  Double_t effSec = 0;
  if (fNAccSec != 0) effSec = Double_t(fNRecSec)  / Double_t(fNAccSec);
  
  Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNEvents);
  
  Double_t rateClones = Double_t(fNClones) / Double_t(fNEvents);
  
  Double_t effGhosts = 0;
  if (fNAccAll != 0) effGhosts = Double_t(fNGhosts) / Double_t(fNAccAll);
  
  Double_t effClones = 0;
  if (fNAccAll != 0) effClones = Double_t(fNClones) / Double_t(fNAccAll);
  
  Double_t rateRecAll  = Double_t(fNRecAll)  / Double_t(fNEvents);
  Double_t rateRecPrim = Double_t(fNRecPrim) / Double_t(fNEvents);
  Double_t rateRecRef  = Double_t(fNRecRef)  / Double_t(fNEvents);
  Double_t rateRecSec  = Double_t(fNRecSec)  / Double_t(fNEvents);
  Double_t rateAccAll  = Double_t(fNAccAll)  / Double_t(fNEvents);
  Double_t rateAccPrim = Double_t(fNAccPrim) / Double_t(fNEvents);
  Double_t rateAccRef  = Double_t(fNAccRef)  / Double_t(fNEvents);
  Double_t rateAccSec  = Double_t(fNAccSec)  / Double_t(fNEvents);
  
  // Run summary to screen
  cout << endl << endl;
  cout << "=======================================================" 
       << endl;
  cout << "            MuchFindTracksQaA: Run summary" << endl << endl;
  cout << "Events analysed : " << fNEvents << endl;
  cout << "Efficiency all tracks       : " << effAll*100 << " % (" 
       << fNRecAll << "/" << fNAccAll <<"), per event (" 
       << rateRecAll << "/" << rateAccAll << ")" << endl;
  cout << "Efficiency vertex tracks    : " << effPrim*100 << " % (" 
       << fNRecPrim << "/" << fNAccPrim <<"), per event (" 
       << rateRecPrim << "/" << rateAccPrim << ")" << endl;
  cout << "Efficiency reference tracks : " << effRef*100 << " % (" 
       << fNRecRef << "/" << fNAccRef <<"), per event (" 
       << rateRecRef << "/" << rateAccRef << ")" << endl;
  cout << "Efficiency secondary tracks : " << effSec*100 << " % (" 
       << fNRecSec << "/" << fNAccSec <<"), per event (" 
       << rateRecSec << "/" << rateAccSec << ")" << endl;
  cout << "Ghost rate " << rateGhosts << " per event, "
       << "ghosts/accepted MC tracks :" << effGhosts*100 << "% ("
       << fNGhosts << "/" << fNAccAll << ")" << endl;
  cout << "Clone rate " << rateClones << " per event, "
       << "clones/accepted MC tracks :" << effClones*100 << "% ("
       << fNClones << "/" << fNAccAll << ")" << endl;
  
  cout << "Number of Much tracks - " << fNAll << ", per event - " 
       << Double_t(fNAll)/ Double_t(fNEvents) << endl;
  cout << "Number of good Much Tracks - " << fNAllGood << ", per event - "
       << Double_t(fNAllGood)/ Double_t(fNEvents) << endl;
  if (fNormType == 2)
       cout << "Number of good Much Tracks w/o matched STS track - " << fNNoSts << ", per event - "
            << Double_t(fNNoSts)/ Double_t(fNEvents) << endl;     
  cout << "Ghosts - " << fNGhosts << ", per event - "
       << Double_t(fNGhosts)/ Double_t(fNEvents) << endl;
  cout << "Clones - " << fNClones << ", per event - "
       << Double_t(fNClones)/ Double_t(fNEvents) << endl;
  
  cout << "Number of Acc Signal Pairs: " << fNAccSignalPairs << endl
       << "Number of Rec Signal Pairs: " << fNRecSignalPairs << endl;     
  
  cout << "Number of rec signal muons PID from STS: " << fNRecMuonsSts << endl
         << "Number of Rec Signal Pairs PID from STS: " << fNRecMuonPairsSts << endl;   
  
  cout << "=======================================================" 
       << endl;
  cout << endl << endl;

  TDirectory *olddir = gDirectory;
  TDirectory *hdir = new TDirectory("MuchFinderQa",
				    "Performance of the track finding in Much");
  hdir->cd();
  TIter next(fHistoList);
  while ( TH1* histo = ((TH1*)next()) ) histo->Write();
  olddir->cd();

}
// -------------------------------------------------------------------------



// -----   Private method GetGeometry   ------------------------------------
InitStatus CbmLitMuchRecQa::GetGeometry() {

  // Get target geometry
  if ( ! fPassGeo ) {
    cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!"
         <<endl;
    fTargetPos.SetXYZ(0., 0., 0.);
    return kERROR;
  }
  TObjArray* passNodes = fPassGeo->GetGeoPassiveNodes();
  if ( ! passNodes ) {
    cout << "-W- " << GetName() << "::GetGeometry: No passive node array" 
         << endl;
    fTargetPos.SetXYZ(0., 0., 0.);
    return kERROR;
  }
  CbmGeoNode* target = (CbmGeoNode*) passNodes->FindObject("targ");
  if ( ! target ) {
    cout << "-E- " << GetName() << "::GetGeometry: No target node" 
         << endl;
    fTargetPos.SetXYZ(0., 0., 0.);
    return kERROR;
  }
  CbmGeoVector targetPos = target->getLabTransform()->getTranslation();
  CbmGeoVector centerPos = target->getCenterPosition().getTranslation();
  Double_t targetX = targetPos.X() + centerPos.X();
  Double_t targetY = targetPos.Y() + centerPos.Y();
  Double_t targetZ = targetPos.Z() + centerPos.Z();
  fTargetPos.SetXYZ(targetX, targetY, targetZ);
  
  return kSUCCESS;

}
// -------------------------------------------------------------------------



// -----   Private method CreateHistos   -----------------------------------
void CbmLitMuchRecQa::CreateHistos() {

  // Histogram list
  fHistoList = new TList();

  // Momentum distributions
  Double_t minMom   =  0.;
  Double_t maxMom   = 10.;
  Int_t    nBinsMom = 10;
  fhMomAccAll  = new TH1F("hMomAccAll", "all accepted tracks",
                           nBinsMom, minMom, maxMom);
  fhMomRecAll  = new TH1F("hMomRecAll", "all reconstructed tracks",
                           nBinsMom, minMom, maxMom);
  fhMomEffAll  = new TH1F("hMomEffAll", "efficiency all tracks",
                           nBinsMom, minMom, maxMom);
  fhMomAccPrim = new TH1F("hMomAccPrim", "accepted vertex tracks",
                           nBinsMom, minMom, maxMom);
  fhMomRecPrim = new TH1F("hMomRecPrim", "reconstructed vertex tracks",
                           nBinsMom, minMom, maxMom);
  fhMomEffPrim = new TH1F("hMomEffPrim", "efficiency vertex tracks",
                           nBinsMom, minMom, maxMom);
  fhMomAccSec  = new TH1F("hMomAccSec", "accepted non-vertex tracks",
                           nBinsMom, minMom, maxMom);
  fhMomRecSec  = new TH1F("hMomRecSec", "reconstructed non-vertex tracks",
                           nBinsMom, minMom, maxMom);
  fhMomEffSec  = new TH1F("hMomEffSec", "efficiency non-vertex tracks",
                           nBinsMom, minMom, maxMom);
  fHistoList->Add(fhMomAccAll);
  fHistoList->Add(fhMomRecAll);
  fHistoList->Add(fhMomEffAll);
  fHistoList->Add(fhMomAccPrim);
  fHistoList->Add(fhMomRecPrim);
  fHistoList->Add(fhMomEffPrim);
  fHistoList->Add(fhMomAccSec);
  fHistoList->Add(fhMomRecSec);
  fHistoList->Add(fhMomEffSec);

  // Number-of-points distributions
  Double_t minNp   = -0.5;
  Double_t maxNp   = 15.5;
  Int_t    nBinsNp = 16;
  fhNpAccAll  = new TH1F("hNpAccAll", "all accepted tracks",
                           nBinsNp, minNp, maxNp);
  fhNpRecAll  = new TH1F("hNpRecAll", "all reconstructed tracks",
                           nBinsNp, minNp, maxNp);
  fhNpEffAll  = new TH1F("hNpEffAll", "efficiency all tracks",
                           nBinsNp, minNp, maxNp);
  fhNpAccPrim = new TH1F("hNpAccPrim", "accepted vertex tracks",
                           nBinsNp, minNp, maxNp);
  fhNpRecPrim = new TH1F("hNpRecPrim", "reconstructed vertex tracks",
                           nBinsNp, minNp, maxNp);
  fhNpEffPrim = new TH1F("hNpEffPrim", "efficiency vertex tracks",
                           nBinsNp, minNp, maxNp);
  fhNpAccSec  = new TH1F("hNpAccSec", "accepted non-vertex tracks",
                           nBinsNp, minNp, maxNp);
  fhNpRecSec  = new TH1F("hNpRecSec", "reconstructed non-vertex tracks",
                           nBinsNp, minNp, maxNp);
  fhNpEffSec  = new TH1F("hNpEffSec", "efficiency non-vertex tracks",
                           nBinsNp, minNp, maxNp);
  fHistoList->Add(fhNpAccAll);
  fHistoList->Add(fhNpRecAll);
  fHistoList->Add(fhNpEffAll);
  fHistoList->Add(fhNpAccPrim);
  fHistoList->Add(fhNpRecPrim);
  fHistoList->Add(fhNpEffPrim);
  fHistoList->Add(fhNpAccSec);
  fHistoList->Add(fhNpRecSec);
  fHistoList->Add(fhNpEffSec);

  // z(vertex) distributions
  Double_t minZ    =  0.;
  Double_t maxZ    = 50.;
  Int_t    nBinsZ  = 50;
  fhZAccSec = new TH1F("hZAccSec", "accepted non-vertex tracks",
                        nBinsZ, minZ, maxZ);
  fhZRecSec = new TH1F("hZRecSecl", "reconstructed non-vertex tracks",
                        nBinsZ, minZ, maxZ);
  fhZEffSec = new TH1F("hZEffRec", "efficiency non-vertex tracks",
                        nBinsZ, minZ, maxZ);
  fHistoList->Add(fhZAccSec);
  fHistoList->Add(fhZRecSec);
  fHistoList->Add(fhZEffSec);

  // Number-of-hit distributions
  fhNhClones  = new TH1F("hNhClones", "number of hits for clones",
                           nBinsNp, minNp, maxNp);
  fhNhGhosts  = new TH1F("hNhGhosts", "number of hits for ghosts",
                           nBinsNp, minNp, maxNp);
  fHistoList->Add(fhNhClones);
  fHistoList->Add(fhNhGhosts);

}
// -------------------------------------------------------------------------



// -----   Private method Reset   ------------------------------------------
void CbmLitMuchRecQa::Reset() {

  TIter next(fHistoList);
  while ( TH1* histo = ((TH1*)next()) ) histo->Reset();

  fNAccAll = fNAccPrim = fNAccRef = fNAccSec = 0;
  fNRecAll = fNRecPrim = fNRecRef = fNRecSec = 0;
  fNGhosts = fNClones = fNEvents = 0;

}
// -------------------------------------------------------------------------




// -----   Private method DivideHistos   -----------------------------------
void CbmLitMuchRecQa::DivideHistos(TH1* histo1, TH1* histo2,
                                        TH1* histo3) {
  
  if ( !histo1 || !histo2 || !histo3 ) {
    cout << "-E- " << GetName() << "::DivideHistos: "
         << "NULL histogram pointer" << endl;
    Fatal("DivideHistos", "Null histo pointer");
  }

  Int_t nBins = histo1->GetNbinsX();
  if ( histo2->GetNbinsX() != nBins || histo3->GetNbinsX() != nBins ) {
    cout << "-E- " << GetName() << "::DivideHistos: "
         << "Different bin numbers in histos" << endl;
    cout << histo1->GetName() << " " << histo1->GetNbinsX() << endl;
    cout << histo2->GetName() << " " << histo2->GetNbinsX() << endl;
    cout << histo3->GetName() << " " << histo3->GetNbinsX() << endl;
   return;
  }

  Double_t c1, c2, c3, ce;
  for (Int_t iBin=0; iBin<nBins; iBin++) {
    c1 = histo1->GetBinContent(iBin);
    c2 = histo2->GetBinContent(iBin);
    if ( c2 ) {
      c3 = c1 / c2;
      ce = TMath::Sqrt( c3 * ( 1. - c3 ) / c2 );
    }
    else {
      c3 = 0.;
      ce = 0.;
    }
    histo3->SetBinContent(iBin, c3);
    histo3->SetBinError(iBin, ce);
  }

}
// -------------------------------------------------------------------------



ClassImp(CbmLitMuchRecQa)