PWGLF/STRANGENESS/LambdaK0PbPb/MainStreamAnalysis/PtMassAna2.C

   1 #include "AliMassFitControl.h"
   2
   3 // Function for fitting back ground away from peak
   4 Float_t quad(Double_t *x, Double_t *par){
   5   if (x[0] > par[3] && x[0] < par[4]){
   6     TF1::RejectPoint();
   7     return 0;
   8   }
   9   return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; //a+bx+cx**2
  10 }
  11
  12 TH1F* PtMassAna2(TH2F *PtMass, Int_t mode,Int_t ihist, const Int_t NControl, TObjArray *ControlArray,Char_t* fulllabel){
  13
  14
  15   //  gStyle->SetOptTitle(0);
  16   // gStyle->SetOptStat(0);
  17   gStyle->SetFrameBorderSize(0);
  18   gStyle->SetFrameBorderMode(0);
  19   gStyle->SetCanvasBorderMode(0);
  20   gStyle->SetCanvasColor(0);
  21
  22   //TString *tLabel = new TString(fulllabel); //not reqd
  23
  24   //Arguments TH2F *PtMass - 2D histogram of some variable versus mass
  25   //          Int_t mode - mode for switching things depending on particle
  26   //          const Int_t NControl - number of bins to project into which must
  27   //          be the same as the number of controller objects (see next arg).
  28   //          TObjArray *ContolArray - holds the controllers objects, 1 per bin
  29
  30   //const Int_t NDefault=27; // Used in default array (change as required)
  31   // Trap user errors
  32   if (ControlArray->GetEntries() != NControl){
  33     cout << "PtMassAna2: Problem: Specified number of bins and number of bins in supplied TObjArray do not match" << endl;
  34     return 0;
  35   }
  36   //Set up mode mapping to particle and names for histograms
  37   if(mode == 0){
  38     const Char_t* part = "K0";
  39     const Char_t* partName = "K^{0}_{s}";
  40   }
  41   else if (mode == 1 ){
  42     const Char_t* part = "LAM";
  43     const Char_t* partName = "#Lambda";
  44   }
  45   else if (mode == 2 ){//since lambda and anti-Lambda have same limits etc.
  46     const Char_t* part = "LAM";
  47     const Char_t* partName = "#bar{#Lambda}";
  48   }
  49   else if (mode ==3) {// This is for combined Lambda + anti-Lambda
  50     const Char_t* part = "LAM";
  51     const Char_t* partName = "#Lambda + #bar{#Lambda}";
  52   }
  53   else if (mode ==4) {
  54     const Char_t* part = "XI";
  55     const Char_t* partName = "#Xi";
  56   }
  57   else{
  58     cout << "Mode not recognized, defaulting to K0" << endl;
  59     const Char_t* part = "K0";
  60   }
  61   cout << "Debug info. Particle: " << part << endl;
  62
  63   Float_t xxlo, xxhi; //Exclusion limits for fitting
  64   // FUTURE - these can also be part of AliMassFitControl to allow different regions
  65   // in different bins
  66   Float_t NSigmaEx=4.0;  // Number of sigma to exclude from bkgd fit
  67   Float_t NSigmaInt=3.5; // Number of sigma to integrate histo over
  68
  69   // Fitting function for initial combined peak + background fit
  70   TF1 *gausQuad = new TF1("gausQuad","(([1]/([2]*pow(6.2831,0.5)))*[6]*exp(-0.5*pow((x-[0])/[2],2)))+ [3] + [4]*x + [5]*x*x",0.3,1.8);
  71   //>>>>>>>>> NB Need to set proper range when it is used <<<<<<<<<<<<<<<<<<<<<
  72   //>>>>>>>>> AND use 'R' in fit option                   <<<<<<<<<<<<<<<<<<<<<
  73   gausQuad->SetLineWidth(1);
  74   gausQuad->SetLineColor(kRed);
  75
  76   // Fitting function for background
  77   TF1* quad = new TF1("quad",quad,0.3,1.8,5);
  78   quad->SetLineWidth(2);
  79   quad->SetLineColor(kBlue);
  80   quad->SetLineStyle(2);
  81
  82   // Function for background under peak - to be integrated and drawn
  83   TF1* qback = new TF1("qback","[0]+[1]*x+[2]*x*x",0.3,1.8);
  84   qback->SetLineWidth(2);
  85   qback->SetLineColor(kRed);
  86   qback->SetLineStyle(2);
  87   qback->SetFillStyle(3013);
  88   qback->SetFillColor(kRed);
  89
  90   // Set ranges and limits according to particle
  91   Float_t xmin,xmax;
  92   Float_t defMean,defArea,defWidth,defa0,defa1,defa2,defBW; //defaults
  93   Float_t defMnLo, defMnHi;
  94   if (part=="K0"){
  95     // FUTURE: xmin, xmax define the range considered in fit so could go into
  96     // AliMassFitControl
  97     //xmin=0.418;
  98     //xmax=0.578;
  99     //    xxlo=0.46;xxhi=0.53;
 100     //    gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001*RBFact);
 101     //gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001);
 102     defMean=0.4977; defArea=1000; defWidth=0.01;
 103     defa0=2; defa1=2; defa2=0;
 104     //gausQuad->SetParLimits(0,0.483,0.523); // Constrain mean
 105     defMnLo = 0.483; defMnHi=0.523;
 106   }
 107   if (part=="LAM"){
 108     // FUTURE: See above
 109     //xmin=1.085;
 110     //xmax=1.17;
 111     //  xxlo=1.10;xxhi=1.135;
 112     //    gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001*RBFact);
 113     //  gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001);
 114     defMean=1.115; defArea=1000; defWidth=0.08;
 115     defa0=100; defa1=-10; defa2=0;
 116     //    gausQuad->SetParLimits(0,1.113,1.123); // Constrain mean
 117     defMnLo=1.113; defMnHi=1.123;
 118   }
 119   if (part=="XI"){
 120     defMean=1.32171; defArea=500; defWidth=0.08;
 121     defa0=100; defa1=-10; defa2=0;
 122     defMnLo=1.301; defMnHi=1.341;
 123   }
 124   //gausQuad->SetRange(xmin,xmax);
 125   gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million
 126   gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV
 127   gausQuad->SetParNames("Mean","Area1","Sigma1","p0","p1","p2","BinWid");
 128
 129   // Deciding how to divide up the canvas
 130   c1->Clear();
 131   c1->SetWindowSize(1024,768);
 132   Int_t NDraw = NControl;
 133   if(NDraw==16){c1->Divide(4,4);}
 134   elseif(NDraw<=36 && NDraw>25){ c1->Divide(6,6);}
 135   elseif(NDraw<=25 && NDraw>20){ c1->Divide(5,5);}
 136   elseif(NDraw<=20 && NDraw>16){  c1->Divide(5,4);}
 137   elseif(NDraw<=15 && NDraw>12){ c1->Divide(5,3);}
 138   elseif(NDraw<=12 && NDraw>8){ c1->Divide(4,3);}
 139   elseif(NDraw<=8 && NDraw>6){ c1->Divide(4,2);}
 140   elseif(NDraw<=6 && NDraw>4){c1->Divide(3,2);}
 141   elseif(NDraw<=4 && NDraw>2){c1->Divide(2,2);}
 142   elseif(NDraw==2){c1->Divide(2,1);}
 143   elseif(NDraw==1){/*do nothing*/;}
 144   else{c1->Divide(7,6);}
 145
 146   //
 147   // Project out the histograms from 2D into 1D 'slices'
 148   //
 149   Char_t id[20], title[80];
 150   //cout << "NControl: " << NControl << endl;
 151   TH1D* hMassSlice[NControl];
 152   //Arrays to store various quantities which can later be histogrammed.
 153   //  const Int_t NBinsArrays = NBins+1-NFirst; // array of 1D projection histograms is counted from 0 but other arrays are used for histograms contents and therefore N+1 are need because element zero is left empty
 154   const Int_t NBinsArrays = NControl+1;
 155   Stat_t sigBkgd[NBinsArrays], errSigBkgd[NBinsArrays]; // Signal+background and error on this
 156   Stat_t bkgd[NBinsArrays], errBkgd[NBinsArrays]; // Background and error on this
 157   Stat_t signal[NBinsArrays], errSignal[NBinsArrays]; // Net signal and error
 158
 159   Stat_t chi2PerDOF1[NBinsArrays], errChi2PerDOF1[NBinsArrays]; // Chi-squared per defree of freedom from the initial fit
 160   Stat_t chi2PerDOF2[NBinsArrays], errChi2PerDOF2[NBinsArrays]; // Chi-squared per defree of freedom from the 2nd (background) fit
 161
 162     Double_t BinCont = 0;
 163     Double_t BinErr = 0;
 164
 165
 166
 167   // REDO won't need this?
 168   // Zero the signal - avoids problems with plotting
 169   for(Int_t j=0;j<NBinsArrays;j++){
 170     signal[j]=0.0;
 171     errSignal[j]=1.0;
 172   } // <<< Not actually in use a the moment <<<
 173
 174   Stat_t mean[NBinsArrays], errMean[NBinsArrays]; // Mean and error
 175   Stat_t sigma[NBinsArrays], errSigma[NBinsArrays]; // Gaussian width and error
 176   Stat_t area[NBinsArrays], errArea[NBinsArrays]; // Peak area from fit and error
 177   Stat_t integLow[NBinsArrays], integUp[NBinsArrays]; // Lower/upper limits for integration
 178
 179   // Float_t DefaultBinPtEdges[NDefault]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
 180   Float_t BinPtEdges[NBinsArrays];
 181
 182   //  Float_t BinPtEdges[NBins+1]={0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
 183
 184   Int_t BinLo, BinHi; //For pt bins
 185   // **** Main loop over the AliMassFitControllers ****
 186
 187   TIter controlIter(ControlArray);
 188   AliMassFitControl *controller;
 189   while (controller=(AliMassFitControl*)controlIter.Next()) {
 190     if(ihist == 0)controller->CalcBinLimits(20); //This BinsPerGeV argument should be calculated from the data
 191     if(ihist == 1)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data
 192     if(ihist == 2)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data
 193     if(ihist == 3)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data
 194     if(ihist == 4)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data
 195     if(ihist == 6)controller->CalcBinLimits(20000); //This BinsPerGeV argument should be calculated from the data
 196     if(ihist == 5)controller->CalcBinLimits(100); //This BinsPerGeV argument should be calculated from the data
 197     //Had to introduce this Nint fn otherwise got inconsistencies after type implicit conversion
 198     //    BinLo=TMath::Nint(1.+BinPtEdges[N]*20.); //cout << "BinLo: " << BinLo << ", " << 1+BinPtEdges[N]*20. << endl;
 199     //    BinHi=TMath::Nint(BinPtEdges[N+1]*20.); //cout << "BinHi: " << BinHi << ", " << BinPtEdges[N+1]*20. << endl;
 200     BinLo=controller->BinLower();
 201     BinHi=controller->BinUpper();
 202     Int_t N=ControlArray->IndexOf(controller) - 1;
 203     sprintf(id,"Mass%d",N);
 204     cout << "Mass histo:" << N << " Firstbin: " << BinLo << " Last bin:" << BinHi << " " << controller->PtLower() << "-" << controller->PtUpper() << " GeV" << endl;
 205     //cout << "About to create mass projection " << N << endl;
 206     if(ihist == 0) hMassSlice[N] = PtMass->ProjectionX(id,BinLo,BinHi);
 207     if(ihist != 0) hMassSlice[N] = PtMass->ProjectionY(id,BinLo,BinHi);
 208     //cout << "Mass projection " << N << " created." << endl;
 209    sprintf(title,"%s Mass, %.2f < p_{t} < %.2f",partName,controller->PtLower(),controller->PtUpper());
 210     hMassSlice[N]->SetTitle(title);
 211     hMassSlice[N]->SetXTitle("GeV/c^{2}");
 212     // cout << "Mass projection " << N << " titles set." << endl;
 213     hMassSlice[N]->Rebin(controller->RebinFactor());
 214     // } // End of creating the projections
 215     Int_t massBinLo,massBinHi; //For mass histogram bins
 216     Int_t NArray;
 217   // Do all the fitting
 218     //REDO this is all one loop (iteration over TObjArray) now
 219   //for(Int_t N=0;N<NBins-NFirst;N++){
 220     NArray=N+1; // Arrays numbered from 1 not 0..
 221     c1->cd(N+1); //Canvas subdivisions numbered from 1 and NOT 0
 222
 223     // Set range - can be different each time
 224     xmin=controller->MinMass();
 225     xmax=controller->MaxMass();
 226     gausQuad->SetRange(xmin,xmax);
 227     // Everytime need to set some sensible parameters for both fits
 228     gausQuad->SetParameters(defMean,defArea,defWidth,defa0,defa1,defa2,0.1);
 229     quad->SetParameters(defa0,defa1,defa2,0.49,0.51);// Last 2 parameters are
 230     // the range to exclude but these are fixed later by FixParameter(3,.. etc.
 231     gausQuad->FixParameter(6,hMassSlice[N]->GetBinWidth(1));
 232     // Release Linear parameter
 233     gausQuad->ReleaseParameter(4);
 234     quad->ReleaseParameter(1);
 235     // Release Quadratic parameter
 236     gausQuad->ReleaseParameter(5);
 237     quad->ReleaseParameter(2);
 238     gausQuad->SetParLimits(0,defMnLo,defMnHi);
 239     gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million
 240     gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV
 241 //     gausQuad->SetParLimits(4,0,0); // Releases linear parameter
 242 //     quad->SetParLimits(1,0,0); // Releases linear parameter
 243 //     gausQuad->SetParLimits(5,0,0); // Releases quadratic parameter
 244 //     quad->SetParLimits(2,0,0); // Releases quadratic parameter
 245     if(controller->FixedLin()){
 246       quad->FixParameter(1,0.0);
 247       gausQuad->FixParameter(4,0.0);
 248       cout << "Info: linear part fixed to zero." << endl;
 249       gausQuad->SetParLimits(3,0.0,1000); // Ensure constant is +ve
 250       // NB documentation indicates that limits only used when option B
 251       // is specified. Also how are they can they be freed later?
 252     }
 253     if(controller->FixedQuad()){
 254       quad->FixParameter(2,0.0);
 255       gausQuad->FixParameter(5,0.0);
 256       cout << "Info: quadratic part fixed to zero." << endl;
 257     }
 258
 259     //setting the errors for histogram
 260     for(int i =1; i<= hMassSlice[N]->GetNbinsX();i++){
 261       BinCont = hMassSlice[N]->GetBinContent(i);
 262       BinErr = TMath::Sqrt(BinCont);
 263       hMassSlice[N]->SetBinError(i,BinErr);
 264     }
 265
 266      hMassSlice[N]->Fit(gausQuad,"LR");
 267      hMassSlice[N]->Draw("E,SAME");
 268      gausQuad->DrawCopy("SAME");
 269     hMassSlice[N]->Clone("hMySlice");
 270     //Fit the background:
 271     xxlo=gausQuad->GetParameter(0)-NSigmaEx*gausQuad->GetParameter(2);
 272     xxhi=gausQuad->GetParameter(0)+NSigmaEx*gausQuad->GetParameter(2);
 273     // Limits have to be adjusted to match bins
 274     massBinLo=hMassSlice[N]->FindBin(xxlo);
 275     xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo);
 276     massBinHi=hMassSlice[N]->FindBin(xxhi)+1;
 277     xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi);
 278     quad->SetParameters(gausQuad->GetParameter(3),gausQuad->GetParameter(4),gausQuad->GetParameter(5),xxlo,xxhi);
 279     quad->FixParameter(3,xxlo);
 280     quad->FixParameter(4,xxhi);
 281     quad->SetRange(xmin,xmax);
 282      hMassSlice[N]->Fit(quad,"LRN+");
 283     quad->SetFillColor(30);
 284     quad->SetFillStyle(3013);
 285     quad->SetRange(xmin,xxlo);
 286       quad->DrawCopy("SAME");
 287     quad->SetRange(xxhi,xmax);
 288          quad->DrawCopy("SAME");
 289
 290     // Draw the background - use a new function since other one not defined everywhere
 291     xxlo=gausQuad->GetParameter(0)-NSigmaInt*gausQuad->GetParameter(2);
 292     xxhi=gausQuad->GetParameter(0)+NSigmaInt*gausQuad->GetParameter(2);
 293     massBinLo=hMassSlice[N]->FindBin(xxlo);
 294     xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo);
 295     massBinHi=hMassSlice[N]->FindBin(xxhi)+1;
 296     xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi);
 297     qback->SetRange(xxlo,xxhi);
 298     qback->SetParameter(0,quad->GetParameter(0));
 299     qback->SetParameter(1,quad->GetParameter(1));
 300     qback->SetParameter(2,quad->GetParameter(2));
 301          qback->DrawCopy("SAME");
 302     if(N==3){
 303 TCanvas *myCan = new TCanvas();
 304 TPad *mypad = new TPad();
 305 myCan->cd();
 306
 307 // TLatex *myTitle = new TLatex(1.04753,20527.5,"ALICE PbPb at  #sqrt{s_{NN}} = 2.76TeV, |y|<0.75");
 308  TLatex *myTitle = new TLatex(0.4753,2527.5,"ALICE PbPb at  #sqrt{s_{NN}} = 2.76TeV, |y|<0.75");
 309  myTitle->SetTextSize(0.04);
 310  myTitle->SetTextColor(1);
 311
 312  // TLatex *myTitle1 = new TLatex(1.14634,17940.8, title);
 313  TLatex *myTitle1 = new TLatex(0.4634,2940.8, title);
 314  myTitle1->SetTextSize(0.04);
 315  myTitle1->SetTextColor(1);
 316
 317 mypad->Draw();
 318
 319
 320 hMassSlice[N]->Draw("E");
 321  hMassSlice[N]->SetXTitle("Inv. Mass [GeV/c^{2}]");
 322  hMassSlice[N]->SetYTitle("N");
 323  myTitle->Draw("same");
 324  myTitle1->Draw("same");
 325
 326 quad->SetRange(xmin,xxlo);
 327 quad->DrawCopy("SAME");
 328 quad->SetRange(xxhi,xmax);
 329 quad->DrawCopy("SAME");
 330 qback->SetRange(xxlo,xxhi);
 331 qback->DrawCopy("SAME");
 332  TFile *hfile = new TFile("06_07MassFitPID.root","RECREATE");
 333  hMassSlice[N]->Write();
 334  hfile->Close();
 335
 336 }
 337
 338
 339     //Integrate the signal+background (i.e. original histo)
 340     sigBkgd[NArray]=hMassSlice[N]->Integral(massBinLo,massBinHi);
 341     errSigBkgd[NArray]=TMath::Sqrt(sigBkgd[NArray]);
 342     //Integrate background - better to do this analytically ?
 343     bkgd[NArray]=qback->Integral(xxlo,xxhi)/gausQuad->GetParameter(6);//Divide by bin width
 344     errBkgd[NArray]=TMath::Sqrt(bkgd[NArray]); //Get errors from fit parameters?
 345
 346     // Fill arrays for diagnostic histograms
 347     mean[NArray]=gausQuad->GetParameter(0);
 348     errMean[NArray]=gausQuad->GetParError(0);
 349     area[NArray]=gausQuad->GetParameter(1);
 350     errArea[NArray]=gausQuad->GetParError(1);
 351     sigma[NArray]=gausQuad->GetParameter(2);
 352     errSigma[NArray]=gausQuad->GetParError(2);
 353     chi2PerDOF1[NArray]=gausQuad->GetChisquare()/(Double_t)gausQuad->GetNDF();
 354     errChi2PerDOF1[NArray]= 0.0; // Don't know how to calc. error for this?
 355     chi2PerDOF2[NArray]=quad->GetChisquare()/(Double_t)quad->GetNDF();
 356     errChi2PerDOF2[NArray]= 0.0; // Don't know how to calc. error for this?
 357
 358     BinPtEdges[N]=controller->PtLower();
 359     //    BinPtEdges(N+1)=controller->PtUpper();
 360   }
 361     BinPtEdges[NBinsArrays-1]=((AliMassFitControl*)ControlArray->Last())->PtUpper();
 362 //     for (Int_t jj=0;jj<NBinsArrays;jj++){
 363 //       cout << "BinPtEdges " << jj << " = " << BinPtEdges[jj] << endl;
 364 //       cout << "Mean " << jj << " = " << mean[jj] << endl;
 365 //       cout << "Sigma " << jj << " = " << sigma[jj] << endl;
 366 //       cout << "Signal + bkgd " << jj << " = " << sigBkgd[jj] << endl;
 367 //       cout << "Background " << jj << " = " << bkgd[jj] << endl;
 368 //     }
 369   // Yields in each bin returned in a histogram
 370   TH1F* hYields= new TH1F("Yield","Raw Particle Yield",NBinsArrays-1,BinPtEdges);
 371   //hYields->Reset("ICE");
 372   hYields->SetXTitle("p_{t} / (GeV/c)");
 373   //hYields->SetContent(sigBkgd);
 374   //hYields->SetError(errSigBkgd);
 375   //c1->cd(4);
 376   //hYields->Draw();
 377
 378   // Fill the diagnostic histograms: clone, set title, contents, errors.
 379   TH1F* hMeans = hYields->Clone("Means");
 380   hMeans->SetTitle("Mean from Gaussian Fit");
 381   hMeans->SetContent(mean);
 382   hMeans->SetError(errMean);
 383   //hMeans->Sumw2();
 384   TH1F* hSigmas = hYields->Clone("Sigmas");
 385   hSigmas->SetTitle("Sigma from Gaussian Fit");
 386   hSigmas->SetContent(sigma);
 387   hSigmas->SetError(errSigma);
 388   TH1F* hAreas = hYields->Clone("Areas");
 389   hAreas->SetTitle("Area from Gaussian Fit");
 390   hAreas->SetContent(area);
 391   hAreas->SetError(errArea);
 392   TH1F* hSigBkgd = hYields->Clone("SigBkgd");
 393   hSigBkgd->SetTitle("Signal + Background (Histo. Integral)");
 394   hSigBkgd->SetContent(sigBkgd);
 395   hSigBkgd->SetError(errSigBkgd);
 396   TH1F* hBkgd = hYields->Clone("Bkgd");
 397   hBkgd->SetTitle("Background (Fn. integral)");
 398   hBkgd->SetContent(bkgd);
 399   hBkgd->SetError(errBkgd);
 400
 401   hYields->Sumw2();
 402   hYields->Add(hSigBkgd,hBkgd,1.0,-1.0);
 403
 404   TH1F* hSoverB = hYields->Clone("SoverB");
 405   hSoverB->SetTitle("Signal to Background ratio");
 406   hSoverB->Sumw2();
 407   hSoverB->Divide(hBkgd);
 408
 409   TH1F* hSoverA = hYields->Clone("SoverA"); // Signal over area
 410   hSoverA->SetTitle("Ratio of Signal to Area from Fit");
 411   hSoverA->Sumw2();
 412   hSoverA->Divide(hAreas);
 413
 414   TH1F* hChi2PerDOF1 = hYields->Clone("Chi2PerDOF1");
 415   hChi2PerDOF1->SetTitle("Chi-squared per d.o.f. - initial fit");
 416   hChi2PerDOF1->SetContent(chi2PerDOF1);
 417   hChi2PerDOF1->SetError(errChi2PerDOF1);
 418
 419   TH1F* hChi2PerDOF2 = hYields->Clone("Chi2PerDOF2");
 420   hChi2PerDOF2->SetTitle("Chi-squared per d.o.f. - background fit");
 421   hChi2PerDOF2->SetContent(chi2PerDOF2);
 422   hChi2PerDOF2->SetError(errChi2PerDOF2);
 423
 424   // Draw the diagnostic histograms on their own canvas
 425   TCanvas* cDiag = new TCanvas("Diag","Diagnostics",600,600);
 426   cDiag->Divide(2,3);
 427
 428   cDiag->cd(1);
 429   Float_t bookmass;
 430   if(part=="LAM"){
 431     hMeans->SetMinimum(1.11);
 432     hMeans->SetMaximum(1.13);
 433     bookmass=1.11563;
 434   } else if (part=="K0"){
 435     hMeans->SetMaximum(0.51);
 436     hMeans->SetMinimum(0.475);
 437     bookmass=0.4976;
 438   } else if (part=="XI") {
 439     hMeans->SetMaximum(1.34);
 440     hMeans->SetMinimum(1.30);
 441     bookmass=1.32171;
 442   }
 443   TLine PDGmass;
 444   PDGmass.SetLineStyle(2);
 445   hMeans->Draw();
 446   PDGmass.DrawLine(0,bookmass,5,bookmass);
 447
 448   cDiag->cd(2);
 449   if (part=="K0"){
 450     //hSigmas->SetMaximum(0.012);
 451     hSigmas->SetMaximum(0.112);
 452   }else{
 453     hSigmas->SetMaximum(0.006);}
 454   hSigmas->SetMinimum(0.0);
 455   hSigmas->Draw();
 456
 457   cDiag->cd(3);
 458   hSoverB->SetMinimum(0.0);
 459   hSoverB->Draw();
 460
 461   cDiag->cd(4);
 462   hSoverA->SetMinimum(0.5);
 463   hSoverA->SetMaximum(1.5);
 464   hSoverA->Draw();
 465
 466   cDiag->cd(5);
 467   hChi2PerDOF1->SetMinimum(0);
 468   hChi2PerDOF1->SetMaximum(6);
 469   hChi2PerDOF1->Draw();
 470
 471   cDiag->cd(6);
 472   hChi2PerDOF2->SetMinimum(0);
 473   hChi2PerDOF2->SetMaximum(60);
 474   hChi2PerDOF2->Draw();
 475
 476   Char_t fileNameBase[80];
 477   sprintf(fileNameBase,"Diagnostics%s_%s",part,fulllabel);
 478   Char_t fileNamePng[80];
 479   sprintf(fileNamePng,"%s.png",fileNameBase);
 480   Char_t fileNameEps[80];
 481   sprintf(fileNameEps,"%s.eps",fileNameBase);
 482   Char_t fileNamePdf[80];
 483   sprintf(fileNamePdf,"%s.pdf",fileNameBase);
 484
 485   cDiag->SaveAs(fileNamePng);
 486   cDiag->SaveAs(fileNameEps);
 487   cDiag->SaveAs(fileNamePdf);
 488
 489   // Scale result by the bin size to get dN/dpT and by bin centre to get 1/pt...
 490   Float_t y, ey;
 491   for(Int_t j=0;j<=hYields->GetNbinsX();j++){
 492     y = hYields->GetBinContent(j)/hYields->GetBinWidth(j);
 493     ey = hYields->GetBinError(j)/hYields->GetBinWidth(j);
 494     hYields->SetBinContent(j,y);
 495     hYields->SetBinError(j,ey);
 496 //    y = hYields->GetBinContent(j)/hYields->GetBinCenter(j);
 497 //    ey = hYields->GetBinError(j)/hYields->GetBinCenter(j);
 498     hYields->SetBinContent(j,y);
 499     hYields->SetBinError(j,ey);
 500   }
 501   return hYields;
 502
 503 }