MergeEvent Class Reference

#include <MergeEvent.h>

Inheritance diagram for MergeEvent:
JobCModule

List of all members.

Public Member Functions

 MergeEvent ()
 ~MergeEvent ()
void EndJob ()
JobCResult Ana (const MomNavigator *mom)
JobCResult Reco (MomNavigator *mom)
void InfileSetUp ()
void OutfileSetUp ()
const RegistryDefaultConfig () const
void Config (const Registry &r)

Private Attributes

std::string fInputFileName
TFile * fInputFile
TTree * fInputTree
Int_t fInputEventNumber
TClonesArray * fInRawDigits
Float_t fInputP4 [4]
Float_t fInputMRM [6]
std::string fOutputFileName
TFile * fOutputFile
TTree * fOutputTree
Int_t fOutputEventNumber
TClonesArray * fOutRawDigits
Float_t fOutputP4 [4]
Float_t fOutputMRM [6]
std::string fMergeAlg
std::string fMergeConfig
std::string fMergeListOut
Int_t fUseTrkForTiming

Detailed Description

Id
MergeEvent.h,v 1.3 2007/11/17 08:11:39 tjyang Exp

(Document me!)

Definition at line 16 of file MergeEvent.h.


Constructor & Destructor Documentation

MergeEvent::MergeEvent (  ) 

Definition at line 66 of file MergeEvent.cxx.

References fInputMRM, fInputP4, fOutputMRM, and fOutputP4.

00066                        :
00067   fInputFileName(""), fInputFile(NULL), fInputTree(NULL),
00068   fInputEventNumber(0), fInRawDigits(NULL),
00069   fOutputFileName(""), fOutputFile(NULL), fOutputTree(NULL),
00070   fOutputEventNumber(0), fOutRawDigits(NULL),
00071   fMergeAlg(""), fMergeConfig(""), fMergeListOut(""), fUseTrkForTiming(1)
00072 {
00073   fInputP4[0] = 0; fInputP4[1] = 0; fInputP4[2] = 0; fInputP4[3] = 0;  
00074   fOutputP4[0] = 0; fOutputP4[1] = 0; fOutputP4[2] = 0; fOutputP4[3] = 0;  
00075   for (int i = 0; i<6; i++){
00076     fInputMRM[i] = 0;
00077     fOutputMRM[i] = 0;
00078   }
00079 }

MergeEvent::~MergeEvent (  ) 

Definition at line 99 of file MergeEvent.cxx.

References fInRawDigits, and fOutRawDigits.

00100 {
00101   if(fInRawDigits) delete fInRawDigits;
00102   if(fOutRawDigits) delete fOutRawDigits;
00103 }


Member Function Documentation

JobCResult MergeEvent::Ana ( const MomNavigator mom  )  [virtual]

Implement this for read only access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 142 of file MergeEvent.cxx.

References RecDataRecord< T >::FindComponent(), RawRecord::FindRawBlock(), fOutputEventNumber, fOutputMRM, fOutputP4, fOutputTree, fOutRawDigits, RawDigit::GetADC(), RawDigit::GetChannel(), RawDigit::GetCrateT0(), RawDigitDataBlock::GetDatumIter(), MomNavigator::GetFragment(), RawRecord::GetRawHeader(), RawDaqHeader::GetRun(), RawDaqSnarlHeader::GetSnarl(), RawDigit::GetTDC(), it, JobCResult::kPassed, Msg::kWarning, MSG, OutfileSetUp(), REROOT_NeuKin::P4Shw(), JobCResult::SetFailed(), and JobCResult::SetWarning().

