// $Id$ // // File: JEventProcessor_pipolarization.cc // Created: Tue Apr 17 08:12:00 EDT 2012 // Creator: davidl (on Darwin genmacbook-wireless.jlab.org 11.3.0 i386) // #include using namespace std; #include "JEventProcessor_pipolarization.h" using namespace jana; #include #include #include // Routine used to create our JEventProcessor #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_pipolarization()); } } // "C" //------------------ // JEventProcessor_pipolarization (Constructor) //------------------ JEventProcessor_pipolarization::JEventProcessor_pipolarization() { pthread_rwlock_init(&root_lock, NULL); } //------------------ // ~JEventProcessor_pipolarization (Destructor) //------------------ JEventProcessor_pipolarization::~JEventProcessor_pipolarization() { } //------------------ // init //------------------ jerror_t JEventProcessor_pipolarization::init(void) { //japp->RootWriteLock(); pthread_rwlock_wrlock(&root_lock); thrown_tree = new TTree("thrown_tree", "Generated values"); thrown_tree->Branch("T", &kinfo, "t/F:phi_pipi:theta:psi:Wpipi"); // Unlock ROOT read-write lock //japp->RootUnLock(); pthread_rwlock_unlock(&root_lock); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_pipolarization::brun(JEventLoop *eventLoop, int runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_pipolarization::evnt(JEventLoop *loop, int eventnumber) { vector throwns; loop->Get(throwns); TLorentzVector pip(0.0, 0.0, 0.0, 0.0), pim(0.0, 0.0, 0.0, 0.0); for(unsigned int i=0; itype == PiPlus)pip = thrown->lorentzMomentum(); if(thrown->type == PiMinus)pim = thrown->lorentzMomentum(); } if(pip.E()==0.0 || pim.E()==0.0)return NOERROR; // Lock ROOT //japp->RootWriteLock(); pthread_rwlock_wrlock(&root_lock); // Assume incident photon is 8.5GeV TLorentzVector beam(0.0, 0.0, 8.5, 8.5); TLorentzVector pip_pim = pip+pim; kinfo.t = (pip_pim - beam).M2(); kinfo.phi_pipi = pip_pim.Phi(); kinfo.Wpipi = pip_pim.M(); // Boost pi+ into pipi frame TVector3 boost_LAB_to_pipi = -pip_pim.BoostVector(); TLorentzVector pip_pipi = pip; pip_pipi.Boost(boost_LAB_to_pipi); //------ Rory's recipe ------ TVector3 k_hat = pip_pim.Vect(); k_hat.SetMag(1.0); // Rotation angle about y to get into y-z plane double theta_yz = atan2(-k_hat.X(), k_hat.Z()); k_hat.RotateY(theta_yz); // pi+ direction in rotated lab frame TVector3 pip_hat = pip.Vect(); pip_hat.RotateY(theta_yz); // Angle to rotate k_hat into z-direction double theta_x = atan2(k_hat.Y(), k_hat.Z()); k_hat.RotateX(theta_x); pip_hat.RotateX(theta_x); double psi = pip_hat.Phi(); if(psi<0.0)psi += TMath::TwoPi(); //---------------------------- // The Z-axis of interest for calculating psi is in the // direction of the pi+/pi- system, but in the pipi frame. // If we boost the pip_pim vector into this frame, the direction // would be undefined (p=0). Instead, take a vector pointing // opposite to pip_pim and boost it into the pipi frame and // then flip it. This should preserve the direction. TVector3 tmpv = -(pip_pim.Vect()); TLorentzVector tmp(tmpv, tmpv.Mag()); tmp.Boost(boost_LAB_to_pipi); TVector3 z_helicity_pipi = -tmp.Vect(); // Find theta and phi angles of z_helicity_pipi and use those // to rotate both the pi+ and polarization vectors such that // the Z-axis is now in the direction of z_helicity_pipi. This // allows us to just use the TVector3::Phi() method to extract // an angle for each of these about z_helicity_pipi. double theta_z_helicity_pipi = z_helicity_pipi.Theta(); double phi_z_helicity_pipi = z_helicity_pipi.Phi(); TVector3 pip_pipi_rotated = pip_pipi.Vect(); pip_pipi_rotated.RotateZ(-phi_z_helicity_pipi); pip_pipi_rotated.RotateY(-theta_z_helicity_pipi); kinfo.theta = pip_pipi_rotated.Theta(); kinfo.psi = psi; thrown_tree->Fill(); // Unlock ROOT //japp->RootUnLock(); pthread_rwlock_unlock(&root_lock); return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_pipolarization::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_pipolarization::fini(void) { // Called before program exit after event processing is finished. return NOERROR; } #if 0 // Dave's broken method // The Z-axis of interest for calculating psi is in the // direction of the pi+/pi- system, but in the pipi frame. // If we boost the pip_pim vector into this frame, the direction // would be undefined (p=0). Instead, take a vector pointing // opposite to pip_pim and boost it into the pipi frame and // then flip it. This should preserve the direction. TVector3 tmpv = -(pip_pim.Vect()); TLorentzVector tmp(tmpv, tmpv.Mag()); tmp.Boost(boost_LAB_to_pipi); TVector3 z_helicity_pipi = -tmp.Vect(); //-------------------------------------------------------------- // 1. Find intial photon direction in pipi frame // 2. Determine rotation matrix needed to transform from initial // photon direction in lab frame (i.e. +z) to direction in // pipi frame // 3. Use rotation matrix to rotate polarization vector to find // its direction in pipi frame // 4. Find angle between polarization vector in pipi frame and // pi+ direction in pipi frame about "z" TLorentzVector beam_pipi = beam; beam_pipi.Boost(boost_LAB_to_pipi); double theta_photon_pipi = beam_pipi.Vect().Theta(); double phi_photon_pipi = beam_pipi.Vect().Phi(); TVector3 pol(0.0, 1.0, 0.0); // in lab frame TVector3 pol_pipi = pol; pol_pipi.RotateY(theta_photon_pipi); pol_pipi.RotateZ(phi_photon_pipi); // Find theta and phi angles of z_helicity_pipi and use those // to rotate both the pi+ and polarization vectors such that // the Z-axis is now in the direction of z_helicity_pipi. This // allows us to just use the TVector3::Phi() method to extract // an angle for each of these about z_helicity_pipi. double theta_z_helicity_pipi = z_helicity_pipi.Theta(); double phi_z_helicity_pipi = z_helicity_pipi.Phi(); TVector3 pol_pipi_rotated = pol_pipi; TVector3 pip_pipi_rotated = pip_pipi.Vect(); pol_pipi_rotated.RotateZ(-phi_z_helicity_pipi); pol_pipi_rotated.RotateY(-theta_z_helicity_pipi); pip_pipi_rotated.RotateZ(-phi_z_helicity_pipi); pip_pipi_rotated.RotateY(-theta_z_helicity_pipi); double psi = pol_pipi_rotated.Phi() - pip_pipi_rotated.Phi(); if(psi<0.0)psi+=TMath::TwoPi(); #endif // Dave's broken method