//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Evaluation of combined FOPI PID information
//      
//
//
// Environment:
//      Software developed for GEM-TPC detector in FOPI.
//
// Author List:
//      Felix Boehmer    TUM            (original author)
//
//
//-----------------------------------------------------------


#include "FopiPidHub.h"
#include "FopiPidProb.h"
#include "FopiPidProbSet.h"

#include "TList.h"
#include "TKey.h"
#include "TFile.h"
#include "TGraphErrors.h"

#include <cmath>
#include <utility>
#include <iostream>
#include <vector>
#include <cassert>

ClassImp(FopiPidHub)
ClassImp(FopiPidMomWindow)



FopiPidMomWindow::FopiPidMomWindow()
  :
  isRange(false),
  isSharp(false),
  momLow(0.),
  momHigh(0.),
  momSharp(0.)
{;}

FopiPidMomWindow::FopiPidMomWindow(double low, double high)
  :
  isRange(true),
  isSharp(false),
  momLow(low),
  momHigh(high),
  momSharp(0.)
{;}

FopiPidMomWindow::FopiPidMomWindow(double sharp)
  :
  isRange(false),
  isSharp(true),
  momLow(0.),
  momHigh(0.),
  momSharp(sharp)
{;}


FopiPidMomWindow&
FopiPidMomWindow::operator=(const FopiPidMomWindow& other) {
  if(this == &other)
    return *this;
  isRange = other.isRange;
  isSharp = other.isSharp;
  momLow = other.momLow;
  momHigh = other.momHigh;
  momSharp = other.momSharp;
  return *this;
}


bool
FopiPidMomWindow::operator==(const FopiPidMomWindow& other) const {
  if(isRange && other.isRange) {
    double dl = std::fabs(momLow - other.momLow);
    double dh = std::fabs(momHigh - other.momHigh);
    if(dl+dh < 1.e-5) 
      return true;
    else 
      return false;
  }
  else if(isRange && other.isSharp) {
    if((other.momSharp >= momLow) && (other.momSharp <= momHigh))
      return true;
    else
      return false;
  }
  
  else if(isSharp && other.isRange) {
    if((momSharp >= other.momLow) && (momSharp <= other.momHigh))
      return true;
    else
      return false;
  }
  else
    return false;
}


bool
FopiPidMomWindow::operator<(const FopiPidMomWindow& other) const {
  if(isRange && other.isRange) {
    return momLow < other.momLow;
  }
  else if(isRange && other.isSharp) {
    return momHigh < other.momSharp;
  }
  else if(isSharp && other.isRange) {
    return momSharp < other.momLow;
  }
  else
    return momLow < other.momLow;
}


//IMPLEMENTATIONS FOR FOPIPIDHUB ---------------------------------------------------

FopiPidHub::FopiPidHub() 
:
  fVerbose(false),
  fInit(false),
  fFractionsMC("")
{
  buildBookMap();
}

void
FopiPidHub::setDetProbFile(const TString& det, const TString& file) {
   fFileMap[det]=file;
}


int
FopiPidHub::setCalibrationFile(const TString& det, const TString& file,
			       const TString& pikey, const TString& protkey,
			       int runLow, int runHigh) {

  assert(runLow<runHigh);
  //obtain valid det names, not using boolean info
  if(fBookMap.find(det) == fBookMap.end()) {
    std::cout<<"FopiPidHub::setCalibrationFile() INVALID DETECTOR NAME"
	     <<std::endl;
    return 1;
  }
  fCalibKeyMap[det] = std::pair<TString,TString>(pikey,protkey);
  fCalibIOMap[det] = file;
  fCalibRangeMap[det] = std::pair<int,int>(runLow,runHigh);
  return 0;
}

int 
FopiPidHub::setUseMCFractions(const TString& file) {
  fFractionsMC = file;
}



TString
FopiPidHub::getDetProbFile(const TString& det) const {
  std::map<TString, TString>::const_iterator it;
  it = fFileMap.find(det);
  if(it==fFileMap.end())
    return TString("");
  else 
    return it->second;
}



//TODO close memory leaks with TFile pointers.
//

