#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 "TPaveStats.h"#include "TLine.h"#include "TPaveLabel.h"#include "root/MTRun.hh"#include "root/MTEvent.hh"#include "root/MTChannel.hh"#include "root/MTMicrorocChip.hh"#include "MicroException.hh"Include dependency graph for fit_scurve_chip.cpp:
Go to the source code of this file.
Defines | |
| #define | PLOT |
| #define | INFORENO |
Functions | |
| void | usage (char *progname) |
| void | set_style_myplain () |
| int | main (int argc, char **argv) |
Variables | |
| TStyle * | style_myplain = new TStyle("myplain","classic style") |
| #define PLOT |
Definition at line 42 of file fit_scurve_chip.cpp.
| #define INFORENO |
Definition at line 46 of file fit_scurve_chip.cpp.
| void usage | ( | char * | progname | ) |
| void set_style_myplain | ( | ) |
Definition at line 402 of file fit_scurve_chip.cpp.
Referenced by main().
00403 { 00404 style_myplain->SetCanvasBorderMode(0); 00405 style_myplain->SetCanvasColor(0); 00406 style_myplain->SetDrawBorder(0); 00407 style_myplain->SetPadBorderMode(0); 00408 style_myplain->SetPadColor(10); 00409 style_myplain->SetFrameLineColor(1); 00410 style_myplain->SetFrameFillColor(5); 00411 style_myplain->SetFrameFillStyle(0); 00412 style_myplain->SetFrameBorderMode(0); 00413 00414 //style_myplain->SetFrameLineWidth(5); //test for setframelinecolor ! 00415 00416 style_myplain->SetLegendBorderSize(1); 00417 00418 style_myplain->SetTextColor(231); 00419 00420 // style_myplain->SetFillColor(0); // If enable : impossible to fill histo with color !!! see note in the end ! 00421 style_myplain->SetLineColor(231); 00422 00423 style_myplain->SetStatColor(10); 00424 style_myplain->SetStatTextColor(1); 00425 style_myplain->SetStatBorderSize(1) ; 00426 style_myplain->SetStatX(0.9); 00427 style_myplain->SetStatY(0.9); 00428 style_myplain->SetStatW(0.1); 00429 style_myplain->SetStatH(0.1); 00430 style_myplain->SetOptStat("");//"emr" by default 00431 style_myplain->SetOptFit(0111); 00432 00433 style_myplain->SetTitleTextColor(1); 00434 style_myplain->SetTitleBorderSize(0); 00435 style_myplain->SetTitleAlign(23); 00436 style_myplain->SetTitleX(0.5); 00437 00438 style_myplain->SetAxisColor(1,"x"); // bleu 00439 style_myplain->SetAxisColor(1,"y"); 00440 style_myplain->SetAxisColor(1,"z"); 00441 style_myplain->SetLabelColor(1,"x"); 00442 style_myplain->SetLabelColor(1,"x"); 00443 style_myplain->SetLabelColor(1,"y"); 00444 style_myplain->SetLabelColor(1,"z"); 00445 style_myplain->SetLabelSize(0.03,"x"); 00446 style_myplain->SetLabelSize(0.03,"y"); 00447 style_myplain->SetLabelSize(0.03,"z"); 00448 style_myplain->SetTitleColor(1,"x"); 00449 style_myplain->SetTitleColor(1,"y"); 00450 style_myplain->SetTitleColor(1,"z"); 00451 style_myplain->SetTitleFillColor(10); 00452 00453 style_myplain->SetHistFillColor(10); 00454 //style_myplain->SetHistFillStyle(0); 00455 style_myplain->SetHistFillStyle(1001); 00456 style_myplain->SetHistLineColor(1); // 231 !!! 00457 00458 style_myplain->SetFuncColor(kOrange+10); // dark orange 00459 style_myplain->SetPalette(1); // temp�rature like 00460 00461 // style_myplain->SetPaperSize(100,100); // for Hi Res plots only ! 00462 00463 style_myplain->cd(); 00464 00465 gROOT->SetStyle("myplain"); 00466 00467 gROOT->ForceStyle(true); 00468 00469 00470 00471 // frame style not applied ! why ???????? 00472 // do it manually : right click on frame then "use current style" 00473 }
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 57 of file fit_scurve_chip.cpp.
00057 { 00058 00059 /*****************************************************************/ 00060 /* ARGUMENTS TO THE PROGRAM **************************************/ 00061 /* 1 - DIRECTORY OF THE ROOT FILE PRODUCED BY THE RECONSTRUCTION */ 00062 /* 2 - NAME OF THE ROOT FILE *************************************/ 00063 /* 3 - THRESHOLD THAT IS MOVED DURING THE CALIBRATION ************/ 00064 /* 0 -> DAC_0 (LOW THRESHOLD) *********************************/ 00065 /* 1 -> DAC_1 (MEDIUM THRESHOLD) ******************************/ 00066 /* 2 -> DAC_2 (HIGH THRESHOLD) ********************************/ 00067 /*****************************************************************/ 00068 if(argc != 5) 00069 { usage(argv[0]) ; 00070 return 1 ; 00071 } 00072 TString root_directory = argv[1]; 00073 TString root_filename = argv[2]; 00074 int dac = atoi(argv[3]); 00075 int chan = atoi(argv[4]); // chip channel 00076 00077 //bool display = 0; 00078 /* #ifdef PLOT 00079 set_style_myplain() ; 00080 #endif 00081 */ 00082 TString root_file = root_directory + root_filename; 00083 TString scurve_file = root_directory + "scurve_" + root_filename; 00084 00085 // cout<<"Wish to write "<<scurve_file<<endl; 00086 00087 int nbChannel = 0; 00088 float nbEvt = 0; 00089 int threshold = 0; 00090 //int chipid = 0; 00091 int hardid = 0; 00092 int value = 0; 00093 //int content = 0; 00094 Double_t content = 0; 00095 UInt_t chipId ; 00096 00097 // TH1I * hvalue = new TH1I("hvalue","",5,0,5); 00098 // TH1I * hnpulse_threshold = new TH1I("hnpulse_threshold","",1024,0,1024); 00099 00100 // TH1I * hscurve[64]; 00101 TString name,title; 00102 00103 // for (int j=0;j<64;j++) 00104 // { name = Form("hscurve_%i",j); 00105 // title = Form("SCurve from channel %i",j); 00106 // hscurve[j] = new TH1I(name,title,1024,0,1024); 00107 // } 00108 00109 name = Form("hscurve_%i", chan); 00110 title = Form("SCurve from channel %i", chan); 00111 TH1I * hscurve = new TH1I(name,"",1025,0,1025); 00112 00113 TFile f(root_file); 00114 TIter nextkey(f.GetListOfKeys()); 00115 TKey *key; 00116 #ifdef INFORENO 00117 cout << "Reading tree" << endl ; 00118 #endif 00119 while (key = (TKey*)nextkey()) 00120 { 00121 TTree *tree = (TTree*)key->ReadObj(); 00122 MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun"); 00123 const MTChip& chip = run->GetOneChip(); 00124 00125 try 00126 { const MTMicrorocChip µChip = dynamic_cast<const MTMicrorocChip &> (chip); 00127 //const MTHardroc2Chip µChip = dynamic_cast<const MTHardroc2Chip &> (chip); 00128 //cout << "threshold dac 0:[" << hr2Chip.GetThresholdDac_0() << endl; 00129 //cout << "threshold dac 1:[" << hr2Chip.GetThresholdDac_1() << endl; 00130 //cout << "threshold dac 2:[" << hr2Chip.GetThresholdDac_2() << endl; 00131 //chipId = chip.GetId() ; 00132 //cout << "Chip identifier: " << chipId << endl ; 00133 if (dac==0){threshold = microChip.GetThresholdDac_0();} 00134 if (dac==1){threshold = microChip.GetThresholdDac_1();} 00135 if (dac==2){threshold = microChip.GetThresholdDac_2();} 00136 } 00137 catch (...) {} 00138 00139 MTEvent *evt = new MTEvent(); 00140 TBranch *branch= tree->GetBranch("MTEvent"); 00141 branch->SetAddress(&evt); 00142 MTChannel* channel = NULL; 00143 nbEvt = tree->GetEntries(); 00144 // cout << "Number of events in the tree: " << nbEvt << endl ; 00145 00146 // hnpulse_threshold->SetBinContent(threshold+1,nbEvt); 00147 00148 for( int evtNum = 0; evtNum < nbEvt ; evtNum++) 00149 { //if (evtNum%1000 == 0) {cout<<evtNum/nbEvt*100.<<" %"<<"\r"<<flush;} 00150 tree->GetEntry(evtNum); 00151 nbChannel = evt->GetNchannel(); 00152 for(int i=0;i<nbChannel ;i++) 00153 { channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i); 00154 chipId = channel->GetChipId(); 00155 hardid = channel->GetHardId(); 00156 value = channel->GetDigitalValue(); 00157 // hvalue->Fill(value); 00158 if (value==(dac+1)) 00159 { hscurve->Fill(threshold); 00160 #ifdef R_DEBUG 00161 cout << "---- digital value:" << value <<endl; 00162 cout << "---- hard id:" << hardid <<endl; 00163 cout << "---- x:"<< channel->GetX() <<endl; 00164 cout << "---- y:" << channel->GetY() <<endl; 00165 cout << "---- chip id:" << chipId <<endl; 00166 #endif 00167 00168 } 00169 } 00170 } 00171 } 00172 #ifdef INFORENO 00173 cout << "Tree has been read !" << endl ; 00174 #endif 00175 00176 //===========================================================================* 00177 // Now analyze scurves 00178 #ifdef INFORENO 00179 cout << "Fitting" << endl ; 00180 #endif 00181 TString t_scurve_file = root_directory + root_filename ;//+ ".root"; 00182 //cout << t_scurve_file<< endl; 00183 00184 int tmin = 0; 00185 int tmed = 0; 00186 int tmax = 0; 00187 int delta = 0; 00188 00189 // TH1I * hscurve[64]; 00190 00191 TH1F * hinflex = new TH1F("hinflex","",1024,0,1024); 00192 TH1F * hsigma = new TH1F("hsigma","",50,-2,2); 00193 00194 TH1F * hinflex_err = new TH1F("hinflex_err","",100,0,10); 00195 TH1F * hsigma_err = new TH1F("hsigma_err","",100,0,2); 00196 00197 TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64); 00198 TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64); 00199 00200 TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64); 00201 TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64); 00202 00203 TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64); 00204 00205 // TString name; 00206 00207 /* for (int j=0;j<64;j++) 00208 { name = Form("hscurve_%i",j); 00209 TH1I * h = (TH1I*)tf.Get(name); 00210 hscurve[j] = h; 00211 } 00212 */ 00213 TString func = "[0]/(1+exp((x-[1])/[2]))"; 00214 TF1 * fermi[64]; 00215 float scurve_max = 0; 00216 float scurve_inflex = 0; 00217 float inflex = 0; 00218 float sigma = 0; 00219 float inflex_err = 0; 00220 float sigma_err = 0; 00221 float content0,content1,content2; 00222 float chi2 = 0; 00223 float ndf = 0; 00224 00225 00226 #ifdef PLOT 00227 set_style_myplain() ; 00228 TApplication theApp("App", &argc, argv); 00229 00230 // TCanvas *c = new TCanvas("c", "The Hello Canvas", 400, 400); 00231 00232 // TPaveLabel *hello = new TPaveLabel(0.2,0.4,0.8,0.6,"Hello World"); 00233 // hello->Draw(); 00234 // TPaveLabel *quit = new TPaveLabel(0.2,0.2,0.8,0.3,"Close via menu File/Quit"); 00235 // quit->Draw(); 00236 // c->Update(); 00237 00238 // Enter event loop, one can now interact with the objects in 00239 // the canvas. Select "Exit ROOT" from Canvas "File" menu to exit 00240 // the event loop and execute the next statements. 00241 // theApp.Run(kTRUE); 00242 00243 // TLine *l = new TLine(0.1,0.2,0.5,0.9); 00244 // l->Draw(); 00245 // quit->SetLabel("Select File/Quit again to exit app"); 00246 // c->Update(); 00247 00248 // Here we don't return from the eventloop. "Exit ROOT" will quit the app. 00249 // theApp.Run(); 00250 00251 00252 00253 00254 00255 00256 // set_style_myplain() ; 00257 TCanvas *c1 = new TCanvas("c1","Test MICROROC");//,0,0,1024,1024); 00258 // TApplication theApp("App", &argc, argv); 00259 // hscurve->Draw("AP"); 00260 // theApp.Run() ; 00261 #endif 00262 // for (int j=0;j<64;j++) 00263 // { 00264 int j = chan ; 00265 name = Form("fermi_%i",j); 00266 fermi[j] = new TF1(name,func,0,1024); 00267 //fermi[j]->SetLineColor(j+1); 00268 fermi[j]->SetLineColor(kRed); 00269 fermi[j]->SetLineWidth(2); 00270 /*If SCurve exist*/ 00271 if (hscurve->GetEntries()!=0) 00272 { /*Find threshold min tmin*/ 00273 for (int b=0;b<1023;b++) 00274 { content = hscurve->GetBinContent(b+1); 00275 if (content != 0) 00276 { tmin = b ; 00277 b = 1023 ; 00278 scurve_max = content; 00279 } 00280 } 00281 /*Find threshold at which scurve = 0 tmax*/ 00282 for (int b=1024;b>tmin;b--) 00283 { if (hscurve->GetBinContent(b+1) != 0) 00284 { tmax = b+1 ; 00285 b = tmin ; 00286 } 00287 } 00288 /*Find threshold at which scurve reaches the max*/ 00289 for (int b=tmax;b>=0;b--) 00290 { if (hscurve->GetBinContent(b+1) == scurve_max) 00291 { tmed = b+1 ; 00292 b = -1 ; 00293 } 00294 } 00295 scurve_inflex = 0.5*(tmax+tmed); 00296 delta = tmax - tmed; 00297 00298 /*If SCurve is a step function, do not fit parameters*/ 00299 if (delta==0) 00300 { inflex = scurve_inflex; 00301 sigma = 0.; 00302 inflex_err = 0.; 00303 sigma_err = 0.; 00304 } 00305 else 00306 { /*Set parameters*/ 00307 fermi[j]->FixParameter(0,scurve_max); 00308 fermi[j]->SetParameter(1,scurve_inflex); 00309 fermi[j]->SetParameter(2,delta/10.); 00310 /*Perform fit*/ 00311 hscurve->Fit(fermi[j],"q","",0,tmax+50); 00312 /*Get parameters*/ 00313 inflex = fermi[j]->GetParameter(1); 00314 sigma = sqrt(pow(fermi[j]->GetParameter(2),2)); 00315 inflex_err = sqrt(pow(fermi[j]->GetParError(1),2)); 00316 sigma_err = sqrt(pow(fermi[j]->GetParError(2),2)); 00317 chi2 = fermi[j]->GetChisquare(); 00318 ndf = fermi[j]->GetNDF(); 00319 } 00320 00321 hc_inflex_sigma->SetBinContent(j+1,inflex); 00322 hc_inflex_sigma->SetBinError(j+1,sigma); 00323 00324 hc_inflex->SetBinContent(j+1,inflex); 00325 hc_sigma->SetBinContent(j+1,sigma); 00326 00327 hc_inflex_err->SetBinContent(j+1,inflex_err); 00328 hc_sigma_err->SetBinContent(j+1,sigma_err); 00329 00330 hinflex->Fill(inflex); 00331 hsigma->Fill(sigma); 00332 00333 #ifdef DISPLAY 00334 cout<<" Channel\t" << j << "\t" << inflex << "\t width\t" << sigma << endl ; 00335 #endif 00336 } 00337 // } 00338 #ifdef WRITE 00339 root_filename.ReplaceAll(".root", ""); 00340 TString outfile = root_directory + "../RESULTS/" + "parameters_" + root_filename + ".txt"; 00341 ofstream out ; 00342 out.open(outfile) ; 00343 for (int j=0;j<64;j++) 00344 { 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; 00345 } 00346 #endif 00347 #ifdef PLOT 00348 #ifdef INFORENO 00349 cout << "Plotting" << endl ; 00350 #endif 00351 00352 TString adj; 00353 if (dac==0){adj = "Low";} 00354 if (dac==1){adj = "Medium";} 00355 if (dac==2){adj = "High";} 00356 c1->Connect("Closed()", "TApplication", &theApp, "Terminate()"); 00357 //c1->Divide(8,8); 00358 c1->cd(); 00359 hscurve->Draw(""); 00360 TString htitle = adj + " threshold"; 00361 /* for (int j=0;j<64;j++) 00362 { //c0->cd(j/8,j%8); 00363 c1->cd(j+1); 00364 hscurve->SetLineWidth(2); 00365 hscurve->Draw(""); 00366 c1->Update() ; 00367 } 00368 */ 00369 hscurve->GetXaxis()->SetTitle("Threshold [DAC]"); 00370 hscurve->GetYaxis()->SetTitle("N"); 00371 TPaveStats *p1 = (TPaveStats*) hscurve->GetListOfFunctions()->FindObject("stats"); 00372 00373 00374 c1->Update() ; 00375 #ifdef INFORENO 00376 cout << "Update..." << endl ; 00377 #endif 00378 theApp.Run() ; 00379 #ifdef INFORENO 00380 cout << "Do you see something ?" << endl ; 00381 #endif 00382 #endif 00383 return 0 ; 00384 }
| TStyle* style_myplain = new TStyle("myplain","classic style") |
1.4.7