// ------------------------------------------------------------------------- // ----- PndMdtDigitization source file ----- // ----- author: Jifeng Hu, hu@to.infn.it // ------------------------------------------------------------------------- #include "PndMdtDigitization.h" #include "PndMdtWaveform.h" #include "PndMdtDigi.h" #include "FairRootManager.h" #include "FairRun.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "TClonesArray.h" #include "TVector3.h" #include "TRandom.h" #include #include "PndMdtIGeometry.h" #include "FairEventHeader.h" using std::map; using std::cout; using std::endl; // ----- Default constructor ------------------------------------------- PndMdtDigitization::PndMdtDigitization(Int_t verbose, Bool_t store) : FairTask(" MDT Digitization") , fPersitence(store), fTimeOrderedDigi(kFALSE) { SetVerbose(verbose); } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- PndMdtDigitization::~PndMdtDigitization() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus PndMdtDigitization::Init() { //FairRun* sim = FairRun::Instance(); //[R.K. 01/2017] unused variable? //FairRuntimeDb* rtdb=sim->GetRuntimeDb(); //[R.K. 01/2017] unused variable? // Get RootManager FairRootManager* ioman = FairRootManager::Instance(); if ( ! ioman ) { cout << "-E- PndMdtDigitization::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // Get input array fWaveformArray = (TClonesArray*) ioman->GetObject("MdtWaveform"); if ( ! fWaveformArray ) { cout << "-W- PndMdtDigitization::Init: " << "No MdtWaveform array!" << endl; return kERROR; } if(fTimeOrderedDigi){ fDigiArray = ioman->Register("MdtDigi","PndMdtDigi","Mdt", fPersitence); }else{ // Create and register output array fDigiBoxArray = ioman->Register("MdtDigiBox","PndMdtDigi","Mdt", fPersitence); fDigiStripArray = ioman->Register("MdtDigiStrip","PndMdtDigi","Mdt", fPersitence); } fWireNoiseSigma = 0.05; fStripNoiseSigma = 0.01; fSamplingInterval = 10.;//nano seconds fGeoIF = PndMdtIGeometry::Instance(); fGeoIF->SetVerbose(3); fGeoIF->AddSensor("MDT"); fGeoIF->AddSensor("GasCell"); Bool_t fGoodGeo = fGeoIF->Init(); HowManyWireWf = HowManyStripWf = 0; HowManyWireDigi = HowManyStripDigi = 0; cout << "-I- PndMdtDigitization: Intialization "<<(fGoodGeo ? "successful." : "failed.") << endl; return fGoodGeo ? kSUCCESS : kERROR; } // ------------------------------------------------------------------------- // ----- Public method Exec -------------------------------------------- void PndMdtDigitization::Exec(Option_t*) { if(fTimeOrderedDigi) exec_t(); else exec_e(); } void PndMdtDigitization::exec_e() { // Reset output array fDigiBoxArray->Delete(); fDigiStripArray->Delete(); Int_t evtNo = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber(); Int_t nWaveform = fWaveformArray->GetEntriesFast(); PndMdtWaveform* theWf(0); Double_t fTimeStamp; Double_t fAmplitude; for (Int_t iW=0; iWAt(iW); if(theWf->IsWire()){ ++HowManyWireWf; if(Digitize(theWf, fTimeStamp, fAmplitude, kTRUE)){ TVector3 fPos; fGeoIF->GetTubeCenter(theWf->GetDetectorID(), fPos); PndMdtDigi* aDigi = AddDigiBox(theWf->GetDetectorID(), fPos, evtNo); aDigi->SetTimeStamp(theWf->GetTimeStamp() + fTimeStamp); ++HowManyWireDigi; } }else{ ++HowManyStripWf; if(Digitize(theWf, fTimeStamp, fAmplitude, kFALSE)){ TVector3 fPos; fGeoIF->GetStripCenter(theWf->GetDetectorID(), fPos); PndMdtDigi* aDigi = AddDigiStrip(theWf->GetDetectorID(), fPos, evtNo); aDigi->SetTimeStamp(theWf->GetTimeStamp() + fTimeStamp); ++HowManyStripDigi; } } } } void PndMdtDigitization::exec_t() { // Reset output array fDigiArray->Delete(); Int_t evtNo = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber(); Int_t nWaveform = fWaveformArray->GetEntriesFast(); PndMdtWaveform* theWf(0); Double_t fTimeStamp; Double_t fAmplitude; // cout<<"event no "<At(iW); //theWf->print(); if(theWf->IsWire()){ ++HowManyWireWf; if(Digitize(theWf, fTimeStamp, fAmplitude, kTRUE)){ TVector3 fPos; fGeoIF->GetTubeCenter(theWf->GetDetectorID(), fPos); PndMdtDigi* aDigi = AddDigi(theWf->GetDetectorID(), fPos, evtNo); aDigi->SetTimeStamp(theWf->GetTimeStamp() + fTimeStamp); ++HowManyWireDigi; } }else{ ++HowManyStripWf; if(Digitize(theWf, fTimeStamp, fAmplitude, kFALSE)){ TVector3 fPos; fGeoIF->GetStripCenter(theWf->GetDetectorID(), fPos); PndMdtDigi* aDigi = AddDigi(theWf->GetDetectorID(), fPos, evtNo); aDigi->SetTimeStamp(theWf->GetTimeStamp() + fTimeStamp); ++HowManyStripDigi; } } } } // ------------------------------------------------------------------------- // ----- Private method AddDigi -------------------------------------------- PndMdtDigi* PndMdtDigitization::AddDigiBox(Int_t detID, TVector3& pos, Int_t evtNo) { // It fills the PndMdtDigi category TClonesArray& clref = *fDigiBoxArray; Int_t size = clref.GetEntriesFast(); return new(clref[size]) PndMdtDigi(detID, pos, evtNo); } // ---- // ----- Private method AddDigi -------------------------------------------- PndMdtDigi* PndMdtDigitization::AddDigiStrip(Int_t detID, TVector3& pos, Int_t evtNo) { // It fills the PndMdtDigi category TClonesArray& clref = *fDigiStripArray; Int_t size = clref.GetEntriesFast(); return new(clref[size]) PndMdtDigi(detID, pos, evtNo); } PndMdtDigi* PndMdtDigitization::AddDigi(Int_t detID, TVector3& pos, Int_t evtNo) { // It fills the PndMdtDigi category TClonesArray& clref = *fDigiArray; Int_t size = clref.GetEntriesFast(); return new(clref[size]) PndMdtDigi(detID, pos, evtNo); } // --- threshold digitization --- Bool_t PndMdtDigitization::Digitize(PndMdtWaveform* theWf, Double_t& time, Double_t& amp, Bool_t isWire) { time = -1; amp = -9999; Double_t fPeak = 0.; Bool_t NotFound = kTRUE; //Double_t fNoiseSigma; //[R.K. 01/2017] unused variable? Double_t fThreshold(1e9); if(isWire) fThreshold = 5.*fWireNoiseSigma; else fThreshold = 5.*fStripNoiseSigma; const std::vector& fSignalData = theWf->GetSignal(); for(size_t is=0; is < fSignalData.size(); ++ is){ if(TMath::Abs(fSignalData[is]) > TMath::Abs(fPeak)) fPeak = fSignalData[is]; if(TMath::Abs(fSignalData[is]) > fThreshold && NotFound) { time = is*fSamplingInterval;//nano seconds NotFound = kFALSE; } } amp = fPeak; return time > 0; } // ---- void PndMdtDigitization::FinishTask() { if(fTimeOrderedDigi){ Int_t nWaveform = fWaveformArray->GetEntriesFast(); if(fVerbose>0){ std::cout<<"PndMdtDigitization::FinishTask, fWaveformArray size #"<