Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include <getopt.h>
#include <fstream>
#include <string>
#include "Examples/pbarpToOmegaPiLS/OmegaTo3PiEventListLS.hh"
#include "Event/EventList.hh"
#include "PwaUtils/KinUtils.hh"
#include "Event/Event.hh"
#include "ErrLogger/ErrLogger.hh"
OmegaTo3PiEventListLS::OmegaTo3PiEventListLS(EventList& evtListData, EventList& evtListMc, unsigned jmax, unsigned pbarmom):
AbsOmegaPiEventListLS(evtListData, evtListMc, jmax, pbarmom)
{
read4Vecs(evtListData, _dataList);
read4Vecs(evtListMc, _mcList);
}
OmegaTo3PiEventListLS::~OmegaTo3PiEventListLS()
{
}
void OmegaTo3PiEventListLS::read4Vecs(EventList& evtList, std::vector<OmPiEvtDataLS*>& omPiEvtList)
{
Event* anEvent;
int evtCount = 0;
bool bPbarmomOnce=false;
while ((anEvent = evtList.nextEvent()) != 0){
if ( evtCount%1000 == 0 ) Info << "4vec calculation for event " << evtCount << endmsg;
Vector4<float> cm_4V( *(anEvent->p4(0))+ *(anEvent->p4(1)) + *(anEvent->p4(2)) + *(anEvent->p4(3)));
Vector4<float> pi0Omega_4V;
Vector4<float> pipOmega_4V;
Vector4<float> pimOmega_4V;
Vector4<float> pi0Rec_4V;
pi0Rec_4V=*(anEvent->p4(0));
pi0Omega_4V=*(anEvent->p4(1));
pipOmega_4V=*(anEvent->p4(2));
pimOmega_4V=*(anEvent->p4(3));
Vector4<float> omega_4V(pi0Omega_4V+pipOmega_4V+pimOmega_4V);
Vector4<float> omega_cm_4V(omega_4V);
omega_cm_4V.Boost(cm_4V);
Vector4<float> pi0Rec_cm_4V(pi0Rec_4V);
pi0Rec_cm_4V.Boost(cm_4V);
if (!bPbarmomOnce)
{
float E=sqrt(pow(938.272,2)+pow(double(_pbarmom),2));
Vector4<float> pbarmom;
pbarmom.SetP4(E,0.,0.,_pbarmom);
pbarmom.Boost(cm_4V);
cout << endl << "pbarmom-cmc =" << pbarmom.P() << " MeV/c" << endl;
bPbarmomOnce=true;
}
Vector4<float> omega_4V2(pi0Rec_4V+pipOmega_4V+pimOmega_4V);
Vector4<float> omega_cm_4V2(omega_4V2);
omega_cm_4V2.Boost(cm_4V);
Vector4<float> pi0HeliOmega4V = helicityVec(cm_4V, omega_4V, pi0Omega_4V);
Vector4<float> pipHeliOmega4V = helicityVec(cm_4V, omega_4V, pipOmega_4V);
Vector4<float> pimHeliOmega4V = helicityVec(cm_4V, omega_4V, pimOmega_4V);
Vector4<float> pi0HeliOmega4V2 = helicityVec(cm_4V, omega_4V2, pi0Rec_4V);
Vector4<float> pipHeliOmega4V2 = helicityVec(cm_4V, omega_4V2, pipOmega_4V);
Vector4<float> pimHeliOmega4V2 = helicityVec(cm_4V, omega_4V2, pimOmega_4V);
Vector4<float> normOmegaDecHeliOmega_4V(0.5*(pi0HeliOmega4V.T()+pipHeliOmega4V.T()+pimHeliOmega4V.T()),
pimHeliOmega4V.Y()*pipHeliOmega4V.Z()-pimHeliOmega4V.Z()*pipHeliOmega4V.Y(),
pimHeliOmega4V.Z()*pipHeliOmega4V.X()-pimHeliOmega4V.X()*pipHeliOmega4V.Z(),
pimHeliOmega4V.X()*pipHeliOmega4V.Y()-pimHeliOmega4V.Y()*pipHeliOmega4V.X());
Vector4<float> normOmegaDecHeliOmega_4V2(0.5*(pi0HeliOmega4V2.T()+pipHeliOmega4V2.T()+pimHeliOmega4V2.T()),
pimHeliOmega4V2.Y()*pipHeliOmega4V2.Z()-pimHeliOmega4V2.Z()*pipHeliOmega4V2.Y(),
pimHeliOmega4V2.Z()*pipHeliOmega4V2.X()-pimHeliOmega4V2.X()*pipHeliOmega4V2.Z(),
pimHeliOmega4V2.X()*pipHeliOmega4V2.Y()-pimHeliOmega4V2.Y()*pipHeliOmega4V2.X());
float theQ=pipHeliOmega4V.E()-pipHeliOmega4V.M()+pimHeliOmega4V.E()-pimHeliOmega4V.M()+pi0HeliOmega4V.E()-pi0HeliOmega4V.M();
float lambdaNorm=theQ*theQ*(theQ*theQ/108.+pimHeliOmega4V.M()*theQ/9.+pimHeliOmega4V.M()*pimHeliOmega4V.M()/3.);
float lambda=normOmegaDecHeliOmega_4V.P()*normOmegaDecHeliOmega_4V.P()/lambdaNorm;
OmPiEvtDataLS* theOmPiEvtData=new OmPiEvtDataLS();
theOmPiEvtData->cm_4Vec=cm_4V;
theOmPiEvtData->omegaHeliCm4Vec=omega_cm_4V;
theOmPiEvtData->omegaHeliCm4Vec2=omega_cm_4V2;
theOmPiEvtData->pi0RecHeliCm4Vec=pi0Rec_cm_4V;
theOmPiEvtData->pi0HeliOmega4Vec=normOmegaDecHeliOmega_4V;
theOmPiEvtData->pi0HeliOmega4Vec2=normOmegaDecHeliOmega_4V2;
theOmPiEvtData->cosPi0HeliOmega4Vec=normOmegaDecHeliOmega_4V.CosTheta();
theOmPiEvtData->pipHeliOmega4Vec=pipHeliOmega4V;
theOmPiEvtData->pimHeliOmega4Vec=pimHeliOmega4V;
theOmPiEvtData->pippimCrossHeliOmega4Vec=normOmegaDecHeliOmega_4V;
theOmPiEvtData->lambda=lambda;
theOmPiEvtData->eventWeight = anEvent->Weight();
for (Spin j=0; j<=_jmax; j++){
for (Spin M=-1; M<=1; M++){
if (fabs(M)>j) continue;
for (Spin lam=-1; lam<=1; lam++){
if (fabs(lam)>j) continue;
Bertram Kopf
committed
theOmPiEvtData->Dfp[j][M][lam]=Wigner_D(0.,omega_cm_4V.Theta(),0,j,M,lam);
theOmPiEvtData->Dfp2[j][M][lam]=Wigner_D(0.,omega_cm_4V2.Theta(),0,j,M,lam);
}
}
}
Spin omegaSpin=1;
Wigner_D(omegaSpin, normOmegaDecHeliOmega_4V.Phi(), normOmegaDecHeliOmega_4V.Theta(),0, theOmPiEvtData->Dfd);
Wigner_D(omegaSpin, normOmegaDecHeliOmega_4V2.Phi(), normOmegaDecHeliOmega_4V2.Theta(),0, theOmPiEvtData->Dfd2);