StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcDAQ2Ped.cxx
1 #include <TF1.h>
2 
3 #include "StEEmcDAQ2Ped.h"
4 #include "StEEmcUtil/database/StEEmcDb.h"
5 #include "StEEmcUtil/database/EEmcDbItem.h"
6 
7 ClassImp(StEEmcDAQ2Ped)
8 
9 //---------------------------
10 //---------------------------
11 //---------------------------
12 StEEmcDAQ2Ped::StEEmcDAQ2Ped(const char* name, TFile* file)
13  : StMaker(name) {
14  mHList=0;
15  mDAQHistos=file;
16  cout <<"StEEmcDAQ2Ped::StEEmcDAQ2Ped()"<<endl;
17 }
18 
19 //---------------------------
20 //---------------------------
21 //---------------------------
22 
23 StEEmcDAQ2Ped::~StEEmcDAQ2Ped()
24 {
25  cout <<"StEEmcDAQ2Ped::~StEEmcDAQ2Ped()"<<endl;
26 }
27 
28 //---------------------------
29 //---------------------------
30 //---------------------------
31 
32 Int_t
33 StEEmcDAQ2Ped::Init() {
34  mEeDb = (StEEmcDb*)this->GetDataSet("StEEmcDb");
35  assert(mEeDb); // eemcDB must be in the chain, fix it
36  return kStOk;
37 }
38 
39 //------------------------------------
40 //------------------------------------
41 
42 Int_t StEEmcDAQ2Ped::InitRun(int runNo){
43 
44  initHistos(); //initialize esmd histos
45 
46  assert(mEeDb);
47  assert(mDAQHistos);
48 
49  // ................... histo for individual pixels ...
50  int icrate,ichan=0;
51 
52  // loop over tower crates
53  for(icrate=0;icrate<MaxTwCrateID;icrate++) {
54  char tt0[100];
55  sprintf(tt0,"ETOW_%d",icrate+1);
56  h2D=(TH2F*)mDAQHistos->Get(tt0); assert(h2D);
57  for (ichan=0;ichan<=MaxTwCrateCh;ichan++) {
58  const EEmcDbItem *x=mEeDb->getByCrate(icrate+1,ichan);
59  if(x==0) continue;
60  // initialize histos only for pixels acquired from DB
61  char tt1[100],tt2[200];
62  sprintf(tt1,"a%s",x->name);
63  sprintf(tt2,"ADC for %s, cr/chan=%3.3d/%3.3d, tube=%s; ADC",x->name,x->crate,x->chan,x->tube);
64  //cout<<tt2<<endl;
65  //x->print();
66  TH1F* h=new TH1F(tt1,tt2,5000,-20.5,4979.5);
67  // Fill 1D histos
68  TAxis* axY=h2D->GetYaxis();
69  int nbY=axY->GetNbins();
70  //if(ichan == 0)
71  // cout<<tt2<<" crate/chan "<<icrate<<"/"<<ichan<<endl;
72  for(int i=1;i<=nbY;i++) //weight entries in 1D histos
73  h->SetBinContent(i+20,h2D->GetBinContent(x->chan+1,i));
74  mHList->Add(h);
75 
76  //now find slope and plot in 2D histo
77  if(x->fail) continue;
78  float entries=h2D->GetEntries()/160;
79  float xInt1=0; float xInt2=0; int xbin=-1;
80  int xHigh = -1; int xLow = -1;
81  for(xbin=h2D->GetNbinsY();xbin>0;xbin--){
82  xInt1 += h2D->GetBinContent(ichan+1,xbin);
83  //cout<<"xbin="<<xbin<<" xInt1="<<xInt1<<endl;
84  if(xInt1 > 0.00036*entries){
85  xHigh = xbin; break;}
86  }
87  for(xbin=h2D->GetNbinsY();xbin>0;xbin--){
88  xInt2 += h2D->GetBinContent(ichan+1,xbin);
89  if(xInt2 > 0.0027*entries){
90  xLow = xbin; break;}
91  }
92  float ped=x->ped; float sigPed=x->sigPed;
93  if(xLow < (ped+(5*sigPed)))
94  cout<<"ichan= "<<ichan<<" x->name= "<<x->name<<" x->tube= "<<x->tube<<endl;
95  xLowEtow->Fill(xLow-ped);
96  xHighEtow->Fill(xHigh-ped);
97  xDiffEtow->Fill(xHigh-xLow);
98  xCorrEtow->Fill(xLow-ped,xHigh-ped);
99  TF1* f1=new TF1("f1","exp([0]+[1]*x)",xLow,xHigh);
100  f1->SetLineColor(2);
101  h->Fit(f1,"Q","",xLow,xHigh);
102  float slope = f1->GetParameter(1);
103  float slopeErr = f1->GetParError(1);
104  //if(slope<0 && slope>-0.1){ //don't plot outliers
105  etow[icrate]->SetBinContent(ichan+1,slope);
106  etow[icrate]->SetBinError(ichan+1,slopeErr);
107 
108  //if(icrate==0) cout<<"ichan= "<<ichan<<" x->name= "<<x->name<<" x->tube= "<<x->tube<<endl;
109  }
110  }
111 
112  //ascii output file for mapping
113  ofstream outfile(mappingFile.Data());
114 
115  // loop over Mapmt crates
116  for(icrate=MinMapmtCrateID;icrate<=MaxMapmtCrateID;icrate++) {
117  char tt0[100];
118  sprintf(tt0,"ESMD_%d",icrate-63);//daq counts crates from 1
119  h2D=(TH2F*)mDAQHistos->Get(tt0); assert(h2D);
120  for (ichan=0;ichan<MaxMapmtCrateCh;ichan++) {
121  const EEmcDbItem *x=mEeDb->getByCrate(icrate,ichan);
122  if(x==0) continue;
123  // initialize histos only for pixels acquired from DB
124  char tt1[100],tt2[200];
125  sprintf(tt1,"a%s",x->name);
126  sprintf(tt2,"ADC for %s, cr/chan=%3.3d/%3.3d, tube=%s; ADC",x->name,x->crate,x->chan,x->tube);
127  TH1F* h=new TH1F(tt1,tt2,5000,-20.5,4979.5);
128  // Fill 1D histos
129  TAxis* axY=h2D->GetYaxis();
130  int nbY=axY->GetNbins();
131  for(int i=1;i<=nbY;i++) //weight entries in 1D histos
132  h->SetBinContent(i+20,h2D->GetBinContent(x->chan+1,i));
133  mHList->Add(h);
134 
135  outfile<<icrate<<","<<ichan<<","<<x->name<<endl;
136  //now find slope and plot in 2D histo
137  if(x->fail) continue;
138  float entries=h2D->GetEntries()/192;
139  float xInt1=0; float xInt2=0; int xbin=-1;
140  int xHigh = -1; int xLow = -1;
141  for(xbin=h2D->GetNbinsY();xbin>0;xbin--){
142  xInt1 += h2D->GetBinContent(ichan+1,xbin);
143  //cout<<"xbin="<<xbin<<" xInt1="<<xInt1<<endl;
144  if(xInt1 > 0.0003*entries){
145  xHigh = xbin; break;}
146  }
147  for(xbin=h2D->GetNbinsY();xbin>0;xbin--){
148  xInt2 += h2D->GetBinContent(ichan+1,xbin);
149  if(xInt2 > 0.0015*entries){
150  xLow = xbin; break;}
151  }
152  float ped=x->ped; //float sigPed=x->sigPed;
153  float stopper=ped; stopper +=25;
154  if(xLow < stopper){
155  xTestEsmd->Fill(xLow-ped-stopper);
156  //cout<<"ichan= "<<ichan<<" x->name= "<<x->name<<" [xLow,xHigh] ["<<xLow<<","<<xHigh<<"] ped+25 = "<<stopper;
157  xLow=(int)stopper+1;
158  //cout<<" new xLow "<<xLow<<endl;
159  //continue;
160  }
161  xDiffEsmd->Fill(xHigh-xLow);
162  if(xHigh - xLow < 10) {
163  //cout<<"Not enough channels fit ichan="<<ichan<<" x->name="<<x->name<<" # chan's fit="<<xHigh-xLow<<" # counts past xLow="<<h2D->Integral(ichan+1,ichan+1,xLow,1024)/entries<<endl;
164  continue;
165  }
166  xLowEsmd->Fill(xLow-ped);
167  xHighEsmd->Fill(xHigh-ped);
168  xCorrEsmd->Fill(xLow-ped,xHigh-ped);
169  TF1* f1=new TF1("f1","exp([0]+[1]*x)",xLow,xHigh);
170  f1->SetLineColor(2);
171  h->Fit(f1,"Q","",xLow,xHigh);
172  float slope = f1->GetParameter(1);
173  float slopeErr = f1->GetParError(1);
174  //fill histos by sector
175  int sec=x->sec; int plane=3; int strip=x->strip;
176  if(x->plane=='U') plane=0;
177  if(x->plane=='V') plane=1;
178  //printf("Sec %d plane %c=%d strip %d : name %s\n",sec,x->plane,plane,strip,x->name);
179  if(plane==0 || plane==1){
180  esmdSec[sec-1][plane]->SetBinContent(strip,slope);
181  esmdSec[sec-1][plane]->SetBinError(strip,slopeErr);
182  }
183  //fill histos by crate
184  if(plane!=0 && plane !=1){//eprs first
185  esmd[icrate-63]->SetBinContent(ichan+1,slope);
186  esmd[icrate-63]->SetBinError(ichan+1,slopeErr);
187  }
188  else if(strip<221 && strip>10){//then esmd for long strips
189  esmd[icrate-63]->SetBinContent(ichan+1,slope);
190  esmd[icrate-63]->SetBinError(ichan+1,slopeErr);
191  }
192  //end gain stuff
193 
194  }
195  }
196 
197  printf("StEEmcDAQ2Ped::InitRun() \n");
198 
199  return kStOk;
200 }
201 
202 //------------------------------------
203 //------------------------------------
204 void
205 StEEmcDAQ2Ped::initHistos()
206 {
207 
208  int icrate=0;
209  for(icrate=MinMapmtCrateID;icrate<=MaxMapmtCrateID;icrate++){
210  //need to initialize slope histos for each crate
211  char name[100]; char title[100];
212  sprintf(name,"ESMD_%d",icrate-63);
213  sprintf(title,"ESMD FEE %d",icrate);
214  esmd[icrate-63] = new TH1F(name,title,192,-0.5,192-0.5);
215  esmd[icrate-63]->SetMarkerStyle(20+mSet);//TCD25=21 and TCD63=22
216  esmd[icrate-63]->SetMarkerColor(mSet+1);//TCD25=2 and TCD63=3
217  mHList->Add(esmd[icrate-63]);
218  }
219  for(icrate=0;icrate<MaxTwCrateID;icrate++) {
220  //need to initialize slope histos for each crate
221  char nameT[100]; char titleT[100];
222  sprintf(nameT,"ETOW_%d",icrate+1);
223  sprintf(titleT,"ETOW FEE %d",icrate+1);
224  etow[icrate] = new TH1F(nameT,titleT,120,-0.5,119.5);
225  etow[icrate]->SetMarkerStyle(20+mSet);//TCD25=21 and TCD63=22
226  etow[icrate]->SetMarkerColor(mSet+1);//TCD25=2 and TCD63=3
227  mHList->Add(etow[icrate]);
228  }
229  for(int isector=0;isector<12;isector++){
230  for(int iuv=0;iuv<2;iuv++){
231  //need to initialize slope histos for each sector + plane
232  char name[100]; char title[100];
233  sprintf(name,"ESMD_sector%d_%cplane",isector+1,'U'+iuv);
234  sprintf(title,"ESMD Sector%d %c-Plane",isector+1,'U'+iuv);
235  esmdSec[isector][iuv] =new TH1F(name,title,288,0.5,288.5);
236  esmdSec[isector][iuv]->SetMarkerStyle(20+mSet);//TCD25=21 and TCD63=22
237  esmdSec[isector][iuv]->SetMarkerColor(mSet+1);//TCD25=2 and TCD63=3
238  mHList->Add(esmdSec[isector][iuv]);
239  }
240  }
241  xLowEtow = new TH1F("xLowEtow","Fit range minimum ETOW;xLow",60,0.0,60.0); mHList->Add(xLowEtow);
242  xHighEtow = new TH1F("xHighEtow","Fit range maximum ETOW;xHigh",100,0.0,100.0); mHList->Add(xHighEtow);
243  xDiffEtow = new TH1F("xDiffEtow","Fit range difference ETOW;xHigh-xLow",80,0.0,80.0); mHList->Add(xDiffEtow);
244  xLowEsmd = new TH1F("xLowEsmd","Fit range minimum ESMD;xLow",150,0.0,150.0); mHList->Add(xLowEsmd);
245  xHighEsmd = new TH1F("xHighEsmd","Fit range maximum ESMD;xHigh",300,0.0,300.0); mHList->Add(xHighEsmd);
246  xDiffEsmd = new TH1F("xDiffEsmd","Fit range difference ESMD;xHigh-xLow",200,0.0,200.0); mHList->Add(xDiffEsmd);
247  xCorrEtow = new TH2F("xCorrEtow","ETOW xHigh vs xLow;xLow;xHigh",60,0.0,60.0,100,0.0,100.0); mHList->Add(xCorrEtow);
248  xCorrEsmd = new TH2F("xCorrEsmd","ESMD xHigh vs xLow;xLow;xHigh",50,0.0,200.0,100,0.0,400.0); mHList->Add(xCorrEsmd);
249  xTestEsmd = new TH1F("xTestEsmd","ESMD xLow - (ped+5*sigPed)",60,-3.0,3.0); mHList->Add(xTestEsmd);
250 
251 }
252 
253 
254 
255 //------------------------------------
256 //------------------------------------
257 Int_t
258 StEEmcDAQ2Ped::Finish()
259 {
260 
261  TH1F* hSlope = new TH1F("hSlope","Average slope by crate;Crate;Slope",48,63.5,111.5);
262  hSlope->SetMarkerStyle(20+mSet);//TCD25=23 and TCD63=22
263  hSlope->SetMarkerColor(mSet+1);//TCD25=4 and TCD63=5
264  int icrate=0;
265  for(icrate=MinMapmtCrateID;icrate<=MaxMapmtCrateID;icrate++){
266  TF1* f=new TF1("f","[0]");
267  f->SetLineColor(5+mSet); //TCD25=4 and TCD63=5
268  if((icrate-63)%4 != 0) //SMD crates
269  esmd[icrate-63]->Fit(f,"Q","",0,192);
270  else {//PRS crates
271  esmd[icrate-63]->Fit(f,"Q","",60,192);
272  continue;
273  }
274  float constant = f->GetParameter(0);
275  float constantErr = f->GetParError(0);
276  hSlope->SetBinContent(icrate-63,constant);
277  hSlope->SetBinError(icrate-63,constantErr);
278  }
279  mHList->Add(hSlope);
280 
281  return kStOk;
282 }
283 
284 
285 //------------------------------------
286 //------------------------------------
287 Int_t
288 StEEmcDAQ2Ped::Make(){
289 
290  return kStOk;
291 }
292 
293 
EEmcDbItem::name
char name[StEEmcNameLen]
ASCII name of the channel, see Readme.
Definition: EEmcDbItem.h:20
EEmcDbItem::tube
char tube[StEEmcNameLen]
name of PMT or MAPMT pixel
Definition: EEmcDbItem.h:21
StEEmcDAQ2Ped::Make
virtual Int_t Make()
Definition: StEEmcDAQ2Ped.cxx:288
StEEmcDAQ2Ped::Finish
virtual Int_t Finish()
Definition: StEEmcDAQ2Ped.cxx:258
EEmcDbItem::chan
int chan
hardware channel
Definition: EEmcDbItem.h:28
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5