TCanvas *c; double par[6]; double *epar; TFile *fdA[12]; FILE *wfd=stdout; // output web page FILE *gfd; //====================================== fitTower() { gStyle->SetStatW(0.22); gStyle->SetStatH(0.22); // f=new TFile("june18.hist.root"); // f=new TFile("/star/data05/scratch/balewski/2005-eemcCal/day49-hist/iter2-out/sum-sect5.hist.root"); TH1F *h1=new TH1F("mpv","MPV gated w/ MIP ; MPV of ADC-ped",40,-5,35); TH1F *h2=new TH1F("mpvE","relative error of MPV , MIP gated; err(MPV)/MPV ",50,0,0.5); // stupid root tricks: hDum=new TH1F("aa","bb",10,1,9); hDum->Fill(5); char cT='T'; openAll(cT); int sec=5; char core[100]; int eta; char sub='C'; for(eta=1;eta<=12;eta++) { int nErr=0, nOK=0; float mpvL=999, mpvH=0; h1->Reset(); h2->Reset(); fprintf(wfd," %d \n",eta); gStyle->SetOptStat(1001111); for(sec=1; sec<=12;sec++) { TFile *f=fdA[sec-1]; for(sub='A';sub<='E';sub++) { sprintf(core,"%02d%c%c%02d",sec,cT,sub,eta); TString coreT=core; ha=(TH1F*)f->Get("a"+coreT); hd=(TH1F*)f->Get("d"+coreT); // specical case of stuck low bits if(eta==7 && sub=='B' && (sec==4 || sec==6 || sec==8) ){ ha->Rebin(4); hd->Rebin(4); printf("tower=%s rebinned\n",core); } c=new TCanvas("aa","aa",400,400); c->Divide(1,2); c->cd(1); hDum->Draw(); gPad->SetLogy(); float mpv, mpvEr; TString errS=plotOne(ha,hd, mpv, mpvEr); printf("errS=%s=\n",errS.Data()); bool isBad=errS.Sizeof()>1; if(isBad) { // report error channel fprintf(wfd," %s ,\n",(coreT+"-"+errS).Data()); nErr++; } if(errS.Contains("mask")) continue; nOK++; //return; if(mpvL>mpv) mpvL=mpv; if(mpvHFill(mpv); if(mpv>0) { h2->Fill(mpvEr/mpv); } } // end of eta bin }// end of sector loop fprintf(wfd," %d %d\n",nOK,nErr); char txt[100],sumN[100], pdfN[100]; sprintf(pdfN,"%cfitEta%02d.pdf",cT,eta); sprintf(txt,"cat *%02d.ps | ps2pdf - %s",eta,pdfN); printf("%s\n",txt); system(txt); sprintf(txt,"mv %s /star/u/balewski/WWW/tmp/",pdfN); printf("%s\n",txt); system(txt); fprintf(wfd," PDF \n",pdfN); fprintf(wfd," %.1f to %.1f \n",mpvL,mpvH); gStyle->SetOptStat(1111111); sprintf(sumN,"mpv%c-eta%02d",cT,eta); c=new TCanvas(sumN,sumN,600,600); c->Divide(1,3); c->cd(1); h1->Draw(); c->cd(2); h2->Draw(); c->Print(); sprintf(txt,"ps2pdf %s.ps %s.pdf",sumN,sumN); printf("%s\n",txt); system(txt); sprintf(txt,"mv %s.pdf /star/u/balewski/WWW/tmp/",sumN); printf("%s\n",txt); system(txt); fprintf(wfd," PDF \n",sumN); fflush(wfd); }// end of loop over eta bins if(wfd!=stdout)fclose(wfd); fclose(gfd); return; } //================================= TString plotOne(TH1F *ha, TH1F *hd, float &MPV, float &MPVerr) { assert(ha); assert(hd); const float feta[]= {1.95,1.855,1.765,1.675,1.59,1.51,1.435,1.365,1.3,1.235,1.17,1.115}; char *core =ha->GetName()+1; TString coreT=core; float xMax=120; ha->SetAxisRange(-10,xMax); hd->SetAxisRange(-10,xMax/2.); int maxbina=ha->GetMaximumBin(); float xcenta=ha->GetBinCenter(maxbina); printf("\n============================\n working on %s\n",ha->GetName()); c->cd(1); ha->SetLineColor(4); ha->SetMinimum(0.9); char *func="gaus"; ha->Fit(func,"RQI","",xcenta-5,xcenta+5); gPad->SetLogy(); TF1* gausa=ha->GetFunction(func); gausa->SetLineWidth(1); gausa->SetLineColor(2); float meanA=gausa->GetParameter(1); float errorA=gausa->GetParError(1); hx=hd->Clone(); // memory leak hx->Draw("same"); c->cd(2); TF1 *f1 = new TF1("myfunc",myfunction,-10,100,5); f1->SetParNames("x0","aL","aG","sigL","sigG"); f1->SetLineColor(kRed); f1->SetLineWidth(1); hd->Fit("gaus","RQ","",0,30); TF1 *ff=hd->GetFunction("gaus"); f1->SetParameter(0,ff->GetParameter(1)); f1->SetParameter(1,ff->GetParameter(0)); f1->FixParameter(2,0); fitGausLand(hd); // do fitting c->Print(coreT+".ps"); // draw before QA on results changing axis limits //....... retrieve components f1->GetParameters(par); epar= f1->GetParErrors(); ha->SetAxisRange(-3,3); // printf("RR=%f\n",ha->Integral()/ha->GetEntries()); if(ha->Integral()/ha->GetEntries()>0.999) return "mask"; float meanD=par[0]; float errorD=epar[0]; //printf("ss=%s=\n",ha->GetName()+5); int ieta=atoi(core+4) -1; // printf("ssieta=%d\n",ieta); MPV=(meanD-meanA); MPVerr=sqrt(errorA*errorA+errorD*errorD+0.09); float gain=2.89*TMath::TanH(feta[ieta])*MPV; float sig=2.89*TMath::TanH(feta[ieta])*MPVerr; fprintf(gfd,"%s %.2f %.2f %5.1f %5.2f %5.2f %5.2f\n",core,gain,sig,meanD,errorD,meanA,errorA); if(fabs(meanA)>1.) return "pedOff"; if(ha->GetEntries() <=0) return "empty"; // deviation of gain from the goal TListIter it(hd->GetListOfFunctions()); TLine *w=(TLine*)it.Next() ; w->Print(); float del=MPV-w->GetX1(); printf(" MPV=%f, del=%f\n",MPV,del); float eps=0.30; if(MPV<(1-eps)*w->GetX1()) return "lowG"; if(MPV>(1+eps)*w->GetX1()) return "highG"; if(par[1]<0 || par[2]<0) return "amplNeg"; if(gain<10.) return "gainNeg"; if(fabs(epar[0])<0.05) return "singF"; // single bin pathology if(par[4]/par[0]>0.60) return "wideG"; hd->SetAxisRange(5,50); float sum=hd->Integral(); if(sum<150) return "lowS"; return ""; } //======================== void openAll(char cT) { int i; char txt[200]; for(i=0;i<12;i++) { sprintf(txt,"/star/data05/scratch/balewski/2005-eemcCal/day49-hist/iter2-out/sum-sect%d.hist.root",i+1); fdA[i]=new TFile(txt); assert(fdA[i]->IsOpen()); } sprintf(txt,"/star/u/balewski/WWW/tmp/gains%c-allSect.dat",cT); gfd=fopen(txt,"w"); assert(gfd); fprintf(gfd,"# final %c-gains from MIP using UxV\n",cT); fprintf(gfd,"# format: gain (ch/GeV) & errGain, GausLandauFit: MPV(ADC) & errMPV, pedFit: means(ADC) & err \n"); sprintf(txt,"/star/u/balewski/WWW/tmp/gains%c.html",cT); wfd=fopen(txt,"w"); assert(wfd); fprintf(wfd,"MIP -- > final %c-gains table \n",cT,cT); fprintf(wfd,"\n"); fprintf(wfd,"
etaBin problems # of good # of errors
+warn
spectra
(60 tiles)
MPV range
(ADC)
summary
plot\n"); } //--------------------------------------- Double_t myfunction(Double_t *x, Double_t *par) { Float_t xx =x[0]; float mpv=par[0]; float aL=par[1]*fabs(par[3]); float aG=par[1]; float sigL=2*par[3]; float sigG=par[4]; Double_t fland = TMath::Landau(xx,mpv,sigL); if(sigL!=0) fland/=sigL; Double_t fgaus = TMath::Gaus(xx,mpv,sigG); double f=aL*fland + aG*fgaus; return f; } // float fitGausLand(TH1F *h) { char *funcd="myfunc"; int maxbind=h->GetMaximumBin(); float sig=h->GetRMS(); float sum=h->Integral(); float xcentd=h->GetBinCenter(maxbind); float eta=atoi(h->GetName()+5); //printf("eta=%f sum=%f\n",eta,sum); TF1 *f1=(TF1 *)gROOT->GetFunction(funcd); f1->SetParameters(xcentd,0,0,sig,sig); float x1=xcentd-10; float x2=xcentd*2.; if(x1<0) x1=0; if(x2Fit(funcd,"R","",x1,x2); //h->Fit(funcd,"R","",20,70); }