/////////////////////////////////////////////////////////////////////
// PixModuleGroup.cxx
// version 1.0.2
/////////////////////////////////////////////////////////////////////
//
// 08/04/03  Version 1.0 (PM)
//           Initial release
// 14/04/03  Version 1.0.1 (CS)
//           Added Configuration DataBase
// 16/04/03  Version 1.0.2 (PM)
//           TPLL interface
//


#include "PixModuleGroup/PixModuleGroup.h"
#include "PixConfDBInterface/PixConfDBInterface.h"

#include "PixController/RodPixController.h"
#include "PixController/TpllPixController.h"
#include "PixController/USBPixController.h"

#include "PixController/PixScan.h"
#include "PixModule/PixModule.h"
#include "Config/Config.h"
#include "PixBoc/PixBoc.h"
#include "PixFe/PixFe.h"
#include "PixMcc/PixMcc.h"
#include "Histo/Histo.h"
#include "Config/Config.h"
#include "PixFe/PixFeI4A.h"
#include "PixFe/PixFeI4B.h"
#include "PixDcs/PixDcs.h"
#include "PixDcs/SleepWrapped.h"
#include <FitClass.h>  

#ifndef WIN32
#include <unistd.h>
#endif
#include <stdlib.h>
#include "cmath"

#define PMG_DEBUG false

using namespace SctPixelRod;
using namespace PixLib;

PixModuleGroup::PixModuleGroup(PixConfDBInterface *dbx, DBInquire *dbi, SctPixelRod::VmeInterface &vme) :
  m_db(dbx), m_vmeIf(vme), m_name(""), m_dbInquire(dbi) {

  // Read database fields
  readDbInquire();
   // Look for Pixel Controller
  fieldIterator f;
  recordIterator r=m_dbInquire->recordBegin();
  for( ; (r!=m_dbInquire->recordEnd()) && ((*r)->getName()!="PixController"); r++) {}
  if(r!=m_dbInquire->recordEnd()) {

    // Create Pixel Controller
    f = (*r)->findField("ActualClassName");

    m_db->DBProcess(f,READ, m_rodName);
    m_pixCtrl = PixController::make(*this, *r, m_rodName);

    // Create Boc if the controller is a ROD
    m_pixBoc = 0;
    RodPixController *rod = dynamic_cast<RodPixController *>(m_pixCtrl);
    if (rod) {
      cout << "searching for BOC in DB" << endl;
      for(recordIterator b = m_dbInquire->recordBegin(); b!=m_dbInquire->recordEnd(); b++){
	if ((*b)->getName() == "PixBoc"){
	  std::cout << "creating BOC" << std::endl;
	  m_pixBoc = new PixBoc(*this, *b);
	}
      }
    }

    // Create modules
    for(recordIterator it = m_dbInquire->recordBegin(); it != m_dbInquire->recordEnd(); it++){
      if((*it)->getName() == "PixModule"){
	PixModule *tempModule = new PixModule(*it, this, (*it)->getDecName().c_str());
	m_modules.push_back(tempModule);
	m_dcsChans.push_back(0);
	tempModule->loadConfig("PHYSICS");
      }
    }
  }
}

PixModuleGroup::PixModuleGroup(PixConfDBInterface *dbx, SctPixelRod::VmeInterface &vme) :
  m_db(dbx), m_vmeIf(vme), m_name(""), m_dbInquire(NULL) {

  // usare conf default  readDbInquire();

  m_pixCtrl = PixController::make(*this, "RodPixController");

}

PixModuleGroup::~PixModuleGroup() {
  // Delete Pixel Controller
  if (m_pixCtrl != NULL) delete m_pixCtrl;

  // Delete modules
  moduleIterator m, mEnd=m_modules.end();
  for(m=m_modules.begin(); m!=mEnd; m++) delete *m;
  m_modules.clear();
  m_dcsChans.clear();
}

void PixModuleGroup::initHW() {
  // Initialize the PixController
  m_pixCtrl->initHW();
  // Initialize the Boc
  if (m_pixBoc != NULL) m_pixBoc->BocConfigure();
}

void PixModuleGroup::testHW() {
}

void PixModuleGroup::init() {
  initHW();
  downloadConfig();
  configure();
}

void PixModuleGroup::readConfig(DBInquire *dbi) {  //! load configuration from the database
  m_dbInquire = dbi;
  // Decode group name
  fieldIterator f = m_dbInquire->findField("ModuleGroupName");
  m_db->DBProcess(f,READ,m_rodName);
  readDbInquire();
  // Load module config
  moduleIterator m, mEnd=m_modules.end();
  for(m=m_modules.begin(); m!=mEnd; m++) {
    (*m)->loadConfig("");
  }
}

void PixModuleGroup::writeConfig(DBInquire *dbi){ // save current configuration into the database
}

void PixModuleGroup::downloadConfig() { // write current configuration into the (possibly present) PixController
  moduleIterator m, mEnd=m_modules.end();
  for(m=m_modules.begin(); m!=mEnd; m++) {
    (*m)->writeConfig();
  }
}

void PixModuleGroup::configure(){ // set current configuration into the actual phisical modules
}                                 //JW: nützlich!

void PixModuleGroup::storeConfig(std::string cfgName) { // give a name to the current config
}

void PixModuleGroup::recallConfig(std::string cfgName){ // recall a named config
}

PixModule* PixModuleGroup::module(int im) {
  for (unsigned int i=0; i<m_modules.size(); i++) {
    if (m_modules[i]->moduleId() == im) return m_modules[i];
  }
  return NULL;
}
void PixModuleGroup::setDcsChan(int modID, PixDcsChan *dcsChan){
  for (unsigned int i=0; i<m_modules.size(); i++) {
    if (m_modules[i]->moduleId() == modID) m_dcsChans[i] = dcsChan;
  }
}

void PixModuleGroup::readDbInquire(){
  // Create the config object
  m_config = new Config("PixModuleGroup");
  Config &conf = *m_config;

  // Group Ggeneral
  conf.addGroup("general");

  conf["general"].addString("ModuleGroupName", m_name, "DEF",
                         "Module Group Name", true);
  conf["general"].addInt("TriggerDelay", m_triggerDelay, 0,
                         "Global Trigger Delay for in-time threshold scan", true);

  // Read from DB
  conf.reset();
  conf.read(m_dbInquire);

  if (m_name == "DEF") {
    fieldIterator f = m_dbInquire->findField("ModuleGroupName");
    m_db->DBProcess(f,READ,m_name);
  }
}

void PixModuleGroup::prepareTDACTuning(int nloop, PixScan *scn) {
	// Module loop
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
		if (m_modules[pmod]->m_readoutActive) {
			unsigned int mod = m_modules[pmod]->m_moduleId;
			int trg = scn->getThresholdTargetValue();
			Histo *hm, *hThr, *hTdac;
			if (scn->scanIndex(nloop) == 0) {
				hThr = new Histo("Threshold", "Threshold", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
				hTdac = new Histo("TDAC", "Tdac Tuning", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
				// Associate different parameters if the tcad tuning is done in loop 1 or 2
				if (nloop == 1) {
					scn->addHisto(*hThr, PixScan::TDAC_THR, mod, scn->scanIndex(2), -1, -1);
					scn->addHisto(*hTdac, PixScan::TDAC_T, mod, scn->scanIndex(2), -1, -1);
				} else if (nloop == 2) {
					scn->addHisto(*hThr, PixScan::TDAC_THR, mod, -1, -1, -1);
					scn->addHisto(*hTdac, PixScan::TDAC_T, mod, -1, -1, -1);
				}
			} else {
				if (nloop == 1) {
					hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1)-1, 0);
					hThr = &scn->getHisto(PixScan::TDAC_THR, mod, scn->scanIndex(2), 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, scn->scanIndex(2), 0, 0);
				} else if (nloop == 2) {
					hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(nloop)-1, 0, 0);
					hThr = &scn->getHisto(PixScan::TDAC_THR, mod, 0, 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, 0, 0, 0);
				}
			}
			// FE loop
			for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
				// check both FE-I4 flavours - fei4-pointer is only used to check if it's NULL
				PixFe* fei4 = dynamic_cast<PixFeI4A*>(*fe);
				if(fei4==0) fei4 = dynamic_cast<PixFeI4B*>(*fe);
				int maxtrim;
				(*fe)->getTrimMax("TDAC", maxtrim);
				// Bisection method
				ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");
				for (unsigned int col=0; col<tdacs.get().size(); col++) {
					for (unsigned int row=0; row<tdacs[0].size(); row++) {
						unsigned int colmod, rowmod;
						int ife = (*fe)->number();
						// 	    if (ife < 8) {
						rowmod = row;
						colmod = col + NCOLS*ife;
						// 	    } else {
						// 	      rowmod = 319 - row;
						// 	      colmod = 143 - col - 18*(ife-8);
						// 	    }
						int newTDAC;
						if (scn->scanIndex(nloop) == 0) {
							newTDAC = (int)(scn->getLoopVarValues(nloop))[0];
						} else {
							int corr = 0;
							// Correct bisection for pixels in the tails
							if (scn->scanIndex(nloop) == 2 && fei4==0){ // correction makes no sense for FE-I4
								int th0 = (int)((*hThr)(colmod, rowmod));
								int th1 = (int)((*hm)(colmod, rowmod));
								int td0 = (int)((scn->getLoopVarValues(nloop))[0]);
								int td1;
								if (th0 > trg) {
									td1 = td0 -(int)((scn->getLoopVarValues(nloop))[1]);
								} else {
									td1 = td0 + (int)((scn->getLoopVarValues(nloop))[1]);
								}
								int delta = (int)((scn->getLoopVarValues(nloop))[2]);
								if (th1 != th0) {
									int ttt = td1 + (trg-th1)*(td1-td0)/(th1-th0);
									if (ttt < 25) {
										ttt = 35;
									} else if (ttt > 117) {
										ttt = 107;
									} else if (ttt < 35) {
										ttt += 10;
									} else if (ttt > 100) {
										ttt -= 6;
									}
									if (td1 < td0 && ttt < td1) {
										corr = ttt - td1 + delta;
									} else if (td1 > td0 && ttt > td1) {
										corr = ttt - td1 - delta;
									}
								}
							}
							// Update best-mach histo
							if (scn->scanIndex(nloop) == 1) {
								hThr->set(colmod, rowmod, (*hm)(colmod,rowmod));
								hTdac->set(colmod, rowmod, tdacs[col][row]);
							} else {
								if (abs((*hm)(colmod,rowmod)-trg) < abs((*hThr)(colmod,rowmod)-trg)) {
									hThr->set(colmod, rowmod, (*hm)(colmod,rowmod));
									hTdac->set(colmod, rowmod, tdacs[col][row]);
								}
							}
							// Compute new TDAC value
							int delta = (int)(scn->getLoopVarValues(nloop))[scn->scanIndex(nloop)];
							if(fei4!=0) delta *= -1;
							if(PMG_DEBUG && col==40 && row==79) cout << "meas. thresh.: " << ((*hm)(colmod,rowmod)) <<
								", delta = " << delta << endl;
							if ((*hm)(colmod,rowmod) > trg) {
								newTDAC = tdacs[col][row] + corr - delta;
							} else {
								newTDAC = tdacs[col][row] + corr + delta;
							}
						}
						if(PMG_DEBUG && col==40 && row==79) cout << "new TDAC before cut-off:" << newTDAC <<
							" in step " << scn->scanIndex(nloop) << endl;
						if (newTDAC > maxtrim)                         newTDAC = maxtrim;
						if (newTDAC < 0 && scn->scanIndex(nloop) != 0) newTDAC = 0;
						if(PMG_DEBUG && col==40 && row==79) cout << "final new TDAC is " << newTDAC << " in step " << scn->scanIndex(nloop) << endl;
						if (newTDAC>=0) tdacs[col][row] = (unsigned int) newTDAC;
					}
				}
			}
			m_modules[pmod]->writeConfig(); //download the config again
		}
	}
}

void PixModuleGroup::prepareTDACFastTuning(int nloop, PixScan *scn) {
	if(PMG_DEBUG) cout << endl << endl << "**********************************\n"
		<< "start prepareTDACFastTuning for nloop == " << nloop << "\n******************************************" << endl;
	if (scn->scanIndex(nloop) > 5){
		cout << " \t !!! ABORTING !!!" << endl;
	} else {
		// Module loop
		for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
			if (m_modules[pmod]->m_readoutActive) {
				unsigned int mod = m_modules[pmod]->m_moduleId;
				int trg = scn->getRepetitions() * scn->getLoopVarNSteps(nloop-1) / 2;
				if(PMG_DEBUG) cout << "Total target Occ: " << trg << endl;
				Histo *hn, *hOcc, /**hSumOcc, */*hTdac;
				if (scn->scanIndex(nloop) == 0) {
					hOcc = new Histo("Occupancy", "Occupancy", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
					hTdac = new Histo("TDAC", "Tdac Tuning", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
					//hSumOcc = new Histo("SumOccupancy", "Sum of Occupancies", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
					// Associate different parameters if the tdac tuning is done in loop 0 or 1 or 2. Most likely 0 in fast tuning...
					if (nloop == 1) {
						scn->addHisto(*hOcc, PixScan::TDAC_OCC, mod, scn->scanIndex(2), -1, -1);
						scn->addHisto(*hTdac, PixScan::TDAC_T, mod, scn->scanIndex(2), -1, -1);
						//scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), -1, -1);
					} else if (nloop == 2) {
						scn->addHisto(*hOcc, PixScan::TDAC_OCC, mod, -1, -1, -1);
						scn->addHisto(*hTdac, PixScan::TDAC_T, mod, -1, -1, -1);
						//	scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), -1, -1);
					}
				} else {
					if (nloop == 1) {
						hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(nloop)-1, 0);
						hOcc = &scn->getHisto(PixScan::TDAC_OCC, mod, scn->scanIndex(2), 0, 0);
						hTdac = &scn->getHisto(PixScan::TDAC_T, mod, scn->scanIndex(2), 0, 0);
					} else if (nloop == 2) {
						hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(nloop)-1, 0, 0);
						hOcc = &scn->getHisto(PixScan::TDAC_OCC, mod, 0, 0, 0);
						hTdac = &scn->getHisto(PixScan::TDAC_T, mod, 0, 0, 0);
					}
				}
				// FE loop
				for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
					PixFe* fei4 = dynamic_cast<PixFeI4A*>(*fe);
					if(fei4==0) fei4 = dynamic_cast<PixFeI4B*>(*fe);
					int maxtrim;
					(*fe)->getTrimMax("TDAC", maxtrim);
					ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");

					int debugfe = (*fe)->number();
					if(PMG_DEBUG) cout << "\n\n############# Start FE loop. FE number: " << debugfe << " ####################" << endl;
					float vcal_a = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient3"])).value();
					float vcal_b = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient2"])).value();
					float vcal_c = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient1"])).value();
					float vcal_d = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient0"])).value();
					std::string capLabels[3]={"CInjLo", "CInjMed", "CInjHi"};
					float cInj     = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"][capLabels[scn->getChargeInjCap()]])).value();
					if(PMG_DEBUG) cout << "PixModuleGroup::prepareTDACTuning : using inj. capacitance of " << cInj << endl;
					if(PMG_DEBUG) cout << "PixModuleGroup::prepareTDACTuning : using PulserDAC slope of " << vcal_c << endl;
					//std::cout << std::dec << (*fe)->readGlobRegister("DAC_VCAL") << std::endl;
					float qBest=0;
					int vBest=0;
					for (int v=0; v<1024; v++) {
						float q = (vcal_a*v*v*v + vcal_b*v*v + vcal_c*v + vcal_d)*cInj/0.160218;
						if (abs(q-scn->getThresholdTargetValue()) < abs(qBest-scn->getThresholdTargetValue())) {
							qBest = q;
							vBest = v;
						}
					}
					//PixFe* fei4 = dynamic_cast<PixFeI4A*>(*fe);
					//if(fei4==0) fei4 = dynamic_cast<PixFeI4B*>(*fe);
					if(fei4!=0)
						(*fe)->writeGlobRegister("PlsrDAC", vBest);
					else
						(*fe)->writeGlobRegister("DAC_VCAL", vBest);
					if(PMG_DEBUG) cout << "best PulserDAC value is " << vBest << endl;
					// Bisection method
					for (unsigned int col=0; col<tdacs.get().size(); col++) {
						for (unsigned int row=0; row<tdacs[0].size(); row++) {
							if(PMG_DEBUG && col == 0 && row == 0) cout << "Col row loop: col low = " << col << " col high = "<< tdacs.get().size() << "row low = " << row << "row high = " << tdacs[0].size() << endl;

							unsigned int colmod, rowmod;
							int ife = (*fe)->number();
							// 	    if (ife < 8) {
							rowmod = row;
							colmod = col + NCOLS*ife;
							// 	    } else {
							// 	      rowmod = 319 - row;
							// 	      colmod = 143 - col - 18*(ife-8);
							// 	    }
							int newTDAC;
							if (scn->scanIndex(nloop) == 0) {
								newTDAC = (int)(scn->getLoopVarValues(nloop))[0];
							} else {
								int corr = 0;
								// Correct bisection for pixels in the tails
								if (scn->scanIndex(nloop) == 2 && fei4==0){ // correction makes no sense for FE-I4
									int occ0 = (int)((*hOcc)(colmod, rowmod));
									int occ1 = (int)((*hn)(colmod, rowmod));
									int td0 = (int)((scn->getLoopVarValues(nloop))[0]);
									int td1;
									if (occ0 < trg) {
										td1 = td0 -(int)((scn->getLoopVarValues(nloop))[1]);
									} else {
										td1 = td0 + (int)((scn->getLoopVarValues(nloop))[1]);
									}
									int delta = (int)((scn->getLoopVarValues(nloop))[2]);
									if (occ1 != occ0) {
										int ttt = td1 + (trg-occ1)*(td1-td0)/(occ1-occ0);
										if (ttt < 25) {
											ttt = 35;
										} else if (ttt > 117) {
											ttt = 107;
										} else if (ttt < 35) {
											ttt += 10;
										} else if (ttt > 100) {
											ttt -= 6;
										}
										if (td1 < td0 && ttt < td1) {
											corr = ttt - td1 + delta;
										} else if (td1 > td0 && ttt > td1) {
											corr = ttt - td1 - delta;
										}
									}
								}
								// Update best-mach histo
								if (scn->scanIndex(nloop) == 1) {
									if(PMG_DEBUG && col == 0 && rowmod == 0) cout << " !!! Setting best mach histo first time !!!" << endl;
									hOcc->set(colmod, rowmod, (*hn)(colmod,rowmod));
									hTdac->set(colmod, rowmod, tdacs[col][row]);
								} else {
									if (abs((*hn)(colmod,rowmod)-trg) < abs((*hOcc)(colmod,rowmod)-trg)) {
										hOcc->set(colmod, rowmod, (*hn)(colmod,rowmod));
										hTdac->set(colmod, rowmod, tdacs[col][row]);
									}
								}
								if(PMG_DEBUG && col == 39 && rowmod == 139) cout << "old tdac before computing: " << tdacs[col][row]<< endl;
								if(PMG_DEBUG && col == 39 && rowmod == 139) cout << "colmod: " << colmod << ", rowmod: " << rowmod << endl;
								// Compute new TDAC value
								if(PMG_DEBUG && col == 0 && rowmod == 0) cout << "---------- Starting to compute new TDAC values ----------------------" << endl;

								int delta = (int)(scn->getLoopVarValues(nloop))[scn->scanIndex(nloop)];

								if(fei4!=0) delta *= -1;

								if(PMG_DEBUG && col == 39 && rowmod == 139) cout << "Measured Occ is = " << (*hn)(colmod,rowmod) << endl;

								if ((*hn)(colmod,rowmod) < trg) {
									if(PMG_DEBUG && col == 39 && rowmod == 139) cout << "incrementing TDAC, delta = " << delta*-1 << ", old TDAC is " << tdacs[col][row] <<endl;
									newTDAC = tdacs[col][row] + corr - delta;
								} else {
									if(PMG_DEBUG && col == 39 && rowmod == 139) cout << "decrementing TDAC, delta = " << delta*-1 << ", old TDAC is " << tdacs[col][row] << endl;
									newTDAC = tdacs[col][row] + corr + delta;
								}
								if(PMG_DEBUG && colmod == 39 && rowmod == 139) cout << "New TDAC value before cutoff is " << newTDAC << endl;
							}
							if (newTDAC > maxtrim)                         newTDAC = maxtrim;
							if (newTDAC < 0 && scn->scanIndex(nloop) != 0) newTDAC = 0;
							if (newTDAC>=0){
								tdacs[col][row] = (unsigned int) newTDAC;
								if(PMG_DEBUG && col == 39 && rowmod == 139) cout << " --------------- Final new TDAC value is " << newTDAC << " --------------- " << endl;
							}
						}
					}
					if(PMG_DEBUG) cout << "Done with col, row, loop" << endl;
				}
				m_modules[pmod]->writeConfig(); //download the config again
			}
		}
		scn->setFeVCal(0x1fff);
	}
	if(PMG_DEBUG) cout << "ending prepareTDACFastTuning for nloop == " << nloop << endl << endl;
}

