#include "DSelector_p2g.h" #include "particleType.h" //for topo analysis for MC void DSelector_p2g::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 = "myfile.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(); /******************************************** EXAMPLE USER INITIALIZATION *******************************************/ //DO WHATEVER YOU WANT HERE //EXAMPLE HISTOGRAM ACTIONS: dHistComboKinematics = new DHistogramAction_ParticleComboKinematics(dComboWrapper, false); //false: use measured data dHistComboPID = new DHistogramAction_ParticleID(dComboWrapper, false); //false: use measured data //change binning here dHistComboKinematics->Initialize(); dHistComboPID->Initialize(); //EXAMPLE CUT ACTIONS: //below: false: measured data, value: +/- N ns, Unknown: All PIDs, SYS_NULL: all timing systems dCutPIDDeltaT = new DCutAction_PIDDeltaT(dComboWrapper, false, 2.0, Unknown, SYS_NULL); dCutPIDDeltaT->Initialize(); //EXAMPLE MANUAL HISTOGRAMS: dHist_MissingMassSquared = new TH1I("MissingMassSquared", ";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.06, 0.06); dHist_BeamEnergy = new TH1I("BeamEnergy", ";Beam Energy (GeV)", 600, 0.0, 12.0); globEvent = 0; //For MC //Histogram for True Signal dHist_True_Signal = new TH1I("True_Signal", " ;True Signal", 3, -1, 2); //check for if there are true signal events //_cut6 means with cuts 0-6 //_cutBE means with cuts 0-6 and the coherent Beam Energy cut //For Data and MC //Delta T between Beam and RF time DeltaTBeamRF = new TH1D("DeltaTBeamRF", "Delta T between Beam and RF time; T_{beam}-T_{RF}; Counts/0.052ns", 1000, -26.0, 26.0); DeltaTBeamRF_pi0_cut6 = new TH1D("DeltaTBeamRF_pi0_cut6", "Delta T between Beam and RF time after cuts 0-6 for pi0; T_{beam}-T_{RF}; Counts/0.052ns", 1000, -26.0, 26.0); DeltaTBeamRF_pi0_cutBE = new TH1D("DeltaTBeamRF_pi0_cutBE", "Delta T between Beam and RF time after BE cut for pi0; T_{beam}-T_{RF}; Counts/0.052ns", 1000, -26.0, 26.0); DeltaTBeamRF_eta_cut6 = new TH1D("DeltaTBeamRF_eta_cut6", "Delta T between Beam and RF time after cuts 0-6 for eta; T_{beam}-T_{RF}; Counts/0.052ns", 1000, -26.0, 26.0); DeltaTBeamRF_eta_cutBE = new TH1D("DeltaTBeamRF_eta_cutBE", "Delta T between Beam and RF time after BE cut for eta; T_{beam}-T_{RF}; Counts/0.052ns", 1000, -26.0, 26.0); //Some distributions in the sideband of Delta_phi 10.0SetParameters(5.2, 3.0, 0.9); //dMinKinFitCL = 0.001; //dMinBeamEnergy = 8.4; //dMaxBeamEnergy = 9.0; //dMinPi0Mass = 0.11; //dMaxPi0Mass = 0.16; //For MC // Histogram of background decomposition locReactionLabels.push_back("TrueSignal"); locReactionLabels.push_back("Accidental"); locReactionLabels.push_back("p#pi^{+}#pi^{-}#pi^{0}"); locReactionLabels.push_back("n#pi^{+}#pi^{0}"); locReactionLabels.push_back("p#pi^{0}#pi^{0}"); locReactionLabels.push_back("p#gamma#pi^{0}"); locReactionLabels.push_back("p#pi^{+}#pi^{-}#eta"); locReactionLabels.push_back("n#pi^{+}#eta"); locReactionLabels.push_back("p#pi^{0}#eta"); locReactionLabels.push_back("p#gamma#eta"); locReactionLabels.push_back("Other"); hMass2gamma_Reaction = new TH2F("hMass2gamma_Reaction", "hMass2gamma_Reaction", locReactionLabels.size(), 0.5, float(locReactionLabels.size()) + 0.5, 100, 0., 1.5); for(size_t loc_i = 0; loc_i < locReactionLabels.size(); ++loc_i) hMass2gamma_Reaction->GetXaxis()->SetBinLabel(1 + loc_i, locReactionLabels[loc_i].c_str()); /***************************************** 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_p2g::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; globEvent++; if(globEvent % 100000 == 0) cout<0)//Just for MC for background study { for (UInt_t i=0;iSet_ArrayIndex(i); if(dThrownWrapper->Get_PID()==Proton) NumProton++; else if(dThrownWrapper->Get_PID()==Neutron) NumNeutron++; else if(dThrownWrapper->Get_PID()==Gamma) NumPhotons++; else if(dThrownWrapper->Get_PID()==Pi0) NumPi0++; else if(dThrownWrapper->Get_PID()==PiPlus) NumPiPlus++; else if(dThrownWrapper->Get_PID()==PiMinus) NumPiMinus++; else if(dThrownWrapper->Get_PID()==Eta) NumEta++; else if(dThrownWrapper->Get_PID()==EtaPrime) NumEtaPrime++; else NumOther++; } // find index of thrown proton (if it exists) for (UInt_t i=0; iSet_ArrayIndex(i); if(dThrownWrapper->Get_PID()==Proton){ Proton_Index = i; break; } } } /**************************************** SETUP AUTOMATIC UNIQUENESS TRACKING ***************************************/ //Reset uniqueness tracking for each action dHistComboKinematics->Reset_NewEvent(); dHistComboPID->Reset_NewEvent(); //INSERT OTHER USER ACTIONS HERE /***************************************** SETUP MANUAL UNIQUENESS TRACKING *****************************************/ //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 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 /* // Reject true pi0 events from sim1 (ie. bggen) since adding signal from AmpTools generator if(Get_IsThrownTopology() && dOption.Contains("sim1")) { dComboWrapper->Set_IsComboCut(true); continue; } */ /********************************************** Topo of background information *********************************************/ UInt_t reaction = 1; if(Get_NumThrown()>0)//Just for MC { Int_t Photon1__MatchID = -1; Int_t Photon2__MatchID = -1; Int_t Thrown__PID_Photon1 = -1; Int_t Thrown__PID_Photon2 = -1; Photon1__MatchID = dPhoton1Wrapper->Get_ThrownIndex(); Photon2__MatchID = dPhoton2Wrapper->Get_ThrownIndex(); Int_t Proton__MatchID = dProtonWrapper->Get_ThrownIndex(); // check all requirements for trueSignal if(Photon1__MatchID>=0 && Photon2__MatchID>=0 && Proton__MatchID>=0 && Proton_Index>=0){ dThrownWrapper->Set_ArrayIndex(Photon1__MatchID); Thrown__PID_Photon1=dThrownWrapper->Get_PID(); TLorentzVector locThrownPhoton1__P4 = dThrownWrapper->Get_P4(); dThrownWrapper->Set_ArrayIndex(Photon2__MatchID); Thrown__PID_Photon2=dThrownWrapper->Get_PID(); TLorentzVector locThrownPhoton2__P4 = dThrownWrapper->Get_P4(); dThrownWrapper->Set_ArrayIndex(Proton__MatchID); TLorentzVector locThrownProton__P4 = dThrownWrapper->Get_P4(); if(Thrown__PID_Photon1==Gamma && Thrown__PID_Photon2==Gamma){ Bool_t TrueBeamParticle = dThrownBeam->Get_IsGenerator(); TLorentzVector thrownBeam__P4 = dComboBeamWrapper->Get_P4(); // identify exclusive p2gamma channel TLorentzVector thrownMissingP4(0,0,thrownBeam__P4.Energy(),thrownBeam__P4.Energy()+0.938); thrownMissingP4 -= TLorentzVector(locThrownPhoton1__P4.Vect(),locThrownPhoton1__P4.E()); thrownMissingP4 -= TLorentzVector(locThrownPhoton2__P4.Vect(),locThrownPhoton2__P4.E()); thrownMissingP4 -= TLorentzVector(locThrownProton__P4.Vect(),locThrownProton__P4.E()); if(NumPhotons == 2 && fabs(thrownMissingP4.M2())<0.01 && TrueBeamParticle) { trueSignal = true; } } } //Fill the histogram for true signal dHist_True_Signal->Fill(trueSignal); // identify reactions with significant background contributions if(trueSignal == false){ if(NumProton==1 && NumNeutron==0 && NumPhotons==2 && (NumPi0==1 || NumEta==1 || NumEtaPrime==1) && NumPiPlus==0 && NumPiMinus==0 && NumOther==0) reaction = 2; // Accidentals with correct topology else if(NumProton==1 && NumNeutron==0 && NumPhotons==2 && NumPi0==1 && NumPiPlus==1 && NumPiMinus==1) reaction = 3; // p pi+ pi- pi0 else if(NumProton==0 && NumNeutron==1 && NumPhotons==2 && NumPi0==1 && NumPiPlus==1 && NumPiMinus==0) reaction = 4; // n pi+ pi0 else if(NumProton==1 && NumNeutron==0 && NumPhotons==4 && NumPi0==2 && NumPiPlus==0 && NumPiMinus==0) reaction = 5; // p pi0 pi0 else if(NumProton==1 && NumNeutron==0 && NumPhotons==3 && NumPi0==1 && NumPiPlus==0 && NumPiMinus==0) reaction = 6; // p pi0 gamma else if(NumProton==1 && NumNeutron==0 && NumPhotons==2 && NumPi0==0 && NumEta==1 && NumPiPlus==1 && NumPiMinus==1) reaction = 7; // p pi+ pi- eta else if(NumProton==0 && NumNeutron==1 && NumPhotons==2 && NumPi0==1 && NumPiPlus==1 && NumPiMinus==0) reaction = 8; // n pi+ eta else if(NumProton==1 && NumNeutron==0 && NumPhotons==4 && NumPi0==1 && NumEta==1 && NumPiPlus==0 && NumPiMinus==0) reaction = 9; // p pi0 eta else if(NumProton==1 && NumNeutron==0 && NumPhotons==3 && NumPi0==0 && NumEta==1 && NumPiPlus==0 && NumPiMinus==0) reaction = 10; // p eta gamma else reaction = 11; // other } } /********************************************** GET PARTICLE INDICES *********************************************/ //cout<<"before step 0"<Get_BeamID(); Int_t locPhoton1NeutralID = dPhoton1Wrapper->Get_NeutralID(); Int_t locPhoton2NeutralID = dPhoton2Wrapper->Get_NeutralID(); 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 locPhoton1P4 = dPhoton1Wrapper->Get_P4(); TLorentzVector locPhoton2P4 = dPhoton2Wrapper->Get_P4(); TLorentzVector locProtonP4 = dProtonWrapper->Get_P4(); // Get Measured P4's: //Step 0 TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured(); TLorentzVector locPhoton1P4_Measured = dPhoton1Wrapper->Get_P4_Measured(); TLorentzVector locPhoton2P4_Measured = dPhoton2Wrapper->Get_P4_Measured(); TLorentzVector locProtonP4_Measured = dProtonWrapper->Get_P4_Measured(); /*********************************************** GET Space-Time **********************************************/ TLorentzVector locBeam_X4_Measured = dComboBeamWrapper->Get_X4_Measured(); TLorentzVector locProton_X4_Measured = dProtonWrapper->Get_X4_Measured(); //The space-time information when the photons hit the detector TLorentzVector locPhoton1_X4_Measured = dPhoton1Wrapper->Get_X4_Shower(); TLorentzVector locPhoton2_X4_Measured = dPhoton2Wrapper->Get_X4_Shower(); //The space-time information when the photons hit the detector TLorentzVector X4_Photon1 = locPhoton1_X4_Measured; double x_g1 = X4_Photon1.X(); double y_g1 = X4_Photon1.Y(); double z_g1 = X4_Photon1.Z(); double r_g1 = X4_Photon1.Perp(); double t_g1 = X4_Photon1.T(); TLorentzVector X4_Photon2 = locPhoton2_X4_Measured; double x_g2 = X4_Photon2.X(); double y_g2 = X4_Photon2.Y(); double z_g2 = X4_Photon2.Z(); double r_g2 = X4_Photon2.Perp(); double t_g2 = X4_Photon2.T(); //the distant between two photons in XY plane, used for splits cut in FCAL double dis_g1_g2 = sqrt((x_g1-x_g2)*(x_g1-x_g2)+(y_g1-y_g2)*(y_g1-y_g2)); //define the detector photon angle TVector3 X3_Photon1 = X4_Photon1.Vect(); TVector3 X3_Photon2 = X4_Photon2.Vect(); TVector3 X3_target(0.0, 0.0, 65.0); double theta_det_g1 = (X3_Photon1-X3_target).Theta(); double theta_det_g2 = (X3_Photon2-X3_target).Theta(); double theta_det_gmin = theta_det_g1; if(theta_det_g1>theta_det_g2) theta_det_gmin = theta_det_g2; double theta_det_g_min = theta_det_gmin*180/TMath::Pi(); /********************************************* RF time and Beam time **********************************************/ double tRF = 0.0; tRF = dComboWrapper->Get_RFTime_Measured(); TLorentzVector P4_Beam = locBeamP4_Measured; TLorentzVector X4_Beam = locBeam_X4_Measured; double tbeam = 0.0; tbeam = X4_Beam.T(); /*********************************************** BCAL & FCAL Energy **********************************************/ double Ephoton1BCAL = dPhoton1Wrapper->Get_Energy_BCAL(); double Ephoton2BCAL = dPhoton2Wrapper->Get_Energy_BCAL(); double Ephoton1FCAL = dPhoton1Wrapper->Get_Energy_FCAL(); double Ephoton2FCAL = dPhoton2Wrapper->Get_Energy_FCAL(); //set the minimum photon energy cut: 0.3 GeV in the BCAL & 0.5 GeV in the FCAL if(Ephoton1BCAL>0) { if(Ephoton1BCAL<0.3) continue; } if(Ephoton1FCAL>0) { if(Ephoton1FCAL<0.5) continue; } if(Ephoton2BCAL>0) { if(Ephoton2BCAL<0.3) continue; } if(Ephoton2FCAL>0) { if(Ephoton2FCAL<0.5) continue; } /****************************************** Shift FCAL energy for photon *****************************************/ // DO YOUR STUFF HERE TLorentzVector P4_Proton = locProtonP4_Measured; TLorentzVector X4_Proton = locProton_X4_Measured; TLorentzVector P4_gamma1 = locPhoton1P4_Measured; TLorentzVector P4_gamma2 = locPhoton2P4_Measured; //Make an FCAL energy shift for photons by a factor of 1.0/0.992 for data double scale = 1.0/0.992; if(Get_NumThrown()==0)//only for Data { double E_g1 = locPhoton1P4_Measured.E(); double px_g1 = locPhoton1P4_Measured.Px(); double py_g1 = locPhoton1P4_Measured.Py(); double pz_g1 = locPhoton1P4_Measured.Pz(); double E_g2 = locPhoton2P4_Measured.E(); double px_g2 = locPhoton2P4_Measured.Px(); double py_g2 = locPhoton2P4_Measured.Py(); double pz_g2 = locPhoton2P4_Measured.Pz(); if(Ephoton1FCAL>0) { P4_gamma1.SetPxPyPzE(scale*px_g1, scale*py_g1, scale*pz_g1, scale*E_g1); } if(Ephoton2FCAL>0) { P4_gamma2.SetPxPyPzE(scale*px_g2, scale*py_g2, scale*pz_g2, scale*E_g2); } } /********************************************* COMBINE FOUR-MOMENTUM ********************************************/ // Combine 4-vectors for Missing momentum for gamma p -> gamma gamma p TLorentzVector locMissingP4_Measured = locBeamP4_Measured + dTargetP4; locMissingP4_Measured -= P4_gamma1 + P4_gamma2 + locProtonP4_Measured; TLorentzVector P4_missing = locMissingP4_Measured; //Combine 4-vectors for 2gamma TLorentzVector P4_2gamma = P4_gamma1 + P4_gamma2; // Combine 4-vectors for Missing momentum for gamma p -> Xp TLorentzVector P4_missing_omega = locBeamP4_Measured + dTargetP4; P4_missing_omega -= locProtonP4_Measured; //The measured 4-vectors of momentum transform TLorentzVector P4_Trans = P4_Proton; P4_Trans -= dTargetP4; /********************************************* Kinematic varias ********************************************/ //vertex information double z = 0.0; double r = 0.0; z = X4_Proton.Z(); r = X4_Proton.Perp(); //theta of proton double theta_proton = P4_Proton.Theta(); //phi of proton and photons double phi_proton = P4_Proton.Phi(); double phi_gamma1 = P4_gamma1.Phi(); double phi_gamma2 = P4_gamma2.Phi(); //energy of photons double energy_gamma1 = P4_gamma1.E(); double energy_gamma2 = P4_gamma2.E(); //dEdx in CDC and the momentum of the proton double mp = 0.938272; double pp = 0.0; double dEdx = 0.0; pp = P4_Proton.P(); dEdx = dProtonWrapper->Get_dEdx_CDC()*1e6; //Missing Mass of Omega double mmo = P4_missing_omega.M(); //Missing Mass double mm = P4_missing.M(); //Missing Mass squared double mmsq = P4_missing.M2(); //Invariant Mass of 2 gamma double im_2g = P4_2gamma.M(); // calculate unused energy //cout<Set_ArrayIndex(loc_i); if(dNeutralHypoWrapper->Get_NeutralID() != locPhoton1NeutralID && dNeutralHypoWrapper->Get_NeutralID() != locPhoton2NeutralID && dNeutralHypoWrapper->Get_PID()==1){ unusedEnergyBCAL += dNeutralHypoWrapper->Get_Energy_BCAL(); unusedEnergyFCAL += dNeutralHypoWrapper->Get_Energy_FCAL(); unusedEnergy = unusedEnergyBCAL + unusedEnergyFCAL; } //Do stuff with the wrapper here ... } //Beam Energy double be = P4_Beam.E(); //Missing Energy double me = P4_missing.E(); //phi_2gamma double phi_2g = (180.0/TMath::Pi())*(P4_2gamma.Phi()); //phi_proton double phi_p = (180/TMath::Pi())*(P4_Proton.Phi()); //Delta_phi double delta_phi = phi_2g - phi_p; if(delta_phi>360.) delta_phi += -360.; if(delta_phi<0.) delta_phi += 360.; //t double t = P4_Trans.M2(); /**************************************** PreSelection ******************************************/ if(pp<0.25) continue; //cut in momentum of proton if(z<49.0 || z>76.0 || r>1.0) continue; //cut in vertex if(dis_g1_g2<12.5) continue; //cut for FCAL splits if((180.0/TMath::Pi())*theta_det_g1<2.5 || (180.0/TMath::Pi())*theta_det_g2<2.5) continue; //cut for beam hole if(((180.0/TMath::Pi())*theta_det_g1<11.5 && (180.0/TMath::Pi())*theta_det_g1>10.3) || ((180.0/TMath::Pi())*theta_det_g2<11.5 && (180.0/TMath::Pi())*theta_det_g2>10.3)) continue; //cut for FCAL-BCAL gap /**************************************** cut condition ******************************************/ bool Cuts[ncuts] = {true}; bool CutsAdd[ncuts] = {true}; for (int i=1; ifMinProton_dEdx->Eval(pp/mp)) ? true : false; //dEdx cut (cut0) Cuts[2]=(abs(delta_phi-180.0)<5.0) ? true : false; //delta_phi cut (cut1) Cuts[3]=(mmsq<0.01&&mmsq>-0.01) ? true : false; //MMsq cut (cut2) Cuts[4]=(me>-0.50&&me<0.70) ? true : false; //ME cut (cut3) Cuts[5]=(be>4) ? true : false; //BE cut (cut4) Cuts[6]=(unusedEnergy<0.01) ? true : false; //Unused Energy cut (cut5) Cuts[7]=(mmo<0.5) ? true : false; //MMXp cut (cut6) for pi0 Cuts[8]=(be>8.4&&be<9.0) ? true : false; //BE cut (cutBE) Cuts[9]=(abs(im_2g-0.135)<0.021) ? true : false; // IM cut (cutpi) for pi0 Cuts[10]=(mmo<0.7) ? true : false; //MMXp cut (cut6) for eta Cuts[11]=(be>8.4&&be<9.0) ? true : false; //BE cut (cutBE) Cuts[12]=(abs(im_2g-0.546)<0.064) ? true : false; // IM cut (cuteta) for eta if(Get_NumThrown()>0)//Only for MC { Cuts[4]=(me>-0.60&&me<0.60) ? true : false; //ME cut (cut3) } for (int i=1; i<10; i++) //i=1-6 for both pi0 and eta; i=7-9 for pi0 { CutsAdd[i]=CutsAdd[i-1]&&Cuts[i]; } CutsAdd[10] = CutsAdd[6]&&Cuts[10]; for (int i=11; i-0.05) ? true : false; //Open MMsq cut (cut2) OpenCuts[2][0] = (me>-1.0&&me<1.0) ? true : false; //Open ME cut (cut3) OpenCuts[3][0] = true; //Open MMXp cut (cut6) OpenCuts[4][0] = true; //Open M_2g cut (pi0/eta) int indexpi0cuts[5] = {2, 3, 4, 7, 9}; //opened cut index for pi0 int indexetacuts[5] = {2, 3, 4, 10, 12}; //opened cut index for eta bool OpenCutsadd[5][2] = {true}; for(int i=0; i<5; i++) { OpenCuts[i][1] = OpenCuts[i][0]; OpenCutsadd[i][0] = OpenCuts[i][0]; OpenCutsadd[i][1] = OpenCuts[i][0]; } //for pi0 for(int j=0; j<5; j++) { for(int i=1; i<10; i++) { if(i!=indexpi0cuts[j]) OpenCutsadd[j][0] = OpenCutsadd[j][0] && Cuts[i]; } } //for eta for(int j=0; j<5; j++) { for(int i=1; i9)) OpenCutsadd[j][1] = OpenCutsadd[j][1] && Cuts[i]; } } bool TopoCuts[7] = {true}; //topo cut for background study if(Get_NumThrown()>0)//Only for MC { TopoCuts[0]=(reaction == 1) ? true : false; //True Signal TopoCuts[1]=(reaction == 2) ? true : false; //Accidentals with correct topology TopoCuts[2]=(reaction == 3) ? true : false; //p pi+ pi- pi0 TopoCuts[3]=(reaction == 5) ? true : false; //p pi0 pi0 TopoCuts[4]=(reaction == 6) ? true : false; //p pi0 gamma TopoCuts[5]=(reaction == 7) ? true : false; // p pi+ pi- eta TopoCuts[6]=(reaction == 9) ? true : false; // p pi0 eta } bool NoUEpi0Cut = CutsAdd[5]&&Cuts[7]; //for pi0 mmo<0.5 bool NoUEetaCut = CutsAdd[5]&&Cuts[10]; //for eta mmo<0.7 bool NoBEpi0Cut = CutsAdd[7]&&Cuts[9]; //with cuts 0-6 and pi0 mass cut bool NoMMXpBECut = CutsAdd[6]&&Cuts[8]; //with cuts 0-5 and coherent Beam Energy cut bool NoMMXppi0Cut = CutsAdd[6]&&Cuts[8]&&Cuts[9];//without MMXp cut for pi0 bool NoMMXpetaCut = CutsAdd[6]&&Cuts[8]&&Cuts[12];//without MMXp cut for eta bool DPSBcut = (abs(delta_phi-180.0)<60.0&&abs(delta_phi-180.0)>10.0) ? true : false; //delta_phi sideband cut (cut1) bool DPSBcut5 = Cuts[1]&&DPSBcut; //cuts 0-6 for pi0 in the delta_phi sideband bool DPSBpi0cut, DPSBetacut; //cuts 0-6 for pi0 or eta in the delta_phi sideband for (int i=3; i<7; i++) { DPSBcut5 = DPSBcut5&&Cuts[i]; } DPSBpi0cut = DPSBcut5&&Cuts[7]; DPSBetacut = DPSBcut5&&Cuts[10]; /**************************************** HISTOGRAM ******************************************/ //Fill Delta T between Beam and RF time with cuts 0-6 (and the coherent Beam Energy cut) for pi0 and eta DeltaTBeamRF->Fill(tbeam-tRF); if(CutsAdd[7]) { DeltaTBeamRF_pi0_cut6->Fill(tbeam-tRF); if(Cuts[8]) { DeltaTBeamRF_pi0_cutBE->Fill(tbeam-tRF); } } if(CutsAdd[10]) { DeltaTBeamRF_eta_cut6->Fill(tbeam-tRF); if(Cuts[11]) { DeltaTBeamRF_eta_cutBE->Fill(tbeam-tRF); } } //Fill the histograms(Delta_phi, MMsq, ME, MMXp, IM_2gamma) with one perticular cut opened or loosed if(abs(tbeam-tRF)<(0.5*4.008)) { for(int j=0; j<2; j++)//j=0 for pi0; j=1 for eta { if(OpenCutsadd[0][j]) //open Delta_phi cut { S_Delta_phi[j]->Fill(delta_phi); } if(OpenCutsadd[1][j]) //open MMsq cut { S_MMsq[j]->Fill(mmsq); } if(OpenCutsadd[2][j]) //open ME cut { S_ME[j]->Fill(me); } if(OpenCutsadd[3][j]) //open MMXp cut { S_MMXp[j]->Fill(mmo); } if(OpenCutsadd[4][j]) //open IM_2g cut { S_IM_2gamma[j]->Fill(im_2g); } if(Get_NumThrown()>0)//For MC (the background study) { //k[7]={True Signal, Accidentals, p pi+ pi- pi0, p pi0 pi0, p pi0 gamma, p pi+ pi- eta, p pi0 eta} for(int k=0; k<7; k++)//for 7 topos { if(TopoCuts[k]) { if(OpenCutsadd[0][j])//open Delta_phi cut { T_Delta_phi[j][k]->Fill(delta_phi); } if(OpenCutsadd[1][j])//open MMsq cut { T_MMsq[j][k]->Fill(mmsq); } if(OpenCutsadd[2][j])//open ME cut { T_ME[j][k]->Fill(me); } if(OpenCutsadd[3][j])//open MMXp cut { T_MMXp[j][k]->Fill(mmo); } if(OpenCutsadd[4][j])//open IM_2g cut { T_IM_2gamma[j][k]->Fill(im_2g); } } } } } } //Fill some distributions in the central peak of tbeam-tRF with perticular cuts if(abs(tbeam-tRF)<(0.5*4.008)) { if(NoUEpi0Cut)//without unused energy cut for pi0 { NU_Unused_Energy_pi0->Fill(unusedEnergy); NU_Unused_Energy_1_pi0->Fill(unusedEnergy); } if(NoUEetaCut)//without unused energy cut for eta { NU_Unused_Energy_eta->Fill(unusedEnergy); NU_Unused_Energy_1_eta->Fill(unusedEnergy); } if(NoMMXppi0Cut)//without MMXp cut for pi0 { NM_MMXp_pi0->Fill(mmo); } if(NoMMXpetaCut)//without MMXp cut for eta { NM_MMXp_eta->Fill(mmo); } if(DPSBpi0cut)//with cuts 0-6 for pi0 in the delta_phi sideband { DPSB_Delta_phi_pi0->Fill(delta_phi); if(Cuts[8]&&Cuts[9])//with cuts 0-6 and BE&&M_2g cuts for pi0 in the delta_phi sideband { DPSB_t_pi0->Fill(-t); DPSB_Phi_Proton_pi0->Fill(phi_proton); DPSB_Phi_Proton_vs_t_pi0->Fill(-t, phi_proton); } } if(DPSBetacut)//with cuts 0-6 for eta in the delta_phi sideband { DPSB_Delta_phi_eta->Fill(delta_phi); if(Cuts[11]&&Cuts[12])//with cuts 0-6 and BE&&M_2g cuts for eta in the delta_phi sideband { DPSB_t_eta->Fill(-t); DPSB_Phi_Proton_eta->Fill(phi_proton); DPSB_Phi_Proton_vs_t_eta->Fill(-t, phi_proton); } } } if(NoBEpi0Cut)// with cuts 0-6 and pi0 mass cut { if(abs(tbeam-tRF)<(0.5*4.008)) { Phi_Proton_vs_BE_cut6->Fill(be, phi_proton); } else { Acci_Phi_Proton_vs_BE_cut6->Fill(be, phi_proton); } } if(NoMMXpBECut)//with cuts 0-5 and coherent Beam Energy cut { if(abs(tbeam-tRF)<(0.5*4.008)) { Phi_Proton_vs_IM2g_cutBE->Fill(im_2g, phi_proton); } else { Acci_Phi_Proton_vs_IM2g_cutBE->Fill(im_2g, phi_proton); } } //Fill the histograms for Missing Energy(ME), Invariant Mass of 2gamma(IM_2gamma), Missing Mass squared(MMsq), Missing Mass of Xp process(MMXp), Unused Energy(UE), Beam Energy(BE), Delta phi of proton(Delta_phi) for(int i=0; iFill(me); IM_2gamma[i]->Fill(im_2g); MMsq[i]->Fill(mmsq); MMXp[i]->Fill(mmo); UE[i]->Fill(unusedEnergy); BE[i]->Fill(be); Delta_phi[i]->Fill(delta_phi); dEdx_vs_p_Proton[i]->Fill(pp/mp, dEdx); dEdx_vs_p_Proton_vs_Phi_Proton[i]->Fill(phi_proton, pp/mp, dEdx); } } else { if(CutsAdd[i]) { Acci_ME[i]->Fill(me); Acci_IM_2gamma[i]->Fill(im_2g); Acci_MMsq[i]->Fill(mmsq); Acci_MMXp[i]->Fill(mmo); Acci_UE[i]->Fill(unusedEnergy); Acci_BE[i]->Fill(be); Acci_Delta_phi[i]->Fill(delta_phi); Acci_dEdx_vs_p_Proton[i]->Fill(pp/mp, dEdx); Acci_dEdx_vs_p_Proton_vs_Phi_Proton[i]->Fill(phi_proton, pp/mp, dEdx); } } } if(abs(tbeam-tRF)<(0.5*4.008)) { //Fill some histograms(theta_Photon_vs_phi_Photon, theta_Proton_vs_phi_Proton, Phi_Proton_vs_t) with cuts for cuts0-6, without BE and with BE cut for pi0 and eta for(int i=0; iFill(phi_gamma1, theta_det_g1); theta_Photon2_vs_phi_Photon2[i]->Fill(phi_gamma2, theta_det_g2); theta_Proton_vs_phi_Proton[i]->Fill(phi_proton, theta_proton); Phi_Proton_vs_t[i]->Fill(-t, phi_proton); Phi_Proton_vs_t_5[i]->Fill(-t, phi_proton); Phi_Proton_vs_t_15[i]->Fill(-t, phi_proton); } } //Fill some distributions (ME_vs_Phi_Proton, ME_vs_theta_Photon, X_Y_FCAL_g1, X_Y_FCAL_g2, E_g1_vs_E_g2) with cuts 0-6 for pi0 if(CutsAdd[7]) //cuts 0-6 for pi0 { if(Get_NumThrown()>0)//Only for MC { // fill histogram of reaction type for background study hMass2gamma_Reaction->Fill(reaction, im_2g); } ME_vs_Phi_Proton_cut6->Fill(phi_proton, me); ME_vs_theta_Photon_min_all_cut6->Fill(theta_det_g_min, me); if(theta_det_g_min<20.0&&theta_det_g_min>0.0) ME_vs_theta_Photon_min_20_cut6->Fill(theta_det_g_min, me); if(theta_det_g_min<4.0&&theta_det_g_min>0.0) ME_vs_theta_Photon_min_cut6->Fill(theta_det_g_min, me); if(Ephoton1FCAL>0 && Ephoton2FCAL>0) { X_Y_FCAL_g1->Fill(y_g1, x_g1); X_Y_FCAL_g2->Fill(y_g2, x_g2); E_g1_vs_E_g2->Fill(energy_gamma2, energy_gamma1); } } //Fill the -t, Phi_Proton distributions with all cuts for pi0 if(CutsAdd[9]) //all cuts for pi0 { t_pi0->Fill(-t); Phi_Proton_pi0->Fill(phi_proton); } //Fill the -t, Phi_Proton distributions with all cuts for eta if(CutsAdd[12]) //all cuts for eta { t_eta->Fill(-t); Phi_Proton_eta->Fill(phi_proton); } } else { //Fill some histograms(theta_Photon_vs_phi_Photon, theta_Proton_vs_phi_Proton, Phi_Proton_vs_t) with cuts for cuts0-6, without BE and with BE cut for pi0 and eta for accidentals for(int i=0; iFill(phi_gamma1, theta_det_g1); Acci_theta_Photon2_vs_phi_Photon2[i]->Fill(phi_gamma2, theta_det_g2); Acci_theta_Proton_vs_phi_Proton[i]->Fill(phi_proton, theta_proton); Acci_Phi_Proton_vs_t[i]->Fill(-t, phi_proton); Acci_Phi_Proton_vs_t_5[i]->Fill(-t, phi_proton); Acci_Phi_Proton_vs_t_15[i]->Fill(-t, phi_proton); } } //Fill some distributions (ME_vs_Phi_Proton, ME_vs_theta_Photon, X_Y_FCAL_g1, X_Y_FCAL_g2, E_g1_vs_E_g2) with cuts 0-6 for pi0 for accidentals if(CutsAdd[7]) //cuts 0-6 for pi0 { Acci_ME_vs_Phi_Proton_cut6->Fill(phi_proton, me); Acci_ME_vs_theta_Photon_min_all_cut6->Fill(theta_det_g_min, me); if(theta_det_g_min<20.0&&theta_det_g_min>0.0) Acci_ME_vs_theta_Photon_min_20_cut6->Fill(theta_det_g_min, me); if(theta_det_g_min<4.0&&theta_det_g_min>0.0) Acci_ME_vs_theta_Photon_min_cut6->Fill(theta_det_g_min, me); if(Ephoton1FCAL>0 && Ephoton2FCAL>0) { Acci_X_Y_FCAL_g1->Fill(y_g1, x_g1); Acci_X_Y_FCAL_g2->Fill(y_g2, x_g2); Acci_E_g1_vs_E_g2->Fill(energy_gamma2, energy_gamma1); } } //Fill the -t, Phi_Proton distributions with all cuts for pi0 for accidentals if(CutsAdd[9]) //all cuts for pi0 { Acci_t_pi0->Fill(-t); Acci_Phi_Proton_pi0->Fill(phi_proton); } //Fill the -t, Phi_Proton distributions with all cuts for eta for accidentals if(CutsAdd[12]) //all cuts for eta { Acci_t_eta->Fill(-t); Acci_Phi_Proton_eta->Fill(phi_proton); } } /**************************************** EXAMPLE: HISTOGRAM KINEMATICS ******************************************/ dHistComboKinematics->Perform_Action(); dHistComboPID->Perform_Action(); /**************************************** EXAMPLE: PID CUT ACTION ************************************************/ /* if(!dCutPIDDeltaT->Perform_Action()) { dComboWrapper->Set_IsComboCut(true); continue; } */ /**************************************** 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[Gamma].insert(locPhoton1NeutralID); locUsedThisCombo_MissingMass[Gamma].insert(locPhoton2NeutralID); locUsedThisCombo_MissingMass[Proton].insert(locProtonTrackID); //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); } //E.g. Cut //if((locMissingMassSquared < -0.04) || (locMissingMassSquared > 0.04)) // continue; //could also mark combo as cut, then save cut results to a new TTree } /******************************************* 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 (all are present, even those not in any combos) 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 (all are present, even those not in any combos) 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 != "") FillOutputTree(); */ return kTRUE; } void DSelector_p2g::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 }