PWGLF/SPECTRA/PiKaPr/TestAOD/HighLevelQA/CheckNSigmaStability.C

   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