void PixModuleGroup::prepareGDACTuning(int nloop, PixScan *scn){
  // already covered by scan processing
//   // Read step value
//   int GDAC_FE = (int)(scn->getLoopVarValues(nloop)[scn->scanIndex(nloop)]);
//   // Modules loop
//   for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
//     if (m_modules[pmod]->m_readoutActive) {
//       // FEs loop
//       for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
// 	(*fe)->writeGlobRegister("GLOBAL_DAC",GDAC_FE);
//       }
//       m_modules[pmod]->writeConfig(); //download the config again

//     }
//   }
}

void PixModuleGroup::prepareGDACFastTuning(int nloop, PixScan *scn) {
  if(PMG_DEBUG) cout << endl << "start prepareGDACFastTuning for nloop == " << nloop << endl;
  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      int trg = scn->getRepetitions() * scn->getLoopVarNSteps(nloop-1) / 2;
      if(PMG_DEBUG) cout << "Total target Occ: " << trg << endl;
      Histo *hn, *hOcc, *hGdac;
      if(PMG_DEBUG) cout << "Scan index is " << scn->scanIndex(nloop) << endl;
      if (scn->scanIndex(nloop) == 0) {
	hOcc = new Histo("Occupancy", "Occupancy", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	hGdac = new Histo("GDAC", "Gdac Tuning", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	//hSumOcc = new Histo("SumOccupancy", "Sum of Occupancies", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	// Associate different parameters if the tdac tuning is done in loop 0 or 1 or 2. Most likely 0 in fast tuning...
	if (nloop == 1) {
	  if(PMG_DEBUG) cout << " adding histograms in step 0" << endl;
	  scn->addHisto(*hOcc, PixScan::GDAC_OCC, mod, scn->scanIndex(2), -1, -1);
	  if(PMG_DEBUG) cout << "added histo hOcc" << endl;
	  scn->addHisto(*hGdac, PixScan::GDAC_T, mod, scn->scanIndex(2), -1, -1);
	  if(PMG_DEBUG) cout << "added histo hGdac" << endl;
	} else if (nloop == 2) {
	  scn->addHisto(*hOcc, PixScan::GDAC_OCC, mod, -1, -1, -1);
	  scn->addHisto(*hGdac, PixScan::GDAC_T, mod, -1, -1, -1);
	  //	scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), -1, -1);
	}
      } else {
	if (nloop == 1) {
	  hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(nloop)-1, 0);
	  hOcc = &scn->getHisto(PixScan::GDAC_OCC, mod, scn->scanIndex(2), 0, 0);
	  hGdac = &scn->getHisto(PixScan::GDAC_T, mod, scn->scanIndex(2), 0, 0);
	} else if (nloop == 2) {
	  hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(nloop)-1, 0, 0);
	  hOcc = &scn->getHisto(PixScan::GDAC_OCC, mod, 0, 0, 0);
	  hGdac = &scn->getHisto(PixScan::GDAC_T, mod, 0, 0, 0);
	}
      }
      // FE loop
      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	PixFe* fei4a = dynamic_cast<PixFeI4A*>(*fe);
	PixFe* fei4b = dynamic_cast<PixFeI4B*>(*fe);
	ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC"); // M.B.: read tdacs to calculate max col and row number later. More elegant method vailable?
	// Get previous GDAC
	int gdac = 0;
	if(fei4a!=0){
	  gdac = (*fe)->readGlobRegister("Vthin_AltFine");
	  gdac += (*fe)->readGlobRegister("Vthin_AltCoarse")<<8;
	}else if(fei4b!=0){
	  gdac = (*fe)->readGlobRegister("Vthin_AltFine");
	  gdac += ((*fe)->readGlobRegister("Vthin_AltCoarse")&0xfe)<<7;
	}else
	  gdac = (*fe)->readGlobRegister("GLOBAL_DAC");
	
	if(PMG_DEBUG) cout << "\t read old gdac value: " << gdac << endl;
	
	float vcal_a = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient3"])).value();
	float vcal_b = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient2"])).value();
	float vcal_c = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient1"])).value();
	float vcal_d = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient0"])).value();
	std::string capLabels[3]={"CInjLo", "CInjMed", "CInjHi"};
	float cInj     = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"][capLabels[scn->getChargeInjCap()]])).value();
	if(PMG_DEBUG) cout << "PixModuleGroup::prepareGDACFastTuning : using inj. capacitance of " << cInj << endl;
	if(PMG_DEBUG) cout << "PixModuleGroup::prepareGDACFastTuning : using PulserDAC slope of " << vcal_c << endl;
	//std::cout << std::dec << (*fe)->readGlobRegister("DAC_VCAL") << std::endl;
	float qBest=0;
	int vBest=0;
	for (int v=0; v<1024; v++) {
	  float q = (vcal_a*v*v*v + vcal_b*v*v + vcal_c*v + vcal_d)*cInj/0.160218;
	  if (abs(q-scn->getThresholdTargetValue()) < abs(qBest-scn->getThresholdTargetValue())) {
	    qBest = q;
	    vBest = v;
	  }
	}
	if(fei4a!=0 || fei4b!=0)
	  (*fe)->writeGlobRegister("PlsrDAC", vBest);
	else
	  (*fe)->writeGlobRegister("DAC_VCAL", vBest);
	if(PMG_DEBUG) cout << "best PulserDAC value is " << vBest << endl;
	
	int newGDAC;
	// get new GDAC from Scan parameter, if loop is executed first time
	if (scn->scanIndex(nloop) == 0) {
	  newGDAC = (int)(scn->getLoopVarValues(nloop))[0];
	} else {
	  
	  // Calculate mean of Occ
	  int hnMean = 0;
	  int nrpixels = 0;
	  for (unsigned int col=0; col<tdacs.get().size(); col++) {
	    for (unsigned int row=0; row<tdacs[0].size(); row++) {
	      unsigned int colmod, rowmod;
	      int ife = (*fe)->number();
	      // 	    if (ife < 8) {
	      rowmod = row;
	      colmod = col + NCOLS*ife;
	      hnMean += (int)(*hn)(colmod, rowmod);
	      nrpixels++;
	    }
	  }
	  if (nrpixels != 0) {
	    hnMean = (int) hnMean / nrpixels;
	    if(PMG_DEBUG) cout << " \t\t MeanOcc = " << hnMean << endl;
	  } else {
	    if(PMG_DEBUG) cout << "!!! Error, nrpixels is 0, set OccMean to 0 !!!" << endl;
	    hnMean = 0;
	  }
	  for (unsigned int col=0; col<tdacs.get().size(); col++) {
	    for (unsigned int row=0; row<tdacs[0].size(); row++) {
	      unsigned int colmod, rowmod;
	      int ife = (*fe)->number();
	      // 	    if (ife < 8) {
	      rowmod = row;
	      colmod = col + NCOLS*ife;
	      // Update best-mach histo
	      if (scn->scanIndex(nloop) == 1) {
		hOcc->set(colmod, rowmod, hnMean);
		hGdac->set(colmod, rowmod, gdac);
	      } else {
		if (abs(hnMean-trg) < abs((*hOcc)(colmod,rowmod)-trg)) {
		  hOcc->set(colmod, rowmod, hnMean);
		  hGdac->set(colmod, rowmod, gdac);
		}
	      }
	    }
	  }
	  if(PMG_DEBUG) cout << "old gdac before computing: " << gdac << endl;
	  // Compute new GDAC value
	  if(PMG_DEBUG) cout << "---------- Starting to compute new GDAC value ----------------------" << endl;
	  int delta = (int)(scn->getLoopVarValues(nloop))[scn->scanIndex(nloop)];
	  if(PMG_DEBUG) cout << "Measured Mean Occ is = " << hnMean << endl;
	  
	  if (hnMean > trg) {
	    if(PMG_DEBUG) cout << "incrementing GDAC, delta = " << delta << ", old GDAC is " << gdac << endl;
	    newGDAC = gdac + delta;
	  } else {
	    if(PMG_DEBUG) cout << "decrementing GDAC, delta = " << delta << ", old GDAC is " << gdac << endl;
	    newGDAC = gdac - delta;
	  }
	  if(PMG_DEBUG) cout << "New GDAC value before cutoff is " << newGDAC << endl;
	}
	
	if (newGDAC > 255)                         newGDAC = 255;
	if (newGDAC < 0 && scn->scanIndex(nloop) != 0) newGDAC = 0;
	if (newGDAC>=0){
	  gdac = (unsigned int) newGDAC;
	  if(PMG_DEBUG) cout << " --------------- Final new GDAC value is " << newGDAC << " --------------- " << endl;
	}
	if(fei4a!=0){
	  (*fe)->writeGlobRegister("Vthin_AltFine", gdac&0xff);
	  (*fe)->writeGlobRegister("Vthin_AltCoarse", (gdac&0xff00)>>8);
	}else if(fei4b!=0){
	  (*fe)->writeGlobRegister("Vthin_AltFine", gdac&0xff);
	  (*fe)->writeGlobRegister("Vthin_AltCoarse", (gdac&0x7f00)>>7);
	}else
	  (*fe)->writeGlobRegister("GLOBAL_DAC", gdac);
      }
    }
    m_modules[pmod]->writeConfig(); //download the config again
  }
  if(PMG_DEBUG) cout << "ending prepareTDACFastTuning for nloop == " << nloop << endl << endl;
}

void PixModuleGroup::prepareFDACTuning(int nloop, PixScan *scn){
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	float vcal_a = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient3"])).value();
	float vcal_b = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient2"])).value();
	float vcal_c = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient1"])).value();
	float vcal_d = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient0"])).value();
	std::string capLabels[3]={"CInjLo", "CInjMed", "CInjHi"};
	float cInj     = (dynamic_cast<ConfFloat &>((*fe)->config()["Misc"][capLabels[scn->getChargeInjCap()]])).value();
	if(PMG_DEBUG) cout << "PixModuleGroup::prepareFDACTuning : using inj. capacitance of " << cInj << endl;
	//std::cout << std::dec << (*fe)->readGlobRegister("DAC_VCAL") << std::endl;
	float qBest=0;
        int vBest=0;
	for (int v=0; v<1024; v++) {
	  float q = (vcal_a*v*v*v + vcal_b*v*v + vcal_c*v + vcal_d)*cInj/0.160218;
	  if (abs(q-scn->getTotTargetCharge()) < abs(qBest-scn->getTotTargetCharge())) {
	    qBest = q;
            vBest = v;
	  }
	}
	PixFe* fei4 = dynamic_cast<PixFeI4A*>(*fe);
	if(fei4==0) fei4 = dynamic_cast<PixFeI4B*>(*fe);
	if(fei4!=0)
	  (*fe)->writeGlobRegister("PlsrDAC", vBest);
	else
	  (*fe)->writeGlobRegister("DAC_VCAL", vBest);
	//std::cout << vBest << " " << (*fe)->readGlobRegister("DAC_VCAL") << std::endl;
      }
      m_modules[pmod]->writeConfig(); //download the config again
    }
  }
  scn->setFeVCal(0x1fff);
}

void PixModuleGroup::prepareIFTuning(int nloop, PixScan *scn){
  prepareFDACTuning(nloop, scn);
}

void PixModuleGroup::prepareT0Set(int nloop, PixScan *scn){
  //prepareFDACTuning(nloop, scn);
}

void PixModuleGroup::prepareIncrTdac(int nloop, PixScan *scn)
{
  bool debug=false;
  int vidx=0;
  if (!scn->getLoopVarValuesFree(nloop)) vidx = scn->scanIndex(nloop);

  int incr = (int)(scn->getLoopVarValues(nloop))[vidx];
  incr = abs(incr);

  if(debug) cout << "read increment " << incr << " at scan step " << scn->scanIndex(nloop) << endl;

  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++)
    {
      if (m_modules[pmod]->m_readoutActive)
	{
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  Histo *hm, *hs, *hc, *hTdac;
	  int maxtrim;
	  (*(m_modules[pmod]->feBegin()))->getTrimMax("TDAC", maxtrim);

	  if (scn->scanIndex(nloop) == 0)
	    {
	      if(debug) cout << "scanIndex ist hier " << scn->scanIndex(nloop) << endl;

	      m_modules[pmod]->storeConfig("IncrTdacPreviousConfig");
	      hTdac = new Histo("TDAC", "Tdac Tuning", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	      // Attach the result histogram in the correct place if the scan is done in loop 1 or 2
	      if (nloop == 1)
		{
		  scn->addHisto(*hTdac, PixScan::TDAC_T, mod, scn->scanIndex(2), -1, -1);
		}
	      else if (nloop == 2)
		{
		  scn->addHisto(*hTdac, PixScan::TDAC_T, mod, -1, -1, -1);
		}
	    }

	  else
	    {
	      if (nloop == 1)
		{
		  hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1)-1, 0);
		  hs = &scn->getHisto(PixScan::SCURVE_SIGMA, mod, scn->scanIndex(2), scn->scanIndex(1)-1, 0);
		  hc = &scn->getHisto(PixScan::SCURVE_CHI2, mod, scn->scanIndex(2), scn->scanIndex(1)-1, 0);
		  hTdac = &scn->getHisto(PixScan::TDAC_T, mod, scn->scanIndex(2), 0, 0);
		}

	      else if (nloop == 2)
		{
		  hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2)-1, 0, 0);
		  hs = &scn->getHisto(PixScan::SCURVE_SIGMA, mod, scn->scanIndex(2), 0, 0);
		  hc = &scn->getHisto(PixScan::SCURVE_CHI2, mod, scn->scanIndex(2), 0, 0);
		  hTdac = &scn->getHisto(PixScan::TDAC_T, mod, 0, 0, 0);
		}
	    }
	  // FE loop
	  /*
	    if (!scn->getLoopVarValuesFree(nloop)) {
	    for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	    ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");
	    ConfMask<unsigned short int> &tdacsPre = ((*fe)->feConfig("IncrTdacPreviousConfig")).readTrim("TDAC");
	    for (unsigned int col=0; col<tdacs.get().size(); col++) {
	    for (unsigned int row=0; row<tdacs[0].size(); row++) {
	    unsigned int colmod, rowmod;
	    int ife = (*fe)->number();
	    if (ife < 8) {
	    rowmod = row;
	    colmod = col + 18*ife;
	    } else {
	    rowmod = 319 - row;
	    colmod = 143 - col - 18*(ife-8);
	    }
	    if ((*hs)(colmod,rowmod) == 0 || (*hm)(colmod,rowmod) == 0 || (*hc)(colmod,rowmod)>20 || (*hc)(colmod,rowmod)<0) {
	    tdacs[col][row] = tdacsPre[col][row];
	    }
	    }
	    }
	    }
	    m_modules[pmod]->storeConfig("IncrTdacPreviousConfig");
	    }
	  */


	  int modif = 0;
	  int zerofit = 0, tozero = 0, td00 = 0, tdp00 = 0;

	  for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++)
	    {
	      ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");
	      ConfMask<unsigned short int> &tdacsPre = ((*fe)->feConfig("IncrTdacPreviousConfig")).readTrim("TDAC");

	      for (unsigned int col=0; col<tdacs.get().size(); col++)
		{
		  for (unsigned int row=0; row<tdacs[0].size(); row++)
		    {
		      unsigned int colmod, rowmod;
		      int ife = (*fe)->number();

// 		      if (ife < 8)
// 			{
			  rowmod = row;
			  colmod = col + NCOLS*ife;
// 			}
// 		      else
// 			{
// 			  rowmod = 319 - row;
// 			  colmod = 143 - col - 18*(ife-8);
// 			}

		      if (scn->getLoopVarValuesFree(nloop))
			{
			  if (scn->scanIndex(nloop) == 0) //Wird nur beim ersten Mal betreten
			    {
			      if (tdacs[col][row] - incr < 0) //Der Fall sollte nicht eintreten
				{
				  if(col==0 && row==0)
				    {
				      if(debug) cout << "if 1" << endl;
				    }
				  tdacsPre[col][row] = tdacs[col][row];
				  tdacs[col][row] = 0;
				}
			      else // Hier werden die incr Werte eingelsen
				{
				  tdacsPre[col][row] = incr;
				  tdacs[col][row] -= incr;

				  if(col==0 && row==0)
				    {
				      if(debug) cout << "else 1" << endl;
				    }

				}
			      modif++;             //neu

			      if(col==0 && row==0)
				{
				  if(debug) cout << "if 0" << endl;
				}

			    }
			  //} //neu
			  else    //Wird sonst immer eingelesen
			    {

			      if(col==0 && row==0)
				{
				  if(debug) cout << "else 2 " << tdacs[col][row] << " " << tdacsPre[col][row] <<  " chi2 " << (*hc)(colmod,rowmod)<< " hs " << (*hs)(colmod,rowmod) << " hm " << (*hm)(colmod,rowmod) << endl;
				}
			      int ptype=0;
// 			      if(row==159 || row==157 || row==155 || row==153) ptype=2; // ganged pixel
// 			      else if(row==158 || row==156 || row==154)        ptype=3; // inter-ganged pixel
// 			      else if(col==0 || col==17)                       ptype=1; // long pixel
			      float noisecut = scn->getMinThrNoiseCut(ptype);
			      // org. version
			      //			      if ((*hs)(colmod,rowmod) <1 || (*hs)(colmod,rowmod) >noisecut || (*hm)(colmod,rowmod) <1 || (*hm)(colmod,rowmod) >15000 || (*hc)(colmod,rowmod)>50 || (*hc)(colmod,rowmod)<1)  //Abbruchbedingung, wann der TDAC Wert zu niedrig ist
			      // new decrement-stop condition after chi2-fix
			      if ((*hs)(colmod,rowmod) >noisecut || (*hm)(colmod,rowmod) <1 || (*hc)(colmod,rowmod)>scn->getMinThrChi2Cut() || (*hc)(colmod,rowmod)<=0)
				{
				  tdacs[col][row] += tdacsPre[col][row];  //TDACs wieder hoch setzten

				  if(col==0 && row==0)
				    {
				      if(debug) cout << "if 5 " << tdacs[col][row] << " " << tdacsPre[col][row] <<  "chi2 " << (*hc)(colmod,rowmod) <<endl;
				    }

				  if (tdacsPre[col][row] > 0)
				    {
				      tdacsPre[col][row] -= 1; //incr kleiner machen

				      if(col==0 && row==0)
					{
					  if(debug) cout << "if 2 " << tdacs[col][row] << " " << tdacsPre[col][row] <<endl;
					}
				    }
				  zerofit++;  //Zahl der nicht gelungenen Fits
				}
			      /*
				 if ((*hm)(colmod,rowmod) > 0 && (*hm)(colmod,rowmod)-4*(*hs)(colmod,rowmod) < 0) {
				 tdacsPre[col][row] = 0;
				 tozero++;
				 } else {
			      */
			      if (tdacsPre[col][row] > 0) //Sollte so sein
				{
				  if (tdacs[col][row] - tdacsPre[col][row] < 0) // Sollte nicht eintrten
				    {
				      if(col==0 && row==0)
					{
					  if(debug) cout << "if 3 " << tdacs[col][row] << " " << tdacsPre[col][row] << endl;
					}

				      tdacsPre[col][row] = 0;
				      tdacs[col][row] = 0;
				    }

				  else  // Sollte eintreten
				    {

				      if(col==0 && row==0)
					{
					  if(debug) cout << "else 3 " << tdacs[col][row] << " " << tdacsPre[col][row] <<endl;
					}

				      tdacs[col][row] -= tdacsPre[col][row]; //TDACs kleiner machen

				      if(col==0 && row==0)
					{
					  if(debug) cout << "if 4 " << tdacs[col][row] << " " << tdacsPre[col][row] <<endl;
					}

				      modif++;

				      if(col==0 && row==0)
					{
					  if(debug) cout << "if 0 modif=" << modif << endl;
					}
				    }
				}
			    }

			  tdp00 = tdacsPre[0][0];
			  td00 = tdacs[0][0];
			  hTdac->set(colmod, rowmod, tdacs[col][row]);
			}

		      else //Tritt nicht ein
			{

			  if(col==0 && row==0)
			    {
			      if(debug) cout << "else 4 " << tdacs[col][row] << " " << tdacsPre[col][row] <<endl;
			    }

			  if (scn->scanIndex(nloop) == 0)
			    {
			      if (tdacs[col][row] + incr > maxtrim)
				{
				  tdacs[col][row] = maxtrim;
				}
			      else if (tdacs[col][row] + incr < 0)
				{
				  tdacs[col][row] = 0;
				}
			      else
				{
				  tdacs[col][row] += incr;
				}
			    }
			  else
			    {
			      if ((*hm)(colmod,rowmod) > 0 && (*hm)(colmod,rowmod)-4*(*hs)(colmod,rowmod) > 0)
				{
				  if (tdacs[col][row] + incr > maxtrim)
				    {
				      tdacs[col][row] = maxtrim;
				    }
				  else if (tdacs[col][row] + incr < 0)
				    {
				      tdacs[col][row] = 0;
				    }
				  else
				    {
				      tdacs[col][row] += incr;
				    }
				}
			    }

			  hTdac->set(colmod, rowmod, tdacs[col][row]);
			}



		    }
		}
	    }

	  if (scn->getLoopVarValuesFree(nloop))
	    {
	      if (modif == 0) scn->terminate(nloop); // Abbruch, wenn alle tdacPre=0????

	      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++)
		{
		  ConfMask<unsigned short int> &tdacs2 = (*fe)->readTrim("TDAC");

		  if(debug) cout << "tdac " << tdacs2[0][0] << endl;

		}

	      if(debug) std::cout << "+++" << std:: dec << modif << " " << zerofit << " " << tozero << " " << td00 << " " << tdp00 << std::endl;
	    }



	  m_modules[pmod]->writeConfig(); //download the config again
	  // }
	}
    }
}

void PixModuleGroup::endTDACTuning(int nloop, PixScan *scn) {
	if(PMG_DEBUG) cout << "Starting PixModuleGroup::endTDACTuning" << endl;
	// Check if the loop is executed on the host
	if (!scn->getDspLoopAction(nloop)) {
		Histo *hm, *hThr, *hTdac;
		int trg = scn->getThresholdTargetValue();
		// Loop on modules
		for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
			if (m_modules[pmod]->m_readoutActive) {
				unsigned int mod = m_modules[pmod]->m_moduleId;
				if (nloop == 1) {
					hThr = &scn->getHisto(PixScan::TDAC_THR, mod, scn->scanIndex(2), 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, scn->scanIndex(2), 0, 0);
					hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), 0);
				} else if (nloop == 2) {
					hThr = &scn->getHisto(PixScan::TDAC_THR, mod, 0, 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, 0, 0, 0);
					hm = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), 0, 0);
				}
				// FE loop
				for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
					ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");
					for (unsigned int col=0; col<tdacs.get().size(); col++) {
						for (unsigned int row=0; row<tdacs[0].size(); row++) {
							unsigned int colmod, rowmod;
							int ife = (*fe)->number();
							// 	      if (ife < 8) {
							rowmod = row;
							colmod = col + NCOLS*ife;
							// 	      } else {
							// 		rowmod = 319 - row;
							// 		colmod = 143 - col - 18*(ife-8);
							// 	      }
							if (abs((*hm)(colmod,rowmod)-trg) < abs((*hThr)(colmod,rowmod)-trg)) {
								hThr->set(colmod, rowmod, (*hm)(colmod,rowmod));
								hTdac->set(colmod, rowmod, tdacs[col][row]);
							}
							// Save in the FEs TDAC tuned values
							tdacs[col][row] = (int)(*hTdac)(colmod,rowmod);
						}
					}
				}
			}
		}
	}
	if(PMG_DEBUG) cout << "End of PixModuleGroup::endTDACTuning" << endl;
}

void PixModuleGroup::endTDACFastTuning(int nloop, PixScan *scn) {
	cout << "Starting PixModuleGroup::endTDACFastTuning" << endl;
	// Check if the loop is executed on the host
	if (!scn->getDspLoopAction(nloop)) {
		Histo *hn, *hOcc, *hTdac;
		int trg = scn->getRepetitions() * scn->getLoopVarNSteps(nloop-1) / 2;
		if(PMG_DEBUG) cout << "Total target Occ: " << trg << endl;
		// Loop on modules
		for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
			if (m_modules[pmod]->m_readoutActive) {
				unsigned int mod = m_modules[pmod]->m_moduleId;
				if (nloop == 1) {
					hOcc = &scn->getHisto(PixScan::TDAC_OCC, mod, scn->scanIndex(2), 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, scn->scanIndex(2), 0, 0);
					hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(nloop), 0);
				} else if (nloop == 2) {
					hOcc = &scn->getHisto(PixScan::TDAC_OCC, mod, 0, 0, 0);
					hTdac = &scn->getHisto(PixScan::TDAC_T, mod, 0, 0, 0);
					hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(nloop), scn->scanIndex(nloop), 0);
				}
				// FE loop
				for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
					ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC");
					for (unsigned int col=0; col<tdacs.get().size(); col++) {
						for (unsigned int row=0; row<tdacs[0].size(); row++) {
							unsigned int colmod, rowmod;
							int ife = (*fe)->number();
							// 	      if (ife < 8) {
							rowmod = row;
							colmod = col + NCOLS*ife;
							// 	      } else {
							// 		rowmod = 319 - row;
							// 		colmod = 143 - col - 18*(ife-8);
							// 	      }
							if (abs((*hn)(colmod,rowmod)-trg) < abs((*hOcc)(colmod,rowmod)-trg)) {
								hOcc->set(colmod, rowmod, (*hn)(colmod,rowmod));
								hTdac->set(colmod, rowmod, tdacs[col][row]);
							}
							// Save in the FEs TDAC tuned values
							tdacs[col][row] = (int)(*hTdac)(colmod,rowmod);
						}
					}
				}
			}
		}
	}
	cout << "End of PixModuleGroup::endTDACFastTuning" << endl;
}

void PixModuleGroup::endOccSumming(int nloop, PixScan *scn) {
	if(PMG_DEBUG) cout << "Starting PixModuleGroup::endOccSumming for nloop == " << nloop << endl;
	// Check if the loop is executed on the host
	if (!scn->getDspLoopAction(nloop)) {
		Histo *hOcc, *hSumOcc;
		// Loop on modules
		for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
			if (m_modules[pmod]->m_readoutActive) {
				unsigned int mod = m_modules[pmod]->m_moduleId;

				//if (scn->scanIndex(nloop+1) == scn->getLoopVarNSteps(nloop+1)) {
				if(PMG_DEBUG) cout << endl << "Creating OccSum" << endl << endl;
				hSumOcc = new Histo("SumOccupancy", "Sum of Occupancies", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
				// Associate different parameters if the tdac tuning is done in loop 0 or 1 or 2. Most likely 0 in fast tuning...
				if (nloop == 0) {
					scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), scn->scanIndex(1), scn->scanIndex(nloop));
				} else if (nloop == 1) {
					scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), -1, -1);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), scn->scanIndex(nloop), 0);
				} else if (nloop == 2) {
					scn->addHisto(*hSumOcc, PixScan::SUM_OCC, mod, scn->scanIndex(2), -1, -1);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(nloop), 0, 0);
				}
				//} else {
				if (nloop == 0) {
					hSumOcc = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(1), 0);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), scn->scanIndex(1), scn->scanIndex(nloop));
				} else if (nloop == 1) {
					hSumOcc = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), 0, 0);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), scn->scanIndex(nloop), 0);
				} else if (nloop == 2) {
					hSumOcc = &scn->getHisto(PixScan::SUM_OCC, mod, 0, 0, 0);
					hOcc = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(nloop), 0, 0);
				}
				//}

				// Loop on steps (finding the mean threshold in each FE)
				for (int step = 0; step<scn->getLoopVarNSteps(nloop); step++){
					if (nloop == 0){
						hOcc      = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), scn->scanIndex(1), step);
					}else if (nloop == 1){
						hOcc      = &scn->getHisto(PixScan::OCCUPANCY, mod, scn->scanIndex(2), step, 0);
					}else if (nloop == 2){
						hOcc      = &scn->getHisto(PixScan::OCCUPANCY, mod, step, 0, 0);
					}
					// Loop on FEs
					for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
						unsigned int colmod, rowmod, colfe, rowfe;
						double mean = 0, meanlocal = 0;
						int ife, npixel=0;
						// Cols and rows are the pixel coordinates in FE level
						ife = (*fe)->number();
						for (unsigned int col=0; col<NCOLS; col++){
							for (unsigned int row=0; row<NROWS; row++){
								// 		if (ife < 8) {
								rowmod = row;
								colmod = col + NCOLS*ife;
								//sum histograms...
								int oldVal = (*hSumOcc)(colmod,rowmod);
								int newOcc = (*hOcc)(colmod,rowmod);
								int newVal = oldVal + newOcc;
								hSumOcc->set(colmod, rowmod, newVal);
								if(PMG_DEBUG && colmod == 39 && rowmod == 139) cout << "Summed histogram: oldSum = " << oldVal	<< ", newOcc = " << newOcc << ", newSum = " << newVal << endl;
							}
						}
					}
				}
			}
		}
		if(PMG_DEBUG) cout << "Ending PixModuleGroup::endOccSumming" << endl;
	}
}

void PixModuleGroup::endGDACTuning(int nloop, PixScan *scn) {
  // Check if the loop is executed on the host
  if (!scn->getDspLoopAction(nloop)) {
    // Modules loop
    for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
      if (m_modules[pmod]->m_readoutActive) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	Histo *hm, *hThr, *hGdac, *hm1, *hm2;
	hThr = new Histo("Threshold", "Threshold", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	hGdac = new Histo("GDAC", "Gdac Tuning", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	// Associate different parameters if the gdac tuning is done in loop 0, 1 or 2
	if (nloop == 0) {
	  scn->addHisto(*hThr, PixScan::GDAC_THR, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	  scn->addHisto(*hGdac, PixScan::GDAC_T, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	} else if (nloop == 1) {
	  scn->addHisto(*hThr, PixScan::GDAC_THR, mod, scn->scanIndex(2), -1, -1);
	  scn->addHisto(*hGdac, PixScan::GDAC_T, mod, scn->scanIndex(2), -1, -1);
	} else if (nloop == 2) {
	  scn->addHisto(*hThr, PixScan::GDAC_THR, mod, -1, -1, -1);
	  scn->addHisto(*hGdac, PixScan::GDAC_T, mod, -1, -1, -1);
	}
	// Create matrix for storing FE mean threshold values for each step
#ifdef WIN32
	std::vector<float> meanmat[16];
	for(int iv=0; iv<16;iv++)
	  meanmat[iv].resize(scn->getLoopVarNSteps(nloop));
#else
	float meanmat[16][scn->getLoopVarNSteps(nloop)];
#endif
	// Create matrix for storing GDAC values for each step
#ifdef WIN32
	std::vector<int> GDACmat[16];
	for(int iv=0; iv<16;iv++)
	  GDACmat[iv].resize(scn->getLoopVarNSteps(nloop));
#else
	int GDACmat[16][scn->getLoopVarNSteps(nloop)];
#endif
	// Get target threshold
	int trg = scn->getThresholdTargetValue();
	// Create vectors and matrix for storing FE mean H&L threshold values and GDAC H&L
	int thrFEL[16], thrFEH[16], GdacL[16], GdacH[16], H[16], L[16];
	// Load protection values in each FE
	int  b, c, d, e;
	b = 500; // Min Threshold
	c = 20000; // Max Threshold
	d = 0;    // Min GDAC
	e = 400;   // Max GDAC
	for (int i=0; i<16; i++){   // Loop on vectorFE
	  thrFEL[i] = b;
	  thrFEH[i] = c;
	  GdacL[i] = d;
	  GdacH[i] = e;
	}
	// Loop on steps (finding the mean threshold in each FE)
	for (int step = 0; step<scn->getLoopVarNSteps(nloop); step++){
	  // Verse determination of GDAC values
	  int dir = 0; // Increasing GDAC values
	  if (scn->getLoopVarValues(nloop)[0] > scn->getLoopVarValues(nloop)[1]) dir = 1; // Decreasing GDAC values ordination
	  // Check if gdac tuning is done in loop 1 or 2
	  if (nloop == 1){
	    hm      = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step, 0);
	    hThr    = &scn->getHisto(PixScan::GDAC_THR, mod, scn->scanIndex(2), step, 0);
	    hGdac   = &scn->getHisto(PixScan::GDAC_T, mod, scn->scanIndex(2), step, 0);
	  }else if (nloop == 2){
	    hm      = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step, 0, 0);
	    hThr    = &scn->getHisto(PixScan::GDAC_THR, mod, step, 0, 0);
	    hGdac   = &scn->getHisto(PixScan::GDAC_T, mod, step, 0, 0);
	  }
	  // Loop on FEs
	  for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	    PixFe* fei4a = dynamic_cast<PixFeI4A*>(*fe);
	    PixFe* fei4b = dynamic_cast<PixFeI4B*>(*fe);
	    unsigned int colmod, rowmod, colfe, rowfe;
	    double mean = 0, meanlocal = 0;
	    int ife, npixel=0;
	    // Cols and rows are the pixel coordinates in FE level
	    ife = (*fe)->number();
	    for (unsigned int col=0; col<NCOLS; col++){
	      for (unsigned int row=0; row<NROWS; row++){
// 		if (ife < 8) {
		  rowmod = row;
		  colmod = col + NCOLS*ife;
// 		} else {
// 		  rowmod = 319 - row;
// 		  colmod = 143 - col - 18*(ife-8);
// 		}
		meanlocal = (int)(*hm)(colmod,rowmod); // Mean value in each pixel
		if (meanlocal > 100 && meanlocal<c) {
		  npixel++;
		  mean = mean + meanlocal;
		}
	      }
	    }
	    if (npixel != 0) mean = mean/npixel;
	    // Save values in matrix
	    meanmat[ife][step] = (int)mean; // Mean threshold values for each FE in each step
	    GDACmat[ife][step] = (int)(scn->getLoopVarValues(nloop))[step]; // GDAC values for each FE in each step
	    // Fill histo of mean threshold for each FE
// 	    if (ife < 8) {
	      colfe = ife;
	      rowfe = 0;
// 	    } else {
// 	      colfe = 15-ife;
// 	      rowfe = 1;
// 	    }
	    if (mean <= trg && mean>thrFEL[ife]) {  // Update low value
	      thrFEL[ife] = (int)mean;
	      GdacL[ife] = (int)(scn->getLoopVarValues(nloop)[step]);
	      L[ife] = step;
	    } else if (mean>trg && mean<thrFEH[ife]) { // Update high value
	      thrFEH[ife] = (int)mean;
	      GdacH[ife] = (int)(scn->getLoopVarValues(nloop)[step]);
	      H[ife] = step;
	    }
	    // If you are in the last step, you have to interpolate between Low and High values to find the tuned value of Gdac
	    if (step == (scn->getLoopVarNSteps(nloop)-1)) {
	      int GdacI = 0;
	      float FinalThr = 0;
	      int test = 0;
	      // Check FE status
	      for (int i=0; i<scn->getLoopVarNSteps(nloop); i++){
		if (meanmat[ife][i] == 0) test = test+1;
	      }
	      // If FE is working
	      if (test != scn->getLoopVarNSteps(nloop)){
		// If scan has found a value higher than trg and one lower
		if (thrFEL[ife]<=trg && thrFEH[ife]>trg && thrFEL[ife]!=b && thrFEH[ife]!=c ) {
		  GdacI = (int)(GdacL[ife] + (trg - thrFEL[ife]) * (GdacH[ife]-GdacL[ife]) / (thrFEH[ife] - thrFEL[ife]));
		  if (nloop == 1){
		    hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), L[ife], 0);
		    hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), H[ife], 0);
		  } else if (nloop == 2){
		    hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, L[ife], 0, 0);
		    hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, H[ife], 0, 0);
		  }
		  for (unsigned int col=0; col<NCOLS; col++){
		    for (unsigned int row=0; row<NROWS; row++){
// 		      if (ife < 8) {
			rowmod = row;
			colmod = col + NCOLS*ife;
// 		      } else {
// 			rowmod = 319 - row;
// 			colmod = 143 - col - 18*(ife-8);
// 		      }
		      if ((*hm2)(colmod,rowmod) != 0 && (*hm1)(colmod,rowmod) != 0){
			FinalThr = (*hm1)(colmod,rowmod) + (GdacI - GdacL[ife]) * ((*hm2)(colmod,rowmod)-(*hm1)(colmod,rowmod))/(GdacH[ife]-GdacL[ife]);
		      } else FinalThr = 0;
		      hGdac->set(colmod, rowmod, GdacI);
		      hThr->set(colmod, rowmod, FinalThr);
		    }
		  }
		} else if (thrFEL[ife] == b){  // If scan has not found a lower value
		  //Check Gdac values ordination
		  if (dir == 0 && (meanmat[ife][1] - meanmat[ife][0]) != 0) {
		    GdacI = (int)(GDACmat[ife][0] + (trg - meanmat[ife][0])*(GDACmat[ife][1] - GDACmat[ife][0])/(meanmat[ife][1] - meanmat[ife][0]));
		    if (nloop == 1){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), 1, 0);
		    } else if (nloop == 2){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, 0, 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, 1, 0, 0);
		    }
		    for (unsigned int col=0; col<NCOLS; col++){
		      for (unsigned int row=0; row<NROWS; row++){
// 			if (ife < 8) {
			  rowmod = row;
			  colmod = col + NCOLS*ife;
// 			} else {
// 			  rowmod = 319 - row;
// 			  colmod = 143 - col - 18*(ife-8);
// 			}
			if ((*hm2)(colmod,rowmod) != 0 && (*hm1)(colmod,rowmod) != 0){
			  FinalThr = (*hm1)(colmod,rowmod) + (GdacI-GDACmat[ife][0]) * ((*hm2)(colmod,rowmod)-(*hm1)(colmod,rowmod)) / (GDACmat[ife][1] - GDACmat[ife][0]);
			} else if ((*hm2)(colmod,rowmod) == 0 && (*hm1)(colmod,rowmod) != 0) {
			  FinalThr = (*hm1)(colmod,rowmod);
			} else FinalThr = 0;
			hGdac->set(colmod, rowmod, GdacI);
			hThr->set(colmod, rowmod, FinalThr);
		      }
		    }
		  } else if (dir == 1 && (meanmat[ife][step-1] - meanmat[ife][step]) != 0) {
		    GdacI = (int)(GDACmat[ife][step] + (trg - meanmat[ife][step])*(GDACmat[ife][step-1] - GDACmat[ife][step])/(meanmat[ife][step-1] - meanmat[ife][step]));
		    if (nloop == 1){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step-1, 0);
		    } else if (nloop == 2){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step, 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step-1, 0, 0);
		    }
		    for (unsigned int col=0; col<NCOLS; col++){
		      for (unsigned int row=0; row<NROWS; row++){
// 			if (ife < 8) {
			  rowmod = row;
			  colmod = col + NCOLS*ife;
// 			} else {
// 			  rowmod = 319 - row;
// 			  colmod = 143 - col - 18*(ife-8);
// 			}
			if ((*hm2)(colmod,rowmod) != 0 && (*hm1)(colmod,rowmod) != 0){
			  FinalThr = (*hm1)(colmod,rowmod)+ (GdacI-GDACmat[ife][step])*((*hm2)(colmod,rowmod)-(*hm1)(colmod,rowmod))/(GDACmat[ife][step-1]-GDACmat[ife][step]);
			} else if ((*hm2)(colmod,rowmod) == 0 && (*hm1)(colmod,rowmod) != 0) {
			  FinalThr = (*hm1)(colmod,rowmod);
			} else FinalThr = 0;
			hGdac->set(colmod, rowmod, GdacI);
			hThr->set(colmod, rowmod, FinalThr);
		      }
		    }
		  }
		} else if (thrFEH[ife] == c){  // If scan has not found a higher value
		  if (dir == 0 && (meanmat[ife][step] - meanmat[ife][step-1]) != 0) {
		    GdacI = (int)(GDACmat[ife][step] + (trg - meanmat[ife][step])*(GDACmat[ife][step-1] - GDACmat[ife][step])/(meanmat[ife][step-1]-meanmat[ife][step]));
		    if (nloop == 1){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step-1, 0);
		    } else if (nloop == 2){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step, 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step-1, 0, 0);
		    }
		    for (unsigned int col=0; col<NCOLS; col++){
		      for (unsigned int row=0; row<NROWS; row++){
// 			if (ife < 8) {
			  rowmod = row;
			  colmod = col + NCOLS*ife;
// 			} else {
// 			  rowmod = 319 - row;
// 			  colmod = 143 - col - 18*(ife-8);
// 			}
			if ((*hm2)(colmod,rowmod) != 0 && (*hm1)(colmod,rowmod) != 0){
			  FinalThr = (*hm1)(colmod,rowmod)+ (GdacI-GDACmat[ife][step])*((*hm2)(colmod,rowmod)-(*hm1)(colmod,rowmod))/(GDACmat[ife][step-1]-GDACmat[ife][step]);
			} else if ((*hm2)(colmod,rowmod) == 0 && (*hm1)(colmod,rowmod) != 0) {
			  FinalThr = (*hm1)(colmod,rowmod);
			} else FinalThr = 0;
			hGdac->set(colmod, rowmod, GdacI);
			hThr->set(colmod, rowmod, FinalThr);
		      }
		    }
		  } else if (dir == 1 && (meanmat[ife][1] - meanmat[ife][0]) != 0 ) {
		    GdacI = (int)(GDACmat[ife][0] + (trg - meanmat[ife][0])*(GDACmat[ife][1] - GDACmat[ife][0])/(meanmat[ife][1] - meanmat[ife][0]));
		    if (nloop == 1){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), 1, 0);
		    } else if (nloop == 2){
		      hm1 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, 0, 0, 0);
		      hm2 = &scn->getHisto(PixScan::SCURVE_MEAN, mod, 1, 0, 0);
		    }
		    for (unsigned int col=0; col<NCOLS; col++){
		      for (unsigned int row=0; row<NROWS; row++){
// 			if (ife < 8) {
			  rowmod = row;
			  colmod = col + NCOLS*ife;
// 			} else {
// 			  rowmod = 319 - row;
// 			  colmod = 143 - col - 18*(ife-8);
// 			}
			if ((*hm2)(colmod,rowmod) != 0 && (*hm1)(colmod,rowmod) != 0 ){
			  FinalThr = (*hm1)(colmod,rowmod) + (GdacI-GDACmat[ife][0]) * ((*hm2)(colmod,rowmod)-(*hm1)(colmod,rowmod)) / (GDACmat[ife][1] - GDACmat[ife][0]);
			} else if ((*hm2)(colmod,rowmod) == 0 && (*hm1)(colmod,rowmod) != 0) {
			  FinalThr = (*hm1)(colmod,rowmod);
			} else FinalThr = 0;
			hGdac->set(colmod, rowmod, GdacI);
			hThr->set(colmod, rowmod, FinalThr);
		      }
		    }
		  }
		} else {
		  GdacI = 14;
		}
// 		if (GdacI > 31) GdacI = 31; // Protection
// 		if (GdacI < 6)  GdacI = 6;  // Protection
	      } else {
		for (unsigned int col=0; col<NCOLS; col++){
		  for (unsigned int row=0; row<NROWS; row++){
// 		    if (ife < 8) {
		      rowmod = row;
		      colmod = col + NCOLS*ife;
// 		    } else {
// 		      rowmod = 319 - row;
// 		      colmod = 143 - col - 18*(ife-8);
// 		    }
		    GdacI = 14;
		    FinalThr = 0;
		    hGdac->set(colmod, rowmod, GdacI);
		    hThr->set(colmod, rowmod, FinalThr);
		  }
		}
	      }
	      // Save in the GDAC tuned values in the configuration
	      if(fei4a!=0){
		(*fe)->writeGlobRegister("Vthin_AltFine", GdacI&0xff);
		(*fe)->writeGlobRegister("Vthin_AltCoarse", (GdacI&0xff00)>>8);
		if(fei4b!=0){
		  (*fe)->writeGlobRegister("Vthin_AltFine", GdacI&0xff);
		  (*fe)->writeGlobRegister("Vthin_AltCoarse", (GdacI&0x7f00)>>7);
		}else
		  (*fe)->writeGlobRegister("GLOBAL_DAC", GdacI);
	      }

	    }
	  }
	}
      }
    }
  }
}

void PixModuleGroup::endGDACFastTuning(int nloop, PixScan *scn) {
	if(PMG_DEBUG) cout << "Starting PixModuleGroup::endGDACFastTuning" << endl;
	// Check if the loop is executed on the host
	if (!scn->getDspLoopAction(nloop)) {
		Histo *hn, *hOcc, *hGdac;
		int trg = scn->getRepetitions() * scn->getLoopVarNSteps(nloop-1) / 2;
		if(PMG_DEBUG) cout << "Total target Occ: " << trg << endl;
		// Loop on modules
		for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
			if (m_modules[pmod]->m_readoutActive) {
				unsigned int mod = m_modules[pmod]->m_moduleId;
				if (nloop == 1) {
					hOcc = &scn->getHisto(PixScan::GDAC_OCC, mod, scn->scanIndex(2), 0, 0);
					hGdac = &scn->getHisto(PixScan::GDAC_T, mod, scn->scanIndex(2), 0, 0);
					hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(2), scn->scanIndex(nloop), 0);
				} else if (nloop == 2) {
					hOcc = &scn->getHisto(PixScan::GDAC_OCC, mod, 0, 0, 0);
					hGdac = &scn->getHisto(PixScan::GDAC_T, mod, 0, 0, 0);
					hn = &scn->getHisto(PixScan::SUM_OCC, mod, scn->scanIndex(nloop), scn->scanIndex(nloop), 0);
				}
				// FE loop
				for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
				        PixFe* fei4a = dynamic_cast<PixFeI4A*>(*fe);
					PixFe* fei4b = dynamic_cast<PixFeI4B*>(*fe);
					ConfMask<unsigned short int> &tdacs = (*fe)->readTrim("TDAC"); // get tdacs similar to prepareGDACFastTuning
					int gdac = 0;
					if(fei4a!=0){
					  gdac = (*fe)->readGlobRegister("Vthin_AltFine");
					  gdac += (*fe)->readGlobRegister("Vthin_AltCoarse")<<8;
					}else if(fei4b!=0){
					  gdac = (*fe)->readGlobRegister("Vthin_AltFine");
					  gdac += ((*fe)->readGlobRegister("Vthin_AltCoarse")&0xfe)<<7;
					}else
						gdac = (*fe)->readGlobRegister("GLOBAL_DAC");
					if(PMG_DEBUG) cout << "\t read old gdac value: " << gdac << endl;

					// Calculate mean of Occ
					int hnMean = 0;
					int nrpixels = 0;
					for (unsigned int col=0; col<tdacs.get().size(); col++) {
						for (unsigned int row=0; row<tdacs[0].size(); row++) {
							unsigned int colmod, rowmod;
							int ife = (*fe)->number();
							// 	    if (ife < 8) {
							rowmod = row;
							colmod = col + NCOLS*ife;
							hnMean += (int)(*hn)(colmod, rowmod);
							nrpixels++;
						}
					}
					if (nrpixels != 0) {
						hnMean = (int) hnMean / nrpixels;
						if(PMG_DEBUG) cout << " \t\t MeanOcc = " << hnMean << endl;
					} else {
						if(PMG_DEBUG) cout << "!!! Aborting, nrpixels is 0 !!!" << endl;
						return;
					}

					for (unsigned int col=0; col<tdacs.get().size(); col++) {
						for (unsigned int row=0; row<tdacs[0].size(); row++) {
							unsigned int colmod, rowmod;
							int ife = (*fe)->number();
							// 	      if (ife < 8) {
							rowmod = row;
							colmod = col + NCOLS*ife;

							// Update best mach histos and write gdac
							if (abs(hnMean-trg) < abs((*hOcc)(colmod,rowmod)-trg)) {
								hOcc->set(colmod, rowmod, (*hn)(colmod,rowmod));
								hGdac->set(colmod, rowmod, gdac);
							}
							gdac = (int)(*hGdac)(colmod,rowmod);
						}
					}

					if(fei4a!=0){
					  (*fe)->writeGlobRegister("Vthin_AltFine", gdac&0xff);
					  (*fe)->writeGlobRegister("Vthin_AltCoarse", (gdac&0xff00)>>8);
					}else if(fei4b!=0){
					  (*fe)->writeGlobRegister("Vthin_AltFine", gdac&0xff);
					  (*fe)->writeGlobRegister("Vthin_AltCoarse", (gdac&0x7f00)>>7);
					}else
					  (*fe)->writeGlobRegister("GLOBAL_DAC", gdac);
				}
			}
		}
	}
	if(PMG_DEBUG) cout << "End of PixModuleGroup::endGDACFastTuning" << endl;
}

void PixModuleGroup::endFDACTuning(int nloop, PixScan *scn) {
  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      Histo *hTot, *hFdac, *hTrg;
      hTot = new Histo("FDAC_TOT", "Tuned Fdacs ToT", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      hFdac = new Histo("FDAC_T", "Tuned Fdacs", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      // Associate different parameters if the fcad tuning is done in loop 0, 1 or 2
      if (nloop == 0) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
      } else if (nloop == 1) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, scn->scanIndex(2), -1, -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, scn->scanIndex(2), -1, -1);
      } else if (nloop == 2) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, -1, -1, -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, -1, -1, -1);
      }
      int trg = scn->getTotTargetValue();
      if(PMG_DEBUG) cout << "DEBUG PixModuleGroup: FDAC ToT target is " << trg << endl;
      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	ConfMask<unsigned short int> &fdacs = (*fe)->readTrim("FDAC");
	for (unsigned int col=0; col<fdacs.get().size(); col++) {
	  for (unsigned int row=0; row<fdacs[0].size(); row++) {
	    //if(PMG_DEBUG) cout << "DEBUG PixModuleGroup: FDAC tuning processing col-row" << col << "-" << row << endl;
	    unsigned int colmod, rowmod;
	    int ife = (*fe)->number();
// 	    if (ife < 8) {
	      rowmod = row;
	      colmod = col + NCOLS*ife;
// 	    } else {
// 	      rowmod = 319 - row;
// 	      colmod = 143 - col - 18*(ife-8);
// 	    }
	    for (int index=0; index<scn->getLoopVarNSteps(nloop); index++) { // Loop on FDACS values
	      if (nloop == 0) {
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), index);
	      } else if (nloop == 1) {
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), index, 0);
	      } else if (nloop == 2){
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, index, 0, 0);
	      }
	      if (index == 0) {
		hTot->set(colmod, rowmod, (*hTrg)(colmod, rowmod));
		hFdac->set(colmod, rowmod, index);
	      } else {
		if (abs((*hTrg)(colmod,rowmod)-trg) < abs((*hTot)(colmod,rowmod)-trg)) {
		  hTot->set(colmod, rowmod, (*hTrg)(colmod,rowmod));
		  hFdac->set(colmod, rowmod, index);
		}
	      }
	    }
            int ix = (int)(*hFdac)(colmod,rowmod);
	    if(PMG_DEBUG && col==40 && row>100 && row<120) cout << "DEBUG PixModuleGroup: FDAC tuning best point: " << ix << endl;
	    fdacs[col][row] = (int)(scn->getLoopVarValues(nloop))[ix];
	    hFdac->set(colmod, rowmod, (scn->getLoopVarValues(nloop))[ix]);
	  }
	}
      }
    }
  }
}

void PixModuleGroup::endFDACTuningAlt(int nloop, PixScan *scn) {
  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      Histo *hTot, *hFdac, *hTrg;
      hTot = new Histo("FDAC_TOT", "Tuned Fdacs ToT", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      hFdac = new Histo("FDAC_T", "Tuned Fdacs", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      // Associate different parameters if the fcad tuning is done in loop 0, 1 or 2
      if (nloop == 0) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
      } else if (nloop == 1) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, scn->scanIndex(2), -1, -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, scn->scanIndex(2), -1, -1);
      } else if (nloop == 2) {
	scn->addHisto(*hTot, PixScan::FDAC_TOT, mod, -1, -1, -1);
	scn->addHisto(*hFdac, PixScan::FDAC_T, mod, -1, -1, -1);
      }
      // determine target from central scan point
      int index = scn->getLoopVarNSteps(nloop)/2-1;
      if (nloop == 0) {
	hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), index);
      } else if (nloop == 1) {
	hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), index, 0);
      } else if (nloop == 2){
	hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, index, 0, 0);
      }
      double trg = 0.;
      int nvals=0;
      for (unsigned int col=0; col<(unsigned int)hTrg->nBin(0); col++) {
	for (unsigned int row=0; row<(unsigned int)hTrg->nBin(1); row++) {
	  double currval = (*hTrg)(col,row);
	  if(currval>0){
	    trg += currval;
	    nvals++;
	  }
	}
      }
      if(nvals>0) trg /= (double)nvals;
      else        trg = 0.;
      if(PMG_DEBUG) cout << "DEBUG PixModuleGroup: FDAC ToT target (alt.) from scan pt. " << index << " is " << trg << endl;

      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	ConfMask<unsigned short int> &fdacs = (*fe)->readTrim("FDAC");
	for (unsigned int col=0; col<fdacs.get().size(); col++) {
	  for (unsigned int row=0; row<fdacs[0].size(); row++) {
	    unsigned int colmod, rowmod;
	    int ife = (*fe)->number();
	    rowmod = row;
	    colmod = col + NCOLS*ife;
	    for (int index=0; index<scn->getLoopVarNSteps(nloop); index++) { // Loop on FDACS values
	      if (nloop == 0) {
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), index);
	      } else if (nloop == 1) {
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), index, 0);
	      } else if (nloop == 2){
		hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, index, 0, 0);
	      }
	      if (index == 0) {
		hTot->set(colmod, rowmod, (*hTrg)(colmod, rowmod));
		hFdac->set(colmod, rowmod, index);
	      } else {
		if (abs((*hTrg)(colmod,rowmod)-trg) < abs((*hTot)(colmod,rowmod)-trg)) {
		  hTot->set(colmod, rowmod, (*hTrg)(colmod,rowmod));
		  hFdac->set(colmod, rowmod, index);
		}
	      }
	    }
            int ix = (int)(*hFdac)(colmod,rowmod);
	    if(PMG_DEBUG && col==40 && row>100 && row<120) cout << "DEBUG PixModuleGroup: FDAC tuning best point: " << ix << endl;
	    fdacs[col][row] = (int)(scn->getLoopVarValues(nloop))[ix];
	    hFdac->set(colmod, rowmod, (scn->getLoopVarValues(nloop))[ix]);
	  }
	}
      }
    }
  }
}

