StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
calibBprsBox.C
1 // cat ps/bprs* |ps2pdf - BPM12.pmt1x.pdf
2 //cat psV6/soft* | ps2pdf - allTower.pdf
3 //cat psV6/bprsBPM* | ps2pdf - allBPM.pdf
4 
5 #include <vector>
6 const int mxBTile=2;
7 
8 vector <int> softL; // list of Bprs softID's belonging to given PMT
9 TH2F *h2DoutRaw=0; // raw 2D spectrum
10 TH2F *h2DoutMip=0; // mip 2D spectrum
11 TH1F *h1MaP=0,*h1MaT=0; // MIP ADC
12 TH2F *h2Cr=0; // MIP ADC vs. slopw
13 
14 TH2F* h2DinRawP, *h2DinTrP, *h2DinMipP;
15 TH2F* h2DinRawT, *h2DinTrT, *h2DinMipT;
16 TH1F* janDb_mipMean[mxBTile], * janDb_mipSig[mxBTile], * janDb_mipStat[mxBTile];
17 int doPS=2; // postscript : 1=per MPT, 2=per softID
18 TFile *fout=0;
19 
20 calibBprsBox(int box=11, int pmt=5) {
21  gStyle->SetOptStat(1110);
22  gStyle->SetOptFit(1);
23  gStyle->SetPalette(1,0);
24 
25 
26  char *fname="calib-nov-21-2008/sumD43-70V6.hist.root";
27  fd=new TFile(fname); assert(fd->IsOpen());
28  printf("Work with %s\n", fd->GetName()); //fd->ls();
29 
30  char *fnameO="barrelMipSpectV6.hist.root";
31  fout=new TFile(fnameO,"update"); assert(fout->IsOpen());
32  printf("Write output to %s\n", fout->GetName());
33 
34  //... load used calibration ....
35  char txt[100];
36  sprintf(txt,"calib-nov-20-2008/mipGainBprs+Btow_v2.hist.root");
37  TFile* fd5=new TFile(txt); assert(fd5->IsOpen());
38  fd5->ls();
39  char *core2[mxBTile]={"btow","bprs"};
40  for(int ibp=0;ibp<mxBTile;ibp++) {
41  TString tit=core2[ibp]; tit+="MipGain";
42  janDb_mipMean[ibp]=(TH1F *)fd5->Get(tit); assert( janDb_mipMean[ibp]);
43  tit=core2[ibp]; tit+="MipSig";
44  janDb_mipSig[ibp]=(TH1F *)fd5->Get(tit); assert( janDb_mipSig[ibp]);
45  tit=core2[ibp]; tit+="MipStat";
46  janDb_mipStat[ibp]=(TH1F *)fd5->Get(tit); assert( janDb_mipStat[ibp]);
47  }
48 
49 
50  h2DinRawP=(TH2F*)fd->Get("BPRS_c0"); assert(h2DinRawP);
51  h2DinTrP=(TH2F*)fd->Get("mipBprsTr"); assert(h2DinTrP);
52  h2DinMipP=(TH2F*)fd->Get("mipBprsTrBt"); assert(h2DinMipP);
53 
54  h2DinRawT=(TH2F*)fd->Get("BTOW_c0"); assert(h2DinRawT);
55  h2DinTrT=(TH2F*)fd->Get("mipBtowTr"); assert(h2DinTrT);
56  h2DinMipT=(TH2F*)fd->Get("mipBtowTrPr"); assert(h2DinMipT);
57 
58  buildSoftList(box,pmt);
59 
60  char core[1000];
61  sprintf(core, "BPRS PMB=%d pmt=%d",box,pmt);
62  int page=0;
63  can=new TCanvas("aa","aa",800,820);
64  TPad *c=makeTitle(can,core,page);
65  c->Divide(1,3);
66 
67  axY=h2DinRawP->GetYaxis();
68  float y1=axY->GetXmin();
69  float y2=axY->GetXmax();
70  int nbY=axY->GetNbins();
71  printf("Y-axis range --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
72 
73  char tt1[100], tt2[1000];
74  sprintf(tt1, "rawBPM%d_%d",box,pmt);
75  h2DoutRaw=new TH2F(tt1,"aa2;;ADC-ped",16,0.5,16.5,nbY,y1,y2);
76  h2DoutRaw->GetXaxis()->SetTitleSize(0.05);
77  sprintf(tt1, "mipBPM%d_%d",box,pmt);
78  h2DoutMip=new TH2F(tt1,"aa1;;ADC-ped",16,0.5,16.5,nbY,y1,y2);
79  h2DoutMip->GetXaxis()->SetTitleSize(0.05);
80 
81  sprintf(tt1, "gainBPM%d_%d",box,pmt);
82  sprintf(tt2, "BPRS gain BPM-%d_%d, all tiles; average MIP (ADC)",box,pmt);
83  h1MaP=new TH1F(tt1,tt2,40,0,40);
84  h1MaT=new TH1F("aa4","BTOW tower: mean MIP; ADC-ped",25,0,50);
85  h2Cr=new TH2F("aa5","BPRS comparison (both axis); slope (raw) ; average MIP (fit)",20,-0.15,-0.015,20,0,40);
86 
87  TString xtit=core; xtit+="; softID:";
88  for(int i=0;i<softL.size();i++) {
89  int softID=softL[i];
90  xtit+=softID; xtit+=" , ";
91  processOne(softID,core,i+1);
92  // break;
93  //if(i>1) break;
94  }
95  // return;
96  ln00=new TLine(0,0,17,0); ln00->SetLineStyle(2);
97  TString tit;
98  c->cd(3);
99  h2DoutRaw->Draw("colz"); gPad->SetLogz(); h2DoutRaw->SetAxisRange(-10,60,"y");
100  tit=core; tit+=", all values"+xtit;
101  h2DoutRaw->SetTitle(tit);
102  ln00->Draw();
103 
104  c->cd(2);
105  h2DoutMip->Draw("colz");
106  tit=core; tit+=", reqire MIP @ TPC & BTOW"+xtit;
107  h2DoutMip->SetTitle(tit);
108  h2DoutMip->Rebin2D(1,3);
109  ln00->Draw();
110  h2DoutMip->SetAxisRange(-10,60,"y");
111 
112  c->cd(1);
113  pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
114  pad->Draw();
115  pad->Divide(3,1);
116  pad->cd(1);
117  h1MaP->Draw();
118 
119  pad->cd(2);
120  h2Cr->Draw("box");
121 
122  pad->cd(3);
123  h1MaT->Draw(); // h1MaT->SetLineColor(kGreen);
124 
125  if(doPS) {
126  sprintf(core, "ps/bprsPmtPix%02d.ps",0);
127  if(doPS==2) sprintf(core, "ps/bprsBPM%03d_%d.ps", box , pmt);
128  can->Print(core);
129  }
130 
131  //....... write summary histos to output file
132  fout->cd();
133  h2DoutRaw->Write();
134  h2DoutMip->Write();
135  h1MaP->Write();
136  h2Cr->Write();
137  // fout->ls();
138 
139 }
140 
141 //------------------------
142 //------------------------
143 //------------------------
144 void processOne(int softID, char *core0,int pix) {
145  TString ttx=core0; ttx+=pix;
146  can=new TCanvas(ttx,ttx,800,700);
147  char core[1000];
148  sprintf(core, "%s, softID=%d",core0,softID);
149 
150  c=makeTitle(can,core,pix);
151  c->Divide(2,2);
152  c->cd(1);
153 
154  TH1F * hrP= getSlice( h2DinRawP,softID,"bprsR");
155  TH1F * htP= getSlice( h2DinTrP,softID,"bprsT");
156  TH1F * hmP= getSlice( h2DinMipP,softID,"bprsM");
157 
158  TH1F * hrT= getSlice( h2DinRawT,softID,"btowR");
159  TH1F * htT= getSlice( h2DinTrT,softID,"btowT");
160  TH1F * hmT= getSlice( h2DinMipT,softID,"btowM");
161 
162  hrP->SetTitle("BPRS spectrum, all events");
163  hrP->Fit("expo","","R",13,50);
164  TF1 *ff=hrP->GetFunction("expo");
165  float slope=ff->GetParameter(1);
166  hrP->Fit("gaus","+","R",-10,10);
167 
168  hmP->SetLineColor(kRed); hmP->SetFillColor(kRed);
169  hmP->SetTitle("BPRS spectrum, MIP : TPC, TPC+BTOW");
170  hmP->Fit("gaus","","R",-10,50);
171 
172  TF1 *ff=hmP->GetFunction("gaus");
173  ff->SetLineColor(kBlue); ff->SetLineWidth(1);
174  float mean=ff->GetParameter(1);
175  float meanS=ff->GetParameter(2);
176  float cut_adcL=2;
177  float cut_sig=1.8;
178  if(mean>cut_adcL && meanS>cut_sig) { // cut if MIP peak is too close to ped
179  h1MaP->Fill(mean);
180  h2Cr->Fill(slope,mean);
181  }
182 
183  hrT->SetTitle("BTOW spectrum, all events");
184  hrT->Fit("expo","","R",10,40);
185  hmT->SetLineColor(kGreen);hmT->SetFillColor(kGreen);
186  hmT->Fit("gaus","","R",-10,60);
187  TF1 *ffT=hmT->GetFunction("gaus");
188  ffT->SetLineColor(kBlue); ffT->SetLineWidth(1);
189  float meanT=ffT->GetParameter(1);
190  float meanST=ffT->GetParameter(2);
191  h1MaT->Fill(meanT);
192 
193  hmT->SetTitle("BTOW spectrum, MIP : TPC, TPC+BPRS");
194 
195  //..... BPRS .....
196  c->cd(1);
197  hrP->Draw(); gPad->SetLogy(); hrP->SetMinimum(0.7);
198  hmP->Draw("same"); gPad->SetGrid();
199  float xC=mean-meanS,yC=1e7;
200  ln=new TLine(xC,0,xC,yC); ln->SetLineColor(kMagenta); ln->Draw();
201 
202  c->cd(3);
203  hmP->Draw(); hmP->SetMaximum(htP->GetMaximum()/3);
204  htP->Draw("same");
205 
206  //..... BTOW.....
207  c->cd(2);
208  hrT->Draw(); gPad->SetLogy(); hrT->SetMinimum(0.7);
209  hmT->Draw("same");gPad->SetGrid();
210 
211  c->cd(4);
212  hmT->Draw(); hmT->SetMaximum(htT->GetMaximum()/3);
213  htT->Draw("same");
214 
215  ln=new TLine(10,-3,25,-3); ln->SetLineColor(kMagenta); ln->Draw();ln->SetLineWidth(2);
216  int nb=hmP->GetXaxis()->GetNbins();
217  for(int b=1;b<=nb;b++) {
218  h2DoutRaw->SetBinContent(pix,b,hrP->GetBinContent(b));
219  h2DoutMip->SetBinContent(pix,b,hmP->GetBinContent(b));
220  }
221 
222  FILE *fcsv=fopen("bprs+btowMipGain.csv","a");
223  // fprintf(fcsv,"average MIP ADC for BPRS tiles and BTOW towers; ver=1.6\n");
224  //fprintf(fcsv,"softID,PMB_pmt, pmt pix #, , bprs nMIP, bprs avrMIP (adc), bprs sigMIP (adc), , btow nMIP, btow avrMIP (adc), btow sigMIP (adc)\n");
225  fprintf(fcsv,"%d, %s, %d, , %d, %.2f, %.2f, , %d, %.2f, %.2f\n",
226  softID,core0+9,pix,(int)hmP->GetEntries(),mean,meanS,
227  (int)hmT->GetEntries(),meanT,meanST );
228  fclose(fcsv);
229  //..... compute ADC limmits
230  for(int ibp=0;ibp<mxBTile;ibp++) {
231  float adcL=5., adcH=40.; // default for _BTOW_
232  if(ibp==1) adcL=3.5; // BPRS
233  int statGain=(int)janDb_mipStat[ibp]->GetBinContent(softID);
234  if(!statGain){ //valid gain
235  float mean=janDb_mipMean[ibp]->GetBinContent(softID);
236  float sig=janDb_mipSig[ibp]->GetBinContent(softID);
237  adcL=mean-sig;
238  adcH=mean+sig;
239  if(ibp==0 && adcL<5) adcL=5; // BTOW
240  if(ibp==1 && adcL<3.5) adcL=3.5;// BPRS
241  if(adcH>2*mean) adcH=2*mean;
242  printf("ibp=%d id=%4d MIP mean=%.1f sig=%.1f adcL=%.1f H=%.1f\n",ibp,softID,mean,sig,adcL,adcH);
243  }
244  float y=1e5;
245  TH1F *h1=hmT; if (ibp==1) h1=hmP;
246  Lx=h1->GetListOfFunctions(); assert(Lx);
247  lnL=new TLine( adcL,0, adcL,y); lnL->SetLineColor(kBlue);lnL->SetLineWidth(2);lnL->SetLineStyle(2);
248  Lx->Add(lnL);
249  lnH=new TLine( adcH,0, adcH,y); lnH->SetLineColor(kBlue);lnH->SetLineWidth(2);lnH->SetLineStyle(2);
250  Lx->Add(lnH);
251  } // Loop over B/P
252 
253  //....... write 1D histos to output file
254  fout->cd();
255  char tt1[100], tt2[1000];
256  sprintf(tt2, "id=%d %s pix=%d, ",softID,core0,pix);
257  sprintf(tt1, "bprs%d",softID);
258 
259  TString Tt=tt1, T2=tt2;
260  hrP->SetName(Tt+"a"); hrP->SetTitle(T2+"inclusive"); hrP->Write();
261  htP->SetName(Tt+"t"); htP->SetTitle(T2+"MIP in TPC"); htP->Write();
262  hmP->SetName(Tt+"m"); hmP->SetTitle(T2+"MIP in TPC+BTOW");hmP->Write();
263 
264  sprintf(tt2, "id=%d BTOW , ",softID);
265  sprintf(tt1, "btow%d",softID);
266  Tt=tt1, T2=tt2;
267  hrT->SetName(Tt+"a"); hrT->SetTitle(T2+"inclusive"); hrT->Write();
268  htT->SetName(Tt+"t"); htT->SetTitle(T2+"MIP in TPC"); htT->Write();
269  hmT->SetName(Tt+"m"); hmT->SetTitle(T2+"MIP in TPC+BPRS");hmT->Write();
270  //...... FINISH .....
271  sprintf(core, "softID=%d",softID);
272  can->SetTitle(core);
273  if(doPS) {
274  sprintf(core, "ps/bprsPmtPix%02d.ps",pix);
275  if(doPS==2) sprintf(core, "ps/softID%04d.ps",softID);
276  can->Print(core);
277  }
278 
279 
280  }
281 
282 
283 //=================
284 TH1F * getSlice(TH2F * h2, int id, char *ctit) {
285  axX=h2->GetXaxis();
286  float x1=axX->GetXmin();
287  float x2=axX->GetXmax();
288  int nbX=axX->GetNbins();
289  printf("X-axis range --> [%.1f, %.1f], nb=%d %s\n",x1,x2,nbX,axX->GetTitle());
290 
291  axY=h2->GetYaxis();
292  float y1=axY->GetXmin();
293  float y2=axY->GetXmax();
294  int nbY=axY->GetNbins();
295  printf("Y-axis range --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
296 
297  assert(id>=1 && id<=nbX);
298  // do projections
299  char txt1[100], txt2[1000];
300  sprintf(txt1,"%s_id%d",ctit,id);
301  sprintf(txt2,"%s soft id=%d;%s ",ctit, id,axY->GetTitle());
302  TH1F*h=new TH1F(txt1,txt2,nbY,y1,y2); // working histo for 1-D spectrum
303 
304  int i;
305  for(i=1;i<=nbY;i++) h->SetBinContent(i,h2->GetBinContent(id,i));
306  float x1=-20, x2=70;
307  h->SetAxisRange(x1,x2);
308  h->SetEntries(h->Integral());
309 
310  return h;
311 }
312 
313 
314 //------------------------
315 void buildSoftList(int box=11, int pmt=1) {
316  fd1=new TFile("calib-nov-8-2008/map-softID-bprsPmt-Rory.root"); assert(fd1->IsOpen());
317  printf("Opened: %s\n",fd1->GetName()); fd1->ls();
318  TH1I * mapBprsPmt=(TH1I*) fd1->Get("bprsPmt"); assert(mapBprsPmt);
319  //mapBprsPmt->Draw();
320 
321  int y= box*10+pmt;
322  for(int b=1;b<=4800;b++) {
323  if(y!=(int)mapBprsPmt->GetBinContent(b)) continue;
324  softL.push_back(b);
325  }
326  printf("found pmb=%d pmt=%d Bprs tiles =%d\n",box,pmt,softL.size());
327 }
328 
329 //------------------------
330 TPad *makeTitle(TCanvas *c,char *core, int page) {
331 
332  c->Range(0,0,1,1);
333  TPad *pad0 = new TPad("pad0", "apd0",0.0,0.95,1.,1.);
334  pad0->Draw();
335  pad0->cd();
336 
337  TPaveText *pt = new TPaveText(0,0.,1,1,"br");
338  pt->Draw();
339  TDatime dt;
340  TString txt2=core;
341  txt2+=", pixel=";
342  txt2+=page;
343  txt2+=", ";
344  txt2+=dt.AsString();
345  pt->AddText(txt2);
346  txt2="--";
347  pt->AddText(txt2);
348 
349  c->cd();
350  pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
351  pad->Draw();
352  return pad;
353 }
354 
Generated by   doxygen 1.8.5