[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / HighLevelQA / CheckNSigmaStability.C

///////////////////////////////////////////////////////////////////////////////////////////
// CheckNSigmaStability.C (called by AODQAChecks.C)                                      //
//                                                                                       //
// Written by John Groh                                                                  //
///////////////////////////////////////////////////////////////////////////////////////////

void CheckNSigmaStability(AliSpectraAODHistoManager * hman,
			  TH1F*& TPCnsigMeanTrendPion,
			  TH1F*& TPCnsigMeanTrendKaon,
			  TH1F*& TPCnsigMeanTrendProton,
			  TH1F*& TPCnsigSigmaTrendPion,
			  TH1F*& TPCnsigSigmaTrendKaon,
			  TH1F*& TPCnsigSigmaTrendProton,
			  TH1F*& TOFnsigMeanTrendPion,
			  TH1F*& TOFnsigMeanTrendKaon,
			  TH1F*& TOFnsigMeanTrendProton,
			  TH1F*& TOFnsigSigmaTrendPion,
			  TH1F*& TOFnsigSigmaTrendKaon,
			  TH1F*& TOFnsigSigmaTrendProton,
			  Int_t runs[],
			  Int_t nRuns,
			  Int_t irun,
			  Bool_t useMC)
{
  // ranges for the projections
  const Double_t ProjTPC[2] = {0.4, 0.5};
  const Double_t ProjTOF[2] = {0.9, 1.0};

  // ranges for the gaussian fits
  const Double_t fitRangeTPC[nPart][2] = {{-3.5, 3.5},
					  {-1.4, 1.4},
					  {-2.5, 2.5}};
  const Double_t fitRangeTOF[nPart][2] = {{-2.3, 2.3},
					  {-2.0, 2.0},
					  {-2.0, 2.0}};

  //------------------------------------------------------
  //                        TPC                          -
  //------------------------------------------------------

  // canvas for printing the projections and fits to a pdf file (once per run)
  TCanvas * cNSigProjFits = new TCanvas("cNSigProjFits","cNSigProjFits");
  cNSigProjFits->Divide(3,2);

  // define the 2D nsigma histograms, the projections, and the fitting functions
  // (one each for pions, kaons, and protons)
  TH2F * nsig_TPC[nPart];
  TH1F * nsig_TPC_proj[nPart];
  TF1 * fitFuncTPC[nPart];
  for (Int_t ipart; ipart<nPart; ipart++)
    {
      nsig_TPC[ipart] = new TH2F;
      nsig_TPC_proj[ipart] = new TH1F;
    }

  // for plotting data/fit
  TGraph gDataOverFitTPC[nPart];
  TGraph gDataOverFitTOF[nPart];
  TCanvas * cDataOverFit = new TCanvas("cDataOverFit","cDataOverFit");
  cDataOverFit->Divide(2,1);
   
  for (Int_t ipart=0; ipart<nPart; ipart++)
    {
      // create projection
      nsig_TPC[ipart] = (TH2F*)((TH2F*)hman->GetNSigHistogram(Form("hHistNSig%sTPC",Particle[ipart].Data())))->Clone();
      if (useMC)
	nsig_TPC_proj[ipart] = (TH1F*)nsig_TPC[ipart]->ProjectionY(Form("TPC NsigProjection %s, mc [%.1f,%.1f]",Particle[ipart].Data(),ProjTPC[0],ProjTPC[1]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[0]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[1]));
      else
	nsig_TPC_proj[ipart] = (TH1F*)nsig_TPC[ipart]->ProjectionY(Form("TPC NsigProjection %s, data [%.1f,%.1f]",Particle[ipart].Data(),ProjTPC[0],ProjTPC[1]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[0]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[1]));

      // fit the peak of interest with a gaussian
      fitFuncTPC[ipart] = new TF1("fitFuncTPC","gaus",fitRangeTPC[ipart][0],fitRangeTPC[ipart][1]);
      nsig_TPC_proj[ipart]->Fit(fitFuncTPC[ipart],"NRLQ");
    
      // draw the projections and fits
      cNSigProjFits->cd(ipart+1);
      gPad->SetLogy();
      nsig_TPC_proj[ipart]->GetXaxis()->SetTitle(Form("TPC nsigma for %.1f GeV/c < p < %.1f GeV/c",ProjTPC[0],ProjTPC[1]));
      nsig_TPC_proj[ipart]->GetXaxis()->SetRangeUser(-10,10);
      nsig_TPC_proj[ipart]->SetStats(kFALSE);
      nsig_TPC_proj[ipart]->DrawCopy();
      //fitFuncTPC[ipart]->SetLineWidth(1);
      fitFuncTPC[ipart]->DrawCopy("same");
      TLegend * lNSigTPC = new TLegend(0.59,0.59,0.99,0.99);
      lNSigTPC->SetFillColor(0);
      lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("%ss, TPC",Particle[ipart].Data()),"");
      lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("Run %i",runs[irun]),"");
      if (useMC) lNSigTPC->AddEntry(nsig_TPC_proj[ipart],"MC","");
      else lNSigTPC->AddEntry(nsig_TPC_proj[ipart],"DATA","");
      //      lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("#chi^{2}/nDOF = %.2f",(Float_t)resultTPC->Chi2() / (Float_t)resultTPC->Ndf()),"");
      lNSigTPC->DrawClone();

      // fill gDataOverFitTPC[] graphs
      Int_t nGraphPoints = 0;
      for (Float_t iPt = fitRangeTPC[ipart][0]; iPt <= fitRangeTPC[ipart][1]; iPt += (fitRangeTPC[ipart][1] - fitRangeTPC[ipart][0])/15)
	{
	  gDataOverFitTPC[ipart].SetPoint(nGraphPoints,iPt, nsig_TPC_proj[ipart]->GetBinContent(nsig_TPC_proj[ipart]->FindBin(iPt)) / fitFuncTPC[ipart]->Eval(iPt));
	  nGraphPoints++;
	}

      // fill the histograms with the fit parameters and run numbers
      switch (ipart)
	{
	case 0: // pion
	  TPCnsigMeanTrendPion->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
	  TPCnsigSigmaTrendPion->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));	 
	  TPCnsigMeanTrendPion->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
	  TPCnsigSigmaTrendPion->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
	  break;
	case 1: // kaon
	  TPCnsigMeanTrendKaon->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
	  TPCnsigSigmaTrendKaon->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));
	  TPCnsigMeanTrendKaon->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
	  TPCnsigSigmaTrendKaon->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
	  break;
	case 2: // proton
	  TPCnsigMeanTrendProton->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
	  TPCnsigSigmaTrendProton->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));
	  TPCnsigMeanTrendProton->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
	  TPCnsigSigmaTrendProton->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
	  break;
	default:
	  Printf("\n!!! ERROR in TPC switch control structure in CheckNSigmaStability.C !!!\n");
	  break;
	  }
    } // end loop over ipart

  // draw the data/fit
  cDataOverFit->cd(1);
  TH2F * hAxesDataOverFitTPC = new TH2F("hAxesDataOverFitTPC","",100,-5,5,100,-3,5);
  hAxesDataOverFitTPC->SetStats(kFALSE);
  hAxesDataOverFitTPC->SetTitle(";nsigma for 0.4 GeV/c < p < 0.5 GeV/c;Data/Fit (TPC)");
  hAxesDataOverFitTPC->DrawCopy();
  TLegend * lDataOverFitTPC = new TLegend(.59,.69,.99,.99);
  lDataOverFitTPC->SetFillColor(0);
  lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],Form("Run %i",runs[irun]),"");
  if (useMC) lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],"MC","");
  else lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],"DATA","");
  for (Int_t ipart=0; ipart<nPart; ipart++)
    {
      gDataOverFitTPC[ipart].SetMarkerStyle(Marker[ipart]);
      gDataOverFitTPC[ipart].SetMarkerColor(Color[ipart]);
      gDataOverFitTPC[ipart].DrawClone("Psame");
      lDataOverFitTPC->AddEntry(&gDataOverFitTPC[ipart],Particle[ipart].Data(),"p");
    }
  lDataOverFitTPC->DrawClone();

  //------------------------------------------------------
  //                        TOF                          -
  //------------------------------------------------------

  // define the 2D nsigma histograms, the projections, and the fitting functions
  // (one each for pions, kaons, and protons)
  TH2F * nsig_TOF[nPart];
  TH1F * nsig_TOF_proj[nPart];
  TF1 * fitFuncTOF[nPart];
  for (Int_t ipart; ipart<nPart; ipart++)
    {
      nsig_TOF[ipart] = new TH2F;
      nsig_TOF_proj[ipart] = new TH1F;
    }
  
  for (Int_t ipart=0; ipart<nPart; ipart++)
    {
      // create projection
      nsig_TOF[ipart] = (TH2F*)((TH2F*)hman->GetNSigHistogram(Form("hHistNSig%sTOF",Particle[ipart].Data())))->Clone();
      if (useMC)
	nsig_TOF_proj[ipart] = (TH1F*)nsig_TOF[ipart]->ProjectionY(Form("TOF NsigProjection %s, mc [%.1f,%.1f]",Particle[ipart].Data(),ProjTOF[0],ProjTOF[1]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[0]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[1]));
      else
	nsig_TOF_proj[ipart] = (TH1F*)nsig_TOF[ipart]->ProjectionY(Form("TOF NsigProjection %s, data [%.1f,%.1f]",Particle[ipart].Data(),ProjTOF[0],ProjTOF[1]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[0]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[1]));
	
      // fit the peak of interest with a gaussian
      fitFuncTOF[ipart] = new TF1("fitFuncTOF","gaus",fitRangeTOF[ipart][0],fitRangeTOF[ipart][1]);
      TFitResultPtr resultTOF = nsig_TOF_proj[ipart]->Fit(fitFuncTOF[ipart],"NRSLQ");

      // draw the projections and fits
      cNSigProjFits->cd(ipart+1+nPart);
      gPad->SetLogy();

      nsig_TOF_proj[ipart]->GetXaxis()->SetTitle(Form("TOF nsigma for %.1f GeV/c < p < %.1f GeV/c",ProjTOF[0],ProjTOF[1]));
      nsig_TOF_proj[ipart]->SetStats(kFALSE);
      nsig_TOF_proj[ipart]->DrawCopy();
      fitFuncTOF[ipart]->DrawCopy("same");
      TLegend * lNSigTOF = new TLegend(0.59,0.59,0.99,0.99);
      lNSigTOF->SetFillColor(0);
      lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("%ss, TOF",Particle[ipart].Data()),"");
      lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("Run %i",runs[irun]),"");
      if (useMC) lNSigTOF->AddEntry(nsig_TOF_proj[ipart],"MC","");
      else lNSigTOF->AddEntry(nsig_TOF_proj[ipart],"DATA","");
      //lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("#chi^{2}/nDOF = %.2f",fitFuncTOF[ipart]->GetChisquare()/fitFuncTOF[ipart]->GetNDF()),"");
      lNSigTOF->DrawClone();

      // fill gDataOverFitTOF[] graphs
      Int_t nGraphPoints = 0;
      for (Float_t iPt = fitRangeTOF[ipart][0]; iPt <= fitRangeTOF[ipart][1]; iPt += (fitRangeTOF[ipart][1] - fitRangeTOF[ipart][0])/15)
	{
	  gDataOverFitTOF[ipart].SetPoint(nGraphPoints,iPt, nsig_TOF_proj[ipart]->GetBinContent(nsig_TOF_proj[ipart]->FindBin(iPt)) / fitFuncTOF[ipart]->Eval(iPt));
	  nGraphPoints++;
	}

    
      // fill the histograms with the fit parameters and run numbers
      switch (ipart)
	{
	case 0: // pion
	  TOFnsigMeanTrendPion->SetBinContent(irun+1, resultTOF->Parameter(1));
	  TOFnsigMeanTrendPion->SetBinError(irun+1, resultTOF->ParError(1));
	  TOFnsigSigmaTrendPion->SetBinContent(irun+1, resultTOF->Parameter(2));
	  TOFnsigSigmaTrendPion->SetBinError(irun+1, resultTOF->ParError(2));
	  break;
	case 1: // kaon
	  TOFnsigMeanTrendKaon->SetBinContent(irun+1, resultTOF->Parameter(1));
	  TOFnsigMeanTrendKaon->SetBinError(irun+1, resultTOF->ParError(1));
	  TOFnsigSigmaTrendKaon->SetBinContent(irun+1, resultTOF->Parameter(2));
	  TOFnsigSigmaTrendKaon->SetBinError(irun+1, resultTOF->ParError(2));
	  break;
	case 2: // proton
	  TOFnsigMeanTrendProton->SetBinContent(irun+1, resultTOF->Parameter(1));
	  TOFnsigMeanTrendProton->SetBinError(irun+1, resultTOF->ParError(1));
	  TOFnsigSigmaTrendProton->SetBinContent(irun+1, resultTOF->Parameter(2));
	  TOFnsigSigmaTrendProton->SetBinError(irun+1, resultTOF->ParError(2));
	  break;
	default:
	  Printf("\n!!! ERROR in TOF switch control structure in CheckNSigmaStability.C !!!\n");
	  break;
	}
    } // end loop over ipart

  // draw the data/fit
  cDataOverFit->cd(2);
  TH2F * hAxesDataOverFitTOF = new TH2F("hAxesDataOverFitTOF","",100,-5,5,100,-3,5);
  hAxesDataOverFitTOF->SetStats(kFALSE);
  hAxesDataOverFitTOF->SetTitle(";nsigma for 0.9 GeV/c < p < 1.0 GeV/c;Data/Fit (TOF)");
  hAxesDataOverFitTOF->DrawCopy();
  TLegend * lDataOverFitTOF = new TLegend(.59,.69,.99,.99);
  lDataOverFitTOF->SetFillColor(0);
  lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],Form("Run %i",runs[irun]),"");
  if (useMC) lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],"MC","");
  else lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],"DATA","");
  for (Int_t ipart=0; ipart<nPart; ipart++)
    {
      gDataOverFitTOF[ipart].SetMarkerStyle(Marker[ipart]);
      gDataOverFitTOF[ipart].SetMarkerColor(Color[ipart]);
      gDataOverFitTOF[ipart].DrawClone("Psame");
      lDataOverFitTOF->AddEntry(&gDataOverFitTOF[ipart],Particle[ipart].Data(),"p");
    }
  lDataOverFitTOF->DrawClone();


  // save the projections and fits to a pdf file once per run
  if (useMC)
    {
      if (irun == 0) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf(","pdf");
      else if (irun < nRuns-1) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf","pdf");
      else if (irun == nRuns-1) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf)","pdf");
    } 
  else
    {
      if (irun == 0) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf(","pdf");
      else if (irun < nRuns-1) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf","pdf");
      else if (irun == nRuns-1) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf)","pdf");
    }
  cNSigProjFits->Close();


  // also save the data/fit plots to a pdf file once per run
  if (useMC)
    {
      if (irun == 0) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf(","pdf");
      else if (irun < nRuns-1) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf","pdf");
      else if (irun == nRuns-1) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf)","pdf");
    }
  else
    {
      if (irun == 0) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf(","pdf");
      else if (irun < nRuns-1) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf","pdf");
      else if (irun == nRuns-1) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf)","pdf");
    }
  cDataOverFit->Close();
}
Commit	Line	Data
22d1b439	1	///////////////////////////////////////////////////////////////////////////////////////////
	2	// CheckNSigmaStability.C (called by AODQAChecks.C) //
	3	// //
	4	// Written by John Groh //
	5	///////////////////////////////////////////////////////////////////////////////////////////
	6
	7	void CheckNSigmaStability(AliSpectraAODHistoManager * hman,
	8	TH1F*& TPCnsigMeanTrendPion,
	9	TH1F*& TPCnsigMeanTrendKaon,
	10	TH1F*& TPCnsigMeanTrendProton,
	11	TH1F*& TPCnsigSigmaTrendPion,
	12	TH1F*& TPCnsigSigmaTrendKaon,
	13	TH1F*& TPCnsigSigmaTrendProton,
	14	TH1F*& TOFnsigMeanTrendPion,
	15	TH1F*& TOFnsigMeanTrendKaon,
	16	TH1F*& TOFnsigMeanTrendProton,
	17	TH1F*& TOFnsigSigmaTrendPion,
	18	TH1F*& TOFnsigSigmaTrendKaon,
	19	TH1F*& TOFnsigSigmaTrendProton,
	20	Int_t runs[],
	21	Int_t nRuns,
	22	Int_t irun,
	23	Bool_t useMC)
	24	{
	25	// ranges for the projections
	26	const Double_t ProjTPC[2] = {0.4, 0.5};
	27	const Double_t ProjTOF[2] = {0.9, 1.0};
	28
	29	// ranges for the gaussian fits
	30	const Double_t fitRangeTPC[nPart][2] = {{-3.5, 3.5},
	31	{-1.4, 1.4},
	32	{-2.5, 2.5}};
	33	const Double_t fitRangeTOF[nPart][2] = {{-2.3, 2.3},
	34	{-2.0, 2.0},
	35	{-2.0, 2.0}};
	36
	37	//------------------------------------------------------
	38	// TPC -
	39	//------------------------------------------------------
	40
	41	// canvas for printing the projections and fits to a pdf file (once per run)
	42	TCanvas * cNSigProjFits = new TCanvas("cNSigProjFits","cNSigProjFits");
	43	cNSigProjFits->Divide(3,2);
	44
	45	// define the 2D nsigma histograms, the projections, and the fitting functions
	46	// (one each for pions, kaons, and protons)
	47	TH2F * nsig_TPC[nPart];
	48	TH1F * nsig_TPC_proj[nPart];
	49	TF1 * fitFuncTPC[nPart];
	50	for (Int_t ipart; ipart<nPart; ipart++)
	51	{
	52	nsig_TPC[ipart] = new TH2F;
	53	nsig_TPC_proj[ipart] = new TH1F;
	54	}
	55
	56	// for plotting data/fit
	57	TGraph gDataOverFitTPC[nPart];
	58	TGraph gDataOverFitTOF[nPart];
	59	TCanvas * cDataOverFit = new TCanvas("cDataOverFit","cDataOverFit");
	60	cDataOverFit->Divide(2,1);
	61
	62	for (Int_t ipart=0; ipart<nPart; ipart++)
	63	{
	64	// create projection
65	nsig_TPC[ipart] = (TH2F)((TH2F)hman->GetNSigHistogram(Form("hHistNSig%sTPC",Particle[ipart].Data())))->Clone();
66	if (useMC)
67	nsig_TPC_proj[ipart] = (TH1F*)nsig_TPC[ipart]->ProjectionY(Form("TPC NsigProjection %s, mc [%.1f,%.1f]",Particle[ipart].Data(),ProjTPC[0],ProjTPC[1]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[0]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[1]));
68	else
69	nsig_TPC_proj[ipart] = (TH1F*)nsig_TPC[ipart]->ProjectionY(Form("TPC NsigProjection %s, data [%.1f,%.1f]",Particle[ipart].Data(),ProjTPC[0],ProjTPC[1]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[0]), nsig_TPC[ipart]->GetXaxis()->FindBin(ProjTPC[1]));
70
71	// fit the peak of interest with a gaussian
72	fitFuncTPC[ipart] = new TF1("fitFuncTPC","gaus",fitRangeTPC[ipart][0],fitRangeTPC[ipart][1]);
73	nsig_TPC_proj[ipart]->Fit(fitFuncTPC[ipart],"NRLQ");
74
75	// draw the projections and fits
76	cNSigProjFits->cd(ipart+1);
77	gPad->SetLogy();
78	nsig_TPC_proj[ipart]->GetXaxis()->SetTitle(Form("TPC nsigma for %.1f GeV/c < p < %.1f GeV/c",ProjTPC[0],ProjTPC[1]));
79	nsig_TPC_proj[ipart]->GetXaxis()->SetRangeUser(-10,10);
80	nsig_TPC_proj[ipart]->SetStats(kFALSE);
81	nsig_TPC_proj[ipart]->DrawCopy();
82	//fitFuncTPC[ipart]->SetLineWidth(1);
83	fitFuncTPC[ipart]->DrawCopy("same");
84	TLegend * lNSigTPC = new TLegend(0.59,0.59,0.99,0.99);
85	lNSigTPC->SetFillColor(0);
86	lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("%ss, TPC",Particle[ipart].Data()),"");
87	lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("Run %i",runs[irun]),"");
88	if (useMC) lNSigTPC->AddEntry(nsig_TPC_proj[ipart],"MC","");
89	else lNSigTPC->AddEntry(nsig_TPC_proj[ipart],"DATA","");
90	// lNSigTPC->AddEntry(nsig_TPC_proj[ipart],Form("#chi^{2}/nDOF = %.2f",(Float_t)resultTPC->Chi2() / (Float_t)resultTPC->Ndf()),"");
91	lNSigTPC->DrawClone();
92
93	// fill gDataOverFitTPC[] graphs
94	Int_t nGraphPoints = 0;
95	for (Float_t iPt = fitRangeTPC[ipart][0]; iPt <= fitRangeTPC[ipart][1]; iPt += (fitRangeTPC[ipart][1] - fitRangeTPC[ipart][0])/15)
96	{
97	gDataOverFitTPC[ipart].SetPoint(nGraphPoints,iPt, nsig_TPC_proj[ipart]->GetBinContent(nsig_TPC_proj[ipart]->FindBin(iPt)) / fitFuncTPC[ipart]->Eval(iPt));
98	nGraphPoints++;
99	}
100
101	// fill the histograms with the fit parameters and run numbers
102	switch (ipart)
103	{
104	case 0: // pion
105	TPCnsigMeanTrendPion->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
106	TPCnsigSigmaTrendPion->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));
107	TPCnsigMeanTrendPion->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
108	TPCnsigSigmaTrendPion->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
109	break;
110	case 1: // kaon
111	TPCnsigMeanTrendKaon->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
112	TPCnsigSigmaTrendKaon->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));
113	TPCnsigMeanTrendKaon->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
114	TPCnsigSigmaTrendKaon->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
115	break;
116	case 2: // proton
117	TPCnsigMeanTrendProton->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(1));
118	TPCnsigSigmaTrendProton->SetBinContent(irun+1, fitFuncTPC[ipart]->GetParameter(2));
119	TPCnsigMeanTrendProton->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(1));
120	TPCnsigSigmaTrendProton->SetBinError(irun+1, fitFuncTPC[ipart]->GetParError(2));
121	break;
122	default:
123	Printf("\n!!! ERROR in TPC switch control structure in CheckNSigmaStability.C !!!\n");
124	break;
125	}
126	} // end loop over ipart
127
128	// draw the data/fit
129	cDataOverFit->cd(1);
130	TH2F * hAxesDataOverFitTPC = new TH2F("hAxesDataOverFitTPC","",100,-5,5,100,-3,5);
131	hAxesDataOverFitTPC->SetStats(kFALSE);
132	hAxesDataOverFitTPC->SetTitle(";nsigma for 0.4 GeV/c < p < 0.5 GeV/c;Data/Fit (TPC)");
133	hAxesDataOverFitTPC->DrawCopy();
134	TLegend * lDataOverFitTPC = new TLegend(.59,.69,.99,.99);
135	lDataOverFitTPC->SetFillColor(0);
136	lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],Form("Run %i",runs[irun]),"");
137	if (useMC) lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],"MC","");
138	else lDataOverFitTPC->AddEntry(&gDataOverFitTPC[0],"DATA","");
139	for (Int_t ipart=0; ipart<nPart; ipart++)
140	{
141	gDataOverFitTPC[ipart].SetMarkerStyle(Marker[ipart]);
142	gDataOverFitTPC[ipart].SetMarkerColor(Color[ipart]);
143	gDataOverFitTPC[ipart].DrawClone("Psame");
144	lDataOverFitTPC->AddEntry(&gDataOverFitTPC[ipart],Particle[ipart].Data(),"p");
145	}
146	lDataOverFitTPC->DrawClone();
147
148	//------------------------------------------------------
149	// TOF -
150	//------------------------------------------------------
151
152	// define the 2D nsigma histograms, the projections, and the fitting functions
153	// (one each for pions, kaons, and protons)
154	TH2F * nsig_TOF[nPart];
155	TH1F * nsig_TOF_proj[nPart];
156	TF1 * fitFuncTOF[nPart];
157	for (Int_t ipart; ipart<nPart; ipart++)
158	{
159	nsig_TOF[ipart] = new TH2F;
160	nsig_TOF_proj[ipart] = new TH1F;
161	}
162
163	for (Int_t ipart=0; ipart<nPart; ipart++)
164	{
165	// create projection
166	nsig_TOF[ipart] = (TH2F)((TH2F)hman->GetNSigHistogram(Form("hHistNSig%sTOF",Particle[ipart].Data())))->Clone();
167	if (useMC)
168	nsig_TOF_proj[ipart] = (TH1F*)nsig_TOF[ipart]->ProjectionY(Form("TOF NsigProjection %s, mc [%.1f,%.1f]",Particle[ipart].Data(),ProjTOF[0],ProjTOF[1]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[0]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[1]));
169	else
170	nsig_TOF_proj[ipart] = (TH1F*)nsig_TOF[ipart]->ProjectionY(Form("TOF NsigProjection %s, data [%.1f,%.1f]",Particle[ipart].Data(),ProjTOF[0],ProjTOF[1]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[0]), nsig_TOF[ipart]->GetXaxis()->FindBin(ProjTOF[1]));
171
172	// fit the peak of interest with a gaussian
173	fitFuncTOF[ipart] = new TF1("fitFuncTOF","gaus",fitRangeTOF[ipart][0],fitRangeTOF[ipart][1]);
174	TFitResultPtr resultTOF = nsig_TOF_proj[ipart]->Fit(fitFuncTOF[ipart],"NRSLQ");
175
176	// draw the projections and fits
177	cNSigProjFits->cd(ipart+1+nPart);
178	gPad->SetLogy();
179
180	nsig_TOF_proj[ipart]->GetXaxis()->SetTitle(Form("TOF nsigma for %.1f GeV/c < p < %.1f GeV/c",ProjTOF[0],ProjTOF[1]));
181	nsig_TOF_proj[ipart]->SetStats(kFALSE);
182	nsig_TOF_proj[ipart]->DrawCopy();
183	fitFuncTOF[ipart]->DrawCopy("same");
184	TLegend * lNSigTOF = new TLegend(0.59,0.59,0.99,0.99);
185	lNSigTOF->SetFillColor(0);
186	lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("%ss, TOF",Particle[ipart].Data()),"");
187	lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("Run %i",runs[irun]),"");
188	if (useMC) lNSigTOF->AddEntry(nsig_TOF_proj[ipart],"MC","");
189	else lNSigTOF->AddEntry(nsig_TOF_proj[ipart],"DATA","");
190	//lNSigTOF->AddEntry(nsig_TOF_proj[ipart],Form("#chi^{2}/nDOF = %.2f",fitFuncTOF[ipart]->GetChisquare()/fitFuncTOF[ipart]->GetNDF()),"");
191	lNSigTOF->DrawClone();
192
193	// fill gDataOverFitTOF[] graphs
194	Int_t nGraphPoints = 0;
195	for (Float_t iPt = fitRangeTOF[ipart][0]; iPt <= fitRangeTOF[ipart][1]; iPt += (fitRangeTOF[ipart][1] - fitRangeTOF[ipart][0])/15)
196	{
197	gDataOverFitTOF[ipart].SetPoint(nGraphPoints,iPt, nsig_TOF_proj[ipart]->GetBinContent(nsig_TOF_proj[ipart]->FindBin(iPt)) / fitFuncTOF[ipart]->Eval(iPt));
198	nGraphPoints++;
199	}
200
201
202	// fill the histograms with the fit parameters and run numbers
203	switch (ipart)
204	{
205	case 0: // pion
206	TOFnsigMeanTrendPion->SetBinContent(irun+1, resultTOF->Parameter(1));
207	TOFnsigMeanTrendPion->SetBinError(irun+1, resultTOF->ParError(1));
208	TOFnsigSigmaTrendPion->SetBinContent(irun+1, resultTOF->Parameter(2));
209	TOFnsigSigmaTrendPion->SetBinError(irun+1, resultTOF->ParError(2));
210	break;
211	case 1: // kaon
212	TOFnsigMeanTrendKaon->SetBinContent(irun+1, resultTOF->Parameter(1));
213	TOFnsigMeanTrendKaon->SetBinError(irun+1, resultTOF->ParError(1));
214	TOFnsigSigmaTrendKaon->SetBinContent(irun+1, resultTOF->Parameter(2));
215	TOFnsigSigmaTrendKaon->SetBinError(irun+1, resultTOF->ParError(2));
216	break;
217	case 2: // proton
218	TOFnsigMeanTrendProton->SetBinContent(irun+1, resultTOF->Parameter(1));
219	TOFnsigMeanTrendProton->SetBinError(irun+1, resultTOF->ParError(1));
220	TOFnsigSigmaTrendProton->SetBinContent(irun+1, resultTOF->Parameter(2));
221	TOFnsigSigmaTrendProton->SetBinError(irun+1, resultTOF->ParError(2));
222	break;
223	default:
224	Printf("\n!!! ERROR in TOF switch control structure in CheckNSigmaStability.C !!!\n");
225	break;
226	}
227	} // end loop over ipart
228
229	// draw the data/fit
230	cDataOverFit->cd(2);
231	TH2F * hAxesDataOverFitTOF = new TH2F("hAxesDataOverFitTOF","",100,-5,5,100,-3,5);
232	hAxesDataOverFitTOF->SetStats(kFALSE);
233	hAxesDataOverFitTOF->SetTitle(";nsigma for 0.9 GeV/c < p < 1.0 GeV/c;Data/Fit (TOF)");
234	hAxesDataOverFitTOF->DrawCopy();
235	TLegend * lDataOverFitTOF = new TLegend(.59,.69,.99,.99);
236	lDataOverFitTOF->SetFillColor(0);
237	lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],Form("Run %i",runs[irun]),"");
238	if (useMC) lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],"MC","");
239	else lDataOverFitTOF->AddEntry(&gDataOverFitTOF[0],"DATA","");
240	for (Int_t ipart=0; ipart<nPart; ipart++)
241	{
242	gDataOverFitTOF[ipart].SetMarkerStyle(Marker[ipart]);
243	gDataOverFitTOF[ipart].SetMarkerColor(Color[ipart]);
244	gDataOverFitTOF[ipart].DrawClone("Psame");
245	lDataOverFitTOF->AddEntry(&gDataOverFitTOF[ipart],Particle[ipart].Data(),"p");
246	}
247	lDataOverFitTOF->DrawClone();
248
249
250
251
252	// save the projections and fits to a pdf file once per run
253	if (useMC)
254	{
255	if (irun == 0) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf(","pdf");
256	else if (irun < nRuns-1) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf","pdf");
257	else if (irun == nRuns-1) cNSigProjFits->Print("Plots/MC/AODnSigmaProjFits.pdf)","pdf");
258	}
259	else
260	{
261	if (irun == 0) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf(","pdf");
262	else if (irun < nRuns-1) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf","pdf");
263	else if (irun == nRuns-1) cNSigProjFits->Print("Plots/DATA/AODnSigmaProjFits.pdf)","pdf");
264	}
265	cNSigProjFits->Close();
266
267
268
269	// also save the data/fit plots to a pdf file once per run
270	if (useMC)
271	{
272	if (irun == 0) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf(","pdf");
273	else if (irun < nRuns-1) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf","pdf");
274	else if (irun == nRuns-1) cDataOverFit->Print("Plots/MC/NSigDataOverFit.pdf)","pdf");
275	}
276	else
277	{
278	if (irun == 0) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf(","pdf");
279	else if (irun < nRuns-1) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf","pdf");
280	else if (irun == nRuns-1) cDataOverFit->Print("Plots/DATA/NSigDataOverFit.pdf)","pdf");
281	}
282	cDataOverFit->Close();
283	}
284
285
286
287