void PixModuleGroup::endIFTuning(int nloop, PixScan *scn) {
  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      Histo *hTot, *hIFdac, *hTrg;
      hTot = new Histo("IF_TOT", "Tuned IF ToT", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      hIFdac = new Histo("IF_T", "Tuned IF", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      // Associate different parameters if the fcad tuning is done in loop 0, 1 or 2
      if (nloop == 0) {
	scn->addHisto(*hTot, PixScan::IF_TOT, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	scn->addHisto(*hIFdac, PixScan::IF_T, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
      } else if (nloop == 1) {
	scn->addHisto(*hTot, PixScan::IF_TOT, mod, scn->scanIndex(2), -1, -1);
	scn->addHisto(*hIFdac, PixScan::IF_T, mod, scn->scanIndex(2), -1, -1);
      } else if (nloop == 2) {
	scn->addHisto(*hTot, PixScan::IF_TOT, mod, -1, -1, -1);
	scn->addHisto(*hIFdac, PixScan::IF_T, mod, -1, -1, -1);
      }
      int trg = scn->getTotTargetValue();
      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++) {
	int ife = (*fe)->number();
        int bestIF;
        double bestMean;
	for (int index=0; index<scn->getLoopVarNSteps(nloop); index++) { // Loop on FDACS values
	  if (nloop == 0) {
	    hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), index);
	  } else if (nloop == 1) {
	    hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), index, 0);
	  } else if (nloop == 2){
	    hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, index, 0, 0);
	  }
          int count = 0;
          double mean = 0;
	  for (unsigned int col=0; col<NCOLS; col++) {
	    for (unsigned int row=0; row<NROWS; row++) {
	      unsigned int colmod, rowmod;
// 	      if (ife < 8) {
		rowmod = row;
		colmod = col + NCOLS*ife;
// 	      } else {
// 		rowmod = 319 - row;
// 		colmod = 143 - col - 18*(ife-8);
// 	      }
	      if (((*hTrg)(colmod, rowmod)) > 0) {
		mean += (*hTrg)(colmod, rowmod); count++;
	      }
	    }
	  }
	  mean /= count;
	  if (index == 0) {
	    bestMean = mean;
	    bestIF = index;
	  } else {
	    if ( abs(mean-trg) < abs(bestMean-trg) ) {
	      bestMean = mean;
	      bestIF = index;
	    }
	  }
	}
	if (nloop == 0) {
	  hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), scn->scanIndex(1), bestIF);
	} else if (nloop == 1) {
	  hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, scn->scanIndex(2), bestIF, 0);
	} else if (nloop == 2){
	  hTrg = &scn->getHisto(PixScan::TOT_MEAN, mod, bestIF, 0, 0);
	}
	for (unsigned int col=0; col<NCOLS; col++) {
	  for (unsigned int row=0; row<NROWS; row++) {
	    unsigned int colmod, rowmod;
// 	    if (ife < 8) {
	      rowmod = row;
	      colmod = col + NCOLS*ife;
// 	    } else {
// 	      rowmod = 319 - row;
// 	      colmod = 143 - col - 18*(ife-8);
// 	    }
	    hIFdac->set(colmod, rowmod, (scn->getLoopVarValues(nloop))[bestIF]);
	    hTot->set(colmod, rowmod, (*hTrg)(colmod, rowmod));
	  }
	}
	//	(*fe)->writeGlobRegister("DAC_IF",(int)(scn->getLoopVarValues(nloop))[bestIF]);
	(*fe)->writeGlobRegister("PrmpVbpf",(int)(scn->getLoopVarValues(nloop))[bestIF]);
      }
    }
  }
}

void PixModuleGroup::endT0Set(int nloop, PixScan *scn) {
  const int min_nok = 200;
  // We assume loop0 is on strobe delay and loop1 on trigger delay
  if (nloop != 1) return; //throw an exception?

  std::vector<float> mean_del[32][NFES]; // first index: module, second index: FE
  std::vector<int> npix_ok, n_del[32][NFES];
  int il;
  for (il=0; il<scn->getLoopVarNSteps(nloop); il++) {
    npix_ok.push_back(0);
  }
  Histo *hTime;
  // Loop on modules
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    // if the module is active
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      // Loop over LVL1 slices
      for (il=0; il<scn->getLoopVarNSteps(nloop); il++) {
	hTime  = &scn->getHisto(PixScan::TIMEWALK, mod, scn->scanIndex(2), il, 0);
	unsigned int colmod, rowmod;
	for(unsigned int iFE=0;iFE<NFES;iFE++){
	  mean_del[mod][iFE].push_back(0);
	  n_del[mod][iFE].push_back(0);
	  for (unsigned int col=0; col<NCOLS; col++) {
	    for (unsigned int row=0; row<NROWS; row++) {
	      rowmod = row;
	      colmod = col + NCOLS*iFE;
	      double a = (*hTime)(colmod, rowmod);
	      if(a>0){
		mean_del[mod][iFE][il] += a;
		(n_del[mod][iFE][il])++;
		npix_ok[il]++;
	      }
	    }
	  }
	}
      }
      for (il=0; il<scn->getLoopVarNSteps(nloop); il++) {
	for(unsigned int iFE=0;iFE<NFES;iFE++){
	  if(n_del[mod][iFE][il]>min_nok){
	    mean_del[mod][iFE][il]/=(float)n_del[mod][iFE][il];
	    if(PMG_DEBUG) cout << "PixModuleGroup::endT0Set: module " << mod << " FE " << iFE << " : mean delay for step " << il
			       << " is OK and " << mean_del[mod][iFE][il] << endl;
	  }else
	    mean_del[mod][iFE][il] = 0;
	}
      }
    }
  }
  // Determination of optimal strobe delays
  int max_ok = 0, il_ok = -1;
  for (il=0; il<scn->getLoopVarNSteps(nloop); il++) {
    if (npix_ok[il] > max_ok) {
      max_ok = npix_ok[il];
      il_ok = il;
    }
  }
  if (il_ok >= 0) {
    m_triggerDelay = (int)((scn->getLoopVarValues(nloop))[il_ok]);
    if(PMG_DEBUG) cout << "PixModuleGroup::endT0Set: setting grp. trg. delay to " << m_triggerDelay << endl;
    int il_ok2=-1;
    if(il_ok > 0 && npix_ok[il_ok-1]>min_nok &&
       !(il_ok<(scn->getLoopVarNSteps(nloop)-1) && npix_ok[il_ok-1]<npix_ok[il_ok+1]))
      il_ok2 = il_ok-1;
    if(il_ok<(scn->getLoopVarNSteps(nloop)-1)  && npix_ok[il_ok+1]>min_nok &&
       !(il_ok>0 && npix_ok[il_ok+1]<npix_ok[il_ok-1]))
      il_ok2 = il_ok+1;
    for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
      // if the module is active
      if (m_modules[pmod]->m_readoutActive) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	hTime  = &scn->getHisto(PixScan::TIMEWALK, mod, scn->scanIndex(2), il, 0);
	for(unsigned int iFE=0;iFE<NFES;iFE++){
	  PixFe* fei4 = dynamic_cast<PixFeI4A*>(m_modules[pmod]->pixFE(iFE));
	  if(fei4==0) fei4 = dynamic_cast<PixFeI4B*>(m_modules[pmod]->pixFE(iFE));
	  // select best mean value above or below chosen trigger delay point
	  double oldconv = 1., conv = 1.;
	  if(il_ok2>=0 && fabs(mean_del[mod][iFE][il_ok2]-mean_del[mod][iFE][il_ok])>0.)
	    conv = 25./(fabs(mean_del[mod][iFE][il_ok2]-mean_del[mod][iFE][il_ok]));
	  if(PMG_DEBUG) cout << "PixModuleGroup::endT0Set: conversion input: " << ((il_ok2>=0)?mean_del[mod][iFE][il_ok2]:-1.)
			     << " -> conv = " << conv << endl;
	  if(fei4!=0){
	    oldconv = fei4->getDelayCalib();
	    fei4->setDelayCalib(conv * oldconv);
	  }
	  // re-calibrate histograms
	  unsigned int colmod, rowmod;
	  for (il=0; il<scn->getLoopVarNSteps(nloop); il++) {
	    for (unsigned int col=0; col<NCOLS; col++) {
	      for (unsigned int row=0; row<NROWS; row++) {
		rowmod = row;
		colmod = col + NCOLS*iFE;
		double a = (*hTime)(colmod, rowmod);
		hTime->set(colmod,rowmod,a*conv);
	      }
	    }
	  }
	  // subtact 5ns from mean
	  if(mean_del[mod][iFE][il_ok]>0.) mean_del[mod][iFE][il_ok]=(mean_del[mod][iFE][il_ok]+5./conv)/oldconv+0.5;
	  if(mean_del[mod][iFE][il_ok]>0. && mean_del[mod][iFE][il_ok]<64.){
	    if(PMG_DEBUG) cout << "PixModuleGroup::endT0Set: setting FE strb. delay to "
			       << ((int)(mean_del[mod][iFE][il_ok])) << " in mod. " << mod << ", FE " << iFE << endl;
	    if(fei4!=0) fei4->writeGlobRegister("PlsrDelay", (int)(mean_del[mod][iFE][il_ok]));
	  } else
	    cerr << "PixModuleGroup::endT0Set: FE strb. delay is out of range: "
		 << ((int)(mean_del[mod][iFE][il_ok])) << " in mod. " << mod << ", FE " << iFE << endl;
	}
      }
    }
  }
}

void PixModuleGroup::mccDelFit(int nloop, PixScan *scn) {
  // check if this action makes sense at all
  if(!scn->getHistogramFilled(PixScan::TIMEWALK) || !scn->getHistogramFilled(PixScan::OCCUPANCY)) return;
  // Loop on modules
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    // if the module is active
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      // create the histos
      Histo *hTime, *hSigma, *hChi;
      // Conversion from MCC range to ns
      std::stringstream a;
      a << scn->getStrobeMCCDelayRange();
      float conv = 1;
      if(m_modules[pmod]->pixMCC()!=0)
	conv = ((ConfFloat&)m_modules[pmod]->pixMCC()->config()["Strobe"]["DELAY_"+a.str()]).m_value;
      // fit MCC delay
      hTime = new Histo("TIMEWALK", "Timewalk", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      hSigma = new Histo("SIGMA", "Sigma", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
      hChi = new Histo("CHI2", "Chi2", NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
//       if(!scn->getHistogramKept(PixScan::OCCUPANCY)) // temporarily retrieve histos
// 	scn->downloadHisto(m_pixCtrl, mod, PixScan::OCCUPANCY);
      PixScanHisto &sh = scn->getHisto(PixScan::OCCUPANCY);
      PixScanHisto *sc = 0;
      if (nloop == 0) {
	if(scn->getHistogramKept(PixScan::TIMEWALK))
	  scn->addHisto(*hTime, PixScan::TIMEWALK, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
	sc = &(sh[mod][scn->scanIndex(2)][scn->scanIndex(1)]);
      }else if (nloop == 1) {
	if(scn->getHistogramKept(PixScan::TIMEWALK))
	  scn->addHisto(*hTime, PixScan::TIMEWALK, mod, scn->scanIndex(2), -1, -1);
	sc = &(sh[mod][scn->scanIndex(2)]);
      }else if (nloop == 2){
	if(scn->getHistogramKept(PixScan::TIMEWALK))
	  scn->addHisto(*hTime, PixScan::TIMEWALK, mod, -1, -1, -1);
	sc = &(sh[mod]);
      }
      // S-curve fit
      scn->fitSCurve(*sc, *hTime, *hSigma, *hChi, 0, scn->getRepetitions());
      // clear occupancy histos if requested by user
//       if(!scn->getHistogramKept(PixScan::OCCUPANCY)){ //sc->clear();
// 	// just clear histograms and leave scan entry intact
// 	// otherwise will fail at higher loops...
// 	for(unsigned int isc=1;isc<sc->size();isc++) ((*sc)[isc]).clear();
//       }
      // Convert all histograms (MCC-del. -> ns)
      for (unsigned int col=0; col<NCOLM; col++) {
	for (unsigned int row=0; row<NROWM; row++) {
	  double val = (*hTime)(col, row)*(double)conv;
	  hTime->set(col, row, val);
	}
      }
      delete hSigma;
      delete hChi;
    }
  }
}

void PixModuleGroup::fitCalib(int nloop, PixScan *scn){
  // check if this action makes sense at all
  if(!scn->getHistogramFilled(PixScan::SCURVE_MEAN)) return;
  // Loop on modules
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    // if the module is active
    if (m_modules[pmod]->m_readoutActive) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      Histo *hm;
      // Loop on FEs
      for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	int ife = (*fe)->number();
	int nreg=0;
	float sx=0., sy=0., sx2=0., sxy=0.;
	// Loop over steps (finding the mean threshold in each FE)
	for (int step = 0; step<scn->getLoopVarNSteps(nloop); step++){
	  // Check if cap. scan was done in loop 1 or 2
	  if (nloop == 1){
	    hm      = &scn->getHisto(PixScan::SCURVE_MEAN, mod, scn->scanIndex(2), step, 0);
	  }else if (nloop == 2){
	    hm      = &scn->getHisto(PixScan::SCURVE_MEAN, mod, step, 0, 0);
	  }
	  unsigned int colmod, rowmod;
	  double mean = 0, meanlocal = 0;
	  int npixel=0;
	  // Cols and rows are the pixel coordinates in FE level
	  for (unsigned int col=0; col<NCOLS; col++){
	    for (unsigned int row=0; row<NROWS; row++){
	      rowmod = row;
	      colmod = col + NCOLS*ife;
	      meanlocal = (int)(*hm)(colmod,rowmod); // Mean value in each pixel
	      if (meanlocal > 100) {
		npixel++;
		mean = mean + meanlocal;
	      }
	    }
	  }
	  if (npixel != 0) mean = mean/npixel;
	  std::string capLabels[3]={"CInjLo", "CInjMed", "CInjHi"};
	  int scanVal = (int) scn->getLoopVarValues(nloop)[step];
	  if(scanVal>=0 || scanVal<3){
	    float cInj     = ((dynamic_cast<ConfFloat &>((*fe)->config()["Misc"][capLabels[scanVal]])).value())/0.160218;
	    // to do: check avg. threshold vs capacity
	    if(PMG_DEBUG) cout << "PixModuleGroup::fitCalib : FE " << ife << " has avg. thresh. of " << mean << " at cap. of " << cInj << endl;
	    // fill linear regression variables
	    sx  += cInj;
	    sy  += mean;
	    sx2 += cInj*cInj;
	    sxy += cInj*mean;
	    nreg++;
	  }
	}
	// get (-1)*slope from straight line and fill into FE calib
	float VCALoffset = 0.;
	if((sx2*(float)nreg-sx*sx)!=0.) VCALoffset = (sx*sy - sxy*(float)nreg)/(sx2*(float)nreg-sx*sx);
	ConfFloat &offsCfg = dynamic_cast<ConfFloat &>((*fe)->config()["Misc"]["VcalGradient0"]);
	VCALoffset += offsCfg.value();
	if(PMG_DEBUG) cout << "PixModuleGroup::fitCalib : FE " << ife << " has VCAL offset of " << VCALoffset << endl;
	offsCfg.m_value = VCALoffset;
      }
    }
  }
}

void PixModuleGroup::endIncrTdac(int nloop, PixScan *scn) {
  // Module loop
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
    if (m_modules[pmod]->m_readoutActive) {
      m_modules[pmod]->deleteConfig("IncrTdacPreviousConfig");
    }
  }
}

void PixModuleGroup::initScan(PixScan *scn){

  //This method has to prepare the module group for the scan.

  //m_dcs->setLVOn(); //switch on voltages
  //m_dcs->setHVOn();

  //check controller status;
  m_pixCtrl->testHW();

  // Save module configurations if restore is needed at the end of the scan
  if (scn->getRestoreModuleConfig()) {
    for (unsigned int pmod=0; pmod<m_modules.size(); pmod++){
      for (int gr=0; gr<MAX_GROUPS;gr++){
	if ((scn->getModuleMask(gr)) & ((0x1)<<(m_modules[pmod]->m_moduleId))) {
	  m_modules[pmod]->storeConfig("PreScanConfig");
	}
      }
    }
  }

  if (scn->getRunType() != PixScan::RAW_PATTERN && scn->getRunType() != PixScan::RAW_EVENT && scn->getConfigAtStart()){
    //configure modules - only for non-BOC scans
    m_pixCtrl->setConfigurationMode(); //preparo il rod a compiere l`azione di caricamento dei file di config dei moduli
    m_pixCtrl->sendModuleConfig(0x0); // manda a tutti i moduli materialmente
  }

  //before scanning put all modules properties false
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++){
    m_modules[pmod]->m_configActive = false;
    m_modules[pmod]->m_triggerActive = false;
    m_modules[pmod]->m_strobeActive = false;
    m_modules[pmod]->m_readoutActive = false;
    //if the module in the configuration is the same in the group that I`m going to scan, switch on its properties
    for (int gr=0; gr<MAX_GROUPS;gr++){
      if ((scn->getModuleMask(gr)) & ((0x1)<<(m_modules[pmod]->m_moduleId))) {
	if (scn->getConfigEnabled(gr)) m_modules[pmod]->m_configActive = true;
	if (scn->getTriggerEnabled(gr)) m_modules[pmod]->m_triggerActive = true;
	if (scn->getStrobeEnabled(gr)) m_modules[pmod]->m_strobeActive = true;
	if (scn->getReadoutEnabled(gr)) m_modules[pmod]->m_readoutActive = true;
      }
    }
  }
  if (scn->getRunType() == PixScan::RAW_PATTERN || scn->getRunType() == PixScan::RAW_EVENT) {
    RodPixController *rod = dynamic_cast<RodPixController*>(m_pixCtrl);
    if (rod) {
      // Raw pattern or event
      rod->prepareTestPattern(scn);
    }
    //JW: Sollte sicher sein. Es ist rod=NULL wenn m_pixCtrl!=RodPixController.
  }
}

void PixModuleGroup::scanLoopStart(int nloop, PixScan *scn) {
  // This method prepares the module group for the beginning of a particular loop.
  // The typical operation to perform  in this phase is the setting of the loop scan
  // variable to the appropriate initial value.

  // Backward comaptibility for apps without scanTerminate
  if (nloop == 0) m_execToTerminate = false;
  // Do nothing else if loop is disabled or handled by the ROD
  if (scn->getLoopActive(nloop)) {
    switch (scn->getLoopAction(nloop)) {
    case PixScan::TOTCAL_FIT:
      break;
    case PixScan::SCURVE_FIT:
      break;
    case PixScan::TDAC_TUNING:
      break;
    case PixScan::TDAC_FAST_TUNING:
      break;
	case PixScan::GDAC_FAST_TUNING:
      break;
    case PixScan::GDAC_TUNING:
      break;
    case PixScan::FDAC_TUNING:
    case PixScan::FDAC_TUNING_ALT:
      prepareFDACTuning(nloop, scn);
      break;
    case PixScan::IF_TUNING:
      prepareIFTuning(nloop, scn);
      break;
    case PixScan::T0_SET:
      break;
    case PixScan::MCCDEL_FIT:
      break;
    case PixScan::OFFSET_CALIB:
      break;
    case PixScan::NO_ACTION:
    default:
      break;
    }
  }
}

void PixModuleGroup::prepareStep(int nloop, PixScan *scn) {
  // Backward comaptibility for apps without scanTerminate
  if (nloop == 0 && m_execToTerminate) scanTerminate(scn);
  // If you are in the innermost loop & if mask staging is under host control
  if (nloop == 0 && !scn->getDspMaskStaging() && scn->newMaskStep()) setupMasks(0, scn); //setup masks
  // Do nothing else if loop is disabled or handled by the ROD
  if (scn->getLoopActive(nloop) && !scn->getDspProcessing(nloop)) {
    if(!scn->getLoopVarValuesFree(nloop)){ // don't try to change scan variable if there isn't any in this loop
      if (nloop == 0) {
	if (scn->newScanStep()) setupScanVariable(0, scn); //loop on variable
	// If you are in the other loops
      } else {
	setupScanVariable(nloop, scn);
      }
    }

    // Loop action specific settings
    if(PMG_DEBUG) cout << "PixModuleGroup::prepareStep: loop " << nloop << " action is " << scn->getLoopAction(nloop) << " - NO_ACTION: " << PixScan::NO_ACTION << endl;
    switch (scn->getLoopAction(nloop)) {
    case PixScan::TDAC_TUNING:
      prepareTDACTuning(nloop, scn);
      break;
    case PixScan::GDAC_TUNING:
      //      prepareGDACTuning(nloop, scn);
      break;
    case PixScan::FDAC_TUNING:
    case PixScan::FDAC_TUNING_ALT:
      break;
    case PixScan::IF_TUNING:
      break;
    case PixScan::T0_SET:
      prepareT0Set(nloop, scn);
      break;
    case PixScan::MIN_THRESHOLD:
      prepareIncrTdac(nloop, scn);
      break;
    case PixScan::MCCDEL_FIT:
      break;
	case PixScan::TDAC_FAST_TUNING:
		prepareTDACFastTuning(nloop, scn);
		break;
	case PixScan::OCC_SUM:
		//cout << "call prepareOccSumming" << endl;
		//prepareOccSumming(nloop, scn);
		break;
	case PixScan::GDAC_FAST_TUNING:
		prepareGDACFastTuning(nloop, scn);
		break;
    case PixScan::NO_ACTION:
    default:
      break;
    }
  }
}

