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