TruthCompareAna Class Reference

#include <TruthCompareAna.h>

List of all members.

Public Member Functions

 TruthCompareAna ()
 ~TruthCompareAna ()
void ana (ParticleObjectHolder *poh, TruthCompare *t)
void Reset ()

Private Member Functions

std::vector< int > MatchRecoParticle (Particle3D *p, ParticleObjectHolder *poh)
void ProcessTrueParticles (ParticleObjectHolder *poh)
std::vector< int > MatchRecoElectron (Particle3D *p, ParticleObjectHolder *poh)
std::vector< int > MatchRecoMuon (Particle3D *p, ParticleObjectHolder *poh)
std::vector< int > MatchRecoProton (Particle3D *p, ParticleObjectHolder *poh)
void ComputeStatistics (ParticleObjectHolder *poh, TruthCompare *t)

Private Attributes

std::vector< int > true_idx
std::vector< int > matched
double reco_matched_visible_e
double reco_visible_e
double matchangle
double meupergev
double sigcorpermip

Detailed Description

Definition at line 9 of file TruthCompareAna.h.


Constructor & Destructor Documentation

TruthCompareAna::TruthCompareAna (  ) 

Definition at line 10 of file TruthCompareAna.cxx.

References meupergev, and Reset().

00011 {
00012         Reset();
00013         
00014         meupergev = 25.;//should use something more precise here!
00015 }

TruthCompareAna::~TruthCompareAna (  ) 

Definition at line 17 of file TruthCompareAna.cxx.

References Reset().

00018 {
00019         Reset();
00020 }


Member Function Documentation

void TruthCompareAna::ana ( ParticleObjectHolder poh,
TruthCompare t 
)

Definition at line 33 of file TruthCompareAna.cxx.

