CbmStsSimulationQa.cxx

Go to the documentation of this file.

// -------------------------------------------------------------------------
// -----                  CbmStsSimulationQa source file               -----
// -----                  Created 02/02/07  by R. Karabowicz           -----
// -------------------------------------------------------------------------
 
#include "CbmStsSimulationQa.h"
 
#include "CbmStsPoint.h"
// #include "CbmStsTrack.h"
// #include "CbmStsTrackMatch.h"
#include "BmnDetectorList.h"
#include "CbmGeoStsPar.h"
#include "CbmMCTrack.h"
#include "FairGeoNode.h"
#include "FairGeoPassivePar.h"
#include "FairGeoTransform.h"
#include "FairGeoVector.h"
#include "FairMCPoint.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TLegend.h"
#include "TList.h"
#include "TPad.h"
#include "TPaveText.h"
#include "TVector3.h"
 
#include <iomanip>
#include <iostream>
 
using std::cout;
using std::endl;
using std::flush;
using std::setprecision;
 
// -----   Default constructor   -------------------------------------------
CbmStsSimulationQa::CbmStsSimulationQa()
    : fMCTracks()
    , fSTSPoints()
    , fPassGeo()
    , fStsGeo()
    , fTargetPos()
    , fNStations(0)
    , fhMomAll()
    , fhYPtMapAll()
    , fhPdgCodeAll()
    , fhStsPointsAll()
    , fhMomRec()
    , fhYPtMapRec()
    , fhPdgCodeRec()
    , fhStsPointsRec()
    , fhMomStsPoints()
    , fhStsPointsPosition()
    , fhNofEvents()
    , fhNofStsStations()
    , fHistoList()
    , fNEvents(0)
    , fOnlineAnalysis(kFALSE)
    , fOnlineCanvas()
{}
 
// -----   Standard constructor   ------------------------------------------
CbmStsSimulationQa::CbmStsSimulationQa(Bool_t visualizeBool, Int_t iVerbose)
    : FairTask("STS Simulation QA", iVerbose)
    , fMCTracks()
    , fSTSPoints()
    , fPassGeo()
    , fStsGeo()
    , fTargetPos()
    , fNStations(0)
    , fhMomAll()
    , fhYPtMapAll()
    , fhPdgCodeAll()
    , fhStsPointsAll()
    , fhMomRec()
    , fhYPtMapRec()
    , fhPdgCodeRec()
    , fhStsPointsRec()
    , fhMomStsPoints()
    , fhStsPointsPosition()
    , fhNofEvents()
    , fhNofStsStations()
    , fHistoList()
    , fNEvents(0)
    , fOnlineAnalysis(kFALSE)
    , fOnlineCanvas()
{
    fOnlineAnalysis = visualizeBool;
}
 
// -----   Destructor   ----------------------------------------------------
CbmStsSimulationQa::~CbmStsSimulationQa()
{
 
    fHistoList->Delete();
    delete fHistoList;
}
 
// -----   Public method SetParContainers   --------------------------------
void CbmStsSimulationQa::SetParContainers()
{
 
    // Get Run
    FairRunAna* run = FairRunAna::Instance();
    if (!run) {
        cout << "-E- " << GetName() << "::SetParContainers: No FairRunAna!" << endl;
        return;
    }
 
    // Get Runtime Database
    FairRuntimeDb* runDb = run->GetRuntimeDb();
    if (!run) {
        cout << "-E- " << GetName() << "::SetParContainers: No runtime database!" << endl;
        return;
    }
 
    // Get passive geometry parameters
    fPassGeo = (FairGeoPassivePar*)runDb->getContainer("FairGeoPassivePar");
    if (!fPassGeo) {
        cout << "-E- " << GetName() << "::SetParContainers: "
             << "No passive geometry parameters!" << endl;
        return;
    }
 
    // Get STS geometry parameters
    fStsGeo = (CbmGeoStsPar*)runDb->getContainer("CbmGeoStsPar");
    if (!fStsGeo) {
        cout << "-E- " << GetName() << "::SetParContainers: "
             << "No Sts geometry parameters!" << endl;
        return;
    }
}
// -------------------------------------------------------------------------
 
