Skip to content
Snippets Groups Projects
Commit 2c4eb8c8 authored by Julian Pychy's avatar Julian Pychy
Browse files

added angle histogramming for 3 particle decay

parent 3a2783ee
No related branches found
No related tags found
No related merge requests found
......@@ -124,17 +124,26 @@ void pbarpEnv::setup(pbarpParser* thePbarpParser){
std::string tmpName;
std::vector<std::string> currentStringDecVec;
std::vector<std::string> currentStringDecVec2;
std::vector<std::string> currentStringMotherVec;
bool isDecParticle=true;
bool isFirstDecParticle=true;
short nBodyDecay=2;
while(stringStr >> tmpName){
if(tmpName=="from") {
isDecParticle=false;
continue;
}
if(isDecParticle) currentStringDecVec.push_back(tmpName);
else if(tmpName=="and") {
isFirstDecParticle=false;
nBodyDecay=3;
continue;
}
if(isDecParticle && isFirstDecParticle) currentStringDecVec.push_back(tmpName);
else if(isDecParticle && !isFirstDecParticle) currentStringDecVec2.push_back(tmpName);
else currentStringMotherVec.push_back(tmpName);
}
boost::shared_ptr<angleHistData> currentAngleHistData(new angleHistData(currentStringMotherVec, currentStringDecVec));
boost::shared_ptr<angleHistData> currentAngleHistData(new angleHistData(currentStringMotherVec, currentStringDecVec, currentStringDecVec2, nBodyDecay));
_angleHistDataVec.push_back(currentAngleHistData);
}
......
......@@ -212,16 +212,25 @@ void pbarpHist::fillAngleHists(EvtData* theData, double weight, std::map<boost::
std::map<boost::shared_ptr<angleHistData>, std::pair<TH1F*, TH1F*>, pawian::Collection::SharedPtrLess >::iterator it;
for(it= toFill.begin(); it!= toFill.end(); ++it){
short nBodyDecay = it->first->_nBodyDecay;
Vector4<double> combinedDec4Vec(0.,0.,0.,0.);
Vector4<double> combinedDec4Vec2(0.,0.,0.,0.);
Vector4<double> combinedMother4Vec(0.,0.,0.,0.);
Vector4<double> all4Vec=theData->FourVecsString["all"];
std::vector<std::string> decNames=it->first->_decPNames;
std::vector<std::string> decNames2=it->first->_decPNames2;
std::vector<std::string>::iterator itStr;
for(itStr=decNames.begin(); itStr!=decNames.end(); ++itStr){
Vector4<double> tmp4vec=theData->FourVecsString[*itStr];
combinedDec4Vec+=tmp4vec;
}
for(itStr=decNames2.begin(); itStr!=decNames2.end(); ++itStr){
Vector4<double> tmp4vec=theData->FourVecsString[*itStr];
combinedDec4Vec2+=tmp4vec;
}
std::vector<std::string> motherNames=it->first->_motherPNames;
......@@ -231,6 +240,8 @@ void pbarpHist::fillAngleHists(EvtData* theData, double weight, std::map<boost::
}
Vector4<double> result4Vec(0.,0.,0.,0.);
Vector4<double> result4Vec2(0.,0.,0.,0.);
if( fabs(all4Vec.E()-combinedMother4Vec.E()) < 1e-5
&& fabs(all4Vec.Px()-combinedMother4Vec.Px()) < 1e-5
&& fabs(all4Vec.Py()-combinedMother4Vec.Py()) < 1e-5
......@@ -238,9 +249,26 @@ void pbarpHist::fillAngleHists(EvtData* theData, double weight, std::map<boost::
result4Vec=combinedDec4Vec;
result4Vec.Boost(all4Vec);
}
else result4Vec=helicityVec(all4Vec, combinedMother4Vec, combinedDec4Vec);
else{
result4Vec=helicityVec(all4Vec, combinedMother4Vec, combinedDec4Vec);
if(nBodyDecay==3)
result4Vec2=helicityVec(all4Vec, combinedMother4Vec, combinedDec4Vec2);
}
it->second.first->Fill( result4Vec.CosTheta(), weight);
it->second.second->Fill( result4Vec.Phi(), weight);
if(nBodyDecay == 2){
it->second.first->Fill( result4Vec.CosTheta(), weight);
it->second.second->Fill( result4Vec.Phi(), weight);
}
else if(nBodyDecay == 3){
Vector4<double> result4VecPart1 = result4Vec;
Vector4<double> result4VecPart2 = result4Vec2;
Vector4<float> normVec(0,
result4VecPart1.Y()*result4VecPart2.Z()-result4VecPart1.Z()*result4VecPart2.Y(),
result4VecPart1.Z()*result4VecPart2.X()-result4VecPart1.X()*result4VecPart2.Z(),
result4VecPart1.X()*result4VecPart2.Y()-result4VecPart1.Y()*result4VecPart2.X());
it->second.first->Fill( normVec.CosTheta(), weight);
it->second.second->Fill( normVec.Phi(), weight);
}
}
}
......@@ -59,10 +59,12 @@ struct massHistData {
};
struct angleHistData {
angleHistData(std::vector<std::string>& motherPNames, std::vector<std::string>& decPNames) :
_name("")
angleHistData(std::vector<std::string>& motherPNames, std::vector<std::string>& decPNames, std::vector<std::string>& decPNames2, short nBodyDecay) :
_nBodyDecay(nBodyDecay)
,_name("")
,_motherPNames(motherPNames)
,_decPNames(decPNames)
,_decPNames2(decPNames2)
{
std::vector<std::string>::iterator it;
......@@ -70,6 +72,13 @@ struct angleHistData {
_name+=(*it);
}
if(nBodyDecay == 3)
_name+="AND";
for(it=decPNames2.begin(); it!=decPNames2.end(); ++it){
_name+=(*it);
}
_name+="_Heli";
for(it=motherPNames.begin(); it!=motherPNames.end(); ++it){
......@@ -77,9 +86,11 @@ struct angleHistData {
}
}
short _nBodyDecay;
std::string _name;
std::vector<std::string> _motherPNames;
std::vector<std::string> _decPNames;
std::vector<std::string> _motherPNames;
std::vector<std::string> _decPNames;
std::vector<std::string> _decPNames2;
virtual bool operator==(const angleHistData& compare) const {
bool result=false;
......
......@@ -25,12 +25,12 @@ spinDensityHist::spinDensityHist(boost::shared_ptr<AbsLh> theLh, fitParams& theF
,_maxEvents(2000)
,_theLh(theLh)
{
_dataList=_theLh->getEventList()->getMcVecs();
_dataList=_theLh->getEventList()->getMcVecs();
boost::dynamic_pointer_cast<pbarpBaseLh>(_theLh)->updateFitParams(theFitParams);
std::stringstream spinDensityRootFileName;
spinDensityRootFileName << "./spinDensity" << pbarpEnv::instance()->outputFileNameSuffix() << ".root";
_spinDensityRootFile = new TFile(spinDensityRootFileName.str().c_str(), "recreate");
std::vector<string> spinDensityParticles = pbarpEnv::instance()->spinDensityNames();
......@@ -39,7 +39,7 @@ spinDensityHist::spinDensityHist(boost::shared_ptr<AbsLh> theLh, fitParams& theF
for(it=spinDensityParticles.begin(); it!=spinDensityParticles.end(); ++it){
Info << "Calculating spin density matrix for particle " << (*it) << endmsg;
calcSpinDensityMatrix(*it);
}
}
}
......@@ -52,7 +52,7 @@ spinDensityHist::~spinDensityHist(){
void spinDensityHist::calcSpinDensityMatrix(std::string& particleName){
int J = pbarpEnv::instance()->particleTable()->particle(particleName)->J();
if(J!=1){
......@@ -79,7 +79,7 @@ void spinDensityHist::calcSpinDensityMatrix(std::string& particleName){
Info << "Calculating Element (" << M1 << ", " << M2 << ")" << endmsg;
calcSpinDensityMatrixElement(particleName, M1, M2);
}
}
}
}
......@@ -101,7 +101,7 @@ void spinDensityHist::calcSpinDensityMatrixElement(std::string& particleName, Sp
int eventsRead=0;
std::vector<EvtData*>::iterator it;
for(it = _dataList.begin(); it != _dataList.end(); ++it){
complex<double> tempSpinDensity =
complex<double> tempSpinDensity =
boost::dynamic_pointer_cast<pbarpBaseLh>(_theLh)->calcSpinDensity(M1, M2, particleName, *it);
fillHistogram(particleName, newHistoReal, *it, tempSpinDensity.real());
......@@ -123,18 +123,18 @@ void spinDensityHist::calcSpinDensityMatrixElement(std::string& particleName, Sp
void spinDensityHist::fillHistogram(std::string& particleName, TH1F* theHisto,
void spinDensityHist::fillHistogram(std::string& particleName, TH1F* theHisto,
EvtData* theData, double spinDensityValue){
boost::shared_ptr<AbsDecayList> absDecayList = pbarpEnv::instance()->absDecayList();
std::vector<boost::shared_ptr<AbsDecay> > absDecays = absDecayList->getList();
std::vector<boost::shared_ptr<AbsDecay> >::iterator it;
Vector4<double> all4Vec=theData->FourVecsString["all"];
Vector4<double> all4Vec=theData->FourVecsString["all"];
Vector4<double> particle4Vec(0.,0.,0.,0.);
for(it=absDecays.begin(); it!=absDecays.end(); ++it){
if( (*it)->motherPart()->name() == particleName){
std::vector<Particle*> fsParticles = (*it)->finalStateParticles();
std::vector<Particle*>::iterator it2;
......@@ -148,8 +148,8 @@ void spinDensityHist::fillHistogram(std::string& particleName, TH1F* theHisto,
particle4Vec.Boost(all4Vec);
break;
}
}
}
theHisto->Fill(particle4Vec.CosTheta(), spinDensityValue * theData->evtWeight);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment