loon: RemoveMuon.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 #include "MuonRemoval/RemoveMuon.h"
00008 #include "MessageService/MsgService.h"
00009 #include "MinosObjectMap/MomNavigator.h"
00010 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00011 
00012 #include <CandData/CandHeader.h>
00013 #include "CandData/CandRecord.h"
00014 #include <RecoBase/CandTrackListHandle.h>
00015 #include <RecoBase/CandFitTrackHandle.h>
00016 #include <RecoBase/CandTrackHandle.h>
00017 #include <RecoBase/CandShowerHandle.h>
00018 #include <RecoBase/CandShowerListHandle.h>
00019 #include <RecoBase/CandStripHandle.h>
00020 #include <RecoBase/CandStripListHandle.h>
00021 #include <RecoBase/CandEventListHandle.h>
00022 #include <RecoBase/CandEventHandle.h>
00023 #include "CandDigit/CandDigitList.h"
00024 #include "CandDigit/CandDigitListHandle.h"
00025 
00026 #include "Algorithm/AlgConfig.h"
00027 #include "Algorithm/AlgFactory.h"
00028 #include "Algorithm/AlgHandle.h"
00029 
00030 #include "Candidate/CandContext.h"
00031 
00032 #include "TFile.h"
00033 #include "TTree.h"
00034 #include "TObjArray.h"
00035 
00036 #include "MuonRemoval/SelectEvent.h"
00037 #include "MuonRemoval/CandRmMuListHandle.h"
00038 #include "MuonRemoval/CandRmMuList.h"
00039 
00040 #include <cmath>
00041 
00042 JOBMODULE(RemoveMuon, "RemoveMuon",
00043           "removes the muon");
00044 CVSID("$Id: RemoveMuon.cxx,v 1.5 2007/10/22 18:25:22 boehm Exp $");
00045 //......................................................................
00046 
00047 RemoveMuon::RemoveMuon() : 
00048   fAlgName(""), 
00049   fAlgConfig(""), 
00050   fNameListOut(""),
00051   fAlgOutputName(""), 
00052   fAlgOutputConfig(""), 
00053   fMakeNtp(0),
00054   fNtpFileName("muon_removal.root"), fNtpTreeName("ms"),
00055   fTrackListName("CandFitTrackCamList"), 
00056   fNtpFile(NULL), fNtpTree(NULL), 
00057   fNtpRun(0), fNtpSnarl(0),
00058   fNtpVtxX(0), fNtpVtxY(0), fNtpVtxZ(0), fNtpVtxP(0), fNtpNPln(0),
00059   fNtpPRng(0), fNtpPCrv(0), fNtpVDrX(0), fNtpVDrY(0), fNtpVDrZ(0),
00060   fNtpFitP(0), fNtpEndC(0),
00061   fNtpPass(0), fNtpPMuX(0), fNtpPMuY(0), fNtpPMuZ(0) 
00062 {
00063 }
00064 //......................................................................
00065 
00066 RemoveMuon::~RemoveMuon()
00067 {
00068 }
00069 
00070 //......................................................................
00071 
00072 void RemoveMuon::BeginJob()
00073 {
00074   if(fMakeNtp){
00075     TDirectory* current_directory = gDirectory;
00076     fNtpFile = TFile::Open(fNtpFileName.c_str(), "recreate");
00077     fNtpFile->cd();
00078     fNtpTree = new TTree(fNtpTreeName.c_str(), "");
00079     fNtpTree->SetDirectory(fNtpFile);
00080     fNtpTree->Branch("run", &fNtpRun, "run/I");
00081     fNtpTree->Branch("snarl", &fNtpSnarl, "snarl/I");
00082     fNtpTree->Branch("vtxx", &fNtpVtxX, "vtxx/F");
00083     fNtpTree->Branch("vtxy", &fNtpVtxY, "vtxy/F");
00084     fNtpTree->Branch("vtxz", &fNtpVtxZ, "vtxz/F");
00085     fNtpTree->Branch("vtxp", &fNtpVtxP, "vtxp/I");
00086     fNtpTree->Branch("npln", &fNtpNPln, "npln/I");
00087     fNtpTree->Branch("prng", &fNtpPRng, "prng/F");
00088     fNtpTree->Branch("pcrv", &fNtpPCrv, "pcrv/F");
00089     fNtpTree->Branch("pvdx", &fNtpVDrX, "pvdx/F");
00090     fNtpTree->Branch("pvdy", &fNtpVDrY, "pvdy/F");
00091     fNtpTree->Branch("pvdz", &fNtpVDrZ, "pvdz/F");
00092     fNtpTree->Branch("fitp", &fNtpFitP, "fitp/I");
00093     fNtpTree->Branch("endc", &fNtpEndC, "endc/I");
00094     fNtpTree->Branch("pass", &fNtpPass, "pass/I");
00095     fNtpTree->Branch("pmux", &fNtpPMuX, "pmux/F");
00096     fNtpTree->Branch("pmuy", &fNtpPMuY, "pmuy/F");
00097     fNtpTree->Branch("pmuz", &fNtpPMuZ, "pmuz/F");
00098 
00099     current_directory->cd();
00100   }
00101 }
00102 
00103 //......................................................................
00104 
00105 void RemoveMuon::EndJob()
00106 {
00107   if(fMakeNtp){
00108     TDirectory* current_directory = gDirectory;
00109     fNtpFile->cd();
00110     fNtpTree->Write();
00111     current_directory->cd();
00112   }  
00113 }
00114 
00115 //......................................................................
00116 
00117 void RemoveMuon::BeginFile()
00118 {
00119 }
00120 
00121 //......................................................................
00122 
00123 void RemoveMuon::EndFile()
00124 {
00125 }
00126 
00127 //......................................................................
00128 
00129 JobCResult RemoveMuon::Ana(const MomNavigator* mom)
00130 {
00135   JobCResult result(JobCResult::kPassed);
00136   
00137   CandRecord* record = dynamic_cast<CandRecord*>(mom->GetFragment("CandRecord"));
00138   if(!record){
00139     MSG("RmMu",Msg::kError) << " Unable to find CandRecord in mom !!! " << endl;
00140     result.SetFailed();
00141     return result;
00142   }
00143   const CandHeader* header  = dynamic_cast<const CandHeader*>(record->GetHeader());
00144   
00145   CandEventListHandle * eventlist = dynamic_cast<CandEventListHandle*>
00146     (record->FindCandHandle("CandEventListHandle"));
00147   if(!eventlist){
00148     MSG("RmMun",Msg::kDebug) << " Rejecting event as it has no events " << endl;
00149     result.SetFailed();
00150     return result;
00151   }
00152 
00153   //bool has_good_track = 0;
00154 
00155   TIter event_iter(eventlist->GetDaughterIterator());
00156   while( const CandEventHandle* event = dynamic_cast<const CandEventHandle*>(event_iter()) ){
00157     cout << " Event has " << event->GetLastTrack()+1 << " tracks " <<endl;
00158     if(SelectEvent(event)){      
00159       if(fMakeNtp){
00160         const CandTrackHandle* track_tmp = GetRemovableTrack(event);
00161         const CandFitTrackHandle* track  = dynamic_cast<const CandFitTrackHandle*>(track_tmp);
00162         if(track==NULL){
00163           MSG("RmMu",Msg::kError) << " GetRemovableTrack failed to return a good track " << endl;
00164           continue;
00165         }
00166         fNtpRun = header->GetRun();
00167         fNtpSnarl = header->GetSnarl();
00168         fNtpVtxX = .707106781*(track->GetVtxU() - track->GetVtxV());
00169         fNtpVtxY = .707106781*(track->GetVtxU() + track->GetVtxV());
00170         fNtpVtxZ = track->GetVtxZ();
00171         fNtpVtxP = track->GetVtxPlane();
00172         fNtpNPln = abs(track->GetEndPlane() - track->GetVtxPlane()) +1;
00173         fNtpPRng = track->GetMomentumRange();
00174         fNtpPCrv = track->GetMomentumCurve();
00175         fNtpVDrX = 0.707106781*(track->GetVtxDirCosU() - track->GetVtxDirCosV()) ;
00176         fNtpVDrY= 0.707106781*(track->GetVtxDirCosU() + track->GetVtxDirCosV()) ;
00177         fNtpVDrZ =  track->GetVtxDirCosZ();      
00178         fNtpFitP =  track->GetPass();
00179         fNtpEndC = TrackEndContained(track);
00180         fNtpPass = 0;
00181         Float_t momentum = fNtpPRng;
00182         fNtpPass = 1;
00183         if(!fNtpEndC){
00184           fNtpPass =  fNtpFitP;
00185           momentum = fNtpPCrv;
00186         }
00187         fNtpPMuX = fNtpVDrX*fabs(momentum);
00188         fNtpPMuY = fNtpVDrY*fabs(momentum);
00189         fNtpPMuZ = fNtpVDrZ*fabs(momentum);      
00190         fNtpTree->Fill();
00191       }    
00192     }else{    
00193       result.SetFailed();
00194     }
00195   } //iterate over events.
00196   return result;
00197 
00198 }
00199 
00200 //......................................................................
00201 
00202 JobCResult RemoveMuon::Reco(MomNavigator* mom)
00203 {
00204   
00205   //
00206   //This function trys to remove the muon from the event and writes out a new cand digit list 
00207   //  NB: it does no event selection, i.e. it doesnt decide which events are suitable for 
00208   //      event removal.  To run the event selection, run the Ana function befor ethe reco funtion.
00209   //
00210   
00211   JobCResult result(JobCResult::kPassed);
00212   MSG("RmMu", Msg::kDebug) << " Starting RemoveMuon::Reco() " <<endl; 
00213   MSG("RmMu", Msg::kDebug) << "  Alg        : " << fAlgName<< endl;
00214   MSG("RmMu", Msg::kDebug) << "  AlgConfig  : " << fAlgConfig<< endl;
00215   MSG("RmMu", Msg::kDebug) << "  NameListOut: " << fNameListOut<< endl;
00216   MSG("RmMu", Msg::kDebug) << "  TrackListName: " << fTrackListName<< endl;
00217     
00218   //
00219   // Find PrimaryCandidateRecord fragment in MOM.
00220   //
00221   CandRecord *record = dynamic_cast<CandRecord *>(mom->GetFragment("CandRecord", 
00222                                                                    "PrimaryCandidateRecord"));
00223   if (record == 0) {
00224     MSG("RmMu", Msg::kError) << "No PrimaryCandidateRecord in MOM."<< endl;
00225     result.SetFatal().SetFailed();
00226     return result;
00227   }
00228 
00229   //
00230   //Remove the muon 
00231   //
00232   CandContext cx(NULL, mom);  
00233   TObjArray alg_input;
00234   alg_input.Add(record);
00235   cx.SetDataIn(&alg_input);
00236   cx.SetCandRecord(record);
00237   AlgHandle algorithm = AlgFactory::GetInstance().GetAlgHandle(fAlgName.c_str(), 
00238                                                                fAlgConfig.c_str());
00239   AlgConfig& algconfig = algorithm.GetAlgConfig();
00240   Int_t cMaxTrackLikePlanes = 0;
00241   if(!algconfig.Get("cMaxTrackLikePlanes",cMaxTrackLikePlanes)) 
00242     cMaxTrackLikePlanes = 6;
00243   
00244   CandDigitListHandle digitlist = CandDigitList::MakeCandidate(algorithm, cx);
00245   digitlist.SetName(fNameListOut.c_str());
00246   digitlist.SetTitle(fNameListOut.c_str());  
00247   record->SecureCandHandle(digitlist);
00248 
00249   // make a candidate to hold the removed muon info 
00250   // along with handles to the remaining digits
00251   CandEventListHandle * eventlist = dynamic_cast<CandEventListHandle*>
00252     (record->FindCandHandle("CandEventListHandle"));
00253   CandTrackListHandle * tracklist = dynamic_cast<CandTrackListHandle*>
00254     (record->FindCandHandle("CandTrackListHandle", fTrackListName.c_str()));
00255   CandDigitListHandle * digitlist2 = dynamic_cast<CandDigitListHandle*>
00256     (record->FindCandHandle("CandDigitListHandle"));
00257 
00258   CandContext cx2(NULL, mom);
00259   TObjArray alg_input2;
00260   alg_input2.Add(digitlist2);
00261   alg_input2.Add(eventlist);
00262   alg_input2.Add(tracklist);
00263   cx2.SetDataIn(&alg_input2);
00264   cx2.SetCandRecord(record);
00265   AlgHandle algorithm2 = AlgFactory::GetInstance().
00266     GetAlgHandle(fAlgOutputName.c_str(),fAlgOutputConfig.c_str());
00267   AlgConfig& algconfig2 = algorithm2.GetAlgConfig();
00268   algconfig2.UnLockKeys();
00269   algconfig2.UnLockValues();
00270   algconfig2.Set("MaxTrackLikePlanes",cMaxTrackLikePlanes);
00271   algconfig2.LockKeys();
00272   algconfig2.LockValues();
00273   CandRmMuListHandle rmmulist = CandRmMuList::MakeCandidate(algorithm2, cx2);
00274   rmmulist.SetName("candrmmulist");
00275   rmmulist.SetTitle("candrmmulist");
00276   record->SecureCandHandle(rmmulist);
00277   
00278   MSG("RmMu", Msg::kDebug) << " Leaving RemoveMuon::Reco() " <<endl; 
00279   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00280 }
00281 
00282 //......................................................................
00283 
00284 const Registry& RemoveMuon::DefaultConfig() const
00285 {
00289   static Registry r; // Default configuration for module
00290 
00291   // Set name of config
00292   std::string name = this->GetName();
00293   name += ".config.default";
00294   r.SetName(name.c_str());
00295 
00296   // Set values in configuration
00297   r.UnLockValues();
00298   r.Set("AlgName","AlgMuonRemoval");
00299   r.Set("AlgConfig","default");
00300   r.Set("StripListOut", "stripdigitlist");
00301   r.Set("AlgOutputName","AlgRmMuList");
00302   r.Set("AlgOutputConfig","default");
00303   r.Set("MakeNtp", 0);
00304   r.Set("NtpFileName", "muon_summary.root");
00305   r.Set("NtpTreeName", "muonparams");
00306 
00307   r.LockValues();
00308   
00309   return r;
00310 }
00311 
00312 //......................................................................
00313 
00314 void RemoveMuon::Config(const Registry& r)
00315 {
00319   int tmpi = 0;
00320   
00321   fAlgName = r.GetCharString("AlgName");
00322   fAlgConfig = r.GetCharString("AlgConfig");
00323   fNameListOut = r.GetCharString("StripListOut");
00324   fAlgOutputName = r.GetCharString("AlgOutputName");
00325   fAlgOutputConfig = r.GetCharString("AlgOutputConfig");
00326   if(r.Get("MakeNtp", tmpi)) fMakeNtp = tmpi;
00327   fNtpFileName = r.GetCharString("NtpFileName");
00328   fNtpTreeName = r.GetCharString("NtpTreeName");
00329   if(r.Get("TrackListName", tmpi)) fTrackListName = tmpi;
00330   
00331 }
00332