StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BprsCapPolygraph.cxx
1 #include <stdio.h>
2 #include <TH2F.h>
3 #include <TObjArray.h>
4 #include <TGraph.h>
5 
6 #include "BprsCapPolygraph.h"
7 #include "StJanBarrelDbMaker.h"
8 #include "JanBarrelEvent.h"
9 
10 
11 //________________________________________________
12 //________________________________________________
13 BprsCapPolygraph::BprsCapPolygraph( TObjArray *HList, StJanBarrelDbMaker* jdb, int pedFlag) {
14  par_pedFlag=pedFlag;
15  mJanDbMaker=jdb; assert(mJanDbMaker);
16 
17  const char *fname="all"; //tmp
18  //QA histos
19 
20  hChiB=new TH1F("bppg_chb","BprsPoly logN(chi2/dof), nominal capID ; logN(chi2/dof)",200,-1.,9.);
21  hChiGap=new TH1F("bppg_chgI","BprsPoly chi2/dof (nominal - best_capID), any minimum; logN( nominal_chi2 - best_chi2 ) ",200,-6,10);
22  hChiGap2=new TH1F("bppg_chgA","BprsPoly chi2/dof (nominal - best_capID), acceted; logN( nominal_chi2 - best_chi2 ) ",200,-6,10);
23 
24  char tt2[1000];
25  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);
26  sprintf(tt2,"BprsPoly caps in sync, %s; nominal capID; BPRS crate ID",fname);
27  hCapGood=new TH2F("bppg_capSyn",tt2, 128,-0.5,127.5,mxBprsCrate,-0.5,mxBprsCrate-0.5);
28  sprintf(tt2,"BprsPoly caps NOT sync, %s; nominal capID; BPRS crateID",fname);
29  hCapCorr=new TH2F("bppg_capDeSyn",tt2, 128,-0.5,127.5,mxBprsCrate,-0.5,mxBprsCrate-0.5);
30 
31  HList->Add( hChiB);
32  HList->Add( hChiGap);
33  HList->Add( hChiGap2);
34  HList->Add( hCh2D);
35  HList->Add(hCapGood);
36  HList->Add(hCapCorr);
37 
38  int nb=200;
39  float adc1=-50;
40  float adc2=adc1+nb;
41 
42  hAdcGood=new TH1F("bppg_adcGd","BprsPoly ADC for capID in sync; rawAdc - ped(nominal capID)",nb,adc1,adc2);
43  hAdcCorr=new TH1F("bppg_adcCor","BprsPoly ADC for capID NOT in sync,fixed; rawAdc - ped(best capID)",nb,adc1,adc2);
44  hAdcCorr->SetLineColor(kRed);
45 
46  HList->Add(hAdcGood);
47  HList->Add( hAdcCorr);
48 
49  }
50 
51 //________________________________________________
52 //________________________________________________
53 void BprsCapPolygraph::doBaseline( JanBprsEveA & bprsEve, JanBarrelEvent &fullEve){
54 
55  int crateID=bprsEve.crateID;
56  int capID=bprsEve.capID;
57  printf("doBaseline: bprsEve: %d rawADC, doCrate=%d nominalCapID=%d\n", bprsEve.raw.GetN(),crateID,capID);
58  TGraph res; //residua: X=ADC-ped, y=index of raw;
59 
60  int ibp=kBPrs;
61  for(int i=0;i< bprsEve.raw.GetN();i++) {
62  double rawAdc, yid;
63  bprsEve.raw.GetPoint(i,rawAdc,yid);
64  int id= (int) yid;
65  // printf("K i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
66  int stat =mJanDbMaker->statTile(ibp,id);
67  if(stat) continue;
68  float ped =mJanDbMaker->pedTile(ibp,id,capID);
69  float sig =mJanDbMaker->sigPedTile(ibp,id,capID);
70  double del=rawAdc-ped;
71  if(del> cut_adcMax) continue;
72  // if(i<20)
73  //printf("i=%d id=%d del=%f sig=%f\n",i,id,del,sig);
74  assert(sig!=0);
75  del/=sig;
76  res.SetPoint(res.GetN(),del,i);
77  } // end of loop over raw data
78 
79  res.Sort();
80  int i2=(int)(res.GetN()*cut_fracSkip);
81  int i1=res.GetN()-i2;
82  printf("res size=%d , use=[%d,%d]\n", res.GetN(),i1,i2);
83  // res.Print();
84 
85  double sum=0;
86  for(int i=i1; i<i2;i++) {
87  double rdel, yj;
88  res.GetPoint(i,rdel,yj);
89  sum+=rdel*rdel; // compute base chi2/dof
90  //if(i<20) printf("CC i=%d rdel=%f sig=%f sum=%f\n",i,rdel,sig,sum);
91  // prepare array of index for chi2 minimalization
92  int j= (int) yj; // index in 'raw' array
93  double rawAdc, yid;
94  bprsEve.raw.GetPoint(j,rawAdc,yid);
95  bprsEve.addGoodValue(yid,rawAdc);
96  }// base chi2 computed
97  int dof=i2-i1;
98  float chi2dof=sum/dof;
99  printf(" base chi2=%f DOF=%d, chi/DOF=%f nGood=%d\n",sum,dof,chi2dof,bprsEve.good.GetN());
100  bprsEve.chi2=sum;
101  bprsEve.chi2dof=chi2dof;
102  hChiB->Fill(log(chi2dof));
103 }
104 
105 //________________________________________________
106 //________________________________________________
107 void BprsCapPolygraph::findBestCap(JanBprsEveA & bprsEve, JanBarrelEvent &fullEve){
108 
109  printf("find BPRS bestCap: eve: %d good\n", bprsEve.good.GetN());
110  int ibp=kBPrs;
111 
112  float minChi2dof=bprsEve.chi2dof-0.0001;
113  int bestCapID=-1;
114 
115  int cap1=bprsEve.capID-par_mxDelCap+mxBcap;
116  int cap2=bprsEve.capID+par_mxDelCap+mxBcap;
117  for(int ic=cap1; ic<=cap2;ic++) {
118  int cap=ic%mxBcap;
119  // printf("xx cap=%d\n",cap);
120  double sum=0;
121  for(int i=0;i< bprsEve.good.GetN();i++) {
122  double rawAdc, yid;
123  bprsEve.good.GetPoint(i,rawAdc,yid);
124  int id= (int) yid;
125  // printf("bK i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
126  int stat =mJanDbMaker->statTile(ibp,id);
127 
128  assert(stat==0);// bad channels were already discarded
129  float ped =mJanDbMaker->pedTile(ibp,id,cap);
130  float sig =mJanDbMaker->sigPedTile(ibp,id,cap);
131  double del=rawAdc-ped;
132  //if(i<20) printf("i=%d id=%d del=%f sig=%f\n",i,id,del,sig);
133  del/=sig;
134  sum+=del*del;
135  } // end of data loop
136  float chi2dof=sum/bprsEve.good.GetN();
137  printf(" capID=%d chi2=%.1f DOF=%d, chi/DOF=%f\n",cap,sum, bprsEve.good.GetN(),chi2dof);
138  if(minChi2dof<=chi2dof) continue;
139  //.... found new minimum
140  minChi2dof=chi2dof;
141  bestCapID=cap;
142  }// end of loop over caps
143 
144  // float minChi2dof=bprsEve.chi2dof-cut_stepChi2dof;
145 
146  if(bestCapID>=0)hChiGap->Fill(log(bprsEve.chi2dof-minChi2dof));
147 
148  if( minChi2dof<bprsEve.chi2dof-cut_stepChi2dof) { // better minimum was found, event corrupted
149  assert(bestCapID>=0);
150  // printf("ccc %f %f %f\n",minChi2dof,bprsEve.chi2dof ,log(bprsEve.chi2dof-minChi2dof));
151  hChiGap2->Fill(log(bprsEve.chi2dof-minChi2dof));
152  hCh2D->Fill(log(minChi2dof),log(bprsEve.chi2dof));
153  hCapCorr->Fill(bprsEve.capID,bprsEve.crateID);
154  printf("BETTER capID=%d\n",bestCapID);
155  bprsEve.bestCapID=bestCapID;
156  bprsEve.bestChi2dof=minChi2dof;
157  bprsEve.useFix=true;
158  } else { // caps are in cync
159  hCapGood->Fill(bprsEve.capID,bprsEve.crateID);
160  }
161 }
162 
163 //________________________________________________
164 //________________________________________________
165 void BprsCapPolygraph::doPedResidua(JanBprsEveA & bprsEve){
166  printf("doPedRes: eveID=%d: %d good\n", bprsEve.eveID,bprsEve.good.GetN());
167 
168  int capID=bprsEve.capID;
169  int ibp=kBPrs;
170 
171  switch(par_pedFlag){
172  case 0: break; // all w/o correction
173  case 1: if(bprsEve.bestCapID>=0) return; // skip corrupted events
174  break;
175  case 2: if(bprsEve.bestCapID<0) return; // skip good event
176  capID=bprsEve.bestCapID; // fix capID
177  break;
178  case 3: if(bprsEve.bestCapID>=0) capID=bprsEve.bestCapID; // fix & keep all
179  break;
180  }
181 
182  int nX=0;//tmp
183 
184  for(int i=0;i< bprsEve.raw.GetN();i++) {
185  double rawAdc, yid;
186  bprsEve.raw.GetPoint(i,rawAdc,yid);
187  int id= (int) yid;
188  // printf("cK i=%d id=%d rawAdc=%f\n",i,id,rawAdc);
189  int stat =mJanDbMaker->statTile(ibp,id);
190  if(stat) continue;// drop bad channels
191  float ped =mJanDbMaker->pedTile(ibp,id,capID);
192  double del=rawAdc-ped;
193  // hPedRes2D->Fill(id,del);
194  if(bprsEve.bestCapID<0) hAdcGood->Fill(del); // good
195  else hAdcCorr->Fill(del);// fixed
196 
197  if(del<-10) nX++;
198  } // end of data loop
199  printf("nX=%d\n",nX);
200  }
201 
Generated by   doxygen 1.8.5