int 
FopiPidHub::init() {
  //init maps from probability input
  std::map<TString, TString>::const_iterator fit;
  for(fit=fFileMap.begin(); fit!=fFileMap.end(); fit++) {
    TString detName = fit->first;
    TString fileName = fit->second;
    
    if(fileName.Length() > 1) {
      TFile* ifile = new TFile(fileName,"read");
      if(ifile->IsZombie()) {
	std::cerr<<"FopiPidHub::init() Could not read file for detector "
		 <<detName<<std::endl;
	delete ifile;
	return 1;
      }
      
      //now get all keys from input file and fill maps:
      TList* keys = ifile->GetListOfKeys();
      for(unsigned int k=0; k<keys->GetEntries(); k++) {
	TKey* ikey = (TKey*) keys->At(k);
	TString cname = ikey->GetClassName();
	//ADD THE SET:
	if(cname.EqualTo("FopiPidProbSet")) {
	  TString name = ikey->GetName();
	  FopiPidProbSet* iset = (FopiPidProbSet*) ifile->Get(name);
	  if(iset == NULL) {
	    std::cerr<<"FopiPidHub::init() Got invalid ProbSet ... abort init"
		     <<std::endl;
	    return 1;
	  }
	  double mLow, mHigh;
	  iset->getMomWindow(mLow,mHigh);
	  FopiPidMomWindow window(mLow,mHigh);
	  
	  std::map<FopiPidMomWindow, FopiPidProbSet*>* thismap = fBookMap[detName];
	  std::pair<std::map<FopiPidMomWindow, FopiPidProbSet*>::iterator,bool> test;
	  test = thismap->insert(std::pair<FopiPidMomWindow, FopiPidProbSet*>(window,iset));
	  //check if insertion was successful:
	  if(test.second==false) {
	     std::cerr<<"FopiPidHub::init() Could not insert ProbSet!"
		      <<std::endl;
	     return 1;
	  }
	}
	
      } //end loop over keys
      
    } //end check for filename set	
	
  } //end loop over files

  //build calibration maps:  
  int calStat = buildCalibrationMaps();
  if(calStat)
    return calStat;
  
  //build MC fraction map:
  if(!fFractionsMC.EqualTo("")) {
    TFile* fracFile = new TFile(fFractionsMC.Data(), "read");
    if(fracFile->IsZombie())
      std::cout<<"FopiPidHub::init() Could not open MC fraction input file"
	       <<std::endl;

    TList* fracKeys = fracFile->GetListOfKeys();
    for(unsigned int k=0; k<fracKeys->GetEntries(); k++) {
      TKey* ikey = (TKey*) fracKeys->At(k);
      TString cname = ikey->GetClassName();
      //ADD THE SET:
      if(cname.EqualTo("FopiPidFractionSet")) {
	TString name = ikey->GetName();
	FopiPidFractionSet* iset = (FopiPidFractionSet*) fracFile->Get(name.Data());
	if(iset == NULL) {
	  std::cerr<<"FopiPidHub::init() Got invalid FractionSet ... abort init"
		   <<std::endl;
	  fracFile->Close();
	  delete fracFile;
	  return 1;
	}
	double mLow, mHigh;
	iset->getMomWindow(mLow,mHigh);
	FopiPidMomWindow window(mLow,mHigh);
	
	std::pair<std::map<FopiPidMomWindow, FopiPidFractionSet>::iterator,bool> test;
	std::pair<FopiPidMomWindow, FopiPidFractionSet> ins(window,FopiPidFractionSet(*iset));
	test = fFracMCMap.insert(ins);
	 //check if insertion was successful:
	  if(test.second==false) {
	     std::cerr<<"FopiPidHub::init() Could not insert FractionSet!"
		      <<std::endl;
	     fracFile->Close();
	     delete fracFile;
	     return 1;
	  }
      }
    }
    fracFile->Close();
    delete fracFile;
  }

  fInit = true;
  return 0;
}


