DetailedParticle Class Reference

#include <DetailedParticle.h>

List of all members.

Public Member Functions

 DetailedParticle (const char *fname)
 ~DetailedParticle ()
void Run ()
int PassesPreselec ()

Public Attributes

TChain * input
ParticleObjectHolderpoh
TChain * inputPA
PRecordprec

Detailed Description

Definition at line 10 of file DetailedParticle.h.


Constructor & Destructor Documentation

DetailedParticle::DetailedParticle ( const char *  fname  ) 

Definition at line 14 of file DetailedParticle.cxx.

References input, inputPA, poh, and prec.

00015 {
00016         input=new TChain("PO");
00017         inputPA=new TChain("PA");
00018         input->Add(fname);
00019         inputPA->Add(fname);
00020         poh=new ParticleObjectHolder();
00021         prec=new PRecord();
00022         input->SetBranchAddress("ParticleObjectHolder",&poh);
00023         inputPA->SetBranchAddress("PRecord",&prec);
00024 }

DetailedParticle::~DetailedParticle (  ) 

Definition at line 26 of file DetailedParticle.cxx.

00027 {}


Member Function Documentation

int DetailedParticle::PassesPreselec (  ) 

Definition at line 30 of file DetailedParticle.cxx.

References Particles::elec_vise, Particles::emfrac, PRecord::event, Particles::largest_particle_e, Particles::longest_s_particle_s, Event::nstrips, PRecord::particles, prec, Particles::primary_long_e, Particles::totvise, Event::vtx_u, Event::vtx_v, and Event::vtx_z.

Referenced by Run().

00031 {
00032         int pass=1;
00033         
00034 //      pass = pass && prec->event.inFiducial==1;
00035 
00036 
00037         pass = pass && ( (prec->event.vtx_z>0.5 && prec->event.vtx_z<14.5) || (prec->event.vtx_z>16.5 && prec->event.vtx_z<29.5) );
00038 
00039         double r = sqrt(prec->event.vtx_u*prec->event.vtx_u+prec->event.vtx_v*prec->event.vtx_v);
00040         
00041         pass = pass && (r>0.25 && r<3.9);
00042         
00043         pass = pass && prec->event.nstrips>10;
00044         pass = pass && prec->event.nstrips<50;
00045         pass = pass && prec->particles.longest_s_particle_s<0.9;
00046         pass = pass && prec->particles.emfrac>0.4;
00047         pass = pass && prec->particles.elec_vise/25.>1;
00048         pass = pass && prec->particles.primary_long_e<40;
00049         pass = pass && prec->particles.largest_particle_e/prec->particles.totvise>0.65;
00050 /*      pass = pass && prec->event.visenergy/25.-prec->particles.totvise/25>0;
00051         pass = pass && prec->event.visenergy/25.-prec->particles.totvise/25<0.5;
00052 */      
00053         
00054         /*
00055 pass = pass && prec->particles.ntot>0;
00056 pass = pass && prec->particles.longest_s_particle_s<1.1;
00057 pass = pass && prec->particles.elec_vise<200;
00058 pass = pass && prec->particles.nelec>0;
00059 pass = pass && prec->particles.largest_particle_type==11;
00060 pass = pass && prec->event.nstrips>10;
00061 pass = pass && prec->event.nstrips<55;
00062 pass = pass && prec->particles.largest_particle_avg_rms>0.004;
00063 */
00064 
00065 /*
00066         pass = pass && prec->event.nstrips>10 && prec->event.nstrips<50;        
00067         pass = pass && prec->particles.ntot>0;
00068         pass = pass &&  prec->event.vtx_z-prec->event.min_z<0.35;
00069         pass = pass &&  prec->event.nstrips>15 && prec->event.nstrips<45;
00070         pass = pass &&  prec->particles.rms_r<1.1;
00071         pass = pass &&  prec->particles.prim_par_b>0.1&&prec->particles.prim_par_b<1.2;
00072         pass = pass &&  prec->particles.ntot/prec->event.visenergy<0.14;
00073         pass = pass &&  prec->particles.mol_rad_r<0.0425;
00074         pass = pass &&  prec->particles.largest_particle_avg_rms<0.017;
00075         pass = pass &&  prec->particles.largest_particle_avg_rms>0.005;
00076         pass = pass &&  prec->particles.ntot<8 ;
00077         pass = pass &&  prec->particles.longest_z<0.8;
00078         pass = pass &&  prec->particles.primary_long_e<60;
00079         
00080 */      
00081 
00082         
00083 /*      
00084         pass = pass && prec->particles.rough_primary_theta_z<0.7;
00085         pass = pass && prec->particles.primary_long_e<15;
00086         
00087         pass = pass && sqrt((prec->event.max_u-prec->event.min_u)*(prec->event.max_u-prec->event.min_u)
00088                         +(prec->event.max_v-prec->event.min_v)*(prec->event.max_v-prec->event.min_v))>0.3;
00089                         
00090         pass = pass && prec->particles.largest_particle_avg_rms>0.0045;
00091         pass = pass && prec->particles.largest_particle_e/prec->event.visenergy>0.85;
00092 
00093         pass = pass && prec->particles.total_long_e<25;
00094         pass = pass && prec->particles.mol_rad_r>0.001;
00095         pass = pass && prec->particles.mol_rad_r<0.035;
00096 
00097         pass = pass && prec->particles.nelec>0;
00098         pass = pass && prec->particles.ntot<4;
00099         pass = pass && prec->particles.prim_par_b>0.1;
00100         pass = pass && prec->particles.prim_par_b<1;
00101 
00102         pass = pass && prec->particles.longest_z<1;
00103         pass = pass && prec->particles.longest_z>0.1;
00104         pass = pass && prec->particles.longest_particle_type<14;
00105 */
00106         return pass;
00107 }