References ComputeStatistics(), TruthCompare::emenergy, RecCandHeader::GetEvent(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, it, Msg::kDebug, Msg::kError, NtpMCStdHep::mass, TruthCompare::matchangle, matchangle, matched, MatchRecoParticle(), ParticleObjectHolder::mctrue, meupergev, MSG, NtpMCStdHep::p4, NtpMCStdHep::parent, ParticleObjectHolder::particles3d, Particle3D::particletype, ProcessTrueParticles(), TruthCompare::reco_matched, reco_matched_visible_e, TruthCompare::reco_not_matched, TruthCompare::reco_to_true_match, reco_visible_e, TruthCompare::Reset(), TruthCompare::stage, Particle3D::start_u, Particle3D::start_v, Particle3D::start_z, MCTrue::stdhep, Particle3D::sum_e, true_idx, TruthCompare::truelepE, and TruthCompare::truelepType.

Referenced by PRecordAna::ana(), and ParticleCheck::Run().

00034 {
00035         
00036         t->Reset();
00037         
00038         MSG("TruthCompareAna",Msg::kDebug)<<"----------------------------\n";
00039         MSG("TruthCompareAna",Msg::kDebug)<<"Run: "<<poh->GetHeader().GetRun()<<" Snarl: "<<poh->GetHeader().GetSnarl()<<" Event: "<<poh->GetHeader().GetEvent()<<"\n"; 
00040         ProcessTrueParticles(poh);
00041          MSG("TruthCompareAna",Msg::kDebug)<<matched.size()<<" true particles to be considered\n";
00042 
00043 
00044         //find lep energy from CC
00045         t->truelepE=0;
00046         t->truelepType=0;
00047         for(unsigned int i=0;i<poh->mctrue.stdhep.size();i++)
00048         {
00049                 NtpMCStdHep *mc = &(poh->mctrue.stdhep[i]);
00050                 
00051                 
00052                 if(mc->parent[0]!=0)continue;
00053                                 
00054                 if(TMath::Abs(mc->IdHEP)==11 || TMath::Abs(mc->IdHEP)==13 || TMath::Abs(mc->IdHEP)==15)
00055                 {
00056                         
00057                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00058                         if(trueE2<0)
00059                         {
00060                                 trueE2=-trueE2;
00061                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00062                         }
00063                         double trueE=sqrt(trueE2);
00064 
00065                 
00066                         t->truelepE=trueE;
00067                         t->truelepType=mc->IdHEP;
00068                         break;
00069                 }
00070         }
00071                 
00072                 
00073 
00074         //find em energy
00075         t->emenergy=0;
00076         for(unsigned int i=0;i<poh->mctrue.stdhep.size();i++)
00077         {
00078                 NtpMCStdHep *mc = &(poh->mctrue.stdhep[i]);
00079                 if(mc->IstHEP!=1)continue;
00080                 
00081                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00082                         if(trueE2<0)
00083                         {
00084                                 trueE2=-trueE2;
00085                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00086                         }
00087                         double trueE=sqrt(trueE2);
00088 
00089                 
00090 
00091                 if(TMath::Abs(mc->IdHEP)==11 || TMath::Abs(mc->IdHEP)==111 || TMath::Abs(mc->IdHEP)==22)
00092                 {
00093                         t->emenergy+=trueE;
00094                 }
00095         }
00096                 
00097                 
00098 
00099 
00100 
00101         
00102         std::map<double,int> particle_e_index;
00103         for(unsigned int i=0;i<poh->particles3d.size();i++)
00104         {
00105                 Particle3D * p = &poh->particles3d[i];
00106 
00107                 particle_e_index.insert(std::pair<double,int>(p->sum_e,i));
00108         
00109         }
00110         
00111         std::map<double,int>::reverse_iterator it;
00112         for(it=particle_e_index.rbegin();it!=particle_e_index.rend();it++)
00113         {
00114                 Particle3D * p = &poh->particles3d[it->second];
00115                 MSG("TruthCompareAna",Msg::kDebug)<<"looking for match for recoed particle type "<<p->particletype<<" with e "<<p->sum_e/meupergev<<" ("<<p->start_u<<", "<<p->start_v<<", "<<p->start_z<<")\n";
00116                 std::vector<int> truematch = MatchRecoParticle(p,poh);
00117                 
00118                 reco_visible_e+=p->sum_e/meupergev;
00119                 
00120                 if(truematch.size()>0)
00121                 {
00122                         t->reco_to_true_match.push_back(truematch);
00123                         t->reco_matched.push_back(it->second);
00124                         t->matchangle=matchangle;
00125                 
00126                         MSG("TruthCompareAna",Msg::kDebug)<<"reco particle "<<p->particletype<<" with e "<<p->sum_e/meupergev<<" matched to:\n";
00127                         for(unsigned int i=0;i<truematch.size();i++)
00128                         {
00129                                 NtpMCStdHep *mc = &(poh->mctrue.stdhep[truematch[i]]);
00130 
00131                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00132                         if(trueE2<0)
00133                         {
00134                                 trueE2=-trueE2;
00135                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00136                         }
00137                         double trueE=sqrt(trueE2);
00138 
00139 
00140 
00141                                 MSG("TruthCompareAna",Msg::kDebug)<<"\t true particle "<<mc->IdHEP<<" with e "<<trueE<<"\n";
00142                         }
00143                 
00144                 
00145                         //mark the true particles as matched...
00146                         for(unsigned int k=0;k<truematch.size();k++)
00147                                 for(unsigned int p=0;p<matched.size();p++)
00148                                         if(true_idx[p]==truematch[k])matched[p]=truematch[k];
00149                         reco_matched_visible_e=p->sum_e/meupergev;
00150                         
00151                         int left=0;for(unsigned int i=0;i<matched.size();i++)if(!matched[i])left++;
00152                         MSG("TruthCompareAna",Msg::kDebug)<<"true particles left "<<left<<"\n";
00153                 }else{
00154                         t->reco_not_matched.push_back(it->second);
00155                 }
00156                 
00157                 
00158         }       
00159         
00160 
00161                 for(unsigned int i=0;i<true_idx.size();i++)
00162                 {
00163                         if(matched[i])continue;//already used!
00164                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_idx[i]]);
00165 
00166 
00167                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00168                         if(trueE2<0)
00169                         {
00170                                 trueE2=-trueE2;
00171                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00172                         }
00173                         double trueE=sqrt(trueE2);
00174 
00175 
00176                         MSG("TruthCompareAna",Msg::kDebug)<<"Did NOT match true type "<<mc->IdHEP<< " with e "<<trueE<<"\n";
00177                 }                       
00178 
00179                 
00180         
00181         
00182         //we've done the matching...
00183         //compute some statistical values about the matching now...
00184         ComputeStatistics(poh,t);
00185         
00186         
00187         t->stage=1;
00188         
00189         MSG("TruthCompareAna",Msg::kDebug)<<"\n";
00190 }

void TruthCompareAna::ComputeStatistics ( ParticleObjectHolder poh,
TruthCompare t 
) [private]

Definition at line 194 of file TruthCompareAna.cxx.

References TruthCompare::frac_e_found, TruthCompare::frac_particles_found, NtpMCStdHep::IdHEP, Msg::kDebug, Msg::kError, NtpMCStdHep::mass, matched, ParticleObjectHolder::mctrue, meupergev, MSG, NtpMCStdHep::p4, ParticleObjectHolder::particles3d, TruthCompare::particles_matched_to_true, Particle3D::particletype, TruthCompare::possible_true_particles, TruthCompare::reco_matched, reco_matched_visible_e, TruthCompare::reco_matched_visible_e, TruthCompare::reco_to_true_match, TruthCompare::reco_visible_e, reco_visible_e, MCTrue::stdhep, Particle3D::sum_e, true_idx, TruthCompare::true_not_recoed, TruthCompare::true_visible_e, and MCTrue::visible_energy.

Referenced by ana().