// -----   Public method Init   --------------------------------------------
InitStatus CbmStsSimulationQa::Init()
{
 
    cout << "===========================================================" << endl;
    ;
    cout << GetName() << ": Initialising..." << endl;
 
    // Get RootManager
    FairRootManager* ioman = FairRootManager::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 kERROR;
    }
 
    // Get StsPoints array
    fSTSPoints = (TClonesArray*)ioman->GetObject("StsPoint");
    if (!fSTSPoints) {
        // cout << "-E- " << GetName() << "::Init: No StsPoint array!" << endl;
        return kERROR;
    }
 
    // Get the geometry of target and Sts
    InitStatus geoStatus = GetGeometry();
    if (geoStatus != kSUCCESS) {
        cout << "-E- " << GetName() << "::Init: Error in reading geometry!" << endl;
        return geoStatus;
    }
 
    // Create histograms
    CreateHistos();
    Reset();
    if (fOnlineAnalysis) {
        fOnlineCanvas = new TCanvas("StsSimOnline", "Sts simulation online", 10, 10, 600, 600);
        fOnlinePad[0] = new TPad("titlePad", "Title pad", 0.00, 0.90, 1.00, 1.00);
        fOnlinePad[1] = new TPad("momentumPad", "Momentum pad", 0.00, 0.35, 0.50, 0.90);
        fOnlinePad[2] = new TPad("printoutPad", "Print information pad    ", 0.10, 0.10, 0.35, 0.35);
        fOnlinePad[3] = new TPad("pointPad", "Points per track pad", 0.50, 0.50, 1.00, 0.90);
        fOnlinePad[4] = new TPad("stationPad1", "Points on 1st station pad", 0.50, 0.25, 0.75, 0.50);
        fOnlinePad[5] = new TPad("stationPad2", "Points on 3rd station pad", 0.75, 0.25, 1.00, 0.50);
        fOnlinePad[6] = new TPad("stationPad3", "Points on 5th station pad", 0.50, 0.00, 0.75, 0.25);
        fOnlinePad[7] = new TPad("stationPad4", "Points on 7th station pad", 0.75, 0.00, 1.00, 0.25);
        fOnlinePad[1]->SetLogy();
        for (Int_t ipad = 0; ipad < 8; ipad++) {
            fOnlinePad[ipad]->SetFillColor(0);
            fOnlinePad[ipad]->SetBorderMode(0);
            fOnlinePad[ipad]->Draw();
        }
 
        fOnlinePad[0]->cd();
        TLegend* brp = new TLegend(0.1, 0.1, 0.9, 0.9, "Online Sts simulation");
        brp->SetTextAlign(22);
        brp->SetTextSize(0.6);
        brp->SetTextColor(1);
        brp->SetBorderSize(0);
        brp->SetFillColor(0);
        brp->Draw();
        fOnlinePad[0]->Update();
    }
 
    fhNofStsStations->SetBinContent(1, fNStations);
 
    // Output
    cout << "   Number of Sts stations : " << fNStations << endl;
    if (fActive)
        cout << "   *****   Task is ACTIVE   *****" << endl;
    cout << "===========================================================" << endl << endl;
 
    return geoStatus;
}
// -------------------------------------------------------------------------
 
// -----   Public method ReInit   ------------------------------------------
InitStatus CbmStsSimulationQa::ReInit()
{
 
    cout << "===========================================================" << endl;
    ;
    cout << GetName() << ": Reinitialising..." << endl;
 
    // Get the geometry of target and Sts
    InitStatus geoStatus = GetGeometry();
    if (geoStatus != kSUCCESS) {
        cout << "-E- " << GetName() << "::ReInit: Error in reading geometry!" << endl;
        return geoStatus;
    }
 
    // Output
    cout << "   Number of Sts stations : " << fNStations << endl;
    if (fActive)
        cout << "   *****   Task is ACTIVE   *****" << endl;
    cout << "===========================================================" << endl << endl;
 
    return geoStatus;
}
// -------------------------------------------------------------------------
 
