Added B(z) for two superconducting layers insulated by a barrier layer to libTFitPofB

src/external/TFitPofB-lib/classes/TBofZCalc.cpp

@@ -52,7 +52,7 @@ void TBofZCalc::Calculate()
 
 #pragma omp parallel for default(shared) private(j,ZZ) schedule(dynamic)
   for (j=0; j<fSteps; j++) {
-    ZZ = fParam[1] + (double)j*fDZ;
+    ZZ = fParam[1] + static_cast<double>(j)*fDZ;
     fZ[j] = ZZ;
     fBZ[j] = GetBofZ(ZZ);
   }
@@ -227,7 +227,7 @@ void TLondon1D_1L::SetBmin()
   double minZ;
   // check if the minimum is in the first layer
   minZ=-0.5*fParam[3]*log(b_a);
-  if (minZ > fParam[1] && minZ <= fParam[2]) {
+  if (minZ > fParam[1] && minZ <= fParam[1]+fParam[2]) {
     fMinZ = minZ;
     fMinB = GetBofZ(minZ);
     return;
@@ -1088,6 +1088,120 @@ vector< pair<double, double> > TLondon1D_3LS::GetInverseAndDerivative(double BB)
   return inv;
 }
 
+//------------------
+// Constructor of the TLondon1D_3LwInsulator class
+// 1D-London screening in a two thin superconducting layers fully insulated by a buffer layer
+// Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], lambda1[nm], lambda2[nm] 
+//------------------
+
+TLondon1D_3LwInsulator::TLondon1D_3LwInsulator(const vector<double> &param, unsigned int steps)
+{
+  fSteps = steps;
+  fDZ = (param[2]+param[3]+param[4])/static_cast<double>(steps);
+  fParam = param;
+  fMinZ = -1.0;
+  fMinB = -1.0;
+
+// thicknesses have to be greater or equal to zero
+  for(unsigned int i(1); i<5; i++) {
+    if(param[i] < 0.){
+      fParam[i] = 0.;
+    }
+  }
+
+// lambdas have to be greater than zero
+  for(unsigned int i(5); i<7; i++) {
+    if(param[i] < 0.1){
+      fParam[i] = 0.1;
+    }
+  }
+
+// Calculate the coefficients of the exponentials
+  double N0(fParam[0]/(1.0+exp(fParam[2]/fParam[5])));
+  double N1(fParam[0]/(1.0+exp(fParam[4]/fParam[6])));
+
+  fCoeff[0]=N0*exp((fParam[1]+fParam[2])/fParam[5]);
+  fCoeff[1]=N0*exp(-fParam[1]/fParam[5]);
+  fCoeff[2]=N1*exp((fParam[1]+fParam[2]+fParam[3]+fParam[4])/fParam[6]);
+  fCoeff[3]=N1*exp(-(fParam[1]+fParam[2]+fParam[3])/fParam[6]);
+
+// none of the coefficients should be zero
+  for(unsigned int i(0); i<4; i++)
+    assert(fCoeff[i]);
+
+  SetBmin();
+}
+
+double TLondon1D_3LwInsulator::GetBofZ(double ZZ) const
+{
+  if(ZZ < 0. || ZZ < fParam[1] || ZZ > fParam[1]+fParam[2]+fParam[3]+fParam[4] || \
+   (ZZ > fParam[1]+fParam[2] && ZZ < fParam[1]+fParam[2]+fParam[3]))
+    return fParam[0];
+
+  if(ZZ <= fParam[1]+fParam[2])
+    return fCoeff[0]*exp(-ZZ/fParam[5])+fCoeff[1]*exp(ZZ/fParam[5]);
+  else
+    return fCoeff[2]*exp(-ZZ/fParam[6])+fCoeff[3]*exp(ZZ/fParam[6]);
+}
+
+double TLondon1D_3LwInsulator::GetBmax() const
+{
+  // return applied field
+  return fParam[0];
+}
+
+double TLondon1D_3LwInsulator::GetBmin() const
+{
+  // return field minimum
+  return fMinB;
+}
+
+void TLondon1D_3LwInsulator::SetBmin()
+{
+  double b_a(fCoeff[1]/fCoeff[0]);
+  double d_c(fCoeff[3]/fCoeff[2]);
+
+  if(b_a<1E-7) {
+    b_a = 1E-7;
+  }
+  if(d_c<1E-7) {
+    d_c = 1E-7;
+  }
+
+  double minZ1, minZ2, temp;
+  // check if the minimum is in the first or third layer
+  minZ1=-0.5*fParam[5]*log(b_a);
+  minZ2=-0.5*fParam[6]*log(d_c);
+
+  if (minZ1 > fParam[1] && minZ1 <= fParam[1]+fParam[2]) {
+    fMinZ = minZ1;
+    fMinB = GetBofZ(minZ1);
+    if (minZ2 > fParam[1]+fParam[2]+fParam[3] && minZ2 <= fParam[1]+fParam[2]+fParam[3]+fParam[4]) {
+      temp = GetBofZ(minZ2);
+      if (temp < fMinB) {
+        fMinZ = minZ2;
+        fMinB = temp;
+      }
+    }
+    return;
+  } else if (minZ2 > fParam[1]+fParam[2]+fParam[3] && minZ2 <= fParam[1]+fParam[2]+fParam[3]+fParam[4]) {
+    fMinZ = minZ2;
+    fMinB = GetBofZ(minZ2);
+    return;
+  }
+
+//  assert(fMinZ > 0. && fMinB > 0.);
+  if(fMinZ <= 0.){
+    fMinZ = 0.;
+  }
+  if(fMinB <= 0.){
+    fMinB = 0.;
+  }
+
+  return;
+}
+
+
 
 // //------------------
 // // Constructor of the TLondon1D_4L class

src/external/TFitPofB-lib/include/TBofZCalc.h

@@ -239,16 +239,25 @@ private:
   double fCoeff[6];
 };
 
