// $Id$ // // File: JEventProcessor_TrackStudies.cc // Created: Thu Nov 5 13:15:00 EST 2015 // Creator: staylor (on Linux gluon110.jlab.org 2.6.32-358.23.2.el6.x86_64 x86_64) // #include "JEventProcessor_TrackStudies.h" using namespace jana; #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "TMatrixFSym.h" #include #include #include #include #define PHI_CUT 0.012 #define LAMBDA_CUT 0.012 #define SIGMA0_CUT 0.015 #define PI0_CUT 0.018 #define ETA_CUT 0.018 #define LAM1520_MASS 1.5195 #define LAM1520_CUT 0.012 // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_TrackStudies()); } } // "C" //------------------ // JEventProcessor_TrackStudies (Constructor) //------------------ JEventProcessor_TrackStudies::JEventProcessor_TrackStudies() { } //------------------ // ~JEventProcessor_TrackStudies (Destructor) //------------------ JEventProcessor_TrackStudies::~JEventProcessor_TrackStudies() { } //------------------ // init //------------------ jerror_t JEventProcessor_TrackStudies::init(void) { japp->CreateLock("myTSlock"); DEBUG=0; gPARMS->SetDefaultParameter("TRACKSTUDIES:DEBUG",DEBUG); CL_CUT=0.001; gPARMS->SetDefaultParameter("TRACKSTUDIES:CL_CUT",CL_CUT); NUM_SIG_P=3.; gPARMS->SetDefaultParameter("TRACKSTUDIES:NUM_SIG_P",NUM_SIG_P); NUM_SIG_TANL=3.; gPARMS->SetDefaultParameter("TRACKSTUDIES:NUM_SIG_TANL",NUM_SIG_TANL); NUM_SIG_PHI=3.; gPARMS->SetDefaultParameter("TRACKSTUDIES:NUM_SIG_PHI",NUM_SIG_PHI); HYPOTHESIS_TO_STUDY=9; gPARMS->SetDefaultParameter("TRACKSTUDIES:HYPOTHESIS_TO_STUDY", HYPOTHESIS_TO_STUDY); RING_TO_SKIP=0; gPARMS->SetDefaultParameter("TRKFIT:RING_TO_SKIP",RING_TO_SKIP); PLANE_TO_STUDY=0; gPARMS->SetDefaultParameter("TRKFIT:PLANE_TO_SKIP",PLANE_TO_STUDY); MC_DATA_MATCH_CUT=3; gPARMS->SetDefaultParameter("TRACKSTUDIES:MC_DATA_MATCH_CUT",MC_DATA_MATCH_CUT); KINFIT_CL_CUT=0.001; gPARMS->SetDefaultParameter("TRACKSTUDIES:KINFIT_CL_CUT",KINFIT_CL_CUT); PIM_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:PIM_CL_CUT",PIM_CL_CUT); PIP_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:PIP_CL_CUT",PIP_CL_CUT); PROTON_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:PROTON_CL_CUT",PROTON_CL_CUT); ANTIPROTON_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:ANTIPROTON_CL_CUT",ANTIPROTON_CL_CUT); KM_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:KM_CL_CUT",KM_CL_CUT); KP_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:KP_CL_CUT",KP_CL_CUT); ELECTRON_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:ELECTRON_CL_CUT",ELECTRON_CL_CUT); POSITRON_CL_CUT=0.0; gPARMS->SetDefaultParameter("TRACKSTUDIES:POSITRON_CL_CUT",POSITRON_CL_CUT); // This is called once at program startup. If you are creating // and filling historgrams in this plugin, you should lock the // ROOT mutex like this: // japp->RootWriteLock(); { gDirectory->mkdir("P4Only")->cd(); { P4ConfidenceLevel=new TH1F("P4ConfidenceLevel","CL",1000,0,1); P4ChiSq=new TH1F("P4ChiSq","#chi^{2}/ndf",1000,0,100); for (int i=1;i<7;i++){ char hname[80]; sprintf(hname,"P4PiPlusPull%d",i); P4PiPlusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4PiMinusPull%d",i); P4PiMinusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4ProtonPull%d",i); P4ProtonPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4PiPlusPullVsTheta%d",i); P4PiPlusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4PiMinusPullVsTheta%d",i); P4PiMinusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4ProtonPullVsTheta%d",i); P4ProtonPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); } } gDirectory->cd("../"); gDirectory->mkdir("P4Vertex")->cd(); { P4VertexConfidenceLevel=new TH1F("P4VertexConfidenceLevel","CL",1000,0,1); P4VertexChiSq=new TH1F("P4VertexChiSq","#chi^{2}/ndf",1000,0,100); for (int i=1;i<7;i++){ char hname[80]; sprintf(hname,"P4VertexPiPlusPull%d",i); P4VertexPiPlusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexPiMinusPull%d",i); P4VertexPiMinusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexProtonPull%d",i); P4VertexProtonPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexPiPlusPullVsTheta%d",i); P4VertexPiPlusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4VertexPiMinusPullVsTheta%d",i); P4VertexPiMinusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4VertexProtonPullVsTheta%d",i); P4VertexProtonPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4VertexPiPlusPullECut%d",i); P4VertexPiPlusPullECut[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexPiMinusPullECut%d",i); P4VertexPiMinusPullECut[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexProtonPullECut%d",i); P4VertexProtonPullECut[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"P4VertexPiPlusPullECutVsTheta%d",i); P4VertexPiPlusPullECutVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4VertexPiMinusPullECutVsTheta%d",i); P4VertexPiMinusPullECutVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"P4VertexProtonPullECutVsTheta%d",i); P4VertexProtonPullECutVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); } } gDirectory->cd("../"); gDirectory->mkdir("VertexOnly")->cd(); { VertexConfidenceLevel=new TH1F("VertexConfidenceLevel","CL",1000,0,1); VertexChiSq=new TH1F("VertexChiSq","#chi^{2}/ndf",1000,0,100); for (int i=1;i<7;i++){ char hname[80]; sprintf(hname,"PiPlusPull%d",i); PiPlusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"PiMinusPull%d",i); PiMinusPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"ProtonPull%d",i); ProtonPull[i]=new TH2F(hname,hname,100,-M_PI,M_PI,100,-5,5); sprintf(hname,"PiPlusPullVsTheta%d",i); PiPlusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"PiMinusPullVsTheta%d",i); PiMinusPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); sprintf(hname,"ProtonPullVsTheta%d",i); ProtonPullVsTheta[i]=new TH2F(hname,hname,120,0,120,100,-5,5); } } gDirectory->cd("../"); gDirectory->mkdir("MC")->cd(); { gDirectory->mkdir("Pulls")->cd(); { gDirectory->mkdir("PiMinus")->cd(); { Pim_PullX_vs_phi=new TH2F("Pim_PullX_vs_phi","#Deltax/#sigmax", 180,-M_PI,M_PI,100,-5,5); Pim_PullY_vs_phi=new TH2F("Pim_PullY_vs_phi","#Deltay/#sigmay", 180,-M_PI,M_PI,100,-5,5); Pim_PullZ_vs_phi=new TH2F("Pim_PullZ_vs_phi","#Deltaz/#sigmaz", 180,-M_PI,M_PI,100,-5,5); Pim_PullPx_vs_phi=new TH2F("Pim_PullPx_vs_phi","#DeltaPx/#sigmaPx", 180,-M_PI,M_PI,100,-5,5); Pim_PullPy_vs_phi=new TH2F("Pim_PullPy_vs_phi","#DeltaPy/#sigmaPy", 180,-M_PI,M_PI,100,-5,5); Pim_PullPz_vs_phi=new TH2F("Pim_PullPz_vs_phi","#DeltaPz/#sigmaPz", 180,-M_PI,M_PI,100,-5,5); Pim_PullX_vs_theta=new TH2F("Pim_PullX_vs_theta","#Deltax/#sigmax", 120,0,120,100,-5,5); Pim_PullY_vs_theta=new TH2F("Pim_PullY_vs_theta","#Deltay/#sigmay", 120,0,120,100,-5,5); Pim_PullZ_vs_theta=new TH2F("Pim_PullZ_vs_theta","#Deltaz/#sigmaz", 120,0,120,100,-5,5); Pim_PullPx_vs_theta=new TH2F("Pim_PullPx_vs_theta", "#DeltaPx/#sigmaPx", 120,0,120,100,-5,5); Pim_PullPy_vs_theta=new TH2F("Pim_PullPy_vs_theta", "#DeltaPy/#sigmaPy", 120,0,120,100,-5,5); Pim_PullPz_vs_theta=new TH2F("Pim_PullPz_vs_theta", "#DeltaPz/#sigmaPz", 120,0,120,100,-5,5); } gDirectory->cd("../"); } gDirectory->cd("../"); } gDirectory->cd("../"); } japp->RootUnLock(); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_TrackStudies::brun(JEventLoop *eventLoop, int32_t runnumber) { DApplication* dapp=dynamic_cast(eventLoop->GetJApplication()); dgeom = dapp->GetDGeometry(runnumber); bfield = dapp->GetBfield(runnumber); JCalibration *jcalib = dapp->GetJCalibration(runnumber); // Get the particle ID algorithms // Get the CDC wire table from the XML //jout<< "Getting map of cdc wires from the XML" <GetCDCWires(cdcwires); /* double dphi=0.0001; for (unsigned int i=0;iorigin.Phi(); DVector3 wirepos=cdcwires[i][j]->origin+(75./cdcwires[i][j]->udir.z())*cdcwires[i][j]->udir; double dxd=wirepos.Perp()*(cos(phi+dphi)-cos(phi)); double dyd=wirepos.Perp()*(sin(phi+dphi)-sin(phi)); DVector3 wirepos2=cdcwires[i][j]->origin+(-75./cdcwires[i][j]->udir.z())*cdcwires[i][j]->udir; double dxu=wirepos2.Perp()*(cos(phi+dphi)-cos(phi)); double dyu=wirepos2.Perp()*(sin(phi+dphi)-sin(phi)); cout << dxu << " " << dyu << " " << dxd << " " << dyd << endl; } } */ // Beam position and direction map beam_vals; jcalib->Get("PHOTON_BEAM/beam_spot",beam_vals); beam_center.Set(beam_vals["x"],beam_vals["y"]); beam_dir.Set(beam_vals["dxdz"],beam_vals["dydz"]); beam_z0=beam_vals["z"]; // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_TrackStudies::evnt(JEventLoop *loop, uint64_t eventnumber) { // This is called for every event. Use of common resources like writing // to a file or filling a histogram should be mutex protected. Using // loop->Get(...) to get reconstructed objects (and thereby activating the // reconstruction algorithm) should be done outside of any mutex lock // since multiple threads may call this method at the same time. // Here's an example: // // vector mydataclasses; // loop->Get(mydataclasses); // // japp->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); vectorthrowns; loop->Get(throwns); vectorcands; loop->Get(cands); vectorbeamphotons; loop->Get(beamphotons); /* for (unsigned int i=0;icharge(),cands[i]->position(), cands[i]->momentum(),origin,dir,pos1,pos2); printf("d %f rc %f\n",d,cands[i]->rc); vectorcdchits; cands[i]->GetT(cdchits); for (unsigned int j=0;jcharge(),cands[i]->position(),cands[i]->momentum(), cdchits[j]->wire->origin,cdchits[j]->wire->udir,pos1,pos2); printf("ring %d straw %d t %f d %f\n",cdchits[j]->wire->ring, cdchits[j]->wire->straw,cdchits[j]->tdrift,d); } } */ vectortracks; //vectortracks; loop->Get(tracks); vectorhtracks; //vectortracks; loop->Get(htracks); vectorctracks; loop->Get(ctracks); vectortofpoints; loop->Get(tofpoints); vectortofmcs; loop->Get(tofmcs); vectorpseudos; loop->Get(pseudos); vectorcdchits; loop->Get(cdchits); vectorbcals; loop->Get(bcals); vectorfcals; loop->Get(fcals); vectorcdcdigis; loop->Get(cdcdigis); vectorneutrals; loop->Get(neutrals); vector pid_algorithms; loop->Get(pid_algorithms); if(pid_algorithms.size()<1){ _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<mcpips; vectormcpims; vectormcprotons; for (unsigned int i=0;imech==0){ switch(throwns[i]->type){ case Proton: mcprotons.push_back(throwns[i]); break; case PiPlus: mcpips.push_back(throwns[i]); break; case PiMinus: mcpims.push_back(throwns[i]); break; default: break; } } } if(ctracks.size()>=3){ japp->WriteLock("myTSlock"); if (ctracks.size()==3){ unsigned int neg_index=0; unsigned int num_neg=0; for (int i=0;i<3;i++){ if (ctracks[i]->Get_BestTrackingFOM()->charge()<0){ num_neg++; neg_index=i; } } if (num_neg==1){ unsigned int posindexes[2][2]; if (neg_index==0){ posindexes[0][0]=1; posindexes[0][1]=2; posindexes[1][0]=2; posindexes[1][1]=1; } if (neg_index==1){ posindexes[0][0]=0; posindexes[0][1]=2; posindexes[1][0]=2; posindexes[1][1]=0; } if (neg_index==2){ posindexes[0][0]=1; posindexes[0][1]=0; posindexes[1][0]=0; posindexes[1][1]=1; } DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop); DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils); dKinFitUtils->Reset_NewEvent(); dKinFitter->Reset_NewEvent(); const DChargedTrackHypothesis *piminus_hyp=ctracks[neg_index]->Get_Hypothesis(PiMinus); if (piminus_hyp!=NULL){ const DTrackTimeBased *piminus_track=piminus_hyp->Get_TrackTimeBased(); vector >positive_hyps; for (int k=0;k<2;k++){ const DChargedTrackHypothesis *piplus_hyp=ctracks[posindexes[k][0]]->Get_Hypothesis(PiPlus); if (piplus_hyp==NULL) continue; const DChargedTrackHypothesis *proton_hyp=ctracks[posindexes[k][1]]->Get_Hypothesis(Proton); if (proton_hyp==NULL) continue; positive_hyps.push_back(make_pair(piplus_hyp,proton_hyp)); } if (positive_hyps.size()==2){ // MC comparisons if (throwns.size()==3 && mcpims.size()==1 && mcpips.size()==1 && mcprotons.size()==1){ for (unsigned int k=0;k<2;k++){ MCDataComparison(mcpims[0],mcpips[0],mcprotons[0],piminus_track, positive_hyps[k].first->Get_TrackTimeBased(), positive_hyps[k].second->Get_TrackTimeBased()); } } const DTrackTimeBased *proton_track=NULL; const DTrackTimeBased *piplus_track=NULL; double prob1=positive_hyps[0].second->Get_FOM(); double prob2=positive_hyps[1].second->Get_FOM(); if (prob1>prob2 && prob1>0.01){ proton_track=positive_hyps[0].second->Get_TrackTimeBased(); piplus_track=positive_hyps[0].first->Get_TrackTimeBased(); } else if (prob2>prob1 && prob2>0.01){ proton_track=positive_hyps[1].second->Get_TrackTimeBased(); piplus_track=positive_hyps[1].first->Get_TrackTimeBased(); } if (proton_track!=NULL && piplus_track!=NULL){ VertexOnlyFit(dKinFitter,dKinFitUtils,piminus_track,piplus_track, proton_track); } } double t0_rf=piminus_hyp->t0(); for (unsigned int i=0;itime()-t0_rf; double weight=1.; bool got_beam_photon=false; if (fabs(dt_st)>6.&&fabs(dt_st)<18.){ weight=-1./6.; got_beam_photon=true; } if (fabs(dt_st)<2.){ got_beam_photon=true; } if (got_beam_photon==false) continue; for (unsigned int k=0;kGet_TrackTimeBased(); const DTrackTimeBased *proton_track=positive_hyps[k].second->Get_TrackTimeBased(); P4Fit(dKinFitter,dKinFitUtils,piminus_track,piplus_track, proton_track,beamphotons[i],weight); P4VertexFit(dKinFitter,dKinFitUtils,piminus_track,piplus_track, proton_track,beamphotons[i],weight); } } } delete dKinFitUtils; delete dKinFitter; } } japp->Unlock("myTSlock"); } return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_TrackStudies::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_TrackStudies::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } double JEventProcessor_TrackStudies::CalcDoca(double q,const DVector3 &trackpos, const DVector3 &trackmom, const DVector3 &origin, const DVector3 &dir, DVector3 &pos1, DVector3 &pos2) const{ double Bz=fabs(bfield->GetBz(trackpos.x(),trackpos.y(),trackpos.z())); double twok=-0.003*q*Bz/trackmom.Perp(); double phi=trackmom.Phi(); double cosphi=cos(phi); double cosphi_sq=cosphi*cosphi; double sinphi=sin(phi); double sinphi_sq=sinphi*sinphi; double tanl=tan(M_PI_2-trackmom.Theta()); double tanl_sq=tanl*tanl; double tx=dir.x(); double ty=dir.y(); double tz=dir.z(); double tx_sq=tx*tx,ty_sq=ty*ty,tz_sq=tz*tz; double delta_x=origin.x()-trackpos.x(); double delta_y=origin.y()-trackpos.y(); double delta_z=origin.z()-trackpos.z(); if (fabs(tx)<1e-8 && fabs(ty)<1e-8){ double strack=(delta_x*cosphi+delta_y*sinphi) /(1.+twok*(sinphi*delta_x-cosphi*delta_y)); double sline=tanl*strack-delta_z; pos1.SetXYZ(origin.x()+tx*sline,origin.y()+ty*sline,origin.z()+tz*sline); pos2.SetXYZ(trackpos.x()+cosphi*strack,trackpos.y()+sinphi*strack, trackpos.z()+tanl*strack); //printf("strack %f\n",strack); //pos1.Print(); //pos2.Print(); return (pos1-pos2).Mag(); } double b=4*twok*(-(sinphi*tx) + cosphi*ty) *(cosphi*tx + sinphi*ty + tanl*tz) *(tanl*tz*(tx*delta_x + ty*delta_y) - tanl*(tx_sq + ty_sq)*delta_z + cosphi*(-(ty_sq*delta_x) - tz_sq*delta_x + tx*ty*delta_y + tx*tz*delta_z) + sinphi*(tx*ty*delta_x - tx_sq*delta_y + tz*(-(tz*delta_y) + ty*delta_z))) + pow(tanl_sq*(tx_sq + ty_sq) - 2*cosphi*tanl*tx*tz + sinphi_sq*(tx_sq + tz_sq) - sinphi*(2*cosphi*tx*ty + 2*tanl*ty*tz + twok*(-(ty_sq*delta_x) - tz_sq*delta_x + tx*ty*delta_y + tx*tz*delta_z)) + cosphi*(cosphi*(ty_sq + tz_sq) + twok*(tx*ty*delta_x - tx_sq*delta_y + tz*(-(tz*delta_y) + ty*delta_z))),2); if (b<0) return 1e6; double a=tanl_sq*tx_sq + cosphi_sq*ty_sq + tanl_sq*ty_sq - 2*cosphi*tanl*tx*tz + cosphi_sq*tz_sq + sinphi_sq*(tx_sq + tz_sq)- cosphi*twok*tx_sq*delta_y - cosphi*twok*tz_sq*delta_y; double sline=-(a - cosphi*twok*tx*ty*delta_x - 2*cosphi*twok*ty_sq*delta_y - cosphi*twok*ty*tz*delta_z + sinphi*(-2*cosphi*tx*ty - 2*tanl*ty*tz + twok*(2*tx_sq*delta_x + ty_sq*delta_x + tz_sq*delta_x + tx*ty*delta_y + tx*tz*delta_z)) - sqrt(b) )/(2.*twok*(sinphi*tx - cosphi*ty)*(tx_sq + ty_sq + tz_sq)); double strack=-(a + cosphi*twok*tx*ty*delta_x + cosphi*twok*ty*tz*delta_z - sinphi*(2*cosphi*tx*ty + 2*tanl*ty*tz + twok*(-(ty_sq*delta_x)- tz_sq*delta_x + tx*ty*delta_y + tx*tz*delta_z)) - sqrt(b) )/(2.*twok*(sinphi*tx - cosphi*ty)*(cosphi*tx + sinphi*ty + tanl*tz)); //printf("strack %f\n",strack); pos1.SetXYZ(origin.x()+tx*sline,origin.y()+ty*sline,origin.z()+tz*sline); pos2.SetXYZ(trackpos.x()+cosphi*strack,trackpos.y()+sinphi*strack, trackpos.z()+tanl*strack); return (pos1-pos2).Mag(); } void JEventProcessor_TrackStudies::VertexOnlyFit(DKinFitter *dKinFitter, DKinFitUtils_GlueX *dKinFitUtils, const DTrackTimeBased *piminus_track, const DTrackTimeBased *piplus_track, const DTrackTimeBased *proton_track){ //-------------------------------- // Kinematic fit //-------------------------------- dKinFitter->Reset_NewFit(); set> myParticles; auto myProton=dKinFitUtils->Make_DetectedParticle(proton_track); auto myPiPlus=dKinFitUtils->Make_DetectedParticle(piplus_track); auto myPiMinus=dKinFitUtils->Make_DetectedParticle(piminus_track); myParticles.insert(myProton); myParticles.insert(myPiPlus); myParticles.insert(myPiMinus); auto vertexConstraint=dKinFitUtils->Make_VertexConstraint(myParticles,{}, proton_track->position()); dKinFitter->Add_Constraint(vertexConstraint); dKinFitter->Fit_Reaction(); // GET THE FIT RESULTS double confidenceLevel = dKinFitter->Get_ConfidenceLevel(); double reducedChiSq=dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF()); VertexConfidenceLevel->Fill(confidenceLevel); VertexChiSq->Fill(reducedChiSq); if (confidenceLevel>KINFIT_CL_CUT){ map, map > locParticlePulls; dKinFitter->Get_Pulls(locParticlePulls); auto locMapIterator = locParticlePulls.begin(); for(; locMapIterator != locParticlePulls.end(); ++locMapIterator){ auto myParticle = locMapIterator->first; maplocPulls=locMapIterator->second; auto locMapIterator2 = locPulls.begin(); for(; locMapIterator2 != locPulls.end(); ++locMapIterator2){ TVector3 mom=myParticle->Get_Momentum(); double phi=mom.Phi(); double theta=180./M_PI*mom.Theta(); switch(PDGtoPType(myParticle->Get_PID())){ case PiMinus: PiMinusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second); PiMinusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second); break; case PiPlus: PiPlusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second); PiPlusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second); break; case Proton: ProtonPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second); ProtonPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second); break; default: break; } } } } } void JEventProcessor_TrackStudies::P4Fit(DKinFitter *dKinFitter, DKinFitUtils_GlueX *dKinFitUtils, const DTrackTimeBased *piminus_track, const DTrackTimeBased *piplus_track, const DTrackTimeBased *proton_track, const DBeamPhoton *beam,double weight ){ //-------------------------------- // Kinematic fit //-------------------------------- dKinFitter->Reset_NewFit(); set> initParticles, finalParticles; auto myProton=dKinFitUtils->Make_DetectedParticle(proton_track); auto myPiPlus=dKinFitUtils->Make_DetectedParticle(piplus_track); auto myPiMinus=dKinFitUtils->Make_DetectedParticle(piminus_track); finalParticles.insert(myProton); finalParticles.insert(myPiPlus); finalParticles.insert(myPiMinus); auto myBeam=dKinFitUtils->Make_BeamParticle(beam); auto myTarget=dKinFitUtils->Make_TargetParticle(Proton); initParticles.insert(myBeam); initParticles.insert(myTarget); auto locP4Constraint = dKinFitUtils->Make_P4Constraint(initParticles, finalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); dKinFitter->Fit_Reaction(); // GET THE FIT RESULTS double confidenceLevel = dKinFitter->Get_ConfidenceLevel(); double reducedChiSq=dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF()); P4ConfidenceLevel->Fill(confidenceLevel); P4ChiSq->Fill(reducedChiSq); if (confidenceLevel>KINFIT_CL_CUT){ map, map > locParticlePulls; dKinFitter->Get_Pulls(locParticlePulls); auto locMapIterator = locParticlePulls.begin(); for(; locMapIterator != locParticlePulls.end(); ++locMapIterator){ auto myParticle = locMapIterator->first; maplocPulls=locMapIterator->second; auto locMapIterator2 = locPulls.begin(); for(; locMapIterator2 != locPulls.end(); ++locMapIterator2){ TVector3 mom=myParticle->Get_Momentum(); double phi=mom.Phi(); double theta=180./M_PI*mom.Theta(); switch(PDGtoPType(myParticle->Get_PID())){ case PiMinus: P4PiMinusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4PiMinusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); break; case PiPlus: P4PiPlusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4PiPlusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); break; case Proton: P4ProtonPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4ProtonPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); break; default: break; } } } } } void JEventProcessor_TrackStudies::P4VertexFit(DKinFitter *dKinFitter, DKinFitUtils_GlueX *dKinFitUtils, const DTrackTimeBased *piminus_track, const DTrackTimeBased *piplus_track, const DTrackTimeBased *proton_track, const DBeamPhoton *beam,double weight ){ //-------------------------------- // Kinematic fit //-------------------------------- dKinFitter->Reset_NewFit(); set> initParticles, finalParticles; auto myProton=dKinFitUtils->Make_DetectedParticle(proton_track); auto myPiPlus=dKinFitUtils->Make_DetectedParticle(piplus_track); auto myPiMinus=dKinFitUtils->Make_DetectedParticle(piminus_track); finalParticles.insert(myProton); finalParticles.insert(myPiPlus); finalParticles.insert(myPiMinus); auto myBeam=dKinFitUtils->Make_BeamParticle(beam); auto myTarget=dKinFitUtils->Make_TargetParticle(Proton); initParticles.insert(myBeam); initParticles.insert(myTarget); auto vertexConstraint=dKinFitUtils->Make_VertexConstraint(finalParticles, initParticles, proton_track->position()); dKinFitter->Add_Constraint(vertexConstraint); auto locP4Constraint = dKinFitUtils->Make_P4Constraint(initParticles, finalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); dKinFitter->Fit_Reaction(); // GET THE FIT RESULTS double confidenceLevel = dKinFitter->Get_ConfidenceLevel(); double reducedChiSq=dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF()); P4VertexConfidenceLevel->Fill(confidenceLevel); P4VertexChiSq->Fill(reducedChiSq); if (confidenceLevel>KINFIT_CL_CUT){ double EBeam=beam->lorentzMomentum().E(); map, map > locParticlePulls; dKinFitter->Get_Pulls(locParticlePulls); auto locMapIterator = locParticlePulls.begin(); for(; locMapIterator != locParticlePulls.end(); ++locMapIterator){ auto myParticle = locMapIterator->first; maplocPulls=locMapIterator->second; auto locMapIterator2 = locPulls.begin(); for(; locMapIterator2 != locPulls.end(); ++locMapIterator2){ TVector3 mom=myParticle->Get_Momentum(); double phi=mom.Phi(); double theta=180./M_PI*mom.Theta(); switch(PDGtoPType(myParticle->Get_PID())){ case PiMinus: P4VertexPiMinusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexPiMinusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); if (EBeam>8.2&&EBeam<8.8){ P4VertexPiMinusPullECut[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexPiMinusPullECutVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); } break; case PiPlus: P4VertexPiPlusPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexPiPlusPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); if (EBeam>8.2&&EBeam<8.8){ P4VertexPiPlusPullECut[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexPiPlusPullECutVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); } break; case Proton: P4VertexProtonPull[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexProtonPullVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); if (EBeam>8.2&&EBeam<8.8){ P4VertexProtonPullECut[locMapIterator2->first]->Fill(phi,locMapIterator2->second,weight); P4VertexProtonPullECutVsTheta[locMapIterator2->first]->Fill(theta,locMapIterator2->second,weight); } break; default: break; } } } } } void JEventProcessor_TrackStudies::MCDataComparison(const DMCThrown *mcpim, const DMCThrown *mcpip, const DMCThrown *mcproton, const DTrackTimeBased *pim, const DTrackTimeBased *pip, const DTrackTimeBased *proton){ shared_ptrpim_terrm=pim->TrackingErrorMatrix(); shared_ptrpim_errm=pim->errorMatrix(); double mc_pim_px=mcpim->momentum().x(); double mc_pim_py=mcpim->momentum().y(); double mc_pim_pz=mcpim->momentum().z(); double mc_pim_phi=mcpim->momentum().Phi(); double mc_pim_theta=mcpim->momentum().Theta(); double mc_pim_theta_deg=180./M_PI*mc_pim_theta; double pim_px=pim->momentum().x(); double pim_py=pim->momentum().y(); double pim_pz=pim->momentum().z(); double pim_phi=pim->momentum().Phi(); double pim_theta=pim->momentum().Theta(); double pim_var_px=pim_errm->operator()(0,0); double pim_var_py=pim_errm->operator()(1,1); double pim_var_pz=pim_errm->operator()(2,2); double pim_pull_px=(pim_px-mc_pim_px)/sqrt(pim_var_px); double pim_pull_py=(pim_py-mc_pim_py)/sqrt(pim_var_py); double pim_pull_pz=(pim_pz-mc_pim_pz)/sqrt(pim_var_pz); Pim_PullPx_vs_theta->Fill(mc_pim_theta_deg,pim_pull_px); Pim_PullPy_vs_theta->Fill(mc_pim_theta_deg,pim_pull_py); Pim_PullPz_vs_theta->Fill(mc_pim_theta_deg,pim_pull_pz); Pim_PullPx_vs_phi->Fill(mc_pim_phi,pim_pull_px); Pim_PullPy_vs_phi->Fill(mc_pim_phi,pim_pull_py); Pim_PullPz_vs_phi->Fill(mc_pim_phi,pim_pull_pz); double mc_pim_x=mcpim->position().x(); double mc_pim_y=mcpim->position().y(); double mc_pim_z=mcpim->position().z(); double pim_x=pim->position().x(); double pim_y=pim->position().y(); double pim_z=pim->position().z(); double pim_var_x=pim_errm->operator()(3,3); double pim_var_y=pim_errm->operator()(4,4); double pim_var_z=pim_errm->operator()(5,5); double pim_pull_x=(pim_x-mc_pim_x)/sqrt(pim_var_x); double pim_pull_y=(pim_y-mc_pim_y)/sqrt(pim_var_y); double pim_pull_z=(pim_z-mc_pim_z)/sqrt(pim_var_z); Pim_PullX_vs_theta->Fill(mc_pim_theta_deg,pim_pull_x); Pim_PullY_vs_theta->Fill(mc_pim_theta_deg,pim_pull_y); Pim_PullZ_vs_theta->Fill(mc_pim_theta_deg,pim_pull_z); Pim_PullX_vs_phi->Fill(mc_pim_phi,pim_pull_x); Pim_PullY_vs_phi->Fill(mc_pim_phi,pim_pull_y); Pim_PullZ_vs_phi->Fill(mc_pim_phi,pim_pull_z); }