#include "DSelector_2kp.h" void DSelector_2kp::Init(TTree *locTree) { // The Init() function is called when the selector needs to initialize a new tree or chain. // Typically here the branch addresses and branch pointers of the tree will be set. // Init() will be called many times when running on PROOF (once per file to be processed). //SET OUTPUT FILE NAME //can be overriden by user in PROOF dOutputFileName = "2kp.root"; //"" for none dOutputTreeFileName = ""; //"" for none //DO THIS NEXT //Because this function gets called for each TTree in the TChain, we must be careful: //We need to re-initialize the tree interface & branch wrappers, but don't want to recreate histograms bool locInitializedPriorFlag = dInitializedFlag; //save whether have been initialized previously DSelector::Init(locTree); //This must be called to initialize wrappers for each new TTree //gDirectory now points to the output file with name dOutputFileName (if any) if(locInitializedPriorFlag) return; //have already created histograms, etc. below: exit //THEN THIS Get_ComboWrappers(); dPreviousRunNumber = 0; /*********************************** EXAMPLE USER INITIALIZATION: ANALYSIS ACTIONS **********************************/ //ANALYSIS ACTIONS: //Executed in order if added to dAnalysisActions //false/true below: use measured/kinfit data //PID dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false)); //below: value: +/- N ns, Unknown: All PIDs, SYS_NULL: all timing systems dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.5, KPlus, SYS_BCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.5, KPlus, SYS_FCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.35, KPlus, SYS_TOF)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.5, KMinus, SYS_BCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.5, KMinus, SYS_FCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.35, KMinus, SYS_TOF)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.8, Proton, SYS_BCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.5, Proton, SYS_FCAL)); dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.5, Proton, SYS_TOF)); dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false, "PostPID")); //MASSES deque locPhiPIDs; locPhiPIDs.push_back(KPlus); locPhiPIDs.push_back(KMinus); deque locYStarPIDs; locYStarPIDs.push_back(Proton); locYStarPIDs.push_back(KMinus); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi_KinFit")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar_KinFit")); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, false, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0)); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, true, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0, "KinFit")); dAnalysisActions.push_back(new DHistogramAction_MissingMassSquared(dComboWrapper, false, 1000, -0.1, 0.1)); //KINFIT RESULTS dAnalysisActions.push_back(new DHistogramAction_KinFitResults(dComboWrapper)); dAnalysisActions.push_back(new DCutAction_KinFitFOM(dComboWrapper, 1.0E-4)); //MASSES dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi_PostKinFitCut")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi_PostKinFitCut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar_PostKinFitCut")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar_PostKinFitCut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, false, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0, "PostKinFitCut")); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, true, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0, "PostKinFitCut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_MissingMassSquared(dComboWrapper, false, 1000, -0.1, 0.1, "PostKinFitCut")); //BEAM ENERGY dAnalysisActions.push_back(new DHistogramAction_BeamEnergy(dComboWrapper, false)); dAnalysisActions.push_back(new DCutAction_BeamEnergy(dComboWrapper, false, 8.4, 9.05)); //MASSES dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi_PostBeamECut")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locPhiPIDs, 1000, 0.9, 1.4, "Phi_PostBeamECut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar_PostBeamECut")); dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, true, 0, locYStarPIDs, 1800, 1.3, 3.1, "YStar_PostBeamECut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, false, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0, "PostBeamECut")); dAnalysisActions.push_back(new DHistogramAction_Dalitz(dComboWrapper, true, 0, locPhiPIDs, locYStarPIDs, 760, 0.8, 16.0, 1440, 1.6, 16.0, "PostBeamECut_KinFit")); dAnalysisActions.push_back(new DHistogramAction_MissingMassSquared(dComboWrapper, false, 1000, -0.1, 0.1, "PostBeamECut")); //KINEMATICS dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false)); dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, true, "KinFit")); //INITIALIZE ACTIONS //If you create any actions that you want to run manually (i.e. don't add to dAnalysisActions), be sure to initialize them here as well Initialize_Actions(); /******************************** EXAMPLE USER INITIALIZATION: STAND-ALONE HISTOGRAMS *******************************/ //BEST COMBO HISTOGRAMS gDirectory->mkdir("BestCombo", "BestCombo")->cd(); { dHist_YStarMass_BestCombo = new TH1I("YStarMass", ";pK^{#minus} Invariant Mass (GeV/c^{2})", 800, 1.3, 2.1); dHist_MandelstamT = new TH1I("MandelstamT", ";t (GeV^{2})", 200, -5.0, 0.0); dHist_YStarMassVsT = new TH2I("YStarMassVsT", ";t (GeV^{2});pK^{#minus} Invariant Mass (GeV/c^{2})", 200, -5.0, 0.0, 800, 1.3, 2.1); //BEST COMBO: POLARIZATION dHist_ProdPlanePhi_PARA = new TH1I("ProdPlanePhi_PARA", ";Production Plane #phi#circ", 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Sideband_PARA = new TH2I("ProdPlanePhiVsT_Sideband_PARA", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Signal_PARA = new TH2I("ProdPlanePhiVsT_Signal_PARA", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); dHist_ProdPlanePhi_PERP = new TH1I("ProdPlanePhi_PERP", ";Production Plane #phi#circ", 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Sideband_PERP = new TH2I("ProdPlanePhiVsT_Sideband_PERP", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Signal_PERP = new TH2I("ProdPlanePhiVsT_Signal_PERP", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); dHist_ProdPlanePhi_AMO = new TH1I("ProdPlanePhi_AMO", ";Production Plane #phi#circ", 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Sideband_AMO = new TH2I("ProdPlanePhiVsT_Sideband_AMO", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); dHist_ProdPlanePhiVsT_Signal_AMO = new TH2I("ProdPlanePhiVsT_Signal_AMO", ";t (GeV^{2});Production Plane #phi#circ", 200, -5.0, 0.0, 360, -180.0, 180.0); } gDirectory->cd(".."); /***************************************** ADVANCED: CHOOSE BRANCHES TO READ ****************************************/ //TO SAVE PROCESSING TIME //If you know you don't need all of the branches/data, but just a subset of it, you can speed things up //By default, for each event, the data is retrieved for all branches //If you know you only need data for some branches, you can skip grabbing data from the branches you don't need //Do this by doing something similar to the commented code below //dTreeInterface->Clear_GetEntryBranches(); //now get none //dTreeInterface->Register_GetEntryBranch("Proton__P4"); //manually set the branches you want } Bool_t DSelector_2kp::Process(Long64_t locEntry) { // The Process() function is called for each entry in the tree. The entry argument // specifies which entry in the currently loaded tree is to be processed. // // This function should contain the "body" of the analysis. It can contain // simple or elaborate selection criteria, run algorithms on the data // of the event and typically fill histograms. // // The processing can be stopped by calling Abort(). // Use fStatus to set the return value of TTree::Process(). // The return value is currently not used. //CALL THIS FIRST DSelector::Process(locEntry); //Gets the data from the tree for the entry //cout << "RUN " << Get_RunNumber() << ", EVENT " << Get_EventNumber() << endl; /******************************************** GET POLARIZATION ORIENTATION ******************************************/ //Only if the run number changes //RCDB environment must be setup in order for this to work! (Will return false otherwise) UInt_t locRunNumber = Get_RunNumber(); if(locRunNumber != dPreviousRunNumber) { dIsPolarizedFlag = dAnalysisUtilities.Get_IsPolarizedBeam(locRunNumber, dIsPARAFlag); dPreviousRunNumber = locRunNumber; } /********************************************* SETUP UNIQUENESS TRACKING ********************************************/ //ANALYSIS ACTIONS: Reset uniqueness tracking for each action //For any actions that you are executing manually, be sure to call Reset_NewEvent() on them here Reset_Actions_NewEvent(); //PREVENT-DOUBLE COUNTING WHEN HISTOGRAMMING //Sometimes, some content is the exact same between one combo and the next //e.g. maybe two combos have different beam particles, but the same data for the final-state //When histogramming, you don't want to double-count when this happens: artificially inflates your signal (or background) //So, for each quantity you histogram, keep track of what particles you used (for a given combo) //Then for each combo, just compare to what you used before, and make sure it's unique //EXAMPLE 1: Particle-specific info: set locUsedSoFar_BeamEnergy; //Int_t: Unique ID for beam particles. set: easy to use, fast to search //EXAMPLE 2: Combo-specific info: //In general: Could have multiple particles with the same PID: Use a set of Int_t's //In general: Multiple PIDs, so multiple sets: Contain within a map //Multiple combos: Contain maps within a set (easier, faster to search) set > > locUsedSoFar_MissingMass; //INSERT USER ANALYSIS UNIQUENESS TRACKING HERE /************************************************* LOOP OVER COMBOS *************************************************/ //Loop over combos UInt_t locNumSurvivingCombos = 0; set locSurvivingBeamIDs; double locBestMissingMassSquared = 9.9E9; double locT_BestCombo = 0.0, locYStarMass_BestCombo = 0.0, locProdPlanePhi_BestCombo = 0.0; bool locIsYStarSignal_BestCombo = false; for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) { //Set branch array indices for combo and all combo particles dComboWrapper->Set_ComboIndex(loc_i); // Is used to indicate when combos have been cut if(dComboWrapper->Get_IsComboCut()) // Is false when tree originally created continue; // Combo has been cut previously /********************************************** GET PARTICLE INDICES *********************************************/ //Used for tracking uniqueness when filling histograms, and for determining unused particles //Step 0 Int_t locBeamID = dComboBeamWrapper->Get_BeamID(); Int_t locKPlusTrackID = dKPlusWrapper->Get_TrackID(); Int_t locKMinusTrackID = dKMinusWrapper->Get_TrackID(); Int_t locProtonTrackID = dProtonWrapper->Get_TrackID(); /*********************************************** GET FOUR-MOMENTUM **********************************************/ // Get P4's: //is kinfit if kinfit performed, else is measured //dTargetP4 is target p4 //Step 0 TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4(); TLorentzVector locKPlusP4 = dKPlusWrapper->Get_P4(); TLorentzVector locKMinusP4 = dKMinusWrapper->Get_P4(); TLorentzVector locProtonP4 = dProtonWrapper->Get_P4(); // Get Measured P4's: //Step 0 TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured(); TLorentzVector locKPlusP4_Measured = dKPlusWrapper->Get_P4_Measured(); TLorentzVector locKMinusP4_Measured = dKMinusWrapper->Get_P4_Measured(); TLorentzVector locProtonP4_Measured = dProtonWrapper->Get_P4_Measured(); /********************************************* COMBINE FOUR-MOMENTUM ********************************************/ // DO YOUR STUFF HERE // Combine 4-vectors TLorentzVector locMissingP4_Measured = locBeamP4_Measured + dTargetP4; locMissingP4_Measured -= locKPlusP4_Measured + locKMinusP4_Measured + locProtonP4_Measured; TLorentzVector locYStarP4_Measured = locProtonP4_Measured + locKMinusP4_Measured; TLorentzVector locYStarP4_KinFit = locProtonP4 + locKMinusP4; /******************************************** EXECUTE ANALYSIS ACTIONS *******************************************/ // Loop through the analysis actions, executing them in order for the active particle combo if(!Execute_Actions()) //if the active combo fails a cut, IsComboCutFlag automatically set continue; //if you manually execute any actions, and it fails a cut, be sure to call: //dComboWrapper->Set_IsComboCut(true); /************************************************* MISCELLANEOUS *************************************************/ double locYStarMass = locYStarP4_Measured.M(); double locYStarMass_KinFit = locYStarP4_KinFit.M(); ++locNumSurvivingCombos; locSurvivingBeamIDs.insert(locBeamID); bool locIsYStarSignalFlag = ((locYStarMass > 1.495) && (locYStarMass <= 1.545)); double locProdPlanePhi = dAnalysisUtilities.Calc_ProdPlanePhi_Pseudoscalar(locBeamP4.E(), dTargetPID, locKPlusP4); double locT = (locBeamP4 - locKPlusP4).Mag2(); /**************************************************** BEST COMBO ****************************************************/ // Signal selection and background subtraction are performed on the YStar mass peak. // However, a given YStar (proton, pi-, photon combo) can be matched to multiple beam photons (or kaon candidates) // This is OK in a 1D plot (check for duplicates), but we are quoting results in 2D: In bins of Mandelstam-t // We cannot put a signal YStar mass value into two separate t-bins: We would count that signal event in both t-bins. // Therefore, we must choose a t-bin. This is, in effect, a choice of "Best" combo. // We must choose this "best" combo in a way that doesn't bias the background subtraction. // Since the background subtraction is done on the YStar mass peak, we must choose something else. // Ideally we would pick the combo with the best kinematic-fit confidence level (or BDT response when it's ready). // However, at the time of this writing, the tracking errors are bad and this quantity cannot be trusted. // Therefore, we choose the combo that has a missing mass squared closest to zero: It's calculation involves all particles. // As far as the code is concerned, it is easiest to keep track of the relevant quantities for the "Best" combo // and update them as new "Best" combos are found. Then histogram the "Best" results at the end. double locMissingMassSquared = locMissingP4_Measured.M2(); if(fabs(locMissingMassSquared) > fabs(locBestMissingMassSquared)) continue; locBestMissingMassSquared = locMissingMassSquared; locT_BestCombo = locT; locYStarMass_BestCombo = locYStarMass; locProdPlanePhi_BestCombo = locProdPlanePhi; locIsYStarSignal_BestCombo = locIsYStarSignalFlag; } // end of combo loop //FILL HISTOGRAMS: Num combos / events surviving actions Fill_NumCombosSurvivedHists(); //FILL ASYMMETRY INFO FOR BEST COMBOS if(locNumSurvivingCombos > 0) { dHist_YStarMass_BestCombo->Fill(locYStarMass_BestCombo); dHist_YStarMassVsT->Fill(locT_BestCombo, locYStarMass_BestCombo); dHist_MandelstamT->Fill(locT_BestCombo); if((locYStarMass_BestCombo > 1.48) && (locYStarMass_BestCombo <= 1.58)) { if(!dIsPolarizedFlag) { dHist_ProdPlanePhiVsT_Signal_AMO->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); dHist_ProdPlanePhi_AMO->Fill(locProdPlanePhi_BestCombo); } else if(dIsPARAFlag) { dHist_ProdPlanePhiVsT_Signal_PARA->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); dHist_ProdPlanePhi_PARA->Fill(locProdPlanePhi_BestCombo); } else //PERP { dHist_ProdPlanePhiVsT_Signal_PERP->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); dHist_ProdPlanePhi_PERP->Fill(locProdPlanePhi_BestCombo); } } else if((locYStarMass_BestCombo > 1.6) && (locYStarMass_BestCombo <= 1.68)) //sideband { if(!dIsPolarizedFlag) dHist_ProdPlanePhiVsT_Sideband_AMO->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); else if(dIsPARAFlag) dHist_ProdPlanePhiVsT_Sideband_PARA->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); else //PERP dHist_ProdPlanePhiVsT_Sideband_PERP->Fill(locT_BestCombo, locProdPlanePhi_BestCombo); } } /************************************ EXAMPLE: FILL CLONE OF TTREE HERE WITH CUTS APPLIED ************************************/ /* Bool_t locIsEventCut = true; for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) { //Set branch array indices for combo and all combo particles dComboWrapper->Set_ComboIndex(loc_i); // Is used to indicate when combos have been cut if(dComboWrapper->Get_IsComboCut()) continue; locIsEventCut = false; // At least one combo succeeded break; } if(!locIsEventCut && dOutputTreeFileName != "") FillOutputTree(); */ return kTRUE; } void DSelector_2kp::Finalize(void) { //Save anything to output here that you do not want to be in the default DSelector output ROOT file. //Otherwise, don't do anything else (especially if you are using PROOF). //If you are using PROOF, this function is called on each thread, //so anything you do will not have the combined information from the various threads. //Besides, it is best-practice to do post-processing (e.g. fitting) separately, in case there is a problem. //DO YOUR STUFF HERE //CALL THIS LAST DSelector::Finalize(); //Saves results to the output file }