Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventProcessor_BCAL_Shower.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DEventProcessor_BCAL_Shower.cc
4 // Created: Fri Oct 10 16:41:18 EDT 2014
5 // Creator: wmcginle (on Linux ifarm1101 2.6.32-220.7.1.el6.x86_64 x86_64)
6 //
7 
8 #include "DEventProcessor_BCAL_Shower.h"
9 
10 #include <TLorentzVector.h>
11 #include "TMath.h"
12 
13 #include "DANA/DApplication.h"
14 #include "BCAL/DBCALShower.h"
15 #include "BCAL/DBCALTruthShower.h"
16 #include "BCAL/DBCALCluster.h"
17 #include "BCAL/DBCALPoint.h"
18 #include "BCAL/DBCALHit.h"
19 #include "FCAL/DFCALShower.h"
20 #include "TRACKING/DMCThrown.h"
21 #include "ANALYSIS/DAnalysisUtilities.h"
22 //#include "TRACKING/DTrackFinder.h"
23 
24 #include <vector>
25 #include <deque>
26 #include <string>
27 #include <iostream>
28 
29 // Routine used to create our DEventProcessor
30 
31 extern "C"
32 {
33  void InitPlugin(JApplication *locApplication)
34  {
35  InitJANAPlugin(locApplication);
36  locApplication->AddProcessor(new DEventProcessor_BCAL_Shower()); //register this plugin
37  }
38 } // "C"
39 
40 #define FCAL_Z_OFFSET 640.0-65.0
41 //------------------
42 // init
43 //------------------
44 jerror_t DEventProcessor_BCAL_Shower::init(void)
45 {
46  TDirectory *dir = new TDirectoryFile("BCAL","BCAL");
47  dir->cd();
48 
49  two_gamma_mass = new TH1F("two_gamma_mass","two_gamma_mass",2000, 0.0, 3.00);
50  bcal_diphoton_mass = new TH1F("bcal_diphoton_mass","bcal_diphoton_mass",2000,0.0,3.0);
51  bcal_diphoton_mass_highE = new TH1F("bcal_diphoton_mass_highE","bcal_diphoton_mass_highE",2000,0.0,3.0);
52  //two_gamma_mass_corr = new TH1F("two_gamma_mass_corr","two_gamma_mass_corr",100, 0.0, 1.00);
53  //two_gamma_mass_cut = new TH1F("two_gamma_mass_cut","two_gamma_mass_cut",100, 0.0, 0.300);
54  bcal_fcal_two_gamma_mass = new TH1F("one_fcal_one_bcal_gamma_mass","bcal_fcal_two_gamma_mass",1000, 0.0, 1.00);
55  //bcal_fcal_two_gamma_mass_cut = new TH1F("fcal_two_gamma_mass_cut","bcal_fcal_two_gamma_mass_cut",100, 0.0, 0.300);
56  two_fcal_gamma_mass = new TH1F("two_fcal_gamma_mass","two_fcal_gamma_mass",1000,0.0,1.0);
57  two_fcal_gamma_mass_test = new TH1F("two_fcal_gamma_mass_test","two_fcal_gamma_mass_test",1000,0.0,1.0);
58  mass_summed = new TH1F("inv_mass","inv_mass",100,0.0,1.0);
59  //xy_shower = new TH2F("xy_shower","xy_shower",100, -100.0, 100., 100 , -100.0, 100.0);
60  z_shower = new TH1F("z_shower","z_shower",1600, 0.0, 6.0);
61  E_shower = new TH1F("E_shower","E_shower", 1600, 0.0, 6.0);
62  Neutral_E = new TH1F("Neutral_E","Neutral_E",1000,0.0, 400.0);
63  VertexZ = new TH1F("VertexZ","VertexZ",1000,-200,300);
64  goodVertexZ = new TH1F("goodVertexZ","goodVertexZ",1000,-200,300);
65  Theta_Hist_Both_BCAL = new TH1F("Theta_Both_BCAL","Theta_Both_BCAL",300,0.0,150.0);
66  Phi_Hist_Both_BCAL = new TH1F("Phi Both BCAL","Phi Both BCAL",720,-180.0,180.0);
67  Psi_Hist_Both_BCAL = new TH1F("Psi Both BCAL","Psi Both BCAL",1000,0.0,150.0);
68  Theta_Hist_Split_Gammas = new TH1F("Theta Split Gammas","Theta Split Gammas",300,0.0,150.0);
69  Phi_Hist_Split_Gammas = new TH1F("Phi Split Gammas","Phi Split Gammas",720,-180.0,180.0);
70  Psi_Hist_Split_Gammas = new TH1F("Psi Split Gammas","Psi Split Gammas",1000,0.0,150.0);
71  Theta_Hist_Both_FCAL = new TH1F("Theta Both FCAL","Theta Both FCAL",300,0.0,150.0);
72  Phi_Hist_Both_FCAL = new TH1F("Phi Both FCAL","Phi Both FCAL",720,-180.0,180.0);
73  Psi_Hist_Both_FCAL = new TH1F("Psi Both FCAL","Psi Both FCAL",1000,0.0,150.0);
74  E1_Vs_E2CosTheta = new TH2F("E1 Vs E2(1-cosTheta)","E1 Vs E2(1-cosTheta)",2000,0.0,1.5,2000.0,0.0,5.0);
75  Thrown_Gamma_Theta = new TH1F("Thrown Gamma Theta","Thrown Gamma Theta",300,0.0,150.0);
76  Thrown_Inv_Mass = new TH1F("Thrown Gamma Inv Mass","Thrown Gamma Inv Mass",1000,0.0,1.50);
77  Thrown_Vertex = new TH1F("Thrown Gamma Vertex","Thrown Gamma Vertex",2000,-300.0,400.0);
78  Thrown_Vertex_Frequency = new TH1F("Vertex Counting","Vertex Counting",100,0.0,10.0);
79  Pi0_velocity = new TH1F("v","v",1000,0.0,2.0);
80  Cluster_Energy = new TH1F("Cluster E","Cluster E",1000,0.0,2.0);
81  Point_Energy = new TH1F("Point Energy","Point Energy",1000,0.0,2.0);
82  Point_Module = new TH1F("Module","Module",48,0,48.0);
83  Point_z = new TH1F("Point z","Point z",1000,-100.0,600.0);
84  Assoc_E = new TH1F("Assoc E","Assoc E",1000,0,2.0);
85  Layer1_Energy_vs_Channel = new TH2F("Layer1 Energy vs channel","Layer1 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
86  Layer2_Energy_vs_Channel = new TH2F("Layer2 Energy vs channel","Layer2 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
87  Layer3_Energy_vs_Channel = new TH2F("Layer3 Energy vs channel","Layer3 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
88  Layer4_Energy_vs_Channel = new TH2F("Layer4 Energy vs channel","Layer4 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
89  Layer1_Energy = new TH1F("Layer1 Energy","Layer1 Energy",1500,0.0,1.5);
90  Layer2_Energy = new TH1F("Layer2 Energy","Layer2 Energy",1500,0.0,1.5);
91  Layer3_Energy = new TH1F("Layer3 Energy","Layer3 Energy",1500,0.0,1.5);
92  Layer4_Energy = new TH1F("Layer4 Energy","Layer4 Energy",1500,0.0,1.5);
93  Layer1_Energy_v2 = new TH1F("Layer1 Energy v2","Layer1 Energy v2",1500,0.0,1.5);
94  Layer2_Energy_v2 = new TH1F("Layer2 Energy v2","Layer2 Energy v2",1500,0.0,1.5);
95  Layer3_Energy_v2 = new TH1F("Layer3 Energy v2","Layer3 Energy v2",1500,0.0,1.5);
96  Layer4_Energy_v2 = new TH1F("Layer4 Energy v2","Layer4 Energy v2",1500,0.0,1.5);
97  Layer1_Energy_pos = new TH1F("Layer1 Energy pos","Layer1 Energy pos",1500,0.0,0.1);
98  Layer2_Energy_pos = new TH1F("Layer2 Energy pos","Layer2 Energy pos",1500,0.0,0.1);
99  Layer3_Energy_pos = new TH1F("Layer3 Energy pos","Layer3 Energy pos",1500,0.0,0.15);
100  Layer4_Energy_pos = new TH1F("Layer4 Energy pos","Layer4 Energy pos",1500,0.0,0.15);
101  Point_E_M1S1L1 = new TH1F("M1S1L1 E","M1S1L1 E",1500,0.0,1.0);
102  Point_E_M12S2L2 = new TH1F("M12S2L2 E","M12S2L2 E",1500,0.0,1.0);
103  Point_E_M25S3L3 = new TH1F("M25S3L3 E","M25S3L3 E",1500,0.0,1.0);
104  Point_E_M37S4L4 = new TH1F("M37S4L4 E","M37S4L4 E",1500,0.0,1.0);
105  Time_Diff = new TH1F("Time Diff", "Time Diff",3000,0.0,500);
106  Eoverp_plus_nocuts = new TH1F(" E over p plus no cuts "," E over p plus no cuts ", 1000,0.0,5.0);
107  Eoverp_minus_nocuts = new TH1F(" E over p minus no cuts ", " E over p minus no cuts ",1000,0.0,5.0);
108  Eoverp_plus_cuts = new TH1F(" E over p plus cuts "," E over p plus cuts ", 1000,0.0,5.0);
109  Eoverp_minus_cuts = new TH1F(" E over p minus cuts ", " E over p minus cuts ",1000,0.0,5.0);
110  Evsp_plus = new TH2F(" E vs p plus ", " E vs p plus", 1000,0.0,10.0, 1000, 0.0, 10.0);
111  Evsp_minus = new TH2F(" E vs p minus ", " E vs p minus", 1000,0.0,10.0, 1000, 0.0, 10.0);
112  Eoverpvsp_plus = new TH2F(" E over p vs p plus " , " E over p vs p plus " , 1000, 0.0, 10.0, 1000, 0.0, 10.0);
113  Eoverpvsp_minus = new TH2F(" E over p vs p minus " , " E over p vs p minus " , 1000,0.0,10.0,1000,0.0,10.0);
114  FCAL_Evsp = new TH2F(" FCAL E vs p " , " FCAL E vs p " ,500,0.0,5.0,250,0.0,2.5);
115  FCAL_Eoverpvsp = new TH2F( " FCAL Eoverp vs p " , " FCAL Eoverp vs p " , 500, 0.0, 5.0, 120, 0.0, 1.2);
116  FCAL_Eoverp_cuts = new TH1F ( " FCAL E over p cuts " , " FCAL E over p cuts " , 1000, 0.0, 1.0);
117  FCAL_Eoverp_nocuts = new TH1F( " FCAL E over p no cuts " , " FCAL E over p no cuts " , 1000,0.0,1.0);
118 
119  BCALShowerTrack_Energy = new TH1F("Charged energy","Charged energy",1000,0.0,2.0);
120  All_Layer_Energy = new TH1F("All Layer Point Energy","All Layer Point Energy",1500,0.0,1.0);
121 
122  //Erec_over_Ethrown_vs_z = new TH2F("Erec_over_Ethrown_vs_z","Erec_over_Ethrown_vs_z", 200, -50.0, 600.0, 200, 0.0, 2.0);
123  //Ecorr_over_Erec_vs_z = new TH2F("Ecorr_over_Erec_vs_z","Ecorr_over_Erec_vs_z", 200, -50.0, 600.0, 200, 0.0, 4.0);
124  //Ereconstructed_vs_Ethrown = new TH2F("Ereconstructed_vs_Ethrown","BCAL total reconstructed E to total thrown E", 200, 0.0, 6.0, 200, 0.0, 6.0);
125 
126  //Etot_truth = new TH1F("Etot_truth", "Sum of all truth showers (GeV)", 200, 0.0, 6.0);
127  Etot_hits = new TH1F("Etot_hits", "Sum of all hits (GeV)", 200, 0.0, 6.0);
128  //Etruth_over_Ethrown_vs_z = new TH2F("Etruth_over_Ethrown_vs_z","Etruth_over_Ethrown_vs_z", 200, -50.0, 600.0, 200, 0.0, 2.0);
129 
130  //Edeposited_over_Ethrown_vs_z = new TH2F("Edeposited_over_Ethrown_vs_z","Edeposited_over_Ethrown_vs_z", 80, -50.0, 450.0, 200, 0.0, 2.0);
131 
132  // Go back up to the parent directory
133  dir->cd("../");
134 
135 
136  BCALPoint_Charged_neg = new TTree("BCALPointChargedNeg"," from DBCALPoint object");
137  //BCAL_Energy->Branch("shower_energy",&shower_energy);
138  BCALPoint_Charged_neg->Branch("energy_point",&energy_point);
139  BCALPoint_Charged_neg->Branch("energy_shower",&energy_shower);
140  BCALPoint_Charged_neg->Branch("energy_raw_shower",&energy_raw_shower);
141  BCALPoint_Charged_neg->Branch("track_momentum",&track_momentum);
142  BCALPoint_Charged_neg->Branch("layer1_energysum",&layer1_energysum);
143  BCALPoint_Charged_neg->Branch("layer2_energysum",&layer2_energysum);
144  BCALPoint_Charged_neg->Branch("layer3_energysum",&layer3_energysum);
145  BCALPoint_Charged_neg->Branch("layer4_energysum",&layer4_energysum);
146  BCALPoint_Charged_neg->Branch("module",&module,"module/i");
147  BCALPoint_Charged_neg->Branch("sector",&sector,"sector/i");
148  BCALPoint_Charged_neg->Branch("layer",&layer,"layer/i");
149  BCALPoint_Charged_neg->Branch("eventnum",&eventnum,"eventnum/i");
150  BCALPoint_Charged_neg->Branch("theta",&theta);
151  BCALPoint_Charged_neg->Branch("channel",&channel,"channel/i");
152  //BCALPoint_Charged_neg->Branch("charge",&charge);
153  BCALPoint_Charged_neg->Branch("z", &z);
154  BCALPoint_Charged_neg->Branch("r", &r);
155  BCALPoint_Charged_neg->Branch("phi", &phi);
156  BCALPoint_Charged_neg->Branch("L1_pathlength", &L1_pathlength);
157  BCALPoint_Charged_neg->Branch("L2_pathlength", &L2_pathlength);
158  BCALPoint_Charged_neg->Branch("L3_pathlength", &L3_pathlength);
159  BCALPoint_Charged_neg->Branch("L4_pathlength", &L4_pathlength);
160 
161 
162  BCALPoint_Charged_pos = new TTree("BCALPointChargedPos"," from DBCALPoint object pos");
163  //BCAL_Energy->Branch("shower_energy",&shower_energy);
164  BCALPoint_Charged_pos->Branch("energy_point",&energy_point);
165  BCALPoint_Charged_pos->Branch("energy_shower",&energy_shower);
166  BCALPoint_Charged_pos->Branch("energy_raw_shower",&energy_raw_shower);
167  BCALPoint_Charged_pos->Branch("track_momentum",&track_momentum);
168  BCALPoint_Charged_pos->Branch("layer1_energysum",&layer1_energysum);
169  BCALPoint_Charged_pos->Branch("layer2_energysum",&layer2_energysum);
170  BCALPoint_Charged_pos->Branch("layer3_energysum",&layer3_energysum);
171  BCALPoint_Charged_pos->Branch("layer4_energysum",&layer4_energysum);
172  BCALPoint_Charged_pos->Branch("module",&module,"module/i");
173  BCALPoint_Charged_pos->Branch("sector",&sector,"sector/i");
174  BCALPoint_Charged_pos->Branch("layer",&layer,"layer/i");
175  BCALPoint_Charged_pos->Branch("eventnum",&eventnum,"eventnum/i");
176  BCALPoint_Charged_pos->Branch("theta",&theta);
177  BCALPoint_Charged_pos->Branch("channel",&channel,"channel/i");
178  //BCALPoint_Charged_neg->Branch("charge",&charge);
179  BCALPoint_Charged_pos->Branch("z", &z);
180  BCALPoint_Charged_pos->Branch("r", &r);
181  BCALPoint_Charged_pos->Branch("phi", &phi);
182  BCALPoint_Charged_pos->Branch("L1_pathlength", &L1_pathlength);
183  BCALPoint_Charged_pos->Branch("L2_pathlength", &L2_pathlength);
184  BCALPoint_Charged_pos->Branch("L3_pathlength", &L3_pathlength);
185  BCALPoint_Charged_pos->Branch("L4_pathlength", &L4_pathlength);
186 
187  BCALPoint_Neutral = new TTree("BCAL Point Neutral", " neutral bcal points");
188  // BCALPoint_Neutral->Branch("energy_point",&energy_point);
189  BCALPoint_Neutral->Branch("module",&module,"module/i");
190  BCALPoint_Neutral->Branch("sector",&sector,"sector/i");
191  BCALPoint_Neutral->Branch("layer",&layer,"layer/i");
192  BCALPoint_Neutral->Branch("eventnum",&eventnum,"eventnum/i");
193  BCALPoint_Neutral->Branch("theta",&theta);
194  BCALPoint_Neutral->Branch("channel",&channel,"channel/i");
195 
196  BCAL_Neutrals = new TTree("BCALNeutrals","BCALNeutrals");
197  BCAL_Neutrals->Branch("eventnum",&eventnum,"eventnum/i");
198  BCAL_Neutrals->Branch("E1",&E1);
199  BCAL_Neutrals->Branch("E1_raw",&E1_raw);
200  BCAL_Neutrals->Branch("E2",&E2);
201  BCAL_Neutrals->Branch("E2_raw",&E2_raw);
202  BCAL_Neutrals->Branch("p1_mag",&p1_mag);
203  BCAL_Neutrals->Branch("p1_raw_mag",&p1_raw_mag);
204  BCAL_Neutrals->Branch("p2_mag",&p2_mag);
205  BCAL_Neutrals->Branch("p2_raw_mag",&p2_raw_mag);
206  BCAL_Neutrals->Branch("inv_mass",&inv_mass);
207  BCAL_Neutrals->Branch("inv_mass_raw",&inv_mass_raw);
208  BCAL_Neutrals->Branch("x1",&x1);
209  BCAL_Neutrals->Branch("y1",&y1);
210  BCAL_Neutrals->Branch("x2",&x2);
211  BCAL_Neutrals->Branch("y2",&y2);
212  BCAL_Neutrals->Branch("z1",&z1);
213  BCAL_Neutrals->Branch("z2",&z2);
214  BCAL_Neutrals->Branch("phi1",&phi1);
215  BCAL_Neutrals->Branch("phi2",&phi2);
216  BCAL_Neutrals->Branch("t1",&t1);
217  BCAL_Neutrals->Branch("t2",&t2);
218  BCAL_Neutrals->Branch("vertexZ",&vertexZ);
219  BCAL_Neutrals->Branch("vertexX",&vertexX);
220  BCAL_Neutrals->Branch("vertexY",&vertexY);
221  BCAL_Neutrals->Branch("BCALShowers_per_event",&BCALShowers_per_event,"BCALShowers_per_event/i");
222  BCAL_Neutrals->Branch("channel",&channel);
223 
224  FCAL_Neutrals = new TTree("FCALNeutrals","FCALNeutrals");
225  FCAL_Neutrals->Branch("eventnum",&eventnum,"eventnum/i");
226  FCAL_Neutrals->Branch("E1",&E1);
227  FCAL_Neutrals->Branch("E2",&E2);
228  FCAL_Neutrals->Branch("p1_mag",&p1_mag);
229  FCAL_Neutrals->Branch("p2_mag",&p2_mag);
230  FCAL_Neutrals->Branch("inv_mass",&inv_mass);
231  FCAL_Neutrals->Branch("x1",&x1);
232  FCAL_Neutrals->Branch("x2",&x2);
233  FCAL_Neutrals->Branch("y1",&y1);
234  FCAL_Neutrals->Branch("y2",&y2);
235  FCAL_Neutrals->Branch("z1",&z1);
236  FCAL_Neutrals->Branch("z2",&z2);
237  FCAL_Neutrals->Branch("t1",&t1);
238  FCAL_Neutrals->Branch("t2",&t2);
239  FCAL_Neutrals->Branch("vertexZ",&vertexZ);
240  FCAL_Neutrals->Branch("vertexX",&vertexX);
241  FCAL_Neutrals->Branch("vertexY",&vertexY);
242  FCAL_Neutrals->Branch("FCALShowers_per_event",&FCALShowers_per_event,"FCALShowers_per_event/i");
243 
244  FCALClusterNeutrals = new TTree("FCALClusterNeutrals","FCALClusterNeutrals");
245  FCALClusterNeutrals->Branch("eventnum",&eventnum,"eventnum/i");
246  FCALClusterNeutrals->Branch("E1",&E1);
247  FCALClusterNeutrals->Branch("E2",&E2);
248  FCALClusterNeutrals->Branch("p1_mag",&p1_mag);
249  FCALClusterNeutrals->Branch("p2_mag",&p2_mag);
250  FCALClusterNeutrals->Branch("inv_mass",&inv_mass);
251  FCALClusterNeutrals->Branch("x1",&x1);
252  FCALClusterNeutrals->Branch("x2",&x2);
253  FCALClusterNeutrals->Branch("y1",&y1);
254  FCALClusterNeutrals->Branch("y2",&y2);
255  FCALClusterNeutrals->Branch("z1",&z1);
256  FCALClusterNeutrals->Branch("z2",&z2);
257  FCALClusterNeutrals->Branch("t1",&t1);
258  FCALClusterNeutrals->Branch("t2",&t2);
259  FCALClusterNeutrals->Branch("vertexZ",&vertexZ);
260  FCALClusterNeutrals->Branch("vertexX",&vertexX);
261  FCALClusterNeutrals->Branch("vertexY",&vertexY);
262  FCALClusterNeutrals->Branch("FCALClusters_per_event",&FCALShowers_per_event,"FCALClusters_per_event/i");
263 
264 
265  Triple_FCAL_Neutrals = new TTree("TripleFCALNeutrals","TripleFCALNeutrals");
266  Triple_FCAL_Neutrals->Branch("eventnum",&eventnum,"eventnum/i");
267  Triple_FCAL_Neutrals->Branch("E1",&E1);
268  Triple_FCAL_Neutrals->Branch("E2",&E2);
269  Triple_FCAL_Neutrals->Branch("E3",&E3);
270  Triple_FCAL_Neutrals->Branch("p1_mag",&p1_mag);
271  Triple_FCAL_Neutrals->Branch("p2_mag",&p2_mag);
272  Triple_FCAL_Neutrals->Branch("p3_mag",&p3_mag);
273  Triple_FCAL_Neutrals->Branch("inv_mass",&inv_mass);
274  Triple_FCAL_Neutrals->Branch("x1",&x1);
275  Triple_FCAL_Neutrals->Branch("x2",&x2);
276  Triple_FCAL_Neutrals->Branch("x3",&x3);
277  Triple_FCAL_Neutrals->Branch("y1",&y1);
278  Triple_FCAL_Neutrals->Branch("y2",&y2);
279  Triple_FCAL_Neutrals->Branch("y3",&y3);
280  Triple_FCAL_Neutrals->Branch("z1",&z1);
281  Triple_FCAL_Neutrals->Branch("z2",&z2);
282  Triple_FCAL_Neutrals->Branch("z3",&z3);
283  Triple_FCAL_Neutrals->Branch("t1",&t1);
284  Triple_FCAL_Neutrals->Branch("t2",&t2);
285  Triple_FCAL_Neutrals->Branch("t3",&t3);
286  Triple_FCAL_Neutrals->Branch("vertexZ",&vertexZ);
287  Triple_FCAL_Neutrals->Branch("vertexX",&vertexX);
288  Triple_FCAL_Neutrals->Branch("vertexY",&vertexY);
289  Triple_FCAL_Neutrals->Branch("FCALShowers_per_event",&FCALShowers_per_event,"FCALShowers_per_event/i");
290 
291 
292  Split_Gamma_Neutrals = new TTree ("SplitGammasNeutrals","SplitGammasNeutrals");
293  Split_Gamma_Neutrals->Branch("eventnum",&eventnum,"eventnum/i");
294  Split_Gamma_Neutrals->Branch("bcal_E",&bcal_E);
295  Split_Gamma_Neutrals->Branch("bcal_x",&bcal_x);
296  Split_Gamma_Neutrals->Branch("bcal_y",&bcal_y);
297  Split_Gamma_Neutrals->Branch("bcal_z",&bcal_z);
298  Split_Gamma_Neutrals->Branch("bcal_phi",&bcal_phi);
299  Split_Gamma_Neutrals->Branch("bcal_t",&bcal_t);
300  Split_Gamma_Neutrals->Branch("fcal_E",&fcal_E);
301  Split_Gamma_Neutrals->Branch("fcal_x",&fcal_x);
302  Split_Gamma_Neutrals->Branch("fcal_y",&fcal_y);
303  Split_Gamma_Neutrals->Branch("fcal_z",&fcal_z);
304  Split_Gamma_Neutrals->Branch("fcal_t",&fcal_t);
305  Split_Gamma_Neutrals->Branch("bcal_p",&bcal_p);
306  Split_Gamma_Neutrals->Branch("fcal_p",&fcal_p);
307  Split_Gamma_Neutrals->Branch("inv_mass",&inv_mass);
308  Split_Gamma_Neutrals->Branch("BCALShowers_per_event",&BCALShowers_per_event);
309  Split_Gamma_Neutrals->Branch("FCALShowers_per_event",&FCALShowers_per_event);
310  Split_Gamma_Neutrals->Branch("vertexZ",&vertexZ);
311  Split_Gamma_Neutrals->Branch("vertexX",&vertexX);
312  Split_Gamma_Neutrals->Branch("vertexY",&vertexY);
313 
314 
315  Split_Gamma_Neutrals_raw = new TTree ("SplitGammasNeutralsRaw","SplitGammasNeutralsRaw");
316  Split_Gamma_Neutrals_raw->Branch("eventnum",&eventnum,"eventnum/i");
317  Split_Gamma_Neutrals_raw->Branch("bcal_E",&bcal_E);
318  Split_Gamma_Neutrals_raw->Branch("bcal_x",&bcal_x);
319  Split_Gamma_Neutrals_raw->Branch("bcal_y",&bcal_y);
320  Split_Gamma_Neutrals_raw->Branch("bcal_z",&bcal_z);
321  Split_Gamma_Neutrals_raw->Branch("bcal_phi",&bcal_phi);
322  Split_Gamma_Neutrals_raw->Branch("bcal_t",&bcal_t);
323  Split_Gamma_Neutrals_raw->Branch("fcal_E",&fcal_E);
324  Split_Gamma_Neutrals_raw->Branch("fcal_x",&fcal_x);
325  Split_Gamma_Neutrals_raw->Branch("fcal_y",&fcal_y);
326  Split_Gamma_Neutrals_raw->Branch("fcal_z",&fcal_z);
327  Split_Gamma_Neutrals_raw->Branch("fcal_t",&fcal_t);
328  Split_Gamma_Neutrals_raw->Branch("bcal_p",&bcal_p);
329  Split_Gamma_Neutrals_raw->Branch("fcal_p",&fcal_p);
330  Split_Gamma_Neutrals_raw->Branch("inv_mass",&inv_mass);
331  Split_Gamma_Neutrals_raw->Branch("BCALShowers_per_event",&BCALShowers_per_event);
332  Split_Gamma_Neutrals_raw->Branch("FCALClusters_per_event",&FCALShowers_per_event);
333  Split_Gamma_Neutrals_raw->Branch("vertexZ",&vertexZ);
334  Split_Gamma_Neutrals_raw->Branch("vertexX",&vertexX);
335  Split_Gamma_Neutrals_raw->Branch("vertexY",&vertexY);
336 
337  dEventWriterROOT = NULL;
338  dEventWriterREST = NULL;
339 
340  return NOERROR;
341 }
342 
343 
344 //------------------
345 // brun
346 //------------------
347 jerror_t DEventProcessor_BCAL_Shower::brun(jana::JEventLoop* locEventLoop, int locRunNumber)
348 {
349  // This is called whenever the run number changes
350 
351  /*
352  //Optional: Retrieve REST writer for writing out skims
353  locEventLoop->GetSingle(dEventWriterREST);
354  */
355 
356  //vector<const DTrackFinder *> finders;
357  //locEventLoop->Get(finders);
358  //finder = const_cast<DTrackFinder*>(finders[0]);
359 
360  /*
361  //Recommeded: Create output ROOT TTrees (nothing is done if already created)
362  locEventLoop->GetSingle(dEventWriterROOT);
363  dEventWriterROOT->Create_DataTrees(locEventLoop);
364  */
365 
366  return NOERROR;
367 }
368 
369 //------------------
370 // evnt
371 //------------------
372 
373 
374 
375 
376 jerror_t DEventProcessor_BCAL_Shower::evnt(jana::JEventLoop* locEventLoop, uint64_t locEventNumber)
377 {
378  eventnum = locEventNumber;
379  // This is called for every event. Use of common resources like writing
380  // to a file or filling a histogram should be mutex protected. Using
381  // locEventLoop->Get(...) to get reconstructed objects (and thereby activating the
382  // reconstruction algorithm) should be done outside of any mutex lock
383  // since multiple threads may call this method at the same time.
384  //
385  // Here's an example:
386  //
387  // vector<const MyDataClass*> mydataclasses;
388  // locEventLoop->Get(mydataclasses);
389  //
390  // japp->RootWriteLock();
391  // ... fill historgrams or trees ...
392  // japp->RootUnLock();
393 
394  // DOCUMENTATION:
395  // ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software
396 
397 
398 
399 
400 
401 
402  vector<const DMCThrown*> locMCThrowns;
403  locEventLoop->Get(locMCThrowns);
404 
405  vector<const DBCALShower*> locBCALShowers;
406  vector<const DFCALShower*> locFCALShowers;
407  vector<const DBCALTruthShower*> truthshowers;
408  vector<const DMCThrown*> mcthrowns;
409  vector<const DBCALHit*> bcalhits;
410  vector<const DBCALCluster*> locBCALClusters;
411  vector<const DFCALCluster*> locFCALClusters;
412  vector<const DBCALPoint*> locBCALPoints;
413  vector<const DVertex*> kinfitVertex;
414  const DDetectorMatches* locDetectorMatches = NULL;
415  locEventLoop->GetSingle(locDetectorMatches);
416  locEventLoop->Get(locBCALShowers);
417  locEventLoop->Get(locFCALShowers);
418  locEventLoop->Get(bcalhits);
419  locEventLoop->Get(mcthrowns);
420  locEventLoop->Get(truthshowers);
421  locEventLoop->Get(locBCALClusters);
422  locEventLoop->Get(locFCALClusters);
423  locEventLoop->Get(locBCALPoints);
424  locEventLoop->Get(kinfitVertex);
425 
426  vector<const DTrackTimeBased*> locTrackTimeBased;
427  locEventLoop->Get(locTrackTimeBased);
428 
429  const DKinematicData* locKinematicData;
430  //locEventLoop->Get(locKinematicData);
431 
432 
433  DVector3 crude_vertex;
434  vector <const DChargedTrackHypothesis*> locChargedTrackHypothesis;
435  locEventLoop->Get(locChargedTrackHypothesis);
436 
437 // vector <const DTrackCandidate *> locTrackCandidate;
438 // locEventLoop->Get(locTrackCandidate,"StraightLine");
439 
440 // = static_cast<const DChargedTrackHypothesis* (locKinematicData);
441  //DChargedTrackHypothesis * locChargedTrack = const_cast<const DChargedTrackHypothesis *> (locChargedTrackHypothesis);
442  //const DChargedTrackHypothesis* locChargedTrack = static_cast<const DChargedTrackHypothesis* (locKinematicData);
443  //deque <const DChargedTrackHypothesis *> locChargedTrack = static_cast <const DChargedTrackHypothesis*>(locKinematicData);
444 
445 /* vector <const DNeutralParticleHypothesis *> locNeutralParticleHypothesis;
446  locEventLoop->Get(locNeutralParticleHypothesis);
447  double locTheta;
448  locTheta = locNeutralParticleHypothesis->momentum().Theta*180.0/TMath::Pi();
449  Theta->Fill(locTheta);*/
450 
451 
452  vector <const DChargedTrackHypothesis*> locParticles;
453  locEventLoop->Get(locParticles);
454  //vector <const DKinematicData*> locParticles;
455  //locEventLoop->Get(locParticles);
456 
457  // Although we are only filling objects local to this plugin, TTree::Fill() periodically writes to file: Global ROOT lock
458  japp->RootWriteLock(); //ACQUIRE ROOT LOCK
459 
460  vector <const DBCALShower *> matchedShowers;
461  vector < const DBCALShower *> matchedShowersneg;
462  vector < const DBCALShower *> matchedShowerspos;
463  vector <const DTrackTimeBased* > matchedTracks;
464  vector <const DFCALShower *> matchedFCALShowers;
465  vector < const DFCALCluster *> matchedFCALClusters;
466  DVector3 mypos(0.0,0.0,0.0);
467  DVector3 myposL1(0.0,0.0,0.0);
468  DVector3 myposL2(0.0,0.0,0.0);
469  DVector3 myposL3(0.0,0.0,0.0);
470  DVector3 myposL4(0.0,0.0,0.0);
471  DVector3 myposL5(0.0,0.0,0.0);
472  DVector3 myposL41(0.0,0.0,0.0);
473  DVector3 myposL42(0.0,0.0,0.0);
474  DVector3 myposL43(0.0,0.0,0.0);
475  double p;
476  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
477  for (unsigned int j=0; j< locBCALShowers.size(); ++j){
478 
479  double x = locBCALShowers[j]->x;
480  double y = locBCALShowers[j]->y;
481  double z = locBCALShowers[j]->z;
482  double E = locBCALShowers[j]->E;
483  DVector3 pos_bcal(x,y,z);
484  double R = pos_bcal.Perp();
485  double phi = pos_bcal.Phi();
486  double L2 = 0.81*2.54+65.0;
487  double L3 = L2 + 0.81*2.54*2;
488  double L4 = L3 + 0.81*2.54*3;
489  double L5 = L4 + 0.97*2.54*4;
490  double L41 = L4 + 0.97*2.54*4*1/4;
491  double L42 = L4 + 0.97*2.54*4*2/4;
492  double L43 = L4 + 0.97*2.54*4*3/4;
493  double FOM = TMath::Prob(locTrackTimeBased[i]->chisq, locTrackTimeBased[i]->Ndof);
494  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(R, mypos);
495  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(65.0,myposL1);
496  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L2,myposL2);
497  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L3,myposL3);
498  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L4,myposL4);
499  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L5,myposL5);
500  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L41,myposL41);
501  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L42,myposL42);
502  locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L43,myposL43);
503 
504  double q = locTrackTimeBased[i]->rt->q;
505  p = locTrackTimeBased[i]->momentum().Mag();
506  double dPhi = TMath::Abs(mypos.Phi()-pos_bcal.Phi());
507  double dZ = TMath::Abs(mypos.Z() - z);
508  // cout << " energy before matching = " << E << endl;
509  if( mypos.Perp() == R && p > 0.8 && dZ < 30.0 && dPhi < 0.18 && q == 1.0) Eoverp_plus_cuts->Fill(E/p);
510  if( mypos.Perp() == R && p > 0.8 && dZ < 30.0 && dPhi < 0.18 && q == -1.0) Eoverp_minus_cuts->Fill(E/p);
511  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0) Eoverp_plus_nocuts->Fill(E/p);
512  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0) Eoverp_minus_nocuts->Fill(E/p);
513  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0) Evsp_plus->Fill(p,E);
514  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0) Evsp_minus->Fill(p,E);
515  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0) Eoverpvsp_plus->Fill(p,E/p);
516  if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0) Eoverpvsp_minus->Fill(p,E/p);
517 
518  if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R) {
519  matchedShowers.push_back(locBCALShowers[j]);
520  matchedTracks.push_back(locTrackTimeBased[i]);
521  z_shower->Fill(E);
522  //cout << " energy after matching = " << E << " q = " << q << endl;
523  }
524  if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R && q == -1.0){
525  matchedShowersneg.push_back(locBCALShowers[j]);
526  }
527  if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R && q == 1.0){
528  matchedShowerspos.push_back(locBCALShowers[j]);
529  }
530  }
531  }
532 
533  for(unsigned int i = 0 ; i < locTrackTimeBased.size() ; ++i)
534  {
535  for(unsigned int j = 0 ; j < locFCALShowers.size(); ++j)
536  {
537  const DFCALShower *s2 = locFCALShowers[j];
538  //double dx, dy, dz, R, thetarad1, phirad1,
539  double x = s2->getPosition().X();
540  double y = s2->getPosition().Y();
541  double z = s2->getPosition().Z();
542  DVector3 fcalpos(x,y,z);
543  //cout << " x = " << x << " y = " << y << endl;
544  DVector3 norm(0.0,0.0,-1);
545  DVector3 pos;
546  locTrackTimeBased[i]->rt->GetIntersectionWithPlane(fcalpos,norm,pos);
547 
548  double diffX = TMath::Abs(x - pos.X());
549  double diffY = TMath::Abs(y - pos.Y());
550  if(diffX < 3.0 && diffY < 3.0)
551  {
552  matchedFCALShowers.push_back(locFCALShowers[j]);
553  }
554 
555  }
556  }
557 
558  for(unsigned int i = 0 ; i < locTrackTimeBased.size(); ++i)
559  {
560  for(unsigned int j = 0 ; j < locFCALClusters.size(); ++j)
561  {
562  const DFCALCluster *c1 = locFCALClusters[j];
563  double x = c1->getCentroid().X();
564  double y = c1->getCentroid().Y();
565  double z = c1->getCentroid().Z();
566  DVector3 fcalpos(x,y,z);
567  //cout << " x = " << x << " y = " << y << endl;
568  DVector3 norm(0.0,0.0,-1);
569  DVector3 pos;
570  locTrackTimeBased[i]->rt->GetIntersectionWithPlane(fcalpos,norm,pos);
571 
572  double diffX = TMath::Abs(x - pos.X());
573  double diffY = TMath::Abs(y - pos.Y());
574  if(diffX < 3.0 && diffY < 3.0)
575  {
576  matchedFCALClusters.push_back(locFCALClusters[j]);
577  }
578 
579  }
580  }
581 
582 
583  //What to do later
584  // if (matchedShowes.find( <<< cost BCALShower * >>> ) != matchedShower.end()) continue;
585 
586 
587  // DVector3 Calc_CrudeVertex(const deque<const DKinematicData*>& locParticles) const
588  //{
589  DVector3 locVertex(0.0,0.0,dTargetZCenter);
590  double locVertexZ, locVertexY, locVertexX;
591  //cout << " vertex before = " << locVertexZ << " dTargetZcenter = " << dTargetZCenter << endl;
592  //locVertexZ = locVertex.Z();
593 
594  if(locParticles.size() == 0)
595  return NOERROR;
596 
597  double locDOCA, locDOCA2, locSmallestDOCA, locTime0;
598  double locAverageTime = 0.0;
599  DVector3 locTempVertex;
600  DVector3 locPOCA;
601  DVector3 locDeltaVertex;
602  DVector3 locMomentum;
603 
604  locSmallestDOCA = 9.9E9;
605  if(locParticles.size()>1){
606  //cout << " locParticles.size() = " << locParticles.size() << endl;
607  for(int j = 0; j < (int(locParticles.size())-1); ++j)
608  {
609  //cout << " please be this far lol " << endl;
610  for(size_t k= j+1; k < locParticles.size(); ++k)
611  {
612  locDOCA = dAnalysisUtilities->Calc_DOCAVertex(locParticles[j],locParticles[k], locTempVertex);
613 
614  locSmallestDOCA = locDOCA;
615  locVertex = locTempVertex;
616  locVertexZ = locVertex.Z();
617  locVertexY = locVertex.Y();
618  locVertexX = locVertex.X();
619 
620  }
621  }
622  }
623 
624  double kinfitVertexX, kinfitVertexY, kinfitVertexZ, kinfitVertexT;
625  for (int i = 0 ; i < kinfitVertex.size(); i++)
626  {
627  kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
628  kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
629  kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
630  kinfitVertexT = kinfitVertex[i]->dSpacetimeVertex.T();
631  goodVertexZ->Fill(kinfitVertexZ);
632  }
633 
634 
635  //DVector3 poca;
636  //double d = finder->FindDoca(position1,momentum1,position2,momentum2, &poca);
637 
638 
639 /* for(int j = 0; j < locParticles.size(); ++j)
640  {
641  locDOCA2 = dAnalysisUtilities->Calc_DOCAToVertex(locParticles[j], locVertex, locPOCA);
642  locDeltaVertex = locPOCA - locParticles[j]->position();
643  locMomentum = locParticles[j]->momentum();
644  double locTime = locParticles[j]->time() + locDeltaVertex.Dot(locMomentum)*locParticles[j]->energy()/(29.9792458*locMomentum.Mag2());
645  locAverageTime += locTime;
646  }
647  locTime0 = locAverageTime/(double(locParticles.size())); */
648  //cout << " vertex after = " << locVertexZ << endl;
649  VertexZ->Fill(locVertexZ);
650  //if (FUCKTHISDEQUE.size() > 0){
651  //crude_vertex = dAnalysisUtilities->Calc_CrudeVertex(FUCKTHISDEQUE);
652  double VertexZ;
653  //VertexZ = crude_vertex.Z();
654  //cout << " vertex = " << locVertex << endl;
655 // }
656 // cout << "check " << endl;
657  //single photon details
658  double Etot_reconstructed = 0.0;
659  for(unsigned int i=0; i<locBCALShowers.size(); i++)
660  {
661  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
662  // continue;
663  const DBCALShower *s = locBCALShowers[i];
664  //xy_shower->Fill(s->x, s->y);
665  //if (locVertexZ > 50 && locVertexZ < 80)
666  //z_shower->Fill(s->E);
667  E_shower->Fill(s->E);
668  Etot_reconstructed += s->E;
669  }
670 
671  // test to reproduce neutral shower energy histos
672 
673 
674 
675  for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
676  {
677  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[loc_i]))
678  // continue;
679  if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[loc_i]) != matchedShowers.end()) continue;
680  const DBCALShower *s1 = locBCALShowers[loc_i];
681  const DBCALShower* a_shower = locBCALShowers[loc_i];
682  vector<const DBCALCluster*> associated_clusters;
683  a_shower->Get(associated_clusters);
684  double E1 = s1->E;
685  double z_one = s1->z;
686  Neutral_E->Fill(z_one);
687 
688  double E_shower = a_shower->E;
689  double E_shower_raw = a_shower->E_raw;
690 // if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
691  //if(associated_clusters.size() != 0) cout << " assoc cluster is not empty " << endl;
692 
693  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
694  {
695 
696  const DBCALCluster* a_cluster = associated_clusters[loc_j];
697  vector<const DBCALPoint*> associated_points;
698  a_cluster->Get(associated_points);
699 
700  if(associated_points.size() == 0) cout << " assoc points is empty " << endl;
701 
702  for(unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
703  {
704  int module2 = associated_points[loc_k]->module();
705  int layer = associated_points[loc_k]->layer();
706  double point_energy = associated_points[loc_k]->E();
707 
708  //cout << " associated point vector size = " << associated_points.size() << " module = " << module2 << " layer = " << associated_points[loc_k]->layer() << " sector = " << associated_points[loc_k]->sector() << " shower energy = " << E_shower << " raw shower energy = " << E_shower_raw << " point energy = " << point_energy << " associated cluster size = " << associated_clusters.size() << " shower sizse = " << locBCALShowers.size() << endl;
709 
710  //double theta = associated_clusters[loc_j]->theta();
711 
712  //double theta = associated_clusters.rho();
713  //if(layer==1) Layer1_Energy->Fill(point_energy);
714  //if(layer==2) Layer2_Energy->Fill(point_energy);
715  //if(layer==3) Layer3_Energy->Fill(point_energy);
716  //if(layer==4) Layer4_Energy->Fill(point_energy);
717 
718 
719  //Assoc_E->Fill(point_energy);
720  }
721  }
722  }
723 
724 // investigating DBCALPoint members
725 
726  for(unsigned int i=0; i< locBCALPoints.size(); i++)
727  {
728  const DBCALPoint *p1 = locBCALPoints[i];
729  double E = p1->E();
730  double r = p1->r();
731  double phi = p1->phi();
732  double theta = p1->theta();
733  double z = p1->z();
734  int module = p1->module();
735  int layer = p1->layer();
736  int sector = p1->sector();
737  Point_Energy->Fill(E);
738  Point_Module->Fill(module);
739  Point_z->Fill(z);
740  //cout << " Point vector size = " << locBCALPoints.size() << " module = " << module << endl;
741  }
742  for(unsigned int i=0; i< locBCALClusters.size(); i++)
743  {
744  const DBCALCluster *c1 = locBCALClusters[i];
745  double E = c1->E();
746  Cluster_Energy->Fill(E);
747  }
748 
749 
750 /*
751 
752 
753  for(unsigned int i=0; i<locBCALShowers.size(); i++)
754  {
755  if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
756  continue;
757  const DBCALShower *s1 = locBCALShowers[i];
758  vector<const DBCALCluster*> associated_clusters1;
759  s1->Get(associated_clusters1);
760  double dx1 = s1->x - locVertexX;
761  double dy1 = s1->y - locVertexY;
762  double dz1 = s1->z - locVertexZ;
763  double dt1 = s1->t;
764  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
765  //double path1 = sqrt((dx1-locVertexX)*(dx1-locVertexX)+(dy1-locVertexY)*(dy1-locVertexY)+(dz1)*(dz1));
766  double E1 = s1->E;
767  double ECalc = s1->E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
768  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
769  double thetadeg1, thetarad1, phideg1, phirad1;
770  thetadeg1 = p1.Theta()*180.0/TMath::Pi();
771  thetarad1 = p1.Theta();
772  phideg1 = p1.Phi()*180.0/TMath::Pi();
773  phirad1 = p1.Phi();
774  TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
775  for(unsigned int loc_j = 0; loc_j < associated_clusters1.size(); loc_j++)
776  {
777  const DBCALCluster* a_cluster1 = associated_clusters1[loc_j];
778  vector<const DBCALPoint*> associated_points1;
779  a_cluster1->Get(associated_points1);
780  for(unsigned int loc_k = 0; loc_k < associated_points1.size(); loc_k++)
781  {
782  int module1 = associated_points1[loc_k]->module();
783  }
784  }
785  for(unsigned int j=i+1; j<locBCALShowers.size(); j++){
786  if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[j]))
787  continue;
788  const DBCALShower *s2 = locBCALShowers[j];
789  vector<const DBCALCluster*> associated_clusters2;
790  s2->Get(associated_clusters2);
791  double dx2 = s2->x - locVertexX;
792  double dy2 = s2->y - locVertexY;
793  double dz2 = s2->z - locVertexZ; // shift to coordinate relative to center of target
794  double dt2 = s2->t;
795  double R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
796  //double path2 = sqrt((dx2-locVertexX)*(dx2-locVertexX)+(dy2-locVertexY)*(dy2-locVertexY)+(dz2)*(dz2));
797  double E2 = s2->E;
798  double ECalc = s2->E*(1.106+(dz2+65.0-208.4)*(dz2+65.0-208.4)*6.851E-6);
799  TLorentzVector p2(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
800  TLorentzVector p2Calc(ECalc*dx2/R2, ECalc*dy2/R2, ECalc*dz2/R2, ECalc);
801  double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
802  thetarad2 = p2.Theta();
803  phirad2 = p2.Phi();
804  thetadeg2 = p2.Theta()*180.0/TMath::Pi();
805  phideg2 = p2.Phi()*180.0/TMath::Pi();
806  cospsi = sin(thetarad1)*sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
807  psi=acos(cospsi)*180/TMath::Pi();
808  TLorentzVector ptot = p1+p2;
809 
810  Theta_Hist_Both_BCAL->Fill(thetadeg1);
811  Theta_Hist_Both_BCAL->Fill(thetadeg2);
812  Phi_Hist_Both_BCAL->Fill(phideg2);
813  Phi_Hist_Both_BCAL->Fill(phideg2);
814  Psi_Hist_Both_BCAL->Fill(psi);
815  double makes_sense = 0;
816  if (R1 > R2 && dt1 > dt2) makes_sense = 1;
817  if (R1 < R2 && dt1 < dt2) makes_sense = 1;
818  if (R1 < R2 && dt1 > dt2) makes_sense = 0;
819  if (R2 < R1 && dt1 > dt1) makes_sense = 0;
820 
821  for(unsigned int loc_jj = 0; loc_jj < associated_clusters1.size(); loc_jj++)
822  {
823  const DBCALCluster* a_cluster2 = associated_clusters2[loc_jj];
824  vector<const DBCALPoint*> associated_points2;
825  a_cluster2->Get(associated_points2);
826 
827  for(unsigned int loc_kk = 0; loc_kk < associated_points2.size(); loc_kk++)
828  {
829  //int module2 = associated_points2[loc_kk]->module();
830 
831  if ( locVertexZ > 55.0 && locVertexZ < 75.0 && E1 > 0.2 && E2 > 0.2 && makes_sense==1 && locVertexX < 15.0 && locVertexX > -15.0 && locVertexY > -15.0 && locVertexY < 15.0) two_gamma_mass->Fill(ptot.M());
832  }
833  }
834 
835 
836 
837  }
838  }
839 
840 */
841 
842 
843 
844 
845 
846 
847 
848 // MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM MIN IONIZE REAL DATA MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
849 
850 
851 
852  for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
853  {
854  const DBCALShower* locBCALShower = locBCALShowers[loc_i];
855  double E_before = locBCALShower->E;
856  //cout << " E before = " << E_before << endl;
857  if(find(matchedShowers.begin(), matchedShowers.end(), locBCALShower) == matchedShowers.end()) continue;
858 
859  double x = locBCALShower->x;
860  double y = locBCALShower->y;
861  double z = locBCALShower->z;
862  double E1 = locBCALShower->E;
863  double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
864  double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
865  double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
866  double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
867  double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
868  double fourth_pos2_test = TMath::Sqrt( TMath::Power(myposL4.X()-myposL41.X(),2)+TMath::Power(myposL4.Y()-myposL41.Y(),2)+TMath::Power(myposL4.Z()-myposL41.Z(),2));
869  double fourth_pos3_test = TMath::Sqrt( TMath::Power(myposL41.X()-myposL42.X(),2)+TMath::Power(myposL41.Y()-myposL42.Y(),2)+TMath::Power(myposL41.Z()-myposL42.Z(),2));
870  double fourth_pos4_test = TMath::Sqrt( TMath::Power(myposL42.X()-myposL43.X(),2)+TMath::Power(myposL42.Y()-myposL43.Y(),2)+TMath::Power(myposL42.Z()-myposL43.Z(),2));
871  double fourth_pos5_test = TMath::Sqrt( TMath::Power(myposL43.X()-myposL5.X(),2)+TMath::Power(myposL43.Y()-myposL5.Y(),2)+TMath::Power(myposL43.Z()-myposL5.Z(),2));
872 
873  double less_approx_dist = fourth_pos2_test + fourth_pos3_test + fourth_pos4_test + fourth_pos5_test;
874 
875  double L1_pathlength = second_pos - first_pos;
876  double L2_pathlength = third_pos - second_pos;
877  double L3_pathlength = fourth_pos - third_pos;
878  double L4_pathlength = fifth_pos - fourth_pos;
879 
880 
881  //cout << " first pos = " << first_pos << " second_pos = " << second_pos << " third pos = " << third_pos << " fourth pos = " << fourth_pos << " fifth pos = " << fifth_pos << " event num = " << eventnum << endl;
882  //cout << " E after = " << E1 << endl;
883  double E_shower = locBCALShower->E;
884  double E_shower_raw = locBCALShower->E_raw;
885  double x1 = locBCALShower->x - kinfitVertexX;
886  double y1 = locBCALShower->y - kinfitVertexY;
887  double z1 = locBCALShower->z - kinfitVertexZ;
888  double t1 = locBCALShower->t;
889  double R1 = sqrt(x1*x1 + y1*y1 + z1*z1);
890  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
891  TLorentzVector p1(E1*x1/R1, E1*y1/R1, E1*z1/R1, E1);
892  vector<const DBCALCluster*> associated_clusters;
893  locBCALShower->Get(associated_clusters);
894 
895  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
896 
897  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
898  {
899 
900  const DBCALCluster* a_cluster = associated_clusters[loc_j];
901  vector<const DBCALPoint*> associated_points;
902  a_cluster->Get(associated_points);
903 
904  if(associated_points.size() == 0) cout << " assoc points is empty " << endl;
905  double Layer1_Energy_Sum=0.0;
906  double Layer2_Energy_Sum=0.0;
907  double Layer3_Energy_Sum=0.0;
908  double Layer4_Energy_Sum=0.0;
909 
910  double Layer1_Energy_Sum2=0.0;
911  double Layer2_Energy_Sum2=0.0;
912  double Layer3_Energy_Sum2=0.0;
913  double Layer4_Energy_Sum2=0.0;
914 
915  double Layer1_Energy_Sumpos = 0.0;
916  double Layer2_Energy_Sumpos = 0.0;
917  double Layer3_Energy_Sumpos = 0.0;
918  double Layer4_Energy_Sumpos = 0.0;
919  for(unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
920  {
921  const DBCALPoint* a_point = associated_points[loc_k];
922  vector<const DBCALUnifiedHit*> associated_unifiedhits;
923  a_point->Get(associated_unifiedhits);
924  int point_size = associated_points.size();
925  int module = associated_points[loc_k]->module();
926  int layer = associated_points[loc_k]->layer();
927  int sector = associated_points[loc_k]->sector();
928  double theta = associated_points[loc_k]->theta();
929  double sintheta = sin(theta);
930  double point_energy = associated_points[loc_k]->E();
931  int channel_per_module;
932  // int charge = locChargedTrackHypothesis->charge();
933  double z = associated_points[loc_k]->z();
934  double r = associated_points[loc_k]->r();
935  double theta_wrt_vertex = atan(r/(z-kinfitVertexZ));
936  double point_energy_sum;
937  double point_energy_over_angle = TMath::Abs(point_energy/sin(theta_wrt_vertex));
938  point_energy_sum += point_energy_over_angle;
939  int logical = 1;
940  // if( layer == 1 && point_energy/sintheta/point_energy_sum > 0.13) logical = 0;
941  // if( layer == 2 && point_energy/sintheta/point_energy_sum > 2*1.5/associated_points.size()) logical = 0;
942  // if( layer == 3 && point_energy/sintheta/point_energy_sum > 3*1.5/associated_points.size()) logical = 0;
943  // if( layer == 4 && point_energy/sintheta/point_energy_sum > 4*1.5/associated_points.size()) logical = 0;
944  //if(E_shower>0.3) logical = 0;
945  if(associated_points.size()<4) logical = 0;
946  //if(layer==1 && point_energy/sintheta >0.4) logical=0;
947  //if(layer==2 && (point_energy/sintheta < 0.2 || point_energy/sintheta > 0.6)) logical=0;
948  //if(layer==3 && (point_energy/sintheta < 0.4 || point_energy/sintheta > 0.8)) logical=0;
949  //if(layer==4 && (point_energy/sintheta < 0.6 || point_energy/sintheta > 1.0)) logical=0;
950  if (layer == 1 && sector == 1) channel_per_module = 0;
951  if (layer == 1 && sector == 2) channel_per_module = 1;
952  if (layer == 1 && sector == 3) channel_per_module = 2;
953  if (layer == 1 && sector == 4) channel_per_module = 3;
954  if (layer == 2 && sector == 1) channel_per_module = 4;
955  if (layer == 2 && sector == 2) channel_per_module = 5;
956  if (layer == 2 && sector == 3) channel_per_module = 6;
957  if (layer == 2 && sector == 4) channel_per_module = 7;
958  if (layer == 3 && sector == 1) channel_per_module = 8;
959  if (layer == 3 && sector == 2) channel_per_module = 9;
960  if (layer == 3 && sector == 3) channel_per_module = 10;
961  if (layer == 3 && sector == 4) channel_per_module = 11;
962  if (layer == 4 && sector == 1) channel_per_module = 12;
963  if (layer == 4 && sector == 2) channel_per_module = 13;
964  if (layer == 4 && sector == 3) channel_per_module = 14;
965  if (layer == 4 && sector == 4) channel_per_module = 15;
966  int channel = channel_per_module + (module-1)*16;
967 
968  //cout << " point energy = " << point_energy << " E shower = " << E1 << " logical = " << logical << " sine = " << sin(theta_wrt_vertex) << " points size = " << associated_points.size() << endl;
969  if(associated_points.size()>3)
970  {
971 
972  if(layer==1)Layer1_Energy_Sum += point_energy;
973  if(layer==2)Layer2_Energy_Sum += point_energy;
974  if(layer==3)Layer3_Energy_Sum += point_energy;
975  if(layer==4)Layer4_Energy_Sum += point_energy;
976  //cout << " Layer 1 Energy Sum = " << Layer1_Energy_Sum << " layer 2 e sum = " << Layer2_Energy_Sum << " Layer3_Energy_Sum = " << Layer3_Energy_Sum << " Layer 4 e sum = " << Layer4_Energy_Sum << " point size = " << associated_points.size() << " event num = " << eventnum << " iteration = " << loc_k << endl;
977 
978  if(layer==1)Layer1_Energy_Sum2 += point_energy/sintheta;
979  if(layer==2)Layer2_Energy_Sum2 += point_energy/sintheta;
980  if(layer==3)Layer3_Energy_Sum2 += point_energy/sintheta;
981  if(layer==4)Layer4_Energy_Sum2 += point_energy/sintheta;
982 
983  if(layer==1)Layer1_Energy_Sumpos += point_energy/L1_pathlength;
984  if(layer==2)Layer2_Energy_Sumpos += point_energy/L2_pathlength;
985  if(layer==3)Layer3_Energy_Sumpos += point_energy/L3_pathlength;
986  if(layer==4)Layer4_Energy_Sumpos += point_energy/L4_pathlength;
987  All_Layer_Energy->Fill(point_energy/sin(theta_wrt_vertex));
988 
989  if(layer==1) Layer1_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
990  if(layer==2) Layer2_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
991  if(layer==3) Layer3_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
992  if(layer==4) Layer4_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
993  if(loc_k+1 == associated_points.size() ){
994  if(Layer1_Energy_Sum < Layer4_Energy_Sum){
995  if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001) Layer1_Energy_v2->Fill(Layer1_Energy_Sum);
996  if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001) Layer2_Energy_v2->Fill(Layer2_Energy_Sum);
997  if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001) Layer3_Energy_v2->Fill(Layer3_Energy_Sum);
998  if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001) Layer4_Energy_v2->Fill(Layer4_Energy_Sum);
999  }
1000  if(Layer1_Energy_Sum2 < Layer4_Energy_Sum2 && Layer1_Energy_Sum2 > 0.005){
1001  Layer1_Energy->Fill(Layer1_Energy_Sum2);
1002  Layer2_Energy->Fill(Layer2_Energy_Sum2);
1003  Layer3_Energy->Fill(Layer3_Energy_Sum2);
1004  Layer4_Energy->Fill(Layer4_Energy_Sum2);
1005  }
1006  if(Layer1_Energy_Sumpos < Layer4_Energy_Sumpos && Layer1_Energy_Sumpos > 0.005 && Layer4_Energy_Sumpos > 0.005){
1007  Layer1_Energy_pos->Fill(Layer1_Energy_Sumpos);
1008  Layer2_Energy_pos->Fill(Layer2_Energy_Sumpos);
1009  Layer3_Energy_pos->Fill(Layer3_Energy_Sumpos);
1010  Layer4_Energy_pos->Fill(Layer4_Energy_Sumpos);
1011  //cout << " approx dist = " << L4_pathlength << " less approx dist = " << less_approx_dist << " ratio = " << less_approx_dist/L4_pathlength << endl;
1012  }
1013  //cout << " layer 1 e sum = " << Layer1_Energy_Sum2 << " layer 2 e sum = " << Layer2_Energy_Sum2 << " layer 3 e sum = " << Layer3_Energy_Sum2 << " layer 4 e sum = " << Layer4_Energy_Sum2 << " layer 1 path = " << L1_pathlength << " layer 2 path = " << L2_pathlength << " layer 3 path = " << L3_pathlength << " layer 4 path = " << L4_pathlength << " event num = " << eventnum << endl;
1014  }
1015  }
1016  //cout << " charge = " << locChargedTrackHypothesis->charge() << " point energy/sintheta= " << point_energy/sintheta1 << " 'shower energy' = " << E1 << " shower momentum = " << p1.M() << " channel = " << channel << " module = " << module1 << " layer = " << layer1 << " point size = " << associated_points.size() << endl;
1017  for(unsigned int loc_m = 0; loc_m < associated_unifiedhits.size(); loc_m++)
1018  {
1019  int modulehit = associated_unifiedhits[loc_m]->module;
1020  int layerhit = associated_unifiedhits[loc_m]->layer;
1021  int sectorhit = associated_unifiedhits[loc_m]->sector;
1022  int end = associated_unifiedhits[loc_m]->end;
1023  double unifiedhit_energy = associated_unifiedhits[loc_m]->E;
1024  //cout << " point E = " << point_energy << " hit E = " << unifiedhit_energy << " module hit = " << modulehit << " module point = " << module1 << " layer hit = " << layerhit << " layer point = " << layer1 << " end = " << end << " point size = " << associated_points.size() << " hit size = " << associated_unifiedhits.size() << endl;
1025 
1026 
1027  }
1028 
1029 
1030 
1031  }
1032  }
1033  //}
1034 
1035  }
1036  //}
1037  //}
1038 
1039 
1040 /*
1041  vector<const DChargedTrack*> locChargedTracks;
1042  locEventLoop->Get(locChargedTracks);
1043 
1044  for(size_t i = 0 ; i < locChargedTracks.size(); ++i)
1045  {
1046 
1047  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[i]->Get_BestFOM();
1048  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
1049  const DShowerMatchParams& locBCALShowerMatchParams = locChargedTrackHypothesis->dBCALShowerMatchParams;
1050 
1051  if(locBCALShowerMatchParams.dTrackTimeBased != NULL)
1052  {
1053  const DBCALShower* locBCALShower = dynamic_cast <const DBCALShower*>(locBCALShowerMatchParams.dShowerObject);
1054  if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower)){
1055  BCALShowerTrack_Energy->Fill(locBCALShower->E);
1056 
1057  //for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1058  //{
1059  //const DBCALShower* s1 = locBCALShowers;
1060  //const DBCALShower* a_shower = locBCALShowers;
1061 
1062  double E1 = locBCALShower->E;
1063  double E_shower = locBCALShower->E;
1064  double E_shower_raw = locBCALShower->E_raw;
1065  double x1 = locBCALShower->x - kinfitVertexX;
1066  double y1 = locBCALShower->y - kinfitVertexY;
1067  double z1 = locBCALShower->z - kinfitVertexZ;
1068  double t1 = locBCALShower->t;
1069  double R1 = sqrt(x1*x1 + y1*y1 + z1*z1);
1070  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
1071  TLorentzVector p1(E1*x1/R1, E1*y1/R1, E1*z1/R1, E1);
1072  vector<const DBCALCluster*> associated_clusters;
1073  locBCALShower->Get(associated_clusters);
1074 
1075  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
1076 
1077  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1078  {
1079 
1080  const DBCALCluster* a_cluster = associated_clusters[loc_j];
1081  vector<const DBCALPoint*> associated_points;
1082  a_cluster->Get(associated_points);
1083 
1084  if(associated_points.size() == 0) cout << " assoc points is empty " << endl;
1085 
1086  for(unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
1087  {
1088  const DBCALPoint* a_point = associated_points[loc_k];
1089  vector<const DBCALUnifiedHit*> associated_unifiedhits;
1090  a_point->Get(associated_unifiedhits);
1091  int point_size = associated_points.size();
1092  int module = associated_points[loc_k]->module();
1093  int layer = associated_points[loc_k]->layer();
1094  int sector = associated_points[loc_k]->sector();
1095  double theta = associated_points[loc_k]->theta();
1096  double sintheta = sin(theta);
1097  double point_energy = associated_points[loc_k]->E();
1098  int channel_per_module;
1099  int charge = locChargedTrackHypothesis->charge();
1100  double z = associated_points[loc_k]->z();
1101  double r = associated_points[loc_k]->r();
1102  double theta_wrt_vertex = atan(r/(z-kinfitVertexZ));
1103  double point_energy_sum;
1104  point_energy_sum += point_energy/sin(theta_wrt_vertex);
1105  int logical = 1;
1106  // if( layer == 1 && point_energy/sintheta/point_energy_sum > 0.13) logical = 0;
1107  // if( layer == 2 && point_energy/sintheta/point_energy_sum > 2*1.5/associated_points.size()) logical = 0;
1108  // if( layer == 3 && point_energy/sintheta/point_energy_sum > 3*1.5/associated_points.size()) logical = 0;
1109  // if( layer == 4 && point_energy/sintheta/point_energy_sum > 4*1.5/associated_points.size()) logical = 0;
1110  if(E_shower>0.3) logical = 0;
1111  if(associated_points.size()<4) logical = 0;
1112  //if(layer==1 && point_energy/sintheta >0.4) logical=0;
1113  //if(layer==2 && (point_energy/sintheta < 0.2 || point_energy/sintheta > 0.6)) logical=0;
1114  //if(layer==3 && (point_energy/sintheta < 0.4 || point_energy/sintheta > 0.8)) logical=0;
1115  //if(layer==4 && (point_energy/sintheta < 0.6 || point_energy/sintheta > 1.0)) logical=0;
1116  if (layer == 1 && sector == 1) channel_per_module = 0;
1117  if (layer == 1 && sector == 2) channel_per_module = 1;
1118  if (layer == 1 && sector == 3) channel_per_module = 2;
1119  if (layer == 1 && sector == 4) channel_per_module = 3;
1120  if (layer == 2 && sector == 1) channel_per_module = 4;
1121  if (layer == 2 && sector == 2) channel_per_module = 5;
1122  if (layer == 2 && sector == 3) channel_per_module = 6;
1123  if (layer == 2 && sector == 4) channel_per_module = 7;
1124  if (layer == 3 && sector == 1) channel_per_module = 8;
1125  if (layer == 3 && sector == 2) channel_per_module = 9;
1126  if (layer == 3 && sector == 3) channel_per_module = 10;
1127  if (layer == 3 && sector == 4) channel_per_module = 11;
1128  if (layer == 4 && sector == 1) channel_per_module = 12;
1129  if (layer == 4 && sector == 2) channel_per_module = 13;
1130  if (layer == 4 && sector == 3) channel_per_module = 14;
1131  if (layer == 4 && sector == 4) channel_per_module = 15;
1132  int channel = channel_per_module + (module-1)*16;
1133  if(charge < 0 && logical == 1)
1134  {
1135  double Layer1_Energy_Sum;
1136  double Layer2_Energy_Sum;
1137  double Layer3_Energy_Sum;
1138  double Layer4_Energy_Sum;
1139  if(layer==1)Layer1_Energy_Sum += point_energy/sin(theta_wrt_vertex);
1140  if(layer==2)Layer2_Energy_Sum += point_energy/sin(theta_wrt_vertex);
1141  if(layer==3)Layer3_Energy_Sum += point_energy/sin(theta_wrt_vertex);
1142  if(layer==4)Layer4_Energy_Sum += point_energy/sin(theta_wrt_vertex);
1143  double Layer1_Energy_Sum2;
1144  double Layer2_Energy_Sum2;
1145  double Layer3_Energy_Sum2;
1146  double Layer4_Energy_Sum2;
1147  if(layer==1)Layer1_Energy_Sum2 += point_energy/sintheta;
1148  if(layer==2)Layer2_Energy_Sum2 += point_energy/sintheta;
1149  if(layer==3)Layer3_Energy_Sum2 += point_energy/sintheta;
1150  if(layer==4)Layer4_Energy_Sum2 += point_energy/sintheta;
1151  All_Layer_Energy->Fill(point_energy/sin(theta_wrt_vertex));
1152 
1153  if(layer==1) Layer1_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
1154  if(layer==2) Layer2_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
1155  if(layer==3) Layer3_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
1156  if(layer==4) Layer4_Energy_vs_Channel->Fill(channel, point_energy/sin(theta_wrt_vertex));
1157  if(loc_k+1 == associated_points.size() ){
1158  if(Layer1_Energy_Sum < Layer4_Energy_Sum){
1159  if(Layer1_Energy_Sum > 0.001) Layer1_Energy_v2->Fill(Layer1_Energy_Sum);
1160  if(Layer2_Energy_Sum > 0.001) Layer2_Energy_v2->Fill(Layer2_Energy_Sum);
1161  if(Layer3_Energy_Sum > 0.001) Layer3_Energy_v2->Fill(Layer3_Energy_Sum);
1162  if(Layer4_Energy_Sum > 0.001) Layer4_Energy_v2->Fill(Layer4_Energy_Sum);
1163  }
1164  if(Layer1_Energy_Sum2 < Layer4_Energy_Sum2){
1165  Layer1_Energy->Fill(Layer1_Energy_Sum2);
1166  Layer2_Energy->Fill(Layer2_Energy_Sum2);
1167  Layer3_Energy->Fill(Layer3_Energy_Sum2);
1168  Layer4_Energy->Fill(Layer4_Energy_Sum2);
1169  }
1170  }
1171  cout << " point energy/sintheta= " << point_energy/sintheta << " layer 1 e sum = " << Layer1_Energy_Sum2 << " layer 2 e sum = " << Layer2_Energy_Sum2 << " layer 3 e sum = " << Layer3_Energy_Sum2 << " layer 4 e sum = " << Layer4_Energy_Sum2 << "evenr num = " << eventnum << endl;
1172  }
1173  //cout << " charge = " << locChargedTrackHypothesis->charge() << " point energy/sintheta= " << point_energy/sintheta1 << " 'shower energy' = " << E1 << " shower momentum = " << p1.M() << " channel = " << channel << " module = " << module1 << " layer = " << layer1 << " point size = " << associated_points.size() << endl;
1174  for(unsigned int loc_m = 0; loc_m < associated_unifiedhits.size(); loc_m++)
1175  {
1176  int modulehit = associated_unifiedhits[loc_m]->module;
1177  int layerhit = associated_unifiedhits[loc_m]->layer;
1178  int sectorhit = associated_unifiedhits[loc_m]->sector;
1179  int end = associated_unifiedhits[loc_m]->end;
1180  double unifiedhit_energy = associated_unifiedhits[loc_m]->E;
1181  //cout << " point E = " << point_energy << " hit E = " << unifiedhit_energy << " module hit = " << modulehit << " module point = " << module1 << " layer hit = " << layerhit << " layer point = " << layer1 << " end = " << end << " point size = " << associated_points.size() << " hit size = " << associated_unifiedhits.size() << endl;
1182  if(locChargedTrackHypothesis->charge() < 0.0 && module == 1 && sector == 1 && layer == 1 && logical ==1) Point_E_M1S1L1->Fill(point_energy/sintheta);
1183  if(locChargedTrackHypothesis->charge() < 0.0 && module == 12 && sector == 2 && layer == 2 && logical ==1) Point_E_M12S2L2->Fill(point_energy/sintheta);
1184  if(locChargedTrackHypothesis->charge() < 0.0 && module == 25 && sector == 3 && layer == 3 && logical == 1) Point_E_M25S3L3->Fill(point_energy/sintheta);
1185  if(locChargedTrackHypothesis->charge() < 0.0 && module == 37 && sector == 4 && layer == 4 && logical == 1) Point_E_M37S4L4->Fill(point_energy/sintheta);
1186 
1187  }
1188 
1189 
1190 
1191  }
1192  }
1193  //}
1194 
1195  }
1196  }
1197  }
1198 
1199 
1200 */
1201 
1202 
1203 
1204 
1205 
1206 
1207 /* for(size_t i = 0 ; i < locChargedTracks.size(); ++i)
1208  {
1209 
1210  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[i]->Get_BestFOM();
1211  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
1212  const DShowerMatchParams& locBCALShowerMatchParams = locChargedTrackHypothesis->dBCALShowerMatchParams;
1213 
1214  if(locBCALShowerMatchParams.dTrackTimeBased != NULL)
1215  {
1216  const DBCALShower* locBCALShower = dynamic_cast <const DBCALShower*>(locBCALShowerMatchParams.dShowerObject);
1217  if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower)){
1218  BCALShowerTrack_Energy->Fill(locBCALShower->E);
1219 
1220  //for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1221  //{
1222  //const DBCALShower* s1 = locBCALShowers;
1223  //const DBCALShower* a_shower = locBCALShowers;
1224 
1225  double E1 = locBCALShower->E;
1226  double E_shower = locBCALShower->E;
1227  double E_shower_raw = locBCALShower->E_raw;
1228  double x1 = locBCALShower->x - kinfitVertexX;
1229  double y1 = locBCALShower->y - kinfitVertexY;
1230  double z1 = locBCALShower->z - kinfitVertexZ;
1231  double t1 = locBCALShower->t;
1232  double R1 = sqrt(x1*x1 + y1*y1 + z1*z1);
1233  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
1234  TLorentzVector p1(E1*x1/R1, E1*y1/R1, E1*z1/R1, E1);
1235  vector<const DBCALCluster*> associated_clusters;
1236  locBCALShower->Get(associated_clusters);
1237 
1238  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
1239 
1240  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1241  {
1242 
1243  const DBCALCluster* a_cluster = associated_clusters[loc_j];
1244  vector<const DBCALPoint*> associated_points;
1245  a_cluster->Get(associated_points);
1246 
1247  if(associated_points.size() == 0) cout << " assoc points is empty " << endl;
1248 
1249  for(unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
1250  {
1251  const DBCALPoint* a_point = associated_points[loc_k];
1252  vector<const DBCALUnifiedHit*> associated_unifiedhits;
1253  a_point->Get(associated_unifiedhits);
1254  int module = associated_points[loc_k]->module();
1255  int layer = associated_points[loc_k]->layer();
1256  int sector = associated_points[loc_k]->sector();
1257  double theta = associated_points[loc_k]->theta();
1258  double sintheta = sin(theta);
1259  double point_energy = associated_points[loc_k]->E();
1260  int channel_per_module;
1261  int charge = locChargedTrackHypothesis->charge();
1262  float point_energy_sum;
1263  point_energy_sum += point_energy;
1264  int logical = 1;
1265  if( layer == 1 && point_energy/sintheta/point_energy_sum > 1.5/associated_points.size()) logical = 0;
1266  if( layer == 2 && point_energy/sintheta/point_energy_sum > 2*1.5/associated_points.size()) logical = 0;
1267  if( layer == 3 && point_energy/sintheta/point_energy_sum > 3*1.5/associated_points.size()) logical = 0;
1268  if( layer == 4 && point_energy/sintheta/point_energy_sum > 4*1.5/associated_points.size()) logical = 0;
1269  if(E_shower>0.3) logical = 0;
1270  if(associated_points.size()<4) logical = 0;
1271  if(layer==1 && point_energy/sintheta >0.4) logical=0;
1272  if(layer==2 && point_energy/sintheta < 0.2) logical=0;
1273  if(layer==2 && point_energy/sintheta > 0.6) logical=0;
1274  if(layer==3 && point_energy/sintheta < 0.4) logical=0;
1275  if(layer==3 && point_energy/sintheta > 0.8) logical=0;
1276  if(layer==4 && point_energy/sintheta < 0.6) logical=0;
1277  if(layer==4 && point_energy/sintheta > 1.0) logical=0;
1278  if (layer == 1 && sector == 1) channel_per_module = 0;
1279  if (layer == 1 && sector == 2) channel_per_module = 1;
1280  if (layer == 1 && sector == 3) channel_per_module = 2;
1281  if (layer == 1 && sector == 4) channel_per_module = 3;
1282  if (layer == 2 && sector == 1) channel_per_module = 4;
1283  if (layer == 2 && sector == 2) channel_per_module = 5;
1284  if (layer == 2 && sector == 3) channel_per_module = 6;
1285  if (layer == 2 && sector == 4) channel_per_module = 7;
1286  if (layer == 3 && sector == 1) channel_per_module = 8;
1287  if (layer == 3 && sector == 2) channel_per_module = 9;
1288  if (layer == 3 && sector == 3) channel_per_module = 10;
1289  if (layer == 3 && sector == 4) channel_per_module = 11;
1290  if (layer == 4 && sector == 1) channel_per_module = 12;
1291  if (layer == 4 && sector == 2) channel_per_module = 13;
1292  if (layer == 4 && sector == 3) channel_per_module = 14;
1293  if (layer == 4 && sector == 4) channel_per_module = 15;
1294  int channel = channel_per_module + (module-1)*16;
1295  if(charge < 0 && logical == 1)
1296  {
1297 
1298  //if(layer==1) Layer1_Energy_v2->Fill(point_energy/sintheta);
1299  //if(layer==2) Layer2_Energy_v2->Fill(point_energy/sintheta);
1300  //if(layer==3) Layer3_Energy_v2->Fill(point_energy/sintheta);
1301  //if(layer==4) Layer4_Energy_v2->Fill(point_energy/sintheta);
1302  //cout << " charge = " << locChargedTrackHypothesis->charge() << " point energy/sintheta= " << point_energy/sintheta << " 'shower energy' = " << E1 << " energy sum = " << point_energy_sum << " point energy = " << point_energy << " theta = " << theta << " channel = " << channel << " module = " << module << " layer = " << layer << " point size = " << associated_points.size() << endl;
1303  }
1304  //cout << " charge = " << locChargedTrackHypothesis->charge() << " point energy/sintheta= " << point_energy/sintheta1 << " 'shower energy' = " << E1 << " shower momentum = " << p1.M() << " channel = " << channel << " module = " << module1 << " layer = " << layer1 << " point size = " << associated_points.size() << endl;
1305 
1306 
1307 
1308 
1309 
1310 
1311  }
1312  }
1313  //}
1314 
1315  }
1316  }
1317  }
1318 
1319 
1320 */
1321 
1322 
1323  for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1324  {
1325  const DBCALShower* locBCALShower = locBCALShowers[loc_i];
1326  if(find(matchedShowersneg.begin(), matchedShowersneg.end(), locBCALShower) == matchedShowersneg.end()) continue;
1327  energy_shower = locBCALShower->E;
1328  energy_raw_shower = locBCALShower->E_raw;
1329  energy_raw_shower = locBCALShower->E_raw;
1330  track_momentum = p;
1331  double shower_x = locBCALShower->x;
1332  double shower_y = locBCALShower->y;
1333 
1334  r = TMath::Sqrt(TMath::Power(shower_x,2)+TMath::Power(shower_y,2));
1335  z = locBCALShower->z;
1336  //phi = locBCALShower->phi;
1337  //theta = locBCALShower->theta;
1338  vector<const DBCALCluster*> associated_clusters;
1339  locBCALShower->Get(associated_clusters);
1340 
1341  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
1342 
1343  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1344  {
1345 
1346  const DBCALCluster* a_cluster = associated_clusters[loc_j];
1347  vector<const DBCALPoint*> associated_points;
1348  a_cluster->Get(associated_points);
1349  layer1_energysum = 0.0;
1350  layer2_energysum = 0.0;
1351  layer3_energysum = 0.0;
1352  layer4_energysum = 0.0;
1353  double layer1_energysum_inter = 0.0;
1354  double layer2_energysum_inter = 0.0;
1355  double layer3_energysum_inter = 0.0;
1356  double layer4_energysum_inter = 0.0;
1357  phi = associated_clusters[loc_j]->phi();
1358  theta = associated_clusters[loc_j]->theta();
1359 
1360  for(unsigned int loc_k = 0 ; loc_k < associated_points.size(); loc_k++)
1361  {
1362  energy_point = associated_points[loc_k]->E();
1363  module = associated_points[loc_k]->module();
1364  layer = associated_points[loc_k]->layer();
1365  sector = associated_points[loc_k]->sector();
1366  //theta = associated_points[loc_k]->theta();
1367  // z = associated_points[loc_k]->z();
1368  // phi = associated_points[loc_k]->phi();
1369  // r = associated_points[loc_k]->r();
1370  double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
1371  double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
1372  double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
1373  double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
1374  double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
1375  L1_pathlength = second_pos - first_pos;
1376  L2_pathlength = third_pos - second_pos;
1377  L3_pathlength = fourth_pos - third_pos;
1378  L4_pathlength = fifth_pos - fourth_pos;
1379 
1380  if(layer==1) layer1_energysum_inter += energy_point;
1381  if(layer==2) layer2_energysum_inter += energy_point;
1382  if(layer==3) layer3_energysum_inter += energy_point;
1383  if(layer==4) layer4_energysum_inter += energy_point;
1384  if(loc_k+1 == associated_points.size()){
1385  layer1_energysum = layer1_energysum_inter;
1386  layer2_energysum = layer2_energysum_inter;
1387  layer3_energysum = layer3_energysum_inter;
1388  layer4_energysum = layer4_energysum_inter;
1389  }
1390  //cout << " shower E = " << energy_shower << " raw shower E = " << energy_raw_shower << " layer1_energysum = " << layer1_energysum << " layer 2 energy sum = " << layer2_energysum << " layer 3 energy sum = " << layer3_energysum << " layer4 energy sum = " << layer4_energysum << " p = " << track_momentum << endl;
1391  //if(layer1_energysum > 0.0 || layer2_energysum > 0.0 || layer3_energysum > 0.0 || layer4_energysum > 0.0 ) cout << " shower E = " << energy_shower << " raw shower E = " << energy_raw_shower << " layer1_energysum = " << layer1_energysum << " layer 2 energy sum = " << layer2_energysum << " layer 3 energy sum = " << layer3_energysum << " layer4 energy sum = " << layer4_energysum << " p = " << track_momentum << " REPEATER " << endl;
1392  int channel_per_module;
1393  if (layer == 1 && sector == 1) channel_per_module = 0;
1394  if (layer == 1 && sector == 2) channel_per_module = 1;
1395  if (layer == 1 && sector == 3) channel_per_module = 2;
1396  if (layer == 1 && sector == 4) channel_per_module = 3;
1397  if (layer == 2 && sector == 1) channel_per_module = 4;
1398  if (layer == 2 && sector == 2) channel_per_module = 5;
1399  if (layer == 2 && sector == 3) channel_per_module = 6;
1400  if (layer == 2 && sector == 4) channel_per_module = 7;
1401  if (layer == 3 && sector == 1) channel_per_module = 8;
1402  if (layer == 3 && sector == 2) channel_per_module = 9;
1403  if (layer == 3 && sector == 3) channel_per_module = 10;
1404  if (layer == 3 && sector == 4) channel_per_module = 11;
1405  if (layer == 4 && sector == 1) channel_per_module = 12;
1406  if (layer == 4 && sector == 2) channel_per_module = 13;
1407  if (layer == 4 && sector == 3) channel_per_module = 14;
1408  if (layer == 4 && sector == 4) channel_per_module = 15;
1409  channel = channel_per_module + (module-1)*16;
1410  //double sintheta = sin(theta);
1411 
1412  // cout << " energy = " << energy << " charge = " << charge << " module = " << module << " layer = " << layer << " point size = " << associated_points.size() << " event num = " << eventnum << " theta = " << theta << " z = " << z << " r = " << r << " vertexZ = " << kinfitVertexZ << " theta wrt vertex = " << theta_wrt_vertex << endl;
1413 
1414 
1415  if(layer1_energysum > 0.0 || layer2_energysum > 0.0 || layer3_energysum > 0.0 || layer4_energysum > 0.0 ) BCALPoint_Charged_neg->Fill();
1416  }
1417  }
1418  }
1419 
1420 
1421 
1422  for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1423  {
1424  const DBCALShower* locBCALShower = locBCALShowers[loc_i];
1425  if(find(matchedShowerspos.begin(), matchedShowerspos.end(), locBCALShower) == matchedShowerspos.end()) continue;
1426  energy_shower = locBCALShower->E;
1427  energy_raw_shower = locBCALShower->E_raw;
1428  double shower_x = locBCALShower->x;
1429  double shower_y = locBCALShower->y;
1430  track_momentum = p;
1431  z = locBCALShower->z;
1432  //phi = locBCALShower->phi;
1433  //theta = locBCALShower->theta;
1434  r = TMath::Sqrt(TMath::Power(shower_x,2)+TMath::Power(shower_y,2));
1435  vector<const DBCALCluster*> associated_clusters;
1436  locBCALShower->Get(associated_clusters);
1437 
1438  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
1439 
1440  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1441  {
1442 
1443  const DBCALCluster* a_cluster = associated_clusters[loc_j];
1444  vector<const DBCALPoint*> associated_points;
1445  a_cluster->Get(associated_points);
1446 
1447  phi = associated_clusters[loc_j]->phi();
1448  theta = associated_clusters[loc_j]->theta();
1449 
1450  layer1_energysum = 0.0;
1451  layer2_energysum = 0.0;
1452  layer3_energysum = 0.0;
1453  layer4_energysum = 0.0;
1454  double layer1_energysum_inter = 0.0;
1455  double layer2_energysum_inter = 0.0;
1456  double layer3_energysum_inter = 0.0;
1457  double layer4_energysum_inter = 0.0;
1458 
1459  for(unsigned int loc_k = 0 ; loc_k < associated_points.size(); loc_k++)
1460  {
1461  energy_point = associated_points[loc_k]->E();
1462  module = associated_points[loc_k]->module();
1463  layer = associated_points[loc_k]->layer();
1464  sector = associated_points[loc_k]->sector();
1465  //theta = associated_points[loc_k]->theta();
1466  // z = associated_points[loc_k]->z();
1467  // phi = associated_points[loc_k]->phi();
1468  // r = associated_points[loc_k]->r();
1469  double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
1470  double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
1471  double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
1472  double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
1473  double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
1474  L1_pathlength = second_pos - first_pos;
1475  L2_pathlength = third_pos - second_pos;
1476  L3_pathlength = fourth_pos - third_pos;
1477  L4_pathlength = fifth_pos - fourth_pos;
1478 
1479  int channel_per_module;
1480  if (layer == 1 && sector == 1) channel_per_module = 0;
1481  if (layer == 1 && sector == 2) channel_per_module = 1;
1482  if (layer == 1 && sector == 3) channel_per_module = 2;
1483  if (layer == 1 && sector == 4) channel_per_module = 3;
1484  if (layer == 2 && sector == 1) channel_per_module = 4;
1485  if (layer == 2 && sector == 2) channel_per_module = 5;
1486  if (layer == 2 && sector == 3) channel_per_module = 6;
1487  if (layer == 2 && sector == 4) channel_per_module = 7;
1488  if (layer == 3 && sector == 1) channel_per_module = 8;
1489  if (layer == 3 && sector == 2) channel_per_module = 9;
1490  if (layer == 3 && sector == 3) channel_per_module = 10;
1491  if (layer == 3 && sector == 4) channel_per_module = 11;
1492  if (layer == 4 && sector == 1) channel_per_module = 12;
1493  if (layer == 4 && sector == 2) channel_per_module = 13;
1494  if (layer == 4 && sector == 3) channel_per_module = 14;
1495  if (layer == 4 && sector == 4) channel_per_module = 15;
1496  channel = channel_per_module + (module-1)*16;
1497  //double sintheta = sin(theta);
1498 
1499  if(layer==1) layer1_energysum_inter += energy_point;
1500  if(layer==2) layer2_energysum_inter += energy_point;
1501  if(layer==3) layer3_energysum_inter += energy_point;
1502  if(layer==4) layer4_energysum_inter += energy_point;
1503  if(loc_k+1 == associated_points.size()){
1504  layer1_energysum = layer1_energysum_inter;
1505  layer2_energysum = layer2_energysum_inter;
1506  layer3_energysum = layer3_energysum_inter;
1507  layer4_energysum = layer4_energysum_inter;
1508  }
1509  // cout << " energy = " << energy << " charge = " << charge << " module = " << module << " layer = " << layer << " point size = " << associated_points.size() << " event num = " << eventnum << " theta = " << theta << " z = " << z << " r = " << r << " vertexZ = " << kinfitVertexZ << " theta wrt vertex = " << theta_wrt_vertex << endl;
1510 
1511 
1512  if(layer1_energysum > 0.0 || layer2_energysum > 0.0 || layer3_energysum > 0.0 || layer4_energysum > 0.0 ) BCALPoint_Charged_pos->Fill();
1513  }
1514  }
1515  }
1516 
1517 
1518 
1519 
1520 //LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL LOOKING FOR E OVER P LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
1521 
1522 
1523 
1524 
1525  for (unsigned int i=0; i < locTrackTimeBased.size() ; i++)
1526  {
1527  if (find(matchedTracks.begin(), matchedTracks.end(), locTrackTimeBased[i]) == matchedTracks.end()) continue;
1528  double charge = locTrackTimeBased[i]->rt->q;
1529  double p = locTrackTimeBased[i]->momentum().Mag();
1530  //cout << " p = " << p << endl;
1531 
1532  for(unsigned int i=0; i<locBCALShowers.size(); i++)
1533  {
1534  const DBCALShower *s1 = locBCALShowers[i];
1535  if (find(matchedShowers.begin(), matchedShowers.end(),s1) == matchedShowers.end()) continue;
1536  double E_MatchedShower = s1->E;
1537  //cout << " shower E = " << E_MatchedShower << " track p = " << p << endl;
1538  //if(charge==1.0)Eoverp_plus->Fill(E_MatchedShower/p);
1539  //if(charge==-1.0)Eoverp_minus->Fill(E_MatchedShower/p);
1540  //if(charge==1.0)Evsp_plus->Fill(E_MatchedShower,p);
1541  //if(charge==-1.0)Evsp_minus->Fill(E_MatchedShower,p);
1542 
1543  }
1544 
1545  }
1546 
1547 
1548 
1549 
1550 
1551 
1552 
1553 
1554 
1555 
1556 
1557 
1558 // MMMMMMMMMMMMMMMMMMMMMMMMMMMM MATHCING TEST MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
1559 /*
1560 double locDistance=99999;
1561 const DBCALShower* locBCALShowerz;
1562 for(size_t i = 0 ; i < locBCALShowers.size() ; ++i)
1563 {
1564 locDetectorMatches->Get_DistanceToNearestTrack(locBCALShowerz, locDistance);
1565 cout << " shower track doca = " << locDistance << endl;
1566 }
1567 */
1568 
1569 
1570 
1571  //cons locTrackTimeBased[j]->DReferenceTrajectory;
1572  //if(rt == NULL) continue;
1573 
1574  // DVector3 mypos(0.0,0.0,0.0);
1575 
1576  for(size_t loc_i = 0; loc_i < locBCALShowers.size(); ++loc_i)
1577  {
1578  double x = locBCALShowers[loc_i]->x;
1579  double y = locBCALShowers[loc_i]->y;
1580  double z = locBCALShowers[loc_i]->z;
1581  double E = locBCALShowers[loc_i]->E;
1582  DVector3 pos_bcal(x,y,z);
1583  double R = pos_bcal.Perp();
1584  double R1 = sqrt((x-kinfitVertexX)*(x-kinfitVertexX)+(y-kinfitVertexY)*(y-kinfitVertexY)+(z-kinfitVertexZ)*(z-kinfitVertexZ));
1585  double phi = pos_bcal.Phi();
1586  for(size_t j = 0 ; j < locTrackTimeBased.size(); ++j)
1587  {
1588 
1589  locTrackTimeBased[j]->rt->GetIntersectionWithRadius(R, mypos);
1590  //double dPhi = sqrt((mypos.Phi() - pos_bcal.Phi())*(mypos.Phi() - pos_bcal().Phi()));
1591  double dPhi = TMath::Abs(mypos.Phi()-pos_bcal.Phi());
1592  double dZ = TMath::Abs(mypos.Z() - z);
1593  TLorentzVector p1(E*(x-kinfitVertexX)/R, E*(y-kinfitVertexY)/R1, E*(z-kinfitVertexZ)/R1, E);
1594  // if(dZ < 10.0 && dPhi < 1.0) EoverP->Fill(E/p1.M());
1595  // if(R==mypos.Perp())cout << " radius = " << mypos.Perp() << " R = " << R << " phi track = " << mypos.Phi() << " bcal phi = " << pos_bcal.Phi() << " dPhi = " << dPhi << " dZ = " << dZ << endl;
1596  //cout << " pos_bcal.perp = " << pos_bcal.Perp() << " pos_bcal.phi = " << pos_bcal.Phi() << endl;
1597  }
1598 }
1599 
1600 //double x = locBCALShowers->x;
1601 //double y = locBCALShowers->y;
1602 //double z = locBCALShowers->z;
1603 //DVector3 bcal_pos(x, y, z);
1604 //DVector3 proj_pos = rt->GetLastDOCAPoint();
1605 //cout << " z pos = " << proj_pos.Perp() << endl;
1606 
1607 
1608 
1609 for(unsigned int j = 0 ; j < locTrackTimeBased.size(); ++j)
1610  {
1611  for(unsigned int i = 0 ; i < locFCALShowers.size(); ++i)
1612  {
1613  double E = locFCALShowers[i]->getEnergy();
1614  double x = locFCALShowers[i]->getPosition().X();
1615  double y = locFCALShowers[i]->getPosition().Y();
1616  double z = locFCALShowers[i]->getPosition().Z();
1617  DVector3 origin(0.0,0.0,z);
1618  DVector3 norm(0.0,0.0,-1.0);
1619  DVector3 mypos2(0.0,0.0,0.0);
1620  DVector3 mymom(0.0,0.0,0.0);
1621  DVector3 pos_fcal(x,y,z);
1622  double myorigin = pos_fcal.Perp();
1623 
1624 
1625  locTrackTimeBased[j]->rt->GetIntersectionWithPlane(pos_fcal,norm,mypos2,mymom);
1626  double pmag = mymom.Mag();
1627  //cout << " p mag = " << pmag << " my pos perp = " << mypos2.Perp() << " my pos x = " << mypos2.x() << " my pos y = " << mypos2.y() << "my pos z = " << mypos2.z() << " pos fcal x = " << pos_fcal.x() << " pos fcal y = " << pos_fcal.y() << " pos fcal z = " << pos_fcal.z() << " pos fcal perp = " << pos_fcal.Perp() << " energy = " << E << " shower size = " << locFCALShowers.size() << " track size = " << locTrackTimeBased.size () << endl;
1628 
1629 
1630  double dX = TMath::Abs(mypos2.X() - x);
1631  double dY = TMath::Abs(mypos2.Y() - y);
1632  double dZ = TMath::Abs(mypos2.Z() - z);
1633  double dR = TMath::Abs(mypos2.Perp() - pos_fcal.Perp());
1634  //if(pos_fcal.z() == mypos2.z() && pos_fcal.Perp() > 20 && dR < 3.0) FCAL_Evsp->Fill(pmag,E);
1635  //if(pos_fcal.z() == mypos2.z() && pos_fcal.Perp() > 20 && dR < 3.0) FCAL_Eoverpvsp->Fill(pmag,E/pmag);
1636  //if(pos_fcal.z() == mypos2.z() && pos_fcal.Perp() > 20 && dR < 3.0 && pmag > 1.3 && pos_fcal.Perp() > 60.0) FCAL_Eoverp_nocuts->Fill(E/pmag);
1637  //if(pos_fcal.z() == mypos2.z() && dR < 3.0 && pmag > 1.3 && pos_fcal.Perp() > 20.0 && pos_fcal.Perp() < 60.0) FCAL_Eoverp_cuts->Fill(E/pmag);
1638  //cout << " E = " << E << " p mag = " << pmag << " my pos 2 perp = " << mypos2.Perp() << " dx = " << dX << " dy = " << dY << " eventnum = " << eventnum << endl;
1639  }
1640 
1641 }
1642 
1643 
1644 
1645 
1646 
1647 //vector<double> evsp_vec, p_vec;
1648 for(size_t j = 0 ; j < locTrackTimeBased.size(); ++j)
1649 {
1650  for(size_t loc_i = 0; loc_i < locFCALClusters.size(); ++loc_i)
1651 {
1652 //double x = locFCALShowers[loc_i]->getPosition().X();
1653 //double y = locFCALShowers[loc_i]->getPosition().Y();
1654 //double z = locFCALShowers[loc_i]->getPosition().Z();
1655 DVector3 fcalpos=locFCALClusters[loc_i]->getCentroid();
1656 //cout << " cluster x = " << fcalpos.X() << " cluster y = " << fcalpos.Y() << endl;
1657  // fcalpos.SetZ( fcalpos.Z() + m_targetZ );
1658 //DVector3 fcalpos(x,y,z);
1659 DVector3 norm(0.0,0.0,-1);
1660 DVector3 pos,mom;
1661 
1662 //locTrackTimeBased[j]->rt->GetIntersectionWithPlane(origin,norm,pos);
1663 if(locTrackTimeBased[j]->rt->GetIntersectionWithPlane(fcalpos,norm,pos,mom)==NOERROR)
1664 {
1665 double diffX = TMath::Abs(fcalpos.X() - pos.X());
1666 double diffY = TMath::Abs(fcalpos.Y() - pos.Y());
1667 double track_mom = TMath::Sqrt( TMath::Power(mom.X(),2) +
1668 TMath::Power(mom.Y(),2) + TMath::Power(mom.Z(),2));
1669 double cluster_energy = locFCALClusters[loc_i]->getEnergy();
1670 double EvsP = cluster_energy/track_mom;
1671 double cluster_radius = TMath::Sqrt( TMath::Power(fcalpos.X(),2) +
1672 TMath::Power(fcalpos.Y(),2));
1673 //cout << " cluster E = " << cluster_energy << " cluster radius = " << cluster_radius << " track momentum = " << track_mom << " diff x = " << diffX << " diff y = " << diffY << " event num = " << eventnum << endl;
1674 //m_evsp->Fill( EvsP );
1675 //if(diffX < 10 && diffY < 10) cout << " clusrer E = " << cluster_energy << " track p = " << track_mom << " cluster radius = " << cluster_radius << endl;
1676 if(diffX < 10 && diffY < 10){
1677  FCAL_Eoverp_nocuts->Fill( EvsP );
1678  FCAL_Evsp->Fill(track_mom,cluster_energy);
1679  FCAL_Eoverpvsp->Fill(track_mom,cluster_energy/track_mom);
1680 }
1681 //std::cout<<fcalpos.X()<<"\t"<<fcalpos.Y()<<"\t"<<fcalpos.Z()<<endl;
1682 //std::cout<<pos.X()<<"\t"<<pos.Y()<<"\t"<<pos.Z()<<endl;
1683 //std::cout<<endl;
1684 if (diffX < 5 && diffY < 5 && track_mom >1.3 && cluster_radius > 20.0 && cluster_radius < 60.0)
1685 {
1686 //m_posX->Fill( fcalpos.X() - pos.X() );
1687 //m_posY->Fill( fcalpos.Y() - pos.Y() );
1688 //m_energy->Fill(cluster_energy);
1689 //m_mom->Fill(track_mom);
1690 //m_evspcut->Fill( EvsP);
1691 FCAL_Eoverp_cuts->Fill( EvsP );
1692 //m_radius->Fill(cluster_radius);
1693 //evsp_vec.push_back(EvsP);
1694 //p_vec.push_back(track_mom);
1695 //m_ebypvsp->Fill( EvsP, track_mom );
1696 //std::cout<<"Real Tracks:"<<endl;
1697 //std::cout<<fcalpos.X()<<"\t"<<fcalpos.Y()<<"\t"<<fcalpos.Z()<<endl;
1698 //std::cout<<pos.X()<<"\t"<<pos.Y()<<"\t"<<pos.Z()<<endl;
1699 //std::cout<<endl;
1700 }
1701  }
1702 }
1703 }
1704 
1705 
1706 
1707 
1708 
1709 
1710 
1711 
1712 //WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWw Fill ROOT Tree Test WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1713 
1714 /*
1715 for(size_t i = 0 ; i < locChargedTracks.size(); ++i)
1716  {
1717 
1718  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[i]->Get_BestFOM();
1719  // cout << " Best FOM = " << locChargedTracks[i]->Get_BestFOM() << endl;
1720  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
1721  const DShowerMatchParams& locBCALShowerMatchParams = locChargedTrackHypothesis->dBCALShowerMatchParams;
1722 
1723  // if(locBCALShowerMatchParams.dTrackTimeBased != NULL)
1724  // {
1725  double d_min = locBCALShowerMatchParams.dDOCAToShower;
1726  //if(d_min != 0.0) cout << " d min = " << d_min << endl;
1727  const DBCALShower* locBCALShower = dynamic_cast <const DBCALShower*>(locBCALShowerMatchParams.dShowerObject);
1728  if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShower)){
1729  BCALShowerTrack_Energy->Fill(locBCALShower->E);
1730  charge = locChargedTrackHypothesis->charge();
1731  //cout << " shower track E = " << BCALShowerTrack_Energy << endl;
1732  //for(unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1733  //{
1734  //const DBCALShower* s1 = locBCALShowers;
1735  //const DBCALShower* a_shower = locBCALShowers;
1736 
1737  //double shower_energy = locBCALShower->E;
1738  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
1739  vector<const DBCALCluster*> associated_clusters;
1740  locBCALShower->Get(associated_clusters);
1741 
1742  if(associated_clusters.size() == 0) cout << " assoc cluster is empty " << endl;
1743 
1744  for(unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1745  {
1746 
1747  const DBCALCluster* a_cluster = associated_clusters[loc_j];
1748  vector<const DBCALPoint*> associated_points;
1749  a_cluster->Get(associated_points);
1750 
1751  if(associated_points.size() == 0) cout << " assoc points is empty " << endl;
1752 
1753  for(unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
1754  {
1755  //cout << " charged point size = " << associated_points.size() << endl;
1756  energy = associated_points[loc_k]->E();
1757  module = associated_points[loc_k]->module();
1758  layer = associated_points[loc_k]->layer();
1759  sector = associated_points[loc_k]->sector();
1760  theta = associated_points[loc_k]->theta();
1761  z = associated_points[loc_k]->z();
1762  phi = associated_points[loc_k]->phi();
1763  r = associated_points[loc_k]->r();
1764  theta_wrt_vertex =atan(r/(z-kinfitVertexZ));
1765 
1766  int channel_per_module;
1767  if (layer == 1 && sector == 1) channel_per_module = 0;
1768  if (layer == 1 && sector == 2) channel_per_module = 1;
1769  if (layer == 1 && sector == 3) channel_per_module = 2;
1770  if (layer == 1 && sector == 4) channel_per_module = 3;
1771  if (layer == 2 && sector == 1) channel_per_module = 4;
1772  if (layer == 2 && sector == 2) channel_per_module = 5;
1773  if (layer == 2 && sector == 3) channel_per_module = 6;
1774  if (layer == 2 && sector == 4) channel_per_module = 7;
1775  if (layer == 3 && sector == 1) channel_per_module = 8;
1776  if (layer == 3 && sector == 2) channel_per_module = 9;
1777  if (layer == 3 && sector == 3) channel_per_module = 10;
1778  if (layer == 3 && sector == 4) channel_per_module = 11;
1779  if (layer == 4 && sector == 1) channel_per_module = 12;
1780  if (layer == 4 && sector == 2) channel_per_module = 13;
1781  if (layer == 4 && sector == 3) channel_per_module = 14;
1782  if (layer == 4 && sector == 4) channel_per_module = 15;
1783  channel = channel_per_module + (module-1)*16;
1784  //double sintheta = sin(theta);
1785 
1786  // cout << " energy = " << energy << " charge = " << charge << " module = " << module << " layer = " << layer << " point size = " << associated_points.size() << " event num = " << eventnum << " theta = " << theta << " z = " << z << " r = " << r << " vertexZ = " << kinfitVertexZ << " theta wrt vertex = " << theta_wrt_vertex << endl;
1787 
1788 
1789  BCALPoint_Charged->Fill();
1790 
1791  }
1792  }
1793  //}
1794 
1795  }
1796  // }
1797  }
1798 
1799 
1800 */
1801 
1802 
1803 
1804  for(unsigned int i=0; i<locBCALShowers.size(); i++)
1805  {
1806  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
1807  // continue;
1808  if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end()) continue;
1809  const DBCALShower *s1 = locBCALShowers[i];
1810  vector<const DBCALCluster*> associated_clusters1;
1811 
1812  s1->Get(associated_clusters1);
1813  double shower_energy = s1->E;
1814  double raw_shower_energy = s1->E_raw;
1815  //cout << " shower E = " << shower_energy << endl;
1816  //cout << " shower raw E = " << raw_shower_energy << endl;
1817  for(unsigned int loc_j = 0; loc_j < associated_clusters1.size(); loc_j++)
1818  {
1819  double cluster_energy = associated_clusters1[loc_j]->E();
1820  //cout << " cluster E = " << cluster_energy << endl;
1821  const DBCALCluster* a_cluster1 = associated_clusters1[loc_j];
1822  vector<const DBCALPoint*> associated_points1;
1823  a_cluster1->Get(associated_points1);
1824  for(unsigned int loc_k = 0; loc_k < associated_points1.size(); loc_k++)
1825  {
1826  module = associated_points1[loc_k]->module();
1827  layer = associated_points1[loc_k]->layer();
1828  sector = associated_points1[loc_k]->sector();
1829  int channel_per_module;
1830  if (layer == 1 && sector == 1) channel_per_module = 0;
1831  if (layer == 1 && sector == 2) channel_per_module = 1;
1832  if (layer == 1 && sector == 3) channel_per_module = 2;
1833  if (layer == 1 && sector == 4) channel_per_module = 3;
1834  if (layer == 2 && sector == 1) channel_per_module = 4;
1835  if (layer == 2 && sector == 2) channel_per_module = 5;
1836  if (layer == 2 && sector == 3) channel_per_module = 6;
1837  if (layer == 2 && sector == 4) channel_per_module = 7;
1838  if (layer == 3 && sector == 1) channel_per_module = 8;
1839  if (layer == 3 && sector == 2) channel_per_module = 9;
1840  if (layer == 3 && sector == 3) channel_per_module = 10;
1841  if (layer == 3 && sector == 4) channel_per_module = 11;
1842  if (layer == 4 && sector == 1) channel_per_module = 12;
1843  if (layer == 4 && sector == 2) channel_per_module = 13;
1844  if (layer == 4 && sector == 3) channel_per_module = 14;
1845  if (layer == 4 && sector == 4) channel_per_module = 15;
1846  channel = channel_per_module + (module-1)*16;
1847  energy_point = associated_points1[loc_k]->E();
1848  theta = associated_points1[loc_k]->theta();
1849  // cout << " point E = " << energy << " point size = " << associated_points1.size() << " shower size = " << locBCALShowers.size() << " channel = " << channel << " event num = " << eventnum << endl;
1850  BCALPoint_Neutral->Fill();
1851 
1852  }
1853  }
1854  }
1855 
1856 
1857 
1858 
1859 
1860 
1861 /*
1862  vector<const DBCALPoint* > bcalpoints;
1863  locEventLoop->Get(bcalpoints);
1864  for(unsigned int i = 0; i < bcalpoints.size(); i++)
1865  {
1866  if(locDetectorMatches->Get_IsMatchedToTrack(bcalpoints[i]))
1867  continue;
1868 
1869  theta = bcalpoints[i]->theta();
1870  energy = bcalpoints[i]->E();
1871  module = bcalpoints[i]->module();
1872  layer = bcalpoints[i]->layer();
1873  sector = bcalpoints[i]->sector();
1874  int channel_per_module;
1875  if (layer == 1 && sector == 1) channel_per_module = 0;
1876  if (layer == 1 && sector == 2) channel_per_module = 1;
1877  if (layer == 1 && sector == 3) channel_per_module = 2;
1878  if (layer == 1 && sector == 4) channel_per_module = 3;
1879  if (layer == 2 && sector == 1) channel_per_module = 4;
1880  if (layer == 2 && sector == 2) channel_per_module = 5;
1881  if (layer == 2 && sector == 3) channel_per_module = 6;
1882  if (layer == 2 && sector == 4) channel_per_module = 7;
1883  if (layer == 3 && sector == 1) channel_per_module = 8;
1884  if (layer == 3 && sector == 2) channel_per_module = 9;
1885  if (layer == 3 && sector == 3) channel_per_module = 10;
1886  if (layer == 3 && sector == 4) channel_per_module = 11;
1887  if (layer == 4 && sector == 1) channel_per_module = 12;
1888  if (layer == 4 && sector == 2) channel_per_module = 13;
1889  if (layer == 4 && sector == 3) channel_per_module = 14;
1890  if (layer == 4 && sector == 4) channel_per_module = 15;
1891  channel = channel_per_module + (module-1)*16;
1892  BCALPoint_Neutral->Fill();
1893 
1894  }
1895 
1896 
1897 
1898 
1899 
1900 
1901 */
1902 
1903 
1904 
1905 
1906 
1907 
1908 
1909 
1910 
1911 
1912 //OoOoOoOoOoOoOoOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO CALCULATING THE PI0 INVARIANT MASS PORTION O0O0O0O0O0O0O0O0OOooooooooooOooOoOoooooOOooooOOo
1913 
1914 
1915 
1916 
1917 
1918  //2 gamma inv mass
1919 
1920 
1921 
1922 
1923  for(unsigned int i = 0; i < matchedShowers.size(); i++)
1924  {
1925  double E = matchedShowers[i]->E;
1926  //cout << " E = " << E << endl;
1927  const DBCALShower *ms = matchedShowers[i];
1928  vector<const DBCALCluster*> associated_clusters;
1929  ms->Get(associated_clusters);
1930  for(unsigned int j = 0 ; j < associated_clusters.size(); j++)
1931  {
1932  double cluster_E = associated_clusters[j]->E();
1933  //cout << " cluster E = " << cluster_E << endl;
1934  }
1935 
1936  }
1937 
1938 
1939 
1940  for(unsigned int i=0; i<locBCALShowers.size(); i++)
1941  {
1942  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
1943  // continue;
1944  double E_before = locBCALShowers[i]->E;
1945  BCALShowers_per_event = locBCALShowers.size();
1946  //cout << " E before = " << E_before << endl;
1947  if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end()) continue;
1948  const DBCALShower *s1 = locBCALShowers[i];
1949  vector<const DBCALCluster*> associated_clusters1;
1950  s1->Get(associated_clusters1);
1951  double dx1 = s1->x - kinfitVertexX;
1952  double dy1 = s1->y - kinfitVertexY;
1953  double dz1 = s1->z - kinfitVertexZ;
1954  double dt1 = s1->t;
1955  t1 = dt1;
1956  z1 = s1->z;
1957  x1 = s1->x;
1958  y1 = s1->y;
1959 
1960  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
1961  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
1962  E1 = s1->E;
1963  E1_raw = s1->E_raw;
1964 
1965  //cout << " E after = " << E1 << endl;
1966  double ECalc = s1->E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
1967  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
1968  TLorentzVector p1_raw(E1_raw*dx1/R1, E1_raw*dy1/R1, E1_raw*dz1/R1, E1_raw);
1969  double thetadeg1, thetarad1, phideg1, phirad1;
1970  thetadeg1 = p1.Theta()*180.0/TMath::Pi();
1971  thetarad1 = p1.Theta();
1972  phideg1 = p1.Phi()*180.0/TMath::Pi();
1973  phirad1 = p1.Phi();
1974  phi1 = phirad1;
1975  p1_mag = p1.M();
1976  p1_raw_mag = p1_raw.M();
1977  vertexZ = kinfitVertexZ;
1978  vertexX = kinfitVertexX;
1979  vertexY = kinfitVertexY;
1980  TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
1981 
1982 
1983  for(unsigned int j=i+1; j<locBCALShowers.size(); j++){
1984  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[j]))
1985  // continue;
1986  const DBCALShower *s2 = locBCALShowers[j];
1987  //cout << " energy before matching = " << s2->E << endl;
1988  if (find(matchedShowers.begin(), matchedShowers.end(),s2) != matchedShowers.end()) continue;
1989  //cout << " energy after matching = " << s2->E << endl;
1990  vector<const DBCALCluster*> associated_clusters2;
1991  s2->Get(associated_clusters2);
1992  double dx2 = s2->x - kinfitVertexX;
1993  double dy2 = s2->y - kinfitVertexY;
1994  double dz2 = s2->z - kinfitVertexZ; // shift to coordinate relative to center of target
1995  double dt2 = s2->t;
1996  t2=dt2;
1997  z2 = s2->z;
1998  x2 = s2->x;
1999  y2 = s2->y;
2000  double deltaT = TMath::Abs(dt1-dt2);
2001  Time_Diff->Fill(deltaT);
2002  double R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2003  //double path2 = sqrt((dx2-kinfitVertexX)*(dx2-kinfitVertexX)+(dy2-kinfitVertexY)*(dy2-kinfitVertexY)+(dz2)*(dz2));
2004  E2 = s2->E;
2005  E2_raw = s2->E_raw;
2006  double ECalc = s2->E*(1.106+(dz2+65.0-208.4)*(dz2+65.0-208.4)*6.851E-6);
2007  TLorentzVector p2(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
2008  TLorentzVector p2_raw(E2_raw*dx2/R2, E2_raw*dy2/R2, E2_raw*dz2/R2, E2_raw);
2009  p2_mag = p2.M();
2010  p2_raw_mag = p2_raw.M();
2011  TLorentzVector p2Calc(ECalc*dx2/R2, ECalc*dy2/R2, ECalc*dz2/R2, ECalc);
2012  double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
2013  thetarad2 = p2.Theta();
2014  phirad2 = p2.Phi();
2015  phi2 = phirad2;
2016  thetadeg2 = p2.Theta()*180.0/TMath::Pi();
2017  phideg2 = p2.Phi()*180.0/TMath::Pi();
2018  cospsi = sin(thetarad1)*sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
2019  psi=acos(cospsi)*180/TMath::Pi();
2020  TLorentzVector ptot = p1+p2;
2021  TLorentzVector ptot_raw = p1_raw + p2_raw ;
2022  inv_mass = ptot.M();
2023  inv_mass_raw = ptot_raw.M();
2024 
2025  Theta_Hist_Both_BCAL->Fill(thetadeg1);
2026  Theta_Hist_Both_BCAL->Fill(thetadeg2);
2027  Phi_Hist_Both_BCAL->Fill(phideg2);
2028  Phi_Hist_Both_BCAL->Fill(phideg2);
2029  Psi_Hist_Both_BCAL->Fill(psi);
2030 
2031  double makes_sense = 0;
2032  if (R1 > R2 && dt1 > dt2) makes_sense = 1;
2033  if (R1 < R2 && dt1 < dt2) makes_sense = 1;
2034  if (R1 < R2 && dt1 > dt2) makes_sense = 0;
2035  if (R2 < R1 && dt1 > dt1) makes_sense = 0;
2036  /* for(unsigned int loc_j = 0; loc_j < associated_clusters1.size(); loc_j++)
2037  {
2038  const DBCALCluster* a_cluster1 = associated_clusters1[loc_j];
2039  vector<const DBCALPoint*> associated_points1;
2040  a_cluster1->Get(associated_points1);
2041  for(unsigned int loc_jj = 0; loc_jj < associated_clusters2.size(); loc_jj++)
2042  {
2043  const DBCALCluster* a_cluster2 = associated_clusters2[loc_jj];
2044  vector<const DBCALPoint*> associated_points2;
2045  a_cluster2->Get(associated_points2);
2046  for(unsigned int loc_k = 0; loc_k < associated_points1.size(); loc_k++)
2047  {
2048  int module1 = associated_points1[loc_k]->module();
2049  int layer1 = associated_points1[loc_k]->layer();
2050  int sector1 = associated_points1[loc_k]->sector();
2051  double energy1 = associated_points1[loc_k]->E();
2052  int channel_per_module1;
2053  if (layer1 == 1 && sector1 == 1) channel_per_module1 = 0;
2054  if (layer1 == 1 && sector1 == 2) channel_per_module1 = 1;
2055  if (layer1 == 1 && sector1 == 3) channel_per_module1 = 2;
2056  if (layer1 == 1 && sector1 == 4) channel_per_module1 = 3;
2057  if (layer1 == 2 && sector1 == 1) channel_per_module1 = 4;
2058  if (layer1 == 2 && sector1 == 2) channel_per_module1 = 5;
2059  if (layer1 == 2 && sector1 == 3) channel_per_module1 = 6;
2060  if (layer1 == 2 && sector1 == 4) channel_per_module1 = 7;
2061  if (layer1 == 3 && sector1 == 1) channel_per_module1 = 8;
2062  if (layer1 == 3 && sector1 == 2) channel_per_module1 = 9;
2063  if (layer1 == 3 && sector1 == 3) channel_per_module1 = 10;
2064  if (layer1 == 3 && sector1 == 4) channel_per_module1 = 11;
2065  if (layer1 == 4 && sector1 == 1) channel_per_module1 = 12;
2066  if (layer1 == 4 && sector1 == 2) channel_per_module1 = 13;
2067  if (layer1 == 4 && sector1 == 3) channel_per_module1 = 14;
2068  if (layer1 == 4 && sector1 == 4) channel_per_module1 = 15;
2069  channel = channel_per_module1 + (module1-1)*16;
2070 
2071  //cout << " moodule = " << module1 << " layer = " << layer1 << " sector = " << sector1 << " channel = " << channel << endl;
2072 
2073  //if( ptot.M() > 0.12 && ptot.M() < 0.18) cout << "shower1 energy = " << E1 << " shower2 energy = " << E2 << " shower1 time = " << dt1 << " shower 2 time = " << dt2 << " shower 1 z = " << dz1 << " shower 2 z = " << dz2 << " ptot = " << ptot.M() << " event num = " << eventnum << endl;
2074 
2075 
2076  for(unsigned int loc_kk = 0; loc_kk < associated_points1.size(); loc_kk++)
2077  {
2078  // int module2 = associated_points2[loc_kk]->module();
2079  // int layer2 = associated_points2[loc_kk]->layer();
2080  // int sector2 = associated_points2[loc_kk]->sector();
2081  // int energy2 = associated_points2[loc_kk]->E();
2082  // cout << "shower1 energy = " << E1 << " shower2 energy = " << E2 <<" point1 energy = " << energy1 << " point2 energy = " << energy2 << " channel = " << channel << " point size = " << associated_points1.size() << " shower size = " << locBCALShowers.size() << " ptot = " << ptot.M() << " event num = " << eventnum << endl;
2083 
2084  }
2085  }
2086  }
2087  } */
2088  //cout << " z1 = " << z1 << " z2 = " << z2 << " shower size = " << locBCALShowers.size() << " eventnum = " << eventnum << endl;
2089  BCAL_Neutrals->Fill();
2090  if (E1 > 0.5 && E2 > 0.5 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && kinfitVertexX > -15.0 && kinfitVertexX < 15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0) two_gamma_mass->Fill(ptot.M());
2091  if (E1 > 0.2 && E2 > 0.2 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0) bcal_diphoton_mass->Fill(ptot.M());
2092  if (E1 > 0.7 && E2 > 0.7 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && kinfitVertexX > -15.0 && kinfitVertexX < 15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0) bcal_diphoton_mass_highE->Fill(ptot.M());
2093  }
2094  }
2095 
2096 
2097  /* for(unsigned int j=0; j<locFCALShowers.size(); j++)
2098  {
2099  const DFCALShower *s3 = locFCALShowers[j];
2100 
2101  double dx3 = s3->getPosition().X() - kinfitVertexX;
2102  double dy3 = s3->getPosition().Y() - kinfitVertexY;
2103  double dz3 = 620.0 - kinfitVertexZ; // shift to coordinate relative to center of target
2104  double dt3 = s3->getTime();
2105  double R3 = sqrt(dx3*dx3 + dy3*dy3 + dz3*dz3);
2106  double E3 = s3->getEnergy();
2107  //double path = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));
2108  TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2109  int makes_sense2 = 0;
2110  if(R3 > R1 && dt3 > dt1) makes_sense2=1;
2111  if(R3 < R1 && dt3 < dt1) makes_sense2=1;
2112  if(R3 > R1 && dt3 < dt1) makes_sense2=0;
2113  if(R3 < R1 && dt3 > dt1) makes_sense2=0;
2114  double thetadeg3, thetarad3, phideg3, phirad3, cospsi3, psi3;
2115  thetarad3 = p3.Theta();
2116  phirad3 = p3.Phi();
2117  thetadeg3 = p3.Theta()*180.0/TMath::Pi();
2118  phideg3 = p3.Phi()*180.0/TMath::Pi();
2119  cospsi3 = sin(thetarad1)*sin(thetarad3)*cos(phirad1-phirad3)+cos(thetarad1)*cos(thetarad3);
2120  psi3=acos(cospsi3)*180/TMath::Pi();
2121 
2122  TLorentzVector ptot = p1+p3;
2123  //if(E1>0.3 && E2>0.3 && kinfitVertexZ > 50 && kinfitVertexZ < 80)
2124  if( kinfitVertexZ > 62.0 && kinfitVertexZ < 68.0 && E1 > 0.4 && E3 > 0.4 && makes_sense2==1 && kinfitVertexX < 15.0 && kinfitVertexX > -15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2125  //if(E1>0.2 && E>0.2) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2126 
2127  Theta_Hist_Split_Gammas->Fill(thetadeg1);
2128  Theta_Hist_Split_Gammas->Fill(thetadeg3);
2129  Phi_Hist_Split_Gammas->Fill(phideg1);
2130  Phi_Hist_Split_Gammas->Fill(phideg3);
2131  Psi_Hist_Split_Gammas->Fill(psi3);
2132 
2133  if(locBCALShowers.size()==1 && locFCALShowers.size()==1){
2134  //bcal_fcal_two_gamma_mass_cut->Fill(ptot.M());
2135  }
2136 
2137 
2138  } */
2139 
2140 
2141 
2142  for(unsigned int j=0; j<locFCALShowers.size(); j++)
2143  {
2144  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[j]))
2145  // continue;
2146 
2147  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[j]) != matchedFCALShowers.end()) continue;
2148  const DFCALShower *s3 = locFCALShowers[j];
2149 
2150  double dx3 = s3->getPosition().X() - kinfitVertexX;
2151  double dy3 = s3->getPosition().Y() - kinfitVertexY;
2152  double dz3 = s3->getPosition().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2153  double dt3 = s3->getTime();
2154  double R3 = sqrt(dx3*dx3 + dy3*dy3 + dz3*dz3);
2155  double E3 = s3->getEnergy();
2156  fcal_x = s3->getPosition().X();
2157  fcal_y = s3->getPosition().Y();
2158  fcal_z = s3->getPosition().Z();
2159  fcal_E = s3->getEnergy();
2160  fcal_t = s3->getTime();
2161  FCALShowers_per_event = locFCALShowers.size();
2162  //double path = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));
2163  TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2164  fcal_p=p3.M();
2165  double thetadeg3, thetarad3, phideg3, phirad3, cospsi3, psi3;
2166  thetarad3 = p3.Theta();
2167  phirad3 = p3.Phi();
2168  thetadeg3 = p3.Theta()*180.0/TMath::Pi();
2169  phideg3 = p3.Phi()*180.0/TMath::Pi();
2170 
2171 
2172 
2173  for( unsigned int i=0; i<locBCALShowers.size(); i++)
2174  {
2175  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
2176  // continue;
2177 
2178  if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end()) continue;
2179  const DBCALShower *s1 = locBCALShowers[i];
2180  double dx1 = s1->x - kinfitVertexX;
2181  double dy1 = s1->y - kinfitVertexY;
2182  double dz1 = s1->z - kinfitVertexZ;
2183  double dt1 = s1->t;
2184  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2185  bcal_x = s1->x;
2186  bcal_y = s1->y;
2187  bcal_z = s1->z;
2188  bcal_E = s1->E;
2189  bcal_t = s1->t;
2190  vertexZ = kinfitVertexZ;
2191  vertexX = kinfitVertexX;
2192  vertexY = kinfitVertexY;
2193  BCALShowers_per_event = locBCALShowers.size();
2194  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
2195  double E1 = s1->E;
2196  double ECalc = s1->E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
2197  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
2198  bcal_p = p1.M();
2199  double thetadeg1, thetarad1, phideg1, phirad1;
2200  thetadeg1 = p1.Theta()*180.0/TMath::Pi();
2201  thetarad1 = p1.Theta();
2202  phideg1 = p1.Phi()*180.0/TMath::Pi();
2203  phirad1 = p1.Phi();
2204  bcal_phi = phirad1;
2205  TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
2206  TLorentzVector ptot = p1+p3;
2207  inv_mass = ptot.M();
2208 
2209  cospsi3 = sin(thetarad1)*sin(thetarad3)*cos(phirad1-phirad3)+cos(thetarad1)*cos(thetarad3);
2210  psi3=acos(cospsi3)*180/TMath::Pi();
2211  //if(E1>0.3 && E2>0.3 && kinfitVertexZ > 50 && kinfitVertexZ < 80)
2212  if( kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && E1 > 0.5 && E3 > 0.5) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2213  //if(E1>0.2 && E>0.2) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2214 
2215  Theta_Hist_Split_Gammas->Fill(thetadeg1);
2216  Theta_Hist_Split_Gammas->Fill(thetadeg3);
2217  Phi_Hist_Split_Gammas->Fill(phideg1);
2218  Phi_Hist_Split_Gammas->Fill(phideg3);
2219  Psi_Hist_Split_Gammas->Fill(psi3);
2220 
2221 
2222  Split_Gamma_Neutrals->Fill();
2223  }
2224 
2225  }
2226 
2227 
2228 
2229  for(unsigned int j=0; j<locFCALClusters.size(); j++)
2230  {
2231  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[j]))
2232  // continue;
2233 
2234  if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[j]) != matchedFCALClusters.end()) continue;
2235  const DFCALCluster *s3 = locFCALClusters[j];
2236 
2237  double dx3 = s3->getCentroid().X() - kinfitVertexX;
2238  double dy3 = s3->getCentroid().Y() - kinfitVertexY;
2239  double dz3 = s3->getCentroid().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2240  double dt3 = s3->getTime();
2241  double R3 = sqrt(dx3*dx3 + dy3*dy3 + dz3*dz3);
2242  double E3 = s3->getEnergy();
2243  fcal_x = s3->getCentroid().X();
2244  fcal_y = s3->getCentroid().Y();
2245  fcal_z = s3->getCentroid().Z();
2246  fcal_E = s3->getEnergy();
2247  fcal_t = s3->getTime();
2248  FCALClusters_per_event = locFCALClusters.size();
2249  //double path = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));
2250  TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2251  fcal_p=p3.M();
2252  double thetadeg3, thetarad3, phideg3, phirad3, cospsi3, psi3;
2253  thetarad3 = p3.Theta();
2254  phirad3 = p3.Phi();
2255  thetadeg3 = p3.Theta()*180.0/TMath::Pi();
2256  phideg3 = p3.Phi()*180.0/TMath::Pi();
2257 
2258 
2259 
2260  for( unsigned int i=0; i<locBCALShowers.size(); i++)
2261  {
2262  // if(locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i]))
2263  // continue;
2264 
2265  if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end()) continue;
2266  const DBCALShower *s1 = locBCALShowers[i];
2267  double dx1 = s1->x - kinfitVertexX;
2268  double dy1 = s1->y - kinfitVertexY;
2269  double dz1 = s1->z - kinfitVertexZ;
2270  double dt1 = s1->t;
2271  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2272  bcal_x = s1->x;
2273  bcal_y = s1->y;
2274  bcal_z = s1->z;
2275  bcal_E = s1->E_raw;
2276  bcal_t = s1->t;
2277  vertexZ = kinfitVertexZ;
2278  vertexX = kinfitVertexX;
2279  vertexY = kinfitVertexY;
2280  BCALShowers_per_event = locBCALShowers.size();
2281  //double path1 = sqrt((dx1-kinfitVertexX)*(dx1-kinfitVertexX)+(dy1-kinfitVertexY)*(dy1-kinfitVertexY)+(dz1)*(dz1));
2282  double E1_raw = s1->E_raw;
2283  double ECalc = s1->E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
2284  TLorentzVector p1(E1_raw*dx1/R1, E1_raw*dy1/R1, E1_raw*dz1/R1, E1_raw);
2285  bcal_p = p1.M();
2286  double thetadeg1, thetarad1, phideg1, phirad1;
2287  thetadeg1 = p1.Theta()*180.0/TMath::Pi();
2288  thetarad1 = p1.Theta();
2289  phideg1 = p1.Phi()*180.0/TMath::Pi();
2290  phirad1 = p1.Phi();
2291  bcal_phi = phirad1;
2292  TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
2293  TLorentzVector ptot = p1+p3;
2294  inv_mass = ptot.M();
2295 
2296  cospsi3 = sin(thetarad1)*sin(thetarad3)*cos(phirad1-phirad3)+cos(thetarad1)*cos(thetarad3);
2297  psi3=acos(cospsi3)*180/TMath::Pi();
2298  //if(E1>0.3 && E2>0.3 && kinfitVertexZ > 50 && kinfitVertexZ < 80)
2299  //if( kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && E1 > 0.5 && E3 > 0.5) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2300  //if(E1>0.2 && E>0.2) bcal_fcal_two_gamma_mass->Fill(ptot.M());
2301 
2302  Theta_Hist_Split_Gammas->Fill(thetadeg1);
2303  Theta_Hist_Split_Gammas->Fill(thetadeg3);
2304  Phi_Hist_Split_Gammas->Fill(phideg1);
2305  Phi_Hist_Split_Gammas->Fill(phideg3);
2306  Psi_Hist_Split_Gammas->Fill(psi3);
2307 
2308 
2309  Split_Gamma_Neutrals_raw->Fill();
2310  }
2311 
2312  }
2313 
2314 
2315 
2316 
2317 
2318 
2319  for(unsigned int ii=0; ii<locFCALShowers.size(); ii++)
2320  {
2321 
2322  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[ii]))
2323  // continue;
2324  FCALShowers_per_event = locFCALShowers.size();
2325  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[ii]) != matchedFCALShowers.end()) continue;
2326  const DFCALShower *s2 = locFCALShowers[ii];
2327 
2328  //double dx, dy, dz, R, thetarad1, phirad1,
2329  double dx1 = s2->getPosition().X() - kinfitVertexX;
2330  double dy1 = s2->getPosition().Y() - kinfitVertexY;
2331  double dz1 = s2->getPosition().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2332  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2333  E1 = s2->getEnergy();
2334  double dt1 = s2->getTime();
2335  //double path1 = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));g
2336  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
2337  x1 = s2->getPosition().X();
2338  y1 = s2->getPosition().Y();
2339  z1 = s2->getPosition().Z();
2340  t1 = s2->getTime();
2341  p1_mag = p1.M();
2342 
2343  vertexZ = kinfitVertexZ;
2344  vertexX = kinfitVertexX;
2345  vertexY = kinfitVertexY;
2346  //double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
2347 
2348  double thetarad2 = p1.Theta();
2349  double phirad2 = p1.Phi();
2350  double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2351  double phideg2 = p1.Phi()*180.0/TMath::Pi();
2352  //double cospsi = sin(thetarad1)*sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
2353  //double psi=acos(cospsi)*180/TMath::Pi();
2354 
2355 
2356 
2357  for(unsigned int k=ii+1; k<locFCALShowers.size(); k++)
2358  {
2359  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[k]))
2360  // continue;
2361  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[k]) != matchedFCALShowers.end()) continue;
2362  const DFCALShower *s3 = locFCALShowers[k];
2363  double dx2 = s3->getPosition().X() - kinfitVertexX;
2364  double dy2 = s3->getPosition().Y() - kinfitVertexY;
2365  double dz2 = s3->getPosition().Z()-kinfitVertexZ; // shift to coordinate relative to center of target
2366  double R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2367  E2 = s3->getEnergy();
2368  double dt2 = s3->getTime();
2369  x2 = s3->getPosition().X();
2370  y2 = s3->getPosition().Y();
2371  z2 = s3->getPosition().Z();
2372  t2 = s3->getTime();
2373 
2374 
2375  TLorentzVector p2(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
2376  p2_mag = p2.M();
2377  double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2378  int makes_sense3 = 0;
2379  if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2380  if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2381  if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2382  if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2383  TLorentzVector ptot = p2+p1;
2384  inv_mass = ptot.M();
2385  thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2386  thetarad3 = p2.Theta();
2387  phideg3 = p2.Phi()*180.0/TMath::Pi();
2388  phirad3 = p2.Phi();
2389  cospsi3 = sin(thetarad2)*sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2390  psi3=acos(cospsi3)*180/TMath::Pi();
2391  Theta_Hist_Both_FCAL->Fill(thetadeg3);
2392  Theta_Hist_Both_FCAL->Fill(thetadeg2);
2393  Phi_Hist_Both_FCAL->Fill(phideg3);
2394  Phi_Hist_Both_FCAL->Fill(phideg2);
2395  Psi_Hist_Both_FCAL->Fill(psi3);
2396  //if(E2>0.3 && E3>0.3 && kinfitVertexZ > 50 && kinfitVertexZ < 80)
2397  //double deltaT = TMath::Abs(dt1-dt2);
2398  //if(ptot.M() > 0.05 && ptot.M() < 1.0)Time_Diff->Fill(deltaT);
2399  // cout << " E1 = " << E1 << " E2 = " << E2 << " d1 = " << dt1 << " t2 = " << dt2 << " ptot = " << ptot.M() << endl;
2400  FCAL_Neutrals->Fill();
2401  if( kinfitVertexZ > 62.0 && kinfitVertexZ < 68.0 && E2 > 0.5 && E1 > 0.5 && kinfitVertexX < 15.0 && kinfitVertexX > -15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0) two_fcal_gamma_mass->Fill(ptot.M());
2402  //cout << " shower1 energy = " << E1 << " shower2 energy = " << E2 << " ptot = " << ptot.M() << " fcal shower size = " << locFCALShowers.size() << endl;
2403  //if(E2>0.2 && E1>0.2) two_fcal_gamma_mass->Fill(ptot.M());
2404 
2405  }
2406  }
2407 
2408 
2409  for(unsigned int ii=0; ii<locFCALShowers.size(); ii++)
2410  {
2411 
2412  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[ii]))
2413  // continue;
2414  FCALShowers_per_event = locFCALShowers.size();
2415  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[ii]) != matchedFCALShowers.end()) continue;
2416  const DFCALShower *s2 = locFCALShowers[ii];
2417 
2418  //double dx, dy, dz, R, thetarad1, phirad1,
2419  double dx1 = s2->getPosition().X() - kinfitVertexX;
2420  double dy1 = s2->getPosition().Y() - kinfitVertexY;
2421  double dz1 = s2->getPosition().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2422  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2423  E1 = s2->getEnergy();
2424  double dt1 = s2->getTime();
2425  //double path1 = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));g
2426  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
2427  x1 = s2->getPosition().X();
2428  y1 = s2->getPosition().Y();
2429  z1 = s2->getPosition().Z();
2430  t1 = s2->getTime();
2431  p1_mag = p1.M();
2432 
2433  vertexZ = kinfitVertexZ;
2434  vertexX = kinfitVertexX;
2435  vertexY = kinfitVertexY;
2436  //double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
2437 
2438  double thetarad2 = p1.Theta();
2439  double phirad2 = p1.Phi();
2440  double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2441  double phideg2 = p1.Phi()*180.0/TMath::Pi();
2442  //double cospsi = sin(thetarad1)*sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
2443  //double psi=acos(cospsi)*180/TMath::Pi();
2444 
2445 
2446 
2447  for(unsigned int k=ii+1; k<locFCALShowers.size(); k++)
2448  {
2449  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[k]))
2450  // continue;
2451  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[k]) != matchedFCALShowers.end()) continue;
2452  const DFCALShower *s3 = locFCALShowers[k];
2453  double dx2 = s3->getPosition().X() - kinfitVertexX;
2454  double dy2 = s3->getPosition().Y() - kinfitVertexY;
2455  double dz2 = s3->getPosition().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2456  double R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2457  E2 = s3->getEnergy();
2458  double dt2 = s3->getTime();
2459  x2 = s3->getPosition().X();
2460  y2 = s3->getPosition().Y();
2461  z2 = s3->getPosition().Z();
2462  t2 = s3->getTime();
2463 
2464 
2465  TLorentzVector p2(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
2466  p2_mag = p2.M();
2467  double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2468  int makes_sense3 = 0;
2469  if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2470  if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2471  if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2472  if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2473  thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2474  thetarad3 = p2.Theta();
2475  phideg3 = p2.Phi()*180.0/TMath::Pi();
2476  phirad3 = p2.Phi();
2477  cospsi3 = sin(thetarad2)*sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2478  psi3=acos(cospsi3)*180/TMath::Pi();
2479  Theta_Hist_Both_FCAL->Fill(thetadeg3);
2480  Theta_Hist_Both_FCAL->Fill(thetadeg2);
2481  Phi_Hist_Both_FCAL->Fill(phideg3);
2482  Phi_Hist_Both_FCAL->Fill(phideg2);
2483  Psi_Hist_Both_FCAL->Fill(psi3);
2484  for(unsigned int m = ii+2; m < locFCALShowers.size(); m ++)
2485  {
2486  if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[m]) != matchedFCALShowers.end()) continue;
2487  double dx3 = locFCALShowers[m]->getPosition().X() - kinfitVertexX;
2488  double dy3 = locFCALShowers[m]->getPosition().Y() - kinfitVertexY;
2489  double dz3 = locFCALShowers[m]->getPosition().Z() - kinfitVertexZ;
2490  double R3 = sqrt(dx3*dx3+dy3*dy3+dz3*dz3);
2491  E3 = locFCALShowers[m]->getEnergy();
2492  double dt3 = locFCALShowers[m]->getTime();
2493  x3 = locFCALShowers[m]->getPosition().X();
2494  y3 = locFCALShowers[m]->getPosition().Y();
2495  z3 = locFCALShowers[m]->getPosition().Z();
2496  t3 = locFCALShowers[m]->getTime();
2497  TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2498  p3_mag = p3.M();
2499  TLorentzVector ptot = p1 + p2 + p3;
2500  inv_mass = ptot.M();
2501  Triple_FCAL_Neutrals->Fill();
2502  }
2503 
2504  }
2505  }
2506 
2507 
2508 
2509 
2510 
2511 
2512  // time to look at clusters lmao
2513 
2514  for(unsigned int ii=0; ii<locFCALClusters.size(); ii++)
2515  {
2516 
2517  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[ii]))
2518  // continue;
2519  FCALClusters_per_event = locFCALClusters.size();
2520  if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[ii]) != matchedFCALClusters.end()) continue;
2521  const DFCALCluster *s2 = locFCALClusters[ii];
2522 
2523  //double dx, dy, dz, R, thetarad1, phirad1,
2524  double dx1 = s2->getCentroid().X() - kinfitVertexX;
2525  double dy1 = s2->getCentroid().Y() - kinfitVertexY;
2526  double dz1 = s2->getCentroid().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2527  double R1 = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2528  E1 = s2->getEnergy();
2529  double dt1 = s2->getTime();
2530  //double path1 = sqrt((dx-kinfitVertexX)*(dx-kinfitVertexX)+(dy-kinfitVertexY)*(dy-kinfitVertexY)+(dz)*(dz));g
2531  TLorentzVector p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
2532  x1 = s2->getCentroid().X();
2533  y1 = s2->getCentroid().Y();
2534  z1 = s2->getCentroid().Z();
2535  t1 = s2->getTime();
2536  p1_mag = p1.M();
2537 
2538  vertexZ = kinfitVertexZ;
2539  vertexX = kinfitVertexX;
2540  vertexY = kinfitVertexY;
2541  //double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
2542 
2543  double thetarad2 = p1.Theta();
2544  double phirad2 = p1.Phi();
2545  double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2546  double phideg2 = p1.Phi()*180.0/TMath::Pi();
2547  //double cospsi = sin(thetarad1)*sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
2548  //double psi=acos(cospsi)*180/TMath::Pi();
2549 
2550 
2551 
2552  for(unsigned int k=ii+1; k<locFCALClusters.size(); k++)
2553  {
2554  // if(locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[k]))
2555  // continue;
2556  if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[k]) != matchedFCALClusters.end()) continue;
2557  const DFCALCluster *s3 = locFCALClusters[k];
2558  double dx2 = s3->getCentroid().X() - kinfitVertexX;
2559  double dy2 = s3->getCentroid().Y() - kinfitVertexY;
2560  double dz2 = s3->getCentroid().Z() - kinfitVertexZ; // shift to coordinate relative to center of target
2561  double R2 = sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2562  E2 = s3->getEnergy();
2563  double dt2 = s3->getTime();
2564  x2 = s3->getCentroid().X();
2565  y2 = s3->getCentroid().Y();
2566  z2 = s3->getCentroid().Z();
2567  t2 = s3->getTime();
2568 
2569 
2570  TLorentzVector p2(E2*dx2/R2, E2*dy2/R2, E2*dz2/R2, E2);
2571  p2_mag = p2.M();
2572  double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2573  int makes_sense3 = 0;
2574  if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2575  if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2576  if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2577  if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2578  thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2579  thetarad3 = p2.Theta();
2580  phideg3 = p2.Phi()*180.0/TMath::Pi();
2581  phirad3 = p2.Phi();
2582  cospsi3 = sin(thetarad2)*sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2583  psi3=acos(cospsi3)*180/TMath::Pi();
2584  Theta_Hist_Both_FCAL->Fill(thetadeg3);
2585  Theta_Hist_Both_FCAL->Fill(thetadeg2);
2586  Phi_Hist_Both_FCAL->Fill(phideg3);
2587  Phi_Hist_Both_FCAL->Fill(phideg2);
2588  Psi_Hist_Both_FCAL->Fill(psi3);
2589  TLorentzVector ptot = p1 + p2;
2590  inv_mass = ptot.M();
2591  FCALClusterNeutrals->Fill();
2592  }
2593  }
2594 
2595 
2596 
2597 
2598 
2599 
2600 
2601 
2602 
2603 
2604 
2605 
2606  //const DMCThrown* locMCThrown;
2607  //const DMCThrown* locMCThrown2
2608 
2609 
2610 
2611 
2612 /* for(unsigned int i = 0; i < locMCThrowns.size(); i++)
2613  {
2614  const DMCThrown *locMCThrown = locMCThrowns[i];
2615  Particle_t locPID;
2616  locPID = (Particle_t) locMCThrown->type;
2617  if (locPID != Gamma) continue;
2618  int numZ = 0;
2619  double dz1;
2620  //if (locPID == Gamma)
2621  // {
2622  dz1 = locMCThrown->position().Z();
2623  //Thrown_Vertex->Fill(dz1);
2624  // }
2625  for(unsigned int j=i+1; j < locMCThrowns.size(); j++)
2626  {
2627  Particle_t locPID2;
2628  const DMCThrown *locMCThrown2 = locMCThrowns[j];
2629  locPID2 = (Particle_t) locMCThrown2->type;
2630  if (locPID2 != Gamma) continue;
2631  //if(locPID2==Gamma)
2632  //{
2633  double dz2 = locMCThrown2->position().Z();
2634  if(dz1 < (dz2 + 1) && dz1 > (dz2 - 1)) numZ = 1;
2635  else numZ = 2;
2636  //if(dz1 != dz2) numZ = 2;
2637  Thrown_Vertex_Frequency->Fill(numZ);
2638  // cout << " thrown size = " << locMCThrowns.size() << " vertex of gamma 1 = " << dz1 << " vertex of gamma 2 = " << dz2 << " numz = " << numZ << endl;
2639  //}
2640  }
2641 
2642 
2643 
2644  }
2645 
2646 
2647 
2648  Particle_t locPID2;
2649 //=(Particle_t) locMCThrown->type;
2650  for(unsigned int i = 0; i < locMCThrowns.size(); i++)
2651  {
2652  const DMCThrown *locMCThrown = locMCThrowns[i];
2653  Particle_t locPID;
2654  locPID = (Particle_t) locMCThrown->type;
2655  DVector3 locMomentum1;
2656  DLorentzVector lorentzmomentum1;
2657  TLorentzVector p1;
2658  if (locPID == Gamma)
2659  {
2660  DVector3 locMomentum1 = locMCThrown->momentum();
2661  DLorentzVector lorentzmomentum1 = locMCThrown->lorentzMomentum();
2662  double dx = locMCThrown->position().X();
2663  double dy = locMCThrown->position().Y();
2664  double dz = locMCThrown->position().Z();
2665  double R = sqrt(dx*dx + dy*dy + dz*dz);
2666  double E1 = locMCThrown->energy();
2667  TLorentzVector p1(E1*dx/R, E1*dy/R, E1*dz/R,E1);
2668  double locTheta = locMomentum1.Theta()*180.0/TMath::Pi();
2669  //Thrown_Vertex->Fill(dz);
2670  //cout << " vertex in z = " << dz << " locMCThrownSize = " << locMCThrowns.size() << endl;
2671 
2672 
2673  //}
2674  for(unsigned int j =i+1; j < locMCThrowns.size(); j++)
2675  {
2676  const DMCThrown *locMCThrown2 = locMCThrowns[j];
2677  locPID2 = (Particle_t) locMCThrown2->type;
2678  //if (locPID2 == 1)
2679  //{
2680  double dx = locMCThrown2->position().X();
2681  double dy = locMCThrown2->position().Y();
2682  double dz = locMCThrown2->position().Z();
2683  double R = sqrt(dx*dx + dy*dy + dz*dz);
2684  double E2 = locMCThrown2->energy();
2685  TLorentzVector p2(E2*dx/R, E2*dy/R, E2*dz/R,E2);
2686  DVector3 locMomentum2 = locMCThrown2->momentum();
2687  double momentum2Mag = locMomentum2.Mag();
2688  DLorentzVector lorentzmomentum2 = locMCThrown2->lorentzMomentum();
2689  DVector3 locPSum = locMomentum1 + locMomentum2;
2690  DLorentzVector LorentzPSum = lorentzmomentum1 + lorentzmomentum2;
2691  TLorentzVector ptot = p1+p2;
2692 
2693  Thrown_Gamma_Theta->Fill(lorentzmomentum1.Theta()*180.0/TMath::Pi());
2694  Thrown_Gamma_Theta->Fill(lorentzmomentum2.Theta()*180.0/TMath::Pi());
2695  Thrown_Inv_Mass->Fill(LorentzPSum.M());
2696 
2697  }
2698  }
2699  }
2700 */
2701  // Total energy in hits
2702  double Ehit_tot = 0.0;
2703  for(unsigned int i=0; i<bcalhits.size(); i++){
2704  Ehit_tot += bcalhits[i]->E;
2705  }
2706  Etot_hits->Fill(Ehit_tot);
2707 
2708  // Truth values
2709  double Etruth_tot = 0.0;
2710  double z_truth = 0.0;
2711  for(unsigned int i=0; i<truthshowers.size(); i++){
2712  Etruth_tot += truthshowers[i]->E;
2713  z_truth += truthshowers[i]->E*truthshowers[i]->z;
2714  }
2715  z_truth/=Etruth_tot;
2716  //Etot_truth->Fill(Etruth_tot);
2717 
2718  // Compare to thrown values
2719  double Etot_thrown=0.0;
2720  for(unsigned int i=0; i<mcthrowns.size(); i++){
2721  Etot_thrown += mcthrowns[i]->energy();
2722  for(unsigned int j=0; j<locBCALShowers.size(); j++){
2723  double z = locBCALShowers[j]->z;
2724  //Erec_over_Ethrown_vs_z->Fill(z, locBCALShowers[j]->E/mcthrowns[i]->energy());
2725 
2726  double E = locBCALShowers[j]->E*(1.106+(z-208.4)*(z-208.4)*6.851E-6);
2727  //E_over_Erec_vs_z->Fill(z, mcthrowns[i]->energy()/E);
2728  }
2729  }
2730 
2731  //Ereconstructed_vs_Ethrown->Fill(Etot_thrown, Etot_reconstructed);
2732  //Etruth_over_Ethrown_vs_z->Fill(z_truth, Etruth_tot/Etot_thrown);
2733 
2734  // Single thrown particle
2735  if(mcthrowns.size()==1){
2736  const DMCThrown* mcthrown = mcthrowns[0];
2737  if(mcthrown->momentum().Theta()>0.0001){
2738  double z = mcthrown->position().Z() + 65.0/tan(mcthrown->momentum().Theta());
2739  double Ethrown = 1.0; // for some reason, mcthrown->E is zero right now.
2740  // I fudge this for now since I know all of the events thrw 1.0GeV
2741  //Edeposited_over_Ethrown_vs_z->Fill(z, Ehit_tot/Ethrown);
2742  }
2743  }
2744 
2745  japp->RootUnLock(); //RELEASE ROOT LOCK
2746 
2747  /*
2748  //Optional: Save event to output REST file. Use this to create skims.
2749  dEventWriterREST->Write_RESTEvent(locEventLoop, "BCAL_Shower"); //string is part of output file name
2750  */
2751 
2752  return NOERROR;
2753 }
2754 
2755 //------------------
2756 // erun
2757 //------------------
2758 jerror_t DEventProcessor_BCAL_Shower::erun(void)
2759 {
2760  // This is called whenever the run number changes, before it is
2761  // changed to give you a chance to clean up before processing
2762  // events from the next run number.
2763  return NOERROR;
2764 }
2765 
2766 //------------------
2767 // fini
2768 //------------------
2769 jerror_t DEventProcessor_BCAL_Shower::fini(void)
2770 {
2771  // Called before program exit after event processing is finished.
2772  return NOERROR;
2773 }
2774 
DEventProcessor_BCAL_Shower::Layer1_Energy_vs_Channel
TH2F * Layer1_Energy_vs_Channel
Definition: DEventProcessor_BCAL_Shower.h:159
DEventProcessor_BCAL_Shower::layer3_energysum
Float_t layer3_energysum
Definition: DEventProcessor_BCAL_Shower.h:49
DEventProcessor_BCAL_Shower::y1
Float_t y1
Definition: DEventProcessor_BCAL_Shower.h:79
DFCALCluster::getCentroid
DVector3 getCentroid() const
Definition: DFCALCluster.h:183
DEventProcessor_BCAL_Shower::dAnalysisUtilities
const DAnalysisUtilities * dAnalysisUtilities
Definition: DEventProcessor_BCAL_Shower.h:116
DFCALShower::getEnergy
double getEnergy() const
Definition: DFCALShower.h:156
DEventProcessor_BCAL_Shower::mass_summed
TH1F * mass_summed
Definition: DEventProcessor_BCAL_Shower.h:136
DEventProcessor_BCAL_Shower::init
jerror_t init(void)
Called once at program start.
Definition: DEventProcessor_BCAL_Shower.cc:44
DEventProcessor_BCAL_Shower::FCALClusters_per_event
uint32_t FCALClusters_per_event
Definition: DEventProcessor_BCAL_Shower.h:89
DEventProcessor_BCAL_Shower::inv_mass_raw
Float_t inv_mass_raw
Definition: DEventProcessor_BCAL_Shower.h:96
DBCALShower.h
DBCALPoint::module
int module() const
Definition: DBCALPoint.h:64
DEventProcessor_BCAL_Shower::vertexZ
Float_t vertexZ
Definition: DEventProcessor_BCAL_Shower.h:84
DEventProcessor_BCAL_Shower::Split_Gamma_Neutrals_raw
TTree * Split_Gamma_Neutrals_raw
Definition: DEventProcessor_BCAL_Shower.h:41
DEventProcessor_BCAL_Shower::E1_Vs_E2CosTheta
TH2F * E1_Vs_E2CosTheta
Definition: DEventProcessor_BCAL_Shower.h:153
DEventProcessor_BCAL_Shower::layer1_energysum
Float_t layer1_energysum
Definition: DEventProcessor_BCAL_Shower.h:47
DEventProcessor_BCAL_Shower::Eoverp_plus_cuts
TH1F * Eoverp_plus_cuts
Definition: DEventProcessor_BCAL_Shower.h:178
DEventProcessor_BCAL_Shower::Eoverp_minus_nocuts
TH1F * Eoverp_minus_nocuts
Definition: DEventProcessor_BCAL_Shower.h:181
DEventProcessor_BCAL_Shower::vertexY
Float_t vertexY
Definition: DEventProcessor_BCAL_Shower.h:86
DEventProcessor_BCAL_Shower::dEventWriterREST
const DEventWriterREST * dEventWriterREST
Definition: DEventProcessor_BCAL_Shower.h:207
DEventProcessor_BCAL_Shower::VertexZ
TH1F * VertexZ
Definition: DEventProcessor_BCAL_Shower.h:147
DEventProcessor_BCAL_Shower::Phi_Hist_Both_BCAL
TH1F * Phi_Hist_Both_BCAL
Definition: DEventProcessor_BCAL_Shower.h:138
DEventProcessor_BCAL_Shower::brun
jerror_t brun(jana::JEventLoop *locEventLoop, int locRunNumber)
Called every time a new run number is detected.
Definition: DEventProcessor_BCAL_Shower.cc:347
DEventProcessor_BCAL_Shower::Theta_Hist_Both_FCAL
TH1F * Theta_Hist_Both_FCAL
Definition: DEventProcessor_BCAL_Shower.h:143
DEventProcessor_BCAL_Shower::goodVertexZ
TH1F * goodVertexZ
Definition: DEventProcessor_BCAL_Shower.h:148
DEventProcessor_BCAL_Shower::BCALPoint_Neutral
TTree * BCALPoint_Neutral
Definition: DEventProcessor_BCAL_Shower.h:35
DEventProcessor_BCAL_Shower::FCALClusterNeutrals
TTree * FCALClusterNeutrals
Definition: DEventProcessor_BCAL_Shower.h:38
DVector3
TVector3 DVector3
Definition: DVector3.h:14
DEventProcessor_BCAL_Shower::fcal_y
Float_t fcal_y
Definition: DEventProcessor_BCAL_Shower.h:106
DEventProcessor_BCAL_Shower::FCALShowers_per_event
uint32_t FCALShowers_per_event
Definition: DEventProcessor_BCAL_Shower.h:88
DEventProcessor_BCAL_Shower::dTargetZCenter
double dTargetZCenter
Definition: DEventProcessor_BCAL_Shower.h:202
DEventProcessor_BCAL_Shower::Layer3_Energy_pos
TH1F * Layer3_Energy_pos
Definition: DEventProcessor_BCAL_Shower.h:173
x
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
DEventProcessor_BCAL_Shower::Layer3_Energy
TH1F * Layer3_Energy
Definition: DEventProcessor_BCAL_Shower.h:165
DEventProcessor_BCAL_Shower::phi2
Float_t phi2
Definition: DEventProcessor_BCAL_Shower.h:83
DEventProcessor_BCAL_Shower::Evsp_minus
TH2F * Evsp_minus
Definition: DEventProcessor_BCAL_Shower.h:183
DEventProcessor_BCAL_Shower::p2_mag
Float_t p2_mag
Definition: DEventProcessor_BCAL_Shower.h:92
DFCALShower::getTime
double getTime() const
Definition: DFCALShower.h:161
DEventProcessor_BCAL_Shower::L1_pathlength
Float_t L1_pathlength
Definition: DEventProcessor_BCAL_Shower.h:61
DEventProcessor_BCAL_Shower::z2
Float_t z2
Definition: DEventProcessor_BCAL_Shower.h:74
p1
double p1
Definition: JEventProcessor_PS_E_calib.cc:39
DEventProcessor_BCAL_Shower::x1
Float_t x1
Definition: DEventProcessor_BCAL_Shower.h:76
y
#define y
DEventProcessor_BCAL_Shower::fcal_z
Float_t fcal_z
Definition: DEventProcessor_BCAL_Shower.h:107
DEventProcessor_BCAL_Shower::Psi_Hist_Both_BCAL
TH1F * Psi_Hist_Both_BCAL
Definition: DEventProcessor_BCAL_Shower.h:139
DEventProcessor_BCAL_Shower::t3
Float_t t3
Definition: DEventProcessor_BCAL_Shower.h:72
DBCALPoint::phi
float phi() const
Definition: DBCALPoint.h:56
DBCALPoint::sector
int sector() const
Definition: DBCALPoint.h:66
DEventProcessor_BCAL_Shower::bcal_y
Float_t bcal_y
Definition: DEventProcessor_BCAL_Shower.h:99
DEventProcessor_BCAL_Shower::Eoverp_minus_cuts
TH1F * Eoverp_minus_cuts
Definition: DEventProcessor_BCAL_Shower.h:179
DKinematicData::position
const DVector3 & position(void) const
Definition: DKinematicData.h:47
DEventProcessor_BCAL_Shower::sector
uint32_t sector
Definition: DEventProcessor_BCAL_Shower.h:54
DEventProcessor_BCAL_Shower::fcal_x
Float_t fcal_x
Definition: DEventProcessor_BCAL_Shower.h:105
DEventProcessor_BCAL_Shower::two_fcal_gamma_mass
TH1F * two_fcal_gamma_mass
Definition: DEventProcessor_BCAL_Shower.h:135
DEventProcessor_BCAL_Shower::Evsp_plus
TH2F * Evsp_plus
Definition: DEventProcessor_BCAL_Shower.h:182
DEventProcessor_BCAL_Shower::layer2_energysum
Float_t layer2_energysum
Definition: DEventProcessor_BCAL_Shower.h:48
DEventProcessor_BCAL_Shower::bcal_p
Float_t bcal_p
Definition: DEventProcessor_BCAL_Shower.h:103
DApplication.h
DEventProcessor_BCAL_Shower::bcal_E
Float_t bcal_E
Definition: DEventProcessor_BCAL_Shower.h:97
DEventProcessor_BCAL_Shower::BCALShowers_per_event
uint32_t BCALShowers_per_event
Definition: DEventProcessor_BCAL_Shower.h:87
DEventProcessor_BCAL_Shower::channel
uint32_t channel
Definition: DEventProcessor_BCAL_Shower.h:57
DEventProcessor_BCAL_Shower::erun
jerror_t erun(void)
Called every time run number changes, provided brun has been called.
Definition: DEventProcessor_BCAL_Shower.cc:2758
DFCALCluster::getTime
double getTime() const
Definition: DFCALCluster.h:165
DEventProcessor_BCAL_Shower::Phi_Hist_Both_FCAL
TH1F * Phi_Hist_Both_FCAL
Definition: DEventProcessor_BCAL_Shower.h:144
DEventProcessor_BCAL_Shower
Definition: DEventProcessor_BCAL_Shower.h:26
DEventProcessor_BCAL_Shower::BCALPoint_Charged_pos
TTree * BCALPoint_Charged_pos
Definition: DEventProcessor_BCAL_Shower.h:34
DEventProcessor_BCAL_Shower::vertexX
Float_t vertexX
Definition: DEventProcessor_BCAL_Shower.h:85
DEventProcessor_BCAL_Shower::Point_E_M12S2L2
TH1F * Point_E_M12S2L2
Definition: DEventProcessor_BCAL_Shower.h:191
DEventProcessor_BCAL_Shower::Eoverpvsp_plus
TH2F * Eoverpvsp_plus
Definition: DEventProcessor_BCAL_Shower.h:184
DEventProcessor_BCAL_Shower::E_shower
TH1F * E_shower
Definition: DEventProcessor_BCAL_Shower.h:132
DEventProcessor_BCAL_Shower::phi
Float_t phi
Definition: DEventProcessor_BCAL_Shower.h:58
DEventProcessor_BCAL_Shower::track_momentum
Float_t track_momentum
Definition: DEventProcessor_BCAL_Shower.h:46
DEventProcessor_BCAL_Shower::bcal_diphoton_mass
TH1F * bcal_diphoton_mass
Definition: DEventProcessor_BCAL_Shower.h:124
DEventProcessor_BCAL_Shower::p3_mag
Float_t p3_mag
Definition: DEventProcessor_BCAL_Shower.h:94
DAnalysisUtilities::Calc_DOCAVertex
double Calc_DOCAVertex(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2, DVector3 &locDOCAPoint) const
Definition: DAnalysisUtilities.cc:1071
DEventProcessor_BCAL_Shower::t1
Float_t t1
Definition: DEventProcessor_BCAL_Shower.h:70
DBCALPoint::layer
int layer() const
Definition: DBCALPoint.h:65
DBCALTruthShower.h
c1
Double_t c1[2][NMODULES]
Definition: tw_corr.C:68
DEventProcessor_BCAL_Shower::Eoverpvsp_minus
TH2F * Eoverpvsp_minus
Definition: DEventProcessor_BCAL_Shower.h:185
japp
JApplication * japp
DEventProcessor_BCAL_Shower::FCAL_Evsp
TH2F * FCAL_Evsp
Definition: DEventProcessor_BCAL_Shower.h:186
DEventProcessor_BCAL_Shower::Cluster_Energy
TH1F * Cluster_Energy
Definition: DEventProcessor_BCAL_Shower.h:155
DEventProcessor_BCAL_Shower::Psi_Hist_Split_Gammas
TH1F * Psi_Hist_Split_Gammas
Definition: DEventProcessor_BCAL_Shower.h:142
DEventProcessor_BCAL_Shower::Theta_Hist_Both_BCAL
TH1F * Theta_Hist_Both_BCAL
Definition: DEventProcessor_BCAL_Shower.h:137
DEventProcessor_BCAL_Shower::Layer4_Energy_vs_Channel
TH2F * Layer4_Energy_vs_Channel
Definition: DEventProcessor_BCAL_Shower.h:162
fcal_t
Definition: fcal_t.h:19
DEventProcessor_BCAL_Shower::BCALPoint_Charged_neg
TTree * BCALPoint_Charged_neg
Definition: DEventProcessor_BCAL_Shower.h:33
DEventProcessor_BCAL_Shower::Layer2_Energy_vs_Channel
TH2F * Layer2_Energy_vs_Channel
Definition: DEventProcessor_BCAL_Shower.h:160
DEventProcessor_BCAL_Shower::bcal_phi
Float_t bcal_phi
Definition: DEventProcessor_BCAL_Shower.h:101
DEventProcessor_BCAL_Shower::L4_pathlength
Float_t L4_pathlength
Definition: DEventProcessor_BCAL_Shower.h:64
DEventProcessor_BCAL_Shower::z1
Float_t z1
Definition: DEventProcessor_BCAL_Shower.h:73
DEventProcessor_BCAL_Shower::layer
uint32_t layer
Definition: DEventProcessor_BCAL_Shower.h:53
DBCALPoint::theta
float theta() const
Definition: DBCALPoint.h:53
DBCALPoint.h
DEventProcessor_BCAL_Shower::r
Float_t r
Definition: DEventProcessor_BCAL_Shower.h:59
DEventProcessor_BCAL_Shower::Eoverp_plus_nocuts
TH1F * Eoverp_plus_nocuts
Definition: DEventProcessor_BCAL_Shower.h:180
DEventProcessor_BCAL_Shower::two_gamma_mass
TH1F * two_gamma_mass
Definition: DEventProcessor_BCAL_Shower.h:123
DFCALCluster::getEnergy
double getEnergy() const
Definition: DFCALCluster.h:155
DBCALPoint::E
float E() const
Definition: DBCALPoint.h:38
DEventProcessor_BCAL_Shower::Layer2_Energy
TH1F * Layer2_Energy
Definition: DEventProcessor_BCAL_Shower.h:164
DEventProcessor_BCAL_Shower::fcal_p
Float_t fcal_p
Definition: DEventProcessor_BCAL_Shower.h:109
DEventProcessor_BCAL_Shower::FCAL_Eoverp_nocuts
TH1F * FCAL_Eoverp_nocuts
Definition: DEventProcessor_BCAL_Shower.h:189
InitPlugin
InitPlugin_t InitPlugin
DEventProcessor_BCAL_Shower::fini
jerror_t fini(void)
Called after last event of last event source has been processed.
Definition: DEventProcessor_BCAL_Shower.cc:2769
DEventProcessor_BCAL_Shower::Point_Module
TH1F * Point_Module
Definition: DEventProcessor_BCAL_Shower.h:157
DEventProcessor_BCAL_Shower::FCAL_Eoverp_cuts
TH1F * FCAL_Eoverp_cuts
Definition: DEventProcessor_BCAL_Shower.h:188
DEventProcessor_BCAL_Shower::x2
Float_t x2
Definition: DEventProcessor_BCAL_Shower.h:77
DEventProcessor_BCAL_Shower::Layer1_Energy_v2
TH1F * Layer1_Energy_v2
Definition: DEventProcessor_BCAL_Shower.h:167
bcal_t
Definition: bcal_t.h:19
DEventProcessor_BCAL_Shower::bcal_z
Float_t bcal_z
Definition: DEventProcessor_BCAL_Shower.h:100
DEventProcessor_BCAL_Shower::Thrown_Vertex
TH1F * Thrown_Vertex
Definition: DEventProcessor_BCAL_Shower.h:151
DBCALHit.h
DEventProcessor_BCAL_Shower::Point_E_M37S4L4
TH1F * Point_E_M37S4L4
Definition: DEventProcessor_BCAL_Shower.h:193
DEventProcessor_BCAL_Shower::E2_raw
Float_t E2_raw
Definition: DEventProcessor_BCAL_Shower.h:68
DEventProcessor_BCAL_Shower::E1_raw
Float_t E1_raw
Definition: DEventProcessor_BCAL_Shower.h:66
DEventProcessor_BCAL_Shower::Etot_hits
TH1F * Etot_hits
Definition: DEventProcessor_BCAL_Shower.h:198
DEventProcessor_BCAL_Shower::Phi_Hist_Split_Gammas
TH1F * Phi_Hist_Split_Gammas
Definition: DEventProcessor_BCAL_Shower.h:141
DEventProcessor_BCAL_Shower::y2
Float_t y2
Definition: DEventProcessor_BCAL_Shower.h:80
DEventProcessor_BCAL_Shower::L3_pathlength
Float_t L3_pathlength
Definition: DEventProcessor_BCAL_Shower.h:63
DEventProcessor_BCAL_Shower::Layer1_Energy_pos
TH1F * Layer1_Energy_pos
Definition: DEventProcessor_BCAL_Shower.h:171
DBCALPoint::z
float z() const
Definition: DBCALPoint.h:60
DEventProcessor_BCAL_Shower::All_Layer_Energy
TH1F * All_Layer_Energy
Definition: DEventProcessor_BCAL_Shower.h:176
DEventProcessor_BCAL_Shower::module
uint32_t module
Definition: DEventProcessor_BCAL_Shower.h:52
DEventProcessor_BCAL_Shower::E1
Float_t E1
Definition: DEventProcessor_BCAL_Shower.h:65
DEventProcessor_BCAL_Shower::Layer1_Energy
TH1F * Layer1_Energy
Definition: DEventProcessor_BCAL_Shower.h:163
DEventProcessor_BCAL_Shower::Layer4_Energy
TH1F * Layer4_Energy
Definition: DEventProcessor_BCAL_Shower.h:166
DEventProcessor_BCAL_Shower::Thrown_Inv_Mass
TH1F * Thrown_Inv_Mass
Definition: DEventProcessor_BCAL_Shower.h:150
DEventProcessor_BCAL_Shower::Point_E_M1S1L1
TH1F * Point_E_M1S1L1
Definition: DEventProcessor_BCAL_Shower.h:190
DEventProcessor_BCAL_Shower::Time_Diff
TH1F * Time_Diff
Definition: DEventProcessor_BCAL_Shower.h:177
DEventProcessor_BCAL_Shower::p1_mag
Float_t p1_mag
Definition: DEventProcessor_BCAL_Shower.h:90
sqrt
double sqrt(double)
DEventProcessor_BCAL_Shower::Point_E_M25S3L3
TH1F * Point_E_M25S3L3
Definition: DEventProcessor_BCAL_Shower.h:192
DEventProcessor_BCAL_Shower::Thrown_Vertex_Frequency
TH1F * Thrown_Vertex_Frequency
Definition: DEventProcessor_BCAL_Shower.h:152
DEventProcessor_BCAL_Shower::layer4_energysum
Float_t layer4_energysum
Definition: DEventProcessor_BCAL_Shower.h:50
sin
double sin(double)
DEventProcessor_BCAL_Shower::z3
Float_t z3
Definition: DEventProcessor_BCAL_Shower.h:75
DEventProcessor_BCAL_Shower::FCAL_Neutrals
TTree * FCAL_Neutrals
Definition: DEventProcessor_BCAL_Shower.h:37
DEventProcessor_BCAL_Shower::bcal_x
Float_t bcal_x
Definition: DEventProcessor_BCAL_Shower.h:98
DEventProcessor_BCAL_Shower::z
Float_t z
Definition: DEventProcessor_BCAL_Shower.h:60
DFCALShower.h
DEventProcessor_BCAL_Shower::energy_raw_shower
Float_t energy_raw_shower
Definition: DEventProcessor_BCAL_Shower.h:44
DEventProcessor_BCAL_Shower::E2
Float_t E2
Definition: DEventProcessor_BCAL_Shower.h:67
DEventProcessor_BCAL_Shower::Neutral_E
TH1F * Neutral_E
Definition: DEventProcessor_BCAL_Shower.h:133
DKinematicData::momentum
const DVector3 & momentum(void) const
Definition: DKinematicData.h:46
logical
long int logical
Definition: grkuta.cc:46
DEventProcessor_BCAL_Shower::Split_Gamma_Neutrals
TTree * Split_Gamma_Neutrals
Definition: DEventProcessor_BCAL_Shower.h:40
DEventProcessor_BCAL_Shower::BCAL_Neutrals
TTree * BCAL_Neutrals
Definition: DEventProcessor_BCAL_Shower.h:36
DEventProcessor_BCAL_Shower::Layer4_Energy_v2
TH1F * Layer4_Energy_v2
Definition: DEventProcessor_BCAL_Shower.h:170
DAnalysisUtilities.h
DEventProcessor_BCAL_Shower::Point_z
TH1F * Point_z
Definition: DEventProcessor_BCAL_Shower.h:158
DEventProcessor_BCAL_Shower::Layer3_Energy_v2
TH1F * Layer3_Energy_v2
Definition: DEventProcessor_BCAL_Shower.h:169
DEventProcessor_BCAL_Shower::energy_shower
Float_t energy_shower
Definition: DEventProcessor_BCAL_Shower.h:43
DEventProcessor_BCAL_Shower::Psi_Hist_Both_FCAL
TH1F * Psi_Hist_Both_FCAL
Definition: DEventProcessor_BCAL_Shower.h:145
DEventProcessor_BCAL_Shower::Layer3_Energy_vs_Channel
TH2F * Layer3_Energy_vs_Channel
Definition: DEventProcessor_BCAL_Shower.h:161
DEventProcessor_BCAL_Shower::p2_raw_mag
Float_t p2_raw_mag
Definition: DEventProcessor_BCAL_Shower.h:93
DEventProcessor_BCAL_Shower::bcal_fcal_two_gamma_mass
TH1F * bcal_fcal_two_gamma_mass
Definition: DEventProcessor_BCAL_Shower.h:128
DEventProcessor_BCAL_Shower::Pi0_velocity
TH1F * Pi0_velocity
Definition: DEventProcessor_BCAL_Shower.h:154
DEventProcessor_BCAL_Shower::bcal_diphoton_mass_highE
TH1F * bcal_diphoton_mass_highE
Definition: DEventProcessor_BCAL_Shower.h:125
DEventProcessor_BCAL_Shower::y3
Float_t y3
Definition: DEventProcessor_BCAL_Shower.h:81
DEventProcessor_BCAL_Shower::fcal_E
Float_t fcal_E
Definition: DEventProcessor_BCAL_Shower.h:104
DBCALPoint::r
float r() const
Definition: DBCALPoint.h:62
DEventProcessor_BCAL_Shower::t2
Float_t t2
Definition: DEventProcessor_BCAL_Shower.h:71
DEventProcessor_BCAL_Shower::p1_raw_mag
Float_t p1_raw_mag
Definition: DEventProcessor_BCAL_Shower.h:91
DEventProcessor_BCAL_Shower::Layer2_Energy_v2
TH1F * Layer2_Energy_v2
Definition: DEventProcessor_BCAL_Shower.h:168
DEventProcessor_BCAL_Shower::BCALShowerTrack_Energy
TH1F * BCALShowerTrack_Energy
Definition: DEventProcessor_BCAL_Shower.h:175
DEventProcessor_BCAL_Shower::theta
Float_t theta
Definition: DEventProcessor_BCAL_Shower.h:56
DEventProcessor_BCAL_Shower::energy_point
Float_t energy_point
Definition: DEventProcessor_BCAL_Shower.h:45
DEventProcessor_BCAL_Shower.h
DEventProcessor_BCAL_Shower::Point_Energy
TH1F * Point_Energy
Definition: DEventProcessor_BCAL_Shower.h:156
dir
TDirectory * dir
Definition: bcal_hist_eff.C:31
DEventProcessor_BCAL_Shower::Assoc_E
TH1F * Assoc_E
Definition: DEventProcessor_BCAL_Shower.h:134
DEventProcessor_BCAL_Shower::Layer2_Energy_pos
TH1F * Layer2_Energy_pos
Definition: DEventProcessor_BCAL_Shower.h:172
DEventProcessor_BCAL_Shower::z_shower
TH1F * z_shower
Definition: DEventProcessor_BCAL_Shower.h:131
DEventProcessor_BCAL_Shower::inv_mass
Float_t inv_mass
Definition: DEventProcessor_BCAL_Shower.h:95
DEventProcessor_BCAL_Shower::FCAL_Eoverpvsp
TH2F * FCAL_Eoverpvsp
Definition: DEventProcessor_BCAL_Shower.h:187
DEventProcessor_BCAL_Shower::x3
Float_t x3
Definition: DEventProcessor_BCAL_Shower.h:78
p2
double p2
Definition: JEventProcessor_PS_E_calib.cc:40
DFCALShower::getPosition
DVector3 getPosition() const
Definition: DFCALShower.h:151
DMCThrown.h
DEventProcessor_BCAL_Shower::evnt
jerror_t evnt(jana::JEventLoop *locEventLoop, uint64_t locEventNumber)
Called every event.
Definition: DEventProcessor_BCAL_Shower.cc:376
DEventProcessor_BCAL_Shower::two_fcal_gamma_mass_test
TH1F * two_fcal_gamma_mass_test
Definition: DEventProcessor_BCAL_Shower.h:146
DEventProcessor_BCAL_Shower::Layer4_Energy_pos
TH1F * Layer4_Energy_pos
Definition: DEventProcessor_BCAL_Shower.h:174
DEventProcessor_BCAL_Shower::Triple_FCAL_Neutrals
TTree * Triple_FCAL_Neutrals
Definition: DEventProcessor_BCAL_Shower.h:39
DEventProcessor_BCAL_Shower::eventnum
uint32_t eventnum
Definition: DEventProcessor_BCAL_Shower.h:55
DEventProcessor_BCAL_Shower::Thrown_Gamma_Theta
TH1F * Thrown_Gamma_Theta
Definition: DEventProcessor_BCAL_Shower.h:149
DEventProcessor_BCAL_Shower::dEventWriterROOT
const DEventWriterROOT * dEventWriterROOT
Definition: DEventProcessor_BCAL_Shower.h:206
DEventProcessor_BCAL_Shower::phi1
Float_t phi1
Definition: DEventProcessor_BCAL_Shower.h:82
DBCALShower::E_raw
float E_raw
Definition: DBCALShower.h:18
DEventProcessor_BCAL_Shower::Theta_Hist_Split_Gammas
TH1F * Theta_Hist_Split_Gammas
Definition: DEventProcessor_BCAL_Shower.h:140
DEventProcessor_BCAL_Shower::L2_pathlength
Float_t L2_pathlength
Definition: DEventProcessor_BCAL_Shower.h:62
DBCALCluster::E
float E() const
Definition: DBCALCluster.h:41
DBCALCluster.h
DEventProcessor_BCAL_Shower::E3
Float_t E3
Definition: DEventProcessor_BCAL_Shower.h:69