TPC/groot/FillBlobHist.C

0001 ////////////////////////////////////
0002 //        FillBlobHist.C          //
0003 ////////////////////////////////////
0004
0005 // This makes a BH (Blob Bistogram) that can show the raw hits.
0006
0007 #include "FillBlobHist.h"
0008 #include "TH1D.h"
0009 #include "TH2D.h"
0010 #include "ABlob.h"
0011 #include "AZigzag.h"
0012
0013 #include "groot.h"
0014
0015 #include <iostream>
0016 #include <cmath>
0017 #include <string>
0018
0019 TH1D* BH_Q = 0;       // How many blobs of a certain charge
0020 TH1D* BH_Size = 0;    // How many blobs made of a certain # of zig-zag pads
0021 TH2D* BH_QvsSize = 0; // Charge of a blob vs. how many pads it's made of
0022
0023 TH1D* BH_PhiTic = 0;  // Just the phi location of each pad
0024 TH2D* BH_PhivsR = 0;  // Phi vs. Radius
0025
0026 TH1D* BH_Tmax = 0; // This is the time of the zig with the most charge
0027 TH1D* BH_Tdiff = 0;
0028
0029 TH1D* BH_5Residual = 0; // Minor cheating based on arc-length for a quick preliminary resolution
0030 TH2D* BH_5Residual_vsPhi = 0;
0031  TH1D* BH_5Residual_PC = 0;
0032  TH2D* BH_5Residual_vsPhi_PC = 0;
0033   TH1D* BH_5Efficiency = 0; // Checks if a blob even exists on the 5th radii if 4 & 6 are hit
0034
0035
0036 // Histograming ALL Radii through an array
0037 char name[100];
0038 TH1D* BH_Residual_[16]; // Initialize rows 0 and 15 (in case they're useful later)
0039 TH2D* BH_Residual_vsPhi_[16];
0040  TH1D* BH_Residual_PC_[16]; // PC = PhiCut (efficiency/track attempt before Hough)
0041  TH2D* BH_Residual_vsPhi_PC_[16];
0042   TH1D* BH_Efficiency_[16];
0043
0044
0045 using namespace std;
0046
0047 void FillBlobHist()
0048 {
0049   groot* Tree = groot::instance();
0050
0051   // phi=90 and phi=91 are instrumented (connected to hybrids)
0052   double phi90 = Tree->ZigzagMap2[0][90]->PCenter();
0053   double phi91 = Tree->ZigzagMap2[0][91]->PCenter();
0054
0055   double phi0   = 91.0*phi90 - 90.0*phi91;
0056   double phi127 = phi0 + 127.0*(phi91-phi90);
0057
0058   double BIN_FACTOR = 10.0;
0059   double dPhi_BIN   = (phi91-phi90)/BIN_FACTOR;
0060
0061   double min = phi0 - dPhi_BIN/2.0;
0062   double max = phi127 + dPhi_BIN/2.0;
0063   int  N_BIN = 127*BIN_FACTOR;
0064
0065   if (!BH_Q)
0066     {
0067       BH_Q       = new TH1D("BH_Q", "BH_Q",              256, -0.5, 2047.5);
0068       BH_Size    = new TH1D("BH_Size", "BH_Size",         21, -0.5, 20.5);
0069       BH_QvsSize = new TH2D("BH_QvsSize", "BH_QvsSize",  256, -0.5, 2047.5, 21, -0.5, 20.5);
0070       BH_PhiTic  = new TH1D("BH_PhiTic", "BH_PhiTic",    2*N_BIN, min, max); // "normalized" tics
0071       BH_PhivsR  = new TH2D("BH_PhivsR", "BH_PhivsR",      N_BIN, min, max, 16, -0.5, 15.5);
0072       BH_Tmax    = new TH1D("BH_Tmax", "BH_Tmax",        300, -5.0, 30.0); // "normalized" tics
0073       BH_Tdiff   = new TH1D("BH_Tdiff", "BH_Tdiff",      20000, -30.0, 30.0); // "normalized" tic difference
0074
0075       // Resolution-study Histograms
0076       BH_5Residual          = new TH1D("BH_5Residual", "BH_5Residual", 200, -1000.0, 1000.0); //arc length in microns
0077       BH_5Residual_vsPhi    = new TH2D("BH_5Residual_vsPhi", "BH_5Residual_vsPhi", (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
0078       BH_5Residual_PC       = new TH1D("BH_5Residual_PC", "BH_5Residual_PC", 200, -1000.0, 1000.0); //arc length in microns
0079       BH_5Residual_vsPhi_PC = new TH2D("BH_5Residual_vsPhi_PC", "BH_5Residual_vsPhi_PC", (int)(0.2/dPhi_BIN), 1.53, 1.63, 200, -1000.0, 1000.0);
0080       BH_5Efficiency        = new TH1D("BH_5Efficiency", "BH_5Efficiency", 2, 0, 2);
0081
0082
0083       // Histograming ALL Radii through an array
0084       for (int i=1; i<15; i++) // Can't do radii 0 or 15 since we don't have radii -1 or 16 for comparison
0085     {
0086       sprintf(name, "BH_Residual_%d", i);
0087                      BH_Residual_[i]          = new TH1D(name, name, 200,-1000.0, 1000.0);
0088       sprintf(name, "BH_Residual_vsPhi_%d", i);
0089                      BH_Residual_vsPhi_[i]    = new TH2D(name, name, (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
0090       sprintf(name, "BH_Residual_PC_%d", i);
0091                      BH_Residual_PC_[i]       = new TH1D(name, name, 200,-1000.0, 1000.0);
0092       sprintf(name, "BH_Residual_vsPhi_PC_%d", i);
0093                      BH_Residual_vsPhi_PC_[i] = new TH2D(name, name, (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
0094       sprintf(name, "BH_Efficiency_%d", i);
0095                      BH_Efficiency_[i]        = new TH1D(name, name, 2, 0, 2);
0096     }
0097
0098
0099       // Used for Condor
0100       Tree->theHistograms.push_back(BH_Q);
0101       Tree->theHistograms.push_back(BH_Size);
0102       Tree->theHistograms.push_back(BH_QvsSize);
0103       Tree->theHistograms.push_back(BH_PhiTic );
0104       Tree->theHistograms.push_back(BH_PhivsR);
0105       Tree->theHistograms.push_back(BH_Tmax);
0106       Tree->theHistograms.push_back(BH_Tdiff);
0107       Tree->theHistograms.push_back(BH_5Residual);
0108       Tree->theHistograms.push_back(BH_5Residual_PC);
0109       Tree->theHistograms.push_back(BH_5Residual_vsPhi);
0110       Tree->theHistograms.push_back(BH_5Residual_vsPhi_PC);
0111       Tree->theHistograms.push_back(BH_5Efficiency);
0112       for (int i=1; i<15; i++)
0113     {
0114       Tree->theHistograms.push_back(BH_Residual_[i]);
0115       Tree->theHistograms.push_back(BH_Residual_vsPhi_[i]);
0116       Tree->theHistograms.push_back(BH_Residual_PC_[i]);
0117       Tree->theHistograms.push_back(BH_Residual_vsPhi_PC_[i]);
0118       Tree->theHistograms.push_back(BH_Efficiency_[i]);
0119     }
0120     } // End of creating histograms
0121
0122
0123   if ( BH_PhiTic->GetMaximum() != 1) //Gets the positions of the pads. Hist hight 0 to give better BH_->Scale(int) control
0124     {
0125       for (int i=0; i<Nphi; i++)
0126     {
0127       if (Tree->ZigzagMap2[0][i] ) { BH_PhiTic->Fill(Tree->ZigzagMap2[0][i]->PCenter()); }
0128     }
0129     }
0130
0131   // All radii are in mm. Only convert to microns in resolution since all x,y values are also in mm
0132   double RadiusValue[16];
0133   for (int k=0; k<16; k++) { RadiusValue[k] = Tree->theZigzags[k]->RCenter(); }
0134
0135
0136   for (int i=0; i<Nr; i++)
0137     {
0138       for (int j=0; j< Tree->theBlobs[i].size(); j++)
0139     {
0140       double Q   = Tree->theBlobs[i][j]->Q();
0141       double Phi = Tree->theBlobs[i][j]->CentroidP();
0142       int Radius = i;
0143
0144       int numZigs = Tree->theBlobs[i][j]->numZigs();
0145
0146       if ( Tree->theBlobs[i].size()==1 ) // Only 1 blob per radius
0147         {
0148           BH_Q->Fill(Q);
0149           BH_Size->Fill(numZigs);
0150           BH_QvsSize->Fill(Q, numZigs);
0151
0152           BH_PhivsR->Fill(Phi, Radius);
0153
0154           double TMAX = Tree->theBlobs[i][0]->maxT();
0155           BH_Tmax->Fill(TMAX);
0156           for (int k=0; k< Tree->theBlobs[i][0]->manyZigs.size(); k++)
0157         {
0158           BH_Tdiff->Fill( Tree->theBlobs[i][0]->manyZigs[k]->T() - TMAX );
0159         }
0160         }
0161     }
0162     } // End of priliminary histograms
0163
0164
0165
0166   // Original full Resolution studies (over all angles)
0167   if ( Tree->theBlobs[3].size()==1 && Tree->theBlobs[4].size()==1 && Tree->theBlobs[5].size()==1 )
0168     {
0169       // Calculation to handle polar coordinates correctly
0170       double X4 = Tree->theBlobs[3][0]->CentroidX(), Y4 = Tree->theBlobs[3][0]->CentroidY();
0171       double X6 = Tree->theBlobs[5][0]->CentroidX(), Y6 = Tree->theBlobs[5][0]->CentroidY();
0172
0173       double x5a, y5a, x5b, y5b; //line must intersect the circle in 2 points
0174       int found = line_circle(X4, Y4, X6, Y6, RadiusValue[4], x5a, y5a, x5b, y5b); //takes 9 arguments
0175
0176       if (found == 2) // Nonesense unless both are found
0177     {
0178       double X5p = x5a, Y5p = y5a;
0179       if (y5b > 0)
0180         {
0181           X5p = x5b; Y5p = y5b;
0182         }
0183
0184       double Phi5p = atan2(Y5p, X5p);                   // predicted/expected angle, should be around 1.58 [rad]
0185       double Phi5a = Tree->theBlobs[4][0]->CentroidP(); // actual angle
0186       double Residual5 = Phi5a - Phi5p;                 // difference in angles
0187
0188       BH_5Residual->Fill(RadiusValue[4]*Residual5*1000.0); //angle * radius * conversion = arc length [microns]
0189       BH_5Residual_vsPhi->Fill(Phi5p, RadiusValue[4]*Residual5*1000.0); //Fill (predicted/expected angle, arc length[microns]);
0190     } // End of filling actual histographs
0191     } // End of making sure only one blob per radii
0192
0193
0194   // New efficiency histograms with a Phi cut
0195   if ( Tree->theBlobs[3].size()==1 && Tree->theBlobs[5].size()==1 )
0196     {
0197       if ( (Tree->theBlobs[3][0]->CentroidP() >= 1.53) && (Tree->theBlobs[3][0]->CentroidP() <= 1.63) &&
0198        (Tree->theBlobs[5][0]->CentroidP() >= 1.53) && (Tree->theBlobs[5][0]->CentroidP() <= 1.63) )
0199     {
0200       if ( Tree->theBlobs[4].size()==1 )
0201         {
0202           BH_5Efficiency->Fill(1.5);
0203
0204           // Calculation to handle polar coordinates correctly
0205           double X4 = Tree->theBlobs[3][0]->CentroidX(), Y4 = Tree->theBlobs[3][0]->CentroidY();
0206           double X6 = Tree->theBlobs[5][0]->CentroidX(), Y6 = Tree->theBlobs[5][0]->CentroidY();
0207
0208           double x5a, y5a, x5b, y5b; //line must intersect the circle in 2 points
0209           int found = line_circle(X4, Y4, X6, Y6, RadiusValue[4], x5a, y5a, x5b, y5b); //takes 9 arguments
0210
0211           if (found == 2) // Nonesense unless both are found
0212         {
0213           double X5p = x5a, Y5p = y5a;
0214           if (y5b > 0)
0215             {
0216               X5p = x5b; Y5p = y5b;
0217             }
0218
0219           double Phi5p = atan2(Y5p, X5p);                   // predicted/expected angle, should be around 1.58 [rad]
0220           double Phi5a = Tree->theBlobs[4][0]->CentroidP(); // actual angle
0221           double Residual5 = Phi5a - Phi5p;                 // difference in angles
0222
0223           BH_5Residual_PC->Fill(RadiusValue[4]*Residual5*1000.0); //angle * radius * conversion = arc length [microns]
0224           BH_5Residual_vsPhi_PC->Fill(Phi5p, RadiusValue[4]*Residual5*1000.0); //Fill (predicted/expected angle, arc length[microns]);
0225         } // End of filling actual histographs
0226         } // End of one blob on 5
0227       else { BH_5Efficiency->Fill(0.5); }
0228     } // End of one blob on 4 & 6 within angle cut
0229     } // End of one blob per adjacent radii
0230
0231
0232
0233
0234
0235
0236
0237
0238 //---------------------------------------------------------------------------------------------------------------------------------------
0239 //------------------------------------------------|ALL RESOLUTIONS|----------------------------------------------------------------------
0240 //---------------------------------------------------------------------------------------------------------------------------------------
0241
0242   for (int i=1; i<15; i++)
0243     {
0244       // Original full Resolution studies (over all angles)
0245       if ( Tree->theBlobs[i-1].size()==1 && Tree->theBlobs[i].size()==1 && Tree->theBlobs[i+1].size()==1 )
0246     {
0247       // Calculation to handle polar coordinates correctly
0248       double Xpre  = Tree->theBlobs[i-1][0]->CentroidX(), Ypre  = Tree->theBlobs[i-1][0]->CentroidY();
0249       double Xpost = Tree->theBlobs[i+1][0]->CentroidX(), Ypost = Tree->theBlobs[i+1][0]->CentroidY();
0250
0251       double xia, yia, xib, yib; //line must intersect the circle in 2 points
0252       int found = line_circle(Xpre, Ypre, Xpost, Ypost, RadiusValue[i], xia, yia, xib, yib); //takes 9 arguments
0253
0254       if (found == 2) // Nonesense unless both are found
0255         {
0256           double Xip = xia, Yip = yia;
0257           if (yib > 0) { Xip = xib; Yip = yib; }
0258
0259           double Phi_ip = atan2(Yip, Xip);                   // predicted/expected angle, should be around 1.58 [rad]
0260           double Phi_ia = Tree->theBlobs[i][0]->CentroidP(); // actual angle
0261           double Residuali = Phi_ia - Phi_ip;                // difference in angles
0262
0263           BH_Residual_[i]->Fill(RadiusValue[i]*Residuali*1000.0); //angle * radius * conversion = arc length [microns]
0264           BH_Residual_vsPhi_[i]->Fill(Phi_ip, RadiusValue[i]*Residuali*1000.0); //Fill (predicted/expected angle, arc length[microns]);
0265         } // End of filling actual histographs
0266     } // End of making sure only one blob per radii
0267
0268
0269       // New efficiency histograms with a Phi cut
0270       if ( Tree->theBlobs[i-1].size()==1 && Tree->theBlobs[i+1].size()==1 )
0271     {
0272       double PC_LOW  = 1.53;
0273       double PC_HIGH = 1.63;
0274       if ( (Tree->theBlobs[i-1][0]->CentroidP() >= PC_LOW) && (Tree->theBlobs[i-1][0]->CentroidP() <= PC_HIGH) &&
0275            (Tree->theBlobs[i+1][0]->CentroidP() >= PC_LOW) && (Tree->theBlobs[i+1][0]->CentroidP() <= PC_HIGH) )
0276         {
0277           if ( Tree->theBlobs[i].size()==1 )
0278         {
0279           BH_Efficiency_[i]->Fill(1.5);
0280
0281           // Calculation to handle polar coordinates correctly
0282           double Xpre  = Tree->theBlobs[i-1][0]->CentroidX(), Ypre  = Tree->theBlobs[i-1][0]->CentroidY();
0283           double Xpost = Tree->theBlobs[i+1][0]->CentroidX(), Ypost = Tree->theBlobs[i+1][0]->CentroidY();
0284
0285           double xia, yia, xib, yib; //line must intersect the circle in 2 points
0286           int found = line_circle(Xpre, Ypre, Xpost, Ypost, RadiusValue[i], xia, yia, xib, yib); //takes 9 arguments
0287
0288           if (found == 2) // Nonesense unless both are found
0289             {
0290               double Xip = xia, Yip = yia;
0291               if (yib > 0) { Xip = xib; Yip = yib; }
0292
0293               double Phi_ip = atan2(Yip, Xip);                   // predicted/expected angle, should be around 1.58 [rad]
0294               double Phi_ia = Tree->theBlobs[i][0]->CentroidP(); // actual angle
0295               double Residuali = Phi_ia - Phi_ip;                // difference in angles
0296
0297               BH_Residual_PC_[i]->Fill(RadiusValue[i]*Residuali*1000.0); //angle * radius * conversion = arc length [microns]
0298               BH_Residual_vsPhi_PC_[i]->Fill(Phi_ip, RadiusValue[i]*Residuali*1000.0); //Fill (predicted/expected angle, arc length[microns]);
0299             } // End of filling actual histographs
0300         } // End of one blob on 5
0301           else { BH_Efficiency_[i]->Fill(0.5); }
0302         } // End of one blob on 4 & 6 within angle cut
0303     } // End of one blob per adjacent radii
0304
0305     } // end of ARRAY of histograms
0306
0307 //---------------------------------------------------------------------------------------------------------------------------------------
0308 //------------------------------------------------|ALL RESOLUTIONS END|------------------------------------------------------------------
0309 //---------------------------------------------------------------------------------------------------------------------------------------
0310
0311 }
0312
0313
0314
0315
0316
0317
0318 int line_circle(double m, double b, double r, double& Xi1, double& Yi1, double& Xi2, double& Yi2)
0319 {
0320   return line_circle( 0.0, b, 1.0, (m+b), r, Xi1, Yi1, Xi2, Yi2);
0321 }
0322
0323
0324 int line_circle(double X1, double Y1, double X2, double Y2, double r,  double& Xi1, double& Yi1, double& Xi2, double& Yi2) //takes 9 arguments
0325 {
0326   double dx = X2 - X1;
0327   double dy = Y2 - Y1;
0328   double dr = sqrt(dx*dx + dy*dy);
0329   double D  = X1*Y2 - X2*Y1;
0330
0331   double DELTA = r*r*dr*dr - D*D;
0332   if (DELTA < 0) return 0;
0333
0334   if (DELTA == 0)
0335     {
0336       double sign = (dy > 0) - (dy < 0); // Implements the sign() function
0337       Xi1 =  D*dy/(dr*dr);
0338       Yi1 = -D*dx/(dr*dr);
0339       return 1;
0340     }
0341
0342   double sign = (dy > 0) - (dy < 0); // Implements the sign() function
0343   Xi1 = ( D*dy + sign*dx*sqrt(DELTA))/(dr*dr);
0344   Yi1 = (-D*dx + abs(dy)*sqrt(DELTA))/(dr*dr);
0345   Xi2 = ( D*dy - sign*dx*sqrt(DELTA))/(dr*dr);
0346   Yi2 = (-D*dx - abs(dy)*sqrt(DELTA))/(dr*dr);
0347   return 2;
0348 }