loon: DetailedParticle.cxx Source File

00001 #include "NueAna/ParticlePID/ParticleAna/DetailedParticle.h"
00002 #include "TH3F.h"
00003 #include "NueAna/ParticlePID/ParticleFinder/Particle3D.h"
00004 
00005 #include "TH2F.h"
00006 #include "TFile.h"
00007 #include "TMath.h"
00008 
00009 #include "MCNtuple/NtpMCStdHep.h"
00010 
00011 #include "NueAna/ParticlePID/ParticleAna/TruthCompareAna.h"
00012 #include "MessageService/MsgService.h"
00013 
00014 DetailedParticle::DetailedParticle(const char*fname)
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 }
00025 
00026 DetailedParticle::~DetailedParticle()
00027 {}
00028 
00029 
00030 int DetailedParticle::PassesPreselec()
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 }
00108 
00109 void DetailedParticle::Run()
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 }
00758 
00759 
00760 
00761 /*
00762 
00763 TH3F * d_e_x_elec = new TH3F("d_e_x_elec","d_e_x_elec",100,0,1000,100,0,1000,100,0,10);
00764         TH3F * d_e_x_muon = new TH3F("d_e_x_muon","d_e_x_muon",100,0,1000,100,0,1000,100,0,10);
00765         TH3F * d_e_x_prot = new TH3F("d_e_x_prot","d_e_x_prot",100,0,1000,100,0,1000,100,0,10);
00766         TH3F * d_e_x_neut = new TH3F("d_e_x_neut","d_e_x_neut",100,0,1000,100,0,1000,100,0,10);         
00767         TH3F * d_e_x_other = new TH3F("d_e_x_other","d_e_x_other",100,0,1000,100,0,1000,100,0,10);
00768         
00769         TH2F * par_b_elec = new TH2F("par_b_elec","par_b_elec",100,0,10,100,0,100);
00770         TH2F * par_b_other = new TH2F("par_b_other","par_b_other",100,0,10,100,0,100);
00771         TH2F * par_b_else = new TH2F("par_b_else","par_b_else",100,0,10,100,0,100);
00772         
00773         
00774         TH2F * par_b_elec_rms = new TH2F("par_b_elec_rms","par_b_elec_rms",100,0,10,1000,0,0.02);
00775         TH2F * par_b_other_rms = new TH2F("par_b_other_rms","par_b_other_rms",100,0,10,1000,0,0.02);
00776 
00777         TH2F * tmax_vs_e = new TH2F("tmax_vs_e","tmax_vs_e",100,0,10,1000,0,10);
00778         
00779         TH2F * fitvscale = new TH2F("fitvscale","fitvscale",100,0,10000,100,0,10000);
00780         
00781         TH1F * fitvscaleres = new TH1F("fitvscaleres","fitvscaleres",100,-0.2,0.2);
00782         
00783         //MINE
00784         TH2F * neut_tail_z = new TH2F("neut_tail_z", "neut_tail_z", 50, -1.5,0.5,100,0,200);
00785         TH3F * neut_tail_uv = new TH3F("neut_tail_uv", "neut_tail_uv",50, -3,2,50,-1,1,100,0,200);
00786         int event_n_elec, event_n_neut;
00787         double n_tail_z, n_tail_u, n_tail_v, nvtxZ, evtxZ, nvtxU, evtxU, nvtxV, evtxV;
00788         //end MINE
00789         
00790         int ent =input->GetEntries();
00791         printf("%d total entries\n",ent);
00792         for(int i=0;i<ent;i++)
00793         {
00794                 input->GetEntry(i);
00795                 //printf("snarl %d\n",poh->GetHeader().GetSnarl());
00796                 
00797                 
00798                 int type = (poh->mctrue.iaction>0)*(abs(poh->mctrue.inu)/2-5);
00799                 
00800                 if(truetype>-1)if(truetype!=type)continue;
00801                 
00802 //              if(poh->mctrue.iaction!=0)continue;
00803 //              if(poh->mctrue.iresonance!=1001)continue;
00804 //              if(poh->mctrue.inu!=12)continue;
00805         
00806 //              if(poh->particles3d1->GetEntries()>1)continue;
00807                 
00808                 Particle3D * maxParticle=0;
00809 
00810                 //MINE
00811                 event_n_elec = 0;
00812                 event_n_neut = 0;
00813                 //end MINE
00814                 
00815                                 
00816                 for(int j=0;j<poh->particles3d.size();j++)
00817                 {
00818 
00819                         Particle3D * p = &poh->particles3d[j];
00820                         if(!p)continue;
00821                         
00822                         if(maxParticle==0 || maxParticle->sum_e<p->sum_e)maxParticle=p;
00823                         
00824                         double dex=p->sum_e / (p->end_z-p->start_z);
00825                         //double dex=p->sum_e / (double)(p->entries);
00826                         
00827                         switch(p->particletype)
00828                         {
00829                                 case Particle3D::neutron:
00830                                         d_e_x_neut->Fill(dex,p->sum_e,(p->end_z-p->start_z));
00831                                         
00832                                         //MINE
00833                                         event_n_neut++;
00834                                         nvtxZ = p->start_z;
00835                                         nvtxU = p->start_u;
00836                                         nvtxV = p->start_v;
00837                                         //end MINE                                      
00838                                         break;
00839                                 case Particle3D::proton:
00840                                         d_e_x_prot->Fill(dex,p->sum_e,(p->end_z-p->start_z));
00841                                         break;
00842                                 case Particle3D::muon:
00843                                         d_e_x_muon->Fill(dex,p->sum_e,(p->end_z-p->start_z));
00844                                         break;
00845                                 case Particle3D::electron:
00846                                         d_e_x_elec->Fill(dex,p->sum_e,(p->end_z-p->start_z));
00847                                         if(p->emfit_b>0)
00848                                         {
00849                                                 par_b_elec->Fill(p->emfit_b,p->sum_e);  
00850                                                 par_b_elec_rms->Fill(p->emfit_b,p->avg_rms_t);
00851                                         }
00852                                         
00853                                         //MINE
00854                                         event_n_elec++;
00855                                         evtxZ = p->end_z;
00856                                         evtxU = p->end_u;
00857                                         evtxV = p->end_v;
00858                                         //end MINE
00859                                         break;
00860                                 case Particle3D::other:
00861                                         d_e_x_other->Fill(dex,p->sum_e,(p->end_z-p->start_z));
00862                                         if(p->emfit_b>0)
00863                                         {
00864                                                 par_b_other->Fill(p->emfit_b,p->sum_e);
00865                                                 par_b_other_rms->Fill(p->emfit_b,p->avg_rms_t); 
00866                                         }
00867                                         break;
00868                                         
00869                                 //else
00870                                         if(p->emfit_b>0)par_b_else->Fill(p->emfit_b,p->sum_e);  
00871                                         
00872                                         
00873                         
00874                         }
00875                         
00876                         if(p->emfit_e0>0)
00877                         {
00878                                 fitvscale->Fill(p->emfit_e0,p->calibrated_energy);
00879                                 fitvscaleres->Fill((p->emfit_e0-p->calibrated_energy)/p->calibrated_energy);
00880                         }
00881                         
00882                         //MINE
00883                         if(event_n_neut ==1 && event_n_elec ==1){
00884                                         n_tail_z = nvtxZ - evtxZ;
00885                                         n_tail_u = nvtxU - evtxZ;
00886                                         n_tail_v = nvtxV - evtxV;
00887                                         neut_tail_z->Fill(n_tail_z,p->sum_e);
00888                                         neut_tail_uv->Fill(n_tail_u, n_tail_v,p->sum_e);
00889                         }
00890                         //end MINE
00891                 }
00892 
00893                 if(maxParticle!=0)      if(maxParticle->emfit_b>0 && maxParticle->sum_e>0)tmax_vs_e->Fill((maxParticle->emfit_a-1)/maxParticle->emfit_b,TMath::Log(maxParticle->sum_e));//,poh->event.oscprob);
00894 
00895 
00896                 if(i%10000==0)printf("Entry %d\n",i);
00897         }
00898 
00899         TFile *f = new TFile("detailplots.root","RECREATE");
00900         f->cd();
00901         d_e_x_elec->Write();
00902         d_e_x_muon->Write();
00903         d_e_x_prot->Write();
00904         d_e_x_neut->Write();
00905         d_e_x_other->Write();
00906         par_b_other->Write();
00907         par_b_elec->Write();
00908         par_b_else->Write();
00909         par_b_other_rms->Write();
00910         par_b_elec_rms->Write();
00911         tmax_vs_e->Write();
00912         //MINE
00913         neut_tail_z->Write();
00914         neut_tail_uv->Write();
00915         fitvscale->Write();
00916         //end MINE
00917         f->Close();
00918 
00919 */
00920 
00921 
00922