StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcFgtCorrelatorA.cxx
1 /***************************************************************************
2  *
3  * $Id: StEEmcFgtCorrelatorA.cxx,v 1.2 2012/05/09 21:11:58 sgliske Exp $
4  * Author: S. Gliske, May 2012
5  *
6  ***************************************************************************
7  *
8  * Description: see header
9  *
10  ***************************************************************************
11  *
12  * $Log: StEEmcFgtCorrelatorA.cxx,v $
13  * Revision 1.2 2012/05/09 21:11:58 sgliske
14  * updates
15  *
16  * Revision 1.1 2012/05/09 17:26:26 sgliske
17  * creation
18  *
19  *
20  **************************************************************************/
21 
22 #include "StEEmcFgtCorrelatorA.h"
23 
24 #include "StEEmcRawMapMaker.h"
25 #include "StThreeVectorF.hh"
26 
27 #include "StRoot/StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
28 #include "StRoot/StEEmcPool/StEEmcPointMap/StEEmcPointMap.h"
29 #include "StRoot/StEvent/StEvent.h"
30 #include "StRoot/StEvent/StPrimaryVertex.h"
31 #include "StRoot/StMuDSTMaker/COMMON/StMuDst.h"
32 #include "StRoot/StMuDSTMaker/COMMON/StMuEvent.h"
33 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
34 
35 
36 #define DEBUG
37 #define CLUSTER_WIDTH 4
38 #define CLUSTER_SIZE 9
39 #define FAIL_FLAG -1234
40 
41 // constructors
42 StEEmcFgtCorrelatorA::StEEmcFgtCorrelatorA( const Char_t* name, const Char_t* rawMapMkrName ) :
43  StMaker( name ), mInputType(1), mInputName("MuDst"), mEEmcRawMapMkr(0),
44  mMipMin(5), mMipMax(50), mSigThres(3) {
45 
46  mEEmcRawMapMkr = static_cast< StEEmcRawMapMaker* >( GetMaker( rawMapMkrName ) );
47  assert( mEEmcRawMapMkr );
48  assert( CLUSTER_SIZE == CLUSTER_WIDTH*2+1 ); // fix your consts
49  assert( FAIL_FLAG < 0 );
50 };
51 
52 // deconstructor
53 StEEmcFgtCorrelatorA::~StEEmcFgtCorrelatorA(){ /* */ };
54 
55 // init
56 Int_t StEEmcFgtCorrelatorA::Init(){
57  return kStOk;
58 };
59 
60 Int_t StEEmcFgtCorrelatorA::Make(){
61  Int_t ierr = loadVertex();
62 
63  const StEEmcRawMap& eemcMap = mEEmcRawMapMkr->getMap( 4 );
64 
65  std::vector< Int_t > seedIdxVec;
66  std::vector< Float_t > mipClusPosVec[12]; // still in u/v indexing, but with fraction position between strips
67  std::vector< TVector3 > mipPosVec; // 3D space points
68 
69  StEEmcRawMap::const_iterator mapIter;
70  if( !eemcMap.empty() ){
71 
72  // first, find all with signal in range
73  for( mapIter = eemcMap.begin(); mapIter != eemcMap.end(); ++mapIter ){
74  Int_t idx = mapIter->first;
75  const StEEmcRawMapData& data = mapIter->second;
76 
77  Float_t resp = data.rawAdc - data.ped;
78  if( data.fail || data.stat )
79  resp = 0;
80  if( resp && data.pedSigma )
81  resp /= data.pedSigma;
82 
83  if( resp < mMipMax && resp > mMipMin )
84  seedIdxVec.push_back( idx );
85  };
86  };
87 
88 #ifdef DEBUG2
89  LOG_INFO << "zzz " << GetEventNumber() << " found " << seedIdxVec.size() << " MIP seeds out of " << eemcMap.size() << " strips" << endm;
90 #endif
91 
92  Int_t nMipClus = 0;
93  Int_t nFailBothAdj = 0, nFailFail = 0, nFailSideHigh = 0, nFailNonIso = 0;
94 
95  if( !seedIdxVec.empty() ){
96  std::vector< Int_t >::const_iterator seedIter;
97 
98  Float_t clusterShape[CLUSTER_SIZE];
99 
100  for( seedIter = seedIdxVec.begin(); seedIter != seedIdxVec.end(); ++seedIter ){
101 
102  // clear cluster shape
103  for( Float_t *p = clusterShape; p != &clusterShape[CLUSTER_SIZE]; ++p )
104  *p = 0;
105 
106  // compute min and max index for the specific u/v and sector
107  Int_t first = ((*seedIter)/kEEmcNumStrips)*kEEmcNumStrips;
108  Int_t last = first + kEEmcNumStrips;
109 
110  // compute range of the cluster
111  Int_t offset = *seedIter - CLUSTER_WIDTH;
112  Int_t low = std::max( first, *seedIter - CLUSTER_WIDTH );
113  Int_t high = std::min( last, *seedIter + CLUSTER_WIDTH + 1 );
114 
115  // copy response values to the cluster shape array
116  for( Int_t i = offset; i<offset+CLUSTER_SIZE; ++i ){
117  if( i >= low && i < high ){
118  mapIter = eemcMap.find( i );
119 
120  if( mapIter != eemcMap.end() ){
121  const StEEmcRawMapData& data = mapIter->second;
122 
123  Float_t resp = data.rawAdc - data.ped;
124  if( resp && data.pedSigma )
125  resp /= data.pedSigma;
126  if( data.fail || data.stat )
127  resp = FAIL_FLAG; // magic value to flag bad strip
128 
129  clusterShape[i-offset] = resp;
130  };
131  };
132  };
133 
134 #ifdef DEBUG_CLUS_SHAPE
135  LOG_INFO << "zzz cluster shape: ";
136  for( Int_t i=0; i<CLUSTER_SIZE; ++i )
137  LOG_INFO << clusterShape[i] << ' ';
138  LOG_INFO << endm;
139 #endif
140 
141  // now check the cluster shape to see if it qualifies as a MIP
142 
143  // ensure exactly one of the adjacent have signal
144  Bool_t pass = (( clusterShape[CLUSTER_WIDTH-1] > mSigThres ) ^ ( clusterShape[CLUSTER_WIDTH+1] > mSigThres ));
145  if( !pass )
146  ++nFailBothAdj;
147 
148 
149  //cout << "zzz " << clusterShape[CLUSTER_WIDTH-1] << ' ' << clusterShape[CLUSTER_WIDTH+1] << endl;
150 
151  // ensure adjacent are good strips and that they are less than the seed
152  if( pass && (clusterShape[CLUSTER_WIDTH+1] == FAIL_FLAG || clusterShape[CLUSTER_WIDTH-1] == FAIL_FLAG) ){
153  ++nFailFail;
154  pass = 0;
155  };
156  if( pass && ( clusterShape[CLUSTER_WIDTH+1] > clusterShape[CLUSTER_WIDTH] || clusterShape[CLUSTER_WIDTH-1] > clusterShape[CLUSTER_WIDTH] ) ){
157  ++nFailSideHigh;
158  pass = 0;
159  };
160 
161  Bool_t passOlder = pass;
162  // ensure the farther strips have no signal
163  for( Int_t i = 0; i<CLUSTER_WIDTH-1 && pass; ++i )
164  if( clusterShape[i] > mSigThres )
165  pass = 0;
166  for( Int_t i = CLUSTER_WIDTH+2; i<CLUSTER_SIZE && pass; ++i )
167  if( clusterShape[i] > mSigThres )
168  pass = 0;
169 
170  if( passOlder && !pass )
171  ++nFailNonIso;
172 
173  if( pass ){
174  // determine position better
175  Float_t posA = *seedIter;
176  Float_t wA = clusterShape[CLUSTER_WIDTH];
177  Float_t posB = (( clusterShape[CLUSTER_WIDTH+1] > clusterShape[CLUSTER_WIDTH-1] ) ? *seedIter+1 : *seedIter-1 );
178  Float_t wB = std::max( clusterShape[CLUSTER_WIDTH+1], clusterShape[CLUSTER_WIDTH-1] );
179 
180  ++nMipClus;
181  Int_t sec = (*seedIter)/kEEmcNumStrips/2;
182  mipClusPosVec[sec].push_back(( posA*wA + posB*wB ) / ( wA + wB ));
183  };
184  };
185  };
186 
187 #ifdef DEBUG2
188  LOG_INFO << "zzz FOUND " << nMipClus << " possible MIP clusters, "
189  << "failed " << nFailBothAdj << ' ' << nFailFail << ' ' << nFailSideHigh << ' ' << nFailNonIso << endm;
190 #endif
191 
192  if( nMipClus ){
193  std::vector< Float_t >::iterator mipClusPosIter; // still in u/v indexing, but with fraction position between strips
194 
195  // limit to only sectors 6-12, since only have half coverage of FGT in 2012
196  for( Int_t sec = 6; sec<kEEmcNumSectors; ++sec ){
197 #ifdef DEBUG3
198  LOG_INFO << "zzz \t sector " << sec << " has " << mipClusPosVec[sec].size() << " MIP clusters " << endm;
199 
200  for( UInt_t j = 0; j<mipClusPosVec[sec].size(); ++j ){
201  Float_t idx = mipClusPosVec[sec][j];
202  Int_t idx2 = idx;
203  Int_t strip = idx2 % kEEmcNumStrips;
204  Bool_t isV = (idx2 / kEEmcNumStrips)%2;
205  Int_t sec = (idx2 / kEEmcNumStrips)/2;
206 
207  LOG_INFO << "zzz \t\t " << idx << ' ' << sec << ' ' << (isV ? 'v' : 'u') << ' ' << strip << ' ' << strip+idx-idx2 << endm;
208  };
209 #endif
210 
211  if( mipClusPosVec[sec].size() == 2 ){
212  Float_t idx1 = mipClusPosVec[sec][0];
213  Float_t idx2 = mipClusPosVec[sec][1];
214  Bool_t isV1 = ((Int_t)idx1/kEEmcNumStrips)%2;
215  Bool_t isV2 = ((Int_t)idx2/kEEmcNumStrips)%2;
216 
217  if( isV1 ^ isV2 ){
218  Float_t idxU = ( isV1 ? idx2 : idx1 );
219  Float_t idxV = ( isV1 ? idx1 : idx2 );
220  Float_t x = 0, y = 0;
221 
222  StEEmcPointMap_t::instance().convertStripUVtoXY( sec, idxU, idxV, x, y );
223  mipPosVec.push_back( TVector3( x, y, kEEmcZSMD ) );
224  };
225  };
226  };
227  };
228 
229 #ifdef DEBUG
230  if( !mipPosVec.empty() ){
231  LOG_INFO << "zzz EVENT " << GetEventNumber() << " FOUND " << mipPosVec.size() << " MIP points" << endm;
232  };
233 #endif
234 
235  if( !mipPosVec.empty() ){
236 #ifdef DEBUG
237  LOG_INFO << "zzz -> vertex " << mVertex.X() << ' ' << mVertex.Y() << ' ' << mVertex.Z() << endm;
238 #endif
239 
240  for( UInt_t j = 0; j < mipPosVec.size(); ++j ){
241  TVector3& smdPt = mipPosVec[j];
242  TVector3 delta = smdPt - mVertex;
243 
244 #ifdef DEBUG
245  LOG_INFO << "zzz -> ESMD point " << smdPt.X() << ' ' << smdPt.Y() << ' ' << smdPt.Z() << endm;
246 #endif
247  for( Int_t disc = 0; disc<kFgtNumDiscs; ++disc ){
248  Float_t z = StFgtGeom::getDiscZ( disc );
249  Float_t alpha = ( z - mVertex.Z() ) / ( smdPt.Z() - mVertex.Z() );
250  TVector3 pos = alpha*delta + mVertex;
251 
252 #ifdef DEBUG
253  LOG_INFO << "zzz ----> disc " << disc+1 << " line passes through " << pos.X() << ' ' << pos.Y() << ' ' << pos.Z() << endm;
254 #endif
255  };
256  };
257  };
258 
259  return ierr;
260 };
261 
262 void StEEmcFgtCorrelatorA::Clear( Option_t *opt ){
263  // ???
264 };
265 
266 Int_t StEEmcFgtCorrelatorA::setInput( const Char_t *name, Int_t type ){
267  Int_t ierr = kStOk;
268 
269  if( type == 0 || type == 1){
270  mInputType = type;
271  mInputName = name;
272  } else {
273  LOG_ERROR << "Invalid input type" << endm;
274  ierr = kStFatal;
275  };
276  return ierr;
277 };
278 
279 Int_t StEEmcFgtCorrelatorA::loadVertex(){
280  Int_t ierr = kStErr;
281 
282  if( mInputType ){
283  // load from MuDst
284  const StMuDst* muDst = (const StMuDst*)GetInputDS( mInputName.data() );
285  assert( muDst );
286 
287  StMuEvent *event = muDst->event();
288  if( event ){
289  const StThreeVectorF& v = event->primaryVertexPosition();
290 
291  if( v.z() > -300 && v.z() < 300 ){
292  ierr = kStOk;
293  mVertex.SetXYZ( v.x(), v.y(), v.z() );
294  };
295  };
296  } else {
297  // load from StEvent
298  const StEvent *event = (const StEvent*)GetInputDS( mInputName.data() );
299  assert( event );
300  const StPrimaryVertex* vertex = event->primaryVertex();
301 
302  if( vertex ){
303  const StThreeVectorF& v = vertex->position();
304 
305  if( v.z() > -300 && v.z() < 300 ){
306  ierr = kStOk;
307  mVertex.SetXYZ( v.x(), v.y(), v.z() );
308  };
309  };
310  };
311 
312  return ierr;
313 };
314 
315 ClassImp(StEEmcFgtCorrelatorA);
StEEmcFgtCorrelatorA::Clear
virtual void Clear(Option_t *opt="")
User defined functions.
Definition: StEEmcFgtCorrelatorA.cxx:262
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5