// $Id$ // // File: JEventProcessor_AccPi0.cc // Created: Ср. янв. 30 16:59:40 CET 2013 // Creator: alexander (on Linux altro.site 3.4.11-2.16-desktop x86_64) // #include "JEventProcessor_AccPi0.h" using namespace jana; // #include #include // Routine used to create our JEventProcessor #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_AccPi0()); } } // "C" //------------------ // JEventProcessor_AccPi0 (Constructor) //------------------ JEventProcessor_AccPi0::JEventProcessor_AccPi0() { } //------------------ // ~JEventProcessor_AccPi0 (Destructor) //------------------ JEventProcessor_AccPi0::~JEventProcessor_AccPi0() { } //------------------ // init //------------------ jerror_t JEventProcessor_AccPi0::init(void) { // This is called once at program startup // Output tree // fTreeAccPi0 = new TTree("AccPi0", "Reconstruction pi^{0} acceptance"); // All tracks // fTreeAccPi0->Branch("Ntrack", &fNTrack); fHistNtrack = new TH1D("Ntrack", "Number of tracks", 10, 0, 10); fHistTypeGen = new TH1D("Type", "Type of the generated tracks", 10, 0, 10); // Gamma fHistNGamma = new TH1D("NGamma", "Number of gammas", 10, 0, 10); fHistNGammaGen = (TH1D*)fHistNGamma->Clone("NGammaGen"); fHistThetaGammaGen = new TH1D("ThetaGammaGen", "Theta of gammas from Pi0", 400., 0., 15.); fHistThetaGamma = (TH1D*)fHistThetaGammaGen->Clone("ThetaGamma"); fHistThetaGammaGenDet = (TH1D*)fHistThetaGammaGen->Clone("ThetaGammaGenDet"); fHistPhiGammaGen = new TH1D("PhiGammaGen", "Phi of gammas from Pi0", 360., 0., 360.); fHistPhiGamma = (TH1D*)fHistPhiGammaGen->Clone("PhiGamma"); fHistPhiGammaGenDet = (TH1D*)fHistPhiGammaGen->Clone("PhiGammaGenDet"); fHistEGammaGen = new TH1D("EGammaGen", "E of gammas from Pi0", 100., 0., 5.5); fHistEGamma = (TH1D*)fHistEGammaGen->Clone("EGamma"); fHistEGammaGenDet = (TH1D*)fHistEGammaGen->Clone("EGammaGenDet"); fHistFOMgamma = new TH1D("FOMgamma", "FOM for gammas", 1000., 0., 100000000000.); fHistEGammaVSThetaGamma = new TH2D("EGammaVSThetaGamma", "Eg VS Theta",400.,0.,15.,100.,0.,5.5); fHistPhiGammaVSThetaGamma = new TH2D("PhiGammaVSThetaGamma", "Phi VS Theta",400.,0.,15.,360.,0.,360.); fHistPhiGammaGenVSThetaGammaGen = new TH2D("PhiGammaGenVSThetaGammaGen", "PhiGen VS ThetaGen",400.,0.,15.,360.,0.,360.); // Pi0 fHistNPi0 = new TH1D("Npi0", "Number of #pi^0", 10, 0, 10); fHistNPi0Gen = (TH1D*)fHistNPi0->Clone("Npi0Gen"); fHistThetaPi0Gen = new TH1D("ThetaPi0Gen", "Theta of Pi0", 100., 0., 2.); fHistThetaPi0 = (TH1D*)fHistThetaPi0Gen->Clone("ThetaPi0"); fHistThetaPi0GenDet = (TH1D*)fHistThetaPi0Gen->Clone("ThetaPi0GenDet"); fHistPhiPi0Gen = new TH1D("PhiPi0Gen", "Phi of Pi0", 360., 0., 360.); fHistPhiPi0 = (TH1D*)fHistPhiPi0Gen->Clone("PhiPi0"); fHistPhiPi0GenDet = (TH1D*)fHistPhiPi0Gen->Clone("PhiPi0GenDet"); fHistEPi0Gen = new TH1D("EPi0Gen", "E of Pi0", 300., 3.5, 6.5); fHistEPi0 = (TH1D*)fHistEPi0Gen->Clone("EPi0"); fHistEPi0GenDet = (TH1D*)fHistEPi0Gen->Clone("EPi0GenDet"); fHistMgg = new TH1D("Mgg", "gg invariant mass", 200., 0., 0.2); fHistMggNoCut = (TH1D*)fHistMgg->Clone("MggNoCut"); fHistMggNoCutAll = (TH1D*)fHistMgg->Clone("MggNoCutAll"); fHistMggNoCutVsThetaPi0 = new TH2D("MggNoCutVsThetaPi0", "gg invariant mass vs theta of pi0", 20., 0., 2., 200., 0., 0.2); fHistPhiPi0VSThetaPi0 = new TH2D("PhiPi0VSThetaPi0", "Phi VS Theta",20.,0.,2.,360.,0.,360.); fHistPhiPi0ShiftVSThetaPi0 = (TH2D*)fHistPhiPi0VSThetaPi0->Clone("PhiPi0ShiftVSThetaPi0"); fHistPhiPi0GenVSThetaPi0Gen = (TH2D*)fHistPhiPi0VSThetaPi0->Clone("PhiPi0GenVSThetaPi0Gen"); fHistPhiPi0GenDetVSThetaPi0GenDet = (TH2D*)fHistPhiPi0VSThetaPi0->Clone("PhiPi0GenDetVSThetaPi0GenDet"); // Other fHistNOther = new TH1D("NOther", "Number of other particles", 10, 0, 10); fHistNOtherGen = (TH1D*)fHistNOther->Clone("NOtherGen"); fHistThetaOtherGen = new TH1D("ThetaOtherGen", "Theta of others", 180., 0., 180.); fHistThetaOther = (TH1D*)fHistThetaOtherGen->Clone("ThetaOther"); fHistPhiOtherGen = new TH1D("PhiOtherGen", "Phi of others", 360., 0., 360.); fHistPhiOther = (TH1D*)fHistPhiOtherGen->Clone("PhiOther"); fHistEOtherGen = new TH1D("EOtherGen", "E of others", 100., 0., 6.); fHistEOther = (TH1D*)fHistEOtherGen->Clone("EOther"); // return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_AccPi0::brun(JEventLoop *eventLoop, int runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_AccPi0::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. // 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. // ResetEvnt(); // // static const Double_t kMpi0 = 0.13498; static const Double_t kMpi0 = 0.14204; // // Generated // vector throwns; loop->Get(throwns); DLorentzVector q4pi0Gen, q4gGen[2]; for (unsigned int i=0; iFill(throwns[i]->type); DLorentzVector q4 = throwns[i]->lorentzMomentum(); TH1D *hTheta, *hPhi, *hE; switch (throwns[i]->type) { case 1: hTheta = fHistThetaGammaGen; hPhi = fHistPhiGammaGen; hE = fHistEGammaGen; q4gGen[fNGammaGen] = q4; fHistPhiGammaGenVSThetaGammaGen->Fill(q4.Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4.Phi())*TMath::RadToDeg()); ++fNGammaGen; break; case 7: ++fNPi0Gen; hTheta = fHistThetaPi0Gen; hPhi = fHistPhiPi0Gen; hE = fHistEPi0Gen; q4pi0Gen = q4; fHistPhiPi0GenVSThetaPi0Gen->Fill(q4.Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4.Phi())*TMath::RadToDeg()); break; default: ++fNOtherGen; hTheta = fHistThetaOtherGen; hPhi = fHistPhiOtherGen; hE = fHistEOtherGen; } hTheta->Fill(q4.Theta()*TMath::RadToDeg()); hPhi->Fill(TVector2::Phi_0_2pi(q4.Phi())*TMath::RadToDeg()); hE->Fill(q4.T()); } fHistNGammaGen->Fill(fNGammaGen); fHistNPi0Gen->Fill(fNPi0Gen); fHistNOtherGen->Fill(fNOtherGen); // Select only pi0->g+g if (fNGammaGen != 2 || fNPi0Gen != 1) return NOERROR; // // Reconstructed // vector tracks; loop->Get(tracks); fNtrack = tracks.size(); fHistNtrack->Fill(fNtrack); DLorentzVector q4pi0, q4g[2]; Double_t diffMin = 1000.; for (Int_t i=0; iGet_Hypothesis(Gamma); if (!g1) continue; fHistFOMgamma->Fill(g1->dFOM); for (Int_t j=i+1; jGet_Hypothesis(Gamma); if (!g2) continue; fHistFOMgamma->Fill(g2->dFOM); DLorentzVector q = g1->lorentzMomentum() + g2->lorentzMomentum(); Double_t diff = TMath::Abs(q.M()-kMpi0); if (diff < diffMin) { q4pi0 = q; q4g[0] = g1->lorentzMomentum(); q4g[1] = g2->lorentzMomentum(); if (q4pi0.M() >= 0.125 && q4pi0.M() <= 0.17) ++fNPi0; diffMin = diff; } fHistMggNoCutAll->Fill(q.M()); } } if (diffMin != 1000.) { fHistMggNoCut->Fill(q4pi0.M()); fHistMggNoCutVsThetaPi0->Fill(q4pi0.Theta()*TMath::RadToDeg(),q4pi0.M()); } if (fNPi0 != 0) { for (Int_t i=0; i<2; ++i) { fHistEGamma->Fill(q4g[i].E()); fHistEGammaGenDet->Fill(q4gGen[i].E()); fHistThetaGamma->Fill(q4g[i].Theta()*TMath::RadToDeg()); fHistThetaGammaGenDet->Fill(q4gGen[i].Theta()*TMath::RadToDeg()); fHistPhiGamma->Fill(TVector2::Phi_0_2pi(q4g[i].Phi())*TMath::RadToDeg()); fHistPhiGammaGenDet->Fill(TVector2::Phi_0_2pi(q4gGen[i].Phi())*TMath::RadToDeg()); fHistEGammaVSThetaGamma->Fill(q4g[i].Theta()*TMath::RadToDeg(),q4g[i].E()); fHistPhiGammaVSThetaGamma->Fill(q4g[i].Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4g[i].Phi())*TMath::RadToDeg()); } fHistEPi0->Fill(q4pi0.E()); fHistEPi0GenDet->Fill(q4pi0Gen.E()); fHistThetaPi0->Fill(q4pi0.Theta()*TMath::RadToDeg()); fHistThetaPi0GenDet->Fill(q4pi0Gen.Theta()*TMath::RadToDeg()); fHistPhiPi0->Fill(TVector2::Phi_0_2pi(q4pi0.Phi())*TMath::RadToDeg()); fHistPhiPi0GenDet->Fill(TVector2::Phi_0_2pi(q4pi0Gen.Phi())*TMath::RadToDeg()); fHistPhiPi0VSThetaPi0->Fill(q4pi0.Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4pi0.Phi())*TMath::RadToDeg()); fHistPhiPi0ShiftVSThetaPi0->Fill(q4pi0.Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4pi0.Phi()+TMath::PiOver2())*TMath::RadToDeg()); fHistPhiPi0GenDetVSThetaPi0GenDet->Fill(q4pi0Gen.Theta()*TMath::RadToDeg(),TVector2::Phi_0_2pi(q4pi0Gen.Phi())*TMath::RadToDeg()); fHistMgg->Fill(q4pi0.M()); fHistNPi0->Fill(fNPi0); } // // fTreeAccPi0->Fill(); return NOERROR; } //___________________________________________________________ void JEventProcessor_AccPi0::ResetEvnt() { fNtrack = 0; // fNGamma = 0; fNGammaGen = 0; // fNPi0 = 0; fNPi0Gen = 0; // fNOther = 0; fNOtherGen = 0; } //------------------ // erun //------------------ jerror_t JEventProcessor_AccPi0::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_AccPi0::fini(void) { // Called before program exit after event processing is finished. return NOERROR; }