#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/MTTrack.hh"#include "root/MTChamber.hh"#include "root/MTDetector.hh"#include "root/MTBoard.hh"#include "root/MTChip.hh"#include "root/MTChannel.hh"#include "root/MTEvent.hh"#include "root/MTDif.hh"#include "root/MTHardroc1Chip.hh"#include "root/MTHardroc2Chip.hh"#include "root/MTDiracChip.hh"#include "root/MTMicrorocChip.hh"#include "MicroException.hh"Include dependency graph for test_temperature.cpp:
Go to the source code of this file.
Defines | |
| #define | WRITE |
| #define | CHIP_ID 33 |
Functions | |
| void | usage (char *progname) |
| int | main (int argc, char **argv) |
| #define WRITE |
Definition at line 50 of file test_temperature.cpp.
| #define CHIP_ID 33 |
| void usage | ( | char * | progname | ) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 67 of file test_temperature.cpp.
00067 { 00068 00069 /*****************************************************************/ 00070 /* ARGUMENTS TO THE PROGRAM **************************************/ 00071 /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */ 00072 /* 2 - NAME OF THE ROOT FILE *************************************/ 00073 /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/ 00074 /* 0 -> DAC_0 (LOW THRESHOLD) *********************************/ 00075 /* 1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/ 00076 /* 2 -> DAC_2 (HIGH THRESHOLD) ********************************/ 00077 /*****************************************************************/ 00078 if(argc != 4) 00079 { usage(argv[0]) ; 00080 return 1 ; 00081 } 00082 TString root_directory = argv[1]; 00083 TString root_filename = argv[2]; 00084 int dac = atoi(argv[3]); 00085 00086 #ifdef PLOT 00087 set_style_myplain() ; 00088 #endif 00089 00090 TCanvas *c1 = new TCanvas("c1","Test MICROROC",0,0,1000,1000); 00091 TString root_file = root_directory + root_filename; 00092 TString scurve_file = root_directory + "scurve_" + root_filename; 00093 int nbChannel = 0; 00094 float nbEvt = 0; 00095 int threshold = 0; 00096 //int chipid = 0; 00097 int hardid = 0; 00098 int value = 0; 00099 int serial = 0; 00100 //int content = 0; 00101 Double_t content = 0; 00102 UInt_t chipId ; 00103 int dac_min = 1025; 00104 int dac_max = 0; 00105 00106 TH1I * hvalue = new TH1I("hvalue","",5,0,5); 00107 TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024); 00108 00109 TH1I * hscurve[64] ; 00110 TH1I * hscurve_deriv[64] ; 00111 TString name, title, name_deriv, title_deriv ; 00112 for (int j=0;j<64;j++) 00113 { name = Form("hscurve_%i",j); 00114 title = Form("SCurve from channel %i",j); 00115 hscurve[j] = new TH1I(name,title,1024,0,1024); 00116 name_deriv = Form("deriv_%i",j); 00117 title_deriv = Form("Derivative from channel %i",j); 00118 hscurve_deriv[j] = new TH1I(name_deriv,title_deriv,1024,0,1024); 00119 } 00120 00121 TFile f(root_file); 00122 00123 int ttiter=0 ; 00124 00125 // from example file ! 00126 TIter nextkey(f.GetListOfKeys()); 00127 TKey *key; 00128 while (key = (TKey*)nextkey()) 00129 { 00130 TTree *tree = (TTree*)key->ReadObj(); 00131 MTRun *run=NULL; 00132 try 00133 { 00134 run=dynamic_cast<MTRun *>(tree->GetUserInfo()->First()); 00135 // cout << "---------------------------- NEW TTREE ---------------------------"<< endl; 00136 // ttiter++ ; 00137 // cout << ttiter << endl ; 00138 // nbEvt =0; 00139 // nbHit = 0; 00140 00141 bool showConfig = true; 00142 if ( run != NULL && showConfig) 00143 { 00144 //run->Info(); 00145 MTDetector* det = run->GetDetector(); 00146 // loop over chamber config 00147 00148 const std::map<UInt_t,MTChamber*>& chambers = det->GetChambers(); 00149 for( std::map<UInt_t,MTChamber*>::const_iterator ii=chambers.begin(); ii!=chambers.end(); ++ii) 00150 { 00151 const MTChamber& chamber= *((*ii).second); 00152 //cout << "chamber id: " << chamber.GetId() << endl; 00153 //cout << "chamber type: " << chamber.GetType() << endl; 00154 00155 const std::map<UInt_t,MTBoard*>& boards = chamber.GetBoards(); 00156 for( std::map<UInt_t,MTBoard*>::const_iterator ii=boards.begin(); ii!=boards.end(); ++ii) 00157 { 00158 const MTBoard& board= *((*ii).second); 00159 //board.Info(); 00160 00161 //cout << "\t\t\t**** Chips configuration ****" << endl ; 00162 00163 map<UInt_t,MTChip*> chips = board.GetChips(); 00164 for( map<UInt_t,MTChip*>::const_iterator iiChip=chips.begin(); iiChip!=chips.end(); ++iiChip) 00165 { 00166 const MTChip& chip = *((*iiChip).second); 00167 if ( chip.IsConfigured() ) 00168 { 00169 //cout << " -- New Chip configuration:" << chip.GetSoftId().toString() << " softId[" << chip.GetSoftId().GetValue() <<"], serial Number[" << chip.GetSerialNumber() << "]"<< endl <<endl <<endl ; 00170 if (chip.GetSerialNumber() == 14) { // serial 14 = id 32 -> see xml file 00171 try{ 00172 00173 const MTMicrorocChip µChip = dynamic_cast<const MTMicrorocChip &> (chip); 00174 00175 if (dac==0){threshold = microChip.GetThresholdDac_0();} 00176 if (dac==1){threshold = microChip.GetThresholdDac_1();} 00177 if (dac==2){threshold = microChip.GetThresholdDac_2();} 00178 //cout << "Chip: " << chip.GetSerialNumber() << " SoftId: " << chip.GetSoftId().toString() << " Threshold:" << threshold << " SoftId:" << endl ; 00179 } 00180 catch (...) {} 00181 00182 if (threshold<dac_min){dac_min = threshold;} 00183 if (threshold>dac_max){dac_max = threshold;} 00184 00185 MTEvent *evt = new MTEvent(); 00186 TBranch *branch= tree->GetBranch("MTEvent"); 00187 branch->SetAddress(&evt); 00188 MTChannel* channel = NULL; 00189 nbEvt = tree->GetEntries(); 00190 //cout << "Number of events in the tree: " << nbEvt << endl ; 00191 00192 // hnpulse_threshold->SetBinContent(threshold+1,nbEvt); 00193 00194 for( int evtNum = 0; evtNum < nbEvt ; evtNum++) 00195 { //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;} 00196 tree->GetEntry(evtNum); 00197 // if (chip.GetSerialNumber() == 14) { 00198 //ttiter++ ; 00199 //cout << "procees " << ttiter << endl ; 00200 nbChannel = evt->GetNchannel(); 00201 for(int i=0;i<nbChannel ;i++) 00202 { channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i); 00203 //if (channel->getChip()->GetSerialNumber() == 14) 00204 //cout << channel->>GetSerialNumber() << endl ; 00205 if (chipId = channel->GetChipId() == CHIP_ID) // numéro du chip 00206 {hardid = channel->GetHardId(); // numéro de la voie 00207 value = channel->GetDigitalValue(); 00208 hvalue->Fill(value); 00209 if (value==(dac+1)) // ie 1 for low, 2 for med and 3 for high .... 00210 { hscurve[hardid]->Fill(threshold); 00211 #ifdef R_DEBUG 00212 cout << "---- digital value:" << value <<endl; 00213 cout << "---- hard id:" << hardid <<endl; 00214 //cout << "---- x:"<< channel->GetX() <<endl; 00215 //cout << "---- y:" << channel->GetY() <<endl; 00216 cout << "---- chip id:" << chipId <<endl; 00217 #endif 00218 }} 00219 } 00220 } 00221 00222 00223 00224 00225 00226 } // fin serial == 14 00227 } 00228 } 00229 } 00230 } 00231 } 00232 00233 00234 00235 00236 00237 00238 00239 } 00240 catch (MicroException e) {} 00241 } 00242 00243 00244 // cout << "ttree iter " << ttiter << endl ; 00245 00246 //cout << " " << endl ; 00247 00248 00249 //===========================================================================* 00250 // Now analyze scurves 00251 //===========================================================================* 00252 00253 TString t_scurve_file = root_directory + root_filename ;//+ ".root"; 00254 //cout << t_scurve_file<< endl; 00255 00256 int tmin = 0; 00257 int tmed = 0; 00258 int tmax = 0; 00259 int delta = 0; 00260 00261 // TH1I * hscurve[64]; 00262 00263 /* TH1F * hinflex = new TH1F("hinflex","",1024,0,1024); 00264 TH1F * hsigma = new TH1F("hsigma","",50,-2,2); 00265 00266 TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10); 00267 TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2); 00268 */ 00269 TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64); 00270 TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64); 00271 00272 TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64); 00273 TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64); 00274 00275 TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64); 00276 00277 // TString name; 00278 /* for (int j=0;j<64;j++) 00279 { name = Form("hscurve_%i",j); 00280 TH1I * h = (TH1I*)tf.Get(name); 00281 hscurve[j] = h; 00282 }*/ 00283 00284 float scurve_max = 0; 00285 float scurve_inflex = 0; 00286 float inflex = 0; 00287 float sigma = 0; 00288 float inflex_err = 0; 00289 float sigma_err = 0; 00290 float chi2 = 0; 00291 float ndf = 0; 00292 Double_t deriv = 0; 00293 00294 bool start_100 ; // test if a channel starts at 100% else warning 00295 bool has_flipped ; // test if each channel flips from 100% to 0% else, warning message 00296 Int_t noisy; // test if each channel is noisy, ie has entry beyond flip point 00297 Int_t noise_limit = 2; // number of entry accepted 00298 Int_t noise_window = 10 ;// number of bin defore dac_max to verify noise 00299 Int_t flip_limit = 30 ; 00300 00301 for (int j=0;j<64;j++) 00302 { noisy = 0 ; 00303 for (int k=dac_min+1;k<=dac_max;k++) // ATTN : bin = dac + 1 !! 00304 { deriv = ( hscurve[j]->GetBinContent(k)) - (hscurve[j]->GetBinContent(k+1)); 00305 //cout << j << " " << k << " " << hscurve[j]->GetBinContent(k) << " " << hscurve[j]->GetBinContent(k+1) << " " << deriv << endl ; 00306 // test if each channel has flipped 00307 if (k==dac_min+1) // first loop 00308 { if (hscurve[j]->GetBinContent(k)>=100) 00309 start_100 = 1 ; 00310 else start_100 = 0 ; 00311 } 00312 00313 if (k>(dac_max - noise_window)) // test noisy channel 00314 { if (hscurve[j]->GetBinContent(k)>=noise_limit) 00315 noisy++ ; 00316 //cout << "Channel " << j << " dac " << k << " noise " << noisy << endl ; 00317 } 00318 00319 00320 00321 if (k==dac_max-1) // last loop has flipped ? 00322 { if (hscurve[j]->GetBinContent(k)<=flip_limit && start_100 == 1 ) 00323 { has_flipped = 1 ; 00324 //if (hscurve[j]->GetBinContent(k)>0) 00325 // noisy++;// = 0 ; 00326 //else noisy = 1 ; 00327 } 00328 else has_flipped = 0 ; 00329 } 00330 if (deriv > 0.) 00331 { hscurve_deriv[j]->SetBinContent(k,deriv); 00332 hscurve_deriv[j]->SetEntries((hscurve_deriv[j]->GetEntries())+deriv); 00333 //hscurve_deriv[j]->Fill(k,deriv); 00334 } 00335 } 00336 if (has_flipped == 0) 00337 { cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " has not flipped" ; 00338 if (start_100 == 1) cout << " (stuck to 100 %)" << endl ; 00339 else cout << " (stuck to 0 %)" << endl ; 00340 } 00341 if (noisy==noise_window) cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << " may be noisy, check root scurve file !" << endl ; 00342 00343 } 00344 00345 // FIT 00346 TF1 * gaussfunc[64]; 00347 for (int j=0;j<64;j++) 00348 { name = Form("g_%i",j); 00349 gaussfunc[j] = new TF1(name,"gaus",0,1024); 00350 } 00351 00352 for (int j=0;j<64;j++) 00353 { //gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; 00354 gaussfunc[j]->SetParameter(0, hscurve_deriv[j]->GetMaximumBin()); 00355 gaussfunc[j]->SetParLimits(0, 1, 100); 00356 gaussfunc[j]->SetParameter(1, hscurve_deriv[j]->GetMean() ); //GetMaximumBin() 00357 //gaussfunc[j]->FixParameter(1, hscurve_deriv[j]->GetMean() ); 00358 gaussfunc[j]->SetParLimits(1, dac_min, dac_max) ; //hscurve_deriv[j]->GetMaximumBin()-10, hscurve_deriv[j]->GetMaximumBin() + 10) ; 00359 gaussfunc[j]->SetParameter(2, hscurve_deriv[j]->GetRMS() ); 00360 gaussfunc[j]->SetParLimits(2, 0.01, hscurve_deriv[j]->GetRMS()) ; 00361 // if (j==48) 00362 // hscurve_deriv[j]->Fit(gaussfunc[j],"BV", "",hscurve_deriv[j]->GetMean()-15, hscurve_deriv[j]->GetMean() + 15 ); 00363 // else 00364 hscurve_deriv[j]->Fit(gaussfunc[j],"WBQ", "",hscurve_deriv[j]->GetMean()-20, hscurve_deriv[j]->GetMean() + 20 ); 00365 //hscurve_deriv[j]->Fit(gaussfunc[j],"Q"); 00366 inflex = gaussfunc[j]->GetParameter(1); 00367 if (inflex<=dac_min || inflex>=dac_max) 00368 { cout << "****** " << argv[0] << ": Warning! Chip " << serial << " Channel " << j << ", fit may be not correct ! Inflexion: " << inflex << " Dac min: " << dac_min << " Dac max: " << dac_max << endl ; 00369 } 00370 sigma = gaussfunc[j]->GetParameter(2); 00371 inflex_err = gaussfunc[j]->GetParError(1); 00372 sigma_err =gaussfunc[j]->GetParError(2); 00373 00374 hc_inflex->SetBinContent(j+1,inflex); 00375 hc_sigma->SetBinContent(j+1,sigma); 00376 00377 hc_inflex_err->SetBinContent(j+1,inflex_err); 00378 hc_sigma_err->SetBinContent(j+1,sigma_err); 00379 00380 //cout << "Channel " << j << " init at " << hscurve_deriv[j]->GetMean() << " " << gaussfunc[j]->GetParameter(1) << endl ; 00381 // cout << j << " " << inflex << " " << sigma << endl ; 00382 } 00383 00384 00385 ////////////////////////////////////////////////////////////////////////////////////////// 00386 // Write Scurves 00387 ////////////////////////////////////////////////////////////////////////////////////////// 00388 TFile * tf = new TFile(scurve_file,"RECREATE"); 00389 hvalue->Write(); 00390 hnpulse_threshold->Write(); 00391 for (int j=0;j<64;j++) 00392 { hscurve[j]->Write(); 00393 hscurve_deriv[j]->Write(); 00394 //gaussfunc[j]->Write(); 00395 } 00396 tf->Close(); 00397 00398 00399 #ifdef WRITE 00400 root_filename.ReplaceAll(".root", ""); 00401 //TString outfile = root_directory + "../RESULTS/" + "parameters_" + root_filename + ".txt"; 00402 TString outfile = root_directory + "../RESULTS/" + root_filename + ".txt"; 00403 ofstream out ; 00404 out.open(outfile) ; 00405 for (int j=0;j<64;j++) 00406 { 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; 00407 } 00408 #endif 00409 #ifdef PLOT 00410 TApplication theApp("App", &argc, argv); 00411 TString adj; 00412 if (dac==0){adj = "Low";} 00413 if (dac==1){adj = "Medium";} 00414 if (dac==2){adj = "High";} 00415 00416 c1->Divide(8,8); 00417 //c1->cd(); 00418 TString htitle = adj + " threshold"; 00419 for (int j=0;j<64;j++) 00420 { //c0->cd(j/8,j%8); 00421 c1->cd(j+1); 00422 hscurve[j]->SetLineWidth(2); 00423 hscurve[j]->Draw(""); 00424 c1->Update() ; 00425 } 00426 theApp.Run(kTRUE) ; 00427 c1->Update() ; 00428 #endif 00429 return 0 ; 00430 }
1.4.7