loon: AlgRmMu.cxx Source File

00001 
00002 // $Id: AlgRmMu.cxx,v 1.10 2009/12/08 23:04:08 rhatcher Exp $
00003 //
00004 // AlgRmMu.cxx
00005 //
00007 
00008 #include <cassert>
00009 
00010 #include "Algorithm/AlgConfig.h"
00011 #include "Candidate/CandContext.h"
00012 #include "MuonRemoval/AlgRmMu.h"
00013 #include "MuonRemoval/CandRmMuHandle.h"
00014 #include "MuonRemoval/SelectEvent.h"
00015 #include "Conventions/PlaneView.h"
00016 #include "Conventions/DetectorType.h"
00017 #include "DataUtil/PlaneOutline.h"
00018 #include "MessageService/MsgService.h"
00019 #include "MinosObjectMap/MomNavigator.h"
00020 #include "RecoBase/CandStripHandle.h"
00021 #include "RecoBase/LinearFit.h"
00022 #include "UgliGeometry/UgliGeomHandle.h"
00023 #include "Validity/VldContext.h"
00024 #include "RecoBase/CandEventHandle.h"
00025 #include "RecoBase/CandTrackHandle.h"
00026 #include "RecoBase/CandFitTrackHandle.h"
00027 #include "RecoBase/CandShowerHandle.h"
00028 #include "CandDigit/CandDigitHandle.h"
00029 #include "CandDigit/CandDigitListHandle.h"
00030 #include "CandDigit/CandDeMuxDigitHandle.h"
00031 #include "CandDigit/CandDeMuxDigitListHandle.h"
00032 #include "MuonRemoval/SelectEvent.h"
00033 #include "DataUtil/Truthifier.h"
00034 #include "DataUtil/TruthHelper.h"
00035 #include "Record/SimSnarlRecord.h"
00036 #include "REROOT_Classes/REROOT_NeuKin.h"
00037 #include "TClonesArray.h"
00038 #include "Plex/PlexHandle.h"
00039 #include "AstroUtil/Ast.h"
00040 #include "AstroUtil/AstCoordinate.h"
00041 
00042 ClassImp(AlgRmMu)
00043 
00044 //______________________________________________________________________
00045 CVSID("$Id: AlgRmMu.cxx,v 1.10 2009/12/08 23:04:08 rhatcher Exp $");
00046 
00047 //______________________________________________________________________
00048 AlgRmMu::AlgRmMu()
00049 {
00050 }
00051 
00052 //______________________________________________________________________
00053 AlgRmMu::~AlgRmMu()
00054 {
00055 }
00056 
00057 //______________________________________________________________________
00058 void AlgRmMu::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
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::kSynopsis) << " 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 }
00557 
00558 //______________________________________________________________________
00559 void AlgRmMu::Trace(const char * /* c */) const
00560 {
00561 }