PETrimmer Class Reference

#include <PETrimmer.h>

Inheritance diagram for PETrimmer:
JobCModule

List of all members.

Public Member Functions

 PETrimmer ()
 ~PETrimmer ()
JobCResult Reco (MomNavigator *mom)
void TrimRecord (NtpStRecord *str)
void EndJob ()
void BeginJob ()
void Config (const Registry &)
const RegistryDefaultConfig () const
void DumpIt (NtpStRecord *str)
int CleanEvent (NtpSREvent *evt, std::vector< NtpSRShower * > shw, std::vector< NtpSRTrack * > trk, std::vector< const NtpSRStrip * > stp, int *pepass, float *stripe_trl, float *stripe_shw, int *goodshower, int *goodtrack)
int CleanShower (NtpSRShower *evt, std::vector< const NtpSRStrip * > stp, int *pepass, float *stripe)
int CleanTrack (NtpSRTrack *evt, std::vector< const NtpSRStrip * > stp, int *pepass, float *stripe)

Private Attributes

float pecut
float updateEventEnergy

Detailed Description

Definition at line 26 of file PETrimmer.h.


Constructor & Destructor Documentation

PETrimmer::PETrimmer (  ) 

Definition at line 50 of file PETrimmer.cxx.

00050                     :
00051   pecut(0.0),
00052   updateEventEnergy(1)
00053 {}

PETrimmer::~PETrimmer (  ) 

Definition at line 57 of file PETrimmer.cxx.

00058 {}


Member Function Documentation

void PETrimmer::BeginJob ( void   )  [virtual]

Implement for notification of begin of job

Reimplemented from JobCModule.

Definition at line 61 of file PETrimmer.cxx.

00062 {
00063 
00064 }

int PETrimmer::CleanEvent ( NtpSREvent evt,
std::vector< NtpSRShower * >  shw,
std::vector< NtpSRTrack * >  trk,
std::vector< const NtpSRStrip * >  stp,
int *  pepass,
float *  stripe_trl,
float *  stripe_shw,
int *  goodshower,
int *  goodtrack 
)

incorrect... we did not change array size.. workaround for other mods...

Definition at line 507 of file PETrimmer.cxx.

References NtpSRPlane::beg, NtpSRPlane::begu, NtpSRPlane::begv, NtpSRPlane::end, NtpSRPlane::endu, NtpSRPlane::endv, NtpSRStripPulseHeight::gev, Msg::kDebug, MSG, NtpSRPlane::n, NtpSREvent::nshower, NtpSREvent::nstpcnt, NtpSREvent::nstrip, NtpSREvent::ntrack, NtpSRPlane::nu, NtpSRPlane::nv, NtpSRPulseHeight::pe, NtpSREvent::ph, NtpSREvent::plane, NtpSRPulseHeight::raw, showergreater(), NtpSREvent::shw, NtpSRPulseHeight::sigcor, NtpSRPulseHeight::siglin, NtpSRStripPulseHeight::sigmap, NtpSREvent::stp, trackgreater(), NtpSREvent::trk, and updateEventEnergy.

Referenced by TrimRecord().

00508 {
00509      int goodevent=1;
00510 
00511                 //remove bad showers from this event....
00512 
00513                 std::map<int,int> shwstpmap;
00514 
00515                 //mark bad showers as bad, record which strips are effected
00516                 for(int i=0;i<evt->nshower;i++)
00517                 {
00518                         if(!goodshower[evt->shw[i]])
00519                         {
00520                             //cout<<"Removing shw: "<<i<<endl;  
00521                                 for(int j=0;j<shw[evt->shw[i]]->nstrip;j++)
00522                                 shwstpmap[shw[evt->shw[i]]->stp[j]]=1;                                  
00523 
00524 
00525                                 evt->shw[i]=-1;
00526                                 
00527                         }       
00528                 }
00529                 
00530                 //remove bad showers from event
00531                 for(int i=0;i<evt->nshower;i++)
00532                 {
00533                         if(evt->shw[i]<0)
00534                         {
00535                                 for(int j=i+1;j<evt->nshower;j++)
00536                                         evt->shw[j-1]=evt->shw[j];
00537                                 evt->shw[evt->nshower-1]=-1;
00538                                 evt->nshower--;
00539                                 i--;
00540                         }
00541                 }
00542 
00543 
00544 
00545 
00546                 //mark bad trackss as bad, record which strips are effected
00547                 for(int i=0;i<evt->ntrack;i++)
00548                 {      
00549                         if(evt->trk[i]<0)continue; 
00550                         if(!goodtrack[evt->trk[i]])
00551                         {       
00552                            //cout<<"removing track"<<i<<endl;
00553                                 for(int j=0;j<trk[evt->trk[i]]->nstrip;j++)
00554                                 shwstpmap[trk[evt->trk[i]]->stp[j]]=1;
00555 
00556                          
00557                                 evt->trk[i]=-1;
00558                          
00559                         }
00560                 }
00561                 
00562                 //remove bad tracks from event
00563                 for(int i=0;i<evt->ntrack;i++)
00564                 {       
00565                         if(evt->trk[i]<0)
00566                         {       
00567                                 for(int j=i+1;j<evt->ntrack;j++)
00568                                         evt->trk[j-1]=evt->trk[j];
00569                                 evt->trk[evt->ntrack-1]=-1; 
00570                                evt->ntrack--;
00571                                 i--;
00572                         }
00573                 }
00574 
00575 
00576 
00577 
00578 
00579 
00580 
00581 
00582 
00583                 //which strips are in the tracks? removing these strips due to removed showers has no effect on event energy
00584                 for(int i=0;i<evt->ntrack;i++)
00585                 {
00586                         for(int j=0;j<trk[evt->trk[i]]->nstrip;j++)
00587                         {
00588                                 shwstpmap[trk[evt->trk[i]]->stp[j]]=0;
00589                         
00590                 
00591                         }
00592                 }
00593 
00594                 //if we have strips shared between good and bad showers, we are not removing those strips from the event
00595                 for(int i=0;i<evt->nshower;i++)
00596                 {
00597                        for(int j=0;j<shw[evt->shw[i]]->nstrip;j++)
00598                         {
00599                                 shwstpmap[shw[evt->shw[i]]->stp[j]]=0;
00600         
00601 
00602                         }
00603                 }
00604                 
00605 
00606 
00607                 std::map<int,int>::iterator iter;
00608                 for(iter=shwstpmap.begin();iter!=shwstpmap.end();iter++)
00609                 {
00610                         
00611                         
00612                         if(iter->second==1)
00613                         {
00614                           int j=evt->nstrip;
00615                           for(int i=0;i<evt->nstrip;i++) //find the start index in evt->stp for this bad strip
00616                             {
00617                               if(evt->stp[i]==iter->first)
00618                                 {
00619                                   j=i;
00620                                   break;
00621                                 }
00622                             }
00623                                 for(int i=j;i<evt->nstrip-1;i++)
00624                                 {
00625                                   evt->stp[i]=evt->stp[i+1]; //remove that strip from the event strip list
00626                                   
00627                                 }
00628                                 //mark the last strip, so we dont get confused
00629                                 evt->stp[evt->nstrip-1]=-1;
00630                                 evt->nstrip--; //adjust both .... we are loosing access to the end of the array
00631                                 evt->nstpcnt--;
00632                         }       
00633 
00634                 }
00635                 //cout<<"Remaining strips: "<<evt->nstrip<<endl;
00636                 if(evt->nstrip<1)return 0; //mark this event for deletion       
00637 
00638 
00639                 int cutcount =0;
00640         
00641                 //temporary vectors for saving good hits - better than erasing single entries in arrays at a time!
00642                 std::vector<float> stpph0sigmap;
00643                 std::vector<float> stpph0mip;
00644                 std::vector<float> stpph0gev;
00645                 std::vector<float> stpph1sigmap;
00646                 std::vector<float> stpph1mip;
00647                 std::vector<float> stpph1gev;
00648                 std::vector<int> stpid;
00649 
00650                 if(updateEventEnergy){
00651                 //we are recomputing the energies... clear out those variables...
00652                 evt->ph.raw =0;
00653                 evt->ph.pe =0;
00654                 evt->ph.sigcor =0;
00655                 evt->ph.siglin =0;
00656                 evt->ph.sigmap =0;
00657                 //evt->ph.mip =0; //don't change this
00658                 evt->ph.gev =0;         
00659                 }
00660 
00661                 int us=0;
00662                 int vs=0;
00663                 for(int j=0;j<evt->nstpcnt;j++)
00664                 {
00665 
00666         
00667                         if(evt->stp[j]<0)continue;
00668 
00669 
00670 
00671                         
00672                         if(pepass[evt->stp[j]])
00673                         {
00674                         
00675 
00676                                  //printf("stp %d plane %d\n",stp[evt->stp[j]],stp[evt->stp[j]]->plane);
00677 
00678 
00679                                 //keep it!
00680                                 stpid.push_back(evt->stp[j]);
00681 /*                              stpph0sigmap.push_back(evt->stpph0sigmap[j]);
00682                                 stpph1sigmap.push_back(evt->stpph1sigmap[j]);
00683                                 stpph0gev.push_back(evt->stpph0gev[j]);
00684                                 stpph1gev.push_back(evt->stpph1gev[j]);
00685                                 stpph0mip.push_back(evt->stpph0mip[j]);
00686                                 stpph1mip.push_back(evt->stpph1mip[j]);
00687 */
00688 
00689 
00690 
00691 
00692                                 if(updateEventEnergy){
00693 
00694                         //recompute the energy!!!!!!!
00695                                 evt->ph.raw += (stp[evt->stp[j]]->ph0.raw + stp[evt->stp[j]]->ph1.raw);
00696                                 evt->ph.pe += (stp[evt->stp[j]]->ph0.pe + stp[evt->stp[j]]->ph1.pe);
00697                                 evt->ph.sigcor += (stp[evt->stp[j]]->ph0.sigcor + stp[evt->stp[j]]->ph1.sigcor);
00698                                 evt->ph.siglin += (stp[evt->stp[j]]->ph0.siglin + stp[evt->stp[j]]->ph1.siglin);
00699 
00700                                 //give preference to track calibrated strips???
00701                                 //this puts preference for all tracks before all showers
00702                                 //should it primary track, primary shower, other stuff, etc...???
00703                                                 
00704                                 evt->ph.sigmap += stripe_trk[3*evt->stp[j]]>0 ? stripe_trk[3*evt->stp[j]] :stripe_shw[3*evt->stp[j]];
00705                         //      evt->ph.mip += stripe_trk[1+3*evt->stp[j]]>0 ? stripe_trk[1+3*evt->stp[j]] :stripe_shw[1+3*evt->stp[j]];
00706 //                              evt->ph.gev += stripe_trk[1+3*evt->stp[j]]>0 ? stripe_trk[1+3*evt->stp[j]] :stripe_shw[1+3*evt->stp[j]];
00707 
00708 //trying a new method of energy recalculation..... store it in gev - that is recalculated anyways!
00709 //give priority to showers, in order of energy, and then to tracks, in order of energy
00710 //this is to get energies for nue like events!
00711                                 std::vector<std::pair<float, NtpSRShower *> >showers;
00712                                 std::vector<std::pair<float, NtpSRTrack *> >tracks;
00713 
00714                                 for(int k=0;k<evt->nshower;k++)
00715                                         showers.push_back(make_pair((float)shw[evt->shw[k]]->ph.mip,shw[evt->shw[k]]));
00716                                 for(int k=0;k<evt->ntrack;k++)  
00717                                         tracks.push_back(make_pair((float)trk[evt->trk[k]]->ph.mip,trk[evt->trk[k]]));
00718 
00719 
00720                                 std::sort(showers.begin(),showers.end(),showergreater);
00721                                 std::sort(tracks.begin(),tracks.end(),trackgreater);            
00722 
00723                                 float ees[100000];
00724                                 for(int k=0;k<100000;k++)ees[k]=0;
00725 
00726                                 for(unsigned int k=0;k<showers.size();k++)
00727                                 for(int j=0;j<showers[k].second->nstrip;j++)
00728                                 {
00729                                         int sidx = showers[k].second->stp[j];
00730                                         if(sidx>0 && sidx < 100000)                             
00731                                         {
00732                                                 if(ees[sidx]==0)ees[sidx]=showers[k].second->stpph0mip[j]+showers[k].second->stpph1mip[j];
00733                                         }
00734                                 }
00735 
00736                                 for(unsigned int k=0;k<tracks.size();k++)
00737                                 for(int j=0;j<tracks[k].second->nstrip;j++)
00738                                 {
00739                                         int sidx = tracks[k].second->stp[j];
00740                                         if(sidx>0 && sidx < 100000)
00741                                         {
00742                                                 if(ees[sidx]==0)ees[sidx]=tracks[k].second->stpph0mip[j]+tracks[k].second->stpph1mip[j];
00743                                         }
00744                                 }
00745 
00746                                 float summeu=0;
00747                                 for(int k=0;k<evt->nstrip;k++)
00748                                         if(ees[evt->stp[k]] > 0 && ees[evt->stp[k]] < 100000)
00749                                                 summeu+=ees[evt->stp[k]];
00750 
00751                                 evt->ph.gev=summeu;
00752 
00753 
00754                                 }
00755 
00756                                  if(stp[evt->stp[j]]->planeview==2)
00757                                         us++;
00758                                 if(stp[evt->stp[j]]->planeview==3)
00759                                         vs++;
00760 
00761 
00762                         }else{
00763                                 cutcount++;
00764                         }
00765                 }
00766 
00767 
00768                 MSG("PETrimmer",Msg::kDebug)<<"Cut "<<cutcount<<" strips from event!"<<endl;
00769                 MSG("PETrimmer",Msg::kDebug)<<"Shower us "<<us << " vs "<< vs <<endl;
00770 
00771                 if(us==0 || vs==0)
00772                 {
00773                         goodevent=0;
00774                         return goodevent;
00775                 }
00776                 
00777 
00778 
00779 
00780 
00781 
00782 
00783 
00784 
00785                 if (cutcount==0)return goodevent;  //we didn't cut anything, so continue        
00786 
00787                 //rewrite event stp vectors
00788                 for(int j=0;j<evt->nstpcnt;j++)
00789                 {
00790                         evt->stp[j]=-1;
00791         /*              evt->stpph0sigmap[j]=0;
00792                         evt->stpph1sigmap[j]=0;
00793                         evt->stpph0gev[j]=0;
00794                         evt->stpph1gev[j]=0;
00795                         evt->stpph0mip[j]=0;
00796                         evt->stpph1mip[j]=0;*/
00797                 }
00798 
00799 
00800                 //deleteing these arrays and putting new ones in works on my mac, but not on minos machines.... why?
00801 
00802 //              delete[] evt->stp;
00803 /*              delete[] evt->stpph0sigmap;
00804                 delete[] evt->stpph1sigmap;
00805                 delete[] evt->stpph0gev;
00806                 delete[] evt->stpph1gev;
00807                 delete[] evt->stpph0mip;
00808                 delete[] evt->stpph1mip;*/
00809                 
00810 //              evt->stp = new Int_t[stpid.size()];
00811 /*              evt->stpph0sigmap = new Float_t[stpid.size()];
00812                 evt->stpph1sigmap = new Float_t[stpid.size()];
00813                 evt->stpph0gev = new Float_t[stpid.size()];
00814                 evt->stpph1gev = new Float_t[stpid.size()];
00815                 evt->stpph0mip = new Float_t[stpid.size()];
00816                 evt->stpph1mip = new Float_t[stpid.size()];     //these arrays dont seem to be in the ntpst file!
00817                 
00818 */              
00819                 int planes[1000];
00820                 for(int i=0;i<1000;i++)planes[i]=0;     
00821                 
00822                 evt->nstpcnt=stpid.size();
00823                 evt->nstrip=stpid.size();  
00824                 for(int j=0;j<(int)stpid.size();j++)
00825                 {
00826                         if(stpid[j]<0 || stpid[j]>(int)stp.size()-1)continue;
00827                         
00828                         int p=stp[stpid[j]]->plane;
00829                         if(p>-1 && p < 1000) planes[p]=stp[stpid[j]]->planeview; //2 or 3 for u or v
00830 
00831 //                      evt->stp[j]=-1;
00832                         if(stpid[j]>-1 && stpid[j] < (int)stp.size())   
00833                                 evt->stp[j]=stpid[j];
00834 /*                      evt->stpph0sigmap[j]=stpph0sigmap[j];
00835                         evt->stpph1sigmap[j]=stpph1sigmap[j];
00836                         evt->stpph0gev[j]=stpph0gev[j];
00837                         evt->stpph1gev[j]=stpph1gev[j];
00838                         evt->stpph0mip[j]=stpph0mip[j];
00839                         evt->stpph1mip[j]=stpph1mip[j];*/
00840                 }
00841 
00842 
00843                 int totplanes =0;
00844                 int begu=1000;
00845                 int begv=1000;
00846                 int endu=0;
00847                 int endv=0; 
00848                 int nu=0;
00849                 int nv=0;
00850                 for(int i=0;i<1000;i++)
00851                 {
00852                         if(planes[i]>0)totplanes++;
00853                         if(planes[i]==2)
00854                         {
00855                                 nu++;
00856                                 begu=begu<i? begu:i;
00857                                 endu=endu<i? i :endu;
00858                         }else if(planes[i]==3)
00859                         {
00860                                 nv++;
00861                                 begv=begv<i? begv:i;
00862                                 endv=endv<i?i:endv;
00863                         }
00864                 }
00865                 evt->plane.n=totplanes;
00866                 evt->plane.nu=nu;
00867                 evt->plane.nv=nv;
00868                 evt->plane.beg=begu<begv?begu:begv;
00869                 evt->plane.begu=begu;
00870                 evt->plane.begv=begv;
00871                 evt->plane.end=endu<endv?endv:endu;
00872                 evt->plane.endu=endu;
00873                 evt->plane.endv=endv;
00874 
00875 
00876         return goodevent;
00877 
00878 }

