p2pi_hists/DCustomAction_p2pi_hists.cc

Go to the documentation of this file.

// $Id$
//
//    File: DCustomAction_p2pi_hists.cc
// Created: Wed Jan 21 16:53:41 EST 2015
// Creator: jrsteven (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64)
//

#include "DCustomAction_p2pi_hists.h"

void DCustomAction_p2pi_hists::Initialize(JEventLoop* locEventLoop)
{
   DApplication* dapp=dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
   JCalibration *jcalib = dapp->GetJCalibration((locEventLoop->GetJEvent()).GetRunNumber());

   // Parameters for event selection to fill histograms
   endpoint_energy = 12.;
   map<string, double> photon_endpoint_energy;
   if(jcalib->Get("/PHOTON_BEAM/endpoint_energy", photon_endpoint_energy) == false) {
      endpoint_energy = photon_endpoint_energy["PHOTON_BEAM_ENDPOINT_ENERGY"];
   }
   else {
      jout<<"No /PHOTON_BEAM/endpoint_energy for this run number: using default of 12 GeV"<<endl;
   }
   endpoint_energy_bins = (int)(20*endpoint_energy);

   cohmin_energy = 0.;
   cohedge_energy = 12.;
   map<string, double> photon_beam_param;
   if(jcalib->Get("/ANALYSIS/beam_asymmetry/coherent_energy", photon_beam_param) == false) {
      cohmin_energy = photon_beam_param["cohmin_energy"];
      cohedge_energy = photon_beam_param["cohedge_energy"];
   }
   else {
      jout<<"No /ANALYSIS/beam_asymmetry/coherent_energy for this run number: using default range of 0-12 GeV"<<endl;
   }

   dEdxCut = 2.2;
   minMMCut = 0.8;
   maxMMCut = 1.05;
   minMM2Cut = -0.006;
   maxMM2Cut = 0.004;
   missingEnergyCut = 1.0;
   minRhoMassCut = 0.6;
   maxRhoMassCut = 0.88;

   // get PID algos
        const DParticleID* locParticleID = NULL;
        locEventLoop->GetSingle(locParticleID);
        dParticleID = locParticleID;

   locEventLoop->GetSingle(dAnalysisUtilities);

   // check if a particle is missing
   auto locMissingPIDs = Get_Reaction()->Get_MissingPIDs();
   dMissingPID = (locMissingPIDs.size() == 1) ? locMissingPIDs[0] : Unknown;

   //CREATE THE HISTOGRAMS
   //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
   japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
   {
      //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
         //If another thread has already created the folder, it just changes to it. 
      CreateAndChangeTo_ActionDirectory();
      
      dEgamma = GetOrCreate_Histogram<TH1I>("Egamma", "TAGGER photon energy; E_{#gamma}", endpoint_energy_bins, 0., endpoint_energy);
      
      if(dMissingPID == Proton) {
         dMM_M2pi = GetOrCreate_Histogram<TH2I>("MM_M2pi", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
         dMM_M2pi_noEle = GetOrCreate_Histogram<TH2I>("MM_M2pi_noEle", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
         dt_M2pi_noEle = GetOrCreate_Histogram<TH2I>("t_M2pi_noEle", "|t| vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; t", 200, 0.0, 2.0, 1000, 0., 5.);
         dPiPlusPsi_t = GetOrCreate_Histogram<TH2I>("PiPlusPsi_t","#psi vs |t|; |t|; #pi^{+} #psi",1000,0.,5.,360,-180,180);
      }
      else {
         dMM2_M2pi = GetOrCreate_Histogram<TH2I>("MM2_M2pi", "MM^{2} off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM^{2}", 200, 0.0, 2.0, 200, -0.1, 0.1);
         dProton_dEdx_P = GetOrCreate_Histogram<TH2I>("Proton_dEdx_P","dE/dx vs p; p; dE/dx",200,0,2,500,0,25);
         dProton_P_Theta = GetOrCreate_Histogram<TH2I>("Proton_P_Theta","p vs #theta; #theta; p (GeV/c)",180,0,180,120,0,12);
         dDeltaE_M2pi = GetOrCreate_Histogram<TH2I>("DeltaE_M2pi", "#DeltaE vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; #DeltaE (tagger - tracks)", 200, 0.0, 2.0, 200, -5., 5.);
         dDeltaE_M2pi_ProtonTag = GetOrCreate_Histogram<TH2I>("DeltaE_M2pi_ProtonTag", "#DeltaE vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; #DeltaE (tagger - tracks)", 200, 0.0, 2.0, 200, -5., 5.);

         dDalitz_p2pi = GetOrCreate_Histogram<TH2I>("Dalitz_p2pi", "Dalitz p#pi^{+}#pi^{-}; M_{#pi^{+}#pi^{-}}^{2}; M_{p#pi^{+}}^{2}", 200, 0.0, 5.0, 200, 0., 5.);
         dMppiplus_M2pi = GetOrCreate_Histogram<TH2I>("Mppiplus_M2pi", "Psuedo-Dalitz p#pi^{+}#pi^{-}; M_{#pi^{+}#pi^{-}}; M_{p#pi^{+}}", 200, 0.0, 2.0, 200, 1.0, 3.0);
         dMppiminus_M2pi = GetOrCreate_Histogram<TH2I>("Mppiminus_M2pi", "Psuedo-Dalitz p#pi^{+}#pi^{-}; M_{#pi^{+}#pi^{-}}; M_{p#pi^{-}}", 200, 0.0, 2.0, 200, 1.0, 3.0);

         dProtonPhi_Egamma = GetOrCreate_Histogram<TH2I>("ProtonPhi_Egamma","#phi vs E_{#gamma}; E_{#gamma}; proton #phi",endpoint_energy_bins,0,endpoint_energy,360,-180,180);  
         dPiPlusPsi_Egamma = GetOrCreate_Histogram<TH2I>("PiPlusPsi_Egamma","#psi vs E_{#gamma}; E_{#gamma}; #pi^{+} #phi",endpoint_energy_bins,0,endpoint_energy,360,-180,180);
         dPiPlusPsi_t = GetOrCreate_Histogram<TH2I>("PiPlusPsi_t","#psi vs |t|; |t|; #pi^{+} #psi",1000,0.,5.,360,-180,180);

         dEgamma_M2pi = GetOrCreate_Histogram<TH2I>("Egamma_M2pi", "E_{#gamma} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; E_{#gamma}", 200, 0.0, 2.0, endpoint_energy_bins,0,endpoint_energy);     

         dBaryonM_CosTheta_Egamma1 = GetOrCreate_Histogram<TH2I>("BaryonM_CosTheta_Egamma1", "Baryon M_{p#pi^{+}} vs cos#theta: E_{#gamma} < 2.5; cos#theta; M_{p#pi^{+}}", 100, -1, 1, 100, 1.0, 2.5);
         dBaryonM_CosTheta_Egamma2 = GetOrCreate_Histogram<TH2I>("BaryonM_CosTheta_Egamma2", "Baryon M_{p#pi^{+}} vs cos#theta: 2.5 < E_{#gamma} < 3.0; cos#theta; M_{p#pi^{+}}", 100, -1, 1, 100, 1.0, 2.5);
         dBaryonM_CosTheta_Egamma3 = GetOrCreate_Histogram<TH2I>("BaryonM_CosTheta_Egamma3", "Baryon M_{p#pi^{+}} vs cos#theta: E_{#gamma} > 3.0; cos#theta; M_{p#pi^{+}}", 100, -1, 1, 100, 1.0, 2.5);
      }

      //Return to the base directory
      ChangeTo_BaseDirectory();
   }
   japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

bool DCustomAction_p2pi_hists::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{

   const DDetectorMatches* locDetectorMatches = NULL;
        locEventLoop->GetSingle(locDetectorMatches);

   // should only have one reaction step
   const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(0);

   // get beam photon energy and final state particles
   auto locBeamPhoton = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_InitialParticle() : locParticleComboStep->Get_InitialParticle_Measured();
   auto locParticles = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_FinalParticles() : locParticleComboStep->Get_FinalParticles_Measured();
        double locBeamPhotonEnergy = locBeamPhoton->energy();

   // calculate 2pi P4
   DLorentzVector locP4_2pi;
   for(size_t loc_i = 0; loc_i < 3; ++loc_i) {
      if(locParticles[loc_i] == NULL) continue; // missing particle
                if(locParticles[loc_i]->PID() == PiPlus || locParticles[loc_i]->PID() == PiMinus)
                        locP4_2pi += locParticles[loc_i]->lorentzMomentum();
        }

   // calculate missing P4
   DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo,Get_UseKinFitResultsFlag());

   // reconstructed proton variables
   double dEdx = 0.;
   DLorentzVector locProtonP4;
   if(dMissingPID != Proton) {
      const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(0));
      const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(Proton);
      dEdx = locChargedTrackHypothesis->Get_TrackTimeBased()->dEdx()*1e6; // convert to keV
      locProtonP4 = locChargedTrackHypothesis->lorentzMomentum();
   }
   else locProtonP4 = locMissingP4;

   // production kinematics
   DLorentzVector locSumInitP4; 
   DLorentzVector locProtonP4Init(0,0,0,0.938);
   locSumInitP4 += locProtonP4Init;
   locSumInitP4 += locBeamPhoton->lorentzMomentum();
   TLorentzVector locDelta = (locProtonP4 - locProtonP4Init);
   double t = locDelta.M2();

   TLorentzVector locPiPlus_P4 = locParticles[1]->lorentzMomentum();

   // boost to resonance frame for angular distributions    
   TLorentzRotation resonanceBoost( -locP4_2pi.BoostVector() );
   TLorentzVector recoil_res = resonanceBoost * locProtonP4;
   TLorentzVector p1_res = resonanceBoost * locPiPlus_P4;
   
   // normal to the production plane
   TVector3 y = (locBeamPhoton->lorentzMomentum().Vect().Unit().Cross(-locProtonP4.Vect().Unit())).Unit();
   
   // choose helicity frame: z-axis opposite recoil proton in rho rest frame
   TVector3 z = -1. * recoil_res.Vect().Unit();
   TVector3 x = y.Cross(z).Unit();
   TVector3 angles( (p1_res.Vect()).Dot(x),
          (p1_res.Vect()).Dot(y),
          (p1_res.Vect()).Dot(z) );

   double cosTheta = angles.CosTheta();
   double phi = angles.Phi();
   
   TVector3 eps(1.0, 0.0, 0.0); // beam polarization vector
   double Phi = atan2(y.Dot(eps), locBeamPhoton->lorentzMomentum().Vect().Unit().Dot(eps.Cross(y)));

   double locPsi = phi - Phi;
   if(locPsi < -1*TMath::Pi()) locPsi += 2*TMath::Pi();
        if(locPsi > TMath::Pi()) locPsi -= 2*TMath::Pi();

   //FILL HISTOGRAMS
   //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
   //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
   Lock_Action(); //ACQUIRE ROOT LOCK!!
   {
      // Fill histograms here
      dEgamma->Fill(locBeamPhotonEnergy);

      if(dMissingPID == Proton) {
         if(locBeamPhotonEnergy > cohmin_energy && locBeamPhotonEnergy < cohedge_energy) {
            dMM_M2pi->Fill(locP4_2pi.M(), locMissingP4.M());
            
            if(locParticles[1]->lorentzMomentum().Theta()*180/TMath::Pi() > 2. && locParticles[2]->lorentzMomentum().Theta()*180/TMath::Pi() > 2.) {
               dMM_M2pi_noEle->Fill(locP4_2pi.M(), locMissingP4.M());
               
               if(locMissingP4.M() > minMMCut && locMissingP4.M() < maxMMCut) {
                  dt_M2pi_noEle->Fill(locP4_2pi.M(), fabs(t));
                  
                  if(locP4_2pi.M() > minRhoMassCut && locP4_2pi.M() < maxRhoMassCut){
                     dPiPlusPsi_t->Fill(fabs(t), locPsi*180/TMath::Pi());
                  }  
               }
            }
         }           
      }
      else {
         dMM2_M2pi->Fill(locP4_2pi.M(), locMissingP4.M2());
         dDeltaE_M2pi->Fill(locP4_2pi.M(),locMissingP4.E());

         if(locMissingP4.M2() > minMM2Cut && locMissingP4.M2() < maxMM2Cut && fabs(locMissingP4.E()) < missingEnergyCut){
            dProton_dEdx_P->Fill(locProtonP4.Vect().Mag(), dEdx);
            dProton_P_Theta->Fill(locProtonP4.Vect().Theta()*180/TMath::Pi(), locProtonP4.Vect().Mag());

            if(dEdx > dEdxCut) {
               
               dEgamma_M2pi->Fill(locP4_2pi.M(), locBeamPhotonEnergy);
               dDeltaE_M2pi_ProtonTag->Fill(locP4_2pi.M(),locMissingP4.E());

               DLorentzVector locP4_ppiplus = locProtonP4;
               DLorentzVector locP4_ppiminus = locProtonP4;
               locP4_ppiplus += locParticles[1]->lorentzMomentum();
               locP4_ppiminus += locParticles[2]->lorentzMomentum();
               if(locBeamPhotonEnergy > cohmin_energy && locBeamPhotonEnergy < cohedge_energy){
                  dDalitz_p2pi->Fill(locP4_2pi.M2(), locP4_ppiplus.M2());
                  dMppiplus_M2pi->Fill(locP4_2pi.M(), locP4_ppiplus.M());
                  dMppiminus_M2pi->Fill(locP4_2pi.M(), locP4_ppiminus.M());
               }

               // some p pi+ mass distributions (hunt for Delta++...)
               if(locBeamPhotonEnergy < cohmin_energy)
                  dBaryonM_CosTheta_Egamma1->Fill(cosTheta,locP4_ppiplus.M());
               else if(locBeamPhotonEnergy < cohedge_energy)
                  dBaryonM_CosTheta_Egamma2->Fill(cosTheta,locP4_ppiplus.M());
               else 
                  dBaryonM_CosTheta_Egamma3->Fill(cosTheta,locP4_ppiplus.M());


               if(locP4_2pi.M() > minRhoMassCut && locP4_2pi.M() < maxRhoMassCut){
                  dProtonPhi_Egamma->Fill(locBeamPhotonEnergy, locProtonP4.Phi()*180/TMath::Pi());
                  dPiPlusPsi_Egamma->Fill(locBeamPhotonEnergy, locPsi*180/TMath::Pi());
                  
                  if(locBeamPhotonEnergy > cohmin_energy && locBeamPhotonEnergy < cohedge_energy){
                     dPiPlusPsi_t->Fill(fabs(t), locPsi*180/TMath::Pi());  
                  }
               }  
            }  
         }
      }
   }
   Unlock_Action(); //RELEASE ROOT LOCK!!

   return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
}