-// //--------------------
+//--------------------
-// // Class "for Meissner screening" in a thin superconducting film - four layers with four different lambdas
+// Class "for Meissner screening" in a thin superconducting film - tri-layer with insulating buffer layer, two lambda
-// //--------------------
+//--------------------
-// 
+
-// class TLondon1D_4L : public TBofZCalc {
+class TLondon1D_3LwInsulator : public TBofZCalc {
-// 
+
-// public:
+public:
-// 
+
-//   TLondon1D_4L(unsigned int, const vector<double>& );
+  TLondon1D_3LwInsulator(const vector<double>&, unsigned int steps = 3000);
-// 
+  double GetBofZ(double) const;
-// };
+  double GetBmin() const;
+  double GetBmax() const;
+
+private:
+  void SetBmin();
+
+  double fMinZ;
+  double fMinB;
+  double fCoeff[4];
+};
 
 #endif // _BofZCalc_H_

src/external/libCalcMeanFieldsLEM/TCalcMeanFieldsLEM.cpp

@@ -425,3 +425,107 @@ double TMeanFieldsForScTrilayer::CalcMeanB (double E, const vector<double>& inte
     }
     return meanB;
 }
+
+// Constructor: Read the energies from the xml-file and load the according TRIM.SP-data files
+
+TMeanFieldsForScTrilayerWithInsulator::TMeanFieldsForScTrilayerWithInsulator() {
+
+  // read startup file
+  string startup_path_name("TFitPofB_startup.xml");
+
+  TSAXParser *saxParser = new TSAXParser();
+  TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler();
+  saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler);
+  int status (saxParser->ParseFile(startup_path_name.c_str()));
+  // check for parse errors
+  if (status) { // error
+    cout << endl << "**WARNING** reading/parsing TFitPofB_startup.xml failed." << endl;
+  }
+
+  string rge_path(startupHandler->GetDataPath());
+  vector< pair<double, string> > energy_vec(startupHandler->GetEnergies());
+
+  TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+  fImpProfile = x;
+  x = 0;
+
+  // clean up
+  if (saxParser) {
+    delete saxParser;
+    saxParser = 0;
+  }
+  if (startupHandler) {
+    delete startupHandler;
+    startupHandler = 0;
+  }
+
+}
+
+// Operator-method that returns the mean field for a given implantation energy
+// Parameters: field, deadlayer, layer1, layer2, layer3, lambda1, lambda2, weight1, weight2, weight3, weight4, weight5
+
+double TMeanFieldsForScTrilayerWithInsulator::operator()(double E, const vector<double> &par_vec) const{
+  
+  vector<double> interfaces;
+  interfaces.push_back(par_vec[1]);
+  interfaces.push_back(par_vec[1]+par_vec[2]);
+  interfaces.push_back(par_vec[1]+par_vec[2]+par_vec[3]);
+  interfaces.push_back(par_vec[1]+par_vec[2]+par_vec[3]+par_vec[4]);
+  
+  vector<double> weights;
+  weights.push_back(par_vec[7]);
+  weights.push_back(par_vec[8]);
+  weights.push_back(par_vec[9]);
+  weights.push_back(par_vec[10]);
+  weights.push_back(par_vec[11]);
+  
+  // Calculate field profile
+  vector<double> parForBofZ;
+  for (unsigned int i(0); i<7; i++)
+    parForBofZ.push_back(par_vec[i]);
+
+  TLondon1D_3LwInsulator BofZ(parForBofZ);
+
+  vector<double> energies(fImpProfile->Energy());
+  vector<double>::const_iterator energyIter;
+  energyIter = find(energies.begin(), energies.end(), E);
+
+  if (energyIter != energies.end()) { // implantation profile found - no interpolation needed
+    return CalcMeanB(E, interfaces, weights, BofZ);
+  } else {
+    if (E < *energies.begin())
+      return CalcMeanB(*energies.begin(), interfaces, weights, BofZ);
+    if (E > *(energies.end()-1))
+      return CalcMeanB(*(energies.end()-1), interfaces, weights, BofZ);
+
+    energyIter = find_if(energies.begin(), energies.end(), bind2nd( greater<double>(), E));
+//    cout << *(energyIter - 1) << " " << *(energyIter) << endl;
+
+    double E1(*(energyIter - 1));
+    double E2(*(energyIter));
+
+    double B1(CalcMeanB(E1, interfaces, weights, BofZ));
+    double B2(CalcMeanB(E2, interfaces, weights, BofZ));
+
+    return B1 + (B2-B1)/(E2-E1)*(E-E1);
+  }
+}
+
+double TMeanFieldsForScTrilayerWithInsulator::CalcMeanB
+ (double E, const vector<double>& interfaces, const vector<double>& weights, const TLondon1D_3LwInsulator& BofZ) const {
+  //calcData->UseHighResolution(E);
+    fImpProfile->WeightLayers(E, interfaces, weights);
+    fImpProfile->Normalize(E);
+
+    vector<double> z(fImpProfile->DataZ(E));
+    vector<double> nz(fImpProfile->DataNZ(E));
+
+    // calculate mean field
+
+    double meanB(0.);
+
+    for (unsigned int i(0); i<z.size(); i++) {
+      meanB += (z[1]-z[0])*nz[i]*BofZ.GetBofZ(0.1*z[i]);
+    }
+    return meanB;
+}

src/external/libCalcMeanFieldsLEM/TCalcMeanFieldsLEM.h

@@ -98,4 +98,20 @@ private:
   ClassDef(TMeanFieldsForScTrilayer,1)
 };
 
+class TMeanFieldsForScTrilayerWithInsulator : public PUserFcnBase {
+
+public:
+  // default constructor
+  TMeanFieldsForScTrilayerWithInsulator();
+  ~TMeanFieldsForScTrilayerWithInsulator() {delete fImpProfile; fImpProfile = 0;}
+
+  double operator()(double, const vector<double>&) const;
+  double CalcMeanB (double, const vector<double>&, const vector<double>&, const TLondon1D_3LwInsulator&) const;
+
+private:
+  TTrimSPData *fImpProfile;
+
+  ClassDef(TMeanFieldsForScTrilayerWithInsulator,1)
+};
+
 #endif /* _TCalcMeanFieldsLEM_H_ */

src/external/libCalcMeanFieldsLEM/TCalcMeanFieldsLEMLinkDef.h

@@ -40,6 +40,7 @@
 #pragma link C++ class TMeanFieldsForScSingleLayer+;
 #pragma link C++ class TMeanFieldsForScBilayer+;
 #pragma link C++ class TMeanFieldsForScTrilayer+;
+#pragma link C++ class TMeanFieldsForScTrilayerWithInsulator+;
 
 #endif //__CINT__
 // root dictionary stuff --------------------------------------------------