def calibrateAmplitudeAndOffset(self,f):
   """
   Only used when this pulse Analyser has to be used as real analyser, not when using it to see bifurcation
   """
   rowData=Datacube()
   for phi in arange(0,2*math.pi,math.pi/30):
     print "calibration : phi = %f deg" % (phi/math.pi*180)    
     self._pulseGenerator.clearPulse()
     self.clear()
     self._pulseGenerator.generatePulse(duration=20000, frequency=f, amplitude=0.6, DelayFromZero=0,useCalibration=True, phase=phi)
     self.addFrequency(f=f,useCorrection=False)    	
     self._pulseGenerator.sendPulse()
     time.sleep(0.5)
     (av, co, fr)= self.analyse()
     rowData.set(I=av[0,0], Q=av[1,0],phi=phi)          	
     rowData.commit()
   #I0=2/ptp(rowData['I'])
   #Q0=2/ptp(rowData['Q'])
   (I,Q,phi,dphi)=scipy.optimize.fmin_powell(lambda (I,Q,phi0,dphi): sum((I*rowData['I'] - sin(rowData['phi']+phi0+dphi))**2)+sum((Q*rowData['Q'] - cos(rowData['phi']+phi0))**2),(1,1,0,0))
   print (I,Q,phi,dphi)
   f_c=self._MWSource.frequency()
   df=f-f_c
   index=self._calibration.search(f_sb=df,f_c=f_c)
   if index!=None:
     self._calibration.removeRow(index)
   self._calibration.set(I=I,Q=Q,phi=dphi,f_c=f_c,f_sb=df)
   self._calibration.commit()
   self._calibration.savetxt()
   register['%s Cal'% self._name]=self._calibration.filename()
   return rowData

 def clicks(self):
 
   """
   Returns a datacube that contains the individual detector clicks for all measured samples in binary form.
   """
     
   clicks = Datacube("detector clicks",dtype = uint8)
   
   angles = self.bifurcationMapRotation()
   
   clicks1 = self.trends()[0]*cos(angles[0])+self.trends()[1]*sin(angles[0]) > 0
   clicks2 = self.trends()[2]*cos(angles[1])+self.trends()[3]*sin(angles[1]) > 0
   
   def mapper(t):
     (x1,x2) = t
     if x1 and x2:
       return 3
     elif x1 and not x2:
       return 1
     elif x2 and not x1:
       return 2
     else:
       return 0
       
   clicks.createColumn("clicks",map(mapper,zip(clicks1,clicks2)))
   
   return clicks

 def initCalibrationData(self):
   """
   Initialize the datacubes that contain the IQ calibration data.
   """
   self._offsetCalibrationData = Datacube()
   self._offsetCalibrationData.setName("IQ mixer calibration - Offset Calibration Data")
   self._powerCalibrationData = Datacube()
   self._powerCalibrationData.setName("IQ mixer calibration - Power Calibration Data")
   self._sidebandCalibrationData = Datacube()
   self._sidebandCalibrationData.setName("IQ mixer calibration - Sideband Mixing Calibration Data")

 def _findAmplitude(self,voltages,frequency,shape):
   """
   Only used by function calibrate.
   Measure and find the amplitude where the JBA is bi-evaluated, go to this point, and store this amplitude
   """
   ps = []
   max = 0
   maxVoltage = 0
   self.variances = zeros((len(voltages),2))
   self.variances[:,0] = voltages
   data = Datacube("Variance")
   for i in range(0,len(voltages)):
     if self.stopped():
       self._stopped = False
       raise Exception("Got stopped!")
     v = voltages[i] 
     self.init()
     self.frequency(frequency=frequency,shape=v*shape)
     (av,trends, fr)=self.analyse()
     varsum =cov(trends[0,:,0])+cov(trends[1,:,0])
     data.set(v = v)
     data.set(varsum=varsum)
     data.commit()
     self.notify("variance",(data.column("v"),data.column("varsum")))
     self.notify("status","Variance at %g V: %g" % (v,varsum))
     print "Variance at v = %f : %f" % (v,varsum)
     self.variances[i,1] = varsum
     ps.append(varsum)
     if varsum > max:
       max = varsum
       maxVoltage = voltages[i]
   self.frequency(frequency=frequency,shape=maxVoltage*shape)
   return (ps,max,maxVoltage,data)

 def _attenuatorRangeCheck(self,voltages):
   ps = []
   max = 0
   maxVoltage = 0
   self.variances = zeros((len(voltages),2))
   self.variances[:,0] = voltages
   data = Datacube("Variance")
   for i in range(0,len(voltages)):
     if self.stopped():
       self._stopped = False
       raise Exception("Got stopped!")
     v = voltages[i] 
     self._attenuator.setVoltage(v)
     self._acqiris.bifurcationMap(ntimes = 10)
     trends = self._acqiris.trends()
     varsum =cov(trends[self._params["acqirisChannel"]])+cov(trends[self._params["acqirisChannel"]+1])
     data.set(v = v)
     data.set(varsum=varsum)
     data.commit()
     self.notify("variance",(data.column("v"),data.column("varsum")))
     self.notify("status","Variance at %g V: %g" % (v,varsum))
     print "Variance: %f" % varsum
     self.variances[i,1] = varsum
     ps.append(varsum)
     if varsum > max:
       max = varsum
       maxVoltage = voltages[i]
   return (ps,max,maxVoltage)

  def measureSCurve(self,voltages = None,ntimes = 40,microwaveOff = True):
    self.notify("status","Measuring S curve...")
    def getVoltageBounds(v0,jba,variable,ntimes):
      v = v0
      jba.setVoltage(v)
      jba._acqiris.bifurcationMap(ntimes = ntimes)
      p = jba._acqiris.Psw()[variable]
      
      while p > 0.03 and v < v0*2.0:
        v*=1.05
        jba.setVoltage(v)
        jba._acqiris.bifurcationMap(ntimes = ntimes)
        p = jba._acqiris.Psw()[variable]
      vmax = v
      
      v = v0
      jba.setVoltage(v)
      jba._acqiris.bifurcationMap(ntimes = ntimes)
      p = jba._acqiris.Psw()[variable]
      
      while p < 0.98 and v > v0/2.0:
        v/=1.05
        jba.setVoltage(v)
        jba._acqiris.bifurcationMap(ntimes = ntimes)
        p = jba._acqiris.Psw()[variable]
      vmin = v
      return (vmin*0.95,vmax*1.2)

    try:
      v0 = self.voltage()
      state = self._qubitmwg.output()
      self._attenuator.turnOn()
      data = Datacube("S Curve")
      dataManager = DataManager()
      dataManager.addDatacube(data)
      if microwaveOff:
        self._qubitmwg.turnOff()
      if voltages == None:
        self.notify("status","Searching for proper voltage range...")
        (vmin,vmax) = getVoltageBounds(v0,self,self._params["variable"],ntimes)
        voltages = arange(vmin,vmax,0.005)
        self.notify("status","Measuring S curve in voltage range  [%g - %g]..." % (vmin,vmax))
      for v in voltages:
        self.setVoltage(v)
        self._acqiris.bifurcationMap(ntimes = ntimes)
        data.set(v = v)
        data.set(**(self._acqiris.Psw()))
        data.commit()
        self.notify("sCurve",(data.column("v"),data.column(self._params["variable"])))
    finally:
      self.notify("status","S curve complete.")
      self.setVoltage(v0)
      if state:
        self._qubitmwg.turnOn()

  def initialize(self, name, generator, analyser,magnitudeButton='formAmplitude'):
    """
    Initialize instrument, and define default shape in self.shape
    """
    instrumentManager=Manager()
    self._pulseGenerator=instrumentManager.getInstrument(generator)
    self._pulseAnalyser=instrumentManager.getInstrument(analyser)
    self._params=dict()
    self._params["pulseGenerator"]=generator
    self._params["pulseAnalyser"]=analyser
    self._change=True
    self._magnitudeButton=magnitudeButton
    self.data=Datacube()
    try:
      self._variableAttenuator=instrumentManager.getInstrument('jba_att')
    except:
      pass
    self._shapeParams=dict()
    self._shapeParams["risingTime"]    = 10
    self._shapeParams["plateauLength"] = 200
    self._shapeParams["latchLength"]   = 1000
    self._shapeParams["plateau"]       = 0.85

    self.generateShape()

    self.bit=int(self.name()[-1])-1
    self._phase=0

 def test_savetxt(self):
   """
   Test saving a datacube to a text file
   """
   for key in self.testCubes.keys():
     print "Checking plain text loading of test cube {0!s}".format(key)
     cube = self.testCubes[key]
     filename = os.path.normpath(self.dataPath+"/test_{0!s}.txt".format(key))
     cube.savetxt(filename,overwrite = True)
     
     self.assert_(os.path.exists(filename),"File {0!s} has not been created!".format(filename))
     self.assert_(os.path.isfile(filename))
     
     restoredCube = Datacube()
     restoredCube.loadtxt(filename)
           
     self.assert_(restoredCube.equal(cube),"Error: Restored datacube does not match original one!")

 def measureSCurves(self,ntimes=None,data=None):
   if data==None:
     if name==None:name='sCurves %s'%self.name()
     data=Datacube(name)
     dataManager.addDatacube(data)
   self._pulseGenerator._MWSource.setPower(self._rabiPower)
   self._pulseGenerator._MWSource.setFrequency(self._f01)
   self._pulseGenerator.generatePulse(duration=self._rabiDuration,frequency=self._f01,DelayFromZero=register['repetitionPeriod']/2-self._rabiDuration-10,useCalibration=False)
   self._pulseGenerator.sendPulse()
   off=Datacube('sOff')
   data.addChild(off)
   self._pulseGenerator._MWSource.turnOff()
   self._jba.measureSCurve(data=off,ntimes=10)    
   on=Datacube('sOn')
   data.addChild(on)
   self._pulseGenerator._MWSource.turnOn()
   self._jba.measureSCurve(data=on,ntimes=10)    

 def initCal(self):
   """
   Re-init the calibration when using this instrument as real analyser
   """
   self._calibration=Datacube()
   self._calibration.setName('analyser IQ mixer Calibration')
   self._calibration.savetxt()
   register['%s Cal'% self._name]=self._name.filename()

 def calibrateSidebandMixing(self,frequencyRange = None,sidebandRange = arange(-0.5,0.51,0.1)):
   """
   Calibrate the IQ mixer sideband generation.
   """
   if frequencyRange==None:
     frequencyRange=[self._mwg.frequency()]
   try:
     self.setup()
     params = dict()
     params["power"] = self._mwg.power()
     params["channels"] = self._awgChannels
     params["mwg"] = self._mwg.name()
     params["awg"] = self._awg.name()
     params["fsp"] = self._fsp.name()
     self.sidebandCalibrationData().setParameters(params)
     self._mwg.turnOn()
     channels = self._awgChannels
     self.loadSidebandWaveforms()
     for f_c in frequencyRange:
       #We round the center frequency to an accuracy of 1 MHz
       f_c = round(f_c,3)
       self._mwg.setFrequency(f_c)
       self._awg.setAmplitude(channels[0],4.5)
       self._awg.setAmplitude(channels[1],4.5)
       self._awg.setOffset(channels[0],self.iOffset(f_c))
       self._awg.setOffset(channels[1],self.qOffset(f_c))
       data = Datacube("f_c = %g GHz" % f_c)
       rowsToDelete = []
       try:
         for i in range(0,len(self._sidebandCalibrationData.column("f_c"))):
          if abs(self._sidebandCalibrationData.column("f_c")[i]-f_c) < 0.1:
             rowsToDelete.append(i)
       except:
         pass
       self._sidebandCalibrationData.removeRows(rowsToDelete)
       self._sidebandCalibrationData.addChild(data, f_c=f_c)
       self._sidebandCalibrationData.set(f_c = f_c)
       self._sidebandCalibrationData.commit()
       for f_sb in sidebandRange: 
         print "f_c = %g GHz, f_sb = %g GHz" % (f_c,f_sb)
         self._fsp.write("SENSE1:FREQUENCY:CENTER %f GHZ" % (f_c+f_sb))
         result = scipy.optimize.fmin_powell(lambda x,*args: self.measureSidebandPower(x,*args),[0,0],args = [f_sb],full_output = 1,xtol = 0.00001,ftol = 1e-4,maxiter = 2)
         params = result[0]
         value = result[1]
         print "f_c = %g GHz, f_sb = %g GHz, c = %g, phi = %g rad" % (f_c,f_sb,params[0],params[1])
         self.loadSidebandCalibrationWaveform(f_sb = f_sb,c = params[0],phi = params[1])
         for i in [-3,-2,-1,0,1,2,3]:
           self._fsp.write("SENSE1:FREQUENCY:CENTER %f GHZ" % (f_c+f_sb*i))
           if i < 0:
             suppl = "m"
           else:
             suppl = ""
           data.set(**{"p_sb%s%d" % (suppl,abs(i)) : self.measureAveragePower()})
         data.set(f_c = f_c,f_sb = f_sb,c = params[0],phi = params[1])
         data.commit()
       self._sidebandCalibrationData.sortBy("f_c")
       self._sidebandCalibrationData.savetxt()
   finally:
     self.teardown()
   return self._sidebandCalibrationData.filename()

