TPC/groot/ABlob.C

0001 #include "ABlob.h"
0002 #include "AZigzag.h"
0003 #include <iostream>
0004
0005 #include "groot.h"
0006
0007 #include "TH1D.h"
0008 #include "TF1.h"
0009
0010 #include "TMarker.h"
0011
0012 using namespace std;
0013
0014 bool ABlob::RecalibrateCharge = true;
0015 bool ABlob::GaussPosition = false;
0016 TH1* ABlob::BlobPos   = 0;
0017 TH1* ABlob::BlobSigma = 0;
0018 TF1* ABlob::BlobFit   = 0;
0019
0020 ABlob::ABlob(std::vector<AZigzag*> MANYZIGZAGS)
0021 {
0022   manyZigs = MANYZIGZAGS;
0023
0024   correctedq = 0;
0025   correctedcentroid = 0;
0026
0027   if (RecalibrateCharge) FixTheCharges();
0028
0029   Precalc_PHI = GetPHI();
0030   Precalc_R   = manyZigs[0]->RCenter();
0031 }
0032
0033
0034 double ABlob::CentroidX() { return Precalc_R * cos(Precalc_PHI); }
0035
0036 double ABlob::CentroidY() { return Precalc_R * sin(Precalc_PHI); }
0037
0038 double ABlob::CentroidR() { return Precalc_R; }
0039
0040 double ABlob::CentroidP() { return Precalc_PHI; }
0041
0042 double ABlob::maxT()      { return manyZigs[0]->T(); } //time of zigzag with the max charge

0043
0044
0045 double ABlob::Q()
0046 {
0047   double den = 0;
0048   for(int i=0; i<manyZigs.size(); i++)
0049     {
0050       den += manyZigs[i]->Q();
0051     }
0052
0053   return den;
0054 }
0055
0056 double ABlob::CentroidZ()
0057 {
0058   if (Q() == 0)
0059     {
0060       cout << "ABlob:: Error: Blob has zero total charge...no Z centroid" << endl;
0061       return -9999;
0062     }
0063
0064   double num = 0;
0065   for(int i=0; i<manyZigs.size(); i++)
0066     {
0067       num += manyZigs[i]->Q() * manyZigs[i]->ZCenter();
0068     }
0069
0070   return num/Q();
0071 }
0072
0073
0074 double ABlob::GetPHI() //Angle Phi