int PETrimmer::CleanShower ( NtpSRShower evt,
std::vector< const NtpSRStrip * >  stp,
int *  pepass,
float *  stripe 
)

incorrect... we did not change array size.. workaround for other mods...

Definition at line 884 of file PETrimmer.cxx.

References NtpSRStripPulseHeight::gev, Msg::kDebug, NtpSRStripPulseHeight::mip, MSG, NtpSRShower::nstpcnt, NtpSRShower::nstrip, NtpSRPulseHeight::pe, NtpSRShower::ph, NtpSRPulseHeight::raw, NtpSRPulseHeight::sigcor, NtpSRPulseHeight::siglin, NtpSRStripPulseHeight::sigmap, NtpSRShower::stp, NtpSRShower::stpph0gev, NtpSRShower::stpph0mip, NtpSRShower::stpph0sigmap, NtpSRShower::stpph1gev, NtpSRShower::stpph1mip, and NtpSRShower::stpph1sigmap.

Referenced by TrimRecord().

00885 {
00886                 int goodshower=1;
00887                 int cutcount =0;
00888         
00889                 //temporary vectors for saving good hits - better than erasing single entries in arrays at a time!
00890                 std::vector<float> stpph0sigmap;
00891                 std::vector<float> stpph0mip;
00892                 std::vector<float> stpph0gev;
00893                 std::vector<float> stpph1sigmap;
00894                 std::vector<float> stpph1mip;
00895                 std::vector<float> stpph1gev;
00896                 std::vector<int> stpid;
00897 
00898                         int us=0;
00899                         int vs=0;       
00900                 for(int j=0;j<shw->nstpcnt;j++)
00901                 {
00902                 
00903 
00904         
00905 
00906                         if(pepass[shw->stp[j]])
00907                         {
00908                         
00909                                 //keep it!
00910                                 stpid.push_back(shw->stp[j]);
00911                                 stpph0sigmap.push_back(shw->stpph0sigmap[j]);
00912                                 stpph1sigmap.push_back(shw->stpph1sigmap[j]);
00913                                 stpph0gev.push_back(shw->stpph0gev[j]);
00914                                 stpph1gev.push_back(shw->stpph1gev[j]);
00915                                 stpph0mip.push_back(shw->stpph0mip[j]);
00916                                 stpph1mip.push_back(shw->stpph1mip[j]);
00917                                 
00918                                 stripe[shw->stp[j]*3]= shw->stpph0sigmap[j]+shw->stpph1sigmap[j];
00919                                 stripe[1+shw->stp[j]*3]= shw->stpph0mip[j]+shw->stpph1mip[j];
00920                                 stripe[2+shw->stp[j]*3]= shw->stpph0gev[j]+shw->stpph1gev[j];
00921 
00922                                 if(stp[shw->stp[j]]->planeview==2)
00923                                         us++;
00924                                 if(stp[shw->stp[j]]->planeview==3)
00925                                         vs++;
00926 
00927 
00928                         }else{
00929                                 cutcount++;
00930                         
00931                                 //adjust total energy of shower
00932                                 shw->ph.raw -= (stp[shw->stp[j]]->ph0.raw + stp[shw->stp[j]]->ph1.raw);
00933                                 shw->ph.pe -= (stp[shw->stp[j]]->ph0.pe + stp[shw->stp[j]]->ph1.pe);
00934                                 shw->ph.sigcor -= (stp[shw->stp[j]]->ph0.sigcor + stp[shw->stp[j]]->ph1.sigcor);
00935                                 shw->ph.siglin -= (stp[shw->stp[j]]->ph0.siglin + stp[shw->stp[j]]->ph1.siglin);
00936                                 shw->ph.sigmap -= (shw->stpph0sigmap[j] + shw->stpph1sigmap[j]);
00937                                 shw->ph.mip -= (shw->stpph0mip[j] + shw->stpph1mip[j]);
00938                                 shw->ph.gev -= (shw->stpph0gev[j] + shw->stpph1gev[j]);
00939         
00940                         }
00941 
00942                         
00943                 }
00944 
00945                 MSG("PETrimmer",Msg::kDebug)<<"Shower us "<<us << " vs "<< vs <<endl;
00946                 if(us ==0 || vs==0)goodshower=0; //require && here as a test
00947 
00948                 MSG("PETrimmer",Msg::kDebug)<<"Cut "<<cutcount<<" strips from shower!"<<endl;
00949 
00950 
00951                 if (cutcount==0)return goodshower;  //we didn't cut anything, so continue       
00952 
00953                 //rewrite event stp vectors
00954                 for(int j=0;j<shw->nstpcnt;j++)
00955                 {
00956                         shw->stp[j]=-1;
00957                         shw->stpph0sigmap[j]=0;
00958                         shw->stpph1sigmap[j]=0;
00959                         shw->stpph0gev[j]=0;
00960                         shw->stpph1gev[j]=0;
00961                         shw->stpph0mip[j]=0;
00962                         shw->stpph1mip[j]=0;
00963                 }
00964 
00965 
00966 /*              delete[] shw->stp;
00967                 delete[] shw->stpph0sigmap;
00968                 delete[] shw->stpph1sigmap;
00969                 delete[] shw->stpph0gev;
00970                 delete[] shw->stpph1gev;
00971                 delete[] shw->stpph0mip;
00972                 delete[] shw->stpph1mip;
00973                 
00974                 shw->stp = new Int_t[stpid.size()];
00975                 shw->stpph0sigmap = new Float_t[stpid.size()];
00976                 shw->stpph1sigmap = new Float_t[stpid.size()];
00977                 shw->stpph0gev = new Float_t[stpid.size()];
00978                 shw->stpph1gev = new Float_t[stpid.size()];
00979                 shw->stpph0mip = new Float_t[stpid.size()];
00980                 shw->stpph1mip = new Float_t[stpid.size()];     //these arrays dont seem to be in the ntpst file!
00981                 
00982 */              
00983                 
00984                 
00985                 shw->nstpcnt=stpid.size();
00986                 shw->nstrip=stpid.size();  
00987 
00988                 for(int j=0;j<(int)stpid.size();j++)
00989                 {
00990                         shw->stp[j]=-1;
00991                         if(stpid[j]>-1 && stpid[j] < (int)stp.size())   
00992                                 shw->stp[j]=stpid[j];
00993                         shw->stpph0sigmap[j]=stpph0sigmap[j];
00994                         shw->stpph1sigmap[j]=stpph1sigmap[j];
00995                         shw->stpph0gev[j]=stpph0gev[j];
00996                         shw->stpph1gev[j]=stpph1gev[j];
00997                         shw->stpph0mip[j]=stpph0mip[j];
00998                         shw->stpph1mip[j]=stpph1mip[j];
00999                 }
01000 
01001         return goodshower;
01002 
01003 }

int PETrimmer::CleanTrack ( NtpSRTrack evt,
std::vector< const NtpSRStrip * >  stp,
int *  pepass,
float *  stripe 
)

incorrect... we did not change array size.. workaround for other mods...

Definition at line 1009 of file PETrimmer.cxx.

References NtpSRStripPulseHeight::gev, Msg::kDebug, NtpSRStripPulseHeight::mip, MSG, NtpSRTrack::nstpcnt, NtpSRTrack::nstrip, NtpSRPulseHeight::pe, NtpSRTrack::ph, NtpSRPulseHeight::raw, NtpSRPulseHeight::sigcor, NtpSRPulseHeight::siglin, NtpSRStripPulseHeight::sigmap, NtpSRTrack::stp, NtpSRTrack::stpph0gev, NtpSRTrack::stpph0mip, NtpSRTrack::stpph0sigmap, NtpSRTrack::stpph1gev, NtpSRTrack::stpph1mip, and NtpSRTrack::stpph1sigmap.

Referenced by TrimRecord().

01010 {
01011                 int cutcount =0;
01012 int goodtrack=1;        
01013                 //temporary vectors for saving good hits - better than erasing single entries in arrays at a time!
01014                 std::vector<float> stpph0sigmap;
01015                 std::vector<float> stpph0mip;
01016                 std::vector<float> stpph0gev;
01017                 std::vector<float> stpph1sigmap;
01018                 std::vector<float> stpph1mip;
01019                 std::vector<float> stpph1gev;
01020                 std::vector<int> stpid;
01021 
01022 
01023                           int us=0;
01024                         int vs=0;
01025                                 
01026                 for(int j=0;j<trk->nstpcnt;j++)
01027                 {
01028                 
01029 
01030                         if(trk->stp[j]>0)       
01031                         if(pepass[trk->stp[j]])
01032                         {
01033                         
01034                                 //keep it!
01035                                 stpid.push_back(trk->stp[j]);
01036                                 stpph0sigmap.push_back(trk->stpph0sigmap[j]);
01037                                 stpph1sigmap.push_back(trk->stpph1sigmap[j]);
01038                                 stpph0gev.push_back(trk->stpph0gev[j]);
01039                                 stpph1gev.push_back(trk->stpph1gev[j]);
01040                                 stpph0mip.push_back(trk->stpph0mip[j]);
01041                                 stpph1mip.push_back(trk->stpph1mip[j]);
01042                                 
01043                                 
01044                                 stripe[trk->stp[j]*3]= trk->stpph0sigmap[j]+trk->stpph1sigmap[j];
01045                                 stripe[1+trk->stp[j]*3]= trk->stpph0mip[j]+trk->stpph1mip[j];
01046                                 stripe[2+trk->stp[j]*3]= trk->stpph0gev[j]+trk->stpph1gev[j];
01047 
01048 
01049                                 if(stp[trk->stp[j]]->planeview==2)
01050                                         us++;
01051                                 if(stp[trk->stp[j]]->planeview==3)
01052                                         vs++;
01053 
01054 
01055                                 
01056 
01057                         }else{
01058                                 cutcount++;     
01059 
01060                                 //adjust total energy of track // we dont use these values of the track (usually...)
01061                                 trk->ph.raw -= (stp[trk->stp[j]]->ph0.raw + stp[trk->stp[j]]->ph1.raw);
01062                                 trk->ph.pe -= (stp[trk->stp[j]]->ph0.pe + stp[trk->stp[j]]->ph1.pe);
01063                                 trk->ph.sigcor -= (stp[trk->stp[j]]->ph0.sigcor + stp[trk->stp[j]]->ph1.sigcor);
01064                                 trk->ph.siglin -= (stp[trk->stp[j]]->ph0.siglin + stp[trk->stp[j]]->ph1.siglin);
01065                                 trk->ph.sigmap -= (trk->stpph0sigmap[j] + trk->stpph1sigmap[j]);
01066                                 trk->ph.mip -= (trk->stpph0mip[j] + trk->stpph1mip[j]);
01067                                 trk->ph.gev -= (trk->stpph0gev[j] + trk->stpph1gev[j]);
01068 
01069 
01070                         }
01071                 }
01072 
01073                 MSG("PETrimmer",Msg::kDebug)<<"Track us "<<us << " vs "<< vs <<endl;
01074                 if(us ==0 || vs==0)goodtrack=0;
01075            
01076                 MSG("PETrimmer",Msg::kDebug)<<"Cut "<<cutcount<<" strips from track!"<<endl;
01077                 if(!goodtrack)MSG("PETrimmer",Msg::kDebug)<<"Track marked as bad!"<<endl;
01078 
01079                 if (cutcount==0 )return goodtrack;  //we didn't cut anything, so continue       
01080                 
01081                 //rewrite event stp vectors
01082                 for(int j=0;j<trk->nstpcnt;j++)
01083                 {
01084                         trk->stp[j]=-1;
01085                         trk->stpph0sigmap[j]=0;
01086                         trk->stpph1sigmap[j]=0;
01087                         trk->stpph0gev[j]=0;
01088                         trk->stpph1gev[j]=0;
01089                         trk->stpph0mip[j]=0;
01090                         trk->stpph1mip[j]=0;
01091                 }
01092 
01093 /*
01094                 delete[] trk->stp;
01095                 delete[] trk->stpph0sigmap;
01096                 delete[] trk->stpph1sigmap;
01097                 delete[] trk->stpph0gev;
01098                 delete[] trk->stpph1gev;
01099                 delete[] trk->stpph0mip;
01100                 delete[] trk->stpph1mip;
01101                 
01102                 trk->stp = new Int_t[stpid.size()];
01103                 trk->stpph0sigmap = new Float_t[stpid.size()];
01104                 trk->stpph1sigmap = new Float_t[stpid.size()];
01105                 trk->stpph0gev = new Float_t[stpid.size()];
01106                 trk->stpph1gev = new Float_t[stpid.size()];
01107                 trk->stpph0mip = new Float_t[stpid.size()];
01108                 trk->stpph1mip = new Float_t[stpid.size()];     //these arrays dont seem to be in the ntpst file!
01109 */              
01110                 
01111                 
01112                 
01113                 trk->nstpcnt=stpid.size();
01114                 trk->nstrip=stpid.size();  
01115 
01116                 for(int j=0;j<(int)stpid.size();j++)
01117                 {
01118                         trk->stp[j]=-1;
01119                         if(stpid[j]>-1 && stpid[j] < (int)stp.size())           
01120                                 trk->stp[j]=stpid[j];
01121                         trk->stpph0sigmap[j]=stpph0sigmap[j];
01122                         trk->stpph1sigmap[j]=stpph1sigmap[j];
01123                         trk->stpph0gev[j]=stpph0gev[j];
01124                         trk->stpph1gev[j]=stpph1gev[j];
01125                         trk->stpph0mip[j]=stpph0mip[j];
01126                         trk->stpph1mip[j]=stpph1mip[j];
01127                 }
01128 
01129 
01130 
01131         return goodtrack; //track is bad only if it is empty
01132 
01133 
01134 }

