123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238 |
- //------------------------------------------------------------------------------------------------------------------------
- /// \class MpdTofHitProducer
- ///
- /// \brief
- /// \author Sergei Lobastov (LHE, JINR, Dubna)
- //------------------------------------------------------------------------------------------------------------------------
- #include <assert.h>
- #include <TMath.h>
- #include <TH1D.h>
- #include <TH2D.h>
- #include <TEfficiency.h>
- #include <TRandom2.h>
- #include <TClonesArray.h>
- #include "MpdMCTrack.h"
- #include "FairLogger.h"
- #include "MpdTofUtils.h"
- #include "MpdTofPoint.h"
- #include "MpdTofHit.h"
- #include "MpdTofGeoUtils.h"
- #include "MpdTofHitProducerQA.h"
- #include "MpdTofHitProducer.h"
- using namespace std;
- ClassImp(MpdTofHitProducer)
- //------------------------------------------------------------------------------------------------------------------------
- MpdTofHitProducer::MpdTofHitProducer(const char *name, Bool_t useMCdata, Int_t verbose, Bool_t test, const char *flnm, double stripLength)
- : MpdTofHitProducerIdeal(name, useMCdata, verbose, test, true, flnm, false), fStripLength(stripLength)
- {
- pRandom = new TRandom2;
- }
- //------------------------------------------------------------------------------------------------------------------------
- MpdTofHitProducer::~MpdTofHitProducer()
- {
- delete pRandom;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdTofHitProducer::AddParameters(TString& buf)const
- {
- MpdTofHitProducerIdeal::AddParameters(buf);
- buf += TString::Format(", fTimeSigma=%.5g ns", fTimeSigma);
- buf += TString::Format(", fErrX=%.4g cm", fErrX);
- buf += TString::Format(", fErrZ=%.4g cm", fErrZ);
- }
- //------------------------------------------------------------------------------------------------------------------------
- InitStatus MpdTofHitProducer::Init()
- {
- InitStatus status = MpdTofHitProducerIdeal::Init();
- MpdTofGeoUtils::Instance()->ParseTGeoManager(pQA, true);//, "geometryTree.txt"); // forced
-
- MpdTofGeoUtils::Instance()->FindNeighborStrips(0.8, pQA, // 0.8 [cm] <--- thresh. distance between neighbor strips
- true); // forced
- LOG(INFO)<<"[MpdTofHitProducer::Init] Initialization finished succesfully.";
- return status;
- }
- //------------------------------------------------------------------------------------------------------------------------
- Bool_t MpdTofHitProducer::IsHitCreated(Double_t value, Int_t gap) // value - distance to the strip edge [cm]
- {
- constexpr Double_t slope = (0.98 - 0.95)/0.2;
- Double_t efficiency = (value > 0.2) ? 0.98 : ( 0.95 + slope*value);
-
- //-------------------------------------
- // 99% ---------
- // \.
- // \.
- // \.
- // 95% \.
- // <-----------|--|
- // 0.2 0. value
- //-------------------------------------
- if(gap == 1 || gap == 3) efficiency /= 2.; // reduce efficiency on outer strip gap
- if(pRandom->Rndm() < efficiency) return true;
- return false;
- }
- //------------------------------------------------------------------------------------------------------------------------
- Bool_t MpdTofHitProducer::IsCrossHitCreated(Double_t value, Int_t gap) // value - distance to the strip edge [cm]
- {
- constexpr Double_t slope = (0.3 - 0.0)/0.5;
- Double_t efficiency = (value > 0.5) ? 0. : (0.3 - slope*value);
-
- //-------------------------------------
- // 30% /
- // /
- // /
- // /
- // 0% /
- // <-----------|----|
- // 0.5 0. value
- //-------------------------------------
-
- if(efficiency == 0.) return false;
- if(gap == 1 || gap == 3) efficiency /= 2.; // reduce efficiency on outer strip gap
- if(pRandom->Rndm() < efficiency) return true;
- return false;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdTofHitProducer::Exec(Option_t *option)
- {
- aTofHits->Clear();
- TVector3 mcPosition, hitPosition, hitPosError;
- int nSingleHits = 0, nCrossHits = 0;
- auto SmearingAlongStrip = [this](const TGeoCombiTrans& gap, const TGeoCombiTrans& centralGap, const TVector3& point, TVector3& hit, TVector3& error)
- {
- // rotate stip to origin LRS (dX = fStripLength, dY = 0, dZ = strip step)
- Double_t local[3], master[3] = {point.X(), point.Y(), point.Z()};
- gap.MasterToLocal(master, local);
- // smearing the hit position along x axis
- double Xsmeared, x = local[0]; size_t nTries = 0;
- do
- {
- Xsmeared = pRandom->Gaus(x, fErrX);
- if(++nTries > 100) return false;
- }
- while(fabs(Xsmeared) > fStripLength/2.); // check strip boundary exit
- // shift and rotate back to MRS(shift by Y from gap to central gap)
- local[0] = Xsmeared, local[2] = 0.; // smearing along X axis, shift to the strip center along Z axis(at LRS)
- centralGap.LocalToMaster(local, master);
- hit.SetXYZ(master[0], master[1], master[2]);
- // rotate error vector
- local[0] = fErrX; local[1] = 0.; local[2] = fErrZ;
- centralGap.LocalToMaster(local, master);
- error.SetXYZ(master[0], master[1], master[2]);
- return true;
- };
- if(fUseMCData)
- {
- for(Int_t pointIndex = 0, nTofPoint = aMcPoints->GetEntriesFast(); pointIndex < nTofPoint; pointIndex++ ) // cycle by TOF points
- {
- auto pPoint = (MpdTofPoint*) aMcPoints->UncheckedAt(pointIndex);
-
- if(fVerbose > 2) pPoint->Print("");
-
- Int_t tid = pPoint->GetTrackID();
- Int_t suid = pPoint->GetDetectorID();
- Int_t gap = pPoint->GetGap();
- Double_t time = pRandom->Gaus(pPoint->GetTime(), fTimeSigma); // default 100 ps
- pPoint->Position(mcPosition);
- auto strip = MpdTofGeoUtils::Instance()->FindStrip(suid);
- auto centralStrip = MpdTofGeoUtils::Instance()->FindStrip(MpdTofPoint::SetCentralGap(suid));
- if(! SmearingAlongStrip(strip->fMatrix, centralStrip->fMatrix, mcPosition, hitPosition, hitPosError))
- {
- strip->Dump(" [MpdTofHitProducer::Exec] -E- Invalid Rotation matrix.");
- continue;
- }
- LStrip::Side_t side;
- Double_t distance = strip->MinDistanceToEdge(&mcPosition, side); // [cm]
- bool stripFired;
- if( (stripFired = IsHitCreated(distance, gap)) ) // simulate hit efficiency (add flag MpdTofUtils::IsSingle)
- {
- AddHit(MpdTofPoint::ClearGap(suid), hitPosition, hitPosError, pointIndex, tid, time, MpdTofUtils::IsSingle);
- nSingleHits++;
- }
-
- if(pQA)
- {
- if(stripFired) pQA->Point2HitSmearingTest(mcPosition, hitPosition);
- pQA->HitGapEfficiencyTest(stripFired, distance, gap);
- pQA->PositionInsideStripTest(strip->center, mcPosition, hitPosition);
- pQA->RotationToOriginTest(strip->fMatrix, mcPosition, hitPosition);
- pQA->CentralDetectorTest(suid, hitPosition);
- }
-
- if( (stripFired = IsCrossHitCreated(distance, gap)) ) // simulate cross hit efficiency (add flag MpdTofUtils::IsDouble)
- {
- Int_t crossSuid = (side == LStrip::kRight) ? strip->neighboring[LStrip::kRight] : strip->neighboring[LStrip::kLeft];
-
- if(LStrip::kInvalid == crossSuid) continue; // edge strip on detector
-
- strip = MpdTofGeoUtils::Instance()->FindStrip(crossSuid);
- centralStrip = MpdTofGeoUtils::Instance()->FindStrip(MpdTofPoint::SetCentralGap(crossSuid));
- if(! SmearingAlongStrip(strip->fMatrix, centralStrip->fMatrix, mcPosition, hitPosition, hitPosError))
- {
- strip->Dump(" [MpdTofHitProducer::Exec] -E- Invalid Rotation matrix.");
- continue;
- }
- AddHit(MpdTofPoint::ClearGap(crossSuid), hitPosition, hitPosError, pointIndex, tid, time, MpdTofUtils::IsDouble);
- nCrossHits++;
- }
- if(pQA) pQA->FillCrossHitEfficiency(stripFired, distance, gap, mcPosition, hitPosition);
- } // cycle by the TOF points
-
- }
- else // exp. data used
- {
- // FIXME: now not realized
- //AddHit(Int_t detUID, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
- assert(false);
- }
-
- MergeHitsOnStrip(); // save only the fastest hit in the strip
- int Nhits = CompressHits(); // remove blank slotes
- if(fUseMCData && pQA)
- {
- pQA->FillNPointsHits(aMcPoints->GetEntriesFast(), Nhits);
- pQA->PointDistanceTest(aMcPoints);
- }
- LOG(DEBUG1)<<"[MpdTofHitProducer::Exec] single hits= "<<nSingleHits<<", cross hits= "<<nCrossHits<<", final hits= "<<Nhits;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdTofHitProducer::Finish()
- {
- if(pQA) pQA->Finish();
- }
- //--------------------------------------------------------------------------------------------------------------------------------------
- void MpdTofHitProducer::SetSeed(UInt_t seed)
- {
- pRandom->SetSeed(seed);
- }
- //------------------------------------------------------------------------------------------------------------------------
|