[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / Train / TwoTrackQA / FemtoQAPlots.C

void preparepad()
{
  gPad->SetFillColor(0);
  gPad->SetFillStyle(4000);

  gPad->SetTopMargin(0.02);
  //  gPad->SetRightMargin(0.02);
}

void preparehist(TH1D *hist)
{
  hist->SetMarkerSize(0.8);
  hist->SetMarkerStyle(8);
  hist->SetMarkerColor(2);
}

TH2D *plot2ddtpc(TH2D *numdtpc, TH2D *dendtpc, const char *suff)
{
  // Create a 2D CF vs minimum separation distance in the TPC


  int nslices = numdtpc->GetNbinsY();
  int nrange = numdtpc->GetNbinsX();
  char hname[200];
  
  sprintf(hname,"numdtpc%s",suff);
  TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]", 
			numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(), 
			numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
  sprintf(hname,"dendtpc%s",suff);
  TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]", 
			numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(), 
			numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
  sprintf(hname,"ratdtpc%s",suff);
  TH2D *rats = new TH2D(hname,  ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]", 
			numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(), 
			numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
  
  sprintf(hname,"chi2toone%s",suff);
  TH1D *chi2toone = new TH1D(hname,";nominal TPC entrance sep - lower cut-off [cm];\\chi^2",
			     numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());

  char bufname[100];
  for (int iter=0; iter<nslices; iter++) {
    sprintf(bufname, "numbuf%i", iter);
    TH1D *numbuf = numdtpc->ProjectionX(bufname, (iter+1),  nslices, "e");
    sprintf(bufname, "denbuf%i", iter);
    TH1D *denbuf = dendtpc->ProjectionX(bufname, (iter+1),  nslices, "e");
    
    TH1D *ratbuf = new TH1D(*numbuf);
    ratbuf->Divide(denbuf);
    Double_t scale = ratbuf->Integral(nrange-9,nrange)/10.0;

    if (scale > 0.0)
      ratbuf->Scale(1.0/scale);
    
    double chi2 = 0.0;
    for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
      rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
      if (ratbuf->GetBinError(ibin) > 0.0) {
	chi2 += (TMath::Power(ratbuf->GetBinContent(ibin)-1.0,2)/
		 TMath::Power(ratbuf->GetBinError(ibin),2));
      }
    }
    chi2toone->SetBinContent(iter+1, chi2);
  }

  rats->GetXaxis()->SetTitleOffset(0.9);
  rats->GetYaxis()->SetTitleOffset(0.95);

  return rats;
}

TH2D * plot2dshare(TH2D *numshare, TH2D *denshare, const char *suff)
{
  // Create a 2D CF vs minimum separation distance in the TPC
  int nrange = numshare->GetNbinsX();
  char hname[200];
    
  sprintf(hname,"numshare%s",suff);
  TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off", 
			numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
			numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());

  sprintf(hname,"denshare%s",suff);
  TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off", 
			numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
			numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());

  sprintf(hname,"ratshare%s",suff);
  TH2D *rats = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off", 
			numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
			numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());
  
  char bufname[100];
  for (int iter=0; iter<numshare->GetNbinsY(); iter++) {
    sprintf(bufname, "numbuf%i", iter);
    TH1D *numbuf = numshare->ProjectionX(bufname, 1, iter+1, "e");
    sprintf(bufname, "denbuf%i", iter);
    TH1D *denbuf = denshare->ProjectionX(bufname, 1, iter+1, "e");
    
    ratbuf = new TH1D(*numbuf);
    ratbuf->Divide(denbuf);
    Double_t scale = ratbuf->Integral(nrange-24,nrange)/25.0;
    ratbuf->Scale(1.0/scale);

    for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
      rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
    }
  }
  
  rats->GetXaxis()->SetTitleOffset(0.95);
  rats->GetYaxis()->SetTitleOffset(0.9);

  return rats;
}

