Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DSelector_Z2pi_trees_gen.C
Go to the documentation of this file.
1 #include "DSelector_Z2pi_trees.h"
2 
3 void DSelector_Z2pi_trees::Init(TTree *locTree)
4 {
5 
6  // Note: this script is modified from DSelector_Z2pi_trees.C in order to generate an output file for "tagged" events.
7  // It will eliminate events that do not have a hit in the tagger. Otherwise the event is kept.
8  // Elton 5/9/2018
9 
10 
11  // USERS: IN THIS FUNCTION, ONLY MODIFY SECTIONS WITH A "USER" OR "EXAMPLE" LABEL. LEAVE THE REST ALONE.
12 
13  // The Init() function is called when the selector needs to initialize a new tree or chain.
14  // Typically here the branch addresses and branch pointers of the tree will be set.
15  // Init() will be called many times when running on PROOF (once per file to be processed).
16 
17  //USERS: SET OUTPUT FILE NAME //can be overriden by user in PROOF
18  dOutputFileName = "DSelector_Z2pi_trees.root"; //"" for none
19  dOutputTreeFileName = "tree_DSelector_Z2pi_trees.root"; //"" for none
20  dFlatTreeFileName = ""; //output flat tree (one combo per tree entry), "" for none
21  dFlatTreeName = ""; //if blank, default name will be chosen
22 
23  //Because this function gets called for each TTree in the TChain, we must be careful:
24  //We need to re-initialize the tree interface & branch wrappers, but don't want to recreate histograms
25  bool locInitializedPriorFlag = dInitializedFlag; //save whether have been initialized previously
26  DSelector::Init(locTree); //This must be called to initialize wrappers for each new TTree
27  //gDirectory now points to the output file with name dOutputFileName (if any)
28  if(locInitializedPriorFlag)
29  return; //have already created histograms, etc. below: exit
30 
33 
34 
35  // EXAMPLE CUT PARAMETERS:
36  // fMinProton_dEdx = new TF1("fMinProton_dEdx", "exp(-1.*[0]*x + [1]) + [2]", 0., 10.);
37  // fMinProton_dEdx->SetParameters(4.0, 2.5, 1.25);
38  fMaxPion_dEdx = new TF1("fMaxPion_dEdx", "exp(-1.*[0]*x + [1]) + [2]", 0., 10.);
39  fMaxPion_dEdx->SetParameters(4.0, 2.0, 2.5);
40  dMinKinFitCL = 5.73303e-7; //5.73303e-7;
41  dMaxKinFitChiSq = 5.0;
42  dMinBeamEnergy = 5.5;
43  dMaxBeamEnergy = 6.0;
44  dMin2piMass = 0.2;
45  dMax2piMass = 0.6;
48 
49  /*********************************** EXAMPLE USER INITIALIZATION: ANALYSIS ACTIONS **********************************/
50 
51  //ANALYSIS ACTIONS: //Executed in order if added to dAnalysisActions
52  //false/true below: use measured/kinfit data
53 
54  //PID
55  dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false));//false: use measured data
56  dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, true, "KinFit")); //true: use kinfit data;
57  //below: value: +/- N ns, Unknown: All PIDs, SYS_NULL: all timing systems
58  //dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.5, KPlus, SYS_BCAL));
59 
60  //MASSES
61  //dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, Lambda, 1000, 1.0, 1.2, "Lambda"));
62  //dAnalysisActions.push_back(new DHistogramAction_MissingMassSquared(dComboWrapper, false, 1000, -0.1, 0.1));
63 
64  //KINFIT RESULTS
65  dAnalysisActions.push_back(new DHistogramAction_KinFitResults(dComboWrapper));
66 
67  //CUT MISSING MASS
68  //dAnalysisActions.push_back(new DCutAction_MissingMassSquared(dComboWrapper, false, -0.03, 0.02));
69 
70  //BEAM ENERGY
71  //dAnalysisActions.push_back(new DHistogramAction_BeamEnergy(dComboWrapper, false));
72  dAnalysisActions.push_back(new DCutAction_BeamEnergy(dComboWrapper, false, dMinBeamEnergy ,dMaxBeamEnergy ));
73 
74  //KINEMATICS
75  dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false));
76 
77  //INITIALIZE ACTIONS
78  //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
79  Initialize_Actions();
80 
81  /******************************** EXAMPLE USER INITIALIZATION: STAND-ALONE HISTOGRAMS *******************************/
82 
83  dHist_MissingMassSquared = new TH1I("MissingMassSquared", ";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.24, 0.24);
84  dHist_BeamEnergy = new TH1I("BeamEnergy", ";Beam Energy (GeV)", 600, 0.0, 12.0);
85  // dHist_pMomentumMeasured = new TH1I("pMomentumMeasured", ";p Momentum Measured (GeV)", 100, 0.0, 2);
86  dHist_piplusMomentumMeasured = new TH1I("piplusMomentumMeasured", ";#pi^{+} Momentum Measured (GeV)", 600, 0.0, 12);
87  dHist_piminusMomentumMeasured = new TH1I("piminusMomentumMeasured", ";#pi^{-} Momentum Measured (GeV)", 600, 0.0, 12);
88 
89  // dHist_Proton_dEdx_P = new TH2I("Proton_dEdx_P", " ;p_{proton} GeV/c; dE/dx (keV/cm)", 250, 0.0, 5.0, 250, 0.0, 25.);
90  dHist_KinFitChiSq = new TH1I("KinFitChiSq", ";Kinematic Fit #chi^{2}/NDF", 250, 0., 25.);
91  dHist_KinFitCL = new TH1I("KinFitCL", ";Kinematic Fit Confidence Level", 100, 0., 1.);
92 
93  dHist_M2pigen = new TH1I("M2pigen", ";M_{#pi^{+}#pi^{-}} Gen (GeV/c^{2})", 400, 0.2, 0.6);
94  dHist_M2pikin = new TH1I("M2pikin", ";M_{#pi^{+}#pi^{-}} Kin (GeV/c^{2})", 400, 0.2, 0.6);
95  dHist_M2pidiff = new TH1I("M2pidiff", ";M_{#pi^{+}#pi^{-}} Kin - Gen (GeV/c^{2})", 400, -0.05, 0.05);
96  dHist_tgen = new TH1I("tgen", ";|t| Gen (GeV/c)^{2}", 100, 0.0, 0.01);
97  dHist_tkin = new TH1I("tkin", ";|t| Kin (GeV/c)^{2}", 100, 0.0, 0.01);
98  dHist_tdiff = new TH1I("tdiff", ";|t| Kin - Gen (GeV/c)^{2}", 100, -0.01, 0.01);
99  dHist_tkin_tgen = new TH2I("tkin_tgen", "; |t| Gen ; |t| Kin (GeV/c)^{2}", 50, 0, 0.002, 50, 0, 0.002);
100  dHist_CosTheta_psi = new TH2I("CosTheta_psi", "; #psi; Cos#Theta", 90, -180., 180, 200, -1., 1.);
101  dHist_CosThetakin_CosThetagen = new TH2I("CosThetakin_CosThetagen", "; Cos#Theta Gen; Cos#Theta Kin", 50, -1, 1, 50, -1., 1.);
102  dHist_phigen_Phigen = new TH2I("phigen_Phigen", "; #Phi Gen; #phi Gen", 90, -180., 180, 90,-180,180);
103  dHist_phikin_Phikin = new TH2I("phikin_Phikin", "; #Phi Kin; #phi Kin", 90, -180., 180, 90,-180,180);
104  dHist_phimeas_phigen = new TH2I("phimeas_phigen", "; #phi Gen ; #phi Meas", 90, -180., 180, 90,-180,180);
105  dHist_phikin_phigen = new TH2I("phikin_phigen", "; #phi Gen ; #phi Kin", 90, -180., 180, 90,-180,180);
106  dHist_Phikin_Phigen = new TH2I("Phikin_Phigen", "; #Phi Gen ; #Phi Kin", 90, -180., 180, 90,-180,180);
107  dHist_Phimeas_Phigen = new TH2I("Phimeas_Phigen", "; #Phi Gen ; #Phi Meas", 90, -180., 180, 90,-180,180);
108  dHist_Delta_phi = new TH2I("Delta_phi", "; #phi ; #Delta #phi", 90, -180., 180, 90,-180,180);
109  dHist_Delta_Phi = new TH2I("Delta_Phi", "; #Phi ; #Delta #Phi", 90, -180., 180, 90,-180,180);
110  dHist_Delta_phimeas = new TH2I("Delta_phimeas", "; #phi ; #Delta #phi Meas", 90, -180., 180, 90,-180,180);
111  dHist_Delta_Phimeas = new TH2I("Delta_Phimeas", "; #Phi ; #Delta #Phi Meas", 90, -180., 180, 90,-180,180);
112  dHist_CosTheta = new TH1I("CosTheta", "; Cos#Theta", 100, -1., 1.);
113  dHist_CosThetadiff = new TH1I("CosThetadiff", "; Cos#Theta diff Kin-Gen", 50, -0.5, 0.5);
114  dHist_Phigen = new TH1I("Phigen", ";Phigen (degrees)", 360,-180,180);
115  dHist_phigen = new TH1I("phigen", ";phigen (degrees)", 360,-180,180);
116  dHist_Phikin = new TH1I("Phikin", ";Phikin (degrees)", 360,-180,180);
117  dHist_phikin = new TH1I("phikin", ";phikin (degrees)", 360,-180,180);
118  dHist_Phimeas = new TH1I("Phimeas", ";Phimeas (degrees)", 360,-180,180);
119  dHist_phimeas = new TH1I("phimeas", ";phimeas (degrees)", 360,-180,180);
120  dHist_psikin = new TH1I("psikin", ";psi Kin (degrees)", 360,-180,180);
121  dHist_psigen = new TH1I("psigen", ";psi Gen (degrees)", 360,-180,180);
122  dHist_Phidiff = new TH1I("Phidiff", ";Phi Kin - Gen (degrees)", 100,-50,50);
123  dHist_phidiff = new TH1I("phidiff", ";phi Kin - Gen (degrees)", 100,-50,50);
124  dHist_psidiff = new TH1I("psidiff", ";psi Kin - Gen (degrees)", 100,-50,50);
125 
126  dHist_pipDeltap = new TH1I("pipDeltap","; #pi^{+}: Thrown p - KinFit p/Thrown p",100,-0.2,0.2);
127  dHist_pimDeltap = new TH1I("pimDeltap","; #pi^{-}: Thrown p - KinFit p/ Thrown p",100,-0.2,0.2);
128 
129  dHist_pipDeltap_Measured = new TH1I("pipDeltap_Measured","; #pi^{+}: Thrown p - Measured p/Thrown p",100,-0.2,0.2);
130  dHist_pimDeltap_Measured = new TH1I("pimDeltap_Measured","; #pi^{-}: Thrown p - Measured p/ Thrown p",100,-0.2,0.2);
131 
132 
133 
134  /************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - MAIN TREE *************************/
135 
136  //EXAMPLE MAIN TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)):
137  //The type for the branch must be included in the brackets
138  //1st function argument is the name of the branch
139  //2nd function argument is the name of the branch that contains the size of the array (for fundamentals only)
140  /*
141  dTreeInterface->Create_Branch_Fundamental<Int_t>("my_int"); //fundamental = char, int, float, double, etc.
142  dTreeInterface->Create_Branch_FundamentalArray<Int_t>("my_int_array", "my_int");
143  dTreeInterface->Create_Branch_FundamentalArray<Float_t>("my_combo_array", "NumCombos");
144  dTreeInterface->Create_Branch_NoSplitTObject<TLorentzVector>("my_p4");
145  dTreeInterface->Create_Branch_ClonesArray<TLorentzVector>("my_p4_array");
146  */
147 
148  /************************** EXAMPLE USER INITIALIZATION: CUSTOM OUTPUT BRANCHES - FLAT TREE *************************/
149 
150  //EXAMPLE FLAT TREE CUSTOM BRANCHES (OUTPUT ROOT FILE NAME MUST FIRST BE GIVEN!!!! (ABOVE: TOP)):
151  //The type for the branch must be included in the brackets
152  //1st function argument is the name of the branch
153  //2nd function argument is the name of the branch that contains the size of the array (for fundamentals only)
154  /*
155  dFlatTreeInterface->Create_Branch_Fundamental<Int_t>("flat_my_int"); //fundamental = char, int, float, double, etc.
156  dFlatTreeInterface->Create_Branch_FundamentalArray<Int_t>("flat_my_int_array", "flat_my_int");
157  dFlatTreeInterface->Create_Branch_NoSplitTObject<TLorentzVector>("flat_my_p4");
158  dFlatTreeInterface->Create_Branch_ClonesArray<TLorentzVector>("flat_my_p4_array");
159  */
160 
161  /************************************* ADVANCED EXAMPLE: CHOOSE BRANCHES TO READ ************************************/
162 
163  //TO SAVE PROCESSING TIME
164  //If you know you don't need all of the branches/data, but just a subset of it, you can speed things up
165  //By default, for each event, the data is retrieved for all branches
166  //If you know you only need data for some branches, you can skip grabbing data from the branches you don't need
167  //Do this by doing something similar to the commented code below
168 
169  //dTreeInterface->Clear_GetEntryBranches(); //now get none
170  //dTreeInterface->Register_GetEntryBranch("Proton__P4"); //manually set the branches you want
171 }
172 
173 Bool_t DSelector_Z2pi_trees::Process(Long64_t locEntry)
174 {
175  // The Process() function is called for each entry in the tree. The entry argument
176  // specifies which entry in the currently loaded tree is to be processed.
177  //
178  // This function should contain the "body" of the analysis. It can contain
179  // simple or elaborate selection criteria, run algorithms on the data
180  // of the event and typically fill histograms.
181  //
182  // The processing can be stopped by calling Abort().
183  // Use fStatus to set the return value of TTree::Process().
184  // The return value is currently not used.
185 
186  double PI = 3.14159;
187 
188  //CALL THIS FIRST
189  DSelector::Process(locEntry); //Gets the data from the tree for the entry
190  //cout << "RUN " << Get_RunNumber() << ", EVENT " << Get_EventNumber() << endl;
191  //TLorentzVector locProductionX4 = Get_X4_Production();
192 
193  /******************************************** GET POLARIZATION ORIENTATION ******************************************/
194 
195  //Only if the run number changes
196  //RCDB environment must be setup in order for this to work! (Will return false otherwise)
197  UInt_t locRunNumber = Get_RunNumber();
198  if(locRunNumber != dPreviousRunNumber)
199  {
200  dIsPolarizedFlag = dAnalysisUtilities.Get_IsPolarizedBeam(locRunNumber, dIsPARAFlag);
201  dPreviousRunNumber = locRunNumber;
202  }
203 
204  /********************************************* SETUP UNIQUENESS TRACKING ********************************************/
205 
206  //ANALYSIS ACTIONS: Reset uniqueness tracking for each action
207  //For any actions that you are executing manually, be sure to call Reset_NewEvent() on them here
208  Reset_Actions_NewEvent();
209 
210  //PREVENT-DOUBLE COUNTING WHEN HISTOGRAMMING
211  //Sometimes, some content is the exact same between one combo and the next
212  //e.g. maybe two combos have different beam particles, but the same data for the final-state
213  //When histogramming, you don't want to double-count when this happens: artificially inflates your signal (or background)
214  //So, for each quantity you histogram, keep track of what particles you used (for a given combo)
215  //Then for each combo, just compare to what you used before, and make sure it's unique
216 
217  //EXAMPLE 1: Particle-specific info:
218  set<Int_t> locUsedSoFar_BeamEnergy; //Int_t: Unique ID for beam particles. set: easy to use, fast to search
219  // set<Int_t> locUsedSoFar_Proton, locUsedSoFar_PiPlus, locUsedSoFar_PiMinus;
220  set<Int_t> locUsedSoFar_PiPlus, locUsedSoFar_PiMinus;
221 
222  //EXAMPLE 2: Combo-specific info:
223  //In general: Could have multiple particles with the same PID: Use a set of Int_t's
224  //In general: Multiple PIDs, so multiple sets: Contain within a map
225  //Multiple combos: Contain maps within a set (easier, faster to search)
226  set<map<Particle_t, set<Int_t> > > locUsedSoFar_MissingMass, locUsedSoFar_2pi, locUsedSoFar_Angles;
227 
228  //INSERT USER ANALYSIS UNIQUENESS TRACKING HERE
229 
230  /**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/
231 
232  /*
233  Int_t locMyInt = 7;
234  dTreeInterface->Fill_Fundamental<Int_t>("my_int", locMyInt);
235 
236  TLorentzVector locMyP4(4.0, 3.0, 2.0, 1.0);
237  dTreeInterface->Fill_TObject<TLorentzVector>("my_p4", locMyP4);
238 
239  for(int loc_i = 0; loc_i < locMyInt; ++loc_i)
240  dTreeInterface->Fill_Fundamental<Int_t>("my_int_array", 3*loc_i, loc_i); //2nd argument = value, 3rd = array index
241  */
242 
243 
244  /******************************************* LOOP OVER THROWN DATA (OPTIONAL) ***************************************/
245 
246 
247  //Thrown beam: just use directly
248  double locEbeam_Thrown = 0;
249  if(dThrownBeam != NULL)
250  locEbeam_Thrown = dThrownBeam->Get_P4().E();
251 
252  //Loop over throwns
253  TLorentzVector locPb208P4_Thrown;
254  TLorentzVector locPiPlusP4_Thrown;
255  TLorentzVector locPiMinusP4_Thrown;
256 
257  for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
258  {
259  //Set branch array indices corresponding to this particle
260  dThrownWrapper->Set_ArrayIndex(loc_i);
261  Particle_t thrown_pid = dThrownWrapper->Get_PID();
262  // cout << " loc_i=" << loc_i << " thrown_pid=" << thrown_pid << endl;
263  TLorentzVector locP4_Thrown = dThrownWrapper->Get_P4();
264  if (loc_i == 2) locPb208P4_Thrown = locP4_Thrown; // assume order of particles as PID is zero at the moment
265  if (loc_i == 0) locPiPlusP4_Thrown = locP4_Thrown;
266  if (loc_i == 1) locPiMinusP4_Thrown = locP4_Thrown;
267 
268  }
269  cout << endl << "Thrown" << endl;
270  cout << " dThrownBeam="; dThrownBeam->Get_P4().Print();
271  cout << " locPb208P4="; locPb208P4_Thrown.Print();
272  cout << " locPiPlusP4="; locPiPlusP4_Thrown.Print();
273  cout << " locPiMinusP4="; locPiMinusP4_Thrown.Print();
274  TLorentzVector loc2piP4_Thrown = locPiPlusP4_Thrown + locPiMinusP4_Thrown;
275  double tgen = (dThrownBeam->Get_P4() - locPiPlusP4_Thrown - locPiMinusP4_Thrown).M2(); // use beam and 2pi momenta
276 
277 
278 
279 
280  /************************************************* LOOP OVER COMBOS *************************************************/
281 
282  //Loop over combos
283  for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i)
284  {
285  //Set branch array indices for combo and all combo particles
286  dComboWrapper->Set_ComboIndex(loc_i);
287 
288  // Is used to indicate when combos have been cut
289  if(dComboWrapper->Get_IsComboCut()) // Is false when tree originally created
290  continue; // Combo has been cut previously
291 
292  /********************************************** GET PARTICLE INDICES *********************************************/
293 
294  //Used for tracking uniqueness when filling histograms, and for determining unused particles
295 
296  //Step 0
297  Int_t locBeamID = dComboBeamWrapper->Get_BeamID();
298  Int_t locPiPlusTrackID = dPiPlusWrapper->Get_TrackID();
299  Int_t locPiMinusTrackID = dPiMinusWrapper->Get_TrackID();
300 
301  /*********************************************** GET FOUR-MOMENTUM **********************************************/
302 
303  // Get P4's: //is kinfit if kinfit performed, else is measured
304  //dTargetP4 is target p4
305  //Step 0
306  TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4();
307  TLorentzVector locMissingPb208P4 = dMissingPb208Wrapper->Get_P4();
308  TLorentzVector locPiPlusP4 = dPiPlusWrapper->Get_P4();
309  TLorentzVector locPiMinusP4 = dPiMinusWrapper->Get_P4();
310 
311  TLorentzVector locMissingP4 = locBeamP4; // Ignore target mass/recoil
312  locMissingP4 -= locPiPlusP4 + locPiMinusP4;
313  TLorentzVector loc2piP4 = locPiPlusP4 + locPiMinusP4;
314 
315  cout << "Kin Fit" << endl;
316  cout << " locBeamP4="; locBeamP4.Print();
317  cout << " locMissingPb208P4="; locMissingPb208P4.Print();
318  cout << " locPiPlusP4="; locPiPlusP4.Print();
319  cout << " locPiMinusP4="; locPiMinusP4.Print();
320  cout << " locMissingP4="; locMissingP4.Print();
321  cout << " loc2piP4="; loc2piP4.Print();
322 
323 
324  cout << "KIN: MM2 =" << locMissingP4.M2() << " DeltaE=" << locBeamP4.E()-loc2piP4.E() << " GenDeltaE=" << locEbeam_Thrown-locBeamP4.E() << endl;
325 
326  // Get Measured P4's:
327  //Step 0
328  TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured();
329  TLorentzVector locMissingPb208P4_Measured (0,0,0,193.750748);
330  TLorentzVector locMissingPb208P4_Measured_input = dMissingPb208Wrapper->Get_P4_Measured();
331  TLorentzVector locPiPlusP4_Measured = dPiPlusWrapper->Get_P4_Measured();
332  TLorentzVector locPiMinusP4_Measured = dPiMinusWrapper->Get_P4_Measured();
333 
334  TLorentzVector locMissingP4_Measured = locBeamP4_Measured; // Ignore target mass/recoil
335  locMissingP4_Measured -= locPiPlusP4_Measured + locPiMinusP4_Measured;
336  locMissingPb208P4_Measured += locMissingP4_Measured;
337  TLorentzVector loc2piP4_Measured = locPiPlusP4_Measured + locPiMinusP4_Measured;
338 
339  cout << "Measured" << endl;
340  cout << " locBeamP4_Measured="; locBeamP4_Measured.Print();
341  cout << " locMissingPb208P4_Measured="; locMissingPb208P4_Measured.Print();
342  cout << " locPiPlusP4_Measured="; locPiPlusP4_Measured.Print();
343  cout << " locPiMinusP4_Measured="; locPiMinusP4_Measured.Print();
344  cout << " locMissingP4_Measured="; locMissingP4_Measured.Print();
345  cout << " loc2piP4_Measured="; loc2piP4_Measured.Print();
346 
347  /********************************************* COMBINE FOUR-MOMENTUM ********************************************/
348 
349  // DO YOUR STUFF HERE
350 
351  // Combine 4-vectors
352 
353  /******************************************** EXECUTE ANALYSIS ACTIONS *******************************************/
354 
355  // Loop through the analysis actions, executing them in order for the active particle combo
356  if(!Execute_Actions()) //if the active combo fails a cut, IsComboCutFlag automatically set
357  continue;
358 
359  //if you manually execute any actions, and it fails a cut, be sure to call:
360  //dComboWrapper->Set_IsComboCut(true);
361 
362  /**************************************** EXAMPLE: FILL CUSTOM OUTPUT BRANCHES **************************************/
363 
364  /*
365  TLorentzVector locMyComboP4(8.0, 7.0, 6.0, 5.0);
366  //for arrays below: 2nd argument is value, 3rd is array index
367  //NOTE: By filling here, AFTER the cuts above, some indices won't be updated (and will be whatever they were from the last event)
368  //So, when you draw the branch, be sure to cut on "IsComboCut" to avoid these.
369  dTreeInterface->Fill_Fundamental<Float_t>("my_combo_array", -2*loc_i, loc_i);
370  dTreeInterface->Fill_TObject<TLorentzVector>("my_p4_array", locMyComboP4, loc_i);
371  */
372 
373  /**************************************** EXAMPLE: PID CUT ACTION ************************************************/
374 
375 
376  // Proton CDC dE/dx histogram and cut
377  /*double locProton_dEdx_CDC = dProtonWrapper->Get_dEdx_CDC()*1e6;
378  if(locProton_dEdx_CDC < fMinProton_dEdx->Eval(locProtonP4.P())) {
379  dComboWrapper->Set_IsComboCut(true);
380  continue;
381  }*/
382 
383  // Pi+/- CDC dE/dx histogram cut (histograms in HistComboPID action)
384  double locPiPlus_dEdx_CDC = dPiPlusWrapper->Get_dEdx_CDC()*1e6;
385  double locPiMinus_dEdx_CDC = dPiMinusWrapper->Get_dEdx_CDC()*1e6;
386  if(locPiPlus_dEdx_CDC > fMaxPion_dEdx->Eval(locPiPlusP4.P()) || locPiMinus_dEdx_CDC > fMaxPion_dEdx->Eval(locPiMinusP4.P())) {
387  dComboWrapper->Set_IsComboCut(true);
388  continue;
389  }
390 
391  cout << " Passed CDC dE/dX cut " << endl;
392 
393 
394 
395  /************************************ EXAMPLE: SELECTION CUTS AND HISTOGRAMS ************************************/
396 
397  //Missing Mass Squared
398  double locMissingMassSquared = locMissingP4_Measured.M2();
399 
400  //Uniqueness tracking: Build the map of particles used for the missing mass
401  //For beam: Don't want to group with final-state photons. Instead use "Unknown" PID (not ideal, but it's easy).
402  map<Particle_t, set<Int_t> > locUsedThisCombo_MissingMass;
403  locUsedThisCombo_MissingMass[Unknown].insert(locBeamID); //beam
404  locUsedThisCombo_MissingMass[PiPlus].insert(locPiPlusTrackID);
405  locUsedThisCombo_MissingMass[PiMinus].insert(locPiMinusTrackID);
406 
407  // beam energy cut for SDME
408  if(locBeamP4.E() < dMinBeamEnergy || locBeamP4.E() > dMaxBeamEnergy) {
409  dComboWrapper->Set_IsComboCut(true);
410  continue;
411  }
412 
413  cout << " Passed beam energy cut " << endl;
414 
415  // Cut on Missing mass
416  if((locMissingMassSquared < dMinMissingMassSquared ) || (locMissingMassSquared > dMaxMissingMassSquared)){ // measured
417  // if((locMissingP4.M2() < dMinMissingMassSquared ) || (locMissingP4.M2() > dMaxMissingMassSquared)){ // kinfit
418  dComboWrapper->Set_IsComboCut(true);
419  continue;
420  }
421 
422  cout << " Passed Missing mass cut " << endl;
423 
424 
425  // kinematic fit CL cut
426  dHist_KinFitChiSq->Fill(dComboWrapper->Get_ChiSq_KinFit()/dComboWrapper->Get_NDF_KinFit());
427  dHist_KinFitCL->Fill(dComboWrapper->Get_ConfidenceLevel_KinFit());
428  if(dComboWrapper->Get_ConfidenceLevel_KinFit() <= dMinKinFitCL) {
429  dComboWrapper->Set_IsComboCut(true);
430  continue;
431  }
432  cout << " Passed kinematic Fit CL cut " << endl;
433 
434  // 2pi mass histogram and cut
435  map<Particle_t, set<Int_t> > locUsedThisCombo_2piMass;
436  locUsedThisCombo_2piMass[PiPlus].insert(locPiPlusTrackID);
437  locUsedThisCombo_2piMass[PiMinus].insert(locPiMinusTrackID);
438  if(loc2piP4.M() < dMin2piMass || loc2piP4.M() > dMax2piMass) {
439  dComboWrapper->Set_IsComboCut(true);
440  continue;
441  }
442  cout << " Passed 2pi mass cut " << endl;
443 
444 
445  if(locUsedSoFar_2pi.find(locUsedThisCombo_2piMass) == locUsedSoFar_2pi.end())
446  {
447  dHist_M2pikin->Fill(loc2piP4.M());
448  dHist_M2pigen->Fill(loc2piP4_Thrown.M());
449  dHist_M2pidiff->Fill(loc2piP4.M()-loc2piP4_Thrown.M());
450  locUsedSoFar_2pi.insert(locUsedThisCombo_2piMass);
451  }
452 
453 
454 
455  /**************************************** EXAMPLE: HISTOGRAM energies and momenta *****************************************/
456 
457  //Histogram beam energy (if haven't already)
458  if(locUsedSoFar_BeamEnergy.find(locBeamID) == locUsedSoFar_BeamEnergy.end())
459  {
460  dHist_BeamEnergy->Fill(locBeamP4.E());
461  locUsedSoFar_BeamEnergy.insert(locBeamID);
462  }
463 
464  //compare to what's been used so far
465  if(locUsedSoFar_MissingMass.find(locUsedThisCombo_MissingMass) == locUsedSoFar_MissingMass.end())
466  {
467  //unique missing mass combo: histogram it, and register this combo of particles
468  dHist_MissingMassSquared->Fill(locMissingMassSquared);
469  locUsedSoFar_MissingMass.insert(locUsedThisCombo_MissingMass);
470  }
471 
472 
473  /*if(locUsedSoFar_Proton.find(locProtonTrackID) == locUsedSoFar_Proton.end())
474  {
475  dHist_pMomentumMeasured->Fill(locProtonP4.Vect().Mag());
476  dHist_Proton_dEdx_P->Fill(locProtonP4.P(), locProton_dEdx_CDC);
477  dHist_pDeltap->Fill(locProtonP4_Thrown.Vect().Mag() > 0? (locProtonP4_Thrown.Vect().Mag()-locProtonP4.Vect().Mag())/locProtonP4_Thrown.Vect().Mag() : 0);
478  dHist_pDeltap_Measured->Fill(locProtonP4_Thrown.Vect().Mag() > 0? (locProtonP4_Thrown.Vect().Mag()-locProtonP4_Measured.Vect().Mag())/locProtonP4_Thrown.Vect().Mag() : 0);
479  locUsedSoFar_Proton.insert(locProtonTrackID);
480  }*/
481 
482 
483  if(locUsedSoFar_PiPlus.find(locPiPlusTrackID) == locUsedSoFar_PiPlus.end())
484  {
485  dHist_piplusMomentumMeasured->Fill(locPiPlusP4.Vect().Mag());
486  dHist_pipDeltap->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
487  dHist_pipDeltap_Measured->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4_Measured.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
488  locUsedSoFar_PiPlus.insert(locPiPlusTrackID);
489  }
490 
491  if(locUsedSoFar_PiMinus.find(locPiMinusTrackID) == locUsedSoFar_PiMinus.end())
492  {
493  dHist_piminusMomentumMeasured->Fill(locPiMinusP4.Vect().Mag());
494  dHist_pimDeltap->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
495  dHist_pimDeltap_Measured->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4_Measured.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
496  locUsedSoFar_PiMinus.insert(locPiMinusTrackID);
497  }
498 
499  cout << "Ebeam Thrown P4.E=" << dThrownBeam->Get_P4().E() << " Kinfit P4.E=" << locBeamP4.E() << " P4 Measured =" << locBeamP4_Measured.E() << endl;
500  // cout << "Proton Thrown P4.E=" << locPb208P4_Thrown.E() << " Kinfit P4.E=" << locProtonP4.E() << " P4 Measured =" << locProtonP4_Measured.E() << " CL=" << dComboWrapper->Get_ConfidenceLevel_KinFit() << endl;
501  cout << "PiPlus Thrown P4.E=" << locPiPlusP4_Thrown.E() << " Kinfit P4.E=" << locPiPlusP4.E() << " P4 Measured =" << locPiPlusP4_Measured.E() << endl;
502  cout << "PiMinus Thrown P4.E=" << locPiMinusP4_Thrown.E() << " Kinfit P4.E=" << locPiMinusP4.E() << " P4 Measured =" << locPiMinusP4_Measured.E() << endl << endl;
503 
504 
505  if ( locBeamP4_Measured.Z() < 0 || locPiPlusP4_Measured.Z() < 0 || locPiMinusP4_Measured.Z() < 0 ||
506  locPiPlusP4.Z() < 0 || locPiMinusP4.Z() < 0 || (locPiPlusP4_Measured.Z() + locPiMinusP4_Measured.Z() + locMissingPb208P4_Measured_input.Z()) < 0 ||
507  (locPiPlusP4.Z() + locPiMinusP4.Z() + locMissingPb208P4.Z()) < 0 ) {
508  dComboWrapper->Set_IsComboCut(true);
509  cout << "*** Failed Negative pz cut ***" << endl;
510  continue;
511  }
512  cout << " Passed Negative pz cut " << endl;
513 
514 
515  // Thrown (generated) variables
516  // calculate kinematic and angular variables
517  double tgen = (dThrownBeam->Get_P4() - loc2piP4_Thrown).M2(); // use beam and 2pi momenta
518  TLorentzRotation resonanceBoost( -loc2piP4_Thrown.BoostVector() ); // boost into 2pi frame
519  TLorentzVector beam_res = resonanceBoost * dThrownBeam->Get_P4();
520  TLorentzVector recoil_res = resonanceBoost * locPb208P4_Thrown;
521  TLorentzVector p1_res = resonanceBoost * locPiPlusP4_Thrown;
522  TLorentzVector p2_res = resonanceBoost * locPiMinusP4_Thrown;
523 
524  double phipol = 0; // *** Note assumes horizontal polarization plane.
525  TVector3 eps(cos(phipol), sin(phipol), 0.0); // beam polarization vector in lab
526 
527  // choose helicity frame: z-axis opposite recoil target in rho rest frame. Note that for Primakoff recoil is never measured.
528  TVector3 y = (dThrownBeam->Get_P4().Vect().Unit().Cross(-locPb208P4_Thrown.Vect().Unit())).Unit();
529 
530  // choose helicity frame: z-axis opposite recoil proton in rho rest frame
531  TVector3 z = -1. * recoil_res.Vect().Unit();
532  TVector3 x = y.Cross(z).Unit();
533  TVector3 anglesgen( (p1_res.Vect()).Dot(x),
534  (p1_res.Vect()).Dot(y),
535  (p1_res.Vect()).Dot(z) );
536 
537  double CosThetagen = anglesgen.CosTheta();
538  double phigen = anglesgen.Phi();
539 
540  double Phigen = atan2(y.Dot(eps), dThrownBeam->Get_P4().Vect().Unit().Dot(eps.Cross(y)));
541 
542  double psigen = Phigen - phigen;
543  if(psigen < -3.14159) psigen += 2*3.14159;
544  if(psigen > 3.14159) psigen -= 2*3.14159;
545 
546 
547  // Repeat for kinematically fit variables.
548  // calculate kinematic and angular variables
549  double tkin = (locBeamP4 - loc2piP4).M2(); // use beam and 2pi momenta
550  TLorentzRotation resonanceBoost2( -loc2piP4.BoostVector() ); // boost into 2pi frame
551  beam_res = resonanceBoost2 * locBeamP4;
552  recoil_res = resonanceBoost2 * locMissingPb208P4;
553  p1_res = resonanceBoost2 * locPiPlusP4;
554  p2_res = resonanceBoost2 * locPiMinusP4;
555 
556  // choose helicity frame: z-axis opposite recoil target in rho rest frame. Note that for Primakoff recoil is missing P4, including target.
557  y = (locBeamP4.Vect().Unit().Cross(-locMissingPb208P4.Vect().Unit())).Unit();
558 
559  // choose helicity frame: z-axis opposite recoil proton in rho rest frame
560  z = -1. * recoil_res.Vect().Unit();
561  x = y.Cross(z).Unit();
562  TVector3 angleskin( (p1_res.Vect()).Dot(x),
563  (p1_res.Vect()).Dot(y),
564  (p1_res.Vect()).Dot(z) );
565 
566  double CosThetakin = angleskin.CosTheta();
567  double phikin = angleskin.Phi();
568 
569  double Phikin = atan2(y.Dot(eps), locBeamP4.Vect().Unit().Dot(eps.Cross(y)));
570 
571  double psikin = Phikin - phikin;
572  if(psikin < -3.14159) psikin += 2*3.14159;
573  if(psikin > 3.14159) psikin -= 2*3.14159;
574 
575  // Repeat for measured variables.
576  // calculate measured and angular variables
577  double tmeas = (locBeamP4_Measured - loc2piP4_Measured).M2(); // use beam and 2pi momenta
578  TLorentzRotation resonanceBoost3( -loc2piP4_Measured.BoostVector() ); // boost into 2pi frame
579  beam_res = resonanceBoost3 * locBeamP4_Measured;
580  recoil_res = resonanceBoost3 * locMissingPb208P4_Measured;
581  p1_res = resonanceBoost3 * locPiPlusP4_Measured;
582  p2_res = resonanceBoost3 * locPiMinusP4_Measured;
583 
584  // choose helicity frame: z-axis opposite recoil target in rho rest frame. Note that for Primakoff recoil is missing P4, including target.
585  y = (locBeamP4_Measured.Vect().Unit().Cross(-locMissingPb208P4_Measured.Vect().Unit())).Unit();
586 
587  // choose helicity frame: z-axis opposite recoil proton in rho rest frame
588  z = -1. * recoil_res.Vect().Unit();
589  x = y.Cross(z).Unit();
590  TVector3 anglesmeas( (p1_res.Vect()).Dot(x),
591  (p1_res.Vect()).Dot(y),
592  (p1_res.Vect()).Dot(z) );
593 
594  double CosThetameas = anglesmeas.CosTheta();
595  double phimeas = anglesmeas.Phi();
596 
597  double Phimeas = atan2(y.Dot(eps), locBeamP4_Measured.Vect().Unit().Dot(eps.Cross(y)));
598 
599  double psimeas = Phimeas - phimeas;
600  if(psimeas < -3.14159) psimeas += 2*3.14159;
601  if(psimeas > 3.14159) psimeas -= 2*3.14159;
602 
603  double Delta_phi = phikin - phigen;
604  double Delta_Phi = Phikin - Phigen;
605  double Delta_phimeas = phimeas - phigen;
606  double Delta_Phimeas = Phimeas - Phigen;
607 
608  map<Particle_t, set<Int_t> > locUsedThisCombo_Angles;
609  locUsedThisCombo_Angles[Unknown].insert(locBeamID); //beam
610  // locUsedThisCombo_Angles[Proton].insert(locProtonTrackID);
611  locUsedThisCombo_Angles[PiPlus].insert(locPiPlusTrackID);
612  locUsedThisCombo_Angles[PiMinus].insert(locPiMinusTrackID);
613  if(locUsedSoFar_Angles.find(locUsedThisCombo_Angles) == locUsedSoFar_Angles.end())
614  {
615  dHist_tgen->Fill(fabs(tgen));
616  dHist_tkin->Fill(fabs(tkin));
617  dHist_tdiff->Fill(fabs(tkin)-fabs(tgen));
618  dHist_tkin_tgen->Fill(fabs(tgen),fabs(tkin));
619  dHist_CosTheta_psi->Fill(psikin*180./3.14159, CosThetakin);
620  dHist_CosTheta->Fill(CosThetakin);
621  dHist_CosThetadiff->Fill(CosThetakin-CosThetagen);
622  // dHist_phi->Fill(phi*180./3.14159);
623  dHist_CosThetakin_CosThetagen->Fill(CosThetagen, CosThetakin);
624  dHist_phikin_Phikin->Fill(Phikin*180./3.14159,phikin*180./3.14159);
625  dHist_phigen_Phigen->Fill(Phigen*180./3.14159,phigen*180./3.14159);
626  dHist_phimeas_phigen->Fill(phigen*180./3.14159,phimeas*180./3.14159);
627  dHist_phikin_phigen->Fill(phigen*180./3.14159,phikin*180./3.14159);
628  dHist_Phimeas_Phigen->Fill(Phigen*180./3.14159,Phimeas*180./3.14159);
629  dHist_Phikin_Phigen->Fill(Phigen*180./3.14159,Phikin*180./3.14159);
630  dHist_Delta_phi->Fill(phigen*180./3.14159,Delta_phi*180./3.14159);
631  dHist_Delta_Phi->Fill(Phigen*180./3.14159,Delta_Phi*180./3.14159);
632  dHist_Delta_phimeas->Fill(phigen*180./3.14159,Delta_phimeas*180./3.14159);
633  dHist_Delta_Phimeas->Fill(Phigen*180./3.14159,Delta_Phimeas*180./3.14159);
634  dHist_phigen->Fill(phigen*180./3.14159);
635  dHist_Phigen->Fill(Phigen*180./3.14159);
636  dHist_phikin->Fill(phikin*180./3.14159);
637  dHist_Phikin->Fill(Phikin*180./3.14159);
638  dHist_Phimeas->Fill(Phimeas*180./3.14159);
639  dHist_phimeas->Fill(phimeas*180./3.14159);
640  dHist_Phidiff->Fill((Phikin-Phigen)*180./3.14159);
641  dHist_phidiff->Fill((phikin-phigen)*180./3.14159);
642  dHist_psigen->Fill(psigen*180./3.14159);
643  dHist_psikin->Fill(psikin*180./3.14159);
644  dHist_psidiff->Fill((psikin-psigen)*180./3.14159);
645  // dHist_psikin->Fill(psikin*180./3.14159);
646  locUsedSoFar_Angles.insert(locUsedThisCombo_Angles);
647  }
648 
649  /****************************************** FILL FLAT TREE (IF DESIRED) ******************************************/
650 
651 
652  /*
653  //FILL ANY CUSTOM BRANCHES FIRST!!
654  Int_t locMyInt_Flat = 7;
655  dFlatTreeInterface->Fill_Fundamental<Int_t>("flat_my_int", locMyInt_Flat);
656 
657  TLorentzVector locMyP4_Flat(4.0, 3.0, 2.0, 1.0);
658  dFlatTreeInterface->Fill_TObject<TLorentzVector>("flat_my_p4", locMyP4_Flat);
659 
660  for(int loc_j = 0; loc_j < locMyInt_Flat; ++loc_j)
661  {
662  dFlatTreeInterface->Fill_Fundamental<Int_t>("flat_my_int_array", 3*loc_j, loc_j); //2nd argument = value, 3rd = array index
663  TLorentzVector locMyComboP4_Flat(8.0, 7.0, 6.0, 5.0);
664  dFlatTreeInterface->Fill_TObject<TLorentzVector>("flat_my_p4_array", locMyComboP4_Flat, loc_j);
665  }
666  */
667 
668  //FILL FLAT TREE
669  //Fill_FlatTree(); //for the active combo
670  } // end of combo loop
671 
672  //FILL HISTOGRAMS: Num combos / events surviving actions
673  Fill_NumCombosSurvivedHists();
674 
675  /******************************************* LOOP OVER THROWN DATA (OPTIONAL) ***************************************/
676 /*
677  //Thrown beam: just use directly
678  if(dThrownBeam != NULL)
679  double locEnergy = dThrownBeam->Get_P4().E();
680 
681  //Loop over throwns
682  for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
683  {
684  //Set branch array indices corresponding to this particle
685  dThrownWrapper->Set_ArrayIndex(loc_i);
686 
687  //Do stuff with the wrapper here ...
688  }
689 */
690  /****************************************** LOOP OVER OTHER ARRAYS (OPTIONAL) ***************************************/
691 /*
692  //Loop over beam particles (note, only those appearing in combos are present)
693  for(UInt_t loc_i = 0; loc_i < Get_NumBeam(); ++loc_i)
694  {
695  //Set branch array indices corresponding to this particle
696  dBeamWrapper->Set_ArrayIndex(loc_i);
697 
698  //Do stuff with the wrapper here ...
699  }
700 
701  //Loop over charged track hypotheses
702  for(UInt_t loc_i = 0; loc_i < Get_NumChargedHypos(); ++loc_i)
703  {
704  //Set branch array indices corresponding to this particle
705  dChargedHypoWrapper->Set_ArrayIndex(loc_i);
706 
707  //Do stuff with the wrapper here ...
708  }
709 
710  //Loop over neutral particle hypotheses
711  for(UInt_t loc_i = 0; loc_i < Get_NumNeutralHypos(); ++loc_i)
712  {
713  //Set branch array indices corresponding to this particle
714  dNeutralHypoWrapper->Set_ArrayIndex(loc_i);
715 
716  //Do stuff with the wrapper here ...
717  }
718 */
719 
720  /************************************ EXAMPLE: FILL CLONE OF TTREE HERE WITH CUTS APPLIED ************************************/
721 
722  Bool_t locIsEventCut = true;
723  for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) {
724  //Set branch array indices for combo and all combo particles
725  dComboWrapper->Set_ComboIndex(loc_i);
726  // Is used to indicate when combos have been cut
727  if(dComboWrapper->Get_IsComboCut())
728  continue;
729  locIsEventCut = false; // At least one combo succeeded
730  break;
731  }
732  if(!locIsEventCut && dOutputTreeFileName != "")
733  Fill_OutputTree();
734 
735 
736  return kTRUE;
737 }
738 
740 {
741  //Save anything to output here that you do not want to be in the default DSelector output ROOT file.
742 
743  //Otherwise, don't do anything else (especially if you are using PROOF).
744  //If you are using PROOF, this function is called on each thread,
745  //so anything you do will not have the combined information from the various threads.
746  //Besides, it is best-practice to do post-processing (e.g. fitting) separately, in case there is a problem.
747 
748  //DO YOUR STUFF HERE
749 
750  //CALL THIS LAST
751  DSelector::Finalize(); //Saves results to the output file
752 }
void Process(unsigned int &NEvents, unsigned int &NEvents_read)
DKinematicData * dMissingPb208Wrapper
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define y
DBeamParticle * dComboBeamWrapper
#define PI
Definition: kinematics.h:18
DChargedTrackHypothesis * dPiMinusWrapper
Bool_t Process(Long64_t entry)
void Init(TTree *tree)
DChargedTrackHypothesis * dPiPlusWrapper
double sin(double)
Particle_t
Definition: particleType.h:12