00143 {
00144   //set up output tree if first time:
00145   if(!fOutputTree) this->OutfileSetUp();
00146   
00147   JobCResult result(JobCResult::kPassed);
00148   
00149   // make sure raw record is available
00150   RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
00151   if (rr == 0) {
00152     MSG("EventSR", Msg::kWarning) << "No RawRecord in MOM." << endl;
00153     return result;
00154   }
00155 
00156   const RawDaqSnarlHeader* snarlHdr =
00157     dynamic_cast<const RawDaqSnarlHeader*>(rr->GetRawHeader());
00158   if (snarlHdr) {
00159     if (snarlHdr->GetSnarl()%1000==0) cout<<"Run "<<snarlHdr->GetRun()<<" Snarl "<<snarlHdr->GetSnarl()<<endl;
00160   }
00161   else cout<<"No snarl info"<<endl;  
00162 
00163   //fill raw digit array
00164   fOutRawDigits->Clear();
00165   TClonesArray &rawdigitl = *fOutRawDigits;
00166   const RawDigitDataBlock *rddb=dynamic_cast<const RawDigitDataBlock *>
00167     (rr->FindRawBlock("RawDigitDataBlock"));
00168   Int_t rawDigitCounter = 0;
00169   if (rddb) {
00170     TIter it=rddb->GetDatumIter();
00171     while (TObject *obj=it.Next()) {
00172       RawDigit *rd=dynamic_cast<RawDigit *>(obj);
00173       new(rawdigitl[rawDigitCounter]) RawDigitInfo(rd->GetChannel(),
00174                                                    rd->GetADC(),
00175                                                    rd->GetTDC(),
00176                                                    rd->GetCrateT0());
00177       rawDigitCounter++;
00178     }
00179     rawdigitl.Compress();
00180   }
00181 
00182   //get truth:
00183   SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord*>
00184     (mom->GetFragment("SimSnarlRecord"));
00185   if ( !simrec ) {
00186     MSG("NtpMC",Msg::kWarning) << "No SimSnarlRecord in Mom" << endl;
00187     result.SetWarning().SetFailed();
00188     return result;
00189   }
00190   const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>
00191     (simrec->FindComponent("TClonesArray","NeuKinList"));
00192   REROOT_NeuKin* rneukin = dynamic_cast<REROOT_NeuKin*>(neukinarray->At(0));
00193   if (rneukin){
00194     for(int i=0;i<4;i++) fOutputP4[i] = rneukin->P4Shw()[i];
00195   }
00196   else {
00197     const TClonesArray* parts = 
00198       (const TClonesArray*)simrec->FindComponent("TClonesArray","StdHep");
00199     
00200     TIter next(parts);
00201     TObject* obj = 0;
00202     while ( ( obj = next() ) ) {
00203       TParticle* part = (TParticle*)obj;
00204       if (part->GetStatusCode()==1){
00205         fOutputMRM[0] = part->Px();
00206         fOutputMRM[1] = part->Py();
00207         fOutputMRM[2] = part->Pz();
00208         fOutputMRM[3] = part->Energy();
00209       }
00210       if (part->GetStatusCode()==800){
00211         fOutputP4[0] = part->Px();
00212         fOutputP4[1] = part->Py();
00213         fOutputP4[2] = part->Pz();
00214         fOutputP4[3] = part->Energy();
00215         fOutputMRM[4] = part->Vx();
00216         fOutputMRM[5] = part->Vy();
00217         fOutputTree->Fill();
00218         fOutputEventNumber++;
00219         return result;
00220       }
00221     }
00222   }
00223   //fill tree
00224   fOutputTree->Fill();
00225   fOutputEventNumber++;
00226   return result;
00227 }

void MergeEvent::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Reimplemented from JobCModule.

Definition at line 723 of file MergeEvent.cxx.

References fInputFileName, fMergeAlg, fMergeConfig, fMergeListOut, fOutputFileName, fUseTrkForTiming, Registry::GetCharString(), Registry::GetInt(), Msg::kInfo, and MSG.

00724 {
00725   fInputFileName  = r.GetCharString("InputFileName");
00726   fOutputFileName = r.GetCharString("OutputFileName");
00727   fMergeAlg       = r.GetCharString("MergeAlg");
00728   fMergeConfig    = r.GetCharString("MergeConfig");
00729   fMergeListOut   = r.GetCharString("MergeListOut");
00730   fUseTrkForTiming   = r.GetInt("UseTrkForTiming");
00731   MSG("EvtMrg",Msg::kInfo) << " **** " << fInputFileName << endl;
00732   MSG("EvtMrg",Msg::kInfo) << " **** " << fOutputFileName << endl;
00733   MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeAlg << endl;
00734   MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeConfig << endl;
00735   MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeListOut << endl;
00736   if (fUseTrkForTiming) {
00737     MSG("EvtMrg",Msg::kInfo) << " **** Using Trk For Timing" << endl;
00738   }
00739 }

const Registry & MergeEvent::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Reimplemented from JobCModule.

Definition at line 701 of file MergeEvent.cxx.

References JobCModule::GetName(), Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00702 {
00703   // Default configuration for module
00704   static Registry r;
00705   // Set name of config
00706   std::string name = this->GetName();
00707   name += ".config.default";
00708   r.SetName(name.c_str());
00709   // Set values in configuration
00710   r.UnLockValues();
00711   // Config for merging events
00712   r.Set("InputFileName","input_events.root");
00713   r.Set("OutputFileName","output_events.root");
00714   r.Set("MergeAlg","AlgMergeEvent");
00715   r.Set("MergeConfig","default");
00716   r.Set("MergeListOut","mergedigitlist");  
00717   r.Set("UseTrkForTiming",1);
00718   r.LockValues();
00719   return r;
00720 }

void MergeEvent::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 81 of file MergeEvent.cxx.

References fInputEventNumber, fInputFile, fInputFileName, fOutputEventNumber, fOutputFile, fOutputFileName, fOutputTree, Msg::kInfo, and MSG.

00082 {
00083   if(fOutputEventNumber>0) {
00084     fOutputFile->cd();
00085     fOutputTree->Write();
00086     fOutputFile->Close();
00087     MSG("EvtMrg",Msg::kInfo) << " Wrote digits from "<< fOutputEventNumber 
00088                              << " electron events to summary file: "
00089                              << fOutputFileName.c_str() << endl;
00090   }
00091   if(fInputEventNumber>0) {
00092     fInputFile->Close();
00093     MSG("EvtMrg",Msg::kInfo) << " Merged digits from "<< fInputEventNumber 
00094                              << " electron events from summary file: "
00095                              << fInputFileName.c_str() << endl;
00096   }
00097 }

void MergeEvent::InfileSetUp (  ) 

Definition at line 119 of file MergeEvent.cxx.

References fInputFile, fInputFileName, fInputMRM, fInputP4, fInputTree, fInRawDigits, Msg::kDebug, Msg::kFatal, and MSG.

Referenced by Reco().