TH2D *plot2dquality(TH2D *numqual, TH2D *denqual, const char *suff)
{
  // Create a 2D CF vs pair quality 
  int nrange = numqual->GetNbinsX();
  char hname[200];
    
  sprintf(hname,"numqual%s",suff);
  TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off", 
			numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
			numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
  sprintf(hname,"denqual%s",suff);
  TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off", 
			numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
			numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
  sprintf(hname,"ratqual%s",suff);
  TH2D *rats = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off", 
			numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
			numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
  
  char bufname[100];
  for (int iter=0; iter<numqual->GetNbinsY(); iter++) {
    sprintf(bufname, "numbuf%i", iter);
    TH1D *numbuf = numqual->ProjectionX(bufname, 1, iter+1, "e");
    sprintf(bufname, "denbuf%i", iter);
    TH1D *denbuf = denqual->ProjectionX(bufname, 1, iter+1, "e");
    
    ratbuf = new TH1D(*numbuf);
    ratbuf->Divide(denbuf);
    Double_t scale = ratbuf->Integral(nrange-24,nrange)/25.0;
    if (scale > 0.0)
      ratbuf->Scale(1.0/scale);

    for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
      rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
    }
  }
  
  rats->GetXaxis()->SetTitleOffset(0.95);
  rats->GetYaxis()->SetTitleOffset(0.9);

  return rats;
}

TH1D *getcf(TH1D *numq, TH1D *denq, const char *suff)
{
  char hname[200];

  TH1D *rats = new TH1D(*numq);
  rats->Reset("ICE");
  rats->Divide(numq, denq, 1.0*denq->Integral()/numq->Integral(), 1.0);

  sprintf(hname, "ratqinv%s", suff);
  rats->SetName(hname);
  rats->SetTitle(";q_{inv} [GeV/c];C(q_{inv})");
  preparehist(rats);

  return rats;
}

