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