FitNdNonlinQuad.h File Reference

#include "Validity/VldContext.h"
#include "Plex/PlexStripEndId.h"
#include "TClonesArray.h"

Go to the source code of this file.

Classes

class  NDQStripPulsePoint
class  NDQStripfit

Functions

void FitNdNonlinQuad (VldContext inputContext, const char *outputDirectory="ndfits/", bool writeToDb=false, bool useHighGainPins=false, bool drawChannelPlots=false)

Function Documentation

void FitNdNonlinQuad ( VldContext  inputContext,
const char *  outputDirectory = "ndfits/",
bool  writeToDb = false,
bool  useHighGainPins = false,
bool  drawChannelPlots = false 
)

Definition at line 62 of file FitNdNonlinQuad.cxx.

References MuELoss::a, NDQStripfit::a, NDQStripfit::adcmax, PlexPinDiodeId::AsString(), PlexStripEndId::AsString(), PlexPixelSpotId::AsString(), RawChannelId::AsString(), PlexLedId::AsString(), NDQStripfit::b, NDQStripfit::bad, NDQStripfit::beta, PlexStripEndId::BuildPlnStripEndKey(), NDQStripfit::c, NDQStripPulsePoint::caly, NDQStripfit::chisq, DbiWriter< T >::Close(), NDQStripfit::d, NDQStripfit::da, NDQStripfit::data, NDQStripfit::db, db_version, NDQStripfit::dbeta, NDQStripfit::dc, NDQStripPulsePoint::dcaly, NDQStripfit::dd, NDQStripfit::de, NDQStripfit::dof, NDQStripPulsePoint::dx, NDQStripPulsePoint::dy, NDQStripPulsePoint::dz, MuELoss::e, NDQStripfit::e, NDQStripfit::end, NDQStripPulsePoint::f, fcn(), Form(), NDQStripPulsePoint::fres, NDQStripPulsePoint::fxres, NDQStripPulsePoint::gain, PlexHandle::GetAllStripEnds(), PlexPinDiodeId::GetEncoded(), PlexStripEndId::GetEnd(), PulserGainPin::GetError(), PulserGain::GetError(), PlexPinDiodeId::GetGain(), PlexPinDiodeId::GetInBox(), PlexMuxBoxId::GetInRack(), PlexHandle::GetLedId(), PlexLedId::GetLedInBox(), PulserGain::GetMean(), PulserGainPin::GetMean(), PulserGain::GetNumEntries(), PulserGainPin::GetNumEntries(), PulserGainPin::GetNumPoints(), DbiResultPtr< T >::GetNumRows(), PulserGainPin::GetNumTriggers(), PulserGain::GetNumTriggers(), PlexHandle::GetPinDiodeIds(), PlexHandle::GetPixelSpotId(), PlexPlaneId::GetPlane(), PlexLedId::GetPulserBox(), PlexHandle::GetRawChannelId(), DbiResultPtr< T >::GetRow(), DbiResultPtr< T >::GetRowByIndex(), VldTimeStamp::GetSec(), PlexStripEndId::GetStrip(), DbiResultPtr< T >::GetValidityRec(), gSystem(), NDQStripPulsePoint::i, PlexPinDiodeId::IsLowGain(), SimFlag::kData, SimFlag::kMC, NDQStripfit::led, NDQStripfit::maxres, NDQStripfit::maxres2500, NDQStripfit::maxresf, NDQStripfit::maxresf2500, NDQStripfit::meanres2500, n, NDQStripfit::n, NDQStripfit::ok, NDQStripfit::pb, NDQStripPulsePoint::pe, NDQStripfit::pin, NDQStripfit::plane, NDQStripPulsePoint::res, NDQStripfit::Reset(), NDQStripfit::seid, CalPulserFits::SetAdcMax(), CalPulserFits::SetChisq(), CalPulserFits::SetFitType(), CalPulserFits::SetMaxRes(), CalPulserFits::SetMeanRes(), CalPulserFits::SetNPFit(), CalPulserFits::SetPar1(), CalPulserFits::SetPar2(), CalPulserFits::SetSlope(), CalPulserFits::SetSlopeErr(), CalPulserFits::SetXOffset(), NDQStripfit::strip, NDQStripPulsePoint::x, NDQStripPulsePoint::xres, NDQStripPulsePoint::y, NDQStripPulsePoint::z, and NDQStripPulsePoint::zerocorr.