// -----   Public method Exec   --------------------------------------------
void CbmStsSimulationQa::Exec(Option_t* opt)
{
    //  cout << "EXEC" << endl;
 
    Int_t nofMCTracks = fMCTracks->GetEntriesFast();
    Int_t nofSTSPoints = fSTSPoints->GetEntriesFast();
 
    for (Int_t itr = 0; itr < nofMCTracks; itr++) {
        CbmMCTrack* mctrack = (CbmMCTrack*)fMCTracks->At(itr);
        Int_t pdgCode = mctrack->GetPdgCode();
 
        // reject funny particles
        if ((pdgCode == 10010020) || (pdgCode == 10010030) || (pdgCode == 50000050) || (pdgCode == 50010051)
            || (pdgCode == 10020040))
            continue;
 
        TVector3 mom;
        mctrack->GetMomentum(mom);
        Float_t pT = mom.Pt();
        Float_t p = mom.Mag();
 
        TVector3 startvtx;
        mctrack->GetStartVertex(startvtx);
        Float_t vertexZ = startvtx.z();
 
        TLorentzVector mom4;
        mctrack->Get4Momentum(mom4);
        Float_t rapidity = mom4.Rapidity();
 
        Int_t stsPoints = mctrack->GetNPoints(kGEM);
 
        if (stsPoints > 0 && vertexZ <= 100) {
            fhMomAll->Fill(p);
            fhYPtMapAll->Fill(rapidity, pT);
            fhPdgCodeAll->Fill(pdgCode);
            fhStsPointsAll->Fill(stsPoints);
            fhMomStsPoints->Fill(stsPoints, p);
        }
        if (stsPoints > 3 && vertexZ <= 100) {
            fhMomRec->Fill(p);
            fhYPtMapRec->Fill(rapidity, pT);
            fhPdgCodeRec->Fill(pdgCode);
            fhStsPointsRec->Fill(stsPoints);
        }
    }
 
    //  cout << "track loop done" << endl;
 
    for (Int_t ipnt = 0; ipnt < nofSTSPoints; ipnt++) {
        CbmStsPoint* stsPoint = (CbmStsPoint*)fSTSPoints->At(ipnt);
        Float_t z = stsPoint->GetZ();    // [cm]
        Float_t x = stsPoint->GetX(z);   // [cm]
        Float_t y = stsPoint->GetY(z);   // [cm]
 
        fhStsPointsPosition->Fill(z, x, y);
        //     cout << "filled 3d,    MCId = " << stsPoint->GetDetectorID()
        //   << "             nr = " << fStationNrFromMcId[stsPoint->GetDetectorID()] << endl;
        fhStationPoints[fStationNrFromMcId[stsPoint->GetDetectorID()]]->Fill(x, y);
        //     cout << "filled 2d for MCId = " << stsPoint->GetDetectorID()
        //   << " and station nr = " << fStationNrFromMcId[stsPoint->GetDetectorID()] << endl;
    }
 
    //  cout << "hello" << endl;
 
    Double_t tracksPerEvent = (Double_t)(fhMomAll->GetEntries()) / ((Double_t)fNEvents + 1);
    Double_t pointsPerEvent = (Double_t)(fhStsPointsPosition->GetEntries()) / ((Double_t)fNEvents + 1);
 
    if (fOnlineAnalysis) {
        fOnlinePad[1]->cd();
        fhMomRec->Draw();
        fOnlinePad[1]->Update();
        fOnlinePad[2]->cd();
        TPaveText* printoutPave = new TPaveText(0.1, 0.1, 0.9, 0.9);
        printoutPave->SetTextAlign(22);
        printoutPave->SetTextSize(0.1);
        printoutPave->SetTextColor(1);
        printoutPave->SetBorderSize(0);
        printoutPave->SetFillColor(0);
        printoutPave->AddText(Form("%i events", fNEvents + 1));
        printoutPave->AddText(Form("tracks/event = %3.2f", tracksPerEvent));
        printoutPave->AddText(Form("points/event = %3.2f", pointsPerEvent));
        fOnlinePad[2]->Clear();
        printoutPave->Draw();
        fOnlinePad[2]->Update();
        fOnlinePad[3]->cd();
        fhStsPointsRec->Draw();
        fOnlinePad[3]->Update();
 
        if (fNStations) {
            fOnlinePad[4]->cd();
            fhStationPoints[0]->Draw("colz");
            fOnlinePad[4]->Update();
            if (fNStations > 2) {
                fOnlinePad[5]->cd();
                fhStationPoints[2]->Draw("colz");
                fOnlinePad[5]->Update();
                if (fNStations > 4) {
                    fOnlinePad[6]->cd();
                    fhStationPoints[4]->Draw("colz");
                    fOnlinePad[6]->Update();
                    if (fNStations > 6) {
                        fOnlinePad[7]->cd();
                        fhStationPoints[6]->Draw("colz");
                        fOnlinePad[7]->Update();
                    }
                }
            }
        }
    }
 
    //  cout << "\rEvent #" << fNEvents+1 << flush;
    cout << endl << endl;
    cout << "=======================================================" << endl;
    cout << "===== StsSimulationQa: " << endl;
    cout << "===== \rEvent #" << fNEvents + 1 << endl;
    cout << "===== " << setprecision(6) << tracksPerEvent << " tracks/event" << endl;
    cout << "===== " << setprecision(7) << pointsPerEvent << " points/event" << endl;
    cout << "=======================================================" << endl;
 
    fNEvents++;
 
    fhNofEvents->SetBinContent(1, fNEvents);
}
// -------------------------------------------------------------------------
 
