StRoot: StEmcPool/StBarrelCalib2008/macros/calibBprsBox.C Source File

00001 // cat ps/bprs* |ps2pdf - BPM12.pmt1x.pdf
00002 //cat psV6/soft* | ps2pdf - allTower.pdf 
00003 //cat psV6/bprsBPM* | ps2pdf - allBPM.pdf 
00004 
00005 #include <vector>
00006 const int mxBTile=2;
00007 
00008 vector <int> softL; // list of Bprs softID's belonging to given PMT
00009 TH2F *h2DoutRaw=0; // raw 2D spectrum
00010 TH2F *h2DoutMip=0; // mip 2D spectrum
00011 TH1F *h1MaP=0,*h1MaT=0; // MIP ADC
00012 TH2F *h2Cr=0; // MIP ADC vs. slopw
00013 
00014 TH2F* h2DinRawP, *h2DinTrP, *h2DinMipP;
00015 TH2F* h2DinRawT, *h2DinTrT, *h2DinMipT;
00016 TH1F* janDb_mipMean[mxBTile], * janDb_mipSig[mxBTile], * janDb_mipStat[mxBTile];
00017 int doPS=2; // postscript : 1=per MPT, 2=per softID
00018 TFile *fout=0;
00019 
00020 calibBprsBox(int box=11, int pmt=5) {
00021   gStyle->SetOptStat(1110);
00022   gStyle->SetOptFit(1);
00023   gStyle->SetPalette(1,0);
00024  
00025 
00026   char *fname="calib-nov-21-2008/sumD43-70V6.hist.root";
00027   fd=new TFile(fname); assert(fd->IsOpen());
00028   printf("Work with %s\n", fd->GetName()); //fd->ls();
00029  
00030   char *fnameO="barrelMipSpectV6.hist.root";
00031   fout=new TFile(fnameO,"update"); assert(fout->IsOpen());
00032   printf("Write output to %s\n", fout->GetName()); 
00033 
00034   //... load used calibration ....
00035   char txt[100];
00036   sprintf(txt,"calib-nov-20-2008/mipGainBprs+Btow_v2.hist.root");
00037   TFile* fd5=new TFile(txt);  assert(fd5->IsOpen());
00038   fd5->ls();   
00039   char *core2[mxBTile]={"btow","bprs"};
00040   for(int ibp=0;ibp<mxBTile;ibp++) {
00041     TString tit=core2[ibp]; tit+="MipGain";
00042     janDb_mipMean[ibp]=(TH1F *)fd5->Get(tit); assert(   janDb_mipMean[ibp]);
00043     tit=core2[ibp]; tit+="MipSig";
00044     janDb_mipSig[ibp]=(TH1F *)fd5->Get(tit); assert(   janDb_mipSig[ibp]);
00045     tit=core2[ibp]; tit+="MipStat";
00046     janDb_mipStat[ibp]=(TH1F *)fd5->Get(tit); assert(   janDb_mipStat[ibp]);
00047   }
00048  
00049 
00050   h2DinRawP=(TH2F*)fd->Get("BPRS_c0"); assert(h2DinRawP);
00051   h2DinTrP=(TH2F*)fd->Get("mipBprsTr"); assert(h2DinTrP);
00052   h2DinMipP=(TH2F*)fd->Get("mipBprsTrBt"); assert(h2DinMipP);
00053 
00054   h2DinRawT=(TH2F*)fd->Get("BTOW_c0"); assert(h2DinRawT);
00055   h2DinTrT=(TH2F*)fd->Get("mipBtowTr"); assert(h2DinTrT);
00056   h2DinMipT=(TH2F*)fd->Get("mipBtowTrPr"); assert(h2DinMipT);
00057 
00058   buildSoftList(box,pmt);
00059   
00060   char core[1000];
00061   sprintf(core, "BPRS PMB=%d pmt=%d",box,pmt);
00062   int page=0;
00063   can=new TCanvas("aa","aa",800,820);   
00064   TPad *c=makeTitle(can,core,page);
00065   c->Divide(1,3);
00066  
00067   axY=h2DinRawP->GetYaxis();
00068   float y1=axY->GetXmin();
00069   float y2=axY->GetXmax();
00070   int nbY=axY->GetNbins();
00071   printf("Y-axis range  --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
00072 
00073   char tt1[100], tt2[1000];
00074   sprintf(tt1, "rawBPM%d_%d",box,pmt);
00075   h2DoutRaw=new TH2F(tt1,"aa2;;ADC-ped",16,0.5,16.5,nbY,y1,y2);
00076   h2DoutRaw->GetXaxis()->SetTitleSize(0.05);
00077   sprintf(tt1, "mipBPM%d_%d",box,pmt);
00078   h2DoutMip=new TH2F(tt1,"aa1;;ADC-ped",16,0.5,16.5,nbY,y1,y2);
00079   h2DoutMip->GetXaxis()->SetTitleSize(0.05);
00080 
00081   sprintf(tt1, "gainBPM%d_%d",box,pmt);
00082   sprintf(tt2, "BPRS gain BPM-%d_%d, all tiles; average MIP (ADC)",box,pmt);
00083   h1MaP=new TH1F(tt1,tt2,40,0,40);
00084   h1MaT=new TH1F("aa4","BTOW tower: mean MIP; ADC-ped",25,0,50);
00085   h2Cr=new TH2F("aa5","BPRS comparison (both axis); slope (raw) ; average MIP (fit)",20,-0.15,-0.015,20,0,40);
00086 
00087   TString xtit=core; xtit+="; softID:";
00088   for(int i=0;i<softL.size();i++) {
00089     int softID=softL[i];
00090     xtit+=softID; xtit+=" ,  ";
00091     processOne(softID,core,i+1);
00092     // break;
00093     //if(i>1) break;
00094   }
00095   // return;
00096   ln00=new TLine(0,0,17,0); ln00->SetLineStyle(2);
00097   TString tit;
00098   c->cd(3); 
00099   h2DoutRaw->Draw("colz"); gPad->SetLogz(); h2DoutRaw->SetAxisRange(-10,60,"y");
00100   tit=core; tit+=", all values"+xtit;
00101   h2DoutRaw->SetTitle(tit);
00102   ln00->Draw();
00103 
00104   c->cd(2); 
00105   h2DoutMip->Draw("colz"); 
00106   tit=core; tit+=", reqire MIP @ TPC & BTOW"+xtit;
00107   h2DoutMip->SetTitle(tit);
00108   h2DoutMip->Rebin2D(1,3);
00109   ln00->Draw();
00110   h2DoutMip->SetAxisRange(-10,60,"y");
00111   
00112   c->cd(1);
00113   pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
00114   pad->Draw();
00115   pad->Divide(3,1);
00116   pad->cd(1);
00117   h1MaP->Draw(); 
00118 
00119   pad->cd(2);
00120   h2Cr->Draw("box");
00121 
00122   pad->cd(3);
00123   h1MaT->Draw(); //  h1MaT->SetLineColor(kGreen);
00124  
00125   if(doPS) {
00126     sprintf(core, "ps/bprsPmtPix%02d.ps",0);
00127     if(doPS==2) sprintf(core, "ps/bprsBPM%03d_%d.ps", box , pmt);
00128    can->Print(core);
00129   }
00130 
00131   //....... write summary histos to output file
00132   fout->cd();
00133   h2DoutRaw->Write();
00134   h2DoutMip->Write();
00135   h1MaP->Write();
00136   h2Cr->Write();
00137   //  fout->ls();
00138 
00139 }
00140 
00141 //------------------------
00142 //------------------------
00143 //------------------------
00144 void processOne(int softID, char *core0,int pix) {
00145   TString ttx=core0; ttx+=pix;
00146   can=new TCanvas(ttx,ttx,800,700);
00147   char core[1000];
00148   sprintf(core, "%s, softID=%d",core0,softID);
00149  
00150   c=makeTitle(can,core,pix);
00151   c->Divide(2,2);
00152   c->cd(1);
00153 
00154   TH1F * hrP= getSlice( h2DinRawP,softID,"bprsR");
00155   TH1F * htP= getSlice( h2DinTrP,softID,"bprsT");
00156   TH1F * hmP= getSlice( h2DinMipP,softID,"bprsM");
00157 
00158   TH1F * hrT= getSlice( h2DinRawT,softID,"btowR");
00159   TH1F * htT= getSlice( h2DinTrT,softID,"btowT");
00160   TH1F * hmT= getSlice( h2DinMipT,softID,"btowM");
00161 
00162   hrP->SetTitle("BPRS spectrum,  all events");
00163   hrP->Fit("expo","","R",13,50);
00164   TF1 *ff=hrP->GetFunction("expo");
00165   float slope=ff->GetParameter(1);
00166   hrP->Fit("gaus","+","R",-10,10);
00167  
00168   hmP->SetLineColor(kRed); hmP->SetFillColor(kRed);
00169   hmP->SetTitle("BPRS spectrum,  MIP : TPC, TPC+BTOW");
00170   hmP->Fit("gaus","","R",-10,50);
00171 
00172   TF1 *ff=hmP->GetFunction("gaus");
00173   ff->SetLineColor(kBlue);   ff->SetLineWidth(1);
00174   float mean=ff->GetParameter(1);  
00175   float meanS=ff->GetParameter(2);
00176   float cut_adcL=2;
00177   float cut_sig=1.8;
00178   if(mean>cut_adcL && meanS>cut_sig) { // cut if MIP peak is too close to ped
00179     h1MaP->Fill(mean);
00180     h2Cr->Fill(slope,mean);
00181   }
00182 
00183   hrT->SetTitle("BTOW spectrum,  all events");
00184   hrT->Fit("expo","","R",10,40);
00185   hmT->SetLineColor(kGreen);hmT->SetFillColor(kGreen);
00186   hmT->Fit("gaus","","R",-10,60);
00187   TF1 *ffT=hmT->GetFunction("gaus");
00188   ffT->SetLineColor(kBlue);   ffT->SetLineWidth(1);
00189   float meanT=ffT->GetParameter(1);
00190   float meanST=ffT->GetParameter(2);
00191   h1MaT->Fill(meanT);
00192  
00193   hmT->SetTitle("BTOW spectrum,  MIP : TPC, TPC+BPRS");
00194 
00195   //..... BPRS .....
00196   c->cd(1); 
00197   hrP->Draw(); gPad->SetLogy(); hrP->SetMinimum(0.7);
00198   hmP->Draw("same");  gPad->SetGrid();
00199   float xC=mean-meanS,yC=1e7;
00200   ln=new TLine(xC,0,xC,yC); ln->SetLineColor(kMagenta); ln->Draw();
00201 
00202   c->cd(3);
00203   hmP->Draw(); hmP->SetMaximum(htP->GetMaximum()/3); 
00204   htP->Draw("same"); 
00205  
00206   //..... BTOW.....
00207   c->cd(2); 
00208   hrT->Draw(); gPad->SetLogy(); hrT->SetMinimum(0.7);
00209   hmT->Draw("same");gPad->SetGrid();
00210 
00211   c->cd(4);
00212   hmT->Draw(); hmT->SetMaximum(htT->GetMaximum()/3); 
00213   htT->Draw("same"); 
00214  
00215   ln=new TLine(10,-3,25,-3); ln->SetLineColor(kMagenta); ln->Draw();ln->SetLineWidth(2);
00216   int nb=hmP->GetXaxis()->GetNbins();
00217   for(int b=1;b<=nb;b++) {
00218     h2DoutRaw->SetBinContent(pix,b,hrP->GetBinContent(b));
00219     h2DoutMip->SetBinContent(pix,b,hmP->GetBinContent(b));
00220   }
00221   
00222   FILE *fcsv=fopen("bprs+btowMipGain.csv","a");
00223   // fprintf(fcsv,"average MIP ADC for BPRS tiles and BTOW towers; ver=1.6\n");
00224   //fprintf(fcsv,"softID,PMB_pmt, pmt pix #, , bprs nMIP,  bprs avrMIP (adc), bprs sigMIP (adc), , btow nMIP,  btow avrMIP (adc), btow sigMIP (adc)\n");
00225   fprintf(fcsv,"%d, %s, %d,  , %d, %.2f, %.2f, ,  %d, %.2f, %.2f\n", 
00226           softID,core0+9,pix,(int)hmP->GetEntries(),mean,meanS,
00227           (int)hmT->GetEntries(),meanT,meanST     );
00228   fclose(fcsv);
00229   //..... compute ADC limmits
00230   for(int ibp=0;ibp<mxBTile;ibp++) {
00231     float adcL=5., adcH=40.; // default for _BTOW_
00232     if(ibp==1) adcL=3.5; // BPRS
00233     int statGain=(int)janDb_mipStat[ibp]->GetBinContent(softID);
00234      if(!statGain){ //valid gain
00235       float mean=janDb_mipMean[ibp]->GetBinContent(softID);
00236       float sig=janDb_mipSig[ibp]->GetBinContent(softID);
00237       adcL=mean-sig;
00238       adcH=mean+sig;
00239       if(ibp==0 && adcL<5) adcL=5; // BTOW
00240       if(ibp==1 && adcL<3.5) adcL=3.5;// BPRS
00241       if(adcH>2*mean) adcH=2*mean;
00242       printf("ibp=%d id=%4d MIP mean=%.1f sig=%.1f adcL=%.1f H=%.1f\n",ibp,softID,mean,sig,adcL,adcH);      
00243     }
00244     float y=1e5;
00245     TH1F *h1=hmT; if (ibp==1) h1=hmP;
00246     Lx=h1->GetListOfFunctions();    assert(Lx);
00247     lnL=new TLine( adcL,0, adcL,y); lnL->SetLineColor(kBlue);lnL->SetLineWidth(2);lnL->SetLineStyle(2);
00248     Lx->Add(lnL);
00249     lnH=new TLine( adcH,0, adcH,y); lnH->SetLineColor(kBlue);lnH->SetLineWidth(2);lnH->SetLineStyle(2);
00250     Lx->Add(lnH);
00251  } // Loop over B/P 
00252 
00253   //....... write 1D histos to output file
00254   fout->cd();
00255   char tt1[100], tt2[1000];
00256   sprintf(tt2, "id=%d %s pix=%d, ",softID,core0,pix);
00257   sprintf(tt1, "bprs%d",softID);
00258 
00259   TString Tt=tt1, T2=tt2;
00260   hrP->SetName(Tt+"a"); hrP->SetTitle(T2+"inclusive"); hrP->Write();
00261   htP->SetName(Tt+"t"); htP->SetTitle(T2+"MIP in TPC"); htP->Write();
00262   hmP->SetName(Tt+"m"); hmP->SetTitle(T2+"MIP in TPC+BTOW");hmP->Write();
00263  
00264   sprintf(tt2, "id=%d BTOW , ",softID);
00265   sprintf(tt1, "btow%d",softID);
00266   Tt=tt1, T2=tt2;
00267   hrT->SetName(Tt+"a"); hrT->SetTitle(T2+"inclusive"); hrT->Write();
00268   htT->SetName(Tt+"t"); htT->SetTitle(T2+"MIP in TPC"); htT->Write();
00269   hmT->SetName(Tt+"m"); hmT->SetTitle(T2+"MIP in TPC+BPRS");hmT->Write();
00270  //...... FINISH .....
00271   sprintf(core, "softID=%d",softID);
00272   can->SetTitle(core);
00273   if(doPS) {
00274     sprintf(core, "ps/bprsPmtPix%02d.ps",pix);
00275     if(doPS==2) sprintf(core, "ps/softID%04d.ps",softID);
00276     can->Print(core);
00277   }
00278 
00279 
00280  }
00281 
00282 
00283 //=================
00284 TH1F * getSlice(TH2F * h2, int id, char *ctit) {
00285   axX=h2->GetXaxis();
00286   float x1=axX->GetXmin();
00287   float x2=axX->GetXmax();
00288   int nbX=axX->GetNbins();
00289   printf("X-axis range  --> [%.1f, %.1f], nb=%d %s\n",x1,x2,nbX,axX->GetTitle());
00290 
00291   axY=h2->GetYaxis();
00292   float y1=axY->GetXmin();
00293   float y2=axY->GetXmax();
00294   int nbY=axY->GetNbins();
00295   printf("Y-axis range  --> [%.1f, %.1f], nb=%d\n",y1,y2,nbY);
00296 
00297   assert(id>=1 && id<=nbX);
00298   // do projections
00299   char txt1[100], txt2[1000];
00300   sprintf(txt1,"%s_id%d",ctit,id);
00301   sprintf(txt2,"%s soft id=%d;%s ",ctit, id,axY->GetTitle());
00302   TH1F*h=new TH1F(txt1,txt2,nbY,y1,y2); // working histo for 1-D spectrum
00303   
00304   int i;
00305   for(i=1;i<=nbY;i++) h->SetBinContent(i,h2->GetBinContent(id,i));
00306   float x1=-20, x2=70;
00307   h->SetAxisRange(x1,x2);
00308   h->SetEntries(h->Integral());
00309  
00310   return h;
00311 }
00312 
00313 
00314 //------------------------
00315 void buildSoftList(int box=11, int pmt=1) {
00316    fd1=new TFile("calib-nov-8-2008/map-softID-bprsPmt-Rory.root"); assert(fd1->IsOpen());  
00317   printf("Opened: %s\n",fd1->GetName()); fd1->ls();
00318   TH1I * mapBprsPmt=(TH1I*) fd1->Get("bprsPmt"); assert(mapBprsPmt);
00319   //mapBprsPmt->Draw();
00320   
00321   int y= box*10+pmt;
00322   for(int b=1;b<=4800;b++) {
00323     if(y!=(int)mapBprsPmt->GetBinContent(b)) continue;
00324     softL.push_back(b);
00325   }
00326   printf("found pmb=%d pmt=%d Bprs tiles =%d\n",box,pmt,softL.size());
00327 }
00328 
00329 //------------------------
00330 TPad *makeTitle(TCanvas *c,char *core, int page) {
00331 
00332   c->Range(0,0,1,1);
00333   TPad *pad0 = new TPad("pad0", "apd0",0.0,0.95,1.,1.);
00334   pad0->Draw();
00335   pad0->cd();
00336 
00337   TPaveText *pt = new TPaveText(0,0.,1,1,"br");
00338   pt->Draw();
00339   TDatime dt;
00340   TString txt2=core;
00341   txt2+=", pixel=";
00342   txt2+=page;
00343   txt2+=",  ";
00344   txt2+=dt.AsString();
00345   pt->AddText(txt2);
00346   txt2="--";
00347   pt->AddText(txt2);
00348 
00349   c->cd();
00350   pad = new TPad("pad1", "apd1",0.0,0.0,1,.95);
00351   pad->Draw();
00352   return pad;
00353 }
00354