StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEmcFilter.cxx
1 //******************************************************************************
2 //
3 // StEmcFilter.cxx
4 //
5 // Authors: Alexandre Suaide & Marcia Maria de Moura
6 //
7 // Initial version: 2001/12/10
8 //
9 //******************************************************************************
10 
11 #include "StEmcFilter.h"
12 #include "StMcEventTypes.hh"
13 #include "StEventTypes.h"
14 #include "StMcEvent.hh"
15 #include "StEvent.h"
16 #include <math.h>
17 #include "TFile.h"
18 #include "StTpcDedxPidAlgorithm.h"
19 #include "StEventUtilities/StuRefMult.hh"
20 #include "StEmcUtil/geometry/StEmcGeom.h"
21 #include "StEmcPosition.h"
22 #include "StarClassLibrary/SystemOfUnits.h"
23 #include "StuProbabilityPidAlgorithm.h"
24 #include "TString.h"
25 #ifndef ST_NO_NAMESPACES
26 using units::tesla;
27 #endif
28 
29 ClassImp(StEmcFilter)
30 
31 //------------------------------------------------------------------------------
38 StEmcFilter::StEmcFilter(Int_t mode):TObject()
39 {
40  mGeo[0] = StEmcGeom::getEmcGeom("bemc");
41  mGeo[1] = StEmcGeom::getEmcGeom("bprs");
42  mGeo[2] = StEmcGeom::getEmcGeom("bsmde");
43  mGeo[3] = StEmcGeom::getEmcGeom("bsmdp");
44 
45  mPion=StPionPlus::instance();
46  mProton=StProton::instance();
47  mKaon=StKaonPlus::instance();
48  mElectron=StElectron::instance();
49 
50  // event cuts
51  mPrintLog=kFALSE;
52  mEmcPresent=kTRUE;
53  mHaveVertex=kTRUE;
54  mHavePrimaries=kTRUE;
55  mZVertexCut=20;
56  mEmcETotalMin=-100000;
57  mEmcETotalMax=100000;
58  mMinMult=0;
59  mMaxMult=100000;
60  mBField=0.5;
61  mNTriggerMask=0;
62 
63  // tracks cuts
64  mDCACut=300000.0;
65  mPtCut=0.0;
66  mPtCutMax=1000.0;
67  mEtaMin=-10000.;
68  mEtaMax=10000.;
69  mFitPointsCut=0;
70  mMustProjectEmc=kTRUE;
71 
72  // dE/dX cuts
73  mdEdXScale=1.0;
74  mPointsdEdX=0;
75  mdEdXPMax=1.0;
76  mdEdXCut=0.0;
77  mdEdXNSigma=2.0;
78 
79  // V0 Vertex cuts
80  mV0Pt = 0;
81  mV0TrackProjectOnEmc = kTRUE;
82 
83  // EMC tower cuts
84  mMaxTracksPerTower = 0;
85  mEMin = -1000;
86  mEMax = 1000;
87  mPtMin = 0;
88  mPtMax = 1000;
89  mPNeighbor = kTRUE;
90  mSNeighbor = kFALSE;
91  mTNeighbor = kFALSE;
92 
93  // MC Tracks cuts
94  mOnlyHadrons = kTRUE;
95  mMcChargeType = kALL;
96  mMcMustProjectEmc = kTRUE;
97  mMcPtCut = 0.0;
98  mMcEtaMin=-10000.;
99  mMcEtaMax=10000.;
100 
101  mBemcRunning = NULL;
102  mBemcRunningOrig=NULL;
103  mBprsRunning = NULL;
104  mBsmdeRunning = NULL;
105  mBsmdpRunning = NULL;
106 
107  if(mode!=0)
108  {
109  mBemcRunning = new St_emcStatus("bemcStatus",1);
110  emcStatus_st *table=mBemcRunning->GetTable();
111  for(Int_t i=0;i<4800;i++)
112  {
113  table[0].Status[i]=0;
114  if(mode==1 && i<2400) table[0].Status[i]=1;
115  if(mode==2) table[0].Status[i]=1;
116  }
117  mBemcRunningOrig=mBemcRunning;
118  }
119 
120  mEmcPosition = new StEmcPosition();
121 
122 }
123 //------------------------------------------------------------------------------
124 StEmcFilter::~StEmcFilter()
125 {
126  if(mBemcRunningOrig) delete mBemcRunningOrig;
127  delete mEmcPosition;
128 }
129 //------------------------------------------------------------------------------
130 void StEmcFilter::calcCentrality(StEvent* event)
131 {
132  // works only for AuAu at 200 GeV
133 
134  if(fabs(mBField)>0.4) //full field
135  {
136  Int_t cent[] = {14,30,56,94,146,217,312,431,510,1000};
137  for(Int_t i=0;i<10;i++) mCentrality[i] = cent[i];
138  }
139  else // half field
140  {
141  Int_t cent[] = {14,32,59,98,149,216,302,409,474,1000};
142  for(Int_t i=0;i<10;i++) mCentrality[i] = cent[i];
143  }
144 
145  Int_t primaries=uncorrectedNumberOfPrimaries(*event);
146 
147  mCentralityBin=9;
148 
149  for(Int_t i=0;i<9;i++)
150  if(primaries < mCentrality[i]) { mCentralityBin=i; return; }
151  return;
152 }
153 //------------------------------------------------------------------------------
154 Bool_t StEmcFilter::accept(StEvent* event)
155 {
156  if(!event) return kFALSE;
157 
158  if(mNTriggerMask!=0)
159  {
160  UInt_t trigger = event->triggerMask();
161  Bool_t accept = kFALSE;
162  for(UInt_t n = 0;n<mNTriggerMask;n++)
163  if(mTriggerMask[n] == trigger) accept = kTRUE;
164  if(!accept) return kFALSE;
165 
166  }
167 
168  if(mHaveVertex)
169  {
170  StPrimaryVertex *vertex = event->primaryVertex();
171  if(!vertex) return kFALSE;
172 
173  StThreeVectorF v = vertex->position();
174  Float_t vz=v.z();
175  if(fabs(vz)>fabs(mZVertexCut))
176  {
177  cout <<"\nVertex = "<<vz<<endl<<endl;
178  return kFALSE;
179  }
180 
181  if(mHavePrimaries)
182  {
183  StSPtrVecPrimaryTrack& tracks = vertex->daughters();
184  if(tracks.size()==0) return kFALSE;
185  }
186  }
187 
188  calcCentrality(event);
189 
190  Float_t EventMultiplicity = uncorrectedNumberOfNegativePrimaries(*event);
191  if(EventMultiplicity<mMinMult || EventMultiplicity>mMaxMult) return kFALSE;
192 
193  if(mEmcPresent)
194  {
195  StEmcCollection *emc=event->emcCollection();
196  if(!emc) return kFALSE;
197 
198  Float_t emcETotal = getEmcETotal(event);
199  if(emcETotal<mEmcETotalMin || emcETotal>mEmcETotalMax) return kFALSE;
200 
201  }
202 
203  return kTRUE;
204 }
205 //------------------------------------------------------------------------------
206 Float_t StEmcFilter::getEmcETotal(StEvent * event)
207 {
208  StEmcCollection *emc=event->emcCollection();
209  if(!emc) return 0;
210 
211  Float_t ETotal = 0;
212  StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
213  StEmcDetector* emcDet = emc->detector(id) ;
214  if(!emcDet) return 0;
215  for(Int_t m=1;m<121;m++)
216  {
217  StEmcModule* module = emcDet->module(m);
218  if(module)
219  {
220  StSPtrVecEmcRawHit& Hits = module->hits();
221  if(Hits.size()>0)
222  for(Int_t k=0;k<(Int_t)Hits.size();k++)
223  ETotal+=Hits[k]->energy();
224  }
225  }
226  return ETotal;
227 }
228 //------------------------------------------------------------------------------/
229 void StEmcFilter::initEmcTowers(StEvent *event,Int_t mode)
230 {
231  for(Int_t i=0;i<NTOWER;i++)
232  {
233  mPtTower[i]=0;
234  mETower[i]=0;
235  mNTracksTower[i]=0;
236  }
237  if(!event) return;
238 
239  StEmcCollection *emc=event->emcCollection();
240  if(!emc) return;
241 
242  // filing emc hits
243  StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
244  StEmcDetector* emcDet = emc->detector(id) ;
245  if(!emcDet) return;
246  for(Int_t m=1;m<121;m++)
247  {
248  StEmcModule* module = emcDet->module(m);
249  if(module)
250  {
251  StSPtrVecEmcRawHit& Hits = module->hits();
252  if(Hits.size()>0)
253  for(Int_t k=0;k<(Int_t)Hits.size();k++)
254  {
255  Int_t m = Hits[k]->module();
256  Int_t e = Hits[k]->eta();
257  Int_t s = abs(Hits[k]->sub());
258  if(abs(m)<=120)
259  {
260  Int_t rid;
261  mGeo[0]->getId(m,e,s,rid);
262  mETower[rid-1]=Hits[k]->energy();
263  }
264  }
265  }
266  }
267 
268  //checking tracks
269  UInt_t ntracks =0;
270  if(mode==0)
271  {
272  StSPtrVecTrackNode& tracks=event->trackNodes();
273  ntracks=tracks.size();
274  }
275  else
276  {
277  if(event->l3Trigger())
278  {
279  StSPtrVecTrackNode& tracks =event->l3Trigger()->trackNodes();
280  ntracks=tracks.size();
281  }
282  }
283 
284  if(ntracks>0) for (UInt_t i = 0; i < ntracks; i++)
285  {
286 
287  StTrack* track=NULL;
288  if(mode==0) // tpt tracks
289  {
290  StSPtrVecTrackNode& tracks=event->trackNodes();
291  track=tracks[i]->track(0);
292  }
293  else //L3 tracks
294  {
295  StSPtrVecTrackNode& tracks =event->l3Trigger()->trackNodes();
296  track=tracks[i]->track(0);
297  }
298  if(track)
299  {
300  StThreeVectorD position, momentum;
301  if (mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius()))
302  {
303  Float_t eta=position.pseudoRapidity();
304  Float_t phi=position.phi();
305  Int_t m,e,s;
306  mGeo[0]->getBin(phi,eta,m,e,s);
307  if(s==-1) s=1;
308  if(m<1 || m>120) goto NEXT;
309  if(s<1 || s>2 ) goto NEXT;
310  if(e<1 || e>20 ) goto NEXT;
311  Int_t rid;
312  mGeo[0]->getId(m,e,s,rid);
313  if(getEmcStatus(1,rid)==kBAD) goto NEXT;
314  mNTracksTower[rid-1]++;
315  mPtTower[rid-1]+=track->geometry()->momentum().perp();
316  }
317  NEXT: continue;
318  }
319  }
320  return;
321 }
322 //------------------------------------------------------------------------------/
323 Bool_t StEmcFilter::accept(StTrack *track)
324 {
325  if(!track) return kFALSE;
326 
327  StThreeVectorD p = track->geometry()->momentum();
328  if(p.mag()==0) return kFALSE;
329 
330  if(p.perp()<mPtCut) return kFALSE;
331  if(p.perp()>mPtCutMax) return kFALSE;
332 
333  if(p.pseudoRapidity()<mEtaMin || p.pseudoRapidity()>mEtaMax) return kFALSE;
334  if(track->fitTraits().numberOfFitPoints()<mFitPointsCut) return kFALSE;
335  if(track->impactParameter()>mDCACut) return kFALSE;
336 
337  if(mMustProjectEmc)
338  {
339  StThreeVectorD position, momentum;
340  if (!mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius())) return kFALSE;
341 
342  //check if it hit a valid EMC tower
343  Float_t eta=position.pseudoRapidity();
344  Float_t phi=position.phi();
345  Int_t m,e,s;
346  mGeo[0]->getBin(phi,eta,m,e,s);
347  if(s==-1) s=1;
348  if(m<1 || m>120) return kFALSE;
349  if(s<1 || s>2 ) return kFALSE;
350  if(e<1 || e>20 ) return kFALSE;
351  Int_t id;
352  mGeo[0]->getId(m,e,s,id);
353  if(getEmcStatus(1,id)==kBAD) return kFALSE;
354  }
355 
356  return kTRUE;
357 
358 }
359 //------------------------------------------------------------------------------/
360 Bool_t StEmcFilter::accept(StMcEvent *event)
361 {
362  if(!event) return kFALSE;
363 
364  StMcVertex* vertex = event->primaryVertex();
365  if(!vertex) return kFALSE;
366 
367  StThreeVectorF position = vertex->position();
368  if(fabs(position.z())>mZVertexCut) return kFALSE;
369 
370  return kTRUE;
371 }
372 //------------------------------------------------------------------------------/
373 Bool_t StEmcFilter::accept(StV0Vertex *vertex)
374 {
375  if(!vertex) return kFALSE;
376  StThreeVectorF p = vertex->momentum();
377  if(p.perp()<mV0Pt) return kFALSE;
378  if(!mV0TrackProjectOnEmc) return kTRUE;
379 
380  UInt_t nd = vertex->numberOfDaughters();
381  if(nd==0) return kFALSE;
382 
383  for(UInt_t i=0;i<nd;i++)
384  {
385  StTrack* track = vertex->daughter(i);
386  if(track)
387  {
388  StThreeVectorD position, momentum;
389  if (mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius()))
390  {
391  //check if it hit a valid EMC tower
392  Float_t eta=position.pseudoRapidity();
393  Float_t phi=position.phi();
394  Int_t m,e,s;
395  mGeo[0]->getBin(phi,eta,m,e,s);
396  if(s==-1) s=1;
397  if(m>=1 && m<=120) if(s>=1 && s<=2 ) if(e>=1 && e<=20 )
398  {
399  Int_t id;
400  mGeo[0]->getId(m,e,s,id);
401  if(getEmcStatus(1,id)==kGOOD) return kTRUE;
402  }
403  }
404  }
405  }
406  return kFALSE;
407 }
408 //------------------------------------------------------------------------------/
409 Bool_t StEmcFilter::accept(StMcTrack *track)
410 {
411  if(!track) return kFALSE;
412 
413  if(mOnlyHadrons)
414  {
415  Int_t geantId = (Int_t)track->geantId();
416  if(geantId<8 || geantId>32) return kFALSE;
417  }
418 
419  if(mMcChargeType != kALL)
420  {
421  Int_t charge = (Int_t)track->particleDefinition()->charge();
422  Int_t signal = 0;
423 
424  if(charge!=0) signal = charge/abs(charge);
425 
426  if(signal == 0 && mMcChargeType != kNEUTRAL) return kFALSE;
427  if(signal == -1 && (mMcChargeType != kCHARGED && mMcChargeType != kNEGATIVE)) return kFALSE;
428  if(signal == +1 && (mMcChargeType != kCHARGED && mMcChargeType != kPOSITIVE)) return kFALSE;
429  }
430 
431  StThreeVectorD p = track->momentum();
432  if(p.mag()==0) return kFALSE;
433  if(p.perp()<mMcPtCut) return kFALSE;
434  if(p.pseudoRapidity()<mMcEtaMin || p.pseudoRapidity()>mMcEtaMax) return kFALSE;
435 
436  if(mMcMustProjectEmc)
437  {
438  // first check if there is a stop vertex
439  StMcVertex *stop=track->stopVertex();
440  if(stop) // check if the stop vertex is before EMC
441  {
442  StThreeVectorF po = stop->position();
443  Float_t stopRadius = ::sqrt(po.x()*po.x()+po.y()*po.y());
444  if(stopRadius<mGeo[0]->Radius()) return kFALSE; // track stopped before EMC. Not accepted
445  }
446 
447  StThreeVectorD position, momentum;
448  if (!mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius())) return kFALSE;
449 
450  //check if it hit a valid EMC tower
451  Float_t eta=position.pseudoRapidity();
452  Float_t phi=position.phi();
453  Int_t m,e,s;
454  mGeo[0]->getBin(phi,eta,m,e,s);
455  if(s==-1) s=1;
456  if(m<1 || m>120) return kFALSE;
457  if(s<1 || s>2 ) return kFALSE;
458  if(e<1 || e>20 ) return kFALSE;
459  Int_t id;
460  mGeo[0]->getId(m,e,s,id);
461  if(getEmcStatus(1,id)==kBAD) return kFALSE;
462  }
463  return kTRUE;
464 }
465 //------------------------------------------------------------------------------/
466 Bool_t StEmcFilter::accept(Int_t rid)
467 {
468  if(rid<1 || rid>4800) return kFALSE;
469  if(getEmcStatus(1,rid)!=kGOOD) return kFALSE;
470  if(mNTracksTower[rid-1]>mMaxTracksPerTower) return kFALSE;
471  if(mETower[rid-1]<mEMin || mETower[rid-1]>mEMax) return kFALSE;
472  if(mPtTower[rid-1]<mPtMin || mPtTower[rid-1]>mPtMax) return kFALSE;
473  Int_t size = 0;
474  if(mTNeighbor) { size = 3; goto DOIT;}
475  if(mSNeighbor) { size = 2; goto DOIT;}
476  if(mPNeighbor) size = 1;
477  DOIT:
478  if(size==0) return kTRUE;
479  Float_t eta,phi;
480  mGeo[0]->getEtaPhi(rid,eta,phi);
481  for(Int_t i = -size; i<= size; i++)
482  for(Int_t j = -size; j<= size; j++)
483  {
484  if(!(i==0 && j==0))
485  {
486  Int_t rid1 = mEmcPosition->getNextTowerId(eta,phi,i,j);
487  if(rid1>0)
488  if(mNTracksTower[rid1-1]>0) return kFALSE;
489  }
490  }
491  return kTRUE;
492 }
493 //------------------------------------------------------------------------------/
494 Bool_t StEmcFilter::accept(Int_t m, Int_t e, Int_t s)
495 {
496  Int_t rid;
497  mGeo[0]->getId(m,e,s,rid);
498  return accept(rid);
499 }
500 //------------------------------------------------------------------------------/
501 Int_t StEmcFilter::getNTracksTower(Int_t rid)
502 {
503  if(rid<1 || rid>4800) return 0;
504  return mNTracksTower[rid-1];
505 }
506 //------------------------------------------------------------------------------/
507 Int_t StEmcFilter::getNTracksTower(Int_t m, Int_t e, Int_t s)
508 {
509  Int_t rid;
510  mGeo[0]->getId(m,e,s,rid);
511  return getNTracksTower(rid);
512 }
513 //------------------------------------------------------------------------------/
514 Float_t StEmcFilter::getPtTower(Int_t rid)
515 {
516  if(rid<1 || rid>4800) return 0;
517  return mPtTower[rid-1];
518 }
519 //------------------------------------------------------------------------------/
520 Float_t StEmcFilter::getPtTower(Int_t m, Int_t e, Int_t s)
521 {
522  Int_t rid;
523  mGeo[0]->getId(m,e,s,rid);
524  return getPtTower(rid);
525 }
526 //------------------------------------------------------------------------------
532 Bool_t StEmcFilter::getTrackId(StTrack *track,Float_t& mass,Int_t& id)
533 {
534  Float_t nSigma[4];
535  Int_t order[4];
536  Float_t dEdX;
537  Int_t npt;
538  return getTrackId(track,npt,dEdX,mass,id,order,nSigma);
539 }
540 //------------------------------------------------------------------------------
553 Bool_t StEmcFilter::getTrackId(StTrack *track,Int_t& nPoints,Float_t& dEdX,Float_t& mass,Int_t& id,Int_t *idOrder,Float_t *nSigmaFinal)
554 {
555  // dE/dx
556  id=8;
557  mass=mPion->mass();
558  if(!track) return kFALSE;
559  if(!track->geometry()) return kFALSE;
560  Int_t charge=track->geometry()->charge();
561  if(charge<0) id=9;
562 
563  Double_t momentum = fabs(track->geometry()->momentum().mag());
564  if(momentum==0) return kFALSE;
565  StPtrVecTrackPidTraits traits = track->pidTraits(kTpcId);
566  UInt_t size = traits.size();
567  if(size==0) return kFALSE;
568 
569  Float_t m[4];
570  m[0] = mPion->mass();
571  m[1] = mProton->mass();
572  m[2] = mKaon->mass();
573  m[3] = mElectron->mass();
574  Int_t kk[4]={0,0,0,1};
575 
576  if (size>0)
577  {
578  StDedxPidTraits* pid=NULL;
579  for (UInt_t i = 0; i < traits.size(); i++)
580  {
581  pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
582  if (pid) if(pid->method() == kTruncatedMeanId) break;
583  }
584  if(!pid) return kFALSE;
585 
586  dEdX = (Float_t)pid->mean();
587  if(dEdX==0) return kFALSE;
588  Double_t npt = (Double_t)pid->numberOfPoints();
589  nPoints=(Int_t) npt;
590  if(nPoints==0) return kFALSE;
591  Double_t dedx_expected;
592  Double_t dedx_resolution = (Double_t)pid->errorOnMean();
593  if(dedx_resolution<=0) dedx_resolution=npt > 0 ? 0.45/::sqrt(npt) : 1000.;
594  Double_t z;
595  Float_t nSigma[4],nSigmaTmp[4];
596  Float_t length = (Float_t)pid->length();
597  if(length<=0) length = 60.;
598  for(Int_t i=0;i<4;i++)
599  {
600  //dedx_expected=mBB(momentum/m[i])*mdEdXScale;
601  dedx_expected = 1.0e-6*mBB.Sirrf(momentum/m[i],length,kk[i])*mdEdXScale;
602  z = ::log((Double_t)dEdX/dedx_expected);
603  nSigmaTmp[i]=(Float_t) z/dedx_resolution;
604  nSigma[i]=fabs(nSigmaTmp[i]) ;
605  }
606 
607  Float_t SigmaOrder[4];
608  Int_t type[4];
609  for(Int_t i=0;i<4;i++)
610  {
611  SigmaOrder[i]=999999;
612  type[i]=0;
613  for(Int_t k=0;k<4;k++)
614  {
615  if(nSigma[k]<SigmaOrder[i]) {SigmaOrder[i]=nSigma[k]; nSigmaFinal[i]=nSigmaTmp[k]; type[i]=k;}
616  }
617  nSigma[type[i]]=9999;
618  }
619 
620  for(Int_t i=0;i<4;i++)
621  {
622  if(type[i]==0 && charge>0) idOrder[i] = 8;
623  if(type[i]==0 && charge<0) idOrder[i] = 9;
624  if(type[i]==1 && charge>0) idOrder[i] = 14;
625  if(type[i]==1 && charge<0) idOrder[i] = 15;
626  if(type[i]==2 && charge>0) idOrder[i] = 11;
627  if(type[i]==2 && charge<0) idOrder[i] = 12;
628  if(type[i]==3 && charge>0) idOrder[i] = 2;
629  if(type[i]==3 && charge<0) idOrder[i] = 3;
630  }
631 
632  if(momentum>mdEdXPMax) return kFALSE;
633  if(npt<mPointsdEdX) return kFALSE;
634  if(dEdX<mdEdXCut) return kFALSE;
635 
636  if(SigmaOrder[0]<=mdEdXNSigma)
637  {
638  if(type[0]==0 && charge>0) {mass=mPion->mass(); id = 8; return kTRUE;}
639  if(type[0]==0 && charge<0) {mass=mPion->mass(); id = 9; return kTRUE;}
640 
641  if(type[0]==1 && charge>0) {mass=mProton->mass(); id = 14; return kTRUE;}
642  if(type[0]==1 && charge<0) {mass=mProton->mass(); id = 15; return kTRUE;}
643 
644  if(type[0]==2 && charge>0) {mass=mKaon->mass(); id = 11; return kTRUE;}
645  if(type[0]==2 && charge<0) {mass=mKaon->mass(); id = 12; return kTRUE;}
646 
647  if(type[0]==3 && charge>0) {mass=mElectron->mass();id = 2; return kTRUE;}
648  if(type[0]==3 && charge<0) {mass=mElectron->mass();id = 3; return kTRUE;}
649  }
650  //cout <<"final id = "<<id<<endl;
651  return kFALSE;
652  }
653  return kFALSE;
654 }
655 //------------------------------------------------------------------------------
662 EmcStatus StEmcFilter::getEmcStatus(Int_t det, Int_t id)
663 {
664  switch (det)
665  {
666  case 1: // bemc
667  if(mBemcRunning)
668  {
669  emcStatus_st *st=mBemcRunning->GetTable();
670  if(id<1 || id >4800) return kBAD;
671  else
672  {
673  Int_t status = (Int_t) st[0].Status[id-1];
674  if(status == 1) return kGOOD; else return kBAD;
675  return kOTHER;
676  }
677  }
678  else
679  {
680  if ((id >= 1 && id <= 340) || (id >= 1861 && id <= 2340) ) return kGOOD;
681  else return kBAD;
682  }
683  break;
684  case 2: // bprs
685  if(mBprsRunning)
686  {
687  emcStatus_st *st=mBprsRunning->GetTable();
688  if(id<1 || id >4800) return kBAD;
689  else
690  {
691  Int_t status = (Int_t) st[0].Status[id-1];
692  if(status == 1) return kGOOD;else return kBAD;
693  return kOTHER;
694  }
695  }
696  else return kBAD;
697  break;
698  case 3: // bsmde
699  if(mBsmdeRunning)
700  {
701  smdStatus_st *st=mBsmdeRunning->GetTable();
702  if(id<1 || id >18000) return kBAD;
703  else
704  {
705  Int_t status = (Int_t) st[0].Status[id-1];
706  if(status == 1) return kGOOD;else return kBAD;
707  return kOTHER;
708  }
709  }
710  else return kBAD;
711  break;
712  case 4: // bsmdp
713  if(mBsmdpRunning)
714  {
715  smdStatus_st *st=mBsmdpRunning->GetTable();
716  if(id<1 || id >18000) return kBAD;
717  else
718  {
719  Int_t status = (Int_t) st[0].Status[id-1];
720  if(status == 1) return kGOOD;else return kBAD;
721  return kOTHER;
722  }
723  }
724  else return kBAD;
725  break;
726  }
727  return kBAD;
728 }
729 //------------------------------------------------------------------------------
735 Float_t StEmcFilter::getFraction(Int_t det,Int_t etabin,Int_t side)
736 {
737  if(det<1 || det>4) return 0.;
738 
739  Float_t ntowers=0;
740  Float_t TotalTowers=60*mGeo[det-1]->NSub(); //60 modules x sub/module
741 
742  if(etabin<1 || etabin>mGeo[det-1]->NEta()) return 0.;
743 
744  Int_t mi=1,mf=60;
745  if(side==1) {mi=61; mf=120;}
746 
747  for(Int_t m=mi;m<=mf;m++)
748  for(Int_t s=1;s<=mGeo[det-1]->NSub();s++)
749  {
750  Int_t id;
751  mGeo[det-1]->getId(m,etabin,s,id);
752  if(getEmcStatus(det,id)==kGOOD) ntowers++;
753  }
754 
755  return ntowers/TotalTowers;
756 }
757 //------------------------------------------------------------------------------
761 Float_t StEmcFilter::getWestFraction(Int_t det)
762 {
763  if(det<1 || det>4) return 0.;
764  Int_t nEta = mGeo[det-1]->NEta();
765 
766  Float_t phiFr=0;
767  for(Int_t i=1;i<=nEta;i++)
768  {
769  phiFr+=getFraction(det,i,0);
770  }
771  return phiFr/(Float_t)nEta;
772 }
773 //------------------------------------------------------------------------------
777 Float_t StEmcFilter::getEastFraction(Int_t det)
778 {
779  if(det<1 || det>4) return 0.;
780  Int_t nEta = mGeo[det-1]->NEta();
781 
782  Float_t phiFr=0;
783  for(Int_t i=1;i<=nEta;i++)
784  {
785  phiFr+=getFraction(det,i,1);
786  }
787  return phiFr/(Float_t)nEta;
788 }
789 //------------------------------------------------------------------------------
793 Float_t StEmcFilter::getTotalFraction(Int_t det)
794 {
795  if(det<1 || det>4) return 0.;
796  Float_t WPhiFr=getWestFraction(det);
797  Float_t EPhiFr=getEastFraction(det);
798 
799  return (WPhiFr+EPhiFr)/2.;
800 }
801 //------------------------------------------------------------------------------
802 void StEmcFilter::printCuts()
803 {
804  TString TF[]={"kFALSE","kTRUE"};
805  TString CH[]={"kNEGATIVE", "kNEUTRAL", "kPOSITIVE", "kCHARGED", "kALL"};
806  cout <<"EMC Filter ---------------------------------------------------\n\n";
807  cout <<" 1. Event Cuts \n";
808  cout <<" EmcPresent = "<<TF[(Int_t)mEmcPresent].Data()<<endl;
809  cout <<" HaveVertex = "<<TF[(Int_t)mHaveVertex].Data()<<endl;
810  cout <<" HavePrimaries = "<<TF[(Int_t)mHavePrimaries].Data()<<endl;
811  cout <<" ZVertexCut = "<<mZVertexCut<<" cm\n";
812  cout <<" EmcETotalMin = "<<mEmcETotalMin<<" GeV\n";
813  cout <<" EmcETotalMax = "<<mEmcETotalMax<<" GeV\n";
814  cout <<" MinMult = "<<mMinMult<<endl;
815  cout <<" MaxMult = "<<mMaxMult<<endl;
816  cout <<" BField = "<<mBField<<" Tesla \n";
817  cout <<endl;
818  cout <<" 2. Tracks cuts\n";
819  cout <<" DCACut = "<<mDCACut<<" cm\n";
820  cout <<" PtCut = "<<mPtCut<<" GeV/c\n";
821  cout <<" EtaMin = "<<mEtaMin<<endl;
822  cout <<" EtaMax = "<<mEtaMax<<endl;
823  cout <<" FitPointsCut = "<<mFitPointsCut<<endl;
824  cout <<" MusProjectEmc = "<<TF[(Int_t)mMustProjectEmc].Data()<<endl;
825  cout <<endl;
826  cout <<" 3. dE/dX cuts for StTrack Id\n";
827  cout <<" dEdXScale = "<<mdEdXScale<<endl;
828  cout <<" PointsdEdX = "<<mPointsdEdX<<endl;
829  cout <<" dEdXPMax = "<<mdEdXPMax<<" GeV/c\n";
830  cout <<" dEdXCut = "<<mdEdXCut<<" GeV/cm\n";
831  cout <<" dEdXNSigma = "<<mdEdXNSigma<<endl;
832  cout <<endl;
833  cout <<" 4. MC Tracks cuts\n";
834  cout <<" OnlyHadrons = "<<TF[(Int_t)mOnlyHadrons].Data()<<endl;
835  cout <<" McChargeType = "<<CH[(Int_t)mMcChargeType+1].Data()<<endl;
836  cout <<" McMustProjectEmc = "<<TF[(Int_t)mMcMustProjectEmc].Data()<<endl;
837  cout <<" McEtaMin = "<<mMcEtaMin<<endl;
838  cout <<" McEtaMac = "<<mMcEtaMax<<endl;
839  cout <<"--------------------------------------------------------------\n";
840  return;
841 }
842 
843 
844 
845 
StEmcPosition::trackOnEmc
Bool_t trackOnEmc(StThreeVectorD *position, StThreeVectorD *momentum, const StTrack *const track, Double_t magField, Double_t emcRadius=225.405) const
Track projection utility.
Definition: StEmcPosition.cxx:270
StEmcFilter::getEmcStatus
EmcStatus getEmcStatus(Int_t, Int_t)
Return EMC status (kGOOD, kBAD, kOTHER) for a given detector and bin.
Definition: StEmcFilter.cxx:662
StEmcFilter::getFraction
Float_t getFraction(Int_t, Int_t, Int_t=0)
Return fraction of EMC live on a given detector and eta bin.
Definition: StEmcFilter.cxx:735
StEmcFilter::printCuts
void printCuts()
Print filter parameters.
Definition: StEmcFilter.cxx:802
StEmcFilter::getEmcETotal
Float_t getEmcETotal(StEvent *)
Return total energy on EMC.
Definition: StEmcFilter.cxx:206
StEmcFilter::getNTracksTower
Int_t getNTracksTower(Int_t)
Returns number of tracks pointing to the tower with given id.
Definition: StEmcFilter.cxx:501
StEmcFilter::initEmcTowers
void initEmcTowers(StEvent *, Int_t=0)
Use this function before using accept(), getNTracksTower() and getPtTower() methods for towers...
Definition: StEmcFilter.cxx:229
StEmcFilter::~StEmcFilter
~StEmcFilter()
StEmcFilter destructor.
Definition: StEmcFilter.cxx:124
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StEmcPosition::getNextTowerId
Int_t getNextTowerId(Float_t eta, Float_t phi, Int_t nTowersdEta, Int_t nTowersdPhi) const
Return neighbor tower id&#39;s.
Definition: StEmcPosition.cxx:356
StEmcFilter::accept
Bool_t accept(StEvent *)
Returns kTRUE/kFALSE if StEvent is accepted.
Definition: StEmcFilter.cxx:154
StEmcFilter::getTrackId
Bool_t getTrackId(StTrack *, Float_t &, Int_t &)
Return track id based on dE/dX.
Definition: StEmcFilter.cxx:532
StEmcFilter::getEastFraction
Float_t getEastFraction(Int_t=1)
Return fraction of EMC live on east side of STAR.
Definition: StEmcFilter.cxx:777
StEmcGeom::getBin
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Definition: StEmcGeom.h:321
StEmcFilter::getPtTower
Float_t getPtTower(Int_t)
Returns the total pt from projected tracks in tower with id.
Definition: StEmcFilter.cxx:514
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
StEmcFilter::getWestFraction
Float_t getWestFraction(Int_t=1)
Return fraction of EMC live on west side of STAR.
Definition: StEmcFilter.cxx:761
StEmcFilter::getTotalFraction
Float_t getTotalFraction(Int_t=1)
Return fraction of EMC live on STAR.
Definition: StEmcFilter.cxx:793
Generated by   doxygen 1.8.5