diff --git a/reco/detectors/psd/CbmPsdHitProducer.cxx b/reco/detectors/psd/CbmPsdHitProducer.cxx
index f646df175355e7817b80bd0b965bec36dd3f52c0..46e174427be0a255fded646b0ab41a1c861a098e 100644
--- a/reco/detectors/psd/CbmPsdHitProducer.cxx
+++ b/reco/detectors/psd/CbmPsdHitProducer.cxx
@@ -5,40 +5,40 @@
#include "CbmPsdHitProducer.h"
#include "CbmDigiManager.h"
+#include "CbmEvent.h"
#include "CbmPsdDigi.h"
#include "CbmPsdHit.h"
-#include "FairRootManager.h"
+#include <FairRootManager.h>
#include <Logger.h>
-#include "TClonesArray.h"
-#include "TMath.h"
+#include <TClonesArray.h>
#include <TFile.h>
+#include <TMath.h>
+#include <TStopwatch.h>
#include <fstream>
+#include <iomanip>
#include <iostream>
#include <map>
using std::cout;
using std::endl;
+using std::fixed;
+using std::left;
+using std::pair;
+using std::setprecision;
+using std::setw;
// ----- Default constructor -------------------------------------------
-CbmPsdHitProducer::CbmPsdHitProducer()
- : FairTask("Ideal Psd Hit Producer", 1)
- , fNHits(0)
- , fHitArray(NULL)
- , fDigiMan(nullptr)
- , fXi()
- , fYi()
- , fhModXNewEn(NULL) {
- // Reset();
-}
+CbmPsdHitProducer::CbmPsdHitProducer() : FairTask("Ideal Psd Hit Producer", 1), fXi(), fYi() {}
// -------------------------------------------------------------------------
// ----- Destructor ----------------------------------------------------
-CbmPsdHitProducer::~CbmPsdHitProducer() {
+CbmPsdHitProducer::~CbmPsdHitProducer()
+{
if (fHitArray) {
fHitArray->Delete();
delete fHitArray;
@@ -48,139 +48,201 @@ CbmPsdHitProducer::~CbmPsdHitProducer() {
// ----- Public method Init --------------------------------------------
-InitStatus CbmPsdHitProducer::Init() {
-
- //ifstream fxypos("psd_geo_xy.txt");
- //for (Int_t ii=1; ii<50; ii++) { //SELIM: 49 modules, including central & corner modules (rejected in analysis/flow/eventPlane.cxx)
- // fxypos>>fXi[ii]>>fYi[ii];
- // cout<<ii<<" "<<fXi[ii]<<" "<<fYi[ii]<<endl;
- //}
- //fxypos.close();
+InitStatus CbmPsdHitProducer::Init()
+{
fhModXNewEn = new TH1F("hModXNewEn", "X distr, En", 300, -150., 150.);
fhModXNewEn->Print();
- // Get RootManager
+ std::cout << std::endl << std::endl;
+ LOG(info) << "======= " << GetName() << ": Init =====================";
+
+ // --- Get RootManager
FairRootManager* ioman = FairRootManager::Instance();
if (!ioman) {
- LOG(fatal)
- << "-W- CbmPsdHitProducer::Init: RootManager not instantised!"; //FLORIAN & SELIM
+ LOG(fatal) << "-W- CbmPsdHitProducer::Init: RootManager not instantised!"; //FLORIAN & SELIM
return kFATAL;
}
+ // --- Get event array, if present
+ fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("Event"));
+ if (fEvents) LOG(info) << GetName() << ": found Event branch, run event-by-event";
+ else {
+ fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
+ if (fEvents) LOG(info) << GetName() << ": found CbmEvent branch, run event-by-event";
+ }
+ if (!fEvents) {
+ LOG(info) << GetName() << ": no event branch found; run in time-based mode";
+ LOG(warn) << "*** Note that the current PSD hit producer is not suited for time-based reconstruction!";
+ LOG(warn) << "*** Unless you have run the simulation event-by-event, it will not produce sensible results.";
+ }
+
// --- Initialise digi manager
CbmDigiManager* digiMan = CbmDigiManager::Instance();
digiMan->Init();
+
// --- Check input branch (PsdDigi). If not present, set task inactive.
if (!digiMan->IsPresent(ECbmModuleId::kPsd)) {
LOG(error) << GetName() << ": No PsdDigi input array present; "
<< "task will be inactive.";
return kERROR;
}
+ LOG(info) << GetName() << ": found PsdDigi branch";
-
- // Create and register output array
+ // --- Create and register output array
fHitArray = new TClonesArray("CbmPsdHit", 1000);
- ioman->Register(
- "PsdHit", "PSD", fHitArray, IsOutputBranchPersistent("PsdHit"));
-
+ ioman->Register("PsdHit", "PSD", fHitArray, IsOutputBranchPersistent("PsdHit"));
fHitArray->Dump();
- cout << "-I- CbmPsdHitProducer: Intialisation successfull " << kSUCCESS
- << endl;
+ LOG(info) << GetName() << ": registered branch PsdHit";
+
+ LOG(info) << GetName() << ": intialisation successfull";
+ LOG(info) << "======================================================\n";
return kSUCCESS;
}
-
-
// -------------------------------------------------------------------------
// ----- Public method Exec --------------------------------------------
-void CbmPsdHitProducer::Exec(Option_t* /*opt*/) {
+void CbmPsdHitProducer::Exec(Option_t* /*opt*/)
+{
- cout << " CbmPsdHitProducer::Exec(Option_t* /*opt*/) " << endl;
- fhModXNewEn->Print();
+ // --- Timer
+ TStopwatch timer;
+ timer.Start();
- // Reset output array
+ // --- Reset output array
Reset();
- // Declare some variables
- const CbmPsdDigi* dig = NULL;
+ // --- Local variables
+ pair<Int_t, Int_t> result;
+ Int_t nEvents = 0;
+ Int_t nDigis = 0;
+ Int_t nHits = 0;
+ Int_t nDigisAll = fDigiMan->GetNofDigis(ECbmModuleId::kPsd);
+ CbmEvent* event = nullptr;
+
+ // --- Time-slice mode: process entire digi array
+ if (!fEvents) {
+ result = ProcessData(nullptr);
+ nDigis = result.first;
+ nHits = result.second;
+ }
+
+ // --- Event mode: loop over events and process digis within
+ else {
+ assert(fEvents);
+ nEvents = fEvents->GetEntriesFast();
+ for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
+ event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
+ assert(event);
+ result = ProcessData(event);
+ LOG(debug) << GetName() << ": Event " << iEvent << " / " << nEvents << ", digis " << result.first << ", hits "
+ << result.second;
+ nDigis += result.first;
+ nHits += result.second;
+ } //# events
+ } //? event mode
+
+ // --- Timeslice log and statistics
+ timer.Stop();
+ if (fEvents)
+ LOG(info) << std::setw(20) << std::left << GetName() << fixed << setprecision(2) << " [" << timer.RealTime() * 1000.
+ << " ms] TS " << fNofTs << ", events " << nEvents << ", digis " << nDigis << " / " << nDigisAll
+ << ", hits " << fNofHits;
+ else
+ LOG(info) << std::setw(20) << std::left << GetName() << fixed << setprecision(2) << " [" << timer.RealTime() * 1000.
+ << " ms] TS " << fNofTs << ", digis " << nDigis << " / " << nDigisAll << ", hits " << fNofHits;
+ fNofTs++;
+ fNofEvents += nEvents;
+ fNofDigis += nDigis;
+ fNofHits += nHits;
+ fTimeTot += timer.RealTime();
+}
+// -------------------------------------------------------------------------
+
+
+// ----- End-of-timeslice action ---------------------------------------
+void CbmPsdHitProducer::Finish()
+{
+
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Time slices : " << fNofTs;
+ LOG(info) << "Digis / TS : " << fixed << setprecision(2) << Double_t(fNofDigis) / Double_t(fNofTs);
+ LOG(info) << "Hits / TS : " << fixed << setprecision(2) << Double_t(fNofHits) / Double_t(fNofTs);
+ LOG(info) << "Time / TS : " << fixed << setprecision(2) << 1000. * fTimeTot / Double_t(fNofTs) << " ms";
+ if (fEvents) {
+ LOG(info) << "Events : " << fNofEvents;
+ LOG(info) << "Events / TS : " << fixed << setprecision(2) << Double_t(fNofEvents) / Double_t(fNofTs);
+ LOG(info) << "Digis / event : " << fixed << setprecision(2) << Double_t(fNofDigis) / Double_t(fNofEvents);
+ LOG(info) << "Hits / event : " << fixed << setprecision(2) << Double_t(fNofHits) / Double_t(fNofEvents);
+ }
+
+ // --- Write energy deposition histogram
+ TFile* oldFile = gFile;
+ TDirectory* oldDir = gDirectory;
+ TFile* outfile = new TFile("EdepHistos.root", "RECREATE");
+ outfile->cd();
+ fhModXNewEn->Write();
+ outfile->Close();
+ gFile = oldFile;
+ gDirectory = oldDir;
+ LOG(info) << "Histograms written to EdepHistos.root";
+ LOG(info) << "=====================================\n";
+}
+// -------------------------------------------------------------------------
+
- /*
- Double_t edep[NPsdMod];//marina
- for (Int_t imod=0; imod<NPsdMod; imod++) { edep[imod]=0.; }//marina
-*/
+// ----- Data processing -----------------------------------------------
+pair<Int_t, Int_t> CbmPsdHitProducer::ProcessData(CbmEvent* event)
+{
- std::map<int, Double_t> edepmap;
+ // --- Local variables
+ std::map<int, Double_t> edepmap; // energy deposition per module, key is module ID
+ Int_t nHits = 0;
+ const CbmPsdDigi* digi = nullptr;
+ Int_t nDigis = (event ? event->GetNofData(ECbmDataType::kPsdDigi) : fDigiMan->GetNofDigis(ECbmModuleId::kPsd));
- // Loop over PsdDigits
- Int_t nDigi = fDigiMan->GetNofDigis(ECbmModuleId::kPsd);
- cout << " nDigits " << nDigi << endl;
+ // --- Loop over PsdDigis. Distribute energy deposition into the modules.
+ Int_t digiIndex = -1;
+ for (Int_t iDigi = 0; iDigi < nDigis; iDigi++) {
- for (Int_t idig = 0; idig < nDigi; idig++) {
- dig = fDigiMan->Get<CbmPsdDigi>(idig);
- if (!dig) continue;
- Int_t mod = dig->GetModuleID();
- // Int_t sec = dig->GetSectionID();
- Double_t eDep = (Double_t) dig->GetEdep();
- //edep[mod-1] += (Double_t) dig->GetEdep(); //DEBUG: SELIM
+ digiIndex = (event ? event->GetIndex(ECbmDataType::kPsdDigi, iDigi) : iDigi);
+ digi = fDigiMan->Get<const CbmPsdDigi>(digiIndex);
+ assert(digi);
+ Int_t mod = digi->GetModuleID();
+ Double_t eDep = (Double_t) digi->GetEdep();
auto insert_result = edepmap.insert(std::make_pair(mod, eDep));
if (!insert_result.second) { // entry was here before
(*insert_result.first).second += eDep;
}
- } // Loop over MCPoints
+ } // # digis
- fNHits = 0;
+ // --- Create hits, one per activated module
+ UInt_t hitIndex = -1;
for (auto edep_entry : edepmap) {
int modID = edep_entry.first;
Double_t eDep = edep_entry.second;
- new ((*fHitArray)[fNHits]) CbmPsdHit(modID, eDep);
- fNHits++;
- }
-
- /*
- for (Int_t imod=0; imod<NPsdMod; imod++) //marina
- {
- if (edep[imod]>0.)
- {
- new ((*fHitArray)[fNHits]) CbmPsdHit(imod+1, edep[imod]);
- fNHits++;
- //cout<<"MARINA CbmPsdHitProducer " <<fNHits <<" " <<imod+1 <<" " <<edep[imod] <<endl;
- //fhModXNewEn->Fill(fXi[imod],TMath::Sqrt(edep[imod]) );
- //cout<<"CbmPsdHitProducer "<<fNHits<<" "<<imod<<" "<<edep[imod]<<endl;
- }
+ hitIndex = fHitArray->GetEntriesFast();
+ new ((*fHitArray)[hitIndex]) CbmPsdHit(modID, eDep);
+ if (event) event->AddData(ECbmDataType::kPsdHit, hitIndex);
+ nHits++;
}
- */
- // Event summary
- cout << "-I- CbmPsdHitProducer: " << fNHits << " CbmPsdHits created." << endl;
+ return std::make_pair(nDigis, nHits);
}
// -------------------------------------------------------------------------
-void CbmPsdHitProducer::Finish() {
- cout << " CbmPsdHitProducer::Finish() " << endl;
- /// Save old global file and folder pointer to avoid messing with FairRoot
- TFile* oldFile = gFile;
- TDirectory* oldDir = gDirectory;
-
- TFile* outfile = new TFile("EdepHistos.root", "RECREATE");
- outfile->cd();
- fhModXNewEn->Write();
- /// Restore old global file and folder pointer to avoid messing with FairRoot
- gFile = oldFile;
- gDirectory = oldDir;
-
- //outfile->Close(); //SELIM
-}
-
-// ----- Private method Reset ------------------------------------------
-void CbmPsdHitProducer::Reset() {
- fNHits = 0;
+// ----- Reset: clear output array -------------------------------------
+void CbmPsdHitProducer::Reset()
+{
if (fHitArray) fHitArray->Delete();
}
+// -------------------------------------------------------------------------
ClassImp(CbmPsdHitProducer)
diff --git a/reco/detectors/psd/CbmPsdHitProducer.h b/reco/detectors/psd/CbmPsdHitProducer.h
index 05fd5a5ae292cb9ba9e06a9f7e2f3581c5539ffe..03b4707d7a6bd6ceb149820c7109c4d7ec12b438 100644
--- a/reco/detectors/psd/CbmPsdHitProducer.h
+++ b/reco/detectors/psd/CbmPsdHitProducer.h
@@ -17,10 +17,13 @@
#include "CbmDigiManager.h"
-#include "FairTask.h"
-#include "TH1F.h"
+
+#include <FairTask.h>
+
+#include <TH1F.h>
class TClonesArray;
+class CbmEvent;
const Int_t NPsdMod = 44; //with 4 central mods
@@ -45,13 +48,28 @@ public:
private:
- Int_t fNHits;
+ Int_t fNHits = 0;
/** Output array of CbmPsdHit **/
- TClonesArray* fHitArray;
+ TClonesArray* fHitArray = nullptr;
/** Digi Manager for input **/
- CbmDigiManager* fDigiMan; //!
+ CbmDigiManager* fDigiMan = nullptr; //!
+
+ /** Event array **/
+ TClonesArray* fEvents = nullptr; //!
+
+
+ /** @brief Processing of digis **/
+ std::pair<Int_t, Int_t> ProcessData(CbmEvent* event);
+
+ /** Counters **/
+ ULong64_t fNofTs = 0;
+ ULong64_t fNofEvents = 0;
+ ULong64_t fNofDigis = 0;
+ ULong64_t fNofHits = 0;
+ Double_t fTimeTot = 0.;
+
CbmPsdHitProducer(const CbmPsdHitProducer&);
CbmPsdHitProducer operator=(const CbmPsdHitProducer&);
@@ -61,7 +79,7 @@ private:
Float_t fXi[NPsdMod]; //X coordinate of center of module
Float_t fYi[NPsdMod]; //X coordinate of center of module
- TH1F* fhModXNewEn; //edep in each module for Marina
+ TH1F* fhModXNewEn = nullptr; //edep in each module for Marina
ClassDef(CbmPsdHitProducer, 2);