loon: PreFilter.cxx Source File

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