void PETrimmer::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Reimplemented from JobCModule.

Definition at line 1166 of file PETrimmer.cxx.

References Registry::Get(), JobCModule::GetConfig(), Msg::kDebug, Registry::LockValues(), MSG, pecut, Registry::UnLockValues(), updateEventEnergy, and Registry::ValuesLocked().

01167 {
01168 //======================================================================
01169 // Configure the module given the Registry r
01170 //======================================================================
01171   MSG("PETrimmer",Msg::kDebug)<<"In Trimmer::Config"<<endl;
01172                                                                                
01173   bool islocked = this->GetConfig().ValuesLocked();
01174   if (islocked) this->GetConfig().UnLockValues();
01175  
01176  
01177   double fmps;
01178   if(r.Get("PECut", fmps)) {pecut = fmps;}
01179 
01180   int imps;
01181   if(r.Get("updateEventEnergy", imps)) {updateEventEnergy = imps;}
01182   if (islocked) this->GetConfig().LockValues();
01183 
01184 }

const Registry & PETrimmer::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Reimplemented from JobCModule.

Definition at line 1143 of file PETrimmer.cxx.

References JobCModule::GetName(), Msg::kDebug, Registry::LockValues(), MSG, Registry::Set(), and Registry::UnLockValues().

01144 {
01145 //======================================================================
01146 // Supply the default configuration for the module
01147 //======================================================================
01148    MSG("PETrimmer",Msg::kDebug)<<"In Trimmer::DefaultConfig"<<endl;
01149 
01150   
01151   static Registry r;
01152   std::string name = this->JobCModule::GetName();
01153   name += ".config.default";
01154   r.SetName(name.c_str());
01155 
01156   r.UnLockValues();
01157   r.Set("PECut", 0.0);
01158   r.Set("updateEventEnergy",1);               
01159 
01160   r.LockValues();             
01161                
01162                                                                               
01163   return r;
01164 }

void PETrimmer::DumpIt ( NtpStRecord str  ) 

Definition at line 1188 of file PETrimmer.cxx.

References NtpStRecord::evt, NtpStRecord::GetStrips(), Msg::kDebug, Msg::kError, MSG, NtpStRecord::shw, tc, and NtpStRecord::trk.

Referenced by TrimRecord().

