// $Id$ // // File: JEventProcessor_PIDCheck.cc // Created: Mon Jan 28 08:50:55 EST 2019 // Creator: staylor (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64) // #include "JEventProcessor_PIDCheck.h" using namespace jana; #include #include #include //#include #include #include #include #include #include #include #include #include #include // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_PIDCheck()); } } // "C" //------------------ // JEventProcessor_PIDCheck (Constructor) //------------------ JEventProcessor_PIDCheck::JEventProcessor_PIDCheck() { } //------------------ // ~JEventProcessor_PIDCheck (Destructor) //------------------ JEventProcessor_PIDCheck::~JEventProcessor_PIDCheck() { } //------------------ // init //------------------ jerror_t JEventProcessor_PIDCheck::init(void) { japp->CreateLock("mylock"); // This is called once at program startup. PIM_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:PIM_CL_CUT",PIM_CL_CUT); PIP_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:PIP_CL_CUT",PIM_CL_CUT); KM_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:KM_CL_CUT",PIM_CL_CUT); KP_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:KP_CL_CUT",PIM_CL_CUT); PROTON_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:PROTON_CL_CUT",PROTON_CL_CUT); ANTIPROTON_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:ANTIPROTON_CL_CUT",PROTON_CL_CUT); ELECTRON_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:ELECTRON_CL_CUT",ELECTRON_CL_CUT); POSITRON_CL_CUT=0.01; gPARMS->SetDefaultParameter("PIDCheck:POSITRON_CL_CUT",POSITRON_CL_CUT); CHECK_TOF=false; gPARMS->SetDefaultParameter("PIDCheck:CHECK_TOF",CHECK_TOF); CHECK_FMWPC=false; gPARMS->SetDefaultParameter("PIDCheck:CHECK_FMWPC",CHECK_FMWPC); FCALdt_vs_E=new TH2F("FCALdt_vs_E","dt vs E",100,0,2,100,-10,10); CtofTs=new TH2F("CtofTs",";t (ns);",2000,0.,1000.,8,0.5,8.5); CtofdEs=new TH2F("CtofdEs",";#deltaE (MeV);",100,0.,10.,8,0.5,8.5); CtofT=new TH2F("CtofT",";t (ns);bar",1600,-100.,100.,4,0.5,4.5); CtofdE=new TH2F("CtofdE",";#deltaE (MeV);bar",200,0.,10.,4,0.5,4.5); Ctofdy=new TH2F("Ctofdy",";#deltay (cm);bar",500,-50.,50.,4,0.5,4.5); Ctofdxdy=new TH2F("Ctofdxdy",";#deltax (cm);#deltay (cm)",500,-50.,50., 500,-50.,50.); CtofExtrapXY=new TH2F("CtofExtrapXY",";x (cm);y (cm)",1000,-100.,100., 1000,-100.,100.); CtofExtrapXYwithHit=new TH2F("CtofExtrapXYwithHit",";x (cm);y (cm)",1000,-100.,100., 1000,-100.,100.); CtofYCorr=new TH2F("CtofYCorr",";y(track) (cm);y(ctof) (cm)",1000,-100.,100., 1000,-100.,100.); Ctofdt=new TH2F("Ctofdt",";#deltat (ns);bar",400,-80.,20.,4,0.5,4.5); ITOFxy=new TH2F("ITOFxy",";x[cm];y[cm]",900,-45,45,900,-45,45); gDirectory->mkdir("Positron")->cd(); { PositrondEdxCL=new TH1F("PositrondEdxCL","e+ dEdx CL",100,0,1.); PositronBCALdtCL=new TH1F("PositronBCALdtCL","e+ BCAL dt CL",100,0,1.); PositronFCALdtCL=new TH1F("PositronFCALdtCL","e+ FCAL dt CL",100,0,1.); PositronBCALEOverPCL=new TH1F("PositronBCALEOverPCL","e+ BCAL E/p CL",100,0,1.); PositronEOverPCL=new TH1F("PositronEOverPCL","e+ FCAL E/p CL",100,0,1.); PositrondEdxFDCCL=new TH1F("PositrondEdxFDCCL","e+ dEdx CL",100,0,1.); PositronXdiff_vs_p=new TH2F("PositronXdiff_vs_p",";p [GeV/c];#deltax [cm]", 100,0,10,100,-20,20); PositronYdiff_vs_p=new TH2F("PositronYdiff_vs_p",";p [GeV/c];#deltay [cm]", 100,0,10,100,-20,20); PositronDiff_vs_theta=new TH2F("PositronDiff_vs_theta",";#theta [degrees];#deltay [cm]", 100,0,10,100,0,20); DPositrondEdxFDC=new TH2F("DPositrondEdxFDC","#delta(dEdx) vs p/M, e+ cands.",200,0,20000,1000,-25.0,25.0); PositronTOFdtCL=new TH1F("PositronTOFdtCL","e+ TOF dt CL",100,0,1.); PositronPIDCL=new TH1F("PositronPIDCL","e+ PID CL",100,0,1.); DPositrondEdx=new TH2F("DPositrondEdx","#Delta(dEdx) vs p/M, e+ cands.",200,0,20000,500,-25.,25.0); DEOverPPlus=new TH2F("DEOverPPlus","d(E/p) vs p",100,0,10,200,-1.5,1.5); DEOverPPlus_BCAL=new TH2F("DEOverPPlus_BCAL","d(E/p) vs p",100,0,10,200,-1.5,1.5); EOverPPlus=new TH2F("EOverPPlus","E/p vs p",100,0,10,100,0,1.5); EOverPPlus_BCAL=new TH2F("EOverPPlus_BCAL","E/p vs p",100,0,10,100,0,1.5); EOverPPlus_cut=new TH2F("EOverPPlus_cut","E/p vs p",100,0,10,100,0,1.5); EOverPPlus_BCAL_cut=new TH2F("EOverPPlus_BCAL_cut","E/p vs p",100,0,10,100,0,1.5); PositronTOF=new TH2F("PositronTOF","#deltat vs p for e+",200,0,10,1000,-20,20) ; PositronBCALdt=new TH2F("PositronBCALdt","#deltat vs p for e+",200,0,10,1000,-20,20); PositronFCALdt=new TH2F("PositronFCALdt","#deltat vs p for e+",200,0,10,1000,-20,20); PositrondEdxFDC=new TH2F("PositrondEdxFDC","dEdx vs p/M, e+ cands.", 200,0,20000,500,0,25.0); PositrondEdx=new TH2F("PositrondEdx","dEdx vs p/M, e+ cands.",200,0,20000,500,0,25.0); PositronITOFdt=new TH2F("PositronITOFdt","#deltat vs p for e+",100,0,10,1000,-20,20); PositronITOFdXdY=new TH2F("PositronITOFdXdY","Track matching to ITOF;#deltaX[cm];#deltaY[cm]",400,-20,20,400,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("Electron")->cd(); { ElectrondEdxCL=new TH1F("ElectrondEdxCL","e- dEdx CL",100,0,1.); ElectronBCALdtCL=new TH1F("ElectronBCALdtCL","e- BCAL dt CL",100,0,1.); ElectronFCALdtCL=new TH1F("ElectronFCALdtCL","e- FCAL dt CL",100,0,1.); ElectronTOFdtCL=new TH1F("ElectronTOFdtCL","e- TOF dt CL",100,0,1.); ElectronPIDCL=new TH1F("ElectronPIDCL","e- PID CL",100,0,1.); DElectrondEdx=new TH2F("DElectrondEdx","#Delta(dEdx) vs p/M, e- cands.",200,0,20000,500,-25.,25.0); ElectronBCALEOverPCL=new TH1F("ElectronBCALEOverPCL","e- BCAL E/p CL",100,0,1.); ElectronEOverPCL=new TH1F("ElectronEOverPCL","e- FCAL E/p CL",100,0,1.); ElectronXdiff_vs_p=new TH2F("ElectronXdiff_vs_p",";p [GeV/c];#deltax [cm]", 100,0,10,100,-20,20); ElectronYdiff_vs_p=new TH2F("ElectronYdiff_vs_p",";p [GeV/c];#deltay [cm]", 100,0,10,100,-20,20); ElectronDiff_vs_theta=new TH2F("ElectronDiff_vs_theta",";#theta [degrees];#deltay [cm]", 100,0,10,100,0,20); ElectrondEdxFDCCL=new TH1F("ElectrondEdxFDCCL","e- dEdx CL",100,0,1.); DElectrondEdxFDC=new TH2F("DElectrondEdxFDC","#delta(dEdx) vs p/M, e- cands.",200,0,20000,1000,-25.0,25.0); EOverPMinus_cut=new TH2F("EOverPMinus_cut","E/p vs p",100,0,10,100,0,1.5); EOverPMinus_BCAL_cut=new TH2F("EOverPMinus_BCAL_cut","E/p vs p",100,0,10,100,0,1.5); DEOverPMinus=new TH2F("DEOverPMinus","d(E/p) vs p",100,0,10,200,-1.5,1.5); DEOverPMinus_BCAL=new TH2F("DEOverPMinus_BCAL","d(E/p) vs p",100,0,10,200,-1.5,1.5); EOverPMinus=new TH2F("EOverPMinus","E/p vs p",100,0,10,100,0,1.5); EOverPMinus_BCAL=new TH2F("EOverPMinus_BCAL","E/p vs p",100,0,10,100,0,1.5); ElectronTOF=new TH2F("ElectronTOF","#deltat vs p for e-",200,0,10,1000,-20,20); ElectronBCALdt=new TH2F("ElectronBCALdt","#deltat vs p for e-",200,0,10,1000,-20,20); ElectronBCALdtCut=new TH2F("ElectronBCALdtCut","#deltat vs p for e-",200,0,10,1000,-20,20); ElectronFCALdt=new TH2F("ElectronFCALdt","#deltat vs p for e-",200,0,10,1000,-20,20); ElectrondEdxFDC=new TH2F("ElectrondEdxFDC","dEdx vs p/M, e- cands.",200,0,20000,500,0,25.0); ElectrondEdx=new TH2F("ElectrondEdx","dEdx vs p/M, e- cands.",200,0,20000,500,0,25.0); ElectrondEdxCut=new TH2F("ElectrondEdxCut","dEdx vs p/M, e- cands., cut",200,0,20000,500,0,25.0); ElectrondEdxCut->SetXTitle("p/M"); ElectrondEdxCut->SetYTitle("dEdx [keV/cm]"); ElectrondEdxFDCCut=new TH2F("ElectrondEdxFDCCut","dEdx vs p/M, e- cands., cut",200,0,20000,500,0,25.0); ElectrondEdxFDCCut->SetXTitle("p/M"); ElectrondEdxFDCCut->SetYTitle("dEdx [keV/cm]"); ElectronITOFdt=new TH2F("ElectronITOFdt","#deltat vs p for e-",100,0,10,1000,-20,20); ElectronITOFdXdY=new TH2F("ElectronITOFdXdY","Track matching to ITOF;#deltaX[cm];#deltaY[cm]",400,-20,20,400,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("PiPlus")->cd(); { ExtrapToTOF_PiPlus=new TH2F("ExtrapToTOF_PiPlus","y vs x at TOF, no TOF hit",1000,-100,100,1000,-100,100); ExtrapToFMWPC_PiPlus=new TH2F("ExtrapToFMWPC_PiPlus","y vs x at FMWPC",1000,-100,100,1000,-100,100); FMWPCxy=new TH2F("FMWPCxy","y vs x at FMWPC",1000,-100,100,1000,-100,100); FMWPCdX_vs_p_Pip=new TH2F("FMWPCdX_vs_Pip","Track matching to plane 1;p [GeV/c];#deltax [cm]",100,0,10,500,-50,50); FMWPCdY_vs_p_Pip=new TH2F("FMWPCdY_vs_Pip","Track matching to plane 2;p [GeV]/c]; #deltay [cm]",100,0,10,500,-50,50); FMWPCdXY_vs_layer_Pip=new TH2F("FMWPCdXY_vs_layer_Pip",";layer;Residual [cm]",6,0.5,6.5,500,-50,50); FMWPCTime_PiPlus=new TH2F("FMWPCTime_PiPlus",";layer;t [ns]",6,0.5,6.5,1000,0,1000); FMWPCQ_PiPlus=new TH2F("FMWPCQ_PiPlus",";layer;q [pC]",6,0.5,6.5,500,0,1000); FMWPCAmp_PiPlus=new TH2F("FMWPCAmp_PiPlus",";layer;amplitude [counts]",6,0.5,6.5,500,0,10000); ExtrapToFMWPC_PiPlus_with_FMWPC=new TH2F("ExtrapToFMWPC_PiPlus_with_FMWPC","y vs x at FMWPC",1000,-100,100,1000,-100,100); ExtrapToTOF_NoHoriz=new TH2F("ExtrapToTOF_NoHoriz","y vs x at TOF, no Horiz hit",1000,-100,100,1000,-100,100); ExtrapToTOF_NoVerti=new TH2F("ExtrapToTOF_NoVerti","y vs x at TOF, no Verti hit",1000,-100,100,1000,-100,100); PiPlus_MatchFCALHit_d_vs_p=new TH2F("PiPlus_MatchFCALHit_d_vs_p","d vs p",100,0,10,100,0,20); PiPlus_MatchFCALHit_E_vs_p=new TH2F("PiPlus_MatchFCALHit_E_vs_p","E vs p",100,0,10,100,0,1); PiPlus_MatchFCALHit_E_over_p_vs_p=new TH2F("PiPlus_MatchFCALHit_E_over_p_vs_p","E/p vs p",100,0,10,100,0,2); PiPlus_MatchFCALHit_dt_vs_p=new TH2F("PiPlus_MatchFCALHit_dt_vs_p","dt vs p",100,0,10,100,-10,10); PiPlus_MatchFCALHit_dt_vs_E=new TH2F("PiPlus_MatchFCALHit_dt_vs_E","dt vs E",100,0,2,100,-10,10); PiPlusdEdxCL=new TH1F("PiPlusdEdxCL","#pi+ dEdx CL",100,0,1.); PiPlusBCALdtCL=new TH1F("PiPlusBCALdtCL","#pi+ BCAL dt CL",100,0,1.); PiPlusFCALdtCL=new TH1F("PiPlusFCALdtCL","#pi+ FCAL dt CL",100,0,1.); PiPlusTOFdtCL=new TH1F("PiPlusTOFdtCL","#pi+ TOF dt CL",100,0,1.); PiPlusPIDCL=new TH1F("PiPlusPIDCL","#pi+ PID CL",100,0,1.); PiPlusdEdxFDCCL=new TH1F("PiPlusdEdxFDCCL","pi dEdx CL",100,0,1.); DPiPlusdEdxFDC=new TH2F("DPiPlusEdxFDC","#delta(dEdx) vs p/M, pi cands.",2000,0,80,1000,-25.0,25.0); DPiPlusdEdx=new TH2F("DPiPlusdEdx","#delta(dEdx) vs p/M, pi+ cands.", 2000,0,80,1000,-25.0,25.0); PiPlusdEdx=new TH2F("PiPlusdEdx","dEdx vs p/M, pi+ cands.",2000,0,80,500,0,25.0); PiPlusdEdx->SetXTitle("p/M"); PiPlusdEdx->SetYTitle("dEdx [keV/cm]"); PiPlusdEdxFDC=new TH2F("PiPlusdEdxFDC","dEdx vs p/M, pi+ cands.",2000,0,80,500,0,25.0); PiPlusdEdxCut=new TH2F("PiPlusdEdxCut","dEdx vs p/M, pi+ cands., cut", 2000,0,80,500,0,25.0); PiPlusdEdxCut->SetXTitle("p"); PiPlusdEdxCut->SetYTitle("dEdx [keV/cm]"); PiPlusdEdxFDCCut=new TH2F("PiPlusdEdxFDCCut","dEdx vs p/M, pi+ cands., cut",2000,0,80,500,0,25.0); PiPlusdEdxFDCCut->SetXTitle("p"); PiPlusdEdxFDCCut->SetYTitle("dEdx [keV/cm]"); PiPlusTOFdEdx=new TH2F("PiPlusTOFdEdx","dEdx vs p/M, pi+ cands., cut", 2000,0,80,100,0,10.0e-3); PiPlusTOFdEdx->SetXTitle("p/M"); PiPlusTOFdEdx->SetYTitle("dEdx [GeV/cm]"); PiPlusTOFdEdxCut=new TH2F("PiPlusTOFdEdxCut","dEdx vs p/M, pi+ cands., cut", 2000,0,80,100,0,10.0e-3); PiPlusTOFdEdxCut->SetXTitle("p/M"); PiPlusTOFdEdxCut->SetYTitle("dEdx [GeV/cm]"); PiPlusBCALdt=new TH2F("PiPlusBCALdt","#deltat vs #beta#gamma for pi+",200,0,10,1000,-20,20); PiPlusTOFdt=new TH2F("PiPlusTOFdt","#deltat vs p for pi+",200,0,10,1000,-20,20); PipRvsZ=new TH2F("PipRvsZ","r vs z for pi+ POCA to beam line",400,0,200,100,0,20); PipRvsZ->SetXTitle("z [cm]"); PipRvsZ->SetYTitle("r [cm]"); PiPlusTOFdtCut=new TH2F("PiPlusTOFdtCut","#deltat vs p for pi+",200,0,10,1000,-20,20); PiPlusFCALdt=new TH2F("PiPlusFCALdt","#deltat vs p for pi+",200,0,10,1000,-20,20); PiPlusSCdt=new TH2F("PiPlusSCdt","#deltat vs p for pi+",100,0,5,1000,-20,20); PiPlusBCALdtCut=new TH2F("PiPlusBCALdtCut","#deltat vs p for pi+",100,0,5,1000,-20,20); PiPlusCTOFMatchdXdY=new TH2F("PiPlusCTOFMatchdXdY","Track matching to CTOF;#deltaX[cm];#deltaY[cm]",200,-50,50,200,-50,50); PiPlusCTOFdt=new TH1F("PiPlusCTOFdt","CTOF-track time difference;#deltat [ns]",400,-100,100); PipITOFdXdY=new TH2F("PipITOFdXdY","Track matching to ITOF;#deltaX[cm];#deltaY[cm]",400,-20,20,400,-20,20); PiPlusITOFdt=new TH2F("PiPlusITOFdt","#deltat vs p for pi+",100,0,10,1000,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("PiMinus")->cd(); { ExtrapToTOF_PiMinus=new TH2F("ExtrapToTOF_PiMinus","y vs x at TOF, no TOF hit",1000,-100,100,1000,-100,100); ExtrapToFMWPC_PiMinus=new TH2F("ExtrapToFMWPC_PiMinus","y vs x at FMWPC",1000,-100,100,1000,-100,100); FMWPCdX_vs_p_Pim=new TH2F("FMWPCdX_vs_Pim","dx vs p at FMWPC",100,0,10,200,-25,25); FMWPCdY_vs_p_Pim=new TH2F("FMWPCdY_vs_Pim","dy vs p at FMWPC",100,0,10,200,-25,25); FMWPCdXY_vs_layer_Pim=new TH2F("FMWPCdXY_vs_layer_Pim","dx vs layer at FMWPC",6,0.5,6.5,500,-50,50); FMWPCTime_PiMinus=new TH2F("FMWPCTime_PiMinus",";layer;t [ns]",6,0.5,6.5,1000,0,1000); FMWPCQ_PiMinus=new TH2F("FMWPCQ_PiMinus",";layer;q [pC]",6,0.5,6.5,500,0,1000); FMWPCAmp_PiMinus=new TH2F("FMWPCAmp_PiMinus",";layer;amplitude [counts]",6,0.5,6.5,500,0,10000); ExtrapToFMWPC_PiMinus_with_FMWPC=new TH2F("ExtrapToFMWPC_PiMinus_with_FMWPC","y vs x at FMWPC",1000,-100,100,1000,-100,100); PiMinus_MatchFCALHit_d_vs_p=new TH2F("PiMinus_MatchFCALHit_d_vs_p","d vs p",100,0,10,100,0,20); PiMinus_MatchFCALHit_E_vs_p=new TH2F("PiMinus_MatchFCALHit_E_vs_p","E vs p",100,0,10,100,0,1); PiMinus_MatchFCALHit_E_over_p_vs_p=new TH2F("PiMinus_MatchFCALHit_E_over_p_vs_p","E/p vs p",100,0,10,100,0,2); PiMinus_MatchFCALHit_dt_vs_p=new TH2F("PiMinus_MatchFCALHit_dt_vs_p","dt vs p",100,0,10,100,-10,10); PiMinus_MatchFCALHit_dt_vs_E=new TH2F("PiMinus_MatchFCALHit_dt_vs_E","dt vs E",100,0,2,100,-10,10); PiMinusdEdxCL=new TH1F("PiMinusdEdxCL","#pi- dEdx CL",100,0,1.); PiMinusBCALdtCL=new TH1F("PiMinusBCALdtCL","#pi- BCAL dt CL",100,0,1.); PiMinusFCALdtCL=new TH1F("PiMinusFCALdtCL","#pi- FCAL dt CL",100,0,1.); PiMinusTOFdtCL=new TH1F("PiMinusTOFdtCL","#pi- TOF dt CL",100,0,1.); PiMinusPIDCL=new TH1F("PiMinusPIDCL","#pi- PID CL",100,0,1.); PiMinusdEdxFDCCL=new TH1F("PiMinusdEdxFDCCL","pi dEdx CL",100,0,1.); DPiMinusdEdxFDCVsNall=new TH2F("DPiMinusdEdxFDCVsNall","#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-5.0,5.0); DPiMinusdEdxFDC=new TH2F("DPiMinusdEdxFDC","#delta(dEdx) vs p/M, pi cands.",2000,0,80,1000,-25.0,25.0); for (int i=0;i<40;i++){ char mytitle[80]; sprintf(mytitle,"DPiMinusdEdxFDCVsN%d",i); DPiMinusdEdxFDCVsN[i]=new TH2F(mytitle,"#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-10.0,10.0); } DPiMinusdEdxFDCVsNcorr=new TH2F("DPiMinusdEdxFDCVsNcorr","#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-5.0,5.0); DPiMinusdEdx=new TH2F("DPiMinusdEdx","#delta(dEdx) vs p/M, pi- cands.", 2000,0,80,1000,-25.0,25.0); DPiMinusdEdxVsNall=new TH2F("DPiMinusdEdxVsNall","#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-5.0,5.0); for (int i=0;i<40;i++){ char mytitle[80]; sprintf(mytitle,"DPiMinusdEdxVsN%d",i); DPiMinusdEdxVsN[i]=new TH2F(mytitle,"#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-10.0,10.0); } DPiMinusdEdxVsNcorr=new TH2F("DPiMinusdEdxVsNcorr","#delta(dEdx) vs N, pi- cands.", 28,0.5,28.5,1000,-5.0,5.0); PiMinusTOFdEdx=new TH2F("PiMinusTOFdEdx","dEdx vs p/M, pi- cands., cut", 2000,0,80,100,0,10.0e-3); PiMinusTOFdEdx->SetXTitle("p/M"); PiMinusTOFdEdx->SetYTitle("dEdx [GeV/cm]"); PiMinusTOFdEdxCut=new TH2F("PiMinusTOFdEdxCut","dEdx vs p/M, pi- cands., cut", 2000,0,80,100,0,10.0e-3); PiMinusTOFdEdxCut->SetXTitle("p/M"); PiMinusTOFdEdxCut->SetYTitle("dEdx [GeV/cm]"); PiMinusTOFdt=new TH2F("PiMinusTOFdt","#deltat vs #beta#gamma for pi-",200,0,10,1000,-20,20); PiMinusBCALdt=new TH2F("PiMinusBCALdt","#deltat vs #beta#gamma for pi-",200,0,10,1000,-20,20); PiMinusBCALdtVsTheta=new TH2F("PiMinusBCALdtVsTheta","#deltat vs #theta for pi-",180,0,180,1000,-20,20); PiMinusFCALdt=new TH2F("PiMinusFCALdt","#deltat vs p for pi-",200,0,10,1000,-20,20); PiMinusFCALdtNoTOF=new TH2F("PiMinusFCALdtNoTOF","#deltat vs p for pi-",200,0,10,1000,-20,20); PiMinusdEdxFDC=new TH2F("PiMinusdEdxFDC","dEdx vs p/M, pi- cands.",2000,0,80,500,0,25.0); PiMinusdEdx=new TH2F("PiMinusdEdx","dEdx vs p/M, pi- cands.",2000,0,80,500,0,25.0); PiMinusdEdx->SetXTitle("p/M"); PiMinusdEdx->SetYTitle("dEdx [keV/cm]"); PiMinusdEdxCut=new TH2F("PiMinusdEdxCut","dEdx vs p/M, pi- cands., cut", 2000,0,80,500,0,25.0); PiMinusdEdxCut->SetXTitle("p"); PiMinusdEdxCut->SetYTitle("dEdx [keV/cm]"); PiMinusdEdxFDCCut=new TH2F("PiMinusdEdxFDCCut","dEdx vs p/M, pi- cands., cut",2000,0,80,500,0,25.0); PiMinusdEdxFDCCut->SetXTitle("p"); PiMinusdEdxFDCCut->SetYTitle("dEdx [keV/cm]"); // EOverPMinus_Pim_cut=new TH2F("EOverPMinus_Pim_cut","E/p vs p",100,0,10,100,0,1.5); PimRvsZ=new TH2F("PimRvsZ","r vs z for pi- POCA to beam line",400,0,200,100,0,20); PimRvsZ->SetXTitle("z [cm]"); PimRvsZ->SetYTitle("r [cm]"); PiMinusCTOFMatchdXdY=new TH2F("PiMinusCTOFMatchdXdY","Track matching to CTOF;#deltaX[cm];#deltaY[cm]",200,-50,50,200,-50,50); PiMinusCTOFdt=new TH1F("PiMinusCTOFdt","CTOF-track time difference;#deltat [ns]",400,-100,100); PimITOFdXdY=new TH2F("PimITOFdXdY","Track matching to ITOF;#deltaX[cm];#deltaY[cm]",400,-20,20,400,-20,20); PiMinusITOFdt=new TH2F("PiMinusITOFdt","#deltat vs p for pi-",100,0,10,1000,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("KPlus")->cd(); { KPlusdEdxCL=new TH1F("KplusdEdxCL","K+ dEdx CL",100,0,1.); KPlusBCALdtCL=new TH1F("KPlusBCALdtCL","K+ BCAL dt CL",100,0,1.); KPlusFCALdtCL=new TH1F("KPlusFCALdtCL","K+ FCAL dt CL",100,0,1.); KPlusTOFdtCL=new TH1F("KPlusTOFdtCL","K+ TOF dt CL",100,0,1.); KPlusPIDCL=new TH1F("KPlusPIDCL","K+ PID CL",100,0,1.); KPlusdEdxFDCCL=new TH1F("KPlusdEdxFDCCL","K+ dEdx CL",100,0,1.); DKPlusdEdxFDC=new TH2F("DKPlusdEdxFDC","#delta(dEdx) vs p/M, K- cands.",200,0,24,1000,-25.0,25.0); DKPlusdEdx=new TH2F("DKPlusdEdx","#delta(dEdx) vs p/M, K+ cands.",200,0,4,1000,-25.0,25.0); KPlusdEdxCut=new TH2F("KPlusdEdxCut","dEdx vs p/M, K+ cands.",200,0,4,500,0,25.0); KPlusdEdx=new TH2F("KPlusdEdx","dEdx vs p/M, K+ cands.",200,0,4,500,0,25.0); KPlusdEdx->SetXTitle("p/M"); KPlusdEdx->SetYTitle("dEdx [keV/cm]"); KPlusdEdxFDC=new TH2F("KPlusdEdxFDC","dEdx vs p/M, K+ cands.",200,0,24,500,0,25.0); KPlusdEdxFDC->SetXTitle("p/M"); KPlusdEdxFDC->SetYTitle("dEdx [keV/cm]"); KPlusdEdxFDCCut=new TH2F("KPlusdEdxFDCCut","dEdx vs p/M, K+ cands.",200,0,4,500,0,25.0); KPlusTOFdt=new TH2F("KPlusTOFdt","#deltat vs p for K+",200,0,10,1000,-20,20); KPlusTOFdtCut=new TH2F("KPlusTOFdtCut","#deltat vs p for K+",200,0,10,1000,-20,20); KPlusBCALdt=new TH2F("KPlusBCALdt","#deltat vs p for K+",100,0,5,1000,-20,20); KPlusFCALdt=new TH2F("KPlusFCALdt","#deltat vs p for K+",100,0,5,1000,-20,20); KPlusBCALdtCut=new TH2F("KPlusBCALdtCut","#deltat vs p for K+",100,0,5,1000,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("KMinus")->cd(); { KMinusdEdxCL=new TH1F("KMinusdEdxCL","K- dEdx CL",100,0,1.); KMinusdEdxFDCCL=new TH1F("KMinusdEdxFDCCL","K- dEdx CL",100,0,1.); KMinusBCALdtCL=new TH1F("KMinusBCALdtCL","K- BCAL dt CL",100,0,1.); KMinusFCALdtCL=new TH1F("KMinusFCALdtCL","K- FCAL dt CL",100,0,1.); KMinusTOFdtCL=new TH1F("KMinusTOFdtCL","K- TOF dt CL",100,0,1.); KMinusPIDCL=new TH1F("KMinusPIDCL","K- PID CL",100,0,1.); DKMinusdEdxFDC=new TH2F("DKMinusdEdxFDC","#delta(dEdx) vs p/M, K- cands.",200,0,24,1000,-25.0,25.0); DKMinusdEdxVsN=new TH2F("DKMinusdEdxVsN","#delta(dEdx) vs N, K- cands.", 28,0.5,28.5,1000,-5.0,5.0); DKMinusdEdx=new TH2F("DKMinusdEdx","#delta(dEdx) vs p/M, K- cands.",200,0,24,1000,-25.0,25.0); KMinusdEdx=new TH2F("KMinusdEdx","dEdx vs p/M, K- cands.",200,0,24,500,0,25.0); KMinusdEdx->SetXTitle("p/M"); KMinusdEdx->SetYTitle("dEdx [keV/cm]"); KMinusdEdxFDC=new TH2F("KMinusdEdxFDC","dEdx vs p/M, K- cands.",200,0,24,500,0,25.0); KMinusdEdxFDC->SetXTitle("p/M"); KMinusdEdxFDC->SetYTitle("dEdx [keV/cm]"); KMinusTOFdt=new TH2F("KMinusTOFdt","#deltat vs p for K-",200,0,10,1000,-20,20); KMinusBCALdt=new TH2F("KMinusBCALdt","#deltat vs p for K-",100,0,5,1000,-20,20); KMinusFCALdt=new TH2F("KMinusFCALdt","#deltat vs p for K-",100,0,5,1000,-20,20); KMinusdEdxCut=new TH2F("KMinusdEdxCut","dEdx vs p/M, K- cands.",200,0,4,500,0,25.0); KMinusdEdxFDCCut=new TH2F("KMinusdEdxFDCCut","dEdx vs p/M, K- cands.",200,0,4,500,0,25.0); KMinusTOFdtCut=new TH2F("KMinusTOFdtCut","#deltat vs p for K-",200,0,10,1000,-20,20); KMinusBCALdtCut=new TH2F("KMinusBCALdtCut","#deltat vs p for K-",100,0,5,1000,-20,20); } gDirectory->cd("../"); gDirectory->mkdir("Proton")->cd(); { ProtondEdxCL=new TH1F("ProtondEdxCL","proton dEdx CL",100,0,1.); ProtondEdxFDCCL=new TH1F("ProtondEdxFDCCL","proton dEdx CL",100,0,1.); ProtonBCALdtCL=new TH1F("ProtonBCALdtCL","proton BCAL dt CL",100,0,1.); ProtonFCALdtCL=new TH1F("ProtonFCALdtCL","proton FCAL dt CL",100,0,1.); ProtonTOFdtCL=new TH1F("ProtonTOFdtCL","proton TOF dt CL",100,0,1.); ProtonPIDCL=new TH1F("ProtonPIDCL","proton PID CL",100,0,1.); ProtondEdxSCCL=new TH1F("ProtondEdxSCCL","proton SC dEdx CL",100,0,1.); DProtondEdxSC=new TH2F("DProtondEdxSC","d(dEdx) vs p/M, proton cands.",200,0,4,500,-25.,25.0); ProtondEdxSC=new TH2F("ProtondEdxSC","dEdx vs p/M, proton cands.",200,0,4, 1000,0,50.0); ProtondEdxSCvsTheta=new TH2F("ProtondEdxSCvsTheta","dEdx vs #theta, proton cands.",180,0,180, 1000,0,50.0); ProtondEdxFDC=new TH2F("ProtondEdxFDC","dEdx vs p/M, proton cands.",200,0,10,500,0,25.0); ProtondEdxFDC->SetXTitle("p/M"); ProtondEdxFDC->SetYTitle("dEdx [keV/cm]"); ProtondEdxFDCwithTOFcut=new TH2F("ProtondEdxFDCwithTOFcut","dEdx vs p/M, proton cands.",200,0,10,500,0,25.0); DProtondEdxVsN=new TH2F("DProtondEdxVsN","#delta(dEdx) vs N, proton cands.", 28,0.5,28.5,1000,-5.0,5.0); ProtondEdxAmp=new TH2F("ProtondEdxAmp","dEdx vs p/M, proton cands.",200,0,4,1000,0,50.0); ProtondEdxVsTheta=new TH2F("ProtondEdxVsTheta","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta1=new TH2F("ProtondEdxVsTheta1","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta2=new TH2F("ProtondEdxVsTheta2","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta3=new TH2F("ProtondEdxVsTheta3","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta4=new TH2F("ProtondEdxVsTheta4","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta5=new TH2F("ProtondEdxVsTheta5","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta6=new TH2F("ProtondEdxVsTheta6","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta7=new TH2F("ProtondEdxVsTheta7","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta8=new TH2F("ProtondEdxVsTheta8","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdxVsTheta9=new TH2F("ProtondEdxVsTheta9","dEdx vs theta",140,0,140,1000,0,50.0); ProtondEdx=new TH2F("ProtondEdx","dEdx vs p/M, proton cands.",200,0,4,500,0,25.0); ProtondEdx->SetXTitle("p/M"); ProtondEdx->SetYTitle("dEdx [keV/cm]"); DProtondEdx=new TH2F("DProtondEdx","#delta(dEdx) vs p/M, proton cands.",200,0,4,1000,-25.0,25.0); DProtondEdxFDC=new TH2F("DProtondEdxFDC","#delta(dEdx) vs p/M, proton cands.",200,0,10,1000,-25.0,25.0); ProtondEdxCut=new TH2F("ProtondEdxCut","dEdx vs p/M, proton cands., cut",200,0,4,1000,0,50.0); ProtondEdxCut->SetXTitle("p/M"); ProtondEdxCut->SetYTitle("dEdx [keV/cm]"); ProtondEdxCutFDC=new TH2F("ProtondEdxCutFDC","dEdx vs p/M, proton cands., cut",200,0,10,500,0,25.0); ProtondEdxCutFDC->SetXTitle("p/M"); ProtondEdxCutFDC->SetYTitle("dEdx [keV/cm]"); ProtonTOFdEdx=new TH2F("ProtonTOFdEdx","dEdx vs p/M, proton cands., cut", 200,0,10,100,0,10.0e-2); ProtonTOFdEdx->SetXTitle("p/M"); ProtonTOFdEdx->SetYTitle("dEdx [GeV/cm]"); ProtonTOFdEdx1=new TH2F("ProtonTOFdEdx1","dEdx vs p/M, proton cands., plane 1", 200,0,10,100,0,10.0e-2); ProtonTOFdEdx1->SetXTitle("p/M"); ProtonTOFdEdx1->SetYTitle("dEdx [GeV/cm]"); \ ProtonTOFdEdx2=new TH2F("ProtonTOFdEdx2","dEdx vs p/M, proton cands., plane 2", 200,0,10,100,0,10.0e-2); ProtonTOFdEdx1->SetXTitle("p/M"); ProtonTOFdEdx1->SetYTitle("dEdx [GeV/cm]"); ProtonTOFdEdxCut=new TH2F("ProtonTOFdEdxCut","dEdx vs p/M, proton cands., cut", 200,0,4,100,0,10.0e-2); ProtonTOFdEdxCut->SetXTitle("p/M"); ProtonTOFdEdxCut->SetYTitle("dEdx [GeV/cm]"); DProtonTOFdEdx1=new TH2F("DProtonTOFdEdx1","#delta(dEdx) vs p/M, proton cands., cut", 200,0,10,100,-5e-2,5.0e-2); DProtonTOFdEdx2=new TH2F("DProtonTOFdEdx2","#delta(dEdx) vs p/M, proton cands., cut", 200,0,10,100,-5e-2,5.0e-2); ProtonTOFdt=new TH2F("ProtonTOFdt","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonTOFdtCut=new TH2F("ProtonTOFdtCut","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonSCdt=new TH2F("ProtonSCdt","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonFCALdt=new TH2F("ProtonFCALdt","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonBCALdt=new TH2F("ProtonBCALdt","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonBCALdtCut=new TH2F("ProtonBCALdtCut","#deltat vs p for proton",100,0,5,1000,-20,20); ProtonBCALdtVsTheta=new TH2F("ProtonBCALdtVsTheta","#deltat vs #theta for proton",180,0,180,1000,-20,20); ProtonThetaVsP=new TH2F("ProtonThetaVsP","#theta vs p for proton",200,0,10,120,0,120); ProtonThetaVsP->SetXTitle("p [GeV/c]"); ProtonThetaVsP->SetYTitle("#theta [degrees]"); ProtonRvsZ=new TH2F("ProtonRvsZ","r vs z for proton POCA to beam line",400,0,200,100,0,20); ProtonRvsZ->SetXTitle("z [cm]"); ProtonRvsZ->SetYTitle("r [cm]"); } gDirectory->cd("../"); gDirectory->mkdir("AntiProton")->cd(); { AntiProtondEdxCL=new TH1F("AntiProtondEdxCL","p-bar dEdx CL",100,0,1.); AntiProtonBCALdtCL=new TH1F("AntiProtonBCALdtCL","p-bar BCAL dt CL",100,0,1.); AntiProtonFCALdtCL=new TH1F("AntiProtonFCALdtCL","p-bar FCAL dt CL",100,0,1.); AntiProtonTOFdtCL=new TH1F("AntiProtonTOFdtCL","p-bar TOF dt CL",100,0,1.); AntiProtonPIDCL=new TH1F("AntiProtonPIDCL","p-bar PID CL",100,0,1.); AntiProtondEdxSCCL=new TH1F("AntiProtondEdxSCCL","anti-proton SC dEdx CL",100,0,1.); AntiProtondEdxFDCCL=new TH1F("AntiProtondEdxFDCCL","p-bar dEdx CL",100,0,1.); DAntiProtondEdxFDC=new TH2F("DAntiProtondEdxFDC","#delta(dEdx) vs p/M, p-bar cands.",200,0,4,1000,-25.0,25.0); DAntiProtondEdxSC=new TH2F("DAntiProtondEdxSC","d(dEdx) vs p/M, antiproton cands.",200,0,4,500,-25.,25.0); DAntiProtondEdx=new TH2F("DAntiProtondEdx","#delta(dEdx) vs p/M, p-bar cands.",200,0,4,1000,-25.0,25.0); AntiProtondEdxSC=new TH2F("AntiProtondEdxSC","dEdx vs p/M, antiproton cands.",200,0,4,500,0,25.0); AntiProtondEdxFDC=new TH2F("AntiProtondEdxFDC","dEdx vs p/M, antiproton cands.",200,0,4,500,0,25.0); AntiProtondEdxFDC->SetXTitle("p/M"); AntiProtondEdxFDC->SetYTitle("dEdx [keV/cm]"); AntiProtondEdxAmp=new TH2F("AntiProtondEdxAmp","dEdx vs p/M, antiproton cands.",200,0,4,500,0,25.0); AntiProtondEdx=new TH2F("AntiProtondEdx","dEdx vs p/M, antiproton cands.",200,0,4,500,0,25.0); AntiProtondEdx->SetXTitle("p/M"); AntiProtondEdx->SetYTitle("dEdx [keV/cm]"); AntiProtondEdxCut=new TH2F("AntiProtondEdxCut","dEdx vs p/M, proton cands., cut",200,0,4,500,0,25.0); AntiProtondEdxCut->SetXTitle("p/M"); AntiProtondEdxCut->SetYTitle("dEdx [keV/cm]"); AntiProtondEdxCutFDC=new TH2F("AntiProtondEdxCutFDC","dEdx vs p/M, proton cands., cut",200,0,4,500,0,25.0); AntiProtondEdxCutFDC->SetXTitle("p/M"); AntiProtondEdxCutFDC->SetYTitle("dEdx [keV/cm]"); AntiProtonTOFdt=new TH2F("AntiProtonTOFdt","#deltat vs p for proton",100,0,5,1000,-20,20); AntiProtonBCALdt=new TH2F("AntiProtonBCALdt","#deltat vs p for proton",100,0,5,1000,-20,20); AntiProtonFCALdt=new TH2F("AntiProtonFCALdt","#deltat vs p for anti-proton",100,0,5,1000,-20,20); AntiProtonThetaVsP=new TH2F("AntiProtonThetaVsP","#theta vs p for proton",200,0,2,90,0,90); AntiProtonThetaVsP->SetXTitle("p [GeV/c]"); AntiProtonThetaVsP->SetYTitle("#theta [degrees]"); AntiProtonRvsZ=new TH2F("AntiProtonRvsZ","r vs z for proton POCA to beam line",400,0,200,100,0,20); AntiProtonRvsZ->SetXTitle("z [cm]"); AntiProtonRvsZ->SetYTitle("r [cm]"); } gDirectory->cd("../"); gDirectory->mkdir("KpKm")->cd(); { HKpKmMass=new TH1F("HKpKmMass",";m(K^{+}K^{-}) [GeV]",1000,0.9,1.9); HKpKmMassCut=new TH1F("HKpKmMassCut",";m(K^{+}K^{-}) [GeV]",1000,0.9,1.9); HKpKmMassKinFitCut=new TH1F("HKpKmMassKinFitCut",";m(K^{+}K^{-}) [GeV]",1000,0.9,1.9); HMissingMass=new TH1F("HMissingMass",";M_X(#gamma d #rightarrow K^{+}K^{-}X) [GeV]", 1000,0,10); KinFitCLMissingDeuteron=new TH1F("KinFitCLMissingDeuteron",";CL",100,0,1); HKpKmMassWithProton=new TH1F("HKpKmMassWithProton",";m(K^{+}K^{-}) [GeV]",1000,0.9,1.9); HKpKmPMissingMass=new TH1F("HKpKmPMissingMass",";M_X(#gamma d #rightarrow K^{+}K^{-}pX) [GeV]", 1000,0,10); } gDirectory->cd("../"); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_PIDCheck::brun(JEventLoop *eventLoop, int32_t runnumber) { DApplication* dapp=dynamic_cast(eventLoop->GetJApplication()); JCalibration * jcalib = dapp->GetJCalibration(runnumber); vector >vals; if (jcalib->Get("FDC/dEdxNdep",vals)==false){ map &row = vals[0]; FDC_N.push_back(row["p1"]); FDC_N.push_back(row["p2"]); FDC_N.push_back(row["p3"]); cout << "FDC dEdx N dependence:"; for (unsigned int i=0;iGet("FDC/dEdxSigma",vals)==false){ for(unsigned int i=0; i &row = vals[i]; switch(int(row["PID"])){ case 8: FDC_P.push_back(row["s1"]); FDC_P.push_back(row["s2"]); FDC_P.push_back(row["s3"]); FDC_P.push_back(row["s4"]); break; default: break; } } */ if (jcalib->Get("CDC/dEdxNdep",vals)==false){ map &row = vals[0]; CDC_N.push_back(row["p1"]); CDC_N.push_back(row["p2"]); CDC_N.push_back(row["p3"]); cout << "CDC dEdx N dependence:"; for (unsigned int i=0;itracks; loop->Get(tracks); if (tracks.size()==0) return NOERROR; vectorbeamphotons; loop->Get(beamphotons); if (beamphotons.size()==0) return RESOURCE_UNAVAILABLE; vectorfcalhits; loop->Get(fcalhits); vectortofpaddles; loop->Get(tofpaddles); //vectoritofhits; //loop->Get(itofhits); vectorfmwpchits; loop->Get(fmwpchits); vectorfmwpcmatches; if (CHECK_FMWPC) loop->Get(fmwpcmatches); vectorctofhits; if (CHECK_FMWPC) loop->Get(ctofhits); vectorctofpoints; if (CHECK_FMWPC) loop->Get(ctofpoints); const DParticleID *mypid_algorithm; loop->GetSingle(mypid_algorithm); // Find t0 (at "vertex") for event and assign list of neutral particles double t0_rf=tracks[0]->Get_BestTrackingFOM()->t0(); // Fill vector of beam photons to pass on the analysis along with weights // for in-time/ out-of-time bool in_time=false; vector>beamPhotons; 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; in_time=true; } if (got_beam_photon){ beamPhotons.push_back(make_pair(beamphotons[i],weight)); } } if (in_time==false) return NOERROR; japp->WriteLock("mylock"); // Setup for performing kinematic fits map>chargedParticles; DKinFitUtils_GlueX *dKinFitUtils = new DKinFitUtils_GlueX(loop); DKinFitter *dKinFitter = new DKinFitter(dKinFitUtils); dKinFitUtils->Reset_NewEvent(); dKinFitter->Reset_NewEvent(); if (tracks.size()==3){ double q1=tracks[0]->Get_Charge(); double q2=tracks[1]->Get_Charge(); double q3=tracks[2]->Get_Charge(); unsigned int ip1=0,ip2=0,in=0; bool got_2p1n=false; if (q1<0 && q2>0 && q3>0){ got_2p1n=true; in=0, ip1=1, ip2=2; } else if (q1>0 && q2<0 && q3>0){ got_2p1n=true; in=1, ip1=0, ip2=2; } else if (q1>0 && q2>0 && q3<0){ got_2p1n=true; in=2, ip1=1, ip2=0; } if (got_2p1n){ const DChargedTrackHypothesis *Kmhyp=tracks[in]->Get_Hypothesis(KMinus); if (Kmhyp!=NULL){ DLorentzVector km=Kmhyp->lorentzMomentum(); const DChargedTrackHypothesis *Kphyp=tracks[ip1]->Get_Hypothesis(KPlus); const DChargedTrackHypothesis *Phyp=tracks[ip2]->Get_Hypothesis(Proton); if (Phyp!=NULL && Kphyp!=NULL){ chargedParticles[KPlus].push_back(Kphyp->Get_TrackTimeBased()); chargedParticles[KMinus].push_back(Kmhyp->Get_TrackTimeBased()); chargedParticles[Proton].push_back(Phyp->Get_TrackTimeBased()); DLorentzVector proton=Phyp->lorentzMomentum(); DLorentzVector kp=Kphyp->lorentzMomentum(); DLorentzVector kp_km=kp+km; HKpKmMassWithProton->Fill(kp_km.M()); for (unsigned int i=0;ilorentzMomentum() +target-kp_km-proton; HKpKmPMissingMass->Fill(missing.M(),weight); } } Kphyp=tracks[ip2]->Get_Hypothesis(KPlus); Phyp=tracks[ip1]->Get_Hypothesis(Proton); if (Phyp!=NULL && Kphyp!=NULL){ chargedParticles.clear(); chargedParticles[KPlus].push_back(Kphyp->Get_TrackTimeBased()); chargedParticles[KMinus].push_back(Kmhyp->Get_TrackTimeBased()); chargedParticles[Proton].push_back(Phyp->Get_TrackTimeBased()); DLorentzVector proton=Phyp->lorentzMomentum(); DLorentzVector kp=Kphyp->lorentzMomentum(); DLorentzVector kp_km=kp+km; HKpKmMassWithProton->Fill(kp_km.M()); for (unsigned int i=0;ilorentzMomentum() +target-kp_km-proton; HKpKmPMissingMass->Fill(missing.M(),weight); } } } } } if (tracks.size()==2){ chargedParticles.clear(); double q1=tracks[0]->Get_Charge(); double q2=tracks[1]->Get_Charge(); if (q1*q2<0){ unsigned int ip=(q1>q2)?0:1; unsigned int in=(q1>q2)?1:0; const DChargedTrackHypothesis *Kphyp=tracks[ip]->Get_Hypothesis(KPlus); const DChargedTrackHypothesis *Kmhyp=tracks[in]->Get_Hypothesis(KMinus); if (Kphyp!=NULL && Kmhyp!=NULL){ DLorentzVector kp=Kphyp->lorentzMomentum(); DLorentzVector km=Kmhyp->lorentzMomentum(); chargedParticles[KPlus].push_back(Kphyp->Get_TrackTimeBased()); chargedParticles[KMinus].push_back(Kmhyp->Get_TrackTimeBased()); DLorentzVector kp_km=kp+km; HKpKmMass->Fill(kp_km.M()); for (unsigned int i=0;ilorentzMomentum() +target-kp_km; HMissingMass->Fill(missing.M(),weight); if (missing.M()>1.5 && missing.M()<2.1){ HKpKmMassCut->Fill(kp_km.M(),weight); if (DoKinematicFit(t0_rf,beamPhotons[i].first,chargedParticles, dKinFitUtils,dKinFitter,Deuteron)){ double CL=dKinFitter->Get_ConfidenceLevel(); KinFitCLMissingDeuteron->Fill(CL,weight); if (CL>0.01){ HKpKmMassKinFitCut->Fill(kp_km.M(),weight); } } } } } } } /* for (unsigned int i=0;iFill(itofhits[i]->x,itofhits[i]->y); } */ for (unsigned int i=0;ibar << " " << ctofhits[i]->end << " " // << ctofhits[i]->dE << " " << ctofhits[i]->t << endl; int index=2.*ctofhits[i]->bar-1+ctofhits[i]->end; CtofTs->Fill(ctofhits[i]->t,index); CtofdEs->Fill(1000.*ctofhits[i]->dE,index); } for (unsigned int i=0;iFill(point->t,point->bar); CtofdE->Fill(1000*point->dE,point->bar); } if (CHECK_FMWPC){ for (unsigned int j=0;jtbt << " " // << fmwpcmatches[j]->FMWPC_closest_wire[k] << endl; if (fmwpcmatches[j]->FMWPC_closest_wire[k]>0){ //cout << fmwpcmatches[j]->FMWPC_dist_closest_wire[k] << endl; } } } bool got_x=false; bool got_y=false; double x=0,y=0; for (unsigned int j=0;jlayer; if (mylayer==1){ if (got_x==false){ x=-71.5+(fmwpchits[j]->wire-1); got_x=true; } } if (mylayer==2){ if (got_y==false){ y=-71.5+(fmwpchits[j]->wire-1); got_y=true; } } } if (got_x && got_y){ FMWPCxy->Fill(x,y); } } double my_t0=-1000.; /* if (tracks.size()>0){ const DChargedTrackHypothesis *hyp=tracks[0]->Get_BestTrackingFOM(); if (hyp){ DVector3 pos=hyp->position(); shared_ptrmyErrorMatrix=hyp->errorMatrix(); pos.Print(); myErrorMatrix->Print(); my_t0=hyp->t0(); } } */ if (tracks.size()==1){ const DChargedTrackHypothesis *hyp=tracks[0]->Get_Hypothesis(Proton); if (hyp){ double t0=hyp->t0(); DVector3 pos=hyp->position(); if (pos.z()>50.&&pos.z()<80 && pos.Perp()<1.){ for (unsigned int i=0;ix,fcalhits[i]->y,625.); double dt=fcalhits[i]->t-(fpos-pos).Mag()/29.98-t0; FCALdt_vs_E->Fill(fcalhits[i]->E,dt); } } } } Particle_t myPIDS[8]={Electron,Positron,PiPlus,PiMinus,KPlus,KMinus,Proton, AntiProton}; for (unsigned int j=0;jGet_Hypothesis(PiPlus); if (hyp!=NULL){ const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double p=track->momentum().Mag(); shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); shared_ptrfcalhitparms=hyp->Get_FCALSingleHitMatchParams(); shared_ptrctofparms=hyp->Get_CTOFHitMatchParams(); double t0=hyp->t0(); if (ctofparms!=NULL){ double t=ctofparms->dCTOFPoint->t; double tflight=ctofparms->dFlightTime; double dt=t-tflight-t0; PiPlusCTOFdt->Fill(dt); PiPlusCTOFMatchdXdY->Fill(ctofparms->dDeltaXToHit, ctofparms->dDeltaYToHit); } if (fcalhitparms!=NULL){ PiPlus_MatchFCALHit_d_vs_p->Fill(p,fcalhitparms->dDOCAToHit); PiPlus_MatchFCALHit_E_vs_p->Fill(p,fcalhitparms->dEHit); PiPlus_MatchFCALHit_E_over_p_vs_p->Fill(p,fcalhitparms->dEHit/p); double dt=fcalhitparms->dTHit-fcalhitparms->dFlightTime-t0; PiPlus_MatchFCALHit_dt_vs_p->Fill(p,dt); PiPlus_MatchFCALHit_dt_vs_E->Fill(fcalhitparms->dEHit,dt); } if (CHECK_FMWPC){ vectorfmwpc_extraps=track->extrapolations.at(SYS_FMWPC); if (fmwpc_extraps.size()>0){ DVector3 pos=fmwpc_extraps[0].position; ExtrapToFMWPC_PiPlus->Fill(pos.x(),pos.y()); if (fmwpchits.size()>0){ //cout << fmwpc_extraps.size() << endl; ExtrapToFMWPC_PiPlus_with_FMWPC->Fill(pos.x(),pos.y()); for (unsigned int j=0;jlayer; FMWPCTime_PiPlus->Fill(mylayer,fmwpchits[j]->t); FMWPCQ_PiPlus->Fill(mylayer,fmwpchits[j]->q); FMWPCAmp_PiPlus->Fill(mylayer,fmwpchits[j]->q); if (mylayer<=fmwpc_extraps.size()){ pos=fmwpc_extraps[mylayer-1].position; //cout << " layer " << mylayer << endl; //pos.Print(); if (mylayer==1 || mylayer==3 || mylayer==5){ double dx=-72.644+1.016*(fmwpchits[j]->wire-1)-pos.x(); //cout << " dx " << dx << endl; FMWPCdXY_vs_layer_Pip->Fill(mylayer,dx); } else{ double dy=-72.644+1.016*(fmwpchits[j]->wire-1)-pos.y(); //cout << " dy " << dy << endl; FMWPCdXY_vs_layer_Pip->Fill(mylayer,dy); } if (mylayer==1){ double dx=-72.644+1.016*(fmwpchits[j]->wire-1)-pos.x(); FMWPCdX_vs_p_Pip->Fill(fmwpc_extraps[0].momentum.Mag(),dx); } if (mylayer==2){ double dy=-72.644+1.016*(fmwpchits[j]->wire-1)-pos.y(); FMWPCdY_vs_p_Pip->Fill(fmwpc_extraps[1].momentum.Mag(),dy); } } } } } } if (CHECK_TOF&&p>1.&&track->extrapolations.size()>0){ vectortof_extraps=track->extrapolations.at(SYS_TOF); if (tof_extraps.size()>0){ DVector3 pos=tof_extraps[0].position; vectorverti; vectorhoriz; for (unsigned int k=0;korientation==0) verti.push_back(tofpaddles[k]); else horiz.push_back(tofpaddles[k]); } if (horiz.size()==0 && verti.size()>0){ ExtrapToTOF_NoHoriz->Fill(pos.x(),pos.y()); } if (horiz.size()>0 && verti.size()==0){ ExtrapToTOF_NoVerti->Fill(pos.x(),pos.y()); } if (tofparms==NULL && fcalparms!=NULL){ ExtrapToTOF_PiPlus->Fill(pos.x(),pos.y()); } } } } if ((hyp=tracks[j]->Get_Hypothesis(PiMinus))!=NULL){ const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double p=track->momentum().Mag(); if (CHECK_FMWPC){ vectorctof_extraps=track->extrapolations.at(SYS_CTOF); if (ctof_extraps.size()>0){ double projx=ctof_extraps[0].position.x(); double projy=ctof_extraps[0].position.y(); CtofExtrapXY->Fill(projx,projy); if (ctofpoints.size()>0&&ctof_extraps[0].momentum.Mag()>1.){ CtofExtrapXYwithHit->Fill(projx,projy); } for (unsigned int i=0;ipos.y()-projy; double dx=point->pos.x()-projx; double dt=point->t-ctof_extraps[0].t-my_t0; Ctofdy->Fill(dy,point->bar); Ctofdxdy->Fill(dx,dy); CtofYCorr->Fill(projy,point->pos.y()); Ctofdt->Fill(dt,point->bar); } } vectorfmwpc_extraps=track->extrapolations.at(SYS_FMWPC); if (fmwpc_extraps.size()>0){ DVector3 pos=fmwpc_extraps[0].position; ExtrapToFMWPC_PiMinus->Fill(pos.x(),pos.y()); if (fmwpchits.size()>0){ ExtrapToFMWPC_PiMinus_with_FMWPC->Fill(pos.x(),pos.y()); for (unsigned int j=0;jlayer; FMWPCTime_PiMinus->Fill(mylayer,fmwpchits[j]->t); FMWPCQ_PiMinus->Fill(mylayer,fmwpchits[j]->q); FMWPCAmp_PiMinus->Fill(mylayer,fmwpchits[j]->q); if (mylayer<=fmwpc_extraps.size()){ // cout << mylayer << " x/y " << -72.644+1.016*double(fmwpchits[j]->wire-1) << " x " << pos.x()<<" Y " << pos.y() << endl; if (mylayer==1 || mylayer==3 || mylayer==5 ){ double dx=-72.644+1.016*double(fmwpchits[j]->wire-1)-pos.x(); FMWPCdXY_vs_layer_Pim->Fill(mylayer,dx); } else{ double dy=-72.644+1.016*double(fmwpchits[j]->wire-1)-pos.y(); FMWPCdXY_vs_layer_Pim->Fill(mylayer,dy); } if (mylayer==1){ double dx=-72.644+1.016*double(fmwpchits[j]->wire-1)-pos.x(); // cout << "x " << -71.5+double(fmwpchits[j]->wire-1) << " " << pos.x() << endl; FMWPCdX_vs_p_Pim->Fill(fmwpc_extraps[0].momentum.Mag(),dx); } if (mylayer==2){ double dy=-72.644+1.016*double(fmwpchits[j]->wire-1)-pos.y(); FMWPCdY_vs_p_Pim->Fill(fmwpc_extraps[1].momentum.Mag(),dy); } } } } } } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); shared_ptrfcalhitparms=hyp->Get_FCALSingleHitMatchParams(); shared_ptrctofparms=hyp->Get_CTOFHitMatchParams(); double t0=hyp->t0(); if (ctofparms!=NULL){ double t=ctofparms->dCTOFPoint->t; double tflight=ctofparms->dFlightTime; double dt=t-tflight-t0; PiMinusCTOFdt->Fill(dt); PiMinusCTOFMatchdXdY->Fill(ctofparms->dDeltaXToHit, ctofparms->dDeltaYToHit); } if (fcalhitparms!=NULL){ PiMinus_MatchFCALHit_d_vs_p->Fill(p,fcalhitparms->dDOCAToHit); PiMinus_MatchFCALHit_E_vs_p->Fill(p,fcalhitparms->dEHit); PiMinus_MatchFCALHit_E_over_p_vs_p->Fill(p,fcalhitparms->dEHit/p); double dt=fcalhitparms->dTHit-fcalhitparms->dFlightTime-t0; PiMinus_MatchFCALHit_dt_vs_p->Fill(p,dt); PiMinus_MatchFCALHit_dt_vs_E->Fill(fcalhitparms->dEHit,dt); } if (tofparms==NULL && fcalparms!=NULL){ if (CHECK_TOF&&p>1.&&track->extrapolations.size()>0){ vectortof_extraps=track->extrapolations.at(SYS_TOF); if (tof_extraps.size()>0){ DVector3 pos=tof_extraps[0].position; ExtrapToTOF_PiMinus->Fill(pos.x(),pos.y()); } } } } if ((hyp=tracks[j]->Get_Hypothesis(Positron))!=NULL){ const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); double theta=track->momentum().Theta()*180/M_PI; if (false && theta>15.&& bcalparms!=NULL && fcalparms!=NULL){ cout <<"Event: "<< eventnumber << " positron theta " << theta << endl; cout << bcalparms->dBCALShower->x<< " " << bcalparms->dBCALShower->y << " " << bcalparms->dBCALShower->z << endl; fcalparms->dFCALShower->getPosition().Print(); cout << " doca to fcal " << fcalparms->dDOCAToShower << endl; } } if ((hyp=tracks[j]->Get_Hypothesis(Electron))!=NULL){ const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); double theta=track->momentum().Theta()*180/M_PI; if (false && theta>15 && bcalparms!=NULL && fcalparms!=NULL){ cout << "Event: " <dBCALShower->x<< " " << bcalparms->dBCALShower->y << " " << bcalparms->dBCALShower->z << endl; fcalparms->dFCALShower->getPosition().Print(); cout << " doca to fcal " << fcalparms->dDOCAToShower << endl; } } vector >probabilities; for (unsigned int i=0;i<8;i++){ if ((hyp=tracks[j]->Get_Hypothesis(myPIDS[i]))!=NULL){ probabilities.push_back(make_pair(hyp->Get_FOM(),myPIDS[i])); FillPIDHistos(mypid_algorithm,hyp); } } if (probabilities.size()>0){ sort(probabilities.begin(),probabilities.end(),SortParticleProbability); switch(probabilities[0].second){ case Positron: PositronPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>POSITRON_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(Positron); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double p=track->momentum().Mag(); double t0=hyp->t0(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); // shared_ptritofparms=hyp->Get_ITOFHitMatchParams(); /* if (itofparms!=NULL){ double t=itofparms->dITOFHit->t; double tflight=itofparms->dFlightTime; double dt=t-tflight-t0; PositronITOFdt->Fill(p,dt); PositronITOFdXdY->Fill(itofparms->dDeltaXToHit,itofparms->dDeltaYToHit); } */ if (bcalparms!=NULL){ double E=bcalparms->dBCALShower->E; EOverPPlus_BCAL_cut->Fill(p,E/p); } if (fcalparms!=NULL){ double E=fcalparms->dFCALShower->getEnergy(); EOverPPlus_cut->Fill(p,E/p); /* if (track->extrapolations.size()>0){ vectorfcal_extraps=track->extrapolations.at(SYS_FCAL); if (fcal_extraps.size()>0){ DVector3 fcal_pos=fcalparms->dFCALShower->getPosition(); DVector3 track_pos=fcal_extraps[0].position; DVector3 track_mom=fcal_extraps[0].momentum; track_pos+=((fcal_pos.z()-track_pos.z())/track_mom.z())*track_mom; //cout <<"-----------"<Fill(track_mom.Mag(),diff.x()); PositronYdiff_vs_p->Fill(track_mom.Mag(),diff.y()); PositronDiff_vs_theta->Fill(180./M_PI*track_mom.Theta(), diff.Perp()); } } */ } } break; case Electron: ElectronPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>ELECTRON_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(Electron); double t0=hyp->t0(); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double p=track->momentum().Mag(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); /* shared_ptritofparms=hyp->Get_ITOFHitMatchParams(); if (itofparms!=NULL){ double t=itofparms->dITOFHit->t; double tflight=itofparms->dFlightTime; double dt=t-tflight-t0; ElectronITOFdt->Fill(p,dt); ElectronITOFdXdY->Fill(itofparms->dDeltaXToHit,itofparms->dDeltaYToHit); } */ double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; double betagamma=p/track->mass(); if (dEdx>0.){ ElectrondEdxCut->Fill(betagamma,dEdx); } if (dEdx_fdc>0){ ElectrondEdxFDCCut->Fill(betagamma,dEdx_fdc); } if (bcalparms!=NULL){ double E=bcalparms->dBCALShower->E; EOverPMinus_BCAL_cut->Fill(p,E/p); /* double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-hyp->t0()-pid_algorithm->GetTimeMean(SYS_BCAL,Positron,p); ElectronBCALdtCut->Fill(p,dt_bcal); */ } if (fcalparms!=NULL){ double E=fcalparms->dFCALShower->getEnergy(); EOverPMinus_cut->Fill(p,E/p); /* if (track->extrapolations.size()>0){ vectorfcal_extraps=track->extrapolations.at(SYS_FCAL); if (fcal_extraps.size()>0){ DVector3 fcal_pos=fcalparms->dFCALShower->getPosition(); DVector3 track_pos=fcal_extraps[0].position; DVector3 track_mom=fcal_extraps[0].momentum; track_pos+=((fcal_pos.z()-track_pos.z())/track_mom.z())*track_mom; //cout <<"-----------"<Fill(track_mom.Mag(),diff.x()); ElectronYdiff_vs_p->Fill(track_mom.Mag(),diff.y()); ElectronDiff_vs_theta->Fill(180./M_PI*track_mom.Theta(), diff.Perp()); } } */ } } break; case PiPlus: PiPlusPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>PIP_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(PiPlus); double t0=hyp->t0(); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; DVector3 mom=track->momentum(); DVector3 pos=track->position(); double p=mom.Mag(); double betagamma=mom.Mag()/track->mass(); PipRvsZ->Fill(pos.z(),pos.Perp()); if (dEdx>0.){ PiPlusdEdxCut->Fill(betagamma,dEdx); } if (dEdx_fdc>0){ PiPlusdEdxFDCCut->Fill(betagamma,dEdx_fdc); } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); /* shared_ptritofparms=hyp->Get_ITOFHitMatchParams(); if (itofparms!=NULL){ double t=itofparms->dITOFHit->t; double tflight=itofparms->dFlightTime; double dt=t-tflight-t0; PiPlusITOFdt->Fill(p,dt); PipITOFdXdY->Fill(itofparms->dDeltaXToHit,itofparms->dDeltaYToHit); } */ if (bcalparms!=NULL){ /* double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-hyp->t0()-pid_algorithm->GetTimeMean(SYS_BCAL,Positron,p); PiPlusBCALdtCut->Fill(p,dt_bcal); */ } if (tofparms!=NULL){ //cout << " TOF match " << endl; if (tofparms->dEdx1>0.){ PiPlusTOFdEdxCut->Fill(betagamma,tofparms->dEdx1); } if (tofparms->dEdx2>0.){ PiPlusTOFdEdxCut->Fill(betagamma,tofparms->dEdx2); } double dt=tofparms->dHitTime-tofparms->dFlightTime-hyp->t0(); PiPlusTOFdtCut->Fill(p,dt); } } break; case PiMinus: PiMinusPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>PIM_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(PiMinus); double t0=hyp->t0(); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; DVector3 mom=track->momentum(); double p=mom.Mag(); DVector3 pos=track->position(); double betagamma=mom.Mag()/track->mass(); PimRvsZ->Fill(pos.z(),pos.Perp()); if (dEdx>0.){ PiMinusdEdxCut->Fill(betagamma,dEdx); } if (dEdx_fdc>0){ PiMinusdEdxFDCCut->Fill(betagamma,dEdx_fdc); } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); /* shared_ptritofparms=hyp->Get_ITOFHitMatchParams(); if (itofparms!=NULL){ double t=itofparms->dITOFHit->t; double tflight=itofparms->dFlightTime; double dt=t-tflight-t0; PiMinusITOFdt->Fill(p,dt); PimITOFdXdY->Fill(itofparms->dDeltaXToHit,itofparms->dDeltaYToHit); } */ if (tofparms!=NULL){ if (tofparms->dEdx1>0.){ PiMinusTOFdEdxCut->Fill(betagamma,tofparms->dEdx1); } if (tofparms->dEdx2>0.){ PiMinusTOFdEdxCut->Fill(betagamma,tofparms->dEdx2); } } } break; case KPlus: KPlusPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>KP_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(KPlus); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); DVector3 mom=track->momentum(); DVector3 pos=track->position(); double p=mom.Mag(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; double betagamma=mom.Mag()/track->mass(); if (dEdx_fdc>0){ KPlusdEdxFDCCut->Fill(betagamma,dEdx_fdc); } if (dEdx>0.){ KPlusdEdxCut->Fill(betagamma,dEdx); } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); if (bcalparms!=NULL){ /* double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-hyp->t0()-pid_algorithm->GetTimeMean(SYS_BCAL,Positron,p); KPlusBCALdtCut->Fill(p,dt_bcal); */ } if (tofparms!=NULL){ double dt=tofparms->dHitTime-tofparms->dFlightTime-hyp->t0(); KPlusTOFdtCut->Fill(p,dt); } } break; case KMinus: KMinusPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>PIM_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(KMinus); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); DVector3 mom=track->momentum(); DVector3 pos=track->position(); double p=mom.Mag(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; double betagamma=mom.Mag()/track->mass(); if (dEdx_fdc>0){ KMinusdEdxFDCCut->Fill(betagamma,dEdx_fdc); } if (dEdx>0.){ KMinusdEdxCut->Fill(betagamma,dEdx); } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); if (bcalparms!=NULL){ /* double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-hyp->t0()-pid_algorithm->GetTimeMean(SYS_BCAL,Positron,p); KMinusBCALdtCut->Fill(p,dt_bcal); */ } if (tofparms!=NULL){ double dt=tofparms->dHitTime-tofparms->dFlightTime-hyp->t0(); KMinusTOFdtCut->Fill(p,dt); } } break; case Proton: ProtonPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>PROTON_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(Proton); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); DVector3 mom=track->momentum(); double p=mom.Mag(); DVector3 pos=track->position(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; double betagamma=mom.Mag()/track->mass(); ProtonRvsZ->Fill(pos.z(),pos.Perp()); ProtonThetaVsP->Fill(mom.Mag(),180./M_PI*mom.Theta()); if (dEdx_fdc>0.){ ProtondEdxCutFDC->Fill(betagamma,dEdx_fdc); } if (dEdx>0.){ ProtondEdxCut->Fill(betagamma,dEdx); } shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); if (bcalparms!=NULL){ /* double dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-hyp->t0()-pid_algorithm->GetTimeMean(SYS_BCAL,Positron,p); ProtonBCALdtCut->Fill(p,dt_bcal); */ } if (tofparms!=NULL){ if (tofparms->dEdx1>0.){ ProtonTOFdEdx1->Fill(betagamma,tofparms->dEdx1); } if (tofparms->dEdx2>0.){ ProtonTOFdEdx2->Fill(betagamma,tofparms->dEdx2); } double dt_tof=tofparms->dHitTime-tofparms->dFlightTime-hyp->t0(); ProtonTOFdtCut->Fill(p,dt_tof); } } break; case AntiProton: AntiProtonPIDCL->Fill(probabilities[0].first); if (probabilities[0].first>ANTIPROTON_CL_CUT){ hyp=tracks[j]->Get_Hypothesis(AntiProton); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double dEdx=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; DVector3 mom=track->momentum(); DVector3 pos=track->position(); double betagamma=mom.Mag()/track->mass(); AntiProtonRvsZ->Fill(pos.z(),pos.Perp()); AntiProtonThetaVsP->Fill(mom.Mag(),180./M_PI*mom.Theta()); if (dEdx>0.) AntiProtondEdxCut->Fill(betagamma,dEdx); if (dEdx_fdc>0.) AntiProtondEdxCutFDC->Fill(betagamma,dEdx_fdc); } break; default: break; } } } delete dKinFitter; delete dKinFitUtils; japp->Unlock("mylock"); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_PIDCheck::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_PIDCheck::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } void JEventProcessor_PIDCheck::FillPIDHistos(const DParticleID* mypid_algorithm,const DChargedTrackHypothesis *hyp) const { Particle_t myPID=hyp->PID(); double t0_rf=hyp->t0(); const DTrackTimeBased *track=hyp->Get_TrackTimeBased(); double p=track->momentum().Mag(); double theta=180./M_PI*track->momentum().Theta(); DVector3 pos=track->position(); double z=pos.z(); shared_ptrtofparms=hyp->Get_TOFHitMatchParams(); shared_ptrbcalparms=hyp->Get_BCALShowerMatchParams(); shared_ptrfcalparms=hyp->Get_FCALShowerMatchParams(); shared_ptrscparms=hyp->Get_SCHitMatchParams(); double dt=1000.,dt_bcal=1000.; double betagamma=p/track->mass(); double beta=p/track->energy(); // double bg2=betagamma*betagamma; double chisq_dEdx=0.; double chisq_bcal=0.; double chisq_tof=0.; // printf("t1 detector %d time %f\n",hyp->t1_detector(),hyp->t1()); double dEdx=1e6*track->ddEdx_CDC; //_DBG_ << dEdx << endl; double dEdxAmp=1e6*hyp->Get_dEdx_CDC_amp(); double dEdx_fdc=1e6*track->ddEdx_FDC; unsigned int n_dedx=track->dNumHitsUsedFordEdx_CDC; unsigned int n_dedx_fdc=track->dNumHitsUsedFordEdx_FDC; if (scparms!=NULL){ //printf("SC\n"); double dEdxSC=1e3*scparms->dEdx; // double dEdxmean=1.25/bg2+0.5/betagamma+1.5+0.1*betagamma;//+0.5; double dEdxmean=1e3*mypid_algorithm->GetProtondEdxMean_SC(beta); double diff=dEdxSC-dEdxmean; //double dEdxsig=0.3/bg2-0.45/betagamma+0.75; double dEdxsig=1e3*mypid_algorithm->GetProtondEdxSigma_SC(beta); double chisq=diff*diff/(dEdxsig*dEdxsig); double scdt = scparms->dHitTime-scparms->dFlightTime-t0_rf; switch(myPID){ case PiPlus: PiPlusSCdt->Fill(p,scdt); break; case Proton: ProtondEdxSC->Fill(betagamma,dEdxSC); ProtondEdxSCvsTheta->Fill(theta,dEdxSC); DProtondEdxSC->Fill(betagamma,diff); ProtondEdxSCCL->Fill(TMath::Prob(chisq,1.)); ProtonSCdt->Fill(p,scdt); break; case AntiProton: AntiProtondEdxSC->Fill(betagamma,dEdxSC); DAntiProtondEdxSC->Fill(betagamma,diff); AntiProtondEdxSCCL->Fill(TMath::Prob(chisq,1.)); break; default: break; } } if (dEdx>0){ double dedx_mean=0.; if (mypid_algorithm->GetdEdxMean_CDC(beta,0,dedx_mean,myPID)==NOERROR){ double dedx_diff=dEdxAmp-1e6*dedx_mean; double dedx_sigma=0.; // printf("%f ", dedx_sigma); mypid_algorithm->GetdEdxSigma_CDC(beta,n_dedx,dedx_sigma,myPID); dedx_sigma*=1e6; // printf("%f\n", dedx_sigma); chisq_dEdx=dedx_diff*dedx_diff/(dedx_sigma*dedx_sigma); double nscale=CDC_N[0]*pow(double(n_dedx),CDC_N[1])+CDC_N[2]; switch(myPID){ case Positron: PositrondEdx->Fill(betagamma,dEdxAmp); DPositrondEdx->Fill(betagamma,dedx_diff); PositrondEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case Electron: ElectrondEdx->Fill(betagamma,dEdxAmp); DElectrondEdx->Fill(betagamma,dedx_diff); ElectrondEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case PiPlus: PiPlusdEdx->Fill(betagamma,dEdxAmp); DPiPlusdEdx->Fill(betagamma,dedx_diff); PiPlusdEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case PiMinus: { int index=int(floor(p/0.25)); PiMinusdEdx->Fill(betagamma,dEdxAmp); DPiMinusdEdxVsNall->Fill(n_dedx,dedx_diff/(dedx_sigma/nscale)); if (index<40){ DPiMinusdEdxVsN[index]->Fill(n_dedx,dedx_diff); } DPiMinusdEdxVsNcorr->Fill(n_dedx,dedx_diff/dedx_sigma); DPiMinusdEdx->Fill(betagamma,dedx_diff); PiMinusdEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); } break; case KPlus: KPlusdEdx->Fill(betagamma,dEdxAmp); DKPlusdEdx->Fill(betagamma,dedx_diff); KPlusdEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case KMinus: KMinusdEdx->Fill(betagamma,dEdxAmp); DKMinusdEdx->Fill(betagamma,dedx_diff); DKMinusdEdxVsN->Fill(n_dedx,dedx_diff/dedx_sigma); KMinusdEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case Proton: ProtondEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); ProtondEdx->Fill(betagamma,dEdx); //if (z>0&&z<2) { if (betagamma>0.35 && betagamma<0.4&&n_dedx>14){ ProtondEdxVsTheta->Fill(theta,dEdxAmp); } if (betagamma>0.4 && betagamma<0.45&&n_dedx>14){ ProtondEdxVsTheta1->Fill(theta,dEdxAmp); } if (betagamma>0.45 && betagamma<0.5&&n_dedx>14){ ProtondEdxVsTheta2->Fill(theta,dEdxAmp); } if (betagamma>0.5 && betagamma<0.55&&n_dedx>14){ ProtondEdxVsTheta3->Fill(theta,dEdxAmp); } if (betagamma>0.55 && betagamma<0.6&&n_dedx>14){ ProtondEdxVsTheta4->Fill(theta,dEdxAmp); } if (betagamma>0.6 && betagamma<0.65&&n_dedx>14){ ProtondEdxVsTheta5->Fill(theta,dEdxAmp); } if (betagamma>0.65 && betagamma<0.7&&n_dedx>14){ ProtondEdxVsTheta6->Fill(theta,dEdxAmp); } if (betagamma>0.7 && betagamma<0.75&&n_dedx>14){ ProtondEdxVsTheta7->Fill(theta,dEdxAmp); } if (betagamma>0.75 && betagamma<0.8&&n_dedx>14){ ProtondEdxVsTheta8->Fill(theta,dEdxAmp); } if (betagamma>0.8 && betagamma<0.85&&n_dedx>14){ ProtondEdxVsTheta9->Fill(theta,dEdxAmp); } DProtondEdxVsN->Fill(n_dedx,dedx_diff/dedx_sigma); DProtondEdx->Fill(betagamma,dedx_diff); ProtondEdxAmp->Fill(betagamma,dEdxAmp); } break; case AntiProton: AntiProtondEdx->Fill(betagamma,dEdx); DAntiProtondEdx->Fill(betagamma,dedx_diff); AntiProtondEdxAmp->Fill(betagamma,dEdxAmp); AntiProtondEdxCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; default: break; } } } if (dEdx_fdc>0.){ double dedx_mean=0.; if (mypid_algorithm->GetdEdxMean_FDC(beta,0,dedx_mean,myPID)==NOERROR){ double dedx_diff=dEdx_fdc-1e6*dedx_mean; double dedx_sigma=0.; // printf("%f ", dedx_sigma); mypid_algorithm->GetdEdxSigma_FDC(beta,n_dedx_fdc,dedx_sigma,myPID); dedx_sigma*=1e6; // printf("%f\n", dedx_sigma); chisq_dEdx=dedx_diff*dedx_diff/(dedx_sigma*dedx_sigma); double nscale=FDC_N[0]*pow(double(n_dedx_fdc),FDC_N[1])+FDC_N[2]; /* double nscale2=FDC_P[0]/(betagamma*betagamm)+FDC_P[1]/betagamma +FDC_P[2]+FDC_P[3]*betagamma; */ switch(myPID){ case Positron: PositrondEdxFDC->Fill(betagamma,dEdx_fdc); DPositrondEdxFDC->Fill(betagamma,dedx_diff); PositrondEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case Electron: ElectrondEdxFDC->Fill(betagamma,dEdx_fdc); DElectrondEdxFDC->Fill(betagamma,dedx_diff); ElectrondEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case PiPlus: PiPlusdEdxFDC->Fill(betagamma,dEdx_fdc); DPiPlusdEdxFDC->Fill(betagamma,dedx_diff); PiPlusdEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case PiMinus: { int index=int(floor(p/0.25)); PiMinusdEdxFDC->Fill(betagamma,dEdx_fdc); DPiMinusdEdxFDC->Fill(betagamma,dedx_diff); DPiMinusdEdxFDCVsNall->Fill(n_dedx_fdc,dedx_diff/(dedx_sigma/nscale)); if (index>=0&&index<40){ DPiMinusdEdxFDCVsN[index]->Fill(n_dedx_fdc,dedx_diff); } DPiMinusdEdxFDCVsNcorr->Fill(n_dedx_fdc,dedx_diff/(dedx_sigma)); PiMinusdEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); } break; case KPlus: KPlusdEdxFDC->Fill(betagamma,dEdx_fdc); DKPlusdEdxFDC->Fill(betagamma,dedx_diff); KPlusdEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case KMinus: KMinusdEdxFDC->Fill(betagamma,dEdx_fdc); DKMinusdEdxFDC->Fill(betagamma,dedx_diff); KMinusdEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; case Proton: ProtondEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); ProtondEdxFDC->Fill(betagamma,dEdx_fdc); DProtondEdxFDC->Fill(betagamma,dedx_diff); break; case AntiProton: AntiProtondEdxFDC->Fill(betagamma,dEdx_fdc); DAntiProtondEdxFDC->Fill(betagamma,dedx_diff); AntiProtondEdxFDCCL->Fill(TMath::Prob(chisq_dEdx,1.)); break; default: break; } } } if (bcalparms!=NULL){ //printf("BCAL\n"); dt_bcal=bcalparms->dBCALShower->t-bcalparms->dFlightTime-t0_rf; //dt_bcal-=pid_algorithm->GetTimeMean(SYS_BCAL,myPID,p); ////cout << pid_algorithm->GetTimeMean(SYS_BCAL,myPID,p) <dBCALShower->E/p;//bcalparms->dMomentum; double E_over_p_mean=mypid_algorithm->GetEOverPMean(SYS_BCAL,p);//bcalparms->dMomentum); double E_over_p_diff=E_over_p-E_over_p_mean; double E_over_p_sigma=mypid_algorithm->GetEOverPSigma(SYS_BCAL,p);//bcalparms->dMomentum); double chisq_E_over_p=E_over_p_diff*E_over_p_diff/(E_over_p_sigma*E_over_p_sigma); //double chisq_E_over_p=hyp->Get_ChiSq_EoverP(); // cout << chisq_E_over_p << endl; switch(myPID){ case Positron: PositronBCALdt->Fill(p,dt_bcal); vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,Positron,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; PositronBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); EOverPPlus_BCAL->Fill(p/*bcalparms->dMomentum*/,E_over_p); DEOverPPlus_BCAL->Fill(p/*bcalparms->dMomentum*/,E_over_p_diff); PositronBCALEOverPCL->Fill(TMath::Prob(chisq_E_over_p,1.)); break; case Electron: ElectronBCALdt->Fill(p,dt_bcal); vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,Electron,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; ElectronBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); EOverPMinus_BCAL->Fill(p/*bcalparms->dMomentum*/,E_over_p); DEOverPMinus_BCAL->Fill(p/*bcalparms->dMomentum*/,E_over_p_diff); ElectronBCALEOverPCL->Fill(TMath::Prob(chisq_E_over_p,1.)); break; case PiPlus: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,PiPlus,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; PiPlusBCALdt->Fill(p,dt_bcal); PiPlusBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); break; case PiMinus: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,PiMinus,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; PiMinusBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); PiMinusBCALdt->Fill(p,dt_bcal); PiMinusBCALdtVsTheta->Fill(theta,dt_bcal); break; case KPlus: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,KPlus,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; KPlusBCALdt->Fill(p,dt_bcal); KPlusBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); break; case KMinus: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,KMinus,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; KMinusBCALdt->Fill(p,dt_bcal); KMinusBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); break; case Proton: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,Proton,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; ProtonBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); ProtonBCALdt->Fill(p,dt_bcal); ProtonBCALdtVsTheta->Fill(theta,dt_bcal); break; case AntiProton: vart_bcal=mypid_algorithm->GetTimeVariance(SYS_BCAL,AntiProton,p); chisq_bcal=dt_bcal*dt_bcal/vart_bcal; AntiProtonBCALdt->Fill(p,dt_bcal); AntiProtonBCALdtCL->Fill(TMath::Prob(chisq_bcal,1.)); break; default: break; } } if (tofparms!=NULL){ double betagamma_tof=p/*tofparms->dMomentum*//track->mass(); //printf("TOF\n"); dt=tofparms->dHitTime-tofparms->dFlightTime-t0_rf; // double sigt_tof=0.01988/p+0.04664; double vart_tof=0.; chisq_tof=0.; switch(myPID){ case Positron: PositronTOF->Fill(p,dt); vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,Positron,p); chisq_tof=dt*dt/vart_tof; PositronTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case Electron: ElectronTOF->Fill(p,dt); vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,Electron,p); chisq_tof=dt*dt/vart_tof; ElectronTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case PiPlus: PiPlusTOFdt->Fill(p,dt); if (tofparms->dEdx1>0.){ PiPlusTOFdEdx->Fill(betagamma_tof,tofparms->dEdx1); } if (tofparms->dEdx2>0.){ PiPlusTOFdEdx->Fill(betagamma_tof,tofparms->dEdx2); } vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,PiPlus,p); chisq_tof=dt*dt/vart_tof; PiPlusTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case PiMinus: PiMinusTOFdt->Fill(p,dt); if (tofparms->dEdx1>0.){ PiMinusTOFdEdx->Fill(betagamma_tof,tofparms->dEdx1); } if (tofparms->dEdx2>0.){ PiMinusTOFdEdx->Fill(betagamma_tof,tofparms->dEdx2); } vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,PiMinus,p); chisq_tof=dt*dt/vart_tof; PiMinusTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case KPlus: KPlusTOFdt->Fill(p,dt); vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,KPlus,p); chisq_tof=dt*dt/vart_tof; KPlusTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case KMinus: KMinusTOFdt->Fill(p,dt); vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,KMinus,p); chisq_tof=dt*dt/vart_tof; KMinusTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; case Proton: vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,Proton,p); chisq_tof=dt*dt/vart_tof; ProtonTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); { double dedx_tof=mypid_algorithm->GetProtondEdxMean_TOF(beta); double dedx_tof_diff1=tofparms->dEdx1-dedx_tof; double dedx_tof_diff2=tofparms->dEdx2-dedx_tof; if (tofparms->dEdx1>0.){ ProtonTOFdEdx1->Fill(betagamma_tof,tofparms->dEdx1); DProtonTOFdEdx1->Fill(betagamma_tof,dedx_tof_diff1); } if (tofparms->dEdx2>0.){ ProtonTOFdEdx2->Fill(betagamma_tof,tofparms->dEdx2); DProtonTOFdEdx2->Fill(betagamma_tof,dedx_tof_diff2); } } ProtonTOFdt->Fill(p,dt); if (dEdx_fdc>0.&&TMath::Prob(chisq_tof,1)>0.01){ ProtondEdxFDCwithTOFcut->Fill(betagamma_tof,dEdx_fdc); } break; case AntiProton: AntiProtonTOFdt->Fill(p,dt); vart_tof=mypid_algorithm->GetTimeVariance(SYS_TOF,AntiProton,p); chisq_tof=dt*dt/vart_tof; AntiProtonTOFdtCL->Fill(TMath::Prob(chisq_tof,1.)); break; default: break; } } if (fcalparms!=NULL){ //printf("FCAL\n"); double E_over_p=fcalparms->dFCALShower->getEnergy()/p;//fcalparms->dMomentum; double E_over_p_mean=mypid_algorithm->GetEOverPMean(SYS_FCAL,p);//fcalparms->dMomentum); double E_over_p_diff=E_over_p-E_over_p_mean; double E_over_p_sigma=mypid_algorithm->GetEOverPSigma(SYS_BCAL,p);//fcalparms->dMomentum); double chisq_E_over_p=E_over_p_diff*E_over_p_diff/(E_over_p_sigma*E_over_p_sigma); // double chisq_E_over_p=hyp->Get_ChiSq_EoverP(); //cout << chisq_E_over_p << endl; if (myPID==Electron){ EOverPMinus->Fill(p/*fcalparms->dMomentum*/,E_over_p); DEOverPMinus->Fill(p/*fcalparms->dMomentum*/,E_over_p_diff); ElectronEOverPCL->Fill(TMath::Prob(chisq_E_over_p,1.)); } else if (myPID==Positron){ EOverPPlus->Fill(p/*fcalparms->dMomentum*/,E_over_p); DEOverPPlus->Fill(p/*fcalparms->dMomentum*/,E_over_p_diff); PositronEOverPCL->Fill(TMath::Prob(chisq_E_over_p,1.)); } double dt_fcal=fcalparms->dFCALShower->getTime()-fcalparms->dFlightTime-t0_rf;//-pid_algorithm->GetTimeMean(SYS_FCAL,Positron,p); //cout << dt_fcal << " " << fcalparms->dFlightTime << endl; double vart_fcal=0.; double chisq_fcal=0.; switch(myPID){ case Positron: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,Positron,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; PositronFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); PositronFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; case Electron: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,Electron,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; ElectronFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); ElectronFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; case PiPlus: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,PiPlus,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; PiPlusFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); PiPlusFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; case PiMinus: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,PiMinus,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; PiMinusFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); PiMinusFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); if (tofparms==NULL){ PiMinusFCALdtNoTOF->Fill(p/*fcalparms->dMomentum*/,dt_fcal); } break; case KPlus: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,KPlus,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; KPlusFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); KPlusFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; case KMinus: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,KMinus,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; KMinusFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); KMinusFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; case Proton: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,Proton,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; ProtonFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); ProtonFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); break; case AntiProton: vart_fcal=mypid_algorithm->GetTimeVariance(SYS_FCAL,AntiProton,p); chisq_fcal=dt_fcal*dt_fcal/vart_fcal; AntiProtonFCALdt->Fill(p/*fcalparms->dMomentum*/,dt_fcal); AntiProtonFCALdtCL->Fill(TMath::Prob(chisq_fcal,1.)); break; default: break; } } } // Run the kinematic fitter requiring energy and momentum conservation bool JEventProcessor_PIDCheck::DoKinematicFit(double t0_rf, const DBeamPhoton *beamphoton, map >&chargedParticles, DKinFitUtils_GlueX *dKinFitUtils, DKinFitter *dKinFitter, Particle_t RecoilPID) const { set> InitialParticles, FinalParticles; dKinFitter->Reset_NewFit(); shared_ptrmyBeam=dKinFitUtils->Make_BeamParticle(beamphoton); shared_ptrmyTarget=dKinFitUtils->Make_TargetParticle(Deuteron); InitialParticles.insert(myBeam); InitialParticles.insert(myTarget); // Vertex info //DLorentzVector vert4; // vert4.SetVect(vertex); //vert4.SetT(t0_rf); // Missing particle shared_ptrmyRecoil=dKinFitUtils->Make_MissingParticle(RecoilPID); FinalParticles.insert(myRecoil); Particle_t particle_list[6]={PiPlus,PiMinus,KPlus,KMinus,Proton,AntiProton}; for (unsigned int k=0;k<6;k++){ vectormyParticles=chargedParticles[particle_list[k]]; for (unsigned int m=0;mmyParticle=dKinFitUtils->Make_DetectedParticle(myParticles[m]); FinalParticles.insert(myParticle); } } // make constraint shared_ptrlocP4Constraint = dKinFitUtils->Make_P4Constraint(InitialParticles, FinalParticles); // set constraint dKinFitter->Add_Constraint(locP4Constraint); // PERFORM THE KINEMATIC FIT return dKinFitter->Fit_Reaction(); }