00069 {
00070   PlexHandle plex(inputContext);
00071   DbiResultPtr<PulserGain>    fPmtTable(inputContext,1);//Pulser::kNearEnd);
00072   DbiResultPtr<PulserGainPin> fPinTable(inputContext,1);    
00073 
00074   if(fPmtTable.GetNumRows()==0) {
00075     cout << "No PMT data available for cx " << inputContext.AsString() << endl;
00076     return;
00077   }
00078   if(fPinTable.GetNumRows()==0) {
00079     cout << "No PIN data available for cx " << inputContext.AsString() << endl;
00080     return;
00081   }
00082 
00083 
00084   // Create output directory if it doesn't exist.
00085   gSystem->MakeDirectory(outputDirectory);
00086   
00087 
00088   DbiWriter<CalPulserFits>* writer = 0;
00089   if(writeToDb) {
00090     // Do this for the first prototype table.
00091     //VldTimeStamp start(1970,1,1,0,0,0,0);
00092     //VldTimeStamp stop(2038,1,1,0,0,0,0); // 1 year later.
00093     //VldTimeStamp creation = start + VldTimeStamp((time_t)db_version,0);
00094 
00095     // Otherwise, use the validity time from the input data.
00096     VldTimeStamp start = fPmtTable.GetValidityRec(fPmtTable.GetRow(0))->GetVldRange().GetTimeStart();
00097     VldTimeStamp stop(start.GetSec()+(3*365*24*3600),0); // 3 years later.
00098     VldTimeStamp creation(start.GetSec() + (int)db_version,0);
00099 
00100 
00101     VldRange range((int)inputContext.GetDetector(),
00102                    SimFlag::kData + SimFlag::kMC,
00103                    start,
00104                    stop,
00105                    "");
00106     writer = new DbiWriter<CalPulserFits>(range,
00107                                           -1, // aggragate
00108                                           44, // task
00109                                           creation,
00110                                           "offline", //db
00111                                           "ND Quadratic fit");                                    
00112     cout << "Starting DB write with range: " << endl;
00113     range.Print();
00114   }
00115 
00116   gStyle->SetPadColor(kWhite);
00117   
00118   TCanvas* c1;
00119   if((c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c1"))) {
00120   }
00121   else   c1 = new TCanvas("c1","c1",800,850);
00122   c1->cd();
00123 
00124   TFile* f = new TFile("ndlin_fits.root","RECREATE");
00125   f->cd();
00126 
00127   NDQStripfit* fitptr = new NDQStripfit(40);
00128   NDQStripfit& fit = *fitptr;
00129   TTree* tree = new TTree("fits","fits");
00130   tree->Branch("fitbranch","NDQStripfit",&fitptr);
00131 
00132   TH1F* hchi = new TH1F("hchi","ChiSq/dof",200,0,50);
00133   TH1F* hres = new TH1F("hres","Fractional Residual",500,-1,1);
00134   TH2F* hres2= new TH2F("hres2","Fractional Residual",
00135                         500,0,35000,
00136                         500,-1,1);
00137 
00138   TH2F* hnon2= new TH2F("hnon2","Approx. Nonlinearity",
00139                         500,0,35000,
00140                         500,-1,0.3);
00141 
00142   TH2F* htotres2= new TH2F("htotres2","TotalResidual",
00143                            500,0,35000,
00144                            500,-1000,1000);
00145 
00146 
00147   TH2F* hres2i= new TH2F("hres2i","TotalResidual",
00148                            40,0,40,
00149                            100,-1,1);
00150 
00151   std::map<PlexPinDiodeId,TH2*> mres2pin;
00152 
00153   TH1F* ha  = new TH1F("ha","Fit Parameter a (offset)",250,-500,500);
00154   ha->SetXTitle("ADC");
00155   ha->SetYTitle("stripends");
00156   TH1F* hb  = new TH1F("hb","Fit Parameter b (linear)",180,-10,80);
00157   hb->SetXTitle("ADC/PIN");
00158   hb->SetYTitle("stripends");
00159   TH1F* hc  = new TH1F("hc","Fit Parameter c (quad)",200,-0.01,0.01);
00160   hc->SetXTitle("ADC/PIN^{2}");
00161   hc->SetYTitle("stripends");
00162     
00163   TH1F* hd  = new TH1F("hd","Constant d: PIN offset",300,-300,300);
00164   hd->SetXTitle("PIN");
00165   hd->SetYTitle("stripends");
00166   TH1F* he  = new TH1F("he","Constant e: PIN/ADC conversion",240,-20,100);
00167   he->SetXTitle("PIN/ADC");
00168   he->SetYTitle("stripends");
00169   TH1F* hbeta  = new TH1F("hbeta","Nonlinearity Term #beta",200,-2e-5,1e-5);
00170   he->SetXTitle("ADC^{-1}");
00171   he->SetYTitle("stripends");
00172 
00173   
00174   TH1F* hbadres   = new TH1F("hbadres","Worst Residual",200,0,10000);
00175   TH1F* hbadresf  = new TH1F("hbadresf","Worst Fractional Residual",200,0,2);
00176   TH1F* hbadres2500   = new TH1F("hbadres2500","Worst Residual above 2500 ADC",200,0,10000);
00177   TH1F* hbadresf2500  = new TH1F("hbadresf2500","Worst Fractional Residual above 2500 ADC",200,0,2);
00178  
00179   ofstream ofit("fitdump.dat");
00180 
00181   //  for(int stp=10;stp<20;stp++) {
00182   //  PlexStripEndId seid(Detector::kNear, 10, stp, StripEnd::kWest);
00183 
00184   TGraphErrors* gr =0;
00185   TGraphErrors* grnon =0;
00186   TGraphErrors* grres =0;
00187   TGraphErrors* grgain =0;
00188   TGraphErrors* grcal  =0;
00189   TGraphErrors* grzoom =0;
00190   
00191   const std::vector<PlexStripEndId>& list = plex.GetAllStripEnds();
00192   for(UInt_t stp=0;
00193       stp<list.size();
00194       //stp<100;
00195       stp++) {
00196     PlexStripEndId seid = list[stp];
00197   
00198     //for(int stp=0;stp<1;stp++) {
00199     //  PlexStripEndId seid(Detector::kNear,1,15,StripEnd::kWest);
00200 
00201     cout << "PlexStripEndId: " << seid.AsString() << endl;
00202     PlexLedId        led = plex.GetLedId(seid);
00203     PlexPinDiodeId pinid = plex.GetPinDiodeIds(led).first;
00204     if(useHighGainPins) pinid = plex.GetPinDiodeIds(led).second;
00205 
00206     // Select only the good HG pin diode.
00207     //if(led.GetPulserBox()!=1) continue;
00208     //if(led.GetLedInBox() <9 ) continue;
00209     //if(led.GetLedInBox() >11) continue;
00210     //if(led.GetLedInBox() != 4) continue;
00211     
00212 
00213     // Fill out id.
00214     fit.Reset(40);
00215     fit.seid = seid;
00216     fit.plane = seid.GetPlane();
00217     fit.strip = seid.GetStrip();
00218     fit.end = seid.GetEnd();
00219     fit.pb  = led.GetPulserBox();
00220     fit.led = led.GetLedInBox();
00221     fit.pin = pinid.GetEncoded();
00222     //fit.pmt = pmtid.GetEncoded();
00223 
00224     fit.ok = true;
00225     
00226     bool isBad = false;
00227 
00228     if(gr)     { delete gr;     gr     =0; }
00229     if(grnon)  { delete grnon;  grnon  =0; }
00230     if(grres)  { delete grres;  grres  =0; }
00231     if(grgain) { delete grgain; grgain =0; }
00232     if(grcal)  { delete grcal; grcal =0; }
00233     if(grzoom) { delete grzoom; grzoom =0; }
00234 
00235  
00236     const PulserGain*    pmt = fPmtTable.GetRowByIndex(seid.BuildPlnStripEndKey());
00237     const PulserGainPin* pin = fPinTable.GetRowByIndex(pinid.GetEncoded());
00238 
00239     if(!pmt){ cout << "No Pmt data. " << endl; fit.ok = false; tree->Fill(); continue; }
00240     if(!pin){ cout << "No pin data. " << endl; fit.ok = false; tree->Fill(); continue; }
00241     // check the pin, which may be bad:
00242     double pintot = 0;
00243     for(int i=0;i<pin->GetNumPoints();i++) pintot += (pin->GetMean())[i];
00244     if(pintot<=0) {
00245       cout << "Pin data is bad." << endl;
00246       fit.ok = false;
00247       tree->Fill();
00248       continue;
00249     }
00250 
00251     int n = pin->GetNumPoints();
00252     fit.n = n;
00253 
00254     // Stripend
00255     const Float_t* pmtmean = pmt->GetMean(); 
00256     const Float_t* pmterr  = pmt->GetError(); 
00257     const Float_t* pmtent  = pmt->GetNumEntries(); 
00258     const Float_t* pmttrg  = pmt->GetNumTriggers();
00259     // Pin
00260     const Float_t* pinmean = pin->GetMean(); 
00261     const Float_t* pinerr  = pin->GetError(); 
00262     const Float_t* pinent  = pin->GetNumEntries(); 
00263     const Float_t* pintrg  = pin->GetNumTriggers();
00264 
00265     Int_t dof = 0;
00266     Float_t lowerbound = 0;
00267     Float_t upperbound = 9e9;
00268 
00269     gr = new TGraphErrors(n);    
00270     for(int i=0;i<n;i++) {
00271       float x = pinmean[i];
00272       float dx = pinerr[i] + 4.; // Increase errors on x-coord to improve Chisq.
00273       float z = pmtmean[i];//*pmtent[i]/pmttrg[i]; // Zero correction.
00274       float dz = pmterr[i];//*pmtent[i]/pmttrg[i]; // Zero correction.
00275 
00276       if(z<300){
00277         // lowerbound = x+1.; // exclude point from fit
00278         // Correct for zeroes, and add a big error to these points:
00279         float zerocorr = pmtent[i]/pmttrg[i];
00280         z  = z*zerocorr;
00281         dz = dz*zerocorr + z*0.5*(zerocorr-1);
00282       }
00283       dof++;             // include point
00284 
00285       gr->SetPoint     (i,x, z );
00286       gr->SetPointError(i,dx,dz); // Increase the error on the PIN, which greatly improves the chisq.      
00287     }
00288     
00289     const char* fcnname = "pol2";
00290     
00291     //gr->Draw("a p e0");
00292     gr->Fit(fcnname,"eq","goff",lowerbound,upperbound);
00293     //gr->Draw("a p e0");
00294     TF1* fcn = gr->GetFunction(fcnname);  
00295     
00296 
00297     double chisq = fcn->GetChisquare();
00298     double a = fcn->GetParameter(0);
00299     double b = fcn->GetParameter(1);
00300     double c = fcn->GetParameter(2);
00301     double da = fcn->GetParError(0);
00302     double db = fcn->GetParError(1);
00303     double dc = fcn->GetParError(2);
00304 
00305     cout << chisq << "/" << dof << endl;
00306 
00307     // Badness based on initial fit:
00308     if( (1.0 - 4*a*c/(b*b)) < 0 ) {
00309       // Failed: no intercept!
00310       isBad=true;
00311     };
00312 
00313     // d is x-intercept of fit. (PIN count offset)
00314     double d = (b/(2.*c)) * (-1.0 + sqrt(1.0 - 4*a*c/(b*b)) );
00315     // e is conversion factor in units of PIN/ADC
00316     double e = 1.0/(b+2.*d*c);
00317     double de = (db/b)*e; // roughly correct.
00318     // so, z = y + y*y*beta
00319     double beta = e*e*c;
00320     double dbeta = (dc/c)*beta;
00321 
00322     fit.dof = dof;
00323     fit.a = a;
00324     fit.b = b;
00325     fit.c = c;
00326     fit.d = d;
00327     fit.e = e;
00328     fit.beta = beta;
00329     fit.da = da;
00330     fit.db = db;
00331     fit.dc = dc;
00332     fit.dd = 0;
00333     fit.de = de;
00334     fit.dbeta = dbeta;
00335 
00336     // Badness based on parameter results:
00337     if(beta > -0) isBad = true;
00338 
00339     hchi->Fill(chisq/(float)n);
00340     ha->Fill(a);
00341     hb->Fill(b);
00342     hc->Fill(c);
00343     hd->Fill(d);
00344     he->Fill(e);
00345     hbeta->Fill(beta);
00346 
00347 
00348     TH2* hres2pin = 0;
00349     std::map<PlexPinDiodeId,TH2*>::iterator hit = mres2pin.find(pinid);
00350     if(hit == mres2pin.end()) {
00351       hres2pin = 
00352         new TH2F(Form("hpin_PB%02d_LED%02d_%c",led.GetPulserBox(),led.GetLedInBox(),(pinid.IsLowGain()?'L':'H')),
00353                  Form("hpin %s (%c)",led.AsString(),(pinid.IsLowGain()?'L':'H')),
00354                  //100,0,35000,
00355                  40,0,40,
00356                  100,-1,1);
00357       mres2pin[pinid] = hres2pin;
00358     } else {             
00359       hres2pin = hit->second;
00360     }
00361 
00362 
00363     float maxresidual = 0;
00364     float maxfracresidual = 0;
00365     float maxresidual2500 = 0;
00366     float maxfracresidual2500 = 0;
00367 
00368     float meanres2500 = 0;
00369     float meanres = 0;   
00370     float nres2500 = 0;
00371     float adcmax = 0;
00372     
00373     if(drawChannelPlots) {
00374       grres = new TGraphErrors(n);
00375       grnon = new TGraphErrors(n);
00376       grgain = new TGraphErrors(n);
00377       grcal = new TGraphErrors(n);
00378       grzoom = new TGraphErrors(n);
00379     }
00380 
00381     for(int i=0;i<n;i++) {
00382       // x is pin counts
00383       double x = pinmean[i];
00384       double dx= pinerr[i] + 4.;
00385       
00386       // y is linearized ADC counts.
00387       double y = (x-d)/e;
00388       double dy= dx/e;
00389 
00390       double zerocorr = pmtent[i]/pmttrg[i]; // Zero correction.
00391 
00392       // z is nonlinear ADC counts:
00393       double z = pmtmean[i]; 
00394       double dz = pmterr[i];
00395       if(z<300){
00396         // Correct for zeroes, and add a big error to these points:
00397         z  = z*zerocorr;
00398         dz = dz*zerocorr + 0.5*z*(1-zerocorr);
00399       }
00400 
00401       // The fit at this point (to be compared with z, the data)
00402       double f = fcn->Eval(x);
00403 
00404       // Do a trial calibration:
00405       double s = 1+4.*z*beta;
00406       if(s<0) s=0;
00407       double  caly =  1.0/(2.*beta) * ( sqrt(s) - 1. );
00408 
00409       double dcaly =  0;
00410       if(s>0) dcaly = 2.0 /  sqrt(s) *  dz;
00411 
00412       double res = z-f;
00413       double fres = (z-f)/f;
00414       double gain = (dz*dz*pmtent[i])/z * 0.8;
00415       double pe   = z/gain;
00416       
00417       // Horizontal residual:
00418       double xres = 9999;
00419       double temp = b*b - 4.*c*(a-z);
00420       if(temp>0) {
00421         temp = sqrt(temp);
00422         double x1 = (-b + temp)/(2.*c);
00423         double x2 = (-b - temp)/(2.*c);
00424         if(fabs(x-x1) < fabs(x-x2))
00425           xres = x-x1;
00426         else
00427           xres = x-x2;
00428       }
00429       double fxres = xres/x;
00430       
00431       NDQStripPulsePoint& datum = *((NDQStripPulsePoint*) fit.data->operator[](i));
00432       datum.i  =i;
00433       datum.x  =x;
00434       datum.dx =dx;
00435       datum.y  =y;
00436       datum.dy =dy;
00437       datum.z  =z;
00438       datum.dz =dz;
00439       datum.f  =f;
00440       datum.caly  = caly;
00441       datum.dcaly =dcaly;
00442       datum.res = res;
00443       datum.fres= fres;
00444       datum.xres = xres;
00445       datum.fxres = fxres;
00446       datum.gain = gain;
00447       datum.pe   = pe;
00448       datum.zerocorr = zerocorr;
00449       
00450       hres->Fill(fres);
00451       hres2->Fill(y,fres);
00452       hnon2->Fill(y,(z-y)/y);
00453       htotres2->Fill(y,res);
00454 
00455       //hres2i->  Fill(i,fres);
00456       //hres2pin->Fill(i,fres);
00457       hres2i->  Fill(i,fxres);
00458       hres2pin->Fill(i,fxres);
00459 
00460 
00461       if( fabs(res) > maxresidual) maxresidual = fabs(res);
00462       if( fabs(fres) > maxfracresidual) maxfracresidual = fabs(fres);
00463       meanres += fabs(res);
00464       if(y>2500) {
00465         if( fabs(res)  > maxresidual2500)       maxresidual2500 = fabs(res);
00466         if( fabs(fres) > maxfracresidual2500) maxfracresidual2500 = fabs(fres);
00467         meanres2500 += fabs(res);
00468         nres2500 += 1;
00469       }
00470       if(y>adcmax) adcmax = y;
00471  
00472       if(drawChannelPlots) {
00473         grres->SetPoint     (i,y,  (z - f)/f );
00474         grres->SetPointError(i,dy, dz/f);
00475 
00476         grnon->SetPoint     (i,y,  (z-y)/y);
00477         grnon->SetPointError(i,dy, dz/y);      
00478         
00479         grgain->SetPoint(i,y, (dz*dz*pmtent[i])/pmtmean[i] * 0.8 );
00480 
00481         grcal->SetPoint     (i,y, caly/y);
00482         grcal->SetPointError(i,dy,dcaly/y);
00483 
00484         grzoom->SetPoint     (i,y, z);
00485         grzoom->SetPointError(i,fabs(1-pinent[i]/pintrg[i])*y,dz);
00486       }
00487       
00488       //cout << x << "\t" << y << "\t" << z << "\t" << caly
00489       //           << "\t" << y+y*y*beta << endl;    
00490     }
00491 
00492     
00493     hbadres->Fill(maxresidual);
00494     hbadresf->Fill(maxfracresidual);
00495     hbadres2500->Fill(maxresidual2500);
00496     hbadresf2500->Fill(maxfracresidual2500);
00497     if(nres2500>0) meanres2500 /= nres2500;
00498     if(n>0) meanres/=n;
00499     if(chisq/(float)(dof)>75) isBad = true;
00500 
00501     fit.chisq = chisq;
00502     fit.maxres = maxresidual;
00503     fit.maxresf= maxfracresidual;
00504     fit.maxres2500  = maxresidual2500;
00505     fit.maxresf2500 = maxfracresidual2500;
00506     fit.meanres2500 = meanres2500;
00507     fit.adcmax = adcmax;
00508     fit.bad = isBad;
00509 
00510     if(isBad) {      
00511       ofit << Form("%s/p%03ds%03d.gif",outputDirectory,seid.GetPlane(),seid.GetStrip()) << endl;
00512     /*
00513       ofit << Form("%3d %3d %5.1f %2d  ",
00514       seid.GetPlane(),
00515       seid.GetStrip(),
00516       fcn->GetChisquare(),
00517       n);
00518       ofit << Form("  %.1f %.1f  %.1f %.1f  %.4f %.4f",
00519                    fcn->GetParameter(0), fcn->GetParError(0),
00520                    fcn->GetParameter(1), fcn->GetParError(1),
00521                    fcn->GetParameter(2), fcn->GetParError(2) );
00522       ofit << Form("  %.1f %.2f",maxresidual, maxfracresidual);
00523       ofit << endl;
00524     */
00525     }
00526 
00527     if(drawChannelPlots) {
00528       c1->cd();
00529       gPad->Clear();
00530       TVirtualPad* pad=gPad;
00531       TVirtualPad* spad1 = new TPad("spad1","spad1",0,  0,0.8,1);
00532       TVirtualPad* spad2 = new TPad("spad2","spad2",0.75,0,1,  1);
00533       spad1->Draw();
00534       spad2->Draw();
00535       
00536       // Text.
00537       spad2->cd();
00538       TPaveText* txt = new TPaveText(0,0.1,0.95,0.9,"NDC T");
00539       txt->SetFillColor(kWhite);
00540       txt->SetTextAlign(12);
00541       txt->SetTextSize(0.05);
00542       txt->AddText(seid.AsString());
00543       txt->AddText(Form("Plane: %d  Strip: %d",seid.GetPlane(),seid.GetStrip()));
00544       txt->AddText(plex.GetPixelSpotId(seid).AsString("p"));
00545       txt->AddText(plex.GetRawChannelId(seid).AsString("e"));
00546       txt->AddText(led.AsString("e"));
00547       txt->AddText(Form("%s (%02d%02d%01d)",pinid.AsString("e"),pinid.GetInRack(),pinid.GetInBox(),pinid.GetGain()));
00548       txt->AddText(Form("#{Chi}^{2}/DoF = %.1f/%d",chisq,dof));
00549       txt->AddText(Form("a: %.1f +/- %.1f ADC",a,da));
00550       txt->AddText(Form("b: %.1f +/- %.1f ADC/PIN",b,db));
00551       txt->AddText(Form("c: %.1e +/- %.1e ADC/PIN^{2}",c,dc));
00552       txt->AddText(" ");
00553       txt->AddText(Form("d: %.2f PIN",d));
00554       txt->AddText(Form("e: %.3f PIN/ADC",e));
00555       txt->AddText(Form("#beta: %.2e ADC^{-1}",beta));
00556       txt->Draw();      
00557       
00558       spad1->Divide(1,4,0.001,0.001,kWhite);
00559       spad1->cd(1); 
00560       gr->Draw("a p e0");
00561       gr->GetHistogram()->SetTitle("Full Curve with fit");
00562       gr->GetHistogram()->SetXTitle("PIN counts");
00563       gr->GetHistogram()->SetYTitle("PMT counts");
00564       
00565       spad1->cd(2);
00566       grnon->Draw("a p e0");
00567       grnon->GetHistogram()->SetTitle("Nonlinearity");
00568       grnon->GetHistogram()->SetXTitle("Linear PMT counts");
00569       grnon->GetHistogram()->SetYTitle("(Nonlinear - Linear)/Linear");
00570       
00571       spad1->cd(3);
00572       grres->Draw("a p e0");
00573       grres->GetHistogram()->SetTitle("Fractional Fit Residual");
00574       grres->GetHistogram()->SetXTitle("Linear PMT counts");
00575       grres->GetHistogram()->SetYTitle("(Data - Fit)/Fit");
00576       
00577       spad1->cd(4);
00578       //grgain->SetMarkerStyle(20);
00579       //grgain->Draw("a p");
00580       //grgain->GetHistogram()->SetTitle("Apparent gain");
00581       //grgain->GetHistogram()->SetXTitle("Linear PMT counts");
00582       //grgain->GetHistogram()->SetYTitle("PEs");
00583       //grcal->Draw("a p e0");
00584       //grcal->GetHistogram()->SetTitle("Calibrated result");
00585       //grcal->GetHistogram()->SetXTitle("Approximate linear ADC");
00586       //grcal->GetHistogram()->SetYTitle("Calibrated ADC / Pin-esimated ADC");
00587       static TH2* hframezoom = 0;
00588       static TLine* lzoom = 0;
00589       if(hframezoom == 0){
00590         hframezoom = new TH2F("hframezoom","Zoom",20,0,6000,20,0,6000);
00591         hframezoom->SetXTitle("Linear PMT counts");
00592         hframezoom->SetYTitle("Nonlinear PMT counts");
00593         lzoom = new TLine(0,0,6000,6000);
00594       }
00595       hframezoom->Draw();
00596       lzoom->Draw();
00597       grzoom->Draw("e0");
00598       
00599     
00600       pad->Update();
00601       pad->cd();
00602       
00603       c1->Print(Form("%s/p%03ds%03d.gif",outputDirectory,seid.GetPlane(),seid.GetStrip()));
00604     }
00605     
00606     tree->Fill();
00607     
00608     
00609     if(writeToDb) {     
00610       // Adjust for bad fits:
00611       if(beta>0) {
00612         dbeta +=fabs(beta);
00613         beta = 0;
00614       }
00615       if( (1.0 - 4*a*c/(b*b)) < 0 ) {
00616         d = 0;
00617       };
00618 
00619       CalPulserFits row(-1, //aggno
00620                         seid.BuildPlnStripEndKey(),
00621                         0, // fit type
00622                         (TF1*)0, 0, 0, 0);
00623       row.SetFitType(44);
00624       row.SetNPFit(dof);
00625       row.SetSlope(e);
00626       row.SetSlopeErr(de);
00627       row.SetXOffset(d);
00628       row.SetPar1(beta);
00629       row.SetPar2(dbeta);
00630       row.SetChisq(chisq);
00631       if(nres2500>0) {
00632         row.SetMeanRes(meanres2500);
00633         row.SetMaxRes(maxresidual2500);
00634       } else {
00635         row.SetMeanRes(meanres);
00636         row.SetMaxRes(maxresidual);
00637       }
00638       row.SetAdcMax(adcmax);
00639       //row.Dump();
00640       (*writer) << row;
00641     }
00642     
00643   }
00644   
00645   if(writeToDb) {
00646     writer->Close();
00647     delete writer;
00648   }
00649 
00650   tree->Write();
00651   hchi->Write();
00652   hres->Write();
00653   htotres2->Write();
00654 
00655   hnon2->Write();
00656 
00657   c1->cd(); c1->Clear();
00658   hres2->Draw("colz");
00659   c1->Update();
00660   c1->Print(Form("%s/_res2.gif",outputDirectory));
00661   hres2->Write();
00662 
00663   hres2i->Draw("colz");
00664   c1->Print(Form("%s/_res2i.gif",outputDirectory));
00665   c1->Update();
00666   hres2i->Write();
00667 
00668   ha->Write();
00669   hb->Write();
00670   hc->Write();
00671   hd->Write();
00672   he->Write();
00673   hbeta->Write();
00674   hbadres->Write();
00675   hbadresf->Write();
00676   hbadres2500->Write();
00677   hbadresf2500->Write();
00678   for(std::map<PlexPinDiodeId,TH2*>::iterator hit = mres2pin.begin();
00679       hit!=mres2pin.end();
00680       hit++) {
00681     TH2* h = hit->second;
00682     c1->cd(); c1->Clear();
00683     h->Draw("colz");
00684     c1->Update();
00685     c1->Print(Form("%s/_%s.gif",outputDirectory,h->GetName()));
00686     h->Write();
00687   }
00688   f->Write();
00689 }


Generated on 2 Nov 2017 for loon by  doxygen 1.6.1