From GlueXWiki
{
gROOT->Reset();
char adcn_sim[50],file_sim[50],file_jlab[50],histN[50];
float adcmean_sim[18];
int t=0;
char drawthis[1000] = "(sqrt((";
char drawthat[500] = " ";
char adccalib[1100] = "Ephoton/";
char adccalib2[20];
char cut1[200];
char adcn[50];
char S1[20],S2[20],S3[20],S4[20],S5[20],S6[20],S7[20],S8[20],S9[20],S10[20],S11[20],S12[20],S13[20],S14[20],S15[20],S16[20],S17[20],S18[20];
char N1[20],N2[20],N3[20],N4[20],N5[20],N6[20],N7[20],N8[20],N9[20],N10[20],N11[20],N12[20],N13[20],N14[20],N15[20],N16[20],N17[20],N18[20];
float adcmeann[18],adcmeans[18],n8ratioN[18],n8ratioS[18],simsclfctrN[18],simsclfctrS[18],adcrmsN[18];
Double_t fitpar[9],chisqr[9],fitndf[9],what;
int k=0;
int m=0;
float gaincor[18],gaincor2[18],epsilon,lastsmall;
ofstream outputfile;
outputfile.open("peakfit.txt"); // text file where the means are output
//// root file
TFile *jlabfile = new TFile("/home/leverin/work/root/dst/bcal_dst02334.root");
//for( int i = 0; i<= 17; i++){
// gaincor[i] = 1.0;
// gaincor2[i] = 1.0;
//}
fitpar[2] = 1.0;
//// The correction factors from the last run. Better place to start from.
gaincor[0] = 0.97;
gaincor[1] = 0.99;
gaincor[2] = 1.02;
gaincor[3] = 0.98;
gaincor[4] = 0.99;
gaincor[5] = 0.99;
gaincor[6] = 0.99;
gaincor[7] = 0.79;
gaincor[8] = 1.05;
gaincor[9] = 0.99;
gaincor[10] = 1.03;
gaincor[11] = 0.99;
gaincor[12] = 1.04;
gaincor[13] = 0.94;
gaincor[14] = 0.95;
gaincor[15] = 1.02;
gaincor[16] = 1.02;
gaincor[17] = 0.99;
gaincor2[0] = 1.02;
gaincor2[1] = 1.02;
gaincor2[2] = 1.02;
gaincor2[3] = 0.99;
gaincor2[4] = 0.98;
gaincor2[5] = 1.02;
gaincor2[6] = 0.98;
gaincor2[7] = 0.68;
gaincor2[8] = 1.04;
gaincor2[9] = 1.11;
gaincor2[10] = 0.97;
gaincor2[11] = 1.05;
gaincor2[12] = 0.99;
gaincor2[13] = 1.04;
gaincor2[14] = 0.99;
gaincor2[15] = 1.03;
gaincor2[16] = 1.03;
gaincor2[17] = 1.02;
for (int j = 1; j<=2; j++){
for( int i = 0; i <= 17; i++ ){
//// initialize some parameters
k = 0;
epsilon = 0.01; // step size for each iteration
fitpar[2] = 1.0;
if (j==1){lastsmall = gaincor[i];}
if (j==2){lastsmall = gaincor2[i];}
if (gaincor[i]<1.0){epsilon=-0.01;} //should go faster if you start in the same direction as last time
if (gaincor2[i]<1.0){epsilon=-0.01;}
//// START of fitting loop
do {
k += 1;
fitpar[5] = fitpar[2];
//// define histograms
TH1F *hist = new TH1F( "hist","z", 400,-1.5,1.5);
TH1F *hist2 = new TH1F( "hist2","Ephoton/GeomMean", 200,0.0,0.4);
//// put together the argument to Draw
sprintf( S1, "adcS.s1*%f",gaincor[0]);strcat(drawthis,S1);
sprintf( S2, "+adcS.s2*%f",gaincor[1]);strcat(drawthis,S2);
sprintf( S3, "+adcS.s3*%f",gaincor[2]);strcat(drawthis,S3);
sprintf( S4, "+adcS.s4*%f",gaincor[3]);strcat(drawthis,S4);
sprintf( S5, "+adcS.s5*%f",gaincor[4]);strcat(drawthis,S5);
sprintf( S6, "+adcS.s6*%f",gaincor[5]);strcat(drawthis,S6);
sprintf( S7, "+adcS.s7*%f",gaincor[6]);strcat(drawthis,S7);
sprintf( S8, "+adcS.s8*%f",gaincor[7]);strcat(drawthis,S8);
sprintf( S9, "+adcS.s9*%f",gaincor[8]);strcat(drawthis,S9);
sprintf( S10, "+adcS.s10*%f",gaincor[9]);strcat(drawthis,S10);
sprintf( S11, "+adcS.s11*%f",gaincor[10]);strcat(drawthis,S11);
sprintf( S12, "+adcS.s12*%f",gaincor[11]);strcat(drawthis,S12);
sprintf( S13, "+adcS.s13*%f",gaincor[12]);strcat(drawthis,S13);
sprintf( S14, "+adcS.s14*%f",gaincor[13]);strcat(drawthis,S14);
sprintf( S15, "+adcS.s15*%f",gaincor[14]);strcat(drawthis,S15);
sprintf( S16, "+adcS.s16*%f",gaincor[15]);strcat(drawthis,S16);
sprintf( S17, "+adcS.s17*%f",gaincor[16]);strcat(drawthis,S17);
sprintf( S18, "+adcS.s18*%f",gaincor[17]);strcat(drawthis,S18);
sprintf( N1, "adcN.n1*%f",gaincor2[0]);strcat(drawthat,N1);
sprintf( N2, "+adcN.n2*%f",gaincor2[1]);strcat(drawthat,N2);
sprintf( N3, "+adcN.n3*%f",gaincor2[2]);strcat(drawthat,N3);
sprintf( N4, "+adcN.n4*%f",gaincor2[3]);strcat(drawthat,N4);
sprintf( N5, "+adcN.n5*%f",gaincor2[4]);strcat(drawthat,N5);
sprintf( N6, "+adcN.n6*%f",gaincor2[5]);strcat(drawthat,N6);
sprintf( N7, "+adcN.n7*%f",gaincor2[6]);strcat(drawthat,N7);
sprintf( N8, "+adcN.n8*%f",gaincor2[7]);strcat(drawthat,N8);
sprintf( N9, "+adcN.n9*%f",gaincor2[8]);strcat(drawthat,N9);
sprintf( N10, "+adcN.n10*%f",gaincor2[9]);strcat(drawthat,N10);
sprintf( N11, "+adcN.n11*%f",gaincor2[10]);strcat(drawthat,N11);
sprintf( N12, "+adcN.n12*%f",gaincor2[11]);strcat(drawthat,N12);
sprintf( N13, "+adcN.n13*%f",gaincor2[12]);strcat(drawthat,N13);
sprintf( N14, "+adcN.n14*%f",gaincor2[13]);strcat(drawthat,N14);
sprintf( N15, "+adcN.n15*%f",gaincor2[14]);strcat(drawthat,N15);
sprintf( N16, "+adcN.n16*%f",gaincor2[15]);strcat(drawthat,N16);
sprintf( N17, "+adcN.n17*%f",gaincor2[16]);strcat(drawthat,N17);
sprintf( N18, "+adcN.n18*%f",gaincor2[17]);strcat(drawthat,N18);
strcat(drawthis,")*(");
strcat(drawthis,drawthat);
strcat(adccalib,drawthis);
strcat(adccalib,")))>>hist2");
// cout << adccalib << "xxxx" << endl;
//// cuts
sprintf(cut1 , "trig\=\=8\&\&Nphotons\=\=1\&\&Ephoton\>200\.0\&\&Ephoton\<400\.0\&\&Ntdctrig\>\=8%s","");
//// draw and find the mean for E_photon / E_bcal(geometric mean)
bcal->Draw(adccalib,cut1);
adcmeann[1] = hist2->GetMean();
adcrmsN[1] = hist2->GetRMS();
g2 = new TF1("g2","gaus",(adcmeann[1]-adcrmsN[1]),(adcmeann[1]+adcrmsN[1]));
hist2->Fit("g2","R");
fitpar[7] = g2->GetParameter(1);
//// add the mean to the "z" Draw arguement
sprintf(adccalib2, "))*%f",fitpar[7]);
strcat(drawthis,adccalib2);
strcat(drawthis,"-Ephoton)/Ephoton>>hist");
// cout << drawthis << "xxxx" << endl;
//// Draw and find the parameters for z for this iteration
bcal->Draw(drawthis,cut1);
adcmeann[0] = hist->GetMean();
adcrmsN[0] = hist->GetRMS();
g1 = new TF1("g1","gaus",(-adcrmsN[0]),(+adcrmsN[0]));
hist->Fit("g1","R");
fitpar[0] = g1->GetParameter(0); //contant
fitpar[1] = g1->GetParameter(1); //mean
fitpar[2] = g1->GetParameter(2); //sigma
chisqr[0] = g1->GetChisquare(); //Chi squared
fitndf[0] = g1->GetNDF(); //degrees of freedom
//// delete histograms and reset Draw arguement char array
hist->Delete();
hist2->Delete();
strcpy(drawthis, "(sqrt((");
strcpy(drawthat, " ");
strcpy(adccalib, "Ephoton/");
//// output paramters on the go
cout << "i="<< i <<" "<< fitpar[0] << " " << fitpar[1] << " "<< fitpar[2] << " "<< chisqr[0] << "\/" << fitndf[0] << " " << k << " " << gaincor[i] <<" "<< gaincor2[i] << endl;
//// determine if the new sigma is larger than the old one, then decide what to do
/// North gaincor
if (j==1) {
if (epsilon >0) {
if (fitpar[2] > fitpar[5]) {
if (k == 2) {
fitpar[5] = fitpar[2]; //dummy step so it does't exit
epsilon = -epsilon;
gaincor[i] += epsilon;
} // if after the first iteration sigma is larger, go the other direction
}
if (fitpar[2] < fitpar[5]) {
lastsmall = gaincor[i];
gaincor[i] += epsilon;}
}// if the sigma is smaller on this iteration keep going
if (epsilon < 0) {
if (fitpar[2] > fitpar[5]) {
if (k == 2) {
fitpar[5] = fitpar[2]; //dummy step so it doesn't exit
epsilon = -epsilon;
gaincor[i] += epsilon;
}// if after the first iteration sigma is larger, go the other direction
}
if (fitpar[2] < fitpar[5]) {
lastsmall = gaincor[i];
gaincor[i] += epsilon;}
}
if (fitpar[2] > fitpar[5] && k>2) {gaincor[i] = lastsmall;} //use the last gaincor[i] that gave the smallest sigma, go to next PMT
}
// South gaincor
if (j==2) {
if (epsilon >0) {
if (fitpar[2] > fitpar[5]) {
if (k == 2) {
fitpar[5] = fitpar[2]; //dummy step so it doesn't exit
epsilon = -epsilon;
gaincor2[i] += epsilon;
}// if after the first iteration sigma is larger, go the other direction
}
if (fitpar[2] < fitpar[5]) {
lastsmall = gaincor2[i];
gaincor2[i] += epsilon;}
}// if the sigma is smaller on this iteration keep going
if (epsilon < 0) {
if (fitpar[2] > fitpar[5]) {
if (k == 2) {
fitpar[5] = fitpar[2]; //dummy step so it doesn't exit
epsilon = -epsilon;
gaincor2[i] += epsilon;
}// if after the first iteration sigma is larger, go the other direction
}
if (fitpar[2] < fitpar[5]) {
lastsmall = gaincor2[i];
gaincor2[i] += epsilon;}
}
if (fitpar[2] > fitpar[5] && k>2) {gaincor2[i] = lastsmall;}
}
}
while (fitpar[2] <= fitpar[5]); //end the loop for this photomultiplier tube
}
}
for ( int l = 0;l<=17;l++){
outputfile << fitpar[0] << " " << fitpar[1] << " "<< fitpar[2] << " "<< chisqr[0] << "\/" << fitndf[0] << " " << gaincor[l] <<" "<<gaincor2[l] << endl;
// outputs the current gain correction factors to "peakfit.txt". Final ones are at the end of the file.
}
}
