// $Id$ // // File: JEventProcessor_single_eta_resolution.cc // Created: Fri Jul 14 16:58:12 EDT 2017 // Creator: dalton (on Linux gluon106.jlab.org 2.6.32-642.3.1.el6.x86_64 x86_64) // #include "JEventProcessor_single_eta_resolution.h" using namespace jana; #include // std::min #include "HistogramTools.h" // Creat and fill histograms tread safe #include "PID/DNeutralShower.h" #include "TRACKING/DMCThrown.h" #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_single_eta_resolution()); } } // "C" //------------------ // JEventProcessor_single_eta_resolution (Constructor) //------------------ JEventProcessor_single_eta_resolution::JEventProcessor_single_eta_resolution() { } //------------------ // ~JEventProcessor_single_eta_resolution (Destructor) //------------------ JEventProcessor_single_eta_resolution::~JEventProcessor_single_eta_resolution() { } //------------------ // init //------------------ jerror_t JEventProcessor_single_eta_resolution::init(void) { // This is called once at program startup. gStyle->SetOptStat(1110); // set style gStyle->SetTitleOffset(1, "Y"); gStyle->SetTitleSize(0.05,"xyz"); gStyle->SetTitleSize(0.06,"h"); gStyle->SetLabelSize(0.05,"xyz"); gStyle->SetTitleX(0); gStyle->SetTitleAlign(13); gStyle->SetNdivisions(505,"xy"); gRandom = new TRandom3(); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_single_eta_resolution::brun(JEventLoop *eventLoop, int32_t runnumber) { // This is called whenever the run number changes return NOERROR; } void Print4vec(TLorentzVector vec) { //printf("theta=%5.1f deg, M=%8.2E\n", printf("theta=%5.1f deg, M=%.3f\n", vec.Theta()*57.2957, vec.M()); } //------------------ // evnt //------------------ jerror_t JEventProcessor_single_eta_resolution::evnt(JEventLoop *loop, uint64_t eventnumber) { // This is called for every event. Use of common resources like writing // to a file or filling a histogram should be mutex protected. Using // loop->Get(...) to get reconstructed objects (and thereby activating the // reconstruction algorithm) should be done outside of any mutex lock // since multiple threads may call this method at the same time. // Here's an example: // // vector mydataclasses; // loop->Get(mydataclasses); // // japp->RootFillLock(this); // ... fill historgrams or trees ... // japp->RootFillUnLock(this); vector myshower; loop->Get(myshower); vector mythrown; loop->Get(mythrown); char title[255]; sprintf(title, "Num thrown;num thrown"); Fill1DHistogram("singleeta", "basic", "num_thrown", mythrown.size(), title, 15, -0.5, 14.5); for (uint i=0; iPID(), title, 15, -0.5, 14.5, 20, -0.5, 19.5); } sprintf(title, "Num showers;num showers"); Fill1DHistogram("singleeta", "basic", "num_showers", myshower.size(), title, 15, -0.5, 14.5); if (mythrown.size() == 3) { if (myshower.size()!=2) return NOERROR; // ** Plot the thrown eta mass TLorentzVector thrownEta_P4 = mythrown[0]->lorentzMomentum(); TLorentzVector thrown1_P4 = mythrown[1]->lorentzMomentum(); TLorentzVector thrown2_P4 = mythrown[2]->lorentzMomentum(); // TLorentzVector thrownEtaP4 = thrown1_P4 + thrown2_P4; // sprintf(title, "thrown eta mass;M_{#eta->#gamma#gamma}"); // Fill1DHistogram("singleeta", "basic", "dHist_EtaMass_thrown", // thrownEtaP4.M(), title, 200, 0.2, 0.8); // **** Use the thrown values and manually broaden them. // **** Need to comment out later part where 4 vectors are created // float m1 = gRandom->Gaus(1, sqrt(0.0025/thrown1_P4.E()+0.000625)); // float m2 = gRandom->Gaus(1, sqrt(0.0025/thrown2_P4.E()+0.000625)); // sprintf(title, ";Photon P (GeV);modifier"); // Fill2DHistogram("singleeta", "basic", "random", // thrown1_P4.E(), m1, title, 100, 0, 10, 100, 0, 2); // Fill2DHistogram("singleeta", "basic", "random", // thrown2_P4.E(), m2, title, 100, 0, 10, 100, 0, 2); // TLorentzVector Photon1P4 = thrown1_P4; // TLorentzVector Photon2P4 = thrown2_P4; // Photon1P4.SetRho(thrown1_P4.E()*m1); // Convert position 3-vector to momentum // Photon1P4.SetE(thrown1_P4.E()*m1); // Set energy // Photon2P4.SetRho(thrown2_P4.E()*m2); // Convert position 3-vector to momentum // Photon2P4.SetE(thrown2_P4.E()*m2); // Set energy // Create detected 4-vectors TLorentzVector Photon1X4 = myshower[0]->dSpacetimeVertex; TLorentzVector Photon2X4 = myshower[1]->dSpacetimeVertex; double Photon1E = myshower[0]->dEnergy; double Photon2E = myshower[1]->dEnergy; TLorentzVector Photon1P4 = Photon1X4; TLorentzVector Photon2P4 = Photon2X4; Photon1P4.SetPz(Photon1P4.Pz() - 59.46); // subtract target position Photon1P4.SetRho(Photon1E); // Convert position 3-vector to momentum Photon1P4.SetE(Photon1E); // Set energy Photon2P4.SetPz(Photon2P4.Pz() - 59.46); Photon2P4.SetRho(Photon2E); Photon2P4.SetE(Photon2E); TLorentzVector locEtaP4 = Photon1P4 + Photon2P4; float etaP = locEtaP4.P(); float etaTheta = locEtaP4.Theta(); sprintf(title, "eta mass;M_{#eta->#gamma#gamma}"); Fill1DHistogram("singleeta", "basic", "dHist_EtaMass", locEtaP4.M(), title, 200, 0.2, 0.8); if (Photon1E<1 && Photon2E<1) { sprintf(title, "eta mass;M_{#eta->#gamma#gamma}"); Fill1DHistogram("singleeta", "basic", "dHist_EtaMass_1GeV", locEtaP4.M(), title, 200, 0.2, 0.8); } float PhotonE_small = min(Photon1E,Photon2E); float PhotonE_large = max(Photon1E,Photon2E); float Photon1Theta = Photon1P4.Theta()*57.2957; float Photon2Theta = Photon2P4.Theta()*57.2957; float PhotonTheta_small = min(Photon1P4.Theta(),Photon2P4.Theta())*57.2957; float PhotonTheta_large = max(Photon1P4.Theta(),Photon2P4.Theta())*57.2957; float PhotonMeanE = (Photon1E+Photon2E)/2; sprintf(title, ";Photon 1 P (GeV);Photon 2 P (GeV)"); Fill2DHistogram("singleeta", "basic", "g1P_g2P", PhotonE_small, PhotonE_large, title, 100, 0, 10, 100, 0, 10); sprintf(title, ";Photon 1 Theta (deg);Photon 2 Theta (deg)"); Fill2DHistogram("singleeta", "basic", "g1Theta_g2Theta", PhotonTheta_small, PhotonTheta_large, title, 130, 0, 130, 130, 0, 130); sprintf(title, ";Photon P (GeV);Photon Theta (deg)"); Fill2DHistogram("singleeta", "basic", "gP_gTheta", Photon1E, Photon1Theta, title, 100, 0, 10, 130, 0, 130); Fill2DHistogram("singleeta", "basic", "gP_gTheta", Photon2E, Photon2Theta, title, 100, 0, 10, 130, 0, 130); sprintf(title, ";P_{#eta} (GeV);#theta_{#eta} (deg)"); Fill2DHistogram("singleeta", "basic", "etaP_etaTheta", etaP, etaTheta, title, 100, 0, 10, 130, 0, 130); // ******************* // eta mass resolution sprintf(title, "eta mass;Photon P (GeV);M_{#eta->#gamma#gamma} (MeV)"); Fill2DHistogram("singleeta", "resolution", "dHist_EtaMass_gP", Photon1E, locEtaP4.M(), title, 120, 0, 6, 50, 0.3, 0.8); Fill2DHistogram("singleeta", "resolution", "dHist_EtaMass_gP", Photon2E, locEtaP4.M(), title, 120, 0, 6, 50, 0.3, 0.8); // E symmetrical photons if (abs(Photon1E-Photon2E)<0.1) { sprintf(title, "eta mass (E_{#gamma_{1}} #approx E_{#gamma_{2}});Photon P (GeV);M_{#eta->#gamma#gamma} (MeV)"); Fill2DHistogram("singleeta", "resolution", "dHist_EtaMass_gPsym100", PhotonMeanE, locEtaP4.M(), title, 120, 0, 6, 50, 0.3, 0.8); } if (abs(Photon1E-Photon2E)<0.5) { sprintf(title, "eta mass (E_{#gamma_{1}} #approx E_{#gamma_{2}});Photon P (GeV);M_{#eta->#gamma#gamma} (MeV)"); Fill2DHistogram("singleeta", "resolution", "dHist_EtaMass_gPsym500", PhotonMeanE, locEtaP4.M(), title, 120, 0, 6, 50, 0.3, 0.8); } // Theta symmetrical photons if (abs(Photon1Theta-Photon2Theta)<10) { sprintf(title, "eta mass;Photon Theta (deg);M_{#eta->#gamma#gamma} (MeV)"); Fill2DHistogram("singleeta", "resolution", "dHist_EtaMass_gThetasym", (Photon1Theta+Photon2Theta)/2, locEtaP4.M(), title, 130, 0, 130, 50, 0.3, 0.8); } // ************************ // photon energy resolution // Match thrown and detected and calculate difference in energy //float dist1 = abs(mythrown[1]->energy() - Photon1E) + abs(mythrown[2]->energy() - Photon2E); //float dist2 = abs(mythrown[2]->energy() - Photon1E) + abs(mythrown[1]->energy() - Photon2E); float dist1 = abs((Photon1P4-thrown1_P4).Mag() + (Photon2P4-thrown2_P4).Mag()); float dist2 = abs((Photon2P4-thrown1_P4).Mag() + (Photon1P4-thrown2_P4).Mag()); if (Photon1E<1 && eventnumber<100) { //Photon1X4.Print(); //Photon2X4.Print(); //Photon1P4.Print(); //Photon2P4.Print(); thrown1_P4.Print(); thrown2_P4.Print(); locEtaP4.Print(); thrownEta_P4.Print(); Print4vec(Photon1P4); Print4vec(Photon2P4); Print4vec(locEtaP4); Print4vec(thrownEta_P4); // printf("%f %f %f \n",PhotonMeanE,dist1,dist2); //printf("%f %f\n",m1,m2); printf("\n"); } float thdist11 = mythrown[1]->momentum().Theta()*57.2957 - Photon1Theta; float thdist22 = mythrown[2]->momentum().Theta()*57.2957 - Photon2Theta; float thdist21 = mythrown[2]->momentum().Theta()*57.2957 - Photon1Theta; float thdist12 = mythrown[1]->momentum().Theta()*57.2957 - Photon2Theta; float Photon1Ediff, Photon2Ediff; float Photon1Thetadiff, Photon2Thetadiff; float Photon1Thetadiff_alt, Photon2Thetadiff_alt; if (dist1 < dist2) { Photon1Ediff = mythrown[1]->energy() - Photon1E; Photon2Ediff = mythrown[2]->energy() - Photon2E; Photon1Thetadiff = thdist11; Photon2Thetadiff = thdist22; Photon1Thetadiff_alt = thdist12; Photon2Thetadiff_alt = thdist21; } else { Photon1Ediff = mythrown[2]->energy() - Photon1E; Photon2Ediff = mythrown[1]->energy() - Photon2E; Photon1Thetadiff = thdist21; Photon2Thetadiff = thdist12; Photon1Thetadiff_alt = thdist11; Photon2Thetadiff_alt = thdist22; } float distsmall = min(dist1,dist2); float distlarge = max(dist1,dist2); sprintf(title, "photon E matching;distsmall;distlarge"); Fill2DHistogram("singleeta", "resolution", "E_match", distsmall, distlarge, title, 100, 0, 2, 100, 0, 2); sprintf(title, "photon E matching;distsmall;distlarge"); Fill2DHistogram("singleeta", "resolution", "match_theta", Photon1Thetadiff, Photon1Thetadiff_alt, title, 130, 0, 130, 130, 0, 130); Fill2DHistogram("singleeta", "resolution", "match_theta", Photon2Thetadiff, Photon2Thetadiff_alt, title, 130, 0, 130, 130, 0, 130); sprintf(title, "photon E resolution;E_{#gamma} (GeV);(E_{thrown} - E_{shower})/E_{shower}"); if (Photon1E>0.05) Fill2DHistogram("singleeta", "resolution", "gPres_gP", Photon1E, Photon1Ediff/Photon1E, title, 120, 0, 6, 200, -0.5, 0.5); if (Photon2E>0.05) Fill2DHistogram("singleeta", "resolution", "gPres_gP", Photon2E, Photon2Ediff/Photon2E, title, 120, 0, 6, 200, -0.5, 0.5); sprintf(title, "photon E resolution;E_{thrown1} - E_{shower1} (GeV);E_{thrown2} - E_{shower2} (GeV)"); Fill2DHistogram("singleeta", "resolution", "gPres_gPres", Photon1Ediff, Photon2Ediff, title, 200, -0.5, 0.5, 200, -0.5, 0.5); sprintf(title, "photon #theta resolution;#theta_{thrown1} - #theta_{shower1} (GeV);#theta_{thrown2} - #theta_{shower2} (GeV)"); Fill2DHistogram("singleeta", "resolution", "gThetares_gThetares", Photon1Thetadiff, Photon2Thetadiff, title, 200, -10, 10, 200, -10, 10); sprintf(title, "photon #theta resolution;#theta_{shower} (GeV);#theta_{thrown} - #theta_{shower} (GeV)"); Fill2DHistogram("singleeta", "resolution", "gThetares_gTheta", Photon1Theta, Photon1Thetadiff, title, 130, 0, 130, 200, -20, 20); Fill2DHistogram("singleeta", "resolution", "gThetares_gTheta", Photon2Theta, Photon2Thetadiff, title, 130, 0, 130, 200, -20, 20); } return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_single_eta_resolution::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_single_eta_resolution::fini(void) { // Called before program exit after event processing is finished. SortDirectories(); // Sort the histogram directories by name return NOERROR; }