// -----   Private method Finish   -----------------------------------------
void CbmStsSimulationQa::Finish()
{
    // Run summary to screen
    //   cout << endl << endl;
    //   cout << "======================================================="<< endl;
    //   cout << "            StsSimulationQa: Run summary" << endl << endl;
    //   cout << "======================================================="
    //        << endl;
    //   cout << endl << endl;
 
    gDirectory->mkdir("STSSimulationQA");
    gDirectory->cd("STSSimulationQA");
    TIter next(fHistoList);
    while (TH1* histo = ((TH1*)next()))
        histo->Write();
    gDirectory->cd("..");
}
// -------------------------------------------------------------------------
 
// -----   Private method GetGeometry   ------------------------------------
InitStatus CbmStsSimulationQa::GetGeometry()
{
 
    cout << "GET GEOMETRY" << endl;
    // Get Sts geometry
    if (!fStsGeo) {
        cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!" << endl;
        fNStations = 0;
        return kERROR;
    }
    TObjArray* stsNodes = fStsGeo->GetGeoSensitiveNodes();
    if (!stsNodes) {
        cout << "-E- " << GetName() << "::GetGeometry: No Sts node array" << endl;
        fNStations = 0;
        return kERROR;
    }
    Int_t tempNofStations = stsNodes->GetEntries();
 
    //   cout << "There are " << tempNofStations << " nodes" << (tempNofStations > 10 ? "!!!" : "" ) << endl;
 
    TString geoNodeName;
    fNStations = 0;
    TString stationNames[10];
    for (Int_t ist = 0; ist < tempNofStations; ist++) {
        FairGeoNode* stsNode = (FairGeoNode*)stsNodes->At(ist);
        if (!stsNode) {
            cout << "-W- CbmStsDigiScheme::Init: station#" << ist << " not found among sensitive nodes " << endl;
            continue;
        }
        geoNodeName = stsNode->GetName();
        TArrayD* params = stsNode->getParameters();
 
        Bool_t stationKnown = kFALSE;
        // check if the node belongs to some station, save the MCId and outer radius
        for (Int_t ikst = 0; ikst < fNStations; ikst++)
            if (geoNodeName.Contains(stationNames[ikst])) {
                FairGeoTransform* transform = stsNode->getLabTransform();
                FairGeoVector translat = transform->getTranslation();
                FairGeoTransform center = stsNode->getCenterPosition();
                FairGeoVector centerV = center.getTranslation();
                Double_t radialSize = TMath::Sqrt(
                    (translat.X() + centerV.X() + params->At(0)) * (translat.X() + centerV.X() + params->At(0))
                    + (translat.Y() + centerV.Y() + params->At(1)) * (translat.Y() + centerV.Y() + params->At(1)));
                if (radialSize > fStationRadius[ikst]) {
                    //    cout << "Found bigger radius for station " << ikst << " -> " << radialSize << endl;
                    fStationRadius[ikst] = radialSize;
                }
                fStationNrFromMcId[stsNode->getMCid()] = ikst;
                stationKnown = kTRUE;
            }
 
        if (stationKnown)
            continue;
 
        // if not known, register it and save MCId
        fStationNrFromMcId[stsNode->getMCid()] = fNStations;
 
        if (geoNodeName.Length() == 12)
            fStationRadius[fNStations] = params->At(1);
        else {
            FairGeoTransform* transform = stsNode->getLabTransform();
            FairGeoVector translat = transform->getTranslation();
            FairGeoTransform center = stsNode->getCenterPosition();
            FairGeoVector centerV = center.getTranslation();
            Double_t radialSize = TMath::Sqrt(
                (translat.X() + centerV.X() + params->At(0)) * (translat.X() + centerV.X() + params->At(0))
                + (translat.Y() + centerV.Y() + params->At(1)) * (translat.Y() + centerV.Y() + params->At(1)));
            fStationRadius[fNStations] = radialSize;
        }
        // it will work only if the node name is organized as:
        // for station name is "stsstationXX", where XX is the station number (f.e. XX=07 for station number 7)
        // for sector  name is "stsstationXXanythingHereToDistinguishDifferentSectors"
        geoNodeName.Remove(12, geoNodeName.Length() - 12);
        stationNames[fNStations] = geoNodeName.Data();
        fNStations++;
 
        //     cout << "station #" << fNStations << " has MCID = " << stsNode->getMCid() << " and name " <<
        //     stsNode->GetName() << endl;
 
        //    fStationsMCId[fNStations] = stsNode->getMCid(); // not used
    }
    //   cout << "There are " << fNStations << " stations" << endl;
 
    return kSUCCESS;
}
// -------------------------------------------------------------------------
 
