STDATAFILTERMAKER: StMtdEvtFilterMaker/StMtdEvtFilterMaker.cxx Source File

00001 #include <iostream>
00002 #include <math.h>
00003 #include <vector>
00004 #include <stdlib.h>
00005 #include <bitset>
00006 
00007 #include "TH1F.h"
00008 
00009 #include "StMessMgr.h"
00010 
00011 #include "StEventTypes.h"
00012 #include "StThreeVectorF.hh"
00013 #include "StEvent.h"
00014 #include "StTrack.h"
00015 #include "StTrackGeometry.h"
00016 #include "StDcaGeometry.h"
00017 #include "StTpcDedxPidAlgorithm.h"
00018 #include "StDedxPidTraits.h"
00019 #include "StTrackPidTraits.h"
00020 
00021 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00022 #include "StMuDSTMaker/COMMON/StMuDst.h"
00023 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00024 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00025 #include "StMuDSTMaker/COMMON/StMuMtdHit.h"
00026 #include "StMuDSTMaker/COMMON/StMuMtdPidTraits.h"
00027 
00028 
00029 #include "StMtdEvtFilterMaker.h"
00030 #include "tables/St_trgOfflineFilter_Table.h"
00031 #include "tables/St_mtdEventFilterCuts_Table.h"
00032 #include "tables/St_MtdTrackFilterTag_Table.h"
00033 
00034 #include "tables/St_HighPtTag_Table.h"
00035 
00036 ClassImp(StMtdEvtFilterMaker)
00037 
00038 //_____________________________________________________________________________
00039 StMtdEvtFilterMaker::StMtdEvtFilterMaker(const Char_t *name) : StMaker(name)
00040 {
00041   // default constructor
00042 
00043   mStEvent               = NULL;
00044   mMuDst                 = NULL;
00045   mIsJpsiEvent           = kFALSE;
00046   mIsDiMuon              = kFALSE;
00047   mIsDiMuonOnly          = kFALSE;
00048 
00049      
00050    mMinTrkPtAll           = 1.0;
00051    mMinTrkPtLead          = 0;
00052    mMinNHitsFit           = 15;
00053    mMinNHitsDedx          = 10;
00054    mMinFitHitsFraction    = 0.52;
00055    mMaxDca                = 1e4;
00056    mMinNsigmaPi           = -1e4;
00057    mMaxNsigmaPi           = 1e4;
00058    mMaxDeltaZ             = 1e4;
00059    nMinMuonCandidates     = 2;
00060 }
00061  
00062 //_____________________________________________________________________________
00063 StMtdEvtFilterMaker::~StMtdEvtFilterMaker()
00064 {
00065   // default destructor
00066 }
00067 
00068 //_____________________________________________________________________________
00069 Int_t StMtdEvtFilterMaker::Init()
00070 {
00071   SetFlavor("MtdEvtFilter","trgOfflineFilter");
00072 
00073   // book histograms
00074   bookHistos();
00075 
00076   return kStOK;
00077 }
00078 
00079 //_____________________________________________________________________________
00080 Int_t StMtdEvtFilterMaker::InitRun(const Int_t runNumber)
00081 {
00082   LOG_INFO << "Retrieving trigger ID from database ..." << endm;
00083   mTriggerIDs.clear();
00084   mOtherTrigIDs.clear();
00085 
00086   St_trgOfflineFilter* flaggedTrgs =
00087     static_cast<St_trgOfflineFilter *>(GetDataBase("Calibrations/trg/trgOfflineFilter"));
00088   if (!flaggedTrgs) 
00089     {
00090       LOG_ERROR << "Could not find Calibrations/trg/trgOfflineFilter in database" << endm;
00091       return kStErr;
00092     }
00093 
00094   trgOfflineFilter_st* flaggedTrg = flaggedTrgs->GetTable();
00095   for (long j = 0; j < flaggedTrgs->GetNRows(); j++, flaggedTrg++) 
00096     {
00097       int trigid = flaggedTrg->trigid;
00098       if(trigid==0) continue;
00099       else if (trigid<1e9) mTriggerIDs.push_back(trigid);
00100       else mOtherTrigIDs.push_back(trigid-1e9);
00101     }
00102 
00103   if(Debug())
00104     {
00105       for(unsigned int i=0; i<mTriggerIDs.size(); i++){
00106         LOG_INFO << "Di-muon trigger ID: " << mTriggerIDs[i] << endm; }
00107 
00108       for(unsigned int i=0; i<mOtherTrigIDs.size(); i++){
00109         LOG_INFO << "Other MTD trigger ID: " << mOtherTrigIDs[i] << endm; }
00110     }
00111 
00112   LOG_INFO << "Retrieving event filtering cuts from database ..." << endm;
00113 
00114   St_mtdEventFilterCuts *mtdEventFilterCuts = 
00115     static_cast<St_mtdEventFilterCuts*>(GetDataBase("Calibrations/mtd/mtdEventFilterCuts"));
00116   if(!mtdEventFilterCuts)
00117     {
00118       LOG_ERROR << "No mtdEventFilterCuts table found in database" << endm;
00119       return kStErr;
00120     }
00121   mtdEventFilterCuts_st *table = static_cast<mtdEventFilterCuts_st*>(mtdEventFilterCuts->GetTable());
00122     
00123   mMinTrkPtAll           = table->minTrkPtAll;
00124   mMinTrkPtLead          = table->minTrkPtLead;
00125   mMinNHitsFit           = (int)table->minNHitsFit;
00126   mMinNHitsDedx          = (int)table->minNHitsDedx;
00127   mMinFitHitsFraction    = table->minFitHitsFrac;
00128   mMaxDca                = table->maxDca;
00129   mMinNsigmaPi           = table->minNsigmaPi;
00130   mMaxNsigmaPi           = table->maxNsigmaPi;
00131   mMaxDeltaZ             = table->maxDeltaZ;
00132   nMinMuonCandidates     = (int)table->minNMuons;
00133   if(Debug())
00134     {
00135       LOG_INFO << "minTrkPtAll    = " << mMinTrkPtAll        << endm;
00136       LOG_INFO << "minTrkPtLead   = " << mMinTrkPtLead       << endm;
00137       LOG_INFO << "minNHitsFit    = " << mMinNHitsFit        << endm;
00138       LOG_INFO << "minNHitsDedx   = " << mMinNHitsDedx       << endm;
00139       LOG_INFO << "minFitHitsFrac = " << mMinFitHitsFraction << endm;
00140       LOG_INFO << "maxDca         = " << mMaxDca             << endm;
00141       LOG_INFO << "minNsigmaPi    = " << mMinNsigmaPi        << endm;
00142       LOG_INFO << "maxNsigmaPi    = " << mMaxNsigmaPi        << endm;
00143       LOG_INFO << "maxDeltaZ      = " << mMaxDeltaZ          << endm;
00144       LOG_INFO << "minNMuons      = " << nMinMuonCandidates  << endm;
00145     }
00146 
00147   
00148   return kStOK;
00149 }
00150 
00151 //_____________________________________________________________________________
00152 Int_t StMtdEvtFilterMaker::Make()
00153 {
00154 
00155   mIsJpsiEvent = kFALSE;
00156   mIsDiMuon = kFALSE;
00157   mIsDiMuonOnly = kFALSE;
00158 
00159 
00160   // Check the availability of input data
00161   Int_t iret = -1;
00162 
00163   mStEvent = (StEvent*) GetInputDS("StEvent");
00164   if(mStEvent)
00165     {
00166       LOG_DEBUG << "Running on StEvent ..." << endm;
00167       iret = processStEvent();
00168     }
00169   else 
00170     {
00171       StMuDstMaker *muDstMaker = (StMuDstMaker*) GetMaker("MuDst");
00172       if(muDstMaker) 
00173         {
00174           mMuDst = muDstMaker->muDst();
00175           iret = processMuDst();
00176         }
00177       else
00178         {
00179           LOG_ERROR << "No muDST is available ... "<< endm;
00180           iret = kStErr;
00181         }
00182     }
00183   
00184   MtdTrackFilterTag_st tagTable;
00185   StMaker* maskMk = GetMakerInheritsFrom("StMtdTrackingMaskMaker");
00186   tagTable.tpcSectors = (maskMk ? maskMk->UAttr("TpcSectorsByMtd") : ~0U);
00187   tagTable.isRejectEvent = (isRejectEvent() ? 1 : 0);
00188   tagTable.shouldHaveRejectEvent = shouldHaveRejectEvent();
00189   St_MtdTrackFilterTag* MtdTrackFilterTag = new St_MtdTrackFilterTag("MtdTrackFilterTag",1);
00190   MtdTrackFilterTag->AddAt(&tagTable,0);
00191   AddData(MtdTrackFilterTag);
00192 
00193   if(mStEvent)
00194     {
00195       LOG_INFO << "Is event rejected:        " << tagTable.isRejectEvent << endm;
00196       LOG_INFO << "Should event be rejected: " << tagTable.shouldHaveRejectEvent << endm;
00197     }
00198 
00199   return iret;
00200 }
00201 
00202 //_____________________________________________________________________________
00203 bool StMtdEvtFilterMaker::isRejectEvent()
00204 {
00209 
00210   return (mIsDiMuonOnly && !mIsJpsiEvent);
00211 }
00212 
00213 //_____________________________________________________________________________
00214 int StMtdEvtFilterMaker::shouldHaveRejectEvent()
00215 {
00220 
00221   if(mIsDiMuon)
00222     {
00223       if(!mIsDiMuonOnly && !mIsJpsiEvent) return 1;
00224       else return 2;
00225     }
00226   else return 0;
00227 }
00228 
00229 //_____________________________________________________________________________
00230 void StMtdEvtFilterMaker::checkTriggerIDs(const vector<unsigned int> triggers)
00231 {
00236 
00237   for(unsigned int i=0; !mIsDiMuon && i<triggers.size(); i++)
00238     {
00239       for(unsigned int j=0; !mIsDiMuon && j<mTriggerIDs.size(); j++)
00240         {
00241           mIsDiMuon = ((int)triggers[i]==mTriggerIDs[j]);
00242         }
00243     }
00244 
00245   mIsDiMuonOnly = mIsDiMuon;
00246   for(unsigned int i=0; mIsDiMuonOnly && i<triggers.size(); i++)
00247     {
00248       for(unsigned int j=0; mIsDiMuonOnly && j<mOtherTrigIDs.size(); j++)
00249         {
00250           mIsDiMuonOnly = !((int)triggers[i]==mOtherTrigIDs[j]);
00251         }
00252     }
00253 }
00254 
00255 //_____________________________________________________________________________
00256 Int_t StMtdEvtFilterMaker::processMuDst()
00257 {
00258   mhEventStat->Fill(0.5);
00259 
00260   // check trigger id
00261   vector<unsigned int> triggers = mMuDst->event()->triggerIdCollection().nominal().triggerIds();
00262   checkTriggerIDs(triggers);
00263   if(!mIsDiMuon) return kStOK;
00264   mhEventStat->Fill(1.5);
00265 
00266   // Global tracks
00267   int nNodes = mMuDst->numberOfGlobalTracks();
00268   LOG_DEBUG << "# of global tracks " << nNodes << endm;
00269 
00270   int nMuonCandidates = 0;
00271   double pt_max = 0;
00272 
00273   for(int i=0; i<nNodes; i++)
00274     {
00275       StMuTrack *gTrack = mMuDst->globalTracks(i);
00276       if(!isValidTrack(gTrack)) continue;
00277       if(!isMuonCandidate(gTrack)) continue;
00278 
00279       nMuonCandidates ++;
00280       double pt = gTrack->pt();
00281       if(pt_max<pt) pt_max = pt;
00282     }
00283 
00284   mhNMuonCandidates->Fill(nMuonCandidates);
00285   if(nMuonCandidates>=nMinMuonCandidates && pt_max>mMinTrkPtLead)
00286     {
00287       mIsJpsiEvent = kTRUE;
00288       mhEventStat->Fill(2.5);
00289     }
00290 
00291   return kStOK;
00292 }
00293 
00294 //_____________________________________________________________________________
00295 Int_t StMtdEvtFilterMaker::processStEvent()
00296 {
00297   mhEventStat->Fill(0.5);
00298 
00299   // check trigger id
00300   vector<unsigned int> triggers = mStEvent->triggerIdCollection()->nominal()->triggerIds();
00301   checkTriggerIDs(triggers);
00302   if(!mIsDiMuon) return kStOK;
00303   mhEventStat->Fill(1.5);
00304 
00305   // Global tracks
00306   StSPtrVecTrackNode& nodes = mStEvent->trackNodes();
00307   int nNodes = nodes.size();
00308   LOG_DEBUG << "# of global tracks " << nNodes << endm;
00309 
00310   int nMuonCandidates = 0;
00311   double pt_max = 0;
00312   for(int i=0; i<nNodes; i++)
00313     {
00314       StTrack *gTrack = nodes[i]->track(global);
00315 
00316       if(!isValidTrack(gTrack)) continue;
00317       if(!isMuonCandidate(gTrack)) continue;
00318 
00319       nMuonCandidates ++;
00320       double pt = gTrack->geometry()->momentum().perp();
00321       if(pt_max<pt) pt_max = pt;
00322     }
00323 
00324   mhNMuonCandidates->Fill(nMuonCandidates);
00325   if(nMuonCandidates>=nMinMuonCandidates && pt_max>mMinTrkPtLead)
00326     {
00327       mIsJpsiEvent = kTRUE;
00328       mhEventStat->Fill(2.5);
00329     }
00330 
00331   return kStOK;
00332 }
00333 
00334 //_____________________________________________________________________________
00335 bool StMtdEvtFilterMaker::isValidTrack(StTrack *track) 
00336 {
00340   if(!track) return kFALSE;
00341   StThreeVectorF mom = track->geometry()->momentum();
00342   Float_t pt      = mom.perp();
00343   int nHitsFit  = track->fitTraits().numberOfFitPoints(kTpcId);
00344   int nHitsPoss = track->numberOfPossiblePoints(kTpcId);
00345   int nHitsDedx = 0;
00346   StTpcDedxPidAlgorithm pidAlgorithm;
00347   const StParticleDefinition *pd = track->pidTraits(pidAlgorithm);
00348   if(pd && pidAlgorithm.traits())
00349     nHitsDedx = pidAlgorithm.traits()->numberOfPoints();
00350 
00351   if(pt < mMinTrkPtAll)                            return kFALSE;
00352   if(nHitsFit<mMinNHitsFit)                        return kFALSE;
00353   if(nHitsDedx<mMinNHitsDedx)                      return kFALSE;
00354   if(nHitsFit/(1.0*nHitsPoss)<mMinFitHitsFraction) return kFALSE;
00355   if(mMaxDca<1e4)
00356     {
00357       // Only check track DCA when the cut is reset
00358       // by the user
00359       
00360       StPrimaryVertex* pvtx = mStEvent->primaryVertex();
00361       if(!pvtx) return kFALSE;
00362       StThreeVectorD vtxPos = pvtx->position();
00363       StDcaGeometry* trDcaGeom = ((StGlobalTrack*) track)->dcaGeometry();
00364       if(!trDcaGeom) return kFALSE;
00365       StPhysicalHelixD dcahh = trDcaGeom->helix();
00366       double dca = dcahh.distance(vtxPos,kFALSE);
00367       if(dca>mMaxDca) return kFALSE;
00368     }
00369 
00370   return kTRUE;
00371 }
00372 
00373 //_____________________________________________________________________________
00374 bool StMtdEvtFilterMaker::isValidTrack(StMuTrack *track) 
00375 {
00379 
00380   if(!track) return kFALSE;
00381   StThreeVectorF mom = track->momentum();
00382   double pt = mom.perp();
00383   int nHitsFit  = track->nHitsFit(kTpcId);
00384   int nHitDedx  = track->nHitsDedx();
00385   int nHitPoss  = track->nHitsPoss(kTpcId);
00386   double dca    = track->dcaGlobal().mag();
00387   
00388   if(pt < mMinTrkPtAll)                           return kFALSE;
00389   if(nHitsFit<mMinNHitsFit)                       return kFALSE;
00390   if(nHitDedx<mMinNHitsDedx)                      return kFALSE;
00391   if(dca>mMaxDca)                                 return kFALSE;
00392   if(nHitsFit/(1.0*nHitPoss)<mMinFitHitsFraction) return kFALSE;
00393   return kTRUE;
00394 }
00395 
00396 //_____________________________________________________________________________
00397 bool StMtdEvtFilterMaker::isMuonCandidate(StTrack *track) 
00398 {
00402 
00403   if(!track) return kFALSE;
00404 
00405   if(mMaxNsigmaPi<1e4)
00406     {
00407       // nSigmaPi cut
00408       double nSigmaPi = -999.;
00409       StTpcDedxPidAlgorithm pidAlgorithm;
00410       const StParticleDefinition *pd = track->pidTraits(pidAlgorithm);
00411       if(pd && pidAlgorithm.traits())
00412         {
00413           static StPionPlus* Pion = StPionPlus::instance();
00414           nSigmaPi = pidAlgorithm.numberOfSigma(Pion);
00415         }
00416       if(nSigmaPi<mMinNsigmaPi || nSigmaPi>mMaxNsigmaPi) return kFALSE;
00417     }
00418 
00419   // dz cut
00420   StMtdPidTraits* mtdpid = 0;
00421   StSPtrVecTrackPidTraits& traits = track->pidTraits();
00422   for(unsigned int it=0; it<traits.size(); it++)
00423     {
00424       if (traits[it]->detector() == kMtdId)
00425         {
00426           mtdpid = dynamic_cast<StMtdPidTraits*>(traits[it]);
00427           break;
00428         }
00429     }
00430   if(!mtdpid) return kFALSE;
00431   if(mMaxDeltaZ<1e4)
00432     {
00433       double dz = mtdpid->deltaZ();
00434       if(dz > mMaxDeltaZ) return kFALSE;
00435     }
00436 
00437   return kTRUE;
00438  }
00439 
00440 //_____________________________________________________________________________
00441 bool StMtdEvtFilterMaker::isMuonCandidate(StMuTrack *track) 
00442 {
00446 
00447   if(!track) return kFALSE;
00448 
00449   if(mMaxNsigmaPi<1e4)
00450     {
00451       // nSigmaPi cut
00452       double nSigmaPi = track->nSigmaPion();
00453       if(nSigmaPi<mMinNsigmaPi || nSigmaPi>mMaxNsigmaPi)  return kFALSE;
00454     }
00455 
00456   // dz cut
00457   const StMuMtdHit *hit = track->mtdHit(); 
00458   if(!hit) return kFALSE;
00459   if(mMaxDeltaZ<1e4)
00460     {
00461       const StMuMtdPidTraits mtdPid = track->mtdPidTraits();
00462       double dz = mtdPid.deltaZ();
00463       if(dz > mMaxDeltaZ) return kFALSE;
00464     }
00465 
00466   return kTRUE;
00467  }
00468 
00469 //_____________________________________________________________________________
00470 void StMtdEvtFilterMaker::bookHistos()
00471 {
00472   mhEventStat = new TH1F("hEventStat","Event statistics",5,0,5);
00473   mhEventStat->GetXaxis()->SetBinLabel(1,"All events");
00474   mhEventStat->GetXaxis()->SetBinLabel(2,"Good trigger");
00475   mhEventStat->GetXaxis()->SetBinLabel(3,"Accepted");
00476 
00477   mhNMuonCandidates = new TH1F("hNMuonCandidates","Number of muon candidates per event;N",10,0,10);
00478 
00479   if(mSaveHistos)
00480     {
00481       AddHist(mhEventStat);
00482       AddHist(mhNMuonCandidates);
00483     }
00484 
00485 }
00486 
00487 // $Id: StMtdEvtFilterMaker.cxx,v 1.2 2015/04/23 21:10:19 marr Exp $
00488 // $Log: StMtdEvtFilterMaker.cxx,v $
00489 // Revision 1.2  2015/04/23 21:10:19  marr
00490 // 1. remove dz and pTlead cuts in the filtering by default
00491 // 2. change the number scheme for shouldHaveRejectEvent()
00492 //
00493 // Revision 1.1  2015/04/07 14:10:37  jeromel
00494 // First version of StMtdEvtFilterMaker - R.Ma - review closed 2015/04/06
00495 //
00496 //