00195 {
00196 
00197          t->frac_particles_found=0;
00198 
00199         for(unsigned int i=0;i<matched.size();i++)
00200                 if(!matched[i])t->true_not_recoed.push_back(true_idx[i]);
00201                 else t->frac_particles_found++;
00202 
00203         t->particles_matched_to_true=t->reco_to_true_match.size();
00204         t->possible_true_particles=matched.size();
00205         
00206         if(t->frac_particles_found)t->frac_particles_found  /= t->possible_true_particles;
00207         
00208         t->reco_matched_visible_e=reco_matched_visible_e;
00209         t->reco_visible_e=reco_visible_e;
00210         
00211         t->true_visible_e=poh->mctrue.visible_energy;
00212 
00213         if(t->true_visible_e)t->frac_e_found=t->reco_matched_visible_e/t->true_visible_e;
00214 
00215 
00216         for(unsigned int i=0;i<t->reco_matched.size();i++)
00217         {
00218                 Particle3D * p = &poh->particles3d[t->reco_matched[i]];
00219                 MSG("TruthCompareAna",Msg::kDebug)<<"Matched reco particle type"<<p->particletype<<" with e "<<p->sum_e/meupergev<<"\n";
00220                 std::vector<int> trues=t->reco_to_true_match[i];
00221                 for(unsigned int j=0;j<trues.size();j++)
00222                 {
00223                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[trues[j]]);
00224 
00225                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00226                         if(trueE2<0)
00227                         {
00228                                 trueE2=-trueE2;
00229                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00230                         }
00231                         double trueE=sqrt(trueE2);
00232 
00233 
00234                         MSG("TruthCompareAna",Msg::kDebug)<<"\tto type: "<<mc->IdHEP<<" e "<<trueE<<"\n";
00235                 }
00236                 
00237         }
00238 
00239 
00240 }

std::vector< int > TruthCompareAna::MatchRecoElectron ( Particle3D p,
ParticleObjectHolder poh 
) [private]

Definition at line 291 of file TruthCompareAna.cxx.

References Anp::angle(), MuELoss::e, Particle3D::end_z, Particle3D::entries, NtpMCStdHep::IdHEP, it, Msg::kDebug, Msg::kError, NtpMCStdHep::mass, matchangle, matched, ParticleObjectHolder::mctrue, meupergev, MSG, n, NtpMCStdHep::p4, Particle3D::start_z, MCTrue::stdhep, Particle3D::sum_e, true_idx, Particle3D::u, Particle3D::v, and Particle3D::z.

Referenced by MatchRecoParticle().

