//************************************************************************// // // // Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) // // Julian Pychy (julian@ep1.rub.de) // // - Ruhr-Universität Bochum // // // // This file is part of Pawian. // // // // Pawian is free software: you can redistribute it and/or modify // // it under the terms of the GNU General Public License as published by // // the Free Software Foundation, either version 3 of the License, or // // (at your option) any later version. // // // // Pawian is distributed in the hope that it will be useful, // // but WITHOUT ANY WARRANTY; without even the implied warranty of // // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // // GNU General Public License for more details. // // // // You should have received a copy of the GNU General Public License // // along with Pawian. If not, see <http://www.gnu.org/licenses/>. // // // //************************************************************************// // EvtDataBaseList class definition file. -*- C++ -*- // Copyright 2012 Bertram Kopf #include <getopt.h> #include "PwaUtils/EvtDataBaseList.hh" #include "PwaUtils/KinUtils.hh" #include "PwaUtils/AbsDecay.hh" #include "PwaUtils/AbsDecayList.hh" #include "PwaUtils/AbsDynamics.hh" #include "PwaUtils/DynRegistry.hh" #include "PwaUtils/GlobalEnv.hh" #include "Event/EventList.hh" #include "Event/Event.hh" #include "Particle/Particle.hh" #include "Particle/ParticleTable.hh" #include "Utils/PawianCollectionUtils.hh" #include "ConfigParser/ParserBase.hh" #include "ErrLogger/ErrLogger.hh" EvtDataBaseList::EvtDataBaseList(ChannelID channelID) : _channelID(channelID) ,_noOfWeightedDataEvts(0.) ,_noOfWeightedMcEvts(0.) ,_alreadyRead(false) ,_finalStateParticles(GlobalEnv::instance()->Channel(_channelID)->finalStateParticles()) { } EvtDataBaseList::~EvtDataBaseList() { } void EvtDataBaseList::read(EventList& evtListData, EventList& evtListMc){ if(_alreadyRead){ Alert << "4 vectors already read " << endmsg; // << endmsg; exit(1); } read4Vecs(evtListData, _evtDataList, _noOfWeightedDataEvts, evtListData.size(), 0 ); read4Vecs(evtListMc, _mcDataList, _noOfWeightedMcEvts, evtListMc.size(), evtListData.size() ); _alreadyRead=true; } void EvtDataBaseList::read4Vecs(EventList& evtList, std::vector<EvtData*>& theEvtList, double& evtWeightSum, int maxEvts, int startNo){ Event* anEvent; int evtCount = 0; while ((anEvent = evtList.nextEvent())){ if (evtCount>= maxEvts) break; if (evtCount%500 == 0) Info << "4vec calculation for event " << evtCount ; // << endmsg; EvtData* currentEvt=convertEvent(anEvent, startNo+evtCount); if (evtCount%10000 == 0){ Vector4<double> V4_all_lab=currentEvt->FourVecsString.at("all"); Info << "4vec all in lab system" << "\n" << " px: " << V4_all_lab.Px() <<"\t" << " py: " << V4_all_lab.Py() <<"\t" << " pz: " << V4_all_lab.Pz() <<"\t" << " e : " << V4_all_lab.E() << "\t" << " m : " << V4_all_lab.M() ; // << endmsg; } theEvtList.push_back(currentEvt); evtWeightSum += anEvent->Weight(); ++evtCount; } } EvtData* EvtDataBaseList::convertEvent(Event* theEvent, int evtNo){ Vector4<double> V4_all_lab(0.,0.,0.,0.); std::vector< Vector4<double> > finalState4Vecs; std::map<std::string, Vector4<double> > particle4VecMap; std::vector<Particle*>::iterator itPart; int counter=0; for (itPart=_finalStateParticles.begin(); itPart != _finalStateParticles.end(); ++itPart){ Vector4<float> current4VecFloat=*(theEvent->p4(counter)); Vector4<double> current4Vec(current4VecFloat.E(), current4VecFloat.Px(), current4VecFloat.Py(), current4VecFloat.Pz()); if(*(*itPart)==*(GlobalEnv::instance()->particleTable()->particle("photon"))){ //set mass=0 current4Vec.SetP4(current4Vec.P(), current4Vec.Px(), current4Vec.Py(), current4Vec.Pz()); } // if( (current4Vec.M() != current4Vec.M()) || current4Vec.M()<1e-4) current4Vec.SetP4(sqrt(current4VecFloat.Px()*current4VecFloat.Px()+current4VecFloat.Py()*current4VecFloat.Py()+current4VecFloat.Pz()*current4VecFloat.Pz()), current4VecFloat.Px(), current4VecFloat.Py(), current4VecFloat.Pz()); finalState4Vecs.push_back(current4Vec); particle4VecMap.insert(std::map<std::string, Vector4<double> >::value_type((*itPart)->name(), current4Vec)); V4_all_lab += current4Vec; counter++; } EvtData* evtData=new EvtData(); evtData->evtNo=evtNo; evtData->evtWeight=theEvent->Weight(); evtData->FourVecsString.insert(mapString4Vec::value_type("all",V4_all_lab)); //cache 4 vectors of inital state particles std::map<std::string, Vector4<double> >::iterator it4VecMap; for (it4VecMap=particle4VecMap.begin(); it4VecMap!=particle4VecMap.end(); ++it4VecMap){ evtData->FourVecsString.insert(mapString4Vec::value_type(it4VecMap->first, it4VecMap->second)); } //fill WignerD functions const std::shared_ptr<AbsChannelEnv> absChannelEnvPtr=GlobalEnv::instance()->Channel(_channelID); std::vector<std::shared_ptr<AbsDecay> > theDecays=absChannelEnvPtr->prodDecayList()->getList(); std::vector<std::shared_ptr<AbsDecay> >::iterator itIso; for (itIso=theDecays.begin(); itIso!=theDecays.end(); ++itIso){ std::string theRefKey=(*itIso)->refKey(); (*itIso)->fillWignerDs(particle4VecMap, V4_all_lab, evtData, theRefKey); } //fill 4Vecs for dynamics std::vector<std::shared_ptr<AbsDynamics> > theDynVec=DynRegistry::instance()->getDynVec(); std::vector<std::shared_ptr<AbsDynamics> >::iterator itDyn; for (itDyn=theDynVec.begin(); itDyn!=theDynVec.end(); ++itDyn){ (*itDyn)->fillMasses(evtData); }; //reset filled map for (itIso=theDecays.begin(); itIso!=theDecays.end(); ++itIso){ (*itIso)->resetFilledMap(); } return evtData; } std::string EvtDataBaseList::getName(std::vector<Particle*>& theVec){ std::string result; std::vector<Particle*>::iterator it; for(it=theVec.begin(); it!=theVec.end(); ++it){ result+=(*it)->name(); } return result; }