Peakfit.C

{
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.
}
}