Skip to content
Snippets Groups Projects
Commit 5fd5b723 authored by Bertram Kopf's avatar Bertram Kopf
Browse files

modified fitServerMode method to get pipi scattering data running in coupled channel mode

parent abec6e18
No related branches found
No related tags found
No related merge requests found
......@@ -149,14 +149,14 @@ void AppBase::readEvents(EventList& theEventList, std::vector<std::string>& file
theEventList.rewind();
}
void AppBase::readScatteringEvents(EventList& theEventList, std::vector<std::string>& fileNames, ChannelID channelID){
void AppBase::readScatteringEvents(EventList& theEventList, std::vector<std::string>& fileNames, ChannelID channelID, int evtStart, int evtStop){
std::vector< std::shared_ptr<MassRangeCut> > massRangeCuts=GlobalEnv::instance()->Channel(channelID)->massRangeCuts();
EventReaderScattering evtScatterReader(fileNames, 2, 0, false);
if(GlobalEnv::instance()->Channel(channelID)->useMassRange()){
evtScatterReader.setMassRange(massRangeCuts);
}
evtScatterReader.fill(theEventList);
evtScatterReader.fill(theEventList, evtStart, evtStop);
}
void AppBase::createLhObjects(){
......@@ -603,6 +603,7 @@ bool AppBase::calcAndSendClientLh(NetworkClient& theClient, std::shared_ptr<AbsP
}
void AppBase::fitServerMode(std::shared_ptr<AbsPawianParameters> upar){
double evtWeightSumData=0;
ChannelEnvList channelEnvs=GlobalEnv::instance()->ChannelEnvs();
std::map<short, std::tuple<long, double, long> > numEventMap;
......@@ -621,35 +622,46 @@ void AppBase::fitServerMode(std::shared_ptr<AbsPawianParameters> upar){
std::vector<std::string> mcFileNames;
mcFileNames.push_back(mcFile);
double noOfMcEvts=0.;
EventList eventsData;
readEvents(eventsData, dataFileNames, (*it).first->channelID(), (*it).first->useDataEvtWeight(), 0, noOfDataEvents);
EventList mcData;
if(GlobalEnv::instance()->Channel((*it).first->channelID())->channelType() == AbsChannelEnv::CHANNEL_PIPISCATTERING){
readScatteringEvents(eventsData, dataFileNames, (*it).first->channelID(), 0, noOfDataEvents);
}
else{
readEvents(eventsData, dataFileNames, (*it).first->channelID(), (*it).first->useDataEvtWeight(), 0, noOfDataEvents);
}
double noOfWeightedDataEvts=EvtDataBaseList::noOfWeightedEvts(eventsData, (*it).first->channelID(), noOfDataEvents, 0);
int noDataEvts=eventsData.size();
double noOfWeightedDataEvts=EvtDataBaseList::noOfWeightedEvts(eventsData, (*it).first->channelID(), noOfDataEvents, 0);
evtWeightSumData+=noOfWeightedDataEvts;
eventsData.removeAndDeleteEvents(0, eventsData.size()-1);
int maxMcEvts=noDataEvts*ratioMcToData;
EventList mcData;
readEvents(mcData, mcFileNames, (*it).first->channelID(), (*it).first->useMCEvtWeight(), 0, maxMcEvts-1);
double noOfMcEvts=EvtDataBaseList::noOfWeightedEvts(mcData, (*it).first->channelID(), maxMcEvts-1, 0);
if(GlobalEnv::instance()->Channel((*it).first->channelID())->channelType() != AbsChannelEnv::CHANNEL_PIPISCATTERING){
int maxMcEvts=noDataEvts*ratioMcToData;
readEvents(mcData, mcFileNames, (*it).first->channelID(), (*it).first->useMCEvtWeight(), 0, maxMcEvts-1);
noOfMcEvts=EvtDataBaseList::noOfWeightedEvts(mcData, (*it).first->channelID(), maxMcEvts-1, 0);
}
evtWeightSumData+=noOfWeightedDataEvts;
numEventMap[(*it).first->channelID()] = std::tuple<long, double,long>(noDataEvts, noOfWeightedDataEvts, mcData.size());
mcData.removeAndDeleteEvents(0, mcData.size()-1);
if(noOfWeightedDataEvts<10.){
Alert << "number of wieghted data events too small: " << noOfWeightedDataEvts << endmsg;
exit(1);
}
if(noOfMcEvts<10.){
Alert << "number of wieghted Monte Carlo events too small: " << noOfMcEvts << endmsg;
Alert << "number of weighted data events too small: " << noOfWeightedDataEvts << endmsg;
exit(1);
}
if(GlobalEnv::instance()->Channel((*it).first->channelID())->channelType() != AbsChannelEnv::CHANNEL_PIPISCATTERING && noOfMcEvts<10.){
if(noOfMcEvts<10.){
Alert << "number of weighted Monte Carlo events too small: " << noOfMcEvts << endmsg;
exit(1);
}
}
}
std::shared_ptr<AbsFcn> absFcn;
std::shared_ptr<NetworkServer> theServer(new NetworkServer(GlobalEnv::instance()->parser()->serverPort(), GlobalEnv::instance()->parser()->noOfClients(), numEventMap, GlobalEnv::instance()->parser()->clientNumberWeights()));
// theServer->WaitForFirstClientLogin();
......
......@@ -56,7 +56,7 @@ public:
virtual void dumpDefaultParams();
virtual void generate(std::shared_ptr<AbsPawianParameters> theParams);
virtual void readEvents(EventList& theEventList, std::vector<std::string>& fileNames, ChannelID channelID, bool withEvtWeight=false, int evtStart=0, int evtStop=1000000);
virtual void readScatteringEvents(EventList& theEventList, std::vector<std::string>& fileNames, ChannelID channelID);
virtual void readScatteringEvents(EventList& theEventList, std::vector<std::string>& fileNames, ChannelID channelID, int evtStart=0, int evtStop=1000000);
virtual void createLhObjects();
virtual void qaMode(std::shared_ptr<AbsPawianParameters> startParams, double evtWeightSumData);
virtual void qaModeSimple(EventList& dataEventList, EventList& mcEventList, std::shared_ptr<AbsPawianParameters> startParams);
......
......@@ -73,7 +73,8 @@ int main(int __argc,char *__argv[]){
std::string currentArgv(__argv[i]);
if(currentArgv !=(char*)"--pbarpFiles"
&& currentArgv !=(char*)"--epemFiles"
&& currentArgv !=(char*)"--resFiles"
&& currentArgv !=(char*)"--resFiles"
&& currentArgv !=(char*)"--pipiScatteringFiles"
&& currentArgv !="-c"
&& currentArgv !="--configFile"){
argvWoCfgFile[argcWoCfgFile]=__argv[i];
......
......@@ -167,33 +167,36 @@ double EvtDataBaseList::noOfWeightedEvts(EventList& evtList, ChannelID channelID
int evtCount = 0;
while ((anEvent = evtList.nextEvent())){
if (evtCount>= maxEvts) break;
if (evtCount%500 == 0) InfoMsg << "4vec calculation for event " << evtCount ; // << endmsg;
Vector4<double> V4_all_lab(0.,0.,0.,0.);
std::vector<Particle*> finalStateParticles=GlobalEnv::instance()->Channel(channelID)->finalStateParticles();
std::vector<Particle*>::iterator itPart;
int counter=0;
for (itPart=finalStateParticles.begin(); itPart != finalStateParticles.end(); ++itPart){
Vector4<float> current4VecFloat=*(anEvent->p4(counter));
Vector4<double> current4Vec(current4VecFloat.E(), current4VecFloat.Px(), current4VecFloat.Py(), current4VecFloat.Pz());
V4_all_lab += current4Vec;
counter++;
}
if (evtCount%10000 == 0){
InfoMsg << "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() <<"\n"
<< "weight: " << anEvent->Weight(); // << endmsg;
}
result += anEvent->Weight();
++evtCount;
if(GlobalEnv::instance()->Channel(channelID)->channelType() != AbsChannelEnv::CHANNEL_PIPISCATTERING){
if (evtCount%500 == 0) InfoMsg << "4vec calculation for event " << evtCount ; // << endmsg;
Vector4<double> V4_all_lab(0.,0.,0.,0.);
std::vector<Particle*> finalStateParticles=GlobalEnv::instance()->Channel(channelID)->finalStateParticles();
std::vector<Particle*>::iterator itPart;
int counter=0;
for (itPart=finalStateParticles.begin(); itPart != finalStateParticles.end(); ++itPart){
Vector4<float> current4VecFloat=*(anEvent->p4(counter));
Vector4<double> current4Vec(current4VecFloat.E(), current4VecFloat.Px(), current4VecFloat.Py(), current4VecFloat.Pz());
V4_all_lab += current4Vec;
counter++;
}
if (evtCount%10000 == 0){
InfoMsg << "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() <<"\n"
<< "weight: " << anEvent->Weight(); // << endmsg;
}
}
}
return result;
}
......
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