//* $Id: */

// -------------------------------------------------------------------------
// -----                    CbmStsFindHitsQa source file               -----
// -----                  Created 01/07/2008  by R. Karabowicz         -----
// -------------------------------------------------------------------------

#include "CbmStsFindHitsQa.h"

#include "CbmGeoStsPar.h"
#include "CbmStsAddress.h"
#include "CbmStsDigi.h"
#include "CbmStsCluster.h"
#include "CbmStsDigiPar.h"
#include "CbmStsDigiScheme.h"
#include "CbmStsPoint.h"
#include "CbmStsHit.h"
#include "legacy/CbmStsSensor_old.h"
#include "CbmStsSector.h"
#include "CbmStsStation_old.h"
#include "CbmMCTrack.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TClonesArray.h"
#include "TMath.h"
#include "TGeoManager.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TLegend.h"
#include "TPaveText.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TList.h"
#include "TVector3.h"

#include <iomanip>

using std::cout;
using std::cerr;
using std::endl;
using std::flush;
using std::fixed;
using std::right;
using std::left;
using std::setw;
using std::setprecision;
using std::set;
using std::map;

// -----   Default constructor   ------------------------------------------
CbmStsFindHitsQa::CbmStsFindHitsQa() 
: FairTask("STS Hit Finder", 1),
  fGeoPar(NULL),
  fDigiPar(NULL),
  fDigiScheme(NULL),
  fMCTracks(NULL),          // MCtrack
  fStsPoints(NULL),        // StsPoints
  fStsHits(NULL),          // StsHits
  fMatches(NULL),          // StsTrackMatch
  fStsDigis(NULL),          // StsDigi
  fStsClusters(NULL),          // StsCluster
  fTimer(),
  fhHitPointCorrelation(),
  fhHitPointCorrelationU(),
  fhHitPointCorrelationB(),
  fhEnergyLoss(),
  fhIncAngle(),
  fhPoints(),
  fhRecoPoints(),
  fhHitFindingEfficiency(NULL),
  fhEffIncAng(NULL),
  fhEffMom(NULL),
  fhEffPdgSec(NULL),
  fhEffPdgPrim(NULL),
  fhHitPointPull(NULL),
  fhHitPointCorr(NULL),
  fhNofHits(),
  fHistoList(NULL),
  fHistoListPS(),
  fNStations(0),
  fNEvents(0),
  fTime1(0.),
  fNofHits(),
  fNofDigisPChip(),
  fNofPoints(),
  fNofRecoPoints(),
  fNofPointsIncAng(),
  fNofRecoPointsIncAng(),
  fNofPointsMom(),
  fNofRecoPointsMom(),
  fNofRecoPdgSec(),
  fNofPointsPdgSec(),
  fNofRecoPdgPrim(),
  fNofPointsPdgPrim(),
  fNofPointsPrim(),
  fNofPointsSec(),
  fNofRecoPrim(),
  fNofRecoSec(),
  fNofPointsMomSum(),
  fNofRecoPointsMomSum(),
  fOnlineAnalysis(kFALSE),
  recoCanvas(NULL),
  recoPad()
{  
  fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Standard constructor   ------------------------------------------
CbmStsFindHitsQa::CbmStsFindHitsQa(Bool_t visualizeBool, Int_t iVerbose)
  : FairTask("STSFindHitsQa", iVerbose), 
    fGeoPar(NULL),
    fDigiPar(NULL),
    fDigiScheme(NULL),
    fMCTracks(NULL),          // MCtrack
    fStsPoints(NULL),          // StsPoints
    fStsHits(NULL),          // StsHits
    fMatches(NULL),          // StsTrackMatch
    fStsDigis(NULL),          // StsDigi
    fStsClusters(NULL),          // StsCluster
    fTimer(),
    fhHitPointCorrelation(),
    fhHitPointCorrelationU(),
    fhHitPointCorrelationB(),
    fhEnergyLoss(),
    fhIncAngle(),
    fhPoints(),
    fhRecoPoints(),
    fhHitFindingEfficiency(NULL),
    fhEffIncAng(NULL),
    fhEffMom(NULL),
    fhEffPdgSec(NULL),
    fhEffPdgPrim(NULL),
    fhHitPointPull(NULL),
    fhHitPointCorr(NULL),
    fhNofHits(),
    fHistoList(NULL),
    fHistoListPS(),
    fNStations(0),
    fNEvents(0),
    fTime1(0.),
    fNofHits(),
    fNofDigisPChip(),
    fNofPoints(),
    fNofRecoPoints(),
    fNofPointsIncAng(),
    fNofRecoPointsIncAng(),
    fNofPointsMom(),
    fNofRecoPointsMom(),
    fNofRecoPdgSec(),
    fNofPointsPdgSec(),
    fNofRecoPdgPrim(),
    fNofPointsPdgPrim(),
    fNofPointsPrim(),
    fNofPointsSec(),
    fNofRecoPrim(),
    fNofRecoSec(),
    fNofPointsMomSum(),
    fNofRecoPointsMomSum(),
    fOnlineAnalysis(visualizeBool),
    recoCanvas(NULL),
    recoPad()
{  
  fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Constructor with name   -----------------------------------------
CbmStsFindHitsQa::CbmStsFindHitsQa(const char* name, Int_t iVerbose) 
  : FairTask(name, iVerbose), 
    fGeoPar(NULL),
    fDigiPar(NULL),
    fDigiScheme(NULL),
    fMCTracks(NULL),          // MCtrack
    fStsPoints(NULL),          // StsPoints
    fStsHits(NULL),          // StsHits
    fMatches(NULL),          // StsTrackMatch
    fStsDigis(NULL),          // StsDigi
    fStsClusters(NULL),          // StsCluster
    fTimer(),
    fhHitPointCorrelation(),
    fhHitPointCorrelationU(),
    fhHitPointCorrelationB(),
    fhEnergyLoss(),
    fhIncAngle(),
    fhPoints(),
    fhRecoPoints(),
    fhHitFindingEfficiency(NULL),
    fhEffIncAng(NULL),
    fhEffMom(NULL),
    fhEffPdgSec(NULL),
    fhEffPdgPrim(NULL),
    fhHitPointPull(NULL),
    fhHitPointCorr(NULL),
    fhNofHits(),
    fHistoList(NULL),
    fHistoListPS(),
    fNStations(0),
    fNEvents(0),
    fTime1(0.),
    fNofHits(),
    fNofDigisPChip(),
    fNofPoints(),
    fNofRecoPoints(),
    fNofPointsIncAng(),
    fNofRecoPointsIncAng(),
    fNofPointsMom(),
    fNofRecoPointsMom(),
    fNofRecoPdgSec(),
    fNofPointsPdgSec(),
    fNofRecoPdgPrim(),
    fNofPointsPdgPrim(),
    fNofPointsPrim(),
    fNofPointsSec(),
    fNofRecoPrim(),
    fNofRecoSec(),
    fNofPointsMomSum(),
    fNofRecoPointsMomSum(),
    fOnlineAnalysis(kFALSE),
    recoCanvas(NULL),
    recoPad()
{
  fDigiScheme = new CbmStsDigiScheme();
}
// -------------------------------------------------------------------------



// -----   Destructor   ----------------------------------------------------
CbmStsFindHitsQa::~CbmStsFindHitsQa() {
  fHistoList->Delete();
  delete fHistoList;
  for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
    fHistoListPS[istat]->Delete();
  }
  delete[] fHistoList;
}
// -------------------------------------------------------------------------



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

  fTimer.Start();
  Bool_t warn = kFALSE;

  // Check for digi scheme
  if ( ! fDigiScheme ) {
    cerr << "-E- " << fName << "::Exec: No digi scheme!" << endl;
    return;
  }

  CbmStsStation_old* station = NULL;
  CbmStsSector*  sector  = NULL;
  
  Int_t nofStsDigis = fStsDigis->GetEntriesFast();
  for ( Int_t idigi = 0 ; idigi < nofStsDigis ; idigi++ ) {
    CbmStsDigi *stsDigi     = (CbmStsDigi*)fStsDigis->At(idigi);
    //  cout << "digi side = " << stsDigi->GetSide() << endl;
//    fNofFiredDigis[stsDigi->GetStationNr()-1][stsDigi->GetSectorNr()-1][stsDigi->GetSide()] += 1;
    Int_t iStation = CbmStsAddress::GetElementId(stsDigi->GetAddress(), kStsStation);
    //Int_t iSector = CbmStsAddress::GetElementId(stsDigi->GetAddress(), kStsModule);
    Int_t iSector = stsDigi->GetSectorNr();
    Int_t iSide = CbmStsAddress::GetElementId(stsDigi->GetAddress(), kStsSide);
    Int_t iChannel = CbmStsAddress::GetElementId(stsDigi->GetAddress(), kStsChannel);
    fNofDigisPChip[iStation-1][iSector-1][iSide][(Int_t)(iChannel/128)] += 1;
  }
  Int_t nofStsHits = fStsHits->GetEntriesFast();
  Int_t nofStsPoints = fStsPoints->GetEntriesFast();
  Int_t hitStationLimits[2][100];

  for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
    hitStationLimits[0][ist] = -1;
    hitStationLimits[1][ist] = -1;
  }

  // check for limits of hit indices on different stations...
  
  for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
    CbmStsHit *stsHit     = (CbmStsHit*)fStsHits->At(ihit);
    fNofHits[CbmStsAddress::GetElementId(stsHit->GetAddress(), kStsStation)][stsHit->GetSectorNr()-1][CbmStsAddress::GetElementId(stsHit->GetAddress(), kStsSensor)] += 1; // count hits per sensor
    if ( hitStationLimits[0][CbmStsAddress::GetElementId(stsHit->GetAddress(), kStsStation)] == -1 ) {
      hitStationLimits[0][CbmStsAddress::GetElementId(stsHit->GetAddress(), kStsStation)] = ihit;
    }
    CbmStsHit *stsHitBack = (CbmStsHit*)fStsHits->At(nofStsHits-ihit-1);
    if ( hitStationLimits[1][CbmStsAddress::GetElementId(stsHitBack->GetAddress(), kStsStation)] == -1 ) {
      hitStationLimits[1][CbmStsAddress::GetElementId(stsHitBack->GetAddress(), kStsStation)] = nofStsHits-ihit;
    }
  }
  for ( Int_t istat = 0 ; istat < fDigiScheme->GetNStations() ; istat++ ) {
    station = fDigiScheme->GetStation(istat);
    for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
      sector = station->GetSector(isect);

      CbmStsSensor_old* sensor = (CbmStsSensor_old*)sector->GetSensor(0);
      for ( Int_t iside = 0 ; iside<2; iside++){
        for (Int_t iChip = 0; iChip<8; iChip++){
          if (fNEvents==999) cout<<" St "<<istat+1<<" sect "<<isect+1<<" nofsens "<<sector->GetNSensors()<<" sectX "<<sensor->GetX0()<<" sectY "<<sensor->GetY0()<<" side "<<iside<<" chip "<<iChip+1<<" fNofDigisPChip "<<fNofDigisPChip[istat][isect][iside][iChip]<<endl;
        }
      }

      for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
	//	cout << "filling   " << istat << " " << isect << " " << isens << endl;
	fhNofHits[istat][isect][isens] -> Fill(fNofHits[istat][isect][isens]);
	//	cout << "---> done " << isens << endl;
      
      }