00120 {
00121   TDirectory* current_directory = gDirectory;
00122   fInputFile = TFile::Open(fInputFileName.c_str(),"READ");
00123   if(!fInputFile || !(fInputFile->IsOpen()) ){
00124     MSG("EvtMrg",Msg::kFatal) << " Unable to open: "<< fInputFileName << endl;
00125     return;
00126   }
00127   fInputTree = dynamic_cast<TTree*>(fInputFile->Get("ElectronTree"));
00128   if(!fInputTree){
00129     MSG("EvtMrg",Msg::kFatal) << " Unable to get tree : ElectronTree from  "
00130                               << fInputFileName << endl;
00131     return;
00132   }
00133   fInRawDigits = new TClonesArray("RawDigitInfo",1000);
00134   fInputTree->SetBranchAddress("true_p4",fInputP4);
00135   if (fInputTree->GetBranchStatus("mrminfo")) fInputTree->SetBranchAddress("mrminfo",fInputMRM);
00136   fInputTree->SetBranchAddress("rawdigits",&fInRawDigits);
00137   current_directory->cd();
00138   MSG("EvtMrg",Msg::kDebug) << " Opened input file: " << fInputFileName <<endl;
00139 }

void MergeEvent::OutfileSetUp (  ) 

Definition at line 106 of file MergeEvent.cxx.

References fOutputFile, fOutputFileName, fOutputMRM, fOutputP4, fOutputTree, fOutRawDigits, Msg::kInfo, and MSG.

Referenced by Ana().

00107 {
00108   fOutputFile = new TFile(fOutputFileName.c_str(),"RECREATE");  
00109   fOutputTree = new TTree("ElectronTree","Electron Tree");
00110   fOutputTree->Branch("true_p4",&fOutputP4,"true_p4[4]/F",32000);
00111   fOutputTree->Branch("mrminfo",&fOutputMRM,"mrminfo[6]/F",32000);
00112   fOutRawDigits = new TClonesArray("RawDigitInfo",1000);
00113   fOutputTree->Branch("rawdigits","TClonesArray",&fOutRawDigits,8000,1);  
00114   MSG("EvtMrg",Msg::kInfo) << " Opened summary file "<< fOutputFileName
00115                            << " for writing electron events" << endl;
00116 }

JobCResult MergeEvent::Reco ( MomNavigator mom  )  [virtual]

Implement this for read-write access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 231 of file MergeEvent.cxx.