def measureSpectroscopy(qubit,frequencies,data = None,ntimes = 20,amplitude = 1,measureAtReadout = False,delay = 0,f_sb = 0,delayAtReadout = 1500,pulseLength = 500,gaussian=False):
  if data == None:
    data = Datacube()
  if measureAtReadout:
    qubit.loadRabiPulse(length = pulseLength,f_sb = f_sb,delay = delay,gaussian=gaussian)
  else:
    qubit.loadRabiPulse(length = pulseLength,f_sb = f_sb,delay = delay,gaussian=gaussian)
  qubit.turnOnDrive()
  data.setParameters(dict(data.parameters(),**instrumentManager.parameters()))
  try:
    for f in frequencies:
      qubit.setDriveFrequency(f+f_sb)
      qubit.setDriveAmplitude(I = amplitude,Q = amplitude)
      acqiris.bifurcationMap(ntimes = ntimes)
      data.set(f = f)
      data.set(**acqiris.Psw())
      data.commit()
  except StopThread:
    return data

def findAnticrossing(voltageRead, voltageWrite, searchRange1, searchRange2, cube=None):
    if cube == None:
        cube = Datacube()
    spectro = Datacube()
    spectro.setName("Qubit 1 Spectroscopy")
    cube.addChild(spectro)
    jba.calibrate()
    params = getQubitFrequency(mwg, initialRange, variable, spectro)

  def getTrace(self,waitFullSweep = False,timeOut = 1600,fromMemory=False):
    """ 
    Get a raw trace in the VNA, without correcting the data, except for internal attenuators.
    Get the memory instead of main trace if fromMemory=True.
    Restart a sweep and wait for its completion if fromMemory=False and waitFullSweep=True.
    FOR INTERNAL USE ONLY.
    USE INSTEAD getFreqMagPhase(waitFullSweep = False,fromMemory=False,timeOut=60,addedAttenuators=0,unwindPhase=False,subtractedSlope=None,deltaPhase=None,offset=None).
    """
    trace = Datacube('Spectrum')
    handle = self.getHandle()
    handle.timeout = timeOut
    if waitFullSweep:
      print "Getting trace...",
      # freqs = self.ask_for_values("HLD;TRS;WFS;fma;msb;OFV;") 2011/12 VS 
      self.write('TRS;WFS;')
    freqs = self.ask_for_values('fma;msb;OFV;')
    data = self.write('fma;msb;')
    if(fromMemory):
      data = self.ask_for_values('MEM;OFD;')
      self.write('DTM;')
    else: 
      data = self.ask_for_values('OFD;')
    if waitFullSweep:
      print "done."
    freqs.pop(0)
    data.pop(0)
    mag = []
    phase = []
    #If y length is twice the x length, we got phase and magnitude.
    if len(data) == 2*len(freqs):
      for i in range(0,len(data)):
        if i%2 == 0:	
          mag.append(data[i])
        else:
				 phase.append(data[i])
    else:
      mag = data
    att=self.attenuation()
    trace.setParameters( {'attenuation':att,'power':self.totalPower()})
    trace.createCol(name='freq',values=freqs)
    trace.createCol(name='mag',values=array(mag)+att)
    if len(phase)!=0: trace.createCol(name='phase',values=phase)
    return trace