//       for ( Int_t iside = 0 ; iside < 2 ; iside++ ) {
// 	fhNofFiredDigis[istat][isect][iside]->Fill(fNofFiredDigis[istat][isect][iside]);
// 	for ( Int_t ichip = 0 ; ichip < 8 ; ichip++ ) {
// 	  if ( fhNofDigisPChip[istat][isect][iside][ichip] )
// 	    fhNofDigisPChip[istat][isect][iside][ichip]->Fill(fNofDigisPChip[istat][isect][iside][ichip]);
// 	}
//       }

    }
  }
  for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
    CbmStsPoint *stsPoint = (CbmStsPoint*)fStsPoints->At(ipnt);

    //    cout << "point " << ipnt << flush;

    Int_t trackPdg  = 0;
    Int_t motherId  = 0;
    Int_t motherPdg = 0;
    TVector3 trackVertex;
    TVector3 motherVertex;
    TVector3 momentum;

    Double_t xin = stsPoint->GetXIn();
    Double_t yin = stsPoint->GetYIn();
    Double_t zin = stsPoint->GetZIn();
    Double_t xout = stsPoint->GetXOut();
    Double_t yout = stsPoint->GetYOut();
    Double_t zout = stsPoint->GetZOut();

    gGeoManager->FindNode(xin,yin,zin);
    TGeoNode* curNode = gGeoManager->GetCurrentNode();
    
    Int_t trackId = stsPoint->GetTrackID();
    CbmMCTrack *mcTrack = (CbmMCTrack*) fMCTracks->At(trackId);
    trackPdg = mcTrack -> GetPdgCode();
    motherId = mcTrack -> GetMotherId();
    mcTrack->GetMomentum(momentum);
    Double_t mom = momentum.Mag();

    if (motherId>=0 ) {
      CbmMCTrack *motherTrack = (CbmMCTrack*) fMCTracks->At(motherId);
      motherPdg = motherTrack->GetPdgCode();
      mcTrack->GetStartVertex(trackVertex);
      motherTrack->GetStartVertex(motherVertex);
      
    }

    //    cout << " in node \"" << curNode->GetName() << "\"" << flush;
    TString curPath = fDigiScheme->GetCurrentPath();
    CbmStsSensor_old* sensor = fDigiScheme->GetSensorByName(curPath);
    Int_t stationNr = sensor->GetStationNr();

    Int_t startHit = hitStationLimits[0][stationNr-1];
    Int_t finalHit = hitStationLimits[1][stationNr-1];

    //    cout << " at station " << stationNr << " ( " << startHit << " - " << finalHit << " )" << endl;
    
    fhEnergyLoss[stationNr-1]->Fill(stsPoint->GetXIn(),stsPoint->GetYIn(),stsPoint->GetEnergyLoss());
 
    fhIncAngle[stationNr-1]->Fill(stsPoint->GetXIn(),stsPoint->GetYIn(),TMath::Abs(TMath::ATan2((stsPoint->GetZOut()-stsPoint->GetZIn()),(stsPoint->GetXOut()-stsPoint->GetXIn()))*180./3.14));
//     Float_t zP = stsPoint->GetZ();
//     Float_t xP = stsPoint->GetX(zP);
//     Float_t yP = stsPoint->GetY(zP);

    Bool_t foundPnt = kFALSE;

    if ( startHit == -1 && finalHit == -1 ) continue;
    
    for ( Int_t ihit = startHit ; ihit < finalHit ; ihit++ ) {
      CbmStsHit *stsHit= (CbmStsHit*)fStsHits->At(ihit);
      if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .01 ) &&
 	   ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .04 ) ) {
 	fhHitPointCorrelation[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
						 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
	fhHitPointCorr->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
						 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
	fhHitPointPull -> Fill( (stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()))/stsHit->GetDx() , (stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()))/stsHit->GetDy());
	
	if ( stsPoint->GetY(stsHit->GetZ()) > 0. )
	  fhHitPointCorrelationU[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
						    stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
	if ( stsPoint->GetY(stsHit->GetZ()) < 0. )
	  fhHitPointCorrelationB[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
						    stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
      }
      
      if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .01 ) &&
	   ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .04 ) ) {
	foundPnt = kTRUE;
      }
//       cout <<"St "<<stationNr<<" Sens "<<sensor->GetSensorNr()<< " chan "<<sensor->GetFrontChannel(stsHit->GetX(),stsHit->GetY(),stsHit->GetZ())<<" X "<<stsHit->GetX();
    }
    Int_t incAngle = (Int_t) (TMath::ATan2((stsPoint->GetZOut()-stsPoint->GetZIn()),(stsPoint->GetXOut()-stsPoint->GetXIn()))*180./3.14);