References RawDigitInfo::adc, RawDigitDataBlock::At(), det, digit(), MuELoss::e, CandRecord::FindCandHandle(), fInputEventNumber, fInputMRM, fInputP4, fInputTree, fInRawDigits, fMergeAlg, fMergeConfig, fMergeListOut, EnergyCorrections::FullyCorrectMomentumFromRange(), EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(), fUseTrkForTiming, RawDigit::GetADC(), AlgFactory::GetAlgHandle(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandRecoHandle::GetEndPlane(), MomNavigator::GetFragment(), RecMinos::GetHeader(), AlgFactory::GetInstance(), RawRecord::GetRawBlockIter(), GetRemovableTrack(), CandHeader::GetRun(), CandHeader::GetSnarl(), RawDigit::GetTDC(), Calibrator::GetTimeFromTDC(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), header, InfileSetUp(), Calibrator::Instance(), ReleaseType::kCedarPhyDaikon, Msg::kDebug, EnergyCorrections::kDefault, Msg::kError, Detector::kFar, Msg::kFatal, Msg::kInfo, JobCResult::kPassed, CalDigitType::kPE, RmMuMask::kRMMU_ISMCELEC_MASK, Msg::kVerbose, Msg::kWarning, CandDigitList::MakeCandidate(), MSG, RawDigitInfo::rcid, Calibrator::ReInitialise(), CandRecord::SecureCandHandle(), SelectEvent(), CandContext::SetCandRecord(), CandContext::SetDataIn(), JobCResult::SetFailed(), JobCResult::SetFatal(), CandHandle::SetName(), CandHandle::SetTitle(), JobCResult::SetWarning(), RawDigitInfo::t0, and RawDigitInfo::tdc.

00232 {
00233   JobCResult result(JobCResult::kPassed);
00234   MSG("EvtMrg", Msg::kDebug) << " Starting RemoveMuon::Reco() " <<endl;
00235   MSG("EvtMrg", Msg::kDebug) << "  Alg        : " << fMergeAlg<< endl;
00236   MSG("EvtMrg", Msg::kDebug) << "  AlgConfig  : " << fMergeConfig<< endl;
00237   MSG("EvtMrg", Msg::kDebug) << "  NameListOut: " << fMergeListOut<< endl;
00238 
00239   //set up input tree if first time:
00240   if(!fInputTree) this->InfileSetUp();
00241   
00242   // Find PrimaryCandidateRecord fragment in MOM.
00243   CandRecord *record = dynamic_cast<CandRecord *>(mom->GetFragment("CandRecord"));
00244   if (record == 0) {
00245     MSG("EvtMrg", Msg::kError) << "No PrimaryCandidateRecord in MOM."<< endl;
00246     result.SetWarning().SetFailed();
00247     return result;
00248   }
00249   const CandHeader* header  = dynamic_cast<const CandHeader*>(record->GetHeader());  
00250   
00251   //Find the raw data in MOM
00252   RawRecord *rawrecord = 
00253     dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord",0,"DaqSnarl"));
00254   assert(rawrecord);
00255   RawDigitDataBlock *input_datablock =  NULL;
00256   if(rawrecord){
00257     TIter rdbit = rawrecord->GetRawBlockIter();
00258     TObject *tob;
00259     while ((tob = rdbit())) {
00260       if (tob->InheritsFrom("RawDigitDataBlock")) {
00261         input_datablock = (RawDigitDataBlock *) tob;
00262       }
00263     }
00264   }
00265 
00266   CandEventListHandle * eventlist = 
00267     dynamic_cast<CandEventListHandle*>(record->FindCandHandle("CandEventListHandle"));
00268   if(eventlist==NULL) { 
00269     MSG("EvtMrg",Msg::kWarning) << " Bailing out of Event eventlist = " << eventlist 
00270                                 << endl; 
00271     result.SetWarning().SetFailed();
00272     return result;
00273   }
00274 
00275   //get muon removed digit list:
00276   CandDigitListHandle* digitlist = dynamic_cast<CandDigitListHandle*>
00277     (record->FindCandHandle("CandDigitListHandle","stripdigitlist"));
00278 
00279   //make a new TClonesArray to hold track time corrected 
00280   //RawDigitInfo objects for whole snarl
00281   TClonesArray *fRawDigits = new TClonesArray("RawDigitInfo",5000);
00282 
00283   //get rmmulist to check that input electrons match with removed muons
00284   CandRmMuListHandle *rmmulist = 
00285     dynamic_cast<CandRmMuListHandle*>(record->FindCandHandle("CandRmMuListHandle"));  
00286 
00287   //map to link MC electron RawDigits to RmMu objects
00288   std::map<Int_t,Int_t> linkRmMuRawDig;
00289 
00290   //Set up calibrator
00291   Calibrator& calibrator = Calibrator::Instance();
00292   calibrator.ReInitialise(*(record->GetVldContext()));
00293 
00294   Detector::Detector_t det = record->GetVldContext()->GetDetector();
00295 
00296   //loop over removed muons:
00297   Int_t rmmuIndex = 0;
00298   TIter rmmuItr(rmmulist->GetDaughterIterator());
00299   while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
00300     //read in first electron:
00301     if(fInputTree->GetEntry(fInputEventNumber++)<=0){
00302       MSG("EvtMrg", Msg::kFatal) << "Unable to read input event number: " 
00303                                  << fInputEventNumber-1 <<endl;
00304       result.SetFatal();
00305       return result;
00306     }
00307     MSG("EvtMrg", Msg::kDebug)  << " Cand:" << header->GetRun() << "/" 
00308                                 << header->GetSnarl() << "/" 
00309                                 << rmmu->GetOrigEvtIndex() << endl;
00310     MSG("EvtMrg",Msg::kDebug) << " Read New event P4: " 
00311                               << fInputP4[0] <<", " 
00312                               << fInputP4[1] <<", "
00313                               << fInputP4[2] <<", "
00314                               << fInputP4[3] <<endl;
00315 
00316     //Check that inputs and rmmu values agree
00317     const Float_t emass  = 0.000511; //GeV
00318     const Float_t mumass = 0.105658; //GeV
00319     
00320     Float_t etot = 0;
00321     Float_t pelec = 0;
00322 
00323     double rangemom = rmmu->GetMomRange();
00324     double curvemom = rmmu->GetMomCurv();
00325 
00326     const VldContext vc = *(record->GetVldContext());
00327     ReleaseType::Release_t release = ReleaseType::kCedarPhyDaikon;
00328 
00329     EnergyCorrections::WhichCorrection_t corrver = EnergyCorrections::kDefault;
00330     if (rangemom>0) rangemom=EnergyCorrections::FullyCorrectMomentumFromRange(rangemom,vc,release,corrver);
00331     curvemom = EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(curvemom,vc,release,corrver);
00332 
00333 
00334     if(rmmu->GetIsCont()==1){
00335       etot =  sqrt(rangemom*rangemom + pow(mumass,2));
00336       pelec  = sqrt(etot*etot - emass*emass);
00337     }
00338     else { //use range:
00339       etot = sqrt(pow(curvemom,2)+mumass*mumass);
00340       pelec  = sqrt(etot*etot - emass*emass);
00341     }
00342 
00343     if(etot > 100) {etot = emass; pelec = 0.0;}
00344     if(fInputP4[3] > 100) {
00345       fInputP4[0] = 0;
00346       fInputP4[1] = 0;
00347       fInputP4[2] = 0;
00348       fInputP4[3] = emass;
00349     } 
00350    float cut = 1e-3;
00351 
00352     if( TMath::Abs(rmmu->GetVtxDCosX()*pelec - fInputP4[0])> cut ||
00353         TMath::Abs(rmmu->GetVtxDCosY()*pelec - fInputP4[1])> cut ||
00354         TMath::Abs(rmmu->GetVtxDCosZ()*pelec - fInputP4[2])> cut ||
00355         TMath::Abs(etot - fInputP4[3])> cut ) {
00356       MSG("EvtMrg", Msg::kInfo) << " Input momentum: " << fInputP4[0] << " "
00357                                 << fInputP4[1] << " " << fInputP4[2] << " " 
00358                                 << fInputP4[3] << endl
00359                                 << " RmMu momentum: " << rmmu->GetVtxDCosX()*pelec 
00360                                 << " " << rmmu->GetVtxDCosY()*pelec << " " 
00361                                 << rmmu->GetVtxDCosZ()*pelec << " " << etot << endl;
00362       MSG("EvtMrg", Msg::kError)
00363         << " Input electron file out of synch with muon removal..." 
00364         << " look for a match nearby." << endl;
00365       
00366       Bool_t gotMatch = false;
00367       for(int ii=1;ii<50;ii++){
00368         for(int jj=-1;jj<=1;jj+=2){       
00369           if(fInputTree->GetEntry(fInputEventNumber+(ii*jj))) {
00370             if( TMath::Abs(rmmu->GetVtxDCosX()*pelec - fInputP4[0])< cut &&
00371                 TMath::Abs(rmmu->GetVtxDCosY()*pelec - fInputP4[1])< cut &&
00372                 TMath::Abs(rmmu->GetVtxDCosZ()*pelec - fInputP4[2])< cut &&
00373                 TMath::Abs(etot - fInputP4[3])<1e-4 ) {
00374               //got a match!
00375               fInputEventNumber+=(ii*jj)+1;
00376               gotMatch = true;
00377               MSG("EvtMrg", Msg::kInfo) 
00378                 << "Found an input momentum match; "
00379                 "ii = " << ii << " jj = " << jj << " "
00380                 << rmmu->GetVtxDCosX()*pelec - fInputP4[0] << " " 
00381                 << rmmu->GetVtxDCosY()*pelec - fInputP4[1] << " " 
00382                 << rmmu->GetVtxDCosZ()*pelec - fInputP4[2] << " " 
00383                 << etot - fInputP4[3] << endl;
00384               break;
00385             }
00386           }       
00387         }
00388         if(gotMatch) break;
00389       }
00390 
00391       if(!gotMatch) {
00392         MSG("EvtMrg", Msg::kError) << "Can't find a match... quitting" << endl;
00393         result.SetWarning().SetFailed();
00394         return result;
00395       }
00396     }
00397 
00398     //We have a rmmu entry and a set of electron hits that match
00399     //check if the MC electron/muon digit list is empty
00400     //This might not be a problem for electrons. But sometimes the muon event has no hits since the MC muon momentum is independent of the original muon momentum and can be very small. TY 11/16/07
00401     if (!fInRawDigits->GetEntries()){
00402       MSG("EvtMrg", Msg::kWarning) << "The lepton has no hits, don't add it to the muon removal event." << endl;
00403       continue;
00404     }
00405     
00406     //save MRM information: pmu, Q2 and Eshw
00407     rmmu->SetMRMX(fInputMRM[0]);
00408     rmmu->SetMRMY(fInputMRM[1]);
00409     rmmu->SetMRMZ(fInputMRM[2]);
00410     rmmu->SetMRMQ2(fInputMRM[4]);
00411     rmmu->SetMRMEshw(fInputMRM[5]);
00412 
00413     //now find original event and track to get the timing offsets
00414     TIter evtItr(eventlist->GetDaughterIterator());
00415     const CandEventHandle *event = NULL;
00416     for(int i=0;i<rmmu->GetOrigEvtIndex()+1;i++) 
00417       event = dynamic_cast<const CandEventHandle*>(evtItr());      
00418     if(!SelectEvent(event)) {
00419       MSG("EvtMrg", Msg::kError) << " This event was selected for muon removal " 
00420                                  << " Now it is not passing... quitting"
00421                                  << endl;
00422       result.SetWarning().SetFailed();
00423       return result;
00424     }
00425     const CandTrackHandle* track = GetRemovableTrack(event);
00426     if(!track) {
00427       MSG("EvtMrg", Msg::kError) << " Cannot find removable track from event" 
00428                                  << " ... quitting"
00429                                  << endl;
00430       result.SetWarning().SetFailed();
00431       return result;
00432     }
00433     //check that track matches rmmu:
00434     if(rmmu->GetVtxPlane()!=track->GetVtxPlane() ||
00435        rmmu->GetNPlane()!= TMath::Abs(track->GetEndPlane()-
00436                                       track->GetVtxPlane())+1) {
00437       MSG("EvtMrg", Msg::kError) << " Removable track properties do not match"
00438                                  << " rmmu properties... quitting"
00439                                  << endl;
00440       result.SetWarning().SetFailed();
00441       return result;
00442     }
00443 
00444     //now calculate mean offsets between track/event and MC electron:
00445     int tracktdc = 0;
00446     int tracktdc_ndigit = 0 ;
00447     int track_offset = -1;
00448     if(fUseTrkForTiming){
00449       MSG("EvtMrg", Msg::kDebug) << " Using Track for timing " << endl;
00450       TIter trkstpIter(track->GetDaughterIterator());
00451       while(const CandStripHandle* strip = 
00452             dynamic_cast<const CandStripHandle*>(trkstpIter())){
00453         TIter trkdigIter(strip->GetDaughterIterator());
00454         if(strip->GetPlane()<track->GetVtxPlane()+10){
00455           while(const CandDigitHandle* digit = 
00456                 dynamic_cast<const CandDigitHandle*>(trkdigIter())){
00457             const RawDigit* raw_digit = input_datablock->At(digit->GetRawDigitIndex());
00458             if(raw_digit){
00459               if(track_offset==-1) track_offset = raw_digit->GetTDC();    
00460               tracktdc += raw_digit->GetTDC() - track_offset;
00461               tracktdc_ndigit++;
00462               MSG("EvtMrg", Msg::kDebug) << " TDC:  " << raw_digit->GetTDC() 
00463                                          << " ADC:  " << raw_digit->GetADC() << endl;
00464             }
00465             else MSG("EvtMrg", Msg::kError) << "Digit with no raw digit"<<endl;   
00466           }
00467         }
00468       }
00469     }
00470     else {
00471       //Use all digts for finding the mean time of the event. 
00472       //  NB: You want to be careful that you dont have outliers screwing it up so 
00473       //      use a 5.0PE cut rather than the usual 3PE cut when doing timing 
00474       //      along tracks. 
00475       TIter digIter(digitlist->GetDaughterIterator());
00476       while(const CandDigitHandle* digit = 
00477             dynamic_cast<const CandDigitHandle*>(digIter())){
00478         if(digit->GetCharge(CalDigitType::kPE)>5.0){
00479           const RawDigit* raw_digit = input_datablock->At(digit->GetRawDigitIndex());   
00480           if(raw_digit){
00481             if(track_offset==-1) track_offset = raw_digit->GetTDC();      
00482             tracktdc+=raw_digit->GetTDC() - track_offset;
00483             tracktdc_ndigit++;
00484             MSG("EvtMrg", Msg::kDebug) << " TDC:  " << raw_digit->GetTDC() 
00485                                        << " ADC:  " << raw_digit->GetADC() << endl;
00486           }
00487           else {
00488             MSG("EvtMrg", Msg::kError) << "Digit with no raw digit"<<endl;
00489           }
00490         }
00491       }    
00492     }
00493 
00494     int meantracktdc = (track_offset!=-1)?((int)(round(((float)tracktdc)/ 
00495                                                        tracktdc_ndigit))+track_offset):0;
00496     MSG("EvtMrg", Msg::kDebug) << " Mean Track TDC:  " << meantracktdc << endl;
00497 
00498     
00499     //Calc mean and mode timing offset for the "MC" events.
00500     int rawtdc = 0;
00501     int rawtdc_ndigit = 0;
00502     TIter rawdigitIter(fInRawDigits);
00503     RawDigitInfo *tmp_digit = 0;
00504     std::map<Int_t,Int_t> freqTDC;
00505     while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))){
00506       rawtdc_ndigit++;
00507       rawtdc += (int)(tmp_digit->tdc);
00508       if(freqTDC.find(int(tmp_digit->tdc))!=freqTDC.end()) 
00509         freqTDC[int(tmp_digit->tdc)] += 1;
00510       else freqTDC[int(tmp_digit->tdc)] = 1;
00511     }
00512     int moderawtdc = 0;
00513     int maxrawtdc = 0;
00514     std::map<Int_t,Int_t>::iterator begTDC = freqTDC.begin();
00515     std::map<Int_t,Int_t>::iterator endTDC = freqTDC.end();
00516     while(begTDC!=endTDC){
00517       if(begTDC->second>maxrawtdc) {
00518         moderawtdc = begTDC->first;
00519         maxrawtdc = begTDC->second;
00520       }
00521       begTDC++;
00522     }
00523     int meanrawtdc = 0;
00524     if(rawtdc_ndigit>0) meanrawtdc = (int)(round(((float)rawtdc)/rawtdc_ndigit));
00525     MSG("EvtMrg", Msg::kDebug) << "MODE, MAX, MEAN: " << moderawtdc << " " 
00526                                << maxrawtdc << " " << meanrawtdc << endl;
00527 
00528     //remove large outliers which skew the mean:
00529     if(maxrawtdc>4) meanrawtdc = moderawtdc; //trust the mode
00530     rawtdc = rawtdc_ndigit = 0;
00531     rawdigitIter.Reset();
00532     while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))){   
00533       if(TMath::Abs(tmp_digit->tdc-meanrawtdc)<=5) { //ND: 5 buckets ~ 100ns
00534         rawtdc_ndigit++;
00535         rawtdc += (int)(tmp_digit->tdc);
00536         MSG("EvtMrg", Msg::kDebug) << " TDC:  " << tmp_digit->tdc 
00537                                    << " ADC:  " << tmp_digit->adc << endl;
00538       }
00539     }
00540     if(rawtdc_ndigit>0) meanrawtdc = (int)(round(((float)rawtdc)/rawtdc_ndigit));
00541     MSG("EvtMrg", Msg::kDebug) << " Mean Electron TDC:  " << meanrawtdc << endl;
00542 
00543     MSG("EvtMrg", Msg::kDebug)<< "Delta T muon and electron events: " << meantracktdc
00544                               <<" - "<< meanrawtdc <<" = " 
00545                               << meantracktdc - meanrawtdc <<endl;     
00546     
00547     //Now add RawDigitInfo objects to new TClonesArray for the snarl
00548     //and correct for the average delta time between the muon and electron
00549     rawdigitIter.Reset();
00550     TClonesArray &rawdigitl = *fRawDigits;
00551     Int_t rawDigitCounter = fRawDigits->GetEntries();
00552     tmp_digit = 0;
00553     while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))) {
00554       Double_t TDC = tmp_digit->tdc - meanrawtdc + meantracktdc;
00555       if(TDC<0) TDC = 0;
00556       linkRmMuRawDig[rawDigitCounter] = rmmuIndex;
00557       new(rawdigitl[rawDigitCounter++]) RawDigitInfo(tmp_digit->rcid,tmp_digit->adc,
00558                                                      TDC,tmp_digit->t0);
00559     }
00560     rmmuIndex += 1;
00561   }
00562   
00563   // Add back in the electrons
00564   MSG("EvtMrg", Msg::kDebug) << " Running Event Merging:  " << fMergeAlg 
00565                              << ":"<< fMergeConfig<< " to produce " 
00566                              << fMergeListOut <<endl;   
00567   
00568   CandContext merge_cx(NULL, mom);
00569   TObjArray merge_input;
00570   merge_input.Add(fRawDigits);
00571   merge_input.Add(digitlist);
00572   merge_input.Add(rawrecord);
00573   merge_cx.SetDataIn(&merge_input);
00574   merge_cx.SetCandRecord(record);
00575   AlgHandle merge_algorithm = 
00576     AlgFactory::GetInstance().GetAlgHandle(fMergeAlg.c_str(), fMergeConfig.c_str());
00577   CandDigitListHandle merge_digitlist = 
00578     CandDigitList::MakeCandidate(merge_algorithm, merge_cx);
00579   merge_digitlist.SetName(fMergeListOut.c_str());
00580   merge_digitlist.SetTitle(fMergeListOut.c_str());
00581   record->SecureCandHandle(merge_digitlist);
00582 
00583   //now need to replace daughters in rmmu objects with those in new digitlist
00584   CandDigitListHandle* mergedlist = dynamic_cast<CandDigitListHandle*>
00585     (record->FindCandHandle("CandDigitListHandle",fMergeListOut.c_str()));
00586   TIter mergedlistIter(mergedlist->GetDaughterIterator());
00587 
00588   rmmuItr.Reset();
00589   while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
00590     std::vector<CandDigitHandle*> mdigs(rmmu->GetNDaughters());
00591     std::vector<CandDigitHandle*> rmdigs(rmmu->GetNDaughters());
00592     std::vector<Int_t> rfk(rmmu->GetNDaughters());
00593     Int_t rmmuDigCounter = 0;
00594     TIter rmmuDigIter(rmmu->GetDaughterIterator());
00595     while(CandDigitHandle *rmmuDig = dynamic_cast<CandDigitHandle*>(rmmuDigIter())){      
00596       rfk[rmmuDigCounter] = rmmu->ReasonForKeeping(rmmuDig);
00597       mergedlistIter.Reset();
00598       while(CandDigitHandle *mergeDig = dynamic_cast<CandDigitHandle*>(mergedlistIter())){
00599         if( ( mergeDig->GetRawDigitIndex()==rmmuDig->GetRawDigitIndex() ) ){
00600           mdigs[rmmuDigCounter] = mergeDig;
00601           rmdigs[rmmuDigCounter] = rmmuDig;
00602           break;
00603         }
00604       }
00605       rmmuDigCounter += 1;
00606     }
00607     for(int i=0;i<rmmu->GetNDaughters();i++) {
00608       if(rmdigs[i]){
00609         rmmu->RemoveDaughter(rmdigs[i]);
00610         rmmu->AddDaughterLink(*mdigs[i]);
00611         rmmu->SetReasonForKeeping(mdigs[i],rfk[i]);
00612       }
00613     }
00614   }
00615 
00616   //now add electron digits to RmMu objects, and set ReasonForKeeping IsMCElec flag
00617   //loop over new raw digits:
00618   TIter rawDigIter(fRawDigits);
00619   Int_t rawDigitCounter = 0;
00620   RawDigitInfo *tmp_digit = 0;
00621   
00622   while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawDigIter()))) {
00623     Bool_t gotDigit = false;
00624     //get rmmu index
00625     Int_t rmmuInd = linkRmMuRawDig[rawDigitCounter];
00626     MSG("EvtMrg", Msg::kDebug) << "RawDigit: " << rawDigitCounter 
00627                                << " rmmu index: " << rmmuInd << endl;
00628     //now loop over merged digit list
00629     mergedlistIter.Reset();
00630     while(CandDigitHandle *mergeDig = dynamic_cast<CandDigitHandle*>(mergedlistIter())){
00631       //look for equivalent channels:
00632       if(tmp_digit->rcid == mergeDig->GetChannelId()) {
00633         MSG("EvtMrg", Msg::kVerbose) << "Got channel match for electron rawdigit " 
00634                                      << rawDigitCounter << endl;
00635         //get raw digit and check that TDC is the same
00636         const RawDigit* raw_digit = input_datablock->At(mergeDig->GetRawDigitIndex());
00637         //also check against digit info since electron-only hits have no associated RawDigit
00638         Double_t rdTime = calibrator.GetTimeFromTDC(int(tmp_digit->tdc), tmp_digit->rcid);
00639 
00640         Double_t rdAdc = tmp_digit->adc - 50.;   // 50 is a qie offset
00641         if(det == Detector::kFar)  rdAdc = tmp_digit->adc;
00642 
00643         if( ( raw_digit!=NULL && tmp_digit->tdc==raw_digit->GetTDC() ) || 
00644             ( TMath::Abs(rdTime - mergeDig->GetTime())<1e-9 && 
00645               TMath::Abs(rdAdc  - mergeDig->GetCharge())<1e-6 ) ) {
00646           //found match; 
00647           MSG("EvtMrg", Msg::kDebug) << "Got TDC match for electron rawdigit " 
00648                                      << rawDigitCounter << endl;
00649           MSG("EvtMrg", Msg::kDebug) << "Time: " << rdTime*1e6 << " " 
00650                                      << mergeDig->GetTime()*1e6
00651                                      << " ADC: " << rdAdc << " " << mergeDig->GetCharge()
00652                                      << endl;
00653           rmmuItr.Reset();
00654           Int_t rmmuIndex = 0;
00655           while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
00656             if(rmmuIndex==rmmuInd) {
00657               MSG("EvtMrg", Msg::kDebug) << "Got RmMu match for electron rawdigit " 
00658                                          << rawDigitCounter << endl;
00659               //got rmmu object, now look for match in rmmu digit list
00660               TIter rmmuDigIter(rmmu->GetDaughterIterator());
00661               while(CandDigitHandle *rmmuDig = dynamic_cast<CandDigitHandle*>(rmmuDigIter())){
00662                 Int_t rfk = rmmu->ReasonForKeeping(rmmuDig);
00663                 if(rmmuDig->IsEquivalent(mergeDig)) {
00664                   MSG("EvtMrg", Msg::kDebug) << "Got RmMu Digit match for electron rawdigit " 
00665                                              << rawDigitCounter << endl;
00666                   //this digit is already in list, set IsMCElec flag
00667                   rfk |= RmMuMask::kRMMU_ISMCELEC_MASK;
00668                   rmmu->SetReasonForKeeping(rmmuDig,rfk);
00669                   gotDigit = true;
00670                   break;
00671                 }
00672               }
00673               if(!gotDigit) {
00674                 //then add the daughter link and set the flag
00675                 MSG("EvtMrg", Msg::kDebug) << "Adding electron digit for electron rawdigit " 
00676                                            << rawDigitCounter << endl;
00677                 Int_t rfk = 0;
00678                 rfk |= RmMuMask::kRMMU_ISMCELEC_MASK;
00679                 rmmu->AddDaughterLink(*mergeDig);
00680                 rmmu->SetReasonForKeeping(mergeDig,rfk);
00681               }
00682             }
00683             if(gotDigit) break;
00684             rmmuIndex++;
00685           } // end of loop over rmmulist
00686         } // found a raw digit match
00687       } // found a channel ID match
00688       if(gotDigit) break;
00689     } // end of loop over new merged digitlist
00690     rawDigitCounter++;
00691   } // end of loop over MC electron raw digits
00692 
00693   fRawDigits->Clear();
00694   delete fRawDigits;
00695   
00696   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00697 }


