NOvA Offline: clusterss::CompareClusters Class Reference

Constructor & Destructor Documentation

CompareClusters::CompareClusters

(

const char *

version

)

Definition at line 38 of file CompareClusters.cxx.

References clustdiff, matchvsize, maxcellmipnomatch, maxmipvsize, mipnototmatch, mucatchproblem, ncellnototmatch, ncells, nclust, nomatchncellvmip, perfectmatchcluster, perfectmatchevent, and jobc::Module::SetCfgVersion().

00038                                                     :
00039   jobc::Module("makeclusterss::CompareClusters"),
00040   nevt(0)
00041 {
00042   nclust = new TH2D("nclust",";Found Clusters Old; Found Clusters New",21,-0.5,20.5,21,-0.5,20.5);
00043   clustdiff = new TH1D("clustdiff",";Diff in Found Clusters",41,-20.5,20.5);
00044 
00045   ncells = new TH2D("ncells",";Found Cells Old; Found Cells New",21,-0.5,20.5,21,-0.5,20.5);
00046 
00047   matchvsize = new TH2D("matchvsize",";Cluster Size; N Matched",21,-0.5,20.5,21,-0.5,20.5);
00048   nomatchncellvmip = new TH2D("nomatchncellvmip",
00049                               ";NCells in Unmatched Cluster; MIP in Unmatched",21,-0.5,20.5,500,0,10);
00050   ncellnototmatch = new TH1D("ncellnototmatch",";NCells in Unmatched Cluster",21,-0.5,20.5);
00051   mipnototmatch = new TH1D("mipnototmatch",";MIP in Unmatched Cluster",500,0,10);
00052   maxcellmipnomatch = new TH1D("maxcellmipnomatch",";Max MIP of unmatched cell",500,0,10);
00053   maxmipvsize = new TH2D("maxmipvsize",";NCells; Max MIP unmatched cell",21,-0.5,20.5,500,0,10);
00054   mucatchproblem = new TH1D("mucatchproblem","",2,-0.5,1.5);
00055   perfectmatchcluster = new TH1D("perfectmatchcluster","",2,-0.5,1.5);
00056   perfectmatchevent = new TH1D("perfectmatchevent","",2,-0.5,1.5);
00057   // Set the configuration version tag
00058   this->SetCfgVersion(version); // Required!
00059 }

CompareClusters::~CompareClusters

(

)

Definition at line 69 of file CompareClusters.cxx.

00070 {
00071 }

Member Function Documentation

jobc::Result CompareClusters::Ana ( const edm::EventHandle & evt ) [virtual]

Reimplemented from jobc::Module.

Definition at line 83 of file CompareClusters.cxx.

References edm::EventHandle::Cal(), clustdiff, edm::EventHandle::DAQ(), geo::Geometry::DetLength(), edm::EventHeader::Event(), recobase::ShowerBase::GetNCell(), edm::EventHandle::Header(), geo::Geometry::Instance(), MatchCells(), matchvsize, maxcellmipnomatch, maxmipvsize, edm::EventHandle::MC(), mipnototmatch, mucatchproblem, recobase::Cluster::NCell(), ncellnototmatch, ncells, nclust, nevt, nomatchncellvmip, perfectmatchcluster, perfectmatchevent, recobase::Cluster::Q(), edm::EventHandle::Reco(), and edm::EventHeader::Run().