def measureSpectroscopy(qubit,frequencies,data = None,ntimes = 20,amplitude = 0.1,measureAtReadout = False,delay = 0,f_sb = 0):
  if data == None:
    data = Datacube()
  if measureAtReadout:
    qubit.loadRabiPulse(length = 500,readout = qubit.parameters()["timing.readout"]+500,f_sb = f_sb,delay = delay)
  else:
    qubit.loadRabiPulse(length = 500,readout = qubit.parameters()["timing.readout"],f_sb = f_sb,delay = delay)
  qubit.turnOnDrive()
  data.setParameters(instrumentManager.parameters())
  try:
    for f in frequencies:
      qubit.setDriveFrequency(f+f_sb)
      qubit.setDriveAmplitude(I = amplitude,Q = amplitude)
      acqiris.bifurcationMap(ntimes = ntimes)
      data.set(f = f)
      data.set(**acqiris.Psw())
      data.commit()
  finally:
    return data

def ramsey(delays,cube = None,ntimes = 20,length = 20):
	if cube == None:
		ramseyData = Datacube()
	else:
		ramseyData = cube
	ramseyData.setParameters(instruments.parameters())
	for delay in delays:
		generateRamseyWaveform(length = length,delay = delay)
		acqiris.bifurcationMap(ntimes = ntimes)
		ramseyData.set(delay = delay,**acqiris.Psw())
		ramseyData.commit()
	return rabiData

 def initialize(self, name, acqiris, MWSource,pulse_generator):
   """
   Initialize the instrument
   """
   instrumentManager=Manager()
   self._name=name
   self._acqiris=instrumentManager.getInstrument(acqiris)
   self._MWSource=instrumentManager.getInstrument(MWSource)
   self._pulse_generator=instrumentManager.getInstrument(pulse_generator)
   self._params=dict()
   self._params["acqiris"]=acqiris
   self._params["MWSource"]=MWSource
   self._frequencies=zeros(0)
   try:
     self._calibration=Datacube()
     self._calibration.setName('analyser IQ mixer Calibration')
     self._calibration.loadtxt(register.parameters()['%s Cal'% self._name])
   except:
     pass
   self._Ilist=[]
   self._Qlist=[]
   self._philist=[]