void PixModuleGroup::scanExecute(PixScan *scn) {

  RodPixController *rod = dynamic_cast<RodPixController*>(m_pixCtrl);
  USBPixController *usb = dynamic_cast<USBPixController*>(m_pixCtrl);

  // Backward comaptibility for apps without scanTerminate
  m_execToTerminate = true;
  // Check if the controller is a ROD
  if (rod) {
    if (scn->getRunType() == PixScan::NORMAL_SCAN) {
      // Temp; configure modules
      m_pixCtrl->setConfigurationMode();
      m_pixCtrl->sendModuleConfig(0x0);
      // Enable trapping
      for (int i=0; i<4; i++) {
	rod->setupTrapping(i);
      }
      rod->writeScanConfig(*scn);
      std::string txt;
      rod->getDiagBuffer(txt);
      rod->setCalibrationMode();
      rod->startScan();
    } else {
      // Raw pattern or event
      std::vector< Histo* > vh;
      rod->runTestPattern(scn, vh);
      for (int mod=0; mod<32; mod++) {
	if (vh[mod] != NULL) {
	  scn->addHisto(*vh[mod], PixScan::RAW_DATA_0, mod, scn->scanIndex(2), scn->scanIndex(1),scn->scanIndex(0));
	}
      }
    }
  }
  else if (usb) {
    if(scn->getRunType() == PixScan::NORMAL_SCAN) {
      // configure modules
      m_pixCtrl->setConfigurationMode();
      //m_pixCtrl->sendModuleConfig(0x0);
      // for FE-I4, be more specific and only send what is needed
      for(int iloop=0;iloop<3;iloop++){
	if(scn->getLoopActive(iloop)){
	  switch(scn->getLoopParam(iloop)){
	  case PixScan::VCAL:
	    m_pixCtrl->sendGlobal(0, "PlsrDAC");
	    break;
	  case PixScan::GDAC:
	    m_pixCtrl->sendGlobal(0, "Vthin_AltCoarse");
	    m_pixCtrl->sendGlobal(0, "Vthin_AltFine");
	    break;
	  case PixScan::IF:
	    m_pixCtrl->sendGlobal(0, "PrmpVbpf");
	    break;
	  case PixScan::FEI4_GR:
	    m_pixCtrl->sendGlobal(0, scn->getLoopFEI4GR(iloop));
	    break;
	  case PixScan::STROBE_DELAY:
	    m_pixCtrl->sendGlobal(0, "PlsrDelay");
	    break;
	  case PixScan::TDACS:
	    m_pixCtrl->sendPixel(0, "TDAC4");
	    m_pixCtrl->sendPixel(0, "TDAC3");
	    m_pixCtrl->sendPixel(0, "TDAC2");
	    m_pixCtrl->sendPixel(0, "TDAC1");
	    m_pixCtrl->sendPixel(0, "TDAC0");
	    break;
	  case PixScan::FDACS:
	    m_pixCtrl->sendPixel(0, "FDAC3");
	    m_pixCtrl->sendPixel(0, "FDAC2");
	    m_pixCtrl->sendPixel(0, "FDAC1");
	    m_pixCtrl->sendPixel(0, "FDAC0");
	    break;
	  default:
	    break; // scan var. is not on the FE
	  }
	}
      }
      // read DCS if requested and fill histogram
      if(scn->getHistogramFilled(PixScan::DCS_DATA) && scn->getHistogramKept(PixScan::DCS_DATA)){
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  std::ostringstream mnum;
	  mnum << mod;
	  if(scn->getMaskStageIndex()==0 && scn->scanIndex(0)==0){
	    if(PMG_DEBUG) cout << "PixModuleGroup::scanExecute : creating DCS histo at scan pts " << scn->scanIndex(2) << " - " << scn->scanIndex(1)
			       << " - " << scn->scanIndex(0) << " and mask step " << scn->getMaskStageIndex() << endl;
	    // prepare histograms
	    std::string nam, tit;
	    nam = "DCS_DATA_" + mnum.str();
	    tit = "Readings from DCS channel "+scn->getDcsChan()+" module "+mnum.str();
	    Histo *hdcs = new Histo(nam,tit, scn->getLoopVarNSteps(0)*scn->getMaskStageSteps(), 0., (double)(scn->getLoopVarNSteps(0)*scn->getMaskStageSteps()-1));
	    scn->addHisto(*hdcs, PixScan::DCS_DATA, mod, scn->scanIndex(2), scn->scanIndex(1),-1);
	  }
 	  Histo *hdcs = &(scn->getHisto(PixScan::DCS_DATA, mod, scn->scanIndex(2), scn->scanIndex(1),-1));
 	  double dcsval=0.;
 	  std::string rtypes = "unknown";
	  switch(scn->getDcsMode()){
	  case PixScan::VOLTAGE:
	    rtypes = "voltage";
	    break;
	  case PixScan::CURRENT:
	    rtypes = "current";
	    break;
	  }
	  if(m_dcsChans[pmod]!=0 && rtypes != "unknown") dcsval = m_dcsChans[pmod]->ReadParam(rtypes);
	  if(PMG_DEBUG) cout << "PixModuleGroup::scanExecute : DCS chan: " << ((m_dcsChans[pmod]!=0)?m_dcsChans[pmod]->name():"NULL") <<
	    " reads type " << rtypes << " of value " << dcsval << endl;
	  hdcs->set(scn->scanIndex(0)*scn->getMaskStageSteps()+scn->getMaskStageIndex(), dcsval);
	}
      }
      // read source rates if requested and fill histogram
      if((scn->getHistogramFilled(PixScan::HIT_RATE) && scn->getHistogramKept(PixScan::HIT_RATE)) ||
	 (scn->getHistogramFilled(PixScan::TRG_RATE) && scn->getHistogramKept(PixScan::TRG_RATE))){
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  std::ostringstream mnum;
	  mnum << mod;
	  if(scn->getMaskStageIndex()==0 && scn->scanIndex(0)==0){
	    // prepare histograms
	    std::string nam, tit;
	    if(scn->getHistogramFilled(PixScan::HIT_RATE) && scn->getHistogramKept(PixScan::HIT_RATE)){
	      nam = "HIT_RATE_" + mnum.str();
	      tit = "Src hit rate module "+mnum.str();
	      Histo *hhrate = new Histo(nam,tit, scn->getLoopVarNSteps(0)*scn->getMaskStageSteps(), 0., (double)(scn->getLoopVarNSteps(0)*scn->getMaskStageSteps()-1));
	      scn->addHisto(*hhrate, PixScan::HIT_RATE, mod, scn->scanIndex(2), scn->scanIndex(1),-1);
	    }
	    if(scn->getHistogramFilled(PixScan::TRG_RATE) && scn->getHistogramKept(PixScan::TRG_RATE)){
	      nam = "TRG_RATE_" + mnum.str();
	      tit = "Src trigger rate module "+mnum.str();
	      Histo *hhrate = new Histo(nam,tit, scn->getLoopVarNSteps(0)*scn->getMaskStageSteps(), 0., (double)(scn->getLoopVarNSteps(0)*scn->getMaskStageSteps()-1));
	      scn->addHisto(*hhrate, PixScan::TRG_RATE, mod, scn->scanIndex(2), scn->scanIndex(1),-1);
	    }
	  }
	  Histo *hrate = 0;
	  double trval, hrval;
	  const int nmeas=100;
	  if(PMG_DEBUG) cout << "PixModuleGroup : starting src meas." << endl;
	  usb->startScan(scn);
	  hrval = 0.;
	  trval = 0.;
	  if(PMG_DEBUG) cout << "PixModuleGroup : starting rate meas." << endl;
	  for(int i=0;i<nmeas;){
	    sleep(1200); // must wait > 1s, otherwise FPGA isn't updated
	    usb->nTrigger(); // update readings
	    if(usb->getSRAMFillLevel()<85) { // only read when SRAM not close to full, otherwise rates might be wrong
	      hrval += (double)usb->getHitRate();
	      trval += (double)usb->getTriggerRate();
	      i++;
	    }
	  }
	  if(PMG_DEBUG) cout << "PixModuleGroup : stopping src meas." << endl;
	  usb->stopScan();
	  if(PMG_DEBUG) cout << "PixModuleGroup : waiting for meas. stop" << endl;
	  int iwait=0;
	  usb->nTrigger();
	  while(usb->runStatus() && iwait<1000){
	    sleep(100);
	    iwait++;
	    usb->nTrigger();
	  }
	  if(PMG_DEBUG) cout << "PixModuleGroup : storing rate meas., waited " << iwait << endl;
	  hrval /= (double)nmeas;
	  trval /= (double)nmeas;
	  if(scn->getHistogramFilled(PixScan::HIT_RATE) && scn->getHistogramKept(PixScan::HIT_RATE)){
	    hrate = &(scn->getHisto(PixScan::HIT_RATE, mod, scn->scanIndex(2), scn->scanIndex(1),-1));
	    hrate->set(scn->scanIndex(0)*scn->getMaskStageSteps()+scn->getMaskStageIndex(), hrval);
	  }
	  if(scn->getHistogramFilled(PixScan::TRG_RATE) && scn->getHistogramKept(PixScan::TRG_RATE)){
	    hrate = &(scn->getHisto(PixScan::TRG_RATE, mod, scn->scanIndex(2), scn->scanIndex(1),-1));
	    hrate->set(scn->scanIndex(0)*scn->getMaskStageSteps()+scn->getMaskStageIndex(), trval);
	  }
	}
      }
      // process remaining scan part
      // check if histograms from controller are requested
      bool histoOnCtrl=false;
      std::map<std::string, int> htv = scn->getHistoTypes();
      for(std::map<std::string, int>::iterator IT = htv.begin(); IT!=htv.end();IT++){
	if(PixLib::PixScan::DCS_DATA != (PixLib::PixScan::HistogramType)IT->second &&
	   PixLib::PixScan::HIT_RATE != (PixLib::PixScan::HistogramType)IT->second &&
	   PixLib::PixScan::TRG_RATE != (PixLib::PixScan::HistogramType)IT->second)
	  histoOnCtrl |= scn->getHistogramFilled((PixLib::PixScan::HistogramType)IT->second);
      }
      if(scn->getDspMaskStaging() && histoOnCtrl){
	if(PMG_DEBUG) cout << "PixModuleGroup::scanExecute : starting PixCOntroller::startScan" << endl;
	usb->startScan(scn);
      }else{
	if(scn->getMaskStageIndex()==0 && scn->scanIndex(0)==0){
	  // prepare histograms
	  Histo *ml;
	  std::string nam, tit;
	  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	    unsigned int mod = m_modules[pmod]->m_moduleId;
	    std::ostringstream mnum;
	    mnum << mod;
	    if (scn->getHistogramFilled(PixScan::MON_LEAK) && scn->getHistogramKept(PixScan::MON_LEAK)){
	      if(PMG_DEBUG) cout << "creating ML histo for module " << mod << endl;
	      nam = "ML_" + mnum.str();
	      tit = "MonLeak mod " + mnum.str();
	      ml = new Histo(nam, tit, NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	      scn->addHisto(*ml, PixScan::MON_LEAK, mod, scn->scanIndex(2), scn->scanIndex(1),scn->scanIndex(0));
	    }
	    if (scn->getHistogramFilled(PixScan::HB_SCALER_CNT) && scn->getHistogramKept(PixScan::HB_SCALER_CNT)){
	      if(PMG_DEBUG) cout << "creating HB histo for module " << mod << endl;
	      nam = "HBS_" + mnum.str();
	      tit = "HitBus status mod " + mnum.str();
	      ml = new Histo(nam, tit, NCOLM, -0.5, NCOLM-0.5, NROWM, -0.5, NROWM-0.5);
	      scn->addHisto(*ml, PixScan::HB_SCALER_CNT, mod, scn->scanIndex(2), scn->scanIndex(1),scn->scanIndex(0));
	    }
	  }
	}
	// fill histograms pixel by pixel
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  for(std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++){
	    int ife = (*fe)->number();
	    if (scn->getHistogramFilled(PixScan::MON_LEAK) && scn->getHistogramKept(PixScan::MON_LEAK)){
	      double mlval = 0.;
	      std::string rtypes = "unknown";
	      switch(scn->getDcsMode()){
	      case PixScan::VOLTAGE:
		rtypes = "voltage";
		break;
	      case PixScan::CURRENT:
		rtypes = "current";
		break;
	      }
	      if(m_dcsChans[pmod]!=0 && rtypes != "unknown") mlval = m_dcsChans[pmod]->ReadParam(rtypes);
	      if(PMG_DEBUG) cout << "DCS chan: " << ((m_dcsChans[pmod]!=0)?m_dcsChans[pmod]->name():"NULL") << " reads type " << rtypes << " of value " << mlval << endl;
	      Histo *hml = &(scn->getHisto(PixScan::MON_LEAK, mod, scn->scanIndex(2), scn->scanIndex(1),0));
	      fillHistosPerDc(*scn, ife, hml, mlval);
	    }
	    if(scn->getHistogramFilled(PixScan::HB_SCALER_CNT) && scn->getHistogramKept(PixScan::HB_SCALER_CNT)){
	      double hocc = (double)usb->readHitBusScaler(mod, ife,scn);
	      Histo *hho = &(scn->getHisto(PixScan::HB_SCALER_CNT, mod, scn->scanIndex(2), scn->scanIndex(1),0));
	      fillHistosPerDc(*scn, ife, hho, hocc);
 	    }
	  }
	}
      }
    } else { // sync scan
      if(PMG_DEBUG) cout << "PixModuleGroup::scanExecute : starting sync scan" << endl;
      // configure global registers on FE, otherwise scan will not work
      usb->sendGlobal(0);
      // run scan and store result in histo
      Histo *h = usb->syncScan();
      scn->addHisto(*h, PixScan::RAW_DATA_0, 0, scn->scanIndex(2), scn->scanIndex(1), scn->scanIndex(0));
    }
  }
}

void PixModuleGroup::scanTerminate(PixScan *scn) {

  RodPixController *rod = dynamic_cast<RodPixController*>(m_pixCtrl);
  USBPixController *usb = dynamic_cast<USBPixController*>(m_pixCtrl);

  // Backward comaptibility for apps without scanTerminate
  m_execToTerminate = false;
  // Check if the controller is a ROD
  if (rod) {
    // Download dsp error histograms
    if (scn->getHistogramFilled(PixScan::DSP_ERRORS) && scn->getHistogramKept(PixScan::DSP_ERRORS)) {
      for (int dsp=0; dsp<4; dsp++) {
	scn->downloadErrorHisto(m_pixCtrl, dsp);
      }
    }
    // Stop histogramming and trapping tasks
    rod->setConfigurationMode();
  }
  else if (usb) {
    if (scn->getHistogramFilled(PixScan::DSP_ERRORS) && scn->getHistogramKept(PixScan::DSP_ERRORS)) {
      scn->downloadErrorHisto(m_pixCtrl, 0);
    }
  }

  // Download histograms with the FILLED and KEPT flags (for loop 0 if executed by the HOST)
  // or if not kept but needed for end-of-loop action later (to be deleted after action in scanLoopEnd)
  if (scn->getLoopActive(0) && !scn->getDspProcessing(0)) {
    std::map<std::string, int> &hTypes = scn->getDspHistoTypes();
    std::map<std::string, int>::iterator it;
    for (it = hTypes.begin(); it != hTypes.end(); ++it)  {
      std::string name = (*it).first;
      PixScan::HistogramType type = (PixScan::HistogramType)((*it).second);
      if (scn->getHistogramFilled(type) && (scn->getHistogramKept(type) || scn->getLoopAction(0)!=PixScan::NO_ACTION)) {
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  if(PMG_DEBUG) cout << "PixModuleGroup::scanTerminate : downloading histo " << type << endl;
	  scn->downloadHisto(m_pixCtrl, mod, type);
	}
      }
    }
  }
}