00292 {
00293         std::vector<int>truematch;
00294 
00295                         //calculate some numbers about this particle....
00296                 MSG("TruthCompareAna",Msg::kDebug)<<"attempting to match elec\n";       
00297 
00298                         double sum_z2=0;
00299                         double sum_z=0;
00300                         double sum_u=0;
00301                         double sum_zu=0;                        
00302                         double sum_v=0;
00303                         double sum_zv=0;
00304                         int n=0;        
00305                                                 
00306                         for(int i=0;i<p->entries;i++)
00307                         {
00308                                 sum_z2+=p->z[i]*p->z[i];
00309                                 sum_z+=p->z[i];
00310                                 sum_zu+=p->z[i]*p->u[i];
00311                                 sum_u+=p->u[i];
00312                                 sum_zv+=p->z[i]*p->v[i];
00313                                 sum_v+=p->v[i]; 
00314                                 n++;                    
00315                         }
00316                 
00317                         if(n<2)return truematch;
00318         
00319                         double off_u=0, slope_u=0, off_v=0, slope_v=0;
00320                         if(n*sum_z2-sum_z*sum_z>1e-10)
00321                         {
00322                                 off_u=  (sum_u*sum_z2 - sum_z*sum_zu)/(n*sum_z2-sum_z*sum_z);
00323                                 slope_u=  (n*sum_zu - sum_z*sum_u)/(n*sum_z2-sum_z*sum_z);
00324                                                 
00325                                 off_v=  (sum_v*sum_z2 - sum_z*sum_zv)/(n*sum_z2-sum_z*sum_z);
00326                                 slope_v=  (n*sum_zv - sum_z*sum_v)/(n*sum_z2-sum_z*sum_z);
00327                         }
00328 
00330                         double dz=p->end_z-p->start_z;
00331                         double du=dz*slope_u;
00332                         double dv=dz*slope_v;                   
00333                         
00334                         double r2 = du*du+dv*dv+dz*dz;
00335                         r2=sqrt(r2);
00336                         double pu= r2 ? p->sum_e/meupergev * du/r2 : 0;
00337                         double pv= r2 ? p->sum_e/meupergev * dv/r2 : 0;
00338                         double pz= r2 ? p->sum_e/meupergev * dz/r2 : 0;
00339                         
00340         
00341         MSG("TruthCompareAna",Msg::kDebug)<<"recoed elec puvz "<<pu<<" "<<pv<<" "<<pz<<"\n";
00342                         
00343                         //look for large true electrons
00344                         std::map<double,int>true_elec;
00345                         std::map<double,int>true_pi0;
00346                         std::map<double,int>true_gamma;
00347                         std::map<double,int>true_proton;
00348                                                                         
00349                 for(unsigned int i=0;i<true_idx.size();i++)
00350                 {
00351                         if(matched[i])continue;//already used!
00352                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_idx[i]]);
00353 
00354                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00355                         if(trueE2<0)
00356                         {
00357                                 trueE2=-trueE2;
00358                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00359                         }
00360                         double trueE=sqrt(trueE2);
00361 
00362 
00363                         if(TMath::Abs(mc->IdHEP)==11)true_elec.insert(std::pair<double,int>(trueE,true_idx[i]));
00364                         if(TMath::Abs(mc->IdHEP)==111)true_pi0.insert(std::pair<double,int>(trueE,true_idx[i]));
00365                         if(TMath::Abs(mc->IdHEP)==22)true_gamma.insert(std::pair<double,int>(trueE,true_idx[i]));
00366                         if(TMath::Abs(mc->IdHEP)==2212)true_proton.insert(std::pair<double,int>(trueE,true_idx[i]));
00367                 }                       
00368 
00369                 MSG("TruthCompareAna",Msg::kDebug)<<"true elec: "<<true_elec.size()<<" pi0: "<<true_pi0.size()<<" gamma: "<<true_gamma.size()<<" prot "<<true_proton.size()<<"\n";
00370                 
00371                 std::map<double,int>::reverse_iterator it;
00372                 
00373                 double trueE_matched=0;
00374                 
00375                 //need to add up true vectors to get potential angle!
00376                 double cur_pu=0;
00377                 double cur_pv=0;
00378                 double cur_pz=0;
00379                 
00380                 for(it=true_elec.rbegin();it!=true_elec.rend();it++)
00381                 {               
00382                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00383 
00384                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pu;
00385                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pv;
00386                         double true_pz=mc->p4[2]+cur_pz;
00387 
00388         MSG("TruthCompareAna",Msg::kDebug)<<"current true puvz "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00389                         
00390                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00391                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00392                                         
00393                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00394                         
00395                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00396                         if(trueE2<0)
00397                         {
00398                                 trueE2=-trueE2;
00399                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00400                         }
00401                         double trueE=sqrt(trueE2);
00402 
00403 
00404                         
00405                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00406 
00407                         //require somewhat the same direction
00408                         if(angle<0.5)
00409                         {
00410                                 //matched
00411                                 truematch.push_back(it->second);
00412                                 trueE_matched+=trueE;
00413                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00414                                 cur_pu=true_pu;
00415                                 cur_pv=true_pv;
00416                                 cur_pz=true_pz;
00417                         }else{
00418                                 MSG("TruthCompareAna",Msg::kDebug)<<"\n";
00419                         }
00420                         
00421                         if( trueE_matched> 0.9 *p->sum_e/meupergev)
00422                         {
00423                                 matchangle=angle;
00424                                 return truematch;
00425                         }
00426                 }
00427                 
00428                 
00429                 for(it=true_pi0.rbegin();it!=true_pi0.rend();it++)
00430                 {               
00431                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00432 
00433                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pu;
00434                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pv;
00435                         double true_pz=mc->p4[2]+cur_pz;
00436 
00437         MSG("TruthCompareAna",Msg::kDebug)<<"current true puvz "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00438                         
00439                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00440                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00441                                         
00442                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00443                         
00444                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00445                         if(trueE2<0)
00446                         {
00447                                 trueE2=-trueE2;
00448                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00449                         }
00450 
00451                         double trueE=sqrt(trueE2);
00452 
00453                         
00454                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00455 
00456                         //require somewhat the same direction
00457                         if(angle<0.5)
00458                         {
00459                                 //matched
00460                                 truematch.push_back(it->second);
00461                                 trueE_matched+=trueE;
00462                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00463                                 cur_pu=true_pu;
00464                                 cur_pv=true_pv;
00465                                 cur_pz=true_pz;
00466                         }
00467                         
00468                         if( trueE_matched> 0.9 *p->sum_e/meupergev)
00469                         {
00470                                 matchangle=angle;
00471                                 return truematch;
00472                         }
00473                 }
00474                                 
00475                 
00476                 
00477                 for(it=true_gamma.rbegin();it!=true_gamma.rend();it++)
00478                 {               
00479                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00480 
00481                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pu;
00482                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pv;
00483                         double true_pz=mc->p4[2]+cur_pz;
00484 
00485         MSG("TruthCompareAna",Msg::kDebug)<<"current true puvz "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00486                         
00487                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00488                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00489                                         
00490                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00491                         
00492 
00493                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00494                         if(trueE2<0)
00495                         {
00496                                 trueE2=-trueE2;
00497                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00498                         }
00499                         double trueE=sqrt(trueE2);
00500                 
00501                         
00502                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00503 
00504                         //require somewhat the same direction
00505                         if(angle<0.5)
00506                         {
00507                                 //matched
00508                                 truematch.push_back(it->second);
00509                                 trueE_matched+=trueE;
00510                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00511                                 cur_pu=true_pu;
00512                                 cur_pv=true_pv;
00513                                 cur_pz=true_pz;
00514                         }
00515                         
00516                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00517                                 matchangle=angle;
00518                                 return truematch;
00519                         }
00520                 }
00521                                                 
00522 
00523                 
00524                 for(it=true_proton.rbegin();it!=true_proton.rend();it++)
00525                 {               
00526                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00527 
00528                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pu;
00529                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0)+cur_pv;
00530                         double true_pz=mc->p4[2]+cur_pz;
00531 
00532         MSG("TruthCompareAna",Msg::kDebug)<<"current true puvz "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00533                         
00534                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00535                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00536                                         
00537                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00538                         
00539                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00540                         if(trueE2<0)
00541                         {
00542                                 trueE2=-trueE2;
00543                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00544                         }
00545                         double trueE=sqrt(trueE2);
00546 
00547                         
00548                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00549 
00550                         //require somewhat the same direction
00551                         if(angle<0.5)
00552                         {
00553                                 //matched
00554                                 truematch.push_back(it->second);
00555                                 trueE_matched+=trueE;
00556                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00557                                 cur_pu=true_pu;
00558                                 cur_pv=true_pv;
00559                                 cur_pz=true_pz;
00560                         }
00561                         
00562                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00563                                 matchangle=angle;
00564                                 return truematch;
00565                         }
00566                 }
00567                                         
00568                         
00569         return truematch;               
00570 }