class Instr(Instrument):

  def saveState(self,name):
    d=dict()
    d['frequencies']=self._frequencies
    return d

  def restoreState(self,state):
    self._frequencies=state['frequencies']

  def addFrequency(self, f,name, useCorrection=False,bit=0):
    """
    Add a new frequency to the analyser
    """
    self._frequencies[name]=[f,True,bit]
    return True
        
  def startAnalyseFrequency(self,name):
    self._frequencies[name][1]=True
        
  def stopAnalyseFrequency(self,name):
    self._frequencies[name][1]=False

  def measureAll(self):
    return self._acqiris.frequenciesAnalyse(frequencies=self._frequencies.values)



  def analyse(self,nLoops=1,fast=False):
    """
    Acquire and analyse the frequencies previously sent and returns (components and probabilities)
    """
    
    maxBit=0
    for v in self._frequencies.values():
    	if v[2]>maxBit and v[1]:
    		maxBit=v[2]
    return self._acqiris.frequenciesAnalyse(frequencies=self._frequencies.values(),nLoops=nLoops,maxBit=maxBit,fast=fast)
    

  def measureBifurcationProbabilities(self):
    """
    Acquire, analyse the frequencies, convert it in clicks, and calculate averages values
    """
    (av,co,fr)=self.analyse()
    r=self._acqiris.multiplexedBifurcationMapAdd(co,fr)
    p=self._acqiris.convertToProbabilities(r)
   
  def clear(self):
    """
    Clear the list of frequencies to be analysed and the calibration paramaters associated
    """
    self._Ilist=[]  
    self._Qlist=[]
    self._philist=[]
    self._frequencies=dict()

  def calibrateAmplitudeAndOffset(self,f):
    """
    Only used when this pulse Analyser has to be used as real analyser, not when using it to see bifurcation
    """
    rowData=Datacube()
    for phi in arange(0,2*math.pi,math.pi/30):
      print "calibration : phi = %f deg" % (phi/math.pi*180)    
      self._pulseGenerator.clearPulse()
      self.clear()
      self._pulseGenerator.generatePulse(duration=20000, frequency=f, amplitude=0.6, DelayFromZero=0,useCalibration=True, phase=phi)
      self.addFrequency(f=f,useCorrection=False)    	
      self._pulseGenerator.sendPulse()
      time.sleep(0.5)
      (av, co, fr)= self.analyse()
      rowData.set(I=av[0,0], Q=av[1,0],phi=phi)          	
      rowData.commit()
    #I0=2/ptp(rowData['I'])
    #Q0=2/ptp(rowData['Q'])
    (I,Q,phi,dphi)=scipy.optimize.fmin_powell(lambda (I,Q,phi0,dphi): sum((I*rowData['I'] - sin(rowData['phi']+phi0+dphi))**2)+sum((Q*rowData['Q'] - cos(rowData['phi']+phi0))**2),(1,1,0,0))
    print (I,Q,phi,dphi)
    f_c=self._MWSource.frequency()
    df=f-f_c
    index=self._calibration.search(f_sb=df,f_c=f_c)
    if index!=None:
      self._calibration.removeRow(index)
    self._calibration.set(I=I,Q=Q,phi=dphi,f_c=f_c,f_sb=df)
    self._calibration.commit()
    self._calibration.savetxt()
    register['%s Cal'% self._name]=self._calibration.filename()
    return rowData
      
  def parameters(self):    
    """
    Returns intrument parameters
    """
    return self._params
  
  def initCal(self):
    """
    Re-init the calibration when using this instrument as real analyser
    """
    self._calibration=Datacube()
    self._calibration.setName('analyser IQ mixer Calibration')
#.........这里部分代码省略.........

def rabi02(qubit,durations,data = None,variable ="p1x",f_sb = -0.1,amplitude = 1.0,averaging = 20,delay = 0,callback = None):

  from instruments.qubit import PulseSequence

  if data == None:
    data = Datacube()
  data.setParameters(instrumentManager.parameters())
  data.setName("Rabi Sequence 02 - %s" % qubit.name())
  f_sb_12 = f_sb-qubit.parameters()["frequencies.f02"]+qubit.parameters()["frequencies.f01"]*2
  qubit.setDriveFrequency(qubit.parameters()["frequencies.f01"]+f_sb)
  qubit.setDriveAmplitude(I = amplitude,Q = amplitude)
  pi12Length = len(qubit.generateRabiPulse(phase = qubit.parameters()["pulses.xy.t_pi12"],f_sb = f_sb))
  try:
    for duration in durations:
      pulseLength = len(qubit.generateRabiPulse(length = duration,f_sb = f_sb_12))
      seq = PulseSequence()
      seq.addPulse(qubit.generateRabiPulse(length = duration,delay = qubit.parameters()["timing.readout"]-pulseLength-pi12Length,f_sb = f_sb))
      seq.addPulse(qubit.generateRabiPulse(length = qubit.parameters()["pulses.xy.t_pi12"],delay = qubit.parameters()["timing.readout"]-pi12Length,f_sb = f_sb_12))
      qubit.loadWaveform(seq.getWaveform(),readout = qubit.parameters()["timing.readout"])
      if callback != None:
        callback(duration)
      acqiris.bifurcationMap(ntimes = averaging)
      data.set(duration = duration)
      data.set(**acqiris.Psw())
      data.commit()
  finally:
    data.savetxt()
  return data

def sCurves(qubit,jba,variable = "p1x",data = None,ntimes = 20,optimize = "v20"):
  """
  Measures the s curves of the JBA. Assumes that the qubit is alread preset to a pi-pulse.
  """
  def getVoltageBounds(v0,jba,variable,ntimes):
    v = v0
    jba.setVoltage(v)
    acqiris.bifurcationMap(ntimes = ntimes)
    p = acqiris.Psw()[variable]
    
    while p > 0.03 and v < v0*2.0:
      v*=1.05
      jba.setVoltage(v)
      acqiris.bifurcationMap()
      p = acqiris.Psw()[variable]
    vmax = v
    
    v = v0
    jba.setVoltage(v)
    acqiris.bifurcationMap(ntimes = ntimes)
    p = acqiris.Psw()[variable]
    
    while p < 0.98 and v > v0/2.0:
      v/=1.05
      jba.setVoltage(v)
      acqiris.bifurcationMap()
      p = acqiris.Psw()[variable]
    vmin = v
    return (vmin*0.95,vmax*1.2)

  try:
    v0 = jba.voltage()
    if data == None:
    	sData = Datacube()
    else:
    	sData = data
    sData.setName("S curves - %s" % qubit.name())
    sData.setParameters(instrumentManager.parameters())
    s0 = Datacube("S0")
    s1 = Datacube("S1")
    s2 = Datacube("S2")
    
    sData.addChild(s0)
    sData.addChild(s1)
    sData.addChild(s2)
    
    qubit.turnOffDrive()
    
    (vmin,vmax) = getVoltageBounds(v0,jba,variable,ntimes)
    measureSingleS(voltages = arange(vmin,vmax,0.005),data = s0,jba = jba,ntimes = ntimes)

    qubit.turnOnDrive()

    loadPi01Pulse(qubit)
    measureSingleS(voltages = arange(vmin,vmax,0.005),data = s1,jba = jba,ntimes = ntimes)

    failed12 = False

    try:
      loadPi012Pulse(qubit)
      measureSingleS(voltages = arange(vmin,vmax,0.005),data = s2,jba = jba,ntimes = ntimes)  
      s1.createColumn("contrast20",s2.column(variable)-s0.column(variable))
      s1.createColumn("contrast21",s2.column(variable)-s1.column(variable))
      qubit.parameters()["readout.v20"] = s1.column("v")[argmax(s1.column("contrast20"))]
      qubit.parameters()["readout.v21"] = s1.column("v")[argmax(s1.column("contrast21"))]
    except:
      failed12 = True
      raise
      
    s1.createColumn("contrast10",s1.column(variable)-s0.column(variable))
    qubit.parameters()["readout.v10"] = s1.column("v")[argmax(s1.column("contrast10"))]
    
    if optimize == "v20" and not failed12:
      imax = argmax(s1.column("contrast20"))
      qubit.parameters()["readout.p11"] = s2.column(variable)[imax]
      v0 = s1.column("v")[imax]
    elif optimize == "v21" and not failed12:
      imax = argmax(s1.column("contrast21"))
      v0 = s1.column("v")[imax]
    else:
      imax = argmax(s1.column("contrast10"))
      qubit.parameters()["readout.p11"] = s1.column(variable)[imax]
      v0 = s1.column("v")[imax]
    
    #To do: Add dialog to ask to which voltage (v10,v20,v21) in 
    
    qubit.parameters()["readout.p00"] = 1.0-s0.column(variable)[imax]

    return (sData,v0)
    
  finally: 
    jba.setVoltage(v0)
    data.savetxt()