// $Id$ // // File: JEventProcessor_PhiStudies.cc // Created: Thu Jun 27 13:35:27 EDT 2019 // Creator: staylor (on Linux ifarm1401.jlab.org 3.10.0-327.el7.x86_64 x86_64) // #include "JEventProcessor_PhiStudies.h" using namespace jana; #include #include #include #include #define PHI_CUT 0.012 #define LAMBDA_CUT 0.012 #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_PhiStudies()); } } // "C" //------------------ // JEventProcessor_PhiStudies (Constructor) //------------------ JEventProcessor_PhiStudies::JEventProcessor_PhiStudies() { } //------------------ // ~JEventProcessor_PhiStudies (Destructor) //------------------ JEventProcessor_PhiStudies::~JEventProcessor_PhiStudies() { } //------------------ // init //------------------ jerror_t JEventProcessor_PhiStudies::init(void) { // This is called once at program startup. japp->CreateLock("myKFlock"); gDirectory->mkdir("PhiAnalysis")->cd(); /* gDirectory->mkdir("PhiTree")->cd(); string leaves = "Eg/F:t:M:prob:weight"; phiTree=new TTree("phiTree","Phi Analysis"); phiTree->Branch("T",treeval,leaves.c_str()); gDirectory->cd("../"); */ KpKmConfidenceLevel=new TH1F("KpKmConfidenceLevel","CL",100,0,1); KpKmChiSq=new TH1F("KpKmChisSq","#chi^{2}/ndf",1000,0,100); KpKmMass= new TH1F("KpKmMass","K+K- mass",2000,0.9,2.9); KpKmMass->SetXTitle("M(K^{+}K^{-}) [GeV]"); KpKmMass->SetYTitle("counts / 0.001 GeV"); KmPMass=new TH1F("KmPMass","pK- mass",2000,1.25,3.25); KmPMass->SetXTitle("M(pK^{-}) [GeV]"); KmPMass->SetYTitle("counts / 0.001 GeV"); NumCombos=new TH1F("NumCombos","Number of combos",2,0.5,2.5); TwoGammaMass=new TH1F("ExtraTwoGammaMass","m(2#gamma) [GeV]",100,0,1); DocaToTrack=new TH2F("DocaToTrack","d(track-fcal) vs r",100,0,100,100,0,100); gDirectory->mkdir("MC")->cd(); EgammaMC_vs_t=new TH2F("EgammaMC_vs_t","E vs t",50,0,2,86,2.9,11.5); gDirectory->cd("../"); gDirectory->mkdir("t-dependence")->cd(); for (unsigned int i=0;i<86;i++){ char myname[80]; sprintf(myname,"KpKm_vs_t_%3.2fGeV_NoPID",2.95+0.1*i); KpKm_vs_t[i]=new TH2F(myname,myname,50,0,2,300,0.9,1.2); } gDirectory->cd("../"); gDirectory->mkdir("PID")->cd(); EBCAL_vs_EFCAL=new TH2F("EBCAL_vs_EFCAL","E(BCAL) vs E(FCAL)",1000,0,10,1000,0,10); EBCAL_vs_EFCAL_PhiCut=new TH2F("EBCAL_vs_EFCAL_PhiCut","E(BCAL) vs E(FCAL)",1000,0,10,1000,0,10); KMinusMomentumVsTheta=new TH2F("KMinusMomentumVsTheta","p vs theta",100,0,100,100,0,10); KMinusMomentumVsThetaMatched=new TH2F("KMinusMomentumVsThetaMatched","p vs theta",100,0,100,100,0,10); KmTOFdt=new TH2F("KmTOFdt","#deltat(TOF)",100,0,10,200,-5,5); KmFCALdt=new TH2F("KmFCALdt","#deltat(FCAL)",100,0,10,200,-5,5); KmBCALdt=new TH2F("KmBCALdt","#deltat(BCAL)",100,0,10,200,-5,5); KmTOFdtPull=new TH2F("KmTOFdtPull","#deltat(TOF)/#sigmat(TOF)",100,0,10,200,-5,5); KmFCALdtPull=new TH2F("KmFCALdtPull","#deltat(FCAL)/#sigmat(FCAL)",100,0,10,200,-5,5); KmBCALdtPull=new TH2F("KmBCALdtPull","#deltat(BCAL)/#sigmat(BCAL)",100,0,10,200,-5,5); KmTOFdtCL=new TH1F("KmTOFdtCL","CL for TOF",100,0,1); KmFCALdtCL=new TH1F("KmFCALdtCL","CL for FCAL",100,0,1); KmBCALdtCL=new TH1F("KmBCALdtCL","CL for BCAL",100,0,1); KmdEdxCDC=new TH2F("KmdEdxCDC","dEdx vs p/M, K- cands.",1000,0,10,500,0,25.0); KmdEdxCDCPull=new TH2F("KmdEdxCDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, K- cands.",1000,0,10,100,-5,5.0); KmdEdxCDCCL=new TH1F("KmdEdxCDCCL","CL for CDC dEdx",100,0,1); KmdEdxFDC=new TH2F("KmdEdxFDC","dEdx vs p/M, K- cands.",1000,0,10,500,0,25.0); KmdEdxFDCPull=new TH2F("KmdEdxFDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, K- cands.",1000,0,10,100,-5,5.0); KmdEdxFDCCL=new TH1F("KmdEdxFDCCL","CL for FDC dEdx",100,0,1); KmEOverPFCAL=new TH2F("KmEOverPFCAL","E/p",100,0,10,125,0,1.25); KmEFCAL=new TH2F("KmEFCAL","E(FCAL)",100,0,10,125,0,1.25); KmEFCALvsR=new TH2F("KmEFCALvsR","E(FCAL) vs R ",120,0,120,125,0,1.25); KPlusMomentumVsTheta=new TH2F("KPlusMomentumVsTheta","p vs theta",100,0,100,100,0,10); KPlusMomentumVsThetaMatched=new TH2F("KPlusMomentumVsThetaMatched","p vs theta",100,0,100,100,0,10); KpTOFdt=new TH2F("KpTOFdt","#deltat(TOF)",100,0,10,200,-5,5); KpFCALdt=new TH2F("KpFCALdt","#deltat(FCAL)",100,0,10,200,-5,5); KpBCALdt=new TH2F("KpBCALdt","#deltat(BCAL)",100,0,10,200,-5,5); KpTOFdtPull=new TH2F("KpTOFdtPull","#deltat(TOF)/#sigmat(TOF)",100,0,10,200,-5,5); KpFCALdtPull=new TH2F("KpFCALdtPull","#deltat(FCAL)/#sigmat(FCAL)",100,0,10,200,-5,5); KpBCALdtPull=new TH2F("KpBCALdtPull","#deltat(BCAL)/#sigmat(BCAL)",100,0,10,200,-5,5); KpTOFdtCL=new TH1F("KpTOFdtCL","CL for TOF",100,0,1); KpFCALdtCL=new TH1F("KpFCALdtCL","CL for FCAL",100,0,1); KpBCALdtCL=new TH1F("KpBCALdtCL","CL for BCAL",100,0,1); KpdEdxCDC=new TH2F("KpdEdxCDC","dEdx vs p/M, K- cands.",1000,0,10,500,0,25.0); KpdEdxCDCPull=new TH2F("KpdEdxCDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, K- cands.",1000,0,10,100,-5,5.0); KpdEdxCDCCL=new TH1F("KpdEdxCDCCL","CL for CDC dEdx",100,0,1); KpdEdxFDC=new TH2F("KpdEdxFDC","dEdx vs p/M, K- cands.",1000,0,10,500,0,25.0); KpdEdxFDCPull=new TH2F("KpdEdxFDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, K- cands.",1000,0,10,100,-5,5.0); KpdEdxFDCCL=new TH1F("KpdEdxFDCCL","CL for FDC dEdx",100,0,1); KpEOverPFCAL=new TH2F("KpEOverPFCAL","E/p",100,0,10,125,0,1.25); KpEFCAL=new TH2F("KpEFCAL","E(FCAL)",100,0,10,125,0,1.25); PTOFdt=new TH2F("PTOFdt","#deltat(TOF)",100,0,10,200,-5,5); PFCALdt=new TH2F("PFCALdt","#deltat(FCAL)",100,0,10,200,-5,5); PBCALdt=new TH2F("PBCALdt","#deltat(BCAL)",100,0,10,200,-5,5); PTOFdtPull=new TH2F("PTOFdtPull","#deltat(TOF)/#sigmat(TOF)",100,0,10,200,-5,5); PFCALdtPull=new TH2F("PFCALdtPull","#deltat(FCAL)/#sigmat(FCAL)",100,0,10,200,-5,5); PBCALdtPull=new TH2F("PBCALdtPull","#deltat(BCAL)/#sigmat(BCAL)",100,0,10,200,-5,5); PTOFdtCL=new TH1F("PTOFdtCL","CL for TOF",100,0,1); PFCALdtCL=new TH1F("PFCALdtCL","CL for FCAL",100,0,1); PBCALdtCL=new TH1F("PBCALdtCL","CL for BCAL",100,0,1); PdEdxCDCPull=new TH2F("PdEdxCDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, proton cands.",1000,0,10,100,-5,5.0); PdEdxCDCCL=new TH1F("PdEdxCDCCL","CL for CDC dEdx",100,0,1); PdEdxFDCPull=new TH2F("PdEdxFDCPull","#delta(dE/dx)/#sigma(dE/dx) vs p/M, proton cands.",1000,0,10,100,-5,5.0); PdEdxFDCCL=new TH1F("PdEdxFDCCL","CL for CDC dEdx",100,0,1); PdEdxCDC=new TH2F("PdEdxCDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0); PdEdxFDC=new TH2F("PdEdxFDC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0); PdEdxSC=new TH2F("PdEdxSC","dEdx vs p/M, proton cands.",1000,0,10,500,0,25.0); gDirectory->mkdir("purity")->cd(); KpDenom=new TH2F("KpDenom","P vs #theta",100,0,100,100,0.05,10.05); KpNumer_Kp=new TH2F("KpNumer_Kp","P vs #theta",100,0,100,100,0.05,10.05); KpNumer_Pip=new TH2F("KpNumer_Pip","P vs #theta",100,0,100,100,0.05,10.05); KpNumer_P=new TH2F("KpNumer_P","P vs #theta",100,0,100,100,0.05,10.05); KpNumer_Ep=new TH2F("KpNumer_Ep","P vs #theta",100,0,100,100,0.05,10.05); KmDenom=new TH2F("KmDenom","P vs #theta",100,0,100,100,0.05,10.05); KmNumer_Km=new TH2F("KmNumer_Km","P vs #theta",100,0,100,100,0.05,10.05); KmNumer_Pim=new TH2F("KmNumer_Pim","P vs #theta",100,0,100,100,0.05,10.05); KmNumer_aP=new TH2F("KmNumer_aP","P vs #theta",100,0,100,100,0.05,10.05); KmNumer_Em=new TH2F("KmNumer_Em","P vs #theta",100,0,100,100,0.05,10.05); PDenom=new TH2F("PDenom","P vs #theta",80,0,80,30,0.05,3.05); PNumer_Kp=new TH2F("PNumer_Kp","P vs #theta",80,0,80,30,0.05,3.05); PNumer_Pip=new TH2F("PNumer_Pip","P vs #theta",80,0,80,30,0.05,3.05); PNumer_P=new TH2F("PNumer_P","P vs #theta",80,0,80,30,0.05,3.05); PNumer_Ep=new TH2F("PNumer_Ep","P vs #theta",80,0,80,30,0.05,3.05); gDirectory->cd("../../../"); CL_CUT=0.01; gPARMS->SetDefaultParameter("PHISTUDIES:CL_CUT",CL_CUT); PROTON_FOM_CUT=0.001; gPARMS->SetDefaultParameter("PHISTUDIES:PROTON_FOM_CUT",PROTON_FOM_CUT); CUT_EVENTS_WITH_UNUSED_ENERGY=false; gPARMS->SetDefaultParameter("PHISTUDIES:CUT_EVENTS_WITH_UNUSED_ENERGY",CUT_EVENTS_WITH_UNUSED_ENERGY); EMULATE_TRIGGER=true; gPARMS->SetDefaultParameter("PHISTUDIES:EMULATE_TRIGGER",EMULATE_TRIGGER); CUT_DECAYS=false; gPARMS->SetDefaultParameter("PHISTUDIES:CUT_DECAYS",CUT_DECAYS); ADD_VERTEX_CONSTRAINT=true; gPARMS->SetDefaultParameter("PHISTUDIES:ADD_VERTEX_CONSTRAINT",ADD_VERTEX_CONSTRAINT); APPLY_SPLIT_CUT=false; gPARMS->SetDefaultParameter("PHISTUDIES:APPLY_SPLIT_CUT",APPLY_SPLIT_CUT); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_PhiStudies::brun(JEventLoop *eventLoop, int32_t runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_PhiStudies::evnt(JEventLoop *loop, uint64_t eventnumber) { if (EMULATE_TRIGGER){ vector triggers; loop->Get(triggers); if (triggers.size()>0){ int trig_mask=triggers[0]->Get_L1TriggerBits(); if (!((trig_mask&0x1)||(trig_mask&0x4))){ //cout << "Not a valid trigger: " << trig_mask <tagm_geom; loop->Get(tagm_geom); vectortagh_geom; loop->Get(tagh_geom); vectorfcal_showers; loop->Get(fcal_showers); vectorbcal_showers; loop->Get(bcal_showers); double bcal_sum=0.; double fcal_sum=0.; for (unsigned int i=0;iE; } for (unsigned int i=0;igetEnergy(); } vectorthrowns; loop->Get(throwns); bool DecayedInDetector=false; if (throwns.size()>0){ japp->WriteLock("myKFlock"); const DMCReaction* reaction=NULL; loop->GetSingle(reaction); if (reaction!=NULL){ double mc_beam_energy=reaction->beam.energy(); double Egam=0.; if (tagm_geom.size()&&tagh_geom.size()){ unsigned int counter=0,column=0; if ( tagm_geom[0]->E_to_column(mc_beam_energy,column)){ Egam=0.5*(tagm_geom[0]->getEhigh(column) +tagm_geom[0]->getElow(column)); } else if (tagh_geom[0]->E_to_counter(mc_beam_energy,counter)){ Egam=0.5*(tagh_geom[0]->getEhigh(counter) +tagh_geom[0]->getElow(counter)); } vectormcprotons,mcKps,mcKms; for (unsigned int i=0;iparentid==0){ switch(throwns[i]->type){ case Proton: mcprotons.push_back(throwns[i]->lorentzMomentum()); break; case KPlus: mcKps.push_back(throwns[i]->lorentzMomentum()); break; case KMinus: mcKms.push_back(throwns[i]->lorentzMomentum()); break; default: break; } } else{ if (throwns[i]->charge()<0&&throwns[i]->position().Perp()<60. && throwns[i]->position().z()<345.){ DecayedInDetector=true; break; } } } if (mcprotons.size()==1 && mcKps.size()==1 && mcKms.size()==1){ DLorentzVector recoil=mcprotons[0]; DLorentzVector target(0,0,0,ParticleMass(Proton)); double t_mc=(recoil-target).M2(); EgammaMC_vs_t->Fill(-t_mc,Egam); } } } japp->Unlock("myKFlock"); } if (DecayedInDetector&& CUT_DECAYS){ //cout << "K- decayed!" << endl; return RESOURCE_UNAVAILABLE; } vectorctracks; loop->Get(ctracks); if (ctracks.size()!=3) return RESOURCE_UNAVAILABLE; 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!"<neutrals; loop->Get(neutrals); japp->RootWriteLock(); EBCAL_vs_EFCAL->Fill(fcal_sum,bcal_sum); for (unsigned int i=0;iGet_Hypothesis(Gamma); const DNeutralShower *shower=gamHyp->Get_NeutralShower(); if (shower->dDetectorSystem==SYS_FCAL){ const DFCALShower *fcal=dynamic_cast(shower->dBCALFCALShower); DVector3 pos=fcal->getPosition(); double doca=fcal->getDocaTrack(); DocaToTrack->Fill(pos.Perp(),doca); } } japp->RootUnLock(); if (neutrals.size()>0 && CUT_EVENTS_WITH_UNUSED_ENERGY){ double t0_rf=ctracks[0]->Get_BestTrackingFOM()->t0(); unsigned int num_gamma_in_time=0; for (unsigned int i=0;iGet_Hypothesis(Gamma); DLorentzVector gam1=gamHyp->lorentzMomentum(); double tdiff=gamHyp->time()-t0_rf; if (fabs(tdiff)<2.){ if (APPLY_SPLIT_CUT){ const DNeutralShower *shower=gamHyp->Get_NeutralShower(); if (shower->dQuality<0.5) continue; } num_gamma_in_time++; } } //cout << num_gamma_in_time << endl; if (num_gamma_in_time) return RESOURCE_UNAVAILABLE; } japp->WriteLock("myKFlock"); 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){ 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 *kminus_hyp=ctracks[neg_index]->Get_Hypothesis(KMinus); if (kminus_hyp!=NULL){ const DTrackTimeBased *kminus_track=kminus_hyp->Get_TrackTimeBased(); if (kminus_track!=NULL){ double t0_rf=kminus_hyp->t0(); for (unsigned int i=0;itime()-t0_rf; double weight=1.; bool got_beam_photon=false; if (fabs(dt_st)>6.012&&fabs(dt_st)<18.036){ weight=-1./6.; got_beam_photon=true; } if (fabs(dt_st)<2.004){ got_beam_photon=true; } if (got_beam_photon==false) continue; int num_combos=0; for (int k=0;k<2;k++){ const DChargedTrackHypothesis *kplus_hyp=ctracks[posindexes[k][0]]->Get_Hypothesis(KPlus); if (kplus_hyp==NULL) continue; const DChargedTrackHypothesis *proton_hyp=ctracks[posindexes[k][1]]->Get_Hypothesis(Proton); if (proton_hyp==NULL) continue; const DTrackTimeBased *kplus_track=kplus_hyp->Get_TrackTimeBased(); if (kplus_track==NULL) continue; const DTrackTimeBased *proton_track=proton_hyp->Get_TrackTimeBased(); if (proton_track==NULL) continue; if (proton_hyp->Get_FOM()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_ptrmyKPlus=dKinFitUtils->Make_DetectedParticle(kplus_track); shared_ptrmyKMinus=dKinFitUtils->Make_DetectedParticle(kminus_track); FinalParticles.insert(myProton); FinalParticles.insert(myKPlus); FinalParticles.insert(myKMinus); // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // Production Vertex constraint if (ADD_VERTEX_CONSTRAINT){ set> locFullConstrainParticles, locNoConstrainParticles; locFullConstrainParticles.insert(myKPlus); locFullConstrainParticles.insert(myKMinus); locFullConstrainParticles.insert(myProton); locNoConstrainParticles.insert(myBeam); locNoConstrainParticles.insert(myTarget); shared_ptr locProductionVertexConstraint = dKinFitUtils->Make_VertexConstraint(locFullConstrainParticles, locNoConstrainParticles, proton_track->position()); dKinFitter->Add_Constraint(locProductionVertexConstraint); } // PERFORM THE KINEMATIC FIT dKinFitter->Fit_Reaction(); // GET THE FIT RESULTS double locConfidenceLevel = dKinFitter->Get_ConfidenceLevel(); KpKmConfidenceLevel->Fill(locConfidenceLevel); KpKmChiSq->Fill(dKinFitter->Get_ChiSq()/double(dKinFitter->Get_NDF())); if (locConfidenceLevel>0.){ num_combos++; unsigned int flag=KpKmAnalysis(dKinFitter, ctracks[posindexes[k][1]],// proton ctracks[posindexes[k][0]], //K+ ctracks[neg_index], // K- locConfidenceLevel, weight); if (locConfidenceLevel>CL_CUT){ EBCAL_vs_EFCAL_PhiCut->Fill(fcal_sum,bcal_sum); for (unsigned int m=0;mGet_Hypothesis(Gamma); DLorentzVector gam1=gamHyp->lorentzMomentum(); for (unsigned int n=m+1;nGet_Hypothesis(Gamma); DLorentzVector gam2=gamHyp->lorentzMomentum(); DLorentzVector twogam=gam1+gam2; TwoGammaMass->Fill(twogam.M(),weight); } } if (flag>0){ FillPIDHistos(t0_rf,proton_hyp,proton_track,kplus_hyp, kplus_track,kminus_hyp,kminus_track, pid_algorithm); } } NumCombos->Fill(num_combos); } } } // loop over beam photons } // check for K- track } delete dKinFitter; delete dKinFitUtils; } japp->Unlock("myKFlock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_PhiStudies::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_PhiStudies::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } unsigned int JEventProcessor_PhiStudies::KpKmAnalysis(DKinFitter *dKinFitter, const DChargedTrack *proton_hyp, const DChargedTrack *kplus_hyp, const DChargedTrack *kminus_hyp, double prob, double weight){ DLorentzVector beam_kf,gamma_kf; DLorentzVector proton_kf,kp_kf,km_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 321: // K+ kp_kf=(*locParticleIterator)->Get_P4(); break; case -321: // K- km_kf=(*locParticleIterator)->Get_P4(); break; default: break; } } } DLorentzVector pKm_kf=proton_kf+km_kf; DLorentzVector KpKm_kf=km_kf+kp_kf; DLorentzVector target(0,0,0,ParticleMass(Proton)); double t_mandelstam=(proton_kf-target).M2(); /* treeval[0]=beam_kf.E(); treeval[1]=-t_mandelstam; treeval[2]=KpKm_kf.M(); treeval[3]=prob; treeval[4]=weight; phiTree->Fill(); */ if (probFill(pKm_kf.M(),weight); KpKmMass->Fill(KpKm_kf.M(),weight); if (beam_kf.E()>2.9){ int index=(int)(floor(((beam_kf.E()-2.9)/0.1))); KpKm_vs_t[index]->Fill(-t_mandelstam,KpKm_kf.M(),weight); } unsigned int flag=0; if (fabs(pKm_kf.M()-LAM1520_MASS)0 && flag>0){ double kp_theta=kp_kf.Theta()*180./M_PI; double kp_momentum=kp_kf.P(); KpDenom->Fill(kp_theta,kp_momentum); const DTrackTimeBased *best_fom_track=kplus_hyp->Get_BestFOM()->Get_TrackTimeBased(); switch(best_fom_track->PID()){ case PiPlus: KpNumer_Pip->Fill(kp_theta,kp_momentum); break; case KPlus: KpNumer_Kp->Fill(kp_theta,kp_momentum); break; case Proton: KpNumer_P->Fill(kp_theta,kp_momentum); break; case Positron: KpNumer_Ep->Fill(kp_theta,kp_momentum); break; default: break; } double km_theta=km_kf.Theta()*180./M_PI; double km_momentum=km_kf.P(); KmDenom->Fill(km_theta,km_momentum); best_fom_track=kminus_hyp->Get_BestFOM()->Get_TrackTimeBased(); switch(best_fom_track->PID()){ case PiMinus: KmNumer_Pim->Fill(km_theta,km_momentum); break; case KMinus: KmNumer_Km->Fill(km_theta,km_momentum); break; case AntiProton: KmNumer_aP->Fill(km_theta,km_momentum); break; case Electron: KmNumer_Em->Fill(km_theta,km_momentum); break; default: break; } double proton_theta=proton_kf.Theta()*180./M_PI; double proton_momentum=proton_kf.P(); PDenom->Fill(proton_theta,proton_momentum); best_fom_track=proton_hyp->Get_BestFOM()->Get_TrackTimeBased(); switch(best_fom_track->PID()){ case PiPlus: PNumer_Pip->Fill(proton_theta,proton_momentum); break; case KPlus: PNumer_Kp->Fill(proton_theta,proton_momentum); break; case Proton: PNumer_P->Fill(proton_theta,proton_momentum); break; case Positron: PNumer_Ep->Fill(proton_theta,proton_momentum); break; default: break; } } return flag; } void JEventProcessor_PhiStudies::FillPIDHistos(double t0_rf, const DChargedTrackHypothesis *proton_hyp, const DTrackTimeBased *proton_track, const DChargedTrackHypothesis *kplus_hyp, const DTrackTimeBased *kplus_track, const DChargedTrackHypothesis *kminus_hyp, const DTrackTimeBased *kminus_track, const DParticleID *pid_algorithm){ double p=kminus_track->momentum().Mag(); double theta=180./M_PI*kminus_track->momentum().Theta(); KMinusMomentumVsTheta->Fill(theta,p); double beta=p/kminus_track->energy(); double gammabeta=p/ParticleMass(KMinus); //double dEdxAmp=1e6*kminus_track->ddEdx_CDC_amp; double dEdxAmp=kminus_hyp->Get_dEdx_CDC_amp()*1e6; double dEdx_fdc=1e6*kminus_track->ddEdx_FDC; if (dEdxAmp>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_CDC(beta,0,dEdxMean,KMinus); pid_algorithm->GetdEdxSigma_CDC(beta,0,dEdxSigma,KMinus); double dEdxDiff=1e-6*dEdxAmp-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; KmdEdxCDC->Fill(gammabeta,dEdxAmp); KmdEdxCDCPull->Fill(gammabeta,dEdxPull); KmdEdxCDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } if (dEdx_fdc>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_FDC(beta,0,dEdxMean,KMinus); pid_algorithm->GetdEdxSigma_FDC(beta,0,dEdxSigma,KMinus); double dEdxDiff=kminus_track->ddEdx_FDC-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; KmdEdxFDCPull->Fill(gammabeta,dEdxPull); KmdEdxFDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); KmdEdxFDC->Fill(gammabeta,dEdx_fdc); } shared_ptrbcalparms=kminus_hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=kminus_hyp->Get_FCALShowerMatchParams(); shared_ptrtofparms=kminus_hyp->Get_TOFHitMatchParams(); if (bcalparms!=NULL){ KMinusMomentumVsThetaMatched->Fill(theta,p); double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf; KmBCALdt->Fill(p,dt_bcal); double var=pid_algorithm->GetTimeVariance(SYS_BCAL,KMinus,p); double pull=dt_bcal/sqrt(var); double chisq=dt_bcal*dt_bcal/var; KmBCALdtPull->Fill(p,pull); KmBCALdtCL->Fill(TMath::Prob(chisq,1.)); } if (tofparms!=NULL){ double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf; KmTOFdt->Fill(p,dt); double var=pid_algorithm->GetTimeVariance(SYS_TOF,KMinus,p); double pull=dt/sqrt(var); double chisq=dt*dt/var; KmTOFdtPull->Fill(p,pull); KmTOFdtCL->Fill(TMath::Prob(chisq,1.)); } if (fcalparms!=NULL){ KMinusMomentumVsThetaMatched->Fill(theta,p); double dt_fcal=fcalparms->dFCALShower->getTime()-fcalparms->dFlightTime-t0_rf; KmFCALdt->Fill(p,dt_fcal); double var=pid_algorithm->GetTimeVariance(SYS_FCAL,KMinus,p); double pull=dt_fcal/sqrt(var); double chisq=dt_fcal*dt_fcal/var; KmFCALdtPull->Fill(p,pull); KmFCALdtCL->Fill(TMath::Prob(chisq,1.)); double E_over_p=fcalparms->dFCALShower->getEnergy()/p; KmEOverPFCAL->Fill(p,E_over_p); KmEFCAL->Fill(p,fcalparms->dFCALShower->getEnergy()); KmEFCALvsR->Fill(fcalparms->dFCALShower->getPosition().Perp(),fcalparms->dFCALShower->getEnergy()); } p=kplus_track->momentum().Mag(); theta=180./M_PI*kplus_track->momentum().Theta(); KPlusMomentumVsTheta->Fill(theta,p); beta=p/kplus_track->energy(); gammabeta=p/ParticleMass(KPlus); //dEdxAmp=1e6*kplus_track->ddEdx_CDC_amp; dEdxAmp=kplus_hyp->Get_dEdx_CDC_amp()*1e6; dEdx_fdc=1e6*kplus_track->ddEdx_FDC; if (dEdxAmp>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_CDC(beta,0,dEdxMean,KPlus); pid_algorithm->GetdEdxSigma_CDC(beta,0,dEdxSigma,KPlus); double dEdxDiff=dEdxAmp-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; KpdEdxCDC->Fill(gammabeta,dEdxAmp); KpdEdxCDCPull->Fill(gammabeta,dEdxPull); KpdEdxCDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } if (dEdx_fdc>0){ double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_FDC(beta,0,dEdxMean,KPlus); pid_algorithm->GetdEdxSigma_FDC(beta,0,dEdxSigma,KPlus); double dEdxDiff=kplus_track->ddEdx_FDC-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; KpdEdxFDCPull->Fill(gammabeta,dEdxPull); KpdEdxFDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); KpdEdxFDC->Fill(gammabeta,dEdx_fdc); } bcalparms=kplus_hyp->Get_BCALShowerMatchParams(); fcalparms=kplus_hyp->Get_FCALShowerMatchParams(); tofparms=kplus_hyp->Get_TOFHitMatchParams(); if (bcalparms!=NULL){ KPlusMomentumVsThetaMatched->Fill(theta,p); double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf; KpBCALdt->Fill(p,dt_bcal); double var=pid_algorithm->GetTimeVariance(SYS_BCAL,KPlus,p); double pull=dt_bcal/sqrt(var); double chisq=dt_bcal*dt_bcal/var; KpBCALdtPull->Fill(p,pull); KpBCALdtCL->Fill(TMath::Prob(chisq,1.)); } if (tofparms!=NULL){ double dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf; KpTOFdt->Fill(p,dt); double var=pid_algorithm->GetTimeVariance(SYS_TOF,KPlus,p); double pull=dt/sqrt(var); double chisq=dt*dt/var; KpTOFdtPull->Fill(p,pull); KpTOFdtCL->Fill(TMath::Prob(chisq,1.)); } if (fcalparms!=NULL){ KPlusMomentumVsThetaMatched->Fill(theta,p); double dt_fcal=fcalparms->dFCALShower->getTime()-fcalparms->dFlightTime-t0_rf; KpFCALdt->Fill(p,dt_fcal); double var=pid_algorithm->GetTimeVariance(SYS_FCAL,KPlus,p); double pull=dt_fcal/sqrt(var); double chisq=dt_fcal*dt_fcal/var; KpFCALdtPull->Fill(p,pull); KpFCALdtCL->Fill(TMath::Prob(chisq,1.)); double E_over_p=fcalparms->dFCALShower->getEnergy()/p; KpEOverPFCAL->Fill(p,E_over_p); KpEFCAL->Fill(p,fcalparms->dFCALShower->getEnergy()); } p=proton_track->momentum().Mag(); beta=p/proton_track->energy(); gammabeta=p/ParticleMass(Proton); //dEdxAmp=1e6*proton_track->ddEdx_CDC_amp; dEdxAmp=proton_hyp->Get_dEdx_CDC_amp()*1e6; dEdx_fdc=1e6*proton_track->ddEdx_FDC; if (dEdxAmp>0){ PdEdxCDC->Fill(gammabeta,dEdxAmp); double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_CDC(beta,0,dEdxMean,Proton); pid_algorithm->GetdEdxSigma_CDC(beta,0,dEdxSigma,Proton); double dEdxDiff=1e-6*dEdxAmp-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; PdEdxCDC->Fill(gammabeta,dEdxAmp); PdEdxCDCPull->Fill(gammabeta,dEdxPull); PdEdxCDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } if (dEdx_fdc>0){ PdEdxFDC->Fill(gammabeta,dEdx_fdc); double dEdxMean=0.,dEdxSigma=0.; pid_algorithm->GetdEdxMean_FDC(beta,0,dEdxMean,Proton); pid_algorithm->GetdEdxSigma_FDC(beta,0,dEdxSigma,Proton); double dEdxDiff=proton_track->ddEdx_FDC-dEdxMean; double dEdxPull=dEdxDiff/dEdxSigma; double dEdxChiSq=dEdxPull*dEdxPull; PdEdxFDCPull->Fill(gammabeta,dEdxPull); PdEdxFDCCL->Fill(TMath::Prob(dEdxChiSq,1.)); } shared_ptrscparms=proton_hyp->Get_SCHitMatchParams(); if (scparms!=NULL){ double dEdxSC=1e3*scparms->dEdx; PdEdxSC->Fill(p/ParticleMass(Proton),dEdxSC); } shared_ptrbcalparms2=proton_hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms2=proton_hyp->Get_FCALShowerMatchParams(); shared_ptrtofparms2=proton_hyp->Get_TOFHitMatchParams(); if (bcalparms2!=NULL){ double dt_bcal=bcalparms2->dBCALShower->t-bcalparms2->dFlightTime-t0_rf; PBCALdt->Fill(p,dt_bcal); double var=pid_algorithm->GetTimeVariance(SYS_BCAL,Proton,p); double pull=dt_bcal/sqrt(var); double chisq=dt_bcal*dt_bcal/var; PBCALdtPull->Fill(p,pull); PBCALdtCL->Fill(TMath::Prob(chisq,1.)); } if (tofparms2!=NULL){ double dt=tofparms2->dHitTime-tofparms2->dFlightTime-t0_rf; PTOFdt->Fill(p,dt); double var=pid_algorithm->GetTimeVariance(SYS_TOF,Proton,p); double pull=dt/sqrt(var); double chisq=dt*dt/var; PTOFdtPull->Fill(p,pull); PTOFdtCL->Fill(TMath::Prob(chisq,1.)); } if (fcalparms2!=NULL){ double dt_fcal=fcalparms2->dFCALShower->getTime() -fcalparms2->dFlightTime-t0_rf; PFCALdt->Fill(p,dt_fcal); double var=pid_algorithm->GetTimeVariance(SYS_FCAL,Proton,p); double pull=dt_fcal/sqrt(var); double chisq=dt_fcal*dt_fcal/var; PFCALdtPull->Fill(p,pull); PFCALdtCL->Fill(TMath::Prob(chisq,1.)); } }