00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "StEmcFilter.h"
00012 #include "StMcEventTypes.hh"
00013 #include "StEventTypes.h"
00014 #include "StMcEvent.hh"
00015 #include "StEvent.h"
00016 #include <math.h>
00017 #include "TFile.h"
00018 #include "StTpcDedxPidAlgorithm.h"
00019 #include "StEventUtilities/StuRefMult.hh"
00020 #include "StEmcUtil/geometry/StEmcGeom.h"
00021 #include "StEmcPosition.h"
00022 #include "StarClassLibrary/SystemOfUnits.h"
00023 #include "StuProbabilityPidAlgorithm.h"
00024 #include "TString.h"
00025 #ifndef ST_NO_NAMESPACES
00026 using units::tesla;
00027 #endif
00028
00029 ClassImp(StEmcFilter)
00030
00031
00038 StEmcFilter::StEmcFilter(Int_t mode):TObject()
00039 {
00040 mGeo[0] = StEmcGeom::getEmcGeom("bemc");
00041 mGeo[1] = StEmcGeom::getEmcGeom("bprs");
00042 mGeo[2] = StEmcGeom::getEmcGeom("bsmde");
00043 mGeo[3] = StEmcGeom::getEmcGeom("bsmdp");
00044
00045 mPion=StPionPlus::instance();
00046 mProton=StProton::instance();
00047 mKaon=StKaonPlus::instance();
00048 mElectron=StElectron::instance();
00049
00050
00051 mPrintLog=kFALSE;
00052 mEmcPresent=kTRUE;
00053 mHaveVertex=kTRUE;
00054 mHavePrimaries=kTRUE;
00055 mZVertexCut=20;
00056 mEmcETotalMin=-100000;
00057 mEmcETotalMax=100000;
00058 mMinMult=0;
00059 mMaxMult=100000;
00060 mBField=0.5;
00061 mNTriggerMask=0;
00062
00063
00064 mDCACut=300000.0;
00065 mPtCut=0.0;
00066 mPtCutMax=1000.0;
00067 mEtaMin=-10000.;
00068 mEtaMax=10000.;
00069 mFitPointsCut=0;
00070 mMustProjectEmc=kTRUE;
00071
00072
00073 mdEdXScale=1.0;
00074 mPointsdEdX=0;
00075 mdEdXPMax=1.0;
00076 mdEdXCut=0.0;
00077 mdEdXNSigma=2.0;
00078
00079
00080 mV0Pt = 0;
00081 mV0TrackProjectOnEmc = kTRUE;
00082
00083
00084 mMaxTracksPerTower = 0;
00085 mEMin = -1000;
00086 mEMax = 1000;
00087 mPtMin = 0;
00088 mPtMax = 1000;
00089 mPNeighbor = kTRUE;
00090 mSNeighbor = kFALSE;
00091 mTNeighbor = kFALSE;
00092
00093
00094 mOnlyHadrons = kTRUE;
00095 mMcChargeType = kALL;
00096 mMcMustProjectEmc = kTRUE;
00097 mMcPtCut = 0.0;
00098 mMcEtaMin=-10000.;
00099 mMcEtaMax=10000.;
00100
00101 mBemcRunning = NULL;
00102 mBemcRunningOrig=NULL;
00103 mBprsRunning = NULL;
00104 mBsmdeRunning = NULL;
00105 mBsmdpRunning = NULL;
00106
00107 if(mode!=0)
00108 {
00109 mBemcRunning = new St_emcStatus("bemcStatus",1);
00110 emcStatus_st *table=mBemcRunning->GetTable();
00111 for(Int_t i=0;i<4800;i++)
00112 {
00113 table[0].Status[i]=0;
00114 if(mode==1 && i<2400) table[0].Status[i]=1;
00115 if(mode==2) table[0].Status[i]=1;
00116 }
00117 mBemcRunningOrig=mBemcRunning;
00118 }
00119
00120 mEmcPosition = new StEmcPosition();
00121
00122 }
00123
00124 StEmcFilter::~StEmcFilter()
00125 {
00126 if(mBemcRunningOrig) delete mBemcRunningOrig;
00127 delete mEmcPosition;
00128 }
00129
00130 void StEmcFilter::calcCentrality(StEvent* event)
00131 {
00132
00133
00134 if(fabs(mBField)>0.4)
00135 {
00136 Int_t cent[] = {14,30,56,94,146,217,312,431,510,1000};
00137 for(Int_t i=0;i<10;i++) mCentrality[i] = cent[i];
00138 }
00139 else
00140 {
00141 Int_t cent[] = {14,32,59,98,149,216,302,409,474,1000};
00142 for(Int_t i=0;i<10;i++) mCentrality[i] = cent[i];
00143 }
00144
00145 Int_t primaries=uncorrectedNumberOfPrimaries(*event);
00146
00147 mCentralityBin=9;
00148
00149 for(Int_t i=0;i<9;i++)
00150 if(primaries < mCentrality[i]) { mCentralityBin=i; return; }
00151 return;
00152 }
00153
00154 Bool_t StEmcFilter::accept(StEvent* event)
00155 {
00156 if(!event) return kFALSE;
00157
00158 if(mNTriggerMask!=0)
00159 {
00160 UInt_t trigger = event->triggerMask();
00161 Bool_t accept = kFALSE;
00162 for(UInt_t n = 0;n<mNTriggerMask;n++)
00163 if(mTriggerMask[n] == trigger) accept = kTRUE;
00164 if(!accept) return kFALSE;
00165
00166 }
00167
00168 if(mHaveVertex)
00169 {
00170 StPrimaryVertex *vertex = event->primaryVertex();
00171 if(!vertex) return kFALSE;
00172
00173 StThreeVectorF v = vertex->position();
00174 Float_t vz=v.z();
00175 if(fabs(vz)>fabs(mZVertexCut))
00176 {
00177 cout <<"\nVertex = "<<vz<<endl<<endl;
00178 return kFALSE;
00179 }
00180
00181 if(mHavePrimaries)
00182 {
00183 StSPtrVecPrimaryTrack& tracks = vertex->daughters();
00184 if(tracks.size()==0) return kFALSE;
00185 }
00186 }
00187
00188 calcCentrality(event);
00189
00190 Float_t EventMultiplicity = uncorrectedNumberOfNegativePrimaries(*event);
00191 if(EventMultiplicity<mMinMult || EventMultiplicity>mMaxMult) return kFALSE;
00192
00193 if(mEmcPresent)
00194 {
00195 StEmcCollection *emc=event->emcCollection();
00196 if(!emc) return kFALSE;
00197
00198 Float_t emcETotal = getEmcETotal(event);
00199 if(emcETotal<mEmcETotalMin || emcETotal>mEmcETotalMax) return kFALSE;
00200
00201 }
00202
00203 return kTRUE;
00204 }
00205
00206 Float_t StEmcFilter::getEmcETotal(StEvent * event)
00207 {
00208 StEmcCollection *emc=event->emcCollection();
00209 if(!emc) return 0;
00210
00211 Float_t ETotal = 0;
00212 StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
00213 StEmcDetector* emcDet = emc->detector(id) ;
00214 if(!emcDet) return 0;
00215 for(Int_t m=1;m<121;m++)
00216 {
00217 StEmcModule* module = emcDet->module(m);
00218 if(module)
00219 {
00220 StSPtrVecEmcRawHit& Hits = module->hits();
00221 if(Hits.size()>0)
00222 for(Int_t k=0;k<(Int_t)Hits.size();k++)
00223 ETotal+=Hits[k]->energy();
00224 }
00225 }
00226 return ETotal;
00227 }
00228
00229 void StEmcFilter::initEmcTowers(StEvent *event,Int_t mode)
00230 {
00231 for(Int_t i=0;i<NTOWER;i++)
00232 {
00233 mPtTower[i]=0;
00234 mETower[i]=0;
00235 mNTracksTower[i]=0;
00236 }
00237 if(!event) return;
00238
00239 StEmcCollection *emc=event->emcCollection();
00240 if(!emc) return;
00241
00242
00243 StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
00244 StEmcDetector* emcDet = emc->detector(id) ;
00245 if(!emcDet) return;
00246 for(Int_t m=1;m<121;m++)
00247 {
00248 StEmcModule* module = emcDet->module(m);
00249 if(module)
00250 {
00251 StSPtrVecEmcRawHit& Hits = module->hits();
00252 if(Hits.size()>0)
00253 for(Int_t k=0;k<(Int_t)Hits.size();k++)
00254 {
00255 Int_t m = Hits[k]->module();
00256 Int_t e = Hits[k]->eta();
00257 Int_t s = abs(Hits[k]->sub());
00258 if(abs(m)<=120)
00259 {
00260 Int_t rid;
00261 mGeo[0]->getId(m,e,s,rid);
00262 mETower[rid-1]=Hits[k]->energy();
00263 }
00264 }
00265 }
00266 }
00267
00268
00269 UInt_t ntracks =0;
00270 if(mode==0)
00271 {
00272 StSPtrVecTrackNode& tracks=event->trackNodes();
00273 ntracks=tracks.size();
00274 }
00275 else
00276 {
00277 if(event->l3Trigger())
00278 {
00279 StSPtrVecTrackNode& tracks =event->l3Trigger()->trackNodes();
00280 ntracks=tracks.size();
00281 }
00282 }
00283
00284 if(ntracks>0) for (UInt_t i = 0; i < ntracks; i++)
00285 {
00286
00287 StTrack* track=NULL;
00288 if(mode==0)
00289 {
00290 StSPtrVecTrackNode& tracks=event->trackNodes();
00291 track=tracks[i]->track(0);
00292 }
00293 else
00294 {
00295 StSPtrVecTrackNode& tracks =event->l3Trigger()->trackNodes();
00296 track=tracks[i]->track(0);
00297 }
00298 if(track)
00299 {
00300 StThreeVectorD position, momentum;
00301 if (mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius()))
00302 {
00303 Float_t eta=position.pseudoRapidity();
00304 Float_t phi=position.phi();
00305 Int_t m,e,s;
00306 mGeo[0]->getBin(phi,eta,m,e,s);
00307 if(s==-1) s=1;
00308 if(m<1 || m>120) goto NEXT;
00309 if(s<1 || s>2 ) goto NEXT;
00310 if(e<1 || e>20 ) goto NEXT;
00311 Int_t rid;
00312 mGeo[0]->getId(m,e,s,rid);
00313 if(getEmcStatus(1,rid)==kBAD) goto NEXT;
00314 mNTracksTower[rid-1]++;
00315 mPtTower[rid-1]+=track->geometry()->momentum().perp();
00316 }
00317 NEXT: continue;
00318 }
00319 }
00320 return;
00321 }
00322
00323 Bool_t StEmcFilter::accept(StTrack *track)
00324 {
00325 if(!track) return kFALSE;
00326
00327 StThreeVectorD p = track->geometry()->momentum();
00328 if(p.mag()==0) return kFALSE;
00329
00330 if(p.perp()<mPtCut) return kFALSE;
00331 if(p.perp()>mPtCutMax) return kFALSE;
00332
00333 if(p.pseudoRapidity()<mEtaMin || p.pseudoRapidity()>mEtaMax) return kFALSE;
00334 if(track->fitTraits().numberOfFitPoints()<mFitPointsCut) return kFALSE;
00335 if(track->impactParameter()>mDCACut) return kFALSE;
00336
00337 if(mMustProjectEmc)
00338 {
00339 StThreeVectorD position, momentum;
00340 if (!mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius())) return kFALSE;
00341
00342
00343 Float_t eta=position.pseudoRapidity();
00344 Float_t phi=position.phi();
00345 Int_t m,e,s;
00346 mGeo[0]->getBin(phi,eta,m,e,s);
00347 if(s==-1) s=1;
00348 if(m<1 || m>120) return kFALSE;
00349 if(s<1 || s>2 ) return kFALSE;
00350 if(e<1 || e>20 ) return kFALSE;
00351 Int_t id;
00352 mGeo[0]->getId(m,e,s,id);
00353 if(getEmcStatus(1,id)==kBAD) return kFALSE;
00354 }
00355
00356 return kTRUE;
00357
00358 }
00359
00360 Bool_t StEmcFilter::accept(StMcEvent *event)
00361 {
00362 if(!event) return kFALSE;
00363
00364 StMcVertex* vertex = event->primaryVertex();
00365 if(!vertex) return kFALSE;
00366
00367 StThreeVectorF position = vertex->position();
00368 if(fabs(position.z())>mZVertexCut) return kFALSE;
00369
00370 return kTRUE;
00371 }
00372
00373 Bool_t StEmcFilter::accept(StV0Vertex *vertex)
00374 {
00375 if(!vertex) return kFALSE;
00376 StThreeVectorF p = vertex->momentum();
00377 if(p.perp()<mV0Pt) return kFALSE;
00378 if(!mV0TrackProjectOnEmc) return kTRUE;
00379
00380 UInt_t nd = vertex->numberOfDaughters();
00381 if(nd==0) return kFALSE;
00382
00383 for(UInt_t i=0;i<nd;i++)
00384 {
00385 StTrack* track = vertex->daughter(i);
00386 if(track)
00387 {
00388 StThreeVectorD position, momentum;
00389 if (mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius()))
00390 {
00391
00392 Float_t eta=position.pseudoRapidity();
00393 Float_t phi=position.phi();
00394 Int_t m,e,s;
00395 mGeo[0]->getBin(phi,eta,m,e,s);
00396 if(s==-1) s=1;
00397 if(m>=1 && m<=120) if(s>=1 && s<=2 ) if(e>=1 && e<=20 )
00398 {
00399 Int_t id;
00400 mGeo[0]->getId(m,e,s,id);
00401 if(getEmcStatus(1,id)==kGOOD) return kTRUE;
00402 }
00403 }
00404 }
00405 }
00406 return kFALSE;
00407 }
00408
00409 Bool_t StEmcFilter::accept(StMcTrack *track)
00410 {
00411 if(!track) return kFALSE;
00412
00413 if(mOnlyHadrons)
00414 {
00415 Int_t geantId = (Int_t)track->geantId();
00416 if(geantId<8 || geantId>32) return kFALSE;
00417 }
00418
00419 if(mMcChargeType != kALL)
00420 {
00421 Int_t charge = (Int_t)track->particleDefinition()->charge();
00422 Int_t signal = 0;
00423
00424 if(charge!=0) signal = charge/abs(charge);
00425
00426 if(signal == 0 && mMcChargeType != kNEUTRAL) return kFALSE;
00427 if(signal == -1 && (mMcChargeType != kCHARGED && mMcChargeType != kNEGATIVE)) return kFALSE;
00428 if(signal == +1 && (mMcChargeType != kCHARGED && mMcChargeType != kPOSITIVE)) return kFALSE;
00429 }
00430
00431 StThreeVectorD p = track->momentum();
00432 if(p.mag()==0) return kFALSE;
00433 if(p.perp()<mMcPtCut) return kFALSE;
00434 if(p.pseudoRapidity()<mMcEtaMin || p.pseudoRapidity()>mMcEtaMax) return kFALSE;
00435
00436 if(mMcMustProjectEmc)
00437 {
00438
00439 StMcVertex *stop=track->stopVertex();
00440 if(stop)
00441 {
00442 StThreeVectorF po = stop->position();
00443 Float_t stopRadius = ::sqrt(po.x()*po.x()+po.y()*po.y());
00444 if(stopRadius<mGeo[0]->Radius()) return kFALSE;
00445 }
00446
00447 StThreeVectorD position, momentum;
00448 if (!mEmcPosition->trackOnEmc(&position, &momentum, track, mBField,mGeo[0]->Radius())) return kFALSE;
00449
00450
00451 Float_t eta=position.pseudoRapidity();
00452 Float_t phi=position.phi();
00453 Int_t m,e,s;
00454 mGeo[0]->getBin(phi,eta,m,e,s);
00455 if(s==-1) s=1;
00456 if(m<1 || m>120) return kFALSE;
00457 if(s<1 || s>2 ) return kFALSE;
00458 if(e<1 || e>20 ) return kFALSE;
00459 Int_t id;
00460 mGeo[0]->getId(m,e,s,id);
00461 if(getEmcStatus(1,id)==kBAD) return kFALSE;
00462 }
00463 return kTRUE;
00464 }
00465
00466 Bool_t StEmcFilter::accept(Int_t rid)
00467 {
00468 if(rid<1 || rid>4800) return kFALSE;
00469 if(getEmcStatus(1,rid)!=kGOOD) return kFALSE;
00470 if(mNTracksTower[rid-1]>mMaxTracksPerTower) return kFALSE;
00471 if(mETower[rid-1]<mEMin || mETower[rid-1]>mEMax) return kFALSE;
00472 if(mPtTower[rid-1]<mPtMin || mPtTower[rid-1]>mPtMax) return kFALSE;
00473 Int_t size = 0;
00474 if(mTNeighbor) { size = 3; goto DOIT;}
00475 if(mSNeighbor) { size = 2; goto DOIT;}
00476 if(mPNeighbor) size = 1;
00477 DOIT:
00478 if(size==0) return kTRUE;
00479 Float_t eta,phi;
00480 mGeo[0]->getEtaPhi(rid,eta,phi);
00481 for(Int_t i = -size; i<= size; i++)
00482 for(Int_t j = -size; j<= size; j++)
00483 {
00484 if(!(i==0 && j==0))
00485 {
00486 Int_t rid1 = mEmcPosition->getNextTowerId(eta,phi,i,j);
00487 if(rid1>0)
00488 if(mNTracksTower[rid1-1]>0) return kFALSE;
00489 }
00490 }
00491 return kTRUE;
00492 }
00493
00494 Bool_t StEmcFilter::accept(Int_t m, Int_t e, Int_t s)
00495 {
00496 Int_t rid;
00497 mGeo[0]->getId(m,e,s,rid);
00498 return accept(rid);
00499 }
00500
00501 Int_t StEmcFilter::getNTracksTower(Int_t rid)
00502 {
00503 if(rid<1 || rid>4800) return 0;
00504 return mNTracksTower[rid-1];
00505 }
00506
00507 Int_t StEmcFilter::getNTracksTower(Int_t m, Int_t e, Int_t s)
00508 {
00509 Int_t rid;
00510 mGeo[0]->getId(m,e,s,rid);
00511 return getNTracksTower(rid);
00512 }
00513
00514 Float_t StEmcFilter::getPtTower(Int_t rid)
00515 {
00516 if(rid<1 || rid>4800) return 0;
00517 return mPtTower[rid-1];
00518 }
00519
00520 Float_t StEmcFilter::getPtTower(Int_t m, Int_t e, Int_t s)
00521 {
00522 Int_t rid;
00523 mGeo[0]->getId(m,e,s,rid);
00524 return getPtTower(rid);
00525 }
00526
00532 Bool_t StEmcFilter::getTrackId(StTrack *track,Float_t& mass,Int_t& id)
00533 {
00534 Float_t nSigma[4];
00535 Int_t order[4];
00536 Float_t dEdX;
00537 Int_t npt;
00538 return getTrackId(track,npt,dEdX,mass,id,order,nSigma);
00539 }
00540
00553 Bool_t StEmcFilter::getTrackId(StTrack *track,Int_t& nPoints,Float_t& dEdX,Float_t& mass,Int_t& id,Int_t *idOrder,Float_t *nSigmaFinal)
00554 {
00555
00556 id=8;
00557 mass=mPion->mass();
00558 if(!track) return kFALSE;
00559 if(!track->geometry()) return kFALSE;
00560 Int_t charge=track->geometry()->charge();
00561 if(charge<0) id=9;
00562
00563 Double_t momentum = fabs(track->geometry()->momentum().mag());
00564 if(momentum==0) return kFALSE;
00565 StPtrVecTrackPidTraits traits = track->pidTraits(kTpcId);
00566 UInt_t size = traits.size();
00567 if(size==0) return kFALSE;
00568
00569 Float_t m[4];
00570 m[0] = mPion->mass();
00571 m[1] = mProton->mass();
00572 m[2] = mKaon->mass();
00573 m[3] = mElectron->mass();
00574 Int_t kk[4]={0,0,0,1};
00575
00576 if (size>0)
00577 {
00578 StDedxPidTraits* pid=NULL;
00579 for (UInt_t i = 0; i < traits.size(); i++)
00580 {
00581 pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
00582 if (pid) if(pid->method() == kTruncatedMeanId) break;
00583 }
00584 if(!pid) return kFALSE;
00585
00586 dEdX = (Float_t)pid->mean();
00587 if(dEdX==0) return kFALSE;
00588 Double_t npt = (Double_t)pid->numberOfPoints();
00589 nPoints=(Int_t) npt;
00590 if(nPoints==0) return kFALSE;
00591 Double_t dedx_expected;
00592 Double_t dedx_resolution = (Double_t)pid->errorOnMean();
00593 if(dedx_resolution<=0) dedx_resolution=npt > 0 ? 0.45/::sqrt(npt) : 1000.;
00594 Double_t z;
00595 Float_t nSigma[4],nSigmaTmp[4];
00596 Float_t length = (Float_t)pid->length();
00597 if(length<=0) length = 60.;
00598 for(Int_t i=0;i<4;i++)
00599 {
00600
00601 dedx_expected = 1.0e-6*mBB.Sirrf(momentum/m[i],length,kk[i])*mdEdXScale;
00602 z = ::log((Double_t)dEdX/dedx_expected);
00603 nSigmaTmp[i]=(Float_t) z/dedx_resolution;
00604 nSigma[i]=fabs(nSigmaTmp[i]) ;
00605 }
00606
00607 Float_t SigmaOrder[4];
00608 Int_t type[4];
00609 for(Int_t i=0;i<4;i++)
00610 {
00611 SigmaOrder[i]=999999;
00612 type[i]=0;
00613 for(Int_t k=0;k<4;k++)
00614 {
00615 if(nSigma[k]<SigmaOrder[i]) {SigmaOrder[i]=nSigma[k]; nSigmaFinal[i]=nSigmaTmp[k]; type[i]=k;}
00616 }
00617 nSigma[type[i]]=9999;
00618 }
00619
00620 for(Int_t i=0;i<4;i++)
00621 {
00622 if(type[i]==0 && charge>0) idOrder[i] = 8;
00623 if(type[i]==0 && charge<0) idOrder[i] = 9;
00624 if(type[i]==1 && charge>0) idOrder[i] = 14;
00625 if(type[i]==1 && charge<0) idOrder[i] = 15;
00626 if(type[i]==2 && charge>0) idOrder[i] = 11;
00627 if(type[i]==2 && charge<0) idOrder[i] = 12;
00628 if(type[i]==3 && charge>0) idOrder[i] = 2;
00629 if(type[i]==3 && charge<0) idOrder[i] = 3;
00630 }
00631
00632 if(momentum>mdEdXPMax) return kFALSE;
00633 if(npt<mPointsdEdX) return kFALSE;
00634 if(dEdX<mdEdXCut) return kFALSE;
00635
00636 if(SigmaOrder[0]<=mdEdXNSigma)
00637 {
00638 if(type[0]==0 && charge>0) {mass=mPion->mass(); id = 8; return kTRUE;}
00639 if(type[0]==0 && charge<0) {mass=mPion->mass(); id = 9; return kTRUE;}
00640
00641 if(type[0]==1 && charge>0) {mass=mProton->mass(); id = 14; return kTRUE;}
00642 if(type[0]==1 && charge<0) {mass=mProton->mass(); id = 15; return kTRUE;}
00643
00644 if(type[0]==2 && charge>0) {mass=mKaon->mass(); id = 11; return kTRUE;}
00645 if(type[0]==2 && charge<0) {mass=mKaon->mass(); id = 12; return kTRUE;}
00646
00647 if(type[0]==3 && charge>0) {mass=mElectron->mass();id = 2; return kTRUE;}
00648 if(type[0]==3 && charge<0) {mass=mElectron->mass();id = 3; return kTRUE;}
00649 }
00650
00651 return kFALSE;
00652 }
00653 return kFALSE;
00654 }
00655
00662 EmcStatus StEmcFilter::getEmcStatus(Int_t det, Int_t id)
00663 {
00664 switch (det)
00665 {
00666 case 1:
00667 if(mBemcRunning)
00668 {
00669 emcStatus_st *st=mBemcRunning->GetTable();
00670 if(id<1 || id >4800) return kBAD;
00671 else
00672 {
00673 Int_t status = (Int_t) st[0].Status[id-1];
00674 if(status == 1) return kGOOD; else return kBAD;
00675 return kOTHER;
00676 }
00677 }
00678 else
00679 {
00680 if ((id >= 1 && id <= 340) || (id >= 1861 && id <= 2340) ) return kGOOD;
00681 else return kBAD;
00682 }
00683 break;
00684 case 2:
00685 if(mBprsRunning)
00686 {
00687 emcStatus_st *st=mBprsRunning->GetTable();
00688 if(id<1 || id >4800) return kBAD;
00689 else
00690 {
00691 Int_t status = (Int_t) st[0].Status[id-1];
00692 if(status == 1) return kGOOD;else return kBAD;
00693 return kOTHER;
00694 }
00695 }
00696 else return kBAD;
00697 break;
00698 case 3:
00699 if(mBsmdeRunning)
00700 {
00701 smdStatus_st *st=mBsmdeRunning->GetTable();
00702 if(id<1 || id >18000) return kBAD;
00703 else
00704 {
00705 Int_t status = (Int_t) st[0].Status[id-1];
00706 if(status == 1) return kGOOD;else return kBAD;
00707 return kOTHER;
00708 }
00709 }
00710 else return kBAD;
00711 break;
00712 case 4:
00713 if(mBsmdpRunning)
00714 {
00715 smdStatus_st *st=mBsmdpRunning->GetTable();
00716 if(id<1 || id >18000) return kBAD;
00717 else
00718 {
00719 Int_t status = (Int_t) st[0].Status[id-1];
00720 if(status == 1) return kGOOD;else return kBAD;
00721 return kOTHER;
00722 }
00723 }
00724 else return kBAD;
00725 break;
00726 }
00727 return kBAD;
00728 }
00729
00735 Float_t StEmcFilter::getFraction(Int_t det,Int_t etabin,Int_t side)
00736 {
00737 if(det<1 || det>4) return 0.;
00738
00739 Float_t ntowers=0;
00740 Float_t TotalTowers=60*mGeo[det-1]->NSub();
00741
00742 if(etabin<1 || etabin>mGeo[det-1]->NEta()) return 0.;
00743
00744 Int_t mi=1,mf=60;
00745 if(side==1) {mi=61; mf=120;}
00746
00747 for(Int_t m=mi;m<=mf;m++)
00748 for(Int_t s=1;s<=mGeo[det-1]->NSub();s++)
00749 {
00750 Int_t id;
00751 mGeo[det-1]->getId(m,etabin,s,id);
00752 if(getEmcStatus(det,id)==kGOOD) ntowers++;
00753 }
00754
00755 return ntowers/TotalTowers;
00756 }
00757
00761 Float_t StEmcFilter::getWestFraction(Int_t det)
00762 {
00763 if(det<1 || det>4) return 0.;
00764 Int_t nEta = mGeo[det-1]->NEta();
00765
00766 Float_t phiFr=0;
00767 for(Int_t i=1;i<=nEta;i++)
00768 {
00769 phiFr+=getFraction(det,i,0);
00770 }
00771 return phiFr/(Float_t)nEta;
00772 }
00773
00777 Float_t StEmcFilter::getEastFraction(Int_t det)
00778 {
00779 if(det<1 || det>4) return 0.;
00780 Int_t nEta = mGeo[det-1]->NEta();
00781
00782 Float_t phiFr=0;
00783 for(Int_t i=1;i<=nEta;i++)
00784 {
00785 phiFr+=getFraction(det,i,1);
00786 }
00787 return phiFr/(Float_t)nEta;
00788 }
00789
00793 Float_t StEmcFilter::getTotalFraction(Int_t det)
00794 {
00795 if(det<1 || det>4) return 0.;
00796 Float_t WPhiFr=getWestFraction(det);
00797 Float_t EPhiFr=getEastFraction(det);
00798
00799 return (WPhiFr+EPhiFr)/2.;
00800 }
00801
00802 void StEmcFilter::printCuts()
00803 {
00804 TString TF[]={"kFALSE","kTRUE"};
00805 TString CH[]={"kNEGATIVE", "kNEUTRAL", "kPOSITIVE", "kCHARGED", "kALL"};
00806 cout <<"EMC Filter ---------------------------------------------------\n\n";
00807 cout <<" 1. Event Cuts \n";
00808 cout <<" EmcPresent = "<<TF[(Int_t)mEmcPresent].Data()<<endl;
00809 cout <<" HaveVertex = "<<TF[(Int_t)mHaveVertex].Data()<<endl;
00810 cout <<" HavePrimaries = "<<TF[(Int_t)mHavePrimaries].Data()<<endl;
00811 cout <<" ZVertexCut = "<<mZVertexCut<<" cm\n";
00812 cout <<" EmcETotalMin = "<<mEmcETotalMin<<" GeV\n";
00813 cout <<" EmcETotalMax = "<<mEmcETotalMax<<" GeV\n";
00814 cout <<" MinMult = "<<mMinMult<<endl;
00815 cout <<" MaxMult = "<<mMaxMult<<endl;
00816 cout <<" BField = "<<mBField<<" Tesla \n";
00817 cout <<endl;
00818 cout <<" 2. Tracks cuts\n";
00819 cout <<" DCACut = "<<mDCACut<<" cm\n";
00820 cout <<" PtCut = "<<mPtCut<<" GeV/c\n";
00821 cout <<" EtaMin = "<<mEtaMin<<endl;
00822 cout <<" EtaMax = "<<mEtaMax<<endl;
00823 cout <<" FitPointsCut = "<<mFitPointsCut<<endl;
00824 cout <<" MusProjectEmc = "<<TF[(Int_t)mMustProjectEmc].Data()<<endl;
00825 cout <<endl;
00826 cout <<" 3. dE/dX cuts for StTrack Id\n";
00827 cout <<" dEdXScale = "<<mdEdXScale<<endl;
00828 cout <<" PointsdEdX = "<<mPointsdEdX<<endl;
00829 cout <<" dEdXPMax = "<<mdEdXPMax<<" GeV/c\n";
00830 cout <<" dEdXCut = "<<mdEdXCut<<" GeV/cm\n";
00831 cout <<" dEdXNSigma = "<<mdEdXNSigma<<endl;
00832 cout <<endl;
00833 cout <<" 4. MC Tracks cuts\n";
00834 cout <<" OnlyHadrons = "<<TF[(Int_t)mOnlyHadrons].Data()<<endl;
00835 cout <<" McChargeType = "<<CH[(Int_t)mMcChargeType+1].Data()<<endl;
00836 cout <<" McMustProjectEmc = "<<TF[(Int_t)mMcMustProjectEmc].Data()<<endl;
00837 cout <<" McEtaMin = "<<mMcEtaMin<<endl;
00838 cout <<" McEtaMac = "<<mMcEtaMax<<endl;
00839 cout <<"--------------------------------------------------------------\n";
00840 return;
00841 }
00842
00843
00844
00845