void PixModuleGroup::scanLoopEnd(int nloop, PixScan *scn) {

  // This method will perform the end-of-loop actions.
  // Typically this method will upload histograms or fit results from the
  // ROD and store them in the appropriate PixScan structures.

   USBPixController *usb = dynamic_cast<USBPixController*>(m_pixCtrl);

  // Backward comaptibility for apps without scanTerminate
  if (nloop == 0 && m_execToTerminate) scanTerminate(scn);

  // Nothing to do if the loop is inactive, unless it's loop 0
  if (!scn->getLoopActive(nloop) && nloop != 0) return;

  // Nothing to do for loop 0 if loop 1 is executed by the ROD
  if (nloop == 0 && scn->getDspProcessing(1)) return;

  // Download histograms with the FILLED and KEPT flags (for loop 0 and 1 if executed by the ROD)
  // or if not kept but needed for end-of-loop action later (to be deleted after action below)
  if ((nloop == 0 && scn->getDspProcessing(0)) ||
      (nloop == 1 && scn->getDspProcessing(1)) ||
      (nloop == 0 && !scn->getLoopActive(0))) {
    // Download dsp histograms
    std::map<std::string, int> &hTypes = scn->getDspHistoTypes();
    std::map<std::string, int>::iterator it;
    for (it = hTypes.begin(); it != hTypes.end(); ++it)  {
      std::string name = (*it).first;
      PixScan::HistogramType type = (PixScan::HistogramType)((*it).second);
      if (scn->getHistogramFilled(type) && (scn->getHistogramKept(type) || (nloop==0 && scn->getLoopAction(0)!=PixScan::NO_ACTION)
					     || (nloop==1 && scn->getLoopAction(1)!=PixScan::NO_ACTION) )) {
	for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	  unsigned int mod = m_modules[pmod]->m_moduleId;
	  if(PMG_DEBUG) cout << "PixModuleGroup::scanLoopEnd : downloading histo " << type << endl;
	  scn->downloadHisto(m_pixCtrl, mod, type);
	}
      }
    }
  }

  // Create loop-summary RAW histograms
  if (nloop == 0) {
    if (scn->getHistogramFilled(PixScan::RAW_DATA_0)) {
      // Module loop
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	if (m_modules[pmod]->m_readoutActive) {
          unsigned int mod = m_modules[pmod]->m_moduleId;
	  if (scn->getHistogramFilled(PixScan::RAW_DATA_1) || scn->getHistogramFilled(PixScan::RAW_DATA_DIFF_1)) {
	    bool first = true;
	    Histo *h1, *hDiff1;
	    int nl = scn->getLoopVarNSteps(0);
	    for (int il=0; il<nl; il++) {
	      Histo& h = scn->getHisto(PixScan::RAW_DATA_0, mod, scn->scanIndex(2), scn->scanIndex(1), il);
	      if (first) {
		if (scn->getHistogramFilled(PixScan::RAW_DATA_1) && scn->getHistogramKept(PixScan::RAW_DATA_1)) {
		  std::ostringstream nam, tit;
		  nam << "RAW_1_" << mod << "_" << scn->scanIndex(1) << "_" << scn->scanIndex(2);
		  tit << "Raw data 1 mod " << mod << " L1=" << scn->scanIndex(1) << " L2=" << scn->scanIndex(2);
		  h1 = new Histo(nam.str(), tit.str(), h.nBin(0), -0.5, h.nBin(0)-0.5, nl, -0.5, nl-0.5);
		  scn->addHisto(*h1, PixScan::RAW_DATA_1, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
		}
		if (scn->getHistogramFilled(PixScan::RAW_DATA_DIFF_1)) {
		  std::ostringstream nam, tit;
		  nam << "RAW_DIFF_1_" << mod << "_" << scn->scanIndex(1) << "_" << scn->scanIndex(2);
		  tit << "Raw data diff 1 mod " << mod << " L1=" << scn->scanIndex(1) << " L2=" << scn->scanIndex(2);
		  hDiff1 = new Histo(nam.str(), tit.str(), h.nBin(0), -0.5, h.nBin(0)-0.5, nl, -0.5, nl-0.5);
		  scn->addHisto(*hDiff1, PixScan::RAW_DATA_DIFF_1, mod, scn->scanIndex(2), scn->scanIndex(1), -1);
		}
		first = false;
	      }
	      for (int id=0; id<h.nBin(0); id++) {
		if (scn->getHistogramFilled(PixScan::RAW_DATA_1) && scn->getHistogramKept(PixScan::RAW_DATA_1)) {
		  h1->set(id, il, h(id, 0));
		}
		if (scn->getHistogramFilled(PixScan::RAW_DATA_DIFF_1)) {
		  hDiff1->set(id, il, h(id, 1));
		}
	      }
	    }
	  }
	  if (!scn->getHistogramKept(PixScan::RAW_DATA_0)) scn->clearHisto(mod, PixScan::RAW_DATA_0);
	}
      }
    }
  } else if (nloop == 1) {
    if (scn->getHistogramFilled(PixScan::RAW_DATA_DIFF_1)) {
      // Module loop
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	if (m_modules[pmod]->m_readoutActive) {
          unsigned int mod = m_modules[pmod]->m_moduleId;
	  if (scn->getHistogramFilled(PixScan::RAW_DATA_DIFF_2) || scn->getHistogramKept(PixScan::RAW_DATA_DIFF_2)) {
	    bool first = true;
	    Histo *hDiff2;
	    int nl = scn->getLoopVarNSteps(1);
	    for (int il=0; il<nl; il++) {
	      Histo& h = scn->getHisto(PixScan::RAW_DATA_DIFF_1, mod, scn->scanIndex(2), il, -1);
	      if (first) {
		std::ostringstream nam, tit;
		nam << "RAW_DIFF_2_" << mod << "_" << scn->scanIndex(2);
		tit << "Raw data diff 2 mod " << mod << " L2=" << scn->scanIndex(2);
		hDiff2 = new Histo(nam.str(), tit.str(), h.nBin(1), -0.5, h.nBin(1)-0.5, nl, -0.5, nl-0.5);
		scn->addHisto(*hDiff2, PixScan::RAW_DATA_DIFF_2, mod, scn->scanIndex(2), -1, -1);
		first = false;
	      }
	      for (int id=0; id<h.nBin(1); id++) {
		hDiff2->set(id, il, h(0, id));
	      }
	    }
	  }
	  if (!scn->getHistogramKept(PixScan::RAW_DATA_DIFF_1)) scn->clearHisto(mod, PixScan::RAW_DATA_DIFF_1);
	}
      }
    }
  }
  // Execute the end-loop action
  switch (scn->getLoopAction(nloop)) {
  case PixScan::SCURVE_FIT:
    if (scn->getDspLoopAction(nloop)) {
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	scn->downloadHisto(m_pixCtrl, mod, PixScan::SCURVE_MEAN);
	scn->downloadHisto(m_pixCtrl, mod, PixScan::SCURVE_SIGMA);
	scn->downloadHisto(m_pixCtrl, mod, PixScan::SCURVE_CHI2);
      }
    } else {
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	scn->calcThr(m_pixCtrl, mod, scn->scanIndex(2), scn->scanIndex(1));
      }
    }
    break;
  case PixScan::TOTCAL_FIT:
    if (scn->getDspLoopAction(nloop)) {
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	scn->downloadHisto(m_pixCtrl, mod, PixScan::TOTCAL_PARA);
	scn->downloadHisto(m_pixCtrl, mod, PixScan::TOTCAL_PARB);
	scn->downloadHisto(m_pixCtrl, mod, PixScan::TOTCAL_PARC);
	scn->downloadHisto(m_pixCtrl, mod, PixScan::TOTCAL_CHI2);
      }
    } else {
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	scn->calcTotCal(m_pixCtrl, mod, scn->scanIndex(2), scn->scanIndex(1));
      }
    }
    break;
  case PixScan::TOTCAL_FEI4:
	TOTcalib_FEI4(nloop, scn, m_pixCtrl);
	break;
  case PixScan::TDAC_TUNING:
    endTDACTuning(nloop, scn);
    break;
  case PixScan::TDAC_FAST_TUNING:
	  endTDACFastTuning(nloop, scn);
	  break;
 case PixScan::OCC_SUM:
	 /*if (scn->scanIndex(nloop) == scn->getLoopVarNSteps(nloop)-1)
	 cout << "call endOccSumming" << endl;*/
	  endOccSumming(nloop, scn);
	  break;
  case PixScan::GDAC_FAST_TUNING:
	  endGDACFastTuning(nloop, scn);
	  break;
  case PixScan::GDAC_TUNING:
    endGDACTuning(nloop, scn);
    break;
  case PixScan::FDAC_TUNING:
    endFDACTuning(nloop, scn);
    break;
  case PixScan::FDAC_TUNING_ALT:
    endFDACTuningAlt(nloop, scn);
    break;
  case PixScan::IF_TUNING:
    endIFTuning(nloop, scn);
    break;
  case PixScan::T0_SET:
    endT0Set(nloop, scn);
    break;
  case PixScan::MCCDEL_FIT:
    for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
      unsigned int mod = m_modules[pmod]->m_moduleId;
      scn->calcThr(m_pixCtrl, mod, scn->scanIndex(2), scn->scanIndex(1), true);
    }
    //    mccDelFit(nloop, scn);
    break;
  case PixScan::OFFSET_CALIB:
    fitCalib(nloop, scn);
    break;
  case PixScan::CLEAR_IOMUX_BITS:
	  if(usb!=0){
		  usb->setIOMUXin(0);
		  usb->updateRegs();
	  }
  case PixScan::NO_ACTION:
  default:
    break;
  }

  //JW: clear SRAM of USB board after loop 0 ends
  if(nloop == 0 && scn->getLoopActive(1) && usb) usb->setRunMode();

  // clear temporary histos
  std::map<std::string, int> &hTypes = scn->getDspHistoTypes();
  std::map<std::string, int>::iterator it;
  for (it = hTypes.begin(); it != hTypes.end(); ++it)  {
    std::string name = (*it).first;
    PixScan::HistogramType type = (PixScan::HistogramType)((*it).second);
    if (scn->getHistogramFilled(type) && !scn->getHistogramKept(type) && ((nloop==0 && scn->getLoopAction(0)!=PixScan::NO_ACTION)
									  || (nloop==1 && scn->getLoopAction(1)!=PixScan::NO_ACTION) )){
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	unsigned int mod = m_modules[pmod]->m_moduleId;
	PixScanHisto *sc = 0;
	PixScanHisto &sh = scn->getHisto(type);
	if (nloop == 0) {
	  sc = &(sh[mod][scn->scanIndex(2)][scn->scanIndex(1)]);
	}else if (nloop == 1) {
	  sc = &(sh[mod][scn->scanIndex(2)]);
	}else if (nloop == 2){
	  sc = &(sh[mod]);
	}
	// just clear histograms and leave scan entry intact
	// otherwise will fail at higher loops...
	for(unsigned int isc=1;isc<sc->size();isc++) ((*sc)[isc]).clear();
      }
    }
  }

}

void PixModuleGroup::terminateScan(PixScan *scn){

  //This method has to execute end-of-loop actions.

  // Restore module configurations if needed
  for (unsigned int pmod=0; pmod<m_modules.size(); pmod++){
    if (scn->getRestoreModuleConfig()) {
      for (int gr=0; gr<MAX_GROUPS;gr++){
	if ((scn->getModuleMask(gr)) & ((0x1)<<(m_modules[pmod]->m_moduleId))) {
	  m_modules[pmod]->restoreConfig("PreScanConfig");
	  m_modules[pmod]->deleteConfig("PreScanConfig");
	}
      }
    }
    m_modules[pmod]->writeConfig(); //download the config again
  }
  USBPixController *usb = dynamic_cast<USBPixController*>(m_pixCtrl);
  if(usb) usb->setRunMode();

}

void PixModuleGroup::setupMasks(int nloop, PixScan *scn){
  // translate mask options
  int shiftmask;
  int capMask[3]={SHIFT_CAP1, SHIFT_CAP0, SHIFT_CAP1+SHIFT_CAP0};
  switch(scn->getMaskStageMode()) {
  case PixScan::SEL:          shiftmask = capMask[scn->getChargeInjCap()]; break;
  case PixScan::ENA:          shiftmask = SHIFT_ENABLE; break;
  case PixScan::SEL_ENA:      shiftmask = SHIFT_ENABLE+capMask[scn->getChargeInjCap()]; break;
  case PixScan::HITB:         shiftmask = SHIFT_HITBUS; break;
  case PixScan::HITB_INV:     shiftmask = 0x10; break;
  case PixScan::ENA_HITB:     shiftmask = SHIFT_ENABLE+SHIFT_HITBUS; break;
  case PixScan::SEL_ENA_HITB: shiftmask = SHIFT_ENABLE+SHIFT_HITBUS+capMask[scn->getChargeInjCap()]; break;
  case PixScan::STATIC:       shiftmask = 0x00; break;
  default:                    shiftmask = 0x00; break; // static
  }
  if(scn->getMaskStageIndex()==0){
    // set up mask for mask step 0
    if(PMG_DEBUG) cout << "PMG DEBUG: setting up masks" <<endl;
    m_pixCtrl->writeScanConfig(*scn);
  } else{
  // shift mask
    if(PMG_DEBUG) cout << "PMG DEBUG: shifting masks "<< shiftmask <<endl;
    m_pixCtrl->shiftPixMask(shiftmask);
  }
}

