loon: AlgCosmicMuonRemoval.cxx Source File

00001 //  Added by J. Ling
00002 //  An algorithm of get rid of cosmic ray muon from the cosmic event
00003 //  in order to analyze the EM shower
00004 
00005 #include <cassert>
00006 #include <cmath>
00007 
00008 #include "Algorithm/AlgConfig.h"
00009 #include "Algorithm/AlgFactory.h"
00010 #include "Algorithm/AlgHandle.h"
00011 #include "CandData/CandRecord.h"
00012 #include "MuonRemoval/AlgCosmicMuonRemoval.h"
00013 
00014 #include "CandDigit/CandDeMuxDigitHandle.h"
00015 #include "CandDigit/CandDeMuxDigitListHandle.h"
00016 
00017 #include "RecoBase/CandSliceHandle.h"
00018 #include "RecoBase/CandTrackListHandle.h"
00019 #include "RecoBase/CandTrackHandle.h"
00020 #include "RecoBase/CandShowerListHandle.h"
00021 #include "RecoBase/CandShowerHandle.h"
00022 #include "RecoBase/CandEventListHandle.h"
00023 #include "RecoBase/CandEventHandle.h"
00024 
00025 #include "RecoBase/CandStripHandle.h"
00026 #include "CandData/CandHeader.h"
00027 #include "Candidate/CandContext.h"
00028 #include "MessageService/MsgService.h"
00029 #include "Plex/PlexSEIdAltL.h"
00030 #include "Plex/PlexSEIdAltLItem.h"
00031 #include "Validity/VldContext.h"
00032 #include "Validity/VldTimeStamp.h"
00033 
00034 #include "Plex/PlexHandle.h"
00035 
00036 #include "RawDigitInfo.h"
00037 #include "TClonesArray.h"
00038 
00039 #include "SelectEvent.h"
00040 
00041 CVSID("$Id: AlgCosmicMuonRemoval.cxx,v 1.1 2008/09/17 02:56:19 tjyang Exp $");
00042  
00043 ClassImp(AlgCosmicMuonRemoval)
00044 
00045 AlgCosmicMuonRemoval::AlgCosmicMuonRemoval()
00046 {
00047   Reset();
00048   cTrackLikePlanes = 5;
00049   cShowerLikePlanes = 2;
00050 }
00051 
00052 
00053 AlgCosmicMuonRemoval::~AlgCosmicMuonRemoval()
00054 {
00055 }
00056 
00057 void AlgCosmicMuonRemoval::Reset()
00058 {
00059   for(int ipln=0; ipln<500; ++ipln){
00060     fTrkQ[ipln] = 0;
00061     fPlnQ[ipln] = 0;
00062     fTrkMIPDcosz[ipln] = 0;
00063     fTrkStrip[ipln] = 0;
00064     fNTrackStrip[ipln] = 0;
00065     fNTotalStrip[ipln] = 0;
00066   }
00067   fTrkMIPCalibOkay = false;
00068   fTrkMIPDcoszMean = 0;
00069 }
00070 
00071 void AlgCosmicMuonRemoval::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
00072 {
00073    //
00074    //set up the sub-algs
00075    //  
00076   AlgFactory &algfactory = AlgFactory::GetInstance();
00077   const char *tmpcs = 0;
00078   Int_t tmpi = 5;
00079   Int_t tmpt = 2;
00080   
00081   const char *pass_algorithm = 0;
00082   const char *pass_algconfig = 0;
00083   
00084   if (ac.Get("PassThruAlgorithm", tmpcs)) pass_algorithm = tmpcs;
00085   if (ac.Get("PassThruAlgConfig", tmpcs)) pass_algconfig = tmpcs;
00086   if (ac.Get("TrackLikePlanes", tmpi)) cTrackLikePlanes = tmpi;
00087   if (ac.Get("ShowerLikePlanes", tmpt)) cShowerLikePlanes = tmpt;
00088 
00089   AlgHandle pass_alg = algfactory.GetAlgHandle(pass_algorithm, pass_algconfig);
00090   
00091   CandContext cxx(this, cx.GetMom());
00092   cxx.SetCandRecord(cx.GetCandRecord());
00093   
00094   MSG("RmMu",Msg::kDebug) << " AlgCosmicMuonRemoval::RunAlg() " <<endl;
00095   //MSG("RmMu",Msg::kDebug) << "   SubAlg         :  " <<pass_algorithm<<endl;
00096   //MSG("RmMu",Msg::kDebug) << "   SubAlg Config  :  " <<pass_algconfig<<endl;
00097  
00098   //
00099   //Do event setup
00100   //
00101   Reset();
00102 
00103   TObjArray *input = (TObjArray *)(cx.GetDataIn());  
00104   const CandRecord* record = dynamic_cast<const CandRecord*>(input->At(0));
00105   assert(record);
00106   const CandHeader* header  = dynamic_cast<const CandHeader*>(record->GetHeader());
00107   if(header){
00108     MSG("RmMu",Msg::kDebug) << " Snarl: " << header->GetRun() << " / "<< header->GetSnarl() <<endl; 
00109   }
00110 
00111   const CandEventListHandle * eventlist = dynamic_cast<CandEventListHandle*>(record->FindCandHandle("CandEventListHandle"));
00112   const CandDigitListHandle* digitlist = dynamic_cast<const CandDigitListHandle*>(record->FindCandHandle("CandDigitListHandle"));
00113 
00114   if(eventlist==NULL || digitlist==NULL){
00115     MSG("RmMu",Msg::kError) << " Bailing out of Event eventlist = " << eventlist<< " digitlist = " << digitlist <<endl; 
00116     return; //bail out   
00117   }
00118 
00119 
00120   //
00121   //A list of digits to chuck/modify
00122   //
00123   std::vector<int> digitidx;
00124   std::vector<float> digitweight;
00125   long ntrackdigits = 0;
00126 
00127   TIter event_iter(eventlist->GetDaughterIterator());
00128   while( const CandEventHandle* event = dynamic_cast<const CandEventHandle*>(event_iter()) ){
00129     MSG("RmMu",Msg::kDebug) << " New event " <<endl; 
00130     if(SelectCosmicEvent(event)){      
00131       MSG("RmMu",Msg::kDebug) << "   Event is selected for cosmic muon removal " <<endl; 
00132       
00133       Reset();
00134       
00135       const CandTrackHandle* track = GetRemovableTrack(event);
00136       if(!track){
00137         MSG("RmMu",Msg::kError) << "   NO removalable track! "<<endl;
00138       }
00139 
00140       const CandShowerHandle* shower = event->GetShower(0);
00141       if(!shower) {
00142         MSG("RmMu",Msg::kError) << "   NO shower in event! " << endl;  
00143       }
00144       //
00145       //Fill some local info
00146       //
00147       FillEvtInfo(event,track,shower);
00148       
00149       MSG("RmMu",Msg::kDebug) << "   fTrkMIPCalibOkay:  "<< fTrkMIPCalibOkay << "   fTrkMIPDcoszMean   : " << fTrkMIPDcoszMean <<endl;
00150       
00151       //TY loop through hits in the slice instead of event
00152       const CandSliceHandle *csh = event->GetCandSlice();
00153                               
00154       TIter stpIter(csh->GetDaughterIterator());
00155       while(const CandStripHandle* strip = dynamic_cast<const CandStripHandle*>(stpIter())){
00156         const int planeno = strip->GetPlane();
00157         const int stripno = strip->GetStrip();
00158         const bool isintrack =(track->FindDaughter(strip)!=NULL);
00159 
00160         bool notmuon = fTrkMIPDcosz[planeno]<= 0.3;     
00161         bool isinshower = (planeno>=fMinShowerPlane && planeno<=fMaxShowerPlane) || (planeno>=fMinShowerPlane-cShowerLikePlanes && planeno<fMinShowerPlane&& fNTotalStrip[planeno]>1) || (planeno>fMaxShowerPlane && planeno<=fMaxShowerPlane+cShowerLikePlanes && fNTotalStrip[planeno]>1);        
00162 
00163         bool doscale = isinshower && isintrack && fTrkMIPCalibOkay && fTrkMIPDcosz[planeno] >= fTrkMIPDcoszMean*1.2;
00164 
00165         const double scale = ((fTrkMIPDcosz[planeno]>0.3)? (fTrkMIPDcosz[planeno]-fTrkMIPDcoszMean)/fTrkMIPDcosz[planeno] : 1.);
00166 
00167         int modifydigit = 0;
00168 
00169         if(!isintrack || stripno<0 || doscale || (isintrack && notmuon && fTrkMIPCalibOkay)){ // we are going to keep this digit.         
00170           if(doscale) modifydigit = 1;
00171         }else{
00172           modifydigit = 2;
00173         }
00174 
00175         if(planeno < fMinShowerPlane-cShowerLikePlanes || planeno > fMaxShowerPlane+cShowerLikePlanes) modifydigit = 2;
00176 
00177         if(modifydigit==0){
00178           MSG("RmMu",Msg::kVerbose) << " KEEP : " << planeno << "/"<< stripno 
00179                                     << " Q: " << strip->GetCharge() 
00180                                     << " TrkMIP: " <<  fTrkMIPDcosz[planeno]
00181                                     << " TrkQ: " <<  fTrkQ[planeno]
00182                                     << " PlnQ: " <<  fPlnQ[planeno]
00183                                     << " intrk: " << isintrack 
00184                                     << " notmuon: " << notmuon 
00185                                     << " doscale: " << doscale 
00186                                     << " scale: " << scale <<endl;
00187         }else if(modifydigit==1){
00188           MSG("RmMu",Msg::kVerbose) << " SCALE : " << planeno << "/"<< stripno 
00189                                     << " Q: " << strip->GetCharge() 
00190                                     << " TrkMIP: " <<  fTrkMIPDcosz[planeno]
00191                                     << " TrkQ: " <<  fTrkQ[planeno]
00192                                     << " PlnQ: " <<  fPlnQ[planeno]
00193                                     << " intrk: " << isintrack 
00194                                     << " notmuon: " << notmuon 
00195                                     << " doscale: " << doscale 
00196                                     << " scale: " << scale <<endl;
00197 
00198         }else if(modifydigit==2){
00199           MSG("RmMu",Msg::kVerbose) << " REJECT : " << planeno << "/"<< stripno 
00200                                     << " Q: " << strip->GetCharge() 
00201                                     << " TrkMIP: " <<  fTrkMIPDcosz[planeno]
00202                                     << " TrkQ: " <<  fTrkQ[planeno]
00203                                     << " PlnQ: " <<  fPlnQ[planeno]
00204                                     << " intrk: " << isintrack 
00205                                     << " notmuon: " << notmuon 
00206                                     << " doscale: " << doscale 
00207                                     << " scale: " << scale <<endl;
00208 
00209         }
00210         
00211         if(modifydigit || isintrack){ //        
00212           TIter digitIter(strip->GetDaughterIterator());
00213           while( CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitIter()) ) {    
00214             if(modifydigit){
00215               digitidx.push_back(digit->GetRawDigitIndex());
00216               if(modifydigit==1) digitweight.push_back(scale);
00217               else digitweight.push_back(-1.);
00218             }
00219             if(isintrack) ntrackdigits++;
00220           }
00221         }
00222       }            
00223     }
00224   }
00225    
00226   MSG("RmMu",Msg::kDebug) << " There are " << digitidx.size()<< " digits to be modified  and " << ntrackdigits <<" track digits"  <<endl; 
00227 
00228   //iterate over base digits and do modifications
00229   TIter digitIter(digitlist->GetDaughterIterator());
00230   const unsigned int nmoddigit = digitidx.size();
00231   while( CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitIter()) ) {    
00232     const int thisidx = digit->GetRawDigitIndex(); 
00233     bool ditch = 0;
00234     bool scale = 0;
00235     float scale_factor = 0;
00236     for(unsigned int idig = 0; idig<nmoddigit; ++idig){
00237       if(digitidx[idig] == thisidx ){
00238         if(digitweight[idig]>0){
00239           scale = 1;
00240           scale_factor = digitweight[idig];
00241         }else{
00242           ditch = 1;
00243         }
00244         digitidx[idig] = -1;
00245       }
00246     }
00247     if(!ditch){
00248       AlgConfig& algconfig = pass_alg.GetAlgConfig();
00249       algconfig.UnLockKeys();
00250       algconfig.UnLockValues();
00251       if(scale){
00252         algconfig.Set("doscale", 1);
00253         algconfig.Set("scalefactor", scale_factor);
00254       }else{
00255         algconfig.Set("doscale", 0);
00256         algconfig.Set("scalefactor", 1.);
00257       }
00258       algconfig.LockKeys();
00259       algconfig.LockValues();
00260       TObjArray digitalg_input;
00261       digitalg_input.Add(digit);
00262       cxx.SetDataIn(&digitalg_input);
00263       CandDigitHandle cdh = CandDigit::MakeCandidate(pass_alg, cxx);
00264       ch.AddDaughterLink(cdh, kFALSE); 
00265     }
00266   }
00267   
00268   //
00269   //global list stuff
00270   //
00271   CandDigitListHandle &cdlh = dynamic_cast<CandDigitListHandle &>(ch);
00272   cdlh.SetAbsTime(digitlist->GetAbsTime());                   
00273   cdlh.SetIsSparse(digitlist->GetIsSparse());   
00274 }
00275 
00276 void AlgCosmicMuonRemoval::FillEvtInfo(const CandEventHandle* event,
00277                                  const CandTrackHandle* track, const CandShowerHandle* shower){
00278   //
00279   //fill info about track direction at each plane
00280   //
00281   const double nplanes  = TMath::Abs(track->GetEndPlane() - track->GetVtxPlane()) + 1;
00282   
00283   const int vtxplane = track->GetVtxPlane();
00284   const int endplane = track->GetEndPlane();
00285   
00286   const double vtxangle = track->GetVtxDirCosZ();
00287   const double endangle = track->GetEndDirCosZ(); 
00288 
00289   if (vtxplane < endplane) {
00290     fMinTrackPlane = vtxplane; 
00291     fMaxTrackPlane = endplane;
00292     fMinTrackAngle = vtxangle;
00293     fMaxTrackAngle = endangle;
00294   } else {
00295     fMinTrackPlane = endplane;
00296     fMaxTrackPlane = vtxplane;
00297     fMinTrackAngle = endangle;
00298     fMaxTrackAngle = vtxangle;
00299   }
00300 
00301   MSG("RmMu",Msg::kDebug) << " Track planes: " << fMinTrackPlane << "  " << fMaxTrackPlane << endl;
00302   MSG("RmMu",Msg::kDebug) << " Track angle: "  << fMinTrackAngle << "  " << fMaxTrackAngle << endl; 
00303  
00304   for(int ipln = fMinTrackPlane; ipln<=fMaxTrackPlane; ++ipln){
00305 
00306     double dcosz = fabs(fMinTrackAngle + ((ipln - fMinTrackPlane)/nplanes)*(fMaxTrackAngle - fMinTrackAngle));
00307  
00308     int ip1 =-1;
00309     int ip2 =-1;
00310     if(track->IsTPosValid(ipln-1) && track->IsTPosValid(ipln+1)){
00311       ip2 = ipln+1; ip1 = ipln-1;
00312     }else if(track->IsTPosValid(ipln) && track->IsTPosValid(ipln+1)){
00313       ip2 = ipln+1; ip1 = ipln;
00314     }else if(track->IsTPosValid(ipln) && track->IsTPosValid(ipln-1)){
00315       ip2 = ipln; ip1 = ipln-1;
00316     }
00317     if(ip1!=-1 && ip2!=-1){
00318       const double l = sqrt( (track->GetU(ip2) - track->GetU(ip1)) * 
00319                              (track->GetU(ip2) - track->GetU(ip1)) 
00320                            + (track->GetV(ip2) - track->GetV(ip1)) * 
00321                              (track->GetV(ip2) - track->GetV(ip1)) 
00322                            + (track->GetZ(ip2) - track->GetZ(ip1)) *
00323                              (track->GetZ(ip2) - track->GetZ(ip1)) ); 
00324       if(l!=0) dcosz = fabs( track->GetZ(ip2) - track->GetZ(ip1) ) / l;
00325       else dcosz = 1;
00326     }
00327     fTrkMIPDcosz[ipln] = track->GetPlaneCharge(ipln, CalStripType::kMIP)*dcosz;
00328     if(track->GetPlaneCharge(ipln, CalStripType::kMIP)<1e-3){
00329       fTrkMIPDcosz[ipln] =1.;
00330     }
00331   }
00332 
00333   TIter trkstpIter(track->GetDaughterIterator());
00334   int ntrkcalibstps = 0;
00335   while(const CandStripHandle* strip = dynamic_cast<const CandStripHandle*>(trkstpIter())){
00336     fTrkQ[strip->GetPlane()]+=strip->GetCharge(CalDigitType::kPE);
00337     if(track->GetStripCharge(strip, CalStripType::kMIP)!=0) ntrkcalibstps++; 
00338     fTrkStrip[strip->GetPlane()] = strip->GetStrip();
00339   }
00340   fTrkMIPCalibOkay = (ntrkcalibstps>4);
00341 
00342 // get the shower infomation
00343   
00344   const int shower_vtxplane = shower->GetVtxPlane();
00345   const int shower_endplane = shower->GetEndPlane();
00346  
00347   fMinShowerPlane = (shower_vtxplane<shower_endplane)?shower_vtxplane:shower_endplane;
00348   fMaxShowerPlane = (shower_vtxplane>shower_endplane)?shower_vtxplane:shower_endplane;
00349  
00350   MSG("RmMu",Msg::kDebug) << " Shower planes: " << fMinShowerPlane << "  " << fMaxShowerPlane << endl;
00351 
00352   // change the min and max track plane for the track
00353   
00354   if (fMinTrackPlane<fMinShowerPlane-cTrackLikePlanes) fMinTrackPlane = fMinShowerPlane - cTrackLikePlanes;
00355   if (fMaxTrackPlane>fMaxShowerPlane+cTrackLikePlanes) fMaxTrackPlane = fMaxShowerPlane + cTrackLikePlanes;
00356 
00357   Float_t total_mips = 0;
00358   int mip_count_planes = 0;
00359   
00360   if ((fMinShowerPlane<=fMinTrackPlane)&&(fMaxShowerPlane>=fMaxTrackPlane)) { // shower is longer than track on both end
00361     for (int ipln = fMinTrackPlane; ipln<=fMaxTrackPlane; ++ipln) {
00362       mip_count_planes++;
00363       total_mips += 1.0;  
00364     }
00365 
00366     MSG("RmMu",Msg::kDebug) << "  Shower is longer than track on both sides ... " 
00367                             << "  planes : " << mip_count_planes << " mips :" << total_mips << endl; 
00368         
00369   } else if (fMinShowerPlane<=fMinTrackPlane) {  // shower start before track
00370     for (int ipln = fMaxShowerPlane+1; ipln<=fMaxTrackPlane; ++ipln) {
00371       mip_count_planes++;
00372       total_mips += fTrkMIPDcosz[ipln];
00373     }
00374 
00375     MSG("RmMu",Msg::kDebug) << "  Shower starts before track ... "
00376                             << "  planes : " << mip_count_planes << " mips :" << total_mips << endl;
00377         
00378   } else if (fMaxShowerPlane>=fMaxTrackPlane) {  // shower end after track
00379     for (int ipln = fMinTrackPlane; ipln<fMinShowerPlane; ++ipln) {
00380       mip_count_planes++;
00381       total_mips += fTrkMIPDcosz[ipln];
00382     }
00383   
00384     MSG("RmMu",Msg::kDebug) << "  Shower ends after track ... "
00385                             << "  planes : " << mip_count_planes << " mips :" << total_mips << endl;
00386         
00387   } else { // shower is inside of track planes
00388     for (int ipln = fMinTrackPlane; ipln<fMinShowerPlane; ++ipln) { 
00389       mip_count_planes++;
00390       total_mips += fTrkMIPDcosz[ipln];
00391     }   
00392 
00393     for (int ipln = fMaxShowerPlane+1; ipln<=fMaxTrackPlane; ++ipln) {
00394       mip_count_planes++;
00395       total_mips += fTrkMIPDcosz[ipln];
00396     }
00397     
00398     MSG("RmMu",Msg::kDebug) << "  Shower is inside of track ... "
00399                             << "  planes : " << mip_count_planes << " mips :" << total_mips << endl;
00400   }
00401 
00402   if(mip_count_planes!=0) fTrkMIPDcoszMean = (total_mips/mip_count_planes);
00403   
00404   if(fTrkMIPDcoszMean<0.3) {
00405     MSG("RmMu",Msg::kWarning) << "total_mips, planes, fTrkMIPDcoszMean : " << total_mips << " , " << mip_count_planes << " , " << fTrkMIPDcoszMean << endl;
00406     MSG("RmMu",Msg::kWarning) << "fTrkMIPDcoszMean = " << fTrkMIPDcoszMean << ", which is less than 0.3, set to be 1.0." << endl;
00407     fTrkMIPDcoszMean = 1.0;     
00408   }
00409 
00410   TIter stpIter(event->GetDaughterIterator());
00411   while(const CandStripHandle* strip = dynamic_cast<const CandStripHandle*>(stpIter())){
00412     const int plane = strip->GetPlane();
00413     fPlnQ[plane]+=strip->GetCharge(CalDigitType::kPE);
00414     if(plane>0 && plane<500) {
00415       fNTotalStrip[plane]++; 
00416       if (TMath::Abs(strip->GetStrip()-fTrkStrip[plane])<=1){
00417         fNTrackStrip[plane]++;
00418       }
00419     }
00420   }
00421 }
00422 
00423 void AlgCosmicMuonRemoval::Trace(const char * /* c */) const{
00424 }