Changeset 293 in Idarraga

Timestamp:

May 1, 2012, 3:07:05 PM (12 years ago)

Author:

idarraga

Message:

 

Location:

mafalda

Files:

• MAFTools/MAFTools.cpp (modified) (2 diffs)
• MAFTools/MAFTools.h (modified) (1 diff)
• MPXAlgo/MediPixAlgo.cpp (modified) (3 diffs)
• MPXAlgo/MediPixAlgo.h (modified) (1 diff)
• MPXStoreGate/MPXStoreGate.cpp (modified) (1 diff)

mafalda/MAFTools/MAFTools.cpp

@@ -539 +539 @@
 }
 
-
-TimePixCalibration::TimePixCalibration(const char * af,
-                const char * bf,
-                const char * cf,
-                const char * tf,
-                int width, int height){
-
-        m_width = width;
-        m_height = height;
-
-        ifstream afs(af, fstream::in);
-        ifstream bfs(bf, fstream::in);
-        ifstream cfs(cf, fstream::in);
-        ifstream tfs(tf, fstream::in);
-
-        ////////////////////////////////////
-        // a
-        m_a = new double[m_width*m_height];
-        double temp = 0.;
-        int i = 0;
-        while ( afs.good() ) {
-                afs >> temp;
-                if(afs.eof()) break;
-                m_a[i++] = temp;
-        }
-        afs.close();
-        if(i != m_width*m_height){
-                cout << "Calibration does not seem to be complete for paremeter a." << endl;
-                cout << "Got " << i << " items.  " << m_width << "*" << m_height
-                                << " were requested.  Giving up." << endl;
-                exit(1);
-        }
-
-        ////////////////////////////////////
-        // b
-        m_b = new double[m_width*m_height];
-        i = 0;
-        while ( bfs.good() ) {
-                bfs >> temp;
-                if(bfs.eof()) break;
-                m_b[i++] = temp;
-        }
-        bfs.close();
-        if(i != m_width*m_height){
-                cout << "Calibration does not seem to be complete for paremeter b." << endl;
-                cout << "Got " << i << " items.  " << m_width << "*" << m_height
-                                << " were requested.  Giving up." << endl;
-                exit(1);
-        }
-
-        ////////////////////////////////////
-        // c
-        m_c = new double[m_width*m_height];
-        i = 0;
-        while ( cfs.good() ) {
-                cfs >> temp;
-                if(cfs.eof()) break;
-                m_c[i++] = temp;
-        }
-        cfs.close();
-        if(i != m_width*m_height){
-                cout << "Calibration does not seem to be complete for paremeter c." << endl;
-                cout << "Got " << i << " items.  " << m_width << "*" << m_height
-                                << " were requested.  Giving up." << endl;
-                exit(1);
-        }
-
-        ////////////////////////////////////
-        // t
-        m_t = new double[m_width*m_height];
-        i = 0;
-        while ( tfs.good() ) {
-                tfs >> temp;
-                if(tfs.eof()) break;
-                m_t[i++] = temp;
-        }
-        tfs.close();
-        if(i != m_width*m_height){
-                cout << "Calibration does not seem to be complete for paremeter t." << endl;
-                cout << "Got " << i << " items.  " << m_width << "*" << m_height
-                                << " were requested.  Giving up." << endl;
-                exit(1);
-        }
-
-}
-
-TimePixCalibration::~TimePixCalibration(){
-}
-
-double TimePixCalibration::GetE(pair<int,int> pix, int tot){
-
-        //TF1 * surr = GetSurrogateTF1(pix);
-
-        int index = MAFTools::XYtoC(pix, m_width);
-        double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] };
-
-        /*
-        if(pix.first == 111 && pix.second == 70) {
-                std::cout << "a = " << par[0] << std::endl;
-                std::cout << "b = " << par[1] << std::endl;
-                std::cout << "c = " << par[2] << std::endl;
-                std::cout << "t = " << par[3] << std::endl;
-        }*/
-
-
-        double evalTOT = 0.;
-        double prev_e = 0.;
-        double dist = 0.;
-        double prev_dist = 10000;
-
-        // search
-        for(double e = 3.0; e <= 10000. ; e+=0.05) { // steps of 0.05 keV, start search
-
-                //evalTOT = surr->Eval(e);
-                evalTOT = SurrogateCalibFunction(e, par);
-                if(evalTOT < 0) continue; // don't consider those values. Not physical.
-                dist = TMath::Abs(evalTOT - tot);
-                //if(pix.first == 111 && pix.second == 70) cout << "dist = " << dist << " | e = " << e << " | evalTOT = " << evalTOT << " | tot = " << tot << endl;
-                if( dist > prev_dist ) break; // if I am going far, break.
-
-                prev_dist = dist;
-                prev_e = e; // keep the last energy
-        }
-
-        return prev_e;
-}
-
-int TimePixCalibration::GetTOT (pair<int, int> pix, double E) {
-
-        int index = MAFTools::XYtoC(pix, m_width);
-        double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] };
-
-        return SurrogateCalibFunction(E, par);
-        //return GetSurrogateTF1(pix)->Eval(E);
-}
-
-double TimePixCalibration::SurrogateCalibFunction(double x, double * par){
-        return par[0]*x + par[1] - (par[2]/(x-par[3]));
-}
-
-// Only if needed one can fetch a few functions as TF1 objects
-TF1 * TimePixCalibration::GetSurrogateTF1(pair<int, int> pix){
-
-        // check first if the function has been defined.  Use the Set.
-        m_surrogateItr = m_surrogateSet.find(pix);
-
-        // If surrogate TF1 does not exist yet, build the function first
-        if( m_surrogateItr == m_surrogateSet.end() ) {
-                TString ftitle = "surrogate_";
-                ftitle += pix.first; ftitle += "_"; ftitle += pix.second;
-                // Usually the parameters come given in keV !!! WARNING !!!
-                m_surrogateMap[pix] = new TF1(ftitle,"[0]*x + [1] - ([2]/(x-[3]))",0,1000); // 1 MeV
-                // Populate parameters.  Look for right index first
-                int index = MAFTools::XYtoC(pix, m_width);
-                m_surrogateMap[pix]->SetParameters(m_a[index], m_b[index], m_c[index], m_t[index]);
-        }
-
-        return m_surrogateMap[pix];
-}
-
-TGraph2D * ConvertClusterToGraph2D(blob b, int div, int id){
-
-        TString funcname = "TGraph2D_";
-        funcname += id;
+// 3D representation of a blob with z being the TOT
+
+#define __minwidth_x_graph2D 3
+#define __minwidth_y_graph2D 3
+
+bool GoodCandidateGraph2D(blob b){
+
+        if(b.bP.width_x < __minwidth_x_graph2D ||
+                        b.bP.width_y < __minwidth_y_graph2D
+        )
+        {
+                return false;
+        }
+
+        return true;
+}
+
+TGraph2D * ConvertClusterToGraph2D(blob b, int div, long frameId){
+
+        TString funcname = "TGraph2D_x";
+        funcname += (int)TMath::Floor(b.bP.geoCenter_x);
+        funcname += "_y";
+        funcname += (int)TMath::Floor(b.bP.geoCenter_y);
+        funcname += "_frame";
+        funcname += frameId;
         return ConvertClusterToGraph2D(b, div, funcname);
 }
@@ -708 +569 @@
 TGraph2D * ConvertClusterToGraph2D(blob b, int div, TString gname){
 
+        // check first if this is a good candidate
+        if( !GoodCandidateGraph2D(b) ) {
+                return 0x0;
+        }
+
+        // Some definitions
         vector<double> x;
         vector<double> y;
         vector<double> z;
-        list< pair < pair<int, int>, int > >::iterator i;
+        list< pair < pair<int, int>, int > >::iterator itr;
         list< pair < pair<int, int>, int > > cld = b.GetClusterDescription();
 
-        if(div == 1 || div == 0) {
-                for ( i = cld.begin() ; i != cld.end() ; i++ ) {
-                        x.push_back((*i).first.first);
-                        y.push_back((*i).first.second);
-                        z.push_back((*i).second); // TOT
-                }
-        }
-        /*
-        else {
-        }
-        */
+        // Whatever happens I need first to copy the cluster in a
+        // matrix of size (width_x)*(width_y)
+
+        int i, j;
+
+        int ** c;
+        int Nx = b.bP.width_x; // size of the original grid x
+        int Ny = b.bP.width_y; // size of the original grid y
+        c = new int * [Nx];
+        for(int i = 0 ; i < Nx ; i++) {
+                c[i] = new int[Ny];
+        }
+        // Fill with zeros
+        for(int i = 0 ; i < Nx ; i++) {
+                for(int j = 0 ; j < Ny ; j++) {
+                        c[i][j] = 0;
+                }
+        }
+        // Find the left-bottom extreme
+        int ext_x = TMath::FloorNint(b.bP.geoCenter_x) - (Nx/2);
+        int ext_y = TMath::FloorNint(b.bP.geoCenter_y) - (Ny/2);
+        // Now copy the old cluster.  If the value is not available we will have 0
+        for ( itr = cld.begin() ; itr != cld.end() ; itr++ ) {
+                // bring it back to i,j indexes in the matrix
+                i = (*itr).first.first - ext_x;
+                j = (*itr).first.second - ext_y;
+                c[i][j] = (*itr).second; // TOT
+                //cout << "i = " << i << ", j = " << j << " : " << c[i][j] << endl;
+        }
+
+        // Print the cluster in the original matrix size width_x * width_y with extra zeros
+        //PrintMatrix(c, Nx, Ny);
+
+        if(div == 0) {
+
+                // Fill the vectors for TGraph2D from the original copy
+                // of the matrix.
+                for(int i = 0; i < Nx ; i++) {
+                        for(int j = 0 ; j < Ny ; j++) {
+                                x.push_back( i );
+                                y.push_back( j );
+                                z.push_back( c[i][j] ); // TOT
+                        }
+                }
+
+
+        } else {
+
+                int Sx = Nx;
+                int Sy = Ny;
+
+                int ** cm = c; // initialize with the old grid
+                int ** c_erase; int sox;
+                for (int gridItr = 0 ; gridItr < div ; gridItr++) {
+                        // store point to the old grid to erase after
+                        c_erase = cm; sox = Sx;
+                        // initially Sx is the size of the old grid
+                        cm = RefineGridOnce(cm, &Sx, &Sy);
+                        // now Sx, Sy is the size of the new grid
+                        //PrintMatrix(cm, Sx, Sy);
+                        // erase old grid
+                        for(int i = 0 ; i < sox ; i++){delete [] c_erase[i];}
+                        delete [] c_erase;
+                }
+
+                // Fill the vectors for TGraph2D
+                for(int i = 0; i < Sx ; i++) {
+                        for(int j = 0 ; j < Sy ; j++) {
+                                x.push_back( i );
+                                y.push_back( j );
+                                z.push_back( cm[i][j] ); // TOT
+                        }
+                }
+        }
+
         TGraph2D * g2 = new TGraph2D((int)x.size(), &x[0], &y[0], &z[0]);
         g2->SetName(gname);
 
+        // erase last matrix
+
         return g2;
 }
 
+int ** RefineGridOnce(int ** c, int * prev_x, int * prev_y){
+
+        // newNx and newNy bring the previous size of the grid Nx * Ny and are
+        // also used to return the new
+
+        // Bring the cluster to a matrix of size (2*width_x - 1)*(2*width_y - 1)
+        // prepare the new matrix
+
+        int ** cm;
+        int Sx = 2*(*prev_x) - 1; // size of the new grid x // new size !
+        int Sy = 2*(*prev_y) - 1; // size of the new grid y
+        cm = new int * [Sx];
+        for(int i = 0 ; i < Sx ; i++) {
+                cm[i] = new int[Sy];
+        }
+        // Fill with zeros
+        for(int i = 0 ; i < Sx ; i++) {
+                for(int j = 0 ; j < Sy ; j++) {
+                        cm[i][j] = 0;
+                }
+        }
+        // 1) Fill the normal nodes
+        for(int i = 0 ; i < Sx ; i+=2) {
+                for(int j = 0 ; j < Sy ; j+=2) {
+                        cm[i][j] = c[i/2][j/2];
+                }
+        }
+        // 2) Fill the inner with only 2 neighbors
+        for(int i = 1; i < Sx ; i+=2) {
+                for(int j = 0 ; j < Sy ; j+=2) {
+                        cm[i][j] = (cm[i-1][j] + cm[i+1][j]) / 2;
+                }
+        }
+        for(int i = 0; i < Sx ; i+=2) {
+                for(int j = 1 ; j < Sy ; j+=2) {
+                        cm[i][j] = (cm[i][j-1] + cm[i][j+1]) / 2;
+                }
+        }
+        // 3) Crossed values with 4 neighbors
+        for(int i = 1; i < Sx ; i+=2) {
+                for(int j = 1 ; j < Sy ; j+=2) {
+                        cm[i][j] = (cm[i-1][j-1] + cm[i+1][j-1] + cm[i-1][j+1] + cm[i+1][j+1]) / 4;
+                }
+        }
+
+        // new grid size
+        *prev_x = Sx;
+        *prev_y = Sy;
+
+        return cm;
+}
+
+void PrintMatrix(int ** c, int Nx, int Ny){
+
+        int val;
+        int nchars;
+
+        // find max value
+        int max = 0;
+        for(int i = 0 ; i < Nx ; i++) {
+                for(int j = 0 ; j < Ny ; j++) {
+                        val = c[i][j];
+                        if(val > max) max = val;
+                }
+        }
+        // according to max value this is the number of
+        // caracters to print in each entry to make the
+        // matrix look good
+        nchars = TMath::FloorNint( TMath::Log10(max) );
+
+        // extract each char and print properly
+        char valC[256];
+        char uc = 'a';
+        int uc_cntr = 0;
+
+        for(int j = Ny-1 ; j >= 0 ; j--) {  // Watch the direction here !
+                for(int i = 0 ; i < Nx ; i++) { // trying to match what you observe in the screen :)
+
+                        val = c[i][j];
+                        sprintf(valC,"%d",val);
+                        uc = 'a'; uc_cntr = 0;
+                        while(uc != '\0'){
+                                uc = valC[uc_cntr++];
+                        }
+                        uc_cntr--; // <-- number of characters in valC
+                        // Now print nchars - uc_cntr blank spaces + valC + space
+                        for(int k = 0 ; k <= nchars-uc_cntr ; k++) { printf(" "); }
+                        printf("%s",valC);
+                        printf(" ");
+                }
+                printf("\n");
+        }
+
+}
+
 }
 #endif

