Commit 8c138465 authored by Dominik Smith's avatar Dominik Smith Committed by Florian Uhlig
Browse files

Cleaned up CbmTrdHitProducerQa and added it to run_qa.C.

parent fa24ad9f
......@@ -172,7 +172,7 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
//run->AddTask(new CbmTrdHitRateQa()); //opens lots of windows
//run->AddTask(new CbmTrdDigitizerPRFQa()); //works put currently doesn't do anything
//run->AddTask(new CbmTrdHitRateFastQa()); //opens lots of windows
run->AddTask(new CbmTrdHitProducerQa()); //Histograms currently don't appear in output file
run->AddTask(new CbmTrdHitProducerQa());
}
// ------------------------------------------------------------------------
......
......@@ -9,6 +9,7 @@
#include "CbmTrdHitProducerQa.h"
#include "CbmDigiManager.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmTrdDigi.h"
......@@ -39,9 +40,8 @@ CbmTrdHitProducerQa::CbmTrdHitProducerQa() : CbmTrdHitProducerQa("TrdHitProducer
// ---- Standard constructor ------------------------------------------------
CbmTrdHitProducerQa::CbmTrdHitProducerQa(const char* name, const char*)
: FairTask(name)
, fOutFolder("TrdHitProducerQA", "TrdHitProducerQA")
, fTrdHitCollection(NULL)
, fTrdDigiCollection(NULL)
, fTrdDigiMatchCollection(NULL)
, fTrdPointCollection(NULL)
, fMCTrackArray(NULL)
, fNoTrdStations(-1)
......@@ -73,6 +73,20 @@ CbmTrdHitProducerQa::~CbmTrdHitProducerQa() {}
// ---- Initialisation ------------------------------------------------------
InitStatus CbmTrdHitProducerQa::Init()
{
fOutFolder.Clear();
histFolder = fOutFolder.AddFolder("hist", "Histogramms");
histFolder->Add(fHitPoolsX);
histFolder->Add(fHitPoolsY);
histFolder->Add(S1L1eTR15);
histFolder->Add(S1L1edEdx15);
histFolder->Add(S1L1edE15);
histFolder->Add(S1L1edEall);
histFolder->Add(S1L1pidE15);
histFolder->Add(S1L1pidEall);
histFolder->Add(S3L4edEall);
histFolder->Add(S3L4pidEall);
// Get pointer to the ROOT I/O manager
FairRootManager* rootMgr = FairRootManager::Instance();
if (NULL == rootMgr) {
......@@ -81,27 +95,20 @@ InitStatus CbmTrdHitProducerQa::Init()
return kFATAL;
}
// Get pointer to TRD point array
fTrdPointCollection = (TClonesArray*) rootMgr->GetObject("TrdPoint");
if (NULL == fTrdPointCollection) {
// Get a pointer to the previous already existing data level
fDigiMan = CbmDigiManager::Instance();
if (NULL == fDigiMan) {
cout << "-W- CbmTrdHitProducerQa::Init : "
<< "no TRD point array !" << endl;
<< "no digi manager found !" << endl;
return kERROR;
}
fDigiMan->Init();
// Get pointer to TRD digi array
fTrdDigiCollection = (TClonesArray*) rootMgr->GetObject("TrdDigi");
if (NULL == fTrdDigiCollection) {
cout << "-W- CbmTrdHitProducerQa::Init : "
<< "no TRD digi array !" << endl;
return kERROR;
}
// Get pointer to TRD digi array
fTrdDigiMatchCollection = (TClonesArray*) rootMgr->GetObject("TrdDigiMatch");
if (NULL == fTrdDigiMatchCollection) {
// Get pointer to TRD point array
fTrdPointCollection = (TClonesArray*) rootMgr->GetObject("TrdPoint");
if (NULL == fTrdPointCollection) {
cout << "-W- CbmTrdHitProducerQa::Init : "
<< "no TRD digi match array !" << endl;
<< "no TRD point array !" << endl;
return kERROR;
}
......@@ -120,6 +127,12 @@ InitStatus CbmTrdHitProducerQa::Init()
return kFATAL;
}
// Get pointer to TRD digi array match
if (!fDigiMan->IsMatchPresent(ECbmModuleId::kTrd)) {
cout << GetName() << ": no TRD match branch in digi manager." << endl;
return kERROR;
}
return kSUCCESS;
}
......@@ -129,58 +142,36 @@ InitStatus CbmTrdHitProducerQa::Init()
// ---- Task execution ------------------------------------------------------
void CbmTrdHitProducerQa::Exec(Option_t*)
{
// Declare variables
CbmTrdHit* trdHit = NULL;
// CbmTrdDigi* trdDigi = NULL;
CbmMatch* trdDigiMatch = NULL;
CbmTrdPoint* trdPoint = NULL;
// CbmMCTrack* mctrack = NULL;
Float_t hitPoolX = 0;
Float_t hitPoolY = 0;
Float_t hitPosX = 0;
Float_t hitPosY = 0;
// Float_t hitErrX = 0;
// Float_t hitErrY = 0;
Float_t pointPosX = 0;
Float_t pointPosY = 0;
Int_t plane = 0;
Int_t partID = 0;
Float_t momentum;
// Event counters
Int_t nTrdHit = fTrdHitCollection->GetEntriesFast();
// Loop over TRD hits
for (Int_t iHit = 0; iHit < nTrdHit; iHit++) {
trdHit = (CbmTrdHit*) fTrdHitCollection->At(iHit);
const CbmTrdHit* trdHit = (CbmTrdHit*) fTrdHitCollection->At(iHit);
if (NULL == trdHit) continue;
// This will have to change in the future, when the creation of the poin
// will not be necessarily connected to existence of tyhe point
// This will have to change in the future, when the creation of the point
// will not be necessarily connected to existence of the point
trdDigiMatch = (CbmMatch*) fTrdDigiMatchCollection->At(trdHit->GetRefId());
const CbmMatch* trdDigiMatch = fDigiMan->GetMatch(ECbmModuleId::kTrd, trdHit->GetRefId());
if (NULL == trdDigiMatch) continue;
trdPoint = (CbmTrdPoint*) fTrdPointCollection->At(trdDigiMatch->GetMatchedLink().GetIndex());
const CbmTrdPoint* trdPoint = (CbmTrdPoint*) fTrdPointCollection->At(trdDigiMatch->GetMatchedLink().GetIndex());
if (NULL == trdPoint) continue;
plane = trdHit->GetPlaneId();
const Int_t plane = trdHit->GetPlaneId();
if (plane == 1) {
partID = (((CbmMCTrack*) fMCTrackArray->At(trdPoint->GetTrackID()))->GetPdgCode());
partID = (((CbmMCTrack*) fMCTrackArray->At(trdPoint->GetTrackID()))->GetPdgCode());
momentum = TMath::Sqrt((trdPoint->GetPx() * trdPoint->GetPx()) + (trdPoint->GetPy() * trdPoint->GetPy())
+ (trdPoint->GetPz() * trdPoint->GetPz()));
if ((TMath::Abs(partID) == 11) && (momentum > 1.25) && (momentum < 1.75)) {
S1L1eTR15->Fill(0); //(trdHit->GetELossTR())*1000000);
S1L1edEdx15->Fill((trdHit->GetELoss()) * 1000000);
S1L1edE15->Fill((trdHit->GetELoss()) * 1000000);
......@@ -193,32 +184,26 @@ void CbmTrdHitProducerQa::Exec(Option_t*)
}
if (plane == 12) {
partID = (((CbmMCTrack*) fMCTrackArray->At(trdPoint->GetTrackID()))->GetPdgCode());
partID = (((CbmMCTrack*) fMCTrackArray->At(trdPoint->GetTrackID()))->GetPdgCode());
momentum = TMath::Sqrt((trdPoint->GetPx() * trdPoint->GetPx()) + (trdPoint->GetPy() * trdPoint->GetPy())
+ (trdPoint->GetPz() * trdPoint->GetPz()));
if (TMath::Abs(partID) == 11) { S3L4edEall->Fill((trdHit->GetELoss()) * 1000000); }
if (TMath::Abs(partID) == 211) { S3L4pidEall->Fill((trdHit->GetELoss()) * 1000000); }
}
if (plane == 1) {
// compute the Hit pools for X and Y coordinate
hitPosX = trdHit->GetX();
const Float_t hitPosX = trdHit->GetX();
// hitErrX = trdHit->GetDx();
// hitErrX /= 10000; // micrometers to centimeters
pointPosX = trdPoint->GetX();
hitPoolX = (hitPosX - pointPosX); ///hitErrX;
const Float_t pointPosX = trdPoint->GetX();
const Float_t hitPoolX = (hitPosX - pointPosX); ///hitErrX;
hitPosY = trdHit->GetY();
const Float_t hitPosY = trdHit->GetY();
// hitErrY = trdHit->GetDy();
// hitErrY /= 10000; // micrometers to centimeters
pointPosY = trdPoint->GetY();
hitPoolY = (hitPosY - pointPosY); ///hitErrY;
const Float_t pointPosY = trdPoint->GetY();
const Float_t hitPoolY = (hitPosY - pointPosY); ///hitErrY;
// fill histograms
fHitPoolsX->Fill(hitPoolX);
......@@ -226,82 +211,18 @@ void CbmTrdHitProducerQa::Exec(Option_t*)
}
}
}
// --------------------------------------------------------------------------
// ---- Finish --------------------------------------------------------------
void CbmTrdHitProducerQa::Finish() { WriteHistograms(); }
// --------------------------------------------------------------------------
// ---- Prepare test histograms ---------------------------------------------
/*
void CbmTrdHitProducerQa::PrepareHistograms()
{
fHitPoolsX=NULL;
fHitPoolsY=NULL;
S1L1eTR15=NULL;
S1L1edEdx15=NULL;
S1L1edE15=NULL;
S1L1edEall=NULL;
S1L1pidE15=NULL;
S1L1pidEall=NULL;
S3L4eTR15=NULL;
S3L4edEdx15=NULL;
S3L4edE15=NULL;
S3L4edEall=NULL;
S3L4pidE15=NULL;
S3L4pidEall=NULL;
fHitPoolsX = new TH1F("fHitPoolsX", "", 500, -50, 50);
fHitPoolsY = new TH1F("fHitPoolsY", "", 500, -50, 50);
S1L1eTR15 = new TH1F("S1L1eTR15","TR of e- for first layer ",600,0.,60.);
S1L1edEdx15 = new TH1F("S1L1edEdx15","dEdx of e- for first layer ",600,0.,60.);
S1L1edE15 = new TH1F("S1L1edE15","dEdx+TR of e- for first layer ",600,0.,60.);
S1L1edEall = new TH1F("S1L1edEall","dEdx+TR of e- for first layer ",600,0.,60.);
S1L1pidE15 = new TH1F("S1L1pidE15","dEdx+TR of pi- for first layer ",600,0.,60.);
S1L1pidEall = new TH1F("S1L1pidEall","dEdx+TR of pi- for first layer ",600,0.,60.);
S3L4edEall = new TH1F("S3L4edEall","dEdx+TR of e- for layer 12",600,0.,60.);
S3L4pidEall = new TH1F("S3L4pidEall","dEdx+TR of pi- for layer 12",600,0.,60.);
}
*/
// --------------------------------------------------------------------------
// ---- Write test histograms ------------------------------------------------
void CbmTrdHitProducerQa::WriteHistograms()
{
gDirectory->mkdir("TrdHitProducerQA");
gDirectory->cd("TrdHitProducerQA");
if (fHitPoolsX) fHitPoolsX->Write();
if (fHitPoolsY) fHitPoolsY->Write();
if (S1L1eTR15) S1L1eTR15->Write();
if (S1L1edEdx15) S1L1edEdx15->Write();
if (S1L1edE15) S1L1edE15->Write();
if (S1L1edEall) S1L1edEall->Write();
if (S1L1pidE15) S1L1pidE15->Write();
if (S1L1pidEall) S1L1pidEall->Write();
if (S3L4eTR15) S3L4eTR15->Write();
if (S3L4edEdx15) S3L4edEdx15->Write();
if (S3L4edE15) S3L4edE15->Write();
if (S3L4edEall) S3L4edEall->Write();
if (S3L4pidE15) S3L4pidE15->Write();
if (S3L4pidEall) S3L4pidEall->Write();
gDirectory->cd("..");
FairSink* sink = FairRootManager::Instance()->GetSink();
sink->WriteObject(&fOutFolder, nullptr);
}
// --------------------------------------------------------------------------
ClassImp(CbmTrdHitProducerQa);
......@@ -23,6 +23,9 @@
#include "FairTask.h"
#include <TFolder.h>
class CbmDigiManager;
class TClonesArray;
class TH1F;
......@@ -46,12 +49,14 @@ public:
/* Finish at the end of each event */
virtual void Finish();
private:
TFolder fOutFolder; /// output folder with histos and canvases
TFolder* histFolder; /// subfolder for histograms
CbmDigiManager* fDigiMan = nullptr;
/* Data branches*/
TClonesArray* fTrdHitCollection;
TClonesArray* fTrdDigiCollection;
TClonesArray* fTrdDigiMatchCollection;
TClonesArray* fTrdPointCollection;
TClonesArray* fMCTrackArray;
......@@ -61,6 +66,8 @@ private:
/** Number of layers per station **/
Int_t fNoTrdPerStation;
/* Write test histograms */
void WriteHistograms();
/* Test histograms*/
TH1F* fHitPoolsX; // = ((Hit - Point) / HitError) in X
......@@ -80,18 +87,7 @@ private:
TH1F* S3L4pidE15; //
TH1F* S3L4pidEall; //
private:
CbmTrdHitProducerQa(const CbmTrdHitProducerQa&);
CbmTrdHitProducerQa& operator=(const CbmTrdHitProducerQa&);
//public:
/* Write test histograms */
void WriteHistograms();
ClassDef(CbmTrdHitProducerQa, 2)
ClassDef(CbmTrdHitProducerQa, 3)
};
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment