00001 #if !defined(__CINT__)
00002
00003 #else
00004 #if !defined(__CINT__) || defined(__MAKECINT__)
00005
00006 #else
00007
00008 #endif
00009 #endif
00010
00011 #if !defined(__CINT__) || defined(__MAKECINT__)
00012 #include "Riostream.h"
00013 #include <stdio.h>
00014 #include "TROOT.h"
00015 #include "TSystem.h"
00016 #include "TMath.h"
00017 #include "TH1.h"
00018 #include "TH2.h"
00019 #include "TH3.h"
00020 #include "TStyle.h"
00021 #include "TF1.h"
00022 #include "TProfile.h"
00023 #include "TTree.h"
00024 #include "TChain.h"
00025 #include "TFile.h"
00026 #include "TNtuple.h"
00027 #include "TCanvas.h"
00028 #include "TFileSet.h"
00029 #include "TDataSetIter.h"
00030 #include "TDataSet.h"
00031 #include "TClassTable.h"
00032
00033 #include "TMinuit.h"
00034 #include "TSpectrum.h"
00035 #include "StBichsel/Bichsel.h"
00036 #include "TString.h"
00037 #include "TLine.h"
00038 #include "TText.h"
00039 #include "TList.h"
00040 #include "TPolyMarker.h"
00041 #include "TKey.h"
00042 #include "TLegend.h"
00043 #else
00044 class TMinuit;
00045 class TF1;
00046 class TH1F;
00047 class TH2F;
00048 class TH3F;
00049 class TProfile;
00050 class TH2D;
00051 class TCanvas;
00052 class TSpectrum;
00053 class TSystem;
00054 class Bichsel;
00055 #endif
00056 #include "Names.h"
00057 const Int_t NH = NHYPS/2;
00058 Int_t N = 0;
00059 const Int_t noPoints = 1000;
00060 Double_t X[noPoints];
00061 Double_t dX[noPoints];
00062 Double_t Nu[noPoints];
00063 Double_t Mu[noPoints];
00064 Double_t dMu[noPoints];
00065 Double_t Sigma[noPoints];
00066 Double_t dSigma[noPoints];
00067 TCanvas *canvas = 0;
00068 Double_t Xlog10bg, Ylog2dx, Z;
00069 TFile *newf = 0;
00070 const Char_t *NAMES[6] = {"e","p","K","pi","mu","d"};
00071
00072
00073
00074 static const Double_t peaks[5] = {0,
00075 1.33136848552729070e+00,
00076 5.14780361811221443e-01,
00077 4.22811868484973097e-01,
00078 2.55316350542050152e+00};
00080 struct peak_t {Double_t peak, sigma, mass; const Char_t *Name;};
00081
00082 #if 0
00083 static const peak_t Peaks[6] = {{0. , 0., 0.13956995},
00084 {1.425420, 0.105841, 0.93827231},
00085 {0.571553, 0.064864, 0.493677 },
00086 {0.419100, 0.004250, 0.51099907e-3},
00087 {2.656584, 0.119078, 1.875613},
00088 {0.004818, 0.002861, 0.1056584}};
00089 #else
00090 static const peak_t Peaks[6] = {
00091
00092
00093
00094
00095
00096
00097
00098 { 0. , 0., 0.13956995, "pion"},
00099 { 1.42574, 0.101741, 0.938272, "proton"},
00100 { 0.565411, 0.0616611, 0.493677, "kaon"},
00101 { 0.424919, 0.00408318, 0.000510999, "e"},
00102 { 2.65548, 0.123809, 1.87561, "deuteron"},
00103 { 0.000717144, 0.00490783, 0.105658, "mu"}};
00104 #endif
00105
00106 Bichsel *gBichsel = 0;
00107
00108 TF1 *func = 0;
00109 TF1 *LandauF = 0;
00110 class St_TpcSecRowCor;
00111
00112 TF1* Landau(){
00113 if (!LandauF)
00114 #if 0
00115 LandauF =
00116 new TF1("LandauF","exp([0]-0.5*((x-[1])/[2])*((x-[1])/[2])+exp([3]-0.5*((x-[4])/[5])*((x-[4])/[5])))",-5,10);
00117 Double_t params[6] = {
00118 -3.93739e+00,
00119 1.98550e+00,
00120 1.56338e+00,
00121 1.44692e+00,
00122 -4.93793e-01,
00123 1.54585e+00
00124 };
00125 Double_t sigma_p[4] = {
00126 6.67647e-01,
00127 -1.58690e-01,
00128 2.79764e-02,
00129 -1.67180e-03
00130 };
00131 Double_t sigmaI = 5.07430e-01;
00132 Double_t sigmaO = 3.80682e-01;
00133 #else
00134 LandauF =
00135 new TF1("LandauF","exp([0]-0.5*((x-[1])/[2])**2+exp([3]-0.5*((x-[4])/[5])**2+exp([6]-0.5*((x-[7])/[8])**2)))",-5,10);
00136
00137
00138 Double_t params[9] = {
00139 -7.74975e+00,
00140 6.53414e+00,
00141 1.21524e+00,
00142 3.31409e+00,
00143 -2.58291e+00,
00144 3.51463e+00,
00145 -3.47755e+00,
00146 3.77698e-02,
00147 6.67913e-01};
00148 #endif
00149 LandauF->SetParameters(params);
00150 return LandauF;
00151 }
00152
00153 Double_t fithfcn(Double_t *x,Double_t *par) {
00154 Double_t z = x[0];
00155 Double_t sigma = par[0];
00156 Double_t norm = 1./TMath::Sqrt(2*TMath::Pi())/sigma;
00157 Double_t value = 0;
00158
00159
00160 for (int k = 0; k < NH; k++) {
00161 Double_t mu = par[k+1];
00162
00163 Double_t dev = (z - mu)/sigma;
00164 value += norm*TMath::Exp(par[k+1+NH]-dev*dev/2.);
00165 }
00166 return value;
00167 }
00168
00169 void FitH(const Char_t *set="z", Int_t Hyp = -1, Int_t Bin=-1) {
00170 if (!gBichsel) {
00171 gSystem->Load("StBichsel");
00172 gBichsel = Bichsel::Instance();
00173 }
00174 TString Set(set);
00175 const Double_t window = 0.4;
00176 const Double_t range[2] = {-2., 4.};
00177 const Double_t LFrMin = -10;
00178 if (! canvas) {
00179 canvas = new TCanvas("FitHC","FitH Canvas");
00180 canvas->SetGrid();
00181 }
00182 const Int_t nHYPS = NHYPS;
00183 TH2 *hists[nHYPS];
00184 TProfile *histp[nHYPS];
00185 TFile *fRootFile = (TFile *) gDirectory->GetFile();
00186 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
00187 else printf("%s found\n",fRootFile->GetName());
00188 TString newfile("FitH");
00189 newfile += Set;
00190 newfile += gSystem->BaseName(fRootFile->GetName());
00191 TString FileN("FitPars");
00192 FileN += Set;
00193 FileN += ".h";
00194 TFile * f = 0;
00195 if (Bin < 0) f = new TFile(newfile.Data(),"recreate");
00196 TH1 *proj = 0;
00197 TF1 *g = new TF1("g",fithfcn,range[0],range[1],2+2*NH);
00198 g->SetParName(0,"Sigma"); g->SetParLimits(0,0.06,0.12);
00199 g->SetParName(1+2*NH,"Hyp");
00200 Int_t Nfit = 0;
00201 Int_t nh1 = 0;
00202 Int_t nh2 = NHYPS-1;
00203 if (Hyp >=0 ) {nh1 = nh2 = Hyp;}
00204 Int_t parstatus[NH];
00205 for (Int_t hyp = nh1; hyp<=nh2; hyp++) {
00206 Int_t kh = hyp%NH;
00207 Char_t *HistN = (Char_t *) HistNames[hyp];
00208 if (Set.Contains("70")) HistN = (Char_t *) HistNames70[hyp];
00209 hists[hyp] = (TH2 *) fRootFile->Get(HistN);
00210 if (!hists[hyp]) {printf("Cannot histogram %s\n",HistNames[hyp]); continue;}
00211 histp[hyp] = (TProfile *) fRootFile->Get(HistNameP[hyp]);
00212 if (!histp[hyp]) {printf("Cannot histogram %s\n",HistNameP[hyp]); continue;;}
00213 Int_t k = 0;
00214 if (hyp > NH) k = NH;
00215 for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
00216 Int_t lh = hypl%NH;
00217 TString name1("LF.");
00218 name1 += Names[hypl];
00219 TString name2("");
00220 name2 += Names[hypl];
00221 g->SetParName(lh+1+NH,name1.Data());
00222 g->SetParName(lh+1,name2.Data());
00223 }
00224 Int_t nx = hists[hyp]->GetNbinsX();
00225 Int_t jbin1 = 1;
00226 Int_t jbin2 = nx;
00227 if (Bin > 0) {jbin1 = jbin2 = Bin;}
00228 for (int jbin=jbin1; jbin<=jbin2; jbin++) {
00229 TString name(Form("%s_%i_%i",hists[hyp]->GetName(),hyp,jbin));
00230 proj = hists[hyp]->ProjectionY(name.Data(),jbin,jbin);
00231 Int_t ix1=proj->GetXaxis()->FindBin(-window);
00232 Int_t ix2=proj->GetXaxis()->FindBin( window);
00233 Double_t dinT = 5*proj->Integral();
00234 Double_t dint = proj->Integral(ix1,ix2);
00235 if (dint < 100.) {
00236 printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hists[hyp]->GetName(),jbin,hyp,dint);
00237 delete proj;
00238 continue;
00239 }
00240 for (Int_t lh = 0; lh < NH; lh++) {
00241 g->ReleaseParameter(lh+1);
00242 g->ReleaseParameter(lh+1+NH);
00243 parstatus[lh] = 0;
00244 }
00245 printf("hist:%s bin %i for hyp %i has %10.0f entries\n",hists[hyp]->GetName(),jbin,hyp,dint);
00246 Double_t bg = TMath::Power(10., hists[hyp]->GetBinCenter(jbin));
00247 Double_t pmom = Masses[hyp]*bg;
00248 if (pmom > 0.1) {parstatus[4] = 1; printf("fix muon\n");}
00249 if (kh == 4 && pmom > 0.1) continue;
00250 Double_t devZ[NH];
00251 Double_t zref = 0;
00252 if (TString(set) == "70") zref = 1.e-6*gBichsel->GetI70(TMath::Log10(bg),1.0);
00253 else zref = 1.e-6*TMath::Exp(gBichsel->GetMostProbableZ(TMath::Log10(bg),1.0));
00254
00255 Int_t iok = 1;
00256 for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
00257 Int_t lh = hypl%NH;
00258 Double_t z = 0;
00259 bg = pmom/Masses[hypl];
00260 if (Set.Contains("70")) z = 1.e-6*gBichsel->GetI70(TMath::Log10(bg),1.0);
00261 else z = 1.e-6*TMath::Exp(gBichsel->GetMostProbableZ(TMath::Log10(bg),1.0));
00262 devZ[lh] = TMath::Log(z/zref);
00263 }
00264 for (int l = 0; l < NH; l++) {
00265 if (l == kh) continue;
00266 if (devZ[l] < -2.0 || devZ[l] > 4.0) {
00267 printf("Fix %s dZ = %f\n", Names[NH*(hyp/NH)+l],devZ[l]);
00268 parstatus[l] = 1;
00269 continue;
00270 }
00271 if (TMath::Abs(devZ[l]) < 0.01 ||
00272 (hyp%NH == 3 && l == 4 && TMath::Abs(devZ[l]) < 0.04)) {
00273 printf("Fix %s dZ = %f\n", Names[NH*(hyp/NH)+l],devZ[l]);
00274 parstatus[l] = 1;
00275 continue;
00276 }
00277 }
00278 for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
00279 AGAIN:
00280 Int_t lh = hypl%NH;
00281 if (parstatus[lh]) continue;
00282 Double_t windowN = devZ[lh]-window;
00283 Double_t windowP = devZ[lh]+window;
00284 for (int m = 0; m < NH; m++) {
00285 if (m != lh && ! parstatus[m]) {
00286 Double_t dev = 0.5*(devZ[m] - devZ[lh]);
00287 if (dev < 0 && windowN < devZ[lh] + dev) windowN = devZ[lh] + dev;
00288 if (dev > 0 && windowP > devZ[lh] + dev) windowP = devZ[lh] + dev;
00289 }
00290 }
00291 #if 0
00292 printf("Set Limits for %s ref = %f in [%f,%f]\n",Names[hypl],devZ[lh],windowN,windowP);
00293 #endif
00294 g->SetParameter(lh+1,devZ[lh]);
00295 g->SetParLimits(lh+1,windowN, windowP);
00296 g->SetParLimits(lh+1+NH,LFrMin, TMath::Log(dinT));
00297 if (hypl == hyp)
00298 g->SetParameter(lh+1+NH,TMath::Log(dint));
00299 else {
00300 ix1=proj->GetXaxis()->FindBin(windowN);
00301 ix2=proj->GetXaxis()->FindBin(windowP);
00302 Double_t din= proj->Integral(ix1,ix2);
00303
00304 if (din > 0) g->SetParameter(lh+1+NH,TMath::Log(din));
00305 else {
00306 printf("Fix %s din = %f\n", Names[hypl],din);
00307 parstatus[lh] = 1;
00308 goto AGAIN;
00309 }
00310 }
00311 }
00312 for (int l = 0; l < NH; l++)
00313 if (parstatus[l]) {
00314 printf("Fix %s\n",Names[NH*(hyp/NH)+l]);
00315 g->FixParameter(l+1+NH, LFrMin);
00316 g->FixParameter(l+1, devZ[l]);
00317 }
00318 g->FixParameter(1+2*NH, hyp);
00319 if (parstatus[kh]) iok = 0;
00320 if (! iok) {printf("Too close\n"); continue;}
00321
00322 proj->Fit(g->GetName(),"R");
00323 canvas->Update();
00324 Int_t lh = (hyp%NH)+1;
00325 Double_t Nu = g->GetParameter(lh);
00326 printf("hyp = %i lh = %i Nu = %f zref =%f %f\n",hyp,lh,Nu,zref,g->GetParameter(lh));
00327 Double_t Mu = 1.e6*zref*TMath::Exp(Nu);
00328 Double_t dMu = g->GetParError(lh);
00329 Int_t NFitPoints = g->GetNumberFitPoints();
00330 Int_t NDF = g->GetNDF();
00331 Double_t prob = g->GetProb();
00332 Double_t chisq = g->GetChisquare();
00333 Double_t X = hists[hyp]->GetXaxis()->GetBinCenter(jbin);
00334 printf("Nu = %f Mu = %f dMu = %f\n",Nu,Mu,dMu);
00335 printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
00336 "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
00337 Names[hyp],hyp,jbin,Nfit,X,pmom,Nu,Mu,dMu,chisq,NFitPoints,NDF,prob);
00338 printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
00339 Names[hyp],hyp,jbin,Nfit,X,pmom,Nu,Mu,dMu,chisq,NFitPoints,NDF,prob,g->GetName());
00340 if (f) {
00341 FILE *fp = fopen(FileN.Data(),"a");
00342 if (fp) {
00343 if (Nfit == 0) {
00344 TDatime time;
00345 fprintf (fp,"dEdxPoint_t dEdxZ[] = {\n");
00346 fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
00347 fprintf(fp,
00348 "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
00349 }
00350 if (Nfit == 0) fprintf(fp," {");
00351 else fprintf(fp,",{");
00352 fprintf(fp,
00353 "\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
00354 Names[hyp],hyp,jbin,Nfit,X,pmom,Nu,Mu,dMu,chisq,NFitPoints,NDF,prob,g->GetName());
00355 fclose(fp);
00356 Nfit++;
00357 }
00358 proj->Write();
00359 }
00360 }
00361 }
00362 if (f) {
00363 FILE *fp = fopen(FileN.Data(),"a");
00364 if (fp) fprintf(fp,"};\n");
00365 fclose(fp);
00366 delete f;
00367 }
00368 }
00369
00370
00371
00372 TH2F *Project(TH3F *hist) {
00373 if (!hist) return 0;
00374 Int_t nx = hist->GetNbinsX();
00375 Int_t ny = hist->GetNbinsY();
00376 TString name(hist->GetName());
00377 name += "_xy";
00378 TH2F *h = new TH2F(name.Data(),hist->GetTitle(),
00379 nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
00380 ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
00381 Double_t params[5];
00382 Double_t error;
00383 TF1 *g = new TF1("g","gaus(0)+pol1(3)");
00384 for (int i=1;i<=nx;i++){
00385 for (int j=1;j<=ny;j++){
00386 TH1D *proj = hist->ProjectionZ("f_11",i,i,j,j);
00387 Double_t mean = proj->GetMean();
00388 Double_t sum = proj->Integral();
00389 Double_t rms = proj->GetRMS();
00390 params[0] = sum;
00391 params[1] = mean;
00392 params[2] = rms;
00393 params[3] = 0;
00394 params[4] = 0;
00395 g->SetParameters(params);
00396
00397 g->SetRange(-1.,1.);
00398 error = 0;
00399 if (sum > 100) {
00400 proj->Fit("g","RQ");
00401 g->GetParameters(params);
00402 error = g->GetParError(1);
00403 }
00404 else error = 0;
00405 h->SetCellContent(i,j,params[1]);
00406 h->SetCellError(i,j,error);
00407 printf("i:%i j:%i sum:%f mean:%f rms:%f mu:%f sigma:%f\n",
00408 i,j,sum,mean,rms,params[1],params[2]);
00409 delete proj;
00410 }
00411 }
00412 return h;
00413 }
00414
00415 TH2D *GetRelYZError(TH3 *hist, const Char_t *Name="_yz") {
00416 if (!hist) return 0;
00417 Int_t nx = hist->GetNbinsX();
00418 Int_t ny = hist->GetNbinsY();
00419 Int_t nz = hist->GetNbinsZ();
00420 TString name(hist->GetName());
00421 name += Name;
00422
00423 TAxis *fYaxis = hist->GetYaxis();
00424 TAxis *fZaxis = hist->GetZaxis();
00425 TH2D *h = new TH2D(name.Data(),"Relative error in Space charge (%)",
00426 ny,fYaxis->GetXmin(),fYaxis->GetXmax(),
00427 nz,fZaxis->GetXmin(),fZaxis->GetXmax());
00428 h->SetXTitle("Row number");
00429 h->SetYTitle("Z (cm)");
00430 Int_t bin;
00431 Double_t Scale = TMath::Sqrt(24.*3726./1153.);
00432 for (int iybin=0;iybin<ny;iybin++){
00433 for (int izbin=0;izbin<nz;izbin++){
00434 Double_t cont = 0, err = 0;
00435 for (int ixbin=1;ixbin<=nx;ixbin++){
00436 bin = hist->GetBin(ixbin,iybin,izbin);
00437 cont += hist->GetBinContent(bin);
00438 err += hist->GetBinError(bin)*hist->GetBinError(bin);
00439 }
00440 if (cont > 0) {
00441 bin = h->GetBin(iybin,izbin);
00442 Double_t val = 100*TMath::Sqrt(err)/cont*Scale;
00443 h->SetBinContent(bin,val);
00444 }
00445 }
00446 }
00447 return h;
00448 }
00449
00450 TH2F *ProjectX(TH3F *hist, const Char_t *Name="_yz",const Int_t binx1=0,const Int_t binx2=10000) {
00451 if (!hist) return 0;
00452 Int_t nx = hist->GetNbinsX();
00453 if (nx > binx2) nx = binx2;
00454 Int_t ny = hist->GetNbinsY();
00455 Int_t nz = hist->GetNbinsZ();
00456 TString name(hist->GetName());
00457 name += Name;
00458
00459 TAxis *fYaxis = hist->GetYaxis();
00460 TAxis *fZaxis = hist->GetZaxis();
00461 TH2F *h = new TH2F(name.Data(),hist->GetTitle(),
00462 ny,fYaxis->GetXmin(),fYaxis->GetXmax(),
00463 nz,fZaxis->GetXmin(),fZaxis->GetXmax());
00464
00465
00466
00467 for (int iybin=0;iybin<ny;iybin++){
00468 for (int izbin=0;izbin<nz;izbin++){
00469 for (int ixbin=binx1;ixbin<nx;ixbin++){
00470 Int_t bin = hist->GetBin(ixbin,iybin,izbin);
00471 Double_t cont = hist->GetBinContent(bin);
00472 if (cont) h->Fill(fYaxis->GetBinCenter(iybin),fZaxis->GetBinCenter(izbin), cont);
00473 }
00474 }
00475 }
00476 return h;
00477 }
00478
00479 TF1 *FitGP(TH1 *proj, Option_t *opt="RQ", Double_t nSigma=3, Int_t pow=3) {
00480 if (! proj) return 0;
00481 TString Opt(opt);
00482
00483 TF1 *g = 0, *g0 = 0;
00484 TF1 *gaus = (TF1*) gROOT->GetFunction("gaus");
00485 if (pow >= 0) g0 = new TF1("g0",Form("gaus(0)+pol%i(3)",pow),-0.2,0.2);
00486 else g0 = new TF1("g0","gaus",-0.2,0.2);
00487 g0->SetParName(0,"Constant");
00488 g0->SetParName(1,"Mean");
00489 g0->SetParName(2,"Sigma");
00490 for (int i=0; i<=pow;i++) g0->SetParName(3+i,Form("a%i",i));
00491 TF1 *g1 = new TF1("g1",Form("gaus(0)+pol%i(3)",pow+1),-0.2,0.2);
00492 g1->SetParName(0,"Constant");
00493 g1->SetParName(1,"Mean");
00494 g1->SetParName(2,"Sigma");
00495 for (int i=0; i<=pow+1;i++) g1->SetParName(3+i,Form("a%i",i));
00496 TF1 *g2 = new TF1("g2",Form("gaus(0)+pol%i(3)",pow+2),-0.2,0.2);
00497 g2->SetParName(0,"Constant");
00498 g2->SetParName(1,"Mean");
00499 g2->SetParName(2,"Sigma");
00500 for (int i=0; i<=pow+2;i++) g2->SetParName(3+i,Form("a%i",i));
00501 Double_t params[9];
00502 #if 0
00503 TSpectrum *spec = new TSpectrum();
00504 Int_t nPeaks = spec->Search(proj,3,"");
00505 if (!quet) {cout << proj->GetName() << "\tfound " << nPeaks << " peaks" << endl;}
00506 if (! nPeaks) return 0;
00507 Float_t *xpeak = spec->GetPositionX();
00508 Int_t peak = -1;
00509 Double_t dist = 9999.;
00510 for (Int_t i = 0; i < nPeaks; i++) {
00511 if (!quet) cout << "\tx \t" << xpeak[i];
00512 if (TMath::Abs(xpeak[i]) < dist) {
00513 dist = TMath::Abs(xpeak[i]); peak = i;
00514 }
00515 }
00516 if (peak < 0) return 0;
00517
00518 GetMaximumBin()
00519 if (!quet) {
00520 cout << endl;
00521 cout << "Take x\t" << xpeak[peak] << endl;
00522 }
00523
00524 Int_t peakbin = proj->GetXaxis()->FindBin(xpeak[peak]);
00525 params[0] = 1.e5;
00526 params[1] = proj->GetBinCenter(peakbin);
00527 params[2] = 0.2;
00528 #else
00529 Int_t peak = proj->GetMaximumBin();
00530 Double_t peakX = proj->GetBinCenter(peak);
00531 params[0] = proj->GetBinContent(peak);
00532 if (peakX > 0.5) {
00533 params[1] = 0;
00534 params[2] = 0.2;
00535 }
00536 else {
00537 params[1] = peakX;
00538 params[2] = proj->GetRMS();
00539 if (params[2] > 0.25) params[2] = 0.25;
00540 }
00541 #endif
00542 params[3] = 0;
00543 params[4] = 0;
00544 params[5] = 0;
00545 params[6] = 0;
00546 params[7] = 0;
00547 params[8] = 0;
00548 if (gaus) {
00549 g = gaus;
00550 g->SetParameters(params);
00551 g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
00552 proj->Fit(g,opt);
00553 g->GetParameters(params);
00554 if (g->GetProb() > 0.01) return g;
00555 params[2] = TMath::Abs(params[2]);
00556 }
00557 g = g0;
00558 g->SetParameters(params);
00559 g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
00560 proj->Fit(g,opt);
00561 #if 1
00562 if (g->GetProb() > 0.01) return g;
00563 g->GetParameters(params);
00564 g = g1;
00565 params[2] = TMath::Abs(params[2]);
00566 g->SetParameters(params);
00567 g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
00568 proj->Fit(g,opt);
00569 if (g->GetProb() > 0.01) return g;
00570 g->GetParameters(params);
00571 g = g2;
00572 params[2] = TMath::Abs(params[2]);
00573 g->SetParameters(params);
00574 g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
00575 proj->Fit(g,opt);
00576 #endif
00577 if (! Opt.Contains("q",TString::kIgnoreCase)) {
00578 g->GetParameters(params);
00579 Double_t X = params[1];
00580 Double_t Y = params[0]/TMath::Sqrt(2*TMath::Pi()*params[2]);
00581 TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
00582 proj->GetListOfFunctions()->Add(pm);
00583 pm->SetMarkerStyle(23);
00584 pm->SetMarkerColor(kRed);
00585 pm->SetMarkerSize(1.3);
00586 proj->Draw();
00587 }
00588 return g;
00589 }
00590
00591 Double_t gfFunc(Double_t *x, Double_t *par) {
00592
00593
00594
00595
00596
00597
00598
00599
00600 Double_t sigma = par[2];
00601 Double_t frac[5];
00602 Int_t i;
00603 frac[0] = 1;
00604 for (i = 1; i < 5; i++) {
00605 frac[i] = TMath::Sin(par[2+i]);
00606 frac[i] *= frac[i];
00607 frac[0] -= frac[i];
00608 }
00609 if (frac[0] < 0.4) return 0;
00610 Double_t Value = 0;
00611 Int_t icase = (Int_t) par[8];
00612 Int_t i1 = 0;
00613 Int_t i2 = 4;
00614 if (icase >= 0) {i1 = i2 = icase;}
00615 for (i = i1; i <= i2; i++) {
00616 Double_t Sigma = TMath::Sqrt(sigma*sigma + Peaks[i].sigma*Peaks[i].sigma);
00617 Value += frac[i]*TMath::Gaus(x[0],par[1]+Peaks[i].peak,Sigma,1);
00618
00619 }
00620 return par[7]*TMath::Exp(par[0])*Value;
00621 }
00622
00623 TF1 *FitGF(TH1 *proj, Option_t *opt="") {
00624
00625 if (! proj) return 0;
00626 TString Opt(opt);
00627
00628 TF1 *g2 = (TF1*) gROOT->GetFunction("GF");
00629 if (! g2) {
00630 g2 = new TF1("GF",gfFunc, -5, 5, 9);
00631 g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
00632 g2->SetParName(1,"mu");
00633 g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.2,0.8);
00634 g2->SetParName(3,"P");
00635 g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.5);
00636 g2->SetParName(5,"e"); g2->SetParLimits(5,0.0,0.5);
00637 g2->SetParName(6,"d"); g2->SetParLimits(6,0.0,0.5);
00638 g2->SetParName(7,"Total");
00639 g2->SetParName(8,"Case");
00640
00641 }
00642
00643 Double_t total = proj->Integral()*proj->GetBinWidth(5);
00644 g2->SetParameters(0, 1e-3, 0.32, 0.4, 0.1, 0.1, 0.1,0.1,-1.);
00645 g2->FixParameter(3,2.86731e-01);
00646 g2->FixParameter(4,4.27788e-01);
00647 g2->FixParameter(5,1e-6);
00648 g2->FixParameter(6,5.77664e-02);
00649 g2->FixParameter(7,total);
00650 g2->FixParameter(8,-1);
00651 proj->Fit(g2,Opt.Data());
00652 g2->ReleaseParameter(3); g2->SetParLimits(3,0.0,TMath::Pi()/2);
00653 g2->ReleaseParameter(4); g2->SetParLimits(4,0.0,TMath::Pi()/2);
00654 g2->ReleaseParameter(5); g2->SetParLimits(5,0.0,TMath::Pi()/2);
00655 g2->ReleaseParameter(6); g2->SetParLimits(6,0.0,TMath::Pi()/2);
00656 Int_t iok = proj->Fit(g2,Opt.Data());
00657 if ( iok ) {
00658 cout << g2->GetName() << " fit has failed with " << iok << " for "
00659 << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
00660 proj->Fit(g2,Opt.Data());
00661 }
00662 Opt += "m";
00663 iok = proj->Fit(g2,Opt.Data());
00664 if (! Opt.Contains("q",TString::kIgnoreCase)) {
00665 Double_t params[10];
00666 g2->GetParameters(params);
00667 Double_t X = params[1];
00668 Double_t Y = TMath::Exp(params[0]);
00669 TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
00670 proj->GetListOfFunctions()->Add(pm);
00671 pm->SetMarkerStyle(23);
00672 pm->SetMarkerColor(kRed);
00673 pm->SetMarkerSize(1.3);
00674 for (int i = 0; i <= 4; i++) {
00675 TF1 *f = new TF1(*g2);
00676 f->SetName(Peaks[i].Name);
00677 f->FixParameter(8,i);
00678 f->SetLineColor(i+2);
00679 proj->GetListOfFunctions()->Add(f);
00680 }
00681 proj->Draw();
00682 }
00683 return g2;
00684 }
00685
00686 Double_t gbFunc(Double_t *x, Double_t *par) {
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 #if 1
00697 if (! LandauF) Landau();
00698 static Double_t sigma_p[3] = {
00699 5.31393e-01,
00700 -1.43277e-01,
00701 2.43800e-02
00702 };
00703 Double_t sigmaC = par[2];
00704 Double_t frac[5];
00705 Int_t i;
00706 frac[0] = 1;
00707 for (i = 1; i < 5; i++) {
00708 frac[i] = TMath::Sin(par[2+i]);
00709 frac[i] *= frac[i];
00710 frac[0] -= frac[i];
00711 }
00712 Double_t Value = 0;
00713 static Double_t pMom = 0.475;
00714 static Double_t Xlog10bg[5];
00715 Double_t Ylog2dx = TMath::Log2(par[8]);
00716 Double_t Sigma = sigma_p[2];
00717 for (int n = 1; n>=0; n--) Sigma = Ylog2dx*Sigma + sigma_p[n];
00718 Double_t zMostProb[5];
00719 for (i = 0; i < 5; i++) {
00720 Xlog10bg[i] = TMath::Log10(pMom/Peaks[i].mass);
00721 zMostProb[i] = gBichsel->GetMostProbableZ(Xlog10bg[i],Ylog2dx);
00722 Double_t sigma = Sigma + sigmaC;
00723
00724
00725 Double_t xi = (x[0] - par[1] - Peaks[i].peak)/sigma;
00726 Double_t Prob = LandauF->Eval(xi);
00727
00728 Value += frac[i]*Prob;
00729
00730 }
00731 return par[7]*TMath::Exp(par[0])*Value;
00732 #else
00733 Double_t sigmaC = par[2];
00734 Double_t frac[5];
00735 Int_t i;
00736 frac[0] = 1;
00737 for (i = 1; i < 5; i++) {
00738 frac[i] = TMath::Sin(par[2+i]);
00739 frac[i] *= frac[i];
00740 frac[0] -= frac[i];
00741 }
00742 Double_t Value = 0;
00743 static Double_t pMom = 0.475;
00744 static Double_t Xlog10bg[5];
00745 Double_t Ylog2dx = TMath::Log2(par[8]);
00746 Double_t zMostProb[5];
00747 for (i = 0; i < 5; i++) {
00748 Xlog10bg[i] = TMath::Log10(pMom/Peaks[i].mass);
00749 zMostProb[i] = gBichsel->GetMostProbableZ(Xlog10bg[i],Ylog2dx);
00750 Double_t sigma = gBichsel->GetRmsZ(Xlog10bg[i],Ylog2dx) + sigmaC;
00751
00752 Double_t xi = (x[0] - par[1] - Peaks[i].peak)/sigma;
00753 Double_t Phi = gBichsel->GetProbability(Xlog10bg[i],Ylog2dx,xi);
00754 Double_t Prob = Phi/sigma;
00755
00756 Value += frac[i]*Prob;
00757
00758 }
00759 return par[7]*TMath::Exp(par[0])*Value;
00760 #endif
00761 }
00762
00763 TF1 *FitGB(TH1 *proj, Option_t *opt="", Double_t dX = 2.364) {
00764 if (!gBichsel) {
00765 gSystem->Load("StBichsel");
00766 gBichsel = Bichsel::Instance();
00767 }
00768
00769 if (! proj) return 0;
00770 TString Opt(opt);
00771
00772 TF1 *g2 = (TF1*) gROOT->GetFunction("GB");
00773 if (! g2) {
00774 g2 = new TF1("GB",gbFunc, -5, 5, 9);
00775 g2->SetParName(0,"norm");
00776 g2->SetParName(1,"mu"); g2->SetParLimits(1,-1.5,1.5);
00777 g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.0,0.8);
00778 g2->SetParName(3,"P");
00779 g2->SetParName(4,"K");
00780 g2->SetParName(5,"e");
00781 g2->SetParName(6,"d");
00782 g2->SetParName(7,"Total");
00783 g2->SetParName(8,"dX");
00784
00785 }
00786
00787 Double_t total = proj->Integral()*proj->GetBinWidth(5);
00788 #if 0
00789 TAxis *xx = proj->GetXaxis();
00790 for (int i = 0; i <5; i++) {
00791 Double_t x = 0;
00792 if (i) x = peaks[i];
00793 Int_t bin = xx->FindBin(x);
00794 Int_t p = 0;
00795 if (i != 0) p = i+3;
00796 Double_t cont = proj->GetBinContent(bin);
00797 g2->ReleaseParameter(p);
00798 if (cont > 0) g2->SetParameter(p,TMath::Log(cont));
00799 else g2->FixParameter(p,-99.);
00800 }
00801 #endif
00802 g2->SetParameters(0, 1e-3, 0.01, 0.4, 0., 0., 0.,0.);
00803 g2->FixParameter(7,total);
00804 g2->FixParameter(8,dX);
00805 Int_t iok = proj->Fit(g2,Opt.Data());
00806 if ( iok ) {
00807 cout << g2->GetName() << " fit has failed with " << iok << " for "
00808 << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
00809 proj->Fit(g2,Opt.Data());
00810 }
00811 Opt += "m";
00812 iok = proj->Fit(g2,Opt.Data());
00813 if (! Opt.Contains("q",TString::kIgnoreCase)) {
00814 Double_t params[10];
00815 g2->GetParameters(params);
00816 Double_t X = params[1];
00817 Double_t Y = TMath::Exp(params[0]);
00818 TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
00819 proj->GetListOfFunctions()->Add(pm);
00820 pm->SetMarkerStyle(23);
00821 pm->SetMarkerColor(kRed);
00822 pm->SetMarkerSize(1.3);
00823 proj->Draw();
00824 }
00825 return g2;
00826 }
00827
00828 TF1 *FitG2(TH1 *proj, Option_t *opt="RQ") {
00829 if (! proj) return 0;
00830 Double_t params[9];
00831 TF1 *gaus = new TF1("gaus","gaus",-5.,5.);
00832 proj->Fit(gaus,opt);
00833 params[0] = gaus->GetParameter(0);
00834 params[1] = gaus->GetParameter(1);
00835 params[2] = gaus->GetParameter(2);
00836 params[3] = 1;
00837 params[4] = 2;
00838 params[5] = gaus->GetParameter(2);
00839 params[6] = 0;
00840 params[7] = 0;
00841 params[8] = 0;
00842 TF1 *g = new TF1("g","gaus(0)+gaus(3)",-5.,5.);
00843 g->SetParameters(params);
00844 g->SetParLimits(0,0,1.e10);
00845 g->SetParLimits(1,-5,5);
00846 g->SetParLimits(2,0.01,5);
00847 g->SetParLimits(3,0,1.e10);
00848 g->SetParLimits(4,-5,5);
00849 g->SetParLimits(5,0.01,5);
00850 proj->Fit(g,opt);
00851 g->GetParameters(params);
00852 if (TMath::Abs(params[1]) > TMath::Abs(params[4])) {
00853 g->SetParameter(0,params[3]);
00854 g->SetParameter(1,params[4]);
00855 g->SetParameter(2,params[5]);
00856 g->SetParameter(3,params[0]);
00857 g->SetParameter(4,params[1]);
00858 g->SetParameter(5,params[2]);
00859 proj->Fit(g,opt);
00860 }
00861 proj->Fit(g,opt);
00862 delete gaus;
00863 return g;
00864 }
00865
00866 TF1 *FitG3(TH1 *proj, Option_t *opt="RQ") {
00867 if (! proj) return 0;
00868 Double_t params[9];
00869 TF1 *gaus = new TF1("gaus","gaus",-5.,5.);
00870 TF1 *g = 0;
00871 proj->Fit(gaus,opt);
00872 params[0] = gaus->GetParameter(0);
00873 params[1] = gaus->GetParameter(1);
00874 params[2] = gaus->GetParameter(2);
00875 params[3] = 1;
00876 params[4] = 2;
00877 params[5] = gaus->GetParameter(2);
00878 params[6] = 0;
00879 params[7] = 0;
00880 params[8] = 0;
00881 delete gaus;
00882 g = new TF1("g","gaus(0)+gaus(3)",-5.,5.);
00883 g->SetParameters(params);
00884 g->SetParLimits(0,0,1.e10);
00885 g->SetParLimits(1,-5,5);
00886 g->SetParLimits(2,0.1,5);
00887 g->SetParLimits(3,0,1.e10);
00888 g->SetParLimits(4,-5,5);
00889 g->SetParLimits(5,0.1,5);
00890 proj->Fit(g,opt);
00891 g->GetParameters(params);
00892 if (TMath::Abs(params[1]) > TMath::Abs(params[4])) {
00893 g->SetParameter(0,params[3]);
00894 g->SetParameter(1,params[4]);
00895 g->SetParameter(2,params[5]);
00896 g->SetParameter(3,params[0]);
00897 g->SetParameter(4,params[1]);
00898 g->SetParameter(5,params[2]);
00899 }
00900 proj->Fit(g,opt);
00901 if (g->GetProb() > 0) return g;
00902 g->GetParameters(params);
00903 delete g;
00904 g = new TF1("g","gaus(0)+gaus(3)+gaus(6)",-5.,5.);
00905 g->SetParLimits(6,0,1.e10);
00906 g->SetParLimits(7,-5,5);
00907 g->SetParLimits(8,0.1,5);
00908 params[6] = 10.;
00909 params[7] = 0.5*(params[1]+params[4]);
00910 params[8] = params[2];
00911 g->SetParameters(params);
00912 proj->Fit(g,opt);
00913 return g;
00914 }
00915
00916 void FitB4G(Int_t icase = 0, Int_t hyp=-1, Int_t bin=0,
00917 Double_t xmin1=-1.0, Double_t xmax1 = 1.0,
00918 #if 0
00919 Double_t Mu2 =-3.0,Double_t xmin2=-.2, Double_t xmax2 = .2,
00920 Double_t Mu3 = 0.0,Double_t xmin3=-2., Double_t xmax3 = 4.,
00921 Double_t Mu4 = 0.0
00922 #else
00923 Double_t Mu2 = -.5,Double_t xmin2=-2., Double_t xmax2 = 4.,
00924 Double_t Mu3 = 0.1,Double_t xmin3=-2., Double_t xmax3 = 4.,
00925 Double_t Mu4 = 0.3
00926 #endif
00927 )
00928 {
00929 Double_t sigmas[2] = {0.06,0.12};
00930 const Int_t nHYPS = NHYPS;
00931 TH2 *hists[nHYPS];
00932 TProfile *histp[nHYPS];
00933 TFile *fRootFile = (TFile *) gDirectory->GetFile();
00934 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
00935 TFile *f = 0;
00936 if (bin == 0) {
00937 TString newfile("FitB4G");
00938 newfile += gSystem->BaseName(fRootFile->GetName());
00939 f = new TFile(newfile.Data(),"update");
00940 }
00941 Int_t i;
00942 for (i = 0; i< nHYPS; i++) {
00943 if (hyp >= 0 && hyp != i) continue;
00944 if (! icase) {
00945 hists[i] = (TH2 *) fRootFile->Get(HistNames[i]);
00946 if (!hists[i]) {printf("Cannot histogram %s\n",HistNames[i]); return;}
00947 }
00948 else {
00949 hists[i] = (TH2 *) fRootFile->Get(HistNames70[i]);
00950 if (!hists[i]) {printf("Cannot histogram %s\n",HistNames70[i]); return;}
00951 }
00952 histp[i] = (TProfile *) fRootFile->Get(HistNameP[i]);
00953 if (!histp[i]) {printf("Cannot histogram %s\n",HistNameP[i]); return;}
00954 }
00955 TF1 *g1 = new TF1("g1","gaus",xmin1,xmax1);
00956 TF1 *g = 0, *g2 = 0, *g3 = 0, *g4 = 0, *ga = 0;
00957 TCanvas *canvas = new TCanvas("canvas","canvas");
00958 if (Mu2 != 0) {
00959 if (Mu2 < -1.) {
00960 if (Mu2 <= -4.) ga = new TF1("ga","gaus(0)+exp(pol3(3))",xmin2,xmax2);
00961 else {
00962 if (Mu2 <= -3.) ga = new TF1("ga","gaus(0)+exp(pol2(3))",xmin2,xmax2);
00963 else ga = new TF1("ga","gaus(0)+exp(pol1(3))",xmin2,xmax2);
00964 }
00965 g = ga;
00966 g->SetParName(0,"Constant1");
00967 g->SetParName(1,"Mean1");
00968 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
00969 g->SetParName(3,"Const");
00970 g->SetParName(4,"Slope1");
00971 g->SetParName(5,"Slope2");
00972 g->SetParName(6,"Slope3");
00973 g->SetParName(7,"Slope4");
00974 }
00975 else {
00976 g2 = new TF1("g2","gaus(0)+gaus(3)",xmin2,xmax2);
00977 g = g2;
00978 g->SetParName(0,"Constant1");
00979 g->SetParName(1,"Mean1");
00980 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
00981 g->SetParName(3,"Constant2");
00982 g->SetParName(4,"Mean2");
00983 g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
00984 if (Mu3 != 0) {
00985 g3= new TF1("g3","gaus(0)+gaus(3)+gaus(6)",xmin3,xmax3);
00986 g = g3;
00987 g->SetParName(0,"Constant1");
00988 g->SetParName(1,"Mean1");
00989 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
00990 g->SetParName(3,"Constant2");
00991 g->SetParName(4,"Mean2");
00992 g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
00993 g->SetParName(6,"Constant3");
00994 g->SetParName(7,"Mean3");
00995 g->SetParName(8,"Sigma3"); g->SetParLimits(8,sigmas[0],sigmas[1]);
00996 if (Mu4 != 0) {
00997 g4= new TF1("g4","gaus(0)+gaus(3)+gaus(6)+gaus(9)",xmin3,xmax3);
00998 g4->SetParName(0,"Constant1");
00999 g4->SetParName(1,"Mean1");
01000 g4->SetParName(2,"Sigma1"); g4->SetParLimits(2,sigmas[0],sigmas[1]);
01001 g4->SetParName(3,"Constant2");
01002 g4->SetParName(4,"Mean2");
01003 g4->SetParName(5,"Sigma2"); g4->SetParLimits(5,sigmas[0],sigmas[1]);
01004 g4->SetParName(6,"Constant3");
01005 g4->SetParName(7,"Mean3");
01006 g4->SetParName(8,"Sigma3"); g4->SetParLimits(8,sigmas[0],sigmas[1]);
01007 g4->SetParName(9,"Constant4");
01008 g4->SetParName(10,"Mean4");
01009 g4->SetParName(11,"Sigma4"); g4->SetParLimits(11,sigmas[0],sigmas[1]);
01010 }
01011 }
01012 }
01013 }
01014 TH1 *proj = 0;
01015 Double_t params[12];
01016 Int_t k;
01017 Int_t Bin = TMath::Abs(bin);
01018 for (k = 0; k<nHYPS; k++) {
01019 if (hyp != -1 && k != hyp) continue;
01020 TH2 *hist = hists[k];
01021 Int_t nx = hist->GetNbinsX();
01022 for (i=1; i<=nx; i++) {
01023 if (bin != 0 && Bin != i) continue;
01024
01025 Char_t line[40];
01026 sprintf(line,"%s_%i",hist->GetName(),i);
01027 TString name(line);
01028 proj = hist->ProjectionY(name.Data(),i,i);
01029 Int_t ix1=proj->GetXaxis()->FindBin(-.1);
01030 Int_t ix2=proj->GetXaxis()->FindBin(0.1);
01031 if (proj->Integral(ix1,ix2) < 100.) {
01032 printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hist->GetName(),Bin,k,proj->Integral());
01033 delete proj;
01034 continue;
01035 }
01036 Int_t NFitPoints = 0;
01037 Double_t chisq;
01038 Double_t xrange = 0.6;
01039 g = g1;
01040 g->SetParLimits(0,0,1.e7);
01041 g->SetParLimits(1,-xrange,xrange);
01042 g->SetParLimits(2,sigmas[0],sigmas[1]);
01043 proj->Fit(g->GetName(),"r");
01044 g->GetParameters(params);
01045 Int_t kcase = -1;
01046 if (xmin1 < -0.5 && xmax1 > 0.5 && g->GetChisquare() < 1.e3
01047 ) goto Done;
01048
01049 if (ga) {
01050 g = ga;
01051 g1->GetParameters(params);
01052 params[1] = 0.;
01053 params[3] = 0.;
01054 params[4] = 0.;
01055 params[5] = 0.;
01056 params[6] = 0.;
01057 params[7] = 0.;
01058 params[8] = 0.;
01059 g->SetParameters(params);
01060 kcase = 0;
01061 if (g->GetChisquare() < 1.e3
01062 ) goto Done;
01063 }
01064 else {
01065 if (g2) {
01066 params[1] = 0.;
01067 params[3] = params[0];
01068 params[4] = Mu2;
01069 params[5] = 2.0*params[2];
01070 g = g2;
01071 g->SetParameters(params);
01072 g->SetParLimits(0,0,1.e7);
01073 g->SetParLimits(1,-xrange,xrange);
01074
01075 g->SetParLimits(2,sigmas[0],sigmas[1]);
01076 g->SetParLimits(3,0,1.e7);
01077 g->SetParLimits(5,sigmas[0],sigmas[1]);
01078 proj->Fit(g->GetName(),"r");
01079 g->GetParameters(params);
01080 kcase = 2;
01081 if (g->GetChisquare() < 1.e3
01082 ) goto Done;
01083 if (g3) {
01084 params[1] = 0.;
01085 params[6] = params[0];
01086 params[7] = Mu3;
01087 params[8] = 2.0*params[2];
01088 g = g3;
01089 g->SetParameters(params);
01090 g->SetParLimits(0,0,1.e7);
01091 g->SetParLimits(1,-xrange,xrange);
01092
01093 g->SetParLimits(2,sigmas[0],sigmas[1]);
01094 g->SetParLimits(3,0,1.e7);
01095 g->SetParLimits(5,sigmas[0],sigmas[1]);
01096 g->SetParLimits(6,0,1.e7);
01097 g->SetParLimits(8,sigmas[0],sigmas[1]);
01098 proj->Fit(g->GetName(),"r");
01099 g->GetParameters(params);
01100 kcase = 3;
01101 if (g->GetChisquare() < 1.e3
01102 ) goto Done;
01103 if (g4) {
01104 params[1] = 0.;
01105 params[9] = params[0];
01106 params[10] = Mu4;
01107 params[11] = 2.0*params[2];
01108 g = g4;
01109 g->SetParameters(params);
01110 g->SetParLimits(0,0,1.e7);
01111 g->SetParLimits(1,-xrange,xrange);
01112
01113 g->SetParLimits(2,sigmas[0],sigmas[1]);
01114 g->SetParLimits(3,0,1.e7);
01115 g->SetParLimits(5,sigmas[0],sigmas[1]);
01116 g->SetParLimits(6,0,1.e7);
01117 g->SetParLimits(8,sigmas[0],sigmas[1]);
01118 g->SetParLimits(9,0,1.e7);
01119 g->SetParLimits(11,sigmas[0],sigmas[1]);
01120 proj->Fit(g->GetName(),"r");
01121 kcase = 4;
01122 g->GetParameters(params);
01123 }
01124 }
01125 }
01126 }
01127 Done:
01128 if (g) {
01129 proj->Fit(g->GetName(),"RM");
01130 canvas->Update();
01131 Int_t l = 1;
01132 Double_t mu = g->GetParameter(l);
01133 for (int m=2;m<=kcase;m++) {
01134 if (TMath::Abs(mu) > TMath::Abs(g->GetParameter(3*m-2))) {
01135 l = 3*m-2; mu = g->GetParameter(l);
01136
01137 }
01138 }
01139 Nu[N] = g->GetParameter(l);
01140 Mu[N] = Nu[N] + histp[k]->GetBinContent(i);
01141 dMu[N] = g->GetParError(l);
01142
01143 Sigma[N] = g->GetParameter(2);
01144 dSigma[N] = g->GetParError(2);
01145 NFitPoints = g->GetNumberFitPoints();
01146 Int_t NDF = g->GetNDF();
01147 Double_t prob = g->GetProb();
01148 chisq = g->GetChisquare();
01149 X[N] = hist->GetXaxis()->GetBinCenter(i);
01150 dX[N] = hist->GetXaxis()->GetBinWidth(i);
01151 Double_t pionM = Masses[k]*pow(10.,X[N]);
01152 printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
01153 "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
01154 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob);
01155
01156
01157
01158 if ( bin >= 0) {
01159 printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
01160 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
01161 TString FileN("FitPars");
01162 if (hyp > -1) FileN += hist->GetName();
01163 FileN += ".h";
01164 FILE *fp = fopen(FileN.Data(),"a");
01165 if (fp) {
01166 if (N == 0) {
01167 TDatime time;
01168 fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
01169 fprintf(fp,
01170 "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
01171 }
01172 fprintf(fp,
01173 "{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
01174 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
01175 fclose(fp);
01176 }
01177 #if 1
01178 proj->Write();
01179 #endif
01180 }
01181 else printf ("================== Skip it\n");
01182 }
01183 N++;
01184 proj->Draw();
01185 }
01186 }
01187 if (f) {delete f;}
01188 }
01189
01190 void dEdxFit(const Char_t *HistName = "Time",const Char_t *FitName = "GP",
01191 Option_t *opt="R",
01192 Int_t mergeX=1, Int_t mergeY=1, Int_t ix = 0, Int_t jy = 0,
01193 Double_t nSigma=3, Int_t pow=1) {
01194 TCanvas *canvas = 0;
01195 TString Opt(opt);
01196 if (! Opt.Contains("Q",TString::kIgnoreCase)) canvas = new TCanvas("Fit","Fit results");
01197 TList *list = (TList *) gROOT->GetListOfFiles();
01198 if (! list) {printf("File list is empty\n"); return;}
01199 TIter next(list);
01200 TFile *fRootFile = (TFile *) next();
01201 if (! fRootFile ) {printf("There is no opened file\n"); return;}
01202 TH1 *hist = 0;
01203 fRootFile->GetObject(HistName,hist);
01204 if (!hist) {printf("Cannot find %s\n", HistName); return;}
01205 TAxis *xax = hist->GetXaxis();
01206 Int_t nx = xax->GetNbins(); printf ("nx = %i",nx);
01207 Axis_t xmin = xax->GetXmin(); printf (" xmin = %f",xmin);
01208 Axis_t xmax = xax->GetXmax(); printf (" xmax = %f\n",xmax);
01209 TAxis *yax = hist->GetYaxis();
01210 Int_t dim = hist->GetDimension();
01211 Int_t ny = yax->GetNbins();
01212 if (dim < 3) ny = 1; printf ("ny = %i",ny);
01213 Axis_t ymin = yax->GetXmin(); printf (" ymin = %f",ymin);
01214 Axis_t ymax = yax->GetXmax(); printf (" ymax = %f\n",ymax);
01215 struct Fit_t {
01216 Float_t i;
01217 Float_t j;
01218 Float_t x;
01219 Float_t y;
01220 Float_t mean;
01221 Float_t rms;
01222 Float_t peak;
01223 Float_t mu;
01224 Float_t sigma;
01225 Float_t entries;
01226 Float_t chisq;
01227 Float_t prob;
01228 Float_t a0;
01229 Float_t a1;
01230 Float_t a2;
01231 Float_t a3;
01232 Float_t a4;
01233 Float_t a5;
01234 Float_t Npar;
01235 Float_t dpeak;
01236 Float_t dmu;
01237 Float_t dsigma;
01238 Float_t da0;
01239 Float_t da1;
01240 Float_t da2;
01241 Float_t da3;
01242 Float_t da4;
01243 Float_t da5;
01244 };
01245 Fit_t Fit;
01246 TString NewRootFile(gSystem->DirName(fRootFile->GetName()));
01247 NewRootFile += "/";
01248 NewRootFile += HistName;
01249 NewRootFile += FitName;
01250 if (mergeX != 1) NewRootFile += Form("_x%i",mergeX);
01251 if (mergeY != 1) NewRootFile += Form("_y%i",mergeY);
01252
01253 NewRootFile += gSystem->BaseName(fRootFile->GetName());
01254 TFile *f = 0;
01255 TNtuple *FitP = 0;
01256 if (ix == 0 && jy == 0) {
01257 f = new TFile(NewRootFile.Data(),"update");
01258 if (f) cout << NewRootFile << " has been opened." << endl;
01259 else {cout << "Failed to open " << NewRootFile << endl; return;}
01260
01261
01262 #if 0
01263 if (TString(FitName) == "GF")
01264 FitP = new TNtuple("FitP","Fit results",
01265 "i:j:x:y:mean:rms:peak:mu:sigma:entries:chisq:prob:pi:P:K:e:d:a5:Npar:dpeak:dmu:dsigma:dpi:dP:dK:de:dd:da5");
01266 else if (TString(FitName) == "GB")
01267 FitP = new TNtuple("FitP","Fit results",
01268 "i:j:x:y:mean:rms:peak:mu:sigma:entries:chisq:prob:pi:P:K:e:d:dX:Npar:dpeak:dmu:dsigma:dpi:dP:dK:de:dd:ddX");
01269 else
01270 #endif
01271 FitP = new TNtuple("FitP","Fit results",
01272 "i:j:x:y:mean:rms:peak:mu:sigma:entries:chisq:prob:a0:a1:a2:a3:a4:a5:Npar:dpeak:dmu:dsigma:da0:da1:da2:da3:da4:da5");
01273 FitP->SetMarkerStyle(20);
01274 FitP->SetLineWidth(2);
01275 }
01276 TH1 *mean;
01277 TH1 *rms;
01278 TH1 *entries;
01279 TH1 *mu;
01280 TH1 *sigma;
01281 TH1 *chisq;
01282 if (dim == 3) {
01283 mean = new TH2D("mean",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01284 rms = new TH2D("rms",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01285 entries = new TH2D("entries",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01286 mu = new TH2D("mu",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01287 sigma = new TH2D("sigma",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01288 chisq = new TH2D("chisq",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
01289 }
01290 else {
01291 if (dim == 2 || dim == 1) {
01292 mean = new TH1D("mean",hist->GetTitle(),nx,xmin,xmax);
01293 rms = new TH1D("rms",hist->GetTitle(),nx,xmin,xmax);
01294 entries = new TH1D("entries",hist->GetTitle(),nx,xmin,xmax);
01295 mu = new TH1D("mu",hist->GetTitle(),nx,xmin,xmax);
01296 sigma = new TH1D("sigma",hist->GetTitle(),nx,xmin,xmax);
01297 chisq = new TH1D("chisq",hist->GetTitle(),nx,xmin,xmax);
01298 }
01299 else {
01300 printf("Histogram %s has wrong dimension %i\n", hist->GetName(),dim);
01301 return;
01302 }
01303 }
01304 Double_t params[9];
01305 TH1 *proj = 0;
01306 TF1 *g = 0;
01307 Int_t ix1 = ix, jy1 = jy;
01308 if (ix > 0) nx = ix;
01309 if (jy > 0) ny = jy;
01310 if (ix1 <= 0) ix1 = 0;
01311 if (jy1 <= 0) jy1 = 0;
01312 for (int i=ix1;i<=nx-mergeX+1;i++){
01313 Int_t ir0 = i;
01314 Int_t ir1=i+mergeX-1;
01315 if (i == 0) {ir0 = 1; ir1 = nx;}
01316 for (int j=jy1;j<=ny-mergeY+1;j++){
01317 Int_t jr0 = j;
01318 Int_t jr1 = j+mergeY-1;
01319 if (j == 0) {jr0 = 1; jr1 = ny;}
01320 if (dim == 3) {
01321 if (ir0 == ir1 && jr0 == jr1)
01322 proj = ((TH3 *) hist)->ProjectionZ(Form("f%i_%i", ir0, jr0 ),ir0,ir1,jr0,jr1);
01323 else
01324 proj = ((TH3 *) hist)->ProjectionZ(Form("f%i_%i_%i_%i",ir0,ir1,jr0,jr1),ir0,ir1,jr0,jr1);
01325 if (! proj) continue;
01326 TString title(proj->GetTitle());
01327 title += Form(" in x [%5.1f,%5.1f] and y [%5.1f,%5.1f] range",
01328 xax->GetBinLowEdge(ir0), xax->GetBinUpEdge(ir1),
01329 yax->GetBinLowEdge(jr0), yax->GetBinUpEdge(jr1));
01330 proj->SetTitle(title.Data());
01331 }
01332 else {
01333 if (ir0 == ir1)
01334 proj = ((TH2 *) hist)->ProjectionY(Form("f%i", ir0), ir0,ir1);
01335 else
01336 proj = ((TH2 *) hist)->ProjectionY(Form("f%i_%i",ir0,ir1),ir0,ir1);
01337 if (! proj) continue;
01338 TString title(proj->GetTitle());
01339 title += Form("in x [%5.1f,%5.1f] range",
01340 xax->GetBinLowEdge(ir0), xax->GetBinUpEdge(ir1));
01341 proj->SetTitle(title.Data());
01342 }
01343 cout << "i/j\t" << i << "/" << j << "\t" << proj->GetName() << "\t" << proj->GetTitle() << "\t" << proj->Integral() << endl;
01344
01345 memset (&Fit, 0, sizeof(Fit_t));
01346 Fit.i = (2.*i+mergeX-1.)/2;
01347 Fit.j = (2.*j+mergeY-1.)/2;
01348 Fit.x = 0.5*(xax->GetBinLowEdge(i) + xax->GetBinUpEdge(i+mergeX-1));
01349 Fit.y = 0.5*(yax->GetBinLowEdge(j) + yax->GetBinUpEdge(j+mergeY-1));
01350 Fit.mean = proj->GetMean();
01351 Fit.rms = proj->GetRMS();
01352 Fit.chisq = -100;
01353 Fit.prob = 0;
01354 Fit.entries = proj->Integral();
01355 if (Fit.entries < 100) {delete proj; continue;}
01356 if (TString(FitName) == "GP") g = FitGP(proj,opt,nSigma,pow);
01357 else if (TString(FitName) == "G2") g = FitG2(proj,opt);
01358 else if (TString(FitName) == "GF") g = FitGF(proj,opt);
01359 else if (TString(FitName) == "GB") {
01360
01361 Double_t dX = 2.364;
01362 if ((nx == 48 && (i > 0 && i <= 24)) ||
01363 (nx == 24 && (j > 0 && j <= 13)) ||
01364 (nx == 45 && (i > 0 && i <= 13)))
01365 dX = 1.372;
01366 g = FitGB(proj,opt,dX);
01367 }
01368 else {cout << FitName << " has not been definded" << endl; break;}
01369 if (! g ) {delete proj; continue;}
01370 g->GetParameters(params);
01371 Fit.Npar = g->GetNpar();
01372 Fit.chisq = g->GetChisquare();
01373 Fit.prob = g->GetProb();
01374 Fit.peak = params[0];
01375 Fit.mu = params[1];
01376 Fit.sigma = TMath::Abs(params[2]);
01377 Fit.a0 = params[3];
01378 Fit.a1 = params[4];
01379 Fit.a2 = params[5];
01380 Fit.a3 = params[6];
01381 Fit.a4 = params[7];
01382 Fit.a5 = params[8];
01383 Fit.dpeak = g->GetParError(0);
01384 Fit.dmu = g->GetParError(1);
01385 Fit.dsigma = g->GetParError(2);
01386 Fit.da0 = g->GetParError(3);
01387 Fit.da1 = g->GetParError(4);
01388 Fit.da2 = g->GetParError(5);
01389 Fit.da3 = g->GetParError(6);
01390 Fit.da4 = g->GetParError(7);
01391 Fit.da5 = g->GetParError(8);
01392
01393
01394 if (dim == 3) {
01395 mean->SetBinContent(i,j,Fit.mean);
01396 rms->SetBinContent(i,j,Fit.rms);
01397 entries->SetBinContent(i,j,Fit.entries);
01398 mu->SetBinContent(i,j,Fit.mu);
01399 mu->SetBinError(i,j,g->GetParError(1));
01400 sigma->SetBinContent(i,j,Fit.sigma);
01401 sigma->SetBinError(i,j,g->GetParError(2));
01402 chisq->SetBinContent(i,j,Fit.chisq);
01403 }
01404 else {
01405 mean->SetBinContent(i,Fit.mean);
01406 rms->SetBinContent(i,Fit.rms);
01407 entries->SetBinContent(i,Fit.entries);
01408 mu->SetBinContent(i,Fit.mu);
01409 mu->SetBinError(i,g->GetParError(1));
01410 sigma->SetBinContent(i,Fit.sigma);
01411 sigma->SetBinError(i,g->GetParError(2));
01412 chisq->SetBinContent(i,Fit.chisq);
01413 }
01414
01415 printf("%i/%i %f/%f mean %f rms = %f entries = %f mu = %f sigma = %f chisq = %f prob = %f\n",
01416 i,j,Fit.x,Fit.y,Fit.mean,Fit.rms,Fit.entries,Fit.mu,Fit.sigma,Fit.chisq,Fit.prob);
01417 if (FitP) FitP->Fill(&Fit.i);
01418 if (canvas) {
01419 canvas->Update();
01420 #if 0
01421 int ii;
01422 cout << "Type i ";
01423 cin >> ii; if (ii <0) return;
01424 #endif
01425 }
01426 if (! canvas) delete proj;
01427
01428 }
01429 }
01430 if (f) {
01431 f->Write();
01432 delete f;
01433 }
01434 }
01435
01436 Int_t FitG(TH1 *proj, TF1 *g, TF1 *ga){
01437 Double_t params[10];
01438 proj->Fit("g","R");
01439 Double_t chisq = g->GetChisquare();
01440 if (chisq <= 0. || chisq > 1.e10) return -1;
01441 g->GetParameters(params);
01442 params[3] = 0;
01443 params[4] = 0;
01444 params[5] = 0;
01445 params[6] = 0;
01446 params[7] = 0;
01447 params[8] = 0;
01448 params[9] = 0;
01449 ga->SetParameters(params);
01450 proj->Fit("ga","r");
01451
01452 chisq = g->GetChisquare();
01453 if (chisq <= 0. || chisq > 1.e10) return -1;
01454 return 0;
01455 }
01456
01457 Int_t FitGG(TH1 *proj, TF1 *g1, TF1 *g2=0, TF1 *ga2=0, Double_t scaleM=-2., Double_t scaleP=2.) {
01458 Double_t params[9];
01459 proj->Fit("g1","R");
01460 g1->GetParameters(params);
01461 Double_t chisq = g1->GetChisquare();
01462 if (chisq <= 0. || chisq > 1.e10) return -1;
01463 params[3] = 1;
01464 params[4] = params[1]+0.1;
01465 params[5] = params[2];
01466 g2->SetParameters(params);
01467 proj->Fit("g2");
01468 chisq = g2->GetChisquare();
01469 if (chisq <= 0. || chisq > 1.e10) return -1;
01470 g2->GetParameters(params);
01471 if (ga2) {
01472 params[6] = 0;
01473 params[7] = 0;
01474 params[8] = 0;
01475 ga2->SetParameters(params);
01476 proj->Fit("ga2","R");
01477 ga2->GetParameters(params);
01478 ga2->SetRange(params[1]+scaleM*params[2],params[1]+scaleP*params[2]);
01479 proj->Fit("ga2","R");
01480 }
01481
01482
01483
01484
01485 return 0;
01486 }
01487
01488 void FitX(TH2 *hist=0, Double_t range=1, Int_t Ibin = 0) {
01489 TFile *fRootFile = (TFile *) gDirectory->GetFile();
01490 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
01491 if (!hist) {
01492 const Char_t *HistName = "FPoints";
01493 hist = (TH2D *) fRootFile->Get(HistName);
01494 if (!hist) {printf("Cannot histogram %s\n",HistName); return;}
01495 }
01496
01497 TF1 *g = new TF1("g","gaus",-range,range);
01498
01499 TString fitN("gaus(0)+exp(pol3(3))");
01500 if (range > 2.) fitN = "gaus(0)+exp(pol5(3))";
01501 TF1 *ga= new TF1("ga",fitN.Data(),-range,range);
01502 ga->SetParName(0,"Area Pi");
01503 ga->SetParName(1,"Mean Pi");
01504 ga->SetParName(2,"Sigma Pi");
01505 ga->SetParName(3,"A0");
01506 ga->SetParName(4,"A1");
01507 ga->SetParName(5,"A2");
01508 ga->SetParName(6,"A3");
01509 ga->SetParName(7,"A4");
01510 ga->SetParName(8,"A5");
01511
01512 TCanvas *c = new TCanvas("Fit","Fit results");
01513 int i;
01514 TString name(hist->GetName());
01515 name += "MuFG";
01516 Int_t nBins = hist->GetXaxis()->GetNbins();
01517 Double_t xlow = hist->GetXaxis()->GetXmin();
01518 Double_t xup = hist->GetXaxis()->GetXmax();
01519 TH1D *MuF = new TH1D(name.Data(),"Avarage shift versus no. of measurement points",
01520 nBins, xlow, xup);
01521 name = hist->GetName();
01522 name += "SigmaFG";
01523 TH1D *SigmaF = new TH1D(name.Data(),"Sigma of z versus no. of measurement points",
01524 nBins, xlow, xup);
01525 Int_t i1 = 1, i2 = nBins;
01526 if (Ibin > 0 && Ibin <= nBins) {i1 = Ibin; i2 = Ibin;}
01527 Double_t XFitP, dXFitP, MuFitP,dMuFitP,SigmaFitP,dSigmaFitP;
01528 TH1 *proj = 0;
01529 for (i=i1; i<=i2; i++) {
01530 if (proj) delete proj;
01531 proj = hist->ProjectionY("proj",i,i);
01532 XFitP = hist->GetXaxis()->GetBinCenter(i);
01533 dXFitP = hist->GetXaxis()->GetBinWidth(i);
01534 Double_t chisq = -999;
01535 if (proj->Integral() < 1000) {
01536 if (N>0) goto NEXT;
01537 continue;
01538 }
01539 if (FitG(proj,g,ga)) goto NEXT;
01540
01541 if (ga->GetParError(1) <= 0 || ga->GetParError(1) > 0.01 ) {
01542 printf("============= REJECT ================\n");
01543 goto NEXT;
01544 }
01545 MuFitP = ga->GetParameter(1);
01546 dMuFitP = ga->GetParError(1);
01547
01548 SigmaFitP = ga->GetParameter(2);
01549 dSigmaFitP = ga->GetParError(2);
01550 chisq = ga->GetChisquare();
01551 if (chisq > 1.e4) goto NEXT;
01552 MuF->SetCellContent(i,0,MuFitP);
01553 MuF->SetCellError(i,0,dMuFitP);
01554 SigmaF->SetCellContent(i,0,SigmaFitP);
01555 SigmaF->SetCellError(i,0,dSigmaFitP);
01556 proj->Draw(); c->Update();
01557 printf("Bin: %i x: %f +/- %f MuF: %f+/-%f Sigma: %f+/-%f chisq: %f \n",
01558 i,XFitP,dXFitP,MuFitP,dMuFitP,SigmaFitP,dSigmaFitP,chisq);
01559 NEXT:
01560 N++;
01561
01562 }
01563 if (! Ibin ) {
01564 c->Clear();
01565
01566
01567
01568
01569 SigmaF->SetMarkerStyle(20);
01570 TF1 *p = new TF1("p","[0]/pow(x,[1])");
01571 p->SetParameter(0,0.64);
01572 p->SetParameter(1,0.50);
01573 SigmaF->Fit("p","R");
01574 SigmaF->SetMinimum(0.06);
01575 SigmaF->Draw();
01576 }
01577 TString NewRootFile(hist->GetName());
01578 NewRootFile += fRootFile->GetName();
01579 TFile *f = new TFile(NewRootFile.Data(),"update");
01580 MuF->Write();
01581 SigmaF->Write();
01582 delete f;
01583 fRootFile->cd();
01584 }
01585 #if 0
01586
01587 void FitC4G(const Char_t *set="z") {
01588 Double_t sigmas[2] = {0.06,0.12};
01589 const Int_t nHYPS = NHYPS;
01590 TH2 *hists[nHYPS];
01591 TProfile *histp[nHYPS];
01592 TFile *fRootFile = (TFile *) gDirectory->GetFile();
01593 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
01594 TString newfile("FitC4G");
01595 newfile += "set";
01596 newfile += gSystem->BaseName(fRootFile->GetName());
01597 TFile * f = new TFile(newfile.Data(),"update");
01598 for (Int_t i = 0; i< nHYPS; i++) {
01599 Char_t *HistN = HistNames[i];
01600 if (TString(set) == "70") HistN = HistNames70[i];
01601 fRootFile->Get(HistN);
01602 if (!hists[i]) {printf("Cannot histogram %s\n",HistNames[i]); return;}
01603 histp[i] = (TProfile *) fRootFile->Get(HistNameP[i]);
01604 if (!histp[i]) {printf("Cannot histogram %s\n",HistNameP[i]); return;}
01605 }
01606 TF1 *g[4];
01607 TF1 *g1 = new TF1("g1","gaus",-0.3,0.3);
01608 TF1 *g = 0, *g2 = 0, *g3 = 0, *g4 = 0, *ga = 0;
01609 TCanvas *canvas = new TCanvas("canvas","canvas");
01610 if (Mu2 != 0) {
01611 if (Mu2 < -1.) {
01612 if (Mu2 <= -4.) ga = new TF1("ga","gaus(0)+exp(pol3(3))",xmin2,xmax2);
01613 else {
01614 if (Mu2 <= -3.) ga = new TF1("ga","gaus(0)+exp(pol2(3))",xmin2,xmax2);
01615 else ga = new TF1("ga","gaus(0)+exp(pol1(3))",xmin2,xmax2);
01616 }
01617 g = ga;
01618 g->SetParName(0,"Constant1");
01619 g->SetParName(1,"Mean1");
01620 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
01621 g->SetParName(3,"Const");
01622 g->SetParName(4,"Slope1");
01623 g->SetParName(5,"Slope2");
01624 g->SetParName(6,"Slope3");
01625 g->SetParName(7,"Slope4");
01626 }
01627 else {
01628 g2 = new TF1("g2","gaus(0)+gaus(3)",xmin2,xmax2);
01629 g = g2;
01630 g->SetParName(0,"Constant1");
01631 g->SetParName(1,"Mean1");
01632 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
01633 g->SetParName(3,"Constant2");
01634 g->SetParName(4,"Mean2");
01635 g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
01636 if (Mu3 != 0) {
01637 g3= new TF1("g3","gaus(0)+gaus(3)+gaus(6)",xmin3,xmax3);
01638 g = g3;
01639 g->SetParName(0,"Constant1");
01640 g->SetParName(1,"Mean1");
01641 g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
01642 g->SetParName(3,"Constant2");
01643 g->SetParName(4,"Mean2");
01644 g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
01645 g->SetParName(6,"Constant3");
01646 g->SetParName(7,"Mean3");
01647 g->SetParName(8,"Sigma3"); g->SetParLimits(8,sigmas[0],sigmas[1]);
01648 if (Mu4 != 0) {
01649 g4= new TF1("g4","gaus(0)+gaus(3)+gaus(6)+gaus(9)",xmin3,xmax3);
01650 g4->SetParName(0,"Constant1");
01651 g4->SetParName(1,"Mean1");
01652 g4->SetParName(2,"Sigma1"); g4->SetParLimits(2,sigmas[0],sigmas[1]);
01653 g4->SetParName(3,"Constant2");
01654 g4->SetParName(4,"Mean2");
01655 g4->SetParName(5,"Sigma2"); g4->SetParLimits(5,sigmas[0],sigmas[1]);
01656 g4->SetParName(6,"Constant3");
01657 g4->SetParName(7,"Mean3");
01658 g4->SetParName(8,"Sigma3"); g4->SetParLimits(8,sigmas[0],sigmas[1]);
01659 g4->SetParName(9,"Constant4");
01660 g4->SetParName(10,"Mean4");
01661 g4->SetParName(11,"Sigma4"); g4->SetParLimits(11,sigmas[0],sigmas[1]);
01662 }
01663 }
01664 }
01665 }
01666 TH1 *proj = 0;
01667 Double_t params[12];
01668 Int_t k;
01669 Int_t Bin = TMath::Abs(bin);
01670 for (k = 0; k<nHYPS; k++) {
01671 if (hyp != -1 && k != hyp) continue;
01672 TH2 *hist = hists[k];
01673 Int_t nx = hist->GetNbinsX();
01674 for (i=1; i<=nx; i++) {
01675 if (bin != 0 && Bin != i) continue;
01676
01677 Char_t line[40];
01678 sprintf(line,"%s_%i",hist->GetName(),i);
01679 TString name(line);
01680 proj = hist->ProjectionY(name.Data(),i,i);
01681 Int_t ix1=proj->GetXaxis()->FindBin(-.1);
01682 Int_t ix2=proj->GetXaxis()->FindBin(0.1);
01683 if (proj->Integral(ix1,ix2) < 100.) {
01684 printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hist->GetName(),Bin,k,proj->Integral());
01685 delete proj;
01686 continue;
01687 }
01688 Int_t NFitPoints = 0;
01689 Double_t chisq;
01690 g = g1;
01691 g->SetParLimits(0,0,1.e7);
01692 g->SetParLimits(1,-.1,.1);
01693 g->SetParLimits(2,sigmas[0],sigmas[1]);
01694 proj->Fit(g->GetName(),"ri");
01695 g->GetParameters(params);
01696 Int_t kcase = -1;
01697 if (xmin1 < -0.5 && xmax1 > 0.5 && g->GetProb() > 1.e-3
01698 && TMath::Abs(g->GetParameter(1)) < 0.05) goto Done;
01699 if (ga) {
01700 g = ga;
01701 g1->GetParameters(params);
01702 params[1] = 0.;
01703 params[3] = 0.;
01704 params[4] = 0.;
01705 params[5] = 0.;
01706 params[6] = 0.;
01707 params[7] = 0.;
01708 params[8] = 0.;
01709 g->SetParameters(params);
01710 kcase = 0;
01711 if (g->GetProb() > 1.e-3) goto Done;
01712 }
01713 else {
01714 if (g2) {
01715 params[1] = 0.;
01716 params[3] = params[0];
01717 params[4] = Mu2;
01718 params[5] = 2.0*params[2];
01719 g = g2;
01720 g->SetParameters(params);
01721 g->SetParLimits(0,0,1.e7);
01722 g->SetParLimits(1,-.1,.1);
01723 g->SetParLimits(2,sigmas[0],sigmas[1]);
01724 g->SetParLimits(3,0,1.e7);
01725 g->SetParLimits(5,sigmas[0],sigmas[1]);
01726 proj->Fit(g->GetName(),"ri");
01727 g->GetParameters(params);
01728 kcase = 2;
01729 if (g->GetProb() > 1.e-3) goto Done;
01730 if (g3) {
01731 params[1] = 0.;
01732 params[6] = params[0];
01733 params[7] = Mu3;
01734 params[8] = 2.0*params[2];
01735 g = g3;
01736 g->SetParameters(params);
01737 g->SetParLimits(0,0,1.e7);
01738 g->SetParLimits(1,-.1,.1);
01739 g->SetParLimits(2,sigmas[0],sigmas[1]);
01740 g->SetParLimits(3,0,1.e7);
01741 g->SetParLimits(5,sigmas[0],sigmas[1]);
01742 g->SetParLimits(6,0,1.e7);
01743 g->SetParLimits(8,sigmas[0],sigmas[1]);
01744 proj->Fit(g->GetName(),"ri");
01745 g->GetParameters(params);
01746 kcase = 3;
01747 if (g->GetProb() > 1.e-3) goto Done;
01748 if (g4) {
01749 params[1] = 0.;
01750 params[9] = params[0];
01751 params[10] = Mu4;
01752 params[11] = 2.0*params[2];
01753 g = g4;
01754 g->SetParameters(params);
01755 g->SetParLimits(0,0,1.e7);
01756 g->SetParLimits(1,-.1,.1);
01757 g->SetParLimits(2,sigmas[0],sigmas[1]);
01758 g->SetParLimits(3,0,1.e7);
01759 g->SetParLimits(5,sigmas[0],sigmas[1]);
01760 g->SetParLimits(6,0,1.e7);
01761 g->SetParLimits(8,sigmas[0],sigmas[1]);
01762 g->SetParLimits(9,0,1.e7);
01763 g->SetParLimits(11,sigmas[0],sigmas[1]);
01764 proj->Fit(g->GetName(),"ri");
01765 kcase = 4;
01766 g->GetParameters(params);
01767 }
01768 }
01769 }
01770 }
01771 Done:
01772 if (g) {
01773 proj->Fit(g->GetName(),"RIM");
01774 canvas->Update();
01775 Int_t l = 1;
01776 Double_t mu = g->GetParameter(l);
01777 for (int m=2;m<=kcase;m++) {
01778 if (TMath::Abs(mu) > TMath::Abs(g->GetParameter(3*m-2))) {
01779 l = 3*m-2; mu = g->GetParameter(l);
01780
01781 }
01782 }
01783 Nu[N] = g->GetParameter(l);
01784 Mu[N] = Nu[N] + histp[k]->GetBinContent(i);
01785 dMu[N] = g->GetParError(l);
01786
01787 Sigma[N] = g->GetParameter(2);
01788 dSigma[N] = g->GetParError(2);
01789 NFitPoints = g->GetNumberFitPoints();
01790 Int_t NDF = g->GetNDF();
01791 Double_t prob = g->GetProb();
01792 chisq = g->GetChisquare();
01793 X[N] = hist->GetXaxis()->GetBinCenter(i);
01794 dX[N] = hist->GetXaxis()->GetBinWidth(i);
01795 Double_t pionM = Masses[k]*pow(10.,X[N]);
01796 printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
01797 "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
01798 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob);
01799
01800
01801 if ( bin >= 0 && TMath::Abs(Nu[N]) < 0.10 || bin < 0) {
01802 printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
01803 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
01804 TString FileN("FitPars");
01805 if (hyp > -1) FileN += HistNames[hyp];
01806 FileN += ".h";
01807 FILE *fp = fopen(FileN.Data(),"a");
01808 if (fp) {
01809 if (N == 0) {
01810 TDatime time;
01811 fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
01812 fprintf(fp,
01813 "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
01814 }
01815 fprintf(fp,
01816 "{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
01817 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
01818 fclose(fp);
01819 }
01820 #if 1
01821 proj->Write();
01822 #endif
01823 }
01824 else printf ("================== Skip it\n");
01825 }
01826 N++;
01827 proj->Draw();
01828 }
01829 }
01830 if (f) {delete f;}
01831 }
01832 #endif
01833
01834 void Fit4G(Int_t ng=2, Int_t hyp=-1, Int_t bin=0,
01835 Double_t xmin1=-1.0, Double_t xmax1 = 1.0,
01836 Double_t Mu2 = -.5,Double_t xmin2=-1., Double_t xmax2 = 1.,
01837 Double_t Mu3 = 0.0,Double_t xmin3=-1., Double_t xmax3 = 1.,
01838 Double_t Mu4 = 0.0){
01839 Double_t sigmas[2] = {0.06,0.12};
01840 const Int_t nHYPS = 4;
01841 TH2 *hists[nHYPS];
01842 TProfile *histp[nHYPS];
01843 TFile *fRootFile = (TFile *) gDirectory->GetFile();
01844 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
01845 TFile *f = 0;
01846 if (bin == 0) {
01847 TString newfile("BBFit");
01848 newfile += gSystem->BaseName(fRootFile->GetName());
01849 f = new TFile(newfile.Data(),"update");
01850 }
01851 for (int i = 0; i< nHYPS; i++) {
01852 if (hyp >= 0 && hyp != i) continue;
01853 hists[i] = (TH2 *) fRootFile->Get(HistNames[i]);
01854 if (!hists[i]) {printf("Cannot histogram %s\n",HistNames[i]); return;}
01855 histp[i] = (TProfile *) fRootFile->Get(HistNameP[i]);
01856 if (!histp[i]) {printf("Cannot histogram %s\n",HistNameP[i]); return;}
01857 }
01858 TF1 *g1 = new TF1("g1","gaus",xmin1,xmax1);
01859 TF1 *g = 0, *g2 = 0, *g3 = 0, *g4 = 0, *ga = 0;
01860 TCanvas *canvas = new TCanvas("canvas","canvas");
01861 if (ng<0) {
01862 ga = new TF1("ga",Form("gaus(0)+exp(pol%i(3))",-ng),xmin2,xmax2);
01863 ga->SetParName(0,"Constant1");
01864 ga->SetParName(1,"Mean1");
01865 ga->SetParName(2,"Sigma1");
01866 ga->SetParName(3,"Const");
01867 for (int m = 0; m <= -ng; m++) ga->SetParName(m+3,Form("Slope%i",m));
01868 }
01869 else {
01870 for (int m = 2; m <= ng; m++) {
01871 if (m == 2) {g2 = new TF1("g2","gaus(0)+gaus(3)",xmin2,xmax2); g = g2;}
01872 if (m == 3) {g3 = new TF1("g3","gaus(0)+gaus(3)+gaus(6)",xmin3,xmax3); g = g3;}
01873 if (m == 4) {g4 = new TF1("g4","gaus(0)+gaus(3)+gaus(6)+gaus(9)",-1.,1.); g = g4;}
01874 for (int k = 0; k < 3*m; k += 3) {
01875 g->SetParName(k ,Form("Constant%i",m));
01876 g->SetParName(k+1,Form("Mean%i",m));
01877 g->SetParName(k+2,Form("Sigma%i",m));
01878 }
01879 }
01880 }
01881 TH1 *proj = 0;
01882 Double_t params[12];
01883 Int_t k,i;
01884 Int_t Bin = TMath::Abs(bin);
01885 for (k = 0; k<nHYPS; k++) {
01886 if (hyp != -1 && k != hyp) continue;
01887 TH2 *hist = hists[k];
01888 Int_t nx = hist->GetNbinsX();
01889 for (i=1; i<=nx; i++) {
01890 if (bin != 0 && Bin != i) continue;
01891
01892 Char_t line[40];
01893 sprintf(line,"%s_%i",hist->GetName(),i);
01894 TString name(line);
01895 proj = hist->ProjectionY(name.Data(),i,i);
01896 Int_t ix1=proj->GetXaxis()->FindBin(-.1);
01897 Int_t ix2=proj->GetXaxis()->FindBin(0.1);
01898 if (proj->Integral(ix1,ix2) < 100) {
01899 printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hist->GetName(),Bin,k,proj->Integral());
01900 delete proj;
01901 continue;
01902 }
01903 Int_t NFitPoints = 0;
01904 Double_t xrange = 0.6;
01905 Double_t chisq;
01906 g = g1;
01907 g->SetParLimits(0,0,1.e7);
01908 g->SetParLimits(1,-xrange,xrange);
01909
01910 g->SetParLimits(2,sigmas[0],sigmas[1]);
01911 proj->Fit(g->GetName(),"ri");
01912 g->GetParameters(params);
01913 Int_t kcase = -1;
01914 if (xmin1 < -0.5 && xmax1 > 0.5 && g->GetProb() > 1.e-7
01915 && TMath::Abs(g->GetParameter(1)) < 0.05) goto Done;
01916 if (ga) {
01917 g = ga;
01918 g1->GetParameters(params);
01919 params[1] = 0.;
01920 params[3] = 0.;
01921 params[4] = 0.;
01922 params[5] = 0.;
01923 params[6] = 0.;
01924 params[7] = 0.;
01925 params[8] = 0.;
01926 g->SetParameters(params);
01927 kcase = 0;
01928 if (g->GetProb() > 1.e-7) goto Done;
01929 }
01930 else {
01931 if (g2) {
01932 params[1] = 0.;
01933 params[3] = params[0];
01934 params[4] = Mu2;
01935 params[5] = 2.0*params[2];
01936 g = g2;
01937 g->SetParameters(params);
01938 g->SetParLimits(0,0,1.e7);
01939 g->SetParLimits(1,-xrange,xrange);
01940
01941 g->SetParLimits(2,sigmas[0],sigmas[1]);
01942 g->SetParLimits(3,0,1.e7);
01943 g->SetParLimits(5,sigmas[0],sigmas[1]);
01944 proj->Fit(g->GetName(),"ri");
01945 g->GetParameters(params);
01946 kcase = 2;
01947 if (g->GetProb() > 1.e-7) goto Done;
01948 if (g3) {
01949 params[1] = 0.;
01950 params[6] = params[0];
01951 params[7] = Mu3;
01952 params[8] = 2.0*params[2];
01953 g = g3;
01954 g->SetParameters(params);
01955 g->SetParLimits(0,0,1.e7);
01956 g->SetParLimits(1,-xrange,xrange);
01957
01958 g->SetParLimits(2,sigmas[0],sigmas[1]);
01959 g->SetParLimits(3,0,1.e7);
01960 g->SetParLimits(5,sigmas[0],sigmas[1]);
01961 g->SetParLimits(6,0,1.e7);
01962 g->SetParLimits(8,sigmas[0],sigmas[1]);
01963 proj->Fit(g->GetName(),"ri");
01964 g->GetParameters(params);
01965 kcase = 3;
01966 if (g->GetProb() > 1.e-7) goto Done;
01967 if (g4) {
01968 params[1] = 0.;
01969 params[9] = params[0];
01970 params[10] = Mu4;
01971 params[11] = 2.0*params[2];
01972 g = g4;
01973 g->SetParameters(params);
01974 g->SetParLimits(0,0,1.e7);
01975 g->SetParLimits(1,-xrange,xrange);
01976
01977 g->SetParLimits(2,sigmas[0],sigmas[1]);
01978 g->SetParLimits(3,0,1.e7);
01979 g->SetParLimits(5,sigmas[0],sigmas[1]);
01980 g->SetParLimits(6,0,1.e7);
01981 g->SetParLimits(8,sigmas[0],sigmas[1]);
01982 g->SetParLimits(9,0,1.e7);
01983 g->SetParLimits(11,sigmas[0],sigmas[1]);
01984 proj->Fit(g->GetName(),"ri");
01985 kcase = 4;
01986 g->GetParameters(params);
01987 }
01988 }
01989 }
01990 }
01991 Done:
01992 if (g) {
01993 proj->Fit(g->GetName(),"RIM");
01994 canvas->Update();
01995 Int_t l = 1;
01996 Double_t mu = g->GetParameter(l);
01997 for (int m=2;m<=kcase;m++) {
01998 if (TMath::Abs(mu) > TMath::Abs(g->GetParameter(3*m-2))) {
01999 l = 3*m-2; mu = g->GetParameter(l);
02000
02001 }
02002 }
02003 Nu[N] = g->GetParameter(l);
02004 Mu[N] = Nu[N] + histp[k]->GetBinContent(i);
02005 dMu[N] = g->GetParError(l);
02006
02007 Sigma[N] = g->GetParameter(2);
02008 dSigma[N] = g->GetParError(2);
02009 NFitPoints = g->GetNumberFitPoints();
02010 Int_t NDF = g->GetNDF();
02011 Double_t prob = g->GetProb();
02012 chisq = g->GetChisquare();
02013 X[N] = hist->GetXaxis()->GetBinCenter(i);
02014 dX[N] = hist->GetXaxis()->GetBinWidth(i);
02015 Double_t pionM = Masses[k]*pow(10.,X[N]);
02016 printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
02017 "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
02018 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob);
02019
02020
02021 if ( (bin >= 0 && TMath::Abs(Nu[N]) < 0.10) || bin < 0) {
02022 printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
02023 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
02024 TString FileN("FitPars");
02025 if (hyp > -1) FileN += HistNames[hyp];
02026 FileN += ".h";
02027 FILE *fp = fopen(FileN.Data(),"a");
02028 if (fp) {
02029 if (N == 0) {
02030 TDatime time;
02031 fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
02032 fprintf(fp,
02033 "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
02034 }
02035 fprintf(fp,
02036 "{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
02037 Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
02038 fclose(fp);
02039 }
02040 #if 1
02041 proj->Write();
02042 #endif
02043 }
02044 else printf ("================== Skip it\n");
02045 }
02046 N++;
02047 proj->Draw();
02048 }
02049 }
02050 if (f) {delete f;}
02051 }
02052 #ifdef PLUSHIKIN
02053
02054 void DrawMonths(TH1 *TimeMuFG,Double_t ymin=-.5, Double_t ymax=2.5) {
02055
02056 TimeMuFG->SetMaximum(ymax);
02057 TimeMuFG->SetMinimum(ymin);
02058 TimeMuFG->SetStats(0);
02059 TimeMuFG->Draw();
02060 TimeMuFG->SetMarkerStyle(20);
02061 TimeMuFG->SetXTitle("seconds since the Epoch");
02062 TimeMuFG->SetYTitle("z Fit");
02063 TimeMuFG->SetTitle("");
02064 TimeMuFG->GetXaxis()->SetTimeDisplay(1);
02065 TimeMuFG->GetXaxis()->SetTimeFormat("%m/%d/%y");
02066 TimeMuFG->Draw();
02067 for (Int_t i=7; i<13; i++) {
02068 Int_t d = 20010001 + 100*i;
02069 TDatime t(d,0);
02070 TString ts(t.AsSQLString());
02071 TString ymd(ts.Data(),ts.Index("-01 "));
02072 Double_t x = t.Convert();
02073 TLine *l = new TLine(x,ymin,x,ymax);
02074 l->SetLigneColor(i-6);
02075 l->Draw();
02076 TText *text = new TText(x, 0.8*(ymin+ymax),ymd.Data());
02077 text->Draw();
02078 }
02079 }
02080
02081 void dEdxSpC(TTree *DeDxTree) {
02082 const Int_t noSectors = 24;
02083 const Int_t noPadrows = 45;
02084 const Int_t noSigns = 2;
02085 Int_t nx = noSectors*noPadrows*noSigns;
02086 TCut Cut("m_flag0points>36 && m_dx>0 && m_dx<10 && m_de>0 && m_de<1.e-4 && m_pmag > 0.5");
02087 Char_t line[80];
02088 sprintf(line,"log(m_de/m_dx/m_zPi):(m_charge+1)/2+%i*(m_padrow-1+%i*(m_sector-1))",noSigns,noPadrows);
02089 TString Plot(line); printf("Plot = %s\n",Plot.Data());
02090 TString Title = Plot;
02091 Title += "versus Sector, Padrow, Charge for ";
02092 Title += Cut.GetTitle();
02093 TH2F *SpC = new TH2F("SpC",Title.Data(),nx,0,nx,200,-3.,7.);
02094 Plot += " >> SpC";
02095 DeDxTree->Draw(Plot.Data(),Cut);
02096 TFile *f = new TFile("SpC.root","RECREATE");
02097 SpC->Write();
02098 delete f;
02099 }
02100
02101 void plotProfSpC(TCanvas *c1,TProfile *prof)
02102 {
02103 if (c1 && prof) {
02104 c1->Clear();
02105 c1->Divide(4,6);
02106 char line[10];
02107 for (Int_t i=0;i<24;i++) {
02108 c1->cd(i+1);
02109 sprintf(line,"Sector %i",i+1);
02110 prof->SetTitle(line);
02111 Int_t i1 = 90*i;
02112 prof->SetAxisRange(i1,i1+89); prof->Draw("R");
02113 }
02114 }
02115 }
02116
02117 Double_t LogLandau(Double_t *x, Double_t *par){
02118 Double_t yy =x[0];
02119 Double_t xx = TMath::Exp(yy);
02120 Double_t f = xx*TMath::Landau(xx,par[0],par[1]);
02121 return f;
02122 }
02123
02124 void Func()
02125 {
02126 TF1 *f1 = new TF1("Func",LogLandau,0,10,2);
02127 f1->SetParameters(0,1);
02128 f1->SetParNames("mean","sigma");
02129 f1->Draw();
02130 }
02131
02132 void SpCfit(TH2F *hist, Int_t i1, Int_t i2)
02133 {
02134 TH1 *h1 = hist->ProjectionY("bin",i1,i2);
02135 TF1 *f1= (TF1 *) gROOT->GetFunction("Func");
02136 f1->SetParameters(0,1);
02137 h1->Fit("Func");
02138 }
02139
02140 void Fee(const Char_t *topDir = "/star/rcf/scratch/fisyak/NTuples2/",
02141 const Char_t *TreeName = "DeDxTree") {
02142 Int_t NoSector = 24;
02143 Int_t NoFee = 82;
02144 TH3F *Fee = new TH3F("Fee","dEdx versus sector and Fee",
02145 NoSector,1,NoSector+1,
02146 NoFee,0,NoFee,
02147 200,-5.,5.);
02148 if (gClassTable->GetID("TDataSet")<0) gSystem->Load("libTable");
02149 TFileSet dirs(topDir);
02150 TDataSetIter next(&dirs,0);
02151 TDataSet *set = 0;
02152 Int_t NFiles = 0;
02153 while ( (set = next()) ) {
02154 if (strcmp(set->GetTitle(),"file") ||
02155 !strstr(set->GetName(),".root")) continue;
02156 TString File(gSystem->ConcatFileName(topDir,set->Path()));
02157
02158 TString HName("F");
02159 HName += set->Path();
02160 HName.ReplaceAll("/","");
02161 HName.ReplaceAll(".root","");
02162 TFile f(File.Data());
02163 cout << " ================== Opened " << File.Data() << endl;
02164 NFiles++;
02165 TTree *DeDxTree = (TTree *) f.Get(TreeName);
02166 if (DeDxTree) {
02167 if (NFiles > 2) break;
02168 TH3F *fee = new TH3F(HName.Data(),"dEdx versus sector and Fee",
02169 NoSector,1,NoSector+1,
02170 NoFee,0,NoFee,
02171 200,-5.,5.);
02172 DeDxTree->Draw("log(m_de/m_dx/m_zPi):m_fee:m_sector>>Fee","m_de>0 && m_de<2e-4");
02173
02174
02175
02176
02177
02178
02179
02180
02181 Fee->Add(fee);
02182 TFile *ff = new TFile("fee.root","UPDATE");
02183 fee->Write();
02184 delete ff;
02185 delete fee;
02186 }
02187 else {
02188 cout << " -------------------- File is empty" << endl;
02189
02190 }
02191 Fee->Draw("colz");
02192 gPad->Update();
02193 TFile *ff = new TFile("fee.root","UPDATE");
02194 Fee->Write();
02195 delete ff;
02196 }
02197 }
02198
02199 void MakeTimeGain(TH1 *hist, const Char_t *TabNam = "TpcTimeGain"){
02200 char filename[80];
02201 Double_t params[2];
02202
02203
02204
02205 if (gClassTable->GetID("TTable") < 0) gSystem->Load("libTable");
02206 if (gClassTable->GetID("St_TpcTimeGain") < 0) gSystem->Load("St_Tables");
02207 params[0] = 0;
02208 St_TpcTimeGain *timegain = new St_TpcTimeGain(TabNam,1);
02209 timegain->SetNRows(1);
02210 TpcTimeGain_st *gain = timegain->GetTable();
02211 Int_t Nbins = hist->GetNbinsX();
02212 for (Int_t i=1; i<=Nbins; i++) {
02213 Double_t scale = hist->GetBinContent(i);
02214 Double_t error = hist->GetBinError(i);
02215 if (error == 0) continue;
02216 if (error > 0.1) continue;
02217
02218 if (TMath::Abs(scale) > 0.50) continue;
02219 gain->ScaleFactor = TMath::Exp(-(scale-params[0]));
02220 gain->ErrorScaleFactor = error;
02221 Double_t dt = hist->GetBinLowEdge(i);
02222 UInt_t date = dt;
02223 TDatime time; time.Set(date);
02224 printf("bin: %i date: %i Date: %i Time: %i\n",i,date,time.GetDate(),time.GetTime());
02225 sprintf(filename,"./StarDb/Calibrations/tpc/%s.%08d.%06d.C",TabNam,time.GetDate(),time.GetTime());
02226 printf("Create %s\n",filename);
02227 TString dirname = gSystem->DirName(filename);
02228 if (gSystem->OpenDirectory(dirname.Data())==0) {
02229 if (gSystem->mkdir(dirname.Data())) {
02230 cout << "Directory " << dirname << " creation failed" << endl;
02231 cout << "Putting " << TabNam << ".C in current directory" << endl;
02232 for (int i=0;i<80;i++){filename[i]=0;}
02233 sprintf(filename,"%s.%08d.%06d.C",TabNam.time.GetDate(),time.GetTime());
02234 }
02235 }
02236 ofstream *out = new ofstream(filename);
02237 timegain->SavePrimitive(*out,"");
02238 delete out;
02239
02240 }
02241 }
02242
02243 void FitF(TH2 *hist=0, TF1 *ga=0) {
02244 if (!hist || !ga) return;
02245 Int_t nx = hist->GetNbinsX();
02246 TH1 *proj = 0;
02247 TCanvas *c = new TCanvas("Fit");
02248 TString name(hist->GetName());
02249 name += "MuFG";
02250 Int_t nBins = hist->GetXaxis()->GetNbins();
02251 Double_t xlow = hist->GetXaxis()->GetXmin();
02252 Double_t xup = hist->GetXaxis()->GetXmax();
02253 TH1D *MuF = new TH1D(name.Data(),"Avarage shift versus no. of measurement points",
02254 nBins, xlow, xup);
02255 name = hist->GetName();
02256 name += "SigmaFG";
02257 TH1D *SigmaF = new TH1D(name.Data(),"Sigma of z versus no. of measurement points",
02258 nBins, xlow, xup);
02259 N = 0;
02260 Double_t chisq = -999;
02261 for (int i=1; i<=nBins; i++) {
02262 if (proj) delete proj;
02263 proj = hist->ProjectionY("proj",i,i);
02264 X[N] = hist->GetXaxis()->GetBinCenter(i);
02265 dX[N] = hist->GetXaxis()->GetBinWidth(i);
02266 if (proj->Integral() < 1000) {
02267 if (N>0) goto NEXT;
02268 continue;
02269 }
02270 ga->SetParameter(0,1.e6);
02271 ga->SetParameter(1,0.);
02272
02273 ga->SetParameter(2,0.);
02274
02275 #if 0
02276 ga->SetParameter(1,1.e-9);
02277 ga->SetParLimits(1,1.e-9,1.e-9);
02278
02279 ga->SetParameter(2,1.);
02280 ga->SetParLimits(2,1.,1.);
02281 #endif
02282 ga->SetParameter(3,X[N]);
02283 ga->SetParLimits(3,X[N],X[N]);
02284 proj->Fit(ga->GetName(),"R");
02285 Mu[N] = ga->GetParameter(1);
02286 dMu[N] = ga->GetParError(1);
02287
02288 Sigma[N] = ga->GetParameter(2);
02289 dSigma[N] = ga->GetParError(2);
02290 chisq = ga->GetChisquare();
02291 MuF->SetCellContent(i,0,Mu[N]);
02292 MuF->SetCellError(i,0,dMu[N]);
02293 SigmaF->SetCellContent(i,0,Sigma[N]);
02294 SigmaF->SetCellError(i,0,dSigma[N]);
02295 proj->Draw(); c->Update();
02296 printf("Bin: %i x: %f +/- %f MuF: %f+/-%f Sigma: %f+/-%f chisq: %f \n",
02297 N,X[N],dX[N],Mu[N],dMu[N],Sigma[N],dSigma[N],chisq);
02298 NEXT:
02299 N++;
02300 }
02301 }
02302
02303 void MakeProjY(TH2 *hist, UInt_t ls=0) {
02304 if (! hist) {printf("Missing hist\n"); return;}
02305 if (ls > 1) {printf("Illegal ls = %i\n",ls); return;}
02306 Int_t nx = hist->GetNbinsX();
02307 Int_t ny = hist->GetNbinsY();
02308 TString name(hist->GetName());
02309 static Double_t ShapeMu[2][3] = {
02310 {-2.32656e-01, 7.63952e-02, 6.59020e-03},
02311 {-3.75038e-01, 1.28860e-01,-6.93063e-02}
02312 };
02313 static Double_t ShapeSigma[2][3] = {
02314 { 6.41952e-01,-2.19547e-01, 9.63439e-03},
02315 { 8.29636e-01,-3.71708e-01, 8.65136e-02}
02316 };
02317 name += "pY";
02318 TH1D *h = new TH1D(name.Data(),hist->GetTitle(),ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
02319 for (int i=1;i<=nx;i++){
02320 Double_t x = hist->GetXaxis()->GetBinCenter(i);
02321 Double_t yy = ShapeMu[ls][0] + x*(ShapeMu[ls][1] + x*ShapeMu[ls][2]);
02322 Double_t dy = ShapeSigma[ls][0] + x*(ShapeSigma[ls][1] + x*ShapeSigma[ls][2]);
02323 for (int j=1;j<=ny;j++){
02324 Double_t y = hist->GetYaxis()->GetBinCenter(j);
02325 h->Fill((y-yy)/dy,hist->GetCellContent(i,j));
02326 }
02327 }
02328 }
02329
02330 void badRun(TH1 *hist, Double_t ymin = -0.2, Double_t ymax = 0.) {
02331 if (! hist) return;
02332 Int_t nx = hist->GetNbinsX();
02333 TDatime t;
02334 for (Int_t i=1; i<=nx; i++) {
02335 Double_t y = hist->GetBinContent(i);
02336 if (y < ymin || y > ymax) {
02337 UInt_t d = (UInt_t) hist->GetBinLowEdge(i);
02338 t.Set(d);
02339 cout << "Wrong scale =" << y << " Date " << t.AsSQLString()
02340 << " Date : " << t.AsString() << endl;
02341 }
02342 }
02343 }
02344
02345 Double_t bdEdx(Double_t *xx, Double_t *par) {
02346 Double_t zz = TMath::Log(xx[0]);
02347 Double_t x = par[0];
02348 Double_t y = par[1];
02349 Double_t zprob = gBichsel->GetMostProbableZ(x,y)-5.07402529167365057e-01;
02350 Double_t sigma = gBichsel->GetRmsZ(x,y);
02351 Double_t z = (zz - zprob)/sigma;
02352 return gBichsel->GetProbability(x,y,z)/xx[0]/sigma;
02353 }
02354
02355 Double_t bFunc(Double_t *xx, Double_t *par) {
02356 Double_t z = xx[0];
02357 Double_t x = par[0];
02358 Double_t y = par[1];
02359 return gBichsel->GetProbability(x,y,z);
02360 }
02361
02362 Double_t bFuncPA(Double_t *xx, Double_t *par) {
02363 Double_t x = xx[0];
02364 Double_t y = par[0];
02365 return TMath::Exp(gBichsel->GetMostProbableZ(x,y));
02366 }
02367
02368 Double_t bFuncP(Double_t *xx, Double_t *par) {
02369 Double_t x = xx[0];
02370 Double_t y = par[0];
02371 return gBichsel->GetMostProbableZ(x,y);
02372 }
02373
02374 Double_t bFuncA(Double_t *xx, Double_t *par) {
02375 Double_t x = xx[0];
02376 Double_t y = par[0];
02377 return TMath::Exp(gBichsel->GetAverageZ(x,y));
02378 }
02379
02380 Double_t bFunc70(Double_t *xx, Double_t *par) {
02381 Double_t x = xx[0];
02382 Double_t y = par[0];
02383 return gBichsel->GetI70(x,y);
02384 }
02385
02386 Double_t bFunc60(Double_t *xx, Double_t *par) {
02387 Double_t x = xx[0];
02388 Double_t y = par[0];
02389 return gBichsel->GetI60(x,y);
02390 }
02391
02392 Double_t fncMip(Double_t *xx, Double_t *par) {
02393 Z = xx[0];
02394 Double_t zMostProb = par[0];
02395 Double_t sigma = gBichsel->GetRmsZ(Xlog10bg,Ylog2dx) + par[1];
02396
02397
02398
02399
02400 Double_t z = (Z - zMostProb)/sigma;
02401
02402
02403 Double_t Value = gBichsel->GetProbability(Xlog10bg,Ylog2dx,z)/sigma;
02404 #ifdef PRINT
02405 cout << "Xlog10bg/Ylog2dx/Z =\t" << Xlog10bg << "/" << Ylog2dx << "/" << Z
02406 << "\tzMostProb =\t" << zMostProb
02407 << "\tsigma = \t" << sigma
02408 << "\tValue = \t" << Value << endl;
02409 #endif
02410 return Value;
02411 }
02412
02413 Double_t fncMipX(Double_t *xx, Double_t *par) {
02414 Z = xx[0];
02415 Double_t dx = pow(2.,Ylog2dx);
02416 Double_t dX = dx + par[0];
02417 if (dX < 0.1) dX = 0.1;
02418 Double_t ylog2dx = TMath::Log2(dX);
02419 Double_t zMostProb =
02420 gBichsel->GetMostProbableZ(Xlog10bg,Ylog2dx) -
02421 gBichsel->GetMostProbableZ(Xlog10bg,ylog2dx);
02422 Double_t sigma = gBichsel->GetRmsZ(Xlog10bg,ylog2dx);
02423 Double_t ZZ = Z - TMath::Log(dX/dx);
02424 Double_t z = (ZZ - zMostProb)/sigma;
02425 Double_t Value = gBichsel->GetProbability(Xlog10bg,ylog2dx,z)/sigma;
02426 #ifdef PRINT
02427 cout << "Xlog10bg/Ylog2dx/Z =\t" << Xlog10bg << "/" << Ylog2dx << "/" << Z
02428 << "\tzMostProb =\t" << zMostProb
02429 << "\tsigma = \t" << sigma
02430 << "\tValue = \t" << Value << endl;
02431 #endif
02432 return Value;
02433 }
02434
02435 void bFitMip(const Int_t iX = 8,const Int_t iY=8) {
02436 TDirectory *dir = gDirectory; cout << "Directory: " << dir->GetName() << endl;
02437 if (! gBichsel) {
02438 gSystem->Load("StBichsel");
02439 gBichsel = Bichsel::Instance();
02440 };
02441 dir->cd(); cout << "Directory: " << gDirectory->GetName() << endl;
02442 TString name3D("SecRow3Mip");
02443 canvas = new TCanvas("BinFit","Fit parameters");
02444 TH3 *hist = (TH3 *) gDirectory->Get(name3D);
02445 if (!hist) return;
02446 Int_t nx = hist->GetNbinsX();
02447 Int_t ny = hist->GetNbinsY();
02448 Int_t hyp = 0;
02449 TString tFName("SecRowMipFit");
02450 tFName += NAMES[hyp];
02451 tFName += gSystem->BaseName(gDirectory->GetName());
02452 if (! newf) newf = new TFile(tFName.Data(),"update");
02453 TNtuple *FitP = (TNtuple *) newf->Get("FitP");
02454 if (! FitP) {
02455 FitP = new TNtuple("FitP","Fit results",
02456 "i:j:x:y:mean:rms:peak:mu:sigma:p2:emu:esigma:ep2:sum:chisq:prob:hyp:chisqGP:probGP:peakGP:muGP:sigmaGP:a0:a1:a2:a3:a4:a5:Npar");
02457 FitP->SetMarkerStyle(20);
02458 FitP->SetLineWidth(2);
02459 }
02460 TH2D *p0 = (TH2D *) newf->Get("p0");
02461 if (! p0) p0 = new TH2D("p0","shift of most probable value",
02462 nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
02463 ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
02464 TH2D *p1 = (TH2D *) newf->Get("p1");
02465 if (! p1) p1 = new TH2D("p1","shift of most probable value",
02466 nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
02467 ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
02468 #if 0
02469 TH2D *p2 = (TH2D *) newf->Get("p2");
02470 if (! p2) p2 = new TH2D("p2","shift of most probable value",
02471 nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
02472 ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
02473 #endif
02474 TH2D *chisq = (TH2D *) newf->Get("chisq");
02475 if (! chisq) chisq = new TH2D("chisq","shift of most probable value",
02476 nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
02477 ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
02478 if (! func) func = new TF1("func",fncMip,hist->GetZaxis()->GetXmin(),hist->GetZaxis()->GetXmax(),2);
02479
02480
02481
02482 struct FitMip_t {
02483 Float_t i;
02484 Float_t j;
02485 Float_t x;
02486 Float_t y;
02487 Float_t mean;
02488 Float_t rms;
02489 Float_t peak;
02490 Float_t p0;
02491 Float_t p1;
02492 Float_t p2;
02493 Float_t ep0;
02494 Float_t ep1;
02495 Float_t ep2;
02496 Float_t sum;
02497 Float_t chisq;
02498 Float_t prob;
02499 Float_t hyp;
02500 #if 0
02501 Float_t chisqGP;
02502 Float_t probGP;
02503 Float_t peakGP;
02504 Float_t muGP;
02505 Float_t sigmaGP;
02506 Float_t a0;
02507 Float_t a1;
02508 Float_t a2;
02509 Float_t a3;
02510 Float_t a4;
02511 Float_t a5;
02512 #endif
02513 };
02514 FitMip_t Fit;
02515 Double_t params[9];
02516 Int_t i1 = 1, i2 = nx; if (iX > 0) {i1 = iX; i2 = iX;}
02517 Int_t j1 = 1, j2 = ny; if (iY > 0) {j1 = iY; j2 = iY;}
02518 TF1 *g = 0;
02519 for (Int_t i=i1; i<=i2; i++) {
02520 for (Int_t j=j1; j<=j2; j++) {
02521 newf->cd();
02522 TString projName(Form("%s_%i_%i",NAMES[hyp],i,j));
02523 TH1 *proj = (TH1 *) newf->Get(projName.Data());
02524 if (! proj) proj = hist->ProjectionZ(projName.Data(),i,i,j,j);
02525 if (! proj) continue;
02526 double xx = hist->GetXaxis()->GetBinCenter(i);
02527 double yy = hist->GetYaxis()->GetBinCenter(j);
02528 Ylog2dx = yy;
02529 #if 0
02530 double dy = hist->GetYaxis()->GetBinWidth(j);
02531 double b = 7.5;
02532 if (xx >=14) b = 4;
02533 Double_t xm = yy + 1./b - 0.5*dy/TMath::TanH(0.5*dy*b);
02534 Ylog2dx = xm;
02535 #endif
02536 Double_t sum = proj->Integral();
02537 memset (&Fit.i, 0, sizeof(Fit));
02538 Fit.sum = sum;
02539 Fit.i = i;
02540 Fit.j = j;
02541 Fit.mean = proj->GetMean();
02542 Fit.rms = proj->GetRMS();
02543 Fit.x = xx;
02544 Fit.y = Ylog2dx;
02545 Int_t Row = (int) xx;
02546
02547 Xlog10bg = TMath::Log10(0.448/0.14);
02548 if (sum < 1.e2) {delete proj; continue;}
02549 cout << "Projection:\t"
02550 << proj->GetName()
02551 << "\ti/j\t" << i << "/" << j
02552 << "\tXlog10bg/Ylog2dx\t" << Row << "/" << Ylog2dx
02553 << "\tRow/dx\t" << Row << "/" << pow(2.,Ylog2dx) << "/" << pow(2.,yy)
02554 << "\tIntegral = \t" << sum << endl;
02555 proj->SetTitle(Form("Row = %i dx = %6.2f", Row, pow(2.,Ylog2dx)));
02556 Double_t bw = proj->GetBinWidth(1);
02557 Int_t nb = proj->GetNbinsX();
02558 Int_t l1 = 999, l2 = 0;
02559 for (int l=1; l<=nb; l++) {
02560 Double_t val = proj->GetBinContent(l);
02561 val = val/sum;
02562 Double_t err = TMath::Sqrt(val*(1.-val)/2./sum);
02563 proj->SetBinContent(l,val/bw);
02564 proj->SetBinError(l,err/bw);
02565 if (val <= 0.0) continue;
02566 if (l < l1 ) l1 = l;
02567 l2 = l;
02568 }
02569 if (l1 < l2) proj->GetXaxis()->SetRange(l1,l2);
02570
02571
02572 Int_t lx = proj->GetMaximumBin();
02573 Fit.peak = proj->GetBinCenter(lx);
02574
02575
02576 #if 0
02577 g = FitGP(proj,"r");
02578 g->GetParameters(params);
02579 Fit.Npar = g->GetNpar();
02580 Fit.chisqGP = g->GetChisquare();
02581 Fit.probGP = g->GetProb();
02582 Fit.peakGP = params[0];
02583 Fit.muGP = params[1];
02584 Fit.sigmaGP = params[2];
02585 Fit.a0 = params[3];
02586 Fit.a1 = params[4];
02587 Fit.a2 = params[5];
02588 Fit.a3 = params[6];
02589 Fit.a4 = params[7];
02590 Fit.a5 = params[8];
02591 #endif
02592 params[0] = 0;
02593 params[1] = 0;
02594 func->SetParameters(params);
02595 proj->Fit("func","+em");
02596
02597 Fit.p0 = func->GetParameter(0);
02598 Fit.ep0 = func->GetParError(0);
02599 p0->SetBinContent(i,j,func->GetParameter(0));
02600 p0->SetBinError(i,j,func->GetParError(0));
02601 Fit.p1 = func->GetParameter(1);
02602 Fit.ep1 = func->GetParError(1);
02603 p1->SetBinContent(i,j,func->GetParameter(1));
02604 p1->SetBinError(i,j,func->GetParError(1));
02605 Fit.Npar = g->GetNpar();
02606 Fit.chisq = func->GetChisquare();
02607 chisq->SetBinContent(i,j,func->GetChisquare());
02608 Fit.prob = func->GetProb();
02609 Fit.hyp = hyp;
02610 FitP->Fill(&Fit.i);
02611 proj->Draw();
02612 if (canvas) canvas->Update();
02613 Double_t zz = 0.1303 + gBichsel->GetMostProbableZ(Xlog10bg,Ylog2dx);
02614 func->Eval(zz);
02615 }
02616 }
02617 if (iX == 0 && iY == 0) { newf->cd(); newf->Write(); delete newf;}
02618 }
02619
02620 void MakeTable(St_tpcFeeGainCor *gain) {
02621 if (!gain) return;
02622 TFile *fRootFile = (TFile *) gDirectory->GetFile();
02623 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
02624
02625 Char_t *Nmean = "mu";
02626 TH2D *mean = (TH2D *) fRootFile->Get(Nmean);
02627 if (!mean) {printf("Cannot histogram %s\n",Nmean); return;}
02628 TH2D *entries = (TH2D *) fRootFile->Get("entries");
02629 if (!entries) {printf("Cannot histogram entries\n"); return;}
02630 TH2D *chisq = (TH2D *) fRootFile->Get("chisq");
02631 if (!chisq) {printf("Cannot histogram chisq\n"); return;}
02632 TH2D *sigma = (TH2D *) fRootFile->Get("sigma");
02633 if (!sigma) {printf("Cannot histogram sigma\n"); return;}
02634
02635 Int_t nx = mean->GetNbinsX(); printf ("nx = %i\n",nx);
02636 Int_t ny = mean->GetNbinsY(); printf ("ny = %i\n",ny);
02637 for (int i=1;i<=nx;i++){
02638 tpcFeeGainCor_st row; memset(&row,0,gain->GetRowSize());
02639 for (int j=1;j<=ny;j++){
02640 Double_t res = -1;
02641 if (entries->GetCellContent(i,j) > 0.5e3) {
02642 if (sigma->GetCellContent(i,j) > 0.2) {
02643 Double_t d = mean->GetCellContent(i,j);
02644 if (d > -0.5 && d < 0.5) res = TMath::Exp(-d);
02645 else printf("d: %f i:%i j:%i\n",d,i,j);
02646 }
02647 else printf("entires:%f sigma:%f i:%i j:%i\n"
02648 ,entries->GetCellContent(i,j),sigma->GetCellContent(i,j),i,j);
02649 }
02650 int fee = (j-1)/2;
02651 int eo = (j-1)%2;
02652 row.Gain[fee][eo] = res;
02653 }
02654 gain->AddAt(&row,i-1);
02655 }
02656 }
02657
02658 void FillTable() {
02659 TFile *fRootFile = (TFile *) gDirectory->GetFile();
02660 if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
02661 Char_t *Nmean = "mu";
02662 TH2D *mean = (TH2D *) fRootFile->Get(Nmean);
02663 if (!mean) {printf("Cannot histogram %s\n",Nmean); return;}
02664 TH2D *entries = (TH2D *) fRootFile->Get("entries");
02665 if (!entries) {printf("Cannot histogram entries\n"); return;}
02666 TH2D *chisq = (TH2D *) fRootFile->Get("chisq");
02667 if (!chisq) {printf("Cannot histogram chisq\n"); return;}
02668
02669 Int_t nx = mean->GetNbinsX(); printf ("nx = %i\n",nx);
02670 Int_t ny = mean->GetNbinsY(); printf ("ny = %i\n",ny);
02671 TString NewFile("Correction_");
02672 NewFile += fRootFile->GetName();
02673 NewFile->ReplaceAll(".root",".h");
02674 FILE *fp = fopen(NewFile.Data(),"w");
02675 fprintf(fp,"static Double_t correction[%i][%i] = {// %s %s\n",
02676 ny,nx,mean->GetName(),mean->GetTitle());
02677 for (int j=1;j<=ny;j++){
02678 fprintf(fp,"{");
02679 for (int i=1;i<=nx;i++){
02680 Double_t res = -1;
02681 if (entries->GetCellContent(i,j) > 0.7e4 &&
02682 chisq->GetCellContent(i,j)>0 &&
02683 chisq->GetCellContent(i,j)<2e3
02684 ) {
02685 res = TMath::Exp(-mean->GetCellContent(i,j));
02686 }
02687 if(i != 1) fprintf(fp,",");
02688 fprintf(fp,"%6.3f",res);
02689 if (i==12) fprintf(fp,"\n");
02690 }
02691 fprintf(fp,"}, // fee%i\n",j);
02692 }
02693 fprintf(fp,"};\n");
02694 fclose(fp);
02695 }
02696
02697 void Make2Dproj(TH2D *h, Int_t N = 100, Double_t xlow = -1., Double_t xup = 1.) {
02698 if (!h) return;
02699 Int_t nx = h->GetNbinsX();
02700 Int_t ny = h->GetNbinsY();
02701 TString name(h->GetName());
02702 name += "_P";
02703 TH1D *hist = new TH1D(name.Data(),h->GetTitle(),N,xlow,xup);
02704 for (int i = 1; i<=nx;i++) {
02705 for (int j = 1; j <= ny; j++){
02706 if (TMath::Abs(h->GetCellContent(i,j)) < 0.01) continue;
02707 hist->Fill(h->GetCellContent(i,j));
02708 }
02709 }
02710 }
02711 #endif
02712
02713
02714 TList *ListOfKeys() {
02715 TList *list = 0;
02716 TSeqCollection *files = gROOT->GetListOfFiles();
02717 TIter next(files);
02718 TFile *f = 0;
02719 while ((f = (TFile *) next())) {
02720 f->cd();
02721 list = f->GetListOfKeys();
02722 break;
02723 }
02724 return list;
02725 }
02726
02727 TH1 *FindHistograms(const Char_t *Name, const Char_t *fit = "") {
02728 TH1 *hist = 0;
02729 #ifdef DEBUG
02730 cout << "FindHistograms(\"" << Name << "\",\"" << fit << "\")" << endl;
02731 #endif
02732 TSeqCollection *files = gROOT->GetListOfFiles();
02733 TIter next(files);
02734 TFile *f = 0;
02735 TString Fit(fit);
02736 if (Fit != "") {
02737 while ((f = (TFile *) next())) {
02738 #ifdef DEBUG
02739 cout << "look in " << f->GetName() << endl;
02740 #endif
02741 TString fName(gSystem->BaseName(f->GetName()));
02742 if (fName.BeginsWith(Name)) {
02743 hist = (TH1 *) f->Get(fit);
02744 #ifdef DEBUG
02745 if (hist) cout << "Found histogram " << hist->GetName() << " in file " << f->GetName() << endl;
02746 #endif
02747 return hist;
02748 }
02749 }
02750 } else {
02751 while ((f = (TFile *) next())) {
02752 f->cd();
02753 #ifdef DEBUG
02754 cout << "look in " << f->GetName() << endl;
02755 #endif
02756 hist = (TH1 *) f->Get(Name);
02757 if (hist) {
02758 #ifdef DEBUG
02759 cout << "Found histogram " << hist->GetName() << " in file " << f->GetName() << endl;
02760 #endif
02761 return hist;
02762 }
02763 }
02764 }
02765 return hist;
02766 }
02767
02768 void DrawSummary0(TFile *f, const Char_t *opt) {
02769 if (! f) return;
02770 gStyle->SetTimeOffset(788936400);
02771 f->cd();
02772 TList *keys = f->GetListOfKeys();
02773 if (! keys) return;
02774 keys->Sort();
02775 TIter next(keys);
02776 TKey *key = 0;
02777 TObjArray hists;
02778 TString Opt(opt);
02779 TString Title;
02780 TF1 *powfit = new TF1("powfit","[0]*pow(x,[1])",40,80);
02781 powfit->SetParameters(0.5,-0.5);
02782 while ((key = (TKey *) next())) {
02783 TH1 *hist = FindHistograms(key->GetName());
02784 if (! hist) continue;
02785 if (hist->GetEntries() < 1) continue;
02786 if (Opt != "") {
02787 TString name(hist->GetName());
02788 if (! name.Contains(Opt.Data())) continue;
02789 }
02790 hists.AddLast(hist);
02791 TString hName(hist->GetName());
02792 if ((( hist->IsA()->InheritsFrom( "TH3C" ) ||
02793 hist->IsA()->InheritsFrom( "TH3S" ) ||
02794 hist->IsA()->InheritsFrom( "TH3I" ) ||
02795 hist->IsA()->InheritsFrom( "TH3F" ) ||
02796 hist->IsA()->InheritsFrom( "TH3D" ) ) ||
02797 ( hist->IsA()->InheritsFrom( "TH2C" ) ||
02798 hist->IsA()->InheritsFrom( "TH2S" ) ||
02799 hist->IsA()->InheritsFrom( "TH2I" ) ||
02800 hist->IsA()->InheritsFrom( "TH2F" ) ||
02801 hist->IsA()->InheritsFrom( "TH2D" ) ))&&
02802 ( ! (hName.BeginsWith("Points") || hName.BeginsWith("TPoints") || hName.BeginsWith("MPoints")) ) ) {
02803 TH1 *mu = FindHistograms(key->GetName(),"mu");
02804 if (! mu) continue;
02805 mu->SetName(Form("%s_mu",hist->GetName()));
02806 hists.AddLast(mu);
02807 TH1 *sigma = FindHistograms(key->GetName(),"sigma");
02808 if (! sigma) continue;
02809 sigma->SetName(Form("%s_sigma",hist->GetName()));
02810 hists.AddLast(sigma);
02811 }
02812 }
02813 Int_t N = hists.GetEntriesFast();
02814 Int_t nx = (Int_t) (TMath::Sqrt(N));
02815 if (nx*nx < N) nx++;
02816 Int_t ny = N/nx;
02817 if (nx*ny < N) ny++;
02818 cout << "N " << N << " nx " << nx << " ny " << ny << endl;
02819 TString Tag(gSystem->BaseName(f->GetName()));
02820 Tag.ReplaceAll(".root","");
02821 Tag += opt;
02822 TCanvas *c1 = new TCanvas(Tag,Tag,200,10,700,780);
02823 Double_t dx = 0.98/nx;
02824 Double_t dy = 0.98/ny;
02825 TObjArray pads(N);
02826 for (Int_t iy = 0; iy < ny; iy++) {
02827 for (Int_t ix = 0; ix < nx; ix++) {
02828 Int_t i = ix + nx*iy;
02829 if (i >= N) continue;
02830 TH1 *hist = (TH1 *) hists[i];
02831 if (! hist) continue;
02832 Double_t x1 = 0.01 + dx* ix;
02833 Double_t x2 = 0.01 + dx*(ix+1);
02834 Double_t y1 = 0.99 - dy*(iy+1);
02835 Double_t y2 = 0.99 - dy* iy;
02836 TPad *pad = new TPad(Form("pad_%s",hist->GetName()),hist->GetTitle(),x1,y1,x2,y2);
02837
02838 pad->Draw();
02839 pads.AddLast(pad);
02840 }
02841 }
02842 for (Int_t iy = 0; iy < ny; iy++) {
02843 for (Int_t ix = 0; ix < nx; ix++) {
02844 Int_t i = ix + nx*iy;
02845 if (i >= N) continue;
02846 TH1 *hist = (TH1 *) hists[i];
02847 if (! hist) continue;
02848 TPad *pad = (TPad *) pads[i];
02849 if (! pad) continue;
02850 pad->cd();
02851 TProfile *prof = 0;
02852 if ( hist->IsA()->InheritsFrom( "TProfile" ) ) prof = (TProfile *) hist;
02853 Int_t NbinsX = hist->GetNbinsX();
02854 Int_t xmin = NbinsX;
02855 Int_t xmax = 0;
02856 Double_t ymin = 1e9;
02857 Double_t ymax = -1e9;
02858
02859
02860 if ( hist->IsA()->InheritsFrom( "TH3C" ) ||
02861 hist->IsA()->InheritsFrom( "TH3S" ) ||
02862 hist->IsA()->InheritsFrom( "TH3I" ) ||
02863 hist->IsA()->InheritsFrom( "TH3F" ) ||
02864 hist->IsA()->InheritsFrom( "TH3D" ) ) {
02865 ((TH3 *)hist)->Project3DProfile("xy")->Draw("colz");
02866 goto TIMEAxis;
02867 }
02868
02869 if ( hist->IsA()->InheritsFrom( "TH2C" ) ||
02870 hist->IsA()->InheritsFrom( "TH2S" ) ||
02871 hist->IsA()->InheritsFrom( "TH2I" ) ||
02872 hist->IsA()->InheritsFrom( "TH2F" ) ||
02873 hist->IsA()->InheritsFrom( "TH2D" ) ) {
02874 if (hist->GetMaximum() > 0 && hist->GetMinimum() >= 0) pad->SetLogz(1);
02875 TString hName(hist->GetName());
02876 if (hName.BeginsWith("Points") || hName.BeginsWith("TPoints") || hName.BeginsWith("MPoints")) {
02877 Title = hist->GetTitle();
02878 Title.ReplaceAll("/sigma","");
02879 hist->SetTitle(Title);
02880 TAxis *y = hist->GetYaxis();
02881 Int_t iy1 = y->FindBin(-0.15);
02882 Int_t iy2 = y->FindBin( 0.50);
02883 y->SetRange(iy1,iy2);
02884 ((TH2 *)hist)->Draw("colz");
02885 TH1 *mu = FindHistograms(hist->GetName(),"mu");
02886 if (! mu) goto TIMEAxis;
02887 mu->SetMarkerColor(1);
02888 mu->SetMarkerStyle(20);
02889 mu->Draw("same");
02890 mu->Fit("pol0","e0r","",10,120);
02891 TLegend *leg = new TLegend(0.25,0.6,0.9,0.9,"");
02892 TF1 *f = mu->GetFunction("pol0");
02893 if (f) {
02894 f->Draw("same");
02895 Title = Form("#mu = %5.2f +/- %5.2f %\%",100*f->GetParameter(0),100*f->GetParError(0));
02896 cout << Title << endl;
02897 leg->AddEntry(mu,Title.Data());
02898 }
02899 TH1 *sigma = FindHistograms(hist->GetName(),"sigma");
02900 if (! sigma) goto TIMEAxis;
02901 sigma->SetMarkerColor(1);
02902 sigma->SetMarkerStyle(21);
02903 sigma->Draw("same");
02904 sigma->Fit(powfit,"r0");
02905 f = sigma->GetFunction("powfit");
02906 if (f) {
02907 f->Draw("same");
02908 Title = Form("#sigma(@76cm) = %5.2f%\%",100*f->Eval(76));
02909 cout << Title << endl;
02910 leg->AddEntry(sigma,Title.Data());
02911 }
02912 leg->Draw();
02913 } else hist->Draw("colz");
02914 goto TIMEAxis;
02915 }
02916
02917 NbinsX = hist->GetNbinsX();
02918 xmin = NbinsX;
02919 xmax = 0;
02920 ymin = 1e9;
02921 ymax = -1e9;
02922 for (Int_t i = 1; i <= NbinsX; i++) {
02923 Double_t y = hist->GetBinContent(i);
02924 if ((prof && prof->GetBinEntries(i)) || y > 0) {
02925 Int_t x = i;
02926 xmin = TMath::Min(xmin,x);
02927 xmax = TMath::Max(xmax,x);
02928 ymin = TMath::Min(ymin,y);
02929 ymax = TMath::Max(ymax,y);
02930 }
02931 }
02932 if (ymin < ymax) {
02933 hist->SetMaximum(1.1*ymax);
02934 hist->SetMinimum(0.9*ymin);
02935
02936 }
02937 if (xmin < xmax) hist->GetXaxis()->SetRange(xmin,xmax);
02938 hist->Draw();
02939 TIMEAxis:
02940 TAxis *xax = hist->GetXaxis();
02941 if (xax->GetBinLowEdge(1) > 1e6) {
02942 xax->SetTimeDisplay(1);
02943 gPad->Modified();
02944 }
02945 }
02946 }
02947 }
02948
02949 void DrawSummary(const Char_t *opt="") {
02950 TSeqCollection *files = gROOT->GetListOfFiles();
02951 TIter next(files);
02952 TFile *f = 0;
02953 TString FName("");
02954 while ((f = (TFile *) next())) {
02955 FName = gSystem->BaseName(f->GetName());
02956 if (FName.BeginsWith("Hist")) break;
02957 }
02958 if (! f) {
02959 cout << "Hist* root file has not been found" << endl;
02960 return;
02961 }
02962 DrawSummary0(f,opt);
02963 }