01189 {
01190         std::vector< NtpSREvent* > evts;
01191         std::vector< NtpSRShower* > shw;
01192         std::vector< NtpSRTrack* > trk;
01193 
01194         TClonesArray* tc = str->evt;    
01195         for(int i=0;i<tc->GetEntries();i++)
01196                 evts.push_back((NtpSREvent*)tc->At(i));
01197                 
01198         TClonesArray* ts = str->shw;    
01199         for(int i=0;i<ts->GetEntries();i++)
01200                 shw.push_back((NtpSRShower*)ts->At(i));
01201                 
01202         TClonesArray* tt = str->trk;    
01203         for(int i=0;i<tt->GetEntries();i++)
01204                 trk.push_back((NtpSRTrack*)tt->At(i));
01205 
01206         std::vector<const NtpSRStrip* > stp = str->GetStrips();
01207                 
01208         
01209   //iterate over showers, print(stripindx, pe)
01210 
01211         
01212         for(int i=0;i<(int)shw.size();i++)
01213           {
01214             ostringstream os;
01215 
01216             os <<"Shower "<<i<<": ";
01217             for(int j=0;j<shw[i]->nstrip;j++)
01218               {
01219 
01220                 if(i<0 || i > (int)shw.size()-1)
01221                 {
01222                         os << "bad shower index!";
01223                         continue;
01224                 }
01225                 int strip = shw[i]->stp[j]; 
01226                 if(strip<0)
01227                 {
01228                         os<<"("<<shw[i]->stp[j]<<",NO STRIP HERE....) ";                
01229                         MSG("PETrimmer",Msg::kError)<<"Shower "<<i<<" is missing a strip at index "<<j<< " strip index "<<shw[i]->stp[j]<<endl;
01230                 }
01231                 else
01232                 os<<"("<<shw[i]->stp[j]<<","<<stp[shw[i]->stp[j]]->ph0.pe+stp[shw[i]->stp[j]]->ph1.pe<<") ";
01233               }
01234             os <<"\n";          
01235             MSG("PETrimmer",Msg::kDebug)<<os.str()<<endl;
01236 
01237         
01238           }
01239 
01240         for(int i=0;i<(int)trk.size();i++)
01241           {
01242             ostringstream os;
01243 
01244             os <<"Track "<<i<<": ";
01245             for(int j=0;j<trk[i]->nstrip;j++)
01246               {
01247                 if(i<0 || i > (int)trk.size()-1)
01248                 {
01249                         os << "bad track index!";
01250                         continue;
01251                 }       
01252                 int strip = trk[i]->stp[j];
01253                 if(strip<0)
01254                 os<<"("<<trk[i]->stp[j]<<",NO STRIP HERE...) ";
01255                 else
01256                 os<<"("<<trk[i]->stp[j]<<","<<stp[trk[i]->stp[j]]->ph0.pe+stp[trk[i]->stp[j]]->ph1.pe<<") ";
01257               }
01258             os <<"\n";          
01259             MSG("PETrimmer",Msg::kDebug)<<os.str()<<endl;
01260           }
01261 
01262 
01263         for(int i=0;i<(int)evts.size();i++)
01264           {
01265             ostringstream os;
01266 
01267             os <<"Event "<<i<<": ";
01268             for(int j=0;j<evts[i]->nstrip;j++)
01269               {
01270 
01271                 if(i<0 || i > (int)evts.size()-1)
01272                 {
01273                         os << "bad event index!";
01274                         continue;
01275                 }
01276                 int strip = evts[i]->stp[j];
01277                 if(strip<0)
01278                  os<<"("<<evts[i]->stp[j]<<",NO STRIP HERE...) ";
01279                 else
01280                 os<<"("<<evts[i]->stp[j]<<","<<stp[evts[i]->stp[j]]->ph0.pe+stp[evts[i]->stp[j]]->ph1.pe<<") ";
01281               }
01282             os <<"\n";          
01283 
01284 
01285             os<<"\tShowers: ";
01286             for(int j=0;j<evts[i]->nshower;j++)
01287               os<<evts[i]->shw[j]<<" ";
01288 
01289             os<<"\tTracks: ";
01290             for(int j=0;j<evts[i]->ntrack;j++)
01291               os<<evts[i]->trk[j]<<" ";
01292 
01293 
01294 
01295 
01296             MSG("PETrimmer",Msg::kDebug)<<os.str()<<endl;
01297           }
01298 
01299 
01300 }

void PETrimmer::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 1139 of file PETrimmer.cxx.

01140 {
01141 }

JobCResult PETrimmer::Reco ( MomNavigator mom  )  [virtual]

Implement this for read-write access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 67 of file PETrimmer.cxx.

References MomNavigator::GetFragmentList(), JobCResult::kPassed, and TrimRecord().

00068 {
00069    bool foundST=false;
00070                                                                                
00071 //   VldContext vc;
00072    //first ask mom for a NtpStRecord
00073 
00074 //iterate over all NtpStRecords in MOM
00075 
00076 std::vector<TObject*> records = mom->GetFragmentList("NtpStRecord");
00077 
00078 
00079 for(unsigned int i=0;i<records.size();i++)
00080 {
00081 
00082    NtpStRecord *str=dynamic_cast<NtpStRecord *>(records[i]); //(mom->GetFragment("NtpStRecord","Primary"));
00083    if(str){
00084      foundST=true;
00085   //   vc=str->GetHeader().GetVldContext();
00086    }else{
00087         continue;
00088    }
00089    //cout<<"Running over "<<str->GetName()<<endl;
00090 
00091    TrimRecord(str);
00092 
00093 }
00094 
00095    return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00096 
00097 }

void PETrimmer::TrimRecord ( NtpStRecord str  ) 

Definition at line 100 of file PETrimmer.cxx.

References NtpSRPlane::beg, CleanEvent(), CleanShower(), CleanTrack(), DumpIt(), MuELoss::e, NtpSRPlane::end, NtpStRecord::evt, NtpVtxFinder::FindVertex(), RecRecordImp< T >::GetHeader(), RecPhysicsHeader::GetSnarl(), NtpStRecord::GetStrips(), NtpSREvent::index, Msg::kDebug, Msg::kError, Msg::kWarning, NtpSRShowerPulseHeight::linNCgev, NtpSRStripPulseHeight::mip, MSG, NtpSREvent::nstpcnt, pecut, NtpSRShower::ph, NtpSRTrack::plane, NtpSRVertex::plane, NtpStRecord::shw, NtpSRShower::shwph, NtpSREvent::stp, NtpSRVertex::t, tc, NtpStRecord::trk, NtpSRVertex::u, NtpSRVertex::v, NtpSRShower::vtx, NtpSREvent::vtx, NtpSRTrack::vtx, NtpVtxFinder::VtxPlane(), NtpVtxFinder::VtxT(), NtpVtxFinder::VtxU(), NtpVtxFinder::VtxV(), NtpVtxFinder::VtxX(), NtpVtxFinder::VtxY(), NtpVtxFinder::VtxZ(), NtpSRVertex::x, NtpSRVertex::y, and NtpSRVertex::z.

Referenced by Reco().

