Changeset 285 in Idarraga

Timestamp:

Apr 20, 2012, 5:11:44 PM (12 years ago)

Author:

idarraga

Message:

 

Location:

mafalda/MAFTools

Files:

• MAFTools.cpp (modified) (5 diffs)
• MAFTools.h (modified) (3 diffs)

mafalda/MAFTools/MAFTools.cpp

@@ -17 +17 @@
 #include <set>
 #include <iostream>
+#include <fstream>
+#include <stdlib.h>
 
 #include <TROOT.h>
@@ -349 +351 @@
 
                         //cout << "next: " <<  nextNeighbor.first << " , " << nextNeighbor.second << " ---> "
-                                        //<< totmap[ nextNeighbor.second * maxx + nextNeighbor.first ] << endl;
+                        //<< totmap[ nextNeighbor.second * maxx + nextNeighbor.first ] << endl;
 
                         // check the Set see if it existed
@@ -390 +392 @@
         else
                 return (2*tolerancePixels + 2)*4 +
-                GetLengthOfSpiral(--tolerancePixels);
+                                GetLengthOfSpiral(--tolerancePixels);
 
 }
@@ -437 +439 @@
         spiral->yCntr = 1;
         spiral->dir = 1;
+
+}
+
+void FillValuesForDisplay(Highlighter * hl, blob bl){
+
+        hl->UploadDisplayValue("Number of Pixels     ", bl.bP.nPixels);
+        hl->UploadDisplayValue("Inner Pixels         ", bl.bP.nInnerPixels);
+        hl->UploadDisplayValue("Total charge         ", bl.bP.totalCharge);
+        hl->UploadDisplayValue("Width x              ", bl.bP.width_x);
+        hl->UploadDisplayValue("Width y              ", bl.bP.width_y);
+        hl->UploadDisplayValue("Geo center x         ", bl.bP.geoCenter_x);
+        hl->UploadDisplayValue("Geo center y         ", bl.bP.geoCenter_y);
+        hl->UploadDisplayValue("Weighted center x    ", bl.bP.weightedCenter_x);
+        hl->UploadDisplayValue("Weighted center y    ", bl.bP.weightedCenter_y);
+        hl->UploadDisplayValue("Box area             ", bl.bP.boxArea);
+        hl->UploadDisplayValue("Circle area          ", bl.bP.circleArea);
+        hl->UploadDisplayValue("Ellipse area         ", bl.bP.ellipseArea);
+        hl->UploadDisplayValue("  circle/elipse      ", bl.bP.circleArea/bl.bP.ellipseArea);
+        hl->UploadDisplayValue("Elipse A             ", bl.bP.ellipseA);
+        hl->UploadDisplayValue("Elipse B             ", bl.bP.ellipseB);
+        hl->UploadDisplayValue("Rotation angle [deg] ", bl.bP.rotAngle * 180 / TMath::Pi());
+        hl->UploadDisplayValue("Chisquare/dof        ", bl.bP.chisquare_OverDof);
+        hl->UploadDisplayValue("Fit slope            ", bl.bP.fitSlope);
+        hl->UploadDisplayValue("Fit cut              ", bl.bP.fitCut);
+        hl->UploadDisplayValue("Balance to minimum   ", bl.bP.balanceToMin);
+        hl->UploadDisplayValue("Balance to maximum   ", bl.bP.balanceToMin);
+        hl->UploadDisplayValue("Min dist to center   ", bl.bP.minToGeoCenter);
+        hl->UploadDisplayValue("Max dist to center   ", bl.bP.maxToGeoCenter);
 
 }
@@ -507 +537 @@
 }
 
-}
-
-
+
+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];
+}
+
+}
 #endif

mafalda/MAFTools/MAFTools.h

@@ -7 +7 @@
 
 #include <set>
+#include <map>
 #include <iostream>
 #include <vector>
@@ -16 +17 @@
 //class blob;
 class MediPixAlgo;
+class Highlighter;
 //class AllBlobsContainer;
 //typedef struct steerSpiral;
