Newer
Older
TB_Chris / TbUT / scripts / CalibAna / .svn / text-base / analyzeCalib.cpp.svn-base
  1. //--------------------------------------------------------------------------------------------------------------------------------------------------------
  2. //
  3. // This code will append tracking information to a ROOT tree that contains cluster information from a UT DUT
  4. //
  5. // To compile
  6. // > SetupProject LHCb
  7. // > make
  8. //
  9. // Example of how to run interactively:
  10. //
  11. // ./combDUTwithTrack -i /data2/pmmannin/BoardA4_redo_v2/Run_Bias_Scan-B1-A-227-8711_Tuple.root -t /data2/sblusk/TB/July2015/TelescopeData/BoardA4/Kepler-tuple.root -o /data2/sblusk/test2.root
  12. //
  13. //
  14. // The Telescope data is on eos at:~/eos/lhcb/testbeam/ut/TelescopeData/July2015/RootFiles/RunXXXX/Output/Kepler-tuple.root
  15. // where XXXX = run number from TPIX DAQ system (Kepler run #)
  16. //
  17. // The DUT data needs to be processed through TbUT and put on eos in a standard area.
  18. //
  19. // The easiest, although ot safest way, to get the data directly accessible from eos is to use:
  20. // eosmount ~/eos
  21. //--------------------------------------------------------------------------------------------------------------------------------------------------------
  22.  
  23.  
  24. #include <iostream>
  25. #include <vector>
  26. #include <fstream>
  27. #include <math.h>
  28. #include <stdlib.h>
  29. #include "TFile.h"
  30. #include "TROOT.h"
  31. #include "TSystem.h"
  32. #include "TTree.h"
  33. #include "TH1F.h"
  34. #include "TH2F.h"
  35. #include "TF1.h"
  36. #include "TCanvas.h"
  37.  
  38. #ifdef __MAKECINT__
  39. #pragma link C++ class vector<float>+;
  40. #endif
  41.  
  42. using namespace std;
  43.  
  44.  
  45.  
  46. int main(int __argc, char *__argv[]){
  47.  
  48. gROOT->ProcessLine("#include <vector>");
  49. const char *calib_filename = "calib.root";
  50. const char *out_filename = "outputfile.root";
  51.  
  52. int c;
  53. extern char* optarg;
  54. extern int optind;
  55.  
  56. int channel;
  57. int type;
  58. cout << endl;
  59.  
  60. cout << "number of arguments: " << __argc << endl;
  61. for(int i=0; i < __argc; i++) {
  62.  
  63. cout << __argv[i] << endl;
  64. }
  65. /*____________________________Parse Command Line___________________________*/
  66. while((c = getopt(__argc,__argv,"i:o:n:t:p")) != -1){
  67. switch(c){
  68. case 'h': // help option
  69. cout << "Example: ./analyzeCalib -i calib.root -o outputfile.root -n channel -t type" << endl;
  70. return 0;
  71. break;
  72. case 'i':
  73. calib_filename = optarg;
  74. cout << "---> Calibration ROOT filename: " << calib_filename << endl;
  75. break;
  76. case 'o':
  77. out_filename = optarg;
  78. cout << "---> Output filename: " << out_filename << endl;
  79. break;
  80. case 'n':
  81. channel = atoi(optarg);
  82. cout << "Channel: " << channel << endl;
  83. break;
  84. case 't':
  85. type = atoi(optarg);
  86. cout << "Type: " << type << endl;
  87. break;
  88. default: // unknown option flag
  89. printf("Error!!!! Unknown option -%c\n",c);
  90. cout << "Example: analyzeCalib -i calib.root -o outputfile.root -n channel -t type" << endl;
  91. return 0;
  92. }
  93. }
  94.  
  95. cout << endl;
  96.  
  97. Int_t clusterNumberPerEvent;
  98. UInt_t clustersTDC;
  99. ULong64_t timestamps;
  100. Double_t clustersPosition[10];
  101. Int_t clustersSeedPosition[10];
  102. Double_t clustersCharge[10];
  103. Int_t clustersSize[10];
  104. Double_t clustersSeedCharge[10];
  105. Double_t clustersCharge2StripLeft[10];
  106. Double_t clustersCharge1StripLeft[10];
  107. Double_t clustersCharge1StripRight[10];
  108. Double_t clustersCharge2StripRight[10];
  109.  
  110. // List of branches
  111. TBranch *b_clusterNumberPerEvent; //!
  112. TBranch *b_clustersTDC; //!
  113. TBranch *b_timestamps; //!
  114. TBranch *b_clustersPosition; //!
  115. TBranch *b_clustersSeedPosition; //!
  116. TBranch *b_clustersCharge; //!
  117. TBranch *b_clustersSize; //!
  118. TBranch *b_clustersSeedCharge; //!
  119. TBranch *b_clustersCharge2StripLeft; //!
  120. TBranch *b_clustersCharge1StripLeft; //!
  121. TBranch *b_clustersCharge1StripRight; //!
  122. TBranch *b_clustersCharge2StripRight; //!
  123.  
  124. //dut_filename = "/data2/pmmannin/BoardA4_redo_v2/Run_Bias_Scan-B1-A-227-8711_Tuple.root";
  125. TFile *f_calib = new TFile(calib_filename,"READONLY");
  126. TTree* t_dut = (TTree*)f_calib->Get("TbUT/Clusters");
  127.  
  128. t_dut->SetBranchAddress("clusterNumberPerEvent", &clusterNumberPerEvent, &b_clusterNumberPerEvent);
  129. t_dut->SetBranchAddress("clustersTDC", &clustersTDC, &b_clustersTDC);
  130. t_dut->SetBranchAddress("timestamps", &timestamps, &b_timestamps);
  131. t_dut->SetBranchAddress("clustersPosition", clustersPosition, &b_clustersPosition);
  132. t_dut->SetBranchAddress("clustersSeedPosition", clustersSeedPosition, &b_clustersSeedPosition);
  133. t_dut->SetBranchAddress("clustersCharge", clustersCharge, &b_clustersCharge);
  134. t_dut->SetBranchAddress("clustersSize", clustersSize, &b_clustersSize);
  135. t_dut->SetBranchAddress("clustersSeedCharge", clustersSeedCharge, &b_clustersSeedCharge);
  136. t_dut->SetBranchAddress("clustersCharge2StripLeft", clustersCharge2StripLeft, &b_clustersCharge2StripLeft);
  137. t_dut->SetBranchAddress("clustersCharge1StripLeft", clustersCharge1StripLeft, &b_clustersCharge1StripLeft);
  138. t_dut->SetBranchAddress("clustersCharge1StripRight", clustersCharge1StripRight, &b_clustersCharge1StripRight);
  139. t_dut->SetBranchAddress("clustersCharge2StripRight", clustersCharge2StripRight, &b_clustersCharge2StripRight);
  140.  
  141. int numEvents = t_dut->GetEntriesFast();
  142. cout << " ------------------------------------------------------" << endl;
  143. cout << " | Number of events found = " << numEvents << endl;
  144. cout << " ------------------------------------------------------" << endl;
  145.  
  146. vector<int> maxStrips;
  147.  
  148. int nfound=0;
  149. //out_filename = "/data2/sblusk/test.root";
  150. TFile *f_out = new TFile(out_filename,"recreate");
  151.  
  152. //TFile *fileout = new TFile("temp.root","RECREATE");
  153. TH1F *h_chg = new TH1F("h_chg","Cluster Charge",200,-2000,2000);
  154. TH1F *h_chgLeft = new TH1F("h_chgLeft","Charge Charge (ch-1)",200,-2000,2000);
  155. TH1F *h_chgLeft2 = new TH1F("h_chgLeft2","Charge Charge (ch-2)",200,-2000,2000);
  156. TH1F *h_chgRight = new TH1F("h_chgRight","Charge Charge (ch+1)",200,-2000,2000);
  157. TH1F *h_chgRight2 = new TH1F("h_chgRight2","Charge Charge (ch+2)",200,-2000,2000);
  158. TH1F *h_chgDiff = new TH1F("h_chgDiff","Charge_{left}-Charge_{right}",200,-2000,2000);
  159. TH1F *h_chgDiff2 = new TH1F("h_chgDiff2","Charge_{s-2}-Charge_{s+2}",200,-2000,2000);
  160. TH1F *h_seed = new TH1F("h_seed","Seed Strip",513,-0.5,512.5);
  161. TH1F *h_pos = new TH1F("h_pos","Strip",513,-0.5,512.5);
  162.  
  163. int strip;
  164. int left;
  165. int right;
  166.  
  167. for(int i=0;i < numEvents; i++){
  168. int nb = t_dut->GetEntry(i);
  169. if(nb <= 0) break;
  170. for(int j=0; j < min(clusterNumberPerEvent,10);j++) {
  171. h_pos->Fill(clustersPosition[j]);
  172. }
  173. }
  174.  
  175.  
  176. /*
  177. for(int i=1; i < 511; i++) {
  178. if(h_pos->GetBinContent(i) > 2500) {
  179. maxStrips.push_back((int)h_pos->GetBinCenter(i));
  180. cout << (int)h_pos->GetBinCenter(i) << endl;
  181. }
  182. }
  183.  
  184. for(int i=0; i < maxStrips.size(); i++) {
  185. cout << "Max Strip " << i << ": " << maxStrips.at(i) << endl;
  186. }
  187. */
  188.  
  189.  
  190. if(type == 0) strip = channel+128;
  191. else if(type == 1) strip = channel;
  192. else { cout << "Please enter a valid type: 0 for A, 1 for D, 2 for Micron Mini, 3 for Ham Mini" << endl; return 0; }
  193.  
  194. left = strip-1;
  195. right = strip+1;
  196.  
  197. for(int i=0;i < numEvents; i++){
  198. int nb = t_dut->GetEntry(i);
  199. if(nb <= 0) break;
  200. double chgLeft = 0.0;
  201. double chgRight = 0.0;
  202. for(int j=0; j < min(clusterNumberPerEvent,10);j++) {
  203. if(fabs( clustersPosition[j]-strip) < 0.5 ) {
  204. h_chg->Fill(clustersCharge[j]);
  205. h_seed->Fill(clustersSeedPosition[j]);
  206. h_chgLeft->Fill(clustersCharge1StripLeft[j]);
  207. h_chgRight->Fill(clustersCharge1StripRight[j]);
  208. h_chgLeft2->Fill(clustersCharge2StripLeft[j]);
  209. h_chgRight2->Fill(clustersCharge2StripRight[j]);
  210. h_chgDiff->Fill(clustersCharge1StripLeft[j]-clustersCharge1StripRight[j]);
  211. h_chgDiff2->Fill(clustersCharge2StripLeft[j]-clustersCharge2StripRight[j]);
  212. }
  213. }
  214. }
  215.  
  216. TCanvas *c1 = new TCanvas("c","c",1200,800);
  217.  
  218. c1->Divide(3,3);
  219. c1->cd(1)->SetLeftMargin(0.13);
  220. h_chg->SetTitle(TString(Form("Seed Channel %d (strip %d)",channel,strip)));
  221. h_chg->GetXaxis()->SetTitle("Cluster Charge");
  222. h_chg->Draw();
  223.  
  224. c1->cd(2)->SetLeftMargin(0.13);
  225. h_chgLeft->SetTitle(TString(Form("Left Channel: %d (strip %d)",channel-1,left)));
  226. h_chgLeft->GetXaxis()->SetTitle("Cluster Charge, N-1");
  227. h_chgLeft->Draw();
  228. c1->cd(3)->SetLeftMargin(0.13);
  229. h_chgRight->SetTitle(TString(Form("Right Channel: %d (strip %d)",channel+1,right)));
  230. h_chgRight->GetXaxis()->SetTitle("Cluster Charge, N+1");
  231. h_chgRight->Draw();
  232.  
  233. c1->cd(4)->SetLeftMargin(0.13);
  234. h_chgLeft2->SetTitle(TString(Form("2^{nd} Left Channel: %d (strip %d)",channel-2,left-1)));
  235. h_chgLeft2->GetXaxis()->SetTitle("Cluster Charge, N-2");
  236. h_chgLeft2->Draw();
  237. c1->cd(5)->SetLeftMargin(0.13);
  238. h_chgRight2->SetTitle(TString(Form("2^{nd} Right Channel: %d (strip %d)",channel+2,right+1)));
  239. h_chgRight2->GetXaxis()->SetTitle("Cluster Charge, N+2");
  240. h_chgRight2->Draw();
  241.  
  242.  
  243. c1->cd(6)->SetLeftMargin(0.13);
  244. h_chgDiff->GetXaxis()->SetTitle("Cluster Charge Difference (s-1,s+1)");
  245. h_chgDiff->Draw();
  246. c1->cd(7)->SetLeftMargin(0.13);
  247. h_chgDiff2->GetXaxis()->SetTitle("Cluster Charge Difference (s-2,s+2)");
  248. h_chgDiff2->Draw();
  249. c1->cd(8)->SetLeftMargin(0.13);
  250. h_seed->GetXaxis()->SetTitle("Seed Cluster Position");
  251. h_seed->GetXaxis()->SetRangeUser(h_seed->GetMaximumBin()-10,h_seed->GetMaximumBin()+10);
  252. h_seed->Draw();
  253.  
  254.  
  255.  
  256. c1->SaveAs(TString(Form("Calibration_Type%d_Channel%d.png",type,channel)));
  257.  
  258. // Write and close ROOT file.
  259. f_out->Write();
  260. f_out->Close();
  261.  
  262. return 0;
  263. }