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 
14 using namespace std;
15 // Limits for particular user histograms
16 int L[6] = { 5000000,5000000,2000000,2000000,1000000,1000000 };
17 
18 // very similar to MC_FillUserHistogram from Generate.cxx
19 inline 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 
33 double 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 
49 double dot_product(double *v1, double *v2)
50 {
51  return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
52 }
53 
54 double dot_product(MC4Vector v1, MC4Vector v2)
55 {
56  return v1.GetX1()*v2.GetX1()+v1.GetX2()*v2.GetX2()+v1.GetX3()*v2.GetX3();
57 }
58 
59 double 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 */
68 int 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 
Generated on Tue Mar 12 2024 23:11:41 for C++InterfacetoPHOTOS by   doxygen 1.8.7