loon: DmxDeMuxModule.cxx Source File

00001 
00002 //$Id: DmxDeMuxModule.cxx,v 1.61 2005/03/26 21:41:49 gmieg Exp $
00003 //
00004 //DmxDeMuxModule.cxx
00005 //
00006 //Module for demultiplexing the far detector
00007 //
00008 //Author:  B. Rebel 10/2000
00010 
00011 #include "DeMux/DmxDeMuxModule.h"
00012 #include "DeMux/DmxChiSqrStat.h"
00013 #include "DeMux/DmxUtilities.h"
00014 #include "JobControl/JobCModuleRegistry.h"
00015 #include "JobControl/JobCommand.h"
00016 #include "MessageService/MsgService.h"
00017 #include "MinosObjectMap/MomNavigator.h"
00018 #include "Algorithm/AlgConfig.h"
00019 #include "CandDigit/CandDeMuxDigitListHandle.h"
00020 #include "CandDigit/CandDigitHandle.h"
00021 #include "CandDigit/CandDeMuxDigitHandle.h"
00022 #include "Algorithm/AlgFactory.h"
00023 #include "Algorithm/AlgHandle.h"
00024 #include "CandData/CandRecord.h"
00025 #include "CandData/CandHeader.h"
00026 #include "Candidate/CandContext.h"
00027 #include "Navigation/NavKey.h"
00028 #include "Navigation/NavSet.h"
00029 #include "Conventions/PlaneView.h"
00030 #include "CandTrackSR/CandTrackSR.h"
00031 #include "Navigation/XxxItr.h"
00032 #include "RecoBase/CandTrackHandle.h"
00033 #include "RecoBase/CandTrackListHandle.h"
00034 #include "DeMux/DmxShowerPlane.h"
00035 #include "DeMux/DmxMuonPlane.h"
00036 #include "DeMux/DmxPlane.h"
00037 #include "DeMux/DmxHypothesis.h"
00038 #include "TH1.h"
00039 #include "TH2.h"
00040 #include "TClonesArray.h"
00041 #include "RawData/RawRecord.h"
00042 #include "RawData/RawDaqSnarlHeader.h"
00043 #include "Validity/VldContext.h"
00044 #include "UgliGeometry/UgliGeomHandle.h"
00045 #include "DataUtil/Truthifier.h"
00046 
00047 #include "TObjArray.h"
00048 
00049 #include <cassert>
00050 
00051 //define the NavKey to sort the planes
00052 static NavKey KeyOnPlaneNum( const TObject *tobj){
00053   const DmxPlane *currentPlane = dynamic_cast<const DmxPlane *>(tobj);
00054   return currentPlane->GetPlaneNumber();
00055 }
00056 
00057 //define the NavKey to sort the digits by plane
00058 static NavKey KeyFrmDigitPlane( const CandDeMuxDigitHandle *cdh){
00059   return cdh->GetPlexSEIdAltL().GetPlane();}
00060 
00061 ClassImp(DmxDeMuxModule)
00062 
00063 CVSID("$Id: DmxDeMuxModule.cxx,v 1.61 2005/03/26 21:41:49 gmieg Exp $");
00064 
00065 // Declare this module to JobControll. Arguments are:
00066 //  (1) The class name 
00067 //  (2) The human-readable name 
00068 //  (3) A short, human-readable description of what the module does
00069 JOBMODULE(DmxDeMuxModule, 
00070           "DeMuxModule",
00071           "A module used for demultiplexing the far detector");
00072 
00073 //......................................................................
00074 DmxDeMuxModule::DmxDeMuxModule() :
00075   fAverageTimingOffset(-5.),
00076   fCrate(-1),
00077   fDataType(1),
00078   fDeMuxedStrip(-1),
00079   fDeMuxedStripW(-1),
00080   fDigitCharge(0.),
00081   fDigitDeMuxedCorrectly(0),
00082   fDigitsInEvent(0),
00083   fDigitSide(0),
00084   fDirection(0),
00085   fDontUseCandDigitMasks(false),
00086   fDoubleEnded(0),
00087   fEndPlane(-1),
00088   fEventCharge(0.),
00089   fEventDeMuxed(0),
00090   fEventLength(false),
00091   fEventNumber(0),
00092   fFilter(false),
00093   fHistoFileName("demux_beam.root"),
00094   fHoughSlope(5.),
00095   fHoughIntercept(10.),
00096   fHoughPeak(0.66),
00097   fHypothesisCoG(-1.),
00098   fHypothesisMatedSignalRatio(-1.),
00099   fHypothesisNumber(-1),
00100   fHypothesisRank(-1),
00101   fHypothesisStat(-1.),
00102   fHypothesisTieBreaker(-1.),
00103   fHypothesisValid(-1),
00104   fMakeTrees(true),
00105   fMated(0),
00106   fMatedSignalForValid(0.5),
00107   fMCDigitCharge(0.),
00108   fMCDigitSide(0),
00109   fMCCoG(-1.),
00110   fMCCoGTPos(-1.),
00111   fMCInitialSlopeU(0.),
00112   fMCInitialSlopeV(0.),
00113   fMCVertexPlane(0),
00114   fMCEndPlane(0),
00115   fMultiple(0),
00116   fMuonStartPlane(-1),
00117   fNonPhysicalFailure(0),
00118   fNonPhysicalPlanes(0),
00119   fNoVertexFailure(0),
00120   fNumberOfPlanesInEvent(0),
00121   fOneSidedUPlanes(0),
00122   fOneSidedVPlanes(0),
00123   fPlane(-1),
00124   fPlaneDeMuxedCoG(-1.),
00125   fPlaneStray(0),
00126   fPlanesInWindow(5),
00127   fPlaneType(-1),
00128   fSignalFractionLimit(0.1),
00129   fSingleModuleDeMux(false),
00130   fSingleModuleCharge(-1.),
00131   fSingleModuleCorrectCharge(-1.),
00132   fSingleModuleDigits(-1.),
00133   fSingleModuleCorrectDigits(-1.),
00134   fStrayDeltaStripLimit(6),
00135   fStrip0(0),
00136   fStrip1(0),
00137   fStrip2(0),
00138   fStrip3(0),
00139   fStrip4(0),
00140   fStrip5(0),
00141   fStrip6(0),
00142   fStrip7(0),
00143   fTime(-1.),
00144   fTruthStripE(-1),
00145   fTruthStripW(-1),
00146   fUStrayPlanes(0),
00147   fUValidPlanes(0),
00148   fVStrayPlanes(0),
00149   fVValidPlanes(0),
00150   fXTalkDigit(1),
00151   fUnsetDigitCheck(false),
00152   fValidPlanesFailure(0),
00153   fVertexPlane(-1),
00154   fWriteHistos(false),
00155   fXTalkSignalRatio(-1.),
00156   fUseStrayPlaneCheck(false),
00157   fStrayPlanesLimit(4),
00158   fDeMuxFile(0),
00159   fDeMuxTree(0),
00160   fDeMuxedDigitTree(0),
00161   fDmxPlaneTree(0),
00162   fTrackDigitsTree(0),
00163   fUnusedDigitsTree(0),
00164   fMCDeMuxTree(0),
00165   fMCDeMuxDigitTree(0),
00166   fTFolder(0)
00167 {
00168 
00169   // Put DmxStatus object in /Loon/DeMux folder
00170   if (fTFolder==0) {     
00171     TFolder *lf = dynamic_cast<TFolder *>
00172                          (gROOT->GetRootFolder()->FindObject("Loon"));
00173     if (lf==0) {
00174       lf = gROOT->GetRootFolder()->AddFolder("Loon", "Loon analysis");
00175       gROOT->GetListOfBrowsables()->Add(lf, "Loon");
00176     }
00177     fTFolder = lf->AddFolder("DeMux", "DeMux analysis");
00178     fTFolder->Add(&fStatus);
00179   }
00180 
00181   AlgFactory &af = AlgFactory::GetInstance();
00182   
00183   //register the defaults for AlgDeMuxBeam
00184   af.Register("AlgDeMuxBeam", "default", "libDeMux.so", "AlgConfig");
00185   AlgHandle ah = af.GetAlgHandle("AlgDeMuxBeam", "default");
00186   
00187   //get the AlgConfigDeMux class and set some variables
00188   AlgConfig &acd = ah.GetAlgConfig();
00189   acd.UnLockValues();
00190   acd.Set("HoughInterceptRMSLimit",10.);
00191   acd.Set("HoughPeakLimit",0.60);
00192   acd.Set("HoughSlopeRMSLimit",10.);
00193   acd.Set("HypothesisSize",24);
00194   acd.Set("MuonTrackSlopeLimit",2.0);
00195   acd.Set("NumberOfHypotheses",169);
00196   acd.Set("NumberOfStrips", 192);
00197   acd.Set("PlanesInSet",fPlanesInWindow);
00198   acd.Set("RatioMatedSignalForValid",fMatedSignalForValid);
00199   acd.Set("UseCandDigitMask", 1);
00200   acd.Set("XTalkFractionLimit", 0.1);
00201   acd.Set("IgnoreSignalLimit", 1.5);
00202   acd.Set("MuonPlanesInSet", fPlanesInWindow);
00203   acd.Set("MaxResidual", 0.1);
00204   acd.Set("MaxBadResidual", 1);
00205   acd.Set("MaxMuonRemovalFraction", 0.4);
00206   acd.Set("MaxStripAdjustment", 3);
00207   acd.Set("MinMuonPlanesRequired", 3);
00208   acd.Set("MinMatedChargeFraction", 0.85);
00209   acd.Set("ShowerMuonTransverseDifference", 0.5); //in meters
00210   if(fDataType == 1) acd.Set("AveragePEGainConversion", 60.);
00211   else acd.Set("AveragePEGainConversion", 60.);
00212   acd.Set("StrayDeltaStripLimit",fStrayDeltaStripLimit);
00213   acd.Set("StrayPlanesLimit",10);
00214   acd.Set("UseStrayPlaneCheck", 0);
00215   acd.Set("ResetDigitCut", 0.125); //meters
00216  
00217   if(fUseStrayPlaneCheck){
00218     acd.Set("UseStrayPlaneCheck",(Int_t)fUseStrayPlaneCheck);
00219     //acd.SetStrayPlanesLimit(fStrayPlanesLimit);
00220   }
00221   
00222   if( fDontUseCandDigitMasks ){ acd.Set("UseCandDigitMask",0); }
00223   //MSG("DMX", Msg::kDebug) << fSignalFractionLimit << endl;
00224   if( fSignalFractionLimit > 0.){ acd.Set("XTalkFractionLimit",fSignalFractionLimit);}
00225   
00226   acd.LockValues();
00227   acd.LockKeys();
00228 
00229   //set the demuxing algorithm to use
00230   AlgHandle ahdl = af.GetAlgHandle("AlgDeMuxDigitList", "default");
00231   AlgConfig &acdl = ahdl.GetAlgConfig();
00232   
00233   acdl.UnLockValues();
00234   acdl.Set("DeMuxAlgorithm", "AlgDeMuxBeam");
00235   acdl.Set("DeMuxAlgConfig", "default");
00236   acdl.Set("NormalizeWeights", 0); // Normalize weights to 1 if non-zero
00237   acdl.LockValues();
00238   acdl.LockKeys();
00239 
00240 
00241   MSG("JobC", Msg::kDebug) << "DmxDeMuxModule::Constructor" << endl;
00242 
00243 }
00244 
00245 //----------------------------------------------------------------------
00246 DmxDeMuxModule::~DmxDeMuxModule()
00247 {
00248 
00249   MSG("JobC", Msg::kDebug) << "DmxDeMuxModule::Destructor" << endl;
00250 
00251 }
00252 
00253 //......................................................................
00254 
00255 void DmxDeMuxModule::BeginJob()
00256 {
00257   if(fMakeTrees && !fDeMuxFile){
00258     // save the current working directory
00259     TDirectory* savedir  = gDirectory;  
00260 
00261     //create the new TFile for demuxing
00262     // this changes the value of gDirectory
00263     fDeMuxFile = new TFile(fHistoFileName,"RECREATE");  
00264     
00265     // create TTree's, these will attach themselves to the current 
00266     //working directory (demux_cosmics.root)
00267    
00268     fDeMuxTree = new TTree("fDeMuxTree", "DeMux Info Tree");
00269     fDeMuxedDigitTree = new TTree("fDeMuxedDigitTree", "DeMuxed Digit Info");
00270     fDmxPlaneTree = new TTree("fDmxPlaneTree", "DeMux Plane Info");
00271     
00272     fDeMuxTree->Branch("Event", &fEventNumber, "EventNumber/I"); 
00273     fDeMuxTree->Branch("DeMuxed", &fEventDeMuxed, "EventDeMuxed/I");
00274     fDeMuxTree->Branch("Vertex", &fVertexPlane, "VertexPlane/I");
00275     fDeMuxTree->Branch("VertexZ", &fVertexPlaneZ, "VertexPlaneZ/F");
00276     fDeMuxTree->Branch("Direction", &fDirection, "Direction/I");
00277     fDeMuxTree->Branch("End", &fEndPlane, "EndPlane/I");
00278     fDeMuxTree->Branch("NumberOfPlanesHit", &fNumberOfPlanesInEvent, "NumberOfPlanesInEvent/I");
00279     fDeMuxTree->Branch("EventCharge", &fEventCharge, "EventCharge/F");
00280     fDeMuxTree->Branch("DigitsInEvent", &fDigitsInEvent, "DigitsInEvent/I");
00281     fDeMuxTree->Branch("Multiple", &fMultiple, "Multiple/I");
00282     fDeMuxTree->Branch("MuonStartPlane", &fMuonStartPlane, "MuonStartPlane/I");
00283     fDeMuxTree->Branch("NonPhysicalPlanes", &fNonPhysicalPlanes, "NonPhysicalPlanes/I");
00284     fDeMuxTree->Branch("ValidPlanesFailure", &fValidPlanesFailure, "ValidPlanesFailure/I");
00285     fDeMuxTree->Branch("NoVertexFailure", &fNoVertexFailure, "NoVertexFailure/I");
00286     fDeMuxTree->Branch("NonPhysicalFailure", &fNonPhysicalFailure, "NonPhysicalFailure/I");
00287     fDeMuxTree->Branch("OneSidedUPlanes", &fOneSidedUPlanes, "OneSidedUPlanes/I");
00288     fDeMuxTree->Branch("OneSidedVPlanes", &fOneSidedVPlanes, "OneSidedVPlanes/I");
00289     fDeMuxTree->Branch("SingleModuleCharge", &fSingleModuleCharge, "SingleModuleCharge/F");
00290     fDeMuxTree->Branch("SingleModuleCorrectCharge", &fSingleModuleCorrectCharge, "SingleModuleCorrectCharge/F");
00291     fDeMuxTree->Branch("SingleModuleDigits", &fSingleModuleDigits, "SingleModuleDigits/F");
00292     fDeMuxTree->Branch("SingleModuleCorrectDigits", &fSingleModuleCorrectDigits, "SingleModuleCorrectDigits/F");
00293     fDeMuxTree->Branch("UStrayPlanes", &fUStrayPlanes, "UStrayPlanes/I");
00294     fDeMuxTree->Branch("VStrayPlanes", &fVStrayPlanes, "VStrayPlanes/I");
00295     fDeMuxTree->Branch("UValidPlanes", &fUValidPlanes, "UValidPlanes/I");
00296     fDeMuxTree->Branch("VValidPlanes", &fVValidPlanes, "VValidPlanes/I");
00297     fDeMuxTree->Branch("UOverlappingMultiple", &fUOverlap, "UOverlap/I");
00298     fDeMuxTree->Branch("VOverlappingMultiple", &fVOverlap, "VOverlap/I");
00299     fDeMuxTree->Branch("AvTimingOffset", &fAverageTimingOffset, "AverageTimingOffset/F");
00300     fDeMuxTree->Branch("UMuonMatedFraction", &fUMuonMatedFraction, "UMuonMatedFraction/F");
00301     fDeMuxTree->Branch("VMuonMatedFraction", &fVMuonMatedFraction, "VMuonMatedFraction/F");
00302     fDeMuxTree->Branch("UShowerMatedFraction", &fUShowerMatedFraction, "UShowerMatedFraction/F");
00303     fDeMuxTree->Branch("VShowerMatedFraction", &fVShowerMatedFraction, "VShowerMatedFraction/F");
00304     //fDeMuxTree->SetAutoSave(10000);
00305     
00306     fDeMuxedDigitTree->Branch("Event", &fEventNumber, "EventNumber/I");
00307     fDeMuxedDigitTree->Branch("Plane", &fPlane, "Plane/I");
00308     fDeMuxedDigitTree->Branch("Crate", &fCrate, "Crate/I");
00309     fDeMuxedDigitTree->Branch("ADC", &fDigitCharge, "DigitCharge/F");
00310     fDeMuxedDigitTree->Branch("Time", &fTime, "Time/D");
00311     fDeMuxedDigitTree->Branch("Strip0", &fStrip0, "Strip0/I");
00312     fDeMuxedDigitTree->Branch("Strip1", &fStrip1, "Strip1/I");
00313     fDeMuxedDigitTree->Branch("Strip2", &fStrip2, "Strip2/I");
00314     fDeMuxedDigitTree->Branch("Strip3", &fStrip3, "Strip3/I");
00315     fDeMuxedDigitTree->Branch("Strip4", &fStrip4, "Strip4/I");
00316     fDeMuxedDigitTree->Branch("Strip5", &fStrip5, "Strip5/I");
00317     fDeMuxedDigitTree->Branch("Strip6", &fStrip6, "Strip6/I");
00318     fDeMuxedDigitTree->Branch("Strip7", &fStrip7, "Strip7/I");
00319     fDeMuxedDigitTree->Branch("DeMuxedStrip", &fDeMuxedStrip, "DeMuxedStrip/I");
00320     fDeMuxedDigitTree->Branch("XTalkDigit", &fXTalkDigit, "XTalkDigit/I");
00321     //fDeMuxedDigitTree->SetAutoSave(100000);
00322     
00323     fDmxPlaneTree->Branch("Event", &fEventNumber, "EventNumber/I");
00324     fDmxPlaneTree->Branch("Plane", &fPlane, "Plane/I");
00325     fDmxPlaneTree->Branch("PlaneZ", &fPlaneZ, "PlaneZ/F");
00326     fDmxPlaneTree->Branch("PlaneType", &fPlaneType, "PlaneType/I");
00327     fDmxPlaneTree->Branch("PlaneCoG", &fPlaneDeMuxedCoG, "PlaneDeMuxedCoG/F");
00328     fDmxPlaneTree->Branch("Stray", &fPlaneStray, "PlaneStray/I");
00329     fDmxPlaneTree->Branch("HypothesisCoG", &fHypothesisCoG, "HypothesisCoG/F");
00330     fDmxPlaneTree->Branch("HypothesisCoGTPos", &fHypothesisCoGTPos, "HypothesisCoGTPos/F");
00331     fDmxPlaneTree->Branch("Hypothesis", &fHypothesisNumber, "HypothesisNumber/I");
00332     fDmxPlaneTree->Branch("MatedSignal", &fHypothesisMatedSignalRatio, "HypothesisMatedSignalRatio/F");
00333     fDmxPlaneTree->Branch("Statistic", &fHypothesisStat, "HypothesisStat/F");
00334     fDmxPlaneTree->Branch("TieBreaker", &fHypothesisTieBreaker, "HypothesisTieBreaker/F");
00335     fDmxPlaneTree->Branch("Rank", &fHypothesisRank, "HypothesisRank/I");
00336     fDmxPlaneTree->Branch("Valid", &fHypothesisValid, "HypothesisValid/I");
00337 
00338     fMCDeMuxDigitTree = new TTree("fMCDeMuxDigitTree", "MC DeMux Digit Tree");
00339     fMCDeMuxDigitTree->Branch("Event", &fEventNumber, "EventNumber/I");
00340     fMCDeMuxDigitTree->Branch("TruthDigitCharge", &fMCDigitCharge, "MCDigitChargeE/F");
00341     fMCDeMuxDigitTree->Branch("TruthDigitSide", &fMCDigitSide, "MCDigitSide/I");
00342     fMCDeMuxDigitTree->Branch("Plane", &fPlane, "Plane/I");
00343     fMCDeMuxDigitTree->Branch("Strip", &fDeMuxedStrip, "DeMuxedStrip/I");
00344         fMCDeMuxDigitTree->Branch("XTalkDigit", &fXTalkDigit, "XTalkDigit/I");
00345         fMCDeMuxDigitTree->Branch("DoubleEnded", &fDoubleEnded, "DoubleEnded/I");
00346     fMCDeMuxDigitTree->Branch("MCCoG", &fMCCoG, "fMCCoG/F");
00347     fMCDeMuxDigitTree->Branch("MCCoGTPos", &fMCCoGTPos, "fMCCoGTPos/F");
00348     
00349     fMCDeMuxTree = new TTree("fMCDeMuxTree", "MC DeMux Tree");
00350     fMCDeMuxTree->Branch("Event", &fEventNumber, "EventNumber/I"); 
00351     fMCDeMuxTree->Branch("Vertex", &fMCVertexPlane, "VertexPlane/I");
00352     fMCDeMuxTree->Branch("End", &fMCEndPlane, "EndPlane/I");
00353     fMCDeMuxTree->Branch("InitialSlopeU", &fMCInitialSlopeU, "InitialSlopeU/D");
00354     fMCDeMuxTree->Branch("InitialSlopeV", &fMCInitialSlopeV, "InitialSlopeV/D");
00355     fMCDeMuxTree->Branch("ChargeCorrect", &fMCChargeCorrect, "MCChargeCorrect/F");
00356     fMCDeMuxTree->Branch("ChargeCorrect1Hyp", &fMCChargeCorrect1Hyp, "MCChargeCorrect1Hyp/F");
00357     fMCDeMuxTree->Branch("TotalCharge", &fMCTotalCharge, "MCTotalCorrect/F");
00358     fMCDeMuxTree->Branch("TotalCharge1Hyp", &fMCTotalCharge1Hyp, "MCTotalCharge1Hyp/F");
00359     fMCDeMuxTree->Branch("StripsCorrect", &fMCStripsCorrect, "MCStripsCorrect/F");
00360        
00361     // change back to current working directory before leaving constructor
00362     savedir->cd();   
00363   }
00364 
00365   if( fUnsetDigitCheck && fMakeTrees){
00366 
00367     // save the current working directory
00368     TDirectory* savedir = gDirectory; 
00369     fDeMuxFile->cd();
00370 
00371     fUnusedDigitsTree = new TTree("fUnusedDigitsTree", "UnsetDigitCheck Test Tree");
00372     fUnusedDigitsTree->Branch("Event", &fEventNumber, "EventNumber/I"); 
00373     fUnusedDigitsTree->Branch("DigitCharge", &fDigitCharge, "DigitCharge/F");
00374     fUnusedDigitsTree->Branch("XTalkSignalRatio", &fXTalkSignalRatio, "XTalkSignalRatio/D");
00375     fUnusedDigitsTree->Branch("DigitSide", &fDigitSide, "DigitSide/F");
00376     
00377     // change back to current working directory
00378     savedir->cd(); 
00379   }
00380   
00381   return;
00382 }
00383 
00384 //......................................................................
00385 
00386 JobCResult DmxDeMuxModule::Reco(MomNavigator *mom)
00387 {
00388   //this method is depricated so have it always return a kFailed to
00389   //alert anyone using it that they shouldn't
00390 
00391    JobCResult result(JobCResult::kFailed);
00392 
00393 //   MSG("JobC", Msg::kDebug) << "DeMuxModule::Reco\n";
00394 
00395    // Check that mom exists.
00396    assert(mom);
00397 
00398 //   fStatus.ClearPlaneArray();
00399 
00400 //   // Find PrimaryCandidateRecord fragment in MOM.
00401 //   CandRecord *candrec = dynamic_cast<CandRecord *>
00402 //     (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00403 //   if (candrec == 0) {
00404 //     result.SetError().SetFailed();
00405 //     return result;
00406 //   }
00407   
00408 //   //*-* Find Raw Data CandDigitList fragment in PrimaryCandidateRecord.
00409 //   CandDigitListHandle *crdlh = dynamic_cast<CandDigitListHandle *>
00410 //     (candrec->FindCandHandle("CandDigitListHandle", "canddigitlist"));
00411   
00412 //   // Check that crdlh exists.
00413 //   assert(crdlh);
00414   
00415 //   Int_t event = fStatus.GetEventNumber();
00416 //   const CandHeader* chead = candrec->GetCandHeader();
00417 //   if (!chead) {
00418 //     MSG("JobC",Msg::kWarning) << "No CandHeader, but continuing" << endl;
00419 //   }
00420 //   event++;// chead->GetEvent();
00421   
00422 //   AlgFactory &af = AlgFactory::GetInstance();
00423 //   AlgHandle ah = af.GetAlgHandle("AlgDeMuxBeam", "default");
00424 //   //AlgConfigDeMux &acd = dynamic_cast<AlgConfigDeMux &>(ah.GetAlgConfig());
00425   
00426 //   CandContext cx(this, mom);
00427 //   cx.SetDataIn(&fStatus);
00428   
00429 //   ah.RunAlg(*crdlh, cx);
00430   
00431 //   //if the event fails demuxing for any reason, dont pass it on
00432 //   if( !fStatus.GetEventDeMuxed() && fFilter ){
00433 //     result.SetError().SetFailed();
00434 //     return result;
00435 //   }
00436 
00437 //   MSG("JobC", Msg::kDebug) << "Event DeMuxed" << endl;
00438   
00439    return result;
00440 }
00441 
00442 //......................................................................
00443 //this method can be used to find a figure of merit for real data or
00444 //use DmxValidate with monte carlo data
00445 
00446 JobCResult DmxDeMuxModule::Ana(const MomNavigator *mom)
00447 {
00448         JobCResult result(JobCResult::kPassed);
00449   
00450         //MSG("JobCDmx", Msg::kDebug) << "DeMuxCosmicsModule::Ana" << endl;
00451 
00452         DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
00453         if (status == 0) {
00454                 MSG("DMX", Msg::kWarning) << "//root/Loon/DeMux/DmxStatus not found." << endl;
00455     
00456                 result.SetError().SetFailed();
00457                 return result;
00458         }
00459   
00460         //if the event fails demuxing for any reason and filtering is turned on, 
00461         //dont pass the event on
00462   
00463         if(fFilter && fUseStrayPlaneCheck && 
00464            (status->GetVOverlappingMultiple() || status->GetUOverlappingMultiple())){
00465                 result.SetError().SetFailed();
00466                 return result;
00467         }
00468         else if( !status->GetEventDeMuxed() && fFilter){
00469                 result.SetError().SetFailed();
00470                 return result;
00471         }
00472   
00473         //only do stuff in this method if you are making trees - otherwise there 
00474         //is no point to running the method.
00475         if( fMakeTrees ){
00476     
00477                 DmxUtilities util;
00478 
00479                 // Check that mom exists.
00480                 assert(mom);
00481     
00482                 // Find PrimaryCandidateRecord fragment in MOM.
00483                 CandRecord *candrec = dynamic_cast<CandRecord *>
00484                         (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00485                 if (candrec == 0) {
00486                         result.SetError().SetFailed();
00487                         return result;
00488                 }
00489     
00490                 //Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
00491                 CandDeMuxDigitListHandle *crdlh = dynamic_cast<CandDeMuxDigitListHandle *>
00492                         (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
00493     
00494                 // Check that crdlh exists.
00495                 assert(crdlh);
00496     
00497                 //get an iterator over the digits and one extra for the single module test
00498                 CandDeMuxDigitHandleItr crdhi(crdlh->GetDaughterIterator());
00499                 CandDeMuxDigitHandleItr smcrdhi(crdlh->GetDaughterIterator());
00500 
00501                 //create the KeyFunc
00502                 CandDeMuxDigitHandleKeyFunc *planeKF = crdhi.CreateKeyFunc();
00503                 CandDeMuxDigitHandleKeyFunc *smPlaneKF = smcrdhi.CreateKeyFunc();
00504     
00505                 //program the KeyFunc with the sort function
00506                 planeKF->SetFun(KeyFrmDigitPlane);
00507                 smPlaneKF->SetFun(KeyFrmDigitPlane);
00508     
00509                 //get the NavSet from the iterator and pass the KeyFunc to it
00510                 crdhi.GetSet()->AdoptSortKeyFunc(planeKF);
00511                 smcrdhi.GetSet()->AdoptSortKeyFunc(smPlaneKF);
00512     
00513                 //clear the KF pointer because we no longer own the KeyFunc
00514                 planeKF = 0;
00515                 smPlaneKF = 0;
00516     
00517                 crdhi.ResetFirst();
00518                 smcrdhi.ResetFirst();
00519 
00520                 //get the TObjArray of DmxPlane objects
00521                 const TObjArray *planeArray = status->GetPlaneArray();
00522     
00523                 //create a TObjectItr
00524                 TObjectItr planeItr(planeArray);
00525     
00526                 //create a KeyFunc to sort planes by number
00527                 TObjectKeyFunc *pnKF = planeItr.CreateKeyFunc();
00528     
00529                 //program the KeyFunc with the sort function
00530                 pnKF->SetFun(KeyOnPlaneNum);
00531     
00532                 //get the NavSet from the iterator and pass the KeyFunc to it
00533                 planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
00534     
00535                 //clear the KF pointer because we no longer own the KeyFunc
00536                 pnKF = 0;
00537     
00538                 fNonPhysicalPlanes = 0;
00539                 fHypothesisCoG = -1.;
00540                 fPlaneStray = 0;
00541                 fPlaneDeMuxedCoG = -1.;
00542                 fHypothesisNumber = -1;
00543                 fHypothesisStat = -1.;
00544                 fHypothesisTieBreaker = -1.;
00545                 fHypothesisMatedSignalRatio = -1.;
00546                 fHypothesisValid = -1;
00547                 fPlaneType = -1;
00548                 fCrate = -1;
00549                 fPlane = -1;
00550                 fDeMuxedStrip = -1;
00551                 fOneSidedUPlanes = 0;
00552                 fOneSidedVPlanes = 0;
00553                 fSingleModuleCharge = -1.;
00554                 fSingleModuleDigits = -1.;
00555                 fSingleModuleCorrectCharge = -1.;
00556                 fSingleModuleCorrectDigits = -1.;
00557                 fEventCharge = 0.;
00558                 fDigitsInEvent = 0;
00559                 fDigitCharge = -1.;
00560                 fXTalkSignalRatio = -1.;
00561                 fDigitDeMuxedCorrectly = 0;
00562                 fDigitSide = 0;
00563                 //get generic info for event
00564                 fEventNumber = status->GetEventNumber();
00565                 fEventDeMuxed = (Int_t)status->GetEventDeMuxed();
00566                 fValidPlanesFailure = (Int_t)status->GetValidPlanesFailure();
00567                 fNoVertexFailure = (Int_t)status->GetVertexPlaneFailure();
00568                 fNonPhysicalFailure = (Int_t)status->GetNonPhysicalFailure();
00569                 fUOverlap = (Int_t)status->GetUOverlappingMultiple();
00570                 fVOverlap = (Int_t)status->GetVOverlappingMultiple();
00571                 fAverageTimingOffset = status->GetAverageTimingOffset();
00572                 fUSlope = status->GetUTrackSlope();
00573                 fUIntercept = status->GetUTrackIntercept();
00574                 fUSlopeRMS = status->GetUSlopeRMS();
00575                 fUStrayPlanes = status->GetUStrayPlanes();
00576                 fUInterceptRMS = -1.;
00577                 fUChiSq = status->GetUTrackChiSq();
00578                 fVSlope = status->GetVTrackSlope();
00579                 fVIntercept = status->GetVTrackIntercept();
00580                 fVSlopeRMS = status->GetVSlopeRMS();
00581                 fVStrayPlanes = status->GetVStrayPlanes();
00582                 fVInterceptRMS = -1.;
00583                 fVChiSq = status->GetVTrackChiSq();
00584                 fUValidPlanes = 0;
00585                 fVValidPlanes = 0;
00586                 fUMuonMatedFraction = status->GetUMuonMatedFraction();
00587                 fVMuonMatedFraction = status->GetVMuonMatedFraction();
00588                 fUShowerMatedFraction = status->GetUShowerMatedFraction();
00589                 fVShowerMatedFraction = status->GetVShowerMatedFraction();
00590                 
00591                 //zero values that need to start over with each event
00592                 fVertexPlane = status->GetVertexPlaneNumber();
00593                 fVertexPlaneZ = status->GetVertexZPosition();
00594                 fDirection = status->GetEventDirection();
00595                 fEndPlane = status->GetEndPlaneNumber();
00596                 fMultiple = (Int_t)status->GetMultipleMuon();
00597                 fMuonStartPlane = status->GetMuonStartPlaneNumber();
00598 
00599                 if(fMultiple == 0){planeItr.GetSet()->Slice(fVertexPlane, fEndPlane);}
00600                 else if(fMultiple == 1){planeItr.GetSet()->Slice(fVertexPlane, 500);}
00601     
00602                 //-----------------------------get plane info--------------------------//
00603     
00604                 //get the total charge and digits in the event between vertex and end planes
00605                 while( planeItr.IsValid() ){
00606                         DmxPlane *currentPlane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
00607                         fEventCharge += currentPlane->GetPlaneCharge();
00608                         fPlane = currentPlane->GetPlaneNumber();
00609                         fPlaneZ = currentPlane->GetZPosition();
00610                         fPlaneDeMuxedCoG = currentPlane->GetStripCoG(); 
00611                         fPlaneStray = (Int_t)currentPlane->IsStray();
00612                         if(currentPlane->GetPlaneView() == PlaneView::kU){
00613                                 if(currentPlane->IsValid()){ ++fUValidPlanes; }
00614                         }
00615                         else if(currentPlane->GetPlaneView() == PlaneView::kV){
00616                                 if(currentPlane->IsValid()){ ++fVValidPlanes; }
00617                         }
00618                         //if(fPlaneFitCoG > 192. || fPlaneFitCoG < 0.){ ++fNonPhysicalPlanes; }
00619                         if(TMath::Abs(fPlaneDeMuxedCoG) > 4.){ ++fNonPhysicalPlanes; }
00620       
00621                         crdhi.GetSet()->Slice(fPlane);
00622                         crdhi.GetSet()->Update();
00623       
00624                         //fill the DmxPlane tree
00625                         if( currentPlane->IsValid() && currentPlane->GetPlaneType() == DmxPlaneTypes::kShower ){
00626                                 fPlaneType = 0;
00627         
00628                                 //loop over the best hypotheses to get their stats
00629                                 if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()){
00630                                         fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetLowerBound();
00631                                         fHypothesisRank = 1;
00632                                         fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetCoG();
00633                                         fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestCoG();
00634                                         fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetStat();
00635                                         fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetTieBreakerStat();
00636                                         fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->IsValid();
00637                                         fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetMatedSignalRatio();
00638 
00639                                         if(status->GetEventDeMuxed())  fDmxPlaneTree->Fill();
00640 
00641                                 }
00642                                 if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis() ){
00643                                         fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetLowerBound();
00644                                         fHypothesisRank = 2;
00645                                         fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetCoG();
00646                                         fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestCoG();
00647                                         fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetStat();
00648                                         fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetTieBreakerStat();
00649                                         fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->IsValid();
00650                                         fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
00651                                         if(status->GetEventDeMuxed())  fDmxPlaneTree->Fill();
00652 
00653                                 }
00654                                 if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis() ){
00655                                         fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetLowerBound();
00656                                         fHypothesisRank = 3;
00657                                         fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetCoG();
00658                                         fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestCoG();
00659                                         fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetStat();
00660                                         fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetTieBreakerStat();
00661                                         fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->IsValid();
00662                                         fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
00663                                         if(status->GetEventDeMuxed())  fDmxPlaneTree->Fill();
00664                                 }
00665         
00666                                 //fDmxPlaneTree->Fill();
00667                                 fHypothesisStat = -1.;
00668                                 fHypothesisMatedSignalRatio = -1.;
00669                                 fHypothesisTieBreaker = -1.;
00670                                 fHypothesisRank = -1;
00671                                 fHypothesisValid = -1;
00672                                 fHypothesisCoG = -1.;
00673                         }
00674                         else if( currentPlane->IsValid() && currentPlane->GetPlaneType() == DmxPlaneTypes::kMuon ){
00675                                 fPlaneType = 1;
00676                                 if((Int_t)currentPlane->GetStripCoG() > 11){
00677                                         fHypothesisNumber = (Int_t)currentPlane->GetStripCoG() - 11;
00678                                 }
00679                                 else fHypothesisNumber = 0; 
00680         
00681                                 fHypothesisTieBreaker = -1;
00682                                 fHypothesisRank = 1;
00683                                 fHypothesisCoGTPos = dynamic_cast<DmxMuonPlane *>(currentPlane)->GetInitialCoG();
00684                                 fHypothesisCoG = dynamic_cast<DmxMuonPlane *>(currentPlane)->GetInitialStripCoG();
00685                                 fHypothesisMatedSignalRatio = 1.;
00686                                 fHypothesisValid = (Int_t)currentPlane->IsValid();
00687                                 if(status->GetEventDeMuxed()){
00688                                         fDmxPlaneTree->Fill();
00689                                         if ( fEventNumber%1000 == 0 ) {
00690                                                 // save current working directory
00691                                                 TDirectory* savedir = gDirectory; 
00692                                                 // move to file directory of tree, oddly SaveSelf 
00693                                                 //doesn't work without it
00694                                                 fDeMuxFile->cd();
00695                                                 // save tree to file  
00696                                                 fDmxPlaneTree->AutoSave();
00697                                                 // save file directory containing this tree
00698                                                 fDeMuxFile->SaveSelf();
00699                                                 // change back to working to directory
00700                                                 savedir->cd();  
00701                                         }
00702                                 }
00703                         }
00704       
00705       //-------------------get digit info----------------------------------//
00706       
00707                         fDigitsInEvent += crdhi.SizeSelect();
00708                         
00709                         //loop over the digits to fill the demuxed digit info
00710                         CandDeMuxDigitHandle *digit = 0;
00711                         while( (digit = crdhi()) ){
00712                                 
00713                                 fStrip0 = 0;
00714                                 fStrip1 = 0;
00715                                 fStrip2 = 0;
00716                                 fStrip3 = 0;
00717                                 fStrip4 = 0;
00718                                 fStrip5 = 0;
00719                                 fStrip6 = 0;
00720                                 fStrip7 = 0;
00721                                 fXTalkDigit = 1;
00722                                 fTime = digit->GetTime();
00723                                 //get the strip, charge, crate for the digit
00724                                 if(digit->GetPlexSEIdAltL().GetBestWeight() == 1.){
00725                                         fDeMuxedStrip = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
00726                                 }
00727                                 else fDeMuxedStrip = -1;
00728         
00729                                 //set the alternatives list to the first one
00730                                 digit->GetPlexSEIdAltL().SetFirst();
00731                                 
00732                                 //loop over the strip end alt list
00733                                 Int_t cnt = 0;
00734                                 while( digit->GetPlexSEIdAltL().IsValid() ){
00735                                         if(cnt == 0){ fStrip0 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00736                                         else if(cnt == 1){ fStrip1 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00737                                         else if(cnt == 2){ fStrip2 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00738                                         else if(cnt == 3){ fStrip3 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00739                                         else if(cnt == 4){ fStrip4 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00740                                         else if(cnt == 5){ fStrip5 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00741                                         else if(cnt == 6){ fStrip6 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00742                                         else if(cnt == 7){ fStrip7 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
00743                                         ++cnt;
00744                                         digit->GetPlexSEIdAltL().Next();
00745                                 }
00746         
00747         
00748                                 fDigitCharge = digit->GetCharge();
00749         
00750                                 //MSG("DMX", Msg::kDebug) << status->GetEventNumber() << "\t" << fPlane
00751                                 //                                              << "\t" << fDeMuxedStrip
00752                                 //                                              << "\t" << fDigitCharge << endl;
00753         
00754                                 if( digit->GetDeMuxDigitFlagWord() == 0) fXTalkDigit = 0;
00755                                 
00756                                 fCrate = digit->GetChannelId().GetCrate();
00757                                 if(status->GetEventDeMuxed()){
00758                                         fDeMuxedDigitTree->Fill();
00759                                         if ( fEventNumber%1000 == 0 ) {
00760                                                 TDirectory* savedir = gDirectory; // save current working directory
00761                                                 fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
00762                                                 fDeMuxedDigitTree->AutoSave();  // save tree to file
00763                                                 fDeMuxFile->SaveSelf();  // save file directory containing this tree
00764                                                 savedir->cd();  // change back to working to directory
00765                                         }
00766                                 }
00767                         }//end loop over digits
00768       
00769                         crdhi.GetSet()->ClearSlice(false);
00770                         crdhi.ResetFirst();
00771       
00772                         planeItr.Next();
00773                 }//end loop over planes
00774     
00775                 //if you arent using the cand digit masks to demux, then set them all to true
00776                 //if(fDontUseCandDigitMasks){ crdhi.GetSet()->GetMasks().SetAllMasks(kTRUE); }
00777     
00778                 planeItr.ResetFirst();
00779                 planeItr.GetSet()->ClearSlice(false);
00780 
00781                 Int_t totalPlanes = planeItr.SizeSelect();
00782                 fNumberOfPlanesInEvent = totalPlanes;
00783     
00784                 //------------------start the efficiency tests----------------------------//
00785     
00786                 if( fDataType == 0){
00787       
00788                         if( fUnsetDigitCheck && status->GetEventDeMuxed() ){
00789         
00790                                 //slice to get planes between the vertex and end
00791                                 planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00792                                 planeItr.GetSet()->Update();
00793                                 planeItr.ResetFirst();
00794         
00795                                 //loop over the planes and print their reconstructions
00796                                 while( planeItr.IsValid() ){
00797                                         DmxPlane *plane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
00798           
00799                                         Int_t planeNumber = plane->GetPlaneNumber();
00800           
00801                                         //slice the digit iterator to the current plane
00802                                         crdhi.GetSet()->Slice(planeNumber);
00803                                         crdhi.GetSet()->Update();
00804 
00805           //array to keep track of which pixels on which planes were hit for
00806           //the set digits
00807           Float_t westSet0[16];
00808           Float_t eastSet0[16];
00809           Float_t westSet1[16];
00810           Float_t eastSet1[16];
00811           Float_t westSet2[16];
00812           Float_t eastSet2[16];
00813           for(Int_t i = 0; i < 16; i++){
00814             westSet0[i] = 0.;
00815             eastSet0[i] = 0.;
00816             westSet1[i] = 0.;
00817             eastSet1[i] = 0.;
00818             westSet2[i] = 0.;
00819             eastSet2[i] = 0.;
00820           }
00821           
00822           util.FillHitPixels(crdhi, 0, eastSet0, westSet0);
00823           crdhi.ResetFirst();
00824           util.FillHitPixels(crdhi, 1, eastSet1, westSet1);
00825           crdhi.ResetFirst();
00826           util.FillHitPixels(crdhi, 2, eastSet2, westSet2);
00827           
00828           crdhi.ResetFirst();
00829           
00830           //loop over the digits and see if they have been set
00831           while( crdhi.IsValid() ){
00832             
00833             //get the digit and it's strip, plane, and side
00834             CandDeMuxDigitHandle *digit = crdhi.Ptr();
00835             fDigitCharge = digit->GetCharge()/55.;
00836             
00837             Int_t chip = digit->GetChannelId().GetVaChip();
00838             Int_t pixel = util.GetDigitPixel(digit->GetChannelId().GetVaChannel());
00839             
00840             //see if the best weight is < 1.0 - if so, the digit had no strip chosen
00841             if(digit->GetPlexSEIdAltL().GetBestWeight() < 1.0){
00842               
00843               Float_t nearestNeighborSignal = -1.;
00844               if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
00845                 if(chip == 0){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet0);}
00846                 if(chip == 1){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet1);}
00847                 if(chip == 2){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet2);}
00848                 fDigitSide = 0;
00849               }
00850               if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
00851                 if(chip == 0){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet0);}
00852                 if(chip == 1){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet1);}
00853                 if(chip == 2){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet2);}
00854                 fDigitSide = 1;
00855               }
00856               
00857               if(nearestNeighborSignal > 0.){fXTalkSignalRatio = fDigitCharge / (nearestNeighborSignal/55.);}
00858               
00859               fUnusedDigitsTree->Fill();
00860               if ( fEventNumber%1000 == 0 ) {
00861                 TDirectory* savedir = gDirectory; // save current working directory
00862                 fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
00863                 fUnusedDigitsTree->AutoSave();  // save tree to file
00864                 fDeMuxFile->SaveSelf();  // save file directory containing this tree
00865                 savedir->cd();  // change back to working to directory
00866               }
00867 
00868               fXTalkSignalRatio = -1.;
00869             }
00870             crdhi.Next();
00871           }
00872           
00873           //clear the digit slice
00874           crdhi.GetSet()->ClearSlice(false);
00875           planeItr.Next();
00876         }
00877         
00878         //reset the digit itr
00879         crdhi.GetSet()->ClearSlice(false);
00880         crdhi.ResetFirst();
00881         
00882         //reset the planeItr
00883         planeItr.GetSet()->ClearSlice(false);
00884         planeItr.ResetFirst();
00885         
00886       } //end UnsetDigitCheck test
00887       
00888       if( fEventLength ){
00889         
00890         if( status->GetValidPlanesFailure() ){ 
00891           //MSG("DmxDigit", Msg::kDebug) << "Event = " << status->GetEventNumber() << endl;
00892           
00893           while( planeItr.IsValid() ){
00894             
00895             crdhi.GetSet()->Slice(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber());
00896             crdhi.GetSet()->Update();
00897 
00898             //declare variables to see if the digits from this plane come from each side
00899             Bool_t westDigits = false;
00900             Bool_t eastDigits = false;
00901             
00902             while( crdhi.IsValid() ){
00903               
00904               if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
00905                 westDigits = true;
00906               }
00907               if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
00908                 eastDigits = true;
00909               }
00910               crdhi.Next();
00911               
00912             } //end loop over digits
00913             
00914             crdhi.GetSet()->ClearSlice(false);
00915             crdhi.ResetFirst();
00916             
00917             
00918             if( dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneView() == PlaneView::kU){
00919               if( !westDigits || !eastDigits){ fOneSidedUPlanes += 1; }
00920             }
00921             if( dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneView() == PlaneView::kV){
00922               if( !westDigits || !eastDigits){ fOneSidedVPlanes += 1; }
00923             }
00924             planeItr.Next();
00925             
00926           } //end loop over planes
00927           
00928           planeItr.ResetFirst();
00929           
00930         }//end valid plane failure
00931         else if( status->GetVertexPlaneFailure() ){ 
00932           if( fNumberOfPlanesInEvent > 3 ){
00933             //MSG("DmxVertex", Msg::kDebug) << "Event = " << status->GetEventNumber() << endl;
00934             
00935             bool westSide = false;
00936             bool eastSide = false;
00937             
00938             while( planeItr.IsValid() ){
00939               
00940               crdhi.GetSet()->Slice(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber());
00941               crdhi.GetSet()->Update();
00942 
00943               while( crdhi.IsValid() ){
00944                 
00945                 if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
00946                   westSide = true;
00947                 }
00948                 if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
00949                                   eastSide = true;
00950                 }
00951                 crdhi.Next();
00952                 
00953               } //end loop over digits
00954               
00955               crdhi.GetSet()->ClearSlice(false);
00956               crdhi.ResetFirst();
00957               
00958               planeItr.Next();
00959               
00960             } //end loop over planes
00961             
00962             planeItr.ResetFirst();
00963           } //if more than 10 planes
00964           
00965         }// end vertex plane failure
00966       }// end if EventLength test
00967       
00968       planeItr.GetSet()->ClearSlice(false);
00969     } //end if FarDetector data
00970     //-------------------------------------------------------------------------------------------------------
00971     else if( fDataType == 1 ){
00972 
00973                 //get the Truthifier so that you know what strip a digit came from
00974                 const Truthifier &truth = Truthifier::Instance(mom);
00975                 
00976                 CandDigitHandle *digit = 0;
00977                         
00978                 //MSG("Dmx", Msg::kDebug) << "event " << fEventNumber << endl;
00979                 crdhi.ResetFirst();
00980       
00981                 //reset the plane iterator
00982                 planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00983                 planeItr.GetSet()->Update();
00984                 planeItr.ResetFirst();
00985 
00986                 //get an UgliGeomHandle
00987                 UgliGeomHandle ugh(*crdhi.Ptr()->GetVldContext());
00988                 UgliStripHandle ush;
00989                                                                                                    
00990                 //make arrays to hold the strip hit information
00991                 Float_t chargeCorrect = 0., stripsCorrect = 0.;
00992                 Float_t totalCharge = 0., totalStrips = 0.;
00993                 Float_t totalCharge1Hyp = 0.;//, chargeCorrect1Hyp = 0.;
00994                 Float_t mcNum = 0.;
00995                 Float_t mcNumTPos = 0.;
00996                 Float_t mcDeNom = 0.;
00997 
00998                 Int_t mcStrip = 0;
00999                 Float_t mcCharge = 0;
01000                 Int_t strip = 0;
01001                 Float_t charge = 0.;
01002                 
01003                 //arrays to see if the truth has double ended strips
01004                 Int_t westStrips[251];
01005                 Int_t eastStrips[251];
01006                 
01007                 //loop over the planes in the event
01008                 while( planeItr.IsValid() ){
01009                         DmxPlane *currentPlane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
01010                         fPlane = currentPlane->GetPlaneNumber();
01011                         mcNum = 0.;
01012                         mcNumTPos = 0.;
01013                         mcDeNom = 0.;
01014                         fMCCoG = 0.;
01015                         mcStrip = -1;
01016                         mcCharge = 0.; 
01017                         strip = -1;
01018                         charge = 0.;
01019 
01020                         for(Int_t i = 0; i<251; ++i){
01021                                 eastStrips[i] = 0;
01022                                 westStrips[i] = 0;
01023                         }
01024                         
01025                         MSG("DmxDeMuxModule", Msg::kDebug) << "ana - event " << fEventNumber 
01026                                                                                           << " plane " << fPlane << endl;
01027                         //-------------------get digit info----------------------------------//
01028         
01029                         //slice for the digits in this plane
01030                         crdhi.GetSet()->Slice(fPlane);
01031                         crdhi.GetSet()->Update();
01032                         
01033                         PlexStripEndId trueSeid = truth.BestSEIdOfDigit(crdhi.Ptr());
01034 
01035                         //loop over the mc digits to get the center of gravity for the plane
01036                         while( (digit = crdhi()) ){
01037                                 trueSeid = truth.BestSEIdOfDigit(digit);
01038                                 
01039                                 //get the true stripend id for the digit
01040                                 if(trueSeid.GetStrip()>0 && trueSeid.GetStrip()<192){
01041                                         ush = ugh.GetStripHandle(trueSeid);
01042                                         mcNum += digit->GetCharge() * trueSeid.GetStrip();
01043                                         mcNumTPos += digit->GetCharge() * ush.GetTPos();
01044                                         mcDeNom += digit->GetCharge();
01045                                         if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast)
01046                                                 ++eastStrips[trueSeid.GetStrip()];
01047                                         if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest)
01048                                                 ++westStrips[trueSeid.GetStrip()];                                      
01049                                 }
01050                                 
01051                                 MSG("DmxDeMuxModule", Msg::kDebug) << "plane = " << digit->GetPlexSEIdAltL().GetPlane()
01052                                                                                                   << " digit charge = " << digit->GetCharge() << endl;
01053                         }
01054 
01055                         if(mcDeNom > 0.){
01056                                 fMCCoG = mcNum/mcDeNom;
01057                                 fMCCoGTPos = mcNumTPos/mcDeNom;
01058                         } 
01059                         crdhi.ResetFirst();
01060 
01061                         //loop over the digits and see if the demuxing is correct
01062                         while( (digit = crdhi()) ){
01063                 
01064                                 //get the digit and the strip info from the truth       
01065                                 mcCharge = digit->GetCharge();
01066                                 trueSeid = truth.BestSEIdOfDigit(digit);
01067                                 mcStrip = trueSeid.GetStrip();
01068 
01069                                 MSG("DeMuxModule", Msg::kDebug) << "digit on plane = " << fPlane
01070                                                                                            << " strip = " << mcStrip 
01071                                                                                            << " charge " << mcCharge 
01072                                                                                            << " is xtalk " << (int)truth.DigitIsOnlyCrosstalk(digit) << endl;
01073 
01074                                 //if the digit is cross-talk, ignore it
01075                                 fXTalkDigit = 0;
01076                                 if(truth.DigitIsOnlyCrosstalk(digit) ) fXTalkDigit=1;
01077                                 
01078                                 //if you are looking for the amount of signal correctly reconstructed
01079                                 //in the shower region, make sure it is really a shower region by
01080                                 //looking to see that there are at least 4 planes before the muon
01081                                 //starts
01082                                 if(!fXTalkDigit){
01083                                         if((fPlane<fMuonStartPlane && fMuonStartPlane != -1
01084                                                 && fMuonStartPlane-fVertexPlane>3) 
01085                                            || fMuonStartPlane == -1){
01086                                                 totalCharge += mcCharge;
01087                                                 if(TMath::Abs(fMCCoG-mcStrip)<=23) totalCharge1Hyp += mcCharge;
01088                                         }
01089                                         else if(fMuonStartPlane > -1 && fPlane>= fMuonStartPlane
01090                                                         && TMath::Abs(fMCCoG-mcStrip<=23)  
01091                                                         ) totalStrips += 1.;
01092                                         
01093                                 }//end if not cross-talk
01094 
01095                                 //now get the reco'ed info - the charge really should just be the same - it is the 
01096                                 //same digit afterall
01097                                 charge = digit->GetCharge();
01098                                 strip = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
01099                
01100                                 //the digit is demuxed correctly if the mc digit and demuxed digit
01101                                 //come from the same strip
01102                                 if(mcStrip == strip&&!fXTalkDigit){
01103               
01104                                         if((fPlane<fMuonStartPlane && fMuonStartPlane != -1
01105                                                 && fMuonStartPlane-fVertexPlane>3) 
01106                                            || fMuonStartPlane == -1){
01107                                                 chargeCorrect += charge;
01108                                         }
01109                                         else if(fMuonStartPlane != -1 && fPlane>= fMuonStartPlane
01110                                                         && TMath::Abs(fMCCoG-mcStrip<=23) ) stripsCorrect += 1.;
01111                                 }//end if the demuxed and mc digits are the same
01112 
01113                                 fDeMuxedStrip = mcStrip;
01114                                 fMCDigitCharge = mcCharge;
01115                                 fDoubleEnded = 0;
01116                                 if(eastStrips[mcStrip]&&westStrips[mcStrip]) fDoubleEnded = 1;
01117                                 if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest) fMCDigitSide = 1;
01118                                 else if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast) fMCDigitSide = 0;
01119           
01120                                 fMCDeMuxDigitTree->Fill();
01121                                 if ( fEventNumber%1000 == 0 ) {
01122                                         // save current working directory
01123                                         TDirectory* savedir = gDirectory;
01124                                         // move to file directory of tree, oddly SaveSelf doesn't work without it 
01125                                         fDeMuxFile->cd();  
01126                                         // save tree to file
01127                                         fMCDeMuxDigitTree->AutoSave();
01128                                         // save file directory containing this tree  
01129                                         fDeMuxFile->SaveSelf();
01130                                         // change back to working to directory
01131                                         savedir->cd();  
01132                                 }
01133           
01134                         }//end loop over digits
01135               
01136                         crdhi.GetSet()->ClearSlice(false);
01137                         crdhi.ResetFirst();
01138         
01139                         planeItr.Next();
01140                 } //end loop over planes in event
01141 
01142                 //get the linear fit for the first 5 planes
01143 //              util.FindLinearFit(planeZU, planeTruthCoGU, planeWeightU, 5, intercept, fMCInitialSlopeU, chiSqr);
01144 //              MSG("DmxDeMuxModule", Msg::kDebug) << "slope = " << fMCInitialSlopeU 
01145 //                                                                                << " intercept = " << intercept
01146 //                                                                                << " chiSqr = " << chiSqr << endl;
01147 //              util.FindLinearFit(planeZV, planeTruthCoGV, planeWeightV, 5, intercept, fMCInitialSlopeV, chiSqr);
01148 //              MSG("DmxDeMuxModule", Msg::kDebug) << "slope = " << fMCInitialSlopeV 
01149 //                                                                                << " intercept = " << intercept
01150 //                                                                                << " chiSqr = " << chiSqr << endl;
01151                 
01152       fMCChargeCorrect = -1.;
01153       fMCChargeCorrect1Hyp = -1.;
01154       fMCStripsCorrect = -1.;
01155      
01156                 if(totalCharge != 0.){
01157                   fMCChargeCorrect = chargeCorrect/totalCharge;
01158                   fMCTotalCharge = totalCharge;
01159                 }
01160                 if(totalCharge1Hyp != 0.){
01161                         fMCChargeCorrect1Hyp = chargeCorrect/totalCharge1Hyp;
01162                         if(fMCChargeCorrect1Hyp>1.) fMCChargeCorrect1Hyp=1.;
01163                                 fMCTotalCharge1Hyp = totalCharge1Hyp;
01164                 }
01165                 if(totalStrips != 0.) fMCStripsCorrect = stripsCorrect/totalStrips;
01166                 if( fMCStripsCorrect>1.) fMCStripsCorrect = 1;
01167 
01168       //cout << "filling mc demux tree" << endl;
01169       fMCDeMuxTree->Fill();
01170       if ( fEventNumber%1000 == 0 ) {
01171                   TDirectory* savedir = gDirectory; // save current working directory
01172                   fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
01173                   fMCDeMuxTree->AutoSave();  // save tree to file
01174                   fDeMuxFile->SaveSelf();  // save file directory containing this tree
01175                   savedir->cd();  // change back to working to directory
01176       }
01177 
01178       planeItr.GetSet()->ClearSlice(false);
01179 
01180     } //end if MonteCarlo data
01181     
01182     //fill the DeMux info Tree
01183 
01184     fDeMuxTree->Fill();
01185     if ( fEventNumber%1000 == 0 ) {
01186       TDirectory* savedir = gDirectory; // save current working directory
01187       fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
01188       fDeMuxTree->AutoSave();  // save tree to file
01189       fDeMuxFile->SaveSelf();  // save file directory containing this tree
01190       savedir->cd();  // change back to working to directory
01191     }
01192 
01193   }
01194   //MSG("JobCDmx", Msg::kDebug) << "finished ana" << endl;
01195 
01196   return result;
01197 }
01198 
01199 //......................................................................
01200 
01201 void DmxDeMuxModule::HandleCommand(JobCommand * /* command */) 
01202 {
01203   //
01204   //  Purpose:  Method to interpret module commands.
01205   //
01206   //  Arguments:
01207   //    command   in    Command to interpret.
01208   //
01209   //  Return:   n/a
01210   //
01211   //  Commands implemented:
01212   //    Set TestToMake  Set the type of test to make - beam, cosmics, or xtalk
01213   //    Set DataType    Set the type of data used - far or mc
01214   //
01215 
01216    MSG("JobC", Msg::kDebug) << "DeMuxModule::HandleCommand" << endl;
01217 
01218   //  TString cmd = command->PopCmd();
01219 //    if(cmd == "Set"){   
01220 //      TString opt = command->PopOpt();
01221 //      if(opt == "TestToMake"){ 
01222 //        TString testopt = command->PopOpt();
01223 //        if(testopt == "MCXTalkEvents"){ fMCXTalkEvents = true;}
01224 //        if(testopt == "MCUpAndDownStream"){ fMCUpAndDownStream = true;}
01225 //        if(testopt == "MatedSignalFraction"){ fMatedSignalFraction = true;}
01226 //      }
01227 //      else if(opt == "DataType"){ fDataType = command->PopOpt();}
01228 //      else if(opt == "DontUseCandDigitMasks"){ fDontUseCandDigitMasks = true; }
01229 //      else if(opt == "Filter"){ fFilter = true; }
01230 //      else {
01231 //          MSG("JobC", Msg::kWarning)<< "DeMuxModule: Unrecognized option " << opt << endl;
01232 //       }
01233 //    }
01234 //    else {
01235 //      MSG("JobC", Msg::kWarning)<< "DeMuxModule: Unrecognized command " << cmd << endl;
01236 //    }
01237 
01238    return;
01239 }
01240 
01241 //......................................................................
01242 
01243 const Registry& DmxDeMuxModule::DefaultConfig() const
01244 {
01245   //  see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html 
01246   //for details
01247    
01248   int itrue = 1;  // work around for lack of bool in registry
01249   int ifalse = 0; // work around for lack of bool in registry
01250   
01251   MSG("JobC", Msg::kDebug) << "DeMuxCosmicsModule::DefaultConfig" << endl;
01252   
01253   static Registry r;
01254   
01255   r.UnLockValues();
01256   
01257   r.Set("DataType",            0);
01258   r.Set("WriteHistos",          itrue);
01259   r.Set("MakeTrees",            ifalse);
01260   r.Set("HistoFileName",        "demux_beam.root");
01261   r.Set("Filter",               ifalse);
01262   r.Set("UseCandDigitMasks",    itrue);
01263   r.Set("WindowSize",           5);
01264   r.Set("StrayDeltaStripLimit", 6);
01265   r.Set("StrayPlaneLimit",      4);
01266   r.Set("UseFigureOfMerit",     ifalse);
01267   r.Set("HoughSlopeRMS",        5.);
01268   r.Set("HoughInterceptRMS",    10.);
01269   r.Set("HoughPeak",            0.66);
01270   r.Set("XTalkFractionLimit",  0.1);
01271   r.Set("RatioMatedSignalForValid",  0.5);
01272   r.LockValues();
01273   
01274   return r;
01275 }
01276 
01277 //......................................................................
01278 
01279 void DmxDeMuxModule::Config(const Registry& r)
01280 {
01281   int         tmpb;  // a temp bool. See comment under DefaultConfig...
01282   int         tmpi;  // a temp int.
01283   double      tmpd;  // a temp double.
01284   const char* tmps;  // a temp string
01285 
01286   if (r.Get("SingleModuleDemux",   tmpb)) fSingleModuleDeMux     =  tmpb;
01287   if (r.Get("EventLength",         tmpb)) fEventLength           =  tmpb;
01288   if (r.Get("UnsetDigitCheck",     tmpb)) fUnsetDigitCheck       =  tmpb;
01289   if (r.Get("UseCandDigitMasks",   tmpb)) fDontUseCandDigitMasks = !tmpb;
01290   if (r.Get("WriteHistos",         tmpb)) fWriteHistos           =  tmpb;
01291   if (r.Get("MakeTrees",           tmpb)) fMakeTrees             =  tmpb;
01292   if (r.Get("HistoFileName",       tmps)) fHistoFileName         =  tmps;
01293   if (r.Get("DataType",            tmpi)) fDataType              =  tmpi;
01294   if (r.Get("WindowSize",          tmpi)) fPlanesInWindow        =  tmpi;
01295   if (r.Get("StrayDeltaStripLimit",tmpi)) fStrayDeltaStripLimit  =  tmpi;
01296   if (r.Get("StrayPlaneLimit",     tmpi)) fStrayPlanesLimit      =  tmpi;
01297   if (r.Get("UseFigureOfMerit",    tmpb)) fUseStrayPlaneCheck    =  tmpb;
01298   if (r.Get("Filter",              tmpb)) fFilter                =  tmpb;
01299   if (r.Get("HoughSlopeRMS",       tmpd)) fHoughSlope            =  tmpd;
01300   if (r.Get("HoughInterceptRMS",   tmpd)) fHoughIntercept        =  tmpd;
01301   if (r.Get("HoughPeak",           tmpd)) fHoughPeak             =  tmpd;
01302   if (r.Get("XTalkFractionLimit",  tmpd)) fSignalFractionLimit   =  tmpd;
01303   if (r.Get("RatioMatedSignalForValid",  tmpd)) fMatedSignalForValid   =  tmpd;
01304 
01305 }
01306 
01307 
01308 //......................................................................
01309 void DmxDeMuxModule::Help() 
01310 {
01311   MSG("JobC", Msg::kDebug) 
01312     << "DeMuxModule::Help\n"
01313     <<"DmxDeMuxModule is a module which demultiplexes events "
01314     <<"in the far detector."
01315     << endl;
01316 }
01317 
01318 void DmxDeMuxModule::EndJob() 
01319 {
01320   //save the trees to a file
01321   if(fMakeTrees && fDeMuxFile){
01322     
01323     fDeMuxFile->Write();
01324     fDeMuxFile->Close();
01325   }
01326   return;
01327 }
01328 
01329 
01330 
01331 
01332 
01333 
01334