TPC/groot/ATrace.C

0001 #include "ATrace.h"  //always include class header
0002 #include "TF1.h"
0003 #include "TH1D.h"
0004 #include "TH2D.h"
0005 #include <iostream>
0006
0007 using namespace std;
0008
0009 int ATrace::numtraces=0;
0010
0011 bool ATrace::FastFits = true;
0012
0013 double ATrace::SingleEmean = 0;
0014 double ATrace::SingleEsigma = 0;
0015 double ATrace::DoubleEmean = 0;
0016 double ATrace::DoubleEsigma = 0;
0017
0018 TH2D* heightCompare=0;
0019 TH2D* timeCompare=0;
0020
0021 TH1D* heightDiff=0;
0022 TH1D* timeDiff=0;
0023
0024
0025 ATrace::ATrace()
0026 {
0027   voltage.clear();
0028   height=0;
0029   time=0;
0030   W=0;
0031   dW=0;
0032
0033   sprintf (fcn_name, "fermi%d", numtraces);
0034   fermi = new TF1(fcn_name, "[0]/(exp(-(x-[1])/[2])+1)+[3]", 0, 1024);
0035
0036   sprintf (hst_name, "trace%d", numtraces);
0037   trace = new TH1D(hst_name, hst_name, 1024, -0.5, 1023.5);
0038
0039   numtraces++;
0040 }
0041
0042 ATrace::~ATrace()
0043 {
0044   fermi->Delete();
0045   trace->Delete();
0046 }
0047
0048 double ATrace::FindMaximum(int n)
0049 {
0050
0051   double MaxSum = 9999;
0052   int MaxMid = -1;
0053   for (int mid=1+n; mid<voltage.size()-n; mid++)
0054     {
0055       // We know roughly where true electron pulses should be in time - only look in that range
0056       if( mid < 420 || mid > 460 ) continue;
0057       double Sum=0;
0058       for (int bin=mid-n; bin<mid+n; bin++)
0059     {
0060       Sum += voltage[bin];
0061     }
0062       if (Sum < MaxSum)
0063     {
0064       MaxSum = Sum;
0065       MaxMid = mid;
0066     }
0067     }
0068
0069   return MaxSum/(2.0*n+1.0);
0070 }
0071
0072 int ATrace::FindMaximumMiddle(int n)
0073 {
0074
0075   double MaxSum = 9999;
0076   int MaxMid = -1;
0077   for (int mid=1+n; mid<voltage.size()-n; mid++)
0078     {
0079       // We know roughly where true electron pulses should be in time -- only look in that range
0080       if( mid < 420 || mid > 460 ) continue;
0081       double Sum=0;
0082       for (int bin=mid-n; bin<mid+n; bin++)
0083     {
0084       Sum += voltage[bin];
0085     }
0086       if (Sum < MaxSum)
0087     {
0088       MaxSum = Sum;
0089       MaxMid = mid;
0090     }
0091     }
0092
0093   return MaxMid;
0094 }
0095
0096
0097 void ATrace::FitLeftEdge()
0098 {
0099   MakeTrace(+1);
0100   int bin = trace->GetMaximumBin();
0101   height  = trace->GetBinContent(bin);
0102   time    = 0;
0103
0104   // 1) get time and height from histogram
0105   if (FastFits)
0106     {
0107       bool found = false;
0108       for(int i=bin;i>0;--i)
0109     {
0110       if(found) break;
0111       if(trace->GetBinContent(i) < height/2.0)
0112         {
0113           found = true;
0114           time = trace->GetBinCenter(i);
0115         }
0116     }
0117     }
0118   // 2) get time and height from fit
0119   else
0120     {
0121       if (!heightCompare)
0122     {
0123       heightCompare = new TH2D("heightCompare", "heightCompare", 1000, -0.5, 999.5, 1000, -0.5, 999.5);
0124       timeCompare = new TH2D("timeCompare", "timeCompare", 1000, -0.5, 999.5, 1000, -0.5, 999.5);
0125       heightDiff = new TH1D("heightDiff", "heightDiff", 1000, -499.5, 500.5);
0126       timeDiff = new TH1D("timeDiff", "timeDiff", 1000, -499.5, 500.5);
0127     }
0128       double X = trace->GetBinCenter(bin);
0129
0130       fermi->SetParameter(0,height);
0131       fermi->SetParameter(1,X-2);
0132       fermi->SetParameter(2,5);
0133       fermi->SetParameter(3,2);
0134       fermi->SetRange(5,X+2);
0135       trace->Fit(fermi,"WWR","",5,X+2);
0136
0137       double ht = trace->GetFunction(fcn_name)->GetParameter(0);
0138       double tm = trace->GetFunction(fcn_name)->GetParameter(1);
0139
0140       heightCompare->Fill(height, ht);
0141       timeCompare->Fill(time, tm);
0142       heightDiff->Fill(height-ht);
0143       timeDiff->Fill(time-tm);
0144
0145       height = trace->GetFunction(fcn_name)->GetParameter(0);
0146       time   = trace->GetFunction(fcn_name)->GetParameter(1);
0147     }
0148 }
0149
0150
0151 int ATrace::NAboveThreshold( double thrup, double thrdown )
0152 {
0153
0154   int nAbove = 0;
0155
0156   bool belowThreshold = true;
0157
0158   for(int i=0; i<voltage.size(); i++)
0159     {
0160       if ( belowThreshold && -voltage[i] >= thrup )
0161     {
0162       nAbove++;
0163       belowThreshold = false;
0164     }
0165
0166       else if ( !belowThreshold && -voltage[i] < thrdown )
0167     {
0168       belowThreshold = true;
0169     }
0170     }
0171
0172   return nAbove;
0173
0174 }
0175
0176 double ATrace::PulseWidth( double thrup, double thrdown )
0177 {
0178
0179   //  The results of this routine are ONLY valid
0180   //  if NAbove is one.
0181
0182   bool belowThreshold = true;
0183
0184   int left = 0;
0185   int right = 0;
0186
0187   for(int i=0; i<voltage.size(); i++)
0188     {
0189       if ( belowThreshold && voltage[i] >= thrup )
0190     {
0191       left = i;
0192       belowThreshold = false;
0193     }
0194
0195       else if ( !belowThreshold && voltage[i] < thrdown )
0196     {
0197       right = i;
0198       belowThreshold = true;
0199     }
0200     }
0201
0202   return right-left;
0203
0204 }
0205
0206 void ATrace::MakeTrace(int sign)
0207 {
0208   for(int i=0; i!=voltage.size(); ++i)
0209     {
0210       double vol = sign*voltage[i];
0211       trace->SetBinContent(i+1,vol);
0212     }
0213 }