00084 {
00085   const edm::EventHeader &eh = evt.Header();
00086   int run = eh.Run();
00087   int event = eh.Event();
00088   nevt++;
00089 
00090   //get the list of CellHits
00091   std::vector<const recobase::CellHit*> hitlist(0);
00092   try{ evt.Cal().Get("./MergeHits",hitlist); }
00093   catch(edm::Exception e){
00094     cerr<<"Couldn't find the MergedCellHit list, trying for the original"<<endl;
00095     try{ evt.Cal().Get("./Hits",hitlist); }
00096     catch(edm::Exception e){
00097       cerr<<"Couldn't find the Original CellHit list either, failing."<<endl;
00098       return jobc::kFailed;
00099     }
00100   }
00101   if(hitlist.size()<5) return jobc::kFailed; //lets just look at things that have some real activity
00102   
00103   //figure out if this crosses into the muon catcher
00104   int begplane=-1;
00105   int endplane=-1;
00106   for(unsigned int i=0;i<hitlist.size();i++){
00107     int plane = hitlist[i]->Plane();
00108     if(plane<begplane||begplane==-1){
00109       begplane=plane;
00110     }
00111     if(plane>endplane){
00112       endplane=plane;
00113     }
00114   }
00115   if(begplane<=186&&endplane>=186){
00116     //    cout<<"Event crosses mu catcher, clustering won't match"<<endl;
00117     mucatchproblem->Fill(1);
00118     return jobc::kFailed;
00119   }
00120   
00121 
00122   std::vector<const sim::MCTruth*> truelist(0);
00123   try { evt.MC().Get(".", truelist); }
00124   catch(edm::Exception e){
00125     //    cerr<<"Couldn't get a truelist"<<endl;
00126     return jobc::kFailed;
00127   }
00128 
00129   //now get a header
00130   std::vector<const hdr::Header*> head(0);
00131   try{ evt.DAQ().Get(".",head); }
00132   catch(edm::Exception e){
00133     cerr<<"Couldn't get a head"<<endl;
00134     return jobc::kFailed;
00135   }
00136   geo::Geometry &geom = geo::Geometry::Instance(1,head[0]->DetectorID());
00137 
00138   int tidx=0;
00139   //get some truthiness
00140   const TParticle &nupart = truelist[tidx]->GetNeutrino()->Nu();
00141   double eenu = 1.*nupart.Energy();
00142   //  if(eenu<0.100) return jobc::kFailed;//let's just look at reasonable energy nus
00143 
00144   double nuvtxx=nupart.Vx();
00145   double nuvtxy=nupart.Vy();
00146   double nuvtxz=nupart.Vz();
00147 
00148   int vcont=1;
00149   
00150   //  if(nuvtxx>-1.*(geom.DetHalfWidth()-20)&&nuvtxx<1.*(geom.DetHalfWidth()-20)){//20 cm in from xedges
00151     
00152   //      if(nuvtxy>-1.*(geom.DetHalfHeight()-20)&&nuvtxy<1.*(geom.DetHalfHeight()-20)){//20 cm in from yedges
00153   //      if(nuvtxz>20&&nuvtxz<(geom.DetLength()-300)){//20 cm from front, and 300 cm from back (muon catcher)
00154   //      vcont=1;
00155   //}
00156   //}
00157   //    }
00158   if(nuvtxz>(geom.DetLength()-300)){// 300 cm from back (muon catcher)
00159     vcont=0;
00160   }
00161     
00162 
00163 
00164   if(vcont==0){ 
00165     return jobc::kFailed;
00166   }
00167 
00168 
00169   std::vector<const recobase::Cluster *>xclusterlist(0);
00170   bool foundx=true;
00171   try{evt.Reco().Get("./SSClusterX",xclusterlist);}
00172   catch(edm::Exception e){
00173     //    cerr<<"Couldn't find an xcluster"<<endl;
00174     foundx=false;
00175   }
00176 
00177   std::vector<const recobase::Cluster *>yclusterlist(0);
00178   bool foundy=true;
00179   try{evt.Reco().Get("./SSClusterY",yclusterlist);}
00180   catch(edm::Exception e){
00181     //    cerr<<"Couldn't find an ycluster"<<endl;
00182     foundy=false;
00183   }
00184 
00185   std::vector<const recobase::ShowerBase *>showerlist(0);
00186   bool foundshowers=true;
00187   try{evt.Reco().Get("./Final2DShw",showerlist);}
00188   catch(edm::Exception e){
00189     //    cerr<<"Couldn't find any showers"<<endl;
00190     foundshowers=false;
00191   }
00192 
00193   //  cout<<"I have the clusters and showers: "<<xclusterlist.size()<<" "<<yclusterlist.size()<<" "<<showerlist.size()<<endl;
00194 
00195   nclust->Fill(showerlist.size(),xclusterlist.size()+yclusterlist.size());
00196   clustdiff->Fill(1.*showerlist.size()-1.*(xclusterlist.size()+yclusterlist.size()));
00197   if(showerlist.size()>xclusterlist.size()+yclusterlist.size()){
00198     cout<<"NS found extra clusters: Run: "<<run<<" event "<<event<<" iteration "<<nevt<<" size "<<showerlist.size()<<" compared to "<<xclusterlist.size()+yclusterlist.size()<<endl;
00199   }
00200   std::vector<const recobase::Cluster *> &clist=xclusterlist;
00201 
00202   for(int i=0;i<2;i++){
00203     bool fc=false;
00204     geo::View_t view;
00205     if(i==0){
00206       clist=xclusterlist;
00207       fc = foundx;
00208       view=geo::kX;
00209     }
00210     else{
00211       clist=yclusterlist;
00212       fc = foundy;
00213       view=geo::kY;
00214     }
00215              
00216     if(fc){
00217       bool pmev=true;
00218       for(unsigned int i=0;i<clist.size();i++){
00219         //      cout<<"NEw tv cell"<<endl;
00220         int mostmatchedcells=-2;
00221         int ncellsbest=0;
00222         float bestunmatched=0.;
00223         recobase::Cluster *c = new recobase::Cluster(*clist[i]);
00224         //      cout<<"is Qfilled? "<<c->Q()<<endl;
00225         if(foundshowers){
00226           for(unsigned int j=0;j<showerlist.size();j++){
00227             recobase::ShowerBase *sb = new recobase::ShowerBase(*showerlist[j]);
00228             float unmatchedmip;
00229             int match=MatchCells(view,c,sb,unmatchedmip);
00230             //      cout<<"matches "<<match<<" most matches "<<mostmatchedcells<<endl;
00231             if(match>mostmatchedcells){
00232               mostmatchedcells=match;
00233               ncellsbest = sb->GetNCell();
00234               bestunmatched = unmatchedmip;
00235             }
00236           }
00237         }
00238         //      cout<<"After all that "<<mostmatchedcells<<endl;
00239 
00240         ncells->Fill(ncellsbest,c->NCell(view));
00241         matchvsize->Fill(c->NCell(view),mostmatchedcells);
00242         maxcellmipnomatch->Fill(bestunmatched);
00243         maxmipvsize->Fill(c->NCell(view),bestunmatched);
00244         if((c->NCell(view)-mostmatchedcells)>2&&mostmatchedcells!=0&&bestunmatched>2){
00245           cout<<"MISMATCH!  Run: "<<run<<" event "<<event<<" iteration "<<nevt<<" this cluster has "<<c->NCell(view)
00246               <<" matched "<<mostmatchedcells<<" maxmip unmatched "<<bestunmatched<<endl;
00247         }
00248         if(mostmatchedcells==0){
00249           ncellnototmatch->Fill(c->NCell(view));
00250           //      cout<<"mipnototmatch "<<clist[i]->Q()<<endl;
00251           mipnototmatch->Fill(c->Q());
00252           nomatchncellvmip->Fill(c->NCell(view),c->Q());
00253           if(c->NCell(view)>2&&c->Q()>2){
00254             cout<<"You found an extra cluster!  Run: "<<run<<" event "<<event<<" iteration "<<nevt<<endl;
00255           }
00256         }
00257         if(mostmatchedcells==c->NCell(view)) perfectmatchcluster->Fill(1);
00258         else{
00259           pmev=false;
00260           perfectmatchcluster->Fill(0);
00261         }
00262         
00263       }
00264       if(pmev) perfectmatchevent->Fill(1);
00265       else perfectmatchevent->Fill(0);
00266     }
00267   }
00268 
00269   return jobc::kPassed; // kFailed if you want to fail the event
00270   
00271 }

int CompareClusters::MatchCells	(	geo::View_t	view,
		recobase::Cluster *	c,
		recobase::ShowerBase *	s,
		float &	unmatchedmip
	)

Definition at line 273 of file CompareClusters.cxx.

References rawdata::RawDigit::ADC(), recobase::CellHit::Cell(), recobase::Cluster::Cell(), recobase::ShowerBase::GetCells(), recobase::Cluster::NCell(), and recobase::CellHit::Plane().

Referenced by Ana().

00274 {
00275   int matches=0;
00276   unmatched=0.;
00277 
00278   //  if(s->GetView()!=view) return matches;//wrongview no matches.
00279 
00280   int NCELLCLUSTER=c->NCell(view);
00281 
00282   const std::map<Int_t,std::multimap<Double_t,recobase::CellHit*> >&cellmap = s->GetCells();      
00283 
00284   int lowplane1=-1;
00285   int lowplane2=-1;
00286   int hiplane1=-1;
00287   int hiplane2=-1;
00288 
00289   for(int j=0;j<NCELLCLUSTER;j++){  
00290     bool foundamatch=false;
00291     const recobase::CellHit *thisclustercell = c->Cell(view,j);
00292     int thisplane = thisclustercell->Plane();
00293     if(j==0){
00294       hiplane1=thisplane;
00295       lowplane1=thisplane;
00296     }
00297     if(thisplane>hiplane1){
00298       hiplane1=thisplane;
00299     }
00300     if(thisplane<lowplane1){
00301       lowplane1=thisplane;
00302     }
00303     int thiscellnumber = thisclustercell->Cell();
00304     std::map<Int_t,std::multimap<Double_t,recobase::CellHit*> >::const_iterator mit = cellmap.begin();      
00305     while(mit!=cellmap.end()){
00306       if(mit->first>hiplane2) hiplane2=mit->first;
00307       if(mit->first<lowplane2||lowplane2==-1) lowplane2=mit->first;      
00308       if(mit->first!=thisplane){
00309         mit++;
00310         continue;
00311       }
00312       std::multimap<Double_t,recobase::CellHit*> cellmultimap = mit->second;
00313       std::multimap<Double_t,recobase::CellHit*>::iterator cmit2=cellmultimap.begin();           
00314       while(cmit2!=cellmultimap.end()){
00315         if((*cmit2->second).Plane()==thisplane&&
00316            (*cmit2->second).Cell()==thiscellnumber){
00317           matches++;
00318           foundamatch=true;
00319           break;
00320         }
00321         cmit2++;
00322       }
00323       if(foundamatch=true){
00324         break;
00325       }
00326       mit++;
00327     }
00328     if(!foundamatch){
00329       float mip;
00330       bool calibrated = thisclustercell->MIP(mip);
00331       if(!calibrated){
00332         unsigned short adc;
00333         thisclustercell->ADC(adc);
00334         mip = adc/1.;
00335       }
00336       if(mip>unmatched){
00337         unmatched=mip;
00338       }
00339     }
00340   }
00341   //  cout<<"TV low: "<<lowplane1<<" hi: "<<hiplane1<<endl;
00342   //  cout<<"NS low: "<<lowplane2<<" hi: "<<hiplane2<<endl;
00343   if(lowplane1<=186&&hiplane1>=186){
00344     //    cout<<"Muon catcher straddler TV"<<endl;
00345     return -1;
00346   }
00347   if(lowplane2<=186&&hiplane2>=186){
00348     //    cout<<"Muon catcher straddler NS"<<endl;
00349     return -1;
00350   }
00351   return matches;
00352 }