00101 {
00102    DumpIt(str);
00103 
00104 
00105 
00106         MSG("PETrimmer",Msg::kDebug)<< "Cutting out strips below "<<pecut<<" pe"<<endl;
00107         
00108         std::vector<const NtpSRStrip* > stp = str->GetStrips();
00109         
00110         if(stp.size()<1)return; //no strips!
00111 
00112         std::vector< NtpSREvent* > evts;
00113         std::vector< NtpSRShower* > shw;
00114         std::vector< NtpSRTrack* > trk;
00115 
00116 
00117 MSG("PETrimmer",Msg::kDebug)<< "Snarl " << str->GetHeader().GetSnarl() << endl;
00118 
00119 
00120 
00121         TClonesArray* tc = str->evt;    
00122         for(int i=0;i<tc->GetEntries();i++)
00123                 evts.push_back((NtpSREvent*)tc->At(i));
00124                 
00125         TClonesArray* ts = str->shw;    
00126         for(int i=0;i<ts->GetEntries();i++)
00127                 shw.push_back((NtpSRShower*)ts->At(i));
00128                 
00129         TClonesArray* tt = str->trk;    
00130         for(int i=0;i<tt->GetEntries();i++)
00131                 trk.push_back((NtpSRTrack*)tt->At(i));
00132                 
00133 
00134    for(int i=0;i<(int)evts.size();i++)
00135         {       
00136           MSG("PETrimmer",Msg::kDebug)<< "Event "<<i<<" strips "<<evts[i]->nstrip<<" idx "<<evts[i]->index<<endl;
00137         };
00138 
00139 
00140 
00141 
00143         for(int i=0;i<(int)shw.size();i++)
00144           {
00145        
00146 
00147 
00148                 //found a missing strip in the mrcc far data...  hopefully this will fix it!
00149                 if(shw[i]->stp[shw[i]->nstrip-1]<0)
00150                 {
00151                         int nstrip = shw[i]->nstrip;
00152                         MSG("PETrimmer",Msg::kError)<<"Shower "<<i<<" is missing the last strip!!!! adjusting strip size from "<< nstrip <<" to "<<(nstrip-1)<<endl;
00153                         shw[i]->nstrip=shw[i]->nstrip-1;
00154                         shw[i]->nstpcnt=shw[i]->nstrip;
00155                 }
00156 
00157           }
00158 
00159         
00160 
00161 
00162 
00163 
00164 
00165 
00166 
00167 
00168 
00169 
00170 
00171 
00172 
00173 
00174 
00175 
00176 
00177 
00178 
00179 
00180 
00181 
00182 
00183 
00184 
00185 
00186 
00187 
00188 
00189 
00190 
00191 
00192 
00193       
00194 
00195                 
00196 
00197                 
00198 //      float *pepass = new float[stp.size()];  
00199 //      float *stripe_trk = new float[3 * stp.size()];
00200 //      float *stripe_shw = new float[3 * stp.size()];
00201 
00202 
00203         int pepass[100000];  //will store if that strip passes the pe cut
00204         float stripe_trk[3 * 100000]; //energy of strip as appearing in track objects - will hold sigcor, mip, gev
00205         float stripe_shw[3 * 100000]; //energy of strip as appearing in shower object
00206         
00207 
00208         //these are arrays that will store 1 if that index in shower or event is to be kept
00209         //hopfully there wont be a snarl with more than 1000 showers or events
00210         int goodshower[1000]; //if its bad (a view with no hits).. we will remove them from the event!
00211         for(int i=0;i<1000;i++)goodshower[i]=1;
00212 
00213         int goodevent[1000];
00214         for(int i=0;i<1000;i++)goodevent[i]=1;
00215 
00216         int goodtrack[1000];
00217         for(int i=0;i<1000;i++)goodtrack[i]=1;
00218 
00219 
00220         for(int i=0;i<(int)stp.size();i++)
00221         {
00222                 pepass[i]=stp[i]->ph0.pe+stp[i]->ph1.pe - pecut > 0.0001 ? 1: 0;
00223 
00224                 for(int j=0;j<3;j++)
00225                 {
00226                         stripe_trk[j + 3*i]=0;
00227                         stripe_shw[j + 3*i]=0;
00228                 }
00229         }
00230                 
00231         
00232 
00233 
00234 
00235         for(int i=0;i<(int)shw.size();i++)
00236         {
00237                 MSG("PETrimmer",Msg::kDebug)<< "Shw "<<i<<" nstrip " << shw[i]->nstrip<< " stpcnt "<<shw[i]->nstpcnt<<endl;
00238                 NtpSRShower *sh = shw[i];
00239                 goodshower[i]=CleanShower(sh,stp,pepass,stripe_shw);
00240         }       
00241 
00242         
00243         
00244         for(int i=0;i<(int)trk.size();i++)
00245         {
00246                 MSG("PETrimmer",Msg::kDebug)<< "Trk "<<i<<" nstrip " << trk[i]->nstrip<< " stpcnt "<<trk[i]->nstpcnt<<endl;
00247                 NtpSRTrack *tr = trk[i];
00248                 goodtrack[i]=CleanTrack(tr,stp,pepass,stripe_trk);
00249         }       
00250 
00251         
00252         for(int i=0;i<(int)evts.size();i++)
00253         {
00254                 MSG("PETrimmer",Msg::kDebug)<< "Event "<<i<<" nstrip " << evts[i]->nstrip<< " stpcnt "<<evts[i]->nstpcnt<<endl;
00255                 NtpSREvent* evt = evts[i];
00256                 goodevent[i]=CleanEvent(evt,shw,trk,stp,pepass,stripe_trk,stripe_shw,goodshower,goodtrack);
00257 //                 MSG("PETrimmer",Msg::kInfo)<<"Status: " << goodevent[i]<<endl;
00258 
00259 
00260                 if(!goodevent[i])continue;
00261 
00262                 NtpVtxFinder finder(i,str);
00263                 if(finder.FindVertex() >0)
00264                 {
00265                         evt->vtx.x=finder.VtxX();
00266                         evt->vtx.y=finder.VtxY();
00267                         evt->vtx.z=finder.VtxZ();
00268                         evt->vtx.u=finder.VtxU();
00269                         evt->vtx.v=finder.VtxV();
00270                         evt->vtx.t=finder.VtxT();
00271                         evt->vtx.plane=finder.VtxPlane();
00272                 }
00273 
00274                 int thisVtxPlane=evt->vtx.plane;
00275                 int invtxwindow=0;
00276                 for(int j=0;j<evt->nstpcnt;j++)
00277                 {
00278                         if(evt->stp[j]<0)continue;
00279                         if(pepass[evt->stp[j]])
00280             if(stp[evt->stp[j]]->plane >= thisVtxPlane && stp[evt->stp[j]]->plane < (thisVtxPlane+5))
00281             {
00282                 invtxwindow=1;
00283                 break;
00284             }   
00285         }              
00286 
00287         if(!invtxwindow)
00288         {
00289                         float evtpes=0;
00290             for(int j=0;j<evt->nstpcnt;j++)
00291             {
00292                 if(evt->stp[j]<0)continue;
00293                 if(pepass[evt->stp[j]])
00294                         evtpes+= (stp[evt->stp[j]]->ph0.pe + stp[evt->stp[j]]->ph1.pe);
00295             }
00296                 MSG("PETrimmer",Msg::kWarning)<<"removing "<<evtpes<<" pe event that has no strips in vtx window!"<<endl;
00297 
00298                 goodevent[i]=0;
00299         }
00300 
00301         }
00302 
00303 
00304       
00305      
00306         //remove bad events from the ntpstrecord..
00307 
00308         std::vector<NtpSREvent *>todelrec;
00309 
00310         //int tot=evts.size();
00311         int off=0;
00312         for(int i=0;i<(int)evts.size();i++)
00313         {
00314                 if(!goodevent[i+off])
00315                 {
00316                 
00317                         MSG("PETrimmer",Msg::kDebug)<< "Event removed from snarl "<<i+off<<endl;
00318         
00319                         todelrec.push_back(evts[i]);
00320                         evts.erase(evts.begin()+i);
00321                         i--;    
00322                         off++;
00323                         //str->evthdr.nevent--; //don't adjust this... we are not resizing the event array due to truth matching
00324                 }
00325 
00326         }
00327 
00328         
00329         //clean up primary tracks, showers... 
00330         //for now, if we deleted a primary, set it to -1 (don't try to find a replacement!)
00331 
00332         for(int i=0;i<(int)evts.size();i++)
00333         {
00334            if(evts[i]->primtrk > -1){
00335               bool isOK = false;
00336               int longest = 0;
00337               int longestTrackIndex = -1;
00338 
00339               for(int j=0;j<(int)evts[i]->ntrack && !isOK; j++)
00340               {
00341                  int index = evts[i]->trk[j];
00342                  if(index == -1) continue;
00343                  if(index == evts[i]->primtrk) isOK = true;
00344 
00345                  NtpSRTrack* track = trk[index];
00346                  int length =  TMath::Abs(track->plane.end - track->plane.beg);
00347                  if(length > longest){
00348                     longest = length;  longestTrackIndex = j;
00349                  }
00350               }
00351               if(!isOK) {  evts[i]->primtrk = longestTrackIndex; 
00352   //               cout<<"Assigning new primary track"<<longest<<endl;
00353               }
00354            }
00355 
00356            //check prim shw
00357            if(evts[i]->primshw>-1) {
00358               bool isOK = false;
00359 
00360               //code ruthlessly borrowed from CandEventHandle::SetPrimaryShower
00361 
00362               double biggestE = -1;
00363               double closeE = -1;
00364               int primShwIndex = -1;
00365               int largestShwIndex = -1;
00366               int closeShwIndex = -1;
00367               double dzLarge = 0;
00368               double dzClose = 1e6;
00369               NtpSRTrack *primaryTrack = 0;
00370               if(evts[i]->primtrk != -1) primaryTrack = trk[evts[i]->trk[evts[i]->primtrk]];
00371 
00372               for(int j=0;j<(int)evts[i]->nshower && !isOK;j++)
00373               {
00374                  int index = evts[i]->shw[j];
00375                  if(index == -1) continue;
00376                  if(index==evts[i]->primshw) isOK = true;
00377                  
00378                  NtpSRShower* shower = shw[index];
00379                  if(shower->ph.mip > biggestE){
00380                     largestShwIndex = j;
00381                     biggestE = shower->ph.mip;
00382                     if(primaryTrack)  dzLarge = fabs(shower->vtx.z - primaryTrack->vtx.z); 
00383                  }
00384                  if(primaryTrack){
00385                      if(fabs(shower->vtx.z - primaryTrack->vtx.z) < dzClose){
00386                          dzClose = fabs(shower->vtx.z - primaryTrack->vtx.z);
00387                          closeShwIndex = j;
00388                          closeE = shower->ph.mip;
00389                       }
00390                  }
00391               }
00392               if(!isOK){  //the primary shower has died, we must elect a new leader
00393                 float eratio = 1;
00394                 if(biggestE > 0) eratio = closeE/biggestE;
00395 
00396 //                cout<<primShwIndex<<"  "<<eratio<<"  "<<dzLarge<<"  "<<dzClose<<endl;
00397                 if(closeShwIndex > -1 && (dzLarge - dzClose) > 1 && eratio > 0.2)
00398                    primShwIndex = closeShwIndex;
00399                 else
00400                    primShwIndex = largestShwIndex;
00401 
00402                 if(primShwIndex > -1 && primaryTrack){
00403                   NtpSRShower *PShw = shw[evts[i]->shw[primShwIndex]];
00404                   if(PShw == 0) cout<<"Failed to load shower - VERY BAD"<<endl;
00405                   if(fabs(PShw->vtx.z - primaryTrack->vtx.z) > 0.5 &&
00406                           PShw->shwph.linNCgev < 2.0){
00407                     primShwIndex = -1;
00408                   }
00409                 }
00410                 evts[i]->primshw = primShwIndex;
00411 //                cout<<"New primary shower assigned!"<<primShwIndex<<"  "<<largestShwIndex
00412 //                    <<"  "<<closeShwIndex<<"  "<<endl;
00413               }
00414 //              cout<<"Event Primary shw: "<<evts[i]->primshw<<"  "<<evts[i]->nshower<<endl;                  
00415            } 
00416         }
00417 
00418         //I am not directly delete from the tclonesarray, as it seems to cause problems....
00419         //instead, clear the tclones array, and put the pointers for the good events back in
00420 //check that events make sense....
00421 //str->evt->Clear();
00422 
00423 
00424 for(int i=0;i<(int)todelrec.size();i++)
00425         str->evt->Remove(todelrec[i]);  //does not shift entires...
00426 
00427 
00429 //.......   str->evt->Compress();
00430 
00431 
00432 
00433 
00434 /*
00435    for(int i=0;i<(int)str->evt->GetEntries();i++)
00436      {
00437        ( (NtpSREvent*)str->evt->At(i))->index=i;
00438      }
00439 */
00440    //check that we saved them
00441    evts.clear();
00442    for(int i=0;i<tc->GetEntries();i++)
00443         {
00444                 NtpSREvent * e = ((NtpSREvent*)str->evt->At(i));
00445                 if(!e)continue;
00446                 e->index=i;
00447                 evts.push_back((NtpSREvent*)str->evt->At(i));
00448         }
00449 
00450 
00451    for(int i=0;i<(int)evts.size();i++)
00452    {      
00453           MSG("PETrimmer",Msg::kDebug)<< "Event "<<i<<" strips "<<evts[i]->nstrip<<" idx "<<evts[i]->index<<endl;
00454    };
00455 
00456 
00457       
00458 
00459     
00460    
00461   
00462  
00463 
00464 
00465 
00466 
00467 
00468 
00469 
00470 
00471 
00472 //for debuggins...
00473 
00474 
00475         for(int i=0;i<(int)evts.size();i++)
00476         {
00477         std::ostringstream s;
00478                 s<<" showers ";
00479                 for(int j=0;j<evts[i]->nshower;j++)
00480                         s<<" " << evts[i]->shw[j] <<" ";
00481                 s<<" tracks ";
00482                 for(int j=0;j<evts[i]->ntrack;j++)
00483                         s<<" " << evts[i]->trk[j] <<" ";
00484 
00485 
00486         MSG("PETrimmer",Msg::kDebug)<< "Event done    strips "<<evts[i]->nstrip<< s.str() << endl;
00487 
00488         }
00489         
00490 //      delete[] pepass;
00491 //      delete[] stripe_trk;
00492 //      delete[] stripe_shw;
00493 
00494 
00495 //      DumpIt(str);
00496 
00497    
00498 
00499 }


Member Data Documentation

float PETrimmer::pecut [private]

Definition at line 53 of file PETrimmer.h.

Referenced by Config(), and TrimRecord().

Definition at line 54 of file PETrimmer.h.

Referenced by CleanEvent(), and Config().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1