new features for pbarp test tools

Examples/pbarp/pbarpEnv.cc

@@ -174,6 +174,27 @@ void pbarpEnv::setup(pbarpParser& thePbarpParser){
     _histMassSystems.push_back(currentStringVec);
   }
 
+  std::vector<std::string> theHistAngleNames=thePbarpParser.histAngleNames();
+  //fill vector histMassSystems
+  for ( itStr = theHistAngleNames.begin(); itStr != theHistAngleNames.end(); ++itStr){
+    std::stringstream stringStr;
+    stringStr << (*itStr);
+    
+    std::string tmpName;
+    std::vector<std::string> currentStringDecVec;
+    std::vector<std::string> currentStringMotherVec;
+    bool isDecParticle=true;
+    while(stringStr >> tmpName){
+      if(tmpName=="from") {
+	isDecParticle=false;
+	continue;
+      }
+      if(isDecParticle) currentStringDecVec.push_back(tmpName);
+      else currentStringMotherVec.push_back(tmpName);
+    }
+    boost::shared_ptr<angleHistData> currentAngleHistData(new angleHistData(currentStringMotherVec, currentStringDecVec)); 
+    _angleHistDataVec.push_back(currentAngleHistData);
+  }
 }
 
 

Examples/pbarp/pbarpEnv.hh

@@ -10,12 +10,14 @@
 
 #include "PwaUtils/DataUtils.hh"
 #include "Examples/pbarp/pbarpParser.hh"
+#include "Examples/pbarp/pbarpHist.hh"
 
 class Particle;
 class ParticleTable;
 class pbarpReaction;
 class IsobarDecayList;
 class pbarpEventList;
+//class pbarpHist;
 
 class pbarpEnv {
 
@@ -34,7 +36,7 @@ public:
   boost::shared_ptr<IsobarDecayList> prodDecayList() {return _prodDecList;}
   boost::shared_ptr<pbarpReaction> reaction() {return _pbarpReaction;} 
   std::vector<std::vector<std::string> >& histMassSystems() {return _histMassSystems;} 
-
+  std::vector<boost::shared_ptr<angleHistData> >& angleHistDataVec() {return _angleHistDataVec;} 
 protected:  
 
   pbarpEnv();
@@ -50,4 +52,5 @@ protected:
   boost::shared_ptr<IsobarDecayList> _prodDecList;
   boost::shared_ptr<pbarpReaction> _pbarpReaction;
   std::vector<std::vector<std::string> > _histMassSystems;
+  std::vector<boost::shared_ptr<angleHistData> > _angleHistDataVec;
 };

Examples/pbarp/pbarpEventList.cc