// -----   Private method CreateHistos   -----------------------------------
void CbmStsSimulationQa::CreateHistos()
{
 
    fHistoList = new TList();
 
    fhMomAll = new TH1F("hMomAll", "Momentum - all in STS", 100, 0, 50);
    fhMomAll->SetXTitle("p [GeV/c]");
    fhMomAll->SetYTitle("yield [a.u.]");
    fhYPtMapAll = new TH2F("hYPtMapAll", "Rapidity - trans. mom. map - all in STS", 100, -3, 7, 100, 0, 5);
    fhYPtMapAll->SetXTitle("rapidity");
    fhYPtMapAll->SetYTitle("p_{t} [GeV/c]");
    fhPdgCodeAll = new TH1F("hPdgCodeAll", "PDG code - all in STS", 1000, -500, 500);
    fhStsPointsAll =
        new TH1F("hStsPointsAll", "STSPoints per track - all in STS", 2 * fNStations, 0.5, 2 * fNStations + 0.5);
    fhStsPointsAll->SetXTitle("nof points");
    fhStsPointsAll->SetYTitle("yield [a.u.]");
    fhMomRec = new TH1F("hMomRec", "Momentum - rec in STS", 100, 0, 50);
    fhMomRec->SetXTitle("p [GeV/c]");
    fhMomRec->SetYTitle("yield [a.u.]");
    fhYPtMapRec = new TH2F("hYPtMapRec", "Rapidity - trans. mom. map - rec in STS", 100, -3, 7, 100, 0, 5);
    fhYPtMapRec->SetXTitle("rapidity");
    fhYPtMapRec->SetYTitle("p_{t} [GeV/c]");
    fhPdgCodeRec = new TH1F("hPdgCodeRec", "PDG code - rec in STS", 1000, -500, 500);
    fhStsPointsRec =
        new TH1F("hStsPointsRec", "STSPoints per track - rec in STS", 2 * fNStations, 0.5, 2 * fNStations + 0.5);
    fhStsPointsRec->SetXTitle("nof points");
    fhStsPointsRec->SetYTitle("yield [a.u.]");
    fhMomStsPoints = new TH2F("hMomStsPoints", "momentum vs STSPoints per track", 1000, 0, 9, 100, 0, 50);
    fhStsPointsPosition = new TH3F("hStsPointsPosition", "STS hits", 100, 0, 100, 100, -50, 50, 100, -50, 50);
    fHistoList->Add(fhMomAll);
    fHistoList->Add(fhYPtMapAll);
    fHistoList->Add(fhPdgCodeAll);
    fHistoList->Add(fhStsPointsAll);
    fHistoList->Add(fhMomRec);
    fHistoList->Add(fhYPtMapRec);
    fHistoList->Add(fhPdgCodeRec);
    fHistoList->Add(fhStsPointsRec);
    fHistoList->Add(fhMomStsPoints);
    fHistoList->Add(fhStsPointsPosition);
 
    for (Int_t ist = 0; ist < fNStations; ist++) {
        Int_t histSize = (Int_t)(1.05 * fStationRadius[ist]) + 1;
        fhStationPoints[ist] = new TH2F(Form("hStationPoints%i", ist + 1), Form("Points at station %i", ist + 1),
                                        40 * histSize, -histSize, histSize, 40 * histSize, -histSize, histSize);
        fhStationPoints[ist]->SetXTitle("x [cm]");
        fhStationPoints[ist]->SetYTitle("y [cm]");
        fHistoList->Add(fhStationPoints[ist]);
    }
 
    fhNofEvents = new TH1F("hNofEvents", "Number of events", 1, 0., 1.);
    fhNofStsStations = new TH1F("hNofStsStations", "Number of stations", 1, 0., 1.);
    fHistoList->Add(fhNofEvents);
    fHistoList->Add(fhNofStsStations);
}
// -------------------------------------------------------------------------
 
// -----   Private method Reset   ------------------------------------------
void CbmStsSimulationQa::Reset()
{
    fNEvents = 0;
}
// -------------------------------------------------------------------------