// $Id$ // // File: JEventProcessor_physics_detcom.cc // Created: Fri Dec 5 13:37:34 EST 2014 // Creator: jrsteven (on Linux ifarm1102 2.6.32-220.7.1.el6.x86_64 x86_64) // #include "JEventProcessor_physics_detcom.h" using namespace jana; #include // DAQ header files #include "DAQ/Df125PulseIntegral.h" #include "DAQ/Df125WindowRawData.h" #include "DAQ/Df250PulseIntegral.h" #include "DAQ/Df250WindowRawData.h" // Tracking header files #include "TRACKING/DTrackCandidate.h" #include "TRACKING/DTrackWireBased.h" #include "CDC/DCDCHit.h" #include "CDC/DCDCDigiHit.h" // ST header files #include "START_COUNTER/DSCHit.h" #include "START_COUNTER/DSCDigiHit.h" #include "START_COUNTER/DSCTDCDigiHit.h" // TOF header files #include "TOF/DTOFPaddleHit.h" #include "TOF/DTOFHit.h" #include "TOF/DTOFDigiHit.h" // FCAL/BCAL header files #include "FCAL/DFCALHit.h" #include "FCAL/DFCALCluster.h" #include "BCAL/DBCALShower.h" // PID #include "PID/DDetectorMatches.h" // TAGGER header files #include "TAGGER/DTAGHHit.h" #include "TAGGER/DTAGHDigiHit.h" #include "TAGGER/DTAGMHit.h" // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_physics_detcom()); } } // "C" //------------------ // JEventProcessor_physics_detcom (Constructor) //------------------ JEventProcessor_physics_detcom::JEventProcessor_physics_detcom() { } //------------------ // ~JEventProcessor_physics_detcom (Destructor) //------------------ JEventProcessor_physics_detcom::~JEventProcessor_physics_detcom() { } //------------------ // init //------------------ jerror_t JEventProcessor_physics_detcom::init(void) { // This is called once at program startup. If you are creating // and filling historgrams in this plugin, you should lock the // ROOT mutex like this: // // japp->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); // japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { Energy_BCAL_FCAL=new TH2F("Energy_BCAL_FCAL","Summed BCAL energy vs FCAL energy; E_{FCAL}; E_{BCAL}",100,0.,10.,100,0.,5.); Event_Stat=new TH1F("Event_Status","Event Status",10,0,10); // track quality histograms Track_Stat=new TH1F("Track_Status","Track Status",10,0,10); Track_Ndof_Theta=new TH2F("Track_Ndof_Theta","Track Ndof vs #theta; #theta; Ndof",180,0,180,50,0,50); Track_WireNdof_Theta=new TH2F("Track_WireNdof_Theta","Wire-based Track Ndof vs #theta; #theta; Ndof",180,0,180,50,0,50); Track_Ndof_Theta_NoMatch=new TH2F("Track_Ndof_Theta_NoMatch","Track Ndof vs #theta; #theta; Ndof",180,0,180,50,0,50); Track_WireNdof_Theta_NoMatch=new TH2F("Track_WireNdof_Theta_NoMatch","Wire-based Track Ndof vs #theta; #theta; Ndof",180,0,180,50,0,50); Track_FOM=new TH1F("Track_FOM","Track FOM; FOM",100,0.,1.); NTrack_WireBased=new TH1F("NTrack_WireBased","# of wire based tracks",40,0,40); NTrack_TimeBased=new TH1F("NTrack_TimeBased","# of time based tracks",40,0,40); // dE/dx histograms Track_dEdxVsP_Pos_CDC=new TH2F("Track_dEdxVsP_Pos_CDC","CDC Positive: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Pos_CDC_70=new TH2F("Track_dEdxVsP_Pos_CDC_70","CDC Positive 70%: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Pos_CDC_50=new TH2F("Track_dEdxVsP_Pos_CDC_50","CDC Positive 50%: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Neg_CDC=new TH2F("Track_dEdxVsP_Neg_CDC","CDC Negative: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Neg_CDC_70=new TH2F("Track_dEdxVsP_Neg_CDC_70","CDC Negative 70%: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Neg_CDC_50=new TH2F("Track_dEdxVsP_Neg_CDC_50","CDC Negative 50%: dE/dx vs p",200,0,2,500,0,5); Track_dEdxVsP_Pos_FDC=new TH2F("Track_dEdxVsP_Pos_FDC","FDC Positive: dE/dx vs p",200,0,2,500,0,50); Track_dEdxVsP_Pos_FDC_70=new TH2F("Track_dEdxVsP_Pos_FDC_70","FDC Positive 70%: dE/dx vs p",200,0,2,500,0,50); Track_dEdxVsP_Pos_FDC_50=new TH2F("Track_dEdxVsP_Pos_FDC_50","FDC Positive 50%: dE/dx vs p",200,0,2,500,0,50); Track_dEdxVsP_Neg_FDC=new TH2F("Track_dEdxVsP_Neg_FDC","FDC Negative: dE/dx vs p",200,0,2,500,0,50); Track_dEdxVsP_Neg_FDC_70=new TH2F("Track_dEdxVsP_Neg_FDC_70","FDC Negative 70%: dE/dx vs p",200,0,2,500,0,50); Track_dEdxVsP_Neg_FDC_50=new TH2F("Track_dEdxVsP_Neg_FDC_50","FDC Negative 50%: dE/dx vs p",200,0,2,500,0,50); // SC-Tagger histograms Track_DeltaTProjected_SCtoTAGM=new TH2F("Track_DeltaTProjected_SCtoTAGM","Corrected SC-TAGM Energy vs #Delta t; #Delta t; Energy", 100,-50,50, 160, 2., 10.); Track_DeltaTProjected_SCtoTAGM_Nominal=new TH2F("Track_DeltaTProjected_SCtoTAGM_Nominal","Nominal SC-TAGM Energy vs #Delta t; #Delta t; Energy", 250,-400,100, 160, 2., 10.); Track_DeltaTProjected_SCtoTAGH=new TH2F("Track_DeltaTProjected_SCtoTAGH","Corrected SC-TAGH Energy vs #Delta t; #Delta t; Energy", 100,-50,50, 160, 2., 10.); Track_DeltaTProjected_SCtoTAGH_Nominal=new TH2F("Track_DeltaTProjected_SCtoTAGH_Nominal","Nominal SC-TAGH Energy vs #Delta t; #Delta t; Energy", 250,-400,100, 160, 2., 10.); // SC-TOF histograms Track_BetaVsP_Pos_SCtoTOF=new TH2F("Track_BetaVsP_Pos_SCtoTOF","Positive SC to TOF: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Neg_SCtoTOF=new TH2F("Track_BetaVsP_Neg_SCtoTOF","Negative SC to TOF: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Pos_SCtoTOF_ProtonFDC=new TH2F("Track_BetaVsP_Pos_SCtoTOF_ProtonFDC","Positive SC to TOF: #beta vs. p",200,0,4,200,-1.0,2.0); Track_fADC_TimeCorrelation_SCtoTOF=new TH2F("Track_fADC_TimeCorrelation_SCtoTOF","fADC hit time SC vs TOF; TOF fADC time; SC fADC time",4000,0,4000,4000,0,4000); Track_DeltaTProjected_SCtoTOF=new TH1F("Track_DeltaTProjected_SCtoTOF","#Delta t", 400, -200, 200); Track_DeltaT_SCtoTOF=new TH1F("Track_DeltaT_SCtoTOF","#Delta t", 400, -200, 200); // SC-FCAL histograms Track_BetaVsP_Pos_SCtoFCAL=new TH2F("Track_BetaVsP_Pos_SCtoFCAL","Positive SC to FCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Neg_SCtoFCAL=new TH2F("Track_BetaVsP_Neg_SCtoFCAL","Negative SC to FCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Pos_SCtoFCAL_ProtonFDC=new TH2F("Track_BetaVsP_Pos_SCtoFCAL_ProtonFDC","Positive SC to FCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_DeltaTProjected_SCtoFCAL=new TH1F("Track_DeltaTProjected_SCtoFCAL","#Delta t", 400, -200, 200); Track_DeltaT_SCtoFCAL=new TH1F("Track_DeltaT_SCtoFCAL","#Delta t", 400, -200, 200); Track_EvsP_Pos_FCAL=new TH2F("Track_EvsP_Pos_FCAL","Positive FCAL E vs track p; track p (GeV); FCAL Shower Energy (GeV)", 100, 0., 10., 100, 0., 10.); Track_E2PvsP_Pos_FCAL=new TH2F("Track_E2PvsP_Pos_FCAL","Positive FCAL track: E/p vs p; track p (GeV); FCAL E/p", 100, 0., 10., 100, 0., 4.); Track_E2PvsTheta_Pos_FCAL=new TH2F("Track_E2PvsTheta_Pos_FCAL","Positive FCAL track: E/p vs #theta; track #theta; FCAL E/p", 150, 0., 15., 100, 0., 4.); Track_EvsP_Neg_FCAL=new TH2F("Track_EvsP_Neg_FCAL","Negative FCAL E vs track p; track p (GeV); FCAL Shower Energy (GeV)", 100, 0., 10., 100, 0., 10.); Track_E2PvsP_Neg_FCAL=new TH2F("Track_E2PvsP_Neg_FCAL","Negative FCAL track: E/p vs p; track p (GeV); FCAL E/p", 100, 0., 10., 100, 0., 4.); Track_E2PvsTheta_Neg_FCAL=new TH2F("Track_E2PvsTheta_Neg_FCAL","Negative FCAL track: E/p vs #theta; track #theta; FCAL E/p", 150, 0., 15., 100, 0., 4.); // SC-BCAL histograms Track_BetaVsP_Pos_SCtoBCAL=new TH2F("Track_BetaVsP_Pos_SCtoBCAL","Positive SC to BCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Neg_SCtoBCAL=new TH2F("Track_BetaVsP_Neg_SCtoBCAL","Negative SC to BCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Pos_SCtoBCAL_ProtonCDC=new TH2F("Track_BetaVsP_Pos_SCtoBCAL_ProtonCDC","Positive SC to BCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_BetaVsP_Pos_SCtoBCAL_ProtonFDC=new TH2F("Track_BetaVsP_Pos_SCtoBCAL_ProtonFDC","Positive SC to BCAL: #beta vs. p",200,0,4,200,-1.0,2.0); Track_DeltaTProjected_SCtoBCAL=new TH1F("Track_DeltaTProjected_SCtoBCAL","#Delta t", 400, -200, 200); Track_DeltaT_SCtoBCAL=new TH1F("Track_DeltaT_SCtoBCAL","#Delta t", 400, -200, 200); // track multiplicity Track_Nproton=new TH1F("Track_Nproton","Track N_{proton}; N_{proton}", 5, 0, 5); Track_Npiplus=new TH1F("Track_Npiplus","Track N_{#pi^{+}}; N_{#pi^{+}}", 5, 0, 5); Track_Npiminus=new TH1F("Track_Npiminus","Track N_{#pi^{-}}; N_{#pi^{-}}", 5, 0, 5); Track_Npi0_FCAL=new TH1F("Track_Npi0_FCAL","FCAL N_{#pi^{0}}; N_{#pi^{0}}", 5, 0, 5); Track_Npi0_BCAL=new TH1F("Track_Npi0_BCAL","BCAL N_{#pi^{0}}; N_{#pi^{0}}", 5, 0, 5); // FCAL di-photon histograms DiPhoton_Mass_FCAL=new TH1F("DiPhoton_Mass_FCAL","FCAL Di-photon mass; M(2#gamma) (GeV/c^{2})",100,0.,1.0); Pi0_FCAL_Energy=new TH1F("Pi0_FCAL_Energy","FCAL #pi^{0} Energy; Energy (GeV)",100, 0., 10.0); // BCAL di-photon histograms DiPhoton_Mass_BCAL=new TH1F("DiPhoton_Mass_BCAL","BCAL Di-photon mass; M(2#gamma) (GeV/c^{2})",100,0.,1.0); DiPhoton_Ncell_Mass_BCAL=new TH2F("DiPhoton_Ncell_Mass_BCAL","BCAL Di-photon mass; M(2#gamma) (GeV/c^{2})l Ncell",100,0.,1.0,100,0,100); Pi0_BCAL_Energy=new TH1F("Pi0_BCAL_Energy","BCAL #pi^{0} Energy; Energy (GeV)",100, 0., 10.0); // missing mass off proton MMproton_DeltaT_TAGM=new TH2F("MMproton_DeltaT_TAGM","MM off proton TAGM vs #Delta t; #Delta t (SC-TAGM); MM", 200,-20,20, 200,0.,4.); MMproton_DeltaT_TAGH=new TH2F("MMproton_DeltaT_TAGH","MM off proton TAGH vs #Delta t; #Delta t (SC-TAGH); MM", 200,-20,20, 200,0.,4.); // pi+ pi- final state Mass_2pi=new TH1F("Mass_2pi","#pi^{+} #pi^{-} Mass; Mass (GeV/c^{2})",200,0.,2.); MM_M2pi_TAGM=new TH2F("MM_M2pi_TAGM","MM off #pi^{+}#pi^{-} TAGM vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200,0.,2., 200,0.,4.); MM_M2k_TAGM=new TH2F("MM_M2k_TAGM","MM off K^{+}K^{-} TAGM vs M_{K^{+}K^{-}}; M_{K^{+}K^{-}}; MM", 200,0.,2., 200,0.,4.); DeltaPz_MM_TAGM=new TH2F("DeltaPz_MM_TAGM","TAGM #Delta P_{z} vs MM; MM; p^{#pi^{+}#pi^{-}}_{z} - E_{#gamma} (TAGM)", 200,0.,4., 200,-4.,4.); MM2pi_DeltaT_TAGM=new TH2F("MM2pi_DeltaT_TAGM","MM off #pi^{+}#pi^{-} TAGM vs #Delta t; #Delta t (SC-TAGM); MM", 200,-20,20, 200,0.,4.); M2pi_DeltaT_TAGM=new TH2F("M2pi_DeltaT_TAGM","M_{#pi^{+}#pi^{-}} vs #Delta t TAGM; #Delta t (SC-TAGM); M_{#pi^{+}#pi^{-}}", 200,-20,20, 200,0.,2.); MM_M2pi_TAGH=new TH2F("MM_M2pi_TAGH","MM off #pi^{+}#pi^{-} TAGH vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200,0.,2., 200,0.,4.); MM_M2k_TAGH=new TH2F("MM_M2k_TAGH","MM off K^{+}K^{-} TAGH vs M_{K^{+}K^{-}}; M_{K^{+}K^{-}}; MM", 200,0.,2., 200,0.,4.); DeltaPz_MM_TAGH=new TH2F("DeltaPz_MM_TAGH","TAGH #Delta P_{z} vs MM; MM; p^{#pi^{+}#pi^{-}}_{z} - E_{#gamma} (TAGH)", 200,0.,4., 200,-4.,4.); MM2pi_DeltaT_TAGH=new TH2F("MM2pi_DeltaT_TAGH","MM off #pi^{+}#pi^{-} TAGH vs #Delta t; #Delta t (SC-TAGH); MM", 200,-20,20, 200,0.,4.); M2pi_DeltaT_TAGH=new TH2F("M2pi_DeltaT_TAGH","M_{#pi^{+}#pi^{-}} vs #Delta t TAGH; #Delta t (SC-TAGH); M_{#pi^{+}#pi^{-}}", 200,-20,20, 200,0.,2.); DalitzPlot_p2pi=new TH2F("DalitzPlot_p2pi","Dalitz plot p #pi^{+} #pi^{-}; M^{2}_{p#pi^{+}}; M^{2}_{p#pi^{-}}",200,0.,5.,200,0.,5.); MM2_p2pi_DeltaT_TAGM=new TH2F("MM2_p2pi_DeltaT_TAGM","MM^{2} off p#pi^{+}#pi^{-} TAGM vs #Delta t (SC-TAGM); #Delta t; MM^{2}", 200,-20,20, 200,-4.,4.); Exclusive_M2pi_DeltaT_TAGM=new TH2F("Exclusive_M2pi_DeltaT_TAGM","M_{#pi^{+}#pi^{-}} vs #Delta t TAGM; #Delta t (SC-TAGM); M_{#pi^{+}#pi^{-}}", 200,-20,20, 200,0.,2.); Exclusive_MM2_M2pi_TAGM=new TH2F("Exclusive_MM2_M2pi_TAGM","TAGM MM^{2} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM^{2}", 200,0.,2., 200,-4.,4.); Exclusive_DalitzPlot_p2pi_TAGM=new TH2F("Exclusive_DalitzPlot_p2pi_TAGM","TAGM: Exclusive Dalitz plot p #pi^{+} #pi^{-}; M^{2}_{#pi^{+}#pi^{-}}; M^{2}_{p#pi^{+}}",200,0.,5.,200,0.,5.); MM2_p2pi_DeltaT_TAGH=new TH2F("MM2_p2pi_DeltaT_TAGH","MM^{2} off p#pi^{+}#pi^{-} TAGH vs #Delta t; #Delta t (SC-TAGH); MM^{2}", 200,-20,20, 200,-4.,4.); Exclusive_M2pi_DeltaT_TAGH=new TH2F("Exclusive_M2pi_DeltaT_TAGH","M_{#pi^{+}#pi^{-}} vs #Delta t TAGH; #Delta t (SC-TAGH); M_{#pi^{+}#pi^{-}}", 200,-20,20, 200,0.,2.); Exclusive_MM2_M2pi_TAGH=new TH2F("Exclusive_MM2_M2pi_TAGH","TAGH MM^{2} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM^{2}", 200,0.,2., 200,-4.,4.); Exclusive_DalitzPlot_p2pi_TAGH=new TH2F("Exclusive_DalitzPlot_p2pi_TAGH","TAGH: Exclusive Dalitz plot p #pi^{+} #pi^{-}; M^{2}_{#pi^{+}#pi^{-}}; M^{2}_{p#pi^{+}}",200,0.,5.,200,0.,5.); // recoil proton from pi+ pi- p final state string nameTag[2] = {"TAGM","TAGH"}; for(int itag = 0; itag < 2; itag++){ RecoilP_PhiTrack_PhiMissing[itag] = new TH2F(Form("RecoilP_PhiTrack_PhiMissing_%s",nameTag[itag].data()),Form("Recoil Proton %s: #phi_{track} vs #phi_{missing}; #phi_{missing}; #phi_{track}",nameTag[itag].data()),360,-180,180,360,-180,180); RecoilP_DeltaPhi_DeltaTheta[itag] = new TH2F(Form("RecoilP_DeltaPhi_DeltaTheta_%s",nameTag[itag].data()),Form("Recoil Proton %s: #Delta#phi vs #Delta#theta; #Delta#theta; #Delta#phi",nameTag[itag].data()),360,-180,180,360,-180,180); RecoilP_dEdxVsP_CDC[itag] = new TH2F(Form("RecoilP_dEdxVsP_CDC_%s",nameTag[itag].data()),Form("Recoil Proton %s: CDC dE/dx vs p",nameTag[itag].data()),500,0,5,500,0,5); RecoilP_dEdxVsP_FDC[itag] = new TH2F(Form("RecoilP_dEdxVsP_FDC_%s",nameTag[itag].data()),Form("Recoil Proton %s: FDC dE/dx vs p",nameTag[itag].data()),500,0,5,500,0,50); RecoilP_BetaVsP_SCtoBCAL[itag] = new TH2F(Form("RecoilP_BetaVsP_SCtoBCAL_%s",nameTag[itag].data()),Form("Recoil Proton %s: SC to BCAL #beta vs. p",nameTag[itag].data()),500,0,5,200,-1.0,2.0); } // pi+ pi- pi0 final state DiPhoton_Mass_FCAL_3pi=new TH1F("DiPhoton_Mass_FCAL_3pi","FCAL Di-photon mass; M(2#gamma) (GeV/c^{2})",100,0.,1.0); DiPhoton_Mass_BCAL_3pi=new TH1F("DiPhoton_Mass_BCAL_3pi","BCAL Di-photon mass; M(2#gamma) (GeV/c^{2})",100,0.,1.0); string name[2] = {"FCAL","BCAL"}; for(int iCal = 0; iCal < 2; iCal++){ Mass_3pi[iCal]=new TH1F(Form("Mass_3pi_%s",name[iCal].data()),Form("%s: #pi^{+} #pi^{-} #pi^{0} Mass; Mass (GeV/c^{2})",name[iCal].data()),200,0.,2.); MM_M3pi_TAGM[iCal]=new TH2F(Form("MM_M3pi_TAGM_%s",name[iCal].data()),Form("%s: MM off #pi^{0}#pi^{+}#pi^{-} TAGM vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; MM",name[iCal].data()), 200,0.,2., 200,0.,4.); MM_M3pi_TAGH[iCal]=new TH2F(Form("MM_M3pi_TAGH_%s",name[iCal].data()),Form("%s: MM off #pi^{0}#pi^{+}#pi^{-} TAGH vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; MM",name[iCal].data()), 200,0.,2., 200,0.,4.); DeltaT_M3pi_TAGM[iCal]=new TH2F(Form("DeltaT_M3pi_TAGM_%s",name[iCal].data()),Form("%s: (#pi^{0} - TAGM) #Delta t vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; #Delta t (#pi^{0} - TAGM)",name[iCal].data()), 200,0.,2., 200,-20,20); DeltaTtrack_M3pi_TAGM[iCal]=new TH2F(Form("DeltaTtrack_M3pi_TAGM_%s",name[iCal].data()),Form("%s: TAGM (#pi^{0} - #pi^{-}) #Delta t vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; #Delta t (#pi^{0} - #pi^{-})",name[iCal].data()), 200,0.,2., 200,-20,20); DeltaT_M3pi_TAGH[iCal]=new TH2F(Form("DeltaT_M3pi_TAGH_%s",name[iCal].data()),Form("%s: (#pi^{0} - TAGH) #Delta t vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; #Delta t (#pi^{0} - TAGH)",name[iCal].data()), 200,0.,2., 200,-20,20); DeltaTtrack_M3pi_TAGH[iCal]=new TH2F(Form("DeltaTtrack_M3pi_TAGH_%s",name[iCal].data()),Form("%s: TAGH (#pi^{0} - #pi^{-}) #Delta t vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; #Delta t (#pi^{0} - #pi^{-})",name[iCal].data()), 200,0.,2., 200,-20,20); MM2_p3pi_DeltaT_TAGM[iCal]=new TH2F(Form("MM2_p3pi_DeltaT_TAGM_%s",name[iCal].data()),Form("%s: MM^{2} off p#pi^{0}#pi^{+}#pi^{-} TAGM vs #Delta t; #Delta t (SC-TAGM); MM^{2}",name[iCal].data()), 200,-20,20, 200,-4.,4.); M3pi_DeltaT_TAGM[iCal]=new TH2F(Form("M3pi_DeltaT_TAGM_%s",name[iCal].data()),Form("%s: M_{#pi^{0}#pi^{+}#pi^{-}} vs #Delta t TAGM; #Delta t (SC-TAGM); M_{#pi^{0}#pi^{+}#pi^{-}}",name[iCal].data()), 200,-20,20, 200,0.,2.); MM2_M3pi_TAGM[iCal]=new TH2F(Form("MM2_M3pi_TAGM_%s",name[iCal].data()),Form("%s: MM^{2} off p#pi^{0}#pi^{+}#pi^{-} TAGM vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; MM^{2}",name[iCal].data()), 200,0.,2., 200,-4.,4.); MM2_p3pi_DeltaT_TAGH[iCal]=new TH2F(Form("MM2_p3pi_DeltaT_TAGH_%s",name[iCal].data()),Form("%s: MM^{2} off p#pi^{0}#pi^{+}#pi^{-} TAGH vs #Delta t; #Delta t (SC-TAGH); MM^{2}",name[iCal].data()), 200,-20,20, 200,-4.,4.); M3pi_DeltaT_TAGH[iCal]=new TH2F(Form("M3pi_DeltaT_TAGH_%s",name[iCal].data()),Form("%s: M_{#pi^{0}#pi^{+}#pi^{-}} vs #Delta t TAGH; #Delta t (SC-TAGH); M_{#pi^{0}#pi^{+}#pi^{-}}",name[iCal].data()), 200,-20,20, 200,0.,2.); MM2_M3pi_TAGH[iCal]=new TH2F(Form("MM2_M3pi_TAGH_%s",name[iCal].data()),Form("%s: MM^{2} off p#pi^{0}#pi^{+}#pi^{-} TAGH vs M_{#pi^{0}#pi^{+}#pi^{-}}; M_{#pi^{0}#pi^{+}#pi^{-}}; MM^{2}",name[iCal].data()), 200,0.,2., 200,-4.,4.); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_physics_detcom::brun(JEventLoop *eventLoop, int runnumber) { // This is called whenever the run number changes vector dParticleID_algos; eventLoop->Get(dParticleID_algos); if(dParticleID_algos.size()<1){ _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"< 1500 && runnumber < 1810){ // early running uncalibrated offset_tagm = -159.; offset_tagh = -102.; if(runnumber > 1500 && runnumber < 1530){ // FCAL trigger runs offset_runid = 15.5; } else if(runnumber > 1540 && runnumber < 1610){ // early BCAL trigger runs offset_runid = -44.5; } else if(runnumber > 1779 && runnumber < 1790){ // later BCAL trigger runs offset_runid = -39.5; } else if(runnumber > 1801 && runnumber < 1810){ // later BCAL trigger runs (18NN) offset_runid = 1.0; } } #if 0 else if(runnumber > 2139 && runnumber < 2181){ // calibrated by Simon (2140 to 2180) offset_tagm = -7.5; offset_tagh = -31.; } else if(runnumber > 2180 && runnumber < 2500){ // later running (not sure what's in ccdb?) offset_tagm = -28.; offset_tagh = +28.5; } #endif // scale factor for MC cdcMCscale = 1.; scMCscale = 1.; if(runnumber > 9000) { cdcMCscale = 1e3; scMCscale = 1e3; } return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_physics_detcom::evnt(JEventLoop *loop, int eventnumber) { // This is called for every event. Use of common resources like writing // to a file or filling a histogram should be mutex protected. Using // loop->Get(...) to get reconstructed objects (and thereby activating the // reconstruction algorithm) should be done outside of any mutex lock // since multiple threads may call this method at the same time. // Here's an example: // // vector mydataclasses; // loop->Get(mydataclasses); // // japp->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); vector locFCALShowers; loop->Get(locFCALShowers); double locFCALEnergy = 0; for(uint ifcal=0; ifcalgetEnergy(); vector locBCALShowers; loop->Get(locBCALShowers); double locBCALEnergy = 0; for(uint ibcal=0; ibcalE; japp->RootWriteLock(); if(locFCALEnergy!=0. && locBCALEnergy!=0) Energy_BCAL_FCAL->Fill(locFCALEnergy, locBCALEnergy); japp->RootUnLock(); Event_Stat->Fill(0); // get tagger hits to correlate with downstream events vector< pair > locGoodTAGMHits_Nominal; vector< pair > locGoodTAGHHits_Nominal; get_tagger_hits(loop, locGoodTAGMHits_Nominal, locGoodTAGHHits_Nominal); // correct timing for tagger hits (ad-hoc correction by eye) goodTAGMHits.clear(); goodTAGHHits.clear(); for(uint itagm=0; itagm > locFCALpi0_p4; get_fcal_pi0(loop, locFCALpi0_p4); vector< pair > locBCALpi0_p4; get_bcal_pi0(loop, locBCALpi0_p4); // get tracks and detector matches vector locTrackWireBasedVector; loop->Get(locTrackWireBasedVector); vector locTrackTimeBasedVector; loop->Get(locTrackTimeBasedVector); const DDetectorMatches* locDetectorMatches = NULL; loop->GetSingle(locDetectorMatches); NTrack_WireBased->Fill(locTrackWireBasedVector.size()); NTrack_TimeBased->Fill(locTrackTimeBasedVector.size()); // tracking cuts uint minHitsCDC = 5; uint minHitsFDC = 5; double protonCutCDC = 1.; double protonCutFDC = 15.; double locE2Pcut = 1.0; double locMaxBetaProton_TOF = 0.7; protonP4.clear(); piplusP4.clear(); piplus_chisq.clear(); piminusP4.clear(); piminus_chisq.clear(); kplusP4.clear(); kminusP4.clear(); electronP4.clear(); positronP4.clear(); if(locTrackTimeBasedVector.size() > 0) Event_Stat->Fill(1); if(locTrackTimeBasedVector.size() > 1) Event_Stat->Fill(2); //////////////////////////// // loop time-based tracks // //////////////////////////// for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i){ const DTrackWireBased* locTrackWireBased = NULL; locTrackTimeBasedVector[loc_i]->GetSingle(locTrackWireBased); //if(!locTrackWireBased) continue; bool protonTagCDC = false; bool protonTagFDC = false; bool protonTagTOF = false; bool electronTag = false; bool positronTag = false; double locTimeProjected = -999.; Track_Stat->Fill(0); // cut to get tracks from target double locPositionZ = locTrackTimeBasedVector[loc_i]->position().Z(); if(locPositionZ < 50 || locPositionZ > 80) continue; Track_Stat->Fill(1); int Ndof = locTrackTimeBasedVector[loc_i]->Ndof; int wireNdof = 0; //locTrackWireBased->Ndof; Track_Ndof_Theta_NoMatch->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),Ndof); Track_WireNdof_Theta_NoMatch->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),wireNdof); // cut on match to some detector if(!locDetectorMatches->Get_IsMatchedToHit(locTrackTimeBasedVector[loc_i])) continue; Track_Stat->Fill(2); // cut on track nHits Track_Ndof_Theta->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),Ndof); Track_WireNdof_Theta->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),wireNdof); if(Ndof < 3) continue; Track_Stat->Fill(3); // cut on track quality double FOM = TMath::Prob(locTrackTimeBasedVector[loc_i]->chisq,locTrackTimeBasedVector[loc_i]->Ndof); Track_FOM->Fill(FOM); if(FOM < 1e-3) continue; Track_Stat->Fill(4); int locCharge = locTrackTimeBasedVector[loc_i]->charge(); double locP = locTrackTimeBasedVector[loc_i]->momentum().Mag(); double locTheta = locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(); ///////////////////////////////////////////////////////// // calculate dE/dx different truncated mean thresholds // ///////////////////////////////////////////////////////// double locdEdx_FDC_70, locdEdx_CDC_70; double locdEdx_FDC_50, locdEdx_CDC_50; vector locdEdxHits_CDC, locdEdxHits_FDC; dParticleID->GetDCdEdxHits(locTrackTimeBasedVector[loc_i], locdEdxHits_CDC, locdEdxHits_FDC); // remove 30% of hits calcdEdx(locTrackTimeBasedVector[loc_i], locdEdxHits_CDC, locdEdxHits_FDC, locdEdx_FDC_70, locdEdx_CDC_70, 0.7); // remove 50% of hits (default algo) calcdEdx(locTrackTimeBasedVector[loc_i], locdEdxHits_CDC, locdEdxHits_FDC, locdEdx_FDC_50, locdEdx_CDC_50, 0.5); /////////////////////////// // fill dE/dx in the DCs // /////////////////////////// japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { if(locCharge > 0) { if(locTrackTimeBasedVector[loc_i]->dNumHitsUsedFordEdx_CDC > minHitsCDC) { Track_Stat->Fill(5); Track_dEdxVsP_Pos_CDC->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_CDC*1e3*cdcMCscale); Track_dEdxVsP_Pos_CDC_70->Fill(locP,locdEdx_CDC_70*1e3*cdcMCscale); Track_dEdxVsP_Pos_CDC_50->Fill(locP,locdEdx_CDC_50*1e3*cdcMCscale); if(locTrackTimeBasedVector[loc_i]->ddEdx_CDC*1e3*cdcMCscale > protonCutCDC){ // tag protons protonTagCDC = true; // dump some waveforms for Dave M. //dumpWaveform(locTrackTimeBasedVector[loc_i], loc_i, eventnumber); } } else if(locTrackTimeBasedVector[loc_i]->dNumHitsUsedFordEdx_FDC > minHitsFDC) { Track_Stat->Fill(6); Track_dEdxVsP_Pos_FDC->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_FDC*1e6); Track_dEdxVsP_Pos_FDC_70->Fill(locP,locdEdx_FDC_70*1e6); Track_dEdxVsP_Pos_FDC_50->Fill(locP,locdEdx_FDC_50*1e6); if(locTrackTimeBasedVector[loc_i]->ddEdx_FDC*1e6 > protonCutFDC){ // tag protons protonTagFDC = true; } } } if(locCharge < 0) { if(locTrackTimeBasedVector[loc_i]->dNumHitsUsedFordEdx_CDC > minHitsCDC) { Track_Stat->Fill(5); Track_dEdxVsP_Neg_CDC->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_CDC*1e3*cdcMCscale); Track_dEdxVsP_Neg_CDC_70->Fill(locP,locdEdx_CDC_70*1e3*cdcMCscale); Track_dEdxVsP_Neg_CDC_50->Fill(locP,locdEdx_CDC_50*1e3*cdcMCscale); } else if(locTrackTimeBasedVector[loc_i]->dNumHitsUsedFordEdx_FDC > minHitsFDC) { Track_Stat->Fill(6); Track_dEdxVsP_Neg_FDC->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_FDC*1e6); Track_dEdxVsP_Neg_FDC_70->Fill(locP,locdEdx_FDC_70*1e6); Track_dEdxVsP_Neg_FDC_50->Fill(locP,locdEdx_FDC_50*1e6); } } } japp->RootUnLock(); //RELEASE ROOT LOCK!! ////////////////////////////////////////// // get best matches to SC/TOF/FCAL/BCAL // ////////////////////////////////////////// DSCHitMatchParams locSCHitMatchParams; DTOFHitMatchParams locTOFHitMatchParams; //DShowerMatchParams locFCALShowerMatchParams; //DShowerMatchParams locBCALShowerMatchParams; DFCALShowerMatchParams locFCALShowerMatchParams; DBCALShowerMatchParams locBCALShowerMatchParams; bool foundSC = dParticleID->Get_BestSCMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locSCHitMatchParams); bool foundTOF = dParticleID->Get_BestTOFMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locTOFHitMatchParams); bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locFCALShowerMatchParams); bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locBCALShowerMatchParams); ///////////////////////////// // get E/p for FCAL tracks // ///////////////////////////// double locE2P = 999.; if(foundFCAL){ const DFCALShower* locFCALShower = dynamic_cast(locFCALShowerMatchParams.dFCALShower); //dShowerObject double locFCALEnergy = locFCALShower->getEnergy(); locE2P = locFCALEnergy/locP; japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { if(locCharge > 0) { Track_EvsP_Pos_FCAL->Fill(locP,locFCALEnergy); Track_E2PvsP_Pos_FCAL->Fill(locP,locFCALEnergy/locP); Track_E2PvsTheta_Pos_FCAL->Fill(locTheta,locFCALEnergy/locP); if(locFCALEnergy/locP > locE2Pcut) positronTag = true; } if(locCharge < 0) { Track_EvsP_Neg_FCAL->Fill(locP,locFCALEnergy); Track_E2PvsP_Neg_FCAL->Fill(locP,locFCALEnergy/locP); Track_E2PvsTheta_Neg_FCAL->Fill(locTheta,locFCALEnergy/locP); if(locFCALEnergy/locP > locE2Pcut) electronTag = true; } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } /////////////////////////////// // get SC hit if found match // /////////////////////////////// const DSCHit* locSCHit = NULL; const DSCDigiHit* locSCDigiHit = NULL; const DSCTDCDigiHit* locSCTDCDigiHit = NULL; double locSCTime = -999.; double locSCPathLength = -999.; if(foundSC){ locSCHit = locSCHitMatchParams.dSCHit; locSCHit->GetSingle(locSCDigiHit); locSCHit->GetSingle(locSCTDCDigiHit); locSCTime = locSCHitMatchParams.dHitTime - locSCHitMatchParams.dFlightTime; locSCPathLength = locSCHitMatchParams.dPathLength; locTimeProjected = locSCTime; // find tagger-SC coincidence time for(uint itag=0; itagFill(locDeltaT_TAGM_Nominal, locGoodTAGMHits_Nominal[itag].first); double locDeltaT_TAGM = locSCTime - goodTAGMHits[itag].second; Track_DeltaTProjected_SCtoTAGM->Fill(locDeltaT_TAGM, goodTAGMHits[itag].first); } for(uint itag=0; itagFill(locDeltaT_TAGH_Nominal, locGoodTAGHHits_Nominal[itag].first); double locDeltaT_TAGH = locSCTime - goodTAGHHits[itag].second; Track_DeltaTProjected_SCtoTAGH->Fill(locDeltaT_TAGH, goodTAGHHits[itag].first); } } ////////////////////////// // SC + TOF coincidence // ////////////////////////// if(foundSC && foundTOF) { const DTOFPoint* locTOFPoint = locTOFHitMatchParams.dTOFPoint; ////////////////////////////////////////////////////// // get digiHits for fADC pulse_time and f1TDC time // ////////////////////////////////////////////////////// vector locTOFPaddleHitVector; locTOFPoint->Get(locTOFPaddleHitVector); for(size_t loc_j=0; loc_j locTOFHitVector; locTOFPaddleHitVector[loc_j]->Get(locTOFHitVector); for(size_t loc_k=0; loc_kGetSingle(locTOFDigiHit); if(!locTOFDigiHit || !locSCDigiHit) continue; // shouldn't happen japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { Track_fADC_TimeCorrelation_SCtoTOF->Fill(locTOFDigiHit->pulse_time, locSCDigiHit->pulse_time); } japp->RootUnLock(); //RELEASE ROOT LOCK!! } } double locTOFTime = locTOFPoint->t; double locTOFPathLength = locTOFHitMatchParams.dPathLength; double deltaPathLength = locTOFPathLength - locSCPathLength; double deltaT = locTOFTime - locSCTime; double deltaTProjected = (locTOFPoint->t - locTOFHitMatchParams.dFlightTime) - locSCTime; double beta = deltaPathLength/(29.9792458*deltaT); //cout << "TOF pathlength = " << locTOFPathLength << " SC pathLength = " << locSCPathLength << endl; //cout << "SC time = " << locSCTime << " TOF time = " << locTOFTime << " beta = " << beta <RootWriteLock(); //ACQUIRE ROOT LOCK!! { Track_DeltaTProjected_SCtoTOF->Fill(deltaTProjected); Track_DeltaT_SCtoTOF->Fill(deltaT); // fill beta vs p histograms if(locCharge > 0) { Track_BetaVsP_Pos_SCtoTOF->Fill(locP, beta); if(protonTagFDC) Track_BetaVsP_Pos_SCtoTOF_ProtonFDC->Fill(locP, beta); //if(beta < locMaxBetaProton_TOF && beta > 0) protonTagTOF = true; } if(locCharge < 0) { Track_BetaVsP_Neg_SCtoTOF->Fill(locP, beta); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } }// end SCtoTOF ////////////////////////// // SC + FCAL coincidence // ////////////////////////// if(foundSC && foundFCAL) { const DFCALShower* locFCALShower = dynamic_cast(locFCALShowerMatchParams.dFCALShower); //dShowerObject double locFCALTime = locFCALShower->getTime(); double locFCALPathLength = locFCALShowerMatchParams.dPathLength; double deltaPathLength = locFCALPathLength - locSCPathLength; double deltaT = locFCALTime - locSCTime; double deltaTProjected = (locFCALShower->getTime() - locFCALShowerMatchParams.dFlightTime) - locSCTime; double beta = deltaPathLength/(29.9792458*deltaT); //cout << "FCAL pathlength = " << locFCALPathLength << " SC pathLength = " << locSCPathLength << endl; //cout << "SC time = " << locSCTime << " FCAL time = " << locFCALTime << " beta = " << beta <RootWriteLock(); //ACQUIRE ROOT LOCK!! { Track_DeltaTProjected_SCtoFCAL->Fill(deltaTProjected); Track_DeltaT_SCtoFCAL->Fill(deltaT); // fill beta vs p histograms if(locCharge > 0) { Track_BetaVsP_Pos_SCtoFCAL->Fill(locP, beta); if(protonTagFDC) Track_BetaVsP_Pos_SCtoFCAL_ProtonFDC->Fill(locP, beta); } if(locCharge < 0) { Track_BetaVsP_Neg_SCtoFCAL->Fill(locP, beta); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! }// end SCtoFCAL /////////////////////////// // SC + BCAL coincidence // /////////////////////////// if(foundSC && foundBCAL) { const DBCALShower* locBCALShower = dynamic_cast(locBCALShowerMatchParams.dBCALShower); //dShowerObject double locBCALTime = locBCALShower->t; double locBCALPathLength = locBCALShowerMatchParams.dPathLength; double deltaPathLength = locBCALPathLength - locSCPathLength; double deltaTProjected = (locBCALShower->t - locBCALShowerMatchParams.dFlightTime) - locSCTime; double deltaT = locBCALTime - locSCTime; double beta = deltaPathLength/(29.9792458*deltaT); //cout << "BCAL pathlength = " << locBCALPathLength << " SC pathLength = " << locSCPathLength << endl; //cout << "SC time = " << locSCTime << " BCAL time = " << locBCALTime << " beta = " << beta <RootWriteLock(); //ACQUIRE ROOT LOCK!! { Track_DeltaTProjected_SCtoBCAL->Fill(deltaTProjected); Track_DeltaT_SCtoBCAL->Fill(deltaT); // fill beta vs p histograms if(locCharge > 0) { Track_BetaVsP_Pos_SCtoBCAL->Fill(locP, beta); if(protonTagCDC) Track_BetaVsP_Pos_SCtoBCAL_ProtonCDC->Fill(locP, beta); if(protonTagFDC) Track_BetaVsP_Pos_SCtoBCAL_ProtonFDC->Fill(locP, beta); //if(beta < locMaxBetaProton_BCAL && beta > 0) protonTagBCAL = true; } if(locCharge < 0) { Track_BetaVsP_Neg_SCtoBCAL->Fill(locP, beta); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } // end SCtoBCAL // remove tracks not matched to SC (need for matching to TAGM and TAGH) if(locTimeProjected < -998.) continue; // assign tagged tracks to proper PID (assign projected time at target from SC) // if not tagged as proton/e+/e- then assume is pion (some proton/kaon leakage) if(protonTagCDC || protonTagFDC || protonTagTOF){ TLorentzVector protonTrueMassP4; protonTrueMassP4.SetVectM(locTrackTimeBasedVector[loc_i]->lorentzMomentum().Vect(), 0.938); protonP4.push_back(make_pair(protonTrueMassP4, locTimeProjected)); } else if(positronTag){ positronP4.push_back(make_pair(locTrackTimeBasedVector[loc_i]->lorentzMomentum(), locTimeProjected)); } else if(electronTag){ electronP4.push_back(make_pair(locTrackTimeBasedVector[loc_i]->lorentzMomentum(), locTimeProjected)); } else if(locCharge > 0) { piplusP4.push_back(make_pair(locTrackTimeBasedVector[loc_i]->lorentzMomentum(), locTimeProjected)); piplus_chisq.push_back(locTrackTimeBasedVector[loc_i]->chisq); TLorentzVector kaonTrueMassP4; kaonTrueMassP4.SetVectM(locTrackTimeBasedVector[loc_i]->lorentzMomentum().Vect(), 0.494); kplusP4.push_back(make_pair(kaonTrueMassP4, locTimeProjected)); } else if(locCharge < 0) { piminusP4.push_back(make_pair(locTrackTimeBasedVector[loc_i]->lorentzMomentum(), locTimeProjected)); piminus_chisq.push_back(locTrackTimeBasedVector[loc_i]->chisq); TLorentzVector kaonTrueMassP4; kaonTrueMassP4.SetVectM(locTrackTimeBasedVector[loc_i]->lorentzMomentum().Vect(), 0.494); kminusP4.push_back(make_pair(kaonTrueMassP4, locTimeProjected)); } } // end time-based track loop // fill multiplicity histograms japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { Track_Nproton->Fill(protonP4.size()); Track_Npiplus->Fill(piplusP4.size()); Track_Npiminus->Fill(piminusP4.size()); Track_Npi0_FCAL->Fill(locFCALpi0_p4.size()); Track_Npi0_BCAL->Fill(locBCALpi0_p4.size()); } japp->RootUnLock(); //RELEASE ROOT LOCK!! ///////////////////////////////// // missing mass off the proton // ///////////////////////////////// for(uint loc_i = 0; loc_i < protonP4.size(); loc_i++){ japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { for(uint itagm = 0; itagm < goodTAGMHits.size(); itagm++){ TLorentzVector missingP4(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4 -= protonP4[loc_i].first; MMproton_DeltaT_TAGM->Fill(protonP4[loc_i].second - goodTAGMHits[itagm].second, missingP4.M()); } for(uint itagh = 0; itagh < goodTAGHHits.size(); itagh++){ TLorentzVector missingP4(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4 -= protonP4[loc_i].first; MMproton_DeltaT_TAGH->Fill(protonP4[loc_i].second - goodTAGHHits[itagh].second, missingP4.M()); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } ////////////// // 2pi mass // ////////////// for(uint loc_j = 0; loc_j < piplusP4.size(); loc_j++){ for(uint loc_k = 0; loc_k < piminusP4.size(); loc_k++){ japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { TLorentzVector piplus_piminus = piplusP4[loc_j].first; piplus_piminus += piminusP4[loc_k].first; Mass_2pi->Fill(piplus_piminus.M()); TLorentzVector kplus_kminus = kplusP4[loc_j].first; kplus_kminus += kminusP4[loc_k].first; for(uint itagm = 0; itagm < goodTAGMHits.size(); itagm++){ TLorentzVector missingP4(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4 -= piplus_piminus; MM2pi_DeltaT_TAGM->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, missingP4.M()); M2pi_DeltaT_TAGM->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, piplus_piminus.M()); TLorentzVector missingP4_2k(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4_2k -= kplus_kminus; double deltaT = piminusP4[loc_k].second - goodTAGMHits[itagm].second; if(fabs(deltaT) < 3) { MM_M2pi_TAGM->Fill(piplus_piminus.M(), missingP4.M()); MM_M2k_TAGM->Fill(kplus_kminus.M(), missingP4_2k.M()); // 2pi mass close to rho if(piplus_piminus.M() > 0.6 && piplus_piminus.M() < 1.0){ DeltaPz_MM_TAGM->Fill(missingP4.M(), piplus_piminus.Vect().Z() - goodTAGMHits[itagm].first); // missing mass close to proton if(missingP4.M() > 0.5 && missingP4.M() < 1.2) { cout << "Event number = " << eventnumber << " found exclusive p2pi (M=" << piplus_piminus.M() << ") with TAGM" << endl; //look for awayside protons get_proton_recoil(loop, missingP4, 0); } } } } for(uint itagh = 0; itagh < goodTAGHHits.size(); itagh++){ TLorentzVector missingP4(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4 -= piplus_piminus; MM2pi_DeltaT_TAGH->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, missingP4.M()); M2pi_DeltaT_TAGH->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, piplus_piminus.M()); TLorentzVector missingP4_2k(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4_2k -= kplus_kminus; double deltaT = piminusP4[loc_k].second - goodTAGHHits[itagh].second; if(fabs(deltaT) < 5) { MM_M2pi_TAGH->Fill(piplus_piminus.M(), missingP4.M()); MM_M2k_TAGH->Fill(kplus_kminus.M(), missingP4_2k.M()); // 2pi mass close to rho if(piplus_piminus.M() > 0.6 && piplus_piminus.M() < 1.0){ DeltaPz_MM_TAGH->Fill(missingP4.M(), piplus_piminus.Vect().Z() - goodTAGHHits[itagh].first); // missing mass close to proton if(missingP4.M() > 0.5 && missingP4.M() < 1.2) { cout << "Event number = " << eventnumber << " found exclusive p2pi (M=" << piplus_piminus.M() << ") with TAGH" << endl; //look for awayside protons get_proton_recoil(loop, missingP4, 1); } } } } // exclusive proton reconstruction if(protonP4.size() == 1) { // only events with 1 tagged proton TLorentzVector ppiplus = protonP4[0].first; ppiplus += piplusP4[loc_j].first; TLorentzVector ppiminus = protonP4[0].first; ppiminus += piminusP4[loc_k].first; TLorentzVector p_piplus_piminus = piplus_piminus; p_piplus_piminus += protonP4[0].first; DalitzPlot_p2pi->Fill(piplus_piminus.M2(), ppiplus.M2()); for(uint itagm = 0; itagm < goodTAGMHits.size(); itagm++){ TLorentzVector missingP4(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4 -= p_piplus_piminus; MM2_p2pi_DeltaT_TAGM->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, missingP4.M2()); Exclusive_M2pi_DeltaT_TAGM->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, piplus_piminus.M()); double deltaT = piminusP4[loc_k].second - goodTAGMHits[itagm].second; if(fabs(deltaT) < 3) { Exclusive_MM2_M2pi_TAGM->Fill(piplus_piminus.M(), missingP4.M2()); if(fabs(missingP4.M2()) < 1) Exclusive_DalitzPlot_p2pi_TAGM->Fill(piplus_piminus.M2(), ppiplus.M2()); } } for(uint itagh = 0; itagh < goodTAGHHits.size(); itagh++){ TLorentzVector missingP4(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4 -= p_piplus_piminus; MM2_p2pi_DeltaT_TAGH->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, missingP4.M2()); Exclusive_M2pi_DeltaT_TAGH->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, piplus_piminus.M()); double deltaT = piminusP4[loc_k].second - goodTAGHHits[itagh].second; if(fabs(deltaT) < 5) { Exclusive_MM2_M2pi_TAGH->Fill(piplus_piminus.M(), missingP4.M2()); if(fabs(missingP4.M2()) < 1) Exclusive_DalitzPlot_p2pi_TAGH->Fill(piplus_piminus.M2(), ppiplus.M2()); } } } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } } ////////////// // 3pi mass // ////////////// for(uint loc_i = 0; loc_i < locFCALpi0_p4.size(); loc_i++){ for(uint loc_j = 0; loc_j < piplusP4.size(); loc_j++){ for(uint loc_k = 0; loc_k < piminusP4.size(); loc_k++){ if(loc_j == 0 && loc_k == 0) { // only fill pi0 mass once per event japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { DiPhoton_Mass_FCAL_3pi->Fill(locFCALpi0_p4[loc_i].first.M()); } japp->RootUnLock(); //RELEASE ROOT LOCK!! } if(locFCALpi0_p4[loc_i].first.M() > 0.1 && locFCALpi0_p4[loc_i].first.M() < 0.25) get_3pi(locFCALpi0_p4[loc_i],loc_j,loc_k,0,eventnumber); } } } for(uint loc_i = 0; loc_i < locBCALpi0_p4.size(); loc_i++){ for(uint loc_j = 0; loc_j < piplusP4.size(); loc_j++){ for(uint loc_k = 0; loc_k < piminusP4.size(); loc_k++){ if(loc_j == 0 && loc_k == 0) { // only fill pi0 mass once per event japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { DiPhoton_Mass_BCAL_3pi->Fill(locBCALpi0_p4[loc_i].first.M()); } japp->RootUnLock(); //RELEASE ROOT LOCK!! } if(locBCALpi0_p4[loc_i].first.M() > 0.1 && locBCALpi0_p4[loc_i].first.M() < 0.25) get_3pi(locBCALpi0_p4[loc_i],loc_j,loc_k,1,eventnumber); } } } return NOERROR; } //------------------ // tag and probe proton recoil //------------------ jerror_t JEventProcessor_physics_detcom::get_proton_recoil(JEventLoop *loop, TLorentzVector locMissing_p4, int itag){ // get tracks and detector matches vector locTrackTimeBasedVector; loop->Get(locTrackTimeBasedVector); const DDetectorMatches* locDetectorMatches = NULL; loop->GetSingle(locDetectorMatches); //////////////////////////// // loop time-based tracks // //////////////////////////// for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i){ for(size_t loc_ipiplus = 0; loc_ipiplus < piplus_chisq.size(); loc_ipiplus++) if(locTrackTimeBasedVector[loc_i]->chisq == piplus_chisq[loc_ipiplus]) continue; for(size_t loc_ipiminus = 0; loc_ipiminus < piminus_chisq.size(); loc_ipiminus++) if(locTrackTimeBasedVector[loc_i]->chisq == piminus_chisq[loc_ipiminus]) continue; const DTrackWireBased* locTrackWireBased = NULL; locTrackTimeBasedVector[loc_i]->GetSingle(locTrackWireBased); if(!locTrackWireBased) continue; // cut to get tracks from target double locPositionZ = locTrackTimeBasedVector[loc_i]->position().Z(); if(locPositionZ < 50 || locPositionZ > 80) continue; int Ndof = locTrackTimeBasedVector[loc_i]->Ndof; int wireNdof = locTrackWireBased->Ndof; Track_Ndof_Theta_NoMatch->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),Ndof); Track_WireNdof_Theta_NoMatch->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),wireNdof); // cut on match to some detector if(!locDetectorMatches->Get_IsMatchedToHit(locTrackTimeBasedVector[loc_i])) continue; Track_Stat->Fill(2); // cut on track nHits Track_Ndof_Theta->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),Ndof); Track_WireNdof_Theta->Fill(locTrackTimeBasedVector[loc_i]->momentum().Theta()*180./TMath::Pi(),wireNdof); if(wireNdof < 8 && Ndof < 8) continue; Track_Stat->Fill(3); // cut on track quality double FOM = TMath::Prob(locTrackTimeBasedVector[loc_i]->chisq,locTrackTimeBasedVector[loc_i]->Ndof); Track_FOM->Fill(FOM); if(FOM < 1e-3) continue; Track_Stat->Fill(4); // only positive tracks int locCharge = locTrackTimeBasedVector[loc_i]->charge(); if(locCharge < 1) continue; // calculate deltaPhi double locPhi = locTrackTimeBasedVector[loc_i]->momentum().Phi()*180./TMath::Pi(); double locMissingPhi = locMissing_p4.Vect().Phi()*180./TMath::Pi(); RecoilP_PhiTrack_PhiMissing[itag]->Fill(locMissingPhi, locPhi); double deltaTheta = (locTrackTimeBasedVector[loc_i]->momentum().Theta() - locMissing_p4.Vect().Theta())*180./TMath::Pi(); double deltaPhi = locPhi - locMissingPhi; if(deltaPhi > 180.) deltaPhi -= 360.; else if(deltaPhi < -180.) deltaPhi += 360.; RecoilP_DeltaPhi_DeltaTheta[itag]->Fill(deltaTheta, deltaPhi); // select deltaPhi matching recoils if(fabs(deltaPhi) > 30. || fabs(deltaTheta) > 20.) continue; // plot dE/dx to look at protons double locP = locTrackTimeBasedVector[loc_i]->momentum().Mag(); RecoilP_dEdxVsP_CDC[itag]->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_CDC*1e3*cdcMCscale); RecoilP_dEdxVsP_FDC[itag]->Fill(locP,locTrackTimeBasedVector[loc_i]->ddEdx_FDC*1e6); // plot beta from SC-to-BCAL to look at protons DSCHitMatchParams locSCHitMatchParams; DBCALShowerMatchParams locBCALShowerMatchParams; bool foundSC = dParticleID->Get_BestSCMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locSCHitMatchParams); bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locBCALShowerMatchParams); if(foundSC && foundBCAL) { double locSCTime = locSCHitMatchParams.dHitTime - locSCHitMatchParams.dFlightTime; double locSCPathLength = locSCHitMatchParams.dPathLength; const DBCALShower* locBCALShower = dynamic_cast(locBCALShowerMatchParams.dBCALShower); //dShowerObject double locBCALTime = locBCALShower->t; double locBCALPathLength = locBCALShowerMatchParams.dPathLength; double deltaPathLength = locBCALPathLength - locSCPathLength; //double deltaTProjected = (locBCALShower->t - locBCALShowerMatchParams.dFlightTime) - locSCTime; double deltaT = locBCALTime - locSCTime; double beta = deltaPathLength/(29.9792458*deltaT); RecoilP_BetaVsP_SCtoBCAL[itag]->Fill(locP, beta); } } return NOERROR; } //------------------ // 3pi algo //------------------ jerror_t JEventProcessor_physics_detcom::get_3pi(pair locPi0_p4, int loc_j, int loc_k, int locCal, int eventnumber){ TLorentzVector threepi_p4 = locPi0_p4.first; threepi_p4 += piplusP4[loc_j].first; threepi_p4 += piminusP4[loc_k].first; japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { Mass_3pi[locCal]->Fill(threepi_p4.M()); for(uint itagm = 0; itagm < goodTAGMHits.size(); itagm++){ TLorentzVector missingP4(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4 -= threepi_p4; double deltaT = piminusP4[loc_k].second - goodTAGMHits[itagm].second; if(fabs(deltaT) < 3) { MM_M3pi_TAGM[locCal]->Fill(threepi_p4.M(), missingP4.M()); double deltaT_Pi0ToTag = locPi0_p4.second - goodTAGMHits[itagm].second; double deltaT_Pi0ToTrack = locPi0_p4.second - piminusP4[loc_k].second; if(missingP4.M() > 0.5 && missingP4.M() < 1.2) { DeltaT_M3pi_TAGM[locCal]->Fill(threepi_p4.M(), deltaT_Pi0ToTag); DeltaTtrack_M3pi_TAGM[locCal]->Fill(threepi_p4.M(), deltaT_Pi0ToTrack); } // close to exclusive p omega if(missingP4.M() > 0.5 && missingP4.M() < 1.2 && threepi_p4.M() > 0.7 && threepi_p4.M() < 0.9 && locCal==0) { cout << "Event number = " << eventnumber << " found exclusive p3pi (M=" << threepi_p4.M() << ") with TAGM" << endl; } } } for(uint itagh = 0; itagh < goodTAGHHits.size(); itagh++){ TLorentzVector missingP4(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4 -= threepi_p4; double deltaT = piminusP4[loc_k].second - goodTAGHHits[itagh].second; if(fabs(deltaT) < 5) { MM_M3pi_TAGH[locCal]->Fill(threepi_p4.M(), missingP4.M()); double deltaT_Pi0ToTag = locPi0_p4.second - goodTAGHHits[itagh].second; double deltaT_Pi0ToTrack = locPi0_p4.second - piminusP4[loc_k].second; if(missingP4.M() > 0.5 && missingP4.M() < 1.2) { DeltaT_M3pi_TAGH[locCal]->Fill(threepi_p4.M(), deltaT_Pi0ToTag); DeltaTtrack_M3pi_TAGH[locCal]->Fill(threepi_p4.M(), deltaT_Pi0ToTrack); } // close to exclusive p omega if(missingP4.M() > 0.5 && missingP4.M() < 1.2 && threepi_p4.M() > 0.7 && threepi_p4.M() < 0.9 && locCal==0) { cout << "Event number = " << eventnumber << " found exclusive p3pi (M=" << threepi_p4.M() << ") with TAGH" << endl; } } } // exclusive proton reconstruction if(protonP4.size() == 1) { // only events with 1 tagged proton TLorentzVector p3pi_p4 = protonP4[0].first; p3pi_p4 += threepi_p4; for(uint itagm = 0; itagm < goodTAGMHits.size(); itagm++){ TLorentzVector missingP4(0,0,goodTAGMHits[itagm].first,goodTAGMHits[itagm].first+0.938); missingP4 -= p3pi_p4; MM2_p3pi_DeltaT_TAGM[locCal]->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, missingP4.M2()); M3pi_DeltaT_TAGM[locCal]->Fill(piminusP4[loc_k].second - goodTAGMHits[itagm].second, threepi_p4.M()); double deltaT = piminusP4[loc_k].second - goodTAGMHits[itagm].second; if(fabs(deltaT) < 3) MM2_M3pi_TAGM[locCal]->Fill(threepi_p4.M(), missingP4.M2()); } for(uint itagh = 0; itagh < goodTAGHHits.size(); itagh++){ TLorentzVector missingP4(0,0,goodTAGHHits[itagh].first,goodTAGHHits[itagh].first+0.938); missingP4 -= p3pi_p4; MM2_p3pi_DeltaT_TAGH[locCal]->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, missingP4.M2()); M3pi_DeltaT_TAGH[locCal]->Fill(piminusP4[loc_k].second - goodTAGHHits[itagh].second, threepi_p4.M()); double deltaT = piminusP4[loc_k].second - goodTAGHHits[itagh].second; if(fabs(deltaT) < 5) MM2_M3pi_TAGH[locCal]->Fill(threepi_p4.M(), missingP4.M2()); } } } japp->RootUnLock(); //RELEASE ROOT LOCK!! return NOERROR; } //------------------ // get tagger hits to correlate with downstream events //------------------ jerror_t JEventProcessor_physics_detcom::get_tagger_hits(JEventLoop *loop, vector< pair > &locGoodTAGMHits, vector< pair > &locGoodTAGHHits) { // TAGH hits vector locTAGHHitVector; loop->Get(locTAGHHitVector); for(uint i = 0; i < locTAGHHitVector.size(); i++) { double locTAGHTime = locTAGHHitVector[i]->t; //time_fadc; //if(locTAGHTime == 0.) continue; double locTAGHE = locTAGHHitVector[i]->E; #if 0 // check quality of tagger hit (from Nathan) const DTAGHDigiHit* locTAGHDigiHit; locTAGHHitVector[i]->GetSingle(locTAGHDigiHit); if(!locTAGHDigiHit) continue; const Df250PulseIntegral* locPulseIntegral; locTAGHDigiHit->GetSingle(locPulseIntegral); if(!locPulseIntegral) continue; double ped = locPulseIntegral->pedestal; if(ped == 0.) continue; // some problem with this double pulse_integral = locTAGHDigiHit->pulse_integral - locTAGHDigiHit->nsamples_integral*ped; #endif double pulse_integral = 10000; if(pulse_integral > 1000.) { locGoodTAGHHits.push_back(make_pair(locTAGHE, locTAGHTime)); } } // TAGM hits vector locTAGMHitVector; loop->Get(locTAGMHitVector); for(size_t k = 0; k < locTAGMHitVector.size(); k++){ double locTAGMTime = locTAGMHitVector[k]->t; //time_fadc; //if(locTAGMTime == 0.) continue; double locTAGME = locTAGMHitVector[k]->E; locGoodTAGMHits.push_back(make_pair(locTAGME, locTAGMTime)); } return NOERROR; } //------------------ // Calculate dE/dx (copy of DParticleID::CalcDCdEdx with different truncated mean settings) //------------------ jerror_t JEventProcessor_physics_detcom::calcdEdx(const DTrackTimeBased *locTrackTimeBased, const vector& locdEdxHits_CDC, const vector& locdEdxHits_FDC, double& locdEdx_FDC, double& locdEdx_CDC, double truncMeanThresh) { locdEdx_CDC = 0.0; double locdx_CDC = 0.0; unsigned int locNumHitsUsedFordEdx_CDC = int(locdEdxHits_CDC.size()*truncMeanThresh); if(locNumHitsUsedFordEdx_CDC > 0){ for(unsigned int loc_i = 0; loc_i < locNumHitsUsedFordEdx_CDC; ++loc_i){ locdEdx_CDC += locdEdxHits_CDC[loc_i].dE; //weight is ~ #e- (scattering sites): dx! locdx_CDC += locdEdxHits_CDC[loc_i].dx; } locdEdx_CDC /= locdx_CDC; } //cout << "in calcdEdx locdEdx_CDC = " << locdEdx_CDC << " locdx_CDC = " << locdx_CDC << endl; locdEdx_FDC = 0.0; double locdx_FDC = 0.0; unsigned int locNumHitsUsedFordEdx_FDC = int(locdEdxHits_FDC.size()*truncMeanThresh); if(locNumHitsUsedFordEdx_FDC > 0){ for(unsigned int loc_i = 0; loc_i < locNumHitsUsedFordEdx_FDC; ++loc_i){ locdEdx_FDC += locdEdxHits_FDC[loc_i].dE; //weight is ~ #e- (scattering sites): dx! locdx_FDC += locdEdxHits_FDC[loc_i].dx; } locdEdx_FDC /= locdx_FDC; //weight is dx/dx_total } //cout << "in calcdEdx locdEdx_FDC = " << locdEdx_FDC << " locdx_FDC = " << locdx_FDC << endl; return NOERROR; } //------------------ // BCAL pi0 algo //------------------ jerror_t JEventProcessor_physics_detcom::get_bcal_pi0(JEventLoop *loop, vector< pair > &locBCALpi0_p4) { const DDetectorMatches* locDetectorMatches = NULL; loop->GetSingle(locDetectorMatches); /////////////////////////////// // BCAL pi0 finder from Will // /////////////////////////////// vector locBCALShowerVector; loop->Get(locBCALShowerVector); double minClusterE = 0.5; // loop over BCAL showers for(int i = 0; i < (int)locBCALShowerVector.size()-1; ++i){ const DBCALShower *locBCALShower1 = locBCALShowerVector[i]; if(locBCALShower1->E < minClusterE) continue; if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower1)) continue; for(int j = i+1; j < (int)locBCALShowerVector.size(); ++j){ const DBCALShower *locBCALShower2 = locBCALShowerVector[j]; if(locBCALShower2->E < minClusterE) continue; if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower2)) continue; // cut time difference double deltaT = fabs(locBCALShower1->t - locBCALShower2->t); if(deltaT > 3) continue; // calculate invariant mass TVector3 locGamma1_p3 = TVector3(locBCALShower1->x, locBCALShower1->y, locBCALShower1->z); locGamma1_p3.SetZ(locGamma1_p3.Z() - 65.); double timeShower1 = locBCALShower1->t - locGamma1_p3.Mag()/29.9792458; locGamma1_p3.SetMag(locBCALShower1->E); TLorentzVector locGamma1_p4 = TLorentzVector(locGamma1_p3, locBCALShower1->E); TVector3 locGamma2_p3 = TVector3(locBCALShower2->x, locBCALShower2->y, locBCALShower2->z); locGamma2_p3.SetZ(locGamma2_p3.Z() - 65.); double timeShower2 = locBCALShower2->t - locGamma2_p3.Mag()/29.9792458; locGamma2_p3.SetMag(locBCALShower2->E); TLorentzVector locGamma2_p4 = TLorentzVector(locGamma2_p3, locBCALShower2->E); double locAvgT = (timeShower1 + timeShower2)/2.; TLorentzVector locDiPhoton_p4 = locGamma1_p4 + locGamma2_p4; // fill histos for invariant mass japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { DiPhoton_Mass_BCAL->Fill(locDiPhoton_p4.M()); locBCALpi0_p4.push_back(make_pair(locDiPhoton_p4,locAvgT)); if(locDiPhoton_p4.M() > 0.08 && locDiPhoton_p4.M() < 0.16){ Pi0_BCAL_Energy->Fill(locDiPhoton_p4.E()); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } } return NOERROR; } //------------------ // FCAL pi0 algo //------------------ jerror_t JEventProcessor_physics_detcom::get_fcal_pi0(JEventLoop *loop, vector< pair > &locFCALpi0_p4) { /////////////////////////////// // FCAL pi0 finder from Matt // /////////////////////////////// vector locFCALShowerVector; loop->Get(locFCALShowerVector); // event level cuts double FCAL_Etot = 0.; for(int i = 0; i < (int)locFCALShowerVector.size(); ++i) FCAL_Etot += locFCALShowerVector[i]->getEnergy(); //if(FCAL_Etot > 7.) return NOERROR; //if(locFCALShowerVector.size() >= 8 || locFCALShowerVector.size() < 2) return NOERROR; // loop over FCAL clusters for(int i = 0; i < (int)locFCALShowerVector.size()-1; ++i){ const DFCALShower *locFCALShower1 = locFCALShowerVector[i]; if(locFCALShower1->getPosition().Perp() < 20.) continue; // reject electrons near beamline for(int j = i+1; j < (int)locFCALShowerVector.size(); ++j){ const DFCALShower *locFCALShower2 = locFCALShowerVector[j]; if(locFCALShower2->getPosition().Perp() < 20.) continue; // reject electrons near beamline // cut time difference double deltaT = fabs(locFCALShower1->getTime() - locFCALShower2->getTime()); if(deltaT > 10) continue; // cut on lowest cluster energy double lowShowerEnergy = TMath::Min(locFCALShower1->getEnergy(), locFCALShower2->getEnergy()); if(lowShowerEnergy < 0.5) continue; // calculate invariant mass TVector3 locGamma1_p3 = locFCALShower1->getPosition(); double timeGamma1 = locFCALShower1->getTime() - locGamma1_p3.Mag()/29.9792458; locGamma1_p3.SetMag(locFCALShower1->getEnergy()); TLorentzVector locGamma1_p4 = TLorentzVector(locGamma1_p3, locFCALShower1->getEnergy()); TVector3 locGamma2_p3 = locFCALShower2->getPosition(); double timeGamma2 = locFCALShower2->getTime() - locGamma2_p3.Mag()/29.9792458; locGamma2_p3.SetMag(locFCALShower2->getEnergy()); TLorentzVector locGamma2_p4 = TLorentzVector(locGamma2_p3, locFCALShower2->getEnergy()); double locAvgT = (timeGamma1 + timeGamma2)/2.; TLorentzVector locDiPhoton_p4 = locGamma1_p4 + locGamma2_p4; // fill histos for invariant mass japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { DiPhoton_Mass_FCAL->Fill(locDiPhoton_p4.M()); locFCALpi0_p4.push_back(make_pair(locDiPhoton_p4, locAvgT)); if(locDiPhoton_p4.M() > 0.05 && locDiPhoton_p4.M() < 0.15){ Pi0_FCAL_Energy->Fill(locDiPhoton_p4.E()); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } } #if 0 /////////////////////////////// // FCAL pi0 finder from Matt // /////////////////////////////// vector locFCALClusterVector; loop->Get(locFCALClusterVector); // event level cuts double FCAL_Etot = 0.; for(int i = 0; i < (int)locFCALClusterVector.size(); ++i) FCAL_Etot += locFCALClusterVector[i]->getEnergy(); //if(FCAL_Etot > 7.) return NOERROR; //if(locFCALClusterVector.size() >= 8 || locFCALClusterVector.size() < 2) return NOERROR; // loop over FCAL clusters for(int i = 0; i < (int)locFCALClusterVector.size()-1; ++i){ const DFCALCluster *locFCALCluster1 = locFCALClusterVector[i]; if(locFCALCluster1->getCentroid().Perp() < 20.) continue; // reject electrons near beamline for(int j = i+1; j < (int)locFCALClusterVector.size(); ++j){ const DFCALCluster *locFCALCluster2 = locFCALClusterVector[j]; if(locFCALCluster2->getCentroid().Perp() < 20.) continue; // reject electrons near beamline // cut time difference double deltaT = fabs(locFCALCluster1->getTime() - locFCALCluster2->getTime()); if(deltaT > 10) continue; // cut on lowest cluster energy double lowClusterEnergy = TMath::Min(locFCALCluster1->getEnergy(), locFCALCluster2->getEnergy()); if(lowClusterEnergy < 0.5) continue; // calculate invariant mass TVector3 locGamma1_p3 = locFCALCluster1->getCentroid(); locGamma1_p3.SetZ(locGamma1_p3.Z() - 65.); double timeGamma1 = locFCALCluster1->getTime() - locGamma1_p3.Mag()/29.9792458; locGamma1_p3.SetMag(locFCALCluster1->getEnergy()); TLorentzVector locGamma1_p4 = TLorentzVector(locGamma1_p3, locFCALCluster1->getEnergy()); TVector3 locGamma2_p3 = locFCALCluster2->getCentroid(); locGamma2_p3.SetZ(locGamma2_p3.Z() - 65.); double timeGamma2 = locFCALCluster2->getTime() - locGamma2_p3.Mag()/29.9792458; locGamma2_p3.SetMag(locFCALCluster2->getEnergy()); TLorentzVector locGamma2_p4 = TLorentzVector(locGamma2_p3, locFCALCluster2->getEnergy()); double locAvgT = (timeGamma1 + timeGamma2)/2.; TLorentzVector locDiPhoton_p4 = locGamma1_p4 + locGamma2_p4; // fill histos for invariant mass japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { DiPhoton_Mass_FCAL->Fill(locDiPhoton_p4.M()); locFCALpi0_p4.push_back(make_pair(locDiPhoton_p4, locAvgT)); if(locDiPhoton_p4.M() > 0.1 && locDiPhoton_p4.M() < 0.25){ Pi0_FCAL_Energy->Fill(locDiPhoton_p4.E()); } } japp->RootUnLock(); //RELEASE ROOT LOCK!! } } #endif return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_physics_detcom::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_physics_detcom::fini(void) { // Called before program exit after event processing is finished. return NOERROR; }