#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_2.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 40 of file test_microroc_2.cpp.
| void usage | ( | char * | progname | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 54 of file test_microroc_2.cpp.
00054 { 00055 00056 /*****************************************************************/ 00057 /* ARGUMENTS TO THE PROGRAM **************************************/ 00058 /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */ 00059 /* 2 - NAME OF THE ROOT FILE *************************************/ 00060 /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/ 00061 /* 0 -> DAC_0 (LOW THRESHOLD) *********************************/ 00062 /* 1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/ 00063 /* 2 -> DAC_2 (HIGH THRESHOLD) ********************************/ 00064 /*****************************************************************/ 00065 if(argc != 4) 00066 { usage(argv[0]) ; 00067 return 1 ; 00068 } 00069 TString root_directory = argv[1]; 00070 TString root_filename = argv[2]; 00071 int dac = atoi(argv[3]); 00072 00073 //bool display = 0; 00074 #ifdef PLOT 00075 set_style_myplain() ; 00076 #endif 00077 TCanvas *c1 = new TCanvas("c1","Test MICROROC",0,0,1000,1000); 00078 TString root_file = root_directory + root_filename; 00079 TString scurve_file = root_directory + "scurve_" + root_filename; 00080 00081 // cout<<"Wish to write "<<scurve_file<<endl; 00082 00083 00084 // cout << "Attention ! This version is modified for temperature testing ! Else, comment line 123 of source !" << endl ; 00085 00086 int nbChannel = 0; 00087 float nbEvt = 0; 00088 int threshold = 0; 00089 //int chipid = 0; 00090 int hardid = 0; 00091 int value = 0; 00092 int serial = 0; 00093 //int content = 0; 00094 Double_t content = 0; 00095 UInt_t chipId ; 00096 int dac_min = 1025; 00097 int dac_max = 0; 00098 00099 TH1I * hvalue = new TH1I("hvalue","",5,0,5); 00100 TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024); 00101 00102 TH1I * hscurve[64] ; 00103 TH1I * hscurve_deriv[64] ; 00104 TString name, title, name_deriv, title_deriv ; 00105 00106 for (int j=0;j<64;j++) 00107 { name = Form("hscurve_%i",j); 00108 title = Form("SCurve from channel %i",j); 00109 hscurve[j] = new TH1I(name,title,1024,0,1024); 00110 name_deriv = Form("deriv_%i",j); 00111 title_deriv = Form("Derivative from channel %i",j); 00112 hscurve_deriv[j] = new TH1I(name_deriv,title_deriv,1024,0,1024); 00113 } 00114 00115 TFile f(root_file); 00116 TIter nextkey(f.GetListOfKeys()); 00117 TKey *key; 00118 while (key = (TKey*)nextkey()) 00119 { TTree *tree = (TTree*)key->ReadObj(); 00120 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00121 const MTChip& chip = run->GetOneChip(); 00122 // cout << "------------ Chip[" << chip.GetId() << "], softId: [" << chip.GetSoftId() <<"] configuration ---------------" << endl ; 00123 00124 serial = chip.GetSerialNumber() ; 00125 // cout << serial << " " ; 00126 // if (serial == 14) { 00127 try 00128 { const MTMicrorocChip µChip = dynamic_cast<const MTMicrorocChip &> (chip); 00129 //chipId = chip.GetId() ; 00130 //cout << "Chip identifier: " << chipId << endl ; 00131 if (dac==0){threshold = microChip.GetThresholdDac_0();} 00132 if (dac==1){threshold = microChip.GetThresholdDac_1();} 00133 if (dac==2){threshold = microChip.GetThresholdDac_2();} 00134 } 00135 catch (...) {} 00136 /*Find DAC range*/ 00137 if (threshold<dac_min){dac_min = threshold;} 00138 if (threshold>dac_max){dac_max = threshold;} 00139 00140 MTEvent *evt = new MTEvent(); 00141 TBranch *branch= tree->GetBranch("MTEvent"); 00142 branch->SetAddress(&evt); 00143 MTChannel* channel = NULL; 00144 nbEvt = tree->GetEntries(); 00145 // cout << "Number of events in the tree: " << nbEvt << endl ; 00146 00147 hnpulse_threshold->SetBinContent(threshold+1,nbEvt); 00148 00149 for( int evtNum = 0; evtNum < nbEvt ; evtNum++) 00150 { //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;} 00151 tree->GetEntry(evtNum); 00152 nbChannel = evt->GetNchannel(); 00153 for(int i=0;i<nbChannel ;i++) 00154 { channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i); 00155 chipId = channel->GetChipId(); 00156 hardid = channel->GetHardId(); 00157 value = channel->GetDigitalValue(); 00158 hvalue->Fill(value); 00159 if (value==(dac+1)) // ie 1 for low, 2 for med and 3 for high .... 00160 { hscurve[hardid]->Fill(threshold); 00161 #ifdef R_DEBUG 00162 cout << "---- digital value:" << value <<endl; 00163 cout << "---- hard id:" << hardid <<endl; 00164 //cout << "---- x:"<< channel->GetX() <<endl; 00165 //cout << "---- y:" << channel->GetY() <<endl; 00166 cout << "---- chip id:" << chipId <<endl; 00167 #endif 00168 00169 } 00170 } 00171 } 00172 } 00173 00174 00175 00176 //===========================================================================* 00177 // Now analyze scurves 00178 //===========================================================================* 00179 00180 TString t_scurve_file = root_directory + root_filename ;//+ ".root"; 00181 //cout << t_scurve_file<< endl; 00182 00183 int tmin = 0; 00184 int tmed = 0; 00185 int tmax = 0; 00186 int delta = 0; 00187 00188 // TH1I * hscurve[64]; 00189 00190 /* TH1F * hinflex = new TH1F("hinflex","",1024,0,1024); 00191 TH1F * hsigma = new TH1F("hsigma","",50,-2,2); 00192 00193 TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10); 00194 TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2); 00195 */ 00196 TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64); 00197 TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64); 00198 00199 TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64); 00200 TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64); 00201 00202 TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64); 00203 00204 // TString name; 00205 /* for (int j=0;j<64;j++) 00206 { name = Form("hscurve_%i",j); 00207 TH1I * h = (TH1I*)tf.Get(name); 00208 hscurve[j] = h; 00209 }*/ 00210 00211 float scurve_max = 0; 00212 float scurve_inflex = 0; 00213 float inflex = 0; 00214 float sigma = 0; 00215 float inflex_err = 0; 00216 float sigma_err = 0; 00217 float chi2 = 0; 00218 float ndf = 0; 00219 Double_t deriv = 0; 00220 00221 bool start_100 ; // test if a channel starts at 100% else warning 00222 bool has_flipped ; // test if each channel flips from 100% to 0% else, warning message 00223 Int_t noisy; // test if each channel is noisy, ie has entry beyond flip point 00224 Int_t noise_limit = 2; // number of entry accepted 00225 Int_t noise_window = 10 ;// number of bin defore dac_max to verify noise 00226 Int_t flip_limit = 30 ; 00227 00228 for (int j=0;j<64;j++) 00229 { noisy = 0 ; 00230 for (int k=dac_min+1;k<=dac_max;k++) // ATTN : bin = dac + 1 !! 00231 { deriv = ( hscurve[j]->GetBinContent(k)) - (hscurve[j]->GetBinContent(k+1)); 00232 //cout << j << " " << k << " " << hscurve[j]->GetBinContent(k) << " " << hscurve[j]->GetBinContent(k+1) << " " << deriv << endl ; 00233 // test if each channel has flipped 00234 if (k==dac_min+1) // first loop 00235 { if (hscurve[j]->GetBinContent(k)>=100) 00236 start_100 = 1 ; 00237 else start_100 = 0 ; 00238 } 00239 00240 if (k>(dac_max - noise_window)) // test noisy channel 00241 { if (hscurve[j]->GetBinContent(k)>=noise_limit) 00242 noisy++ ; 00243 //cout << "Channel " << j << " dac " << k << " noise " << noisy << endl ; 00244 } 00245 00246 00247 00248 if (k==dac_max-1) // last loop has flipped ? 00249 { if (hscurve[j]->GetBinContent(k)<=flip_limit && start_100 == 1 ) 00250 { has_flipped = 1 ; 00251 //if (hscurve[j]->GetBinContent(k)>0) 00252 // noisy++;// = 0 ; 00253 //else noisy = 1 ; 00254 } 00255 else has_flipped = 0 ; 00256 } 00257 if (deriv > 0.) 00258 { hscurve_deriv[j]->SetBinContent(k,deriv); 00259 hscurve_deriv[j]->SetEntries((hscurve_deriv[j]->GetEntries())+deriv); 00260 //hscurve_deriv[j]->Fill(k,deriv); 00261 } 00262 } 00263 if (has_flipped == 0) 00264 { cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " has not flipped" ; 00265 if (start_100 == 1) cout << " (stuck to 100 %)" << endl ; 00266 else cout << " (stuck to 0 %)" << endl ; 00267 } 00268 if (noisy==noise_window) cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " may be noisy, check root scurve file !" << endl ; 00269 00270 } 00271 00272 // FIT 00273 TF1 * gaussfunc[64]; 00274 for (int j=0;j<64;j++) 00275 { name = Form("g_%i",j); 00276 gaussfunc[j] = new TF1(name,"gaus",0,1024); 00277 } 00278 00279 for (int j=0;j<64;j++) 00280 { //gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; 00281 gaussfunc[j]->SetParameter(0, hscurve_deriv[j]->GetMaximumBin()); 00282 gaussfunc[j]->SetParLimits(0, 1, 100); 00283 gaussfunc[j]->SetParameter(1, hscurve_deriv[j]->GetMean() ); //GetMaximumBin() 00284 //gaussfunc[j]->FixParameter(1, hscurve_deriv[j]->GetMean() ); 00285 gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; //hscurve_deriv[j]->GetMaximumBin()-10, hscurve_deriv[j]->GetMaximumBin() + 10) ; 00286 gaussfunc[j]->SetParameter(2, hscurve_deriv[j]->GetRMS() ); 00287 gaussfunc[j]->SetParLimits(2, 0.01, hscurve_deriv[j]->GetRMS()) ; 00288 // if (j==48) 00289 // hscurve_deriv[j]->Fit(gaussfunc[j],"BV", "",hscurve_deriv[j]->GetMean()-15, hscurve_deriv[j]->GetMean() + 15 ); 00290 // else 00291 hscurve_deriv[j]->Fit(gaussfunc[j],"WBQ", "",hscurve_deriv[j]->GetMean()-20, hscurve_deriv[j]->GetMean() + 20 ); 00292 //hscurve_deriv[j]->Fit(gaussfunc[j],"Q"); 00293 inflex = gaussfunc[j]->GetParameter(1); 00294 if (inflex<=dac_min || inflex>=dac_max) 00295 { cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << ", fit may be not correct ! Inflexion: " << inflex << " Dac min: " << dac_min << " Dac max: " << dac_max << endl ; 00296 } 00297 sigma = gaussfunc[j]->GetParameter(2); 00298 inflex_err = gaussfunc[j]->GetParError(1); 00299 sigma_err =gaussfunc[j]->GetParError(2); 00300 00301 hc_inflex->SetBinContent(j+1,inflex); 00302 hc_sigma->SetBinContent(j+1,sigma); 00303 00304 hc_inflex_err->SetBinContent(j+1,inflex_err); 00305 hc_sigma_err->SetBinContent(j+1,sigma_err); 00306 00307 //cout << "Channel " << j << " init at " << hscurve_deriv[j]->GetMean() << " " << gaussfunc[j]->GetParameter(1) << endl ; 00308 // cout << j << " " << inflex << " " << sigma << endl ; 00309 } 00310 00311 00312 ////////////////////////////////////////////////////////////////////////////////////////// 00313 // Write Scurves 00314 ////////////////////////////////////////////////////////////////////////////////////////// 00315 TFile * tf = new TFile(scurve_file,"RECREATE"); 00316 hvalue->Write(); 00317 hnpulse_threshold->Write(); 00318 for (int j=0;j<64;j++) 00319 { hscurve[j]->Write(); 00320 hscurve_deriv[j]->Write(); 00321 //gaussfunc[j]->Write(); 00322 } 00323 tf->Close(); 00324 00325 00326 #ifdef WRITE 00327 root_filename.ReplaceAll(".root", ""); 00328 //TString outfile = root_directory + "../RESULTS/" + "parameters_" + root_filename + ".txt"; 00329 TString outfile = root_directory + "../RESULTS/" + root_filename + ".txt"; 00330 ofstream out ; 00331 out.open(outfile) ; 00332 for (int j=0;j<64;j++) 00333 { 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; 00334 } 00335 #endif 00336 #ifdef PLOT 00337 TApplication theApp("App", &argc, argv); 00338 TString adj; 00339 if (dac==0){adj = "Low";} 00340 if (dac==1){adj = "Medium";} 00341 if (dac==2){adj = "High";} 00342 00343 c1->Divide(8,8); 00344 //c1->cd(); 00345 TString htitle = adj + " threshold"; 00346 for (int j=0;j<64;j++) 00347 { //c0->cd(j/8,j%8); 00348 c1->cd(j+1); 00349 hscurve[j]->SetLineWidth(2); 00350 hscurve[j]->Draw(""); 00351 c1->Update() ; 00352 } 00353 theApp.Run(kTRUE) ; 00354 c1->Update() ; 00355 #endif 00356 return 0 ; 00357 }
1.4.7