void FemtoQAPlots(const char *infname)
{
  gStyle->SetStatX(0.8);
  gStyle->SetStatW(0.25);
  gStyle->SetOptStat(11);

  TFile *ifile = new TFile(infname);
  
  gStyle->SetPalette(1,0);

  // Make plot for TPC entrance separation - positive

  TCanvas *candtpcp = new TCanvas("candtpcp","candtpcp",1200,300);
  preparepad();
  candtpcp->Divide(3,1,0.0001,0.0001);
  
  candtpcp->cd(1);  preparepad();
  TH2D *ratdtpcpipstd = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpipstd"), "pipstd");
  ratdtpcpipstd->Draw("COLZ");

  candtpcp->cd(2);  preparepad();
  TH2D *ratdtpcpipnct = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpipnct"), "pipnct");
  ratdtpcpipnct->Draw("COLZ");

  candtpcp->cd(3);  preparepad();
  TH2D *ratdtpcpiptpc = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpiptpc"), "piptpc");
  ratdtpcpiptpc->Draw("COLZ");

  // Make plot for TPC entrance separation - negative

  TCanvas *candtpcm = new TCanvas("candtpcm","candtpcm",1200,300);
  preparepad();
  candtpcm->Divide(3,1,0.0001,0.0001);
  
  candtpcm->cd(1);  preparepad();
  TH2D *ratdtpcpimstd = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimstd"), "pimstd");
  ratdtpcpimstd->Draw("COLZ");

  candtpcm->cd(2);  preparepad();
  TH2D *ratdtpcpimnct = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimnct"), "pimnct");
  ratdtpcpimnct->Draw("COLZ");

  candtpcm->cd(3);  preparepad();
  TH2D *ratdtpcpimtpc = plot2ddtpc((TH2D *) gDirectory->Get("NumTPCSepsqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimtpc"), "pimtpc");
  ratdtpcpimtpc->Draw("COLZ");

  // Make plot for pair fraction of shared hits - positive

  TCanvas *cansharep = new TCanvas("cansharep","cansharep",1200,300);
  preparepad();
  cansharep->Divide(3,1,0.0001,0.0001);
  
  cansharep->cd(1);  preparepad();
  TH2D *ratsharepipstd = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpipstd"), "pipstd");
  ratsharepipstd->Draw("COLZ");

  cansharep->cd(2);  preparepad();
  TH2D *ratsharepipnct = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpipnct"), "pipnct");
  ratsharepipnct->Draw("COLZ");

  cansharep->cd(3);  preparepad();
  TH2D *ratsharepiptpc = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpiptpc"), "piptpc");
  ratsharepiptpc->Draw("COLZ");

  // Make plot for pair fraction of shared hits - negative

  TCanvas *cansharem = new TCanvas("cansharem","cansharem",1200,300);
  preparepad();
  cansharem->Divide(3,1,0.0001,0.0001);
  
  cansharem->cd(1);  preparepad();
  TH2D *ratsharepimstd = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimstd"), "pimstd");
  ratsharepimstd->Draw("COLZ");

  cansharem->cd(2);  preparepad();
  TH2D *ratsharepimnct = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimnct"), "pimnct");
  ratsharepimnct->Draw("COLZ");

  cansharem->cd(3);  preparepad();
  TH2D *ratsharepimtpc = plot2dshare((TH2D *) gDirectory->Get("NumSharesqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimtpc"), "pimtpc");
  ratsharepimtpc->Draw("COLZ");

  // Make signal pair fraction of shared hits - positive

  TCanvas *cansharesignalp = new TCanvas("cansharesignalp","cansharesignalp",1200,300);
  preparepad();
  cansharesignalp->Divide(3,1,0.0001,0.0001);
  
  cansharesignalp->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpipstd"))->Draw("COLZ");

  cansharesignalp->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpipnct"))->Draw("COLZ");

  cansharesignalp->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpiptpc"))->Draw("COLZ");

  // Make signal pair fraction of shared hits - negative

  TCanvas *cansharesignalm = new TCanvas("cansharesignalm","cansharesignalm",1200,300);
  preparepad();
  cansharesignalm->Divide(3,1,0.0001,0.0001);
  
  cansharesignalm->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpimstd"))->Draw("COLZ");

  cansharesignalm->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpimnct"))->Draw("COLZ");

  cansharesignalm->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("NumSharesqqinvcfpimtpc"))->Draw("COLZ");

  // Make background pair fraction of shared hits - positive

  TCanvas *cansharebackp = new TCanvas("cansharebackp","cansharebackp",1200,300);
  preparepad();
  cansharebackp->Divide(3,1,0.0001,0.0001);
  
  cansharebackp->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpipstd"))->Draw("COLZ");

  cansharebackp->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpipnct"))->Draw("COLZ");

  cansharebackp->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpiptpc"))->Draw("COLZ");

  // Make background pair fraction of shared hits - negative

  TCanvas *cansharebackm = new TCanvas("cansharebackm","cansharebackm",1200,300);
  preparepad();
  cansharebackm->Divide(3,1,0.0001,0.0001);
  
  cansharebackm->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpimstd"))->Draw("COLZ");

  cansharebackm->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpimnct"))->Draw("COLZ");

  cansharebackm->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("DenSharesqqinvcfpimtpc"))->Draw("COLZ");

  // Make plot for pair quality - positive

  TCanvas *canqualp = new TCanvas("canqualp","canqualp",1200,300);
  preparepad();
  canqualp->Divide(3,1,0.0001,0.0001);
  
  canqualp->cd(1);  preparepad();
  TH2D *ratqualpipstd = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpipstd"), "pipstd");
  ratqualpipstd->Draw("COLZ");

  canqualp->cd(2);  preparepad();
  TH2D *ratqualpipnct = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpipnct"), "pipnct");
  ratqualpipnct->Draw("COLZ");

  canqualp->cd(3);  preparepad();
  TH2D *ratqualpiptpc = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpiptpc"), "piptpc");
  ratqualpiptpc->Draw("COLZ");

  // Make plot for pair quality - negative

  TCanvas *canqualm = new TCanvas("canqualm","canqualm",1200,300);
  preparepad();
  canqualm->Divide(3,1,0.0001,0.0001);
  
  canqualm->cd(1);  preparepad();
  TH2D *ratqualpimstd = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimstd"), "pimstd");
  ratqualpimstd->Draw("COLZ");

  canqualm->cd(2);  preparepad();
  TH2D *ratqualpimnct = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimnct"), "pimnct");
  ratqualpimnct->Draw("COLZ");

  canqualm->cd(3);  preparepad();
  TH2D *ratqualpimtpc = plot2dquality((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimtpc"), "pimtpc");
  ratqualpimtpc->Draw("COLZ");

  // Make signal pair quality - positive

  TCanvas *canqualitysignalp = new TCanvas("canqualitysignalp","canqualitysignalp",1200,300);
  preparepad();
  canqualitysignalp->Divide(3,1,0.0001,0.0001);
  
  canqualitysignalp->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipstd"))->Draw("COLZ");

  canqualitysignalp->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipnct"))->Draw("COLZ");

  canqualitysignalp->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpiptpc"))->Draw("COLZ");

  // Make signal pair quality - negative

  TCanvas *canqualitysignalm = new TCanvas("canqualitysignalm","canqualitysignalm",1200,300);
  preparepad();
  canqualitysignalm->Divide(3,1,0.0001,0.0001);
  
  canqualitysignalm->cd(1);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimstd"))->Draw("COLZ");

  canqualitysignalm->cd(2);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimnct"))->Draw("COLZ");

  canqualitysignalm->cd(3);  preparepad();
  ((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimtpc"))->Draw("COLZ");

  // Make plot for pair quality - positive

  TCanvas *canqinvp = new TCanvas("canqinvp","canqinvp",1200,300);
  preparepad();
  canqinvp->Divide(3,1,0.0001,0.0001);
  
  canqinvp->cd(1); preparepad();
  TH1D *ratqinvpipstd = getcf((TH1D *) gDirectory->Get("Numqinvcfpipstd"), (TH1D *) gDirectory->Get("Denqinvcfpipstd"), "pipstd");
  ratqinvpipstd->Draw();

  canqinvp->cd(2); preparepad();
  TH1D *ratqinvpipnct = getcf((TH1D *) gDirectory->Get("Numqinvcfpipnct"), (TH1D *) gDirectory->Get("Denqinvcfpipnct"), "pipnct");
  ratqinvpipnct->Draw("");

  canqinvp->cd(3); preparepad();
  TH1D *ratqinvpiptpc = getcf((TH1D *) gDirectory->Get("Numqinvcfpiptpc"), (TH1D *) gDirectory->Get("Denqinvcfpiptpc"), "piptpc");
  ratqinvpiptpc->Draw("");

  // Make plot for pair quality - negative

  TCanvas *canqinvm = new TCanvas("canqinvm","canqinvm",1200,300);
  preparepad();
  canqinvm->Divide(3,1,0.0001,0.0001);
  
  canqinvm->cd(1); preparepad();
  TH1D *ratqinvpimstd = getcf((TH1D *) gDirectory->Get("Numqinvcfpimstd"), (TH1D *) gDirectory->Get("Denqinvcfpimstd"), "pimstd");
  ratqinvpimstd->Draw("");

  canqinvm->cd(2); preparepad();
  TH1D *ratqinvpimnct = getcf((TH1D *) gDirectory->Get("Numqinvcfpimnct"), (TH1D *) gDirectory->Get("Denqinvcfpimnct"), "pimnct");
  ratqinvpimnct->Draw("");

  canqinvm->cd(3); preparepad();
  TH1D *ratqinvpimtpc = getcf((TH1D *) gDirectory->Get("Numqinvcfpimtpc"), (TH1D *) gDirectory->Get("Denqinvcfpimtpc"), "pimtpc");
  ratqinvpimtpc->Draw("");

  
}
Commit	Line	Data
76ce4b5b	1	void preparepad()
	2	{
	3	gPad->SetFillColor(0);
	4	gPad->SetFillStyle(4000);
	5
	6	gPad->SetTopMargin(0.02);
	7	// gPad->SetRightMargin(0.02);
	8	}
	9
	10	void preparehist(TH1D *hist)
	11	{
	12	hist->SetMarkerSize(0.8);
	13	hist->SetMarkerStyle(8);
	14	hist->SetMarkerColor(2);
	15	}
	16
	17	TH2D plot2ddtpc(TH2D numdtpc, TH2D dendtpc, const char suff)
	18	{
	19	// Create a 2D CF vs minimum separation distance in the TPC
	20
	21
	22
	23	int nslices = numdtpc->GetNbinsY();
	24	int nrange = numdtpc->GetNbinsX();
	25	char hname[200];
	26
	27	sprintf(hname,"numdtpc%s",suff);
	28	TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]",
	29	numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(),
	30	numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
	31	sprintf(hname,"dendtpc%s",suff);
	32	TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]",
	33	numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(),
	34	numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
	35	sprintf(hname,"ratdtpc%s",suff);
	36	TH2D *rats = new TH2D(hname, ";q_{inv} [GeV/c];nominal TPC entrance sep - lower cut-off [cm]",
	37	numdtpc->GetNbinsX(), numdtpc->GetXaxis()->GetXmin(), numdtpc->GetXaxis()->GetXmax(),
	38	numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
	39
	40	sprintf(hname,"chi2toone%s",suff);
	41	TH1D *chi2toone = new TH1D(hname,";nominal TPC entrance sep - lower cut-off [cm];\\chi^2",
	42	numdtpc->GetNbinsY(), numdtpc->GetYaxis()->GetXmin(), numdtpc->GetYaxis()->GetXmax());
	43
	44	char bufname[100];
	45	for (int iter=0; iter<nslices; iter++) {
	46	sprintf(bufname, "numbuf%i", iter);
	47	TH1D *numbuf = numdtpc->ProjectionX(bufname, (iter+1), nslices, "e");
	48	sprintf(bufname, "denbuf%i", iter);
	49	TH1D *denbuf = dendtpc->ProjectionX(bufname, (iter+1), nslices, "e");
	50
	51	TH1D ratbuf = new TH1D(numbuf);
	52	ratbuf->Divide(denbuf);
	53	Double_t scale = ratbuf->Integral(nrange-9,nrange)/10.0;
	54
	55	if (scale > 0.0)
	56	ratbuf->Scale(1.0/scale);
	57
	58	double chi2 = 0.0;
	59	for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
	60	rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
	61	if (ratbuf->GetBinError(ibin) > 0.0) {
	62	chi2 += (TMath::Power(ratbuf->GetBinContent(ibin)-1.0,2)/
	63	TMath::Power(ratbuf->GetBinError(ibin),2));
	64	}
65	}
66	chi2toone->SetBinContent(iter+1, chi2);
67	}
68
69	rats->GetXaxis()->SetTitleOffset(0.9);
70	rats->GetYaxis()->SetTitleOffset(0.95);
71
72	return rats;
73	}
74
75	TH2D * plot2dshare(TH2D numshare, TH2D denshare, const char *suff)
76	{
77	// Create a 2D CF vs minimum separation distance in the TPC
78	int nrange = numshare->GetNbinsX();
79	char hname[200];
80
81	sprintf(hname,"numshare%s",suff);
82	TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off",
83	numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
84	numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());
85
86	sprintf(hname,"denshare%s",suff);
87	TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off",
88	numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
89	numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());
90
91	sprintf(hname,"ratshare%s",suff);
92	TH2D *rats = new TH2D(hname, ";q_{inv} [GeV/c];fraction of shared hits - upper cut-off",
93	numshare->GetNbinsX(), numshare->GetXaxis()->GetXmin(), numshare->GetXaxis()->GetXmax(),
94	numshare->GetNbinsY(), numshare->GetYaxis()->GetXmin(), numshare->GetYaxis()->GetXmax());
95
96	char bufname[100];
97	for (int iter=0; iter<numshare->GetNbinsY(); iter++) {
98	sprintf(bufname, "numbuf%i", iter);
99	TH1D *numbuf = numshare->ProjectionX(bufname, 1, iter+1, "e");
100	sprintf(bufname, "denbuf%i", iter);
101	TH1D *denbuf = denshare->ProjectionX(bufname, 1, iter+1, "e");
102
103	ratbuf = new TH1D(*numbuf);
104	ratbuf->Divide(denbuf);
105	Double_t scale = ratbuf->Integral(nrange-24,nrange)/25.0;
106	ratbuf->Scale(1.0/scale);
107
108	for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
109	rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
110	}
111	}
112
113	rats->GetXaxis()->SetTitleOffset(0.95);
114	rats->GetYaxis()->SetTitleOffset(0.9);
115
116	return rats;
117	}
118
119	TH2D plot2dquality(TH2D numqual, TH2D denqual, const char suff)
120	{
121	// Create a 2D CF vs pair quality
122	int nrange = numqual->GetNbinsX();
123	char hname[200];
124
125	sprintf(hname,"numqual%s",suff);
126	TH2D *nums = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off",
127	numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
128	numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
129	sprintf(hname,"denqual%s",suff);
130	TH2D *dens = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off",
131	numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
132	numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
133	sprintf(hname,"ratqual%s",suff);
134	TH2D *rats = new TH2D(hname, ";q_{inv} [GeV/c];pair quality - upper cut-off",
135	numqual->GetNbinsX(), numqual->GetXaxis()->GetXmin(), numqual->GetXaxis()->GetXmax(),
136	numqual->GetNbinsY(), numqual->GetYaxis()->GetXmin(), numqual->GetYaxis()->GetXmax());
137
138	char bufname[100];
139	for (int iter=0; iter<numqual->GetNbinsY(); iter++) {
140	sprintf(bufname, "numbuf%i", iter);
141	TH1D *numbuf = numqual->ProjectionX(bufname, 1, iter+1, "e");
142	sprintf(bufname, "denbuf%i", iter);
143	TH1D *denbuf = denqual->ProjectionX(bufname, 1, iter+1, "e");
144
145	ratbuf = new TH1D(*numbuf);
146	ratbuf->Divide(denbuf);
147	Double_t scale = ratbuf->Integral(nrange-24,nrange)/25.0;
148	if (scale > 0.0)
149	ratbuf->Scale(1.0/scale);
150
151	for (int ibin = 1; ibin<=ratbuf->GetNbinsX(); ibin++) {
152	rats->SetBinContent(ibin, iter+1, ratbuf->GetBinContent(ibin));
153	}
154	}
155
156	rats->GetXaxis()->SetTitleOffset(0.95);
157	rats->GetYaxis()->SetTitleOffset(0.9);
158
159	return rats;
160	}
161
162	TH1D getcf(TH1D numq, TH1D denq, const char suff)
163	{
164	char hname[200];
165
166	TH1D rats = new TH1D(numq);
167	rats->Reset("ICE");
168	rats->Divide(numq, denq, 1.0*denq->Integral()/numq->Integral(), 1.0);
169
170	sprintf(hname, "ratqinv%s", suff);
171	rats->SetName(hname);
172	rats->SetTitle(";q_{inv} [GeV/c];C(q_{inv})");
173	preparehist(rats);
174
175	return rats;
176	}
177
178	void FemtoQAPlots(const char *infname)
179	{
180	gStyle->SetStatX(0.8);
181	gStyle->SetStatW(0.25);
182	gStyle->SetOptStat(11);
183
184	TFile *ifile = new TFile(infname);
185
186	gStyle->SetPalette(1,0);
187
188	// Make plot for TPC entrance separation - positive
189
190	TCanvas *candtpcp = new TCanvas("candtpcp","candtpcp",1200,300);
191	preparepad();
192	candtpcp->Divide(3,1,0.0001,0.0001);
193
194	candtpcp->cd(1); preparepad();
195	TH2D ratdtpcpipstd = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpipstd"), "pipstd");
196	ratdtpcpipstd->Draw("COLZ");
197
198	candtpcp->cd(2); preparepad();
199	TH2D ratdtpcpipnct = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpipnct"), "pipnct");
200	ratdtpcpipnct->Draw("COLZ");
201
202	candtpcp->cd(3); preparepad();
203	TH2D ratdtpcpiptpc = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpiptpc"), "piptpc");
204	ratdtpcpiptpc->Draw("COLZ");
205
206	// Make plot for TPC entrance separation - negative
207
208	TCanvas *candtpcm = new TCanvas("candtpcm","candtpcm",1200,300);
209	preparepad();
210	candtpcm->Divide(3,1,0.0001,0.0001);
211
212	candtpcm->cd(1); preparepad();
213	TH2D ratdtpcpimstd = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimstd"), "pimstd");
214	ratdtpcpimstd->Draw("COLZ");
215
216	candtpcm->cd(2); preparepad();
217	TH2D ratdtpcpimnct = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimnct"), "pimnct");
218	ratdtpcpimnct->Draw("COLZ");
219
220	candtpcm->cd(3); preparepad();
221	TH2D ratdtpcpimtpc = plot2ddtpc((TH2D ) gDirectory->Get("NumTPCSepsqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenTPCSepsqqinvcfpimtpc"), "pimtpc");
222	ratdtpcpimtpc->Draw("COLZ");
223
224	// Make plot for pair fraction of shared hits - positive
225
226	TCanvas *cansharep = new TCanvas("cansharep","cansharep",1200,300);
227	preparepad();
228	cansharep->Divide(3,1,0.0001,0.0001);
229
230	cansharep->cd(1); preparepad();
231	TH2D ratsharepipstd = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpipstd"), "pipstd");
232	ratsharepipstd->Draw("COLZ");
233
234	cansharep->cd(2); preparepad();
235	TH2D ratsharepipnct = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpipnct"), "pipnct");
236	ratsharepipnct->Draw("COLZ");
237
238	cansharep->cd(3); preparepad();
239	TH2D ratsharepiptpc = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpiptpc"), "piptpc");
240	ratsharepiptpc->Draw("COLZ");
241
242	// Make plot for pair fraction of shared hits - negative
243
244	TCanvas *cansharem = new TCanvas("cansharem","cansharem",1200,300);
245	preparepad();
246	cansharem->Divide(3,1,0.0001,0.0001);
247
248	cansharem->cd(1); preparepad();
249	TH2D ratsharepimstd = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimstd"), "pimstd");
250	ratsharepimstd->Draw("COLZ");
251
252	cansharem->cd(2); preparepad();
253	TH2D ratsharepimnct = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimnct"), "pimnct");
254	ratsharepimnct->Draw("COLZ");
255
256	cansharem->cd(3); preparepad();
257	TH2D ratsharepimtpc = plot2dshare((TH2D ) gDirectory->Get("NumSharesqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenSharesqqinvcfpimtpc"), "pimtpc");
258	ratsharepimtpc->Draw("COLZ");
259
260	// Make signal pair fraction of shared hits - positive
261
262	TCanvas *cansharesignalp = new TCanvas("cansharesignalp","cansharesignalp",1200,300);
263	preparepad();
264	cansharesignalp->Divide(3,1,0.0001,0.0001);
265
266	cansharesignalp->cd(1); preparepad();
267	((TH2D *) gDirectory->Get("NumSharesqqinvcfpipstd"))->Draw("COLZ");
268
269	cansharesignalp->cd(2); preparepad();
270	((TH2D *) gDirectory->Get("NumSharesqqinvcfpipnct"))->Draw("COLZ");
271
272	cansharesignalp->cd(3); preparepad();
273	((TH2D *) gDirectory->Get("NumSharesqqinvcfpiptpc"))->Draw("COLZ");
274
275	// Make signal pair fraction of shared hits - negative
276
277	TCanvas *cansharesignalm = new TCanvas("cansharesignalm","cansharesignalm",1200,300);
278	preparepad();
279	cansharesignalm->Divide(3,1,0.0001,0.0001);
280
281	cansharesignalm->cd(1); preparepad();
282	((TH2D *) gDirectory->Get("NumSharesqqinvcfpimstd"))->Draw("COLZ");
283
284	cansharesignalm->cd(2); preparepad();
285	((TH2D *) gDirectory->Get("NumSharesqqinvcfpimnct"))->Draw("COLZ");
286
287	cansharesignalm->cd(3); preparepad();
288	((TH2D *) gDirectory->Get("NumSharesqqinvcfpimtpc"))->Draw("COLZ");
289
290	// Make background pair fraction of shared hits - positive
291
292	TCanvas *cansharebackp = new TCanvas("cansharebackp","cansharebackp",1200,300);
293	preparepad();
294	cansharebackp->Divide(3,1,0.0001,0.0001);
295
296	cansharebackp->cd(1); preparepad();
297	((TH2D *) gDirectory->Get("DenSharesqqinvcfpipstd"))->Draw("COLZ");
298
299	cansharebackp->cd(2); preparepad();
300	((TH2D *) gDirectory->Get("DenSharesqqinvcfpipnct"))->Draw("COLZ");
301
302	cansharebackp->cd(3); preparepad();
303	((TH2D *) gDirectory->Get("DenSharesqqinvcfpiptpc"))->Draw("COLZ");
304
305	// Make background pair fraction of shared hits - negative
306
307	TCanvas *cansharebackm = new TCanvas("cansharebackm","cansharebackm",1200,300);
308	preparepad();
309	cansharebackm->Divide(3,1,0.0001,0.0001);
310
311	cansharebackm->cd(1); preparepad();
312	((TH2D *) gDirectory->Get("DenSharesqqinvcfpimstd"))->Draw("COLZ");
313
314	cansharebackm->cd(2); preparepad();
315	((TH2D *) gDirectory->Get("DenSharesqqinvcfpimnct"))->Draw("COLZ");
316
317	cansharebackm->cd(3); preparepad();
318	((TH2D *) gDirectory->Get("DenSharesqqinvcfpimtpc"))->Draw("COLZ");
319
320	// Make plot for pair quality - positive
321
322	TCanvas *canqualp = new TCanvas("canqualp","canqualp",1200,300);
323	preparepad();
324	canqualp->Divide(3,1,0.0001,0.0001);
325
326	canqualp->cd(1); preparepad();
327	TH2D ratqualpipstd = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpipstd"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpipstd"), "pipstd");
328	ratqualpipstd->Draw("COLZ");
329
330	canqualp->cd(2); preparepad();
331	TH2D ratqualpipnct = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpipnct"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpipnct"), "pipnct");
332	ratqualpipnct->Draw("COLZ");
333
334	canqualp->cd(3); preparepad();
335	TH2D ratqualpiptpc = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpiptpc"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpiptpc"), "piptpc");
336	ratqualpiptpc->Draw("COLZ");
337
338	// Make plot for pair quality - negative
339
340	TCanvas *canqualm = new TCanvas("canqualm","canqualm",1200,300);
341	preparepad();
342	canqualm->Divide(3,1,0.0001,0.0001);
343
344	canqualm->cd(1); preparepad();
345	TH2D ratqualpimstd = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpimstd"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimstd"), "pimstd");
346	ratqualpimstd->Draw("COLZ");
347
348	canqualm->cd(2); preparepad();
349	TH2D ratqualpimnct = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpimnct"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimnct"), "pimnct");
350	ratqualpimnct->Draw("COLZ");
351
352	canqualm->cd(3); preparepad();
353	TH2D ratqualpimtpc = plot2dquality((TH2D ) gDirectory->Get("NumQualitysqqinvcfpimtpc"), (TH2D *) gDirectory->Get("DenQualitysqqinvcfpimtpc"), "pimtpc");
354	ratqualpimtpc->Draw("COLZ");
355
356	// Make signal pair quality - positive
357
358	TCanvas *canqualitysignalp = new TCanvas("canqualitysignalp","canqualitysignalp",1200,300);
359	preparepad();
360	canqualitysignalp->Divide(3,1,0.0001,0.0001);
361
362	canqualitysignalp->cd(1); preparepad();
363	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipstd"))->Draw("COLZ");
364
365	canqualitysignalp->cd(2); preparepad();
366	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpipnct"))->Draw("COLZ");
367
368	canqualitysignalp->cd(3); preparepad();
369	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpiptpc"))->Draw("COLZ");
370
371	// Make signal pair quality - negative
372
373	TCanvas *canqualitysignalm = new TCanvas("canqualitysignalm","canqualitysignalm",1200,300);
374	preparepad();
375	canqualitysignalm->Divide(3,1,0.0001,0.0001);
376
377	canqualitysignalm->cd(1); preparepad();
378	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimstd"))->Draw("COLZ");
379
380	canqualitysignalm->cd(2); preparepad();
381	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimnct"))->Draw("COLZ");
382
383	canqualitysignalm->cd(3); preparepad();
384	((TH2D *) gDirectory->Get("NumQualitysqqinvcfpimtpc"))->Draw("COLZ");
385
386	// Make plot for pair quality - positive
387
388	TCanvas *canqinvp = new TCanvas("canqinvp","canqinvp",1200,300);
389	preparepad();
390	canqinvp->Divide(3,1,0.0001,0.0001);
391
392	canqinvp->cd(1); preparepad();
393	TH1D ratqinvpipstd = getcf((TH1D ) gDirectory->Get("Numqinvcfpipstd"), (TH1D *) gDirectory->Get("Denqinvcfpipstd"), "pipstd");
394	ratqinvpipstd->Draw();
395
396	canqinvp->cd(2); preparepad();
397	TH1D ratqinvpipnct = getcf((TH1D ) gDirectory->Get("Numqinvcfpipnct"), (TH1D *) gDirectory->Get("Denqinvcfpipnct"), "pipnct");
398	ratqinvpipnct->Draw("");
399
400	canqinvp->cd(3); preparepad();
401	TH1D ratqinvpiptpc = getcf((TH1D ) gDirectory->Get("Numqinvcfpiptpc"), (TH1D *) gDirectory->Get("Denqinvcfpiptpc"), "piptpc");
402	ratqinvpiptpc->Draw("");
403
404	// Make plot for pair quality - negative
405
406	TCanvas *canqinvm = new TCanvas("canqinvm","canqinvm",1200,300);
407	preparepad();
408	canqinvm->Divide(3,1,0.0001,0.0001);
409
410	canqinvm->cd(1); preparepad();
411	TH1D ratqinvpimstd = getcf((TH1D ) gDirectory->Get("Numqinvcfpimstd"), (TH1D *) gDirectory->Get("Denqinvcfpimstd"), "pimstd");
412	ratqinvpimstd->Draw("");
413
414	canqinvm->cd(2); preparepad();
415	TH1D ratqinvpimnct = getcf((TH1D ) gDirectory->Get("Numqinvcfpimnct"), (TH1D *) gDirectory->Get("Denqinvcfpimnct"), "pimnct");
416	ratqinvpimnct->Draw("");
417
418	canqinvm->cd(3); preparepad();
419	TH1D ratqinvpimtpc = getcf((TH1D ) gDirectory->Get("Numqinvcfpimtpc"), (TH1D *) gDirectory->Get("Denqinvcfpimtpc"), "pimtpc");
420	ratqinvpimtpc->Draw("");
421
422
423	}