Changeset 293 in Idarraga
- Timestamp:
- May 1, 2012, 3:07:05 PM (12 years ago)
- Location:
- mafalda
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
mafalda/MAFTools/MAFTools.cpp
r288 r293 539 539 } 540 540 541 542 TimePixCalibration::TimePixCalibration(const char * af, 543 const char * bf, 544 const char * cf, 545 const char * tf, 546 int width, int height){ 547 548 m_width = width; 549 m_height = height; 550 551 ifstream afs(af, fstream::in); 552 ifstream bfs(bf, fstream::in); 553 ifstream cfs(cf, fstream::in); 554 ifstream tfs(tf, fstream::in); 555 556 //////////////////////////////////// 557 // a 558 m_a = new double[m_width*m_height]; 559 double temp = 0.; 560 int i = 0; 561 while ( afs.good() ) { 562 afs >> temp; 563 if(afs.eof()) break; 564 m_a[i++] = temp; 565 } 566 afs.close(); 567 if(i != m_width*m_height){ 568 cout << "Calibration does not seem to be complete for paremeter a." << endl; 569 cout << "Got " << i << " items. " << m_width << "*" << m_height 570 << " were requested. Giving up." << endl; 571 exit(1); 572 } 573 574 //////////////////////////////////// 575 // b 576 m_b = new double[m_width*m_height]; 577 i = 0; 578 while ( bfs.good() ) { 579 bfs >> temp; 580 if(bfs.eof()) break; 581 m_b[i++] = temp; 582 } 583 bfs.close(); 584 if(i != m_width*m_height){ 585 cout << "Calibration does not seem to be complete for paremeter b." << endl; 586 cout << "Got " << i << " items. " << m_width << "*" << m_height 587 << " were requested. Giving up." << endl; 588 exit(1); 589 } 590 591 //////////////////////////////////// 592 // c 593 m_c = new double[m_width*m_height]; 594 i = 0; 595 while ( cfs.good() ) { 596 cfs >> temp; 597 if(cfs.eof()) break; 598 m_c[i++] = temp; 599 } 600 cfs.close(); 601 if(i != m_width*m_height){ 602 cout << "Calibration does not seem to be complete for paremeter c." << endl; 603 cout << "Got " << i << " items. " << m_width << "*" << m_height 604 << " were requested. Giving up." << endl; 605 exit(1); 606 } 607 608 //////////////////////////////////// 609 // t 610 m_t = new double[m_width*m_height]; 611 i = 0; 612 while ( tfs.good() ) { 613 tfs >> temp; 614 if(tfs.eof()) break; 615 m_t[i++] = temp; 616 } 617 tfs.close(); 618 if(i != m_width*m_height){ 619 cout << "Calibration does not seem to be complete for paremeter t." << endl; 620 cout << "Got " << i << " items. " << m_width << "*" << m_height 621 << " were requested. Giving up." << endl; 622 exit(1); 623 } 624 625 } 626 627 TimePixCalibration::~TimePixCalibration(){ 628 } 629 630 double TimePixCalibration::GetE(pair<int,int> pix, int tot){ 631 632 //TF1 * surr = GetSurrogateTF1(pix); 633 634 int index = MAFTools::XYtoC(pix, m_width); 635 double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] }; 636 637 /* 638 if(pix.first == 111 && pix.second == 70) { 639 std::cout << "a = " << par[0] << std::endl; 640 std::cout << "b = " << par[1] << std::endl; 641 std::cout << "c = " << par[2] << std::endl; 642 std::cout << "t = " << par[3] << std::endl; 643 }*/ 644 645 646 double evalTOT = 0.; 647 double prev_e = 0.; 648 double dist = 0.; 649 double prev_dist = 10000; 650 651 // search 652 for(double e = 3.0; e <= 10000. ; e+=0.05) { // steps of 0.05 keV, start search 653 654 //evalTOT = surr->Eval(e); 655 evalTOT = SurrogateCalibFunction(e, par); 656 if(evalTOT < 0) continue; // don't consider those values. Not physical. 657 dist = TMath::Abs(evalTOT - tot); 658 //if(pix.first == 111 && pix.second == 70) cout << "dist = " << dist << " | e = " << e << " | evalTOT = " << evalTOT << " | tot = " << tot << endl; 659 if( dist > prev_dist ) break; // if I am going far, break. 660 661 prev_dist = dist; 662 prev_e = e; // keep the last energy 663 } 664 665 return prev_e; 666 } 667 668 int TimePixCalibration::GetTOT (pair<int, int> pix, double E) { 669 670 int index = MAFTools::XYtoC(pix, m_width); 671 double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] }; 672 673 return SurrogateCalibFunction(E, par); 674 //return GetSurrogateTF1(pix)->Eval(E); 675 } 676 677 double TimePixCalibration::SurrogateCalibFunction(double x, double * par){ 678 return par[0]*x + par[1] - (par[2]/(x-par[3])); 679 } 680 681 // Only if needed one can fetch a few functions as TF1 objects 682 TF1 * TimePixCalibration::GetSurrogateTF1(pair<int, int> pix){ 683 684 // check first if the function has been defined. Use the Set. 685 m_surrogateItr = m_surrogateSet.find(pix); 686 687 // If surrogate TF1 does not exist yet, build the function first 688 if( m_surrogateItr == m_surrogateSet.end() ) { 689 TString ftitle = "surrogate_"; 690 ftitle += pix.first; ftitle += "_"; ftitle += pix.second; 691 // Usually the parameters come given in keV !!! WARNING !!! 692 m_surrogateMap[pix] = new TF1(ftitle,"[0]*x + [1] - ([2]/(x-[3]))",0,1000); // 1 MeV 693 // Populate parameters. Look for right index first 694 int index = MAFTools::XYtoC(pix, m_width); 695 m_surrogateMap[pix]->SetParameters(m_a[index], m_b[index], m_c[index], m_t[index]); 696 } 697 698 return m_surrogateMap[pix]; 699 } 700 701 TGraph2D * ConvertClusterToGraph2D(blob b, int div, int id){ 702 703 TString funcname = "TGraph2D_"; 704 funcname += id; 541 // 3D representation of a blob with z being the TOT 542 543 #define __minwidth_x_graph2D 3 544 #define __minwidth_y_graph2D 3 545 546 bool GoodCandidateGraph2D(blob b){ 547 548 if(b.bP.width_x < __minwidth_x_graph2D || 549 b.bP.width_y < __minwidth_y_graph2D 550 ) 551 { 552 return false; 553 } 554 555 return true; 556 } 557 558 TGraph2D * ConvertClusterToGraph2D(blob b, int div, long frameId){ 559 560 TString funcname = "TGraph2D_x"; 561 funcname += (int)TMath::Floor(b.bP.geoCenter_x); 562 funcname += "_y"; 563 funcname += (int)TMath::Floor(b.bP.geoCenter_y); 564 funcname += "_frame"; 565 funcname += frameId; 705 566 return ConvertClusterToGraph2D(b, div, funcname); 706 567 } … … 708 569 TGraph2D * ConvertClusterToGraph2D(blob b, int div, TString gname){ 709 570 571 // check first if this is a good candidate 572 if( !GoodCandidateGraph2D(b) ) { 573 return 0x0; 574 } 575 576 // Some definitions 710 577 vector<double> x; 711 578 vector<double> y; 712 579 vector<double> z; 713 list< pair < pair<int, int>, int > >::iterator i ;580 list< pair < pair<int, int>, int > >::iterator itr; 714 581 list< pair < pair<int, int>, int > > cld = b.GetClusterDescription(); 715 582 716 if(div == 1 || div == 0) { 717 for ( i = cld.begin() ; i != cld.end() ; i++ ) { 718 x.push_back((*i).first.first); 719 y.push_back((*i).first.second); 720 z.push_back((*i).second); // TOT 721 } 722 } 723 /* 724 else { 725 } 726 */ 583 // Whatever happens I need first to copy the cluster in a 584 // matrix of size (width_x)*(width_y) 585 586 int i, j; 587 588 int ** c; 589 int Nx = b.bP.width_x; // size of the original grid x 590 int Ny = b.bP.width_y; // size of the original grid y 591 c = new int * [Nx]; 592 for(int i = 0 ; i < Nx ; i++) { 593 c[i] = new int[Ny]; 594 } 595 // Fill with zeros 596 for(int i = 0 ; i < Nx ; i++) { 597 for(int j = 0 ; j < Ny ; j++) { 598 c[i][j] = 0; 599 } 600 } 601 // Find the left-bottom extreme 602 int ext_x = TMath::FloorNint(b.bP.geoCenter_x) - (Nx/2); 603 int ext_y = TMath::FloorNint(b.bP.geoCenter_y) - (Ny/2); 604 // Now copy the old cluster. If the value is not available we will have 0 605 for ( itr = cld.begin() ; itr != cld.end() ; itr++ ) { 606 // bring it back to i,j indexes in the matrix 607 i = (*itr).first.first - ext_x; 608 j = (*itr).first.second - ext_y; 609 c[i][j] = (*itr).second; // TOT 610 //cout << "i = " << i << ", j = " << j << " : " << c[i][j] << endl; 611 } 612 613 // Print the cluster in the original matrix size width_x * width_y with extra zeros 614 //PrintMatrix(c, Nx, Ny); 615 616 if(div == 0) { 617 618 // Fill the vectors for TGraph2D from the original copy 619 // of the matrix. 620 for(int i = 0; i < Nx ; i++) { 621 for(int j = 0 ; j < Ny ; j++) { 622 x.push_back( i ); 623 y.push_back( j ); 624 z.push_back( c[i][j] ); // TOT 625 } 626 } 627 628 629 } else { 630 631 int Sx = Nx; 632 int Sy = Ny; 633 634 int ** cm = c; // initialize with the old grid 635 int ** c_erase; int sox; 636 for (int gridItr = 0 ; gridItr < div ; gridItr++) { 637 // store point to the old grid to erase after 638 c_erase = cm; sox = Sx; 639 // initially Sx is the size of the old grid 640 cm = RefineGridOnce(cm, &Sx, &Sy); 641 // now Sx, Sy is the size of the new grid 642 //PrintMatrix(cm, Sx, Sy); 643 // erase old grid 644 for(int i = 0 ; i < sox ; i++){delete [] c_erase[i];} 645 delete [] c_erase; 646 } 647 648 // Fill the vectors for TGraph2D 649 for(int i = 0; i < Sx ; i++) { 650 for(int j = 0 ; j < Sy ; j++) { 651 x.push_back( i ); 652 y.push_back( j ); 653 z.push_back( cm[i][j] ); // TOT 654 } 655 } 656 } 657 727 658 TGraph2D * g2 = new TGraph2D((int)x.size(), &x[0], &y[0], &z[0]); 728 659 g2->SetName(gname); 729 660 661 // erase last matrix 662 730 663 return g2; 731 664 } 732 665 666 int ** RefineGridOnce(int ** c, int * prev_x, int * prev_y){ 667 668 // newNx and newNy bring the previous size of the grid Nx * Ny and are 669 // also used to return the new 670 671 // Bring the cluster to a matrix of size (2*width_x - 1)*(2*width_y - 1) 672 // prepare the new matrix 673 674 int ** cm; 675 int Sx = 2*(*prev_x) - 1; // size of the new grid x // new size ! 676 int Sy = 2*(*prev_y) - 1; // size of the new grid y 677 cm = new int * [Sx]; 678 for(int i = 0 ; i < Sx ; i++) { 679 cm[i] = new int[Sy]; 680 } 681 // Fill with zeros 682 for(int i = 0 ; i < Sx ; i++) { 683 for(int j = 0 ; j < Sy ; j++) { 684 cm[i][j] = 0; 685 } 686 } 687 // 1) Fill the normal nodes 688 for(int i = 0 ; i < Sx ; i+=2) { 689 for(int j = 0 ; j < Sy ; j+=2) { 690 cm[i][j] = c[i/2][j/2]; 691 } 692 } 693 // 2) Fill the inner with only 2 neighbors 694 for(int i = 1; i < Sx ; i+=2) { 695 for(int j = 0 ; j < Sy ; j+=2) { 696 cm[i][j] = (cm[i-1][j] + cm[i+1][j]) / 2; 697 } 698 } 699 for(int i = 0; i < Sx ; i+=2) { 700 for(int j = 1 ; j < Sy ; j+=2) { 701 cm[i][j] = (cm[i][j-1] + cm[i][j+1]) / 2; 702 } 703 } 704 // 3) Crossed values with 4 neighbors 705 for(int i = 1; i < Sx ; i+=2) { 706 for(int j = 1 ; j < Sy ; j+=2) { 707 cm[i][j] = (cm[i-1][j-1] + cm[i+1][j-1] + cm[i-1][j+1] + cm[i+1][j+1]) / 4; 708 } 709 } 710 711 // new grid size 712 *prev_x = Sx; 713 *prev_y = Sy; 714 715 return cm; 716 } 717 718 void PrintMatrix(int ** c, int Nx, int Ny){ 719 720 int val; 721 int nchars; 722 723 // find max value 724 int max = 0; 725 for(int i = 0 ; i < Nx ; i++) { 726 for(int j = 0 ; j < Ny ; j++) { 727 val = c[i][j]; 728 if(val > max) max = val; 729 } 730 } 731 // according to max value this is the number of 732 // caracters to print in each entry to make the 733 // matrix look good 734 nchars = TMath::FloorNint( TMath::Log10(max) ); 735 736 // extract each char and print properly 737 char valC[256]; 738 char uc = 'a'; 739 int uc_cntr = 0; 740 741 for(int j = Ny-1 ; j >= 0 ; j--) { // Watch the direction here ! 742 for(int i = 0 ; i < Nx ; i++) { // trying to match what you observe in the screen :) 743 744 val = c[i][j]; 745 sprintf(valC,"%d",val); 746 uc = 'a'; uc_cntr = 0; 747 while(uc != '\0'){ 748 uc = valC[uc_cntr++]; 749 } 750 uc_cntr--; // <-- number of characters in valC 751 // Now print nchars - uc_cntr blank spaces + valC + space 752 for(int k = 0 ; k <= nchars-uc_cntr ; k++) { printf(" "); } 753 printf("%s",valC); 754 printf(" "); 755 } 756 printf("\n"); 757 } 758 759 } 760 733 761 } 734 762 #endif -
mafalda/MAFTools/MAFTools.h
r288 r293 85 85 // More functions on Clustesr 86 86 TGraph2D * ConvertClusterToGraph2D(blob, int, TString); // an identifier, can be string or int 87 TGraph2D * ConvertClusterToGraph2D(blob, int, int); 88 89 /////////////////////////////////////////////// 90 // TimePix specific 91 92 class TimePixCalibration { 93 94 public: 95 // Provide the 4 calibration files, a,b,c,t 96 TimePixCalibration(){}; 97 TimePixCalibration(const char *, const char *, const char *, const char *, int, int); 98 ~TimePixCalibration(); 99 100 double SurrogateCalibFunction(double x, double *); 101 int GetTOT(pair<int, int>, double); 102 double GetE(pair<int, int>, int); 103 104 TF1 * GetSurrogateTF1(pair<int, int>); 105 106 private: 107 108 double * m_a; 109 double * m_b; 110 double * m_c; 111 double * m_t; 112 113 int m_width; 114 int m_height; 115 116 map<pair<int, int>, TF1* > m_surrogateMap; 117 set< pair<int, int> > m_surrogateSet; // verify content 118 set< pair<int, int> >::iterator m_surrogateItr; 119 120 }; 87 TGraph2D * ConvertClusterToGraph2D(blob, int, long); 88 bool GoodCandidateGraph2D(blob); 89 void PrintMatrix(int **, int, int); 90 int ** RefineGridOnce(int ** c, int * , int *); 121 91 122 92 } -
mafalda/MPXAlgo/MediPixAlgo.cpp
r289 r293 11 11 #include "MPXStoreGate/CandidateContainer.h" 12 12 #include "MAFTools/MAFTools.h" 13 14 #include "TLatex.h" 13 15 14 16 ClassImp(MediPixAlgo) … … 599 601 } 600 602 601 TCanvas * MediPixAlgo::DrawInSeparateWindow(vector<TGraph2D *> g){ 603 TCanvas * MediPixAlgo::DrawInSeparateWindow(vector<TGraph2D *> g , MSG::Level vl){ 604 605 if(g.empty()) { 606 Log << vl << "Nothing to draw ... return NULL pointer" << endreq; 607 return 0x0; 608 } 609 610 if(vl == MSG::DEBUG) { 611 vector<TGraph2D *>::iterator i = g.begin(); 612 Log << vl << "Attemping to draw in a separate canvas the following objects : " << endreq; 613 for( ; i != g.end() ; i++){ 614 Log << vl << (*i)->GetName() << endreq; 615 } 616 } 602 617 603 618 /////////////////////////////////////////////////////////////////// … … 613 628 c11->cd(cntr2+1); 614 629 (*itr)->Draw("tri1"); 630 (*itr)->GetXaxis()->SetLabelSize(0.1); 631 (*itr)->GetYaxis()->SetLabelSize(0.1); 632 (*itr)->GetZaxis()->SetLabelSize(0.08); 633 634 TLatex l1; 635 l1.SetTextSize(0.12); 636 l1.DrawLatex(0,0, (*itr)->GetName() ); 637 615 638 cntr2++; 616 639 } -
mafalda/MPXAlgo/MediPixAlgo.h
r289 r293 196 196 197 197 // Extra drawing 198 TCanvas * DrawInSeparateWindow(vector<TGraph2D *> );198 TCanvas * DrawInSeparateWindow(vector<TGraph2D *>, MSG::Level vl = MSG::INFO); 199 199 200 200 ///////////////////////////// -
mafalda/MPXStoreGate/MPXStoreGate.cpp
r25 r293 38 38 } 39 39 40 // This function is for the user. It delivers the numbers of objects 40 41 Int_t MediPixStoreGate::GetNObjWithAuthor(TString authorName){ 41 42 return (Int_t)gateMap[authorName].size();
Note: See TracChangeset
for help on using the changeset viewer.