@@ -22 +24 @@
 
 namespace MAFTools {
-  
-  Bool_t LinearRegression(Float_t * slope, Float_t * cut,
-                          Float_t * chisquare_OverDof, Float_t * chisquare, Float_t * ndf,
-                          Float_t initslope, Float_t initycross,
-                          std::set< std::pair<Int_t, Int_t> > blobContent);
 
-  Int_t XYtoC(pair<Int_t, Int_t>, Int_t);
-  Int_t XYtoC(Int_t, Int_t, Int_t);
+Bool_t LinearRegression(Float_t * slope, Float_t * cut,
+                Float_t * chisquare_OverDof, Float_t * chisquare, Float_t * ndf,
+                Float_t initslope, Float_t initycross,
+                std::set< std::pair<Int_t, Int_t> > blobContent);
 
-  Int_t ConvertXYPositionToInteger(std::pair<Int_t, Int_t> position, Int_t dimX, Int_t dimY);
+Int_t XYtoC(pair<Int_t, Int_t>, Int_t);
+Int_t XYtoC(Int_t, Int_t, Int_t);
 
-  Int_t ConvertXYPositionToInteger(Int_t positionX, Int_t positionY, Int_t dimX, Int_t dimY);
-  
-  Double_t CalcDistance(std::pair<Double_t, Double_t>,
-                        std::pair<Double_t, Double_t>);
-  Double_t CalcDistance(Double_t, Double_t,
-                        Double_t, Double_t);
-  
-  Double_t CalcPerpDistanceToLine(Float_t, Float_t, pair<Int_t, Int_t>);
-  Double_t CalcPerpDistanceToLine(Float_t, Float_t, Int_t, Int_t);
+Int_t ConvertXYPositionToInteger(std::pair<Int_t, Int_t> position, Int_t dimX, Int_t dimY);
+
+Int_t ConvertXYPositionToInteger(Int_t positionX, Int_t positionY, Int_t dimX, Int_t dimY);
+
+Double_t CalcDistance(std::pair<Double_t, Double_t>,
+                std::pair<Double_t, Double_t>);
+Double_t CalcDistance(Double_t, Double_t,
+                Double_t, Double_t);
+
+Double_t CalcPerpDistanceToLine(Float_t, Float_t, pair<Int_t, Int_t>);
+Double_t CalcPerpDistanceToLine(Float_t, Float_t, Int_t, Int_t);
 
 
-  /* Drawing */
-  void DrawLine(MediPixAlgo *, Float_t, Float_t, Float_t, Float_t, Int_t, Int_t, EColor);
+/* Drawing */
+void FillValuesForDisplay(Highlighter * hl, blob bl);
+void DrawLine(MediPixAlgo *, Float_t, Float_t, Float_t, Float_t, Int_t, Int_t, EColor);
 
-  /* Word handling */
-  std::string DumpWordInBinary(UInt_t word);
+/* Word handling */
+std::string DumpWordInBinary(UInt_t word);
 
-  // Blobs handling
-  vector<pair<Float_t, Float_t> > FindLimitPixels(blob);
+// Blobs handling
+vector<pair<Float_t, Float_t> > FindLimitPixels(blob);
 
-  /*
-   * Check if a blob appears within a certain distance from the edges
-   * Parameters:
-   *  - blob
-   *  - vector containing farthest pixels in the blob (see FindLimitPixels(blob))
-   *  - Width of sensor in units of pixels
-   *  - Height of sensor in units of pixels
-   */
-  Bool_t BlobAtADistanceFromEdges(blob, vector<pair<Float_t, Float_t> >, Int_t, Int_t, Int_t, Int_t);
+/*
+ * Check if a blob appears within a certain distance from the edges
+ * Parameters:
+ *  - blob
+ *  - vector containing farthest pixels in the blob (see FindLimitPixels(blob))
+ *  - Width of sensor in units of pixels
+ *  - Height of sensor in units of pixels
+ */
+Bool_t BlobAtADistanceFromEdges(blob, vector<pair<Float_t, Float_t> >, Int_t, Int_t, Int_t, Int_t);
 
-  /**
-   * Clustering stuff
-   */
-  int DigitalBlobReclustering(int, AllBlobsContainer *, blob, int, int, std::map<int,int>);
-  bool StartBlob(int, int, int, AllBlobsContainer *, blob *, int, int, std::map<int,int>);
-  bool FollowBlob(blob *, steerSpiral *, int, int, int, std::map<int,int>);
-  void RewindSpiral(steerSpiral *);
-  Int_t GetLengthOfSpiral(Int_t);
-  pair<Int_t, Int_t> GetNextPosition(pair<Int_t, Int_t>, steerSpiral *);
-  bool IsUnsafePosition(pair<Int_t, Int_t>, int, int);
-  void ClearOneBlobData(blob *);
+/**
+ * Clustering stuff
+ */
+int DigitalBlobReclustering(int, AllBlobsContainer *, blob, int, int, std::map<int,int>);
+bool StartBlob(int, int, int, AllBlobsContainer *, blob *, int, int, std::map<int,int>);
+bool FollowBlob(blob *, steerSpiral *, int, int, int, std::map<int,int>);
+void RewindSpiral(steerSpiral *);
+Int_t GetLengthOfSpiral(Int_t);
+pair<Int_t, Int_t> GetNextPosition(pair<Int_t, Int_t>, steerSpiral *);
+bool IsUnsafePosition(pair<Int_t, Int_t>, int, int);
+void ClearOneBlobData(blob *);
 
-  /* Dump contents of a blob */
-  void DumpBlobContents(blob);
+/* Dump contents of a blob */
+void DumpBlobContents(blob);
 
-  int FindDiscontinuity(blob);
+int FindDiscontinuity(blob);
+
+
+///////////////////////////////////////////////
+// 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;
+
+};
+
 }