Member Data Documentation

Definition at line 42 of file MergeEvent.h.

Referenced by EndJob(), and Reco().

TFile* MergeEvent::fInputFile [private]

Definition at line 40 of file MergeEvent.h.

Referenced by EndJob(), and InfileSetUp().

std::string MergeEvent::fInputFileName [private]

Definition at line 39 of file MergeEvent.h.

Referenced by Config(), EndJob(), and InfileSetUp().

Float_t MergeEvent::fInputMRM[6] [private]

Definition at line 45 of file MergeEvent.h.

Referenced by InfileSetUp(), MergeEvent(), and Reco().

Float_t MergeEvent::fInputP4[4] [private]

Definition at line 44 of file MergeEvent.h.

Referenced by InfileSetUp(), MergeEvent(), and Reco().

TTree* MergeEvent::fInputTree [private]

Definition at line 41 of file MergeEvent.h.

Referenced by InfileSetUp(), and Reco().

TClonesArray* MergeEvent::fInRawDigits [private]

Definition at line 43 of file MergeEvent.h.

Referenced by InfileSetUp(), Reco(), and ~MergeEvent().

std::string MergeEvent::fMergeAlg [private]

Definition at line 57 of file MergeEvent.h.

Referenced by Config(), and Reco().

std::string MergeEvent::fMergeConfig [private]

