Newer
Older
TB_Chris / TbUT / scripts / AddTrigTracks / .svn / text-base / combDUTwithTrack.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.  
  28. #include "TFile.h"
  29. #include "TROOT.h"
  30. #include "TSystem.h"
  31. #include "TTree.h"
  32. #include "TH1F.h"
  33. #include "TH2F.h"
  34. #include "TF1.h"
  35.  
  36. #ifdef __MAKECINT__
  37. #pragma link C++ class vector<float>+;
  38. #endif
  39.  
  40. using namespace std;
  41.  
  42. int main(int __argc, char *__argv[]){
  43.  
  44. gROOT->ProcessLine("#include <vector>");
  45. const char *dut_filename = "dut.root";
  46. const char *tp3_filename = "tp3.root";
  47. const char *out_filename = "outputfile.root";
  48. const char *noise_filename = "noise.root";
  49.  
  50. int c;
  51. extern char* optarg;
  52. extern int optind;
  53.  
  54. int _triggerPlane(-1);
  55. int _timeOffset = 0;
  56.  
  57. cout << endl;
  58. /*____________________________Parse Command Line___________________________*/
  59. while((c = getopt(__argc,__argv,"i:t:n:o:h")) != -1){
  60. switch(c){
  61. case 'h': // help option
  62. cout << "Example: ./combineDUTwithTrack -i dut.root -t tp3.root -o outputfile.root" << endl;
  63. return 0;
  64. break;
  65. case 'i':
  66. dut_filename = optarg;
  67. cout << "---> DUT ROOT filename: " << dut_filename << endl;
  68. break;
  69. case 't':
  70. tp3_filename = optarg;
  71. cout << "---> TimePix3 ROOT filename: " << tp3_filename << endl;
  72. break;
  73. case 'o':
  74. out_filename = optarg;
  75. cout << "---> Combined output filename: " << out_filename << endl;
  76. break;
  77. case 'n':
  78. noise_filename = optarg;
  79. cout << "---> Noise ROOT filename: " << noise_filename << endl;
  80. break;
  81. default: // unknown option flag
  82. printf("Error!!!! Unknown option -%c\n",c);
  83. cout << "Example: combineDUTwithTrack -i dut.root -t tp3.root -o outputfile.root" << endl;
  84. return 0;
  85. }
  86. }
  87.  
  88. cout << endl;
  89.  
  90. // get info from Kepler file
  91. double trk_x, trk_y, trk_tx, trk_ty, chi2ndf, trk_htime, trg_htime0;
  92. ULong64_t trg_time0, trg_time0_25ns;
  93. ULong64_t trk_time, trk_time_25ns, timestamps;
  94. ULong64_t time_window = 1000;
  95. ULong64_t diff_time;
  96. ULong64_t diff_time_loc;
  97. UInt_t trg_Plane;
  98. TBranch *b_timestamps; //!
  99. cout << "Made it this far!" << endl;
  100. //dut_filename = "/data2/pmmannin/BoardA4_redo_v2/Run_Bias_Scan-B1-A-227-8711_Tuple.root";
  101. TFile *f_ut = new TFile(dut_filename,"READONLY");
  102. TFile *f_noise = new TFile(noise_filename,"READONLY");
  103. TTree* t_dut = (TTree*)f_ut->Get("TbUT/Clusters");
  104. t_dut->SetBranchAddress("timestamps",&timestamps, &b_timestamps);
  105.  
  106. cout << "Made it a little further! " << endl;
  107. int numEvents = t_dut->GetEntriesFast();
  108. cout << " ------------------------------------------------------" << endl;
  109. cout << " | Number of triggers found = " << numEvents << endl;
  110. cout << " ------------------------------------------------------" << endl;
  111.  
  112. //tp3_filename = "/data2/sblusk/TB/July2015/TelescopeData/BoardA4/Kepler-tuple.root";
  113. TFile *f_tpix = new TFile(tp3_filename,"READONLY");
  114. TTree *t_trk = (TTree*)f_tpix->Get("TbTupleWriter/Tracks");
  115. t_trk->SetBranchAddress("TkTime",&trk_time);
  116. t_trk->SetBranchAddress("TkHTime",&trk_htime);
  117. t_trk->SetBranchAddress("TkX",&trk_x);
  118. t_trk->SetBranchAddress("TkY",&trk_y);
  119. t_trk->SetBranchAddress("TkTx",&trk_tx);
  120. t_trk->SetBranchAddress("TkTy",&trk_ty);
  121. t_trk->SetBranchAddress("TkChi2PerNdof",&chi2ndf);
  122.  
  123. int nfound=0;
  124. //out_filename = "/data2/sblusk/test.root";
  125. TFile *f_out = new TFile(out_filename,"recreate");
  126. TTree *t_out = (TTree*)t_dut->CloneTree(0); // grab a copy of the existing DUT tree
  127.  
  128. //TFile *fileout = new TFile("temp.root","RECREATE");
  129. TH1F *h1 = new TH1F("h1","Time Diff",1000,-1000,1000);
  130. TH1F *hMeanNoise = new TH1F("hMeanNoise","Mean of noise, CM subtracted",512,0,512);
  131. TH1F *hWidthNoise = new TH1F("hWidthNoise","Gaussian noise, CM subtracted",512,0,512);
  132. TH1D *px[512];
  133.  
  134. // add the new branches w/ track info
  135. int n_tp3_tracks;
  136. std::vector<double> *vec_trk_x = 0;
  137. std::vector<double> *vec_trk_y = 0;
  138. std::vector<double> *vec_trk_tx = 0;
  139. std::vector<double> *vec_trk_ty = 0;
  140. std::vector<double> *vec_trk_htime = 0;
  141. std::vector<double> *vec_trk_chi2ndf = 0;
  142.  
  143.  
  144. t_out->Branch("n_tp3_tracks",&n_tp3_tracks);
  145. t_out->Branch("vec_trk_x",&vec_trk_x);
  146. t_out->Branch("vec_trk_y",&vec_trk_y);
  147. t_out->Branch("vec_trk_tx",&vec_trk_tx);
  148. t_out->Branch("vec_trk_ty",&vec_trk_ty);
  149. t_out->Branch("vec_trk_chi2ndf",&vec_trk_chi2ndf);
  150. //t_out->Branch("vec_trk_htime",&vec_trk_htime); // don't need this one
  151.  
  152. double dtime;
  153. t_out->Branch("dtime",&dtime);
  154.  
  155.  
  156. int entries = (int)t_trk->GetEntries();
  157. int curr_pos(0), j, stop_now;
  158.  
  159. // First short pass, check trigger timing
  160. int iAli = 0;
  161. for(int i=0;i < min(numEvents,10000); i++){
  162. int nb = t_dut->GetEntry(i);
  163. if(nb <= 0) break;
  164. //t_dut->GetEntry(i);
  165. if(i%1000==0) cout << "\nprocessing " << i << " curr_pos: " << curr_pos << " " << endl;
  166. trg_time0_25ns = timestamps - _timeOffset;
  167.  
  168. stop_now = 0;
  169. j = curr_pos;
  170. while(j < entries && stop_now == 0){
  171. t_trk->GetEntry(j);
  172. trk_time_25ns = (trk_time>>12);
  173. if(trg_time0_25ns > trk_time_25ns){
  174. diff_time_loc = trg_time0_25ns - trk_time_25ns;
  175. }
  176. if(trk_time_25ns > trg_time0_25ns){
  177. diff_time_loc = trk_time_25ns - trg_time0_25ns;
  178. if(diff_time_loc > time_window) {stop_now=1; curr_pos = j-1;}
  179. }
  180. if(trk_time_25ns == trg_time0_25ns){
  181. diff_time_loc = trk_time_25ns - trg_time0_25ns;
  182. }
  183.  
  184. double dt = diff_time_loc * 1.0;
  185. h1->Fill(dt);
  186. j++;
  187. }
  188.  
  189. if(curr_pos < 0) curr_pos = 0;
  190. }
  191.  
  192. double _timeOffsetNew = _timeOffset;
  193. int ib = h1->GetMaximumBin();
  194. if(ib>0 && ib<h1->GetNbinsX()) _timeOffsetNew = _timeOffsetNew + h1->GetBinCenter(ib);
  195. cout << "-----------------------------------------------" << endl;
  196. cout << "Resetting time offset to " << _timeOffsetNew << endl;
  197. cout << "Setting time window to 5 ns " << endl;
  198. cout << "-----------------------------------------------" << endl;
  199.  
  200. time_window = 5;
  201. h1->Reset();
  202.  
  203. nfound=0;
  204. iAli = 0;
  205. curr_pos = 0;
  206. //numEvents = 5000;
  207. for(int i=0;i < numEvents; i++){
  208. int nb = t_dut->GetEntry(i);
  209. if(nb <= 0) break;
  210. if(i%10000==0) cout << "\nprocessing " << i << " curr_pos: " << curr_pos << " " << endl;
  211.  
  212. trg_time0_25ns = timestamps - _timeOffsetNew; // in units of 25ns minus the offset
  213. vec_trk_x->clear();
  214. vec_trk_y->clear();
  215. vec_trk_tx->clear();
  216. vec_trk_ty->clear();
  217. vec_trk_chi2ndf->clear();
  218. //vec_trk_htime->clear();
  219.  
  220. stop_now = 0;
  221. j = curr_pos;
  222. while(j < entries && stop_now == 0){
  223. int nb2 = t_trk->GetEntry(j);
  224. if(nb2 <= 0) break;
  225. trk_time_25ns = (trk_time>>12);
  226. if(trg_time0_25ns > trk_time_25ns){
  227. diff_time_loc = trg_time0_25ns - trk_time_25ns;
  228. }
  229. if(trk_time_25ns > trg_time0_25ns){
  230. diff_time_loc = trk_time_25ns - trg_time0_25ns;
  231. if(diff_time_loc > time_window) {stop_now=1; curr_pos = j-1;}
  232. }
  233. if(trk_time_25ns == trg_time0_25ns){
  234. diff_time_loc = trk_time_25ns - trg_time0_25ns;
  235. }
  236.  
  237. double dt = diff_time_loc * 1.0;
  238. h1->Fill(dt);
  239.  
  240. //if(diff_time_loc - _timeOffsetNew < time_window){
  241. if(dt < time_window){
  242. //vec_trk_htime->push_back((trk_htime - trg_htime0));
  243. //cout << diff_time_loc - _timeOffsetNew << " " << dt << endl;
  244. dtime = dt;
  245. vec_trk_x->push_back(trk_x);
  246. vec_trk_y->push_back(trk_y);
  247. vec_trk_tx->push_back(trk_tx);
  248. vec_trk_ty->push_back(trk_ty);
  249. vec_trk_chi2ndf->push_back(chi2ndf);
  250. }
  251. j++;
  252. }
  253.  
  254. if(curr_pos < 0) curr_pos = 0;
  255. n_tp3_tracks = (int)vec_trk_x->size();
  256.  
  257. if(n_tp3_tracks > 0) nfound++;
  258. if(nfound%1000==0 && nfound > 0) cout << "Found " << nfound << endl;
  259. t_out->Fill();
  260. }
  261.  
  262. //fileout->Write();
  263. //fileout->Close();
  264.  
  265. //-------------- Do noise Plots --------------//
  266. //--------------------------------------------//
  267.  
  268.  
  269. f_noise->cd("TbUT");
  270. TH2D *hNoise = (TH2D*)gDirectory->Get("CMSData_vs_channel");
  271.  
  272. f_out->cd();
  273.  
  274. TF1 *gau = new TF1("gau","gaus(0)",-120,120);
  275. gau->SetParameters(10000,0.0,40.0);
  276.  
  277. int i = 128;
  278. for(int i=0;i<512;i++){
  279. if( i>=512 ) {
  280. hMeanNoise->SetBinContent(i,-1000);
  281. continue;
  282. }
  283. px[i] = hNoise->ProjectionY(Form("px%d",i),i,i);
  284. px[i]->GetXaxis()->SetRangeUser(-200,200);
  285. px[i]->SetTitle(Form("Channel %d", i));
  286. px[i]->SetName(Form("Ch%d", i));
  287. if(px[i]->GetEntries() < 10){
  288. hMeanNoise->SetBinContent(i,-1000);
  289. continue;
  290. }
  291. double p = px[i]->GetMaximum();
  292. double m = px[i]->GetMean();
  293. double e = px[i]->GetRMS();
  294. gau->SetParameters(p, m, e);
  295. gau->SetRange(m-3*e,m+3*e);
  296.  
  297. px[i]->Fit(gau,"RQ");
  298. double mn = gau->GetParameter(1);
  299. double wid = gau->GetParameter(2);
  300. double emn = gau->GetParError(1);
  301. double ewid = gau->GetParError(2);
  302. hMeanNoise->SetBinContent(i,mn);
  303. hMeanNoise->SetBinError(i,emn);
  304. hWidthNoise->SetBinContent(i,wid);
  305. hWidthNoise->SetBinError(i,ewid);
  306. }
  307.  
  308. // Save histograms
  309. hMeanNoise->Write();
  310. hWidthNoise->Write();
  311. for(int i=0;i<512;i++){
  312. px[i]->Write();
  313. }
  314. t_out->AutoSave();
  315. delete f_ut;
  316. delete f_noise;
  317. delete f_tpix;
  318. delete f_out;
  319. // Write and close ROOT file.
  320. //f_out->Write();
  321. //f_out->Close();
  322.  
  323. return 0;
  324. }