#include "DSelector_omegaeta.h" void DSelector_omegaeta::Init(TTree *locTree) { // USERS: IN THIS FUNCTION, ONLY MODIFY SECTIONS WITH A "USER" OR "EXAMPLE" LABEL. LEAVE THE REST ALONE. // 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). //USERS: SET OUTPUT FILE NAME //can be overriden by user in PROOF dOutputFileName = "omegaeta.root"; //"" for none dOutputTreeFileName = ""; //"" for none dFlatTreeFileName = ""; //output flat tree (one combo per tree entry), "" for none dFlatTreeName = ""; //if blank, default name will be chosen //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 Get_ComboWrappers(); dPreviousRunNumber = 0; /*********************************** EXAMPLE USER INITIALIZATION: ANALYSIS ACTIONS **********************************/ // EXAMPLE: Create deque for histogramming particle masses: // // For histogramming the phi mass in phi -> K+ K- // // Be sure to change this and dAnalyzeCutActions to match reaction std::deque MyPhi,MyEta; MyPhi.push_back(PiPlus); MyPhi.push_back(PiMinus); MyPhi.push_back(Pi0); MyEta.push_back(Gamma); MyEta.push_back(Gamma); //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)); //MASSES dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, Eta, 1000, 0.35, 0.75, "Eta")); //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,0.01)); dAnalysisActions.push_back(new DCutAction_Energy_UnusedShowers(dComboWrapper, 0.01)); dAnalysisActions.push_back(new DCutAction_NumUnusedTracks(dComboWrapper,0)); //CUT MISSING MASS //dAnalysisActions.push_back(new DCutAction_MissingMassSquared(dComboWrapper, false, -0.03, 0.02)); //BEAM ENERGY dAnalysisActions.push_back(new DHistogramAction_BeamEnergy(dComboWrapper, false)); //dAnalysisActions.push_back(new DCutAction_BeamEnergy(dComboWrapper, false, 8.4, 9.05)); //KINEMATICS dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false)); // ANALYZE CUT ACTIONS // // Change MyPhi to match reaction dAnalyzeCutActions = new DHistogramAction_AnalyzeCutActions( dAnalysisActions, dComboWrapper, true, 0, MyPhi, 1000, 0., 2.4, "CutActionEffect" ); //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(); dAnalyzeCutActions->Initialize(); // manual action, must call Initialize() /******************************** EXAMPLE USER INITIALIZATION: STAND-ALONE HISTOGRAMS *******************************/ //EXAMPLE MANUAL HISTOGRAMS: dHist_MissingMassSquared = new TH1F("MissingMassSquared", ";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.06, 0.06); dHist_PipPimPi0Mass_vs_2GammaMass = new TH2F("PipPimPi0Mass_vs_2GammaMass", ";M(#pi^{0}#pi^{-}#pi^{+}) [GeV] vs M(2#gamma) [GeV]", 400,0.35,0.75,1000, 0.4, 1.4); dHist_OmegaEtaMass = new TH1F("OmegaEtaMass", ";M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_MandelstamT = new TH2F("MandelstamT", "t vs M(#eta#omega);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25,250,0,2.5); dHist_H2000 = new TH1F("H2000", "H(2000);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2020 = new TH1F("H2020", "H(2020);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0020 = new TH1F("H0020", "H(0020);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0022 = new TH1F("H0022", "H(0022);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2222 = new TH1F("H2222", "H(2222);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2200 = new TH1F("H2200", "H(2200);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2220 = new TH1F("H2220", "H(2220);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2120 = new TH1F("H2120", "H(2120);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2100 = new TH1F("H2100", "H(2100);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0021 = new TH1F("H0021", "H(0021);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2221 = new TH1F("H2221", "H(2221);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2212 = new TH1F("H2212", "H(2212);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2110 = new TH1F("H2110", "H(2110);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2021 = new TH1F("H2021", "H(2021);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2121 = new TH1F("H2121", "H(2121);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2122 = new TH1F("H2122", "H(2122);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0040 = new TH1F("H0040", "H(0040);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2040 = new TH1F("H2040", "H(2040);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2140 = new TH1F("H2140", "H(2140);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2240 = new TH1F("H2240", "H(2240);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0041 = new TH1F("H0041", "H(0041);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2041 = new TH1F("H2041", "H(2041);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2141 = new TH1F("H2141", "H(2141);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2241 = new TH1F("H2241", "H(2241);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0042 = new TH1F("H0042", "H(0042);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2042 = new TH1F("H2042", "H(2042);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2142 = new TH1F("H2142", "H(2142);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2242 = new TH1F("H2242", "H(2242);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0043 = new TH1F("H0043", "H(0043);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2043 = new TH1F("H2043", "H(2043);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2143 = new TH1F("H2143", "H(2143);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2243 = new TH1F("H2243", "H(2243);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H0044 = new TH1F("H0044", "H(0044);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2044 = new TH1F("H2044", "H(2044);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2144 = new TH1F("H2144", "H(2144);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_H2244 = new TH1F("H2244", "H(2244);M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_OmegaEtaMass_SB = new TH1F("OmegaEtaMass_SB", ";M(#pi^{0}#pi^{-}#pi^{+}#eta) [GeV]", 50,1.25,3.25); dHist_BeamEnergy = new TH1F("BeamEnergy", ";Beam Energy (GeV)", 600, 0.0, 12.0); dHist_CosThetaGF = new TH2F("CosThetaGF", "cos(#theta_{GF}) vs M(#omega#eta) [GeV]",50, 1.25, 3.25,50,-1,1); dHist_LambdaOmega = new TH1F("LambdaOmega",";#lambda_{#omega}",100,0,1); dHist_OmegaAngDist=new TH2F("OmegaAngDist","#omega angular distr.",100,-M_PI,M_PI,100,-1,1); dHist_OffCombos=new TH2F("OffCombos","M(pair1) vs M(pair2)",1000,0,1,1000,0,1); dHist_OffCombosCut=new TH2F("OffCombosCut","M(pair1) vs M(pair2)",1000,0,1,1000,0,1); /************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - MAIN TREE *************************/ //EXAMPLE MAIN TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)): //The type for the branch must be included in the brackets //1st function argument is the name of the branch //2nd function argument is the name of the branch that contains the size of the array (for fundamentals only) /* dTreeInterface->Create_Branch_Fundamental("my_int"); //fundamental = char, int, float, double, etc. dTreeInterface->Create_Branch_FundamentalArray("my_int_array", "my_int"); dTreeInterface->Create_Branch_FundamentalArray("my_combo_array", "NumCombos"); dTreeInterface->Create_Branch_NoSplitTObject("my_p4"); dTreeInterface->Create_Branch_ClonesArray("my_p4_array"); */ /************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - FLAT TREE *************************/ //EXAMPLE FLAT TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)): //The type for the branch must be included in the brackets //1st function argument is the name of the branch //2nd function argument is the name of the branch that contains the size of the array (for fundamentals only) /* dFlatTreeInterface->Create_Branch_Fundamental("flat_my_int"); //fundamental = char, int, float, double, etc. dFlatTreeInterface->Create_Branch_FundamentalArray("flat_my_int_array", "flat_my_int"); dFlatTreeInterface->Create_Branch_NoSplitTObject("flat_my_p4"); dFlatTreeInterface->Create_Branch_ClonesArray("flat_my_p4_array"); */ /************************************* ADVANCED EXAMPLE: 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_omegaeta::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; //TLorentzVector locProductionX4 = Get_X4_Production(); /******************************************** 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(); dAnalyzeCutActions->Reset_NewEvent(); // manual action, must call Reset_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 /**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/ /* Int_t locMyInt = 7; dTreeInterface->Fill_Fundamental("my_int", locMyInt); TLorentzVector locMyP4(4.0, 3.0, 2.0, 1.0); dTreeInterface->Fill_TObject("my_p4", locMyP4); for(int loc_i = 0; loc_i < locMyInt; ++loc_i) dTreeInterface->Fill_Fundamental("my_int_array", 3*loc_i, loc_i); //2nd argument = value, 3rd = array index */ /************************************************* LOOP OVER COMBOS *************************************************/ //Loop over combos 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 locPiPlusTrackID = dPiPlusWrapper->Get_TrackID(); Int_t locPiMinusTrackID = dPiMinusWrapper->Get_TrackID(); Int_t locProtonTrackID = dProtonWrapper->Get_TrackID(); //Step 1 Int_t locPhoton1NeutralID = dPhoton1Wrapper->Get_NeutralID(); Int_t locPhoton2NeutralID = dPhoton2Wrapper->Get_NeutralID(); //Step 2 Int_t locPhoton3NeutralID = dPhoton3Wrapper->Get_NeutralID(); Int_t locPhoton4NeutralID = dPhoton4Wrapper->Get_NeutralID(); /*********************************************** 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 locPiPlusP4 = dPiPlusWrapper->Get_P4(); TLorentzVector locPiMinusP4 = dPiMinusWrapper->Get_P4(); TLorentzVector locProtonP4 = dProtonWrapper->Get_P4(); //Step 1 TLorentzVector locDecayingPi0P4 = dDecayingPi0Wrapper->Get_P4(); TLorentzVector locPhoton1P4 = dPhoton1Wrapper->Get_P4(); TLorentzVector locPhoton2P4 = dPhoton2Wrapper->Get_P4(); //Step 2 TLorentzVector locPhoton3P4 = dPhoton3Wrapper->Get_P4(); TLorentzVector locPhoton4P4 = dPhoton4Wrapper->Get_P4(); // Get Measured P4's: //Step 0 TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured(); TLorentzVector locPiPlusP4_Measured = dPiPlusWrapper->Get_P4_Measured(); TLorentzVector locPiMinusP4_Measured = dPiMinusWrapper->Get_P4_Measured(); TLorentzVector locProtonP4_Measured = dProtonWrapper->Get_P4_Measured(); //Step 1 TLorentzVector locPhoton1P4_Measured = dPhoton1Wrapper->Get_P4_Measured(); TLorentzVector locPhoton2P4_Measured = dPhoton2Wrapper->Get_P4_Measured(); //Step 2 TLorentzVector locPhoton3P4_Measured = dPhoton3Wrapper->Get_P4_Measured(); TLorentzVector locPhoton4P4_Measured = dPhoton4Wrapper->Get_P4_Measured(); /********************************************* COMBINE FOUR-MOMENTUM ********************************************/ // DO YOUR STUFF HERE // Combine 4-vectors TLorentzVector locMissingP4_Measured = locBeamP4_Measured + dTargetP4; locMissingP4_Measured -= locPiPlusP4_Measured + locPiMinusP4_Measured + locProtonP4_Measured + locPhoton1P4_Measured + locPhoton2P4_Measured + locPhoton3P4_Measured + locPhoton4P4_Measured; /******************************************** EXECUTE ANALYSIS ACTIONS *******************************************/ // Loop through the analysis actions, executing them in order for the active particle combo dAnalyzeCutActions->Perform_Action(); // Must be executed before Execute_Actions() 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); /**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/ /* TLorentzVector locMyComboP4(8.0, 7.0, 6.0, 5.0); //for arrays below: 2nd argument is value, 3rd is array index //NOTE: By filling here, AFTER the cuts above, some indices won't be updated (and will be whatever they were from the last event) //So, when you draw the branch, be sure to cut on "IsComboCut" to avoid these. dTreeInterface->Fill_Fundamental("my_combo_array", -2*loc_i, loc_i); dTreeInterface->Fill_TObject("my_p4_array", locMyComboP4, loc_i); */ /**************************************** EXAMPLE: HISTOGRAM BEAM ENERGY *****************************************/ //Histogram beam energy (if haven't already) if(locUsedSoFar_BeamEnergy.find(locBeamID) == locUsedSoFar_BeamEnergy.end()) { dHist_BeamEnergy->Fill(locBeamP4.E()); locUsedSoFar_BeamEnergy.insert(locBeamID); } /************************************ EXAMPLE: HISTOGRAM MISSING MASS SQUARED ************************************/ //Missing Mass Squared double locMissingMassSquared = locMissingP4_Measured.M2(); //Uniqueness tracking: Build the map of particles used for the missing mass //For beam: Don't want to group with final-state photons. Instead use "Unknown" PID (not ideal, but it's easy). map > locUsedThisCombo_MissingMass; locUsedThisCombo_MissingMass[Unknown].insert(locBeamID); //beam locUsedThisCombo_MissingMass[PiPlus].insert(locPiPlusTrackID); locUsedThisCombo_MissingMass[PiMinus].insert(locPiMinusTrackID); locUsedThisCombo_MissingMass[Proton].insert(locProtonTrackID); locUsedThisCombo_MissingMass[Gamma].insert(locPhoton1NeutralID); locUsedThisCombo_MissingMass[Gamma].insert(locPhoton2NeutralID); locUsedThisCombo_MissingMass[Gamma].insert(locPhoton3NeutralID); locUsedThisCombo_MissingMass[Gamma].insert(locPhoton4NeutralID); //compare to what's been used so far if(locUsedSoFar_MissingMass.find(locUsedThisCombo_MissingMass) == locUsedSoFar_MissingMass.end()) { //unique missing mass combo: histogram it, and register this combo of particles dHist_MissingMassSquared->Fill(locMissingMassSquared); locUsedSoFar_MissingMass.insert(locUsedThisCombo_MissingMass); TLorentzVector locVertex = dComboBeamWrapper->Get_X4_Measured(); float locRFTime = dComboWrapper->Get_RFTime(); float DT_RF = locVertex.T() - (locRFTime + (locVertex.Z() - dTargetCenter.Z())/29.9792458); // Now an event weight can be assigned: // if DT_RF = +/- 2ns within zero the beam photon is in time // within +-4, +-8, +-12, ... the beam photon is out of time double Weight = 1.; if (TMath::Abs(DT_RF)>2.){ Weight = -0.125; } double HalfWeight=0.25*Weight; TLorentzVector g13=locPhoton1P4_Measured+locPhoton3P4_Measured; TLorentzVector g24=locPhoton2P4_Measured+locPhoton4P4_Measured; dHist_OffCombos->Fill(g13.M(),g24.M(),Weight); TLorentzVector g14=locPhoton1P4_Measured+locPhoton4P4_Measured; TLorentzVector g23=locPhoton2P4_Measured+locPhoton3P4_Measured; dHist_OffCombos->Fill(g14.M(),g23.M(),Weight); TLorentzVector myPipPimPi0=locPiPlusP4+locPiMinusP4 +locDecayingPi0P4; TLorentzVector my2Gamma= locPhoton3P4+locPhoton4P4; double my2GammaMass=my2Gamma.M(); double myPipPimPi0Mass=myPipPimPi0.M(); double MandelstamT=(dTargetP4-locProtonP4).M2(); if (my2GammaMass>0.5 && my2GammaMass<0.6){ dHist_OffCombosCut->Fill(g13.M(),g24.M(),Weight); dHist_OffCombosCut->Fill(g14.M(),g23.M(),Weight); } if (fabs(g13.M()-ParticleMass(Pi0))>0.042 &&fabs(g14.M()-ParticleMass(Pi0))>0.042 &&fabs(g23.M()-ParticleMass(Pi0))>0.042 &&fabs(g24.M()-ParticleMass(Pi0))>0.042 ){ dHist_PipPimPi0Mass_vs_2GammaMass->Fill(my2Gamma.M(),myPipPimPi0.M(),Weight); if (my2GammaMass>0.5 && my2GammaMass<0.6 && myPipPimPi0Mass>0.74 && myPipPimPi0Mass<0.82){ TLorentzVector myPipPimPi0Eta=myPipPimPi0+my2Gamma; TVector3 beta=myPipPimPi0Eta.Vect(); beta*=-1./myPipPimPi0Eta.E(); TLorentzVector myPipPimPi0_cm(myPipPimPi0); TLorentzVector my2Gamma_cm(my2Gamma); TLorentzVector myBeam_cm(locBeamP4); my2Gamma_cm.Boost(beta); myPipPimPi0_cm.Boost(beta); myBeam_cm.Boost(beta); TVector3 pipPimPi0dir=myPipPimPi0_cm.Vect(); pipPimPi0dir.SetMag(1.); TVector3 beamDir=myBeam_cm.Vect(); beamDir.SetMag(1.); TVector3 yHatGF=beamDir.Cross(pipPimPi0dir); TVector3 xHatGF=yHatGF.Cross(beamDir); pipPimPi0dir.RotateUz(beamDir); double cosThetaGF=pipPimPi0dir.z(); double phiGF=atan2(pipPimPi0dir.y(),pipPimPi0dir.x()); dHist_CosThetaGF->Fill(myPipPimPi0Eta.M(),cosThetaGF,Weight); TVector3 omega_beta=myPipPimPi0.Vect(); omega_beta*=-1./myPipPimPi0.E(); TLorentzVector myPipCM(locPiPlusP4); TLorentzVector myPimCM(locPiMinusP4); myPipCM.Boost(omega_beta); myPimCM.Boost(omega_beta); TVector3 myPipv3=myPipCM.Vect(); TVector3 myPimv3=myPimCM.Vect(); TVector3 cross=myPimv3.Cross(myPipv3); double lambda=cross.Mag2()/(0.75*pow(myPipPimPi0.M2()/9.-0.14*.14,2)); cross.SetMag(1.); omega_beta.SetMag(1.); cross.RotateUz(omega_beta); double cosThetaH=cross.z(); double sinThetaGF=sin(acos(cosThetaGF)); double sinThetaGFsq=sinThetaGF*sinThetaGF; double cosThetaHsq=cosThetaH*cosThetaH; double cosThetaGFsq=cosThetaGF*cosThetaGF; double sinThetaHsq=1.-cosThetaHsq; double sinThetaH=sin(acos(cosThetaH)); double phiH=cross.Phi(); dHist_OmegaAngDist->Fill(phiH,cosThetaH,Weight); double Weight2000=1.5*cosThetaHsq-0.5; double Weight0020=1.5*cosThetaGFsq-0.5; double Weight0021=-0.5*(1+cosThetaGF)*sinThetaGF*cos(phiGF); double Weight2020=Weight0020*Weight2000; double Weight2021=Weight0021*Weight2000; double Weight2200=cos(2.*phiH)*sqrt(6.)/4.*sinThetaHsq; double Weight2220=Weight2200*Weight0020; double Weight2221=Weight2200*Weight0021; double Weight2100=cos(phiH)*(-0.5)*(1+cosThetaH)*sinThetaH; double Weight2121=Weight2100*Weight0021; double Weight2120=Weight2100*Weight0020; double Weight2110=Weight2100*cosThetaGF; double Weight2111=Weight2100*(-1./sqrt(2.))*sinThetaGF*cos(phiGF); double Weight0022=cos(2.*phiGF)*sqrt(6.)/4.*sinThetaGFsq; double Weight2222=Weight2200*Weight0022; double Weight2122=Weight2100*cos(2*phiGF) *(-0.5)*(1.+cosThetaGF)*sinThetaGF; double Weight2212=Weight2212=Weight2200*cos(2.*cos(phiGF)) *0.5*(1+cosThetaGF)*sinThetaGF; dHist_H2000->Fill(myPipPimPi0Eta.M(),Weight*Weight2000); dHist_H0020->Fill(myPipPimPi0Eta.M(),Weight*Weight0020); dHist_H0021->Fill(myPipPimPi0Eta.M(),Weight*Weight0021); dHist_H2020->Fill(myPipPimPi0Eta.M(),Weight*Weight2020); dHist_H2021->Fill(myPipPimPi0Eta.M(),Weight*Weight2021); dHist_H2200->Fill(myPipPimPi0Eta.M(),Weight*Weight2200); dHist_H2220->Fill(myPipPimPi0Eta.M(),Weight*Weight2220); dHist_H2221->Fill(myPipPimPi0Eta.M(),Weight*Weight2221); dHist_H2100->Fill(myPipPimPi0Eta.M(),Weight*Weight2100); dHist_H2120->Fill(myPipPimPi0Eta.M(),Weight*Weight2120); dHist_H2121->Fill(myPipPimPi0Eta.M(),Weight*Weight2121); dHist_H2122->Fill(myPipPimPi0Eta.M(),Weight*Weight2122); dHist_H2110->Fill(myPipPimPi0Eta.M(),Weight*Weight2110); dHist_H2111->Fill(myPipPimPi0Eta.M(),Weight*Weight2111); dHist_H0022->Fill(myPipPimPi0Eta.M(),Weight*Weight0022); dHist_H2222->Fill(myPipPimPi0Eta.M(),Weight*Weight2222); dHist_H2212->Fill(myPipPimPi0Eta.M(),Weight*Weight2212); dHist_LambdaOmega->Fill(lambda,Weight); dHist_OmegaEtaMass->Fill(myPipPimPi0Eta.M(),Weight); dHist_MandelstamT->Fill(myPipPimPi0Eta.M(),-MandelstamT,Weight); double Weight0040=0.125*(35.*cosThetaGFsq*cosThetaGFsq -30.*cosThetaGFsq+3.); double Weight2040=Weight2000*Weight0040; double Weight0041=-cos(phiGF)*sinThetaGF*(35*cosThetaGF*cosThetaGFsq-15*cosThetaGF)/(4.*sqrt(5.)); double Weight2041=Weight2000*Weight0041; double Weight0042=-cos(2*phiGF)*sinThetaGFsq*(35.*cosThetaGFsq-5)/(4.*sqrt(5.)); double Weight2042=Weight2000*Weight0042; double Weight0043=-cos(3.*phiGF)*sqrt(35.)/4.*sinThetaGFsq*sinThetaGF*cosThetaGF; double Weight2043=Weight2000*Weight0043; double Weight0044=cos(4.*phiGF)*35.*sinThetaGFsq*sinThetaGFsq/(8.*sqrt(70.)); double Weight2044=Weight2000*Weight0044; double Weight2240=Weight2200*Weight0042; double Weight2241=Weight2200*sqrt(2)/8.*cos(phiGF) *sinThetaGF*(1.+cosThetaGF)*(14*cosThetaGFsq-7*cosThetaGF-1); double Weight2242=Weight2200*cos(2.*phiGF)*0.25 *(1-5*cosThetaGF-6*cosThetaGFsq+7*cosThetaGF*cosThetaGFsq +7.*cosThetaGFsq*cosThetaGFsq); // I am not sure about the signs for the following. There seems // to be a discrepancy between Mathematic's WignerD and PDG. double Weight2243=Weight2200*cos(3.*phiGF)*sqrt(14.)/8. *sinThetaGF*(-1+3.*cosThetaGFsq+2*cosThetaGF*cosThetaGFsq); double Weight2244=Weight2200*cos(4.*phiGF)*sqrt(7.)/8. *sinThetaGFsq*pow(1.+cosThetaGF,2); double Weight2140=Weight2100*sinThetaGF*(35*cosThetaGF*cosThetaGFsq-15*cosThetaGF)/(4.*sqrt(5.)); double Weight2141=Weight2100*cos(phiGF)*0.125*(1.+cosThetaGF) *(3-6*cosThetaGF-21*cosThetaGFsq+28*cosThetaGF*cosThetaGFsq); double Weight2142=-Weight2100*cos(2.*phiGF)*sqrt(2)/8. *sinThetaGF*(1.+cosThetaGF)*(14*cosThetaGFsq-7*cosThetaGF-1); double Weight2143= Weight2100*cos(3.*phiGF)*sqrt(7.0)*0.125 *(1.+cosThetaGF)*(4.*cosThetaGF-1); double Weight2144=Weight2100*cos(4.*phiGF)*sqrt(14)*0.125 *sinThetaGF*sinThetaGFsq*(1.+cosThetaGF); dHist_H0040->Fill(myPipPimPi0Eta.M(),Weight*Weight0040); dHist_H2040->Fill(myPipPimPi0Eta.M(),Weight*Weight2040); dHist_H0041->Fill(myPipPimPi0Eta.M(),Weight*Weight0041); dHist_H2041->Fill(myPipPimPi0Eta.M(),Weight*Weight2041); dHist_H0042->Fill(myPipPimPi0Eta.M(),Weight*Weight0042); dHist_H2042->Fill(myPipPimPi0Eta.M(),Weight*Weight2042); dHist_H0043->Fill(myPipPimPi0Eta.M(),Weight*Weight0043); dHist_H2043->Fill(myPipPimPi0Eta.M(),Weight*Weight2043); dHist_H0044->Fill(myPipPimPi0Eta.M(),Weight*Weight0044); dHist_H2044->Fill(myPipPimPi0Eta.M(),Weight*Weight2044); dHist_H2240->Fill(myPipPimPi0Eta.M(),Weight*Weight2240); dHist_H2241->Fill(myPipPimPi0Eta.M(),Weight*Weight2241); dHist_H2242->Fill(myPipPimPi0Eta.M(),Weight*Weight2242); dHist_H2243->Fill(myPipPimPi0Eta.M(),Weight*Weight2243); dHist_H2244->Fill(myPipPimPi0Eta.M(),Weight*Weight2244); dHist_H2140->Fill(myPipPimPi0Eta.M(),Weight*Weight2140); dHist_H2141->Fill(myPipPimPi0Eta.M(),Weight*Weight2141); dHist_H2142->Fill(myPipPimPi0Eta.M(),Weight*Weight2142); dHist_H2143->Fill(myPipPimPi0Eta.M(),Weight*Weight2143); dHist_H2144->Fill(myPipPimPi0Eta.M(),Weight*Weight2144); } if (my2GammaMass>0.4 && my2GammaMass<0.45 && myPipPimPi0Mass>0.74 && myPipPimPi0Mass<0.82){ double mymass=0.5+(my2GammaMass-0.4); TLorentzVector mySB(my2Gamma.Vect(),sqrt(my2Gamma.Vect().Mag2()+mymass*mymass)); TLorentzVector myPipPimPi0Eta=mySB+myPipPimPi0; dHist_OmegaEtaMass_SB->Fill(mySB.M(),HalfWeight); } if (my2GammaMass>0.65 && my2GammaMass<0.7 && myPipPimPi0Mass>0.74 && myPipPimPi0Mass<0.82){ double mymass=0.5+(my2GammaMass-0.65); TLorentzVector mySB(my2Gamma.Vect(),sqrt(my2Gamma.Vect().Mag2()+mymass*mymass)); TLorentzVector myPipPimPi0Eta=mySB+myPipPimPi0; dHist_OmegaEtaMass_SB->Fill(mySB.M(),HalfWeight); } if (my2GammaMass>0.5 && my2GammaMass<0.6 && myPipPimPi0Mass>0.66 && myPipPimPi0Mass<0.7){ double mymass=0.74+(myPipPimPi0Mass-0.66); TLorentzVector mySB(myPipPimPi0.Vect(),sqrt(myPipPimPi0.Vect().Mag2()+mymass*mymass)); TLorentzVector myPipPimPi0Eta=mySB+my2Gamma; dHist_OmegaEtaMass_SB->Fill(myPipPimPi0Eta.M(),HalfWeight); } if (my2GammaMass>0.5 && my2GammaMass<0.6 && myPipPimPi0Mass>0.86 && myPipPimPi0Mass<0.9){ double mymass=0.74+(myPipPimPi0Mass-0.86); TLorentzVector mySB(myPipPimPi0.Vect(),sqrt(myPipPimPi0.Vect().Mag2()+mymass*mymass)); TLorentzVector myPipPimPi0Eta=mySB+my2Gamma; dHist_OmegaEtaMass_SB->Fill(myPipPimPi0Eta.M(),HalfWeight); } } } //E.g. Cut //if((locMissingMassSquared < -0.04) || (locMissingMassSquared > 0.04)) //{ // dComboWrapper->Set_IsComboCut(true); // continue; //} /****************************************** FILL FLAT TREE (IF DESIRED) ******************************************/ /* //FILL ANY CUSTOM BRANCHES FIRST!! Int_t locMyInt_Flat = 7; dFlatTreeInterface->Fill_Fundamental("flat_my_int", locMyInt_Flat); TLorentzVector locMyP4_Flat(4.0, 3.0, 2.0, 1.0); dFlatTreeInterface->Fill_TObject("flat_my_p4", locMyP4_Flat); for(int loc_j = 0; loc_j < locMyInt_Flat; ++loc_j) { dFlatTreeInterface->Fill_Fundamental("flat_my_int_array", 3*loc_j, loc_j); //2nd argument = value, 3rd = array index TLorentzVector locMyComboP4_Flat(8.0, 7.0, 6.0, 5.0); dFlatTreeInterface->Fill_TObject("flat_my_p4_array", locMyComboP4_Flat, loc_j); } */ //FILL FLAT TREE //Fill_FlatTree(); //for the active combo } // end of combo loop //FILL HISTOGRAMS: Num combos / events surviving actions Fill_NumCombosSurvivedHists(); /******************************************* LOOP OVER THROWN DATA (OPTIONAL) ***************************************/ /* //Thrown beam: just use directly if(dThrownBeam != NULL) double locEnergy = dThrownBeam->Get_P4().E(); //Loop over throwns for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i) { //Set branch array indices corresponding to this particle dThrownWrapper->Set_ArrayIndex(loc_i); //Do stuff with the wrapper here ... } */ /****************************************** LOOP OVER OTHER ARRAYS (OPTIONAL) ***************************************/ /* //Loop over beam particles (note, only those appearing in combos are present) for(UInt_t loc_i = 0; loc_i < Get_NumBeam(); ++loc_i) { //Set branch array indices corresponding to this particle dBeamWrapper->Set_ArrayIndex(loc_i); //Do stuff with the wrapper here ... } //Loop over charged track hypotheses for(UInt_t loc_i = 0; loc_i < Get_NumChargedHypos(); ++loc_i) { //Set branch array indices corresponding to this particle dChargedHypoWrapper->Set_ArrayIndex(loc_i); //Do stuff with the wrapper here ... } //Loop over neutral particle hypotheses for(UInt_t loc_i = 0; loc_i < Get_NumNeutralHypos(); ++loc_i) { //Set branch array indices corresponding to this particle dNeutralHypoWrapper->Set_ArrayIndex(loc_i); //Do stuff with the wrapper here ... } */ /************************************ 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 != "") Fill_OutputTree(); */ return kTRUE; } void DSelector_omegaeta::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 }