StRoot: StEmcPool/StBarrelCalib2008/BprsCapPolygraph.cxx Source File

00001 #include <stdio.h>
00002 #include <TH2F.h>
00003 #include <TObjArray.h>
00004 #include <TGraph.h>
00005 
00006 #include "BprsCapPolygraph.h"
00007 #include "StJanBarrelDbMaker.h"
00008 #include "JanBarrelEvent.h"
00009 
00010 
00011 //________________________________________________
00012 //________________________________________________
00013 BprsCapPolygraph::BprsCapPolygraph( TObjArray *HList, StJanBarrelDbMaker* jdb, int pedFlag) {
00014     par_pedFlag=pedFlag;
00015     mJanDbMaker=jdb; assert(mJanDbMaker);
00016 
00017     const char *fname="all"; //tmp
00018     //QA histos
00019 
00020     hChiB=new TH1F("bppg_chb","BprsPoly logN(chi2/dof), nominal capID ; logN(chi2/dof)",200,-1.,9.);
00021     hChiGap=new TH1F("bppg_chgI","BprsPoly chi2/dof (nominal - best_capID), any minimum; logN( nominal_chi2 - best_chi2 ) ",200,-6,10);
00022     hChiGap2=new TH1F("bppg_chgA","BprsPoly chi2/dof (nominal - best_capID), acceted; logN( nominal_chi2 - best_chi2 ) ",200,-6,10);
00023       
00024     char tt2[1000];
00025     hCh2D=new TH2F("bppg_ch2D","BprsPoly chi2/DOF(best capID.NE.nominalCapID) ; LogN( chi2/DOF), best; logN(chi2/DOF), nominal ",50,-1,10,50,-1,10);
00026     sprintf(tt2,"BprsPoly caps in sync, %s; nominal capID; BPRS crate ID",fname);
00027     hCapGood=new TH2F("bppg_capSyn",tt2, 128,-0.5,127.5,mxBprsCrate,-0.5,mxBprsCrate-0.5);
00028     sprintf(tt2,"BprsPoly caps NOT sync, %s; nominal capID; BPRS crateID",fname);
00029     hCapCorr=new TH2F("bppg_capDeSyn",tt2, 128,-0.5,127.5,mxBprsCrate,-0.5,mxBprsCrate-0.5);
00030 
00031     HList->Add( hChiB);
00032     HList->Add( hChiGap);
00033     HList->Add( hChiGap2);
00034     HList->Add( hCh2D);
00035     HList->Add(hCapGood);
00036     HList->Add(hCapCorr);
00037  
00038    int nb=200;
00039     float adc1=-50;
00040     float adc2=adc1+nb;
00041 
00042     hAdcGood=new TH1F("bppg_adcGd","BprsPoly ADC for capID in sync; rawAdc - ped(nominal capID)",nb,adc1,adc2);
00043     hAdcCorr=new TH1F("bppg_adcCor","BprsPoly ADC for capID NOT in sync,fixed; rawAdc - ped(best capID)",nb,adc1,adc2);
00044     hAdcCorr->SetLineColor(kRed);
00045 
00046     HList->Add(hAdcGood);
00047     HList->Add( hAdcCorr);
00048    
00049   }
00050 
00051 //________________________________________________
00052 //________________________________________________
00053 void BprsCapPolygraph::doBaseline( JanBprsEveA & bprsEve, JanBarrelEvent &fullEve){
00054   
00055   int crateID=bprsEve.crateID;
00056   int   capID=bprsEve.capID;
00057   printf("doBaseline: bprsEve: %d rawADC, doCrate=%d nominalCapID=%d\n", bprsEve.raw.GetN(),crateID,capID);
00058   TGraph res; //residua:  X=ADC-ped, y=index of raw; 
00059   
00060   int ibp=kBPrs;
00061   for(int i=0;i< bprsEve.raw.GetN();i++) {
00062     double rawAdc, yid;
00063     bprsEve.raw.GetPoint(i,rawAdc,yid);
00064     int id= (int) yid;
00065     // printf("K i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
00066     int   stat =mJanDbMaker->statTile(ibp,id);
00067     if(stat) continue;
00068     float ped  =mJanDbMaker->pedTile(ibp,id,capID);
00069     float sig  =mJanDbMaker->sigPedTile(ibp,id,capID);
00070     double del=rawAdc-ped;
00071     if(del> cut_adcMax) continue;
00072     //  if(i<20)
00073     //printf("i=%d id=%d del=%f sig=%f\n",i,id,del,sig);
00074     assert(sig!=0);
00075     del/=sig;
00076     res.SetPoint(res.GetN(),del,i);
00077   } // end of loop over raw data
00078   
00079   res.Sort();
00080   int i2=(int)(res.GetN()*cut_fracSkip);
00081   int i1=res.GetN()-i2;
00082   printf("res size=%d , use=[%d,%d]\n", res.GetN(),i1,i2);
00083   // res.Print();
00084   
00085   double sum=0;
00086   for(int i=i1; i<i2;i++) {
00087     double rdel, yj;
00088     res.GetPoint(i,rdel,yj);
00089     sum+=rdel*rdel; // compute base chi2/dof
00090     //if(i<20) printf("CC i=%d rdel=%f sig=%f  sum=%f\n",i,rdel,sig,sum);    
00091     // prepare array of index for chi2 minimalization
00092     int j= (int) yj; // index in 'raw' array
00093     double rawAdc, yid;
00094     bprsEve.raw.GetPoint(j,rawAdc,yid);
00095     bprsEve.addGoodValue(yid,rawAdc);
00096   }// base chi2 computed 
00097   int dof=i2-i1;
00098   float chi2dof=sum/dof;
00099   printf(" base chi2=%f DOF=%d, chi/DOF=%f nGood=%d\n",sum,dof,chi2dof,bprsEve.good.GetN());
00100   bprsEve.chi2=sum;
00101   bprsEve.chi2dof=chi2dof;
00102   hChiB->Fill(log(chi2dof));
00103 }
00104 
00105 //________________________________________________
00106 //________________________________________________
00107 void BprsCapPolygraph::findBestCap(JanBprsEveA & bprsEve, JanBarrelEvent &fullEve){
00108 
00109    printf("find BPRS bestCap: eve: %d good\n", bprsEve.good.GetN());
00110    int ibp=kBPrs;
00111  
00112    float minChi2dof=bprsEve.chi2dof-0.0001;
00113    int bestCapID=-1;
00114 
00115    int cap1=bprsEve.capID-par_mxDelCap+mxBcap;
00116    int cap2=bprsEve.capID+par_mxDelCap+mxBcap;
00117    for(int ic=cap1; ic<=cap2;ic++) {
00118      int cap=ic%mxBcap;
00119      // printf("xx cap=%d\n",cap);
00120      double sum=0;
00121      for(int i=0;i< bprsEve.good.GetN();i++) {
00122       double rawAdc, yid;
00123       bprsEve.good.GetPoint(i,rawAdc,yid);
00124       int id= (int) yid;
00125       // printf("bK i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
00126       int   stat =mJanDbMaker->statTile(ibp,id);
00127 
00128       assert(stat==0);// bad channels were already discarded
00129       float ped  =mJanDbMaker->pedTile(ibp,id,cap);
00130       float sig  =mJanDbMaker->sigPedTile(ibp,id,cap);
00131       double del=rawAdc-ped;
00132       //if(i<20) printf("i=%d id=%d del=%f sig=%f\n",i,id,del,sig);
00133       del/=sig;
00134       sum+=del*del;
00135      } // end of data loop 
00136      float chi2dof=sum/bprsEve.good.GetN();
00137      printf(" capID=%d chi2=%.1f DOF=%d, chi/DOF=%f\n",cap,sum, bprsEve.good.GetN(),chi2dof);
00138      if(minChi2dof<=chi2dof) continue;
00139      //.... found new minimum
00140      minChi2dof=chi2dof;
00141      bestCapID=cap;
00142    }// end of loop over caps
00143 
00144    //  float minChi2dof=bprsEve.chi2dof-cut_stepChi2dof;
00145 
00146    if(bestCapID>=0)hChiGap->Fill(log(bprsEve.chi2dof-minChi2dof));
00147 
00148    if( minChi2dof<bprsEve.chi2dof-cut_stepChi2dof) { // better minimum was found, event corrupted
00149      assert(bestCapID>=0);
00150      // printf("ccc %f %f %f\n",minChi2dof,bprsEve.chi2dof ,log(bprsEve.chi2dof-minChi2dof));
00151      hChiGap2->Fill(log(bprsEve.chi2dof-minChi2dof));
00152      hCh2D->Fill(log(minChi2dof),log(bprsEve.chi2dof));
00153      hCapCorr->Fill(bprsEve.capID,bprsEve.crateID);
00154      printf("BETTER capID=%d\n",bestCapID);
00155      bprsEve.bestCapID=bestCapID;
00156      bprsEve.bestChi2dof=minChi2dof;
00157      bprsEve.useFix=true;
00158    }  else { // caps are in cync
00159      hCapGood->Fill(bprsEve.capID,bprsEve.crateID);
00160    } 
00161 }
00162 
00163 //________________________________________________
00164 //________________________________________________
00165 void BprsCapPolygraph::doPedResidua(JanBprsEveA & bprsEve){
00166    printf("doPedRes: eveID=%d: %d good\n", bprsEve.eveID,bprsEve.good.GetN());
00167  
00168    int capID=bprsEve.capID;
00169    int ibp=kBPrs;
00170 
00171    switch(par_pedFlag){
00172    case 0: break; // all w/o correction
00173    case 1: if(bprsEve.bestCapID>=0) return; // skip corrupted events
00174      break;
00175    case 2: if(bprsEve.bestCapID<0) return; // skip good event
00176      capID=bprsEve.bestCapID; // fix capID
00177      break;
00178    case 3: if(bprsEve.bestCapID>=0) capID=bprsEve.bestCapID; // fix & keep all
00179      break;
00180    }
00181 
00182    int nX=0;//tmp
00183 
00184    for(int i=0;i< bprsEve.raw.GetN();i++) {
00185      double rawAdc, yid;
00186      bprsEve.raw.GetPoint(i,rawAdc,yid);
00187      int id= (int) yid;
00188      // printf("cK i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
00189      int   stat =mJanDbMaker->statTile(ibp,id);
00190      if(stat) continue;// drop bad channels      
00191      float ped  =mJanDbMaker->pedTile(ibp,id,capID);
00192      double del=rawAdc-ped;
00193      //    hPedRes2D->Fill(id,del);
00194      if(bprsEve.bestCapID<0) hAdcGood->Fill(del); // good
00195      else  hAdcCorr->Fill(del);// fixed
00196 
00197      if(del<-10) nX++;
00198      } // end of data loop 
00199    printf("nX=%d\n",nX);
00200  }
00201