@@ -66,7 +66,7 @@ void pbarpEventList::read4Vecs(EventList& evtList, std::vector<EvtDataNew*>& the
 
 
     EvtDataNew* evtData=new EvtDataNew();
-
+    evtData->FourVecsString["all"]=V4_all_lab;
     //cache 4 vectors of inital state particles
     std::map<std::string, Vector4<float> >::iterator it4VecMap;
     for (it4VecMap=particle4VecMap.begin(); it4VecMap!=particle4VecMap.end(); ++it4VecMap){

Examples/pbarp/pbarpHist.cc

@@ -36,6 +36,7 @@ pbarpHist::pbarpHist(boost::shared_ptr<AbsLhNew> theLh, fitParamsNew& theFitPara
     {
       double weight = (*it)->evtWeight;
       fillMassHists((*it), weight, _massDataHistMap);
+      fillThetaHists((*it), weight, _thetaDataHistMap);
       ++it;
     }
 
@@ -45,9 +46,11 @@ pbarpHist::pbarpHist(boost::shared_ptr<AbsLhNew> theLh, fitParamsNew& theFitPara
     {
       double evtWeight = (*it)->evtWeight;
       fillMassHists((*it), evtWeight, _massMcHistMap);
+      fillThetaHists((*it), evtWeight, _thetaMcHistMap);
 
       double fitWeight= theLh->calcEvtIntensity( (*it), theFitParams );
       fillMassHists((*it), evtWeight*fitWeight, _massFitHistMap);
+      fillThetaHists((*it), evtWeight*fitWeight, _thetaFitHistMap);
       ++it;
     }
 
@@ -69,6 +72,11 @@ pbarpHist::pbarpHist(boost::shared_ptr<AbsLhNew> theLh, fitParamsNew& theFitPara
   for(itMassMap= _massFitHistMap.begin(); itMassMap!= _massFitHistMap.end(); ++itMassMap){
     itMassMap->second->Scale(scaleFactor);
   }
+
+  std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess >::iterator itAngleMap;
+  for(itAngleMap= _thetaFitHistMap.begin(); itAngleMap!=_thetaFitHistMap.end(); ++itAngleMap){
+    itAngleMap->second->Scale(scaleFactor);
+  }
 }
 
 pbarpHist::~pbarpHist(){
@@ -106,6 +114,33 @@ void pbarpHist::initRootStuff(){
     _massFitHistMap[tmpMassHistData]=currentMassFitHist;
   }
 
+  std::vector<boost::shared_ptr<angleHistData> > angleHistDataVec=pbarpEnv::instance()->angleHistDataVec();
+
+  std::vector<boost::shared_ptr<angleHistData> >::iterator itAngleVec;
+  for (itAngleVec=angleHistDataVec.begin(); itAngleVec!=angleHistDataVec.end(); ++itAngleVec){
+    std::string tmpBaseName= (*itAngleVec)->_name;
+    boost::replace_all(tmpBaseName,"+","p");
+    boost::replace_all(tmpBaseName,"-","m");
+    std::string histName="data"+tmpBaseName;
+    std::string histDescription = "cos#Theta(" +(*itAngleVec)->_name + ") (data)";
+    TH1F* currentThetaAngleDataHist=new TH1F(histName.c_str(), histDescription.c_str(), 100., -1., 1.);
+    currentThetaAngleDataHist->Sumw2();
+    _thetaDataHistMap[*itAngleVec]=currentThetaAngleDataHist;
+
+    histName="MC"+tmpBaseName;
+    histDescription = "cos#Theta(" +(*itAngleVec)->_name + ") (MC)";
+    TH1F* currentThetaAngleMcHist=new TH1F(histName.c_str(), histDescription.c_str(), 100., -1., 1.);
+    currentThetaAngleMcHist->Sumw2();
+    _thetaMcHistMap[*itAngleVec]=currentThetaAngleMcHist;
+
+    histName="Fit"+tmpBaseName;
+    histDescription = "cos#Theta(" +(*itAngleVec)->_name + ") (fit)";
+    TH1F* currentThetaAngleFitHist=new TH1F(histName.c_str(), histDescription.c_str(), 100., -1., 1.);
+    currentThetaAngleFitHist->Sumw2();
+    _thetaFitHistMap[*itAngleVec]=currentThetaAngleFitHist;
+  } 
+
+ std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess > _thetaDataHistMap;
 }
 
 void pbarpHist::fillMassHists(EvtDataNew* theData, double weight, std::map<boost::shared_ptr<massHistData>, TH1F*, pawian::Collection::SharedPtrLess >& toFill){
@@ -125,3 +160,36 @@ void pbarpHist::fillMassHists(EvtDataNew* theData, double weight, std::map<boost
   }
   
 }
+
+void pbarpHist::fillThetaHists(EvtDataNew* theData, double weight, std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess >& toFill){
+
+  std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess >::iterator it;
+  for(it= toFill.begin(); it!= toFill.end(); ++it){
+    Vector4<double> combinedDec4Vec(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>::iterator itStr;
+    for(itStr=decNames.begin(); itStr!=decNames.end(); ++itStr){
+      Vector4<double> tmp4vec=theData->FourVecsString[*itStr];
+      combinedDec4Vec+=tmp4vec;
+    }
+
+    std::vector<std::string> motherNames=it->first->_motherPNames;
+
+    for(itStr=motherNames.begin(); itStr!=motherNames.end(); ++itStr){
+      Vector4<double> tmp4vec=theData->FourVecsString[*itStr];
+      combinedMother4Vec+=tmp4vec;
+    }
+
+    Vector4<double>  result4Vec(0.,0.,0.,0.);
+    result4Vec=helicityVec(all4Vec, combinedMother4Vec, combinedDec4Vec);
+    // if(all4Vec==combinedMother4Vec){
+    //   result4Vec=combinedDec4Vec;
+    //   result4Vec.Boost(all4Vec);
+    // } 
+    // else result4Vec=helicityVec(all4Vec, combinedMother4Vec, combinedDec4Vec);
+    it->second->Fill( result4Vec.CosTheta(), weight);     
+  }
+}

Examples/pbarp/pbarpHist.hh

@@ -55,7 +55,41 @@ struct massHistData {
   }  
 };
 
+struct angleHistData {
+  angleHistData(std::vector<std::string>& motherPNames, std::vector<std::string>& decPNames) :
+    _name("")
+    ,_motherPNames(motherPNames)
+    ,_decPNames(decPNames)
+  {
 
+    std::vector<std::string>::iterator it;
+    for(it=decPNames.begin(); it!=decPNames.end(); ++it){
+      _name+=(*it);
+    }
+
+    _name+="_Heli";
+
+    for(it=motherPNames.begin(); it!=motherPNames.end(); ++it){
+      _name+=(*it);
+    }
+  }
+
+  std::string _name;
+  std::vector<std::string> _motherPNames; 
+  std::vector<std::string> _decPNames; 
+
+  virtual bool operator==(const angleHistData& compare) const {
+    bool result=false;
+    if ( _name==compare._name) result=true;
+    return result;
+  }
+
+ virtual bool operator<(const angleHistData& compare) const {
+   bool result=false;
+   if(_name < compare._name) result=true; 
+    return result; 
+  }  
+};
 
 class TFile;
 //class TH2F;
@@ -73,11 +107,18 @@ protected:
 private:
   void fillMassHists(EvtDataNew* theData, double weight, std::map<boost::shared_ptr<massHistData>, TH1F*, pawian::Collection::SharedPtrLess >& toFill);
 
+  void fillThetaHists(EvtDataNew* theData, double weight, std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess >& toFill);
+
  TFile* _theTFile;
 
  std::map<boost::shared_ptr<massHistData>, TH1F*, pawian::Collection::SharedPtrLess > _massDataHistMap;
  std::map<boost::shared_ptr<massHistData>, TH1F*, pawian::Collection::SharedPtrLess > _massMcHistMap;
  std::map<boost::shared_ptr<massHistData>, TH1F*, pawian::Collection::SharedPtrLess > _massFitHistMap;
+
+ std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess > _thetaDataHistMap;
+ std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess > _thetaMcHistMap;
+ std::map<boost::shared_ptr<angleHistData>, TH1F*, pawian::Collection::SharedPtrLess > _thetaFitHistMap;
+
   void initRootStuff();
 };
 

Examples/pbarp/pbarpParser.cc

@@ -28,6 +28,7 @@ pbarpParser::pbarpParser(int argc,char **argv):
     ("decay",po::value< vector<string> >(&_decaySystem),  "decay: mother and pair of decay particles")
     ("replaceParamSuffix",po::value< vector<string> >(&_replaceParSuffix),  "replace suffix for fit parameter name")
     ("histMass",po::value< vector<string> >(&_histMass),  "histograms inv mass for the selected final state paricles")
+    ("histAngles",po::value< vector<string> >(&_histAngles),  "histograms decay angles")
     ;
   _config->add(config);
   
@@ -64,5 +65,10 @@ bool pbarpParser::parseCommandLine(int argc, char **argv)
   for (it=_histMass.begin(); it!=_histMass.end();++it){
     std::cout << (*it) << "\n";
   }
+
+  std::cout << "\nhistograms decay angles for systems" << std::endl;
+  for (it=_histAngles.begin(); it!=_histAngles.end();++it){
+    std::cout << (*it) << "\n";
+  }
   return true;
 }

Examples/pbarp/pbarpParser.hh

@@ -20,6 +20,7 @@ class pbarpParser : public ParserBase
   const std::vector<std::string>& decaySystem() const { return _decaySystem; }
   const std::vector<std::string>& replaceSuffixNames() const { return _replaceParSuffix; }
   const std::vector<std::string>& histMassNames() const { return _histMass; }
+  const std::vector<std::string>& histAngleNames() const { return _histAngles;}
 
 protected:
   virtual bool parseCommandLine(int argc,char **argv); 
@@ -30,6 +31,7 @@ protected:
   std::vector<std::string> _decaySystem;
   std::vector<std::string> _replaceParSuffix;
   std::vector<std::string> _histMass;
+  std::vector<std::string> _histAngles;
 };