jobc::Result CompareClusters::Reco ( edm::EventHandle & evt ) [virtual]

Reimplemented from jobc::Module.

Definition at line 75 of file CompareClusters.cxx.

00076 {
00077   
00078   return jobc::kPassed; // kFailed if you want to fail the event
00079 }

void CompareClusters::Update ( const cfg::Config & c ) [virtual]

Implements cfg::Observer.

Definition at line 63 of file CompareClusters.cxx.

00064 {
00065 }

Member Data Documentation

Generated on Mon Nov 23 04:45:30 2009 for NOvA Offline by

1.3.9.1


Public Member Functions
	CompareClusters (const char *version)
void	Update (const cfg::Config &c)
	~CompareClusters ()
jobc::Result	Reco (edm::EventHandle &evt)
jobc::Result	Ana (const edm::EventHandle &evt)
int	MatchCells (geo::View_t view, recobase::Cluster c, recobase::ShowerBase s, float &unmatchedmip)
Public Attributes
TH2D *	nclust
TH1D *	clustdiff
TH2D *	ncells
TH2D *	matchvsize
TH2D *	nomatchncellvmip
TH1D *	ncellnototmatch
TH1D *	mipnototmatch
TH1D *	maxcellmipnomatch
TH2D *	maxmipvsize
TH1D *	mucatchproblem
TH1D *	perfectmatchcluster
TH1D *	perfectmatchevent
int	nevt

clusterss::CompareClusters Class Reference

Public Member Functions

Public Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation