Thursday, December 23, 2010

I know you don't care...

I know you don't care, but I thought I'd post what I've been working on for the past few days. I guess I've been working on it longer than that, really. This is the code I use to turn the output of my machine into something useful. I use a piece of software called ROOT. You can read about it here. It's used mostly by particle physicists....It's runs off of c++, a programing language I didn't know 2 1/2 weeks ago. I still don't know it, but I'm like the foreigner who knows just enough to figure out what's going on and ask questions like, "Where I see emaciation proclamation?" I can use big words, but not necessarily the right way.

...it looks like blogger is trying to read my code as if it were html even though I entered it in this text portion. Oh well...


#include
#include
#include
#include
#include
#include
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TH2.h"
#include "TH1.h"
#include "TGraph.h"
#include "TMatrixDSym.h"
#include "TF1.h"
#include "TFitResult.h"

using namespace std;

int main()
{
string s1, filename; //This string portion and the following are used to generate the names of the files
string filebase = "Test 0";
string digit1 = "0";
string digit2 = "0";
string suffix = ".csv";
string treebase = "myTree";
char *cstr; //This little pointer is used for converting things from strings to digits
TGraph *DispVsLoad, *DispVsStiffness, *DispVsModulus, *DispVsHardness;
ostringstream ss; //This is used in converting ints to strings or chars

int i, j, k, m, loop, outerloop, cap, counter = 0; //my counters
double displacement = 0;
double load = 0;
double time = 0;
double stiffness = 0;
double hardness = 0;
double modulus = 0;
double modmean, moddev, modvar, modM2, moddelta;
double modpmean[4], modpdev[4], modpvar[4], modpM2[4], modpdelta[4];
double hardmean, harddev, hardvar, hardM2, harddelta;
double hardpmean[4], hardpdev[4], hardpvar[4], hardpM2[4], hardpdelta[4];
double loadmean[4], loaddev[4], loadvar[4], loadM2[4], loaddelta[4];
double stiffmean[4], stiffdev[4], stiffvar[4], stiffM2[4], stiffdelta[4];

for (m = 0; m <4;>
modpmean[m] = modpdev[m] = modpvar[m] = modpM2[m] = modpdelta[m] = 0;
hardpmean[m] = hardpdev[m] = hardpvar[m] = hardpM2[m] = hardpdelta[m] = 0;
loadmean[m] = loaddev[m] = loadvar[m] = loadM2[m] = loaddelta[m] = 0;
stiffmean[m] = stiffdev[m] = stiffvar[m] = stiffM2[m] = stiffdelta[m] = 0;
}

ifstream DataFile;
ofstream OFDispVsLoad, OFDispVsStiffness, OFDispVsModulus, OFDispVsHardness;

TF1 *fa = new TF1("fa","[0]*x*x*x + [1]*x*x+[2]*x+[3]",-100,20000);
fa->SetParameter(0,1);
fa->SetParameter(1,1);
fa->SetParameter(2,1);
fa->SetParameter(3,1);
TFile * myRootFile = new TFile("Indenter.root", "RECREATE", "A File with nano indenter output");
TTree *myTree = new TTree("myTree", "nano indenter TTree");
TTree *myTreeTotal = new TTree("myTreeTotal", "nano indenter TTree");
myTreeTotal->Branch("displacement", &displacement, "displacement/D");
myTreeTotal->Branch("load", &load, "load/D");
myTreeTotal->Branch("time", &time, "time/D");
myTreeTotal->Branch("stiffness", &stiffness, "stiffness/D");
myTreeTotal->Branch("hardness", &hardness, "hardness/D");
myTreeTotal->Branch("modulus", &modulus, "modulus/D");
OFDispVsLoad.open("OFDispVsLoad.txt");
OFDispVsLoad<<"Test \tChi2 \tNDF \tPar0 \t\t+/-Error \tPar1 \t\t+/-Error \tPar2 \t\t+/-Error \tPar3 \t\t+/-Error"<
OFDispVsStiffness.open("OFDispVsStiffness.txt");
OFDispVsStiffness<<"Test \tChi2 \tNDF \tPar0 \t\t+/-Error \tPar1 \t\t+/-Error \tPar2 \t\t+/-Error \tPar3 \t\t+/-Error"<
OFDispVsModulus.open("OFDispVsModulus.txt");
OFDispVsModulus<<"Test \tChi2 \tNDF \tPar0 \t\t+/-Error \tPar1 \t\t+/-Error \tPar2 \t\t+/-Error \tPar3 \t\t+/-Error \t 9k-10k avg \t stddev"<
OFDispVsHardness.open("OFDispVsHardness.txt");
OFDispVsHardness<<"Test \tChi2 \tNDF \tPar0 \t\t+/-Error \tPar1 \t\t+/-Error \tPar2 \t\t+/-Error \tPar3 \t\t+/-Error \t 9k-10k avg \t stddev"<

cap = 10; //this is to run the file reading loop from 001 to 025

for (outerloop = 0; outerloop <3; class="Apple-tab-span" style="white-space:pre"> //This is where I enter the file reading loop, first I establish the first digit of the file which matches outerloop
{
if (outerloop == 0) {
loop = 1;
}
else {
loop = 0;
}
if (outerloop ==2){
cap = 5;
}

ss.str("");
ss <<>
digit1 = ss.str();

for (loop = loop; loop < class="Apple-tab-span" style="white-space:pre"> //Here is where I get the second digit of the file name
{
ss.str("");
ss <<>
digit2 = ss.str();
filename = treebase + digit1 + digit2; //Here I use the same first and second digit to name the TTree
cstr = new char [filename.size() + 1]; //point at the address of a chunk the size of 1+ the data just read in
strcpy (cstr, filename.c_str());
DispVsLoad = new TGraph(); //Set up the graphs I want to match curves to
DispVsStiffness = new TGraph();
DispVsModulus = new TGraph();
DispVsHardness = new TGraph();
/* TH2D *myModDis = new TH2D("myModDis","Modulus versus Displacement",17000,-3000.,14000.,1000,0.0,8.0);
TH1D *myWeighting = new TH1D("myWeighting","Weighting Entries",17000,-3000.,14000.);*/
TTree *myTree = new TTree(cstr, "nano indenter TTree"); //Open a TTree with the treename using the last 2 digits of the filename
delete cstr;
myTree->Branch("displacement", &displacement, "displacement/D"); //Put branches on the TTree that match the output file from the nano indenter
myTree->Branch("load", &load, "load/D");
myTree->Branch("time", &time, "time/D");
myTree->Branch("stiffness", &stiffness, "stiffness/D");
myTree->Branch("hardness", &hardness, "hardness/D");
myTree->Branch("modulus", &modulus, "modulus/D");
filename = filebase + digit1 + digit2 + suffix; //Generate the filename of the files to be opened
cout<< "The file name is: "<
DataFile.open(filename);
//First I need to ignore the 12 headings in the file
for (i = 1; i<13>
getline(DataFile, s1, ',');
}
j = 1; //used to keep only the data in the range that represents when the tip is in contact and measuring
i = 0; //used for the switch to put each item in the column in the right branch and right graph
k = 0; //used for averaging regions, like the machine does
modmean = moddev = modvar = modM2 = moddelta =0;
hardmean = harddev = hardvar = hardM2 = moddelta = 0;

//Now I read all the data
while (DataFile.good()){
getline(DataFile, s1, ','); //Read the file to the next ',' save it as a string in s1
cstr = new char [s1.size() + 1]; //point at the address of a chunk the size of 1+ the data just read in
strcpy (cstr, s1.c_str()); //take the string and turn it into characters, put it in the location cstr is pointing at
double f = 0.0; // (default value is 0.0)
f = atof(cstr); //convert the characters into a double
delete cstr; //free the memory
switch(i%6){ //use the switch to assign each item from the file to the right variables
case 0:
displacement = f;
break;
case 1:
load = f;
break;
case 2:
time = f;
break;
case 3:
stiffness = f;
break;
case 4:
hardness = f;
break;
case 5:
modulus = f;
if ((modulus > 0.5) & (modulus <>
{

DispVsLoad->SetPoint(j,displacement,load); //Fill in the TGraphs one point at a time
DispVsStiffness->SetPoint(j,displacement,stiffness);
DispVsHardness->SetPoint(j,displacement,hardness);
DispVsModulus->SetPoint(j,displacement,modulus);
myTree->Fill();
myTreeTotal->Fill();
// myModDis->Fill(displacement,modulus);
// myWeighting->Fill(displacement,modulus);
j++;
}

if((displacement > 9000) & (displacement < class="Apple-tab-span" style="white-space:pre"> //This is the region I want to average modulus and hardness
k++; //this is an algorithm from wikipedia for on-line variance
moddelta = modulus - modmean;
modmean = modmean + moddelta/k;
modM2 = modM2 + moddelta*(modulus - modmean);
if (k!=1) modvar = modM2/(k-1);

harddelta = hardness - hardmean;
hardmean = hardmean + harddelta/k;
hardM2 = hardM2 + harddelta*(hardness - hardmean);
if (k!=1) hardvar = hardM2/(k-1);
}
break;
default:
break;
}
i++;
}

moddev = sqrt(modvar);
harddev = sqrt(hardvar);

//A copy of each set of commands for each of load, stiffness, hardness, and modulus
//These use the filled graphs to get fit data and write that data to the file
{
//DispVsLoad->Draw("ACP");
TFitResultPtr r = DispVsLoad->Fit("fa", "S");
TMatrixDSym cov = r->GetCovarianceMatrix();
double chi2 = r->Chi2();
int ndf = r->Ndf();
double par[4], err[4];
for (int m = 0; m <>
par[m] = fa->GetParameter(m);
err[m] = fa->GetParError(m);
}

if (chi2 != -1){
counter++;
for (int m = 0; m < class="Apple-tab-span" style="white-space:pre">
loaddelta[m] = par[m] - loadmean[m];
loadmean[m] = loadmean[m] + loaddelta[m]/counter;
loadM2[m] = loadM2[m] + loaddelta[m] * (par[m] - loadmean[m]);
if (counter!=1) loadvar[m] = loadM2[m]/(counter-1);
}
}
// r->Write();
OFDispVsLoad <
// DispVsLoad->Write();
delete DispVsLoad;
}

{
//DispVsStiffness->Draw("ACP");
DispVsStiffness->Write();
TFitResultPtr r = DispVsStiffness->Fit("fa", "S");
TMatrixDSym cov = r->GetCovarianceMatrix();
double chi2 = r->Chi2();
int ndf = r->Ndf();
double par[4], err[4];
for (int m = 0; m <>
par[m] = fa->GetParameter(m);
err[m] = fa->GetParError(m);
}

if (chi2 != -1){
for (int m = 0; m <>
stiffdelta[m] = par[m] - stiffmean[m];
stiffmean[m] = stiffmean[m] + stiffdelta[m]/counter;
stiffM2[m] = stiffM2[m] + stiffdelta[m] * (par[m] - stiffmean[m]);
if (counter!=1) stiffvar[m] = stiffM2[m]/(counter-1);
}
}
// r->Write();
OFDispVsStiffness <
// DispVsStiffness->Write();
delete DispVsStiffness;
}

{
TFitResultPtr r = DispVsHardness->Fit("fa", "S");
TMatrixDSym cov = r->GetCovarianceMatrix();
double chi2 = r->Chi2();
int ndf = r->Ndf();
double par[4], err[4];
for (int m = 0; m <>
par[m] = fa->GetParameter(m);
err[m] = fa->GetParError(m);
}

if (chi2 != -1){
for (int m = 0; m <>
hardpdelta[m] = par[m] - hardpmean[m];
hardpmean[m] = hardpmean[m] + hardpdelta[m]/counter;
hardpM2[m] = hardpM2[m] + hardpdelta[m] * (par[m] - hardpmean[m]);
if (counter!=1) hardpvar[m] = hardpM2[m]/(counter-1);
}
}
// r->Write();
OFDispVsHardness <
// DispVsHardness->Write();
delete DispVsHardness;
}

{
TFitResultPtr r = DispVsModulus->Fit("fa", "S");
TMatrixDSym cov = r->GetCovarianceMatrix();
double chi2 = r->Chi2();
int ndf = r->Ndf();
double par[4], err[4];
for (int m = 0; m <>
par[m] = fa->GetParameter(m);
err[m] = fa->GetParError(m);
}

if (chi2 != -1){
for (int m = 0; m <>
modpdelta[m] = par[m] - modpmean[m];
modpmean[m] = modpmean[m] + modpdelta[m]/counter;
modpM2[m] = modpM2[m] + modpdelta[m] * (par[m] - modpmean[m]);
if (counter!=1) modpvar[m] = modpM2[m]/(counter-1);
}
}
// r->Write();
OFDispVsModulus <
// DispVsModulus->Write();
delete DispVsModulus;
}
/*myWeighting->Write();
myModDis->Write();*/
DataFile.close();
myTree->Write();
delete myTree;
}
} //end of activities for each data file

OFDispVsLoad << "Average & Standard Dev.\t";
OFDispVsStiffness << "Average & Standard Dev.\t";
OFDispVsModulus << "Average & Standard Dev.\t";
OFDispVsHardness << "Average & Standard Dev.\t";

for (int m = 0; m<4;>
OFDispVsLoad <<>
OFDispVsStiffness <<>
OFDispVsModulus <<>
OFDispVsHardness <<>
}
myTreeTotal->Write();
myRootFile->Close();
return 0;
}

4 comments:

Jess and Jen said...

I use Python a lot at work; it's pretty easy and straightforward. C++ isn't something I've ever learned, although I took classes in machine language, C, Java, and Visual Basic in college. There's a reason I'm not a computer scientist.

Having said that, I have a real fondness for trying to automate things using scripting like Python, or AML (Arc Macro Language -- it's GIS stuff). But I leave the C++ to someone else. It looks like I can call you now and have you do it for me...

Good luck in your super intensive few weeks of getting the thesis done. Sorry it fell at such a bad time! -Jess

chelsey said...

Huh? Way over my head. I'm not exactly computer program savvy. No clue. Good luck.

Michelle said...

Ok, I don't get it at all. But I'm sure it makes sense to those it should make sense to.

The Duke said...

Are you kidding?!