mafalda/MAFTools/MAFTools.h

@@ -85 +85 @@
 // More functions on Clustesr
 TGraph2D * ConvertClusterToGraph2D(blob, int, TString); // an identifier, can be string or int
-TGraph2D * ConvertClusterToGraph2D(blob, int, int);
-
-///////////////////////////////////////////////
-// TimePix specific
-
-class TimePixCalibration {
-
-public:
-        // Provide the 4 calibration files, a,b,c,t
-        TimePixCalibration(){};
-        TimePixCalibration(const char *, const char *, const char *, const char *, int, int);
-        ~TimePixCalibration();
-
-        double SurrogateCalibFunction(double x, double *);
-        int GetTOT(pair<int, int>, double);
-        double GetE(pair<int, int>, int);
-
-        TF1 * GetSurrogateTF1(pair<int, int>);
-
-private:
-
-        double * m_a;
-        double * m_b;
-        double * m_c;
-        double * m_t;
-
-        int m_width;
-        int m_height;
-
-        map<pair<int, int>, TF1* > m_surrogateMap;
-        set< pair<int, int> > m_surrogateSet; // verify content
-        set< pair<int, int> >::iterator m_surrogateItr;
-
-};
+TGraph2D * ConvertClusterToGraph2D(blob, int, long);
+bool GoodCandidateGraph2D(blob);
+void PrintMatrix(int **, int, int);
+int ** RefineGridOnce(int ** c, int * , int *);
 
 }

