From 591afa9672d585d548d0d86ef6d956e3bf6fdaf2 Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Mon, 8 Apr 2024 16:49:16 -0500 Subject: [PATCH 01/11] StarDb added --- StRoot/StAnalysisMaker/StAnalysisMaker.cxx | 7 +- StRoot/StBFChain/BigFullChain.h | 13 + StRoot/StBFChain/StBFChain.cxx | 5 + StRoot/StEvent/StContainers.cxx | 6 + StRoot/StEvent/StContainers.h | 6 + StRoot/StEvent/StEvent.cxx | 25 + StRoot/StEvent/StEvent.h | 4 + StRoot/StEvent/StEventClusteringHints.cxx | 2 + StRoot/StEvent/StEventTypes.h | 4 + StRoot/StEvent/StGmtCollection.cxx | 58 + StRoot/StEvent/StGmtCollection.h | 90 ++ StRoot/StEvent/StGmtHit.cxx | 52 + StRoot/StEvent/StGmtHit.h | 136 ++ StRoot/StEvent/StGmtHitCollection.cxx | 40 + StRoot/StEvent/StGmtHitCollection.h | 61 + StRoot/StEvent/StGmtPoint.cxx | 65 + StRoot/StEvent/StGmtPoint.h | 67 + StRoot/StEvent/StGmtPointCollection.cxx | 30 + StRoot/StEvent/StGmtPointCollection.h | 52 + StRoot/StEvent/StGmtStrip.cxx | 107 ++ StRoot/StEvent/StGmtStrip.h | 209 +++ StRoot/StEvent/StGmtStripCollection.cxx | 160 +++ StRoot/StEvent/StGmtStripCollection.h | 155 +++ StRoot/StGmtAligner/EventT.cxx | 551 ++++++++ StRoot/StGmtAligner/EventT.h | 141 ++ StRoot/StGmtAligner/HitT.h | 160 +++ StRoot/StGmtAligner/StGmtAligner.cxx | 150 +++ StRoot/StGmtAligner/StGmtAligner.h | 36 + StRoot/StGmtAligner/TrackT.h | 78 ++ StRoot/StGmtAligner/runGmtTree.C | 23 + .../StGmtClusterMaker/StGmtClusterMaker.cxx | 302 +++++ StRoot/StGmtClusterMaker/StGmtClusterMaker.h | 44 + StRoot/StGmtRawMaker/StGmtRawMaker.cxx | 189 +++ StRoot/StGmtRawMaker/StGmtRawMaker.h | 29 + StRoot/StGmtUtil/StGmtConsts.h | 60 + StRoot/StGmtUtil/geometry/StGmtGeom.cxx | 1131 +++++++++++++++++ StRoot/StGmtUtil/geometry/StGmtGeom.h | 270 ++++ .../gmt/GmtOnModule.20140101.000005.C | 21 + StarDb/Geometry/gmt/GmtOnModule.C | 31 + StarDb/Geometry/gmt/GmtOnTpc.C | 40 + 40 files changed, 4607 insertions(+), 3 deletions(-) create mode 100644 StRoot/StEvent/StGmtCollection.cxx create mode 100644 StRoot/StEvent/StGmtCollection.h create mode 100644 StRoot/StEvent/StGmtHit.cxx create mode 100644 StRoot/StEvent/StGmtHit.h create mode 100644 StRoot/StEvent/StGmtHitCollection.cxx create mode 100644 StRoot/StEvent/StGmtHitCollection.h create mode 100644 StRoot/StEvent/StGmtPoint.cxx create mode 100644 StRoot/StEvent/StGmtPoint.h create mode 100644 StRoot/StEvent/StGmtPointCollection.cxx create mode 100644 StRoot/StEvent/StGmtPointCollection.h create mode 100644 StRoot/StEvent/StGmtStrip.cxx create mode 100644 StRoot/StEvent/StGmtStrip.h create mode 100644 StRoot/StEvent/StGmtStripCollection.cxx create mode 100644 StRoot/StEvent/StGmtStripCollection.h create mode 100644 StRoot/StGmtAligner/EventT.cxx create mode 100644 StRoot/StGmtAligner/EventT.h create mode 100644 StRoot/StGmtAligner/HitT.h create mode 100644 StRoot/StGmtAligner/StGmtAligner.cxx create mode 100644 StRoot/StGmtAligner/StGmtAligner.h create mode 100644 StRoot/StGmtAligner/TrackT.h create mode 100644 StRoot/StGmtAligner/runGmtTree.C create mode 100644 StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx create mode 100644 StRoot/StGmtClusterMaker/StGmtClusterMaker.h create mode 100644 StRoot/StGmtRawMaker/StGmtRawMaker.cxx create mode 100644 StRoot/StGmtRawMaker/StGmtRawMaker.h create mode 100644 StRoot/StGmtUtil/StGmtConsts.h create mode 100644 StRoot/StGmtUtil/geometry/StGmtGeom.cxx create mode 100644 StRoot/StGmtUtil/geometry/StGmtGeom.h create mode 100644 StarDb/Geometry/gmt/GmtOnModule.20140101.000005.C create mode 100644 StarDb/Geometry/gmt/GmtOnModule.C create mode 100644 StarDb/Geometry/gmt/GmtOnTpc.C diff --git a/StRoot/StAnalysisMaker/StAnalysisMaker.cxx b/StRoot/StAnalysisMaker/StAnalysisMaker.cxx index ee97767ca92..5c5ccb1e380 100644 --- a/StRoot/StAnalysisMaker/StAnalysisMaker.cxx +++ b/StRoot/StAnalysisMaker/StAnalysisMaker.cxx @@ -1141,13 +1141,14 @@ void StAnalysisMaker::summarizeEvent(StEvent *event, Int_t mEventCounter) { if (event->numberOfPsds()) { LOG_QA << "# PSDs: " << event->numberOfPsds() << endm; } -#ifdef _ST_GMT_HIT_H_ +//#ifdef _ST_GMT_HIT_H_ if (event->gmtCollection() && event->gmtCollection()->getNumHits()) { LOG_QA << "# GMT hits: " << event->gmtCollection()->getNumHits() - << " points: " << event->gmtCollection()->getNumPoints() + << ", points: " << event->gmtCollection()->getNumPoints() + << ", strips: " << event->gmtCollection()->getNumStrips() << endm; } -#endif /* _ST_GMT_HIT_H_ */ +//#endif /* _ST_GMT_HIT_H_ */ if (event->rpsCollection()) { LOG_QA << "# RPS hits: " << event->rpsCollection()->clusters().size() << endm; } diff --git a/StRoot/StBFChain/BigFullChain.h b/StRoot/StBFChain/BigFullChain.h index 917dd831eaa..dd1c84552a8 100644 --- a/StRoot/StBFChain/BigFullChain.h +++ b/StRoot/StBFChain/BigFullChain.h @@ -1381,6 +1381,13 @@ Bfc_st BFC[] = { // standard chains , "MTD hit maker",kFALSE}, {"mtdTrkMask","","","db","StMtdTrackingMaskMaker" ,"StMtdEvtFilterMaker","MTD track masking",kFALSE}, + // GMT + {"gmt" ,"","","gmtDat,gmtClu" ,"","","Gmt data Chain",kFALSE}, + {"gmtDat" ,"","","event","StGmtRawMaker","StGmtRawMaker", "GMT Data reader",kFALSE}, + {"gmtClu" ,"","","gmtutil","StGmtClusterMaker","Spectrum,StGmtClusterMaker","GMT cluster maker",kFALSE}, + {"gmtCosmics" ,"","","Cosmics,gmt","","", "Save only events with GMT clusters and Cosmic tracks",kFALSE}, + {"gmtClusTree","","","","","", "WARNING *** Option is OBSOLETE ***",kFALSE}, + // EPD {"epdHit", "", "", "epdDb,event", "StEpdHitMaker", "StEpdHitMaker","EPD hit maker", kFALSE}, @@ -1782,6 +1789,10 @@ Bfc_st BFC[] = { // standard chains {"mtdCalib" ,"","","db", "StMtdCalibMaker","StMtdCalibMaker","MTD calibration",kFALSE}, {"mtdEvtFilt" ,"","","db", "StMtdEvtFilterMaker","StMtdEvtFilterMaker","MTD event filter" ,kFALSE}, + // GMT + {"gmtMatch" ,"","","", "","","WARNING *** Option is OBSOLETE ***",kFALSE}, + {"gmtPlotter" ,"","","", "","","WARNING *** Option is OBSOLETE ***",kFALSE}, + {"FindVtxSeed" ,"FindVtxSeed" ,"","globT,MuDSTDeps,picoDst","StVertexSeedMaker" , "StPass0CalibMaker", "Performs vertex seed finding",kFALSE}, {"FindEvtVtxSeed","FindEvtVtxSeed","","MuDSTDeps,picoDst","StEvtVtxSeedMaker" @@ -1913,6 +1924,8 @@ Bfc_st BFC[] = { // standard chains , "StSvtPoolEventT,StSvtPoolSvtMatchedTree","Create SvtMatchedTree",kFALSE}, {"LAna" ,"","","in,detDb,StEvent,tpcDb","StLaserAnalysisMaker" , "StLaserAnalysisMaker","Laser data Analysis",kFALSE}, + // GMT + {"gmtAligner" ,"","","detDb", "StGmtAligner","StGmtAligner","GMT cluster plotting",kFALSE}, {"EandBDir","","","in,StEvent,TpcHitMover,nodefault" , "StEandBDirMaker","MathMore,Spectrum,StEandBDirMaker", "get E&B direction",kFALSE}, {"SpinTag" ,"","","" ,"","",STAR_CHAIN_OBSOLETE,kFALSE}, diff --git a/StRoot/StBFChain/StBFChain.cxx b/StRoot/StBFChain/StBFChain.cxx index 566db07be1f..3bd01c1f138 100644 --- a/StRoot/StBFChain/StBFChain.cxx +++ b/StRoot/StBFChain/StBFChain.cxx @@ -550,6 +550,11 @@ Int_t StBFChain::Instantiate() if (GetOption("Cosmics")) mk->SetAttr("Cosmics" ,kTRUE); mk->PrintAttr(); } + + if (maker== "StGmtClusterMaker") { + if (GetOption("gmtCosmics")) mk->SetAttr("gmtCosmics", kTRUE); + } + if (maker=="StKFVertexMaker" && GetOption("laserIT")) mk->SetAttr("laserIT" ,kTRUE); // Sti(ITTF) end if (maker=="StGenericVertexMaker") { diff --git a/StRoot/StEvent/StContainers.cxx b/StRoot/StEvent/StContainers.cxx index 3afb846cf82..9e40ac0a166 100644 --- a/StRoot/StEvent/StContainers.cxx +++ b/StRoot/StEvent/StContainers.cxx @@ -198,6 +198,9 @@ #include "StTrackPidTraits.h" #include "StV0Vertex.h" #include "StXiVertex.h" +#include "StGmtHit.h" +#include "StGmtPoint.h" +#include "StGmtStrip.h" StCollectionImp(BTofHit) StCollectionImp(BTofRawHit) @@ -277,3 +280,6 @@ StCollectionImp(TrackNode) StCollectionImp(TrackPidTraits) StCollectionImp(V0Vertex) StCollectionImp(XiVertex) +StCollectionImp(GmtHit) +StCollectionImp(GmtStrip) +StCollectionImp(GmtPoint) diff --git a/StRoot/StEvent/StContainers.h b/StRoot/StEvent/StContainers.h index 67e5d60c83c..5bc3470e0ab 100644 --- a/StRoot/StEvent/StContainers.h +++ b/StRoot/StEvent/StContainers.h @@ -214,6 +214,9 @@ class StTrackNode; class StTrackPidTraits; class StV0Vertex; class StXiVertex; +class StGmtHit; +class StGmtStrip; +class StGmtPoint; StCollectionDef(BTofHit) StCollectionDef(BTofRawHit) @@ -293,6 +296,9 @@ StCollectionDef(TrackNode) StCollectionDef(TrackPidTraits) StCollectionDef(V0Vertex) StCollectionDef(XiVertex) +StCollectionDef(GmtHit) +StCollectionDef(GmtStrip) +StCollectionDef(GmtPoint) #endif diff --git a/StRoot/StEvent/StEvent.cxx b/StRoot/StEvent/StEvent.cxx index 2ffeb3913fd..f10d844a69b 100644 --- a/StRoot/StEvent/StEvent.cxx +++ b/StRoot/StEvent/StEvent.cxx @@ -256,6 +256,7 @@ #include "StTrackNode.h" #include "StTrack.h" #include "StFwdTrackCollection.h" +#include "StGmtCollection.h" #ifndef ST_NO_NAMESPACES using std::swap; @@ -908,6 +909,23 @@ StEvent::fgtCollection() const return fgtCollection; } +StGmtCollection* +StEvent::gmtCollection() +{ + StGmtCollection *gmtCollection = 0; + _lookup(gmtCollection, mContent); + return gmtCollection; +} + +const StGmtCollection* +StEvent::gmtCollection() const +{ + StGmtCollection *gmtCollection = 0; + _lookup(gmtCollection, mContent); + return gmtCollection; +} + + StIstHitCollection* StEvent::istHitCollection() { @@ -1456,6 +1474,12 @@ StEvent::setFgtCollection(StFgtCollection* val) _lookupAndSet(val, mContent); } +void +StEvent::setGmtCollection(StGmtCollection* val) +{ + _lookupAndSet(val, mContent); +} + void StEvent::setIstHitCollection(StIstHitCollection* val) { @@ -1629,6 +1653,7 @@ void StEvent::statistics() cout << "\t# of hits in EMC: " << (emcCollection() ? emcCollection()->barrelPoints().size() : 0) << endl; cout << "\t# of hits in EEMC: " << (emcCollection() ? emcCollection()->endcapPoints().size() : 0) << endl; cout << "\t# of hits in FGT: " << (fgtCollection() ? fgtCollection()->getNumHits() : 0) << endl; + cout << "\t# of hits in GMT: " << (gmtCollection() ? gmtCollection()->getNumHits() : 0) << endl; cout << "\t# of hits in RICH: " << (richCollection() ? richCollection()->getRichHits().size() : 0) << endl; cout << "\t# of PSDs: " << numberOfPsds() << endl; } diff --git a/StRoot/StEvent/StEvent.h b/StRoot/StEvent/StEvent.h index b0beda64953..6a667b5e19e 100644 --- a/StRoot/StEvent/StEvent.h +++ b/StRoot/StEvent/StEvent.h @@ -228,6 +228,7 @@ class StIstHitCollection; class StFstEvtCollection; class StFstHitCollection; class StFwdTrackCollection; +class StGmtCollection; class StEvent : public StXRefMain { public: @@ -318,6 +319,8 @@ class StEvent : public StXRefMain { const StTriggerIdCollection* triggerIdCollection() const; StTriggerData* triggerData(); const StTriggerData* triggerData() const; + StGmtCollection* gmtCollection(); + const StGmtCollection* gmtCollection() const; StSPtrVecTrackDetectorInfo& trackDetectorInfo(); const StSPtrVecTrackDetectorInfo& trackDetectorInfo() const; @@ -412,6 +415,7 @@ class StEvent : public StXRefMain { void removePsd(StPsd*); void addHitCollection(StSPtrVecHit* p, const Char_t *name); void removeHitCollection(const Char_t *name); + void setGmtCollection(StGmtCollection*); virtual Bool_t Notify(); diff --git a/StRoot/StEvent/StEventClusteringHints.cxx b/StRoot/StEvent/StEventClusteringHints.cxx index 89313aef489..0304683eb9a 100755 --- a/StRoot/StEvent/StEventClusteringHints.cxx +++ b/StRoot/StEvent/StEventClusteringHints.cxx @@ -202,6 +202,7 @@ StEventClusteringHints::StEventClusteringHints() setBranch("StRpsCollection", "evt_aux", 7); setBranch("StFttCollection", "evt_aux", 7); setBranch("StFstEvtCollection", "evt_aux", 7); + setBranch("StGmtCollection", "evt_aux", 7); setBranch("StSsdHitCollection", "evt_hits", 8); setBranch("StSstHitCollection", "evt_hits", 8); setBranch("StSvtHitCollection", "evt_hits", 8); @@ -211,6 +212,7 @@ StEventClusteringHints::StEventClusteringHints() setBranch("StTpcHitCollection", "evt_hits", 8); setBranch("StFtpcHitCollection", "evt_hits", 8); setBranch("StRnDHitCollection", "evt_hits", 8); + setBranch("StGmtHitCollection", "evt_hits", 8); setBranch("StHltEvent", "evt_hlt", 9); setBranch("StFgtCollection", "evt_fgt", 9); diff --git a/StRoot/StEvent/StEventTypes.h b/StRoot/StEvent/StEventTypes.h index 709c1866a2a..cf160cf0719 100644 --- a/StRoot/StEvent/StEventTypes.h +++ b/StRoot/StEvent/StEventTypes.h @@ -169,6 +169,8 @@ **************************************************************************/ #ifndef StEventTypes_hh #define StEventTypes_hh + +#include "StContainers.h" #include "StBbcTriggerDetector.h" #include "StCalibrationVertex.h" @@ -330,6 +332,8 @@ #include "StFgtHitCollection.h" #include "StFgtPoint.h" #include "StFgtPointCollection.h" +#include "StGmtHit.h" +#include "StGmtCollection.h" #include "StIstHit.h" #include "StIstSensorHitCollection.h" #include "StIstLadderHitCollection.h" diff --git a/StRoot/StEvent/StGmtCollection.cxx b/StRoot/StEvent/StGmtCollection.cxx new file mode 100644 index 00000000000..e65f08f53d2 --- /dev/null +++ b/StRoot/StEvent/StGmtCollection.cxx @@ -0,0 +1,58 @@ +/** + * \class StGmtColection + * \brief Holds collections of GMT data + * + * GMT data collection for StEvent (based on StFgtCollection) + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// Load header of the GMT collection +#include "StGmtCollection.h" + +ClassImp(StGmtCollection) + +//________________ +StGmtCollection::StGmtCollection() : StObject() { + // Set the module field for some of the collections + for( int i=0; igetNumStrips(); + } + return n; +} + +//________________ +size_t StGmtCollection::getNumHits() const { + size_t n = 0; + for( const StGmtHitCollection* ptr = &mHitCollection[0]; + ptr != &mHitCollection[kGmtNumModules]; ++ptr ) { + n += ptr->getNumHits(); + } + return n; +} diff --git a/StRoot/StEvent/StGmtCollection.h b/StRoot/StEvent/StGmtCollection.h new file mode 100644 index 00000000000..680524b3316 --- /dev/null +++ b/StRoot/StEvent/StGmtCollection.h @@ -0,0 +1,90 @@ +/** + * \class StGmtCollection + * \brief Holds collections of GMT data + * + * GMT data collection for StEvent (based on StFgtCollection) + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtCollection_hh +#define StGmtCollection_hh + +// StRoot headers +#include "StObject.h" +#include "StGmtStripCollection.h" +#include "StGmtHitCollection.h" +#include "StGmtPointCollection.h" +#include "StEnumerations.h" + +//________________ +class StGmtCollection : public StObject { + public: + /// Constructor + StGmtCollection(); + // StGmtCollection( const StGmtCollection& other ); ---> use default + // StGmtCollection& operator=( const StGmtCollection& other ); ---> use default + + /// Destructor + ~StGmtCollection(); + + /// Number of modules + size_t getNumModules() const { return kGmtNumModules; } + /// Number total number of strips + size_t getNumStrips() const; + /// Number of strips in the i-th module + size_t getNumStrips( unsigned short moduleIdx) const + { return (moduleIdx < kGmtNumModules ? mStripCollection[moduleIdx].getNumStrips() : 0 ); } + /// Number total number of hits + size_t getNumHits() const; + /// Number of hits in the i-th module + size_t getNumHits( unsigned short moduleIdx ) const + { return (moduleIdx < kGmtNumModules ? mHitCollection[moduleIdx].getNumHits() : 0 ); } + /// Number of points + size_t getNumPoints() const + { return mPointCollection.getNumPoints(); } + + // note: ownership of all pointers is retained by the containers. + // Do not delete any pointers received from this class. + + /// Pointer to the GMT strip collection in the i-th module + StGmtStripCollection* getStripCollection( unsigned short moduleIdx ) + { return (moduleIdx < kGmtNumModules ? &mStripCollection[moduleIdx] : 0 ); } + /// Pointer to the GMT strip collection in the i-th module + const StGmtStripCollection* getStripCollection( unsigned short moduleIdx ) const + { return (moduleIdx < kGmtNumModules ? &mStripCollection[moduleIdx] : 0 ); } + + /// Pointer to the GMT hit collection in the i-th module + StGmtHitCollection* getHitCollection( unsigned short moduleIdx ) + { return (moduleIdx < kGmtNumModules ? &mHitCollection[moduleIdx] : 0 ); } + /// Pointer to the GMT hit collection in the i-th module + const StGmtHitCollection* getHitCollection( unsigned short moduleIdx ) const + { return (moduleIdx < kGmtNumModules ? &mHitCollection[moduleIdx] : 0 ); } + + /// Pointer to the GMT point collection + StGmtPointCollection* getPointCollection() { return &mPointCollection; } + /// Pointer to the GMT point collection + const StGmtPointCollection* getPointCollection() const { return &mPointCollection; } + + /// Clear method + void Clear( Option_t *opts = "" ); + + protected: + + /// Friend class of the StMuDSTMaker + friend class StMuDstMaker; + + /// GMT strip collections for all modules + StGmtStripCollection mStripCollection[kGmtNumModules]; + /// GMT hit collections for all modules + StGmtHitCollection mHitCollection[kGmtNumModules]; + /// GMT point collection + StGmtPointCollection mPointCollection; + + private: + ClassDef(StGmtCollection,1) +}; + +#endif // #define StGmtCollection_hh diff --git a/StRoot/StEvent/StGmtHit.cxx b/StRoot/StEvent/StGmtHit.cxx new file mode 100644 index 00000000000..74c8f6bfd72 --- /dev/null +++ b/StRoot/StEvent/StGmtHit.cxx @@ -0,0 +1,52 @@ +/** + * \class StGmtHit + * \brief Holds data for the hit in GMT + * + * Data for an individual ``hit'' in GMT, i.e. a 1D cluster (based on StFgtHit). + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// C++ headers +#include + +// StEvent headers +#include "StGmtHit.h" + +ClassImp(StGmtHit) + +//________________ +StGmtHit::StGmtHit( Int_t key, Int_t module, Float_t adcX, + Float_t adcY, Float_t dadcX, Float_t dadcY, + Float_t localX, Float_t localY, + Float_t localXErr, Float_t localYErr, + Float_t sigmaX, Float_t sigmaY, + Float_t sigmaXErr, Float_t sigmaYErr) : + StHit( StThreeVectorF(localX,localY,0), StThreeVectorF(localXErr,localYErr,0), + module+1, 0.), mKey(key), mAdcX(adcX), mAdcY(adcY), mdAdcX(dadcX), + mdAdcY(dadcY), mSigmaX(sigmaX), mErrSigmaX(sigmaXErr), mSigmaY(sigmaY), + mErrSigmaY(sigmaYErr) { + /* empty */ +} + +//________________ +StGmtHit::~StGmtHit() { + /* empty */ +} + +//________________ +ostream& operator<<(ostream& os, const StGmtHit& v) { + return os << Form("Gmt m %3i ",v.getModule()) + << *((StHit *)&v) + << Form(" Adc X/Y %5.1f/%5.1f locX = %8.3f +/- %7.3f locY = %8.3f +/- %7.3f", + v.getAdcX(), v.getAdcY(), v.getLocalX(), v.getErrorLocalX(), v.getLocalY(), v.getErrorLocalY()); +} + +//________________ +void StGmtHit::Print(Option_t *option) const { + std::cout << *this << std::endl; +} + + diff --git a/StRoot/StEvent/StGmtHit.h b/StRoot/StEvent/StGmtHit.h new file mode 100644 index 00000000000..4296f6b4492 --- /dev/null +++ b/StRoot/StEvent/StGmtHit.h @@ -0,0 +1,136 @@ +/** + * \class StGmtHit + * \brief Holds data for the hit in GMT + * + * Data for an individual ``hit'' in GMT, i.e. a 1D cluster (based on StFgtHit). + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtHit_hh +#define StGmtHit_hh + +// C++ headers +#include + +// StEvent headers +#include "StHit.h" +#include "StGmtStrip.h" + +//________________ +class StGmtHit : public StHit { + public: + /// Constructor + StGmtHit( Int_t key = -1, Int_t module = -1, Float_t adcX = 0, + Float_t adcY = 0, Float_t dadcX = 0, Float_t dadcY = 0, + Float_t localX = 0, Float_t localY = 0, + Float_t localXErr = 10000, Float_t localYErr = 10000, + Float_t sigmaX = 0, Float_t sigmaY = 0, + Float_t sigmaXErr = 10000, Float_t sigmaYErr = 10000); + /// Destructor + ~StGmtHit(); + + /// Print hit information (parameters) + void Print(Option_t *option="") const; + + // + // Getters + // + + /// Unique detector ID + virtual StDetectorId detector() const { return kGmtId; } + + /// Key + Int_t getKey() const { return mKey; } + /// Module + Int_t getModule() const { return hardwarePosition() - 1; } + /// ADC in X + Float_t getAdcX() const { return mAdcX; } + /// ADC error in X + Float_t getErrorAdcX() const { return mdAdcX; } + /// ADC in Y + Float_t getAdcY() const { return mAdcY; } + /// ADC error in Y + Float_t getErrorAdcY() const { return mdAdcY; } + /// Local X coordinate + Float_t getLocalX() const { return position().x(); } + /// Local X coordinate error + Float_t getErrorLocalX() const { return positionError().x(); } + /// Local Y coordinate + Float_t getLocalY() const { return position().y(); } + /// Local Y coordinate error + Float_t getErrorLocalY() const { return positionError().y(); } + /// Position in X + Float_t getSigmaX() const { return mSigmaX; } + /// Position error in X + Float_t getErrorSigmaX() const { return mErrSigmaX; } + /// Position in Y + Float_t getSigmaY() const { return mSigmaY; } + /// Position error in X + Float_t getErrorSigmaY() const { return mErrSigmaY; } + /// Volume ID + Int_t volumeID() const { return 0; } + + // + // Setters + // + + /// Set ADC in X + void setAdcX( Float_t adc ) { mAdcX = adc; } + /// Set ADC error in X + void setErrorAdcX( Float_t error ) { mdAdcX = error; } + /// Set ADC in Y + void setAdcY( Float_t adc ) { mAdcY = adc; } + /// Set ADC error in Y + void setErrorAdcY( Float_t error ) { mdAdcY = error; } + /// Set local X + void setLocalX( Float_t position ) { mPosition.setX(position); } + /// Set local X error + void setErrorLocalX( Float_t error ) { mPositionError.setX(error); } + /// Set local Y + void setLocalY( Float_t position ) { mPosition.setY(position); } + /// Set local Y error + void setErrorLocalY( Float_t error ) { mPositionError.setY(error); } + /// Set local X + void setSigmaX( Float_t sigma ) { mSigmaX = sigma; } + /// Set local X error + void setErrorSigmaX( Float_t error ) { mErrSigmaX = error; } + /// Set local Y + void setSigmaY( Float_t sigma ) { mSigmaY = sigma; } + /// Set local Y error + void setErrorSigmaY( Float_t error ) { mErrSigmaY = error; } + /// Set module ID + void setModule( Short_t module ) { setHardwarePosition(module+1); } + + protected: + /// Unique label + Int_t mKey; + /// ADC counts in X + Float_t mAdcX; + /// ADC counts in Y + Float_t mAdcY; + /// ADC counts in X + Float_t mdAdcX; + /// ADC counts in Y + Float_t mdAdcY; + /// Position in local X + Float_t mSigmaX; + /// Position error in local X + Float_t mErrSigmaX; + /// Position in local Y + Float_t mSigmaY; + /// Position error in local Y + Float_t mErrSigmaY; + /// Uncertanity on the charge for keeping track of which + /// strips constribute to which cluster (not persistant) + Float_t mChargeUncert; + +private: + ClassDef(StGmtHit,2) +}; + +ostream& operator<<(ostream& os, StGmtHit const & v); + +#endif // #define StGmtHit_hh diff --git a/StRoot/StEvent/StGmtHitCollection.cxx b/StRoot/StEvent/StGmtHitCollection.cxx new file mode 100644 index 00000000000..a9179e3cf35 --- /dev/null +++ b/StRoot/StEvent/StGmtHitCollection.cxx @@ -0,0 +1,40 @@ +/** + * \class StGmtHitColection + * \brief Holds collections of GMT hits + * + * A collection of StGmtHit classes for StEvent. + * Basically a wrapper for an StSPtrVecGmtHit. Note, one instance of + * this class corresponds to one module. + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// StEvent headers +#include "StContainers.h" +#include "StGmtHit.h" +#include "StGmtHitCollection.h" + +ClassImp(StGmtHitCollection) + +//________________ +StGmtHitCollection::StGmtHitCollection( short moduleId ) : StObject(), mModule( moduleId ) { + /* empty*/ +} + +//________________ +StGmtHitCollection::~StGmtHitCollection() { + /* empty */ +} + +//________________ +void StGmtHitCollection::Clear( Option_t *opt ) { + // no need to delete the objects in mStripVec, is done within + // its clear function. + + // clear the vector + mHitVec.clear(); +} + + diff --git a/StRoot/StEvent/StGmtHitCollection.h b/StRoot/StEvent/StGmtHitCollection.h new file mode 100644 index 00000000000..0b883e20f3c --- /dev/null +++ b/StRoot/StEvent/StGmtHitCollection.h @@ -0,0 +1,61 @@ +/** + * \class StGmtHitCollection + * \brief Holds collection of GMT hits in the module + * + * A collection of StGmtHit classes for StEvent. + * Basically a wrapper for an StSPtrVecGmtHit. Note, one instance of + * this class corresponds to one module. + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtHitCollection_hh +#define StGmtHitCollection_hh + +// STAR headers +#include "StObject.h" +#include "StContainers.h" + +// Forward declaration +class StGmtHit; + +//________________ +class StGmtHitCollection : public StObject { + public: + /// Constructor + StGmtHitCollection( short moduleId = -1 ); + // StGmtHitCollection( const StGmtHitCollection& other ); ---> use default + // StGmtHitCollection& operator=( const StGmtHitCollection& other ); ---> use default + + /// Deconstructor + ~StGmtHitCollection(); + + /// Vector of hits that belong to the module + StSPtrVecGmtHit& getHitVec() { return mHitVec; } + /// Vector of hits that belong to the module + const StSPtrVecGmtHit& getHitVec() const { return mHitVec; } + + /// Number of hits in the module + size_t getNumHits() const { return mHitVec.size(); } + /// Module number + short getModule() const { return mModule; } + /// Set module number + void setModule( short moduleId ) { mModule = moduleId; } + + /// Clear + void Clear( Option_t *opt = "" ); + + protected: + /// Module index + Short_t mModule; + /// Vector of hits that belong to the current module + StSPtrVecGmtHit mHitVec; + + private: + ClassDef(StGmtHitCollection,1) +}; + +#endif // #define StGmtHitCollection_hh + diff --git a/StRoot/StEvent/StGmtPoint.cxx b/StRoot/StEvent/StGmtPoint.cxx new file mode 100644 index 00000000000..761f93be6f8 --- /dev/null +++ b/StRoot/StEvent/StGmtPoint.cxx @@ -0,0 +1,65 @@ +/** + * \class StGmtPoint + * \brief Holds data for the point (a.k.a. cluster) in GMT + * + * Description: data for individual ``point'' on the GMT, i.e. a pair + * of 1D clusters. Note, if errors during construction, the key will + * be set to -999. Based on StFgtHit. + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// STAR headers +#include "StGmtHit.h" +#include "StGmtPoint.h" +#include "StRoot/St_base/StMessMgr.h" + +ClassImp(StGmtPoint) + +//________________ +StGmtPoint::StGmtPoint() : StHit(), mHitLocalX(0), mHitLocalY(0) { + /* empty */ +}; + +//________________ +StGmtPoint::StGmtPoint( StGmtHit* hit1, StGmtHit* hit2, int key ) : StHit(), mKey( key ) { + if ( !hit1 || !hit2 ) { + LOG_ERROR << "Passed null pointer into StGmtPoint::StGmtPoint( StGmtHit* hit1, StGmtHit* hit2, int key )" << endm; + mKey = -999; + return; + } + + mHitLocalX = mHitLocalY = 0; + // if ( !(hit1->isY()) && hit2->isY() ) { + // mHitLocalX = hit1; + // mHitLocalY = hit2; + // } + // else if ( !(hit2->isY()) && hit1->isY() ) { + // mHitLocalX = hit2; + // mHitLocalY = hit1; + // } + + if ( !mHitLocalX || !mHitLocalY ) { + LOG_ERROR << "Constructor not provided a (LocalX,LocalY) pair." << endm; + mKey = -999; + return; + } + + if ( mHitLocalX->getModule() != mHitLocalY->getModule() ) { + LOG_ERROR << "Cluster pair are not on the same Module." << endm; + mKey = -999; + return; + } + + mHardwarePosition = mHitLocalX->hardwarePosition(); + mCharge = hit1->charge() + hit2->charge(); +} + +//________________ +StGmtPoint::~StGmtPoint() { + /* empty */ +} + + diff --git a/StRoot/StEvent/StGmtPoint.h b/StRoot/StEvent/StGmtPoint.h new file mode 100644 index 00000000000..4cd8f40604e --- /dev/null +++ b/StRoot/StEvent/StGmtPoint.h @@ -0,0 +1,67 @@ +/** + * \class StGmtPoint + * \brief Holds data for the point (a.k.a. cluster) in GMT + * + * Description: data for individual ``point'' on the GMT, i.e. a pair + * of 1D clusters. Note, if errors during construction, the key will + * be set to -999. Based on StFgtHit. + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtPoint_hh +#define StGmtPoint_hh + +// StEvent headers +#include "StHit.h" +#include "StGmtHit.h" +// Forward declaration +//class StGmtHit; + +//________________ +class StGmtPoint : public StHit { + public: + /// Default consturctor + StGmtPoint(); + /// Parametrized constructor + StGmtPoint( StGmtHit* hit1, StGmtHit* hit2, int key ); + // StGmtPoint(const StGmtPoint&); --> use default + // StGmtPoint& operator=(const StGmtPoint&); --> use default + + /// Destructor + ~StGmtPoint(); + + /// Unique detector ID + virtual StDetectorId detector() const { return kGmtId; } + /// Unique label + Int_t getKey() { return mKey; } + /// Module + Int_t getModule() + { return static_cast< Int_t >(mHardwarePosition/8); /* FIX ME */} + /// Return hit in X axis + const StGmtHit* getHitLocalX() const { return mHitLocalX; } + /// Local Y coordinate + const StGmtHit* getHitLocalY() const { return mHitLocalY; } + + /// Hit coordinate in X axis + Float_t getPositionLocalX() const { return mHitLocalX->getLocalX(); } + /// Hit coordinate in Y axis + Float_t getPositionLocalY() const { return mHitLocalY->getLocalY(); } + /// Volume ID + Int_t volumeID() const { return 0; } + + protected: + /// Unique label + Int_t mKey; + /// Hit in X axis + StGmtHit *mHitLocalX; + /// Hit in Y axis + StGmtHit *mHitLocalY; + + private: + ClassDef(StGmtPoint,1) +}; + +#endif // #define StGmtPoint_hh diff --git a/StRoot/StEvent/StGmtPointCollection.cxx b/StRoot/StEvent/StGmtPointCollection.cxx new file mode 100644 index 00000000000..999dd7efdef --- /dev/null +++ b/StRoot/StEvent/StGmtPointCollection.cxx @@ -0,0 +1,30 @@ +/** + * \class StGmtColection + * \brief Holds collections of GMT points + * + * Collection of GMT points for StEvent. Basically a wrapper + * for an StSPtrVecGmtPoint (based on StFgtPointCollection) + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// STAR headers +#include "StContainers.h" +#include "StGmtPoint.h" +#include "StGmtPointCollection.h" + +ClassImp(StGmtPointCollection) + +//________________ +StGmtPointCollection::StGmtPointCollection() : StObject() { + /* empty */ +}; + +//________________ +StGmtPointCollection::~StGmtPointCollection() { + /* empty */ +} + + diff --git a/StRoot/StEvent/StGmtPointCollection.h b/StRoot/StEvent/StGmtPointCollection.h new file mode 100644 index 00000000000..5b32ab3e5f6 --- /dev/null +++ b/StRoot/StEvent/StGmtPointCollection.h @@ -0,0 +1,52 @@ +/** + * \class StGmtCollection + * \brief Holds collections of GMT points + * + * Collection of GMT points for StEvent. Basically a wrapper + * for an StSPtrVecGmtPoint (based on StFgtPointCollection) + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtPointCollection_hh +#define StGmtPointCollection_hh + +// STAR headers +#include "StObject.h" +#include "StContainers.h" + +// Forward declaration +class StGmtPoint; + +//________________ +class StGmtPointCollection : public StObject { + public: + /// Constructor + StGmtPointCollection(); + // StGmtPointCollection( const StGmtPointCollection& other ); ---> use default + // StGmtPointCollection& operator=( const StGmtPointCollection& other ); ---> use default + + /// Destructor + ~StGmtPointCollection(); + + /// Vector with points + StSPtrVecGmtPoint& getPointVec() { return mPointVec; } + /// Vector with points + const StSPtrVecGmtPoint& getPointVec() const { return mPointVec; } + /// Number of points (size of the vector) + size_t getNumPoints() const { return mPointVec.size(); } + + /// Clear + void Clear( Option_t *opt = "" ) { mPointVec.clear(); } + + protected: + /// Vector of GMT points + StSPtrVecGmtPoint mPointVec; + + private: + ClassDef(StGmtPointCollection,1) +}; + +#endif // #define StGmtPointCollection_hh diff --git a/StRoot/StEvent/StGmtStrip.cxx b/StRoot/StEvent/StGmtStrip.cxx new file mode 100644 index 00000000000..af77edd3771 --- /dev/null +++ b/StRoot/StEvent/StGmtStrip.cxx @@ -0,0 +1,107 @@ +/** + * \class StGmtStrip + * \brief Holds data for the strip in GMT + * + * Data for an individual strip in GMT (based on StFgtStrip). + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +// StGmtStrip header +#include "StGmtStrip.h" + +ClassImp(StGmtStrip) + +// Number of defult time bins. +// Was 2 for the FGT, RW 03/15/13 +int StGmtStrip::mDefaultTimeBin = 7; + +//________________ +bool gmtStripPtrLessThan::operator() (const StGmtStrip* strip1, const StGmtStrip* strip2) const { + return ((strip1 && strip2) ? strip1->getGeoId() < strip2->getGeoId() : 0 ); +} + +//________________ +StGmtStrip::StGmtStrip() : StObject(), mGeoId(-1), + mModule(-1), mCoordNum(-1), mIsY(0), mPosition(0.0), + mMaxAdc(-9999), mMaxPedSubtractedAdc(-9999), + mMaxAdcTB(-1), mMaxPedSubtractedAdcTB(-1), + mCharge(kInvalidChargeValue), + mChargeUncert(kInvalidChargeValue), + mRdo(0), mArm(0), mApv(0), mChan(0), + mPed(0), mPedStdDev(0), mPedErr(0), mIsC(0) { + // Constuctor + for( int i = 0; i < kGmtNumTimeBins; ++i ) { + mAdc[i] = -9999; //not set + mPedSubtractedAdc[i] = -9999; //not set + } +} + +//________________ +StGmtStrip::~StGmtStrip() { + /* empty */ +} + +//________________ +StGmtStrip::StGmtStrip(const StGmtStrip& h) : StObject(), // copy the parent + mGeoId( h.mGeoId ), mModule( h.mModule ), mCoordNum( h.mCoordNum ), + mIsY( h.mIsY ), mIsC( h.mIsC ), mPosition( h.mPosition ), + mMaxAdc(h.mMaxAdc), mMaxPedSubtractedAdc(h.mMaxPedSubtractedAdc), + mMaxAdcTB(h.mMaxAdcTB), mMaxPedSubtractedAdcTB(h.mMaxPedSubtractedAdcTB), + // mClusterSeedType(h.mClusterSeedType), + mCharge( h.mCharge ), mChargeUncert(h.mChargeUncert), + mRdo( h.mRdo ), mArm( h.mArm ), mApv( h.mApv ), mChan( h.mChan ), + mPed(h.mPed), mPedStdDev(h.mPedStdDev), mPedErr(h.mPedErr) { + // Copy constructor + for( int i = 0; i < kGmtNumTimeBins; ++i ) { + mAdc[i] = h.mAdc[i]; + mPedSubtractedAdc[i] = h.mPedSubtractedAdc[i]; + } +} + +//________________ +StGmtStrip& StGmtStrip::operator=( const StGmtStrip& h) { + // Assignment operator + mGeoId = h.mGeoId; + mModule = h.mModule; + mCoordNum = h.mCoordNum; + mIsY = h.mIsY; + mIsC = h.mIsC; + mPosition = h.mPosition; + mMaxAdc = h.mMaxAdc; + mMaxPedSubtractedAdc = h.mMaxPedSubtractedAdc; + mMaxAdcTB = h.mMaxAdcTB; + mMaxPedSubtractedAdcTB = h.mMaxPedSubtractedAdcTB; +// mClusterSeedType = h.mClusterSeedType; + mCharge = h.mCharge; + mChargeUncert = h.mChargeUncert; + mRdo = h.mRdo; + mArm = h.mArm; + mApv = h.mApv; + mChan = h.mChan; + mPed=h.mPed; + mPedStdDev=h.mPedStdDev; + mPedErr=h.mPedErr; + + for( int i = 0; i < kGmtNumTimeBins; ++i ) { + mAdc[i] = h.mAdc[i]; + mPedSubtractedAdc[i] = h.mPedSubtractedAdc[i]; + } + + return *this; +} + +//________________ +ostream& operator<<(ostream& os, const StGmtStrip& v) { + return os << Form("GmtStrip gId %3i m %3i C %3i Y %1i C %1i p %8.3f",v.getGeoId(), v.getModule(), v.getCoordNum(), v.isY(), v.isC(), v.getPosition()) + << Form(" Rdo: %2i,Arm: %2i, Apv: %2i, cha: %3i",v.getRdo(), v.getArm(), v.getApv(), v.getChannel()); +} + +//________________ +void StGmtStrip::Print(Option_t *option) const { + std::cout << *this << std::endl; +} + + diff --git a/StRoot/StEvent/StGmtStrip.h b/StRoot/StEvent/StGmtStrip.h new file mode 100644 index 00000000000..c2012faa525 --- /dev/null +++ b/StRoot/StEvent/StGmtStrip.h @@ -0,0 +1,209 @@ +/** + * \class StGmtStrip + * \brief Holds data for the strip in GMT + * + * Data for an individual strip in GMT (based on StFgtStrip). + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtStrip_hh +#define StGmtStrip_hh + +// STAR headers +#include "StObject.h" +#include "St_base/StMessMgr.h" +#include "StEnumerations.h" + +//________________ +class StGmtStrip : public StObject { + public: + /// Constructor + StGmtStrip(); + /// Copy constructor + StGmtStrip(const StGmtStrip&); + /// Assignment operator + StGmtStrip& operator=( const StGmtStrip& ); + /// Destructor + ~StGmtStrip(); + /// Print strip information (parameters) + void Print(Option_t *option="") const; + + // + // Getters + // + + /// Detector ID (8 modules * 2 APV * 128 channels) + Int_t getGeoId() const { return mGeoId; }; + /// Module ID (8 modules in total) + Int_t getModule() const { return mModule; } + /// Coordinate (0-127) + Int_t getCoordNum() const { return mCoordNum; } + /// Is it a pad? + Int_t isY() const { return mIsY; } + /// Is used in a cluster + Int_t isC() const { return mIsC; } + /// Coordinate position relative to local origin (in module) + Float_t getPosition() const { return mPosition; } + /// ADC in a strip for a given time bin + Short_t getAdc( Int_t tb = -1 ) const + { return mAdc[ (tb < 0 || tb >= kGmtNumTimeBins) ? mDefaultTimeBin : tb ]; } + /// Pedestal subtracted ADC for a give time bin + Short_t getPedSubtractedAdc( Int_t tb = -1 ) const + { return mPedSubtractedAdc[ (tb < 0 || tb >= kGmtNumTimeBins) ? mDefaultTimeBin : tb ]; } + /// Maximal ADC over the time bins + Short_t getMaxAdc() const { return mMaxAdc; } + /// Maximal pedestal subtraced ADC over the time bins + Short_t getMaxPedSubtractedAdc() const { return mMaxPedSubtractedAdc; } + /// Maximal over the time bins + Short_t getMaxAdcTB() const { return mMaxAdcTB; } + /// Maximal over the time bins + Short_t getMaxPedSubtractedAdcTB() const { return mMaxPedSubtractedAdcTB; } + //Short_t getClusterSeedType() const { return mClusterSeedType; } + + /// Charge before GEM (in C) + Float_t getCharge() const { return mCharge; } + /// Charge uncertainty + Float_t getChargeUncert() const { return mChargeUncert; } + /// Coordinates from electronics + void getElecCoords( Int_t& rdo, Int_t& arm, Int_t& apv, Int_t& chan ) + { rdo = mRdo; arm = mArm; apv = mApv; chan = mChan; } + /// Pedestal + Float_t getPed() const { return mPed; } + /// Pedestal standard deviation + Float_t getPedStdDev() const { return mPedStdDev; } + /// Pedestal error + Float_t getPedErr() const { return mPedErr; } + /// Check if charge is valid + Bool_t chargeValid() const + { return mCharge != 0 && mCharge != kInvalidChargeValue; } + /// RDO number + Int_t getRdo() const { return mRdo; } + /// Arm + Int_t getArm() const { return mArm; } + /// Apv + Int_t getApv() const { return mApv; } + /// Channel number + Int_t getChannel() const { return mChan; } + /// Default time bin + static Int_t getDefaultTimeBin() { return mDefaultTimeBin; } + + // + // Setters + // + + /// Set detector GeoId + void setGeoId( Int_t geoId ) { mGeoId = geoId; } + /// Set module + void setModule( Int_t module ) { mModule = module; } + /// Set coordinate + void setCoordNum( Int_t coord ) { mCoordNum = coord; } + /// Set is it a pad + void setIsY( Int_t isY ) { mIsY = isY; } + /// Set is used in a cluster + void setIsC( Int_t isC ) { mIsC = isC; } + /// Set position relative to local origin (in module) + void setPosition( Float_t position ) { mPosition = position; } + /// Set ADC for the given time bucket + void setAdc( Short_t adc, Int_t tb = -1 ) { + mAdc[ (tb < 0 || tb >= kGmtNumTimeBins) ? mDefaultTimeBin : tb ] = adc; + if( adc > mMaxAdc ) { mMaxAdc=adc; mMaxAdcTB=tb; } + } + /// Set pedestal stubtracted ADC for the given time bucket + void setPedSubtractedAdc( Short_t adc, Int_t tb = -1 ) { + mPedSubtractedAdc[ (tb < 0 || tb >= kGmtNumTimeBins) ? mDefaultTimeBin : tb ] = adc; + if( adc > mMaxPedSubtractedAdc ) { mMaxPedSubtractedAdc=adc; mMaxPedSubtractedAdcTB=tb; } + } + /// Set maximal ADC over time buckets + void setMaxAdc(Short_t adc) { mMaxAdc=adc; } + /// Set maximal pedestal subtracked ADC over time buckets + void setMaxPedSubtractedAdc(Short_t adc) { mMaxPedSubtractedAdc = adc; } + // void setClusterSeedType( Short_t seedType ) { mClusterSeedType=seedType; } + + /// Set charge before the GEM (in C) + void setCharge( Float_t charge ) { mCharge = charge; } + /// Set charge uncertainty + void setChargeUncert( Float_t chargeUncert) { mChargeUncert = chargeUncert; } + /// Set coordinates from electronics + void setElecCoords( Int_t rdo, Int_t arm, Int_t apv, Int_t chan ) + { mRdo = rdo; mArm = arm; mApv = apv; mChan = chan; } + /// Set pedestal + void setPed(Float_t ped) { mPed=ped; } + /// Set pedestal standard deviation + void setPedStdDev(Float_t pedStdDev) { mPedStdDev=pedStdDev; } + /// Set pedestal error + void setPedErr(Float_t pedErr) { mPedErr=pedErr; } + /// Set charge to the invalid state + void invalidateCharge() { mCharge = kInvalidChargeValue; }; + /// Set default time bin + static void setDefaultTimeBin( Int_t tb ) { mDefaultTimeBin = tb; } + + protected: + /// Indexing: 8 modules * 2 APV * 128 channels = 2048 + Int_t mGeoId; + /// Indexing: 8 modules + Int_t mModule; + /// 0-127 in each dimension (X and Y) + Int_t mCoordNum; + /// Is it a pad (as opposed to a strip)? + Int_t mIsY; + /// Coordinate position relative to local origin (in module) + Float_t mPosition; + /// ADC in a strip. Note "StRoot/RTS/src/DAQ_GMT/daq_gmt.h" uses UShort_t + Short_t mAdc[kGmtNumTimeBins]; + /// ADC after pedestal subtraction + /// Note "StRoot/RTS/src/DAQ_GMT/daq_gmt.h" uses UShort_t + Short_t mPedSubtractedAdc[kGmtNumTimeBins]; + /// Maximal ADC over the time bins + Short_t mMaxAdc; + /// Maximal pedestal subtracted ADC over the time bins + Short_t mMaxPedSubtractedAdc; + /// Maximal over the time bins + Short_t mMaxAdcTB; + /// Max over the time bins + Short_t mMaxPedSubtractedAdcTB; + // Short_t mClusterSeedType; // See types in StGmtConsts.h + /// Charge before GEM, units (C) + /// relation: ADC = ped + charge*gain(r,phi,disc) + Float_t mCharge; + /// Charge uncertainty + Float_t mChargeUncert; + + // elec coords, straight out of the DAQ file + + /// RDO number + Int_t mRdo; + Int_t mArm; + Int_t mApv; + /// Channel number + Int_t mChan; + + /// Pedestal + Float_t mPed; + /// Pedestal standard deviation + Float_t mPedStdDev; + /// Pedestal RMS + Float_t mPedErr; + // Is used in a cluster ? + Int_t mIsC; + /// Time bin + static Int_t mDefaultTimeBin; + + /// To signify an invalid value of the charge + enum { kInvalidChargeValue = -10000 }; + + private: + ClassDef(StGmtStrip,1) +}; + +ostream& operator<<(ostream& os, StGmtStrip const & v); + +// Functor for sorting the strips in the strip weight map. +struct gmtStripPtrLessThan { + bool operator() (const StGmtStrip* strip1, const StGmtStrip* strip2) const; +}; + + +#endif // #define StGmtStrip_hh diff --git a/StRoot/StEvent/StGmtStripCollection.cxx b/StRoot/StEvent/StGmtStripCollection.cxx new file mode 100644 index 00000000000..79909c74199 --- /dev/null +++ b/StRoot/StEvent/StGmtStripCollection.cxx @@ -0,0 +1,160 @@ +/*************************************************************************** + * + * Authors: K.S. Engle and Richard Witt (witt@usna.edu), Jan 2013 + * based on StFgtStripCollection + * + *************************************************************************** + * + * Description: See header file. + * + ***************************************************************************/ + +#include "St_base/StMessMgr.h" + +#include "StContainers.h" +#include "StGmtStrip.h" +#include "StGmtStripCollection.h" +#include +#include + +using namespace std; + +// deconstructor +StGmtStripCollection::~StGmtStripCollection() {/* no op */} + +// remove all hits with negative geoIds or with clusterSeedType set to +// kGmtDeadStrip +void StGmtStripCollection::removeFlagged(){ + if( !mStripVec.empty() ){ + // container to hold a copy + vector< StGmtStrip* > copy; + copy.reserve( mStripVec.size() ); + sortByGeoId(); + + // iterators + StSPtrVecGmtStripIterator srcIter; + StSPtrVecGmtStripIterator lastCopied=mStripVec.begin()-1; + + // copy all valid events + for( srcIter = mStripVec.begin(); srcIter != mStripVec.end(); ++srcIter ) +// if( (*srcIter) && (*srcIter)->getClusterSeedType() != kGmtDeadStrip && (*srcIter)->getGeoId() >= 0 ) + if( (*srcIter) && (*srcIter)->getGeoId() >= 0 ) + copy.push_back( new StGmtStrip( *(*srcIter) ) ); + + if ( copy.size() != mStripVec.size() ){ + // this deletes the objects + mStripVec.clear(); + // note: ownership of new objects passed to StSPtrVec + vector< StGmtStrip* >::iterator copyIter; + for( copyIter = copy.begin(); copyIter != copy.end(); ++copyIter ) + mStripVec.push_back( *copyIter ); + } + } +} + +bool StGmtStripCollection::hitGeoIdLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ){ + return h1->getGeoId() < h2->getGeoId(); +}; + +bool StGmtStripCollection::hitCoordLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ){ + return h1->getCoordNum() < h2->getCoordNum(); +// if( (h1->isY() && h2->isY()) || ( !(h1->isY()) && !(h2->isY()) ) ) +// { +// // LOG_ERROR << "StGmtStripCollection::hitCoordLessThan sort CASE 1: " << +// // "\n\t h1 " << h1 << +// // "\n\t h1->isY() " << h1->isY() << +// // "\n\t h1->GetCoordNum() " << h1->getCoordNum() << +// // "\n\t h2 " << h2 << +// // "\n\t h2->isY() " << h2->isY() << +// // "\n\t h2->GetCoordNum() " << h2->getCoordNum() << +// // "\n\t will return: " << h1->getCoordNum() << " < " << h2->getCoordNum() << endm; +// return h1->getCoordNum() < h2->getCoordNum(); +// } +// else if( h1->isY() && !(h2->isY()) ) +// { +// // LOG_ERROR << "StGmtStripCollection::hitCoordLessThan sort CASE 2: " << +// // "\n\t h1 " << h1 << +// // "\n\t h1->isY() " << h1->isY() << +// // "\n\t h1->GetCoordNum() " << h1->getCoordNum() << +// // "\n\t h2 " << h2 << +// // "\n\t h2->isY() " << h2->isY() << +// // "\n\t h2->GetCoordNum() " << h2->getCoordNum() << +// // "\n\t will return: " << h1->getCoordNum()+kGmtNumStrips << " < " << h2->getCoordNum() << endm; +// return (h1->getCoordNum()+kGmtNumStrips) < h2->getCoordNum(); // order X first +// } +// else if( !(h1->isY()) && h2->isY() ) +// { +// // LOG_ERROR << "StGmtStripCollection::hitCoordLessThan sort CASE 3: " << +// // "\n\t h1 " << h1 << +// // "\n\t h1->isY() " << h1->isY() << +// // "\n\t h1->GetCoordNum() " << h1->getCoordNum() << +// // "\n\t h2 " << h2 << +// // "\n\t h2->isY() " << h2->isY() << +// // "\n\t h2->GetCoordNum() " << h2->getCoordNum() << +// // "\n\t will return: " << h1->getCoordNum() << " < " << h2->getCoordNum()+kGmtNumStrips << endm; +// return h1->getCoordNum() < (h2->getCoordNum()+kGmtNumStrips); // order X first +// } +// else +// { +// // LOG_ERROR << "StGmtStripCollection::hitCoordLessThan sort FAILED: " << +// // "\n\t h1 " << h1 << +// // "\n\t h1->isY() " << h1->isY() << +// // "\n\t h1->GetCoordNum() " << h1->getCoordNum() << +// // "\n\t h2 " << h2 << +// // "\n\t h2->isY() " << h2->isY() << +// // "\n\t h2->GetCoordNum() " << h2->getCoordNum() << +// // "\n\t will return: " << -1 << endm; +// return -1; +// } +}; + +bool StGmtStripCollection::hitLayerLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ){ + return h1->isY() < h2->isY(); +}; + +void StGmtStripCollection::Clear( Option_t *opt ){ + + // no need to delete the objects in mStripVec, is done within its + // clear function. + + // clear the vector + mStripVec.clear(); + + // clear the vector for alternate lookups + for (unsigned int i=0; i(0); + + // clear the other vector for alternate lookups + for (unsigned int i=0; i(0); + +} + +StGmtStrip* StGmtStripCollection::getStrip( Int_t Id ){ // using geoId now instead of elecId so now using more generic index name +// StGmtStrip* &stripPtr = mStripElecIdVec[Id]; + StGmtStrip* &stripPtr = mStripGeoIdVec[Id]; + if( !stripPtr ){ + stripPtr = new StGmtStrip(); + mStripVec.push_back( stripPtr ); + } + return stripPtr; +} + +// StGmtStrip* StGmtStripCollection::getStrip( Int_t Id ){ // using geoId now instead of elecId so now using more generic index name +// StGmtStrip* &stripPtr = mStripVec[Id]; +// if( !stripPtr ){ +// stripPtr = new StGmtStrip(); +// mStripVec.push_back( stripPtr ); +// } +// return stripPtr; +// } + +StGmtStrip* StGmtStripCollection::getSortedStrip( Int_t Id ){ // using geoId now instead of elecId so now using more generic index name +// StGmtStrip* &stripPtr = mStripElecIdVec[Id]; + StGmtStrip* &stripPtr = mStripVec[Id]; + if( !stripPtr ){ + LOG_ERROR << "StGmtStripCollection::getSortedStrip no such Id: " << Id << endm; + return 0; + } + return stripPtr; +} + +ClassImp(StGmtStripCollection) diff --git a/StRoot/StEvent/StGmtStripCollection.h b/StRoot/StEvent/StGmtStripCollection.h new file mode 100644 index 00000000000..da5decce669 --- /dev/null +++ b/StRoot/StEvent/StGmtStripCollection.h @@ -0,0 +1,155 @@ +/** + * \class StGmtCollection + * \brief Holds collections of GMT strips + * + * Collection of GMT strips for StEvent. Basically a wrapper + * for an StSPtrVecGmtStrip (based on StFgtStripCollection) + * + * \author K.S. Engle, Jan. 2013 + * \author Richard Witt (witt@usna.edu), Jan. 2013 + * \author Grigory Nigmatkulov (nigmatkulov@gmail.com), Dec. 2020 + */ + +#ifndef StGmtStripCollection_hh +#define StGmtStripCollection_hh + +// STAR headers +#include "StObject.h" +#include "StContainers.h" +#include "StGmtStrip.h" + +//________________ +class StGmtStripCollection : public StObject { + public: + /// Constructer + StGmtStripCollection( short module = 0 ); + // StGmtStripCollection( const StGmtStripCollection& other ); ---> use default + // StGmtStripCollection& operator=( const StGmtStripCollection& other ); ---> use default + + /// Destructor + ~StGmtStripCollection(); + + // WARNING: never use getStripVec().push_back() or equivelants. + // Instead use StGmtStripCollection::getStrip to add a new strip. + + /// Access vector with strips + StSPtrVecGmtStrip& getStripVec() { return mStripVec; } + /// Access vector with strips + const StSPtrVecGmtStrip& getStripVec() const { return mStripVec; } + + // sort internal vector by geoId + void sortByGeoId(); + + // sort internal vector by coordinate number + void sortByCoord(); + // sort internal vector by coordinate number + void partialSortByCoord(); + + // sort internal vector by layer (X first then Y) + void sortByLayer(); + + // remove all hits with negative geoIds + void removeFlagged(); + + // size of internal vector + size_t getNumStrips() const; + + // modify/access the moduleId + short getModule() const; + void setModule( short module ); + + // Clear + void Clear( Option_t *opt = "" ); + + // Get pointer to a strip -- note: this is the only way to modify a + // strip. New strip is created if it does not exist, but only + // using StGmtStrip() constructor. Ownership is retained by the + // collection. +// StGmtStrip* getStrip( int elecId ); + StGmtStrip* getStrip( int Id ); + StGmtStrip* getSortedStrip( int Id ); + + protected: + /// Function used for sorting strips by geoId + static bool hitGeoIdLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ); + /// Function used for sorting strips by coordinate number + static bool hitCoordLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ); + /// Function used for sorting strips by X then Y + static bool hitLayerLessThan( const StGmtStrip* h1, const StGmtStrip* h2 ); + + /// Module ID + Short_t mModule; + /// Vector with strips + StSPtrVecGmtStrip mStripVec; + + /// Temporary copy of the pointers, indexed by ElecId + /// used for the addStripInfo class + StPtrVecGmtStrip mStripElecIdVec; + /// Temporary copy of the pointers, indexed by GeoId + /// used for the addStripInfo class + StPtrVecGmtStrip mStripGeoIdVec; + + private: + ClassDef(StGmtStripCollection,1) +}; + + +// inline functions + +// inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { +// mStripElecIdVec.resize( kGmtNumElecIds ); +// for (unsigned int i=0; i(0); +// }; + +inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { + mStripGeoIdVec.resize( kGmtNumGeoIds ); + for (unsigned int i=0; i(0); +}; + +// inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { +// mStripVec.resize( kGmtNumGeoIds ); +// for (unsigned int i=0; i(0); +// }; + + +// sort by geoId +inline void StGmtStripCollection::sortByGeoId(){ + sort( mStripVec.begin(), mStripVec.end(), &StGmtStripCollection::hitGeoIdLessThan ); + return; +}; + +// sort by layer (X first then Y) +inline void StGmtStripCollection::sortByLayer(){ + sort( mStripVec.begin(), mStripVec.end(), &StGmtStripCollection::hitLayerLessThan ); + return; +}; + +// sort by coordinate number +inline void StGmtStripCollection::partialSortByCoord(){ + partial_sort( mStripVec.begin(), mStripVec.begin()+kGmtNumStrips, mStripVec.begin()+kGmtNumStrips, &StGmtStripCollection::hitCoordLessThan ); + return; +}; + +// sort by coordinate number +inline void StGmtStripCollection::sortByCoord(){ + sort( mStripVec.begin(), mStripVec.end(), &StGmtStripCollection::hitCoordLessThan ); + return; +}; + + +inline size_t StGmtStripCollection::getNumStrips() const { + return mStripVec.size(); +}; + +inline void StGmtStripCollection::setModule( short moduleId ) { + mModule = moduleId; +}; + +inline short StGmtStripCollection::getModule() const { + return mModule; +}; + +#endif // #define StGmtStripCollection_hh diff --git a/StRoot/StGmtAligner/EventT.cxx b/StRoot/StGmtAligner/EventT.cxx new file mode 100644 index 00000000000..d6f92daec06 --- /dev/null +++ b/StRoot/StGmtAligner/EventT.cxx @@ -0,0 +1,551 @@ +#include +#include "TH1.h" +#include "TH2.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TF1.h" +#include "StEvent.h" +#include "StPrimaryVertex.h" +#include "StEventInfo.h" +#include "StEventSummary.h" +#include "StTrack.h" +#include "StTrackNode.h" +#include "StPrimaryTrack.h" +#include "StGlobalTrack.h" +#include "StTrackDetectorInfo.h" +#include "StTrackGeometry.h" +#include "StDcaGeometry.h" +#include "TGeoMatrix.h" +#include "StarRoot/THelixTrack.h" +#include "EventT.h" +#include "TrackT.h" +#include "HitT.h" +#include "TKey.h" +#include "TDirectory.h" +#include "TClass.h" +#include "TRVector.h" +#include "TRSymMatrix.h" +#include "StGmtHit.h" +#include "StGmtHitCollection.h" +#include "StGmtCollection.h" +#define __OnlyGlobal__ +ClassImp(EventTHeader); +ClassImp(EventT); +ClassImp(TrackT); +ClassImp(HitT); + +TClonesArray *EventT::fgTracks = 0; +TClonesArray *EventT::fgHits = 0; +THashList *EventT::fRotList = 0; + + +static Int_t _debug = 0; +//______________________________________________________________________________ +EventT::EventT() : fIsValid(kFALSE) +{ + // Create an EventT object. + // When the constructor is invoked for the first time, the class static + // variable fgTracks is 0 and the TClonesArray fgTracks is created. + + if (!fgTracks) fgTracks = new TClonesArray("TrackT", 1000); + fTracks = fgTracks; + fNtrack = 0; + if (!fgHits) fgHits = new TClonesArray("HitT", 1000); + fHits = fgHits; + fNhit = 0; +} + +//______________________________________________________________________________ +EventT::~EventT() +{ + Clear(); + SafeDelete(fRotList); +} + +//______________________________________________________________________________ +Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { + static const Int_t NoFitPointCutForGoodTrackT = 15; + Clear(); + Int_t iok = 1; + fIsValid = kFALSE; + if (! pEventT) return iok; + StGmtCollection* GmtCollection = pEventT->gmtCollection(); + if (! GmtCollection) { cout << "No GMT Collections" << endl; return iok;} + StSPtrVecTrackNode& theNodes = pEventT->trackNodes(); + UInt_t nnodes = theNodes.size(); + if (! nnodes) { cout << "No tracks" << endl; return iok;} + StEventInfo* info = pEventT->info(); + Int_t ev = 0, run = 0, time = 0; + if (info) { + ev = info->id(); + run = info->runId(); + time = info->time(); + } + StEventSummary* summary = pEventT->summary(); + Double32_t field = 0; + if (summary) field = summary->magneticField(); + SetHeader(ev,run,time,field); + SetFlag(1); + Int_t nTotMatch = 0; + // Create and Fill the TrackT objects + for (UInt_t i=0; itrack(global); + if (! Track) continue; + StGlobalTrack *gTrack = (StGlobalTrack *) Track; + StDcaGeometry *dca = gTrack->dcaGeometry(); + if (! dca) continue; + StPrimaryTrack *pTrack = (StPrimaryTrack *) node->track(primary); +#ifndef __OnlyGlobal__ + if (! pTrack) continue; + Track = (StTrack *) pTrack; +#endif + StThreeVectorD g3 = Track->geometry()->momentum(); + Double_t pT = g3.perp(); + Double_t pMom = g3.mag(); + if (pMom < pCut) continue; + TrackT *track = AddTrackT(); + Double_t InvpT = 0; + Double_t TanL = 999999; + if (TMath::Abs(pT) > 1.e-7) { + InvpT = Track->geometry()->charge()/pT; + TanL = g3.z()/pT; + } + track->SetInvpT(InvpT); + track->SetPhi(TMath::ATan2(g3.y(),g3.x())); + track->SetTanL(TanL); + static const Double_t EC = 2.9979251E-4; + Double_t Rho = - EC*InvpT*field; + track->SetRho(Rho); + track->SetLength(Track->length()); + StTrackDetectorInfo* dinfo=Track->detectorInfo(); + track->SetNpoint(dinfo->numberOfPoints()); + track->SetNPpoint(Track->numberOfPossiblePoints()); + track->SetN(0); + StPhysicalHelixD helixO = Track->outerGeometry()->helix(); + Double_t R_tof[2]= {210., 216.}; // inner and outer surfaces + pair shR = helixO.pathLength(0.5*(R_tof[0]+R_tof[1])); + if (TMath::Abs(shR.first) > 200 && TMath::Abs(shR.second) > 200) continue; + Double_t stepR = (shR.first > 0) ? shR.first : shR.second; + StThreeVectorD xyzR = helixO.at(stepR); + Double_t phiR = TMath::RadToDeg()*xyzR.phi(); + if (_debug) + cout << "\t shR " << shR.first << "\t" << shR.second << "\tstepR " << stepR + << "\txyzR\t" << xyzR << "\tphiR\t" << phiR << endl; + Int_t NoHitPerTrack = 0; + THashList *fRotList = RotMatrices(); + if (! fRotList) continue; + TIter next(fRotList); + TGeoHMatrix *comb = 0; + + while ((comb = (TGeoHMatrix *) next())) { + TString combName(comb->GetName()); + if (! combName.BeginsWith("R")) continue; + Int_t Id; + sscanf(comb->GetName()+1,"%i",&Id); + UInt_t module = Id; + StGmtHitCollection* GmtHitCollection = GmtCollection->getHitCollection(module); + if (! GmtHitCollection) { cout << "No GMT HitT Collections for mudule " << module << endl; continue;} + StSPtrVecGmtHit& hitvec = GmtHitCollection->getHitVec(); + UInt_t NoHits = hitvec.size(); + if (! NoHits) continue; + if (_debug) { + cout << comb->GetName() << "\tmodule = " << module << endl; + comb->Print(); + } + static Double_t dz[2] = {50.00, 2.10}; + static Double_t dx[2] = {50.00, 3.65}; + Int_t k = 0; // gmt + Double_t *rot = comb->GetRotationMatrix(); + Double_t *tra = comb->GetTranslation(); + const StThreeVectorD unit(1.,0.,0.); + const StThreeVectorD zero(0.,0.,0.); + StThreeVectorD normal(rot[2], rot[5], rot[8]); + StThreeVectorD middle(tra); + if (_debug) cout << "middle:" << middle << "\tnormal:" << normal << endl; + comb->LocalToMaster(zero.xyz(),middle.xyz()); + comb->LocalToMasterVect(unit.xyz(), normal.xyz()); + if (_debug) cout << "middle:" << middle << "\tnormal:" << normal << endl; + Double_t phiM = TMath::RadToDeg()*middle.phi(); + if (_debug) cout << "phiR = " << phiR << "\tphiM = " << phiM << endl; + Double_t dPhi = phiR - phiM; + if (dPhi > 360) dPhi -= 360; + if (dPhi < -360) dPhi += 360; + if (TMath::Abs(dPhi) > 15) continue; + if (_debug) cout << "zR = " << xyzR.z() << "\tzM = " << tra[2] << endl; + if (TMath::Abs(xyzR.z() - tra[2]) > 20) continue; + Double_t sh = helixO.pathLength(middle, normal); + if (_debug) { + cout << "StHelix sh " << sh + << "\t shR " << shR.first << "\t" << shR.second + << endl; + } + StThreeVectorD xyzG = helixO.at(sh); if (_debug) cout << "StHelix xyzG\t" << xyzG << endl; + StThreeVectorD dR = xyzR - xyzG; if (_debug) cout << "dR\t" << dR << " dist = " << dR.magnitude() << endl; + if (dR.magnitude() > 50) continue; + if (sh < -5e2 || sh > 5e2) continue; + if (_debug) { + StThreeVectorD dX = xyzG - helixO.at(0); + cout << "Qi: " << Track->geometry()->charge() + << "\tQo: " << Track->outerGeometry()->charge() + << "\tdX " << dX << endl; + cout << *dca << endl; + } + Double_t uvPred[3]; + comb->MasterToLocal(xyzG.xyz(),uvPred); + TRVector xyzL(3,uvPred); if (_debug) cout << "StHelix xyzL\t" << xyzL << endl; + Double_t dirGPred[3] = {helixO.cx(sh),helixO.cy(sh),helixO.cz(sh)}; + Double_t dxyzL[3]; + comb->MasterToLocalVect(dirGPred,dxyzL); + Double_t tuvPred[2] = {dxyzL[1]/dxyzL[0], dxyzL[2]/dxyzL[0]}; + if (_debug) cout << "StHelix tU/tV = " << tuvPred[0] << "\t" << tuvPred[1] << endl; + if (TMath::Abs(uvPred[1]) > dx[k] + 1.0) continue; + if (TMath::Abs(uvPred[2]) > dz[k] + 1.0) continue; + + Bool_t mIsCrossingMembrain = kTRUE; + Bool_t mIsPrimary = kFALSE; + if(pTrack!=0) { + mIsPrimary = kTRUE; + } + StThreeVectorF firstPoint = Track->detectorInfo()->firstPoint(); + StThreeVectorF lastPoint = Track->detectorInfo()->lastPoint(); + if( (firstPoint.z()>0 && lastPoint.z()>0 && xyzG.z()>0) || + (firstPoint.z()<0 && lastPoint.z()<0 && xyzG.z()<0) ) { + mIsCrossingMembrain = kFALSE; + } + + for (UInt_t l = 0; l < NoHits; l++) { + StHit *hit = hitvec[l]; + if (hit) { + //if (hit->flag()>=4) continue; + //if (hit->flag()< 0) continu; + // cout << "hitFlag=" << hit->flag() << endl; + + HitT *h = AddHitT(); + h->SetHitLength(sh); + h->SetHitLength(stepR); + h->SetHitdR(dR.magnitude()); + h->SetHitFlag(UInt_t(hit->flag())); + h->SetUVPred (uvPred[1],uvPred[2]); + h->SettUVPred(tuvPred[0],tuvPred[1]); + h->SetXyzG(xyzG.xyz()); + h->SetDirG(dirGPred); + h->SetisPrimary(mIsPrimary); + h->SetisCrossingMembrain(mIsCrossingMembrain); + SetHitT(h, hit, comb, track); + NoHitPerTrack++; + h->SetHitPerTrack(NoHitPerTrack); + // SetHitT(h, hit, comb, track, &TPDeriv); + } + } + } + nTotMatch += NoHitPerTrack; + } + fIsValid = kTRUE; + if (nTotMatch) iok = 0; + return iok; +} +//______________________________________________________________________________ +TrackT *EventT::AddTrackT() +{ + // Add a new track to the list of tracks for this event. + // To avoid calling the very time consuming operator new for each track, + // the standard but not well know C++ operator "new with placement" + // is called. If tracks[i] is 0, a new TrackT object will be created + // otherwise the previous TrackT[i] will be overwritten. + + TClonesArray &tracks = *fTracks; + TrackT *track = new(tracks[fNtrack++]) TrackT(); + //Save reference to last TrackT in the collection of Tracks + return track; +} +//______________________________________________________________________________ +HitT *EventT::AddHitT() +{ + // Add a new hit to the list of hits for this event. + // To avoid calling the very time consuming operator new for each hit, + // the standard but not well know C++ operator "new with placement" + // is called. If hits[i] is 0, a new HitT object will be created + // otherwise the previous HitT[i] will be overwritten. + + TClonesArray &hits = *fHits; + HitT *hit = new(hits[fNhit++]) HitT(); + //Save reference to last HitT in the collection of Hits + return hit; +} + +//______________________________________________________________________________ +void EventT::Clear(Option_t * /*option*/) +{ + fTracks->Clear("C"); //will also call TrackT::Clear + fHits->Clear("C"); //will also call HitT::Clear +} + +//______________________________________________________________________________ +void EventT::Reset(Option_t * /*option*/) +{ + // Static function to reset all static objects for this event + // fgTracks->Delete(option); + + delete fgTracks; fgTracks = 0; + delete fgHits; fgHits = 0; +} + +//______________________________________________________________________________ +void EventT::SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field) +{ + fNtrack = 0; + fNhit = 0; + fEvtHdr.Set(i, run, date, field); +} +//________________________________________________________________________________ +void EventT::Print(Option_t *opt) const { + cout << "Run/EventT\t" << fEvtHdr.GetRun() << "/" << fEvtHdr.GetEvtNum() << "\tDate " << fEvtHdr.GetDate() + << "\tField " << fEvtHdr.GetField() << endl; + cout << "Total no. tracks " << GetTotalNoTracks() << "\tRecorded tracks " << GetNtrack() + << "\tRecorded hits " << GetNhit() << endl; + TRVector vertex(3,GetVertex()); + TRSymMatrix cov(3,GetCovMatrix()); + cout << "Primary vertex " << vertex << endl; + cout << "Its cov. matrix " << cov << endl; + for (UInt_t i = 0; i < GetNtrack(); i++) {cout << i << "\t"; GetTrackT(i)->Print();} + for (UInt_t i = 0; i < GetNhit(); i++) {cout << i << "\t"; GetHitT(i)->Print();} + +} +//________________________________________________________________________________ +HitT *EventT::SetHitT(HitT *h, StHit *hit, TGeoHMatrix *comb, TrackT *track) { + UInt_t B = 0, L = 0, l = 0, W = 0, H = 0; + Int_t rdo = 0; + h->SetRDO(rdo); + if (hit->detector() == kGmtId) { + StGmtHit *ht = (StGmtHit *) hit; + B = ht->getModule(); + h->SetId(B,L,l,W,H); + h->SetuvD(ht->getLocalY(), ht->getLocalX()); + h->SetuvDError(ht->getErrorLocalY(), ht->getErrorLocalX()); + h->SetSigma(ht->getSigmaY(), ht->getSigmaX()); + h->SetSigmaError(ht->getErrorSigmaY(), ht->getErrorSigmaX()); + h->SetAdc(ht->getAdcY(), ht->getAdcX()); + h->SetAdcError(ht->getErrorAdcY(), ht->getErrorAdcX()); + h->SetUsedInFit(hit->usedInFit()); + } + StThreeVectorF position = hit->position(); + Double_t xyzG[3] = {position.x(),position.y(),position.z()}; + h->SetGC(xyzG[0],xyzG[1],xyzG[2]); + Double_t xyzL[3] = {0,0,0}; + comb->MasterToLocal(xyzG,xyzL); + // if (TMath::Abs(xyzL[2]) > 0.1) continue; + Double_t uvw[3] = {h->GetU(),h->GetV(),0}; + comb->LocalToMaster(uvw,xyzG); + h->Set(xyzG,uvw); + Double_t *rot = comb->GetRotationMatrix(); + h->SetWG(rot[2],rot[5],rot[8]); + // Int_t IdH = GetIndexOfHitT(h); + Int_t IdH = fNhit - 1; + track->SetHitTId(IdH); + Double_t invpT = track->GetInvpT(); + if (TMath::Abs(invpT) < 1e-7) invpT = 1e-7; + h->SetpT(1./invpT); + h->SetMom(track->GetMomentum()); + h->SetWG(rot[2],rot[5],rot[8]); + TGeoHMatrix *rotL = (TGeoHMatrix *) RotMatrices()->FindObject(Form("WL%s",comb->GetName()+1)); + Double_t xyzLadder[3] = {0,0,0}; + if (rotL) { + + rotL->LocalToMaster(uvw,xyzLadder); + h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); + Double_t uvwP[3] = {h->GetPredU(),h->GetPredV(),0}; + rotL->LocalToMaster(uvwP,xyzLadder); + h->SetXyzL(xyzLadder); + } else { + + cout << Form("WL%s",comb->GetName()+1) << " has not been found" << endl; + h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); + h->SetXyzL(xyzLadder); + } + return h; +} +//________________________________________________________________________________ +void TrackT::Print(Option_t *opt) const { + cout << "TrackT: InvpT " << fInvpT << "\tTanL " << fTanL + << "\tPhi " << fPhi << "\tRho " << fRho + << "\tNpoint " << fNpoint << "\tNsp " << fNsp << endl; + for (UInt_t i = 0; i < fNsp; i++) cout << "\t" << fIdHitT[i]; + cout << endl; +} +//________________________________________________________________________________ +void HitT::SetId(Int_t B, Int_t L, Int_t l, Int_t W, Int_t H) { + barrel = B; layer = L; ladder = l; wafer = W; hybrid = H; +} +//________________________________________________________________________________ +void HitT::Print(Option_t *opt) const { + cout << "HitT: Id " << Id << "\tpT = " << pT << "\tmomentum " << pMom << endl; + TRVector glob(3,&xG); cout << "Global :" << glob << endl; + cout << "Local u/v/w " << u << "/ " << v << "/ " << w << endl; + cout << "Prediction uP/vP " << uP << "/ " << vP << "\ttuP/tvP " << tuP << "/ " << tvP << endl; +} +//________________________________________________________________________________ +void EventT::RestoreListOfRotations() { + if (fRotList) return; + if (! gDirectory) return; + fRotList = new THashList(100,0); + fRotList->SetOwner(); + TIter nextkey(gDirectory->GetListOfKeys() ); + TKey *key; + while ((key = (TKey*) nextkey())) { + TObject *obj = key->ReadObj(); + if ( obj->IsA()->InheritsFrom( "TGeoHMatrix" ) ) { + fRotList->Add(obj); + } + } +} +//________________________________________________________________________________ +void TBase::Loop(Int_t Nevents) { +#if 1 + struct PlotPar_t { + Char_t *Name; + Char_t *Title; + Int_t nx; + Int_t ny; + Double_t xmin; + Double_t xmax; + Double_t ymin; + Double_t ymax; + }; + const PlotPar_t plotUP = // plots for uP + { "uP","track u", 320, 3, -5.,5., 0.,3. }; + + const PlotPar_t plotDu = // plots for u-uP + { "Du","Du before cut", 250, 3, -2.,2., 0.,3. }; + + const PlotPar_t plotDuv = // plots for du & dv + { "Du","Du cuts", 200, 3, -1.,1., 0.,3. }; + + TFile *fOut = new TFile(fOutFileName,"recreate"); + TString Name; + TString Title; + TString uName; + TString uTitle; + enum {NM = 8}; // no. of modules + // B + TH1F *LocPlots[NM]; + TH1F * uPlots[NM]; + TH1F * hpT = new TH1F( "Pt", "pt", 200, -2., 2.); + TH1F * hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); + TH1F * uPAll = new TH1F("UPall","uPall", plotUP.nx, plotUP.xmin, plotUP.xmax); + TH1F * uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); + TH1F * duB[NM][2]; + TH1F * dvB[NM]; + TH1F * uCuts[NM]; + TH1F * xCuts[NM]; + TH1F * uCut = new TH1F("Ucut","uc", plotDu.nx, plotDu.xmin, plotDu.xmax); + TH1F * vCut = new TH1F("Vcut","vc", 200, -3., 3.); + TH2F * dMin = new TH2F("DMin","vumin",100,-0.75,0.75,100,-0.75,0.75); + TH1F * vMin = new TH1F("VMin","vmin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + TH1F * uMin = new TH1F("UMin","umin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + TH1F * uMinC = new TH1F("UMinC","umC", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + memset(LocPlots,0,NM*sizeof(TH1F *)); + memset( uPlots,0,NM*sizeof(TH1F *)); + for (int M = 0; M < NM; M++) {// over gmt Modules + uName = Form("UModule%i", M+1); + uTitle = Form("du for M%i", M+1); + duB[M][0] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + uName = Form("VModule%i", M+1); + uTitle = Form("dv for M%i", M+1); + dvB[M] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + uName = Form("UModule%iVcut", M+1); + uTitle = Form("du for M%i after Vcut", M+1); + duB[M][1] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + Int_t module = M; + uName = plotUP.Name; + uName += Form("M%i", module); + uTitle = Form("uP for Module %i", module); + uPlots[module] = + new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); + uName = Form("%sM%i", plotDu.Name, module); + uTitle = Form("u-uP for M %i", module); + uCuts[module]= + new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + uName = Form("%sxM%i", plotDu.Name, module); + uTitle = Form("u-uP corr M %i", module); + xCuts[module]= + new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + } + Long64_t nentries = fChain->GetEntriesFast(); + if (Nevents > 0 && nentries > Nevents) nentries = Nevents; + Long64_t nbytes = 0, nb = 0; + Int_t TreeNo = -1; + TString currentFile(""); + + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) { + cout << "Read event \t" << jentry + << " so far, switch to file " << fChain->GetCurrentFile()->GetName() + << endl; + cout << " current TreeNo: " << TreeNo + << " new TreeNo: " << fChain->GetTreeNumber() << endl; + TreeNo = fChain->GetTreeNumber(); + } + // if (VertexZCut > 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue; + UInt_t Ntrack = fEvent->GetTracks()->GetEntriesFast(); + // int k_used[100000] = {0}; + for (UInt_t trk = 0; trk < Ntrack; trk++) { + TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); + if (! track) continue; + Int_t Npoints = track->GetNpoint(); +#if 0 + if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue; + if (UseSsd && Npoints < 1000) continue; + if (UseSvt && Npoints < 100) continue; +#endif + Int_t Nsp = track->GetN(); + double dvmin = 1000.; + double dumin = 1000.; + int kmin; + for (Int_t hit = 0; hit < Nsp; hit++) { + Int_t k = track->GetHitTId(hit) - 1; + // assert(k>=0); + HitT *hitT = (HitT *) fEvent->GetHitT(k); + if ( k < 0) cout <<" k <0:"<array with all tracks + TClonesArray *fHits; //->array with all hits + Bool_t fIsValid; // + + static TClonesArray *fgTracks; + static TClonesArray *fgHits; + static THashList *fRotList; + + public: + EventT(); + virtual ~EventT(); + Int_t Build(StEvent *pEventT, Double_t pCut = 0.2); + void Clear(Option_t *option =""); + Bool_t IsValid() const { return fIsValid; } + static void Reset(Option_t *option =""); + void SetNtrack(UInt_t n) { fNtrack = n; } + void SetNhit(UInt_t n) { fNhit = n; } + void SetFlag(UInt_t f) { fFlag = f; } + void SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field); + TrackT *AddTrackT(); + HitT *AddHitT(); + HitT *SetHitT(HitT *h, StHit *hit, TGeoHMatrix *comb, TrackT *track); + Double32_t GetVertex(UInt_t i=0) {return (i<3)?fVertex[i]:0;} + UInt_t GetTotalNoTracks() const {return fNPTracks;} + UInt_t GetNtrack() const { return fNtrack; } + UInt_t GetNhit() const { return fNhit; } + UInt_t GetFlag() const { return fFlag; } + EventTHeader *GetHeader() { return &fEvtHdr; } + const Double32_t *GetVertex() const {return fVertex;} + const Double32_t *GetCovMatrix() const {return fCovariantMatrix;} + TClonesArray *GetTracks() const {return fTracks;} + TClonesArray *GetHits() const {return fHits;} + TrackT *GetTrackT(UInt_t i=0) const {return fTracks && i < fNtrack ? (TrackT*) fTracks->At(i): 0;} + HitT *GetHitT(UInt_t i=0) const {return fHits && i < fNhit ? (HitT*) fHits->At(i): 0;} + Int_t GetIndexOfTrackT(const TrackT *obj) const {return fgTracks->IndexOf(obj);} + Int_t GetIndexOfHitT(const HitT *obj) const {return fgHits->IndexOf(obj);} + static void SetRotMatrices(THashList *Rot) {fRotList = Rot;} + static void RestoreListOfRotations(); + static THashList *RotMatrices() {return fRotList;} + virtual void Print(Option_t *opt="") const; + ClassDef(EventT,1) //EventT structure +}; +class TBase { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + EventT *fEvent; + TString fOutFileName; + TBase(TTree *tree = 0) : fEvent(0) { + // if parameter tree is not specified (or zero), connect the file + // used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("/star/data09/calib/fisyak/Pass112/TpcSsd/065/Event_6065045_raw_1010001.root"); + if (!f) { + f = new TFile("/star/data09/calib/fisyak/Pass112/TpcSsd/065/Event_6065045_raw_1010001.root"); + } + tree = (TTree*)gDirectory->Get("T"); + } + Init(tree); + } + virtual ~TBase() {if (!fChain) return; delete fChain->GetCurrentFile();} + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry) { if (!fChain) return 0; return fChain->GetEntry(entry);} + virtual Long64_t LoadTree(Long64_t entry) { + // Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (!fChain->InheritsFrom(TChain::Class())) return centry; + TChain *chain = (TChain*)fChain; + if (chain->GetTreeNumber() != fCurrent) { + fCurrent = chain->GetTreeNumber(); + Notify(); + } + return centry; + } + virtual void Init(TTree *tree) { + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + fEvent = new EventT(); + TBranch *branch = fChain->GetBranch("EventT"); + branch->SetAddress(&fEvent); + } + virtual void Loop() {Loop(0);} + virtual void Loop(Int_t Nevents); + virtual void SetOutFileName(const Char_t *name="Out.root") {fOutFileName = name;} + virtual Bool_t Notify() {return kTRUE;} + virtual void Show(Long64_t entry = -1) {if (!fChain) return; fChain->Show(entry);} +}; +#endif diff --git a/StRoot/StGmtAligner/HitT.h b/StRoot/StGmtAligner/HitT.h new file mode 100644 index 00000000000..09f451413d1 --- /dev/null +++ b/StRoot/StGmtAligner/HitT.h @@ -0,0 +1,160 @@ +#ifndef __HIT__ +#define __HIT__ +//#define __USE_GLOBAL__ +#include +#include "TObject.h" +class StHit; +class HitT : public TObject { + private: + Char_t start; + Int_t Id; + Int_t sector, barrel, layer, ladder, wafer, hybrid, rdo; // SSD: barrel = layer = hybrid = 0 + Double32_t xG, yG, zG; // hit Global from StEvent + Double32_t xGC, yGC, zGC; // hit Global from local + Double32_t xL, yL, zL; // hit in Ladder CS + Double32_t u, v, w; // hit in Local (Wafer) xL == u_m, yL == v_m + Double32_t tuP, tvP; // tangs + Double32_t uP, vP; // prediction in Wafer CS + Double32_t pT, pMom; // track + Double32_t xPG, yPG, zPG; // Prediction in Global CS + Double32_t cxPG, cyPG, czPG; // Predicted direction cos in Global + Double32_t wGu, wGv, wGw; // Global direction for detector plane + Double32_t xPL, yPL, zPL; // Ladder + Bool_t isPrimary; // If primary track + Bool_t isCrossingMembrain; //If track crosses membrain +#ifdef __USE_GLOBAL__ + Double32_t uPGl, vPGl; // prediction in Wafer CS + Double32_t tuPGl, tvPGl; // tangs + Double32_t xPGlG, yPGlG, zPGlG; // Prediction in Global CS + Double32_t cxPGlG, cyPGlG, czPGlG; // Predicted direction cos in Global + Double32_t xPGlL, yPGlL, zPGlL; // Ladder +#endif + // Double32_t cxPL, cyPL, czPL; // Ladder +#if 0 + Bool_t fValidDerivatives; + Double32_t duPdxV,duPdyV,duPdzV,duPddip,duPdphi,duPdRho; // derivatives uP wrt xV,yV,zV,1/pT,dip,phi + Double32_t dvPdxV,dvPdyV,dvPdzV,dvPddip,dvPdphi,dvPdRho; // derivatives vP wrt xV,yV,zV,1/pT,dip,phi +#endif + Double32_t uM, vM; + Double32_t anode, timeb; + Int_t NoHitPerTrack; + Double32_t uD, vD; // positions of hits from detector + Double32_t duD, dvD; // errors in positions + Double32_t suD, svD; // sigma + Double32_t dsuD, dsvD; // errors in sigma + Double32_t uAdcD, vAdcD; // Adc + Double32_t duAdcD, dvAdcD; // errors in adc + Double32_t uHat; + Double32_t vHat; + Int_t NofHits; // total no. of hits per wafer + Int_t NofFHits;// total no. of fitted hits per wafer + Int_t isFitted; + Int_t isTrack; + Int_t isUsedInFit; + UInt_t hitFlag; + Double32_t sLength; + Double32_t sLengthR; + Double32_t dR; + Char_t end; + public: + HitT(Int_t B = 0, Int_t L = 0, Int_t l = 0, Int_t W = 0, Int_t H = 0, + Double32_t X = 0, Double32_t Y = 0, Double32_t Z = 0, + Double32_t XL = 0, Double32_t YL = 0, Double32_t ZL = 0) { + memset(&start, 0, &end - &start); + SetId(B,L,l,W,H); Set(X,Y,Z,XL,YL,ZL); + } + + virtual ~HitT() {} + void Set(Double32_t X, Double32_t Y, Double32_t Z, + Double32_t XL = 0, Double32_t YL = 0, Double32_t ZL = 0) { + xG = X; yG = Y; zG = Z; + uM = XL; vM = YL; w = ZL; + } + void SetHitLength(Double_t sh) {sLength = sh;} + void SetHitLengthR(Double_t sh) {sLengthR = sh;} + void SetHitdR(Double_t sh) {dR = sh;} + void SetHitFlag(const UInt_t flag) {hitFlag = flag;} + void SetL(Double32_t X, Double32_t Y, Double32_t Z) {xL = X; yL = Y; zL = Z;} + void SetGC(Double32_t X, Double32_t Y, Double32_t Z) {xGC = X; yGC = Y; zGC = Z;} + void SetLM(Double32_t X, Double32_t Z) {u = X; v = Z;} + void SetAnode(Double32_t p=0) {anode=p;} + void SetTimeB(Double32_t p=0) {timeb=p;} + void SetId(Int_t B = 0, Int_t L = 0, Int_t l = 0, Int_t W = 0, Int_t H = 0); +#if 0 + void SetId(StHit *HitT); +#endif + void Set(Double32_t *xyzG, Double32_t *xyzL) {Set(xyzG[0],xyzG[1],xyzG[2],xyzL[0],xyzL[1],xyzL[2]);} + void SetpT(Double32_t p) {pT = p;} + void SetMom(Double32_t p) {pMom = p;} + void SetUVPred(Double32_t u, Double32_t v) {uP = u; vP = v;} + void SettUVPred(Double32_t tu, Double32_t tv) {tuP = tu; tvP = tv;} + void SetWG(Double32_t wu, Double32_t wv, Double32_t ww) { wGu = wu; wGv = wv; wGw = ww;} + // void SetWL(Double32_t wx, Double32_t wy, Double32_t wz) { wLx = wx; wLy = wy; wLz = wz;} +#if 0 + void SetValidDerivatives(Bool_t p=kTRUE) {fValidDerivatives = p;} + void SetDerivatives(Double32_t *der) {Double32_t *d = &duPdxV; for (Int_t i = 0; i < 12; i++) d[i] = der[i];} +#endif + void SetXyzG(const Double_t *x) {Double32_t *xyzPG = &xPG; for (Int_t i = 0; i < 3; i++) xyzPG[i] = x[i];} + void SetDirG(const Double_t *x) {Double32_t *dirPG = &cxPG;for (Int_t i = 0; i < 3; i++) dirPG[i] = x[i];} + void SetXyzL(const Double_t *x) {Double32_t *xyzPL = &xPL; for (Int_t i = 0; i < 3; i++) xyzPL[i] = x[i];} + void SetRDO(Int_t r) {rdo = r;} + void SetuvD(Double32_t u, Double32_t v) {uD = u; vD = v;} + void SetuvDError(Double32_t du, Double32_t dv) {duD = du; dvD = dv;} + void SetSigma(Double32_t su, Double32_t sv) {suD = su; svD = sv;} + void SetSigmaError(Double32_t dsu, Double32_t dsv) {dsuD = dsu; dsvD = dsv;} + void SetAdc(Double32_t uAdc, Double32_t vAdc) {uAdcD = uAdc; vAdcD = vAdc;} + void SetAdcError(Double32_t duAdc, Double32_t dvAdc) {duAdcD = duAdc; dvAdcD = dvAdc;} +#ifdef __USE_GLOBAL__ + void SetUVPredGl(Double32_t u, Double32_t v) {uPGl = u; vPGl = v;} + void SettUVPredGl(Double32_t tu, Double32_t tv) {tuPGl = tu; tvPGl = tv;} + void SetXyzGl(const Double_t *x) {Double32_t *xyzPG = &xPGlG; for (Int_t i = 0; i < 3; i++) xyzPG[i] = x[i];} + void SetDirGl(const Double_t *x) {Double32_t *dirPG = &cxPGlG;for (Int_t i = 0; i < 3; i++) dirPG[i] = x[i];} + void SetXyzGlL(const Double_t *x) {Double32_t *xyzPL = &xPGlL; for (Int_t i = 0; i < 3; i++) xyzPL[i] = x[i];} +#endif + void SetHitPerTrack(Int_t k) {NoHitPerTrack = k;} + void SetuHat(Double_t u) {uHat = u;} + void SetvHat(Double_t v) {vHat = v;} + void SetNofHits(Int_t n) {NofHits = n;} + void SetNofFHits(Int_t n) {NofFHits = n;} + void SetisFitted(Int_t k=1) {isFitted = k;} + void SetisTrack(Int_t k=1) {isTrack = k;} + void SetUsedInFit(Int_t k=0) {isUsedInFit = k;} + void SetisPrimary(Bool_t k) {isPrimary = k;} + void SetisCrossingMembrain(Bool_t k) {isCrossingMembrain = k;} + UInt_t GetHitFlag() const {return hitFlag;} + Double32_t GetU() const {return u;} + Double32_t GetV() const {return v;} + Double32_t GetuD() const {return uD;} + Double32_t GetvD() const {return vD;} + Double32_t *GetXyzP() {return &xPG;} + Double32_t *GetXyzL() {return &xPL;} + Double32_t *GetXyzW() {return &xPL;} + Double32_t GetPredtU() const {return tuP;} + Double32_t GetPredtV() const {return tvP;} + Double32_t GetPredU() const {return uP;} + Double32_t GetPredV() const {return vP;} + Bool_t GetisPrimary() const {return isPrimary;} + Bool_t GetisCrossingMembrain(){return isCrossingMembrain;} +#ifdef __USE_GLOBAL__ + + Double32_t *GetXyzPGl() {return &xPGlG;} + Double32_t *GetXyzLGl() {return &xPGlL;} + Double32_t *GetXyzWGl() {return &xPGlL;} + Double32_t GetPredGltU() const {return tuPGl;} + Double32_t GetPredGltV() const {return tvPGl;} + Double32_t GetPredGlU() const {return uPGl;} + Double32_t GetPredGlV() const {return vPGl;} +#endif + Int_t Barrel() const {return barrel;} + Int_t Layer() const {return layer;} + Int_t Ladder() const {return ladder;} + Int_t Wafer() const {return wafer;} + Int_t GetId() const {return Id;} +#if 0 + const Double32_t *GetDerivatives() const {return &duPdxV;} + Bool_t ValidDerivatives() const {return fValidDerivatives;} +#endif + virtual void Print(Option_t *opt="") const; + ClassDef(HitT,1) +}; +#endif diff --git a/StRoot/StGmtAligner/StGmtAligner.cxx b/StRoot/StGmtAligner/StGmtAligner.cxx new file mode 100644 index 00000000000..757ca7eeab8 --- /dev/null +++ b/StRoot/StGmtAligner/StGmtAligner.cxx @@ -0,0 +1,150 @@ +#include "StGmtAligner.h" +#include +#include "Riostream.h" +#include "TROOT.h" +#include "TSystem.h" +#include "TFile.h" +#include "TKey.h" +#include "TRandom.h" +#include "TTree.h" +#include "TBranch.h" +#include "TStopwatch.h" +#include "StThreeVectorF.hh" +#include "StMatrixF.hh" +#include "TH1.h" +#include "TH2.h" +#include "TProfile.h" +#include "TMath.h" +#include "TVector3.h" +#include "TProcessID.h" +#include "StEvent.h" +#include "StPrimaryVertex.h" +#include "StBFChain.h" +#include "TGeoMatrix.h" +#include "EventT.h" +#include "StTpcDb/StTpcDb.h" +#include "StDetectorDbMaker/StGmtSurveyC.h" +//________________________________________________________________________________ +StGmtAligner::StGmtAligner(const Char_t *name) : StMaker(name),fFile(0), fTree(0), fEvent(0) { + SetMinNoHits(); + SetpCut(); + SetOut(); +} +//________________________________________________________________________________ +Int_t StGmtAligner::Init() { + SetTree(); + return kStOK; +} +//________________________________________________________________________________ +Int_t StGmtAligner::Finish() { + if (fFile) { + fFile = fTree->GetCurrentFile(); //just in case we switched to a new file + fFile->Write(); + fTree->Print(); + } + return kStOK; +} + +//________________________________________________________________________________ +void StGmtAligner::SetTree() { + StBFChain *chain = (StBFChain *) StMaker::GetChain(); + if (! chain) return; + // root.exe + Int_t split = 99; // by default, split Event in sub branches + Int_t comp = 1; // by default file is compressed + Int_t branchStyle = 1; //new style by default + if (split < 0) {branchStyle = 0; split = -1-split;} + Int_t bufsize; + //Authorize Trees up to 2 Terabytes (if the system can do it) + TTree::SetMaxTreeSize(1000*Long64_t(2000000000)); + TFile *f = GetTFile(); + if (f) f->cd(); + else { + TString FName(fOut); + if (fMinNoHits > 0) FName += Form("_%i_%f2.1_",fMinNoHits,fpCut); + FName += gSystem->BaseName(chain->GetFileIn().Data()); + FName.ReplaceAll("st_physics",""); + FName.ReplaceAll(".event",""); + FName.ReplaceAll(".daq",".root"); + fFile = new TFile(FName,"RECREATE","TTree with SVT + SSD hits and tracks"); + fFile->SetCompressionLevel(comp); + } + // Create a ROOT Tree and one superbranch + fTree = new TTree("T","TTree with SVT + SSD hits and tracks"); + fTree->SetAutoSave(1000000000); // autosave when 1 Gbyte written + bufsize = 64000; + if (split) bufsize /= 4; + fEvent = new EventT(); + TTree::SetBranchStyle(branchStyle); + TBranch *branch = fTree->Branch("EventT", &fEvent, bufsize,split); + branch->SetAutoDelete(kFALSE); +} +//________________________________________________________________________________ +Int_t StGmtAligner::Make() { + if (! EventT::RotMatrices()) MakeListOfRotations(); +// StEvent* pEvent = (StEvent*) GetInputDS("StEvent"); + StEvent* pEvent = (StEvent*) StMaker::GetChain()->GetInputDS("StEvent"); + cout << "ALARM!!!!" << endl; + cout << "WE ARE FEELING THE TREE!!!! " << endl; + cout << "LOOOK!!!!! I'M HERE!!!!" << endl; + if (pEvent && !fEvent->Build(pEvent,fpCut)) fTree->Fill(); //fill the tree + return kStOK; +} +//________________________________________________________________________________ +void StGmtAligner::Print(Option_t *opt) const { + if (! EventT::RotMatrices()) return; + TIter next(EventT::RotMatrices()); + TGeoHMatrix *comb = 0; + while ((comb = (TGeoHMatrix *) next())) { + Int_t Id; + sscanf(comb->GetName()+1,"%04i",&Id); + Int_t Ladder = Id%100; + Int_t Layer = Id/1000; if (Layer > 7) Layer = 7; + Int_t Wafer = (Id - 1000*Layer)/100; + cout << comb->GetName() << "\tLayer/Ladder/Wafer = " << Layer << "/" << Ladder << "/" << Wafer << endl; + comb->Print(); + cout << "=================================" << endl; + } +} +//________________________________________________________________________________ +void StGmtAligner::MakeListOfRotations() { + if (EventT::RotMatrices()) return; + THashList *rotMHash = new THashList(100,0); + EventT::SetRotMatrices(rotMHash); + THashList *hash = 0; + const TGeoHMatrix& tpc2Glob = gStTpcDb->Tpc2GlobalMatrix(); + for(int module=0;module< kGmtNumModules;module++) { + TGeoHMatrix *WL = new TGeoHMatrix(StGmtOnModule::instance()->GetMatrix(module)); + WL->SetName(Form("WL%i",module)); + rotMHash->Add(WL); + TGeoHMatrix GmtOnGlob = tpc2Glob * StGmtOnTpc::instance()->GetMatrix(module) * (*WL); + TGeoHMatrix *R = new TGeoHMatrix(GmtOnGlob); + R->SetName(Form("R%i",module)); + rotMHash->Add(R); + } + TIter next(rotMHash); + TGeoHMatrix *comb; + Int_t fail = 0; + while ((comb = (TGeoHMatrix *) next())) { + TString Name(comb->GetName()); + if (Name.BeginsWith("R")) { + TGeoHMatrix *WL = (TGeoHMatrix *) rotMHash->FindObject(Form("WL%s",Name.Data()+1)); + if (! WL) { + cout << Form("WL%s",Name.Data()+1) << " has not been found" << endl; + fail++; + } + } + } + assert(! fail); +#if 0 + if (fFile) { + TDirectory *g = 0; + if (gDirectory != fFile) { + g = gDirectory; + fFile->cd(); + } + rotMHash->Write(); + if (g) g->cd(); + } +#endif +} diff --git a/StRoot/StGmtAligner/StGmtAligner.h b/StRoot/StGmtAligner/StGmtAligner.h new file mode 100644 index 00000000000..7389cc53926 --- /dev/null +++ b/StRoot/StGmtAligner/StGmtAligner.h @@ -0,0 +1,36 @@ +#ifndef __StGmtAligner_H +#define __StGmtAligner_H +#include "StMaker.h" +#include "TFile.h" +#include "TArrayI.h" +#include "TTree.h" + +//________________ +class EventT; + +//________________ +class StGmtAligner : public StMaker { + public: + StGmtAligner(const Char_t *name="GmtAligner"); + virtual ~StGmtAligner() {} + virtual Int_t Init(); + virtual Int_t Make(); + virtual Int_t Finish(); + void SetTree(); + void Print(Option_t *opt="") const; + void SetMinNoHits(Int_t MinNoHits=0) {fMinNoHits = MinNoHits;} + void SetpCut(Double_t pCut=0.0) {fpCut = pCut;} + void SetOut(Char_t *Out="Event") {fOut = Out;} + void MakeListOfRotations(); + virtual const char *GetCVS() const + {static const char cvs[]="Tag $Name: $ $Id: StGmtClusterMaker.h,v 1.1.1.1 2013/09/02 15:01:02 fisyak Exp $ built " __DATE__ " " __TIME__ ; return cvs;} + private: + TFile *fFile; + TTree *fTree; + EventT *fEvent; + Int_t fMinNoHits; + Double_t fpCut; + Char_t *fOut; + ClassDef(StGmtAligner,1) +}; +#endif diff --git a/StRoot/StGmtAligner/TrackT.h b/StRoot/StGmtAligner/TrackT.h new file mode 100644 index 00000000000..61848cc88a7 --- /dev/null +++ b/StRoot/StGmtAligner/TrackT.h @@ -0,0 +1,78 @@ +#ifndef __TRACK__ +#define __TRACK__ +#define NSP 1000 +#include +#include "TObject.h" +#include "TMath.h" +#include "HitT.h" + +//________________ +class TrackT : public TObject { + + private: + Char_t beg; + Double32_t fInvpT; //signed + Double32_t fTanL; + Double32_t fPhi; + Double32_t fRho; +#ifdef __USE_GLOBAL__ + Double32_t fInvpTGl; //signed + Double32_t fTanLGl; + Double32_t fPhiGl; + Double32_t fRhoGl; +#endif + UInt_t fNpoint; //Number of fitted points for this track + UInt_t fNPpoint; //Number of possible points for this track + Short_t fValid; //Validity criterion + UInt_t fNsp; //Number of points for this track with a special value + UInt_t fIdHitT[NSP]; //Index of HitT in fHitT array + Double32_t fdEdx; + Double32_t fLength; + Char_t end; + public: + TrackT() { Clear(); } + virtual ~TrackT() {Clear();} + void Clear(Option_t *option="") {if (option); memset(&beg, 0, &end - &beg);} + Double32_t GetpX() { return GetpT()*TMath::Cos(fPhi);} + Double32_t GetpY() { return GetpT()*TMath::Sin(fPhi);} + Double32_t GetpZ() { return GetpT()*fTanL;} + Double32_t GetInvpT() { return fInvpT;} + Double32_t GetTanL() { return fTanL;} + Double32_t GetDip() { return TMath::ATan(fTanL);} + Double32_t GetPhi() { return fPhi;} + Double32_t GetRho() { return fRho;} + Double32_t GetpT() { return TMath::Abs(fInvpT) > 1.e-7 ? 1./TMath::Abs(fInvpT): 1e7; } + Double32_t GetMomentum(){ return GetpT()*TMath::Sqrt(1. + fTanL*fTanL);} + UInt_t GetNpoint() { return fNpoint; } + UInt_t GetNPpoint() { return fNPpoint; } + Short_t GetCharge() { return (Short_t) TMath::Sign(1., fInvpT); } + Short_t GetValid() { return fValid; } + UInt_t GetN() { return fNsp; } + const UInt_t *GetIndx() const { return fIdHitT;} + Int_t GetHitTId(UInt_t i=0) {return i < fNsp ? ((Int_t) fIdHitT[i])-1 : -1;} + Double32_t GetdEdx() {return fdEdx;} + Double32_t GetLegth() {return fLength;} + virtual void SetInvpT(Double32_t p) {fInvpT = p; } + virtual void SetDip(Double32_t p) {fTanL = TMath::Tan(p); } + virtual void SetTanL(Double32_t p) {fTanL = p; } + virtual void SetPhi(Double32_t p) {fPhi = p; } + virtual void SetRho(Double32_t p) {fRho = p; } +#ifdef __USE_GLOBAL__ + + virtual void SetInvpTGl(Double32_t p) {fInvpTGl = p; } + virtual void SetDipGl(Double32_t p) {fTanLGl = TMath::Tan(p); } + virtual void SetTanLGl(Double32_t p) {fTanLGl = p; } + virtual void SetPhiGl(Double32_t p) {fPhiGl = p; } + virtual void SetRhoGl(Double32_t p) {fRhoGl = p; } +#endif + virtual void SetNpoint(UInt_t p) {fNpoint = p; } + virtual void SetNPpoint(UInt_t p) {fNPpoint = p; } + virtual void SetValid(Short_t p=1) {fValid = p; } + virtual void SetLength(Double_t L) {fLength = L;} + virtual void SetN(UInt_t n) {if (n <= NSP) fNsp = n; else fNsp = NSP;} + virtual void SetHitTId(UInt_t i) {fIdHitT[fNsp] = i+1; if ( fNsp < NSP) fNsp++;} + virtual void SetdEdx(Double_t I70, Double_t L) {fdEdx = I70; fLength = L;} + virtual void Print(Option_t *opt="") const; + ClassDef(TrackT,2) +}; +#endif diff --git a/StRoot/StGmtAligner/runGmtTree.C b/StRoot/StGmtAligner/runGmtTree.C new file mode 100644 index 00000000000..3ef340e8e85 --- /dev/null +++ b/StRoot/StGmtAligner/runGmtTree.C @@ -0,0 +1,23 @@ +void runGmtTree(const Char_t *input, const Char_t *output=0) { + gROOT->LoadMacro("bfc.C"); + Load(); + TString Input(input); + TString Chain("in,StEvent,mysql,libPhysics,db,StarMagField,MagF,svtDb,ssdDb,GmtAligner,nodefault"); + if (Input.Contains("rcf") || Input.Contains("gstar")) Chain += ",y2005e,simu"; + cout << "Chain used:\t" << Chain << endl; + bfc(-1,Chain.Data(),input,0,output); + StMaker *db = chain->Maker("db"); + if (db) db->SetDebug(1); + StGmtAligner *mtree = (StGmtAligner *) chain->Maker("GmtAligner"); + if (! mtree) *mtree = new StGmtAligner(); + if (! mtree) return; +#if 0 + mtree->SetMinNoHits(2); + mtree->SetpCut(0.2); +#else + mtree->SetMinNoHits(0); + mtree->SetpCut(0); +#endif + chain->Init(); + chain->EventLoop(9999); +} diff --git a/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx new file mode 100644 index 00000000000..180b4455366 --- /dev/null +++ b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx @@ -0,0 +1,302 @@ +// +// First Cluster Maker +// \class StGmtClusterMaker +// \authors K.S. Engle and Richard Witt (witt@usna.edu) +// based on StFgtClusterMaker + +#include "StGmtClusterMaker.h" +#include "StRoot/StEvent/StEvent.h" +#include "StRoot/StEvent/StGmtCollection.h" +//#include "StRoot/StEvent/StTpcHitCollection.h" +#include "StEvent/StGmtHit.h" +#include "StRoot/StGmtUtil/geometry/StGmtGeom.h" +#include "TSystem.h" +#include "TSpectrum.h" +#include "StMessMgr.h" +ClassImp(StGmtClusterMaker); + +Int_t StGmtClusterMaker::gmtStat = 0; +const UInt_t CLUS_BINS = 128; +const Double_t CLUS_MIN = 0.0; +const Double_t CLUS_MAX = 128*0.08; +const UInt_t MAX_PEAKS = 20; + +//________________________________________________________________________________ +inline Double_t MyGaus(Double_t x, Double_t mean, Double_t sigma, Double_t delta) { + return TMath::Freq((x-mean+delta/2)/sigma)-TMath::Freq((x-mean-delta/2)/sigma); +} +//________________________________________________________________________________ +Double_t fpeaks(Double_t *x, Double_t *par) { + Float_t result=0.0; + UInt_t nPar=(UInt_t)par[0]; + for (UInt_t p=0;pFit(&back,"Q"); + UInt_t npx=0; + UInt_t nfound = spect.GetNPeaks(); + for(UInt_t i=0; i < nfound; i++) { + Double_t xp=xpeaks[i]; + int bin=hist->GetXaxis()->FindBin(xp); + Double_t yp=hist->GetBinContent(bin); + Double_t err=hist->GetBinError(bin); + if(err<=0.0) continue; + if(bin<=1) continue; + if((yp-err*3) < back.GetParameter(0)) continue; + Double_t yp_left=hist->GetBinContent(bin-1); + Double_t yp_right=hist->GetBinContent(bin+1); + Double_t err_left=hist->GetBinError(bin-1); + Double_t err_right=hist->GetBinError(bin+1); + Double_t yp_sum=yp+yp_left+yp_right; + Double_t err_sum=TMath::Sqrt(err*err+err_left*err_left+err_right*err_right); + if((yp_sum-3*err_sum) < back.GetParameter(0)) continue; + + par[3*npx+1]=TMath::Log(yp); + par[3*npx+2]=xp; + par[3*npx+3]=3*0.08; // sigma + npx++; + } + if (Debug()) {LOG_INFO << hist->GetName() << " found " << nfound << " Accpeted " << npx << endm;} + if (! npx) return 0; + + TString funcName=Form("Func%s",hist->GetName()); + TF1* fitFunc; + if(fitFunc=(TF1*)gROOT->GetListOfFunctions()->FindObject(funcName)) delete fitFunc; + fitFunc=new TF1(funcName,fpeaks,CLUS_MIN,CLUS_MAX,3*npx+1); + + for(UInt_t i=0; i < npx; i++) {fitFunc->SetParLimits(3*i+3,0.08*0.5,10*0.08);} + fitFunc->SetParameters(par); + fitFunc->FixParameter(0,(double)npx); + fitFunc->SetNpx(1000); + fitFunc->SetLineColor(kGreen); + + TVirtualFitter::SetDefaultFitter("Fumili"); + int isOk=hist->Fit(fitFunc); + if(isOk) isOk=hist->Fit(fitFunc); + if(isOk) return 0; + + return fitFunc; +} +//________________________________________________________________________________ +void StGmtClusterMaker::ClusterBuilder(ULong_t events, UInt_t module, StGmtStripCollection& strips, StGmtHitCollection& hits) { + static TCanvas* canv=0; + static TH1F* histX=0; + static TH1F* histY=0; + static TProfile* profX[8]={0}; + static TProfile* profY[8]={0}; + + StGmtStrip* pStrip; + Float_t position; + UInt_t stripsNum=strips.getNumStrips(); + int adc,adc_buf=0; + TString name, title; +// LOG_INFO << "STart TSpectrum" << endm; +// TSpectrum spectX(MAX_PEAKS); TSpectrum spectY(MAX_PEAKS); + TSpectrum spectX(); TSpectrum spectY(); + LOG_INFO << "Created TSpectrum" << endm; + TH1F* histPointer; + TProfile* profPointer; + + if(!profX[module]) { + name="PedestalX_"; name += module; + profX[module]=new TProfile(name,name,CLUS_BINS,CLUS_MIN,CLUS_MAX,"s"); + } + if(!profY[module]) { + name = "PedestalY_"; name += module; + profY[module]=new TProfile(name,name,CLUS_BINS,CLUS_MIN,CLUS_MAX,"s"); + } + if(!histX) histX=new TH1F("ClusterX","ClusterX",CLUS_BINS,CLUS_MIN,CLUS_MAX); + if(!histY) histY=new TH1F("ClusterY","ClusterY",CLUS_BINS,CLUS_MIN,CLUS_MAX); + + if(Debug()>3) { + canv = (TCanvas *) gROOT->GetListOfCanvases()->FindObject("GmtClusters"); + if(!canv) canv=new TCanvas("GmtClusters","GmtClusters",768,768); + else canv->Clear(); + } + + histX->Reset(); histY->Reset(); + TProfile profXold(*profX[module]); TProfile profYold(*profY[module]); + for(UInt_t iStrip=0;iStripisY()) { profPointer=profX[module]; histPointer=histX; } + else { profPointer=profY[module]; histPointer=histY; } + for(UInt_t iTimeBin=0;iTimeBingetAdc(iTimeBin); + if(adc_buf>-999) adc+=adc_buf; + } + position=pStrip->getPosition(); + int bin=profPointer->Fill(position,adc); + Float_t error=TMath::Sqrt(adc); + histPointer->Fill(position,adc); + histPointer->SetBinError(bin,error); + } + if (events < 5) return; + histX->Add(&profXold,-1.0); histY->Add(&profYold,-1.0); + + TF1 *fitX=0, *fitY=0; + fitX = FindPeaks(histX); fitY = FindPeaks(histY); + UInt_t idx[MAX_PEAKS], idy[MAX_PEAKS]; + UInt_t nClusX=0, nClusY=0; + if(fitX) { + for(UInt_t i=0; iGetParameter(0); i++) { + if (fitX->GetParameter(3*i+3) > 0.4) continue; + if (fitX->GetParameter(3*i+1) < 5.0) continue; + idx[nClusX] = i; + nClusX++; + } + if (nClusX) *profX[module]=profXold; + if (Debug()) {LOG_INFO << "######XPEAKS found =" << fitX->GetParameter(0) << ", Clusters fitted =" << nClusX << endm;} + } else if (Debug()) {LOG_INFO << "######XNULL" << endm; } + if(fitY) { + for(UInt_t i=0; iGetParameter(3*i+3); i++) { + if (fitY->GetParameter(3*i+3) > 0.4) continue; + if (fitY->GetParameter(3*i+1) < 5.0) continue; + idy[nClusY] = i; + nClusY++; + } + if (nClusY) *profY[module]=profYold; + if (Debug()) {LOG_INFO << "######YPEAKS found =" << fitY->GetParameter(0) << ", Clusters fitted =" << nClusY << endm;} + } else if (Debug()) {LOG_INFO << "######YNULL" << endm; } + for(UInt_t i=0; i < nClusX; i++) { + UInt_t nx = idx[i]; + for(UInt_t j = 0; j < nClusY; j++) { + UInt_t ny = idy[j]; + StGmtHit* newCluster = new StGmtHit( + hits.getHitVec().size(), + module, + TMath::Exp(fitX->GetParameter(3*nx+1)), // adcX + TMath::Exp(fitY->GetParameter(3*ny+1)), // adcY + fitX->GetParError(3*nx+1), // error(adcX) + fitY->GetParError(3*ny+1), // error(adcY) + fitX->GetParameter(3*nx+2), // meanX + fitY->GetParameter(3*ny+2), // meanY + fitX->GetParError(3*nx+2), // error(meanX) + fitY->GetParError(3*ny+2), // error(meanY) + fitX->GetParameter(3*nx+3), // sigmaX + fitY->GetParameter(3*ny+3), // sigmaY + fitX->GetParError(3*nx+3), // error(sigmaX) + fitY->GetParError(3*ny+3)); // error(sigmaY) + if (Debug()) newCluster->Print(); + hits.getHitVec().push_back(newCluster); + } + } + + if (Debug()>3) { + canv->Divide(2,2); + + canv->cd(1); + profX[module]->Draw(); + canv->cd(2); + profY[module]->Draw(); + + canv->cd(3); + histX->Draw(); + canv->cd(4); + histY->Draw(); + + canv->Modified(); + canv->Update(); + canv->Draw(); + if (nClusX || nClusY) { + while (!gSystem->ProcessEvents()){gSystem->Sleep(200);} + } + } +} +//________________________________________________________________________________ +StGmtClusterMaker::StGmtClusterMaker( const Char_t* name ) : //StMaker(name), + StRTSBaseMaker( "clustser", name ) { + SetAttr("gmtCosmics" ,kFALSE); +}; +//________________________________________________________________________________ +Int_t StGmtClusterMaker::Make() { + LOG_INFO << "MAKE of StGmtClusterMaker" << endm; + Int_t ierr = kStOk; + static ULong_t nEvents=0; + //StEvent* eventPtr = 0; + //eventPtr = (StEvent*)GetInputDS("StEvent"); + StEvent* eventPtr = (StEvent*) (GetInputDS("StEvent")); + //cout << "LLLLLLOOOOOOOOKKK!!!" << endl; + //cout << "TRACK NODES: " << endl; + //cout << eventPtr->trackNodes().size() << endl; + //cout << "TPC HIT COLLECTIONS!" << endl; + //cout << eventPtr->tpcHitCollection()->numberOfHits() << endl; + if(!eventPtr) { + LOG_ERROR << "Error getting pointer to StEvent from '" << ClassName() << "'" << endm; + return kStErr; + } + StGmtCollection* gmtCollectionPtr = eventPtr->gmtCollection(); + if(!gmtCollectionPtr) { + LOG_WARN << "Error getting pointer to StGmtCollection from '" << ClassName() << "'" << endm; + return kStWarn; + } + UInt_t noModWithGMT = 0; + for(UInt_t moduleIdx=0; moduleIdxgetNumModules(); moduleIdx++) { + if(Debug()) { + LOG_INFO << "module: " << moduleIdx << " has strips: \t" << gmtCollectionPtr->getNumStrips(moduleIdx) << endm; + LOG_INFO << "Collection =\t" << gmtCollectionPtr->getNumStrips() << "\t" + << gmtCollectionPtr->getNumHits() << '\t' << gmtCollectionPtr->getNumPoints() << '\t' << endm; + } + Int_t nelements = gmtCollectionPtr->getNumStrips(moduleIdx); + if(nelements < kGmtNumStripsPerModule) { + if(Debug()) { + LOG_WARN <<"StClusterMaker::Make(): no data for module " << moduleIdx << endm; + } + continue; + } + StGmtStripCollection *stripCollectionPtr = gmtCollectionPtr->getStripCollection(moduleIdx); + StGmtHitCollection *hitCollectionPtr = gmtCollectionPtr->getHitCollection(moduleIdx); + ClusterBuilder(nEvents,moduleIdx,*stripCollectionPtr,*hitCollectionPtr); + noModWithGMT++; + if(stripCollectionPtr && hitCollectionPtr && hitCollectionPtr->getHitVec().size()) { + if(Debug()) { + LOG_INFO << "Cluster " << stripCollectionPtr->getNumStrips() << "strips\tin module" << stripCollectionPtr->getModule() << endm; + } + } + } + if (noModWithGMT) nEvents++; + + if(Debug()) { + LOG_INFO << "End of gmt-clust-maker, print all strips & clusters: " << endm; + LOG_INFO <<" gmtCollnumModule=" << gmtCollectionPtr->getNumModules()<<", tot strip=" <getNumStrips() + <<" totClust=" << gmtCollectionPtr->getNumHits() <getNumHits()) return kStERR; + } + if (Debug()) { + UShort_t NumModules = gmtCollectionPtr->getNumModules(); + for (UShort_t m = 0; m < NumModules; m++) { + const StGmtHitCollection *coll = gmtCollectionPtr->getHitCollection(m); + if (! coll) continue; + const StSPtrVecGmtHit &hits = coll->getHitVec(); + UInt_t NoHits = hits.size(); + for (UInt_t l = 0; l < NoHits; l++) { + const StGmtHit *hit = hits[l]; + if (hit) { + hit->Print(""); + } + } + } + } + return ierr; +} +//________________________________________________________________________________ +Int_t StGmtClusterMaker::Init() { + LOG_INFO << "INTI of StGmtClusterMaker" << endm; + if (IAttr("gmtCosmics")) SetAttr(".Privilege",kTRUE); + return kStOk; +} diff --git a/StRoot/StGmtClusterMaker/StGmtClusterMaker.h b/StRoot/StGmtClusterMaker/StGmtClusterMaker.h new file mode 100644 index 00000000000..688f4a7e7d7 --- /dev/null +++ b/StRoot/StGmtClusterMaker/StGmtClusterMaker.h @@ -0,0 +1,44 @@ +// +// First Cluster Maker +// \class StGmtClusterMaker +// \authors K.S. Engle and Richard Witt (witt@usna.edu) +// based on StFgtClusterMaker +#ifndef STAR_StGmtClusterMaker_HH +#define STAR_StGmtClusterMaker_HH +#include "StMaker.h" +#include "TH1.h" +#include "TProfile.h" +#include "TROOT.h" +#include "TCanvas.h" +#include "TPolyMarker.h" +#include "StRoot/StChain/StRTSBaseMaker.h" +#include "Stypes.h" +#include "TSpectrum.h" +#include "TF1.h" +#include "TMath.h" +#include "TVirtualFitter.h" + +class StGmtStripCollection; +class StGmtHitCollection; + +class StGmtClusterMaker : public StRTSBaseMaker { + //omitted assignment operator and copy constructor on purpose + public: + StGmtClusterMaker( const Char_t* name="GmtCluster"); + ~StGmtClusterMaker() {} + Int_t Init(); + Int_t Make(); + /**sets the clustering algorithm. Currently there is the simple Clustering algorithm and the max cluster algorithm. + The simple cluster algorithm is the default one. The max cluster only selects one hit stip per plane, the one with the highest charge + */ + virtual const char *GetCVS() const + {static const char cvs[]="Tag $Name: $ $Id: StGmtClusterMaker.h,v 1.1.1.1 2013/09/02 15:01:02 fisyak Exp $ built " __DATE__ " " __TIME__ ; return cvs;} + static Int_t gmtStat; + protected: + + void ClusterBuilder(ULong_t events, UInt_t module, StGmtStripCollection& strips, StGmtHitCollection& hits); + TF1* FindPeaks(TH1F* hist); + + ClassDef(StGmtClusterMaker,1) +}; +#endif diff --git a/StRoot/StGmtRawMaker/StGmtRawMaker.cxx b/StRoot/StGmtRawMaker/StGmtRawMaker.cxx new file mode 100644 index 00000000000..b72a2783198 --- /dev/null +++ b/StRoot/StGmtRawMaker/StGmtRawMaker.cxx @@ -0,0 +1,189 @@ +// +// \class StGmtRawMaker +// \authors K.S. Engle and Richard Witt (witt@usna.edu) +// based on StFgtRawMaker +// + +#include "St_base/StMessMgr.h" +#include "St_base/Stypes.h" + +#include "StChain/StRtsTable.h" +#include "StEvent/StEvent.h" +#include "DAQ_FGT/daq_fgt.h" +#include "DAQ_READER/daq_dta.h" +// #include "DAQ_GMT/daq_gmt.h" +// #include "DAQ_READER/daq_dta.h" + +#include "StEvent/StGmtCollection.h" +#include "StEvent/StGmtStripCollection.h" +#include "StEvent/StGmtStrip.h" +// #include "StGmtDbMaker/StGmtDbMaker.h" +#include "StGmtUtil/geometry/StGmtGeom.h" +#include "St_base/StMessMgr.h" +#include "St_base/Stypes.h" + +#include "StGmtRawMaker.h" +ClassImp(StGmtRawMaker) + +const Int_t mChIdToSeqId[128] = { + 0,16,32,48,64,80,96,112,4,20,36,52,68,84,100,116,8, + 24,40,56,72,88,104,120,12,28,44,60,76,92,108,124,1, + 17,33,49,65,81,97,113,5,21,37,53,69,85,101,117,9, + 25,41,57,73,89,105,121,13,29,45,61,77,93,109,125,2, + 18,34,50,66,82,98,114,6,22,38,54,70,86,102,118,10, + 26,42,58,74,90,106,122,14,30,46,62,78,94,110,126,3, + 19,35,51,67,83,99,115,7,23,39,55,71,87,103,119,11, + 27,43,59,75,91,107,123,15,31,47,63,79,95,111,127 +}; + +//________________________________________________________________________________ +/** +Function to get pointer to StEvent datastructures. Creates them if they do not exist already. +*/ +Int_t StGmtRawMaker::prepareEnvironment() { + StEvent* eventPtr = (StEvent*)StRTSBaseMaker::GetInputDS("StEvent"); + if (! eventPtr) return kStFatal; + mGmtCollectionPtr = eventPtr->gmtCollection(); + if(!mGmtCollectionPtr) { + mGmtCollectionPtr = new StGmtCollection(); + if(!mGmtCollectionPtr) { + LOG_DEBUG <<"::prepareEnvironment could not create StGmtCollection" <setGmtCollection(mGmtCollectionPtr); + LOG_DEBUG <<"::prepareEnvironment() has added a non existing StGmtCollection()"<Clear(); + } + return kStOK; +} +//________________________________________________________________________________ +/** +Maker main function. Getting pointer to StEvent and fills the event structure +*/ +Int_t StGmtRawMaker::Make() +{ + + LOG_DEBUG <<"StGmtRawMaker::Make()******************************************************************"<begin();it!=rts_tbl->end();it++) { + fgt_adc_t *mGmtRawData=(fgt_adc_t*)*it; + rdo=rts_tbl->Rdo(); + Int_t chanTmp=mGmtRawData->ch; + channel=mChIdToSeqId[chanTmp]; + //this is different from rts_example + timebin=mGmtRawData->tb; + //look at rts_example for the mapping + adc=mGmtRawData->adc; + arm=rts_tbl->Sector(); + apv=rts_tbl->Pad(); + // the next segment is needed because of a lack of ARM port information + if ( (apv >= 0) && (apv <= 3)) { + port = 0; + } else if ( (apv >= 12) && (apv <= 15) ) { + // cout<<"APV "<getStripCollection( moduleIdx ); + if( stripCollectionPtr ) { + geoId = StGmtGeom::encodeGeoId( rdo, arm, apv, channel ); + Int_t elecId = StGmtGeom::getElectIdFromElecCoord( rdo, arm, apv, channel ); + StGmtStrip* stripPtr = stripCollectionPtr->getStrip( geoId ); + if( coordNum == 999 ) { // these are not connected + stripPtr->setAdc( 0, timebin ); + stripPtr->setCharge( 0 ); // was done in separate maker for FGT (StFgtA2CMaker), assume gain=1 for now + stripPtr->setChargeUncert( 0 ); // was done in separate maker for FGT (StFgtA2CMaker), assume gain=1 for now + stripPtr->setGeoId( geoId ); + stripPtr->setModule( moduleIdx ); + stripPtr->setIsY( layer ); + stripPtr->setPosition( position ); + stripPtr->setElecCoords( rdo, arm, apv, channel ); + stripPtr->setCoordNum( coordNum ); + stripPtr->setPed( 0 ); + stripPtr->setPedStdDev( 0 ); + stripPtr->setPedErr( 0 ); + } else {// these are connected (mapping in StGmtGeom.cxx) + if(layer) {// layer here is just an indicator for either a X ( i.e. strip (=0) ) or Y ( i.e. pad (=1) ) element + stripPtr->setCoordNum( coordNum + kGmtNumStrips ); // map Y into 128-255 + //never returns more certain ids + if(channel==100) {LOG_INFO << "Str.=" << channel << "\tLay0\tgeoid=" << geoId<< "\tposition=" << position << endl;} + } else { + stripPtr->setCoordNum( coordNum ); // map X into 0-127 + if(channel==50) LOG_INFO << "Str.=" << channel << "\tLay0\tgeoid=" << geoId<< "\tposition=" << position << endl; + } + stripPtr->setAdc( adc, timebin ); + stripPtr->setCharge( adc ); // was done in separate maker for FGT (StFgtA2CMaker), assume gain=1 for now + stripPtr->setChargeUncert( sqrt(adc) ); // was done in separate maker for FGT (StFgtA2CMaker), assume gain=1 for now + stripPtr->setGeoId( geoId ); + stripPtr->setModule( moduleIdx ); + stripPtr->setIsY( layer ); + stripPtr->setPosition( position ); + stripPtr->setElecCoords( rdo, arm, apv, channel ); + } + if (Debug()) { + LOG_INFO << "StGmtRawMaker::fillHits() Set: " + << *stripPtr << endm; + } + } else { LOG_WARN << "StGmtRawMaker::Make() -- Could not access module " << moduleIdx << endm; } + } + } + return kStOK; +} diff --git a/StRoot/StGmtRawMaker/StGmtRawMaker.h b/StRoot/StGmtRawMaker/StGmtRawMaker.h new file mode 100644 index 00000000000..170d1b706f3 --- /dev/null +++ b/StRoot/StGmtRawMaker/StGmtRawMaker.h @@ -0,0 +1,29 @@ +// \class StGmtRawMaker +// \authors K.S. Engle and Richard Witt (witt@usna.edu) +// based on StFgtRawMaker + +#ifndef STAR_StGmtRawMaker_HH +#define STAR_StGmtRawMaker_HH + +#include +#include "StRoot/StChain/StRTSBaseMaker.h" + +class StGmtCollection; +/** + This is the raw maker for the GMT data. It makes use of its base class functions to read daq + files into the StGmtEvent Data structure. +*/ +class StGmtRawMaker : public StRTSBaseMaker { + public: + StGmtRawMaker(const Char_t* name="GmtRaw") : StRTSBaseMaker( "adc", name ), mGmtCollectionPtr(0) {} + ~StGmtRawMaker() {} + Int_t Make(); + protected: + Int_t fillHits(); + Int_t prepareEnvironment(); + StGmtCollection *mGmtCollectionPtr; + + private: + ClassDef(StGmtRawMaker,1) +}; +#endif diff --git a/StRoot/StGmtUtil/StGmtConsts.h b/StRoot/StGmtUtil/StGmtConsts.h new file mode 100644 index 00000000000..ef1464b2753 --- /dev/null +++ b/StRoot/StGmtUtil/StGmtConsts.h @@ -0,0 +1,60 @@ +/* + * Some basic constants (enums) the GMT. All GMT related constants + * that are not in the database should be in this file, except those + * which cannot be expressed as an int. These remaining constants are + * in the file StGmtGeom.h. + * + * + */ + +#ifndef _ST_GMT_CONSTS_H_ +#define _ST_GMT_CONSTS_H_ + +#include "StEvent/StEnumerations.h" // as many GMT consts are defined there +#include + +// The following are only used in StGmtGeom and StGmtDb, at present +const Int_t kGmtError = -999; +const Char_t kGmtErrorChar = -1; +const std::string kGmtErrorString = "XXXXXX"; + +// Constants related to clustering and pulses +const Int_t kGmtMaxClusterSize = 11; +const Int_t kGmtNumAdditionalStrips = 5; + +// For some reason, the MuDst directory fails during linking, so a seperate value is needed. +// For now, it is required that kMuGmtNumTimeBins == kGmtNumTimeBins +const Int_t kMuGmtNumTimeBins = 15; + + + +// Jan's definitions for the final 400-800 micron pitch design +// Note: +// using #define instead of const double to avoid requiring a .cpp +// file for the constants + +// #define kGmtRout 38.25 // cm , +// #define kGmtRlast 38.1571 // location of last R strip before Rout +// #define kGmtRmid 19.125 // cm, at Rout/2. +// #define kGmtRin 11.5 // cm, +// #define kGmtRfirst 11.5385 // location of first R strip after Rin +// #define kGmtPfirst 0.0324 // location of first Phi strip +// #define kGmtPlast 1.5384 // location of last Phi strip +// #define kGmtRflat 35.85 // cm, +// #define kGmtPhiflat (31.0/180.*3.1416) // rad +// #define kGmtRadPitch 0.09538 // nominal '800 mu pitch' +// #define kGmtPhiPitch 0.08 // 800 mu, at outer radi or at Rmid +// #define kGmtPhiAnglePitch 0.002095 +// #define kGmtDeadQuadEdge 1.2 // (cm) effective dead area along quadrant edges +// #define kGmtMaxClusterSize 11 //maximum cluster size in strips that a cluster algo should return +// #define kGmtNumAdditionalStrips 5 //strips in addition to the cluster size that are passed up. Mainly for debugging. + +// Constants for GMT (RW 3/15/2013) +#define kGmtXPitch 0.08 // nominal '800 mu pitch' +#define kGmtYPitch 0.08 // 800 mu, at outer radi or at Rmid +#define kGmtSfirst 0.00 // location of first strip (local X) relative to readout plane origin +#define kGmtPfirst 0.00 // location of first pad (local Y) relative to readout plane origin +#define kGmtSlast 10.16 // location of last strip [ (128 - 1 at origin) = 127 ]*pitch (in cm) +#define kGmtPlast 10.16 // location of last pad [ (128 - 1 at origin) = 127 ]*pitch (in cm) + +#endif diff --git a/StRoot/StGmtUtil/geometry/StGmtGeom.cxx b/StRoot/StGmtUtil/geometry/StGmtGeom.cxx new file mode 100644 index 00000000000..37246e64052 --- /dev/null +++ b/StRoot/StGmtUtil/geometry/StGmtGeom.cxx @@ -0,0 +1,1131 @@ +/* StGmtGeom.cxx + * + * GMT geometry class declaration. + * + * \authors K.S. Engle and Richard Witt (witt@usna.edu) + * based on StFgtGeom + */ + +#include +//#include +#include +#include +#include "StGmtGeom.h" +#include "StMessMgr.h" + +//ClassImp(StGmtGeom) + +double StGmtGeom::mPi = TMath::Pi(); +double StGmtGeom::mHalfPi = TMath::PiOver2(); + + +Int_t StGmtGeom::encodeGeoId //// WORKING +( +// Int_t module, Char_t layer, Int_t strip + Int_t rdo, Int_t arm, Int_t apv, Int_t channel +) +{ + Short_t module = getModuleIdFromElecCoord( rdo, arm, apv ); + // locally map apv number into [0,1] + if ( apv <= 3 ) + { + apv = apv % 2; + } + else + { + apv = (apv - 12) % 2; + } + + Char_t testS='S'; + Char_t testP='P'; + + if ( module < 0 || module >= kGmtNumModules ) + { + LOG_DEBUG << "Module " << module << " out of range in StGmtGeom::encodeGeoId." << endm; + return kGmtError; + } + else if ( apv > 1 || apv < 0 ) + { + LOG_DEBUG << "APV " << apv << " out of range in StGmtGeom::encodeGeoId." << endm; + return kGmtError; + } + else if ( channel < 0 || channel >= kGmtNumStrips ) + { + LOG_DEBUG << "Channel " << channel << " out of range in StGmtGeom::encodeGeoId." << endm; + return kGmtError; + } + + return + ( module * kGmtNumLayers + apv ) * kGmtNumStrips + channel; // from 0 to 2047 since channel is from 0-127; +} + +Int_t StGmtGeom::decodeGeoId //// WORKING +( +// Int_t geoId, Short_t & module, Char_t & layer, Short_t & strip + Int_t geoId, Short_t & module, Int_t & layer, Short_t & strip +) +{ + if ( geoId < 0 || geoId >= kGmtNumGeoIds ) + { + LOG_DEBUG << "GeoId " << geoId << " out of range in StGmtGeom::decodeGeoId." << endm; + module = kGmtError; +// layer = kGmtErrorChar; + layer = kGmtError; + strip = kGmtError; + + return kGmtError; + } + + strip = geoId % kGmtNumStrips; + geoId /= kGmtNumStrips; + +// layer = ( geoId % kGmtNumLayers ) ? 'P' : 'S'; + Int_t apv = ( geoId % kGmtNumLayers ) ? 1 : 0; + geoId /= kGmtNumLayers; + + StGmtGeomData stripdata = mStrips[ strip + apv*kGmtNumStrips ]; + layer = stripdata.isY; + + module = geoId; + + return 0; +} + +//////////// FIX ME!!!!!!!!!!!!!!!!!!!! +std::string StGmtGeom::encodeGeoName +( + Int_t module, Char_t layer, Int_t strip +) +{ + Char_t testS='S'; + Char_t testP='P'; + + if ( module < 0 || module >= kGmtNumModules ) + { + LOG_DEBUG << "Module " << module << " out of range in StGmtGeom::encodeGeoName." << endm; + return kGmtErrorString; + } + else if ( + layer != testS + && layer != testP + ) + { + LOG_DEBUG << "Layer " << layer << " out of range in StGmtGeom::encodeGeoName." << endm; + return kGmtErrorString; + } + else if ( strip < 0 || strip >= kGmtNumStrips ) + { + LOG_DEBUG << "Strip " << strip << " out of range in StGmtGeom::encodeGeoName." << endm; + return kGmtErrorString; + } + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! + std::stringstream buff; +// buff << disc+1 << (Char_t)(quadrant+'A') << layer; + buff << module+1 << layer; + + if ( strip < 10 ) + buff << "00"; + else if ( strip < 100 ) + buff << "0"; + + buff << strip; + return buff.str(); +} + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! +Int_t StGmtGeom::decodeGeoName +( +// const std::string & geoName, +// Short_t & module, Char_t & layer, Short_t & strip + const std::string & geoName, + Short_t & module, Int_t & layer, Short_t & strip +) +{ + Char_t testS='S'; + Char_t testP='P'; + + //assert( geoName.size() == 6 ); +// disc = geoName[0] - '1'; +// quadrant = geoName[1] - 'A'; +// layer = geoName[2]; +// strip = std::atoi( (geoName.substr(3)).c_str() ); + module = geoName[0] - '1'; + layer = geoName[2]; + strip = std::atoi( (geoName.substr(3)).c_str() ); + + // This is unlikely to catch all errors with the geoName, but it should + // do fairly well. + if ( + module < 0 + || module >= kGmtNumModules + || ( +// layer != testS +// && layer != testP + layer < 0 || layer > 1 + ) + || strip < 0 + || strip > kGmtNumStrips + ) + { + LOG_DEBUG << "Malformed geoName " << geoName << " in StGmtGeom::decodeGeoName." << endm; + module = kGmtError; + layer = kGmtErrorChar; + strip = kGmtError; + + return kGmtError; + } + + return 0; +} + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! +std::string StGmtGeom::translateGeoIdToGeoName( Int_t geoId ) +{ + Short_t module, strip; +// Char_t layer; + Int_t layer; + + if ( geoId < 0 || geoId >= kGmtNumGeoIds ) + { + LOG_DEBUG << "GeoId " << geoId << " out of range in StGmtGeom::translateGeoIdToGeoName." << endm; + return kGmtErrorString; + } + + decodeGeoId( geoId, module, layer, strip ); + return encodeGeoName( module, layer, strip ); +} + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! +Int_t StGmtGeom::translateGeoNameToGeoId( const std::string & geoName ) +{ + Short_t module, strip; +// Char_t layer; + Int_t layer; + Int_t rdo, arm, apv, channel; + + // Error message already taken care of in decodeGeoName. + if ( decodeGeoName( geoName, module, layer, strip ) < 0 ) + return kGmtError; + +// return encodeGeoId( module, layer, strip ); + return encodeGeoId( rdo, arm, apv, channel ); +} + +Int_t StGmtGeom::getPhysicalCoordinate //// WORKING +( + Int_t geoId, +// Short_t & module, Char_t & layer + Short_t & module, Int_t & layer +) +{ + if ( geoId < 0 || geoId >= kGmtNumGeoIds ) + { + LOG_DEBUG << "GeoId " << geoId << " out of range in StGmtGeom::getPhysicalCoordinate." << endm; + module = kGmtError; + layer = kGmtErrorChar; + + return kGmtError; + } + + Short_t strip; + + decodeGeoId( geoId, module, layer, strip ); + + return 0; +} + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! +Int_t StGmtGeom::getPhysicalCoordinate +( +// const std::string & geoName, +// Short_t & module, Char_t & layer + const std::string & geoName, + Short_t & module, Int_t & layer +) +{ + Short_t strip; + + if ( decodeGeoName( geoName, module, layer, strip ) < 0 ) + { + // Error is mostly handled by the decodeGeoName call. + module = kGmtError; +// layer = kGmtErrorChar; + layer = kGmtError; + + return kGmtError; + } + + return 0; +} + +Short_t StGmtGeom::getModuleIdFromElecCoord +( + Int_t rdo, Int_t arm, Int_t apv +) +{ + if ( (rdo - 1) < 0 || (rdo - 1) >= kGmtNumRdos ) + { + LOG_DEBUG << "RDO " << rdo << " out of range in StGmtGeom::getModuleIdFromElecCoord." << endm; + return kGmtError; + } + else if ( arm < 0 || arm >= kGmtNumArms ) + { + LOG_DEBUG << "ARM " << arm << " out of range in StGmtGeom::getModuleIdFromElecCoord." << endm; + return kGmtError; + } + else if ( apv < 0 || apv > kGmtMaxApvId || (apv > 3 && apv < 12) ) + { + LOG_DEBUG << "APV " << apv << " out of range in StGmtGeom::getModuleIdFromElecCoord." << endm; + return kGmtError; + } + + if ( arm == 0 ) + { + if( apv == 0 || apv == 1 ) + return 0; + else if( apv == 2 || apv == 3 ) + return 1; + else if( apv == 12 || apv == 13 ) + return 2; + else if( apv == 14 || apv == 15 ) + return 3; + else + { + LOG_DEBUG << "Invalid electronics coordinates in StGmtGeom::getModuleIdFromElecCoord." << endm; + return kGmtError; + } + } + else + { + if( apv == 0 || apv == 1 ) + return 4; + else if( apv == 2 || apv == 3 ) + return 5; + else if( apv == 12 || apv == 13 ) + return 6; + else if( apv == 14 || apv == 15 ) + return 7; + else + { + LOG_DEBUG << "Invalid electronics coordinates in StGmtGeom::getModuleIdFromElecCoord." << endm; + return kGmtError; + } + } +} + +Int_t StGmtGeom::getCoordNumFromElecCoord +( + Int_t rdo, Int_t arm, Int_t apv, Int_t channel +) +{ + if ( (rdo - 1) < 0 || (rdo - 1) >= kGmtNumRdos ) + { + LOG_DEBUG << "RDO " << rdo << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( arm < 0 || arm >= kGmtNumArms ) + { + LOG_DEBUG << "ARM " << arm << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( apv < 0 || apv > kGmtMaxApvId || (apv > 3 && apv < 12) ) + { + LOG_DEBUG << "APV " << apv << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( channel < 0 || channel >= kGmtNumChannels ) + { + LOG_DEBUG << "Channel " << channel << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + +// // locally map apv number into [0,1] +// if ( apv == 3 ) +// { +// apv = apv/2; +// } +// else +// { +// apv = (apv - 12)/2; +// } + + + StGmtGeomData stripdata = mStrips[ channel + (apv % 2) * kGmtNumStrips ]; + + return stripdata.coordinate; + +} + +Double_t StGmtGeom::getPositionFromElecCoord//here +( + Int_t rdo, Int_t arm, Int_t apv, Int_t channel + ) +{ + Int_t apvi = apv; + if ( (rdo - 1) < 0 || (rdo - 1) >= kGmtNumRdos ) + { + LOG_DEBUG << "RDO " << rdo << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( arm < 0 || arm >= kGmtNumArms ) + { + LOG_DEBUG << "ARM " << arm << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( apv < 0 || apv > kGmtMaxApvId || (apv > 3 && apv < 12) ) + { + LOG_DEBUG << "APV " << apv << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + else if ( channel < 0 || channel >= kGmtNumChannels ) + { + LOG_DEBUG << "Channel " << channel << " out of range in StGmtGeom::getCoordNumFromElecCoord." << endm; + return kGmtError; + } + + // locally map apv number into [0,1] +// if ( apv <= 3 ) +// { +// apv = apv/2; +// } +// else +// { +// apv = (apv - 12)/2; +// } + + apv = apv%2; + + StGmtGeomData stripdata = mStrips[ channel + apv*kGmtNumStrips ];//0-257 + //LOG_INFO << "rdo=" << rdo << "\tarm="<< arm << "\tapv="<< apvi <<"\tchannel = "<< channel << endm; + //LOG_INFO << "apv2= " << apv << "\tstr=" << kGmtNumStrips << "\t[]=" << channel + apv*kGmtNumStrips << "\t===> " <= kGmtNumGeoIds ) + { + LOG_DEBUG << "GeoId " << geoId << " out of range in StGmtGeom::getPhysicalCoordinate." << endm; + module = kGmtError; +// layer = kGmtErrorChar; + layer = kGmtError; + + return kGmtError; + } + + return computeGlobalPhysicalCoordinate( layer, strip); +} + + //////////// FIX ME!!!!!!!!!!!!!!!!!!!! +// The ordinate, lowerSpan and upperSpan are all in centimeters or +// radians, depending on the layer. +Int_t StGmtGeom::getGlobalPhysicalCoordinate +( + const std::string & geoName, +// Short_t & module, Char_t & layer + Short_t & module, Int_t & layer +) +{ + Short_t strip; + + if ( decodeGeoName( geoName, module, layer, strip ) < 0 ) + { + // Error is mostly handled by the decodeGeoName call. + module = kGmtError; +// layer = kGmtErrorChar; + layer = kGmtError; + + return kGmtError; + } + + return computeGlobalPhysicalCoordinate( layer, strip); + +} + + + + +// Please note that the following functions do NOT access the STAR +// database to find mapping information. They assume the most +// straight-forward mapping scheme and use that. +// For those functions that have them, currently rdo can only be 1, arm +// can be 0-1, apv can be 0-23 (although 4-11 are not +// technically valid) and channel is 0-127. +Int_t StGmtGeom::encodeElectronicId //// WORKING +( + Int_t rdo, Int_t arm, Int_t apv, Int_t channel +) +{ + if ( (rdo - 1) < 0 || (rdo - 1) >= kGmtNumRdos ) + { + LOG_DEBUG << "RDO " << rdo << " out of range in StGmtGeom::encodeElectronicId." << endm; + return kGmtError; + } + else if ( arm < 0 || arm >= kGmtNumArms ) + { + LOG_DEBUG << "ARM " << arm << " out of range in StGmtGeom::encodeElectronicId." << endm; + return kGmtError; + } + else if ( apv < 0 || apv > kGmtMaxApvId || (apv > 3 && apv < 12) ) + { + LOG_DEBUG << "APV " << apv << " out of range in StGmtGeom::encodeElectronicId." << endm; + return kGmtError; + } + else if ( channel < 0 || channel >= kGmtNumChannels ) + { + LOG_DEBUG << "Channel " << channel << " out of range in StGmtGeom::encodeElectronicId." << endm; + return kGmtError; + } + + return channel + kGmtNumStrips*(apv + 4*arm); +} + +Int_t StGmtGeom::decodeElectronicId //// WORKING +( + Int_t elecId, + Int_t &rdo, Int_t &arm, Int_t &apv, Int_t &channel +) +{ + if ( elecId < 0 || elecId >= kGmtNumElecIds ) + { + LOG_DEBUG << "Electronic ID " << elecId << " out of range in StGmtGeom::decodeElectronicId." << endm; + + rdo = kGmtError; + arm = kGmtError; + apv = kGmtError; + channel = kGmtError; + + return kGmtError; + } + + channel = elecId % 128; + elecId /= 128; + + apv = elecId % 4; + elecId /= 4; + + arm = elecId; + rdo = 1; + + return 0; +} + + + +// Whether the reverse map is valid +Bool_t StGmtGeom::mReverseNaiveMappingValid = 0; + +// The reverse map data member +Int_t StGmtGeom::mReverseNaiveMapping[ kGmtNumStripsPerModule ]; + + +// Initialize our physical coordinate database here. These are: +// APV,Chan,Strip(0) or Pad(1),coordinate #,Location (cm),Signal *,Readout Order +// The index corresponds to int(apv/2)+channel (assuming that the apv is in +// [0,12). If apv is in [12,24), then the index is int((apv-12)/2)+channel. +StGmtGeom::StGmtGeomData StGmtGeom::mStrips[] = +{ + {0,0,0,54,4.28,"S54",0}, + {0,1,0,37,2.92,"S37",1}, + {0,2,0,31,2.44,"S31",2}, + {0,3,0,8,0.6,"S8",3}, + {0,4,0,51,4.04,"S51",4}, + {0,5,0,35,2.76,"S35",5}, + {0,6,0,29,2.28,"S29",6}, + {0,7,0,21,1.64,"S21",7}, + {0,8,0,60,4.76,"S60",8}, + {0,9,0,33,2.6,"S33",9}, + {0,10,0,27,2.12,"S27",10}, + {0,11,0,12,0.92,"S12",11}, + {0,12,0,41,3.24,"S41",12}, + {0,13,0,55,4.36,"S55",13}, + {0,14,0,25,1.96,"S25",14}, + {0,15,0,14,1.08,"S14",15}, + {0,16,1,999,-999,"NC",16}, // Not connected + {0,17,1,39,3.08,"P39",17}, + {0,18,1,2,0.12,"P2",18}, + {0,19,1,24,1.88,"P24",19}, + {0,20,1,54,4.28,"P54",20}, + {0,21,1,37,2.92,"P37",21}, + {0,22,1,31,2.44,"P31",22}, + {0,23,1,23,1.8,"P23",23}, + {0,24,1,43,3.4,"P43",24}, + {0,25,1,35,2.76,"P35",25}, + {0,26,1,29,2.28,"P29",26}, + {0,27,1,21,1.64,"P21",27}, + {0,28,1,41,3.24,"P41",28}, + {0,29,0,46,3.64,"S46",29}, + {0,30,1,26,2.04,"P26",30}, + {0,31,1,7,0.52,"P7",31}, + {0,32,0,52,4.12,"S52",32}, + {0,33,1,57,4.52,"P57",33}, + {0,34,0,6,0.44,"S6",34}, + {0,35,0,7,0.52,"S7",35}, + {0,36,0,48,3.8,"S48",36}, + {0,37,0,50,3.96,"S50",37}, + {0,38,0,3,0.2,"S3",38}, + {0,39,0,20,1.56,"S20",39}, + {0,40,0,59,4.68,"S59",40}, + {0,41,1,999,-999,"NC",41}, // Not connected + {0,42,1,1,0.04,"P1",42}, + {0,43,0,11,0.84,"S11",43}, + {0,44,0,42,3.32,"S42",44}, + {0,45,0,56,4.44,"S56",45}, + {0,46,0,13,1,"S13",46}, + {0,47,0,16,1.24,"S16",47}, + {0,48,1,61,4.84,"P61",48}, + {0,49,1,46,3.64,"P46",49}, + {0,50,1,33,2.6,"P33",50}, + {0,51,1,25,1.96,"P25",51}, + {0,52,1,56,4.44,"P56",52}, + {0,53,1,50,3.96,"P50",53}, + {0,54,1,10,0.76,"P10",54}, + {0,55,1,15,1.16,"P15",55}, + {0,56,1,47,3.72,"P47",56}, + {0,57,1,55,4.36,"P55",57}, + {0,58,1,12,0.92,"P12",58}, + {0,59,1,9,0.68,"P9",59}, + {0,60,1,42,3.32,"P42",60}, + {0,61,0,62,4.92,"S62",61}, + {0,62,1,27,2.12,"P27",62}, + {0,63,1,6,0.44,"P6",63}, + {0,64,0,47,3.72,"S47",64}, + {0,65,0,49,3.88,"S49",65}, + {0,66,0,5,0.36,"S5",66}, + {0,67,0,10,0.76,"S10",67}, + {0,68,0,45,3.56,"S45",68}, + {0,69,1,59,4.68,"P59",69}, + {0,70,0,4,0.28,"S4",70}, + {0,71,0,9,0.68,"S9",71}, + {0,72,0,39,3.08,"S39",72}, + {0,73,0,61,4.84,"S61",73}, + {0,74,1,4,0.28,"P4",74}, + {0,75,0,18,1.4,"S18",75}, + {0,76,0,40,3.16,"S40",76}, + {0,77,0,57,4.52,"S57",77}, + {0,78,0,24,1.88,"S24",78}, + {0,79,0,2,0.12,"S2",79}, + {0,80,1,62,4.92,"P62",80}, + {0,81,1,48,3.8,"P48",81}, + {0,82,1,32,2.52,"P32",82}, + {0,83,1,22,1.72,"P22",83}, + {0,84,1,58,4.6,"P58",84}, + {0,85,1,52,4.12,"P52",85}, + {0,86,1,30,2.36,"P30",86}, + {0,87,1,17,1.32,"P17",87}, + {0,88,1,49,3.88,"P49",88}, + {0,89,0,44,3.48,"S44",89}, + {0,90,1,28,2.2,"P28",90}, + {0,91,1,11,0.84,"P11",91}, + {0,92,1,53,4.2,"P53",92}, + {0,93,0,63,5,"S63",93}, + {0,94,1,18,1.4,"P18",94}, + {0,95,1,5,0.36,"P5",95}, + {0,96,1,63,5,"P63",96}, + {0,97,0,36,2.84,"S36",97}, + {0,98,0,30,2.36,"S30",98}, + {0,99,0,22,1.72,"S22",99}, + {0,100,0,43,3.4,"S43",100}, + {0,101,0,34,2.68,"S34",101}, + {0,102,0,28,2.2,"S28",102}, + {0,103,0,19,1.48,"S19",103}, + {0,104,0,53,4.2,"S53",104}, + {0,105,0,58,4.6,"S58",105}, + {0,106,0,26,2.04,"S26",106}, + {0,107,0,17,1.32,"S17",107}, + {0,108,0,38,3,"S38",108}, + {0,109,0,32,2.52,"S32",109}, + {0,110,0,23,1.8,"S23",110}, + {0,111,0,1,0.04,"S1",111}, + {0,112,1,60,4.76,"P60",112}, + {0,113,1,40,3.16,"P40",113}, + {0,114,1,8,0.6,"P8",114}, + {0,115,1,14,1.08,"P14",115}, + {0,116,1,45,3.56,"P45",116}, + {0,117,1,38,3,"P38",117}, + {0,118,0,15,1.16,"S15",118}, + {0,119,1,19,1.48,"P19",119}, + {0,120,1,51,4.04,"P51",120}, + {0,121,1,36,2.84,"P36",121}, + {0,122,1,16,1.24,"P16",122}, + {0,123,1,13,1,"P13",123}, + {0,124,1,44,3.48,"P44",124}, + {0,125,1,34,2.68,"P34",125}, + {0,126,1,20,1.56,"P20",126}, + {0,127,1,3,0.2,"P3",127}, + {1,0,0,117,9.32,"S117",128}, + {1,1,0,100,7.96,"S100",129}, + {1,2,0,94,7.48,"S94",130}, + {1,3,0,71,5.64,"S71",131}, + {1,4,0,114,9.08,"S114",132}, + {1,5,0,98,7.8,"S98",133}, + {1,6,0,92,7.32,"S92",134}, + {1,7,0,84,6.68,"S84",135}, + {1,8,0,123,9.8,"S123",136}, + {1,9,0,96,7.64,"S96",137}, + {1,10,0,90,7.16,"S90",138}, + {1,11,0,75,5.96,"S75",139}, + {1,12,0,104,8.28,"S104",140}, + {1,13,0,118,9.4,"S118",141}, + {1,14,0,88,7,"S88",142}, + {1,15,0,77,6.12,"S77",143}, + {1,16,1,999,-999,"NC",144}, // Not connected + {1,17,1,102,8.12,"P102",145}, + {1,18,1,65,5.16,"P65",146}, + {1,19,1,87,6.92,"P87",147}, + {1,20,1,117,9.32,"P117",148}, + {1,21,1,100,7.96,"P100",149}, + {1,22,1,94,7.48,"P94",150}, + {1,23,1,86,6.84,"P86",151}, + {1,24,1,106,8.44,"P106",152}, + {1,25,1,98,7.8,"P98",153}, + {1,26,1,92,7.32,"P92",154}, + {1,27,1,84,6.68,"P84",155}, + {1,28,1,104,8.28,"P104",156}, + {1,29,0,109,8.68,"S109",157}, + {1,30,1,89,7.08,"P89",158}, + {1,31,1,70,5.56,"P70",159}, + {1,32,0,115,9.16,"S115",160}, + {1,33,1,120,9.56,"P120",161}, + {1,34,0,69,5.48,"S69",162}, + {1,35,0,70,5.56,"S70",163}, + {1,36,0,111,8.84,"S111",164}, + {1,37,0,113,9,"S113",165}, + {1,38,0,66,5.24,"S66",166}, + {1,39,0,83,6.6,"S83",167}, + {1,40,0,122,9.72,"S122",168}, + {1,41,1,999,-999,"NC",169}, // Not connected + {1,42,1,64,5.08,"P64",170}, + {1,43,0,74,5.88,"S74",171}, + {1,44,0,105,8.36,"S105",172}, + {1,45,0,119,9.48,"S119",173}, + {1,46,0,76,6.04,"S76",174}, + {1,47,0,79,6.28,"S79",175}, + {1,48,1,124,9.88,"P124",176}, + {1,49,1,109,8.68,"P109",177}, + {1,50,1,96,7.64,"P96",178}, + {1,51,1,88,7,"P88",179}, + {1,52,1,119,9.48,"P119",180}, + {1,53,1,113,9,"P113",181}, + {1,54,1,73,5.8,"P73",182}, + {1,55,1,78,6.2,"P78",183}, + {1,56,1,110,8.76,"P110",184}, + {1,57,1,118,9.4,"P118",185}, + {1,58,1,75,5.96,"P75",186}, + {1,59,1,72,5.72,"P72",187}, + {1,60,1,105,8.36,"P105",188}, + {1,61,0,125,9.96,"S125",189}, + {1,62,1,90,7.16,"P90",190}, + {1,63,1,69,5.48,"P69",191}, + {1,64,0,110,8.76,"S110",192}, + {1,65,0,112,8.92,"S112",193}, + {1,66,0,68,5.4,"S68",194}, + {1,67,0,73,5.8,"S73",195}, + {1,68,0,108,8.6,"S108",196}, + {1,69,1,122,9.72,"P122",197}, + {1,70,0,67,5.32,"S67",198}, + {1,71,0,72,5.72,"S72",199}, + {1,72,0,102,8.12,"S102",200}, + {1,73,0,124,9.88,"S124",201}, + {1,74,1,67,5.32,"P67",202}, + {1,75,0,81,6.44,"S81",203}, + {1,76,0,103,8.2,"S103",204}, + {1,77,0,120,9.56,"S120",205}, + {1,78,0,87,6.92,"S87",206}, + {1,79,0,65,5.16,"S65",207}, + {1,80,1,125,9.96,"P125",208}, + {1,81,1,111,8.84,"P111",209}, + {1,82,1,95,7.56,"P95",210}, + {1,83,1,85,6.76,"P85",211}, + {1,84,1,121,9.64,"P121",212}, + {1,85,1,115,9.16,"P115",213}, + {1,86,1,93,7.4,"P93",214}, + {1,87,1,80,6.36,"P80",215}, + {1,88,1,112,8.92,"P112",216}, + {1,89,0,107,8.52,"S107",217}, + {1,90,1,91,7.24,"P91",218}, + {1,91,1,74,5.88,"P74",219}, + {1,92,1,116,9.24,"P116",220}, + {1,93,0,126,10.04,"S126",221}, + {1,94,1,81,6.44,"P81",222}, + {1,95,1,68,5.4,"P68",223}, + {1,96,1,999,-999,"NC",224}, // Not connected + {1,97,0,99,7.88,"S99",225}, + {1,98,0,93,7.4,"S93",226}, + {1,99,0,85,6.76,"S85",227}, + {1,100,0,106,8.44,"S106",228}, + {1,101,0,97,7.72,"S97",229}, + {1,102,0,91,7.24,"S91",230}, + {1,103,0,82,6.52,"S82",231}, + {1,104,0,116,9.24,"S116",232}, + {1,105,0,121,9.64,"S121",233}, + {1,106,0,89,7.08,"S89",234}, + {1,107,0,80,6.36,"S80",235}, + {1,108,0,101,8.04,"S101",236}, + {1,109,0,95,7.56,"S95",237}, + {1,110,0,86,6.84,"S86",238}, + {1,111,0,64,5.08,"S64",239}, + {1,112,1,123,9.8,"P123",240}, + {1,113,1,103,8.2,"P103",241}, + {1,114,1,71,5.64,"P71",242}, + {1,115,1,77,6.12,"P77",243}, + {1,116,1,108,8.6,"P108",244}, + {1,117,1,101,8.04,"P101",245}, + {1,118,0,78,6.2,"S78",246}, + {1,119,1,82,6.52,"P82",247}, + {1,120,1,114,9.08,"P114",248}, + {1,121,1,99,7.88,"P99",249}, + {1,122,1,79,6.28,"P79",250}, + {1,123,1,76,6.04,"P76",251}, + {1,124,1,107,8.52,"P107",252}, + {1,125,1,97,7.72,"P97",253}, + {1,126,1,83,6.6,"P83",254}, + {1,127,1,66,5.24,"P66",255} +}; + +// // This initialized an idealized mapping. The index of the array is the coordinate number (X coordinates are first then Y). +// // The value in the array is the readout number of the coordinate (i.e. the index into mStrips[] above). +Int_t StGmtGeom::mNaiveMapping[] = +{ + 111, + 79, + 38, + 70, + 66, + 34, + 35, + 3, + 71, + 67, + 43, + 11, + 46, + 15, + 118, + 47, + 107, + 75, + 103, + 39, + 7, + 99, + 110, + 78, + 14, + 106, + 10, + 102, + 6, + 98, + 2, + 109, + 9, + 101, + 5, + 97, + 1, + 108, + 72, + 76, + 12, + 44, + 100, + 89, + 68, + 29, + 64, + 36, + 65, + 37, + 4, + 32, + 104, + 0, + 13, + 45, + 77, + 105, + 40, + 8, + 73, + 61, + 93, + 239, + 207, + 166, + 198, + 194, + 162, + 163, + 131, + 199, + 195, + 171, + 139, + 174, + 143, + 246, + 175, + 235, + 203, + 231, + 167, + 135, + 227, + 238, + 206, + 142, + 234, + 138, + 230, + 134, + 226, + 130, + 237, + 137, + 229, + 133, + 225, + 129, + 236, + 200, + 204, + 140, + 172, + 228, + 217, + 196, + 157, + 192, + 164, + 193, + 165, + 132, + 160, + 232, + 128, + 141, + 173, + 205, + 233, + 168, + 136, + 201, + 189, + 221, + 42, + 18, + 127, + 74, + 95, + 63, + 31, + 114, + 59, + 54, + 91, + 58, + 123, + 115, + 55, + 122, + 87, + 94, + 119, + 126, + 27, + 83, + 23, + 19, + 51, + 30, + 62, + 90, + 26, + 86, + 22, + 82, + 50, + 125, + 25, + 121, + 21, + 117, + 17, + 113, + 28, + 60, + 24, + 124, + 116, + 49, + 56, + 81, + 88, + 53, + 120, + 85, + 92, + 20, + 57, + 52, + 33, + 84, + 69, + 112, + 48, + 80, + 96, + 170, + 146, + 255, + 202, + 223, + 191, + 159, + 242, + 187, + 182, + 219, + 186, + 251, + 243, + 183, + 250, + 215, + 222, + 247, + 254, + 155, + 211, + 151, + 147, + 179, + 158, + 190, + 218, + 154, + 214, + 150, + 210, + 178, + 253, + 153, + 249, + 149, + 245, + 145, + 241, + 156, + 188, + 152, + 252, + 244, + 177, + 184, + 209, + 216, + 181, + 248, + 213, + 220, + 148, + 185, + 180, + 161, + 212, + 197, + 240, + 176, + 208, + 16, + 41, + 144, + 169, + 224 +}; + +// Module locations at the corner of the GEM, +// per email from W.J. Llope to stargmt-l@lists.bnl.gov +// on 2012/10/31 + +Double_t StGmtGeom::getModuleZ(int iModule) { + switch (iModule) { + case 4: + case 0 : return 77.768 * 2.54; // inches => cm + case 5: + case 1 : return 2.729 * 2.54; // inches => cm + case 7: + case 3 : return -2.729 * 2.54; // inches => cm + case 6: + case 2 : return -77.768 * 2.54; // inches => cm + default : return -999; + } +} + +Double_t StGmtGeom::getModulePhi(int iModule) { + double R = 85.606 * 2.54; // inches => cm + double deltaphi = 5./R; // crude radian conversion + switch (iModule) { + case 0: + case 1 : return TMath::Pi()*(10./6.)-deltaphi; + case 2: + case 3 : return TMath::Pi()*(10./6.)+deltaphi; + case 4: + case 5 : return TMath::Pi()*(1./6.)-deltaphi; + case 6: + case 7 : return TMath::Pi()*(1./6.)+deltaphi; + default : return 0; + } +} + diff --git a/StRoot/StGmtUtil/geometry/StGmtGeom.h b/StRoot/StGmtUtil/geometry/StGmtGeom.h new file mode 100644 index 00000000000..dbfc4479f17 --- /dev/null +++ b/StRoot/StGmtUtil/geometry/StGmtGeom.h @@ -0,0 +1,270 @@ +/* StGmtGeom.h + * + * GMT geometry class declaration. + * + * \authors K.S. Engle and Richard Witt (witt@usna.edu) + * based on StFgtGeom + */ + +#ifndef _ST_GMT_GEOM_H_ +#define _ST_GMT_GEOM_H_ + +//#include +#include +#include +#include +#include +#include +#include +#include +#include "StRoot/StGmtUtil/StGmtConsts.h" + + +// StGmtGeom is a "singleton" class. Only one of it needs to exist in any +// program. However, because the data contained in this class is entirely +// static, the class itself is also entirely static. +class StGmtGeom { + public: + + // For all functions where they appear: Disc can be >= 0 (in theory, + // although only values 0-5 work at the moment, I believe). Quadrant + // is 0-3. Layer is 'P' or 'R'. Strip is 0-720 + + // Location of modules in Z. + static Double_t getModuleZ(int iModule); + + // Location of modules in phi. + static Double_t getModulePhi(int iModule); + + static Int_t getNaiveMapping( Int_t idx ) { return mNaiveMapping[ idx ]; } + + // geoId is a unique number used to identify a specific strip + // on a specific disk/quadrant/layer/strip. Please NOTE: + // The set of geoIds IS NOT CONTINUOUS simply becuase strip + // number is not continuous. On the R plane strips 280-399 + // are not implemented. + // static Int_t encodeGeoId( Int_t module, Char_t layer, Int_t strip ); + static Int_t encodeGeoId( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + // static Int_t decodeGeoId( Int_t geoId, Short_t & module, Char_t & layer, Short_t & strip ); + static Int_t decodeGeoId( Int_t geoId, Short_t & module, Int_t & layer, Short_t & strip ); + + + // Geoname is human readable form of geoId + static std::string encodeGeoName( Int_t module, Char_t layer, Int_t strip ); + // static Int_t decodeGeoName( const std::string & geoName, Short_t & module, Char_t & layer, Short_t & strip ); + static Int_t decodeGeoName( const std::string & geoName, Short_t & module, Int_t & layer, Short_t & strip ); + static std::string translateGeoIdToGeoName( Int_t geoId ); + static Int_t translateGeoNameToGeoId( const std::string & geoName ); + + // Returns range upper and lower range of R or Phi valus depending on geoId. + // NOTE phi values are only local - that is they are the same for each quadrant + // The ordinate, lowerSpan and upperSpan are all in centimeters or radians + // static Int_t getPhysicalCoordinate( Int_t geoId, Short_t & module, Char_t & layer); + static Int_t getPhysicalCoordinate( Int_t geoId, Short_t & module, Int_t & layer); + + + + // Returns range upper and lower range of R or Phi valus depending on geoName. + // NOTE phi values are only local - that is they are the same for each quadrant + // The ordinate, lowerSpan and upperSpan are all in centimeters or radians + // static Int_t getPhysicalCoordinate( const std::string & geoName, Short_t & module, Char_t & layer); + static Int_t getPhysicalCoordinate( const std::string & geoName, Short_t & module, Int_t & layer); + + + // Similar to getPhysicalCoordinate but returns phi in STAR coordinate system + // static Int_t getGlobalPhysicalCoordinate( Int_t geoId, Short_t & module, Char_t & layer); + static Int_t getGlobalPhysicalCoordinate( Int_t geoId, Short_t & module, + Int_t & layer); + + //Similar to getPhysicalCoordinate but returns phi in STAR coordinate system + // static Int_t getGlobalPhysicalCoordinate( const std::string & geoName, + // Short_t & module, + // Char_t & layer); + static Int_t getGlobalPhysicalCoordinate( const std::string & geoName, Short_t & module, Int_t & layer); + + static Int_t getCoordNumFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t channel); + static Double_t getPositionFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t channel); + + // Please note that the following functions do NOT access the STAR + // database to find mapping information. They assume the most + // straight-forward mapping scheme and use that. + // For those functions that have them, currently rdo can be 1-2, arm + // can be 0-5, apv can be 0-23 (although 10, 11, 22, and 23 are not + // technically valid) and channel is 0-127. To access database + // functions please use functions in StGmtDb.h + + // Electronic Id is determined from the electronic devices + // rdo/arm/apv/channel and does form a continuous set of integers. + // The mapping from geoId to electronicId is accessible via the + // database or from "naive" functions below + + static Int_t encodeElectronicId( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + + static Int_t decodeElectronicId( Int_t elecId, Int_t &rdo, Int_t &arm, Int_t &apv, Int_t &channel ); + static Int_t getElectIdFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t ch ) { return encodeElectronicId(rdo,arm,apv,ch); } + static Int_t getElecCoordFromElectId( Int_t eID, Int_t& rdo, Int_t& arm, Int_t& apv, Int_t& ch ) { return decodeElectronicId(eID,rdo,arm,apv,ch); } + + static Short_t getModuleIdFromElecCoord(Int_t rdo, Int_t arm, Int_t apv); + + // get the octant for a given layer and strip + static Char_t getOctant( Char_t layer, Int_t strip ); + + // get the octant given the APV number + static Char_t getOctant( Int_t apv ); + + // get the octant given a phi in radians + // maps to i8: 0=A.L, 1=A.S, 2=B.L, .... 7=D.S + static Int_t getOctant( Double_t phi); + + static Int_t getNaiveGeoIdFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + static Int_t getNaiveElecCoordFromGeoId( Int_t geoId, Int_t& rdo, Int_t& arm, Int_t& apv, Int_t& channel ); + static std::string getNaiveGeoNameFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + static Int_t getNaivePhysCoordFromElecCoord( Int_t rdo, Int_t arm, Int_t apv, Int_t channel, Short_t & module, Char_t & layer); + + // This is similar to the above functions, but it takes electronic + // coordinates and only returns the final ordinate. This is here + // primarily so that it can be used as a drop in replacement for + // older code that has similar functionality. + static Double_t getNaiveMapping( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + static bool isNaiveR( Int_t rdo, Int_t arm, Int_t apv, Int_t channel ); + + // Jan's necessary functions start here. These were written by Jan, + // modified slightly by me. + // Jan: I have adjusted the dimensions to match GMT as build, September, 2011 + // static double rIn() { return kGmtRin; } + // static double rMid() { return kGmtRmid; } + // static double rOut() { return kGmtRout; } + static double sFirst() { return kGmtSfirst; } + static double sLast() { return kGmtSlast;} + static double pFirst() { return kGmtPfirst; } + static double pLast() { return kGmtPlast;} + + // static double radStrip_pitch() { return kGmtRadPitch; } // cm + // static double phiStrip_pitch() { return kGmtPhiAnglePitch; } // rad + static double xStrip_pitch() { return kGmtXPitch; } // cm + static double yStrip_pitch() { return kGmtYPitch; } // cm + + // static double yLimit() { return kGmtRout; } + + // deadQuadEdge is in cm, local ref frame + // static double deadQuadEdge() { return kGmtDeadQuadEdge; } + + // static double radStripOff() { return mRadStripOff; } + // static double phiStripOff() { return mPhiStripOff; } + // static double StripOff() { return mStripOff; } + // static double PadOff() { return mPadOff; } + + // static double phiQuadXaxis(int iquad); + // static bool inDisc( TVector3 rLab ); + // static bool belowFlat( TVector3 rLoc ); + // static int getQuad( double phiLab ); + static bool inModule( TVector3 rLab ); + + // What follows are some functions to help with the + // localXYtoStripID function. These are also written by Jan, modified + // slightly by me. + + // These next two return -1 on error. + static int rad2LocalStripId( double rad, double phi, double *binFrac=0 ); + static int phi2LocalStripId( double rad, double phi, double *binFrac=0 ); + + // static double rStrip_Phi_High(int rindex);//return upper phi range for an r strip + // static double rStrip_Phi_Low(int rindex);//return lower phi range for an r strip + // static double pHistrip_R_Low(int pindex);//return lower r range for a phi strip + // static double pHistrip_R_High(int pindex);//return upper r range for a phi strip + + protected: + + // StGmtGeomData stores data on each ordinate associated with each + // global ID used to index individual geometry elements. + // struct StGmtGeomData + // { + // // Bool_t isPhi; + // // Double_t ordinate; + // // Double_t lowerSpan; + // // Double_t upperSpan; + // Bool_t isX; + // Double_t localX; + // Double_t localY; + // }; + + struct StGmtGeomData + { + Int_t apv; // apv number (arm channel 0-3 or 12-15) + Int_t channel; // apv channel (0-127) + Bool_t isY; // is it a pad (true) or strip (false) + Int_t coordinate; // coordinate number + Double_t location; // in cm + TString signal; // label string + Int_t readout; // number in readout order + }; + + friend class StGmtDbFileMaker; + + + // Various constants used in Jan's conversion functions. + static double mPi; + static double mHalfPi; + // static double mRadStripOff; + // static double mPhiStripOff; + // static int mRadStripLOCId_number; + // static int mPhiStripLOCId_number; + // static double mStripOff; + // static double mPadOff; + + + // ---Private member variables--- + // static StGmtGeomData mStrips[ 2*kGmtNumStrips ]; + static StGmtGeomData mStrips[ kGmtNumStripsPerModule ]; // 128 X strip and 128 Y strips + + // maps from (apv*128 + channel) to ((layer=='P')*kGmtNumStrips + stripID) + // static Int_t mNaiveMapping[ kGmtNumChannels*kGmtApvsPerQuad ]; + static Int_t mNaiveMapping[ kGmtNumStripsPerModule ]; + + // reverse mapping: ((layer=='P')*kGmtNumStrips + stripID) to (apv*128 + channel) + static Bool_t mReverseNaiveMappingValid; + static Int_t mReverseNaiveMapping[ kGmtNumStripsPerModule ]; + static void makeReverseNaiveMappingValid(); + + private: + + // Calculates coordinates of strip in global coordinate system + // Units are in cm andradians, depending on the layer. + // static Int_t computeGlobalPhysicalCoordinate(Char_t & layer, + static Int_t computeGlobalPhysicalCoordinate(Int_t & layer, Short_t & strip); +}; + +//________________ +// get the octant given the phi in radians +inline Int_t StGmtGeom::getOctant( Double_t phi ) { + double phiDeg= 75 - ((phi*180)/mPi); + while ( phiDeg < 0 ) phiDeg+=360; + while ( phiDeg > 360 ) phiDeg-=360; + int i8=phiDeg/45; + return i8; +} + +// inline std::string StGmtGeom::getNaiveGeoNameFromElecCoord +// ( +// Int_t rdo, Int_t arm, Int_t apv, Int_t channel +// ) +// { +// return +// translateGeoIdToGeoName( +// getNaiveGeoIdFromElecCoord( rdo, arm, apv, channel ) +// ); +// } +// +// inline Int_t StGmtGeom::getNaivePhysCoordFromElecCoord +// ( +// Int_t rdo, Int_t arm, Int_t apv, Int_t channel, +// Short_t & module, Char_t & layer +// ) +// { +// return getPhysicalCoordinate( +// getNaiveGeoIdFromElecCoord( rdo, arm, apv, channel ), +// module, layer); +// } + +#endif + diff --git a/StarDb/Geometry/gmt/GmtOnModule.20140101.000005.C b/StarDb/Geometry/gmt/GmtOnModule.20140101.000005.C new file mode 100644 index 00000000000..a5c3f866a92 --- /dev/null +++ b/StarDb/Geometry/gmt/GmtOnModule.20140101.000005.C @@ -0,0 +1,21 @@ +///////////////////////////////////////////////////////// +// This file was generated by MakeGmtOnModule.C macros // +///////////////////////////////////////////////////////// + +TDataSet *CreateTable() { + if (!gROOT->GetClass("St_Survey")) return 0; + Survey_st row[8] = { + {0,1.0000,-0.0089,-0.0041,0.0089,1.0000,0.0038,0.0041,-0.0039,1.0000,-0.1329,-5.0360,-0.0710,0.0004,0.0010,0.0009,0.0051,0.0011,0.0036,""}, + {1,0.9999,0.0067,-0.0154,-0.0067,1.0000,0.0022,0.0154,-0.0021,0.9999,0.0286,-5.0587,-0.1892,0.0003,0.0013,0.0014,0.0037,0.0021,0.0018,""}, + {2,1.0000,-0.0032,0.0003,-0.0032,-1.0000,0.0005,0.0003,-0.0004,-1.0000,-0.1442,4.8446,-0.2043,0.0003,0.0009,0.0009,0.0043,0.0010,0.0026,""}, + {3,1.0000,0.0047,-0.0098,0.0047,-1.0000,-0.0016,-0.0099,0.0015,-1.0000,-0.0037,4.8077,-0.2756,0.0003,0.0013,0.0012,0.0034,0.0012,0.0013,""}, + {4,1.0000,0.0051,-0.0030,-0.0051,1.0000,0.0016,0.0030,-0.0015,1.0000,-0.2385,-5.2747,-0.4397,0.0004,0.0011,0.0011,0.0058,0.0013,0.0036,""}, + {5,0.9998,0.0107,-0.0129,-0.0107,0.9999,-0.0001,0.0129,0.0003,0.9999,-0.2755,-5.3099,-0.5536,0.0003,0.0014,0.0014,0.0045,0.0011,0.0017,""}, + {6,1.0000,0.0000,0.0000,0.0000,-1.0000,0.0000,0.0000,0.0000,-1.0000,0.0000,0.0000,0.0000,-0.0100,-0.0100,-0.0100,-0.0010,-0.0010,-0.0010,""}, + {7,0.9999,-0.0122,-0.0127,-0.0121,-0.9999,-0.0007,-0.0126,0.0008,-0.9999,-0.2584,4.6969,-0.8955,0.0003,0.0015,0.0015,0.0039,0.0010,0.0019,""} + }; + Int_t noModules = 8; + St_Survey *tableSet = new St_Survey("GmtOnModule",noModules); + for (Int_t i = 0; i < noModules; i++) tableSet->AddAt(&row[i].Id, i); + return (TDataSet *)tableSet; +} diff --git a/StarDb/Geometry/gmt/GmtOnModule.C b/StarDb/Geometry/gmt/GmtOnModule.C new file mode 100644 index 00000000000..b3a10b8167d --- /dev/null +++ b/StarDb/Geometry/gmt/GmtOnModule.C @@ -0,0 +1,31 @@ +TDataSet *CreateTable() { + // strip coordinate (x,y,0) => module coordinate (y,0,x) + if (!gROOT->GetClass("St_Survey")) return 0; + Survey_st rows[8] = { + //m u v w + {0, 1,0,0, 0, 1,0, 0,0, 1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {1, 1,0,0, 0, 1,0, 0,0, 1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {2, 1,0,0, 0,-1,0, 0,0,-1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {3, 1,0,0, 0,-1,0, 0,0,-1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {4, 1,0,0, 0, 1,0, 0,0, 1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {5, 1,0,0, 0, 1,0, 0,0, 1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {6, 1,0,0, 0,-1,0, 0,0,-1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}, + {7, 1,0,0, 0,-1,0, 0,0,-1, 0,0,0, 0,0,0,0,0,0,"2014AuAu15GMTGlobI"}}; + Int_t noModules = 8;// + St_Survey *tableSet = new St_Survey("GmtOnModule",noModules); + TString Rot; + Double_t dU[8] = { 5.0020, 5.0513, 4.8432, 4.8517, 5.2337, 5.2333, 0.0000, 4.7607}; + // T->Draw("uD-uP:barrel>>dU(8,-0.5,7.5,100,-2.5,2.5)","(dvAdcD/duAdcD)>0.5&&(dvAdcD/duAdcD)<2.&&abs(pT)>0.5&&sLength<40&&abs(vD-vP)<2.5&&abs(uD-uP)<2.5","colz") + Double_t dU1[8]= { 0.0027, 0.0009, 0.0125, 0.0004, 0.0158, 0.0259, 0.0000, -0.0032}; + Double_t dU2[8]= { 0.0023, -0.0031, -0.0119, -0.0088, -0.0054, -0.0072, 0.0000, -0.0067};// J pT > 0.2 + + Double_t dV[8] = {-0.0589, 0.1424, -0.3106, -0.2046, 0.3160, 0.5626, 0.0000, -0.8306}; + Double_t dV1[8]= {-0.0243, 0.0088, -0.0145, -0.0264, -0.0340, 0.0135, 0.0000, -0.0025}; + Double_t dV2[8]= { 0.0170, -0.0119, 0.0009, 0.0075, 0.0068, -0.0110, 0.0000, -0.0012};// J pT > 0.2 + for (Int_t m = 0; m < noModules; m++) { + rows[m].t1 -= (dU[m]+dU1[m])*rows[m].r11; + rows[m].t2 -= (dV[m]+dV1[m])*rows[m].r22; + tableSet->AddAt(&rows[m].Id); + } + return (TDataSet *) tableSet; +} diff --git a/StarDb/Geometry/gmt/GmtOnTpc.C b/StarDb/Geometry/gmt/GmtOnTpc.C new file mode 100644 index 00000000000..3aad996565e --- /dev/null +++ b/StarDb/Geometry/gmt/GmtOnTpc.C @@ -0,0 +1,40 @@ +TDataSet *CreateTable() { + if (!gROOT->GetClass("St_Survey") || !gROOT->GetClass("TGeoRotation")) return 0; + Survey_st row = {0, 1,0,0, 0,1,0, 0,0,1, 0,0,0, 1e-5,1e-5,1e-4,4e-3,4e-3,4e-3,"Pass 5"}; + // + /* A.Lebedev 07/28/15 + GMT modules installed + West mid of sectors 2 and 5 + East mid of sectors 17 and 22 + So it is a continuous 2 lines Phi = 210 and 300. + I don't know exact Z position + Best regards, + Alexie + */ + double inch = 2.54; + double R = 85.606; // inches => cm + double deltaphi = 5./R; // crude radian conversion + // sector = 5/19 sector = 2/22 + // sector 5 5 19 19 2 2 22 22 + // module 0 1 2 3 4 5 6 7 + Double_t z[8] = {77.768, 2.729, -77.768, -2.729, 77.768, 2.729, -77.768, -2.729}; + Int_t phi[8] = { 300, 300, 300, 300, 30, 30, 30, 30}; + + Int_t noModules = 8;// + St_Survey *tableSet = new St_Survey("GmtOnTpc",noModules); + TString Rot; + TGeoRotation *rotm = 0; + for (Int_t m = 0; m < noModules; m++) { + row.Id = m; + TGeoRotation *rotm = new TGeoRotation(Rot); + // rotm->RotateZ(-phi[m]); // Pass 6 + rotm->RotateZ(phi[m]); // Pass 5, 7 + Double_t *rotaion = rotm->GetRotationMatrix(); + memcpy(&row.r00, rotm->GetRotationMatrix(), 9*sizeof(Double_t)); + Double_t xyz[3] = {R*inch, 0, z[m]*inch}; + rotm->LocalToMaster(xyz, &row.t0); + tableSet->AddAt(&row.Id); + delete rotm; + } + return (TDataSet *) tableSet; +} From 7b302ddd25af38d260b2ecdcf078d1c2b5f3e167 Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 12:34:53 -0500 Subject: [PATCH 02/11] Addressing stylish comments --- StRoot/StAnalysisMaker/StAnalysisMaker.cxx | 4 ++-- StRoot/StEvent/StGmtHit.h | 8 ++++---- StRoot/StEvent/StGmtPoint.cxx | 8 -------- StRoot/StEvent/StGmtStrip.cxx | 2 -- StRoot/StEvent/StGmtStrip.h | 3 --- StRoot/StEvent/StGmtStripCollection.h | 16 +--------------- .../EventT.cxx | 0 .../EventT.h | 13 +++++++++---- .../{StGmtAligner => StGmtAligmentMaker}/HitT.h | 0 .../StGmtAligner.cxx | 0 .../StGmtAligner.h | 0 .../TrackT.h | 0 .../runGmtTree.C | 0 StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx | 2 +- 14 files changed, 17 insertions(+), 39 deletions(-) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/EventT.cxx (100%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/EventT.h (94%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/HitT.h (100%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/StGmtAligner.cxx (100%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/StGmtAligner.h (100%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/TrackT.h (100%) rename StRoot/{StGmtAligner => StGmtAligmentMaker}/runGmtTree.C (100%) diff --git a/StRoot/StAnalysisMaker/StAnalysisMaker.cxx b/StRoot/StAnalysisMaker/StAnalysisMaker.cxx index 5c5ccb1e380..fcc4a460952 100644 --- a/StRoot/StAnalysisMaker/StAnalysisMaker.cxx +++ b/StRoot/StAnalysisMaker/StAnalysisMaker.cxx @@ -1141,14 +1141,14 @@ void StAnalysisMaker::summarizeEvent(StEvent *event, Int_t mEventCounter) { if (event->numberOfPsds()) { LOG_QA << "# PSDs: " << event->numberOfPsds() << endm; } -//#ifdef _ST_GMT_HIT_H_ + if (event->gmtCollection() && event->gmtCollection()->getNumHits()) { LOG_QA << "# GMT hits: " << event->gmtCollection()->getNumHits() << ", points: " << event->gmtCollection()->getNumPoints() << ", strips: " << event->gmtCollection()->getNumStrips() << endm; } -//#endif /* _ST_GMT_HIT_H_ */ + if (event->rpsCollection()) { LOG_QA << "# RPS hits: " << event->rpsCollection()->clusters().size() << endm; } diff --git a/StRoot/StEvent/StGmtHit.h b/StRoot/StEvent/StGmtHit.h index 4296f6b4492..3049f8f06b5 100644 --- a/StRoot/StEvent/StGmtHit.h +++ b/StRoot/StEvent/StGmtHit.h @@ -55,13 +55,13 @@ class StGmtHit : public StHit { /// ADC error in Y Float_t getErrorAdcY() const { return mdAdcY; } /// Local X coordinate - Float_t getLocalX() const { return position().x(); } + Float_t getLocalX() const { return mPosition().x(); } /// Local X coordinate error - Float_t getErrorLocalX() const { return positionError().x(); } + Float_t getErrorLocalX() const { return mPositionError().x(); } /// Local Y coordinate - Float_t getLocalY() const { return position().y(); } + Float_t getLocalY() const { return mPosition().y(); } /// Local Y coordinate error - Float_t getErrorLocalY() const { return positionError().y(); } + Float_t getErrorLocalY() const { return mPositionError().y(); } /// Position in X Float_t getSigmaX() const { return mSigmaX; } /// Position error in X diff --git a/StRoot/StEvent/StGmtPoint.cxx b/StRoot/StEvent/StGmtPoint.cxx index 761f93be6f8..460032dce68 100644 --- a/StRoot/StEvent/StGmtPoint.cxx +++ b/StRoot/StEvent/StGmtPoint.cxx @@ -32,14 +32,6 @@ StGmtPoint::StGmtPoint( StGmtHit* hit1, StGmtHit* hit2, int key ) : StHit(), mKe } mHitLocalX = mHitLocalY = 0; - // if ( !(hit1->isY()) && hit2->isY() ) { - // mHitLocalX = hit1; - // mHitLocalY = hit2; - // } - // else if ( !(hit2->isY()) && hit1->isY() ) { - // mHitLocalX = hit2; - // mHitLocalY = hit1; - // } if ( !mHitLocalX || !mHitLocalY ) { LOG_ERROR << "Constructor not provided a (LocalX,LocalY) pair." << endm; diff --git a/StRoot/StEvent/StGmtStrip.cxx b/StRoot/StEvent/StGmtStrip.cxx index af77edd3771..e280b54356a 100644 --- a/StRoot/StEvent/StGmtStrip.cxx +++ b/StRoot/StEvent/StGmtStrip.cxx @@ -50,7 +50,6 @@ StGmtStrip::StGmtStrip(const StGmtStrip& h) : StObject(), mIsY( h.mIsY ), mIsC( h.mIsC ), mPosition( h.mPosition ), mMaxAdc(h.mMaxAdc), mMaxPedSubtractedAdc(h.mMaxPedSubtractedAdc), mMaxAdcTB(h.mMaxAdcTB), mMaxPedSubtractedAdcTB(h.mMaxPedSubtractedAdcTB), - // mClusterSeedType(h.mClusterSeedType), mCharge( h.mCharge ), mChargeUncert(h.mChargeUncert), mRdo( h.mRdo ), mArm( h.mArm ), mApv( h.mApv ), mChan( h.mChan ), mPed(h.mPed), mPedStdDev(h.mPedStdDev), mPedErr(h.mPedErr) { @@ -74,7 +73,6 @@ StGmtStrip& StGmtStrip::operator=( const StGmtStrip& h) { mMaxPedSubtractedAdc = h.mMaxPedSubtractedAdc; mMaxAdcTB = h.mMaxAdcTB; mMaxPedSubtractedAdcTB = h.mMaxPedSubtractedAdcTB; -// mClusterSeedType = h.mClusterSeedType; mCharge = h.mCharge; mChargeUncert = h.mChargeUncert; mRdo = h.mRdo; diff --git a/StRoot/StEvent/StGmtStrip.h b/StRoot/StEvent/StGmtStrip.h index c2012faa525..ed13f01a486 100644 --- a/StRoot/StEvent/StGmtStrip.h +++ b/StRoot/StEvent/StGmtStrip.h @@ -61,7 +61,6 @@ class StGmtStrip : public StObject { Short_t getMaxAdcTB() const { return mMaxAdcTB; } /// Maximal over the time bins Short_t getMaxPedSubtractedAdcTB() const { return mMaxPedSubtractedAdcTB; } - //Short_t getClusterSeedType() const { return mClusterSeedType; } /// Charge before GEM (in C) Float_t getCharge() const { return mCharge; } @@ -120,7 +119,6 @@ class StGmtStrip : public StObject { void setMaxAdc(Short_t adc) { mMaxAdc=adc; } /// Set maximal pedestal subtracked ADC over time buckets void setMaxPedSubtractedAdc(Short_t adc) { mMaxPedSubtractedAdc = adc; } - // void setClusterSeedType( Short_t seedType ) { mClusterSeedType=seedType; } /// Set charge before the GEM (in C) void setCharge( Float_t charge ) { mCharge = charge; } @@ -164,7 +162,6 @@ class StGmtStrip : public StObject { Short_t mMaxAdcTB; /// Max over the time bins Short_t mMaxPedSubtractedAdcTB; - // Short_t mClusterSeedType; // See types in StGmtConsts.h /// Charge before GEM, units (C) /// relation: ADC = ped + charge*gain(r,phi,disc) Float_t mCharge; diff --git a/StRoot/StEvent/StGmtStripCollection.h b/StRoot/StEvent/StGmtStripCollection.h index da5decce669..136bd7eef5a 100644 --- a/StRoot/StEvent/StGmtStripCollection.h +++ b/StRoot/StEvent/StGmtStripCollection.h @@ -65,7 +65,7 @@ class StGmtStripCollection : public StObject { // strip. New strip is created if it does not exist, but only // using StGmtStrip() constructor. Ownership is retained by the // collection. -// StGmtStrip* getStrip( int elecId ); + StGmtStrip* getStrip( int Id ); StGmtStrip* getSortedStrip( int Id ); @@ -95,26 +95,12 @@ class StGmtStripCollection : public StObject { // inline functions - -// inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { -// mStripElecIdVec.resize( kGmtNumElecIds ); -// for (unsigned int i=0; i(0); -// }; - inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { mStripGeoIdVec.resize( kGmtNumGeoIds ); for (unsigned int i=0; i(0); }; -// inline StGmtStripCollection::StGmtStripCollection( short module ) : StObject(), mModule( module ) { -// mStripVec.resize( kGmtNumGeoIds ); -// for (unsigned int i=0; i(0); -// }; - - // sort by geoId inline void StGmtStripCollection::sortByGeoId(){ sort( mStripVec.begin(), mStripVec.end(), &StGmtStripCollection::hitGeoIdLessThan ); diff --git a/StRoot/StGmtAligner/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx similarity index 100% rename from StRoot/StGmtAligner/EventT.cxx rename to StRoot/StGmtAligmentMaker/EventT.cxx diff --git a/StRoot/StGmtAligner/EventT.h b/StRoot/StGmtAligmentMaker/EventT.h similarity index 94% rename from StRoot/StGmtAligner/EventT.h rename to StRoot/StGmtAligmentMaker/EventT.h index 17cd2c78872..53aff37aad2 100644 --- a/StRoot/StGmtAligner/EventT.h +++ b/StRoot/StGmtAligmentMaker/EventT.h @@ -89,24 +89,29 @@ class EventT : public TObject { virtual void Print(Option_t *opt="") const; ClassDef(EventT,1) //EventT structure }; + class TBase { -public : + + public : TTree *fChain; //!pointer to the analyzed TTree or TChain Int_t fCurrent; //!current Tree number in a TChain EventT *fEvent; TString fOutFileName; - TBase(TTree *tree = 0) : fEvent(0) { + + TBase( TTree *tree = 0, + const Char_t *f_name = "/star/data09/calib/fisyak/Pass112/TpcSsd/065/Event_6065045_raw_1010001.root") : fEvent(0) { // if parameter tree is not specified (or zero), connect the file // used to generate this class and read the Tree. if (tree == 0) { - TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("/star/data09/calib/fisyak/Pass112/TpcSsd/065/Event_6065045_raw_1010001.root"); + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(f_name); if (!f) { - f = new TFile("/star/data09/calib/fisyak/Pass112/TpcSsd/065/Event_6065045_raw_1010001.root"); + f = new TFile(f_name); } tree = (TTree*)gDirectory->Get("T"); } Init(tree); } + virtual ~TBase() {if (!fChain) return; delete fChain->GetCurrentFile();} virtual Int_t Cut(Long64_t entry); virtual Int_t GetEntry(Long64_t entry) { if (!fChain) return 0; return fChain->GetEntry(entry);} diff --git a/StRoot/StGmtAligner/HitT.h b/StRoot/StGmtAligmentMaker/HitT.h similarity index 100% rename from StRoot/StGmtAligner/HitT.h rename to StRoot/StGmtAligmentMaker/HitT.h diff --git a/StRoot/StGmtAligner/StGmtAligner.cxx b/StRoot/StGmtAligmentMaker/StGmtAligner.cxx similarity index 100% rename from StRoot/StGmtAligner/StGmtAligner.cxx rename to StRoot/StGmtAligmentMaker/StGmtAligner.cxx diff --git a/StRoot/StGmtAligner/StGmtAligner.h b/StRoot/StGmtAligmentMaker/StGmtAligner.h similarity index 100% rename from StRoot/StGmtAligner/StGmtAligner.h rename to StRoot/StGmtAligmentMaker/StGmtAligner.h diff --git a/StRoot/StGmtAligner/TrackT.h b/StRoot/StGmtAligmentMaker/TrackT.h similarity index 100% rename from StRoot/StGmtAligner/TrackT.h rename to StRoot/StGmtAligmentMaker/TrackT.h diff --git a/StRoot/StGmtAligner/runGmtTree.C b/StRoot/StGmtAligmentMaker/runGmtTree.C similarity index 100% rename from StRoot/StGmtAligner/runGmtTree.C rename to StRoot/StGmtAligmentMaker/runGmtTree.C diff --git a/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx index 180b4455366..f59b8871453 100644 --- a/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx +++ b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx @@ -15,7 +15,7 @@ #include "StMessMgr.h" ClassImp(StGmtClusterMaker); -Int_t StGmtClusterMaker::gmtStat = 0; +int StGmtClusterMaker::gmtStat = 0; const UInt_t CLUS_BINS = 128; const Double_t CLUS_MIN = 0.0; const Double_t CLUS_MAX = 128*0.08; From 0a155631de216e44d15b1808acad87a766190c2c Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 13:45:52 -0500 Subject: [PATCH 03/11] Revert changes due to the inheritance issue --- StRoot/StEvent/StGmtHit.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/StRoot/StEvent/StGmtHit.h b/StRoot/StEvent/StGmtHit.h index 3049f8f06b5..4296f6b4492 100644 --- a/StRoot/StEvent/StGmtHit.h +++ b/StRoot/StEvent/StGmtHit.h @@ -55,13 +55,13 @@ class StGmtHit : public StHit { /// ADC error in Y Float_t getErrorAdcY() const { return mdAdcY; } /// Local X coordinate - Float_t getLocalX() const { return mPosition().x(); } + Float_t getLocalX() const { return position().x(); } /// Local X coordinate error - Float_t getErrorLocalX() const { return mPositionError().x(); } + Float_t getErrorLocalX() const { return positionError().x(); } /// Local Y coordinate - Float_t getLocalY() const { return mPosition().y(); } + Float_t getLocalY() const { return position().y(); } /// Local Y coordinate error - Float_t getErrorLocalY() const { return mPositionError().y(); } + Float_t getErrorLocalY() const { return positionError().y(); } /// Position in X Float_t getSigmaX() const { return mSigmaX; } /// Position error in X From 875a7e421a26dd4d9a9e2952480385a7bfaf8791 Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 15:35:55 -0500 Subject: [PATCH 04/11] Fixing EventT --- StRoot/StGmtAligmentMaker/EventT.cxx | 32 +++++++++++++++------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/StRoot/StGmtAligmentMaker/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx index d6f92daec06..8a04e2acbda 100644 --- a/StRoot/StGmtAligmentMaker/EventT.cxx +++ b/StRoot/StGmtAligmentMaker/EventT.cxx @@ -64,7 +64,9 @@ EventT::~EventT() //______________________________________________________________________________ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { +#if 0 static const Int_t NoFitPointCutForGoodTrackT = 15; +#endif Clear(); Int_t iok = 1; fIsValid = kFALSE; @@ -404,8 +406,8 @@ void EventT::RestoreListOfRotations() { void TBase::Loop(Int_t Nevents) { #if 1 struct PlotPar_t { - Char_t *Name; - Char_t *Title; + const Char_t *Name; + const Char_t *Title; Int_t nx; Int_t ny; Double_t xmin; @@ -413,14 +415,12 @@ void TBase::Loop(Int_t Nevents) { Double_t ymin; Double_t ymax; }; - const PlotPar_t plotUP = // plots for uP - { "uP","track u", 320, 3, -5.,5., 0.,3. }; - - const PlotPar_t plotDu = // plots for u-uP - { "Du","Du before cut", 250, 3, -2.,2., 0.,3. }; - - const PlotPar_t plotDuv = // plots for du & dv - { "Du","Du cuts", 200, 3, -1.,1., 0.,3. }; + // plots for uP + const PlotPar_t plotUP = { "uP", "track u", 320, 3, -5., 5., 0., 3. }; + // plots for u-uP + const PlotPar_t plotDu = { "Du", "Du before cut", 250, 3, -2., 2., 0., 3. }; + // plots for du & dv + const PlotPar_t plotDuv = { "Du", "Du cuts", 200, 3, -1., 1., 0., 3. }; TFile *fOut = new TFile(fOutFileName,"recreate"); TString Name; @@ -431,10 +431,12 @@ void TBase::Loop(Int_t Nevents) { // B TH1F *LocPlots[NM]; TH1F * uPlots[NM]; - TH1F * hpT = new TH1F( "Pt", "pt", 200, -2., 2.); - TH1F * hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); +#if 0 + TH1F *hpT = new TH1F( "Pt", "pt", 200, -2., 2.); + TH1F *hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); + TH1F *uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); +#endif TH1F * uPAll = new TH1F("UPall","uPall", plotUP.nx, plotUP.xmin, plotUP.xmax); - TH1F * uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); TH1F * duB[NM][2]; TH1F * dvB[NM]; TH1F * uCuts[NM]; @@ -444,7 +446,7 @@ void TBase::Loop(Int_t Nevents) { TH2F * dMin = new TH2F("DMin","vumin",100,-0.75,0.75,100,-0.75,0.75); TH1F * vMin = new TH1F("VMin","vmin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); TH1F * uMin = new TH1F("UMin","umin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); - TH1F * uMinC = new TH1F("UMinC","umC", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + //TH1F * uMinC = new TH1F("UMinC","umC", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); memset(LocPlots,0,NM*sizeof(TH1F *)); memset( uPlots,0,NM*sizeof(TH1F *)); for (int M = 0; M < NM; M++) {// over gmt Modules @@ -496,8 +498,8 @@ void TBase::Loop(Int_t Nevents) { for (UInt_t trk = 0; trk < Ntrack; trk++) { TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); if (! track) continue; - Int_t Npoints = track->GetNpoint(); #if 0 + Int_t Npoints = track->GetNpoint(); if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue; if (UseSsd && Npoints < 1000) continue; if (UseSvt && Npoints < 100) continue; From 20dffdf6447a7170f1a6afe70bb5ada06bd27cff Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 15:40:36 -0500 Subject: [PATCH 05/11] Moving from the ROOT to C++ type declaration --- .../StGmtClusterMaker/StGmtClusterMaker.cxx | 30 +++++++++++-------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx index f59b8871453..a704e9512ea 100644 --- a/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx +++ b/StRoot/StGmtClusterMaker/StGmtClusterMaker.cxx @@ -16,16 +16,17 @@ ClassImp(StGmtClusterMaker); int StGmtClusterMaker::gmtStat = 0; -const UInt_t CLUS_BINS = 128; -const Double_t CLUS_MIN = 0.0; -const Double_t CLUS_MAX = 128*0.08; -const UInt_t MAX_PEAKS = 20; +const unsigned int CLUS_BINS = 128; +const double CLUS_MIN = 0.0; +const double CLUS_MAX = 128*0.08; +const unsigned int MAX_PEAKS = 20; -//________________________________________________________________________________ +//_________________ inline Double_t MyGaus(Double_t x, Double_t mean, Double_t sigma, Double_t delta) { return TMath::Freq((x-mean+delta/2)/sigma)-TMath::Freq((x-mean-delta/2)/sigma); } -//________________________________________________________________________________ + +//_________________ Double_t fpeaks(Double_t *x, Double_t *par) { Float_t result=0.0; UInt_t nPar=(UInt_t)par[0]; @@ -37,7 +38,8 @@ Double_t fpeaks(Double_t *x, Double_t *par) { } return result; } -//________________________________________________________________________________ + +//_________________ TF1* StGmtClusterMaker::FindPeaks(TH1F* hist) { TSpectrum spect(MAX_PEAKS); TF1 back("poly","pol0",CLUS_MIN,CLUS_MAX); @@ -90,7 +92,8 @@ TF1* StGmtClusterMaker::FindPeaks(TH1F* hist) { return fitFunc; } -//________________________________________________________________________________ + +//_________________ void StGmtClusterMaker::ClusterBuilder(ULong_t events, UInt_t module, StGmtStripCollection& strips, StGmtHitCollection& hits) { static TCanvas* canv=0; static TH1F* histX=0; @@ -216,12 +219,14 @@ void StGmtClusterMaker::ClusterBuilder(ULong_t events, UInt_t module, StGmtStrip } } } -//________________________________________________________________________________ + +//_________________ StGmtClusterMaker::StGmtClusterMaker( const Char_t* name ) : //StMaker(name), StRTSBaseMaker( "clustser", name ) { SetAttr("gmtCosmics" ,kFALSE); -}; -//________________________________________________________________________________ +} + +//_________________ Int_t StGmtClusterMaker::Make() { LOG_INFO << "MAKE of StGmtClusterMaker" << endm; Int_t ierr = kStOk; @@ -294,7 +299,8 @@ Int_t StGmtClusterMaker::Make() { } return ierr; } -//________________________________________________________________________________ + +//_________________ Int_t StGmtClusterMaker::Init() { LOG_INFO << "INTI of StGmtClusterMaker" << endm; if (IAttr("gmtCosmics")) SetAttr(".Privilege",kTRUE); From b1fb23be8a25d3279e1dd7dbc9faae37aa9d551d Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 16:21:32 -0500 Subject: [PATCH 06/11] Update missing parts --- StRoot/StGmtAligmentMaker/EventT.cxx | 154 +++++++++++++++------------ 1 file changed, 83 insertions(+), 71 deletions(-) diff --git a/StRoot/StGmtAligmentMaker/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx index 8a04e2acbda..5b2868dbc2b 100644 --- a/StRoot/StGmtAligmentMaker/EventT.cxx +++ b/StRoot/StGmtAligmentMaker/EventT.cxx @@ -38,11 +38,10 @@ TClonesArray *EventT::fgTracks = 0; TClonesArray *EventT::fgHits = 0; THashList *EventT::fRotList = 0; +static int _debug = 0; -static Int_t _debug = 0; -//______________________________________________________________________________ -EventT::EventT() : fIsValid(kFALSE) -{ +//_________________ +EventT::EventT() : fIsValid(kFALSE) { // Create an EventT object. // When the constructor is invoked for the first time, the class static // variable fgTracks is 0 and the TClonesArray fgTracks is created. @@ -55,14 +54,13 @@ EventT::EventT() : fIsValid(kFALSE) fNhit = 0; } -//______________________________________________________________________________ -EventT::~EventT() -{ +//_________________ +EventT::~EventT() { Clear(); SafeDelete(fRotList); } -//______________________________________________________________________________ +//_________________ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { #if 0 static const Int_t NoFitPointCutForGoodTrackT = 15; @@ -248,7 +246,7 @@ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { if (nTotMatch) iok = 0; return iok; } -//______________________________________________________________________________ +//_________________ TrackT *EventT::AddTrackT() { // Add a new track to the list of tracks for this event. @@ -262,7 +260,7 @@ TrackT *EventT::AddTrackT() //Save reference to last TrackT in the collection of Tracks return track; } -//______________________________________________________________________________ +//_________________ HitT *EventT::AddHitT() { // Add a new hit to the list of hits for this event. @@ -277,14 +275,14 @@ HitT *EventT::AddHitT() return hit; } -//______________________________________________________________________________ +//_________________ void EventT::Clear(Option_t * /*option*/) { fTracks->Clear("C"); //will also call TrackT::Clear fHits->Clear("C"); //will also call HitT::Clear } -//______________________________________________________________________________ +//_________________ void EventT::Reset(Option_t * /*option*/) { // Static function to reset all static objects for this event @@ -294,14 +292,14 @@ void EventT::Reset(Option_t * /*option*/) delete fgHits; fgHits = 0; } -//______________________________________________________________________________ +//_________________ void EventT::SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field) { fNtrack = 0; fNhit = 0; fEvtHdr.Set(i, run, date, field); } -//________________________________________________________________________________ +//___________________ void EventT::Print(Option_t *opt) const { cout << "Run/EventT\t" << fEvtHdr.GetRun() << "/" << fEvtHdr.GetEvtNum() << "\tDate " << fEvtHdr.GetDate() << "\tField " << fEvtHdr.GetField() << endl; @@ -315,7 +313,7 @@ void EventT::Print(Option_t *opt) const { for (UInt_t i = 0; i < GetNhit(); i++) {cout << i << "\t"; GetHitT(i)->Print();} } -//________________________________________________________________________________ +//___________________ HitT *EventT::SetHitT(HitT *h, StHit *hit, TGeoHMatrix *comb, TrackT *track) { UInt_t B = 0, L = 0, l = 0, W = 0, H = 0; Int_t rdo = 0; @@ -368,7 +366,7 @@ HitT *EventT::SetHitT(HitT *h, StHit *hit, TGeoHMatrix *comb, TrackT *track) { } return h; } -//________________________________________________________________________________ +//___________________ void TrackT::Print(Option_t *opt) const { cout << "TrackT: InvpT " << fInvpT << "\tTanL " << fTanL << "\tPhi " << fPhi << "\tRho " << fRho @@ -376,18 +374,18 @@ void TrackT::Print(Option_t *opt) const { for (UInt_t i = 0; i < fNsp; i++) cout << "\t" << fIdHitT[i]; cout << endl; } -//________________________________________________________________________________ +//___________________ void HitT::SetId(Int_t B, Int_t L, Int_t l, Int_t W, Int_t H) { barrel = B; layer = L; ladder = l; wafer = W; hybrid = H; } -//________________________________________________________________________________ +//___________________ void HitT::Print(Option_t *opt) const { cout << "HitT: Id " << Id << "\tpT = " << pT << "\tmomentum " << pMom << endl; TRVector glob(3,&xG); cout << "Global :" << glob << endl; cout << "Local u/v/w " << u << "/ " << v << "/ " << w << endl; cout << "Prediction uP/vP " << uP << "/ " << vP << "\ttuP/tvP " << tuP << "/ " << tvP << endl; } -//________________________________________________________________________________ +//___________________ void EventT::RestoreListOfRotations() { if (fRotList) return; if (! gDirectory) return; @@ -402,9 +400,12 @@ void EventT::RestoreListOfRotations() { } } } -//________________________________________________________________________________ -void TBase::Loop(Int_t Nevents) { -#if 1 + +//___________________ +void TBase::Loop(Int_t Nevents) { + +//#if 1 + struct PlotPar_t { const Char_t *Name; const Char_t *Title; @@ -463,19 +464,20 @@ void TBase::Loop(Int_t Nevents) { uName = plotUP.Name; uName += Form("M%i", module); uTitle = Form("uP for Module %i", module); - uPlots[module] = - new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); + uPlots[module] = new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); uName = Form("%sM%i", plotDu.Name, module); uTitle = Form("u-uP for M %i", module); - uCuts[module]= - new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); uName = Form("%sxM%i", plotDu.Name, module); uTitle = Form("u-uP corr M %i", module); - xCuts[module]= - new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); - } + xCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + } // for (int M = 0; M < NM; M++) + Long64_t nentries = fChain->GetEntriesFast(); - if (Nevents > 0 && nentries > Nevents) nentries = Nevents; + if (Nevents > 0 && nentries > Nevents) { + nentries = Nevents; + } + Long64_t nbytes = 0, nb = 0; Int_t TreeNo = -1; TString currentFile(""); @@ -485,16 +487,18 @@ void TBase::Loop(Int_t Nevents) { if (ientry < 0) break; nb = fChain->GetEntry(jentry); nbytes += nb; if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) { - cout << "Read event \t" << jentry + std::cout << "Read event \t" << jentry << " so far, switch to file " << fChain->GetCurrentFile()->GetName() - << endl; - cout << " current TreeNo: " << TreeNo - << " new TreeNo: " << fChain->GetTreeNumber() << endl; + << std::endl; + std::cout << " current TreeNo: " << TreeNo + << " new TreeNo: " << fChain->GetTreeNumber() << std::endl; TreeNo = fChain->GetTreeNumber(); - } + } // for (Long64_t jentry=0; jentry 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue; UInt_t Ntrack = fEvent->GetTracks()->GetEntriesFast(); - // int k_used[100000] = {0}; + // int k_used[100000] = {0}; + for (UInt_t trk = 0; trk < Ntrack; trk++) { TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); if (! track) continue; @@ -507,47 +511,55 @@ void TBase::Loop(Int_t Nevents) { Int_t Nsp = track->GetN(); double dvmin = 1000.; double dumin = 1000.; - int kmin; + int kmin{0}; for (Int_t hit = 0; hit < Nsp; hit++) { - Int_t k = track->GetHitTId(hit) - 1; - // assert(k>=0); - HitT *hitT = (HitT *) fEvent->GetHitT(k); - if ( k < 0) cout <<" k <0:"<GetEntriesFast(); From be2bf1c87170984cbc0a99d1d62a8217ae9079bf Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 22 May 2024 16:42:21 -0500 Subject: [PATCH 08/11] Stylish modifications --- StRoot/StGmtAligmentMaker/EventT.cxx | 944 ++++++++++++++------------- 1 file changed, 488 insertions(+), 456 deletions(-) diff --git a/StRoot/StGmtAligmentMaker/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx index 8240da87f18..88bac0d34d8 100644 --- a/StRoot/StGmtAligmentMaker/EventT.cxx +++ b/StRoot/StGmtAligmentMaker/EventT.cxx @@ -42,363 +42,390 @@ static int _debug = 0; //_________________ EventT::EventT() : fIsValid(kFALSE) { - // Create an EventT object. - // When the constructor is invoked for the first time, the class static - // variable fgTracks is 0 and the TClonesArray fgTracks is created. + // Create an EventT object. + // When the constructor is invoked for the first time, the class static + // variable fgTracks is 0 and the TClonesArray fgTracks is created. - if (!fgTracks) fgTracks = new TClonesArray("TrackT", 1000); - fTracks = fgTracks; - fNtrack = 0; - if (!fgHits) fgHits = new TClonesArray("HitT", 1000); - fHits = fgHits; - fNhit = 0; + if (!fgTracks) fgTracks = new TClonesArray("TrackT", 1000); + fTracks = fgTracks; + fNtrack = 0; + if (!fgHits) fgHits = new TClonesArray("HitT", 1000); + fHits = fgHits; + fNhit = 0; } //_________________ EventT::~EventT() { - Clear(); - SafeDelete(fRotList); + Clear(); + SafeDelete(fRotList); } //_________________ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { + #if 0 - static const Int_t NoFitPointCutForGoodTrackT = 15; + static const Int_t NoFitPointCutForGoodTrackT = 15; #endif - Clear(); - Int_t iok = 1; - fIsValid = kFALSE; - if (! pEventT) return iok; - StGmtCollection* GmtCollection = pEventT->gmtCollection(); - if (! GmtCollection) { cout << "No GMT Collections" << endl; return iok;} - StSPtrVecTrackNode& theNodes = pEventT->trackNodes(); - UInt_t nnodes = theNodes.size(); - if (! nnodes) { cout << "No tracks" << endl; return iok;} - StEventInfo* info = pEventT->info(); - Int_t ev = 0, run = 0, time = 0; - if (info) { - ev = info->id(); - run = info->runId(); - time = info->time(); - } - StEventSummary* summary = pEventT->summary(); - Double32_t field = 0; - if (summary) field = summary->magneticField(); - SetHeader(ev,run,time,field); - SetFlag(1); - Int_t nTotMatch = 0; - // Create and Fill the TrackT objects - for (UInt_t i=0; itrack(global); - if (! Track) continue; - StGlobalTrack *gTrack = (StGlobalTrack *) Track; - StDcaGeometry *dca = gTrack->dcaGeometry(); - if (! dca) continue; - StPrimaryTrack *pTrack = (StPrimaryTrack *) node->track(primary); + Clear(); + Int_t iok = 1; + fIsValid = kFALSE; + if (! pEventT) return iok; + StGmtCollection* GmtCollection = pEventT->gmtCollection(); + if (! GmtCollection) { + std::cout << "No GMT Collections" << std::endl; + return iok; + } + StSPtrVecTrackNode& theNodes = pEventT->trackNodes(); + UInt_t nnodes = theNodes.size(); + if (! nnodes) { + std::cout << "No tracks" << std::endl; + return iok; + } + StEventInfo* info = pEventT->info(); + Int_t ev = 0, run = 0, time = 0; + if (info) { + ev = info->id(); + run = info->runId(); + time = info->time(); + } + StEventSummary* summary = pEventT->summary(); + Double32_t field = 0; + if (summary) field = summary->magneticField(); + SetHeader(ev,run,time,field); + SetFlag(1); + Int_t nTotMatch = 0; + // Create and Fill the TrackT objects + for (UInt_t i=0; itrack(global); + if (! Track) continue; + StGlobalTrack *gTrack = (StGlobalTrack *) Track; + StDcaGeometry *dca = gTrack->dcaGeometry(); + if (! dca) continue; + StPrimaryTrack *pTrack = (StPrimaryTrack *) node->track(primary); #ifndef __OnlyGlobal__ - if (! pTrack) continue; - Track = (StTrack *) pTrack; + if (! pTrack) continue; + Track = (StTrack *) pTrack; #endif - StThreeVectorD g3 = Track->geometry()->momentum(); - Double_t pT = g3.perp(); - Double_t pMom = g3.mag(); - if (pMom < pCut) continue; - TrackT *track = AddTrackT(); - Double_t InvpT = 0; - Double_t TanL = 999999; - if (TMath::Abs(pT) > 1.e-7) { - InvpT = Track->geometry()->charge()/pT; - TanL = g3.z()/pT; - } - track->SetInvpT(InvpT); - track->SetPhi(TMath::ATan2(g3.y(),g3.x())); - track->SetTanL(TanL); - static const Double_t EC = 2.9979251E-4; - Double_t Rho = - EC*InvpT*field; - track->SetRho(Rho); - track->SetLength(Track->length()); - StTrackDetectorInfo* dinfo=Track->detectorInfo(); - track->SetNpoint(dinfo->numberOfPoints()); - track->SetNPpoint(Track->numberOfPossiblePoints()); - track->SetN(0); - StPhysicalHelixD helixO = Track->outerGeometry()->helix(); - Double_t R_tof[2]= {210., 216.}; // inner and outer surfaces - pair shR = helixO.pathLength(0.5*(R_tof[0]+R_tof[1])); - if (TMath::Abs(shR.first) > 200 && TMath::Abs(shR.second) > 200) continue; - Double_t stepR = (shR.first > 0) ? shR.first : shR.second; - StThreeVectorD xyzR = helixO.at(stepR); - Double_t phiR = TMath::RadToDeg()*xyzR.phi(); - if (_debug) - cout << "\t shR " << shR.first << "\t" << shR.second << "\tstepR " << stepR - << "\txyzR\t" << xyzR << "\tphiR\t" << phiR << endl; - Int_t NoHitPerTrack = 0; - THashList *fRotList = RotMatrices(); - if (! fRotList) continue; - TIter next(fRotList); - TGeoHMatrix *comb = 0; + StThreeVectorD g3 = Track->geometry()->momentum(); + Double_t pT = g3.perp(); + Double_t pMom = g3.mag(); + if (pMom < pCut) continue; + TrackT *track = AddTrackT(); + Double_t InvpT = 0; + Double_t TanL = 999999; + if (TMath::Abs(pT) > 1.e-7) { + InvpT = Track->geometry()->charge()/pT; + TanL = g3.z()/pT; + } + track->SetInvpT(InvpT); + track->SetPhi(TMath::ATan2(g3.y(),g3.x())); + track->SetTanL(TanL); + static const Double_t EC = 2.9979251E-4; + Double_t Rho = - EC*InvpT*field; + track->SetRho(Rho); + track->SetLength(Track->length()); + StTrackDetectorInfo* dinfo=Track->detectorInfo(); + track->SetNpoint(dinfo->numberOfPoints()); + track->SetNPpoint(Track->numberOfPossiblePoints()); + track->SetN(0); + StPhysicalHelixD helixO = Track->outerGeometry()->helix(); + Double_t R_tof[2]= {210., 216.}; // inner and outer surfaces + pair shR = helixO.pathLength(0.5*(R_tof[0]+R_tof[1])); - while ((comb = (TGeoHMatrix *) next())) { - TString combName(comb->GetName()); - if (! combName.BeginsWith("R")) continue; - Int_t Id; - sscanf(comb->GetName()+1,"%i",&Id); - UInt_t module = Id; - StGmtHitCollection* GmtHitCollection = GmtCollection->getHitCollection(module); - if (! GmtHitCollection) { cout << "No GMT HitT Collections for mudule " << module << endl; continue;} - StSPtrVecGmtHit& hitvec = GmtHitCollection->getHitVec(); - UInt_t NoHits = hitvec.size(); - if (! NoHits) continue; - if (_debug) { - cout << comb->GetName() << "\tmodule = " << module << endl; - comb->Print(); - } - static Double_t dz[2] = {50.00, 2.10}; - static Double_t dx[2] = {50.00, 3.65}; - Int_t k = 0; // gmt - Double_t *rot = comb->GetRotationMatrix(); - Double_t *tra = comb->GetTranslation(); - const StThreeVectorD unit(1.,0.,0.); - const StThreeVectorD zero(0.,0.,0.); - StThreeVectorD normal(rot[2], rot[5], rot[8]); - StThreeVectorD middle(tra); - if (_debug) cout << "middle:" << middle << "\tnormal:" << normal << endl; - comb->LocalToMaster(zero.xyz(),middle.xyz()); - comb->LocalToMasterVect(unit.xyz(), normal.xyz()); - if (_debug) cout << "middle:" << middle << "\tnormal:" << normal << endl; - Double_t phiM = TMath::RadToDeg()*middle.phi(); - if (_debug) cout << "phiR = " << phiR << "\tphiM = " << phiM << endl; - Double_t dPhi = phiR - phiM; - if (dPhi > 360) dPhi -= 360; - if (dPhi < -360) dPhi += 360; - if (TMath::Abs(dPhi) > 15) continue; - if (_debug) cout << "zR = " << xyzR.z() << "\tzM = " << tra[2] << endl; - if (TMath::Abs(xyzR.z() - tra[2]) > 20) continue; - Double_t sh = helixO.pathLength(middle, normal); - if (_debug) { - cout << "StHelix sh " << sh - << "\t shR " << shR.first << "\t" << shR.second - << endl; - } - StThreeVectorD xyzG = helixO.at(sh); if (_debug) cout << "StHelix xyzG\t" << xyzG << endl; - StThreeVectorD dR = xyzR - xyzG; if (_debug) cout << "dR\t" << dR << " dist = " << dR.magnitude() << endl; - if (dR.magnitude() > 50) continue; - if (sh < -5e2 || sh > 5e2) continue; - if (_debug) { - StThreeVectorD dX = xyzG - helixO.at(0); - cout << "Qi: " << Track->geometry()->charge() - << "\tQo: " << Track->outerGeometry()->charge() - << "\tdX " << dX << endl; - cout << *dca << endl; - } - Double_t uvPred[3]; - comb->MasterToLocal(xyzG.xyz(),uvPred); - TRVector xyzL(3,uvPred); if (_debug) cout << "StHelix xyzL\t" << xyzL << endl; - Double_t dirGPred[3] = {helixO.cx(sh),helixO.cy(sh),helixO.cz(sh)}; - Double_t dxyzL[3]; - comb->MasterToLocalVect(dirGPred,dxyzL); - Double_t tuvPred[2] = {dxyzL[1]/dxyzL[0], dxyzL[2]/dxyzL[0]}; - if (_debug) cout << "StHelix tU/tV = " << tuvPred[0] << "\t" << tuvPred[1] << endl; - if (TMath::Abs(uvPred[1]) > dx[k] + 1.0) continue; - if (TMath::Abs(uvPred[2]) > dz[k] + 1.0) continue; + if (TMath::Abs(shR.first) > 200 && TMath::Abs(shR.second) > 200) continue; - Bool_t mIsCrossingMembrain = kTRUE; - Bool_t mIsPrimary = kFALSE; - if(pTrack!=0) { - mIsPrimary = kTRUE; - } - StThreeVectorF firstPoint = Track->detectorInfo()->firstPoint(); - StThreeVectorF lastPoint = Track->detectorInfo()->lastPoint(); - if( (firstPoint.z()>0 && lastPoint.z()>0 && xyzG.z()>0) || - (firstPoint.z()<0 && lastPoint.z()<0 && xyzG.z()<0) ) { - mIsCrossingMembrain = kFALSE; - } - - for (UInt_t l = 0; l < NoHits; l++) { - StHit *hit = hitvec[l]; - if (hit) { - //if (hit->flag()>=4) continue; - //if (hit->flag()< 0) continu; - // cout << "hitFlag=" << hit->flag() << endl; - - HitT *h = AddHitT(); - h->SetHitLength(sh); - h->SetHitLength(stepR); - h->SetHitdR(dR.magnitude()); - h->SetHitFlag(UInt_t(hit->flag())); - h->SetUVPred (uvPred[1],uvPred[2]); - h->SettUVPred(tuvPred[0],tuvPred[1]); - h->SetXyzG(xyzG.xyz()); - h->SetDirG(dirGPred); - h->SetisPrimary(mIsPrimary); - h->SetisCrossingMembrain(mIsCrossingMembrain); - SetHitT(h, hit, comb, track); - NoHitPerTrack++; - h->SetHitPerTrack(NoHitPerTrack); - // SetHitT(h, hit, comb, track, &TPDeriv); - } - } - } - nTotMatch += NoHitPerTrack; - } - fIsValid = kTRUE; - if (nTotMatch) iok = 0; - return iok; + Double_t stepR = (shR.first > 0) ? shR.first : shR.second; + StThreeVectorD xyzR = helixO.at(stepR); + Double_t phiR = TMath::RadToDeg()*xyzR.phi(); + if (_debug) { + std::cout << "\t shR " << shR.first << "\t" << shR.second << "\tstepR " << stepR + << "\txyzR\t" << xyzR << "\tphiR\t" << phiR << std::endl; + } + Int_t NoHitPerTrack = 0; + THashList *fRotList = RotMatrices(); + if (! fRotList) continue; + TIter next(fRotList); + TGeoHMatrix *comb = 0; + + while ((comb = (TGeoHMatrix *) next())) { + TString combName(comb->GetName()); + if (! combName.BeginsWith("R")) continue; + Int_t Id; + sscanf(comb->GetName()+1,"%i",&Id); + UInt_t module = Id; + StGmtHitCollection* GmtHitCollection = GmtCollection->getHitCollection(module); + if (! GmtHitCollection) { cout << "No GMT HitT Collections for mudule " << module << endl; continue;} + StSPtrVecGmtHit& hitvec = GmtHitCollection->getHitVec(); + UInt_t NoHits = hitvec.size(); + if (! NoHits) continue; + if (_debug) { + std::cout << comb->GetName() << "\tmodule = " << module << std::endl; + comb->Print(); + } + static Double_t dz[2] = {50.00, 2.10}; + static Double_t dx[2] = {50.00, 3.65}; + Int_t k = 0; // gmt + Double_t *rot = comb->GetRotationMatrix(); + Double_t *tra = comb->GetTranslation(); + const StThreeVectorD unit(1.,0.,0.); + const StThreeVectorD zero(0.,0.,0.); + StThreeVectorD normal(rot[2], rot[5], rot[8]); + StThreeVectorD middle(tra); + if (_debug) { + std::cout << "middle:" << middle << "\tnormal:" << normal << std::endl; + } + comb->LocalToMaster(zero.xyz(),middle.xyz()); + comb->LocalToMasterVect(unit.xyz(), normal.xyz()); + if (_debug) { + std::cout << "middle:" << middle << "\tnormal:" << normal << std::endl; + } + Double_t phiM = TMath::RadToDeg()*middle.phi(); + if (_debug) { + std::cout << "phiR = " << phiR << "\tphiM = " << phiM << std::endl; + } + Double_t dPhi = phiR - phiM; + if (dPhi > 360) dPhi -= 360; + if (dPhi < -360) dPhi += 360; + if (TMath::Abs(dPhi) > 15) continue; + if (_debug) { + std::cout << "zR = " << xyzR.z() << "\tzM = " << tra[2] << std::endl; + } + if (TMath::Abs(xyzR.z() - tra[2]) > 20) continue; + Double_t sh = helixO.pathLength(middle, normal); + if (_debug) { + std::cout << "StHelix sh " << sh + << "\t shR " << shR.first << "\t" << shR.second + << std::endl; + } + StThreeVectorD xyzG = helixO.at(sh); if (_debug) cout << "StHelix xyzG\t" << xyzG << endl; + StThreeVectorD dR = xyzR - xyzG; if (_debug) cout << "dR\t" << dR << " dist = " << dR.magnitude() << endl; + if (dR.magnitude() > 50) continue; + if (sh < -5e2 || sh > 5e2) continue; + if (_debug) { + StThreeVectorD dX = xyzG - helixO.at(0); + std::cout << "Qi: " << Track->geometry()->charge() + << "\tQo: " << Track->outerGeometry()->charge() + << "\tdX " << dX << std::endl; + std::cout << *dca << std::endl; + } + Double_t uvPred[3]; + comb->MasterToLocal(xyzG.xyz(),uvPred); + TRVector xyzL(3,uvPred); if (_debug) cout << "StHelix xyzL\t" << xyzL << endl; + Double_t dirGPred[3] = {helixO.cx(sh),helixO.cy(sh),helixO.cz(sh)}; + Double_t dxyzL[3]; + comb->MasterToLocalVect(dirGPred,dxyzL); + Double_t tuvPred[2] = {dxyzL[1]/dxyzL[0], dxyzL[2]/dxyzL[0]}; + if (_debug) { + std::cout << "StHelix tU/tV = " << tuvPred[0] << "\t" << tuvPred[1] << std::endl; + } + if (TMath::Abs(uvPred[1]) > dx[k] + 1.0) continue; + if (TMath::Abs(uvPred[2]) > dz[k] + 1.0) continue; + + Bool_t mIsCrossingMembrain = kTRUE; + Bool_t mIsPrimary = kFALSE; + if(pTrack!=0) { + mIsPrimary = kTRUE; + } + StThreeVectorF firstPoint = Track->detectorInfo()->firstPoint(); + StThreeVectorF lastPoint = Track->detectorInfo()->lastPoint(); + if( (firstPoint.z()>0 && lastPoint.z()>0 && xyzG.z()>0) || + (firstPoint.z()<0 && lastPoint.z()<0 && xyzG.z()<0) ) { + mIsCrossingMembrain = kFALSE; + } + + for (UInt_t l = 0; l < NoHits; l++) { + StHit *hit = hitvec[l]; + if (hit) { + //if (hit->flag()>=4) continue; + //if (hit->flag()< 0) continu; + // std::cout << "hitFlag=" << hit->flag() << std::endl; + + HitT *h = AddHitT(); + h->SetHitLength(sh); + h->SetHitLength(stepR); + h->SetHitdR(dR.magnitude()); + h->SetHitFlag(UInt_t(hit->flag())); + h->SetUVPred (uvPred[1],uvPred[2]); + h->SettUVPred(tuvPred[0],tuvPred[1]); + h->SetXyzG(xyzG.xyz()); + h->SetDirG(dirGPred); + h->SetisPrimary(mIsPrimary); + h->SetisCrossingMembrain(mIsCrossingMembrain); + SetHitT(h, hit, comb, track); + NoHitPerTrack++; + h->SetHitPerTrack(NoHitPerTrack); + // SetHitT(h, hit, comb, track, &TPDeriv); + } // if (hit) + } // for (UInt_t l = 0; l < NoHits; l++) + } // while ((comb = (TGeoHMatrix *) next())) + nTotMatch += NoHitPerTrack; + } // for (UInt_t i=0; iClear("C"); //will also call TrackT::Clear - fHits->Clear("C"); //will also call HitT::Clear +void EventT::Clear(Option_t * /*option*/) { + fTracks->Clear("C"); //will also call TrackT::Clear + fHits->Clear("C"); //will also call HitT::Clear } //_________________ -void EventT::Reset(Option_t * /*option*/) -{ - // Static function to reset all static objects for this event - // fgTracks->Delete(option); - - delete fgTracks; fgTracks = 0; - delete fgHits; fgHits = 0; +void EventT::Reset(Option_t * /*option*/) { + // Static function to reset all static objects for this event + // fgTracks->Delete(option); + + delete fgTracks; fgTracks = 0; + delete fgHits; fgHits = 0; } //_________________ -void EventT::SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field) -{ - fNtrack = 0; - fNhit = 0; - fEvtHdr.Set(i, run, date, field); +void EventT::SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field) { + fNtrack = 0; + fNhit = 0; + fEvtHdr.Set(i, run, date, field); } + //___________________ void EventT::Print(Option_t *opt) const { - cout << "Run/EventT\t" << fEvtHdr.GetRun() << "/" << fEvtHdr.GetEvtNum() << "\tDate " << fEvtHdr.GetDate() - << "\tField " << fEvtHdr.GetField() << endl; - cout << "Total no. tracks " << GetTotalNoTracks() << "\tRecorded tracks " << GetNtrack() - << "\tRecorded hits " << GetNhit() << endl; - TRVector vertex(3,GetVertex()); - TRSymMatrix cov(3,GetCovMatrix()); - cout << "Primary vertex " << vertex << endl; - cout << "Its cov. matrix " << cov << endl; - for (UInt_t i = 0; i < GetNtrack(); i++) {cout << i << "\t"; GetTrackT(i)->Print();} - for (UInt_t i = 0; i < GetNhit(); i++) {cout << i << "\t"; GetHitT(i)->Print();} - + std::cout << "Run/EventT\t" << fEvtHdr.GetRun() << "/" << fEvtHdr.GetEvtNum() << "\tDate " << fEvtHdr.GetDate() + << "\tField " << fEvtHdr.GetField() << std::endl; + std::cout << "Total no. tracks " << GetTotalNoTracks() << "\tRecorded tracks " << GetNtrack() + << "\tRecorded hits " << GetNhit() << std::endl; + TRVector vertex(3,GetVertex()); + TRSymMatrix cov(3,GetCovMatrix()); + std::cout << "Primary vertex " << vertex << std::endl; + std::cout << "Its cov. matrix " << cov << std::endl; + for (UInt_t i = 0; i < GetNtrack(); i++) { + std::cout << i << "\t"; GetTrackT(i)->Print(); + } + for (UInt_t i = 0; i < GetNhit(); i++) { + std::cout << i << "\t"; GetHitT(i)->Print(); + } } + //___________________ HitT *EventT::SetHitT(HitT *h, StHit *hit, TGeoHMatrix *comb, TrackT *track) { - UInt_t B = 0, L = 0, l = 0, W = 0, H = 0; - Int_t rdo = 0; - h->SetRDO(rdo); - if (hit->detector() == kGmtId) { - StGmtHit *ht = (StGmtHit *) hit; - B = ht->getModule(); - h->SetId(B,L,l,W,H); - h->SetuvD(ht->getLocalY(), ht->getLocalX()); - h->SetuvDError(ht->getErrorLocalY(), ht->getErrorLocalX()); - h->SetSigma(ht->getSigmaY(), ht->getSigmaX()); - h->SetSigmaError(ht->getErrorSigmaY(), ht->getErrorSigmaX()); - h->SetAdc(ht->getAdcY(), ht->getAdcX()); - h->SetAdcError(ht->getErrorAdcY(), ht->getErrorAdcX()); - h->SetUsedInFit(hit->usedInFit()); - } - StThreeVectorF position = hit->position(); - Double_t xyzG[3] = {position.x(),position.y(),position.z()}; - h->SetGC(xyzG[0],xyzG[1],xyzG[2]); - Double_t xyzL[3] = {0,0,0}; - comb->MasterToLocal(xyzG,xyzL); - // if (TMath::Abs(xyzL[2]) > 0.1) continue; - Double_t uvw[3] = {h->GetU(),h->GetV(),0}; - comb->LocalToMaster(uvw,xyzG); - h->Set(xyzG,uvw); - Double_t *rot = comb->GetRotationMatrix(); - h->SetWG(rot[2],rot[5],rot[8]); - // Int_t IdH = GetIndexOfHitT(h); - Int_t IdH = fNhit - 1; - track->SetHitTId(IdH); - Double_t invpT = track->GetInvpT(); - if (TMath::Abs(invpT) < 1e-7) invpT = 1e-7; - h->SetpT(1./invpT); - h->SetMom(track->GetMomentum()); - h->SetWG(rot[2],rot[5],rot[8]); - TGeoHMatrix *rotL = (TGeoHMatrix *) RotMatrices()->FindObject(Form("WL%s",comb->GetName()+1)); - Double_t xyzLadder[3] = {0,0,0}; - if (rotL) { - - rotL->LocalToMaster(uvw,xyzLadder); - h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); - Double_t uvwP[3] = {h->GetPredU(),h->GetPredV(),0}; - rotL->LocalToMaster(uvwP,xyzLadder); - h->SetXyzL(xyzLadder); - } else { - - cout << Form("WL%s",comb->GetName()+1) << " has not been found" << endl; - h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); - h->SetXyzL(xyzLadder); - } - return h; + UInt_t B = 0, L = 0, l = 0, W = 0, H = 0; + Int_t rdo = 0; + h->SetRDO(rdo); + if (hit->detector() == kGmtId) { + StGmtHit *ht = (StGmtHit *) hit; + B = ht->getModule(); + h->SetId(B,L,l,W,H); + h->SetuvD(ht->getLocalY(), ht->getLocalX()); + h->SetuvDError(ht->getErrorLocalY(), ht->getErrorLocalX()); + h->SetSigma(ht->getSigmaY(), ht->getSigmaX()); + h->SetSigmaError(ht->getErrorSigmaY(), ht->getErrorSigmaX()); + h->SetAdc(ht->getAdcY(), ht->getAdcX()); + h->SetAdcError(ht->getErrorAdcY(), ht->getErrorAdcX()); + h->SetUsedInFit(hit->usedInFit()); + } // if (hit->detector() == kGmtId) + StThreeVectorF position = hit->position(); + Double_t xyzG[3] = {position.x(),position.y(),position.z()}; + h->SetGC(xyzG[0],xyzG[1],xyzG[2]); + Double_t xyzL[3] = {0,0,0}; + comb->MasterToLocal(xyzG,xyzL); + // if (TMath::Abs(xyzL[2]) > 0.1) continue; + Double_t uvw[3] = {h->GetU(),h->GetV(),0}; + comb->LocalToMaster(uvw,xyzG); + h->Set(xyzG,uvw); + Double_t *rot = comb->GetRotationMatrix(); + h->SetWG(rot[2],rot[5],rot[8]); + // Int_t IdH = GetIndexOfHitT(h); + Int_t IdH = fNhit - 1; + track->SetHitTId(IdH); + Double_t invpT = track->GetInvpT(); + if (TMath::Abs(invpT) < 1e-7) invpT = 1e-7; + h->SetpT(1./invpT); + h->SetMom(track->GetMomentum()); + h->SetWG(rot[2],rot[5],rot[8]); + TGeoHMatrix *rotL = (TGeoHMatrix *) RotMatrices()->FindObject(Form("WL%s",comb->GetName()+1)); + Double_t xyzLadder[3] = {0,0,0}; + if (rotL) { + rotL->LocalToMaster(uvw,xyzLadder); + h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); + Double_t uvwP[3] = {h->GetPredU(),h->GetPredV(),0}; + rotL->LocalToMaster(uvwP,xyzLadder); + h->SetXyzL(xyzLadder); + } + else { + std::cout << Form("WL%s",comb->GetName()+1) << " has not been found" << std::endl; + h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]); + h->SetXyzL(xyzLadder); + } + return h; } + //___________________ void TrackT::Print(Option_t *opt) const { - cout << "TrackT: InvpT " << fInvpT << "\tTanL " << fTanL - << "\tPhi " << fPhi << "\tRho " << fRho - << "\tNpoint " << fNpoint << "\tNsp " << fNsp << endl; - for (UInt_t i = 0; i < fNsp; i++) cout << "\t" << fIdHitT[i]; - cout << endl; + std::cout << "TrackT: InvpT " << fInvpT << "\tTanL " << fTanL + << "\tPhi " << fPhi << "\tRho " << fRho + << "\tNpoint " << fNpoint << "\tNsp " << fNsp << std::endl; + for (UInt_t i = 0; i < fNsp; i++) { + std::cout << "\t" << fIdHitT[i]; + } + std::cout << std::endl; } + //___________________ void HitT::SetId(Int_t B, Int_t L, Int_t l, Int_t W, Int_t H) { - barrel = B; layer = L; ladder = l; wafer = W; hybrid = H; + barrel = B; layer = L; ladder = l; wafer = W; hybrid = H; } + //___________________ void HitT::Print(Option_t *opt) const { - cout << "HitT: Id " << Id << "\tpT = " << pT << "\tmomentum " << pMom << endl; - TRVector glob(3,&xG); cout << "Global :" << glob << endl; - cout << "Local u/v/w " << u << "/ " << v << "/ " << w << endl; - cout << "Prediction uP/vP " << uP << "/ " << vP << "\ttuP/tvP " << tuP << "/ " << tvP << endl; + std::cout << "HitT: Id " << Id << "\tpT = " << pT << "\tmomentum " << pMom << std::endl; + TRVector glob(3,&xG); + std::cout << "Global :" << glob << std::endl; + std::cout << "Local u/v/w " << u << "/ " << v << "/ " << w << std::endl; + std::cout << "Prediction uP/vP " << uP << "/ " << vP << "\ttuP/tvP " << tuP << "/ " << tvP << std::endl; } + //___________________ void EventT::RestoreListOfRotations() { - if (fRotList) return; - if (! gDirectory) return; - fRotList = new THashList(100,0); - fRotList->SetOwner(); - TIter nextkey(gDirectory->GetListOfKeys() ); - TKey *key; - while ((key = (TKey*) nextkey())) { - TObject *obj = key->ReadObj(); - if ( obj->IsA()->InheritsFrom( "TGeoHMatrix" ) ) { - fRotList->Add(obj); + if (fRotList) return; + if (! gDirectory) return; + fRotList = new THashList(100,0); + fRotList->SetOwner(); + TIter nextkey(gDirectory->GetListOfKeys() ); + TKey *key; + while ((key = (TKey*) nextkey())) { + TObject *obj = key->ReadObj(); + if ( obj->IsA()->InheritsFrom( "TGeoHMatrix" ) ) { + fRotList->Add(obj); + } } - } } //___________________ @@ -406,162 +433,167 @@ void TBase::Loop(Int_t Nevents) { //#if 1 - struct PlotPar_t { - const Char_t *Name; - const Char_t *Title; - Int_t nx; - Int_t ny; - Double_t xmin; - Double_t xmax; - Double_t ymin; - Double_t ymax; - }; - // plots for uP - const PlotPar_t plotUP = { "uP", "track u", 320, 3, -5., 5., 0., 3. }; - // plots for u-uP - const PlotPar_t plotDu = { "Du", "Du before cut", 250, 3, -2., 2., 0., 3. }; - // plots for du & dv - const PlotPar_t plotDuv = { "Du", "Du cuts", 200, 3, -1., 1., 0., 3. }; - - TFile *fOut = new TFile(fOutFileName,"recreate"); - TString Name; - TString Title; - TString uName; - TString uTitle; - enum {NM = 8}; // no. of modules - // B - TH1F *LocPlots[NM]; - TH1F * uPlots[NM]; -#if 0 - TH1F *hpT = new TH1F( "Pt", "pt", 200, -2., 2.); - TH1F *hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); - TH1F *uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); - TH1F * xCuts[NM]; -#endif - TH1F * uPAll = new TH1F("UPall","uPall", plotUP.nx, plotUP.xmin, plotUP.xmax); - TH1F * duB[NM][2]; - TH1F * dvB[NM]; - TH1F * uCuts[NM]; - TH1F * uCut = new TH1F("Ucut","uc", plotDu.nx, plotDu.xmin, plotDu.xmax); - TH1F * vCut = new TH1F("Vcut","vc", 200, -3., 3.); - TH2F * dMin = new TH2F("DMin","vumin",100,-0.75,0.75,100,-0.75,0.75); - TH1F * vMin = new TH1F("VMin","vmin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); - TH1F * uMin = new TH1F("UMin","umin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); - //TH1F * uMinC = new TH1F("UMinC","umC", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); - memset(LocPlots,0,NM*sizeof(TH1F *)); - memset( uPlots,0,NM*sizeof(TH1F *)); - for (int M = 0; M < NM; M++) {// over gmt Modules - uName = Form("UModule%i", M+1); - uTitle = Form("du for M%i", M+1); - duB[M][0] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); - uName = Form("VModule%i", M+1); - uTitle = Form("dv for M%i", M+1); - dvB[M] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); - uName = Form("UModule%iVcut", M+1); - uTitle = Form("du for M%i after Vcut", M+1); - duB[M][1] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); - Int_t module = M; - uName = plotUP.Name; - uName += Form("M%i", module); - uTitle = Form("uP for Module %i", module); - uPlots[module] = new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); - uName = Form("%sM%i", plotDu.Name, module); - uTitle = Form("u-uP for M %i", module); - uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); - uName = Form("%sxM%i", plotDu.Name, module); - uTitle = Form("u-uP corr M %i", module); -#ifdef 0 - xCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); -#endif - } // for (int M = 0; M < NM; M++) + struct PlotPar_t { + const Char_t *Name; + const Char_t *Title; + Int_t nx; + Int_t ny; + Double_t xmin; + Double_t xmax; + Double_t ymin; + Double_t ymax; + }; + // plots for uP + const PlotPar_t plotUP = { "uP", "track u", 320, 3, -5., 5., 0., 3. }; + // plots for u-uP + const PlotPar_t plotDu = { "Du", "Du before cut", 250, 3, -2., 2., 0., 3. }; + // plots for du & dv + const PlotPar_t plotDuv = { "Du", "Du cuts", 200, 3, -1., 1., 0., 3. }; + + TFile *fOut = new TFile(fOutFileName,"recreate"); + TString Name; + TString Title; + TString uName; + TString uTitle; + enum {NM = 8}; // no. of modules + // B + TH1F *LocPlots[NM]; + TH1F * uPlots[NM]; + #if 0 + TH1F *hpT = new TH1F( "Pt", "pt", 200, -2., 2.); + TH1F *hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); + TH1F *uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); + TH1F * xCuts[NM]; + #endif + TH1F * uPAll = new TH1F("UPall","uPall", plotUP.nx, plotUP.xmin, plotUP.xmax); + TH1F * duB[NM][2]; + TH1F * dvB[NM]; + TH1F * uCuts[NM]; + TH1F * uCut = new TH1F("Ucut","uc", plotDu.nx, plotDu.xmin, plotDu.xmax); + TH1F * vCut = new TH1F("Vcut","vc", 200, -3., 3.); + TH2F * dMin = new TH2F("DMin","vumin",100,-0.75,0.75,100,-0.75,0.75); + TH1F * vMin = new TH1F("VMin","vmin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + TH1F * uMin = new TH1F("UMin","umin", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + //TH1F * uMinC = new TH1F("UMinC","umC", plotDuv.nx, plotDuv.xmin, plotDuv.xmax); + memset(LocPlots,0,NM*sizeof(TH1F *)); + memset( uPlots,0,NM*sizeof(TH1F *)); + // Loop over gmt modules + for (int M = 0; M < NM; M++) { + uName = Form("UModule%i", M+1); + uTitle = Form("du for M%i", M+1); + duB[M][0] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + uName = Form("VModule%i", M+1); + uTitle = Form("dv for M%i", M+1); + dvB[M] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + uName = Form("UModule%iVcut", M+1); + uTitle = Form("du for M%i after Vcut", M+1); + duB[M][1] = new TH1F(uName, uTitle, plotDuv.nx, plotDuv.xmin, plotDuv.xmax ); + Int_t module = M; + uName = plotUP.Name; + uName += Form("M%i", module); + uTitle = Form("uP for Module %i", module); + uPlots[module] = new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); + uName = Form("%sM%i", plotDu.Name, module); + uTitle = Form("u-uP for M %i", module); + uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + uName = Form("%sxM%i", plotDu.Name, module); + uTitle = Form("u-uP corr M %i", module); + #ifdef 0 + xCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + #endif + } // for (int M = 0; M < NM; M++) + + Long64_t nentries = fChain->GetEntriesFast(); + if (Nevents > 0 && nentries > Nevents) { + nentries = Nevents; + } - Long64_t nentries = fChain->GetEntriesFast(); - if (Nevents > 0 && nentries > Nevents) { - nentries = Nevents; - } + Long64_t nbytes = 0, nb = 0; + Int_t TreeNo = -1; + TString currentFile(""); + + // Loop over events in tree + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; - if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) { - std::cout << "Read event \t" << jentry - << " so far, switch to file " << fChain->GetCurrentFile()->GetName() - << std::endl; - std::cout << " current TreeNo: " << TreeNo - << " new TreeNo: " << fChain->GetTreeNumber() << std::endl; - TreeNo = fChain->GetTreeNumber(); - } // for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + if (! jentry%1000 || TreeNo != fChain->GetTreeNumber()) { + std::cout << "Read event \t" << jentry + << " so far, switch to file " << fChain->GetCurrentFile()->GetName() + << std::endl; + std::cout << " current TreeNo: " << TreeNo + << " new TreeNo: " << fChain->GetTreeNumber() << std::endl; + TreeNo = fChain->GetTreeNumber(); + } // for (Long64_t jentry=0; jentry 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue; - UInt_t Ntrack = fEvent->GetTracks()->GetEntriesFast(); - // int k_used[100000] = {0}; + // if (VertexZCut > 0 && TMath::Abs(fVertex[2]) > VertexZCut) continue; + UInt_t Ntrack = fEvent->GetTracks()->GetEntriesFast(); + // int k_used[100000] = {0}; - for (UInt_t trk = 0; trk < Ntrack; trk++) { - TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); - if (! track) continue; + // Loop over tracks + for (UInt_t trk = 0; trk < Ntrack; trk++) { + TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); + if (! track) continue; #if 0 - Int_t Npoints = track->GetNpoint(); - if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue; - if (UseSsd && Npoints < 1000) continue; - if (UseSvt && Npoints < 100) continue; + Int_t Npoints = track->GetNpoint(); + if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue; + if (UseSsd && Npoints < 1000) continue; + if (UseSvt && Npoints < 100) continue; #endif - Int_t Nsp = track->GetN(); - double dvmin = 1000.; - double dumin = 1000.; - int kmin{0}; - for (Int_t hit = 0; hit < Nsp; hit++) { - Int_t k = track->GetHitTId(hit) - 1; - // assert(k>=0); - HitT *hitT = (HitT *) fEvent->GetHitT(k); - if ( k < 0 ) std::cout <<" k <0:"<GetEntriesFast(); @@ -593,7 +593,7 @@ void TBase::Loop(Int_t Nevents) { } //track loop } //jentry loop (event loop) fOut->Write(); - -//#endif // #ifdef 1 + +//#endif // #if 1 } //end of Loop() From 48a7271600419bb3d1538f03d9d6e4236d794717 Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Thu, 23 May 2024 10:23:41 -0500 Subject: [PATCH 10/11] Make preprocessor directive to hide code parts --- StRoot/StGmtAligmentMaker/EventT.cxx | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/StRoot/StGmtAligmentMaker/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx index a1de19c78eb..68d66df038c 100644 --- a/StRoot/StGmtAligmentMaker/EventT.cxx +++ b/StRoot/StGmtAligmentMaker/EventT.cxx @@ -63,7 +63,7 @@ EventT::~EventT() { //_________________ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { -#if 0 +#ifdef Extra static const Int_t NoFitPointCutForGoodTrackT = 15; #endif Clear(); @@ -266,6 +266,7 @@ Int_t EventT::Build(StEvent *pEventT, Double_t pCut) { if (nTotMatch) iok = 0; return iok; } + //_________________ TrackT *EventT::AddTrackT() { // Add a new track to the list of tracks for this event. @@ -459,7 +460,7 @@ void TBase::Loop(Int_t Nevents) { // B TH1F *LocPlots[NM]; TH1F * uPlots[NM]; -#if 0 +#ifdef Extra TH1F *hpT = new TH1F( "Pt", "pt", 200, -2., 2.); TH1F *hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); TH1F *uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); @@ -498,7 +499,7 @@ void TBase::Loop(Int_t Nevents) { uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); uName = Form("%sxM%i", plotDu.Name, module); uTitle = Form("u-uP corr M %i", module); -#if 0 +#ifdef Extra xCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); #endif } // for (int M = 0; M < NM; M++) @@ -535,7 +536,7 @@ void TBase::Loop(Int_t Nevents) { for (UInt_t trk = 0; trk < Ntrack; trk++) { TrackT *track = (TrackT *)fEvent->GetTracks()->UncheckedAt(trk); if (! track) continue; -#if 0 +#ifdef Extra Int_t Npoints = track->GetNpoint(); if (minNoFitPoints > 0 && Npoints%100 < minNoFitPoints) continue; if (UseSsd && Npoints < 1000) continue; @@ -559,7 +560,7 @@ void TBase::Loop(Int_t Nevents) { Double32_t vP = hitT->GetPredV(); Double32_t du = u - uP; Double32_t dv = v - vP; -#if 0 +#ifdef Extra hpT->Fill(hitT->GetXyzPGl()); hpM->Fill(hitT->GetpMom()); if (TMath::Abs(hitT->GetpT()) < 0.2) continue; From 6728fed2a1573dd670932afb68df39c0094b2435 Mon Sep 17 00:00:00 2001 From: Grigory Nigmatkulov Date: Wed, 29 May 2024 10:15:23 -0500 Subject: [PATCH 11/11] Trying to clean up more --- StRoot/StGmtAligmentMaker/EventT.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/StRoot/StGmtAligmentMaker/EventT.cxx b/StRoot/StGmtAligmentMaker/EventT.cxx index 68d66df038c..7478859513d 100644 --- a/StRoot/StGmtAligmentMaker/EventT.cxx +++ b/StRoot/StGmtAligmentMaker/EventT.cxx @@ -465,11 +465,11 @@ void TBase::Loop(Int_t Nevents) { TH1F *hpM = new TH1F( "Ptot", "ptot", 200, 0., 5.); TH1F *uAll = new TH1F("Uall", "ua", plotUP.nx, plotUP.xmin, plotUP.xmax); TH1F * xCuts[NM]; + TH1F * uCuts[NM]; #endif TH1F * uPAll = new TH1F("UPall","uPall", plotUP.nx, plotUP.xmin, plotUP.xmax); TH1F * duB[NM][2]; TH1F * dvB[NM]; - TH1F * uCuts[NM]; TH1F * uCut = new TH1F("Ucut","uc", plotDu.nx, plotDu.xmin, plotDu.xmax); TH1F * vCut = new TH1F("Vcut","vc", 200, -3., 3.); TH2F * dMin = new TH2F("DMin","vumin",100,-0.75,0.75,100,-0.75,0.75); @@ -496,11 +496,11 @@ void TBase::Loop(Int_t Nevents) { uPlots[module] = new TH1F(uName, uTitle, plotUP.nx, plotUP.xmin, plotUP.xmax ); uName = Form("%sM%i", plotDu.Name, module); uTitle = Form("u-uP for M %i", module); - uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); uName = Form("%sxM%i", plotDu.Name, module); uTitle = Form("u-uP corr M %i", module); #ifdef Extra xCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); + uCuts[module] = new TH1F(uName, uTitle, plotDu.nx, plotDu.xmin, plotDu.xmax ); #endif } // for (int M = 0; M < NM; M++) @@ -545,7 +545,7 @@ void TBase::Loop(Int_t Nevents) { Int_t Nsp = track->GetN(); double dvmin = 1000.; double dumin = 1000.; - int kmin{0}; + //int kmin{0}; // Loop over hits for (Int_t hit = 0; hit < Nsp; hit++) { Int_t k = track->GetHitTId(hit) - 1; @@ -569,7 +569,7 @@ void TBase::Loop(Int_t Nevents) { if ( TMath::Abs(dv) < TMath::Abs(dvmin) ) { dvmin = dv; dumin= du; - kmin = k; + //kmin = k; } uPAll->Fill( uP ); uPlots[module]->Fill( uP );