void PixModuleGroup::setupScanVariable(int nloop, PixScan *scn){
  //This method informs the scan about the variable that must be changed

  if (scn->getLoopActive(nloop) && !scn->getDspProcessing(nloop)) {
    PixScan::ScanParam par = scn->getLoopParam(nloop);
    int level = -1;
    if (par == PixScan::NO_PAR)         level = 0;
    if (par == PixScan::STROBE_DELAY)   level = 2;
    if (par == PixScan::STROBE_DURATION)level = 0;
    if (par == PixScan::TRIGGER_DELAY)  level = 0;
    if (par == PixScan::INCR_LAT_TRGDEL)level = 0;
    if (par == PixScan::CAPSEL)         level = 0;
    if (par == PixScan::VCAL)           level = 2;
    if (par == PixScan::BOC_BPH)        level = 0;
    if (par == PixScan::BOC_BPMPH)      level = 0;
    if (par == PixScan::BOC_VFINE)      level = 0;
    if (par == PixScan::BOC_VPH1)       level = 0;
    if (par == PixScan::BOC_VPH0)       level = 0; // 0 means one settings for all the modules
    if (par == PixScan::CAPMEAS)        level = 0;
    if (par == PixScan::IREF_PAD)       level = 0;
	if (par == PixScan::IOMUX_IN)		level = 0;
    if (par == PixScan::BOC_TX_CURR)    level = 1;
    if (par == PixScan::BOC_TX_MS)      level = 1;
    if (par == PixScan::BOC_RX_THR)     level = 1;
    if (par == PixScan::BOC_RX_DELAY)   level = 1; // 1 means module specific settings
    if (par == PixScan::GDAC)           level = 2;
    if (par == PixScan::ID)             level = 2;
    if (par == PixScan::IP)             level = 2;
    if (par == PixScan::IP2)            level = 2;
    if (par == PixScan::ITH1)           level = 2;
    if (par == PixScan::ITH2)           level = 2;
    if (par == PixScan::IL)             level = 2;
    if (par == PixScan::IL2)            level = 2;
    if (par == PixScan::TRIMT)          level = 2;
    if (par == PixScan::TRIMF)          level = 2;
    if (par == PixScan::IVDD2)          level = 2;
    if (par == PixScan::FEI4_GR)        level = 2;
    if (par == PixScan::IF)             level = 2; // 2 means FE specific settings
    if (par == PixScan::FDACS)          level = 3;
    if (par == PixScan::TDACS)          level = 3; // 3 means pixel specific settings

    bool reloadCfg = false;

    double val = (scn->getLoopVarValues(nloop))[scn->scanIndex(nloop)];
    int ival = (int)val;
    double pval = 0.;
    if(scn->scanIndex(nloop)>0) pval = (scn->getLoopVarValues(nloop))[scn->scanIndex(nloop)-1];
    int ipval = (int)pval;

    RodPixController *rod = dynamic_cast<RodPixController*>(m_pixCtrl);
    USBPixController *upc = dynamic_cast<USBPixController*>(m_pixCtrl);

    if(PMG_DEBUG) cout << "PixModuleGroup::setupScanVariable at loop " << nloop << " running on level " << level << endl;

    if (level == 0) {
      // Global settings
      switch (par) {
//       case PixScan::STROBE_DELAY:
// 	scn->setStrobeMCCDelay(ival);
// 	break;
      case PixScan::CAPMEAS: // re-defined to control aux. clk on USB board
	if(upc!=0){
	  upc->setAuxClkDiv(ival);
	  upc->updateRegs();
	}
	break;
      case PixScan::IREF_PAD:
 	if(upc!=0){
	  upc->setIrefPads(ival);
	  upc->updateRegs();
	}
	  case PixScan::IOMUX_IN:
	  if(upc!=0){
		  upc->setIOMUXin(ival);
		  upc->updateRegs();
	  }
	break;
     case PixScan::STROBE_DURATION:
	scn->setStrobeDuration(ival);
	break;
      case PixScan::TRIGGER_DELAY:
	scn->setStrobeLVL1Delay(ival);
	break;
      case PixScan::INCR_LAT_TRGDEL:
	scn->setStrobeLVL1Delay(scn->getStrobeLVL1Delay()+ival-ipval);
	scn->setLVL1Latency(scn->getLVL1Latency()-ival+ipval);
	break;
      case PixScan::CAPSEL:
	scn->setChargeInjCap(ival);
	break;
      case PixScan::VCAL:
	scn->setFeVCal(ival);
	break;
      case PixScan::BOC_BPH:
	if (rod) rod->setBocRegister("BRegClockPhase",ival);
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	break;
      case PixScan::BOC_BPMPH:
	if (rod) rod->setBocRegister("BpmClockPhase",ival);
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	break;
      case PixScan::BOC_VPH0:
	if (rod) rod->setBocRegister("VernierClockPhase0",ival);
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	break;
      case PixScan::BOC_VPH1:
	if (rod) rod->setBocRegister("VernierClockPhase1",ival);
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	break;
      case PixScan::BOC_VFINE:
	if (rod) rod->setBocRegister("VernierFinePhase",ival);
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	break;
      case PixScan::BOC_BVPH:
	if (rod) {
	  unsigned int bph = rod->getBocRegister("BRegClockPhase");
	  unsigned int vph = rod->getBocRegister("VernierClockPhase0") + rod->getBocRegister("VernierClockPhase1");
	  rod->setBocRegister("BRegClockPhase",ival);
	  vph = ival + (vph-bph);
	  if (vph < 25) {
	    rod->setBocRegister("VernierClockPhase0",vph);
	    rod->setBocRegister("VernierClockPhase1",0);
	  } else {
	    rod->setBocRegister("VernierClockPhase0",24);
	    rod->setBocRegister("VernierClockPhase1",vph-24);
	  }
	}
	else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	  case PixScan::NO_PAR:
	  default:
	break;
	  }
    } else {
      // Module loop
      for (unsigned int pmod=0; pmod<m_modules.size(); pmod++) {
	if (m_modules[pmod]->m_readoutActive) {
	  if (level == 1) {
	    // Module specific settings
	    if (rod) {
	      int fmtOff1[8] = { 2, 6, 38, 42, 50, 54, 86, 90 };
	      int fmtOff2[8] = { 2, 4, 6, 8, 86, 88, 90, 92 };
	      // Compute BOC links in case the record is not found in DB (assuming standard BOC mapping)
	      int il, bil, ol[4], bol[2];
	      bil    = ((ConfInt&)(m_modules[pmod]->config())["general"]["BocInputLink"]).getValue();
	      il     = ((ConfInt&)(m_modules[pmod]->config())["general"]["InputLink"]).getValue();
	      bol[0] = ((ConfInt&)(m_modules[pmod]->config())["general"]["BocOutputLink1"]).getValue();
	      bol[1] = ((ConfInt&)(m_modules[pmod]->config())["general"]["BocOutputLink2"]).getValue();
	      ol[0]  = ((ConfInt&)(m_modules[pmod]->config())["general"]["OutputLink1"]).getValue();
	      ol[1]  = ((ConfInt&)(m_modules[pmod]->config())["general"]["OutputLink2"]).getValue();
	      ol[2]  = ((ConfInt&)(m_modules[pmod]->config())["general"]["OutputLink3"]).getValue();
	      ol[3]  = ((ConfInt&)(m_modules[pmod]->config())["general"]["OutputLink4"]).getValue();
	      if (bil < 0) {
		bil = il;
	      }
	      PixScan::MccBandwidth outSpeed = scn->getMccBandwidth();
	      if (bol[0] < 0) {
		if (outSpeed == PixScan::SINGLE_40 || outSpeed == PixScan::DOUBLE_40) {
		  bol[0] = fmtOff1[ol[0]/16] + ol[0]%16;
		} else {
		  bol[0] = fmtOff2[ol[0]/16] + (ol[0]%16)/2;
		}
	      }
	      if (bol[1] < 0) {
		if (outSpeed == PixScan::SINGLE_40 || outSpeed == PixScan::DOUBLE_40) {
		  bol[1] = fmtOff1[ol[1]/16] + ol[1]%16;
		} else {
		  bol[1] = fmtOff2[ol[2]/16] + (ol[2]%16)/2;
		}
	      }
	      if (outSpeed == PixScan::SINGLE_40 || outSpeed == PixScan::SINGLE_80) {
		bol[1] = -1;
	      }
	      //std::cout << "bil = " << bil << " bol = " << bol[0] << " " << bol[1] << std::endl;
	      for (int il = 0; il < 2; il++) {
		if (bol[il] >= 0) {
		  switch (par) {
		  case PixScan::BOC_RX_THR:
		    rod->setBocRegister("RxThreshold", bol[il], ival);
		    break;
		  case PixScan::BOC_RX_DELAY:
		    rod->setBocRegister("DataDelay", bol[il], ival);
		    break;
		  case PixScan::NO_PAR:
		  default:
		    break;
		  }

		}
	      }
	      switch (par) {
	      case PixScan::BOC_TX_CURR:
		rod->setBocRegister("LaserCurrent", bil, ival);
		break;
	      case PixScan::BOC_TX_MS:
		rod->setBocRegister("BpmMarkSpace", bil, ival);
		break;
	      case PixScan::BOC_TX_BPMF:
		rod->setBocRegister("BpmFineDelay", bil, ival);
		break;
	      case PixScan::BOC_TX_BPM:
		rod->setBocRegister("BpmCoarseDelay", bil, ival);
		break;
	      case PixScan::NO_PAR:
	      default:
		break;
	      }
	    } else if (upc){
	      switch (par) {
	      case PixScan::BOC_RX_DELAY:
		break;
	      case PixScan::NO_PAR:
	      default:
		break;
	      }
	    }
	    else cout<<"INFO: No RodPixController present! Aborting!"<<endl;
	  } else {
	    if(par == PixScan::VCAL)  scn->setFeVCal(ival);
	    // FE loop
	    for (std::vector<PixFe*>::iterator fe = m_modules[pmod]->feBegin(); fe != m_modules[pmod]->feEnd(); fe++) {
	      if (level == 2) {
		// FE specific settings (global reg)
		reloadCfg = true;
		PixFe* fei4a = dynamic_cast<PixFeI4A*>(*fe);
		PixFe* fei4b = dynamic_cast<PixFeI4B*>(*fe);
		if(fei4a!=0 || fei4b!=0){ // FE-I4 register names
		  switch (par) {
		  case PixScan::FEI4_GR:
		    if(PMG_DEBUG) cout << "Writing " << ival << " to register " << scn->getLoopFEI4GR(nloop) << endl;
		    (*fe)->writeGlobRegister(scn->getLoopFEI4GR(nloop), ival);
		    break;
		  case PixScan::VCAL:
		    (*fe)->writeGlobRegister("PlsrDAC", ival);
		    break;
		  case PixScan::GDAC:
		    if (scn->getLoopAction(nloop) != PixScan::GDAC_FAST_TUNING) {
		      if(PMG_DEBUG) cout << "setting GlobalRegister GDAC to " << ival << ", nloop is " << nloop << endl;
		      if(fei4a!=0){
			(*fe)->writeGlobRegister("Vthin_AltFine", ival&0xff);
			(*fe)->writeGlobRegister("Vthin_AltCoarse", (ival&0xff00)>>8);
		      }else if(fei4b!=0){
			(*fe)->writeGlobRegister("Vthin_AltFine", ival&0xff);
			(*fe)->writeGlobRegister("Vthin_AltCoarse", (ival&0x7f00)>>7);
			//(*fe)->writeGlobRegister("Vthin_AltFine", ival&0xff);
			//(*fe)->writeGlobRegister("Vthin_AltCoarse", (ival&0xff00)>>8);
		      }
		    }
		    break;
		  case PixScan::IF:
		    if(PMG_DEBUG) cout << "Setting PrmpVbpf to " << ival << endl;
		    (*fe)->writeGlobRegister("PrmpVbpf", ival);
		    break;
		  case PixScan::STROBE_DELAY:
		    if(PMG_DEBUG) cout << "Setting PlsrDelay to " << ival << endl;
		    (*fe)->writeGlobRegister("PlsrDelay", ival);
		    break;
		  default:
		    break;
		  }
		}else{ // FE-I2/3 register names
		  switch (par) {
		  case PixScan::GDAC:
		    if (scn->getLoopAction(nloop) != PixScan::TDAC_TUNING) {
		      (*fe)->writeGlobRegister("GLOBAL_DAC", ival);
		    }
		    break;
		  case PixScan::IF:
		    (*fe)->writeGlobRegister("DAC_IF", ival);
		    break;
		  case PixScan::ID:
		    (*fe)->writeGlobRegister("DAC_ID", ival);
		    break;
		  case PixScan::IP:
		    (*fe)->writeGlobRegister("DAC_IP", ival);
		    break;
		  case PixScan::IP2:
		    (*fe)->writeGlobRegister("DAC_IP2", ival);
		    break;
		  case PixScan::ITH1:
		    (*fe)->writeGlobRegister("DAC_ITH1", ival);
		    break;
		  case PixScan::ITH2:
		    (*fe)->writeGlobRegister("DAC_ITH2", ival);
		    break;
		  case PixScan::IL:
		    (*fe)->writeGlobRegister("DAC_IL", ival);
		    break;
		  case PixScan::IL2:
		    (*fe)->writeGlobRegister("DAC_IL2", ival);
		    break;
		  case PixScan::TRIMT:
		    (*fe)->writeGlobRegister("DAC_ITRIMTH", ival);
		    break;
		  case PixScan::TRIMF:
		    (*fe)->writeGlobRegister("DAC_ITRIMIF", ival);
		    break;
		  case PixScan::IVDD2:
		    (*fe)->writeGlobRegister("DAC_IVDD2", ival);
		    break;
		  case PixScan::VCAL:
		    (*fe)->writeGlobRegister("DAC_VCAL", ival);
		    break;
		  default:
		    break;
		  }
		}
	      } else if (level == 3) {
		reloadCfg = true;
		// Pixel specific settings (masks)
		ConfMask<unsigned short int> *mask = NULL;
		if (par == PixScan::TDACS) {
		  if ((scn->getLoopAction(nloop) != PixScan::TDAC_TUNING) && (scn->getLoopAction(nloop) != PixScan::TDAC_FAST_TUNING)) mask = &((*fe)->readTrim("TDAC"));
		}
		if (par == PixScan::TDACS_VARIATION) {
		  if (scn->getLoopAction(nloop) != PixScan::MIN_THRESHOLD) mask = &((*fe)->readTrim("TDAC"));
		}
		if (par == PixScan::FDACS) mask = &((*fe)->readTrim("FDAC"));
		if (mask != NULL) {
		  for (unsigned int col=0; col<(*mask).get().size(); col++) {
		    for (unsigned int row=0; row<(*mask)[0].size(); row++) {
		      if (par == PixScan::TDACS_VARIATION) {
			int maxtrim;
			(*fe)->getTrimMax("TDAC", maxtrim);
			if ((*mask)[col][row] + ival > maxtrim) {
			  (*mask)[col][row] = maxtrim;
			} else if ((*mask)[col][row] + ival < 0) {
			  (*mask)[col][row] = 0;
			} else {
			  (*mask)[col][row] += ival;
			}
		      } else {
			(*mask)[col][row] = ival;
		      }
		    }
		  }
		}
	      }
	    }
	  }
	}
	if (reloadCfg) m_modules[pmod]->writeConfig(); //download the config again
      }
    }
  }
}
void PixModuleGroup::TOTcalib_FEI4(int nloop, PixScan *scn, PixController *ctrl){

	int NoOfTOTHistos = PixScan::TOT15 - PixScan::TOT0 +1;
	if(PMG_DEBUG) cout << "JR: number of histograms = " <<  NoOfTOTHistos << endl;

	// Check if TOT histos are present
	//for(int tot=PixScan::TOT0; tot<PixScan::TOT15+1; tot++){
	//	if (m_histo.find((PixLib::PixScan::HistogramType)(tot)) == m_histo.end()) {
	//		throw PixScanExc(PixControllerExc::ERROR, "TOT histos not found");
	//	}
	//}
	std::vector< Histo* > vTOTHisto; //vector contains all TOT histograms in order: TOT0(scanpoint 0), TOT1(scanpoint 0),..., TOT15(scanpoint 0), TOT0(scanpoint 1),...
	
	int nr_hits;

	int steps = scn->getLoopVarNSteps(0);
	double start = scn->getLoopVarMin(0);
	double end = scn->getLoopVarMax(0);
	double stepsize = (end-start)/(steps-1);

	if(PMG_DEBUG) cout << "JR : steps = " << steps << " stepsize " << stepsize  << " start "<< start << " end " << end << endl;
	//double cap = scn->getChargeInjCap();   //Gives only the number of the capacitor

	double VCALsumsqr = 0;
	double VCALsum = 0;
	double hitsum = 0;
	double VCALsqravg = 0;
	double VCALavg = 0;
	double VCALsigma = 0;

	double aVCALsqravg[13];
	double aVCALavg[13];
	double aVCALsigma[13];

	double dtot[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
	double aQavg[13], aQsigma[13];

  // Module loop
	//for (unsigned int pmod=0; pmod<m_modules.size(); pmod++){
	//	if (m_modules[pmod]->m_readoutActive){
	//		unsigned int mod = m_modules[pmod]->m_moduleId;
			int mod = 0;

			Histo *parA, *parB, *parC, *chi2;
			std::string nam, tit;
			std::ostringstream mnum;
			mnum << mod;
			nam = "ParA_" + mnum.str();
			tit = "Parameter A mod " + mnum.str();
			parA = new Histo(nam, tit, N_I4_PIXEL_COLUMNS, -0.5, N_I4_PIXEL_COLUMNS-0.5, N_I4_PIXEL_ROWS, -0.5, N_I4_PIXEL_ROWS-0.5);
			nam = "ParB_" + mnum.str();
			tit = "Parameter B mod " + mnum.str();
			parB = new Histo(nam, tit, N_I4_PIXEL_COLUMNS, -0.5, N_I4_PIXEL_COLUMNS-0.5, N_I4_PIXEL_ROWS, -0.5, N_I4_PIXEL_ROWS-0.5);
			nam = "ParC_" + mnum.str();
			tit = "Parameter C mod " + mnum.str();
			parC = new Histo(nam, tit, N_I4_PIXEL_COLUMNS, -0.5, N_I4_PIXEL_COLUMNS-0.5, N_I4_PIXEL_ROWS, -0.5, N_I4_PIXEL_ROWS-0.5);
			nam = "Chi2_" + mnum.str();
			tit = "Chi2 mod " + mnum.str();
			chi2 = new Histo(nam, tit, N_I4_PIXEL_COLUMNS, -0.5, N_I4_PIXEL_COLUMNS-0.5, N_I4_PIXEL_ROWS, -0.5, N_I4_PIXEL_ROWS-0.5);

			FitClass fc;
			int funcID=fc.getFuncID("ToT-calibration FE-I4 Polynomial 2nd order");
			if(funcID<0)
				throw PixScanExc(PixControllerExc::ERROR, "Can't define fit function");

			double par[7]={100.,1000.,1000.,0.,1.,0.,0.};
			bool fix_par[7]={false,false,false,true,true,true,true};
			double xerr[13] = {0,0,0,0,0,0,0,0,0,0,0,0,0};

			int fe = 0;										//Only for single chips! This has to be changed!
			PixModule *pmod=ctrl->getModGroup().module(mod);

			// get FE calib. and put into par[3...6]
			float vcalG0 = (dynamic_cast<ConfFloat &>(pmod->pixFE(fe)->config()["Misc"]["VcalGradient0"])).value();
			float vcalG1 = (dynamic_cast<ConfFloat &>(pmod->pixFE(fe)->config()["Misc"]["VcalGradient1"])).value();
			float vcalG2 = (dynamic_cast<ConfFloat &>(pmod->pixFE(fe)->config()["Misc"]["VcalGradient2"])).value();
			float vcalG3 = (dynamic_cast<ConfFloat &>(pmod->pixFE(fe)->config()["Misc"]["VcalGradient3"])).value();
			std::string capLabels[3]={"CInjLo", "CInjMed", "CInjHi"};
			int chargeInjCap = scn->getChargeInjCap();
			float cInj     = (dynamic_cast<ConfFloat &>(pmod->pixFE(fe)->config()["Misc"][capLabels[chargeInjCap]])).value();
			if(PMG_DEBUG) cout << "PixModuleGroup::calcTotCal : using inj. capacitance of " << cInj << " (switch was " << chargeInjCap << ")" << endl;
			if(PMG_DEBUG) cout << "JR: VCALGrad0 = " << vcalG0 << " vcalG1 " << vcalG1 << " vcalG2 " << vcalG2 << " vcalG3 " << vcalG3 << endl;
			cInj /= 0.160218;									
			par[3] = vcalG0*cInj;
			par[4] = vcalG1*cInj;
			par[5] = vcalG2*cInj;						
			par[6] = vcalG3*cInj;
			
			for(int i=0; i<steps; i++){					//scn->getLoopVarNSteps(i) gives you the number of steps in loop i
				// storing the TOT histograms in vTOTHisto vector for all scanpoints
				for(int tot=PixScan::TOT0; tot<PixScan::TOT15+1; tot++){
					vTOTHisto.push_back((Histo*)(&scn->getHisto((PixLib::PixScan::HistogramType)(tot), mod, 0, 0, i)));
				}
			}

			//TH1D *ToTHistos[16];
			//TCanvas * KCanvas[16];

			for(int col=0; col<N_I4_PIXEL_COLUMNS; col++){
				for(int row=0; row<N_I4_PIXEL_ROWS; row++){
					for(int tot=0; tot<NoOfTOTHistos-3; tot++){
						VCALsumsqr = 0;
						VCALsum = 0;
						hitsum = 0;
						aVCALsqravg[tot] = 0;
						aVCALavg[tot] = 0;
						aVCALsigma[tot] = 0;
						for(int j=0; j<steps; j++){
							VCALsumsqr += ((*vTOTHisto[tot+j*NoOfTOTHistos])(col, row)) * (start + stepsize * j) * (start + stepsize * j);
							VCALsum += ((*vTOTHisto[tot+j*NoOfTOTHistos])(col, row)) * (start + stepsize * j);
							hitsum += ((*vTOTHisto[tot+j*NoOfTOTHistos])(col, row));
						}
						if(hitsum > 1) {
							//VCALsqravg =  VCALsumsqr / hitsum;
							//VCALavg = VCALsum / hitsum;
							//VCALsigma = (VCALsqravg - VCALavg * VCALavg)/(hitsum-1);

							aVCALsqravg[tot] =  VCALsumsqr / hitsum;
							aVCALavg[tot] = VCALsum / hitsum;
							aVCALsigma[tot] = (aVCALsqravg[tot] - aVCALavg[tot] * aVCALavg[tot])/(hitsum-1.0);
						}else{
							//VCALsqravg = 0;
							//VCALavg = 0;
							//VCALsigma = 0;

							aVCALsqravg[tot] = 0;
							aVCALavg[tot] = 0;
							aVCALsigma[tot] = 0;
						}
						if(PMG_DEBUG){
							if(col == 5 && row == 6){
								//cout << "JR: pixel 5,6 VCALsqravg = "<<VCALsqravg<< " VCALavg "<<VCALavg<< " VCALsigma "<< VCALsigma << endl;
								cout << "JR: pixel 5,6 aVCALsqravg = "<<aVCALsqravg[tot]<< " aVCALavg "<<aVCALavg[tot]<< " aVCALsigma "<< aVCALsigma[tot] << " hitsum " << hitsum << endl;
							}
						}	
					}
					//This is where the fit happens. The function returns chi2. Vcal values are converted into charge before.
					
					for(int tot=0; tot<NoOfTOTHistos-3; tot++){
						aQavg[tot]=par[3]+par[4]*aVCALavg[tot];
						aQsigma[tot]=par[3]+par[4]*aVCALsigma[tot];
						if(PMG_DEBUG && (col == 5 && row == 6)) 
						  cout << "JR: aVCALavg["<<tot<<"] = "<< aVCALavg[tot]<<"\t aQavg["<<tot<<"] = "<< aQavg[tot]<<endl; 
					}
					
					double c2 = fc.runFit(NoOfTOTHistos-3, dtot, aQavg, xerr, aQsigma, par, fix_par,funcID, 0., 0.);
					if(c2<=1e-2 || c2>1e20){
					  if(PMG_DEBUG && (col == 5 && row == 6) || (col == 10 && row == 11))
							cout<<"DEBUG: unreasonable chi^2"<<endl;
						parA->set(col,row,0.);
						parB->set(col,row,0.);
						parC->set(col,row,0.);
					} else{
						if(PMG_DEBUG && (col == 5 && row == 6) || (col == 10 && row == 11))
							cout<<"DEBUG: GOOD chi^2"<<endl;
						parA->set(col,row,par[0]);
						parB->set(col,row,par[1]);
						parC->set(col,row,par[2]);
					}
					chi2->set(col,row,c2);

					if(PMG_DEBUG && (col == 5 && row == 6) || (col == 10 && row == 11)){
					  cout << "JR: pixel 5,6 par0 = "<< par[0]  << " par1 "<< par[1] << " par2 "<< par[2]<< " chi2 " << c2 << endl;
					}
				}
			}

			int ix0 = scn->scanIndex(0);
			int ix1 = scn->scanIndex(1);
			int ix2 = scn->scanIndex(2);

			scn->addHisto(*parA, PixScan::TOTCAL_PARA, mod, ix2, ix1, -1);
			scn->addHisto(*parB, PixScan::TOTCAL_PARB, mod, ix2, ix1, -1);
			scn->addHisto(*parC, PixScan::TOTCAL_PARC, mod, ix2, ix1, -1);
			scn->addHisto(*chi2, PixScan::TOTCAL_CHI2, mod, ix2, ix1, -1);  
			scn->setHistogramFilled(PixScan::TOTCAL_PARA, true);
			scn->setHistogramFilled(PixScan::TOTCAL_PARB, true);
			scn->setHistogramFilled(PixScan::TOTCAL_PARC, true);
			scn->setHistogramFilled(PixScan::TOTCAL_CHI2, true);
	//	}
//	}
}

void PixModuleGroup::fillHistosPerDc(PixScan &scn, int ife, Histo *his, double value){
  int DC = 0;
  std::vector<float> scvals = scn.getLoopVarValues(0);
  int scpt = scn.scanIndex(0);
  if(scpt>=0 && scpt<(int)scvals.size()) DC = (int)scvals[scpt];
  if(scn.getMaskStageTotalSteps()==PixScan::STEPS_1){
    unsigned int col, row;
    for(col = (2*DC); col<(unsigned int)(2*DC+2);col++){
      for(row=0;row<N_I4_PIXEL_ROWS;row++){
	if(row>=0 && row<N_I4_PIXEL_ROWS && col>=0 && col<(N_I4_PIXEL_COLUMNS*N_I4_PIXEL_FE_CHIPS)){
	  if(PMG_DEBUG) cout << "Filling value " << value << " into col,row " << col << ", " << row << endl;
	  his->set(col,row, value);
	}else
	  cerr << "Wrong col,row number: " << col << " - " << row << endl;
      }
    }
  } else if(scn.getMaskStageTotalSteps()==PixScan::STEPS_26880){
    // DCS vector filled in order DC 0 mask 0...N, then DC1 mask 0...N, etc
    int Drow = scn.getMaskStageIndex();
    int col = DC*2+ife*N_I4_PIXEL_COLUMNS;
    int row = N_I4_PIXEL_ROWS-Drow-1;
    if(Drow>=N_I4_PIXEL_ROWS){
      col +=1;
      row = Drow-N_I4_PIXEL_ROWS;
    }
    if(row>=0 && row<N_I4_PIXEL_ROWS && col>=0 && col<(N_I4_PIXEL_COLUMNS*N_I4_PIXEL_FE_CHIPS)){
      if(PMG_DEBUG) cout << "Filling value " << value << " into col,row " << col << ", " << row << endl;
      his->set(col,row, value);
    }else
      cerr << "Wrong col,row number: " << col << " - " << row << endl;
  }
}