std::vector< int > TruthCompareAna::MatchRecoMuon ( Particle3D p,
ParticleObjectHolder poh 
) [private]

Definition at line 574 of file TruthCompareAna.cxx.

References Anp::angle(), Particle3D::end_z, Particle3D::entries, NtpMCStdHep::IdHEP, it, Msg::kDebug, Msg::kError, NtpMCStdHep::mass, matchangle, matched, ParticleObjectHolder::mctrue, meupergev, MSG, n, NtpMCStdHep::p4, Particle3D::start_z, MCTrue::stdhep, Particle3D::sum_e, true_idx, Particle3D::u, Particle3D::v, and Particle3D::z.

Referenced by MatchRecoParticle().

00575 {
00576         std::vector<int>truematch;
00577 
00578                         //calculate some numbers about this particle....
00579                 MSG("TruthCompareAna",Msg::kDebug)<<"attempting to match muon\n";       
00580 
00581                         double sum_z2=0;
00582                         double sum_z=0;
00583                         double sum_u=0;
00584                         double sum_zu=0;                        
00585                         double sum_v=0;
00586                         double sum_zv=0;
00587                         int n=0;        
00588                                                 
00589                         int ent=p->entries<10?p->entries:10;                    
00590                         for(int i=0;i<ent;i++)
00591                         {
00592                                 sum_z2+=p->z[i]*p->z[i];
00593                                 sum_z+=p->z[i];
00594                                 sum_zu+=p->z[i]*p->u[i];
00595                                 sum_u+=p->u[i];
00596                                 sum_zv+=p->z[i]*p->v[i];
00597                                 sum_v+=p->v[i]; 
00598                                 n++;                    
00599                         }
00600                 
00601                         if(n<2)return truematch;
00602 
00603         
00604                         double denom = (n*sum_z2-sum_z*sum_z);
00605 
00606                         double off_u=0, slope_u=0, off_v=0, slope_v =0;
00607 
00608                         if(denom)
00609                         {
00610                                 off_u=  (sum_u*sum_z2 - sum_z*sum_zu)/(n*sum_z2-sum_z*sum_z);
00611                                 slope_u=  (n*sum_zu - sum_z*sum_u)/(n*sum_z2-sum_z*sum_z);
00612 
00613                                 off_v=  (sum_v*sum_z2 - sum_z*sum_zv)/(n*sum_z2-sum_z*sum_z);
00614                                 slope_v=  (n*sum_zv - sum_z*sum_v)/(n*sum_z2-sum_z*sum_z);
00615                         }                       
00616 
00617 
00619                         double dz=p->end_z-p->start_z;
00620                         double du=dz*slope_u;
00621                         double dv=dz*slope_v;                   
00622                         
00623                         double r2 = du*du+dv*dv+dz*dz;
00624                         r2=sqrt(r2);
00625                         double pu= r2 ? p->sum_e/meupergev * du/r2 : 0;
00626                         double pv= r2 ? p->sum_e/meupergev * dv/r2 : 0;
00627                         double pz= r2 ? p->sum_e/meupergev * dz/r2 : 0;
00628                         
00629                         
00630                         //look for large true electrons
00631                         std::map<double,int>true_muon;
00632                         std::map<double,int>true_charged_pion;
00633                                                                         
00634                 for(unsigned int i=0;i<true_idx.size();i++)
00635                 {
00636                         if(matched[i])continue;//already used!
00637                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_idx[i]]);
00638                         
00639 
00640                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00641                         if(trueE2<0)
00642                         {
00643                                 trueE2=-trueE2;
00644                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00645                         }
00646                         double trueE=sqrt(trueE2);
00647 
00648 
00649                         if(TMath::Abs(mc->IdHEP)==13)true_muon.insert(std::pair<double,int>(trueE,true_idx[i]));
00650                         if(TMath::Abs(mc->IdHEP)==211 || TMath::Abs(mc->IdHEP)==-211) 
00651                                 true_charged_pion.insert(std::pair<double,int>(trueE,true_idx[i]));
00652                 }
00653 
00654 
00655                 std::map<double,int>::reverse_iterator it;
00656                 
00657                 double trueE_matched=0;
00658                 
00659                 for(it=true_muon.rbegin();it!=true_muon.rend();it++)
00660                 {               
00661                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00662 
00663                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00664                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00665                         double true_pz=mc->p4[2];
00666                         
00667                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00668                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00669                                         
00670                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00671                         
00672                 
00673                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00674                         if(trueE2<0)
00675                         {
00676                                 trueE2=-trueE2;
00677                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00678                         }
00679                         double trueE=sqrt(trueE2);
00680 
00681 
00682                         
00683                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00684                         MSG("TruthCompareAna",Msg::kDebug)<<"reco p "<<pu<<" "<<pv<<" "<<pz<<" true p "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00685                         
00686                         //require somewhat the same direction
00687                         if(angle<0.5)
00688                         {
00689                                 //matched
00690                                 truematch.push_back(it->second);
00691                                 trueE_matched+=trueE;
00692                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00693                         }
00694                         
00695                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00696                                 matchangle=angle;
00697                                 return truematch;
00698                         }
00699                 }
00700                 
00701                 for(it=true_charged_pion.rbegin();it!=true_charged_pion.rend();it++)
00702                 {               
00703                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00704 
00705                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00706                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00707                         double true_pz=mc->p4[2];
00708                         
00709                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00710                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00711                                         
00712                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00713                         
00714                         
00715                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00716                         if(trueE2<0)
00717                         {
00718                                 trueE2=-trueE2;
00719                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00720                         }
00721                         double trueE=sqrt(trueE2);
00722         
00723                 
00724                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00725                         MSG("TruthCompareAna",Msg::kDebug)<<"reco p "<<pu<<" "<<pv<<" "<<pz<<" true p "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00726                         
00727                         
00728                         //require somewhat the same direction
00729                         if(angle<0.5)
00730                         {
00731                                 //matched
00732                                 truematch.push_back(it->second);
00733                                 trueE_matched+=trueE;
00734                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00735                         }
00736                         
00737                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00738                                 matchangle=angle;
00739                                 return truematch;
00740                         }
00741                 }                               
00742                         
00743         return truematch;               
00744 }