double
FopiPidHub::getFraction(const TString& part, double mom) const {
  if(!fInit) {
    std::cout<<"FopiPidHub::getFraction() I am not initialized! Returning -1."
	       <<std::endl;
    return -1.;
  }
  
  FopiPidMomWindow win(mom);
  
  //if MC fractions are set, use them:
  if(!fFractionsMC.EqualTo("")) {
    std::map<FopiPidMomWindow, FopiPidFractionSet>::const_iterator mcit;
    mcit = fFracMCMap.find(win);
    if(mcit==fFracMCMap.end()) {
      std::cout<<"FopiPidHub::getFraction() Could not find MC fraction for mom "<<mom
	       <<std::endl;
      return -1.;
    }
    double mcfrac = mcit->second.getFraction(part);
    if(mcfrac < 0.) {
      std::cout<<"FopiPidHub::getFraction() Lookup failed for particle "
	       <<part.Data()<<std::endl;
      return -1.;
    }
    return mcfrac;
  }
  
  //otherwise use values form TPC dE/dx (large errors!)
  else {
    std::map<TString, std::map<FopiPidMomWindow,FopiPidProbSet*>* >::const_iterator meh;
    meh = fBookMap.find("tpc");
    std::map<FopiPidMomWindow, FopiPidProbSet*>* tpcMap;
    if(meh == fBookMap.end())
      return -1.;
    else 
      tpcMap = meh->second;
    
    std::map<FopiPidMomWindow, FopiPidProbSet*>::const_iterator it;
    it = tpcMap->find(win);
    if(it==tpcMap->end())
      return -1.;
    else 
      return it->second->getFraction(part);  
  }
}

 





double
FopiPidHub::getFracLikelihood(const TString& part, double mom, 
			      const FopiPidInfo& data,
			      bool weighted) const {

  if(!fInit) {
    std::cout<<"FopiPidHub::getFracLikelihood() I am not initialized! Returning -1."
	       <<std::endl;
    return -1.; 
  }
  
  std::map<TString, std::map<TString, double> > lis;
  //get likelihoods; data is calibrated if calibration data is present
  int stat = this->getLikelihoodMap(mom, data, lis);
  if(stat>0) {
    std::cout<<"FopiPidHub::getFracLikelihood() Failed to get internal "
	     <<"likelihood mapping; Returning -1."<<std::endl;
    return -1.;
  }   
  
  std::map<TString, std::vector<double> > likes;  //partname <- list of det likelihoods
  std::map<TString, double> weights;              //partname <- weight
  
  std::map<TString, std::map<TString, double> >::const_iterator detit; 
  //loop over single detector likelihoods
  for(detit=lis.begin(); detit!=lis.end(); detit++) {
    std::map<TString, double>::const_iterator partit;
    if(fVerbose)
      std::cout<<"FopiPidHub::getFracLikelihood()  Looking up'"
	       <<detit->first.Data()<<"' ..."<<std::endl;
    //loop over particle hypothesis contained in the information
    for(partit=detit->second.begin(); partit!=detit->second.end(); partit++) {
      
      //make entry in the weightmap if non-existent
      if(!weights.count(partit->first)) {
	if(weighted) {
	  if(fVerbose)
	    std::cout<<"FopiPidHub::getFracLikelihood()  Making entry in weightmap for '"
		     <<partit->first<<"': "
		     << this->getFraction(partit->first,mom)<<std::endl;
	  weights[partit->first] = this->getFraction(partit->first,mom);
	}
	else 
	  weights[partit->first] = 1.;
      }
      likes[partit->first].push_back(partit->second);
      if(fVerbose) 
	std::cout<<"FopiPidHub::getFracLikelihood() Setting lh for '"
		 <<partit->first<<"': "<<partit->second<<std::endl;
      
    }//end loop over particles hypotheses

  } //end loop over detector likelihoods
  
  double yea = 0.;
  double nay = 0.;

  std::map<TString, double >::const_iterator it;
  for(it=weights.begin(); it!=weights.end(); it++) {
    double weight = it->second; //for clarity
    double partlike = weight;
    std::vector<double> partlikes = likes[it->first];
    for(unsigned int il=0; il<partlikes.size(); il++)
      partlike *= partlikes[il];
    if(!it->first.EqualTo(part))
      nay += partlike;
    else
      yea += partlike;
  }
  
  if(fVerbose) 
    std::cout<<"FopiPidHub::getFracLikelihood() Yea: "<<yea<<", Nay: "
	     <<nay<<std::endl;

  if(yea+nay < 1.e-10)
    return 0.;
  return yea / (yea + nay);
}




//get fractional likelihood of hypothesis part versus another hypothesis
double 
FopiPidHub::getFracLikelihood(const TString& part, 
			      const TString& contestant,
			      double mom, 
			      const FopiPidInfo& data,
			      bool weighted) const {

  return this->getFracLikelihood(part, contestant, TString("none"),
				 mom, data, weighted);
  
}