Definition at line 58 of file MergeEvent.h.

Referenced by Config(), and Reco().

std::string MergeEvent::fMergeListOut [private]

Definition at line 59 of file MergeEvent.h.

Referenced by Config(), and Reco().

Definition at line 51 of file MergeEvent.h.

Referenced by Ana(), and EndJob().

TFile* MergeEvent::fOutputFile [private]

Definition at line 49 of file MergeEvent.h.

Referenced by EndJob(), and OutfileSetUp().

std::string MergeEvent::fOutputFileName [private]

Definition at line 48 of file MergeEvent.h.

Referenced by Config(), EndJob(), and OutfileSetUp().

Float_t MergeEvent::fOutputMRM[6] [private]

Definition at line 54 of file MergeEvent.h.

Referenced by Ana(), MergeEvent(), and OutfileSetUp().

Float_t MergeEvent::fOutputP4[4] [private]

Definition at line 53 of file MergeEvent.h.

Referenced by Ana(), MergeEvent(), and OutfileSetUp().

TTree* MergeEvent::fOutputTree [private]

Definition at line 50 of file MergeEvent.h.

Referenced by Ana(), EndJob(), and OutfileSetUp().

TClonesArray* MergeEvent::fOutRawDigits [private]

Definition at line 52 of file MergeEvent.h.

Referenced by Ana(), OutfileSetUp(), and ~MergeEvent().

Definition at line 60 of file MergeEvent.h.

Referenced by Config(), and Reco().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1