std::vector< int > TruthCompareAna::MatchRecoParticle ( Particle3D p,
ParticleObjectHolder poh 
) [private]

Definition at line 934 of file TruthCompareAna.cxx.

References Particle3D::electron, MatchRecoElectron(), MatchRecoMuon(), MatchRecoProton(), Particle3D::muon, Particle3D::neutron, Particle3D::other, Particle3D::particletype, and Particle3D::proton.

Referenced by ana().

00935 {
00936 
00937 
00938 
00939 
00940                 //now iterate over true particles for consideration........
00941                 //this may be done differently depending on the type of recoed particle!!!
00942                 std::vector<int>truematch;
00943 
00944                         if(p->particletype==Particle3D::electron ||p->particletype==Particle3D::other)truematch= MatchRecoElectron(p,poh);
00945 
00946                         if(p->particletype==Particle3D::muon)truematch= MatchRecoMuon(p,poh);
00947                 
00948                         if(p->particletype==Particle3D::proton || p->particletype==Particle3D::neutron) truematch= MatchRecoProton(p,poh);
00949                 
00950                 return truematch;
00951                 
00952                 
00953 /*              
00954                 for(unsigned int i=0;i<true_idx.size();i++)
00955                 {
00956                         if(matched[i])continue;//already used!
00957                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_idx[i]]);
00958                         
00959                 
00960                 }
00961 */
00962 
00963 /*
00964                         if(p->particletype==Particle3D::electron)
00965                         {
00966                                 printf("largest elec %d with e %f  (%f %f %f) p3 (%f %f %f)\n",p->particletype,p->sum_e/meupergev,p->start_u,p->start_v,p->start_z,pu,pv,pz);
00967                                 if(true_elec_idx>-1)
00968                                 {
00969                                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_elec_idx]);
00970                         
00971                                         printf("compare to true elec %d with e %f at %d  (%f %f %f)  p3 (%f %f %f)\n",mc->IdHEP,sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass),mc->IstHEP,(mc->vtx[0]+mc->vtx[1])/sqrt(2.0),(-mc->vtx[0]+mc->vtx[1])/sqrt(2.0),mc->vtx[2],(mc->p4[0]+mc->p4[1])/sqrt(2.0),(-mc->p4[0]+mc->p4[1])/sqrt(2.0),mc->p4[2]);             
00972                                         //calculate the angle!!!!
00973                                         
00974                                         //angle between p vectors....
00975                                         
00976                                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00977                                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00978                                         double true_pz=mc->p4[2];
00979                                         
00980                                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00981                                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00982                                         
00983                                         double angle = acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) );
00984                                         printf("angle between: %f\n",angle);
00985                                         
00986                                 //      hist_prim_angle_vs_prim_eres[type]->Fill(angle,(sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass)-p->sum_e/meupergev)/(sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass)),weight);          
00987                                 //      hist_prim_angle_vs_prim_true[type]->Fill(angle,sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass));    
00988                                         
00989                                         
00990                                 }
00991                         }
00992 */                      
00993         
00994 }

std::vector< int > TruthCompareAna::MatchRecoProton ( Particle3D p,
ParticleObjectHolder poh 
) [private]

Definition at line 748 of file TruthCompareAna.cxx.

References Anp::angle(), MuELoss::e, Particle3D::end_z, Particle3D::entries, ParticleObjectHolder::event, NtpMCStdHep::IdHEP, it, Msg::kDebug, Msg::kError, NtpMCStdHep::mass, matchangle, matched, ParticleObjectHolder::mctrue, meupergev, MSG, n, NtpMCStdHep::p4, Particle3D::start_z, MCTrue::stdhep, Particle3D::sum_e, true_idx, Particle3D::u, Particle3D::v, ParticleEvent::vtx_u, ParticleEvent::vtx_v, ParticleEvent::vtx_z, and Particle3D::z.