0075 {
0076   if (Q() == 0)
0077     {
0078       cout << "ABlob:: Error: Blob has zero total charge...no Phi centroid" << endl;
0079       return -9999;
0080     }
0081
0082   if (!GaussPosition)
0083     {
0084       double num = 0;
0085       for(int i=0; i<manyZigs.size(); i++)
0086     {
0087       num += manyZigs[i]->Q() * manyZigs[i]->PCenter();
0088     }
0089
0090       return num/Q();
0091     }
0092   else
0093     {
0094       groot* Tree = groot::instance();
0095       if (!BlobPos)
0096     {
0097       double Minphi =  9999;
0098       double Maxphi = -9999;
0099       for (int i=0; i<Nphi; i++)
0100         {
0101           if (Tree->ZigzagMap2[0][i])
0102         {
0103           double phi = Tree->ZigzagMap2[0][i]->PCenter();
0104           if (phi < Minphi) Minphi = phi;
0105           if (phi > Maxphi) Maxphi = phi;
0106         }
0107         }
0108       double dPhi = abs( Tree->ZigzagMap2[0][90]->PCenter() - Tree->ZigzagMap2[0][91]->PCenter() );
0109       int NBINS = int( (Maxphi - Minphi)/dPhi + 1.5 );
0110       double left  = Minphi - 0.5*dPhi;
0111       double right = Maxphi + 0.5*dPhi;
0112
0113       BlobPos   = new TH1D("BlobPos", "BlobPos", NBINS, left, right);
0114       BlobSigma = new TH1D("BlobSigma", "BlobSigma", 250, 0.0, 8.0E-3);
0115       BlobFit   = new TF1("BlobFit", "[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))", left, right);
0116       Tree->theHistograms.push_back(BlobSigma);
0117     }
0118
0119       BlobPos->Reset();
0120
0121       for(int i=0; i<manyZigs.size(); i++)
0122     {
0123       double phi = manyZigs[i]->PCenter();
0124       double Q   = manyZigs[i]->Q();
0125       BlobPos->Fill(phi,Q);
0126     }
0127
0128       for (int i=1; i<=BlobPos->GetNbinsX(); i++)
0129     {
0130       BlobPos->SetBinError( i, 9.0 );
0131     }
0132
0133       double dPhi  = abs( Tree->ZigzagMap2[0][90]->PCenter() - Tree->ZigzagMap2[0][91]->PCenter() );
0134       double max   = BlobPos->GetMaximum();
0135       double mean  = manyZigs[0]->PCenter();
0136       double sigma = 0.003;
0137       if (sigma < dPhi/3.0) sigma = dPhi/3.0;
0138       if (sigma > 5.0*dPhi) sigma = 5.0*dPhi;
0139
0140       BlobFit->SetParameter(0,max);
0141       BlobFit->SetParameter(1,mean);
0142       BlobFit->SetParameter(2,sigma);
0143
0144       BlobFit->SetParLimits(0, max/2.0, 3.0*max);
0145       BlobFit->SetParLimits(1, mean-1.0*dPhi, mean+1.0*dPhi);
0146       BlobFit->SetParLimits(2, 0.0007, 0.0055);
0147
0148       BlobPos->Fit(BlobFit,"Q0");
0149
0150       BlobSigma->Fill(BlobFit->GetParameter(2));
0151
0152       return BlobFit->GetParameter(1);;
0153     }
0154 }
0155
0156 void ABlob::Draw()
0157 {
0158   // color zigzag according to charge fraction of total deposited charge

0159   int shape = 8;
0160   double X = CentroidX();
0161   double Y = CentroidY();
0162
0163   TMarker* ablob = new TMarker(X,Y,shape);
0164
0165   ablob->SetMarkerSize(1);
0166   ablob->SetMarkerColor(6);
0167
0168   ablob->Draw();
0169 }
0170
0171 void ABlob::Report()
0172 {
0173   //cout << "~~BLOB REPORT~~" << endl;

0174   //cout << "Q:" << Q();

0175   //cout << " \tX:" << CentroidX();

0176   //cout << " \tY:" << CentroidY();

0177   //cout << " \tZ:" << CentroidZ();

0178   cout << " \tR:" << CentroidR();
0179   cout << " \tP:" << CentroidP();
0180   //cout << " \tT:" << maxT();

0181   //cout << " \tsize:" << manyZigs.size();

0182   cout << endl;
0183 }
0184
0185 void ABlob::FixTheCharges()
0186 {
0187   //  This routine is an attempt to improve the values of charge

0188   //  from the "edge" zigzags of the blob.  It is premised upon the

0189   //  following assumptions:

0190   //

0191   //  1) When the charge hit is big, the Zigzag does a fine job

0192   //     of determining both Q and T.

0193   //  2) When the charge hit is small, the fit can be drawn away

0194   //     from the correct Q by following a fluctuation at the wrong

0195   //     time.

0196   //  3) Since every blob should have the same time for all zigzags,

0197   //     we can pull a trick for the outlying pads as follows:

0198   //     a:  Determine a reference time by accessing the time from the

0199   //         pad with the most charge.

0200   //     b:  Perform a NEW FIT of the charge on every pad by restricting

0201   //         the time fit range to be "near" to that of the reference time.

0202   //

0203   //  In this way outlier pads will always use the "in time" charge instead

0204   //  of the "highest found" charge and hopefully improve the overall

0205   //  blob position resolution.

0206   //

0207   //                                       TKH 10/10/2018

0208   ///////////////////////////////////////////////////////////////////////

0209
0210   double Reftime = maxT();
0211   double Mintime = Reftime - 2.0;
0212   double Maxtime = Reftime + 2.0;
0213
0214   for(int i=0; i<manyZigs.size(); i++)
0215     {
0216       manyZigs[i]->DetermineQ(Mintime,Maxtime);
0217     }
0218
0219 }