#include <map>#include "TKey.h"#include <TROOT.h>#include <TRint.h>#include <TTree.h>#include <TFile.h>#include <TSystem.h>#include <iostream>#include <TF1.h>#include <TH2.h>#include <TGraph.h>#include <sstream>#include "Log.hh"#include <fstream>#include "TCanvas.h"#include "TApplication.h"#include "TColor.h"#include "TStyle.h"#include "TAxis.h"#include "TString.h"#include "TLegend.h"#include "Bytes.h"#include "root/MTRun.hh"#include "root/MTEvent.hh"#include "root/MTChannel.hh"#include "root/MTMicrorocChip.hh"#include "MicroException.hh"Include dependency graph for test_microroc.cpp:
Go to the source code of this file.
Defines | |
| #define | WRITE |
Functions | |
| void | usage (char *progname) |
| int | main (int argc, char **argv) |
| #define WRITE |
Definition at line 41 of file test_microroc.cpp.
| void usage | ( | char * | progname | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 51 of file test_microroc.cpp.
00051 { 00052 00053 /*****************************************************************/ 00054 /* ARGUMENTS TO THE PROGRAM **************************************/ 00055 /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */ 00056 /* 2 - NAME OF THE ROOT FILE *************************************/ 00057 /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/ 00058 /* 0 -> DAC_0 (LOW THRESHOLD) *********************************/ 00059 /* 1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/ 00060 /* 2 -> DAC_2 (HIGH THRESHOLD) ********************************/ 00061 /*****************************************************************/ 00062 if(argc != 4) 00063 { usage(argv[0]) ; 00064 return 1 ; 00065 } 00066 TString root_directory = argv[1]; 00067 TString root_filename = argv[2]; 00068 int dac = atoi(argv[3]); 00069 00070 //bool display = 0; 00071 #ifdef PLOT 00072 set_style_myplain() ; 00073 #endif 00074 TCanvas *c1 = new TCanvas("c1","Test MICROROC",0,0,1000,1000); 00075 TString root_file = root_directory + root_filename; 00076 TString scurve_file = root_directory + "scurve_" + root_filename; 00077 00078 // cout<<"Wish to write "<<scurve_file<<endl; 00079 00080 int nbChannel = 0; 00081 float nbEvt = 0; 00082 int threshold = 0; 00083 //int chipid = 0; 00084 int hardid = 0; 00085 int value = 0; 00086 //int content = 0; 00087 Double_t content = 0; 00088 UInt_t chipId ; 00089 00090 TH1I * hvalue = new TH1I("hvalue","",5,0,5); 00091 TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024); 00092 00093 TH1I * hscurve[64]; 00094 TString name,title; 00095 00096 for (int j=0;j<64;j++) 00097 { name = Form("hscurve_%i",j); 00098 title = Form("SCurve from channel %i",j); 00099 hscurve[j] = new TH1I(name,title,1024,0,1024); 00100 } 00101 00102 TFile f(root_file); 00103 TIter nextkey(f.GetListOfKeys()); 00104 TKey *key; 00105 while (key = (TKey*)nextkey()) 00106 { TTree *tree = (TTree*)key->ReadObj(); 00107 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00108 const MTChip& chip = run->GetOneChip(); 00109 // cout << "------------ Chip[" << chip.GetId() << "], softId: [" << chip.GetSoftId() <<"] configuration ---------------" << endl ; 00110 00111 try 00112 { const MTMicrorocChip µChip = dynamic_cast<const MTMicrorocChip &> (chip); 00113 //const MTHardroc2Chip µChip = dynamic_cast<const MTHardroc2Chip &> (chip); 00114 //cout << "threshold dac 0:[" << microChip.GetThresholdDac_0() << endl; 00115 //cout << "threshold dac 1:[" << microChip.GetThresholdDac_1() << endl; 00116 //cout << "threshold dac 2:[" << microChip.GetThresholdDac_2() << endl; 00117 //chipId = chip.GetId() ; 00118 //cout << "Chip identifier: " << chipId << endl ; 00119 if (dac==0){threshold = microChip.GetThresholdDac_0();} 00120 if (dac==1){threshold = microChip.GetThresholdDac_1();} 00121 if (dac==2){threshold = microChip.GetThresholdDac_2();} 00122 } 00123 catch (...) {} 00124 00125 MTEvent *evt = new MTEvent(); 00126 TBranch *branch= tree->GetBranch("MTEvent"); 00127 branch->SetAddress(&evt); 00128 MTChannel* channel = NULL; 00129 nbEvt = tree->GetEntries(); 00130 // cout << "Number of events in the tree: " << nbEvt << endl ; 00131 00132 hnpulse_threshold->SetBinContent(threshold+1,nbEvt); 00133 00134 for( int evtNum = 0; evtNum < nbEvt ; evtNum++) 00135 { //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;} 00136 tree->GetEntry(evtNum); 00137 nbChannel = evt->GetNchannel(); 00138 for(int i=0;i<nbChannel ;i++) 00139 { channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i); 00140 chipId = channel->GetChipId(); 00141 hardid = channel->GetHardId(); 00142 value = channel->GetDigitalValue(); 00143 hvalue->Fill(value); 00144 if (value==(dac+1)) 00145 { hscurve[hardid]->Fill(threshold); 00146 #ifdef R_DEBUG 00147 cout << "---- digital value:" << value <<endl; 00148 cout << "---- hard id:" << hardid <<endl; 00149 //cout << "---- x:"<< channel->GetX() <<endl; 00150 //cout << "---- y:" << channel->GetY() <<endl; 00151 cout << "---- chip id:" << chipId <<endl; 00152 #endif 00153 00154 } 00155 } 00156 } 00157 } 00158 00159 00160 00161 00162 //===========================================================================* 00163 // Now analyze scurves 00164 00165 TString t_scurve_file = root_directory + root_filename ;//+ ".root"; 00166 //cout << t_scurve_file<< endl; 00167 00168 int tmin = 0; 00169 int tmed = 0; 00170 int tmax = 0; 00171 int delta = 0; 00172 00173 // TH1I * hscurve[64]; 00174 00175 TH1F * hinflex = new TH1F("hinflex","",1024,0,1024); 00176 TH1F * hsigma = new TH1F("hsigma","",50,-2,2); 00177 00178 TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10); 00179 TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2); 00180 00181 TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64); 00182 TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64); 00183 00184 TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64); 00185 TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64); 00186 00187 TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64); 00188 00189 // TString name; 00190 00191 /* for (int j=0;j<64;j++) 00192 { name = Form("hscurve_%i",j); 00193 TH1I * h = (TH1I*)tf.Get(name); 00194 hscurve[j] = h; 00195 } 00196 */ 00197 TString func = "[0]/(1+exp((x-[1])/[2]))"; 00198 TF1 * fermi[64]; 00199 float scurve_max = 0; 00200 float scurve_inflex = 0; 00201 float inflex = 0; 00202 float sigma = 0; 00203 float inflex_err = 0; 00204 float sigma_err = 0; 00205 float content0,content1,content2; 00206 float chi2 = 0; 00207 float ndf = 0; 00208 00209 for (int j=0;j<64;j++) 00210 { name = Form("fermi_%i",j); 00211 fermi[j] = new TF1(name,func,0,1024); 00212 //fermi[j]->SetLineColor(j+1); 00213 fermi[j]->SetLineColor(kRed); 00214 fermi[j]->SetLineWidth(2); 00215 /*If SCurve exist*/ 00216 if (hscurve[j]->GetEntries()!=0) 00217 { /*Find threshold min tmin*/ 00218 for (int b=0;b<1023;b++) 00219 { content = hscurve[j]->GetBinContent(b+1); 00220 if (content != 0) 00221 { tmin = b ; 00222 b = 1023 ; 00223 scurve_max = content; 00224 } 00225 } 00226 /*Find threshold at which scurve = 0 tmax*/ 00227 for (int b=1024;b>tmin;b--) 00228 { if (hscurve[j]->GetBinContent(b+1) != 0) 00229 { tmax = b+1 ; 00230 b = tmin ; 00231 } 00232 } 00233 /*Find threshold at which scurve reaches the max*/ 00234 for (int b=tmax;b>=0;b--) 00235 { if (hscurve[j]->GetBinContent(b+1) == scurve_max) 00236 { tmed = b+1 ; 00237 b = -1 ; 00238 } 00239 } 00240 scurve_inflex = 0.5*(tmax+tmed); 00241 delta = tmax - tmed; 00242 00243 /*If SCurve is a step function, do not fit parameters*/ 00244 if (delta==0) 00245 { inflex = scurve_inflex; 00246 sigma = 0.; 00247 inflex_err = 0.; 00248 sigma_err = 0.; 00249 } 00250 else 00251 { /*Set parameters*/ 00252 fermi[j]->FixParameter(0,scurve_max); 00253 fermi[j]->SetParameter(1,scurve_inflex); 00254 fermi[j]->SetParameter(2,delta/10.); 00255 /*Perform fit*/ 00256 hscurve[j]->Fit(fermi[j],"q","",0,tmax+50); 00257 /*Get parameters*/ 00258 inflex = fermi[j]->GetParameter(1); 00259 sigma = sqrt(pow(fermi[j]->GetParameter(2),2)); 00260 inflex_err = sqrt(pow(fermi[j]->GetParError(1),2)); 00261 sigma_err = sqrt(pow(fermi[j]->GetParError(2),2)); 00262 chi2 = fermi[j]->GetChisquare(); 00263 ndf = fermi[j]->GetNDF(); 00264 } 00265 00266 hc_inflex_sigma->SetBinContent(j+1,inflex); 00267 hc_inflex_sigma->SetBinError(j+1,sigma); 00268 00269 hc_inflex->SetBinContent(j+1,inflex); 00270 hc_sigma->SetBinContent(j+1,sigma); 00271 00272 hc_inflex_err->SetBinContent(j+1,inflex_err); 00273 hc_sigma_err->SetBinContent(j+1,sigma_err); 00274 00275 hinflex->Fill(inflex); 00276 hsigma->Fill(sigma); 00277 00278 #ifdef DISPLAY 00279 cout<<" Channel\t" << j << "\t" << inflex << "\t width\t" << sigma << endl ; 00280 #endif 00281 } 00282 } 00283 00284 //////////////////////////////////////////////////////////////////////////////// 00285 // Write Scurves 00286 TFile * tf = new TFile(scurve_file,"RECREATE"); 00287 hvalue->Write(); 00288 hnpulse_threshold->Write(); 00289 for (int j=0;j<64;j++) 00290 { hscurve[j]->Write();} 00291 tf->Close(); 00292 //////////////////////////////////////////////////////////////////////////////// 00293 00294 #ifdef WRITE 00295 root_filename.ReplaceAll(".root", ""); 00296 TString outfile = root_directory + "../RESULTS/" + "parameters_" + root_filename + ".txt"; 00297 ofstream out ; 00298 out.open(outfile) ; 00299 for (int j=0;j<64;j++) 00300 { out<<j<<" "<<hc_inflex->GetBinContent(j+1)<<" "<<hc_inflex_err->GetBinContent(j+1)<<" "<<hc_sigma->GetBinContent(j+1)<<" "<<hc_sigma_err->GetBinContent(j+1)<<endl; 00301 } 00302 #endif 00303 #ifdef PLOT 00304 TApplication theApp("App", &argc, argv); 00305 TString adj; 00306 if (dac==0){adj = "Low";} 00307 if (dac==1){adj = "Medium";} 00308 if (dac==2){adj = "High";} 00309 00310 c1->Divide(8,8); 00311 //c1->cd(); 00312 TString htitle = adj + " threshold"; 00313 for (int j=0;j<64;j++) 00314 { //c0->cd(j/8,j%8); 00315 c1->cd(j+1); 00316 hscurve[j]->SetLineWidth(2); 00317 hscurve[j]->Draw(""); 00318 c1->Update() ; 00319 } 00320 theApp.Run(kTRUE) ; 00321 c1->Update() ; 00322 #endif 00323 return 0 ; 00324 }
1.4.7