Referenced by MatchRecoParticle().

00749 {
00750         std::vector<int>truematch;
00751 
00752                         //calculate some numbers about this particle....
00753                 MSG("TruthCompareAna",Msg::kDebug)<<"attempting to match muon\n";       
00754 
00755                         double sum_z2=0;
00756                         double sum_z=0;
00757                         double sum_u=0;
00758                         double sum_zu=0;                        
00759                         double sum_v=0;
00760                         double sum_zv=0;
00761                         int n=0;        
00762                                                 
00763                         int ent=p->entries<10?p->entries:10;                    
00764                         
00765                         double off_u=  0;
00766                         double slope_u=  0;
00767                         
00768                         
00769                         double off_v=  0;
00770                         double slope_v=0;
00771                         
00772                         if(ent>1)
00773                         {
00774                                 for(int i=0;i<ent;i++)
00775                                 {
00776                                         sum_z2+=p->z[i]*p->z[i];
00777                                         sum_z+=p->z[i];
00778                                         sum_zu+=p->z[i]*p->u[i];
00779                                         sum_u+=p->u[i];
00780                                         sum_zv+=p->z[i]*p->v[i];
00781                                         sum_v+=p->v[i]; 
00782                                         n++;                    
00783                                 }
00784                         
00785                                 if(n*sum_z2-sum_z*sum_z)
00786                                 {
00787                                         off_u=  (sum_u*sum_z2 - sum_z*sum_zu)/(n*sum_z2-sum_z*sum_z);
00788                                         slope_u=  (n*sum_zu - sum_z*sum_u)/(n*sum_z2-sum_z*sum_z);
00789         
00790                                         off_v=  (sum_v*sum_z2 - sum_z*sum_zv)/(n*sum_z2-sum_z*sum_z);
00791                                         slope_v=  (n*sum_zv - sum_z*sum_v)/(n*sum_z2-sum_z*sum_z);
00792                                 }
00793                         }else if(TMath::Abs(p->z[0]-poh->event.vtx_z)>1e-5){
00794                                 slope_u=(p->u[0]-poh->event.vtx_u)/(p->z[0]-poh->event.vtx_z);
00795                                 slope_v=(p->v[0]-poh->event.vtx_v)/(p->z[0]-poh->event.vtx_z);
00796                                 
00797                                 off_u=poh->event.vtx_u-slope_u*poh->event.vtx_z;
00798                                 off_v=poh->event.vtx_v-slope_v*poh->event.vtx_z;
00799                                 
00800                         }
00801                         
00802                         
00803 
00805                         double dz=p->end_z-p->start_z;
00806                         double du=dz*slope_u;
00807                         double dv=dz*slope_v;                   
00808                         
00809                         double r2 = du*du+dv*dv+dz*dz;
00810                         r2=sqrt(r2);
00811                         double pu= r2? p->sum_e/meupergev * du/r2:0;
00812                         double pv= r2? p->sum_e/meupergev * dv/r2:0;
00813                         double pz= r2? p->sum_e/meupergev * dz/r2:0;
00814                         
00815                         
00816                         //look for large true electrons
00817                         std::map<double,int>true_proton;
00818                         std::map<double,int>true_neutron;
00819                                                                         
00820                 for(unsigned int i=0;i<true_idx.size();i++)
00821                 {
00822                         if(matched[i])continue;//already used!
00823                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[true_idx[i]]);
00824                         
00825                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00826                         if(trueE2<0)
00827                         {
00828                                 trueE2=-trueE2;
00829                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0   = "<<trueE2<<"\n";
00830                         }
00831                         double trueE=sqrt(trueE2);
00832 
00833                         if(TMath::Abs(mc->IdHEP)==2212)true_proton.insert(std::pair<double,int>(trueE,true_idx[i]));
00834                         if(TMath::Abs(mc->IdHEP)==2112 ) true_neutron.insert(std::pair<double,int>(trueE,true_idx[i]));
00835                 }
00836 
00837 
00838                 std::map<double,int>::reverse_iterator it;
00839                 
00840                 double trueE_matched=0;
00841                 
00842                 for(it=true_proton.rbegin();it!=true_proton.rend();it++)
00843                 {               
00844                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00845 
00846                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00847                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00848                         double true_pz=mc->p4[2];
00849                         
00850                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00851                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00852                                         
00853                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00854                 
00855 
00856                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00857                         if(trueE2<0)
00858                         {
00859                                 trueE2=-trueE2;
00860                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00861                         }
00862                         double trueE=sqrt(trueE2);
00863         
00864 
00865                         
00866                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00867                         MSG("TruthCompareAna",Msg::kDebug)<<"reco p "<<pu<<" "<<pv<<" "<<pz<<" true p "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00868                         
00869                         //require somewhat the same direction
00870                         if(angle<0.5)
00871                         {
00872                                 //matched
00873                                 truematch.push_back(it->second);
00874                                 trueE_matched+=trueE;
00875                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00876                         }
00877                         
00878                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00879                                 matchangle=angle;
00880                                 return truematch;
00881                         }
00882                 }
00883                 
00884                 for(it=true_neutron.rbegin();it!=true_neutron.rend();it++)
00885                 {               
00886                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[it->second]);
00887 
00888                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00889                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00890                         double true_pz=mc->p4[2];
00891                         
00892                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00893                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00894                                         
00895                         double angle = reco_r ? acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) ) : 0;
00896 
00897                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00898                         if(trueE2<0)
00899                         {
00900                                 trueE2=-trueE2;
00901                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00902                         }
00903                         double trueE=sqrt(trueE2);
00904 
00905                         
00906 
00907                         
00908                         MSG("TruthCompareAna",Msg::kDebug)<<"comparing with angle "<<angle <<" trueE "<<trueE<<" recoE "<<p->sum_e/meupergev<<"\n";
00909                         MSG("TruthCompareAna",Msg::kDebug)<<"reco p "<<pu<<" "<<pv<<" "<<pz<<" true p "<<true_pu<<" "<<true_pv<<" "<<true_pz<<"\n";
00910                         
00911                         
00912                         //require somewhat the same direction
00913                         if(angle<0.5)
00914                         {
00915                                 //matched
00916                                 truematch.push_back(it->second);
00917                                 trueE_matched+=trueE;
00918                                 MSG("TruthCompareAna",Msg::kDebug)<<" matched!\n";
00919                         }
00920                         
00921                         if( trueE_matched> 0.9 *p->sum_e/meupergev){
00922                                 matchangle=angle;
00923                                 return truematch;
00924                         }
00925                 }                               
00926                         
00927         return truematch;               
00928 }

