loon: TruthCompareAna.cxx Source File

00001 #include "NueAna/ParticlePID/ParticleAna/TruthCompareAna.h"
00002 #include "TMath.h"
00003 
00004 #include "MCNtuple/NtpMCStdHep.h"
00005 #include <map>
00006 #include "MessageService/MsgService.h"
00007 
00008 CVSID("$Id: TruthCompareAna.cxx,v 1.3 2009/06/24 23:06:46 gmieg Exp $");
00009 
00010 TruthCompareAna::TruthCompareAna()
00011 {
00012         Reset();
00013         
00014         meupergev = 25.;//should use something more precise here!
00015 }
00016 
00017 TruthCompareAna::~TruthCompareAna()
00018 {
00019         Reset();
00020 }
00021 
00022 void TruthCompareAna::Reset()
00023 {
00024         
00025         
00026         true_idx.clear();
00027         matched.clear();
00028         reco_matched_visible_e=0;
00029         reco_visible_e=0;
00030 
00031 }
00032 
00033 void TruthCompareAna::ana(ParticleObjectHolder * poh, TruthCompare * t)
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 }
00191 
00192 
00193 
00194 void TruthCompareAna::ComputeStatistics(ParticleObjectHolder *poh, TruthCompare *t)
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 }
00241 
00242 
00243 //prepare the lists of the true particles that we should attempt to match
00244 void TruthCompareAna::ProcessTrueParticles(ParticleObjectHolder *poh)
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 }
00290 
00291 std::vector<int> TruthCompareAna::MatchRecoElectron(Particle3D *p,ParticleObjectHolder *poh)
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 }
00571 
00572 
00573 
00574 std::vector<int> TruthCompareAna::MatchRecoMuon(Particle3D *p,ParticleObjectHolder *poh)
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 }
00745 
00746 
00747 
00748 std::vector<int> TruthCompareAna::MatchRecoProton(Particle3D *p,ParticleObjectHolder *poh)
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 }
00929 
00930 
00931 
00932 
00933 
00934 std::vector<int> TruthCompareAna::MatchRecoParticle(Particle3D *p, ParticleObjectHolder *poh)
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 }