loon: FitGC.cxx Source File

00001 
00002 // $Id: FitGC.cxx,v 1.14 2007/01/15 21:05:47 rhatcher Exp $
00003 //
00004 // FILL_IN: [Document your code!!]
00005 //
00006 //S. Murgia, T.J. Yang
00008 #include "FitGC.h"
00009 #include "GCSummary.h"
00010 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00011 #include "HistMan/HistMan.h"
00012 #include "Validity/VldRange.h"
00013 
00014 #include "TFile.h"
00015 #include "TBranch.h"
00016 #include "TCanvas.h"
00017 #include "TF1.h"
00018 #include "TError.h"
00019 #include "TStyle.h"
00020 
00021 #include <iostream>
00022 #include <cmath>
00023 
00024 using namespace std;
00025 
00026 ClassImp(FitGC)
00027 
00028 
00029 //......................................................................;
00030 
00031 FitGC::FitGC()
00032   : fittime(0)   
00033   , fDetector(0)
00034   , DoPlots(0)
00035   , xmin(100.)
00036   , xmax(8000.)
00037   , idebug(0)
00038   , ibadfit(0)
00039 {
00040   gStyle->SetOptFit(0111);
00041   gStyle->SetStatX(0.6);
00042   gStyle->SetStatY(0.9);
00043   gStyle->SetStatH(0.1);
00044   gc = new GCSummary();
00045   //fDetector = 0;
00046   pvp = new TGraphErrors(40);
00047   rvpmt = new TGraphErrors(40);
00048   rvpin = new TGraphErrors(40);
00049   fit_formula = new TF1("fit_formula","[0]*((x)-[1])");
00050 }
00051 
00052 //......................................................................
00053 
00054 FitGC::~FitGC()
00055 {
00056   cout<<"Finished."<<endl;
00057 }
00058 
00059 void FitGC::writehm()
00060 {
00061   char hn[100];
00062   sprintf(hn,"%s_plots.root",filename);
00063   hm->WriteOut(hn);
00064 }  
00065 
00066 //.......................................................................
00067 
00068 void FitGC::OpenFile(const char* name)
00069 {
00070   TFile *f = new TFile(name);
00071   gc_tree = 0;
00072   gc_tree = (TTree *)f->Get("gc_tree");
00073   if (!gc_tree) exit(0);
00074   TBranch *branch = gc_tree->GetBranch("");
00075   branch->SetAddress(&gc);
00076   cout<<gc_tree->GetEntries()<<endl;
00077   gc_tree->GetEntry(0);
00078   fDetector = gc->dettype;
00079   //I don't know what's best way to do this.
00080   char str[] = ".root";
00081   int i = strcspn(name,str);
00082   strncpy(filename,name,i);
00083   filename[i]='\0';
00084   cout<<filename<<endl;
00085   char hmn[100];
00086   sprintf(hmn,"%sall",filename);
00087   hm = new HistMan(hmn);
00088   hm->Book<TH1F>("pmtoffset","pmtoffset",1500,-2000,600);
00089   hm->Book<TH1F>("pinoffset","pinoffset",150,-90,60);
00090   hm->Book<TH1F>("slope","slope",100,-10,30);
00091 
00092 //  char hn[100];
00093 //  for (int i = 0; i<40; i++){
00094 //    sprintf(hn,"pmtoffset,pin0=%f",-40.+i*0.5);
00095 //    hm->Book<TH1F>(hn,hn,100,-1000,1000);
00096 //  }
00097 
00098   hm->Book<TH1F>("fitpoints","fitpoints",41,0,41);
00099 
00100   char htitle1[100];
00101   sprintf(htitle1,"res vs pinadc,problem channels");
00102   char htitle2[100];
00103   sprintf(htitle2,"res vs pmtadc,problem channels");
00104   if (!hm->Get<TH2F>(htitle1)){
00105     hm->Book<TH2F>(htitle1,htitle1,1000,0,5000,100,-20,20);
00106     hm->Book<TH2F>(htitle2,htitle2,1000,0,50000,100,-20,20);
00107   }
00108 
00109   //cout<<gc->numPointsPmt<<endl;
00110 }
00111 
00112 void FitGC::Plot(Int_t plane, Int_t strip)
00113 {
00114   TCanvas *c1 = new TCanvas("c1","c1",750,250);
00115 
00116   for (Int_t i = 0; i<gc_tree->GetEntries(); i++){
00117     gc_tree->GetEntry(i);
00118     if (i%1000==0) cout<<i<<"/"<<gc_tree->GetEntries()<<"\r";
00119     if (gc->highlow==0) continue;
00120     if (gc->plane==plane && gc->strip==strip){
00121       c1->Clear();
00122       c1->Divide(3,1);
00123       this->DoFit();
00124       c1->cd(1);
00125       pvp->Draw("ap");
00126       c1->cd(2);
00127       rvpmt->Draw("ap");
00128       c1->cd(3);
00129       rvpin->Draw("ap");
00130     }
00131   }
00132   cout<<endl;
00133 }   
00134 
00135 void FitGC::PlotPL(Int_t ipb, Int_t iled)
00136 {
00137 
00138   gErrorIgnoreLevel = 2;
00139   gStyle->SetTitleSize(0.06,"xyz");
00140 
00141   TCanvas *c1 = new TCanvas("c1","c1",600,600);
00142   char fn1[100];
00143   char fn2[100];
00144   char fn3[100];
00145   char fnd1[100];
00146   char fnd2[100];
00147   char fnd3[100];
00148 
00149   if (ipb == -1 && iled == -1){ //plot everything
00150     sprintf(fn1,"%s_all.ps[",filename);
00151     sprintf(fn2,"%s_all.ps",filename);
00152     sprintf(fn3,"%s_all.ps]",filename);
00153     sprintf(fnd1,"%sdebug_all.ps[",filename);
00154     sprintf(fnd2,"%sdebug_all.ps", filename);
00155     sprintf(fnd3,"%sdebug_all.ps]",filename);
00156   }
00157   else if(ipb != -1 && iled == -1){//plot one pulserbox
00158     sprintf(fn1,"%s_pb%d.ps[",filename,ipb);
00159     sprintf(fn2,"%s_pb%d.ps",filename,ipb);
00160     sprintf(fn3,"%s_pb%d.ps]",filename,ipb);
00161     sprintf(fnd1,"%sdebug_pb%d.ps[",filename,ipb);
00162     sprintf(fnd2,"%sdebug_pb%d.ps", filename,ipb);
00163     sprintf(fnd3,"%sdebug_pb%d.ps]",filename,ipb);
00164   }    
00165   else if(ipb == -1 && iled != -1){//plot one led
00166     sprintf(fn1,"%s_led%d.ps[",filename,iled);
00167     sprintf(fn2,"%s_led%d.ps",filename,iled);
00168     sprintf(fn3,"%s_led%d.ps]",filename,iled);
00169     sprintf(fnd1,"%sdebug_led%d.ps[",filename,iled);
00170     sprintf(fnd2,"%sdebug_led%d.ps", filename,iled);
00171     sprintf(fnd3,"%sdebug_led%d.ps]",filename,iled);
00172   }    
00173   else if(ipb != -1 && iled != -1){//plot one led in one pulserbox
00174     sprintf(fn1,"%s_pb%dled%d.ps[",filename,ipb,iled);
00175     sprintf(fn2,"%s_pb%dled%d.ps",filename,ipb,iled);
00176     sprintf(fn3,"%s_pb%dled%d.ps]",filename,ipb,iled);
00177     sprintf(fnd1,"%sdebug_pb%dled%d.ps[",filename,ipb,iled);
00178     sprintf(fnd2,"%sdebug_pb%dled%d.ps", filename,ipb,iled);
00179     sprintf(fnd3,"%sdebug_pb%dled%d.ps]",filename,ipb,iled);
00180   } 
00181 
00182   if (DoPlots) {
00183     c1->Print(fn1);
00184     c1->Print(fnd1);
00185   }
00186   for (Int_t i = 0; i<gc_tree->GetEntries(); i++){
00187     gc_tree->GetEntry(i);
00188     if (i%1000==0) cout<<i<<"/"<<gc_tree->GetEntries()<<endl;
00189     
00190     if(gc->highlow==0&&gc->pb!=1) continue;  //only use low gain pin diode
00191     
00192     char htitle1[100];
00193     sprintf(htitle1,"PB%d Led%d Gain%d res vs pinadc",gc->pb,gc->led,gc->highlow);
00194     char htitle2[100];
00195     sprintf(htitle2,"PB%d Led%d Gain%d res vs pmtadc",gc->pb,gc->led,gc->highlow);
00196     
00197     char htitle3[100];
00198     sprintf(htitle3,"PB%d Led%d Gain%d PinOffset",gc->pb,gc->led,gc->highlow);
00199     char htitle4[100];
00200     sprintf(htitle4,"PB%d Led%d Gain%d PmtOffset",gc->pb,gc->led,gc->highlow);
00201 
00202     if (!hm->Get<TH2F>(htitle1)){
00203       cout<<"booking histogram, pb "<<gc->pb<<" led "<<gc->led<<endl;
00204       hm->Book<TH2F>(htitle1,htitle1,1000,0,5000,100,-20,20);
00205       hm->Book<TH2F>(htitle2,htitle2,1000,0,50000,100,-20,20);
00206       hm->Book<TH1F>(htitle3,htitle3,150,-90,60);
00207       hm->Book<TH1F>(htitle4,htitle4,1000,-2000,100);
00208     }
00209 
00210     if ((gc->pb==ipb||ipb==-1) && (gc->led==iled||iled==-1)){
00211       this->DoFit();
00212       if (DoPlots){
00213         c1->Clear();
00214         c1->Divide(2,2);
00215         c1->cd(1);
00216         pvp->Draw("ap");
00217         c1->cd(2);
00218         rvpmt->Draw("ap");
00219         c1->cd(3);
00220         rvpin->Draw("ap");
00221         c1->Print(fn2);
00222         if (idebug) c1->Print(fnd2);
00223       }
00224     }
00225   }
00226   if (DoPlots) {
00227     c1->Print(fn3);
00228     c1->Print(fnd3);
00229   }
00230 }
00231 
00232 void FitGC::DoFit()
00233 {
00234   if (gc->numPointsPin!=gc->numPointsPmt){
00235     cout<<"different points"<<endl;
00236     exit(0);
00237   }
00238   
00239   pvp->Clear();
00240   rvpin->Clear();
00241   rvpmt->Clear();
00242 
00243   idebug = 0;
00244   ibadfit = 0;
00245 
00246   for (int i = 0; i<gc->numPointsPin; i++){
00247     pvp->SetPoint(i,gc->meanPin[i],gc->meanPmt[i]);
00248     pvp->SetPointError(i,gc->errorPin[i],gc->errorPmt[i]);
00249   }
00250 
00251   float xmin_fit = xmin;
00252   float xmax_fit = xmax;
00253   int imax_save = 0;
00254 
00255   for (int i = 0; i<gc->numPointsPin; i++) {
00256     if (gc->meanPmt[i] < xmin && gc->meanPin[i] > xmin) 
00257       xmin_fit =gc->meanPin[i]+1;
00258     else {
00259       if (gc->meanPmt[i] > xmax && gc->meanPin[i] < xmax) {
00260         xmax_fit = gc->meanPin[i]+1;
00261         imax_save = i;
00262         break;
00263       }
00264     }
00265   }
00266 
00267   if (fittime==0||fittime==1){
00268     fit_formula->SetParameter(0,1);
00269     fit_formula->SetParameter(1,20);
00270   }
00271   if (fittime==0||fittime==2){
00272     fit_formula->SetParameter(0,1);
00273   }
00274 
00275   fit_formula->SetRange(xmin_fit,xmax_fit);
00276   //TF1 f("line","[0]*(x)+[1]",xmin_fit,xmax_fit);
00277   //f.FixParameter(1,0);
00278   pvp->Fit("fit_formula","QR");
00279 
00280   Int_t fitpoints = fit_formula->GetNumberFitPoints();
00281   //cout<<"Fit Points "<<fitpoints<<endl;
00282 //  if (xmin_fit>xmax_fit) {
00283 //    cout<<"xmin_fit "<<xmin_fit<<" xmax_fit "<<xmax_fit<<" chi2 "<<fit_formula->GetChisquare()<<" "<<fitpoints<<endl;
00284 //  }
00285 
00286   if (fitpoints==0 || fit_formula->GetChisquare()==0){//fitting was failed, try something else
00287     int ipin_min = 0;
00288     int ipin_max = -1;
00289     for (int i = 0; i<gc->numPointsPin-1; i++) {
00290       if (gc->meanPmt[i]<xmin && gc->meanPmt[i+1]>=xmin)
00291         ipin_min = i+1;
00292       if (gc->meanPmt[i]<=xmax && gc->meanPmt[i+1]>xmax)
00293         ipin_max = i;
00294     }
00295     int npot = 3;
00296     if (ipin_max - ipin_min + 1 >= npot){//we have npot points at least
00297       for (int i = 0; i<ipin_max - ipin_min + 1 - npot; i++){
00298         if (i>1) continue;
00299         if (fittime==0||fittime==1){
00300           fit_formula->SetParameter(0,1);
00301           fit_formula->SetParameter(1,20);
00302         }
00303         if (fittime==0||fittime==2){
00304           fit_formula->SetParameter(0,1);
00305         }
00306         fit_formula->SetRange(gc->meanPin[ipin_max-npot+1-i]-1,gc->meanPin[ipin_max]+1);
00307         pvp->Fit("fit_formula","QR");
00308         //cout<<"i="<<i<<" slope="<<fit_formula->GetParameter(0)<<" offset="<<fit_formula->GetParameter(1)<<endl;
00309       }
00310     }
00311   }//finish fitting
00312 
00313   fitpoints = fit_formula->GetNumberFitPoints();
00314   if (fitpoints<3) idebug = 1;
00315   if (fit_formula->GetChisquare()==0) idebug = 1;
00316 
00317   if (fitpoints<3){//the fitting did not fail, however, there are too few points in the linear region, we better include more (at least 3 for now)
00318     for (int i = imax_save; i>=0; i--){
00319       if (gc->meanPin[i]<200) continue; //won't go below 200 in PIN
00320       if (imax_save -i>2) continue;//no more than 3 points
00321       xmin_fit = gc->meanPin[i]-1;
00322     }
00323     if (fittime==0||fittime==1){
00324       fit_formula->SetParameter(0,1);
00325       fit_formula->SetParameter(1,20);
00326     }
00327     if (fittime==0||fittime==2){
00328       fit_formula->SetParameter(0,1);
00329     }
00330 
00331     fit_formula->SetRange(xmin_fit,xmax_fit);
00332     pvp->Fit("fit_formula","QR");
00333   }//finish fitting
00334 
00335   fitpoints = fit_formula->GetNumberFitPoints();
00336 
00337   if (fittime==0||fittime==2){
00338     if(hm->Get<TH1F>("fitpoints")) hm->Fill1d("fitpoints",fitpoints);
00339   }
00340 
00341   char grtitle[100];
00342   sprintf(grtitle,"PMT vs PIN,PB%d Led%d plane%d strip%d Gain%d",gc->pb,gc->led,gc->plane,gc->strip,gc->highlow);
00343   pvp->SetTitle(grtitle);
00344   pvp->SetName(grtitle);  
00345   pvp->SetLineWidth(1);
00346   pvp->SetMarkerStyle(20);
00347   pvp->SetMarkerColor(2);
00348   pvp->SetMarkerSize(0.8);
00349   pvp->GetXaxis()->SetTitle("PIN");
00350   pvp->GetYaxis()->SetTitle("PMT");
00351   pvp->GetXaxis()->SetTitleSize(0.05);
00352   pvp->GetYaxis()->SetTitleSize(0.05);
00353   pvp->GetXaxis()->SetLabelSize(0.04);
00354   pvp->GetYaxis()->SetLabelSize(0.04);
00355   pvp->GetYaxis()->SetTitleOffset(2);
00356 
00357   float invslope;
00358   float sloperr;
00359   float offset;
00360   float offseterr;
00361 
00362   if (fit_formula->GetParameter(0)!=0) {
00363     invslope = 1.0/fit_formula->GetParameter(0);
00364     sloperr = fit_formula->GetParError(0);
00365     offset = fit_formula->GetParameter(1);
00366     offseterr = fit_formula->GetParError(1);
00367 
00368     if (fittime==0||fittime==2){
00369       if(hm->Get<TH1F>("pmtoffset")) hm->Fill1d("pmtoffset",offset/invslope);
00370       if(hm->Get<TH1F>("pinoffset")) hm->Fill1d("pinoffset",offset);
00371       if(hm->Get<TH1F>("slope")) hm->Fill1d("slope",fit_formula->GetParameter(0));
00372     }
00373   } else {
00374     invslope = -9999.;
00375     sloperr = 0.;
00376     offset = -9999.;
00377     offseterr = 0.;
00378   }
00379 
00380   slope = 1./invslope;
00381   intercept = offset;
00382 
00383   char htitle1[100];
00384   sprintf(htitle1,"PB%d Led%d Gain%d res vs pinadc",gc->pb,gc->led,gc->highlow);
00385   char htitle2[100];
00386   sprintf(htitle2,"PB%d Led%d Gain%d res vs pmtadc",gc->pb,gc->led,gc->highlow);
00387   
00388   char htitle3[100];
00389   sprintf(htitle3,"PB%d Led%d Gain%d PinOffset",gc->pb,gc->led,gc->highlow);
00390   char htitle4[100];
00391   sprintf(htitle4,"PB%d Led%d Gain%d PmtOffset",gc->pb,gc->led,gc->highlow);
00392 
00393   char htitle5[100];
00394   sprintf(htitle5,"res vs pinadc,problem channels");
00395   char htitle6[100];
00396   sprintf(htitle6,"res vs pmtadc,problem channels");
00397 
00398   float residual = 0;
00399   float residualerr = 0;
00400   
00401   for (int i=0; i< gc->numPointsPin; i++) {
00402     float stripendmeani = gc->meanPmt[i];
00403     float stripenderri = gc->errorPmt[i];
00404     float pinmeani =  gc->meanPin[i];
00405     float pinerri =  gc->errorPin[i];
00406     //Get mean from fit = slope*pinmean
00407     float meanfit = (pinmeani-offset)/invslope;
00408     //float meanfit = (pinmeani+9.2)/invslope;
00409     float meanfiterr = pow(pinmeani,2)*pow(sloperr,2)+pow(1./invslope,2)*pow(pinerri,2) + pow(sloperr,2)*pow(offset,2)+pow(1./invslope,2)*pow(offseterr,2);
00410     
00411     residual = 100.*(meanfit-stripendmeani)/meanfit;
00412     if (i>gc->numPointsPin/2 && residual < -100){//intent to remove bad fits,pretty loose
00413       ibadfit = 1;
00414     }
00415 
00416     if (meanfit == 0.) { 
00417       residualerr = 0.;
00418     }   else {
00419       residualerr = 100.*sqrt(pow(stripenderri,2)/pow(meanfit,2) + 
00420                               (pow(stripendmeani,2)/pow(meanfit,4))*meanfiterr);
00421     }
00422     //if (ipoint>3 && residual<-5) cout<<ledid.GetPulserBox()<<","<<ledid.GetLedInBox()<<","<<planeno<<","<<stripno<<endl;
00423 //    if (i>1&&residual<-3) idebug = 1;
00424 //    if (i>1&&stripendmeani<8000&&residual>3) idebug = 1;
00425 //    if (offset*invslope>40||offset*invslope<-10) idebug = 1;
00426     rvpin->SetPoint(i, pinmeani,residual);
00427     rvpin->SetPointError(i, pinerri,residualerr);
00428     rvpmt->SetPoint(i, stripendmeani,residual);
00429     rvpmt->SetPointError(i, stripenderri,residualerr);
00430     if (fittime==0||fittime==2){
00431       if(hm->Get<TH2F>(htitle1)) hm->Fill2d(htitle1,pinmeani,residual);
00432       if(hm->Get<TH2F>(htitle2)) hm->Fill2d(htitle2,stripendmeani,residual);
00433       if (idebug){
00434         if(hm->Get<TH2F>(htitle5)) hm->Fill2d(htitle5,pinmeani,residual);
00435         if(hm->Get<TH2F>(htitle6)) hm->Fill2d(htitle6,stripendmeani,residual);
00436       }
00437     }
00438   }  
00439   if (fittime==0||fittime==2){
00440     if(hm->Get<TH1F>(htitle4)) hm->Fill1d(htitle4,offset/invslope);
00441   }
00442   if (fittime==0||fittime==1){
00443     //cout<<"offset "<<offset<<endl;
00444     if(hm->Get<TH1F>(htitle3)) hm->Fill1d(htitle3,offset);
00445   }
00446 
00447   rvpin->SetPoint(0,0,0);
00448   rvpin->SetPointError(0,0,0);
00449   rvpmt->SetPoint(0,0,0);
00450   rvpmt->SetPointError(0,0,0);
00451 
00452   sprintf(grtitle,"Residual vs PMT,PB%d Led%d plane%d strip%d Gain%d",gc->pb,gc->led,gc->plane,gc->strip,gc->highlow);
00453   rvpmt->SetTitle(grtitle);
00454   rvpmt->SetName(grtitle);  
00455   rvpmt->SetLineWidth(1);
00456   rvpmt->SetMarkerStyle(20);
00457   rvpmt->SetMarkerColor(2);
00458   rvpmt->SetMarkerSize(0.8);
00459   rvpmt->GetYaxis()->SetTitle("Residual");
00460   rvpmt->GetXaxis()->SetTitle("PMT");
00461   rvpmt->GetXaxis()->SetTitleSize(0.05);
00462   rvpmt->GetYaxis()->SetTitleSize(0.05);
00463   rvpmt->GetXaxis()->SetLabelSize(0.025);
00464   //if (gc->meanPmt[19]<20000) rvpmt->GetXaxis()->SetLabelSize(0.03);
00465   rvpmt->GetYaxis()->SetLabelSize(0.04);
00466   rvpmt->GetYaxis()->SetTitleOffset(2);
00467 
00468   sprintf(grtitle,"Residual vs PIN,PB%d Led%d plane%d strip%d Gain%d",gc->pb,gc->led,gc->plane,gc->strip,gc->highlow);
00469   rvpin->SetTitle(grtitle);
00470   rvpin->SetName(grtitle);  
00471   rvpin->SetLineWidth(1);
00472   rvpin->SetMarkerStyle(20);
00473   rvpin->SetMarkerColor(2);
00474   rvpin->SetMarkerSize(0.8);
00475   rvpin->GetYaxis()->SetTitle("Residual");
00476   rvpin->GetXaxis()->SetTitle("PIN");
00477   rvpin->GetXaxis()->SetTitleSize(0.05);
00478   rvpin->GetYaxis()->SetTitleSize(0.05);
00479   rvpin->GetXaxis()->SetLabelSize(0.04);
00480   rvpin->GetYaxis()->SetLabelSize(0.04);
00481   rvpin->GetYaxis()->SetTitleOffset(2);
00482 
00483 
00484 }
00485 
00486 //........................................................................
00487 // made change here
00488 void FitGC::GetSiglin(Float_t pmtadc,Int_t pb, Int_t led) {
00489 
00490   this->DrawPlots(0);
00491   fittime = 1;
00492   static int first = 1;
00493   if (first){
00494     this->PlotPL(pb,led);
00495     first = 0;
00496   }
00497        
00498   //this->DrawPlots(1);
00499   
00500   char htitle1[100];
00501   char htitle2[100];
00502   if (pb==-1&&led==-1){
00503     sprintf(htitle1,"Siglin for rawADC:%f",pmtadc);
00504     sprintf(htitle2,"rawADC:%f siglin:raw",pmtadc);
00505   }
00506   else if (pb==-1&&led!=-1){
00507   sprintf(htitle1,"Led%d Siglin for rawADC:%f",led,pmtadc);
00508   sprintf(htitle2,"Led%d rawADC:%f siglin:raw",led,pmtadc);
00509   }    
00510   else if (pb!=-1&&led==-1){
00511   sprintf(htitle1,"PB%d Siglin for rawADC:%f",pb,pmtadc);
00512   sprintf(htitle2,"PB%d rawADC:%f siglin:raw",pb,pmtadc);
00513   }
00514   else if (pb!=-1&&led!=-1){
00515     sprintf(htitle1,"PB%d Led%d Siglin for rawADC:%f",pb,led,pmtadc);
00516     sprintf(htitle2,"PB%d Led%d rawADC:%f siglin:raw",pb,led,pmtadc);
00517   }
00518   if (!hm->Get<TH1F>(htitle1)){
00519       cout<<"booking histogram, pb "<<pb<<" led "<<led<<endl;
00520       hm->Book<TH1F>(htitle1,htitle1,5000,-1000,17000);
00521       hm->Book<TH1F>(htitle2,htitle2,1000,0,2);
00522   }
00523   
00524   TCanvas *c1 = new TCanvas("c1","c1",600,600);
00525   char fn1[100];
00526   char fn2[100];
00527   char fn3[100];
00528   sprintf(fn1,"%s_pb%d_led%d.ps[",filename,pb,led);
00529   sprintf(fn2,"%s_pb%d_led%d.ps",filename,pb,led);
00530   sprintf(fn3,"%s_pb%d_led%d.ps]",filename,pb,led);
00531   if (DoPlots) c1->Print(fn1);
00532 
00533   TGraph *gra = new TGraph(40);
00534 
00535   fittime = 2;
00536 
00537   cout<<"Calculate siglin for raw adc: "<<pmtadc<<endl;
00538 
00539 
00540   for (Int_t i = 0; i<gc_tree->GetEntries(); i++){//for
00541     gc_tree->GetEntry(i);
00542     if (i%1000==0) cout<<i<<"/"<<gc_tree->GetEntries()<<"\r";
00543     if(gc->highlow==0&&gc->pb!=1) continue;
00544     if (!(gc->led>=9 && gc->pb==1 && gc->highlow==0)) continue;    
00545     if ((gc->pb==pb||pb==-1) && (gc->led==led||led==-1)){
00546       //cout<<gc->stripend<<endl;
00547       char htitle3[100];
00548       sprintf(htitle3,"PB%d Led%d Gain%d PinOffset",gc->pb,gc->led,gc->highlow);
00549       Float_t meanpinoffset=0;
00550       if(hm->Get<TH1F>(htitle3)) 
00551         meanpinoffset = hm->Get<TH1F>(htitle3)->GetMean();
00552       fit_formula->SetParameter(1,meanpinoffset);
00553       fit_formula->SetParLimits(1,40,40);
00554       this->DoFit();
00555 
00556       for(Int_t ipoint=0; ipoint< gc->numPointsPmt; ipoint++ ){
00557         gra->SetPoint(ipoint,gc->meanPmt[ipoint],gc->meanPin[ipoint]);
00558         //if (gc->stripend==1616) cout<<ipoint<<" "<<gc->meanPmt[ipoint]<<" "<<gc->meanPin[ipoint]<<" "<<fit_formula->GetParameter(0)<<" "<<fit_formula->GetParameter(1)<<endl;
00559       }
00560 
00561       if (DoPlots){
00562         c1->Clear();
00563         c1->Divide(2,2);
00564         c1->cd(1);
00565         pvp->Draw("ap");
00566         c1->cd(2);
00567         rvpmt->Draw("ap");
00568         c1->cd(3);
00569         rvpin->Draw("ap");
00570         c1->Print(fn2);
00571       }
00572       
00573       // get siglin
00574       if(hm->Get<TH1F>(htitle1)) 
00575         hm->Fill1d(htitle1, fit_formula->Eval(gra->Eval(pmtadc)));
00576       if(hm->Get<TH1F>(htitle2)&&pmtadc) {
00577         hm->Fill1d(htitle2, fit_formula->Eval(gra->Eval(pmtadc))/pmtadc);
00578         //if (gc->stripend==1616) cout<<gra->Eval(pmtadc)<<" "<<fit_formula->GetParameter(0)<<" "<<fit_formula->GetParameter(1)<<" "<<fit_formula->Eval(gra->Eval(pmtadc))/pmtadc<<endl;
00579       }
00580     }
00581   }
00582   cout<<endl;
00583   if (DoPlots) c1->Print(fn3);
00584 }
00585 
00586 //.....................................................................
00587 void FitGC::WriteDB(const Detector::Detector_t &det, std::string fDBName){
00588   
00589   this->DrawPlots(0);
00590 
00591   fittime = 1;
00592 
00593   this->PlotPL();
00594 
00595   fittime = 2;
00596 
00597   this->DrawPlots(0);
00598 
00599   TCanvas *c1 = new TCanvas("c1","c1",600,600);
00600   char fn1[100];
00601   char fn2[100];
00602   char fn3[100];
00603   sprintf(fn1,"%s.ps[",filename);
00604   sprintf(fn2,"%s.ps",filename);
00605   sprintf(fn3,"%s.ps]",filename);
00606   if (DoPlots) c1->Print(fn1);
00607 
00608   //first look at high gain (pingain=0), then low gain (pingain=1)
00609   for (int pingain=0; pingain<2; pingain++){
00610     int aggre = -1;
00611     for (Int_t i = 0; i<gc_tree->GetEntries(); i++){
00612       gc_tree->GetEntry(i);
00613       if (i%1000==0) cout<<i<<"/"<<gc_tree->GetEntries()<<" pin gain "<<pingain<<endl;
00614       if(gc->highlow!=pingain) continue;
00615       if(gc->highlow==0&&gc->pb==0) continue;//right now only low gain pin for pb0
00616       if(gc->highlow==0&&gc->pb==2) continue;//right now only low gain pin for pb2
00617       char htitle3[100];
00618       sprintf(htitle3,"PB%d Led%d Gain%d PinOffset",gc->pb,gc->led,gc->highlow);
00619       Float_t meanpinoffset=0;
00620       if(hm->Get<TH1F>(htitle3)) meanpinoffset = hm->Get<TH1F>(htitle3)->GetMean();
00621       //cout<<"offset "<<meanpinoffset<<endl;
00622       fit_formula->SetParameter(1,meanpinoffset);
00623       fit_formula->SetParLimits(1,40,40);
00624       
00625       this->DoFit();
00626       if (DoPlots){
00627         c1->Clear();
00628         c1->Divide(2,2);
00629         c1->cd(1);
00630         pvp->Draw("ap");
00631         c1->cd(2);
00632         rvpmt->Draw("ap");
00633         c1->cd(3);
00634         rvpin->Draw("ap");
00635         c1->Print(fn2);
00636       }
00637       if (gc->aggno!=aggre){
00638         aggre = gc->aggno;
00639         cout<<"aggregate no "<<aggre<<endl;
00640         VldRange range(det,SimFlag::kData,gc->timeStart,gc->timeEnd,fDBName);
00641         if (gc->highlow==1){//low gain
00642           //fWriter.Open(range,aggre,1,gc->timeStart,fDBName);
00643           fWriter.Open(range,aggre,1);
00644           fWriter.SetDbName(fDBName);
00645         }
00646         else if (gc->highlow==0){//high gain
00647           //fWriter.Open(range,aggre,3,gc->timeStart,fDBName);
00648           fWriter.Open(range,aggre,3);
00649           fWriter.SetDbName(fDBName);
00650         }
00651         
00652         fWriter.SetOverlayCreationDate();
00653       }
00654       float hslope = fit_formula->GetParameter(0);
00655       float hslopeerr = fit_formula->GetParError(0);
00656       float hxoffset = fit_formula->GetParameter(1);
00657       float hchisq = fit_formula->GetChisquare();
00658       if (ibadfit) hchisq = 0;
00659       //this channel is not behavioring well
00660       //if (gc->stripend==6548) hchisq = 0;
00661       int hnumpoints = fit_formula->GetNumberFitPoints();    
00662       
00663       CalPulserFits hrow(aggre,gc->stripend,hnumpoints,hslope,hslopeerr,hxoffset,hchisq);
00664       fWriter << hrow;
00665     }
00666   }
00667     
00668   fWriter.Close();
00669     
00670   if (DoPlots) c1->Print(fn3);
00671 }