#include <AlgRmMu.h>
Inheritance diagram for AlgRmMu:
Public Member Functions | |
| AlgRmMu () | |
| virtual | ~AlgRmMu () |
| virtual void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | Trace (const char *c) const |
|
|
Definition at line 48 of file AlgRmMu.cxx. 00049 {
00050 }
|
|
|
Definition at line 53 of file AlgRmMu.cxx. 00054 {
00055 }
|
|
||||||||||||||||
|
Implements AlgBase. Definition at line 58 of file AlgRmMu.cxx. References abs(), CandHandle::AddDaughterLink(), digit(), PlaneOutline::DistanceToEdge(), RecDataRecord< T >::FindComponent(), Registry::Get(), PlexSEIdAltL::GetBestSEId(), Truthifier::GetBiggestHit(), CandContext::GetCandRecord(), CandDigitHandle::GetCharge(), CandRecoHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandFitTrackHandle::GetEMCharge(), CandRecoHandle::GetEndPlane(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), CandRecoHandle::GetEndZ(), CandRmMuHandle::GetFitPass(), MomNavigator::GetFragment(), DigiSignal::GetHit(), CandRmMuHandle::GetIsCont(), GetMaxTrackLikePlane(), CandContext::GetMom(), CandRmMuHandle::GetMomCurv(), CandFitTrackHandle::GetMomentumCurve(), CandFitTrackHandle::GetMomentumRange(), CandRmMuHandle::GetMomRange(), CandHandle::GetNDaughters(), DigiSignal::GetNumberOfHits(), CandFitTrackHandle::GetPass(), PlexPlaneId::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), CandDigitHandle::GetQieErrorBits(), CandDigitHandle::GetRawDigitIndex(), CandEventHandle::GetShower(), Truthifier::GetSignal(), VldContext::GetSimFlag(), CandDigitHandle::GetTime(), DigiSignal::GetTruth(), CandDigitHandle::GetVaErrorBits(), RecMinos::GetVldContext(), CandRmMuHandle::GetVtxDCosX(), CandRmMuHandle::GetVtxDCosY(), CandRmMuHandle::GetVtxDCosZ(), CandRecoHandle::GetVtxDirCosU(), CandRecoHandle::GetVtxDirCosV(), CandRecoHandle::GetVtxDirCosZ(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandHandle::IsEquivalent(), TruthHelper::IsMuon(), TruthHelper::IsNeutrino(), Truthifier::IsValid(), AstUtil::LocalToHorizon(), MSG, CandRmMuHandle::SetAzimuth(), CandRmMuHandle::SetEndDistance(), CandRmMuHandle::SetEndPlane(), CandRmMuHandle::SetEndX(), CandRmMuHandle::SetEndY(), CandRmMuHandle::SetEndZ(), CandRmMuHandle::SetFitPass(), CandRmMuHandle::SetIsCont(), CandRmMuHandle::SetMaxTrkPlane(), CandRmMuHandle::SetMomCurv(), CandRmMuHandle::SetMomRange(), CandRmMuHandle::SetMomX(), CandRmMuHandle::SetMomY(), CandRmMuHandle::SetMomZ(), CandRmMuHandle::SetNMuonDig(), CandRmMuHandle::SetNMuonDigRetained(), CandRmMuHandle::SetNPlane(), CandRmMuHandle::SetNRejected(), CandRmMuHandle::SetNRejectedBoth(), CandRmMuHandle::SetNRejectedMuon(), CandRmMuHandle::SetNRejectedShw(), CandRmMuHandle::SetNRejShw(), CandRmMuHandle::SetNRejShwFakeTrk(), CandRmMuHandle::SetNRejShwMaxTrk(), CandRmMuHandle::SetNRejShwMix(), CandRmMuHandle::SetNRetained(), CandRmMuHandle::SetNRetainedBoth(), CandRmMuHandle::SetNRetainedMuon(), CandRmMuHandle::SetNRetainedShw(), CandRmMuHandle::SetNShwDig(), CandRmMuHandle::SetNShwDigAtVtx(), CandRmMuHandle::SetNShwDigRetained(), CandRmMuHandle::SetNShwDigRetainedAtVtx(), CandRmMuHandle::SetNShwPE(), CandRmMuHandle::SetNShwPEAtVtx(), CandRmMuHandle::SetNShwPERetained(), CandRmMuHandle::SetNShwPERetainedAtVtx(), CandRmMuHandle::SetOrigEvtIndex(), CandRmMuHandle::SetPass(), CandRmMuHandle::SetPERetained(), CandRmMuHandle::SetPERetainedBoth(), CandRmMuHandle::SetPERetainedMuon(), CandRmMuHandle::SetPERetainedShw(), CandRmMuHandle::SetQP(), CandRmMuHandle::SetReasonForKeeping(), CandRmMuHandle::SetShwCharge(), CandRmMuHandle::SetShwEndPlane(), CandRmMuHandle::SetShwEndX(), CandRmMuHandle::SetShwEndY(), CandRmMuHandle::SetShwEndZ(), CandRmMuHandle::SetShwNPlane(), CandRmMuHandle::SetShwVtxPlane(), CandRmMuHandle::SetShwVtxX(), CandRmMuHandle::SetShwVtxY(), CandRmMuHandle::SetShwVtxZ(), CandRmMuHandle::SetVtxDCosX(), CandRmMuHandle::SetVtxDCosY(), CandRmMuHandle::SetVtxDCosZ(), CandRmMuHandle::SetVtxDistance(), CandRmMuHandle::SetVtxPlane(), CandRmMuHandle::SetVtxX(), CandRmMuHandle::SetVtxY(), CandRmMuHandle::SetVtxZ(), CandRmMuHandle::SetZenith(), TrackEndContained(), and DigiScintHit::TrackId(). 00059 {
00060
00061 MSG("AlgRmMu", Msg::kDebug) << "Starting AlgRmMu::RunAlg()" << endl;
00062
00063 assert(ch.InheritsFrom("CandRmMuHandle"));
00064 CandRmMuHandle &crmh = dynamic_cast<CandRmMuHandle &>(ch);
00065 assert(cx.GetDataIn());
00066 assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00067
00068 const CandEventHandle *event = 0;
00069 const CandTrackHandle *track = 0;
00070
00071 const CandDigitListHandle *digitlist = 0;
00072 const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00073 for (Int_t i=0; i<=tary->GetLast(); i++) {
00074 TObject *tobj = tary->At(i);
00075 if(tobj->InheritsFrom("CandEventHandle")) {
00076 event = dynamic_cast<CandEventHandle*> (tobj);
00077 MSG("AlgRmMu",Msg::kDebug) << "Got Event" << endl;
00078 }
00079 if(tobj->InheritsFrom("CandTrackHandle")) {
00080 track = dynamic_cast<CandTrackHandle*> (tobj);
00081 MSG("AlgRmMu",Msg::kDebug) << "Got Track" << endl;
00082 }
00083 if(tobj->InheritsFrom("CandDigitListHandle")) {
00084 digitlist = dynamic_cast<CandDigitListHandle*> (tobj);
00085 MSG("AlgRmMu",Msg::kDebug) << "Got Digitlist" << endl;
00086 }
00087 }
00088
00089 if(event==NULL || digitlist==NULL) {
00090 MSG("AlgRmMu",Msg::kError) << " Bailing out of Event event = " << event
00091 << " digitlist = " << digitlist <<endl;
00092 return;
00093 }
00094
00095
00096 if(track==NULL) {
00097 MSG("AlgRmMu",Msg::kDebug) << "No track to remove" << endl;
00098 crmh.SetVtxX(0);
00099 crmh.SetVtxY(0);
00100 crmh.SetVtxZ(0);
00101 crmh.SetVtxDistance(0);
00102 crmh.SetEndX(0);
00103 crmh.SetEndY(0);
00104 crmh.SetEndZ(0);
00105 crmh.SetEndDistance(0);
00106 crmh.SetVtxPlane(0);
00107 crmh.SetEndPlane(0);
00108 crmh.SetNPlane(0);
00109 crmh.SetMomRange(0);
00110 crmh.SetMomCurv(0);
00111 crmh.SetVtxDCosX(0);
00112 crmh.SetVtxDCosY(0);
00113 crmh.SetVtxDCosZ(0);
00114 crmh.SetZenith(0);
00115 crmh.SetAzimuth(0);
00116 crmh.SetFitPass(0);
00117 crmh.SetIsCont(0);
00118 crmh.SetPass(0);
00119 crmh.SetMomX(0);
00120 crmh.SetMomY(0);
00121 crmh.SetMomZ(0);
00122 crmh.SetQP(0);
00123 crmh.SetOrigEvtIndex(-1);
00124 crmh.SetMaxTrkPlane(9999);
00125 crmh.SetNMuonDig(0);
00126 crmh.SetNMuonDigRetained(0);
00127 crmh.SetNShwDig(0);
00128 crmh.SetNShwDigRetained(0);
00129 crmh.SetNShwDigAtVtx(0);
00130 crmh.SetNShwDigRetainedAtVtx(0);
00131 crmh.SetNShwPE(0);
00132 crmh.SetNShwPERetained(0);
00133 crmh.SetNShwPEAtVtx(0);
00134 crmh.SetNShwPERetainedAtVtx(0);
00135 crmh.SetNRetained(0);
00136 crmh.SetNRetainedMuon(0);
00137 crmh.SetNRetainedShw(0);
00138 crmh.SetNRetainedBoth(0);
00139 crmh.SetPERetained(0);
00140 crmh.SetPERetainedMuon(0);
00141 crmh.SetPERetainedShw(0);
00142 crmh.SetPERetainedBoth(0);
00143 crmh.SetNRejected(0);
00144 crmh.SetNRejectedMuon(0);
00145 crmh.SetNRejectedShw(0);
00146 crmh.SetNRejectedBoth(0);
00147 crmh.SetNRejShw(0);
00148 crmh.SetNRejShwMaxTrk(0);
00149 crmh.SetNRejShwFakeTrk(0);
00150 crmh.SetNRejShwMix(0);
00151 return;
00152 }
00153
00154 CandRecord* record = dynamic_cast<CandRecord*>(cx.GetCandRecord());
00155 const VldContext &vldc = *(record->GetVldContext());
00156
00157 // added by J. Ling
00158 // get the first shower in the event
00159
00160 const CandShowerHandle* shower = event->GetShower(0);
00161
00162 if (shower==NULL) {
00163 MSG("AlgRmMu",Msg::kDebug) << "No shower in event" << endl;
00164 crmh.SetShwVtxX(0);
00165 crmh.SetShwVtxY(0);
00166 crmh.SetShwVtxZ(0);
00167 crmh.SetShwEndX(0);
00168 crmh.SetShwEndY(0);
00169 crmh.SetShwEndZ(0);
00170 crmh.SetShwVtxPlane(0);
00171 crmh.SetShwEndPlane(0);
00172 crmh.SetShwNPlane(0);
00173 crmh.SetShwCharge(0);
00174 } else {
00175 crmh.SetShwVtxX(0.707106781*(shower->GetVtxU() - shower->GetVtxV()));
00176 crmh.SetShwVtxY(0.707106781*(shower->GetVtxU() + shower->GetVtxV()));
00177 crmh.SetShwVtxZ(shower->GetVtxZ());
00178 crmh.SetShwEndX(0.707106781*(shower->GetEndU() - shower->GetEndV()));
00179 crmh.SetShwEndY(0.707106781*(shower->GetEndU() + shower->GetEndV()));
00180 crmh.SetShwEndZ(shower->GetEndZ());
00181 crmh.SetShwVtxPlane(shower->GetVtxPlane());
00182 crmh.SetShwEndPlane(shower->GetEndPlane());
00183 crmh.SetShwNPlane(TMath::Abs(shower->GetEndPlane() - shower->GetVtxPlane())+1);
00184 crmh.SetShwCharge(shower->GetCharge(CalStripType::kMIP));
00185 }
00186
00187 const CandFitTrackHandle* fittrack = dynamic_cast<const CandFitTrackHandle*>(track);
00188
00189 // added by J. Ling
00190
00191 float vtxdist = -1;
00192 float enddist = -1;
00193 float x,y;
00194
00195 PlexHandle ph(vldc);
00196 PlaneOutline ppo;
00197 Detector::Detector_t Det = vldc.GetDetector();
00198
00199 PlexPlaneId vtxplnid(Det,fittrack->GetVtxPlane(),false);
00200 if (vtxplnid.IsValid()) {
00201 ppo.DistanceToEdge(0.707106781*(fittrack->GetVtxU()-fittrack->GetVtxV()),0.707106781*(fittrack->GetVtxU()+fittrack->GetVtxV()),vtxplnid.GetPlaneView(),vtxplnid.GetPlaneCoverage(),vtxdist,x,y);
00202 }
00203 if (vtxdist<0.) vtxdist = 0;
00204
00205 PlexPlaneId endplnid(Det,fittrack->GetEndPlane(),false);
00206 if (vtxplnid.IsValid()) {
00207 ppo.DistanceToEdge(0.707106781*(fittrack->GetEndU()-fittrack->GetEndV()),0.707106781*(fittrack->GetEndU()+fittrack->GetEndV()),endplnid.GetPlaneView(),endplnid.GetPlaneCoverage(),enddist,x,y);
00208 }
00209 if (enddist<0.) enddist = 0;
00210
00211 //set muon parameters
00212 crmh.SetVtxX(0.707106781*(fittrack->GetVtxU() - fittrack->GetVtxV()));
00213 crmh.SetVtxY(0.707106781*(fittrack->GetVtxU() + fittrack->GetVtxV()));
00214 crmh.SetVtxZ(fittrack->GetVtxZ());
00215 crmh.SetVtxDistance(vtxdist);
00216
00217 crmh.SetEndX(0.707106781*(fittrack->GetEndU() - fittrack->GetEndV()));
00218 crmh.SetEndY(0.707106781*(fittrack->GetEndU() + fittrack->GetEndV()));
00219 crmh.SetEndZ(fittrack->GetEndZ());
00220 crmh.SetEndDistance(enddist);
00221
00222 crmh.SetVtxPlane(fittrack->GetVtxPlane());
00223 crmh.SetEndPlane(fittrack->GetEndPlane());
00224 crmh.SetNPlane(TMath::Abs(fittrack->GetEndPlane() - fittrack->GetVtxPlane())+1);
00225 crmh.SetMomRange(fittrack->GetMomentumRange());
00226 crmh.SetMomCurv(fittrack->GetMomentumCurve());
00227 crmh.SetVtxDCosX(0.707106781*(fittrack->GetVtxDirCosU() - fittrack->GetVtxDirCosV()));
00228 crmh.SetVtxDCosY(0.707106781*(fittrack->GetVtxDirCosU() + fittrack->GetVtxDirCosV()));
00229 crmh.SetVtxDCosZ(fittrack->GetVtxDirCosZ());
00230
00231 double dcosx = -0.707106781 * (fittrack->GetVtxDirCosU() - fittrack->GetVtxDirCosV());
00232 double dcosy = -0.707106781 * (fittrack->GetVtxDirCosU() + fittrack->GetVtxDirCosV());
00233 double dcosz = -(fittrack->GetVtxDirCosZ());
00234
00235 double altitude, azimuth;
00236 AstUtil::LocalToHorizon(dcosx,dcosy,dcosz,Det,altitude,azimuth);
00237
00238 double zenith = 90. - altitude;
00239
00240 crmh.SetZenith((Float_t)zenith);
00241 crmh.SetAzimuth((Float_t)azimuth);
00242
00243 crmh.SetFitPass(fittrack->GetPass());
00244 crmh.SetIsCont(TrackEndContained(fittrack));
00245 Int_t pass = 1;
00246 Float_t momentum = crmh.GetMomRange();
00247 if(!crmh.GetIsCont()){
00248 pass = crmh.GetFitPass();
00249 momentum = crmh.GetMomCurv();
00250 }
00251 crmh.SetPass(pass);
00252 crmh.SetMomX(crmh.GetVtxDCosX()*TMath::Abs(momentum));
00253 crmh.SetMomY(crmh.GetVtxDCosY()*TMath::Abs(momentum));
00254 crmh.SetMomZ(crmh.GetVtxDCosZ()*TMath::Abs(momentum));
00255 if(fittrack->GetMomentumCurve() > 0.001)
00256 crmh.SetQP(fittrack->GetEMCharge()/fittrack->GetMomentumCurve());
00257
00258 //hold some vectors for use later:
00259 std::vector<Float_t> stripRetained(event->GetNDaughters());
00260 std::vector<Bool_t> stripInTrack(event->GetNDaughters());
00261 for(int i=0;i<event->GetNDaughters();i++) {
00262 stripRetained[i] = 0.0;
00263 stripInTrack[i] = 0;
00264 }
00265
00266 //keep map of info to store about hit removal
00267 //keep track of digit<->strip associations
00268 std::map<const CandDigitHandle*,Int_t> infoToStore;
00269 std::multimap<Int_t,const CandDigitHandle*> eventStripMatch;
00270
00271 Int_t detector = vldc.GetDetector();
00272 //add daughter links to appropriate digits
00273 TIter digitIter(digitlist->GetDaughterIterator());
00274 TIter stpIter(event->GetDaughterIterator());
00275 UInt_t stpCounter = 0;
00276 while(CandDigitHandle* digit =
00277 dynamic_cast<CandDigitHandle*>(digitIter())){
00278 infoToStore[digit] = 0;
00279 stpIter.Reset();
00280 stpCounter = 0;
00281 while(const CandStripHandle* strip =
00282 dynamic_cast<const CandStripHandle*>(stpIter())){
00283 TIter stpdigitIter(strip->GetDaughterIterator());
00284 while(const CandDigitHandle* stpdigit =
00285 dynamic_cast<CandDigitHandle*>(stpdigitIter())){
00286 if( ( digit->GetRawDigitIndex() ==
00287 stpdigit->GetRawDigitIndex()) &&
00288 ( digit->GetTime(CalTimeType::kNone) ==
00289 stpdigit->GetTime(CalTimeType::kNone) ) ){
00290 crmh.AddDaughterLink(*digit);
00291 if( (detector == 1 && digit->GetVaErrorBits()!=0) ||
00292 (detector == 2 && digit->GetQieErrorBits()!=0) ) {
00293 stripRetained[stpCounter] = -( digit->GetCharge() /
00294 stpdigit->GetCharge() );
00295 }
00296 else stripRetained[stpCounter] = 1.0;
00297 eventStripMatch.insert(pair<Int_t,const CandDigitHandle*>(stpCounter,digit));
00298 }
00299 }
00300 stpCounter+=1;
00301 }
00302 }
00303
00304 //find out which event strips are in the track:
00305 TIter trkstpItr(fittrack->GetDaughterIterator());
00306 stpIter.Reset();
00307 stpCounter = 1;
00308 while(const CandStripHandle* strip =
00309 dynamic_cast<const CandStripHandle*>(stpIter())) {
00310 trkstpItr.Reset();
00311 while(const CandStripHandle* trkStrip =
00312 dynamic_cast<const CandStripHandle*>(trkstpItr())) {
00313 if(strip->IsEquivalent(trkStrip)) {
00314 stripInTrack[stpCounter] = 1;
00315 }
00316 }
00317 stpCounter += 1;
00318 }
00319
00320 Int_t cMaxTrackLikePlanes = 6;
00321 if(!ac.Get("MaxTrackLikePlanes",cMaxTrackLikePlanes))
00322 cMaxTrackLikePlanes=6;
00323 int maxtrkplane = GetMaxTrackLikePlane(event,track,cMaxTrackLikePlanes);
00324 crmh.SetMaxTrkPlane(maxtrkplane);
00325
00326 //fill reco-related hit removal info:
00327 //loop over event strips:
00328 for(stpCounter=0;stpCounter<(UInt_t)(event->GetNDaughters());stpCounter++){
00329 bool is_retained = false;
00330 if(stripRetained[stpCounter]>0) is_retained = true;
00331 bool is_track = false;
00332 if(stripInTrack[stpCounter]==1) is_track = true;
00333 bool is_retained_scaled = false;
00334 Float_t scale_factor = 1.0;
00335 if(stripRetained[stpCounter]<0) {
00336 is_retained_scaled = true;
00337 scale_factor = TMath::Abs(stripRetained[stpCounter]);
00338 }
00339 //fill in bit-packed info:
00340 std::multimap<Int_t,const CandDigitHandle*>::iterator esm_beg;
00341 std::multimap<Int_t,const CandDigitHandle*>::iterator esm_end;
00342 esm_beg = eventStripMatch.find(stpCounter);
00343 if(esm_beg != eventStripMatch.end()){
00344 esm_end = eventStripMatch.upper_bound(stpCounter);
00345 for ( ; esm_beg != esm_end; ++esm_beg) {
00346 infoToStore[esm_beg->second] = (int(scale_factor*RmMuMask::kRMMU_SCL_FACT_SCALE)<<RmMuMask::kRMMU_NUM_SHIFT);
00347 infoToStore[esm_beg->second] |= RmMuMask::kRMMU_ISRETAIN_MASK;
00348 if (is_retained_scaled)
00349 infoToStore[esm_beg->second] |= RmMuMask::kRMMU_ISSCALED_MASK;
00350 if (is_track)
00351 infoToStore[esm_beg->second] |= RmMuMask::kRMMU_INRECOTRK_MASK;
00352 }
00353 }
00354 }
00355
00356 //do truth analysis?
00357 if(vldc.GetSimFlag()==SimFlag::kMC){
00358 //do truth analysis
00359 Truthifier* truth = new Truthifier(cx.GetMom());
00360 if(truth->IsValid()){
00361 TruthHelper* helper = new TruthHelper(cx.GetMom());
00362 SimSnarlRecord *simsnarl =
00363 dynamic_cast<SimSnarlRecord*>(cx.GetMom()->GetFragment("SimSnarlRecord"));
00364 if(simsnarl) { //simsnarl
00365 const TClonesArray* ctca =
00366 dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","StdHep"));
00367 // const TClonesArray* neukinarray =
00368 // dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","NeuKinList"));
00369 // const REROOT_NeuKin* neukin =
00370 // dynamic_cast<const REROOT_NeuKin*>(neukinarray->At(0));
00371 // assert(neukin);
00372
00373 //variables to fill:
00374 Int_t nMuonDig = 0;
00375 Int_t nMuonDigRetained = 0;
00376 Int_t nShwDig = 0;
00377 Int_t nShwDigRetained = 0;
00378 Int_t nShwDigAtVtx = 0;
00379 Int_t nShwDigRetainedAtVtx = 0;
00380 Float_t peShw = 0.;
00381 Float_t peShwRetained = 0.;
00382 Float_t peShwAtVtx = 0.;
00383 Float_t peShwRetainedAtVtx = 0.;
00384 Int_t nRetained = 0;
00385 Int_t nRetainedMuon = 0;
00386 Int_t nRetainedShw = 0;
00387 Int_t nRetainedBoth = 0;
00388 Float_t peRetained = 0.;
00389 Float_t peRetainedMuon = 0.;
00390 Float_t peRetainedShw = 0.;
00391 Float_t peRetainedBoth = 0.;
00392 Int_t nRejected = 0;
00393 Int_t nRejectedMuon = 0;
00394 Int_t nRejectedShw = 0;
00395 Int_t nRejectedBoth = 0;
00396 Int_t nRejShw = 0;
00397 Int_t nRejShwMaxTrk = 0 ;
00398 Int_t nRejShwFakeTrk = 0;
00399 Int_t nRejShwMix = 0 ;
00400
00401 //loop over event strips:
00402 stpIter.Reset();
00403 stpCounter = 0;
00404 while(const CandStripHandle* strip =
00405 dynamic_cast<const CandStripHandle*>(stpIter())){ //loop over strips
00406 TIter digitItr(strip->GetDaughterIterator());
00407 while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) { //loop over digits
00408 const int planeno = digit->GetPlexSEIdAltL().GetBestSEId().GetPlane();
00409 //const int stripno = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
00410 const CandDigitHandle tcdh = *digit;
00411 //const CandDeMuxDigitHandle* demuxdigit =
00412 //dynamic_cast<const CandDeMuxDigitHandle*>(digit);
00413
00414 bool is_retained = false;
00415 if(stripRetained[stpCounter]>0) is_retained = true;
00416 bool is_track = false;
00417 if(stripInTrack[stpCounter]==1) is_track = true;
00418 bool is_retained_scaled = false;
00419 if(stripRetained[stpCounter]<0) is_retained_scaled = true;
00420 //bool is_recoed_xtalk = (demuxdigit->GetDeMuxDigitFlagWord()&CandDeMuxDigit::kXTalk);
00421 bool is_shower = 0;
00422 bool is_muon = 0;
00423 bool is_physics = 0;
00424
00425 const DigiSignal * signal = truth->GetSignal(tcdh);
00426 const DigiScintHit* biggest_hit = truth->GetBiggestHit(tcdh);
00427 TParticle* biggest_part = NULL;
00428 if(biggest_hit!=NULL){
00429 Int_t biggest_track = abs(biggest_hit->TrackId());
00430 biggest_part = dynamic_cast<TParticle*>((*ctca)[biggest_track]);
00431 }
00432 if(signal!=NULL){
00433 is_physics = ((signal->GetTruth() & DigiSignal::kGenuine)!=0);
00434 for (UInt_t i=0;i<signal->GetNumberOfHits();i++){
00435 const DigiScintHit * hit = signal->GetHit(i);
00436 if(hit){
00437 Int_t track = abs(hit->TrackId());
00438 TParticle* part = dynamic_cast<TParticle*>((*ctca)[track]);
00439 if(!helper->IsNeutrino(part) && !helper->IsMuon(part)){
00440 is_shower=true;
00441 }
00442 if(!helper->IsNeutrino(part) && helper->IsMuon(part)){
00443 is_muon=true;
00444 }
00445 }
00446 }
00447 }
00448
00449 //fill in truth-related hit removal info:
00450 std::multimap<Int_t,const CandDigitHandle*>::iterator esm_beg;
00451 std::multimap<Int_t,const CandDigitHandle*>::iterator esm_end;
00452 esm_beg = eventStripMatch.find(stpCounter);
00453 if(esm_beg != eventStripMatch.end()){
00454 esm_end = eventStripMatch.upper_bound(stpCounter);
00455 for ( ; esm_beg != esm_end; ++esm_beg) {
00456 if (is_muon) infoToStore[esm_beg->second] |= RmMuMask::kRMMU_ISMU_MASK;
00457 if (is_shower) infoToStore[esm_beg->second] |= RmMuMask::kRMMU_ISSHW_MASK;
00458 if (is_physics) infoToStore[esm_beg->second] |= RmMuMask::kRMMU_ISPHYS_MASK;
00459 }
00460 }
00461
00462 if(is_physics){
00463 if(is_shower){
00464 nShwDig++;
00465 peShw+=digit->GetCharge(CalDigitType::kPE);
00466 if(planeno<maxtrkplane){
00467 nShwDigAtVtx++;
00468 peShwAtVtx+=digit->GetCharge(CalDigitType::kPE);
00469 }
00470 if(is_retained || is_retained_scaled){
00471 nShwDigRetained++;
00472 peShwRetained+=digit->GetCharge(CalDigitType::kPE);
00473 if(planeno<maxtrkplane){
00474 nShwDigRetainedAtVtx++;
00475 peShwRetainedAtVtx+=digit->GetCharge(CalDigitType::kPE);
00476 }
00477 }
00478 else{
00479 nRejShw++;
00480 if(planeno>maxtrkplane) nRejShwMaxTrk++;
00481 if(is_muon) nRejShwMix++;
00482 if(!is_muon && is_track) nRejShwFakeTrk++;
00483 }
00484 } //is shower
00485 if(is_muon){
00486 nMuonDig++;
00487 if(is_retained || is_retained_scaled) nMuonDigRetained++;
00488 } //is muon
00489 if(is_retained || is_retained_scaled){
00490 nRetained++;
00491 peRetained+=digit->GetCharge(CalDigitType::kPE);
00492 if(is_muon){
00493 nRetainedMuon++;
00494 peRetainedMuon+=digit->GetCharge(CalDigitType::kPE);
00495 }
00496 if(is_shower){
00497 nRetainedShw++;
00498 peRetainedShw+=digit->GetCharge(CalDigitType::kPE);
00499 }
00500 if(is_muon && is_shower){
00501 nRetainedBoth++;
00502 peRetainedBoth+=digit->GetCharge(CalDigitType::kPE);
00503 }
00504 }
00505 else{
00506 nRejected++;
00507 if(is_muon && !is_shower)nRejectedMuon++;
00508 if(!is_muon && is_shower) nRejectedShw++;
00509 if(is_muon && is_shower) nRejectedBoth++;
00510 }
00511 }
00512 }
00513 stpCounter += 1;
00514 }
00515
00516 crmh.SetNMuonDig(nMuonDig);
00517 crmh.SetNMuonDigRetained(nMuonDigRetained);
00518 crmh.SetNShwDig(nShwDig);
00519 crmh.SetNShwDigRetained(nShwDigRetained);
00520 crmh.SetNShwDigAtVtx(nShwDigAtVtx);
00521 crmh.SetNShwDigRetainedAtVtx(nShwDigRetainedAtVtx);
00522 crmh.SetNShwPE(peShw);
00523 crmh.SetNShwPERetained(peShwRetained);
00524 crmh.SetNShwPEAtVtx(peShwAtVtx);
00525 crmh.SetNShwPERetainedAtVtx(peShwRetainedAtVtx);
00526 crmh.SetNRetained(nRetained);
00527 crmh.SetNRetainedMuon(nRetainedMuon);
00528 crmh.SetNRetainedShw(nRetainedShw);
00529 crmh.SetNRetainedBoth(nRetainedBoth);
00530 crmh.SetPERetained(peRetained);
00531 crmh.SetPERetainedMuon(peRetainedMuon);
00532 crmh.SetPERetainedShw(peRetainedShw);
00533 crmh.SetPERetainedBoth(peRetainedBoth);
00534 crmh.SetNRejected(nRejected);
00535 crmh.SetNRejectedMuon(nRejectedMuon);
00536 crmh.SetNRejectedShw(nRejectedShw);
00537 crmh.SetNRejectedBoth(nRejectedBoth);
00538 crmh.SetNRejShw(nRejShw);
00539 crmh.SetNRejShwMaxTrk(nRejShwMaxTrk);
00540 crmh.SetNRejShwFakeTrk(nRejShwFakeTrk);
00541 crmh.SetNRejShwMix(nRejShwMix);
00542 }
00543 delete helper;
00544 }
00545 delete truth;
00546 }
00547
00548 //set the reason for keeping thing:
00549 std::map<const CandDigitHandle*,Int_t>::iterator beg = infoToStore.begin();
00550 std::map<const CandDigitHandle*,Int_t>::iterator end = infoToStore.end();
00551 while(beg!=end){
00552 crmh.SetReasonForKeeping(beg->first,beg->second);
00553 beg++;
00554 }
00555
00556 }
|
|
|
Reimplemented from AlgBase. Definition at line 559 of file AlgRmMu.cxx. 00560 {
00561 }
|
1.3.9.1