void TruthCompareAna::ProcessTrueParticles ( ParticleObjectHolder poh  )  [private]

Definition at line 244 of file TruthCompareAna.cxx.

References NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, Msg::kError, NtpMCStdHep::mass, matched, ParticleObjectHolder::mctrue, MSG, NtpMCStdHep::p4, MCTrue::stdhep, and true_idx.

Referenced by ana().

00245 {
00247 //              printf("\nstdhep....run %d snarl %d event %d\n",poh->GetHeader().GetRun() , poh->GetHeader().GetSnarl() , poh->GetHeader().GetEvent());
00248 
00249 
00250                 for(unsigned int i=0;i<poh->mctrue.stdhep.size();i++)
00251                 {
00252                         NtpMCStdHep *mc = &(poh->mctrue.stdhep[i]);
00253                         
00254                         if(mc->IstHEP!=1)continue;
00255 
00256 
00257                         double trueE2 = mc->p4[3]*mc->p4[3]-mc->mass*mc->mass;
00258                         if(trueE2<0)
00259                         {
00260                                 trueE2=-trueE2;
00261                                 MSG("TruthCompareAna",Msg::kError)<<"true E2 < 0 = "<<trueE2<<"\n";
00262                         }
00263                         double trueE=sqrt(trueE2);
00264 
00265                         
00266                         //remove particles too small to consider!
00267                         if(TMath::Abs(mc->IdHEP)==13 || TMath::Abs(mc->IdHEP)==211)
00268                         {
00269                                 if(trueE<0.5)continue;
00270                         }else{
00271                                 if(trueE<0.2)continue;
00272                         }
00273                         
00274                         if(TMath::Abs(mc->IdHEP)==12 || TMath::Abs(mc->IdHEP)==14 || TMath::Abs(mc->IdHEP)==16)continue;//can't match a neutrino!
00275                         
00276                         //we'll consider it!
00277                         true_idx.push_back(i);
00278                         matched.push_back(0);
00279                         
00280         //              MSG("TruthCompareAna",Msg::kDebug)<<"true type "<<mc->IdHEP<<" e "<<trueE<<" from "<<mc->dethit[0].vtx[0]<<" "<<mc->dethit[0].vtx[1]<<" "<<mc->dethit[0].vtx[2]<<" to "<<mc->dethit[1].vtx[0]<<" "<<mc->dethit[1].vtx[1]<<" "<<mc->dethit[1].vtx[2]<<"\n";
00281                 }
00282 //              printf("\n");
00283                 
00284                 
00286                 
00287         
00288 
00289 }

void TruthCompareAna::Reset (  ) 

Definition at line 22 of file TruthCompareAna.cxx.

References matched, reco_matched_visible_e, reco_visible_e, and true_idx.

Referenced by TruthCompareAna(), and ~TruthCompareAna().

00023 {
00024         
00025         
00026         true_idx.clear();
00027         matched.clear();
00028         reco_matched_visible_e=0;
00029         reco_visible_e=0;
00030 
00031 }


Member Data Documentation

double TruthCompareAna::matchangle [private]

Definition at line 38 of file TruthCompareAna.h.

Referenced by ana(), MatchRecoElectron(), MatchRecoMuon(), and MatchRecoProton().

std::vector<int> TruthCompareAna::matched [private]
double TruthCompareAna::meupergev [private]

Definition at line 35 of file TruthCompareAna.h.

Referenced by ana(), ComputeStatistics(), and Reset().

Definition at line 36 of file TruthCompareAna.h.

Referenced by ana(), ComputeStatistics(), and Reset().

Definition at line 41 of file TruthCompareAna.h.

std::vector<int> TruthCompareAna::true_idx [private]

The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1