// $Id$ // // File: JEventProcessor_EpEmStudies.cc // Created: Fri Jun 28 10:04:30 EDT 2019 // Creator: staylor (on Linux ifarm1401.jlab.org 3.10.0-327.el7.x86_64 x86_64) // #include "JEventProcessor_EpEmStudies.h" using namespace jana; // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_EpEmStudies()); } } // "C" //------------------ // JEventProcessor_EpEmStudies (Constructor) //------------------ JEventProcessor_EpEmStudies::JEventProcessor_EpEmStudies() { } //------------------ // ~JEventProcessor_EpEmStudies (Destructor) //------------------ JEventProcessor_EpEmStudies::~JEventProcessor_EpEmStudies() { } //------------------ // init //------------------ jerror_t JEventProcessor_EpEmStudies::init(void) { // This is called once at program startup. japp->CreateLock("myEpEmlock"); gDirectory->mkdir("EpEmAnalysis")->cd(); EpEmConfidenceLevel=new TH1F("EpEmConfidenceLevel","CL",100,0,1); EpEmChiSq=new TH1F("EpEmChisSq","#chi^{2}/ndf",1000,0,100); EpEmMassNoGamma= new TH1F("EpEmMassNoGamma","e+e- mass",3500,0.,3.5); EpEmMassNoGamma->SetXTitle("M(e^{+}e^{-}) [GeV]"); EpEmMassNoGamma->SetYTitle("counts / 0.001 GeV"); gDirectory->cd("../"); gDirectory->mkdir("EpEmGammaAnalysis")->cd(); EpEmGammaConfidenceLevel=new TH1F("EpEmGammaConfidenceLevel","CL",100,0,1); EpEmGammaChiSq=new TH1F("EpEmGammaChisSq","#chi^{2}/ndf",1000,0,100); EpEmGammaMass= new TH1F("EpEmGammaMass","e+e-#gamma mass",3500,0.,3.5); EpEmGammaMass->SetXTitle("M(e^{+}e^{-}#gamma) [GeV]"); EpEmGammaMass->SetYTitle("counts / 0.001 GeV"); EpEmMass= new TH1F("EpEmMass","e+e- mass",3500,0.,3.5); EpEmMass->SetXTitle("M(e^{+}e^{-}) [GeV]"); EpEmMass->SetYTitle("counts / 0.001 GeV"); gDirectory->mkdir("MC")->cd(); gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); ElectronTOFdt=new TH2F("ElectronTOFdt","#deltat(TOF)",100,0,10,200,-5,5); ElectronFCALdt=new TH2F("ElectronFCALdt","#deltat(FCAL)",100,0,10,200,-5,5); ElectronBCALdt=new TH2F("ElectronBCALdt","#deltat(BCAL)",100,0,10,200,-5,5); ElectronTOFdtPull=new TH2F("ElectronTOFdtPull","#deltat(TOF)/#sigmat(TOF)",100,0,10,200,-5,5); ElectronFCALdtPull=new TH2F("ElectronFCALdtPull","#deltat(FCAL)/#sigmat(FCAL)",100,0,10,200,-5,5); ElectronBCALdtPull=new TH2F("ElectronBCALdtPull","#deltat(BCAL)/#sigmat(BCAL)",100,0,10,200,-5,5); ElectronTOFdtCL=new TH1F("ElectronTOFdtCL","CL for TOF",100,0,1); ElectronFCALdtCL=new TH1F("ElectronFCALdtCL","CL for FCAL",100,0,1); ElectronBCALdtCL=new TH1F("ElectronBCALdtCL","CL for BCAL",100,0,1); ElectrondEdxCDC=new TH2F("ElectrondEdxCDC","dEdx vs p/M",1000,0,20000,500,0,25.0); ElectrondEdxCDCPull=new TH2F("ElectrondEdxCDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M",1000,0,20000,100,-5,5.0); ElectrondEdxCDCCL=new TH1F("ElectrondEdxCDCCL","CL for CDC dEdx",100,0,1); ElectrondEdxFDC=new TH2F("ElectrondEdxFDC","dEdx vs p/M",1000,0,20000,500,0,25.0); ElectrondEdxFDCPull=new TH2F("ElectrondEdxFDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M",1000,0,20000,100,-5,5.0); ElectrondEdxFDCCL=new TH1F("ElectrondEdxFDCCL","CL for FDC dEdx",100,0,1); ElectronEOverPFCAL=new TH2F("ElectronEOverPFCAL","E/p",100,0,10,125,0,1.25); ElectronEOverPBCAL=new TH2F("ElectronEOverPBCAL","E/p",100,0,10,125,0,1.25); ElectronEOverPPullFCAL=new TH2F("ElectronEOverPPullFCAL","E/p pull",100,0,10,100,-5,5); ElectronEOverPPullBCAL=new TH2F("ElectronEOverPPullBCAL","E/p pull",100,0,10,100,-5,5); ElectronEOverPFCALCL=new TH1F("ElectronEOverPFCALCL","CL for FCAL E/p",100,0,1); ElectronEOverPBCALCL=new TH1F("ElectronEOverPBCALCL","CL for BCAL E/p",100,0,1); gDirectory->mkdir("purity")->cd(); EpDenom=new TH2F("EpDenom","P vs #theta",100,0,100,100,0.05,10.05); EpNumer_Kp=new TH2F("EpNumer_Kp","P vs #theta",100,0,100,100,0.05,10.05); EpNumer_Pip=new TH2F("EpNumer_Pip","P vs #theta",100,0,100,100,0.05,10.05); EpNumer_P=new TH2F("EpNumer_P","P vs #theta",100,0,100,100,0.05,10.05); EpNumer_Ep=new TH2F("EpNumer_Ep","P vs #theta",100,0,100,100,0.05,10.05); EmDenom=new TH2F("EmDenom","P vs #theta",100,0,100,100,0.05,10.05); EmNumer_Em=new TH2F("EmNumer_Em","P vs #theta",100,0,100,100,0.05,10.05); EmNumer_Pim=new TH2F("EmNumer_Pim","P vs #theta",100,0,100,100,0.05,10.05); EmNumer_aP=new TH2F("EmNumer_aP","P vs #theta",100,0,100,100,0.05,10.05); EmNumer_Km=new TH2F("EmNumer_Km","P vs #theta",100,0,100,100,0.05,10.05); gDirectory->cd("../../../"); CL_CUT=0.001; gPARMS->SetDefaultParameter("PHISTUDIES:CL_CUT",CL_CUT); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_EpEmStudies::brun(JEventLoop *eventLoop, int32_t runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_EpEmStudies::evnt(JEventLoop *loop, uint64_t eventnumber) { vectorctracks; loop->Get(ctracks); if (ctracks.size()!=3) return RESOURCE_UNAVAILABLE; vectorneutrals; loop->Get(neutrals); if (neutrals.size()>1) return VALUE_OUT_OF_RANGE; bool use_photon=false; const DNeutralParticleHypothesis *gamma=NULL; if (neutrals.size()==1){ gamma=neutrals[0]->Get_Hypothesis(Gamma); double E_gam=gamma->lorentzMomentum().E(); const DNeutralShower *shower=gamma->Get_NeutralShower(); if (shower->dDetectorSystem==SYS_FCAL){ if (E_gam>0.1 && shower->dQuality>0.5) use_photon=true; } else if (E_gam>0.05) use_photon=true; if (use_photon==false) gamma=NULL; } vectorbeamphotons; loop->Get(beamphotons); vector pid_algorithms; loop->Get(pid_algorithms); if(pid_algorithms.size()<1){ _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<WriteLock("myEpEmlock"); 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){ DAnalysisUtilities *dAnalysisUtilities=new DAnalysisUtilities(loop); DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop); DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils); dKinFitUtils->Reset_NewEvent(); dKinFitter->Reset_NewEvent(); 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; } const DChargedTrackHypothesis *electron_hyp=ctracks[neg_index]->Get_Hypothesis(Electron); if (electron_hyp!=NULL){ const DTrackTimeBased *electron_track=electron_hyp->Get_TrackTimeBased(); if (electron_track!=NULL){ double t0_rf=electron_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 (int k=0;k<2;k++){ const DChargedTrackHypothesis *positron_hyp=ctracks[posindexes[k][0]]->Get_Hypothesis(Positron); if (positron_hyp==NULL) continue; const DChargedTrackHypothesis *proton_hyp=ctracks[posindexes[k][1]]->Get_Hypothesis(Proton); if (proton_hyp==NULL) continue; const DTrackTimeBased *positron_track=positron_hyp->Get_TrackTimeBased(); if (positron_track==NULL) continue; const DTrackTimeBased *proton_track=proton_hyp->Get_TrackTimeBased(); if (proton_track==NULL) continue; //-------------------------------- // Kinematic fit //-------------------------------- dKinFitter->Reset_NewFit(); set> InitialParticles, FinalParticles; shared_ptrmyBeam=dKinFitUtils->Make_BeamParticle(beamphotons[i]); shared_ptrmyTarget=dKinFitUtils->Make_TargetParticle(Proton); InitialParticles.insert(myBeam); InitialParticles.insert(myTarget); shared_ptrmyProton=dKinFitUtils->Make_DetectedParticle(proton_track); shared_ptrmyPositron=dKinFitUtils->Make_DetectedParticle(positron_track); shared_ptrmyElectron=dKinFitUtils->Make_DetectedParticle(electron_track); FinalParticles.insert(myProton); FinalParticles.insert(myPositron); FinalParticles.insert(myElectron); TLorentzVector vert(proton_track->position(),t0_rf); if (gamma!=NULL){ shared_ptrmyGamma1=dKinFitUtils->Make_DetectedParticle(1,0,0.,vert,gamma->momentum(),0.,gamma->errorMatrix()); FinalParticles.insert(myGamma1); } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // Decay Vertex constraint set> locFullConstrainParticles; locFullConstrainParticles.insert(myElectron); locFullConstrainParticles.insert(myPositron); DVector3 second_vertex; dAnalysisUtilities->Calc_DOCAVertex(myElectron.get(),myPositron.get(),second_vertex); shared_ptr locSecondaryVertexConstraint = dKinFitUtils->Make_VertexConstraint(locFullConstrainParticles,{},second_vertex); dKinFitter->Add_Constraint(locSecondaryVertexConstraint); // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); // GET THE FIT RESULTS double locConfidenceLevel = dKinFitter->Get_ConfidenceLevel(); double reducedChiSq=dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF()); if (gamma!=NULL){ EpEmGammaConfidenceLevel->Fill(locConfidenceLevel); EpEmGammaChiSq->Fill(reducedChiSq); if (locConfidenceLevel>CL_CUT){ if (EpEmGammaAnalysis(dKinFitter, ctracks[posindexes[k][1]],// proton ctracks[posindexes[k][0]], //positron ctracks[neg_index], // electron weight)){ FillPIDHistos(t0_rf,proton_hyp,proton_track,positron_hyp, positron_track,electron_hyp,electron_track, pid_algorithm); } } } else{ EpEmConfidenceLevel->Fill(locConfidenceLevel); EpEmChiSq->Fill(reducedChiSq); if (locConfidenceLevel>CL_CUT){ if (EpEmAnalysis(dKinFitter, ctracks[posindexes[k][1]],// proton ctracks[posindexes[k][0]], //positron ctracks[neg_index], // electron weight)){ /* FillPIDHistos(t0_rf,proton_hyp,proton_track,positron_hyp, positron_track,electron_hyp,electron_track, pid_algorithm); */ } } } } } // loop over beam photons } // check for electron track } delete dAnalysisUtilities; delete dKinFitter; delete dKinFitUtils; } japp->Unlock("myEpEmlock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_EpEmStudies::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_EpEmStudies::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } void JEventProcessor_EpEmStudies::FillPIDHistos(double t0_rf, const DChargedTrackHypothesis *proton_hyp, const DTrackTimeBased *proton_track, const DChargedTrackHypothesis *positron_hyp, const DTrackTimeBased *positron_track, const DChargedTrackHypothesis *electron_hyp, const DTrackTimeBased *electron_track, const DParticleID *pid_algorithm){ double p=electron_track->momentum().Mag(); double beta=p/electron_track->energy(); double gammabeta=p/ParticleMass(Electron); // double dEdxAmp=1e6*electron_track->ddEdx_CDC_amp; double dEdxAmp=electron_hyp->Get_dEdx_CDC_amp()*1e6; double dEdx_fdc=1e6*electron_track->ddEdx_FDC; if (dEdxAmp>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_CDC(beta,0,dEdxMean,Electron); pid_algorithm->GetdEdxSigma_CDC(beta,0,dEdxSigma,Electron); double dEdxDiff=1e-6*dEdxAmp-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; ElectrondEdxCDC->Fill(gammabeta,dEdxAmp); ElectrondEdxCDCPull->Fill(gammabeta,dEdxPull); ElectrondEdxCDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } if (dEdx_fdc>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_FDC(beta,0,dEdxMean,Electron); pid_algorithm->GetdEdxSigma_FDC(beta,0,dEdxSigma,Electron); double dEdxDiff=electron_track->ddEdx_FDC-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; ElectrondEdxFDC->Fill(gammabeta,dEdx_fdc); ElectrondEdxFDCPull->Fill(gammabeta,dEdxPull); ElectrondEdxFDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } shared_ptrbcalparms=electron_hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=electron_hyp->Get_FCALShowerMatchParams(); shared_ptrtofparms=electron_hyp->Get_TOFHitMatchParams(); if (bcalparms!=NULL){ double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf; // -pid_algorithm->GetTimeMean(SYS_BCAL,Electron,p); double E_over_p=bcalparms->dBCALShower->E/p; ElectronEOverPBCAL->Fill(p,E_over_p); ElectronBCALdt->Fill(p,dt_bcal); double sigma_E_over_p=pid_algorithm->GetEOverPSigma(SYS_BCAL,p); double diff_E_over_p=E_over_p-pid_algorithm->GetEOverPMean(SYS_BCAL,p); double pull_E_over_p=diff_E_over_p/sigma_E_over_p; double chisq_E_over_p=pull_E_over_p*pull_E_over_p; ElectronEOverPPullBCAL->Fill(p,pull_E_over_p); ElectronEOverPBCALCL->Fill(TMath::Prob(chisq_E_over_p,1.)); double var=pid_algorithm->GetTimeVariance(SYS_BCAL,Electron,p); double pull=dt_bcal/sqrt(var); double chisq=dt_bcal*dt_bcal/var; ElectronBCALdtPull->Fill(p,pull); ElectronBCALdtCL->Fill(TMath::Prob(chisq,1.)); } if (tofparms!=NULL){ double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf; ElectronTOFdt->Fill(p,dt); double var=pid_algorithm->GetTimeVariance(SYS_TOF,Electron,p); double pull=dt/sqrt(var); double chisq=dt*dt/var; ElectronTOFdtPull->Fill(p,pull); ElectronTOFdtCL->Fill(TMath::Prob(chisq,1.)); } if (fcalparms!=NULL){ double dt_fcal=fcalparms->dFCALShower->getTime()-fcalparms->dFlightTime-t0_rf;//-pid_algorithm->GetTimeMean(SYS_FCAL,Electron,p);; double E_over_p=fcalparms->dFCALShower->getEnergy()/p; ElectronEOverPFCAL->Fill(p,E_over_p); ElectronFCALdt->Fill(p,dt_fcal); double sigma_E_over_p=pid_algorithm->GetEOverPSigma(SYS_FCAL,p); double diff_E_over_p=E_over_p-pid_algorithm->GetEOverPMean(SYS_FCAL,p); double pull_E_over_p=diff_E_over_p/sigma_E_over_p; double chisq_E_over_p=pull_E_over_p*pull_E_over_p; ElectronEOverPPullFCAL->Fill(p,pull_E_over_p); ElectronEOverPFCALCL->Fill(TMath::Prob(chisq_E_over_p,1.)); double var=pid_algorithm->GetTimeVariance(SYS_FCAL,Electron,p); double pull=dt_fcal/sqrt(var); double chisq=dt_fcal*dt_fcal/var; ElectronFCALdtPull->Fill(p,pull); ElectronFCALdtCL->Fill(TMath::Prob(chisq,1.)); } } bool JEventProcessor_EpEmStudies::EpEmGammaAnalysis(DKinFitter *dKinFitter, const DChargedTrack *proton_hyp, const DChargedTrack *eplus_hyp, const DChargedTrack *eminus_hyp, double weight){ DLorentzVector beam_kf,gamma_kf; DLorentzVector proton_kf,ep_kf,em_kf; set>myParticles=dKinFitter->Get_KinFitParticles(); set>::iterator locParticleIterator=myParticles.begin(); for(; locParticleIterator != myParticles.end(); ++locParticleIterator){ if ((*locParticleIterator)->Get_KinFitParticleType()==d_BeamParticle){ beam_kf=(*locParticleIterator)->Get_P4(); } else if ((*locParticleIterator)->Get_KinFitParticleType()==d_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 11: // positron ep_kf=(*locParticleIterator)->Get_P4(); break; case -11: // electron em_kf=(*locParticleIterator)->Get_P4(); break; case 1: // Gamma gamma_kf=(*locParticleIterator)->Get_P4(); break; default: break; } } } DLorentzVector epem_kf=ep_kf+em_kf; DLorentzVector epemg_kf=ep_kf+em_kf+gamma_kf; EpEmMass->Fill(epem_kf.M(),weight); EpEmGammaMass->Fill(epemg_kf.M(),weight); if (weight>0 && fabs(epemg_kf.M()-ParticleMass(Pi0))<0.026){ double Ep_theta=ep_kf.Theta()*180./M_PI; double Ep_momentum=ep_kf.P(); EpDenom->Fill(Ep_theta,Ep_momentum); const DTrackTimeBased *best_fom_track=eplus_hyp->Get_BestFOM()->Get_TrackTimeBased(); switch(best_fom_track->PID()){ case PiPlus: EpNumer_Pip->Fill(Ep_theta,Ep_momentum); break; case KPlus: EpNumer_Kp->Fill(Ep_theta,Ep_momentum); break; case Proton: EpNumer_P->Fill(Ep_theta,Ep_momentum); break; case Positron: EpNumer_Ep->Fill(Ep_theta,Ep_momentum); break; default: break; } double em_theta=em_kf.Theta()*180./M_PI; double em_momentum=em_kf.P(); EmDenom->Fill(em_theta,em_momentum); best_fom_track=eminus_hyp->Get_BestFOM()->Get_TrackTimeBased(); switch(best_fom_track->PID()){ case PiMinus: EmNumer_Pim->Fill(em_theta,em_momentum); break; case KMinus: EmNumer_Km->Fill(em_theta,em_momentum); break; case AntiProton: EmNumer_aP->Fill(em_theta,em_momentum); break; case Electron: EmNumer_Em->Fill(em_theta,em_momentum); break; default: break; } return true; } return false; } bool JEventProcessor_EpEmStudies::EpEmAnalysis(DKinFitter *dKinFitter, const DChargedTrack *proton_hyp, const DChargedTrack *eplus_hyp, const DChargedTrack *eminus_hyp, double weight){ DLorentzVector beam_kf; DLorentzVector proton_kf,ep_kf,em_kf; set>myParticles=dKinFitter->Get_KinFitParticles(); set>::iterator locParticleIterator=myParticles.begin(); for(; locParticleIterator != myParticles.end(); ++locParticleIterator){ if ((*locParticleIterator)->Get_KinFitParticleType()==d_BeamParticle){ beam_kf=(*locParticleIterator)->Get_P4(); } else if ((*locParticleIterator)->Get_KinFitParticleType()==d_DetectedParticle){ switch ((*locParticleIterator)->Get_PID()){ case 2212: // Proton proton_kf=(*locParticleIterator)->Get_P4(); break; case 11: // positron ep_kf=(*locParticleIterator)->Get_P4(); break; case -11: // electron em_kf=(*locParticleIterator)->Get_P4(); break; default: break; } } } DLorentzVector epem_kf=ep_kf+em_kf; EpEmMassNoGamma->Fill(epem_kf.M(),weight); return false; }