//    Int_t incAngle = (Int_t) ( TMath::ATan((stsPoint->GetXOut()-stsPoint->GetXIn())/(stsPoint->GetZOut()-stsPoint->GetZIn()))*180./3.14);
    fNofPoints    [stationNr-1][sensor->GetSectorNr()-1][sensor->GetSensorNr()-1] += 1;
    fNofPointsIncAng[TMath::Abs((Int_t)(incAngle))] +=1;
    if (TMath::Abs((Int_t)(mom))<1) fNofPointsMom [TMath::Abs((Int_t)(mom))] +=1;
    if (TMath::Abs((Int_t)(mom))<1) fNofPointsMomSum+=1;
    //  cout << stationNr-1 << " " << sensor->GetSectorNr()-1 << " " << sensor->GetSensorNr()-1 << "/" << flush;
    fhPoints      [stationNr-1]->Fill(stsPoint->GetX(stsPoint->GetZ()),
				      stsPoint->GetY(stsPoint->GetZ()));
    fNofPointsSec = 1;
    if (trackPdg<10000&&motherId>=0) {
      fNofPointsPdgSec[trackPdg]  += 1;
      fNofPointsSec +=1;
    }
    fNofPointsPrim = 1;
    if (trackPdg<10000&&motherId<0) {
      fNofPointsPdgPrim[trackPdg]  += 1;
      fNofPointsPrim +=1;
    }
    if ( foundPnt ) {
      fNofRecoPoints[stationNr-1][sensor->GetSectorNr()-1][sensor->GetSensorNr()-1] += 1;
      fNofRecoPointsIncAng[TMath::Abs((Int_t)(incAngle))] +=1;
      if (TMath::Abs((Int_t)(mom))<1) fNofRecoPointsMom [TMath::Abs((Int_t)(mom))] +=1;
      if (TMath::Abs((Int_t)(mom))<1) fNofRecoPointsMomSum+=1;
      fhRecoPoints  [stationNr-1]->Fill(stsPoint->GetX(stsPoint->GetZ()),
					stsPoint->GetY(stsPoint->GetZ()));
      fNofRecoPrim = 1;
      if (trackPdg<10000&&motherId>=0) {
	fNofRecoPdgSec[trackPdg] += 1;
	fNofRecoSec +=1;
      }
      fNofRecoPrim = 1;
      if (trackPdg<10000&&motherId<0) {
	fNofRecoPdgPrim[trackPdg] += 1;
	fNofRecoPrim +=1;
      }
    }
    
  }
  for (Int_t ang=0; ang<400; ang++) {
    Double_t factor = (Float_t)(fNofRecoPointsIncAng[ang])/(Float_t)(fNofPointsIncAng[ang]);
    fhEffIncAng->Fill(ang,((Float_t)(fNofRecoPointsIncAng[ang]))/((Float_t)(fNofPointsIncAng[ang])));
  }
    for (Int_t momentum=0; momentum<50; momentum++) {
    fhEffMom->Fill(momentum,((Float_t)(fNofRecoPointsMom[momentum]))/((Float_t)(fNofPointsMom[momentum])));
  }
  
  if ( fOnlineAnalysis ) {
    recoPad[0]->cd();
    TLegend* brp = new TLegend(0.1,0.1,0.9,0.9,"STS hit finding");
    brp->SetTextAlign(22);
    brp->SetTextSize(0.6);
    brp->SetTextColor(1);
    brp->SetBorderSize(0);
    brp->SetFillColor(0);
    brp->Draw();
    recoPad[0]->Update();

    recoPad[1]->cd();
    fhEffIncAng->SetLineWidth(2);
    fhEffIncAng->SetLineColor(1);
    fhEffIncAng->SetMarkerColor(kGreen);
    fhEffIncAng->SetMarkerStyle(2);
    fhEffIncAng->SetTitle("Efficiency");
    fhEffIncAng->Draw();
    TLine* oneLine = new TLine(0.0,1.0,10.0,1.0);
    oneLine->SetLineStyle(2);
    oneLine->Draw();
    recoPad[1]->Update();
    
    TH1F* projX[10];
    TH1F* projY[10];
    TF1*   fitX[10];
    TF1*   fitY[10];
    Double_t resolution[2][10];

    for ( Int_t ist = 0 ; ist < 8 ; ist++ ) {
      projX[ist] = (TH1F*)fhHitPointCorrelation[ist]->ProjectionX(Form("projX%i",ist+1));
      projY[ist] = (TH1F*)fhHitPointCorrelation[ist]->ProjectionY(Form("projY%i",ist+1));
      fitX[ist] = new TF1(Form("fitX%i",ist+1),"gaus",-0.02,0.02);
      fitY[ist] = new TF1(Form("fitY%i",ist+1),"gaus",-0.02,0.02);
      projX[ist]->SetAxisRange(-0.02,0.02,"X");
      projY[ist]->SetAxisRange(-0.02,0.02,"X");
      projX[ist]->Fit(fitX[ist],"QN","",-0.02,0.02);
      projY[ist]->Fit(fitY[ist],"QN","",-0.02,0.02);
      resolution[0][ist] = 10000.*fitX[ist]->GetParameter(2);
      resolution[1][ist] = 10000.*fitY[ist]->GetParameter(2);
    }
    
    TH1F* projXPull;
    TH1F* projYPull;
    TF1*   fitXPull;
    TF1*   fitYPull;
    Double_t resolutionPull[2];


      projXPull = (TH1F*)fhHitPointPull->ProjectionX("projXPull");
      projYPull = (TH1F*)fhHitPointPull->ProjectionY("projYPull");
      fitXPull = new TF1("fitXPull","gaus",-2,2);
      fitYPull = new TF1("fitYPull","gaus",-2,2);
      projXPull->SetAxisRange(-2,2,"X");
      projYPull->SetAxisRange(-2,2,"X");
      projXPull->Fit(fitXPull,"QN","",-2,2);
      projYPull->Fit(fitYPull,"QN","",-2,2);
      resolutionPull[0] = fitXPull->GetParameter(2);
      resolutionPull[1] = fitYPull->GetParameter(2);

    
    recoPad[3]->cd();
    fhHitPointCorrelation[0]->Draw("col");
    recoPad[3]->Update();
    
    recoPad[4]->cd();
    fhHitPointPull->Draw("col");
    recoPad[4]->Update();
    
    recoPad[5]->cd();
    projX[0]->SetLineWidth(2);
    projY[0]->SetLineWidth(2);
    projX[0]->SetLineColor(2);
    projY[0]->SetLineColor(3);
    projY[0]->SetXTitle("#Delta x, y [cm]");
    projX[0]->Draw();
    projY[0]->Draw("same");
    //  fitY[fShowStation1]->Draw("same");
    //  fitX[fShowStation1]->Draw("same");
    TLegend* legend1 = new TLegend(0.55,0.6,1.0,0.94);
    legend1->SetBorderSize(0);
    legend1->SetFillColor(0);
    legend1->AddEntry(projX[0],
		    Form("X,#sigma=%3.2f#mum",resolution[0][0]),"l");
    legend1->AddEntry(projY[0],
		    Form("Y,#sigma=%3.2f#mum",resolution[1][0]),"l");
    legend1->Draw();
    recoPad[5]->Update();
    
    recoPad[6]->cd();
    projXPull->SetLineWidth(2);
    projYPull->SetLineWidth(2);
    projXPull->SetLineColor(2);
    projYPull->SetLineColor(3);
    projXPull->SetXTitle("Pull x, y ");
    projXPull->Draw();
    projYPull->Draw("same");
    //  fitY[fShowStation2]->Draw("same");
    //  fitX[fShowStation2]->Draw("same");
    TLegend* legend2a = new TLegend(0.55,0.6,1.0,0.94);
    legend2a->SetBorderSize(0);
    legend2a->SetFillColor(0);
    legend2a->AddEntry(projXPull,Form("X,#sigma=%3.2f#mum",resolutionPull[0]),"l");
    legend2a->AddEntry(projYPull,Form("Y,#sigma=%3.2f#mum",resolutionPull[1]),"l");
    legend2a->Draw();
    recoPad[6]->Update();
    recoPad[2]->cd();
    fhRecoPoints[0]->SetMarkerColor(kCyan);
    fhRecoPoints[0]->Draw();
    recoPad[2]->Update();
    
    recoPad[7]->cd();
    fhRecoPoints[1]->SetMarkerColor(kMagenta);
    fhRecoPoints[1]->Draw();
    recoPad[7]->Update();
    
    recoPad[8]->cd();
    fhRecoPoints[2]->SetMarkerColor(kOrange);
    fhRecoPoints[2]->Draw();
    recoPad[8]->Update();
    
    recoPad[9]->cd();
    fhRecoPoints[4]->SetMarkerColor(kGreen);
    fhRecoPoints[4]->Draw();
    recoPad[9]->Update();
    
    recoPad[10]->cd();
    fhRecoPoints[5]->SetMarkerColor(kBlue);
    fhRecoPoints[5]->Draw();
    recoPad[10]->Update();
    
    recoPad[11]->cd();
    fhRecoPoints[7]->SetMarkerColor(kRed);
    fhRecoPoints[7]->Draw();
    recoPad[11]->Update();
    
    recoPad[12]->cd();
    TPaveText* printoutPave = new TPaveText(0.0,0.0,1.0,1.0);
    printoutPave->SetTextAlign(20);
    printoutPave->SetTextSize(0.08);
    printoutPave->SetTextColor(1);
    printoutPave->SetBorderSize(0);
    printoutPave->SetFillColor(0);
    printoutPave->AddText(Form("%i events",fNEvents+1));
    printoutPave->AddText(Form("Hits/Points %3.2f ",Double_t (nofStsHits)/Double_t (nofStsPoints)));

    printoutPave->AddText(Form("%3.0f points prim, %3.0f reco prim",
			     Double_t (fNofPointsPrim)/Double_t (fNEvents+1),
			     Double_t (fNofRecoPrim)/Double_t (fNEvents+1)));
    printoutPave->AddText(Form("%3.0f points sec, %3.0f reco sec",
			     Double_t (fNofPointsSec)/Double_t (fNEvents+1),
			     Double_t (fNofRecoSec)/Double_t (fNEvents+1)));
   
    printoutPave->AddText(Form("Hit Finding Efficiency: %3.2f  prim, %3.2f sec",
			     (Double_t (fNofRecoPrim)/Double_t (fNofPointsPrim))*100.,
			     (Double_t (fNofRecoSec) /Double_t (fNofPointsSec))*100.));
//     printoutPave->AddText("Tracking efficiencies (p>1.0 GeV/c):");
//     printoutPave->AddText(Form("all = %2.2f%%(%2.2f%%)",100.*effAll,100.*allEff));
//     printoutPave->AddText(Form("vertex = %2.2f%%(%2.2f%%)",100.*effPrim,100.*primEff));
//     printoutPave->AddText(Form("reference = %2.2f%%",100.*effRef));
//     printoutPave->AddText(Form("non-vertex = %2.2f%%(%2.2f%%)",100.*effSec,100.*secEff));
//     printoutPave->AddText(Form("Momentum resolution = %3.2f%%(%3.2f%%)",momentumResolutionAll,momentumResolutionPrim));
    recoPad[12]->Clear();
    printoutPave->Draw();
    recoPad[12]->Update();
    
  }

  // Event summary
  if ( fVerbose ) {
    cout << "----------   StsFindHitsQa : Event " << fNEvents+1 << " summary   ------------"
	 << endl;
    cout << "-----------------------------------------------------------" 
	 << endl;
    cout << endl;
  }
  else
    cout << "\r+ " << setw(15) << left << fName << ": event " << fNEvents+1 << "  " << setprecision(4) 
	 << setw(8) << fixed << right << fTimer.RealTime() << "s." << endl; 

  //cout << "--------        " << fNofRecoPoints[1][4][0] << " / " << fNofPoints[1][4][0] << endl;

  fNEvents++;
  fTime1 += fTimer.RealTime();
}
// -------------------------------------------------------------------------

    


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


  // Histogram list
  fHistoList = new TList();
  for ( Int_t ist = 0 ; ist < fNStations ; ist++ )
    fHistoListPS[ist] = new TList();

  // Momentum distributions
  Double_t minMom   =  0.;
  Double_t maxMom   = 10.;
  Int_t    nBinsMom = 40;

  CbmStsStation_old* station = NULL;
  CbmStsSector*  sector  = NULL;
  CbmStsSensor_old*  sensor  = NULL;

  fhHitFindingEfficiency = new TH1F("hHitFindingEfficiency","Hit finding efficiency",1000000,
				    (Float_t)(2<<24)-0.5,(Float_t)(2<<24)-0.5+1000000.);
  fhHitPointCorr = new TH2F(Form("hHitPointCorrelation"),
					  Form("Hit vs point correlation "),
					  500,-.1, .1,500,-.1,.1);
  fhEffMom     = new TH2F("hEffMom","hEffMom",500,0,25,100,0,1);
  fhEffIncAng  = new TH2F("hEffIncAng","hEffIncAng",100,0,100,100,0,1);
  fhEffPdgSec  = new TH2F("hEffPdgSec", "hEffPdgSec", 350,-1000,2500,10,0,1);
  fhEffPdgPrim = new TH2F("hEffPdgPrim","hEffPdgPrim",350,-1000,2500,10,0,1);
  fhHitPointPull = new TH2F("hHitPointPull","Hit vs sigma at station %i",500,-3., 3.,500,-3.,3.);

  fHistoList->Add(fhEffPdgSec);
  fHistoList->Add(fhEffPdgPrim);
  fHistoList->Add(fhHitFindingEfficiency);
  fHistoList->Add(fhEffIncAng);
  fHistoList->Add(fhHitPointPull);
  fHistoList->Add(fhEffMom);
  fHistoList->Add(fhHitPointCorr);
  for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
    fhHitPointCorrelation[ist] = new TH2F(Form("hHitPointCorrelation%i",ist+1),
					  Form("Hit vs point correlation at station %i",ist+1),
					  500,-.1, .1,500,-.1,.1);
    fhHitPointCorrelation[ist]->SetXTitle("#Delta x [cm]");  
    fhHitPointCorrelation[ist]->SetYTitle("#Delta y [cm]");
    

    
    
    fHistoListPS[ist]->Add(fhHitPointCorrelation[ist]);
    

    fhHitPointCorrelationU[ist] = new TH2F(Form("hHitPointCorrelationU%i",ist+1),
					   Form("Hit vs point correlation at station %i",ist+1),
					   500,-.1, .1,500,-.1,.1);
    fhHitPointCorrelationB[ist] = new TH2F(Form("hHitPointCorrelationB%i",ist+1),
					   Form("Hit vs point correlation at station %i",ist+1),
					   500,-.1, .1,500,-.1,.1);
    fHistoListPS[ist]->Add(fhHitPointCorrelationU[ist]);
    fHistoListPS[ist]->Add(fhHitPointCorrelationB[ist]);

    fhPoints    [ist] = new TH2F(Form("hPoints%i",ist+1),Form("Points on station %i",ist+1),
				 100,-100.,100.,
				 100,-100.,100.);
    fhRecoPoints[ist] = new TH2F(Form("hRecoPoints%i",ist+1),Form("Reconstructed points on station %i",ist+1),
				 100,-100.,100.,
				 100,-100.,100.);
    fHistoListPS[ist]->Add(fhPoints    [ist]);
    fHistoListPS[ist]->Add(fhRecoPoints[ist]);

  }

  Double_t binningNofHits[501];
  binningNofHits[0] = -0.5;
  for ( Int_t itemp = 0 ; itemp < 100 ; itemp++ ) {
    binningNofHits[itemp+  1] = 0.5+(Double_t)itemp;
    binningNofHits[itemp+101] = 101.5+2.*(Double_t)itemp;
    binningNofHits[itemp+201] = 302.5+3.*(Double_t)itemp;
    binningNofHits[itemp+301] = 602.5+3.*(Double_t)itemp;
    binningNofHits[itemp+401] = 902.5+3.*(Double_t)itemp;
  }

  Double_t fSectorWidth = 0.;
  for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
    station = fDigiScheme->GetStation(istat);
   
    fhEnergyLoss[istat] = new TH3F(Form("hEnergyLossSt%d",istat+1),
				   Form("Energy loss on station %d",istat+1),
				   20,-100.,100.,20,-100.,100.,100,0,0.01);
    fhIncAngle [istat] = new TH3F(Form("hIncAngleSt%d",istat+1),
				   Form("Incident Angle on station %d;X;Y;Angle",istat+1),
				   100,-50,50,100,-40,40,20,0,100.);
    fHistoListPS[istat]->Add(fhEnergyLoss[istat]);
    fHistoListPS[istat]->Add(fhIncAngle[istat]);
    for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
      sector = station->GetSector(isect);
      for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
	sensor = sector->GetSensor(isens);
	fSectorWidth = 10.*sensor->GetLx();
	fhNofHits[istat][isect][isens] = new TH1F(Form("hNofHitsSt%dSect%dSens%d",istat+1,isect+1,isens+1),
						  Form("Number of hits in sector %d, sensor %d of station %d",isect+1,isens+1,istat+1),
						  500,binningNofHits);
	//					 500,-0.5,500.5);
	fHistoListPS[istat]->Add(fhNofHits[istat][isect][isens]); 
      }
    }
  }

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



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


  TIter next(fHistoList);
  while ( TH1* histo = ((TH1*)next()) ) histo->Reset();
  for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
    TIter next0(fHistoListPS[ist]);
    while ( TH1* histo = ((TH1*)next0()) ) histo->Reset();
  }

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


