RhoRhoPHOTOSUserTreeAnalysis.C
1//#include "UserTreeAnalysis.H" // remove if copied to user working directory
2#include <stdio.h>
3#include <assert.h>
4#include <math.h>
5#include <vector>
6#include <iostream>
7#include "MC4Vector.H"
8#include "HEPParticle.H"
9#include "TH1.h"
10#include "Setup.H"
11#include "TObjArray.h"
12#include "TMath.h"
13
14using namespace std;
15// Limits for particular user histograms
16int L[6] = { 5000000,5000000,2000000,2000000,1000000,1000000 };
17
18// very similar to MC_FillUserHistogram from Generate.cxx
19inline void fillUserHisto(char *name,double val, double weight=1.0,
20 double min=Setup::bin_min[0][0],
21 double max=Setup::bin_max[0][0])
22{
23 TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
24 if(!h){
25 h=new TH1D(name,name,Setup::nbins[0][0],min,max);
26 if(!h) return;
27 Setup::user_histograms->Add(h);
28 // printf("user histogram created %s\n", name);
29 }
30 h->Fill(val,weight);
31}
32
33double normalised_cross_product(double * v1, double * v2, double * result)
34{
35 result[0] = v1[1]*v2[2] - v1[2]*v2[1];
36 result[1] = v1[2]*v2[0] - v1[0]*v2[2];
37 result[2] = v1[0]*v2[1] - v1[1]*v2[0];
38
39 double normalisation = sqrt(result[0]*result[0]
40 +result[1]*result[1]
41 +result[2]*result[2]);
42
43 for(int i=0;i<3;i++)
44 result[i]=result[i]/normalisation;
45
46 return normalisation;
47}
48
49double dot_product(double *v1, double *v2)
50{
51 return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
52}
53
54double dot_product(MC4Vector v1, MC4Vector v2)
55{
56 return v1.GetX1()*v2.GetX1()+v1.GetX2()*v2.GetX2()+v1.GetX3()*v2.GetX3();
57}
58
59double magnitude(double *v)
60{
61 return sqrt(dot_product(v,v));
62}
63
64
65//see RhoRhoUserTreeAnalysis in TAUOLA
66/** Main function. This does not take any parameters. It assumes
67 the events are something -> tau+ tau-, then tau -> pi+/- pi0 nu */
68int RhoRhoPHOTOSUserTreeAnalysis(HEPParticle *mother,HEPParticleList *stableDaughters, int nparams, double *params)
69{
70 assert(mother!=0);
71 assert(stableDaughters!=0);
72
73 HEPParticleListIterator daughters(*stableDaughters);
74
75 //make temporary 4 vectors for the pions
76 double pi_plus[4]={0};
77 double pi_minus[4]={0};
78 double pi0_plus[4]={0};
79 double pi0_minus[4]={0};
80
81 MC4Vector d_pi0_plus;
82 MC4Vector d_pi0_minus;
83 MC4Vector d_pi_plus;
84 MC4Vector d_pi_minus;
85 MC4Vector d_gamma;
86 MC4Vector temp;
87
88 //make temporary variables to store the center of mass of
89 //the rho+ rho- pair.
90 double cm_px=0;
91 double cm_py=0;
92 double cm_pz=0;
93 double cm_e=0;
94
95 double photon_e = 0;
96
97 //loop over all daughters and sort them by type, filling the
98 //temporary variables.
99 for (HEPParticle *part=daughters.first(); part!=0; part=daughters.next())
100 {
101 temp = part->GetP4();
102 if(part->GetPDGId()!=16&&part->GetPDGId()!=-16)
103 {
104 //Get the center of mass
105 cm_px+=part->GetPx();
106 cm_py+=part->GetPy();
107 cm_pz+=part->GetPz();
108 cm_e+=part->GetE();
109 }
110
111 //Initialize the particle arrays
112 switch(part->GetPDGId())
113 {
114 case 211:
115 d_pi_plus.Set(&temp);
116 d_pi_plus.SetM(part->GetM());
117 break;
118 case -211:
119 d_pi_minus.Set(&temp);
120 d_pi_minus.SetM(part->GetM());
121 break;
122 case 111: //fill the pi0's randomly for the moment.
123 if(d_pi0_minus.GetX1()==0&&d_pi0_minus.GetX2()==0&&d_pi0_minus.GetX3()==0)
124 {
125 d_pi0_minus.Set(&temp);
126 d_pi0_minus.SetM(part->GetM());
127 }
128 else
129 {
130 d_pi0_plus.Set(&temp);
131 d_pi0_plus.SetM(part->GetM());
132 }
133 break;
134 case 22:
135 d_gamma.Set(&temp);
136 d_gamma.SetM(0);
137 // Only the hardest photon counts
138 if(photon_e>d_gamma.GetX0()) return 0;
139 photon_e=d_gamma.GetX0();
140 // Skip photons with energy < 10MeV
141 if(photon_e<0.01) photon_e=0;
142 break;
143 default:
144 break;
145 }
146 }
147
148 //Now check which pi0 is associated with
149 //which pi+/-. Use the angle to decide.
150 double costheta1 = dot_product(d_pi_plus,d_pi0_plus)/(d_pi_plus.Length()*d_pi0_plus.Length());
151 double costheta2 = dot_product(d_pi_minus,d_pi0_plus)/(d_pi_minus.Length()*d_pi0_plus.Length());
152
153
154 if(costheta1<costheta2) //if necessary, change the pi0 vectors
155 {
156 MC4Vector temp;
157 temp.Set(&d_pi0_plus);
158 temp.SetM(d_pi0_plus.GetM());
159 d_pi0_plus.Set(&d_pi0_minus);
160 d_pi0_plus.SetM(d_pi0_minus.GetM());
161 d_pi0_minus.Set(&temp);
162 d_pi0_minus.SetM(temp.GetM());
163 }
164
165 //Now boost everything into the rho+ rho- center of mass frame.
166 double cm_mass = sqrt(cm_e*cm_e-cm_px*cm_px-cm_py*cm_py-cm_pz*cm_pz);
167
168 d_pi0_plus.Boost(cm_px,cm_py,cm_pz,cm_e,cm_mass);
169 d_pi0_minus.Boost(cm_px,cm_py,cm_pz,cm_e,cm_mass);
170 d_pi_plus.Boost(cm_px,cm_py,cm_pz,cm_e,cm_mass);
171 d_pi_minus.Boost(cm_px,cm_py,cm_pz,cm_e,cm_mass);
172 d_gamma.Boost(cm_px,cm_py,cm_pz,cm_e,cm_mass);
173
174 pi0_plus[0]=d_pi0_plus.GetX1();
175 pi0_plus[1]=d_pi0_plus.GetX2();
176 pi0_plus[2]=d_pi0_plus.GetX3();
177 pi0_plus[3]=d_pi0_plus.GetM();
178
179 pi_plus[0]=d_pi_plus.GetX1();
180 pi_plus[1]=d_pi_plus.GetX2();
181 pi_plus[2]=d_pi_plus.GetX3();
182 pi_plus[3]=d_pi_plus.GetM();
183
184 pi0_minus[0]=d_pi0_minus.GetX1();
185 pi0_minus[1]=d_pi0_minus.GetX2();
186 pi0_minus[2]=d_pi0_minus.GetX3();
187 pi0_minus[3]=d_pi0_minus.GetM();
188
189 pi_minus[0]=d_pi_minus.GetX1();
190 pi_minus[1]=d_pi_minus.GetX2();
191 pi_minus[2]=d_pi_minus.GetX3();
192 pi_minus[3]=d_pi_minus.GetM();
193
194 /******* calculate acoplanarity (theta) *****/
195
196 //normal to the plane spanned by pi+ pi0
197 double n_plus[3];
198 normalised_cross_product(pi_plus,pi0_plus,n_plus);
199
200 //normal to the plane spanned by pi- pi0
201 double n_minus[3];
202 normalised_cross_product(pi_minus,pi0_minus,n_minus);
203
204 //get the angle
205 double theta = acos(dot_product(n_plus,n_minus));
206
207 //make theta go between 0 and 2 pi.
208 double pi_minus_n_plus = dot_product(pi_minus,n_plus)/magnitude(pi_minus);
209 if(pi_minus_n_plus>0)
210 theta=2*M_PI-theta;
211
212 /*********** calculate C/D reco (y1y2 in the paper) ***/
213
214 //boost vectors back to the lab frame
215 d_pi0_plus.Boost(-cm_px,-cm_py,-cm_pz,cm_e,cm_mass);
216 d_pi_plus.Boost(-cm_px,-cm_py,-cm_pz,cm_e,cm_mass);
217 d_pi0_minus.Boost(-cm_px,-cm_py,-cm_pz,cm_e,cm_mass);
218 d_pi_minus.Boost(-cm_px,-cm_py,-cm_pz,cm_e,cm_mass);
219
220 //construct effective tau 4 vectors (without neutrino)
221 double e_tau = 120.0/2.0;
222 double m_tau = 1.78;
223 double p_mag_tau = sqrt(e_tau*e_tau - m_tau*m_tau);
224
225 MC4Vector p_tau_plus = d_pi_plus + d_pi0_plus;
226 MC4Vector p_tau_minus = d_pi_minus + d_pi0_minus;
227
228 double norm_plus = p_mag_tau/p_tau_plus.Length();
229 double norm_minus = p_mag_tau/p_tau_minus.Length();
230
231 p_tau_plus.SetX0(e_tau);
232 p_tau_plus.SetX1(norm_plus*p_tau_plus.GetX1());
233 p_tau_plus.SetX2(norm_plus*p_tau_plus.GetX2());
234 p_tau_plus.SetX3(norm_plus*p_tau_plus.GetX3());
235 p_tau_plus.SetM(m_tau);
236
237 p_tau_minus.SetX0(e_tau);
238 p_tau_minus.SetX1(norm_minus*p_tau_minus.GetX1());
239 p_tau_minus.SetX2(norm_minus*p_tau_minus.GetX2());
240 p_tau_minus.SetX3(norm_minus*p_tau_minus.GetX3());
241 p_tau_minus.SetM(m_tau);
242
243 //boost to the (reco) tau's frames
244 d_pi0_plus.BoostP(p_tau_plus);
245 d_pi_plus.BoostP(p_tau_plus);
246 d_pi0_minus.BoostP(p_tau_minus);
247 d_pi_minus.BoostP(p_tau_minus);
248
249 //calculate y1 and y2
250 double y1=(d_pi_plus.GetX0()-d_pi0_plus.GetX0())/(d_pi_plus.GetX0()+d_pi0_plus.GetX0());
251 double y2=(d_pi_minus.GetX0()-d_pi0_minus.GetX0())/(d_pi_minus.GetX0()+d_pi0_minus.GetX0());
252
253 //make seperate plots for different photon energy ranges
254 if(photon_e==0)
255 {
256 if(y1*y2>0)
257 {
258 if(L[0]<=0) return 0;
259 L[0]--;
260 char name[] = "acoplanarity-angle-C-no-photon";
261 fillUserHisto(name,theta,1.0,0,2*M_PI);
262 }
263 else
264 {
265 if(L[1]<=0) return 0;
266 L[1]--;
267 char name[] = "acoplanarity-angle-D-no-photon";
268 fillUserHisto(name,theta,1.0,0,2*M_PI);
269 }
270
271 }
272 if(photon_e>0&&photon_e<1.0)
273 {
274 if(y1*y2>0)
275 {
276 if(L[2]<=0) return 0;
277 L[2]--;
278 char name[] = "acoplanarity-angle-C-photon-up-to-1-GeV";
279 fillUserHisto(name,theta,1.0,0,2*M_PI);
280 }
281 else
282 {
283 if(L[3]<=0) return 0;
284 L[3]--;
285 char name[] = "acoplanarity-angle-D-photon-up-to-1-GeV";
286 fillUserHisto(name,theta,1.0,0,2*M_PI);
287 }
288 }
289 if(photon_e>1.0)
290 {
291 if(y1*y2>0)
292 {
293 if(L[4]<=0) return 0;
294 L[4]--;
295 char name[] = "acoplanarity-angle-C-photon-over-1-GeV";
296 fillUserHisto(name,theta,1.0,0,2*M_PI);
297 }
298 else
299 {
300 if(L[5]<=0) return 0;
301 L[5]--;
302 char name[] = "acoplanarity-angle-D-photon-over-1-GeV";
303 fillUserHisto(name,theta,1.0,0,2*M_PI);
304 }
305 }
306
307 return 0;
308};
309