mafalda/MPXAlgo/MediPixAlgo.cpp

@@ -11 +11 @@
 #include "MPXStoreGate/CandidateContainer.h"
 #include "MAFTools/MAFTools.h"
+
+#include "TLatex.h"
 
 ClassImp(MediPixAlgo)
@@ -599 +601 @@
 }
 
-TCanvas * MediPixAlgo::DrawInSeparateWindow(vector<TGraph2D *> g){
+TCanvas * MediPixAlgo::DrawInSeparateWindow(vector<TGraph2D *> g , MSG::Level vl){
+
+        if(g.empty()) {
+                Log << vl << "Nothing to draw ... return NULL pointer" << endreq;
+                return 0x0;
+        }
+
+        if(vl == MSG::DEBUG) {
+                vector<TGraph2D *>::iterator i = g.begin();
+                Log << vl << "Attemping to draw in a separate canvas the following objects : " << endreq;
+                for( ; i != g.end() ; i++){
+                        Log << vl << (*i)->GetName() << endreq;
+                }
+        }
 
         ///////////////////////////////////////////////////////////////////
@@ -613 +628 @@
                 c11->cd(cntr2+1);
                 (*itr)->Draw("tri1");
+                (*itr)->GetXaxis()->SetLabelSize(0.1);
+                (*itr)->GetYaxis()->SetLabelSize(0.1);
+                (*itr)->GetZaxis()->SetLabelSize(0.08);
+
+                TLatex l1;
+                l1.SetTextSize(0.12);
+                l1.DrawLatex(0,0, (*itr)->GetName() );
+
                 cntr2++;
         }

mafalda/MPXAlgo/MediPixAlgo.h

@@ -196 +196 @@
 
   // Extra drawing
-  TCanvas * DrawInSeparateWindow(vector<TGraph2D *>);
+  TCanvas * DrawInSeparateWindow(vector<TGraph2D *>,  MSG::Level vl = MSG::INFO);
 
   /////////////////////////////

mafalda/MPXStoreGate/MPXStoreGate.cpp

@@ -38 +38 @@
 }
 
+// This function is for the user.  It delivers the numbers of objects
 Int_t MediPixStoreGate::GetNObjWithAuthor(TString authorName){
   return (Int_t)gateMap[authorName].size();