//get fractional likelihood of hypothesis part versus hypotheses contestant 1 & 2
double 
FopiPidHub::getFracLikelihood(const TString& part, 
			      const TString& contestant1,
			      const TString& contestant2,
			      double mom, 
			      const FopiPidInfo& data,
			      bool weighted) const {

  
  if(!fInit) {
    std::cout<<"FopiPidHub::getFracLikelihood() I am not initialized! "
	     <<"Returning -1."
	     <<std::endl;
    return -1.; 
  }
  
  std::map<TString, std::map<TString, double> > lis;
  //get likelihoods; data is calibrated if calibration data is present
  int stat = this->getLikelihoodMap(mom, data, lis);
  if(stat>0) {
    std::cout<<"FopiPidHub::getFracLikelihood() Failed to get internal "
	     <<"likelihood mapping; Returning -1."<<std::endl;
    return -1.;
  }   
  
  std::vector<TString> partlist;
  partlist.push_back(part);
  if(!contestant1.EqualTo("none"))
    partlist.push_back(contestant1);
  if(!contestant2.EqualTo("none"))  
    partlist.push_back(contestant2);
  
  
  std::map<TString, std::vector<double> > likes;  //partname <- list of det likelihoods
  std::map<TString, double> weights;              //partname <- weight
  
  std::map<TString, std::map<TString, double> >::const_iterator detit; 
  //loop over single detector likelihoods
  for(detit=lis.begin(); detit!=lis.end(); detit++) {
    std::map<TString, double>::const_iterator partit;
    if(fVerbose)
      std::cout<<"FopiPidHub::getFracLikelihood()  Looking up'"
	       <<detit->first.Data()<<"' ..."<<std::endl;
    //loop over particle hypothesis contained in the information
    for(unsigned int ipart=0; ipart<partlist.size(); ipart++ ) {
      TString thispart = partlist[ipart];
      partit = detit->second.find(thispart);
      if(partit == detit->second.end()) 
	continue;
      
      //make entry in the weightmap if non-existent
      if(!weights.count(thispart)) {
	if(weighted) {
	  if(fVerbose)
	    std::cout<<"FopiPidHub::getFracLikelihood()  Making entry in weightmap for '"
		     <<thispart<<"': "
		     << this->getFraction(thispart,mom)<<std::endl;
	  weights[thispart] = this->getFraction(thispart,mom);
	}
	else 
	  weights[thispart] = 1.;
      }
      likes[thispart].push_back(partit->second);
      if(fVerbose) 
	std::cout<<"FopiPidHub::getFracLikelihood() Setting lh for '"
		 <<thispart<<"': "<<partit->second<<std::endl;
      
    }//end loop over particles hypotheses

  } //end loop over detector likelihoods
  
  double yea = 0.;
  double nay = 0.;

  
  for(unsigned int ipart=0; ipart<partlist.size(); ipart++ ) {
    double weight = weights[(partlist[ipart])];
    if(weight < 0.) //protection in case no MC weights are used
      continue;
    std::vector<double> likelihoods = likes[(partlist[ipart])];
    double likelihood = weight;
    for(unsigned int il=0; il<likelihoods.size(); il++) 
      likelihood *= likelihoods[il];
    
    if(partlist[ipart].EqualTo(part)) {
      if(fVerbose) 
 	std::cout<<"      Adding "<<likelihood<<" to 'yea' for particle "
		 <<partlist[ipart].Data()<<std::endl;	  
      yea += likelihood;
    }
    else {
      if(fVerbose) 
	std::cout<<"      Adding "<<likelihood<<" to 'nay' for particle "
		 <<partlist[ipart].Data()<<std::endl;	  
      nay += likelihood;
    }
  }
  
  // std::map<TString, double >::const_iterator it;
  // for(it=weights.begin(); it!=weights.end(); it++) {
  //   double weight = it->second; //for clarity
  //   double partlike = weight;
  //   std::vector<double> partlikes = likes[it->first];
  //   for(unsigned int il=0; il<partlikes.size(); il++)
  //     partlike *= partlikes[il];
  //   if(!it->first.EqualTo(part))
  //     nay += partlike;
  //   else
  //     yea += partlike;
  // }
  
  if(fVerbose) 
    std::cout<<"FopiPidHub::getFracLikelihood() Yea: "<<yea<<", Nay: "
	     <<nay<<std::endl;

  if(yea+nay < 1.e-10)
    return 0.;
  return yea / (yea + nay);
  
}