void DetailedParticle::Run (  ) 

hist weighted area ///////

Definition at line 109 of file DetailedParticle.cxx.

References MuELoss::e, Particle3D::e, Particle3D::end_z, Particle3D::entries, ParticleObjectHolder::event, RecCandHeader::GetEvent(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), input, inputPA, ParticleObjectHolder::mctrue, PRecord::mctrue, n, MCTrue::oscprob, ParticleObjectHolder::particles3d, PassesPreselec(), poh, prec, Particle3D::start_z, Particle3D::sum_e, MCTrue::totbeamweight, MCTrue::type, Particle3D::u, Particle3D::v, Particle3D::view, ParticleEvent::vtx_u, ParticleEvent::vtx_v, ParticleEvent::vtx_z, and Particle3D::z.

00110 {
00111 
00112         
00113         TH2F * hist_sumdist[5];
00114         TH2F * hist_weighted_sumdist[5]; 
00115         TH1F * hist_max_dist[5]; 
00116         TH1F * hist_max_dist_per_gev[5];        
00117 
00118 
00119         TH1F * hist_weighted_area[5];   
00120 
00121 
00122         TH1F * hist_first_peak[5];
00123         TH1F * hist_second_peak[5];
00124         TH1F * hist_twopeak_diff[5];
00125         
00126         TH1F * hist_back_peak_dist[5];
00127         
00128         TH1F * hist_view_e_asym[5];             
00129 
00130 
00131         TH2F * hist_prim_angle_vs_prim_eres[5];         
00132         TH2F * hist_prim_angle_vs_prim_true[5];         
00133         
00134         TH2F * hist_sum_curv_vs_part_e[5];
00135         TH2F * hist_pu_pz[5];
00136         TH2F * hist_pv_pz[5];
00137         TH2F * hist_pt_pz[5];
00138         TH2F * hist_pu_pv[5];
00139         
00140         for(int i=0;i<5;i++)
00141         {
00142                 char tmp[200];
00143                 
00144                 sprintf(tmp,"sumdist_%d",i);
00145                 hist_sumdist[i] = new TH2F(tmp,tmp,100,0, 50,20,0,20);
00146 
00147                 sprintf(tmp,"weighted_sumdist_%d",i);           
00148                 hist_weighted_sumdist[i] = new TH2F(tmp,tmp,100,0, 50,20,0,20);
00149 
00150                 sprintf(tmp,"max_dist_%d",i);                   
00151                 hist_max_dist[i] = new TH1F(tmp,tmp,100,0,20);
00152                 
00153                 sprintf(tmp,"max_dist_per_gev_%d",i);           
00154                 hist_max_dist_per_gev[i] = new TH1F(tmp,tmp,100,0,1);
00155         
00156                 sprintf(tmp,"weighted_area_%d",i);              
00157                 hist_weighted_area[i] = new TH1F(tmp,tmp,100,0,20);
00158 
00159                 sprintf(tmp,"first_peak_%d",i);         
00160                 hist_first_peak[i] = new TH1F(tmp,tmp,100,0,2);
00161 
00162                 sprintf(tmp,"second_peak_%d",i);                
00163                 hist_second_peak[i] = new TH1F(tmp,tmp,100,0,2);
00164                 
00165                 sprintf(tmp,"twopeak_diff_peak_%d",i);          
00166                 hist_twopeak_diff[i] = new TH1F(tmp,tmp,100,0,2);
00167                                 
00168                 sprintf(tmp,"back_peak_dist_%d",i);             
00169                 hist_back_peak_dist[i] = new TH1F(tmp,tmp,100,0,2);     
00170                 
00171                 sprintf(tmp,"view_e_asym_%d",i);                
00172                 hist_view_e_asym[i] = new TH1F(tmp,tmp,100,0,1);        
00173 
00174                 sprintf(tmp,"prim_angle_vs_prim_eres_%d",i);
00175                 hist_prim_angle_vs_prim_eres[i] = new TH2F(tmp,tmp,100,0, 3,100,-1,1);          
00176                 sprintf(tmp,"prim_angle_vs_prim_true_%d",i);
00177                 hist_prim_angle_vs_prim_true[i] = new TH2F(tmp,tmp,100,0, 3,100,0,10);  
00178 
00179 
00180                 sprintf(tmp,"sum_curv_vs_part_e_%d",i);
00181                 hist_sum_curv_vs_part_e[i]=new TH2F(tmp,tmp,1000,0,50,100,0,10);
00182                 sprintf(tmp,"pu_pz_%d",i);
00183                 hist_pu_pz[i]=new TH2F(tmp,tmp,100,-10,10,100,0,10);
00184                 sprintf(tmp,"pv_pz_%d",i);
00185                 hist_pv_pz[i]=new TH2F(tmp,tmp,100,-10,10,100,0,10);
00186                 sprintf(tmp,"pt_pz_%d",i);
00187                 hist_pt_pz[i]=new TH2F(tmp,tmp,100,0,10,100,0,10);
00188 
00189                 sprintf(tmp,"pu_pv_%d",i);
00190                 hist_pu_pv[i]=new TH2F(tmp,tmp,100,-10,10,100,-10,10);
00191                                 
00192         }
00193         
00194 
00195 
00196 
00197 
00198         int ent =input->GetEntries();
00199         int entpa =inputPA->GetEntries();
00200         printf("%d %d total entries\n",ent,entpa);
00201         
00202         
00203         double typecnt[5];
00204         for(int i=0;i<5;i++)typecnt[i]=0;
00205         
00206         
00207         int saved=0;
00208         int ji=0;
00209         input->GetEntry(0);
00210         for(int ent_idx=0;ent_idx<entpa;ent_idx++)
00211         {
00212                 //input->GetEntry(i);
00213                 /*if(i==0)*/inputPA->GetEntry(ent_idx);
00214                 
00215                 while(poh->GetHeader().GetRun() != prec->GetHeader().GetRun() || poh->GetHeader().GetSnarl()!=prec->GetHeader().GetSnarl()  || poh->GetHeader().GetEvent()!=prec->GetHeader().GetEvent())
00216                         input->GetEntry(ji++);
00217                         
00218 /*
00219                 int j=0;
00220                 while(poh->GetHeader().GetRun() != prec->GetHeader().GetRun())
00221                 {
00222                         j++;
00223                         inputPA->GetEntry(i+j);
00224         //                                      printf("sync a %d %d\n",poh->GetHeader().GetSnarl(),prec->GetHeader().GetSnarl());
00225                 }
00226                 
00227 
00228                 while(poh->GetHeader().GetSnarl()<prec->GetHeader().GetSnarl() && poh->GetHeader().GetRun()==prec->GetHeader().GetRun())
00229                 {
00230                         i++;
00231                         input->GetEntry(i);
00232         //                              printf("sync b %d %d\n",poh->GetHeader().GetSnarl(),prec->GetHeader().GetSnarl());
00233                 }
00234                 
00235                 j=0;
00236                 while(poh->GetHeader().GetSnarl()>prec->GetHeader().GetSnarl() && poh->GetHeader().GetRun()==prec->GetHeader().GetRun())
00237                 {
00238                         j++;
00239                         inputPA->GetEntry(i+j);
00240         //                              printf("sync c %d %d\n",poh->GetHeader().GetSnarl(),prec->GetHeader().GetSnarl());
00241                 }
00242                 
00243                 if(poh->GetHeader().GetRun()!=prec->GetHeader().GetRun())continue;
00244                 */
00245                 saved++;
00246                 
00247 /*
00248                 while(poh->GetHeader().GetSnarl()>prec->GetHeader().GetSnarl() && poh->GetHeader().GetRun()==prec->GetHeader().GetRun())
00249                 {
00250                         printf("sync %d %d\n",poh->GetHeader().GetSnarl(),prec->GetHeader().GetSnarl());
00251                         j++;
00252                         if(i+j>inputPA->GetEntries()-1)break;
00253                         
00254                         inputPA->GetEntry(i+j);
00255                 }
00256                 
00257                 
00258                 if(i+j>inputPA->GetEntries()-1)break;
00259                 
00260                 if(poh->GetHeader().GetSnarl()<prec->GetHeader().GetSnarl() || poh->GetHeader().GetRun()!=prec->GetHeader().GetRun())
00261                 {
00262                         printf("sync %d %d\n",poh->GetHeader().GetSnarl(),prec->GetHeader().GetSnarl());
00263 
00264                         continue;
00265                 }
00266 */              
00267                 
00268                 if(!PassesPreselec())continue;
00269         
00270         
00271         
00272                 double weight = poh->mctrue.totbeamweight*poh->mctrue.oscprob*3.25/6.5/5;
00273                 
00274                 double paweight = prec->mctrue.totbeamweight*prec->mctrue.oscprob*3.25/6.5/5;
00275                 
00276                 if(fabs(weight-paweight)>0.001)printf("weight mismatch!\n");
00277                 
00278                 int type=poh->mctrue.type;
00279         
00280                 if(type<0 || type>5)continue;
00281                 
00282                 if(type==0)weight=weight/3.;
00283         
00284         
00285         
00286         //      MsgService * ms = MsgService::Instance();
00287         //      ms->GetStream("TruthCompareAna")->SetLogLevel(Msg::kDebug);
00288         //      TruthCompareAna a;
00289         //      a.ana(poh,&prec->truthcompare);
00290         
00291         /*
00293                 printf("\nstdhep....run %d snarl %d event %d\n",prec->GetHeader().GetRun() , prec->GetHeader().GetSnarl() , prec->GetHeader().GetEvent());
00294 
00295                 double true_elec_p4[4];//largest true election
00296                 true_elec_p4[3]=0;
00297                 int true_elec_idx=-1;
00298                 
00299                 for(unsigned int i=0;i<prec->mctrue.stdhep.size();i++)
00300                 {
00301                         NtpMCStdHep *mc = &(prec->mctrue.stdhep[i]);
00302                         
00303                         if(mc->IstHEP!=1)continue;
00304                         
00305                         printf("found %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]);
00306                         
00307                         //need lepton energy from cc
00308                         //need sum of pi0 energy
00309                         //number of pi0
00310                         
00311                         //number of pip,pim
00312                         
00313                         if(fabs(mc->IdHEP)==11 && true_elec_p4[3]<mc->p4[3])
00314                         {
00315                                 true_elec_p4[0]=mc->p4[0];
00316                                 true_elec_p4[1]=mc->p4[1];
00317                                 true_elec_p4[2]=mc->p4[2];
00318                                 true_elec_p4[3]=mc->p4[3];
00319                                 true_elec_idx=i;
00320                         }
00321                 
00322                 }
00323                 printf("\n");
00324                 
00325                 
00327         */
00328         
00329         
00330                 //calculate view momentums of found particles relative to reco vertex
00331                 
00332                 double pu=0;
00333                 double pv=0;
00334                 double pz_u=0;
00335                 double pz_v=0;
00336                 double pz=0;
00337                 double pt=0;                    
00338                 for(int j=0;j<(int)poh->particles3d.size();j++)
00339                 {
00340 
00341                         Particle3D * p = &poh->particles3d[j];
00342                         if(!p)continue;
00343                         
00344                         for(int i=0;i<p->entries;i++)
00345                         {
00346                                 if(p->view[i]==2)
00347                                 {
00348                                         double dt=p->u[i]-poh->event.vtx_u;
00349                                         double dz=p->z[i]-poh->event.vtx_z;
00350                                         double r = sqrt(dt*dt+dz*dz);
00351                                         if(r>0)pu+= p->e[i]/25.*dt/r;
00352                                         if(r>0)pz_u+= p->e[i]/25.*dz/r;
00353                                 }
00354 
00355                                 if(p->view[i]==2)
00356                                 {
00357                                         double dt=p->v[i]-poh->event.vtx_v;
00358                                         double dz=p->z[i]-poh->event.vtx_z;
00359                                         double r = sqrt(dt*dt+dz*dz);
00360                                         if(r>0)pv+= p->e[i]/25.*dt/r;
00361                                         if(r>0)pz_v+= p->e[i]/25.*dz/r;
00362                                 }                       
00363                         
00364                                 double du=p->u[i]-poh->event.vtx_u;
00365                                 double dv=p->v[i]-poh->event.vtx_v;
00366                                 double dz=p->z[i]-poh->event.vtx_z;
00367                                 double r = sqrt(du*du+dv*dv+dz*dz);
00368                                 double dt=sqrt(du*du+dv*dv);
00369                                 if(r>0)pt+=p->e[i]/25.*dt/r;
00370                                 if(r>0)pz+=p->e[i]/25.*dz/r;
00371                         
00372                         }
00373                         
00374                 }
00375                 
00376                 hist_pu_pz[type]->Fill(pu,pz_u,weight);
00377                 hist_pv_pz[type]->Fill(pv,pz_v,weight);
00378                 hist_pt_pz[type]->Fill(pt,pz,weight);
00379                 hist_pu_pv[type]->Fill(pu,pv,weight);
00380                 
00381                 
00382                 int longest_id=-1;
00383                 double longZ=0;
00384                 
00385                 int largest_id=-1;
00386                 double largeE=0;
00387                                 
00388                 for(int j=0;j<(int)poh->particles3d.size();j++)
00389                 {
00390 
00391                         Particle3D * p = &poh->particles3d[j];
00392                         if(!p)continue;
00393                         
00394 
00395                         
00396                         if(p->sum_e>largeE)
00397                         {
00398                                 largest_id=j;
00399                                 largeE=p->sum_e;
00400                         }
00401                         
00402                         if(p->end_z-p->start_z>longZ)
00403                         {
00404                                 longest_id=j;
00405                                 longZ=p->end_z-p->start_z;
00406                         }
00407                         
00408 
00410                         double sum_z2=0;
00411                         double sum_z=0;
00412                         double sum_u=0;
00413                         double sum_zu=0;                        
00414                         double sum_v=0;
00415                         double sum_zv=0;
00416                         int n=0;        
00417                                                 
00418                         for(int k=0;k<p->entries;k++)
00419                         {
00420                                 sum_z2+=p->z[k]*p->z[k];
00421                                 sum_z+=p->z[k];
00422                                 sum_zu+=p->z[k]*p->u[k];
00423                                 sum_u+=p->u[k];
00424                                 sum_zv+=p->z[k]*p->v[k];
00425                                 sum_v+=p->v[k]; 
00426                                 n++;                    
00427                         }
00428                         
00429                         //double off_u=  (sum_u*sum_z2 - sum_z*sum_zu)/(n*sum_z2-sum_z*sum_z);
00430                         double slope_u=  (n*sum_zu - sum_z*sum_u)/(n*sum_z2-sum_z*sum_z);
00431                         
00432                         
00433                         //double off_v=  (sum_v*sum_z2 - sum_z*sum_zv)/(n*sum_z2-sum_z*sum_z);
00434                         double slope_v=  (n*sum_zv - sum_z*sum_v)/(n*sum_z2-sum_z*sum_z);
00435                         
00436 
00437         
00438 
00439                         double dz=p->end_z-p->start_z;
00440                         double du=dz*slope_u;
00441                         double dv=dz*slope_v;                   
00442                         
00443                         double r2 = du*du+dv*dv+dz*dz;
00444                         r2=sqrt(r2);
00445                         //double pu= p->sum_e/25. * du/r2;
00446                         //double pv= p->sum_e/25. * dv/r2;
00447                         //double pz= p->sum_e/25. * dz/r2;
00448                         
00449                         //printf("particle %d with e %f  (%f %f %f) p3 (%f %f %f)\n",p->particletype,p->sum_e/25.,p->start_u,p->start_v,p->start_z,pu,pv,pz);
00450         
00451         
00452 
00453                 }
00454         
00455         
00456         
00457                 
00458                 if(largest_id>-1)
00459                 {
00460                         Particle3D * p = &poh->particles3d[largest_id];
00461                 
00462 
00463                         //calculate particle curv
00464                         double curv=0;
00465                         for(int i=1;i<p->entries;i++)
00466                         {
00467                                 double du = p->u[i]-p->u[i-1];
00468                                 double dv = p->v[i]-p->v[i-1];
00469                                 double dz = p->z[i]-p->z[i-1];
00470                                 double ds = sqrt(du*du+dv*dv+dz*dz);
00471                                 if(dz>0)curv+=ds/dz;
00472                         }
00473                         hist_sum_curv_vs_part_e[type]->Fill(curv,p->sum_e/25.,weight);
00474         
00475 
00476                         double sum_z2=0;
00477                         double sum_z=0;
00478                         double sum_u=0;
00479                         double sum_zu=0;                        
00480                         double sum_v=0;
00481                         double sum_zv=0;
00482                         int n=0;        
00483                                                 
00484                         for(int i=0;i<p->entries;i++)
00485                         {
00486                                 sum_z2+=p->z[i]*p->z[i];
00487                                 sum_z+=p->z[i];
00488                                 sum_zu+=p->z[i]*p->u[i];
00489                                 sum_u+=p->u[i];
00490                                 sum_zv+=p->z[i]*p->v[i];
00491                                 sum_v+=p->v[i]; 
00492                                 n++;                    
00493                         }
00494                         
00495                         double off_u=  (sum_u*sum_z2 - sum_z*sum_zu)/(n*sum_z2-sum_z*sum_z);
00496                         double slope_u=  (n*sum_zu - sum_z*sum_u)/(n*sum_z2-sum_z*sum_z);
00497                         
00498                         
00499                         double off_v=  (sum_v*sum_z2 - sum_z*sum_zv)/(n*sum_z2-sum_z*sum_z);
00500                         double slope_v=  (n*sum_zv - sum_z*sum_v)/(n*sum_z2-sum_z*sum_z);
00501                         
00502 /*
00504                         double dz=p->end_z-p->start_z;
00505                         double du=dz*slope_u;
00506                         double dv=dz*slope_v;                   
00507                         
00508                         double r2 = du*du+dv*dv+dz*dz;
00509                         r2=sqrt(r2);
00510                         double pu= p->sum_e/25. * du/r2;
00511                         double pv= p->sum_e/25. * dv/r2;
00512                         double pz= p->sum_e/25. * dz/r2;
00513 
00514 
00515                         if(p->particletype==Particle3D::electron)
00516                         {
00517                                 printf("largest elec %d with e %f  (%f %f %f) p3 (%f %f %f)\n",p->particletype,p->sum_e/25.,p->start_u,p->start_v,p->start_z,pu,pv,pz);
00518                                 if(true_elec_idx>-1)
00519                                 {
00520                                         NtpMCStdHep *mc = &(prec->mctrue.stdhep[true_elec_idx]);
00521                         
00522                                         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]);             
00523                                         //calculate the angle!!!!
00524                                         
00525                                         //angle between p vectors....
00526                                         
00527                                         double true_pu=(mc->p4[0]+mc->p4[1])/sqrt(2.0);
00528                                         double true_pv=(-mc->p4[0]+mc->p4[1])/sqrt(2.0);
00529                                         double true_pz=mc->p4[2];
00530                                         
00531                                         double reco_r = sqrt(pu*pu+pv*pv+pz*pz);        
00532                                         double true_r = sqrt(true_pu*true_pu+true_pv*true_pv+true_pz*true_pz);
00533                                         
00534                                         double angle = acos( (pu*true_pu+pv*true_pv+pz*true_pz) / (reco_r*true_r) );
00535                                         printf("angle between: %f\n",angle);
00536                                         
00537                                         hist_prim_angle_vs_prim_eres[type]->Fill(angle,(sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass)-p->sum_e/25.)/(sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass)),weight);                
00538                                         hist_prim_angle_vs_prim_true[type]->Fill(angle,sqrt(mc->p4[3]*mc->p4[3]-mc->mass*mc->mass));    
00539                                         
00540                                         
00541                                 }
00542                         }
00543                         
00545         */              
00546                         //calculate rms
00547                         
00548                         double sumdist=0;
00549                         double weighted_sumdist=0;
00550                         double sume=0;
00551                         
00552                         for(int i=0;i<p->entries;i++)
00553                         {
00554                                 double pu = slope_u*p->z[i]+off_u;      
00555                                 double pv = slope_v*p->z[i]+off_v;
00556                                 
00557                                 sumdist += sqrt((pu*p->u[i])*(pu*p->u[i])+(pv*p->v[i])*(pv*p->v[i]));
00558                                 
00559                                 //printf("pt %f %f true %f %f\n",pu,pv,p->u[i],p->v[i]);
00560                                 
00561                                 weighted_sumdist += p->e[i]*sqrt((pu*p->u[i])*(pu*p->u[i])+(pv*p->v[i])*(pv*p->v[i]));
00562                                 sume+=p->e[i];  
00563                         }       
00564                         weighted_sumdist/=sume;
00565                         sumdist/=n;             
00566                         hist_sumdist[type]->Fill(sumdist,n,weight);
00567                         hist_weighted_sumdist[type]->Fill(weighted_sumdist,n,weight);   
00568                 
00569                 
00570                         double max_e=0;
00571                         double max_e_dist=0;
00572                         for(int i=0;i<p->entries;i++)
00573                         {
00574                                 if(max_e<p->e[i])
00575                                 {
00576                                         max_e=p->e[i];
00577                                         max_e_dist=sqrt((p->u[i]-p->u[0])*(p->u[i]-p->u[0]) + (p->v[i]-p->v[0])*(p->v[i]-p->v[0]) + (p->z[i]-p->z[0])*(p->z[i]-p->z[0])   );
00578                                 
00579                                 
00580                                 }
00581                         }
00582                 
00583                 
00584                         hist_max_dist[type]->Fill(max_e_dist,weight);
00585                         hist_max_dist_per_gev[type]->Fill(max_e_dist/sume*25.0,weight);
00586         
00587         
00588 
00589         
00590         
00591         
00594                         
00595                         //draw a line from the first to the last point and calculation the energy weighted r^2 from that line
00596                         
00597                         int ent=p->entries-1;
00598                         double slopeu = (p->u[ent]-p->u[0])/(p->z[ent]-p->z[0]);
00599                         double slopev = (p->v[ent]-p->v[0])/(p->z[ent]-p->z[0]);                        
00600                         
00601                         double offu = p->u[0]-p->z[0]*slopeu;
00602                         double offv = p->v[0]-p->z[0]*slopev;                   
00603                 
00604                         double sum_dr=0;
00605                         
00606                         for(int i=0;i<p->entries;i++)
00607                         {
00608                                 double eu=p->z[i]*slopeu+offu;
00609                                 double ev=p->z[i]*slopev+offv;
00610                                 
00611                                 double dr = sqrt((eu-p->u[i])*(eu-p->u[i]) + (ev-p->v[i])*(ev-p->v[i]));
00612 
00613                                 sum_dr+=dr*p->e[i];
00614                         }
00615                         
00616                         sum_dr/=sume;
00617                         hist_weighted_area[type]->Fill(sum_dr,weight);
00618                         
00620                 
00621                         int fp=-1;
00622                         double fpe=0;
00623                         int sp=-1;
00624                         //second peak distance
00625                         for(int i=1;i<p->entries;i++)
00626                         {               
00627                         
00628                                 if(i==1 &&p->e[0]>p->e[1])fp=0;
00629                                 
00630                                 if(i<p->entries-1)
00631                                 if(p->e[i-1]<p->e[i] && p->e[i+1]<p->e[i])//peak
00632                                 {
00633                                         if(fp<0){
00634                                                 fp=i;
00635                                                 fpe=p->e[i];
00636                                         }
00637                                         else if(p->e[i] *0.5 > fpe)
00638                                         {
00639                                                 sp=i;
00640                                                 break;
00641                                         }
00642                                 
00643                                 }
00644                         
00645                         
00646                                 if(i==p->entries-1 && p->e[i-1]<p->e[i])
00647                                 {
00648                                         if(fp>-1)sp=i;
00649                                         else fp=i;
00650                                 }
00651                         
00652                         }
00653                         
00654                         if(fp>-1)hist_first_peak[type]->Fill(p->z[fp]-p->z[0],weight);
00655                                 else hist_first_peak[type]->Fill(0.,weight);
00656                         if(sp>-1)hist_second_peak[type]->Fill(p->z[sp]-p->z[0],weight);
00657                                 else hist_second_peak[type]->Fill(0.,weight);
00658                         if(fp>-1 && sp>-1)hist_twopeak_diff[type]->Fill(p->z[sp]-p->z[fp],weight);      
00659                                 else hist_twopeak_diff[type]->Fill(0.,weight);
00660 
00661 
00662 
00664                         double maxe=0;
00665                         for(int i=0;i<p->entries;i++)
00666                         {               
00667                                 if(maxe<p->e[i])maxe=p->e[i];
00668                         }
00669                         
00670                         int pkidx=-1;
00671                         for(int i=p->entries-1;i>-1;i--)
00672                         {               
00673                                 if(maxe*0.7<p->e[i])
00674                                 {
00675                                         pkidx=i;
00676                                         break;
00677                                 }
00678                         }                       
00679                 
00680         //              if(pkidx>-1 && p->z[pkidx]-p->z[0]>0.3)continue;
00681                 
00682                         if(pkidx>-1)
00683                                 hist_back_peak_dist[type]->Fill(p->z[pkidx]-p->z[0],weight);
00684                         else
00685                                 hist_back_peak_dist[type]->Fill(0.,weight);
00687 
00688 
00689                         //view e asym
00690                         
00691                         double eu=0;
00692                         double ev=0;
00693                         for(int i=0;i<p->entries;i++)
00694                         {       
00695                                 if(p->view[i]==2)eu+=p->e[i];
00696                                 if(p->view[i]==3)ev+=p->e[i];
00697                         }
00698                         
00699                         hist_view_e_asym[type]->Fill(fabs(eu-ev)/(eu+ev),weight);
00700 
00701 
00702                 }
00703 
00704 
00705 
00706                         typecnt[type]+=weight;
00707 
00708 
00709 
00710                 if(ent_idx%10000==0)printf("Entry %d\n",ent_idx);
00711         }
00712 
00713 
00714         printf("used %d entries\n",saved);
00715 
00716 
00717 
00718 printf("values\n");
00719 for(int i=0;i<5;i++)
00720         printf("%d   %2.2f\n",i,typecnt[i]);
00721         
00722 printf("FOM %2.2f\n",typecnt[1]/sqrt(typecnt[0]+typecnt[2]+typecnt[3]+typecnt[4]));
00723 
00724 
00725         TFile *f = new TFile("detailplots.root","RECREATE");
00726         f->cd();
00727         for(int i=0;i<5;i++)
00728         {
00729                 hist_sumdist[i]->Write();
00730                 hist_weighted_sumdist[i]->Write();
00731         
00732                 hist_max_dist[i]->Write();
00733                 hist_max_dist_per_gev[i]->Write();      
00734                 
00735                 hist_weighted_area[i]->Write();
00736 
00737                 hist_first_peak[i]->Write();
00738                 hist_second_peak[i]->Write();
00739                 hist_twopeak_diff[i]->Write();
00740                 hist_back_peak_dist[i]->Write();
00741                 hist_view_e_asym[i]->Write();
00742 
00743                 hist_prim_angle_vs_prim_eres[i]->Write();               
00744                 hist_prim_angle_vs_prim_true[i]->Write();       
00745 
00746 
00747                 hist_sum_curv_vs_part_e[i]->Write();    
00748                 hist_pu_pz[i]->Write(); 
00749                 hist_pv_pz[i]->Write(); 
00750                 hist_pt_pz[i]->Write(); 
00751                 hist_pu_pv[i]->Write();
00752 
00753         }
00754         f->Close();
00755 
00756 
00757 }


Member Data Documentation

Definition at line 19 of file DetailedParticle.h.

Referenced by DetailedParticle(), and Run().

Definition at line 22 of file DetailedParticle.h.

Referenced by DetailedParticle(), and Run().

Definition at line 20 of file DetailedParticle.h.

Referenced by DetailedParticle(), and Run().

Definition at line 23 of file DetailedParticle.h.

Referenced by DetailedParticle(), PassesPreselec(), and Run().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1