// -----   Private method SetParContainers   -------------------------------
void CbmStsFindHitsQa::SetParContainers() {


  // Get run and runtime database
  FairRunAna* run = FairRunAna::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

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

  // Get STS geometry parameter container
  fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");

  // Get STS digitisation parameter container
  fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");

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




// -----   Private method Init   -------------------------------------------
InitStatus CbmStsFindHitsQa::Init() {


  // Get input array
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) Fatal("Init", "No FairRootManager");

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

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

  // Get StsClusters array
  fStsClusters = (TClonesArray*) ioman->GetObject("StsCluster");
  if ( ! fStsClusters ) {
    cout << "-w- " << GetName() << "::Init: No StsCluster array!" 
	 << endl;
    //  return kERROR;
  }

  // Build digitisation scheme
  Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
  if ( ! success ) return kERROR;

  fNStations = fDigiScheme->GetNStations();

  CbmStsStation_old* station = NULL;
  CbmStsSector*  sector  = NULL;
  for ( Int_t istat = fNStations ; istat > 0 ; ) {
    station = fDigiScheme->GetStation(--istat);
    for ( Int_t isect = station->GetNSectors() ; isect > 0 ; ) {
      sector = station->GetSector(--isect);
      for ( Int_t isens = sector->GetNSensors() ; isens > 0 ; ) {
	fNofPoints    [istat][isect][--isens] = 0;
	fNofRecoPoints[istat][isect][  isens] = 0;
      }
    }
  }
  for (Int_t iang=0; iang<400; iang++) {
    fNofPointsIncAng[iang] = 0;
    fNofRecoPointsIncAng[iang] = 0;
  }
  for (Int_t ipdg=-10000; ipdg<10000; ipdg++) {
    fNofRecoPdgSec    [ipdg] = 0;
    fNofPointsPdgSec  [ipdg] = 0;
    fNofRecoPdgPrim   [ipdg] = 0;
    fNofPointsPdgPrim [ipdg] = 0;
  }

  CreateHistos();
  Reset();
  if ( fOnlineAnalysis ) {
    recoCanvas = new TCanvas("StsRecoCanvas","Sts reconstruction",10,10,600,900);
//     TPad* recoPad[10];

    recoPad[0] = new TPad("titlePad",   "Title pad"                       ,0.00,0.90,1.00,1.00);
    recoPad[1] = new TPad("efficiencyPad","Efficiency pad"                ,0.00,0.60,0.40,0.90);
    recoPad[3] = new TPad("hitVspoint1Pad","Hit vs point on st.1 pad"     ,0.40,0.65,0.70,0.90);
    recoPad[4] = new TPad("hitPullst1","Hit vs sigma on st.1 pad"         ,0.70,0.65,1.00,0.90);
    recoPad[5] = new TPad("hitVspntXY1Pad","Hit vs point XY on st.1 pad"  ,0.40,0.40,0.70,0.65);
    recoPad[6] = new TPad("hitVssigmaXY6Pad","Hit vs sigma XY on st.1 pad",0.70,0.40,1.00,0.65);
    recoPad[2] = new TPad("Reco Points st1","Reco Points st1"         	  ,0.00,0.20,0.33,0.40);
    recoPad[7] = new TPad("Reco Points st2","Reco Points st2"             ,0.33,0.20,0.66,0.40);
    recoPad[8] = new TPad("Reco Points st3","Reco Points st3"             ,0.66,0.20,1.00,0.40);
    recoPad[9] = new TPad("Reco Points st5","Reco Points st5"             ,0.00,0.00,0.33,0.20);
    recoPad[10] = new TPad("Reco Points st6","Reco Points st6"            ,0.33,0.00,0.66,0.20);
    recoPad[11] = new TPad("Reco Points st8","Reco Points st8"            ,0.66,0.00,1.00,0.20);
    recoPad[12] = new TPad("printoutPad","Print information pad"          ,0.00,0.40,0.40,0.60);
//     recoPad[12]->SetLogz();
    for ( Int_t ipad = 0 ; ipad < 13 ; ipad++ ) {
	recoPad[ipad]->SetFillColor(0);
	recoPad[ipad]->SetBorderMode(0);
	recoPad[ipad]->Draw();
     }
  }
  // Control output
  if      (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
  else if (fVerbose >  2) fDigiScheme->Print(kTRUE);
  cout << "-I- " << fName << "::Init: "	<< "STS digitisation scheme succesfully initialised" << endl;
  cout << "    Stations: " << fDigiScheme->GetNStations() 
       << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
       << fDigiScheme->GetNChannels() << endl;
  
  return kSUCCESS;

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




// -----   Private method ReInit   -----------------------------------------
InitStatus CbmStsFindHitsQa::ReInit() {


  // Clear digitisation scheme
  fDigiScheme->Clear();

  // Build new digitisation scheme
  Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
  if ( ! success ) return kERROR;

  // Control output
  if      (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
  else if (fVerbose >  2) fDigiScheme->Print(kTRUE);
  cout << "-I- " << fName << "::Init: "
       << "STS digitisation scheme succesfully reinitialised" << endl;
  cout << "    Stations: " << fDigiScheme->GetNStations() 
       << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: " 
       << fDigiScheme->GetNChannels() << endl;
  
  return kSUCCESS;

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




// -----   Virtual method Finish   -----------------------------------------
void CbmStsFindHitsQa::Finish() {

  CbmStsStation_old* station = NULL;
  CbmStsSector*  sector  = NULL;
  Int_t NofRecoPoints = 0;
  Int_t HFE_reco = 0;
  Int_t HFE_point = 0;
  for ( Int_t istat = 0 ; istat < fDigiScheme->GetNStations() ; istat++ ) {
    station = fDigiScheme->GetStation(istat);
    for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
      sector = station->GetSector(isect);
      for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
        HFE_reco = 1;
        HFE_point = 1;
	if ( fNofPoints[istat][isect][isens] ) {
	  fhHitFindingEfficiency->Fill(2 << 24 | (istat+1) << 16 | (isect+1) << 4 | (isens+1) << 1,
				       ((Float_t)(fNofRecoPoints[istat][isect][isens]))/
				       ((Float_t)(fNofPoints[istat][isect][isens])));
	   HFE_reco += (fNofRecoPoints[istat][isect][isens]);
	   HFE_point += (fNofPoints[istat][isect][isens]);
	}
      }
    }
  }
  for (Int_t iPdg=-10000; iPdg<10000; iPdg++) {
    fhEffPdgSec  -> Fill(iPdg, ((Float_t)(fNofRecoPdgSec[iPdg]))  / ((Float_t)(fNofPointsPdgSec[iPdg])));
    fhEffPdgPrim -> Fill(iPdg, ((Float_t)(fNofRecoPdgPrim[iPdg])) / ((Float_t)(fNofPointsPdgPrim[iPdg])));
  }
 
  cout << endl;
  cout << "============================================================"
       << endl;
  cout << "===== " << fName << ": Run summary " << endl;
  cout << "===== Hit finding efficiency = "<<(((Float_t)(HFE_reco))/((Float_t)(HFE_point)))*100.<<"%"<< endl;
  cout << "============================================================"
       << endl;
	
  gDirectory->mkdir("STSFindHitsQA");
  gDirectory->cd("STSFindHitsQA");
  TIter next(fHistoList);
  while ( TH1* histo = ((TH1*)next()) ) histo->Write();

  for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
    gDirectory->mkdir(Form("Station%d",istat+1));
    gDirectory->cd(Form("Station%d",istat+1));
    TIter nextO(fHistoListPS[istat]);
    while ( TH1* histo = ((TH1*)nextO()) ) histo->Write();

    gDirectory->cd("..");
  }
  gDirectory->cd("..");
}					       
// -------------------------------------------------------------------------

ClassImp(CbmStsFindHitsQa)