loon: PreFilter.cxx Source File

00001 
00002 // $Id: PreFilter.cxx,v 1.7 2007/02/01 21:50:27 rhatcher Exp $
00003 //
00004 // PreFilter
00005 //
00006 // Purpose: PreFilter is a containment-based data filter.
00007 //
00009 #include "PreFilter.h"
00010 
00011 #include <map>
00012 #include <cmath>
00013 
00014 #include "MinosObjectMap/MomNavigator.h"
00015 #include "JobControl/JobCModuleRegistry.h"
00016 #include "CandData/CandRecord.h"
00017 #include "CandData/CandHeader.h"
00018 #include "CandDigit/CandDigitHandle.h"
00019 #include "CandDigit/CandDigitListHandle.h"
00020 #include "Plex/PlexSEIdAltL.h"
00021 #include "Plex/PlexStripEndId.h"
00022 #include "UgliGeometry/UgliGeomHandle.h"
00023 #include "UgliGeometry/UgliStripHandle.h"
00024 #include "Validity/VldContext.h"
00025 #include "RawData/RawRecord.h"
00026 #include "RawData/RawDaqSnarlHeader.h"
00027 #include "RawData/RawDaqHeaderBlock.h"
00028 #include "Digitization/DigiScintHit.h"
00029 #include "Record/SimSnarlRecord.h"
00030 #include "REROOT_Classes/REROOT_NeuKin.h"
00031 #include "MessageService/MsgService.h"
00032 
00033 #include "TFile.h"
00034 #include "TTree.h"
00035 #include "TClonesArray.h"
00036 
00037 JOBMODULE(PreFilter,"PreFilter","");
00038 CVSID("$Id: PreFilter.cxx,v 1.7 2007/02/01 21:50:27 rhatcher Exp $");
00039 
00040 // Begin the code............................................................
00041 PreFilter::PreFilter() :
00042     cDiagnosticMode(0),
00043     cSnrlChrgCut(-1.),
00044     cUVdistcut(-1.),
00045     cXYUVcut(-1.),
00046     cMinNPlanes(-1),
00047     cListIn("demuxdigitlist"),
00048 
00049     fSnrlChrgDecD(-1),
00050     fContainDecD(-1),
00051     fNPlanesDecD(-1),
00052     fSnrlChrgDecT(-1),
00053     fContainDecT(-1),
00054     fNPlanesDecT(-1),
00055     
00056     fProblems(-1),
00057     
00058     cPEcut(-1.),
00059     cDmxWtcut(-1.),
00060     cClustDistcut(-1.),
00061     cDigFitDcut(-1.),
00062     
00063     fMinPln(-1),
00064     fMaxPln(-1),
00065     fSnrlZMin(-1.),
00066     fSnrlZMax(-1.),
00067     fRun(-1),
00068     fSubrun(-1),
00069     fSnrl(-1),
00070     fNRawDig(-1),
00071     fNCandDig(-1),
00072     fNPlanes(-1),
00073     fNTruPlanes(-1),
00074     fUdigSize(-1),
00075     fVdigSize(-1),
00076     fNUClusters(0),
00077     fNVClusters(0),
00078     fNOLClusters(0),
00079     fTotUClustDig(0),
00080     fTotVClustDig(0),
00081     fNSnarls(0),
00082     
00083     fSnrlChrg(0.),
00084     fUDigChrg(0.),
00085     fVDigChrg(0.),
00086     fTotUClustChrg(0.),
00087     fTotVClustChrg(0.),
00088     fTotTruE(0.)
00089 {
00090     MSG("PreFilter",Msg::kVerbose) << "constructor" << std::endl;
00091 }
00092 //.............................................................................
00093 PreFilter::~PreFilter()
00094 {
00095     // cut variables
00096     if(cDiagnosticMode){
00097       MyTree->Branch("cPEcut",&cPEcut,"cPEcut/F");
00098       MyTree->Branch("cDmxWtcut",&cDmxWtcut,"cDmxWtcut/F");  
00099       MyTree->Branch("cUVdistcut",&cUVdistcut,"cUVdistcut/F");  
00100       MyTree->Branch("cClustDistcut",&cClustDistcut,"cClustDistcut/F");  
00101       MyTree->Branch("cXYUVcut",&cXYUVcut,"cXYUVcut/F");
00102       MyTree->Branch("cDigFitDcut",&cDigFitDcut,"cDigFitDcut/F");  
00103       MyTree->Branch("cSnrlChrgCut",&cSnrlChrgCut,"cSnrlChrgCut/F");  
00104       MyTree->Branch("cMinNPlanes",&cMinNPlanes,"cMinNPlanes/I");  
00105       
00106       MyTree->Fill();
00107       outfile->cd();
00108       MyTree->Write();
00109       MyTree->Delete();
00110       outfile->Close();
00111     }
00112     delete MyTree;
00113     delete outfile;
00114 }
00115 //.............................................................................
00116 void PreFilter::BeginJob(){
00117 
00118     MSG("PreFilter",Msg::kVerbose) << "BeginJob" << std::endl;
00119   
00120     if(cDiagnosticMode){
00121       outfile = new TFile("outfile","RECREATE");
00122       MyTree = new TTree("MyTree","MyTree");
00123       
00124       MyTree->Branch("fSnrlChrgDecD",&fSnrlChrgDecD,"fSnrlChrgDecD/I");
00125       MyTree->Branch("fContainDecD",&fContainDecD,"fContainDecD/I");
00126       MyTree->Branch("fNPlanesDecD",&fNPlanesDecD,"fNPlanesDecD/I");
00127       MyTree->Branch("fSnrlChrgDecT",&fSnrlChrgDecT,"fSnrlChrgDecT/I");
00128       MyTree->Branch("fContainDecT",&fContainDecT,"fContainDecT/I");
00129       MyTree->Branch("fNPlanesDecT",&fNPlanesDecT,"fNPlanesDecT/I");
00130       
00131       MyTree->Branch("fProblems",&fProblems,"fProblems/I");
00132       
00133       MyTree->Branch("fMinPln",&fMinPln,"fMinPln/I");
00134       MyTree->Branch("fMaxPln",&fMaxPln,"fMaxPln/I");
00135       MyTree->Branch("fRun",&fRun,"fRun/I");
00136       MyTree->Branch("fSubrun",&fSubrun,"fSubrun/I");
00137       MyTree->Branch("fSnrl",&fSnrl,"fSnrl/I");
00138       MyTree->Branch("fNRawDig",&fNRawDig,"fNRawDig/I");
00139       MyTree->Branch("fNCandDig",&fNCandDig,"fNCandDig/I");
00140       MyTree->Branch("fNPlanes",&fNPlanes,"fNPlanes/I");
00141       MyTree->Branch("fNTruPlanes",&fNTruPlanes,"fNTruPlanes/I");
00142       MyTree->Branch("fUdigSize",&fUdigSize,"fUdigSize/I");
00143       MyTree->Branch("fVdigSize",&fVdigSize,"fVdigSize/I");
00144       MyTree->Branch("fNUClusters",&fNUClusters,"fNUClusters/I");
00145       MyTree->Branch("fNVClusters",&fNVClusters,"fNVClusters/I");
00146       MyTree->Branch("fNOLClusters",&fNOLClusters,"fNOLClusters/I");
00147       MyTree->Branch("fTotUClustDig",&fTotUClustDig,"fTotUClustDig/I");  
00148       MyTree->Branch("fTotVClustDig",&fTotVClustDig,"fTotVClustDig/I");  
00149       MyTree->Branch("fNSnarls",&fNSnarls,"fNSnarls/I");
00150       MyTree->Branch("fSnrlChrg",&fSnrlChrg,"fSnrlChrg/F");  
00151       MyTree->Branch("fUDigChrg",&fUDigChrg,"fUDigChrg/F");  
00152       MyTree->Branch("fVDigChrg",&fVDigChrg,"fVDigChrg/F");  
00153       MyTree->Branch("fTotUClustChrg",&fTotUClustChrg,"fTotUClustChrg/F");  
00154       MyTree->Branch("fTotVClustChrg",&fTotVClustChrg,"fTotVClustChrg/F");
00155       MyTree->Branch("fTotTruE",&fTotTruE,"fTotTruE/F");
00156       MyTree->Branch("fUClustSize",fUClustSize,"fUClustSize[fNUClusters]/I");
00157       MyTree->Branch("fVClustSize",fVClustSize,"fVClustSize[fNVClusters]/I");
00158       MyTree->Branch("fUClustChrg",fUClustChrg,"fUClustChrg[fNUClusters]/F");
00159       MyTree->Branch("fVClustChrg",fVClustChrg,"fVClustChrg[fNVClusters]/F");
00160       MyTree->Branch("fUClustMinZ",fUClustMinZ,"fUClustMinZ[fNUClusters]/F");
00161       MyTree->Branch("fVClustMinZ",fVClustMinZ,"fVClustMinZ[fNVClusters]/F");
00162       MyTree->Branch("fUClustMaxZ",fUClustMaxZ,"fUClustMaxZ[fNUClusters]/F");
00163       MyTree->Branch("fVClustMaxZ",fVClustMaxZ,"fVClustMaxZ[fNVClusters]/F");
00164       MyTree->Branch("fUSlope",fUSlope,"fUSlope[fNUClusters]/F");
00165       MyTree->Branch("fUIntercept",fUIntercept,"fUIntercept[fNUClusters]/F");
00166       MyTree->Branch("fVSlope",fVSlope,"fVSlope[fNVClusters]/F");
00167       MyTree->Branch("fVIntercept",fVIntercept,"fVIntercept[fNVClusters]/F");
00168       MyTree->Branch("fP4Neu",fP4Neu,"fP4Neu[4]/F");
00169     }
00170 }
00171 //.............................................................................
00172 JobCResult PreFilter::Ana(const MomNavigator *mom)
00173 {
00174     JobCResult result(JobCResult::kFailed);
00175     result.SetWarning();
00176 
00177     if(cDiagnosticMode) fNSnarls++;
00178   
00179     this->BeginSnarl();
00180   
00181     if(cDiagnosticMode) this->SimCheck(mom);
00182     
00183     std::vector<CandDigitHandle*> uDigits;
00184     std::vector<CandDigitHandle*> vDigits;
00185     
00186     const RawRecord *rawrec = dynamic_cast<const RawRecord*>
00187       (mom->GetFragment("RawRecord"));
00188     
00189     if(!rawrec){
00190       MSG("PreFilter",Msg::kError) << "No rawrecord in mom!\n";
00191       return result;
00192       //return JobCResult::kError;
00193     }
00194     
00195     // part to get minimum and maximum planes for run...///////////////////
00196     VldContext snrlvld = rawrec->GetRawHeader()->GetVldContext();
00197     UgliGeomHandle snrlugh(snrlvld);
00198     snrlugh.GetZExtent(fSnrlZMin,fSnrlZMax,-1);
00199     fMinPln = snrlugh.GetPlaneIdFromZ(fSnrlZMin).GetPlane();
00200     fMaxPln = snrlugh.GetPlaneIdFromZ(fSnrlZMax).GetPlane();
00202     
00203     if(cDiagnosticMode){
00204       const RawDaqSnarlHeader* rdsh = dynamic_cast<const RawDaqSnarlHeader*>
00205         (rawrec->GetRawHeader());
00206       
00207       if(rdsh){
00208         fSnrl = rdsh->GetSnarl();
00209         fNRawDig = rdsh->GetNumRawDigits();
00210       }
00211       const RawDaqHeaderBlock* rdhb = dynamic_cast<const RawDaqHeaderBlock*>
00212         (rawrec->FindRawBlock("RawDaqHeaderBlock"));
00213       
00214       if(rdhb){
00215         fRun = rdhb->GetRun();
00216         fSubrun = rdhb->GetSubRun();
00217       }
00218     }
00219 
00220     CandRecord *candrec = dynamic_cast<CandRecord *>
00221       (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00222     
00223     if(!candrec){
00224         MSG("PreFilter",Msg::kError) << "No candrecord in mom" << std::endl;
00225         return result;
00226         //return JobCResult::kError;
00227     }
00228     
00229     const CandDigitListHandle *cdlh = dynamic_cast<CandDigitListHandle*>
00230       (candrec->FindCandHandle("CandDeMuxDigitListHandle",cListIn.c_str()));
00231     
00232     if(!cdlh){
00233       MSG("PreFilter",Msg::kError) << "CandDigitListHandle Empty" << std::endl;
00234       return result;
00235       //return JobCResult::kError;
00236     }
00237 
00238     if(cDiagnosticMode){
00239       fNCandDig = cdlh->GetNDaughters();
00240     }
00241 
00242     CandDigitHandleItr cdhiter(cdlh->GetDaughterIterator());
00243     std::vector<int> tmpplanes;
00244     // loop to sum charge in snarl and count planes
00245     for(;cdhiter.IsValid();cdhiter.Next()){
00246       const PlexSEIdAltL &pseidalt = dynamic_cast<const PlexSEIdAltL&>
00247         ((*cdhiter)->GetPlexSEIdAltL());
00248       
00249       if(!pseidalt.IsVetoShield()){
00250         fSnrlChrg += (*cdhiter)->GetCharge(CalDigitType::kPE);
00251         
00252         if(tmpplanes.size()==0){
00253           tmpplanes.push_back(pseidalt.GetPlane());
00254         }
00255         else{
00256           int pln = pseidalt.GetPlane();
00257           std::vector<int>::iterator plnisfound = 
00258             std::find(tmpplanes.begin(),tmpplanes.end(),pln);
00259           if(plnisfound==tmpplanes.end()) tmpplanes.push_back(pln);
00260         }
00261       }
00262     }
00263     if(fSnrlChrg<cSnrlChrgCut){
00264       if(cDiagnosticMode){
00265           MSG("PreFilter",Msg::kError) << "Snarl Charge too low, non-signal"
00266           << " event" << std::endl;
00267         fSnrlChrgDecD = 0;
00268         MyTree->Fill();
00269       }
00270       else{
00271           MSG("PreFilter",Msg::kError) << "Snarl Charge too low, non-signal"
00272           << " event" << std::endl;
00273           return result;
00274           //return JobCResult::kError;
00275       }
00276     }
00277     else{
00278       if(cDiagnosticMode){
00279         fSnrlChrgDecD = 1;
00280       }
00281     }
00282 
00283     fNPlanes = tmpplanes.size();
00284     if(fNPlanes<cMinNPlanes){
00285       if(cDiagnosticMode){
00286         MSG("PreFilter",Msg::kWarning) << "Not enough planes" << std::endl;
00287         fNPlanesDecD = 0;
00288         MyTree->Fill();
00289       }
00290       else{
00291         MSG("PreFilter",Msg::kWarning) << "Not enough planes" << std::endl;
00292         return result;
00293         //return JobCResult::kError;
00294       }
00295     }
00296     else{
00297       if(cDiagnosticMode){
00298         fNPlanesDecD = 1;
00299       }
00300     }
00301     cdhiter.Reset();
00302     int nbaddmx = 0;
00303     for(;cdhiter.IsValid();cdhiter.Next()){
00304       // check to see if digit is greater than 3pe and kludge for now with
00305       // DeMux veto flag.
00306       // check to see if digit is contained in UV and demuxed well
00307       const PlexSEIdAltL &pseidalt = dynamic_cast<const PlexSEIdAltL&>
00308         ((*cdhiter)->GetPlexSEIdAltL());
00309       if(!pseidalt.IsVetoShield()){
00310         // check demux flag for IU demuxer?
00311         if(pseidalt.GetDemuxVetoFlag()==0){
00312           if((*cdhiter)->GetCharge(CalDigitType::kPE) > cPEcut){
00313             // check to see if digit is in first or last planes
00314             if(pseidalt.GetPlane()<=(fMinPln+2) ||
00315                pseidalt.GetPlane()>=(fMaxPln-2)){
00316               if(cDiagnosticMode){
00317                 MSG("PreFilter",Msg::kError) << "Event at end of detector" <<
00318                 std::endl;
00319                 fContainDecD = 0;
00320                 MyTree->Fill();
00321               }
00322               else{
00323                MSG("PreFilter",Msg::kError) << "Event at end of detector" <<
00324                std::endl;
00325                return result;
00326               }
00327             }
00328             if(pseidalt.GetBestWeight()>cDmxWtcut){
00329               //pseidalt.GetDemuxVetoFlag()==0){
00330               UgliGeomHandle ugh(*(*cdhiter)->GetVldContext());
00331               UgliStripHandle ush = ugh.GetStripHandle(pseidalt.GetBestSEId());
00332               
00333               if( (ush.GetTPos()*ush.GetTPos()) > cUVdistcut*cUVdistcut){
00334                 if(cDiagnosticMode){
00335                   MSG("PreFilter",Msg::kError) << "UV containment not met." 
00336                   << std::endl;
00337                   fContainDecD = 0;
00338                   MyTree->Fill();
00339                 }
00340                 else{
00341                   MSG("PreFilter",Msg::kError) << "UV containment not met." 
00342                   << std::endl;
00343                   return result;
00344                 }
00345               }
00346               // if get to this point, charge and UV containment met. need to 
00347               // put digit in vector for later use.
00348               switch(pseidalt.GetPlaneView(true)){
00349               case PlaneView::kU:
00350                 fUDigChrg += (*cdhiter)->GetCharge(CalDigitType::kPE);
00351                 uDigits.push_back((*cdhiter));
00352                 break;
00353                 
00354               case PlaneView::kV:
00355                 fVDigChrg += (*cdhiter)->GetCharge(CalDigitType::kPE);
00356                 vDigits.push_back((*cdhiter));
00357                 break;
00358                 
00359               default:
00360                 MSG("PreFilter",Msg::kWarning) << "SEIdAltL not U or V.";
00361                 MSG("PreFilter",Msg::kWarning) << std::endl;
00362               }
00363             }
00364           }
00365         }
00366         else{
00367             nbaddmx++;
00368             if( ((float)nbaddmx)/((float)fNRawDig) > 0.3){
00369               if(cDiagnosticMode){
00370                 MSG("PreFilter",Msg::kError) << "TOO MANY DEMUX PROBLEMS"
00371                 << std::endl;
00372                 fProblems = 1;
00373                 MyTree->Fill();
00374             }
00375             else{
00376               MSG("PreFilter",Msg::kError) << "TOO MANY DEMUX PROBLEMS"
00377               << std::endl;
00378               return result;
00379             }
00380           }
00381         }
00382        }
00383       }
00384     
00385     fUdigSize = uDigits.size();
00386     fVdigSize = vDigits.size();
00387     
00388     if(fUdigSize<=1 || fVdigSize<=1){
00389       if(cDiagnosticMode){
00390         MSG("PreFilter",Msg::kError) << "Not enough digits for clustering" <<
00391         std::endl;
00392         fProblems = 2;
00393         MyTree->Fill();
00394         return result;
00395       }
00396       else{
00397         MSG("PreFilter",Msg::kError) << "Not enough digits for clustering" <<
00398         std::endl;
00399         return result;
00400       }
00401     }
00402     
00403     std::vector<std::vector<CandDigitHandle*> >
00404       Vclusters(this->Cluster(vDigits));
00405     std::vector<std::vector<CandDigitHandle*> >
00406       Uclusters(this->Cluster(uDigits));
00407     
00408     fNUClusters = Uclusters.size();
00409     fNVClusters = Vclusters.size();
00410     
00411     if(fNUClusters==0 || fNVClusters==0){
00412       if(cDiagnosticMode){
00413         MSG("PreFilter",Msg::kError) << "Clustering failed" << std::endl;
00414         fProblems = 3;
00415         MyTree->Fill();
00416       }
00417       else{
00418        MSG("PreFilter",Msg::kError) << "Clustering failed" << std::endl;
00419        return result;
00420       }
00421     }
00422     
00423     for(unsigned int i=0;i<Vclusters.size();i++){
00424       fVSlope[i] = this->LSFSlope(Vclusters[i]);
00425       fVIntercept[i] = this->LSFIntercept(Vclusters[i]);
00426     }
00427     for(unsigned int i=0;i<Uclusters.size();i++){
00428       fUSlope[i] = this->LSFSlope(Uclusters[i]);
00429       fUIntercept[i] = this->LSFIntercept(Uclusters[i]);
00430     }
00431     
00432     std::vector<std::pair<int,int> >OLclusters(this->ClusterMatch(Uclusters,
00433                                                                   Vclusters));
00434     fNOLClusters = OLclusters.size();
00435     
00436     if(cDiagnosticMode) this->Diagnostic(Uclusters,Vclusters);
00437     
00438     for(unsigned int i=0;i<OLclusters.size();i++){
00439       std::vector<CandDigitHandle*> Ucluster(Uclusters[OLclusters[i].first]);
00440       std::vector<CandDigitHandle*> Vcluster(Vclusters[OLclusters[i].second]);
00441       
00442       // containment check. loop over U and V clusters and check XY containment
00443       // with other view's fit.
00444       
00445       for(std::vector<CandDigitHandle*>::iterator Vclustiter =
00446             Vcluster.begin(); Vclustiter!=Vcluster.end(); Vclustiter++){
00447         
00448         float vpos = this->GetCDHTPos((*Vclustiter));
00449         float zpos = this->GetCDHZPos((*Vclustiter));
00450         
00451         // check to see whether this particular digit is close to the
00452         // fit in its view. if not, don't bother checking containment with it.
00453         if( fabs(vpos - (fVSlope[OLclusters[i].second]*zpos + 
00454                          fVIntercept[OLclusters[i].second])) < cDigFitDcut){
00455           
00456           float upos = fUSlope[OLclusters[i].first]*zpos + 
00457             fUIntercept[OLclusters[i].first];
00458           float xpos = 0.707107*(upos-vpos);
00459           float ypos = 0.707107*(upos+vpos);
00460           
00461           if( ((upos*upos) > (cXYUVcut*cXYUVcut)) || 
00462               ((vpos*vpos) > (cXYUVcut*cXYUVcut)) ||
00463               ((xpos*xpos) > (cXYUVcut*cXYUVcut)) ||
00464               ((ypos*ypos) > (cXYUVcut*cXYUVcut)) ) {
00465             if(cDiagnosticMode){
00466               MSG("PreFilter",Msg::kError)<<"Containment failed" << std::endl;
00467               fContainDecD = 0;
00468               MyTree->Fill();
00469             }
00470             else{
00471              MSG("PreFilter",Msg::kError)<<"Containment failed" << std::endl;
00472              return result;
00473             }
00474           }
00475         }
00476       }
00477       for(std::vector<CandDigitHandle*>::iterator Uclustiter = 
00478             Ucluster.begin(); Uclustiter!=Ucluster.end(); Uclustiter++){
00479         
00480         float upos = this->GetCDHTPos((*Uclustiter));
00481         float zpos = this->GetCDHZPos((*Uclustiter));
00482         
00483         if( fabs(upos - (fUSlope[OLclusters[i].first]*zpos +
00484                          fUIntercept[OLclusters[i].first])) < cDigFitDcut){
00485           
00486           float vpos = fVSlope[OLclusters[i].second]*zpos + 
00487             fVIntercept[OLclusters[i].second];
00488           float xpos = 0.707107*(upos-vpos);
00489           float ypos = 0.707107*(upos+vpos);
00490           
00491           if( ((upos*upos) > (cXYUVcut*cXYUVcut)) ||
00492               ((vpos*vpos) > (cXYUVcut*cXYUVcut)) ||
00493               ((xpos*xpos) > (cXYUVcut*cXYUVcut)) ||
00494               ((ypos*ypos) > (cXYUVcut*cXYUVcut)) ) {
00495             if(cDiagnosticMode){
00496               MSG("PreFilter",Msg::kError)<<"Containment failed" << std::endl;
00497               fContainDecD = 0;
00498               MyTree->Fill();
00499             }
00500             else{
00501              MSG("PreFilter",Msg::kError)<<"Containment failed" << std::endl;
00502              return result;
00503             }
00504           }
00505         }
00506       }
00507     }
00508     
00509     //fProblems = 0;
00510     if(cDiagnosticMode){
00511         if(fContainDecD==-1) fContainDecD = 1;
00512         MyTree->Fill();
00513     }
00514     result.SetAOK().SetPassed();
00515     return result;
00516 }
00517 //.............................................................................
00518 const Registry& PreFilter::DefaultConfig() const { 
00519     static Registry r;
00520   
00521     std::string name = this->JobCModule::GetName();
00522     name += ".config.default";
00523     r.SetName(name.c_str());
00524     r.UnLockValues();
00525     
00526     // Set values in configuration
00527     r.Set("PEcut",3.0);
00528     r.Set("DmxWtcut",0.4);
00529     r.Set("UVdistcut",3.75);
00530     r.Set("ClustDistcut",0.60);
00531     r.Set("XYUVcut",3.75);
00532     r.Set("DigFitDcut",1.0);
00533     r.Set("SnrlChrgCut",50.);
00534     r.Set("MinNPlanes",6);
00535     r.Set("NameListIn","demuxdigitlist");
00536     r.Set("DiagnosticMode",0);
00537     r.LockValues();
00538     
00539     return r;
00540 }
00541 //.............................................................................
00542 void PreFilter::Config(const Registry &r) {
00543 
00544   MSG("PreFilter",Msg::kVerbose) << "We're in Config now." << std::endl;
00545 
00546     double tmpd;
00547 
00548     if(r.Get("PEcut",tmpd)) {cPEcut = tmpd;}
00549     if(r.Get("DmxWtcut",tmpd)) {cDmxWtcut = tmpd;}
00550     if(r.Get("UVdistcut",tmpd)) {cUVdistcut = tmpd;}
00551     if(r.Get("ClustDistcut",tmpd)) {cClustDistcut = tmpd;}
00552     if(r.Get("XYUVcut",tmpd)) {cXYUVcut = tmpd;}
00553     if(r.Get("DigFitDcut",tmpd)) {cDigFitDcut = tmpd;}
00554     if(r.Get("SnrlChrgCut",tmpd)) {cSnrlChrgCut = tmpd;}
00555     
00556     int tmpi;
00557     if(r.Get("MinNPlanes",tmpi)) {cMinNPlanes = tmpi;}
00558     if(r.Get("DiagnosticMode",tmpi)) {cDiagnosticMode = (bool)tmpi;}
00559 
00560     const char *tmps;
00561     if (r.Get("NameListIn",tmps)) { cListIn = tmps; }
00562 }
00563 //.............................................................................
00564 float PreFilter::LSFSlope(std::vector<CandDigitHandle*>& cluster){
00565   
00566     int size = cluster.size();
00567     
00568     std::vector<float> zpos(size,0.);
00569     std::vector<float> tpos(size,0.);
00570     
00571     for(int i=0;i<size;i++){
00572       tpos[i] = this->GetCDHTPos(cluster[i]);
00573       zpos[i] = this->GetCDHZPos(cluster[i]);
00574     }
00575     
00576     float x = 0.;
00577     float xsqr = 0.;
00578     float y = 0.;
00579     float xy = 0.;
00580     
00581     for(int i=0;i<size;i++){
00582       x+=zpos[i];
00583       xsqr+=zpos[i]*zpos[i];
00584       y+=tpos[i];
00585       xy+=zpos[i]*tpos[i];
00586     }
00587     
00588     float numerator = (x*y - size*xy);
00589     float denominator = (x*x - size*xsqr);
00590     
00591     if(denominator!=0.){
00592       float slope = numerator/denominator;
00593       return slope;
00594     }
00595     else return 0.;
00596 }
00597 //.............................................................................
00598 float PreFilter::LSFIntercept(std::vector<CandDigitHandle*>& cluster){
00599 
00600     int size = cluster.size();
00601 
00602     std::vector<float> zpos(size,0.);
00603     std::vector<float> tpos(size,0.);
00604     
00605     for(int i=0;i<size;i++){
00606       tpos[i] = this->GetCDHTPos(cluster[i]);
00607       zpos[i] = this->GetCDHZPos(cluster[i]);
00608     }
00609     
00610     float x = 0.;
00611     float xsqr = 0.;
00612     float y = 0.;
00613     float xy = 0.;
00614     
00615     for(int i=0;i<size;i++){
00616       x+=zpos[i];
00617       xsqr+=zpos[i]*zpos[i];
00618       y+=tpos[i];
00619       xy+=zpos[i]*tpos[i];
00620     }
00621     
00622     float numerator = (xsqr*y - x*xy);
00623     float denominator = (size*xsqr - x*x);
00624     
00625     if(denominator!=0.){
00626       float intercept = numerator/denominator;
00627       return intercept;
00628     }
00629     else return 0.;
00630 }
00631 //.............................................................................
00632 std::vector<std::vector<CandDigitHandle*> >
00633 PreFilter::Cluster(std::vector<CandDigitHandle*> &digits){
00634 
00635     int ndigits = digits.size();
00636     
00637     // initialize dst matrix
00638     std::vector<std::vector<float> > dst (ndigits);
00639     for(int i=0;i<ndigits;i++){
00640       dst[i] = std::vector<float>(ndigits);
00641     }
00642     // initialize group array. this holds the group assignment for each digit.
00643     // -2 => digit hasn't been assigned group yet. -1 => digit didn't fit into 
00644     // any group 
00645     std::vector<int> group(ndigits,-2);
00646     // initialize map of group number/size of group
00647     std::map<int,int> groupsize;
00648     
00649     // fill dst matrix  
00650     for(int i=0;i<ndigits;i++){
00651       
00652       float itpos = this->GetCDHTPos(digits[i]);
00653       float izpos = this->GetCDHZPos(digits[i]);
00654       
00655       for(int j=0;j<ndigits;j++){
00656         
00657         float jtpos = this->GetCDHTPos(digits[j]);
00658         float jzpos = this->GetCDHZPos(digits[j]);
00659         
00660         dst[i][j] = sqrt((itpos-jtpos)*(itpos-jtpos) + 
00661                          (izpos-jzpos)*(izpos-jzpos) );
00662       }
00663     }
00664     
00665     int igroup = 0;
00666     int i0 = 0;
00667     int j0 = 0;
00668     int nleft = ndigits-1;
00669     group[0]=igroup;
00670     
00671     while(nleft>0){
00672       //MSG("PreFilter",Msg::kWarning) << "in clustering while loop" << std::endl;
00673       int nlk = 0;
00674       float dmin;
00675       do{
00676         dmin = cClustDistcut+1.;
00677         for(int i=0;i<ndigits;i++){
00678           for(int j=0;j<ndigits;j++){
00679             if(group[i]==igroup && group[j]==-2 && dst[i][j]>0. && 
00680                dst[i][j]<dmin) {
00681               dmin = dst[i][j];
00682               i0 = i;
00683               j0 = j;
00684             }
00685           }
00686         }
00687         if(dmin<cClustDistcut){
00688           group[j0] = igroup;
00689           nlk++;
00690         }
00691       }while(dmin<cClustDistcut);
00692       
00693       if(nlk==0) group[i0] = -1;
00694       
00695       if(nlk>0){
00696         groupsize[igroup] = nlk+1;
00697         igroup++;
00698       }
00699       
00700       // calculate nleft
00701       nleft = 0;
00702       for(unsigned int i=0;i<group.size();i++){
00703         if(group[i]==-2) nleft++;
00704       }
00705       if(nleft>0){
00706         int i=0;
00707         while(group[i]!=-2) i++;
00708         i0 = i;
00709         group[i0] = igroup;
00710       }
00711     }
00712     
00713     std::vector<std::vector<CandDigitHandle*> >retclusters(groupsize.size());
00714     
00715     for(std::map<int,int>::iterator gsiter = groupsize.begin();
00716         gsiter!=groupsize.end();gsiter++){
00717       
00718       for(unsigned int i=0;i<group.size();i++){
00719         if((gsiter->first)==group[i]) 
00720           retclusters[gsiter->first].push_back(digits[i]);
00721       }
00722     }
00723     return retclusters;
00724 }
00725 //.............................................................................
00726 std::vector<std::pair<int,int > > 
00727 PreFilter::ClusterMatch(std::vector<std::vector<CandDigitHandle*> > &Uclusters,
00728                         std::vector<std::vector<CandDigitHandle*> > &Vclusters)
00729 {
00730   // want to find the clusters that overlap in z to make meaningful containment
00731   // estimate. this method returns a vector holding pairs of ints. the first
00732   // int in each pair represents a particular U cluster, and the second int in
00733   // the pair represents a cluster in V that overlaps with that U cluster.
00734 
00735     std::vector<std::pair<float,float> >Uzminmax(Uclusters.size());
00736     std::vector<std::pair<float,float> >Vzminmax(Vclusters.size());
00737     
00738     for(unsigned int i=0;i<Uclusters.size();i++){
00739       float zmin = 35.;
00740       float zmax = 0.;
00741       
00742       for(unsigned int j=0;j<Uclusters[i].size();j++){
00743         float z = this->GetCDHZPos(Uclusters[i][j]);
00744         
00745         //MSG("PreFilter",Msg::kWarning) << "U digit z is " << z << std::endl;
00746         
00747         if(z<zmin) zmin = z;
00748         if(z>zmax) zmax = z;
00749       }
00750       Uzminmax[i] = std::make_pair(zmin,zmax);
00751     }
00752     
00753     for(unsigned int i=0;i<Vclusters.size();i++){
00754       float zmin = 35.;
00755       float zmax = 0.;
00756       
00757       for(unsigned int j=0;j<Vclusters[i].size();j++){
00758         float z = this->GetCDHZPos(Vclusters[i][j]);
00759         
00760         //MSG("PreFilter",Msg::kWarning) << "V digit z is " << z << std::endl;
00761         
00762         if(z<zmin) zmin = z;
00763         if(z>zmax) zmax = z;
00764       }
00765       Vzminmax[i] = std::make_pair(zmin,zmax);
00766     }
00767     
00768     std::vector<std::pair<int,int> >OLclusters;
00769     
00770     for(unsigned int i=0;i<Uzminmax.size();i++){
00771       float uzmin = Uzminmax[i].first;
00772       float uzmax = Uzminmax[i].second;
00773       
00774       for(unsigned int j=0;j<Vzminmax.size();j++){
00775         float vzmin = Vzminmax[j].first;
00776         float vzmax = Vzminmax[j].second;
00777         
00778         if( ((uzmin > vzmin) && (uzmin < vzmax)) || 
00779           ((uzmax > vzmin) && (uzmax < vzmax)) )
00780           OLclusters.push_back(std::make_pair(i,j));
00781       }
00782     }
00783     return OLclusters;
00784 }
00785 //.............................................................................
00786 float PreFilter::GetCDHTPos(CandDigitHandle *cdh ){
00787     const PlexSEIdAltL &pseidalt = dynamic_cast<const PlexSEIdAltL&>
00788       (cdh->GetPlexSEIdAltL());
00789     const VldContext *vld = dynamic_cast<const VldContext*> 
00790       (cdh->GetVldContext());
00791     UgliGeomHandle ugh(*vld);
00792     UgliStripHandle ush = ugh.GetStripHandle(pseidalt.GetBestSEId());
00793     return (ush.GetTPos());
00794 }
00795 //.............................................................................
00796 float PreFilter::GetCDHZPos(CandDigitHandle *cdh){
00797     const PlexSEIdAltL &pseidalt = dynamic_cast<const PlexSEIdAltL&>
00798       (cdh->GetPlexSEIdAltL());
00799     const VldContext *vld = dynamic_cast<const VldContext*> 
00800       (cdh->GetVldContext());
00801     UgliGeomHandle ugh(*vld);
00802     UgliStripHandle ush = ugh.GetStripHandle(pseidalt.GetBestSEId());
00803     return (ush.GlobalPos(0).Z());
00804 }
00805 //.............................................................................
00806 void PreFilter::BeginSnarl(){
00807     for(int i=0;i<50;i++){
00808         fUClustSize[i] = 0;
00809         fVClustSize[i] = 0;
00810         fUClustChrg[i] = 0.;
00811         fVClustChrg[i] = 0.;
00812         fUClustMinZ[i] = 100.;
00813         fVClustMinZ[i] = 100.;
00814         fUClustMaxZ[i] = -1.;
00815         fVClustMaxZ[i] = -1.;
00816         fUSlope[i] = 0.;
00817         fUIntercept[i] = 0.;
00818         fVSlope[i] = 0.;
00819         fVIntercept[i] = 0.;
00820     }
00821 
00822     if(cDiagnosticMode){
00823       fSnrlChrgDecD = -1;
00824       fContainDecD = -1;
00825       fNPlanesDecD = -1;
00826       fSnrlChrgDecT = -1;
00827       fContainDecT = -1;
00828       fNPlanesDecT = -1;
00829       fProblems = 0;
00830       
00831       fRun = -1;
00832       fSubrun = -1;
00833       fSnrl = -1;
00834       fNRawDig = -1;
00835       fNPlanes = -1;
00836       fNTruPlanes = 0;
00837       fUdigSize = -1;
00838       fVdigSize = -1;
00839       fNUClusters = 0;
00840       fNVClusters = 0;
00841       fNOLClusters = -1;
00842       fTotUClustDig = 0;
00843       fTotVClustDig = 0;
00844       fSnrlChrg = 0.;
00845       fUDigChrg = 0.;
00846       fVDigChrg = 0.;
00847       fTotUClustChrg = 0.;
00848       fTotVClustChrg = 0.;
00849       fTotTruE = 0.;
00850 
00851       for(int i=0;i<4;i++){
00852         fP4Neu[i] = 0.;
00853       }
00854     }    
00855 }
00856 //.............................................................................
00857 // temporary method to fill variables written to tree to evaluate how filter
00858 // is doing
00859 void PreFilter::Diagnostic(std::vector<std::vector<CandDigitHandle*> >
00860                            &Uclusters,
00861                            std::vector<std::vector<CandDigitHandle*> >
00862                            &Vclusters){
00863   
00864     // fill clustsize, clustchrg, clustminZ, clustmaxZ
00865     for(unsigned int i=0;i<Uclusters.size();i++){
00866       fUClustSize[i] = Uclusters[i].size();
00867       fTotUClustDig += Uclusters[i].size();
00868       
00869       for(std::vector<CandDigitHandle*>::iterator 
00870             uclustiter=Uclusters[i].begin();uclustiter!=Uclusters[i].end();
00871           uclustiter++){
00872         
00873         fUClustChrg[i] += (*uclustiter)->GetCharge(CalDigitType::kPE);
00874         
00875         float z = this->GetCDHZPos((*uclustiter));
00876         if(z>fUClustMaxZ[i]) fUClustMaxZ[i] = z;
00877         if(z<fUClustMinZ[i]) fUClustMinZ[i] = z;
00878       }
00879       fTotUClustChrg += fUClustChrg[i];
00880     }
00881     for(unsigned int i=0;i<Vclusters.size();i++){
00882       fVClustSize[i] = Vclusters[i].size();
00883       fTotVClustDig += Vclusters[i].size();
00884       
00885       for(std::vector<CandDigitHandle*>::iterator 
00886             vclustiter=Vclusters[i].begin();vclustiter!=Vclusters[i].end();
00887           vclustiter++){
00888         
00889         fVClustChrg[i] += (*vclustiter)->GetCharge(CalDigitType::kPE);
00890         
00891         float z = this->GetCDHZPos((*vclustiter));
00892         if(z>fVClustMaxZ[i]) fVClustMaxZ[i] = z;
00893         if(z<fVClustMinZ[i]) fVClustMinZ[i] = z;
00894       }
00895       fTotVClustChrg += fVClustChrg[i];
00896     }
00897 }
00898 //............................................................................
00899 void PreFilter::SimCheck(const MomNavigator *mom)
00900 {
00901     const RawRecord *rawrec = dynamic_cast<const RawRecord*>
00902       (mom->GetFragment("RawRecord"));
00903     const RawDaqSnarlHeader* rdsh = dynamic_cast<const RawDaqSnarlHeader*>
00904       (rawrec->GetRawHeader());
00905     const VldContext &vld = rdsh->GetVldContext();
00906     
00907     const SimSnarlRecord* simrec = dynamic_cast<const SimSnarlRecord*>
00908       (mom->GetFragment("SimSnarlRecord"));
00909     
00910     if(simrec){
00911       // true kinematic info
00912       REROOT_NeuKin* neukin = NULL;
00913       int ctor = 0;
00914       TObjArray objarr = simrec->GetComponents();
00915       
00916       for(int i = 0;i <= objarr.GetLast();i++){
00917         neukin = dynamic_cast<REROOT_NeuKin*>(objarr.At(i));
00918         if(neukin){
00919           ctor = i;
00920         }
00921       }
00922       neukin = dynamic_cast<REROOT_NeuKin*>(objarr.At(ctor));
00923       
00924       const Float_t *tmp = neukin->P4Neu();
00925       fP4Neu[0] = tmp[0];
00926       fP4Neu[1] = tmp[1];
00927       fP4Neu[2] = tmp[2];
00928       fP4Neu[3] = tmp[3];
00929       
00930       // containment stuff
00931       const TClonesArray *digiscinthits = dynamic_cast<const TClonesArray*>
00932         (simrec->FindComponent("TClonesArray","DigiScintHits"));
00933       
00934       std::multimap<int,DigiScintHit*> plnhits;
00935       std::map<int,std::multimap<int,DigiScintHit*> > plnstriphits;
00936       
00937       for(int i=0;i<digiscinthits->GetEntries();i++){
00938         DigiScintHit *hit = dynamic_cast<DigiScintHit*>(digiscinthits->At(i));
00939         plnhits.insert(make_pair(hit->Plane(),hit));
00940       }
00941       
00942       for(int i=0;i<487;i++){
00943         std::pair<std::multimap<int,DigiScintHit*>::const_iterator,
00944           std::multimap<int,DigiScintHit*>::const_iterator> plnhititerpair = 
00945           plnhits.equal_range(i);
00946         
00947         if(plnhititerpair.first!=plnhititerpair.second){
00948           std::multimap<int,DigiScintHit*> striphits;
00949           for(std::multimap<int,DigiScintHit*>::const_iterator plnhititer = 
00950                 plnhititerpair.first; plnhititer!=plnhititerpair.second;
00951               plnhititer++){
00952             striphits.insert(std::make_pair((plnhititer->second)->Strip(),
00953                                             plnhititer->second));
00954           }
00955           plnstriphits.insert(std::make_pair(i,striphits));
00956         }
00957       }
00958       
00959       UgliGeomHandle ugh(vld,Ugli::kUseGlobal);
00960       bool contfail = false;
00961       
00962       typedef std::map<int,std::multimap<int,DigiScintHit*> >::const_iterator
00963         MAPITER;
00964       typedef std::multimap<int,DigiScintHit*>::const_iterator MMAPITER;
00965       
00966       for(int i=0;i<487;i++){
00967         bool realplane = false;
00968         std::pair<MAPITER,MAPITER> plnstriphititerpair = 
00969           plnstriphits.equal_range(i);
00970         if(plnstriphititerpair.first!=plnstriphititerpair.second){
00971           
00972           for(int j=0;j<192;j++){
00973             std::pair<MMAPITER,MMAPITER> striphititerpair = 
00974               ((plnstriphititerpair.first)->second).equal_range(j);
00975             if(striphititerpair.first!=striphititerpair.second){
00976           
00977               float stripE = 0.;
00978               for(MMAPITER hititer = striphititerpair.first; 
00979                   hititer!= striphititerpair.second; hititer++){
00980                 stripE += (hititer->second)->DE();
00981               }
00982               fTotTruE += stripE;
00983               if(stripE>0.0015){
00984                 realplane = true;
00985                 
00986                 DigiScintHit *hit = (striphititerpair.first)->second;
00987                 PlexStripEndId pseid = hit->StripEndId();
00988                 const UgliStripHandle ush = ugh.GetStripHandle(pseid);
00989                 TVector3 localvec(hit->X1(),hit->Y1(),hit->Z1());
00990                 TVector3 globalvec(ush.LocalToGlobal(localvec));
00991                 float x = globalvec.X();
00992                 float y = globalvec.Y();
00993                 float z = globalvec.Z();
00994                 float u = 0.707107*(x+y);
00995                 float v = 0.707107*(y-x);
00996 
00997                 if( ((x*x) > cUVdistcut*cUVdistcut) || 
00998                     ((y*y) > cUVdistcut*cUVdistcut) || 
00999                     ((u*u) > cUVdistcut*cUVdistcut) || 
01000                     ((v*v) > cUVdistcut*cUVdistcut) || (z>(fSnrlZMax-0.13)) || 
01001                     (z<(fSnrlZMin+0.13)) ){
01002                   contfail = true;
01003                   //j = 192;
01004                   //i = 487;
01005                 }
01006               }
01007             }
01008           }
01009         }
01010         if(realplane) fNTruPlanes++; 
01011       }
01012       // make cut decisions here. use gross estimate of 
01013       // 1.5 Gev Deposited : 5000 PE to compare true energy/snarl charge      
01014       fSnrlChrgDecT = ((fTotTruE*5000/1.5)<cSnrlChrgCut)?0:1;
01015       fContainDecT = (contfail)?0:1;
01016       fNPlanesDecT = (fNTruPlanes<cMinNPlanes)?0:1;
01017     }
01018 }