Newer
Older
TB_Chris / Kepler / options / ResStudies / .svn / text-base / res_analysis.C.svn-base
  1. #include "TFile.h"
  2. #include "TH2.h"
  3. #include "TH1.h"
  4. #include "TF1.h"
  5. #include "TGraphErrors.h"
  6. #include "iostream"
  7. #include <stdio.h>
  8. #include "Riostream.h"
  9. #include <math.h>
  10. #include <stdlib.h>
  11. #include <string>
  12. #include <vector>
  13. #include "TRandom.h"
  14. #include "TGraph.h"
  15. #include "TMultiGraph.h"
  16. #include "TCanvas.h"
  17. #include "TMath.h"
  18. #include "TLegend.h"
  19. //script opens the kepler-histograms sent by submitRes.py to the grid, specifically taking the unbiased residual for the DUT (plane 4) for both X and Y.
  20. //Fits a guassian to each of the residuals and stores the mean and rms.
  21. //Plots them as a function of Bias or Angle, where the bias and angle values are stored in runList.txt
  22. //stores them as a .png and .root file format (in LHCb style)
  23. //run as root -b if you dont want canvas printed to screen.
  24. //MAY need to change path in ~L62 to your own!!!!!
  25. //Emma Buchanan May 2015 (emma.buchanan@cern.ch)
  26.  
  27. void res_analysis(){
  28. //-------Change Parameters Here------//
  29. const int job = 9; //job number
  30. const int run =5; //number of subjobs in this job
  31. const int beg= 13; //from runList.txt what is the range of the "Block" begining->end
  32. const int end= 17;
  33. int choice =0; //Is this an Bias or angle scan? (0 for Bias, 1 for angle)
  34. const int iplane =4; //number of the plane you wish to analyse (DUT =4)
  35.  
  36. //saving the final plot as a .png, and/or root file
  37. char plot_bias[200]; //the name of the plot, will be printed to a .png at the end of this script.
  38. sprintf(plot_bias, "%s", "sensor_bias_scan.png");
  39. char plot_angle[200];
  40. sprintf(plot_angle, "%s", "sensor_angle_scan.png");
  41. char root_plotb[100]; //this .root file will contain the final plot/s
  42. sprintf(root_plotb, "%s", "sensor_bias_scan.root");
  43. char root_plota[100]; //this .root file will contain the final plot/s
  44. sprintf(root_plota, "%s", "sensor_angle_scan.root");
  45. //----------end of parameters--------//
  46.  
  47.  
  48. //Reading kepler histograms and fitting gaussian to DUT
  49. float DUT_meanX[run]; float DUT_rmsX[run]; //fit parameters for X
  50. float DUT_meanY[run]; float DUT_rmsY[run]; //fit parameters for Y
  51. TF1 *gauss = new TF1("gauss","gaus");
  52. TH1F *XDUT[run];
  53. TH1F *YDUT[run];
  54. TCanvas *Xcanvas[run];
  55. TCanvas *Ycanvas[run];
  56. //opening the DUT (plane 4) from the kepler-histos and fitting a gaussian to each.
  57. char filename[100];
  58. char locationX[100]; char locationY[100];
  59. TFile *openFile[100];
  60. for (int subjob =0; subjob<run; subjob++){
  61. sprintf(filename, "%s%d%s%d%s","/afs/cern.ch/user/e/ebuchana/gangadir/workspace/ebuchana/LocalXML/",job,"/",subjob,"/output/Kepler-histos.root "); //the location of each root file for each run
  62. sprintf(locationX, "%s%d", "Tb/TbTrackPlots/BiasedResiduals/GlobalX/Plane",iplane); //location of the biased residual of the DUT
  63. sprintf(locationY, "%s%d", "Tb/TbTrackPlots/BiasedResiduals/GlobalY/Plane",iplane); //location of the biased residual of the DUT
  64. cout << filename << endl;
  65. openFile[subjob] = new TFile(filename); //opening the root files
  66.  
  67. //------------------ For X ------------------//
  68. cout << locationX << endl;
  69. Xcanvas[subjob] = new TCanvas(); //creating a canvas for each histogram
  70. Xcanvas[subjob]->SetTitle("X Biased Residal for DUT"); //setting canvas title
  71. XDUT[subjob] = (TH1F*)openFile[subjob]->Get(locationX); //getting each of the histograms from the kepler root files
  72. XDUT[subjob]->Fit("gauss", "R"); //fitting guassian
  73. DUT_meanX[subjob]=XDUT[subjob]->GetFunction("gauss")->GetParameter(1); //getting the mean of the postion
  74. DUT_rmsX[subjob]=XDUT[subjob]->GetFunction("gauss")->GetParameter(2); //this should be the RMS
  75. cout << "DUT mean X\t" << DUT_meanX[subjob] <<endl;
  76. cout << "DUT rms X\t" << DUT_rmsX[subjob] <<endl;
  77.  
  78. //------------------ For Y ------------------//
  79. cout << locationY << endl;
  80. Ycanvas[subjob] = new TCanvas(); //creating a canvas for each histogram
  81. Ycanvas[subjob]->SetTitle("Y Biased Residal for DUT"); //setting canvas title
  82. YDUT[subjob] = (TH1F*)openFile[subjob]->Get(locationY); //getting each of the histograms from the kepler root files
  83. YDUT[subjob]->Fit("gauss", "R"); //fitting guassian
  84. DUT_meanY[subjob]=YDUT[subjob]->GetFunction("gauss")->GetParameter(1); //getting the mean of the postion
  85. DUT_rmsY[subjob]=YDUT[subjob]->GetFunction("gauss")->GetParameter(2); //this should be the RMS
  86. cout << "DUT mean Y\t" << DUT_meanY[subjob] <<endl;
  87. cout << "DUT rms Y\t" << DUT_rmsY[subjob] <<endl;
  88.  
  89. }
  90.  
  91. //Printing to the DUT histograms and gaussian fits to a new root file
  92. char rootFile[100];
  93. sprintf(rootFile, "%s", "root_res_analysis.root");
  94. cout << "creating\t" << rootFile << endl;
  95. TFile *myfile = new TFile(rootFile,"RECREATE"); // creating rootfile
  96. for (int subjob =0; subjob<run; subjob++){
  97. XDUT[subjob]->Write(); //writing the plots to file
  98. YDUT[subjob]->Write(); //writing the plots to file
  99. }
  100. //Reading in the parameters from runList.txt
  101. char runList[100];
  102. char str1[10],str2[10],str3[10],str4[10],str5[10],str6[10];
  103. float a,b,c,f; char d[100], e[100];
  104. float full_Angle[100], full_Bias[100], length[100];
  105. float Angle[run], Bias[run], rmsX[run];
  106. int r =0;
  107. sprintf(runList, "%s", "/afs/cern.ch/user/e/ebuchana/cmtuser/KEPLER/KEPLER_HEAD/Tb/Kepler/options/ResStudies/runList.txt");
  108. FILE * runFile = fopen (runList,"read");
  109. cout << "reading in runList.txt " << endl;
  110. fscanf(runFile, "%s \t %s \t %s \t %s \t %s \t %s", str1,str2,str3,str4,str5,str6);
  111. while (!feof(runFile)){
  112. fscanf(runFile,"%f \t %f \t %f \t %s \t %s \t %f", &a, &b, &c, d, e, &f);
  113. full_Angle[r]=b;// all angles in runList.txt
  114. full_Bias[r]=c; // all bias voltages in runList.txt
  115. length[r]=f; // number of line in runList.txt
  116. r++;
  117. }
  118. //initalisng arrays (may not be needed)
  119. for (int range =beg-1; range <end; range ++){
  120. Angle[range]=0;
  121. Bias[range]=0;
  122. }
  123. //saving the range of bias and angles that you need from the full set of parameters
  124. int in =0;
  125. for (int range =beg-1; range <end; range ++){
  126. Angle[in]=full_Angle[range];
  127. Bias[in]=full_Bias[range];
  128. in++;
  129. }
  130.  
  131. //Plots after this line will be in the LHCb style format.
  132. gROOT-> ProcessLine(".x lhcbStyle.C");
  133. if (choice==0){ //if Bias scan
  134. //--------- for X ----------//
  135. TCanvas *c1 = new TCanvas("c1", "Bias vs Xrms (DUT)");
  136. TGraph *g1 = new TGraph(run, Bias, DUT_rmsX);
  137. g1->Draw("ap");
  138. g1->SetTitle("Bias vs. DUT rms; Bias Voltage; DUT X spatial resolution ");
  139. c1->Print(plot_bias);
  140. c1->Print(root_plotb);
  141. //--------- for Y ----------//
  142. TCanvas *c2 = new TCanvas("c2", "Bias vs Yrms (DUT)");
  143. TGraph *g2 = new TGraph(run, Bias, DUT_rmsY);
  144. g2->Draw("ap");
  145. g2->SetTitle("Bias vs. DUT rms; Bias Voltage; DUT Y spatial resolution ");
  146. // c2->Print(plot_bias);
  147. // c2->Print(root_plotb);
  148. }
  149. else{ //if angle scan
  150. TCanvas *c3 = new TCanvas("c3", "Angle vs Xrms (DUT)");
  151. TGraph *g3 = new TGraph(run, Angle, DUT_rmsX);
  152. g3->Draw("ap");
  153. c3->Print(plot_angle);
  154. c3->Print(root_plota);
  155. TCanvas *c4 = new TCanvas("c4", "Angle vs Yrms (DUT)");
  156. TGraph *g4 = new TGraph(run, Angle, DUT_rmsY);
  157. g4->Draw("ap");
  158. // c4->Print(plot_angle);
  159. // c4->Print(root_plota);
  160. }
  161. }
  162.  
  163.  
  164.