//return a likelihood data for a given momentum and PID information
//format:  <detname, <partname, likelihood> >
//contains only valid entries (detectors, particles)
//for TPC and CDC kaon likelihoods are treated equal to the pion/electron peak 
//if no kaon fit is available. However, when using TPC dEdx fractions, the weighted
//kaon likelihhods will be zero.
//Validity ranges of each FopiPidProbSet are checked. If data lies outside, 
//the returned likelihoods will be zero.
int
FopiPidHub::getLikelihoodMap(double mom, const FopiPidInfo& uc_data,
			     std::map<TString, std::map<TString, double> >& mapret) const{
  
  mapret.clear();

  //CALIBRATE DATA
  FopiPidInfo data = getCalibratedData(uc_data);
    
  if(fVerbose)
    std::cout<<"FopiPidHub::getLikeliHoodMap()"<<std::endl;
  
  //reformatting of input data 
  std::vector<TString> dets;
  if(data.hasTpc())
    dets.push_back("tpc");
  if(data.hasCdc())
    dets.push_back("cdc");
  if(data.hasBar())
    dets.push_back("bar");
  if(data.hasRpc())
    dets.push_back("rpc");

  //not using the boolean information:
  std::map<TString,bool>* partmap = makePartMap("pion");
  if(partmap==NULL) {
    std::cerr<<"FopiPidHub::getLikelihoodMap() INCONSISTENT PARTICLE NAME"
	     <<std::endl;
    return -1.;
  }

  std::map<TString, std::map<FopiPidMomWindow,FopiPidProbSet*>* >::const_iterator det;
  //loop over detectors
  for(unsigned int id=0; id<dets.size(); id++) {
    TString detname = dets[id];
    if(fVerbose) 
      std::cout<<"  getting likelihoods from detector: "<<detname.Data()<<std::endl;

    det = fBookMap.find(detname);
    if(det==fBookMap.end()) {
      std::cerr<<"FopiPidHub::getFracLikelihood() INCONSISTENT DETECTOR NAME"
	       <<std::endl;
      delete partmap;
      return -1.;
    }
    
    
    std::map<FopiPidMomWindow, FopiPidProbSet*>* detMap = det->second;
    //build accessor:
    FopiPidMomWindow win(mom);
    std::map<FopiPidMomWindow, FopiPidProbSet*>::const_iterator setit;
    setit = detMap->find(win);

    std::map<TString, double> thisli;
    //in case we found info at the given momentum:
    if(setit != detMap->end()) {
      FopiPidProbSet* thisset = setit->second;
      double rangelow, rangehigh;
      thisset->getValidityRange(rangelow, rangehigh);
      
      std::map<TString,bool>::const_iterator partit;
      
      //loop over hypotheses:
      for(partit=partmap->begin(); partit!=partmap->end(); partit++) {
	if(fVerbose)
	  std::cout<<"    particle: "<<partit->first.Data()<<std::endl;
	const FopiPidProb* thispart = thisset->getProb(partit->first);
	
	if(thispart==NULL) { 
	  //for TPC, CDC:
	  //if kaon fit not present, assume it has been merged with pion peak and 
	  //return that info
	  if(detname.EqualTo("tpc") || detname.EqualTo("cdc") ) {
	    if(partit->first.EqualTo("kaon")) {
	      if(fVerbose)
		std::cout<<"    fetching pion likelihood instead kaon"
			 <<" for mom "<<mom<<std::endl;
	      thispart = thisset->getProb(TString("pion"));
	    }
	    
	  }
	}
	if(thispart==NULL)
	  continue;
	
	//incoming information for this detector:
	double pidval = data.getVal(detname);
	
	double lh = thispart->getLikelihoodValue(pidval,pidval+0.03);
	if(pidval < rangelow || pidval > rangehigh)
	  lh = 0.;
	
	if(fVerbose) 
	  std::cout<<"      likelihood = "<<lh<<" at "<<pidval<<std::endl;
	
	//store likelihood
	thisli[partit->first] = lh;
	
      }//end loop over hypotheses
      
    }
    else if(fVerbose) 
      std::cout<<" could not get Momentum Window for mom = "<<mom<<std::endl;
    mapret[detname] = thisli;
   
    
  }//end loop over detectors
  delete partmap;
  return 0;
}


// TEST IMPLEMENTATIONS BELOW --------------------------------------------



bool
FopiPidHub::addTestWindow(const FopiPidMomWindow& w, double val) {
  std::pair<std::map<FopiPidMomWindow,double>::iterator,bool> ins;
  ins = fTestMap.insert(std::pair<FopiPidMomWindow,double>(w,val));
  return ins.second;
}

double
FopiPidHub::getDouble(const FopiPidMomWindow& win) const {
  std::map<FopiPidMomWindow, double>::const_iterator it;
  it = fTestMap.find(win);
  if(it!=fTestMap.end())
    return it->second;
  else 
    return NULL;
}


 
void
FopiPidHub::buildBookMap() {
  fBookMap["tpc"] = &fTpcMap;
  fBookMap["cdc"] = &fCdcMap;
  fBookMap["bar"] = &fBarMap;
  fBookMap["rpc"] = &fRpcMap;
  //fFileMap["tpc"] = "";
  //fFileMap["cdc"] = "";
  //fFileMap["bar"] = "";
  //fFileMap["rpc"] = "";
  fCalibValMap["tpc"] = &fTpcCalibMap;
  fCalibValMap["cdc"] = &fCdcCalibMap;
  fCalibValMap["bar"] = &fBarCalibMap;
  fCalibValMap["rpc"] = &fRpcCalibMap;
}


//if created successfully, caller takes ownership
std::map<TString,bool>*
FopiPidHub::makePartMap(const TString& part) const {
  std::map<TString,bool>* map = new std::map<TString,bool>();
  (*map)["pion"] = false;
  (*map)["kaon"] = false;
  (*map)["proton"] = false;
  (*map)["deuteron"] = false;
  
  if(map->count(part)) {
    (*map)[part] = true;
    return map;
  }
  else {
    delete map; 
    return NULL;
  }
}



double
FopiPidHub::getSumFractions(double mom) const {
  //not making use of the bools
  std::map<TString,bool>* parts = makePartMap("pion");

  double sumfrac = 0.;
  std::map<TString,bool>::const_iterator it;
  for(it=parts->begin(); it!=parts->end(); it++) {
    double frac = getFraction(it->first,mom);
    std::cout<<"Momentum: "<<mom<<"  particle: "<<it->first.Data()
	     <<"  fraction: "<<frac<<std::endl;
    if(frac>0.)
      sumfrac+=frac;
  }
  std::cout<<"    TOTAL: "<<sumfrac<<std::endl;
  return sumfrac;
}
  

int
FopiPidHub::buildCalibrationMaps() {
  
  std::map<TString, TString >::const_iterator ioit;
  if(fVerbose) 
    std::cout<<"FopiPidHub::buildCalibrationMaps() ... "<<std::endl;
  //loop over detectors
  for(ioit=fCalibIOMap.begin(); ioit!=fCalibIOMap.end(); ioit++) {
    if(fVerbose) 
      std::cout<<" ... loading calibration data for '"<<ioit->first.Data()
	       <<"' ..."<<std::endl;
    
    TFile* cfile = new TFile(ioit->second,"read");
    if(cfile->IsZombie()) {
      std::cerr<<"    could not open file "<<ioit->second.Data()
	       <<"; aborting"<<std::endl;
      return 1;
    }
    
    std::map<TString, std::pair<TString, TString> >::const_iterator keyit;
    keyit = fCalibKeyMap.find(ioit->first);
    std::pair<TString, TString> keypair = keyit->second;
    TString pionkey = keypair.first;
    TString protkey = keypair.second;
    //get data for this det
    TGraphErrors* calPion = (TGraphErrors*) cfile->Get(pionkey);   
    TGraphErrors* calProt = (TGraphErrors*) cfile->Get(protkey);   
    if(calPion == NULL || calProt == NULL ) {
      std::cerr<<"    could not get calibration data" 
	       <<" from file; aborting"<<std::endl;
      return 1;
    }
    
    //identify stable zone:
    int runLow = (fCalibRangeMap[ioit->first]).first;
    int runHigh = (fCalibRangeMap[ioit->first]).second;

    std::cout<<" ... building calibration lookup map for '"<<ioit->first
	     <<"' and stable range ["<<runLow<<", "<<runHigh
	     <<"]"<<std::endl;
    
    unsigned int nPi = calPion->GetN();
    unsigned int nPr = calProt->GetN();
    double x, y, erry, x2, y2;
    double mean = 0.;
    double totErr = 0.;
    for(unsigned int n=0; n<nPi; n++) {
      calPion->GetPoint(n,x,y);
      if(x < runLow)
	continue;
      if(x > runHigh) 
	break;
      erry = calPion->GetErrorY(n);
      mean += (1./erry) * y;
      totErr += (1./erry);
    }
    double meanPi = mean / totErr;

    mean = 0.;
    totErr = 0.;
    for(unsigned int n=0; n<nPr; n++) {
      calProt->GetPoint(n,x,y);
      if(x < runLow)
	continue;
      if(x > runHigh) 
	break;
      erry = calProt->GetErrorY(n);
      mean += (1./erry) * y;
      totErr += (1./erry);
    }
    double meanPr = mean / totErr;
    
    std::cout<<"     found mean values:  pion: "<<meanPi<<", "
	     <<" proton: "<<meanPr<<std::endl;
    
    //now build lookup map:
    std::map<TString, std::map<int, std::pair<double,double> >* >::const_iterator it;
    it = fCalibValMap.find(ioit->first);
    if(it == fCalibValMap.end()) {
      std::cout<<"   something went wrong ... "<<std::endl;
      return 1;
    }
    std::map<int, std::pair<double,double> >* thisCalMap = it->second;
        
    double val1, val2; 
    for(unsigned int n=0; n<nPi; n++) { 
      calPion->GetPoint(n,x,val1);
      calProt->GetPoint(n,x2,val2);
      //check for run number consistency
      if(std::fabs(x-x2) > 1.e-1) {
	std::cout<<" INCONSISTENT DATA FOR PI AND PROT.... aborting"
		 <<std::endl;
	return 1;
      }
      double v1 = val1-meanPi;
      double v2 = val2-meanPr;
      
      //calculate linear calibration parameters
      double m = (v2-v1)/(meanPr-meanPi);
      double t = v1 - m*meanPi;

      double drun = x+0.1;
      int run = (int) drun;
      // if(run == 3445)
      // 	std::cout<<"RUN 3445: v1="<<v1<<", v2="<<v2
      // 		 <<", m="<<m<<", t="<<t<<std::endl;
      (*thisCalMap)[run] = std::pair<double,double>(m,t);
            
    }      

  } //end loop over detectors
   
  std::cout<<"FopiPidHub::buildCalibrationMaps() success."<<std::endl;
  return 0;
}



FopiPidInfo
FopiPidHub::getCalibratedData(const FopiPidInfo& uc_data) const {
  FopiPidInfo ret(uc_data);
  
  std::map<TString, std::map<int, std::pair<double,double> >* >::const_iterator detit;
  //loop over detectors
  for(detit=fCalibValMap.begin(); detit!=fCalibValMap.end(); detit++) {
    if(!uc_data.hasVal(detit->first))
      continue; //nothing to correct

    int run = uc_data.getRunNr();

    std::map<int, std::pair<double,double> >* thisMap;
    thisMap = detit->second;
    if(!thisMap->count(run)) //empty calib data
      continue;
    
    //check if inside calibration range
    std::pair<int,int> calRange;
    calRange = fCalibRangeMap.find(detit->first)->second;
    if(run >= calRange.first && run <= calRange.second)
      continue;  //nothing to do        
    
    std::pair<double,double> calibData;
    calibData = thisMap->find(run)->second;
    
    double val = uc_data.getVal(detit->first);
    double CO = uc_data.getCO(detit->first);
    

    //if we have our values in a different space than the raw PID data spectra used for computing the
    //calibration parameters, we need to transform into that raw system first, apply the calibration and
    //transform back. Usecase: realtive dE/dx spectra
    double valNorm = val - CO;
    

    double m = calibData.first;
    double t = calibData.second;
    
    double newVal = valNorm - m*valNorm - t;
    ret.setVal(detit->first, newVal+CO,CO);
    
  } //end loop over detectors
  
  return ret;

}