WCSimDataCleaner.cc

Go to the documentation of this file.
00001 #include "WCSimDataCleaner.hh"
00002 
00003 #include "WCSimRecoDigit.hh"
00004 #include "WCSimRecoCluster.hh"
00005 #include "WCSimRecoClusterDigit.hh"
00006 
00007 #include <cmath>
00008 #include <iostream>
00009 #include <cassert>
00010 
00011 ClassImp(WCSimDataCleaner)
00012 
00013 static WCSimDataCleaner* fgDataCleaner = 0;
00014 
00015 WCSimDataCleaner* WCSimDataCleaner::Instance()
00016 {
00017   if( !fgDataCleaner ){
00018     fgDataCleaner = new WCSimDataCleaner();
00019   }
00020 
00021   if( !fgDataCleaner ){
00022     assert(fgDataCleaner);
00023   }
00024 
00025   if( fgDataCleaner ){
00026 
00027   }
00028 
00029   return fgDataCleaner;
00030 }
00031   
00032 void WCSimDataCleaner::Config(Int_t config)
00033 {
00034   WCSimDataCleaner::Instance()->SetConfig(config);
00035 }
00036 
00037 void WCSimDataCleaner::MinPulseHeight(Double_t min)
00038 {
00039   WCSimDataCleaner::Instance()->SetMinPulseHeight(min);
00040 }
00041   
00042 void WCSimDataCleaner::NeighbourRadius(Double_t radius)
00043 {
00044   WCSimDataCleaner::Instance()->SetNeighbourRadius(radius);
00045 }
00046   
00047 void WCSimDataCleaner::NeighbourDigits(Int_t digits)
00048 {
00049   WCSimDataCleaner::Instance()->SetNeighbourDigits(digits);
00050 }
00051   
00052 void WCSimDataCleaner::ClusterRadius(Double_t radius)
00053 {
00054   WCSimDataCleaner::Instance()->SetClusterRadius(radius);
00055 }
00056   
00057 void WCSimDataCleaner::ClusterDigits(Int_t digits)
00058 {
00059   WCSimDataCleaner::Instance()->SetClusterDigits(digits);
00060 }
00061   
00062 void WCSimDataCleaner::TimeWindow(Double_t window)
00063 {
00064   WCSimDataCleaner::Instance()->SetTimeWindow(window);
00065 }
00066 
00067 void WCSimDataCleaner::PrintParameters()
00068 {
00069   WCSimDataCleaner::Instance()->RunPrintParameters();
00070 }
00071 
00072 WCSimDataCleaner::WCSimDataCleaner()
00073 {
00074   // cleaning mode
00075   fConfig = WCSimDataCleaner::kPulseHeightAndClusters;
00076 
00077   // default cleaning parameters
00078   fMinPulseHeight = 1.0;     // minimum pulse height (PEs)
00079   fNeighbourRadius = 200.0;  // clustering window (cm)
00080   fMinNeighbourDigits = 2;   // minimum neighbouring digits
00081   fClusterRadius = 200.0;    // clustering window (cm)
00082   fMinClusterDigits = 50;    // minimum clustered digits
00083   fTimeWindow = 25.0;        // timing window (ns)
00084 
00085   // vector of filtered digits
00086   fFilterAll = new std::vector<WCSimRecoDigit*>;
00087   fFilterByPulseHeight = new std::vector<WCSimRecoDigit*>;
00088   fFilterByNeighbours = new std::vector<WCSimRecoDigit*>;
00089   fFilterByClusters = new std::vector<WCSimRecoDigit*>;
00090   
00091   // vector of clusters
00092   fClusterList = new std::vector<WCSimRecoCluster*>;
00093 }
00094 
00095 WCSimDataCleaner::~WCSimDataCleaner()
00096 {
00097   delete fFilterAll;
00098   delete fFilterByPulseHeight;
00099   delete fFilterByNeighbours;
00100   delete fFilterByClusters;
00101   delete fClusterList;
00102 }
00103 
00104 void WCSimDataCleaner::RunPrintParameters()
00105 {
00106   std::cout << " *** WCSimDataCleaner::PrintParameters() *** " << std::endl;
00107   
00108   std::cout << "  Data Cleaner Parameters: " << std::endl
00109             << "   Config = " << fConfig<< std::endl
00110             << "   MinPulseHeight = " << fMinPulseHeight << std::endl
00111             << "   NeighbourRadius = " << fNeighbourRadius << std::endl
00112             << "   MinNeighbourDigits = " << fMinNeighbourDigits << std::endl
00113             << "   ClusterRadius = " << fClusterRadius << std::endl
00114             << "   MinClusterDigits = " << fMinClusterDigits << std::endl
00115             << "   TimeWindow = " << fTimeWindow << std::endl;
00116 
00117   return;
00118 }
00119 
00120 void WCSimDataCleaner::Reset()
00121 {
00122 
00123   return;
00124 }
00125 
00126 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::Run(std::vector<WCSimRecoDigit*>* myDigitList)
00127 {
00128   std::cout << " *** WCSimDataCleaner::Run(...) *** " << std::endl;
00129 
00130   // input digit list
00131   // ================
00132   std::vector<WCSimRecoDigit*>* myInputList = myDigitList;
00133   std::vector<WCSimRecoDigit*>* myOutputList = myDigitList;
00134 
00135   // filter all digits
00136   // =================
00137   myInputList = ResetDigits(myOutputList);
00138   myOutputList = (std::vector<WCSimRecoDigit*>*)(this->FilterAll(myInputList));
00139   myOutputList = FilterDigits(myOutputList);
00140   if( fConfig==WCSimDataCleaner::kNone ) return myOutputList;
00141 
00142   // filter by pulse height
00143   // ======================
00144   myInputList = ResetDigits(myOutputList);
00145   myOutputList = (std::vector<WCSimRecoDigit*>*)(this->FilterByPulseHeight(myInputList));
00146   myOutputList = FilterDigits(myOutputList);
00147   if( fConfig==WCSimDataCleaner::kPulseHeight ) return myOutputList;
00148 
00149   // filter using neighbouring digits
00150   // ================================
00151   myInputList = ResetDigits(myOutputList);
00152   myOutputList = (std::vector<WCSimRecoDigit*>*)(this->FilterByNeighbours(myInputList));
00153   myOutputList = FilterDigits(myOutputList);
00154   if( fConfig==WCSimDataCleaner::kPulseHeightAndNeighbours ) return myOutputList;
00155 
00156   // filter using clustered digits
00157   // =============================
00158   myInputList = ResetDigits(myOutputList);
00159   myOutputList = (std::vector<WCSimRecoDigit*>*)(this->FilterByClusters(myInputList));
00160   myOutputList = FilterDigits(myOutputList);
00161   if( fConfig==WCSimDataCleaner::kPulseHeightAndClusters ) return myOutputList;
00162 
00163   // return vector of filtered digits
00164   // ================================
00165   return myOutputList;
00166 }
00167 
00168 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::ResetDigits(std::vector<WCSimRecoDigit*>* myDigitList)
00169 {
00170   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00171     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00172     recoDigit->ResetFilter();
00173   }
00174 
00175   return myDigitList;
00176 }
00177 
00178 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::FilterDigits(std::vector<WCSimRecoDigit*>* myDigitList)
00179 {
00180   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00181     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00182     recoDigit->PassFilter();
00183   }
00184 
00185   return myDigitList;
00186 }
00187 
00188 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::FilterAll(std::vector<WCSimRecoDigit*>* myDigitList)
00189 {
00190   // clear vector of filtered digits
00191   // ==============================
00192   fFilterAll->clear();
00193 
00194   // filter all digits
00195   // =================
00196   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00197     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00198     fFilterAll->push_back(recoDigit);
00199   }
00200 
00201   // return vector of filtered digits
00202   // ================================
00203   std::cout << "  filter all: " << fFilterAll->size() << std::endl;
00204   
00205   return fFilterAll;
00206 }
00207 
00208 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::FilterByPulseHeight(std::vector<WCSimRecoDigit*>* myDigitList)
00209 {
00210   // clear vector of filtered digits
00211   // ===============================
00212   fFilterByPulseHeight->clear();
00213 
00214   // filter by pulse height
00215   // ======================
00216   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00217     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00218     if( recoDigit->GetQPEs()>fMinPulseHeight ){
00219       fFilterByPulseHeight->push_back(recoDigit);
00220     }
00221   }
00222 
00223   // return vector of filtered digits
00224   // ================================
00225   std::cout << "  filter by pulse height: " << fFilterByPulseHeight->size() << std::endl;
00226   
00227   return fFilterByPulseHeight;
00228 }
00229 
00230 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::FilterByNeighbours(std::vector<WCSimRecoDigit*>* myDigitList)
00231 {
00232   // clear vector of filtered digits
00233   // ===============================
00234   fFilterByNeighbours->clear();
00235 
00236   // create array of neighbours
00237   // ==========================
00238   Int_t Ndigits = myDigitList->size();
00239 
00240   if( Ndigits<=0 ){
00241     return fFilterByNeighbours;
00242   }
00243 
00244   Int_t* numNeighbours = new Int_t[Ndigits];
00245 
00246   for( Int_t idigit=0; idigit<Ndigits; idigit++ ){
00247     numNeighbours[idigit] = 0;
00248   }
00249 
00250   // count number of neighbours
00251   // ==========================
00252   for( UInt_t idigit1=0; idigit1<myDigitList->size(); idigit1++ ){
00253     for( UInt_t idigit2=idigit1+1; idigit2<myDigitList->size(); idigit2++ ){
00254       WCSimRecoDigit* fdigit1 = (WCSimRecoDigit*)(myDigitList->at(idigit1));
00255       WCSimRecoDigit* fdigit2 = (WCSimRecoDigit*)(myDigitList->at(idigit2));
00256 
00257       Double_t dx = fdigit1->GetX() - fdigit2->GetX();
00258       Double_t dy = fdigit1->GetY() - fdigit2->GetY();
00259       Double_t dz = fdigit1->GetZ() - fdigit2->GetZ();
00260       Double_t dt = fdigit1->GetTime() - fdigit2->GetTime();
00261       Double_t drsq = dx*dx + dy*dy + dz*dz;
00262 
00263       if( drsq>0.0
00264        && drsq<fNeighbourRadius*fNeighbourRadius
00265        && fabs(dt)<fTimeWindow ){
00266         numNeighbours[idigit1]++;
00267         numNeighbours[idigit2]++;
00268       }
00269     }
00270   }
00271 
00272   // filter by number of neighbours
00273   // ==============================
00274   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00275     WCSimRecoDigit* fdigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00276     if( numNeighbours[idigit]>=fMinNeighbourDigits ){
00277       fFilterByNeighbours->push_back(fdigit);
00278     }
00279   }
00280 
00281   // delete array of neighbours
00282   // ==========================
00283   delete [] numNeighbours;
00284 
00285   // return vector of filtered digits
00286   // ================================
00287   std::cout << "  filter by neighbours: " << fFilterByNeighbours->size() << std::endl;
00288   
00289   return fFilterByNeighbours;
00290 }
00291 
00292 std::vector<WCSimRecoDigit*>* WCSimDataCleaner::FilterByClusters(std::vector<WCSimRecoDigit*>* myDigitList)
00293 {
00294   // clear vector of filtered digits
00295   // ===============================
00296   fFilterByClusters->clear();
00297 
00298   // run clustering algorithm
00299   // ========================
00300   std::vector<WCSimRecoCluster*>* myClusterList = (std::vector<WCSimRecoCluster*>*)(this->RecoClusters(myDigitList));
00301 
00302   for( UInt_t icluster=0; icluster<myClusterList->size(); icluster++ ){
00303     WCSimRecoCluster* myCluster = (WCSimRecoCluster*)(myClusterList->at(icluster));
00304     for( Int_t idigit=0; idigit<myCluster->GetNDigits(); idigit++ ){
00305       WCSimRecoDigit* myDigit = (WCSimRecoDigit*)(myCluster->GetDigit(idigit));
00306       fFilterByClusters->push_back(myDigit);
00307     }
00308   }
00309   
00310   // return vector of filtered digits
00311   // ================================
00312   std::cout << "  filter by clusters: " << fFilterByClusters->size() << std::endl;
00313   
00314   return fFilterByClusters;
00315 }
00316 
00317 std::vector<WCSimRecoCluster*>* WCSimDataCleaner::RecoClusters(std::vector<WCSimRecoDigit*>* myDigitList)
00318 {  
00319 
00320   // delete cluster digits
00321   // =====================
00322   for( UInt_t i=0; i<vClusterDigitList.size(); i++ ){
00323     delete (WCSimRecoClusterDigit*)(vClusterDigitList.at(i));
00324   }
00325   vClusterDigitList.clear();
00326 
00327   // delete clusters
00328   // ===============
00329   for( UInt_t i=0; i<vClusterList.size(); i++ ){
00330     delete (WCSimRecoCluster*)(vClusterList.at(i));
00331   }
00332   vClusterList.clear();
00333 
00334   // clear vector clusters
00335   // =====================
00336   fClusterList->clear();
00337 
00338   // make cluster digits
00339   // ===================
00340   for( UInt_t idigit=0; idigit<myDigitList->size(); idigit++ ){
00341     WCSimRecoDigit* recoDigit = (WCSimRecoDigit*)(myDigitList->at(idigit));
00342     WCSimRecoClusterDigit* clusterDigit = new WCSimRecoClusterDigit(recoDigit);
00343     vClusterDigitList.push_back(clusterDigit);
00344   }
00345 
00346   // run clustering algorithm
00347   // ========================
00348   for( UInt_t idigit1=0; idigit1<vClusterDigitList.size(); idigit1++ ){
00349     for( UInt_t idigit2=idigit1+1; idigit2<vClusterDigitList.size(); idigit2++ ){
00350 
00351       WCSimRecoClusterDigit* fdigit1 = (WCSimRecoClusterDigit*)(vClusterDigitList.at(idigit1));
00352       WCSimRecoClusterDigit* fdigit2 = (WCSimRecoClusterDigit*)(vClusterDigitList.at(idigit2));
00353 
00354       Double_t dx = fdigit1->GetX() - fdigit2->GetX();
00355       Double_t dy = fdigit1->GetY() - fdigit2->GetY();
00356       Double_t dz = fdigit1->GetZ() - fdigit2->GetZ();
00357       Double_t dt = fdigit1->GetTime() - fdigit2->GetTime();
00358       Double_t drsq = dx*dx + dy*dy + dz*dz;
00359 
00360       if( drsq>0.0
00361        && drsq<fClusterRadius*fClusterRadius
00362        && fabs(dt)<fTimeWindow ){
00363         fdigit1->AddClusterDigit(fdigit2);
00364         fdigit2->AddClusterDigit(fdigit1);
00365       }
00366     }
00367   }
00368 
00369   // collect up clusters
00370   // ===================
00371   Bool_t carryon = 0;
00372 
00373   for( UInt_t idigit=0; idigit<vClusterDigitList.size(); idigit++ ){
00374     WCSimRecoClusterDigit* fdigit = (WCSimRecoClusterDigit*)(vClusterDigitList.at(idigit));
00375 
00376     if( fdigit->IsClustered()==0
00377      && fdigit->GetNClusterDigits()>0 ){
00378         
00379       vClusterDigitCollection.clear();
00380       vClusterDigitCollection.push_back(fdigit);
00381       fdigit->SetClustered();
00382 
00383       carryon = 1;
00384       while( carryon ){
00385         carryon = 0;
00386 
00387         for( UInt_t jdigit=0; jdigit<vClusterDigitCollection.size(); jdigit++ ){
00388           WCSimRecoClusterDigit* cdigit = (WCSimRecoClusterDigit*)(vClusterDigitCollection.at(jdigit));
00389 
00390           if( cdigit->IsAllClustered()==0 ){
00391             for( Int_t kdigit=0; kdigit<cdigit->GetNClusterDigits(); kdigit++ ){
00392               WCSimRecoClusterDigit* cdigitnew = (WCSimRecoClusterDigit*)(cdigit->GetClusterDigit(kdigit));
00393                 
00394               if( cdigitnew->IsClustered()==0 ){
00395                 vClusterDigitCollection.push_back(cdigitnew);
00396                 cdigitnew->SetClustered();
00397                 carryon = 1;
00398               }
00399             }
00400           }
00401         }
00402       }
00403     
00404       if( (Int_t)vClusterDigitCollection.size()>=fMinClusterDigits ){
00405         WCSimRecoCluster* cluster = new WCSimRecoCluster();
00406         fClusterList->push_back(cluster);
00407         vClusterList.push_back(cluster);
00408 
00409         for( UInt_t jdigit=0; jdigit<vClusterDigitCollection.size(); jdigit++ ){
00410           WCSimRecoClusterDigit* cdigit = (WCSimRecoClusterDigit*)(vClusterDigitCollection.at(jdigit));
00411           WCSimRecoDigit* recodigit = (WCSimRecoDigit*)(cdigit->GetRecoDigit());
00412           cluster->AddDigit(recodigit);        
00413         }
00414       }
00415     }
00416   }
00417 
00418   // return vector of clusters
00419   // =========================
00420   return fClusterList;
00421 }
00422 
00423