PWGGA/GammaConv/macros/Extract_Pi0_Characteristics.C

   1
   2 #include <fstream>
   3 #include <Riostream.h>
   4 /*
   5  *
   6  *
   7  *
   8  */
   9
  10 extern TRandom *gRandom;
  11 extern TBenchmark *gBenchmark;
  12 extern TSystem *gSystem;
  13
  14 void Extract_Pi0_Characteristics(const char *cutSelection="", const char *inputRootFile = "AnalysisResults",const char *path = "./",const char*outputDir = "./Output/",Bool_t makeMappingPlots=kTRUE){
  15
  16   gROOT->Reset();
  17   gROOT->SetStyle("Plain");
  18   gStyle->SetOptFit(0);
  19   // rebinning
  20   Int_t RebinMass = 2;
  21
  22   Double_t BGFit_range[2] = {0.17,0.3};
  23
  24   Double_t Pi0Mass = 0.135;
  25   Double_t Pi0Mass_range[2] = {0.12,0.15};
  26   Double_t Pi0Width = 0.003;
  27   Double_t Pi0Width_range[2] = {0.001,0.007};
  28   Double_t Pi0Slope = 0.007;
  29   Double_t Pi0Slope_range[2] = {0.001,0.016};
  30
  31   Double_t Pi0FitRange[2] = {0.1,0.16};
  32
  33   // pt bins
  34   Int_t NBinsPt = 31;
  35   Double_t BinsPt[] = {0.0, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4., 5., 6.,8.,10.,13.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29.,30.};
  36
  37   // input file
  38   TString filename = Form("%s%s.root",path,inputRootFile);
  39
  40   TFile f(filename.Data());
  41
  42   TFile *outputFile;
  43   if(makeMappingPlots==kTRUE){
  44     outputFile = new TFile(Form("%sPi0Characteristics.root",outputDir),"UPDATE");
  45   }
  46   TDirectory *pwg4dir =(TDirectory*)f.Get(Form("PWGGA_GammaConversion_%s",cutSelection));
  47   TList *fHistosGammaConversion = (TList*)pwg4dir->Get(Form("histogramsAliGammaConversion_%s",cutSelection));
  48   TList *fESDContainer = (TList*)fHistosGammaConversion->FindObject("ESD histograms");
  49   TList *fMCContainer = (TList*)fHistosGammaConversion->FindObject("MC histograms");
  50   TList *fBackgroundContainer = (TList*)fHistosGammaConversion->FindObject("Background histograms");
  51
  52   TH1F *ESD_Mother_InvMass =fESDContainer->FindObject("ESD_Mother_InvMass");
  53   ESD_Mother_InvMass->SetName(Form("ESD_Mother_InvMass_%s",cutSelection));
  54   ESD_Mother_InvMass->Write();
  55
  56   TH1F *ESD_Background_InvMass =fBackgroundContainer->FindObject("ESD_Background_InvMass");
  57   ESD_Background_InvMass->SetName(Form("ESD_Background_InvMass_%s",cutSelection));
  58   ESD_Background_InvMass->Write();
  59
  60   TH1F *ESD_dedx =fESDContainer->FindObject("ESD_ConvGamma_E_dEdxP");
  61   ESD_dedx->SetName(Form("ESD_ConvGamma_E_dEdxP_%s",cutSelection));
  62   ESD_dedx->Write();
  63
  64   TH2F *ESD_alfaqt = fESDContainer->FindObject("ESD_ConvGamma_alfa_qt");
  65   ESD_alfaqt->SetName(Form("ESD_ConvGamma_alfa_qt_%s",cutSelection));
  66   ESD_alfaqt->Write();
  67
  68
  69   TH2F *ESD_Mother_InvMass_vs_Pt = fESDContainer->FindObject("ESD_Mother_InvMass_vs_Pt");
  70   ESD_Mother_InvMass_vs_Pt->Sumw2();
  71   TH2F *ESD_Background_InvMass_vs_Pt = fBackgroundContainer->FindObject("ESD_Background_InvMass_vs_Pt");
  72   ESD_Background_InvMass_vs_Pt->Sumw2();
  73
  74   //histos for comparison
  75   TH1F *histoNormYield_Pi0 = new TH1F(Form("Norm_Yield_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  76   TH1F *histoRawYield_Pi0 = new TH1F(Form("Raw_Yield_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  77   TH1F *histoMass_Pi0 = new TH1F(Form("Mass_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  78   TH1F *histoFWHM_Pi0 = new TH1F(Form("FWHM_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  79   TH1F *histoSignificance_Pi0 = new TH1F(Form("Significance_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  80   TH1F *histoSB_Pi0 = new TH1F(Form("SB_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  81   TH1F *histoSBRatio_Pi0 = new TH1F(Form("SBRatio_Pi0_%s",cutSelection),"",NBinsPt,BinsPt);
  82
  83   // for yield correction in pt bins
  84   TH1F *deltaPt = new TH1F("deltaPt","",NBinsPt,BinsPt);
  85
  86   TH1F * ESD_NumberOfContributorsVtx = fESDContainer->FindObject("ESD_NumberOfContributorsVtx");
  87   ESD_NumberOfContributorsVtx->SetAxisRange(1.,100.);
  88   Float_t nGoodEvents = ESD_NumberOfContributorsVtx->Integral();
  89   cout <<"# good events: "<< nGoodEvents << endl;
  90
  91
  92   // binning / max pt from analysis
  93   Double_t maxPt = ESD_Mother_InvMass_vs_Pt->GetYaxis()->GetBinCenter(ESD_Mother_InvMass_vs_Pt->GetNbinsY()) + ESD_Mother_InvMass_vs_Pt->GetYaxis()->GetBinWidth(ESD_Mother_InvMass_vs_Pt->GetNbinsY())/2.;
  94   Double_t binningPt = ESD_Mother_InvMass_vs_Pt->GetNbinsY();
  95
  96
  97   // ---------------------------------------------------------------------------------------------
  98
  99   // extract yield, mass, significance and FWHM in pt bins
 100
 101   TH1D* Mapping_Reco_InvMass_PtBin[32];    // aray of histos for pt slices
 102   TH1D* Mapping_Back_InvMass_PtBin[32];
 103   TH1D* Mapping_Signal_InvMass_PtBin[32];
 104   TH1D* Mapping_SBRatio_InvMass_PtBin[32];
 105
 106   TCanvas *canvas = new TCanvas("canvas","",200,10,600,600);  // gives the page size
 107   canvas->SetFillColor(0);
 108
 109   for(Int_t iPt = 2; iPt < NBinsPt+1; iPt++){
 110     TString histonameReco = Form("Mapping_Reco_InvMass_in_Pt_Bin%s%02d",cutSelection ,iPt);
 111     TString histonameBack = Form("Mapping_Back_InvMass_in_Pt_Bin%s%02d",cutSelection, iPt);
 112
 113     Mapping_Reco_InvMass_PtBin[iPt]=new TH1D(histonameReco.Data(),histonameReco.Data(),ESD_Mother_InvMass_vs_Pt->GetNbinsX(),0.,1.);
 114     Mapping_Back_InvMass_PtBin[iPt]=new TH1D(histonameBack.Data(),histonameBack.Data(),ESD_Background_InvMass_vs_Pt->GetNbinsX(),0.,1.);
 115
 116     Int_t startBin = ESD_Mother_InvMass_vs_Pt->GetYaxis()->FindBin(BinsPt[iPt-1])-1;
 117     Int_t endBin = ESD_Mother_InvMass_vs_Pt->GetYaxis()->FindBin(BinsPt[iPt])-1;
 118
 119     Double_t startpt = startBin*maxPt/binningPt;
 120     Double_t endpt = endBin*maxPt/binningPt;
 121
 122     ESD_Mother_InvMass_vs_Pt->ProjectionX(histonameReco.Data(),startBin,endBin);
 123     ESD_Background_InvMass_vs_Pt->ProjectionX(histonameBack.Data(),startBin,endBin);
 124
 125     Mapping_Reco_InvMass_PtBin[iPt]=(TH1D*)gDirectory->Get(histonameReco.Data());
 126     Mapping_Reco_InvMass_PtBin[iPt]->Rebin(RebinMass);
 127     Mapping_Back_InvMass_PtBin[iPt]=(TH1D*)gDirectory->Get(histonameBack.Data());
 128     Mapping_Back_InvMass_PtBin[iPt]->Rebin(RebinMass);
 129
 130
 131     // normalisation of background
 132     TAxis *xaxis_reco = Mapping_Reco_InvMass_PtBin[iPt]->GetXaxis();
 133     Int_t r_1 = xaxis_reco->FindBin(BGFit_range[0]);
 134     Int_t r_2 = xaxis_reco->FindBin(BGFit_range[1]);
 135     Double_t r =  Mapping_Reco_InvMass_PtBin[iPt]->Integral(r_1,r_2); // Integral(75,125)
 136     TAxis *xaxis_back = Mapping_Back_InvMass_PtBin[iPt]->GetXaxis();
 137     Int_t b_1 = xaxis_back->FindBin(BGFit_range[0]);
 138     Int_t b_2 = xaxis_back->FindBin(BGFit_range[1]);
 139     Double_t b =  Mapping_Back_InvMass_PtBin[iPt]->Integral(b_1,b_2);
 140     Double_t norm = 1;
 141     if(b != 0) norm = r/b;
 142     Mapping_Back_InvMass_PtBin[iPt]->Sumw2();
 143     Mapping_Back_InvMass_PtBin[iPt]->Scale(norm);
 144
 145     TString histonameSignal = Form("Mapping_Signal_InvMass_in_Pt_Bin%s%02d",cutSelection, iPt);
 146     Mapping_Signal_InvMass_PtBin[iPt] = (TH1D*)Mapping_Reco_InvMass_PtBin[iPt]->Clone();
 147     Mapping_Signal_InvMass_PtBin[iPt]->SetName(histonameSignal.Data());
 148     Mapping_Signal_InvMass_PtBin[iPt]->SetTitle(histonameSignal.Data());
 149     Mapping_Signal_InvMass_PtBin[iPt]->Sumw2();
 150     Mapping_Signal_InvMass_PtBin[iPt]->Add(Mapping_Back_InvMass_PtBin[iPt],-1.);
 151
 152     /*    TString histonameSBRatio= Form("Mapping_SBRatio_InvMass_in_Pt_Bin%02d", iPt);
 153     Mapping_SBRatio_InvMass_PtBin[iPt]=->ProjectionX(histonameSBRatio.Data(),startBin,endBin);;
 154     Mapping_SBRatio_InvMass_PtBin[iPt]->Add(Mapping_Reco_InvMass_PtBin[iPt],1);
 155     Mapping_SBRatio_InvMass_PtBin[iPt]->Divide(Mapping_Back_InvMass_PtBin[iPt]);
 156     Mapping_SBRatio_InvMass_PtBin[iPt]->Write();
 157     */
 158     Int_t bins =  Mapping_Signal_InvMass_PtBin[iPt]->GetNbinsX();
 159     for(Int_t ib = 1; ib < bins; ib++){
 160       Double_t error = Mapping_Signal_InvMass_PtBin[iPt]->GetBinError(ib);
 161       if(abs(error) < 0.1) error = Mapping_Signal_InvMass_PtBin[iPt]->GetBinError(ib) +1;
 162       Mapping_Signal_InvMass_PtBin[iPt]->SetBinError(ib,error);
 163     }
 164
 165     Double_t Pi0Amplitude = Mapping_Signal_InvMass_PtBin[iPt]->GetMaximum();
 166     Double_t Pi0Amplitude_Min = Pi0Amplitude*80/100;
 167     Double_t Pi0Amplitude_Max = Pi0Amplitude*100/100;
 168
 169     TAxis *xaxis = Mapping_Signal_InvMass_PtBin[iPt]->GetXaxis();
 170     Int_t ilowmassPi0_bin = xaxis->FindBin(Pi0FitRange[0]);
 171     Int_t ihighmassPi0_bin =xaxis->FindBin(Pi0FitRange[1]);
 172     if(makeMappingPlots){
 173       Mapping_Signal_InvMass_PtBin[iPt]->Write();
 174     }
 175     TF1 *fReco = new TF1(Form("Exp+Gauss PtBin[%02d]",iPt),"(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
 176
 177     // fit the peak around pi0 mass and draw
 178     fReco->SetParameter(0,Pi0Amplitude);
 179     fReco->SetParameter(1,Pi0Mass);
 180     fReco->SetParameter(2,Pi0Width);
 181     fReco->SetParameter(3,Pi0Slope);
 182
 183     fReco->SetParLimits(0,Pi0Amplitude_Min,Pi0Amplitude_Max);
 184     fReco->SetParLimits(1,Pi0Mass_range[0],Pi0Mass_range[1]);
 185     fReco->SetParLimits(2,Pi0Width_range[0],Pi0Width_range[1]);
 186     fReco->SetParLimits(3,Pi0Slope_range[0],Pi0Slope_range[1]);
 187
 188     Mapping_Signal_InvMass_PtBin[iPt]->Fit(fReco,"RME");
 189     //    Mapping_Signal_InvMass_PtBin[iPt]->Draw();
 190     fReco->SetLineColor(3);
 191     fReco->SetLineWidth(0.7);
 192     fReco->SetLineStyle(4);
 193     fReco->Draw("same");
 194     if(makeMappingPlots){
 195       fReco->Write();
 196     }
 197     cout<<"Kenneth: Setting Bin content for bin: "<<iPt<<"  to "<<fReco->GetParameter(1)<<endl;
 198     histoMass_Pi0->SetBinContent(iPt, fReco->GetParameter(1));
 199     cout<<"Kenneth: Setting Bin content for bin: "<<iPt<<"  to "<<fReco->GetParError(1)<<endl;
 200     histoMass_Pi0->SetBinError(iPt, fReco->GetParError(1));
 201
 202     TAxis *xaxis = Mapping_Signal_InvMass_PtBin[iPt]->GetXaxis();
 203     Int_t ilowmassPi0_bin = xaxis->FindBin(Pi0FitRange[0]);
 204     Int_t ihighmassPi0_bin =xaxis->FindBin(Pi0FitRange[1]);
 205
 206     Double_t Reco_error;
 207     Double_t Background_error;
 208
 209     Double_t IntReco = Mapping_Reco_InvMass_PtBin[iPt]->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Reco_error);
 210     Double_t IntBackground = Mapping_Back_InvMass_PtBin[iPt]->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Background_error);
 211     if(IntBackground>0){
 212       Double_t SB= IntReco/IntBackground;
 213       Double_t SB_error = sqrt( pow(Reco_error/IntBackground,2) + pow(IntReco/IntBackground/IntBackground*Background_error,2));
 214
 215       histoSB_Pi0->SetBinContent(iPt,SB);
 216       histoSB_Pi0->SetBinError(iPt, SB_error);
 217     }
 218     //    canvas->Update();
 219     //    canvas->Clear();
 220
 221     Double_t parameter[4];
 222     fReco->GetParameters(parameter);
 223
 224     // calculation of FWHM + error
 225     Int_t Pi0Amplitude_fit = fReco->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 226     Int_t Pi0Amplitude_bin = fReco->GetXaxis()->FindBin(Pi0Amplitude_fit);
 227     Double_t MassFWHM1 = fReco->GetX(Pi0Amplitude_fit*0.5,0.1,0.135);
 228     Double_t MassFWHM2 = fReco->GetX(Pi0Amplitude_fit*0.5,0.135,0.16);
 229     Double_t FWHM = MassFWHM2 - MassFWHM1;
 230
 231     // + error
 232     TF1 *fReco_EPlus = new TF1("Exp+Gauss","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
 233     fReco_EPlus->SetParameter(0,parameter[0] + fReco->GetParError(0));
 234     fReco_EPlus->SetParameter(1,parameter[1] + fReco->GetParError(1));
 235     fReco_EPlus->SetParameter(2,parameter[2] + fReco->GetParError(2));
 236     fReco_EPlus->SetParameter(3,parameter[3] + fReco->GetParError(3));
 237
 238     Int_t Pi0Amplitude_EPlus_fit = fReco_EPlus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 239     Int_t Pi0Amplitude_EPlus_bin = fReco_EPlus->GetXaxis()->FindBin(Pi0Amplitude_EPlus_fit);
 240
 241     Double_t MassFWHM1_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.1,0.135);
 242     Double_t MassFWHM2_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.135,0.16);
 243     Double_t FWHM_EPlus = MassFWHM2_EPlus - MassFWHM1_EPlus;
 244
 245     // - error
 246     TF1 *fReco_EMinus = new TF1("Exp+Gauss","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",Pi0FitRange[0],Pi0FitRange[1]);
 247     fReco_EMinus->SetParameter(0,parameter[0] - fReco->GetParError(0));
 248     fReco_EMinus->SetParameter(1,parameter[1] - fReco->GetParError(1));
 249     fReco_EMinus->SetParameter(2,parameter[2] - fReco->GetParError(2));
 250     fReco_EMinus->SetParameter(3,parameter[3] - fReco->GetParError(3));
 251
 252     Int_t Pi0Amplitude_EMinus_fit = fReco_EMinus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 253     Int_t Pi0Amplitude_EMinus_bin = fReco_EMinus->GetXaxis()->FindBin(Pi0Amplitude_EMinus_fit);
 254     Double_t MassFWHM1_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.1,0.135);
 255     Double_t MassFWHM2_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.135,0.16);
 256     Double_t FWHM_EMinus = MassFWHM2_EMinus - MassFWHM1_EMinus;
 257
 258     Double_t Error1 = TMath::Abs(FWHM-FWHM_EPlus);
 259     Double_t Error2 = TMath::Abs(FWHM-FWHM_EMinus);
 260
 261     Double_t FWHM_error;
 262     if(Error1>=Error2) FWHM_error = Error1;
 263     if(Error1 < Error2) FWHM_error = Error2;
 264
 265     histoFWHM_Pi0->SetBinContent(iPt, FWHM);
 266     histoFWHM_Pi0->SetBinError(iPt,FWHM_error);
 267
 268     Double_t Reco_error;
 269     Double_t Background_error;
 270
 271     Double_t Reco = Mapping_Reco_InvMass_PtBin[iPt]->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Reco_error);
 272     Double_t Background = Mapping_Back_InvMass_PtBin[iPt]->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Background_error);
 273
 274     if(Background>0){
 275       Double_t Significance = Reco/sqrt(Background);
 276       Double_t Significance_error = sqrt( pow(Reco_error/sqrt(Background),2) + pow( -0.5 * Reco/pow(Background,1.5) * Background_error,2) );
 277
 278       histoSignificance_Pi0->SetBinContent(iPt,Significance);
 279       histoSignificance_Pi0->SetBinError(iPt, Significance_error);
 280     }
 281     histoNormYield_Pi0->SetBinContent(iPt,(Reco-Background)/nGoodEvents);
 282     histoNormYield_Pi0->SetBinError(iPt,sqrt(Reco_error*Reco_error + Background_error*Background_error)/nGoodEvents);
 283
 284     histoRawYield_Pi0->SetBinContent(iPt,Reco-Background);
 285     histoRawYield_Pi0->SetBinError(iPt,sqrt(Reco_error*Reco_error + Background_error*Background_error));
 286
 287
 288     deltaPt->SetBinContent(iPt, endpt-startpt);
 289     deltaPt->SetBinError(iPt,0);
 290     if(iPt == 2){
 291       deltaPt->SetBinContent(1,startpt-0);
 292       deltaPt->SetBinError(1,0);
 293     }
 294   }
 295
 296   delete [] Mapping_Reco_InvMass_PtBin;
 297   delete [] Mapping_Back_InvMass_PtBin;
 298   delete [] Mapping_Signal_InvMass_PtBin;
 299   delete [] Mapping_SBRatio_InvMass_PtBin;
 300   // end of extract yield, mass .....
 301
 302   // write the created histos into the output root file
 303
 304   if(makeMappingPlots == kFALSE){// open the root file here so the mapping histograms is not included
 305     outputFile = new TFile(Form("%sPi0Characteristics.root",outputDir),"UPDATE");
 306   }
 307
 308
 309   //    ESD_Mother_InvMass_vs_Pt->ProjectionX(histonameReco.Data(),startBin,endBin);
 310   //    ESD_Background_InvMass_vs_Pt->ProjectionX(histonameBack.Data(),startBin,endBin);
 311   /*
 312   histoSBRatio_Pi0->Add(ESD_Mother_InvMass_vs_Pt,1);
 313   histoSBRatio_Pi0->Divide(ESD_Background_InvMass_vs_Pt);
 314   histoSBRatio_Pi0->Write();
 315   */
 316   histoNormYield_Pi0->SetXTitle("p_{t} (GeV/c)");
 317   histoNormYield_Pi0->Divide(deltaPt);
 318   histoNormYield_Pi0->SetYTitle("#pi^{0} raw yield/per event");
 319   histoNormYield_Pi0->DrawCopy("e1");
 320   histoNormYield_Pi0->Write();
 321   //  outputFile->Write();
 322   canvas->Update();
 323   canvas->Clear();
 324
 325   histoRawYield_Pi0->SetXTitle("p_{t} (GeV/c)");
 326   histoRawYield_Pi0->Divide(deltaPt);
 327   histoRawYield_Pi0->SetYTitle("#pi^{0} raw yield");
 328   histoRawYield_Pi0->DrawCopy("e1");
 329   histoRawYield_Pi0->Write();
 330   //  outputFile->Write();
 331   canvas->Update();
 332   canvas->Clear();
 333
 334   histoMass_Pi0->SetXTitle("p_{t} (GeV/c)");
 335   //  histoMass_Pi0->Divide(deltaPt);
 336   histoMass_Pi0->SetYTitle("#pi^{0} mass");
 337   histoMass_Pi0->DrawCopy("e1");
 338   histoMass_Pi0->Write();
 339   canvas->Update();
 340   canvas->Clear();
 341
 342
 343   histoSB_Pi0->SetXTitle("p_{t} (GeV/c)");
 344   //  histoMass_Pi0->Divide(deltaPt);
 345   histoSB_Pi0->SetYTitle("S/B");
 346   histoSB_Pi0->DrawCopy("e1");
 347   histoSB_Pi0->Write();
 348   canvas->Update();
 349   canvas->Clear();
 350
 351
 352   histoFWHM_Pi0->SetXTitle("p_{t} (GeV/c)");
 353   //  histoFWHM_Pi0->Divide(deltaPt);
 354   histoFWHM_Pi0->SetYTitle("#pi^{0} FWHM/2.36");
 355   histoFWHM_Pi0->DrawCopy("e1");
 356   histoFWHM_Pi0->Write();
 357   canvas->Update();
 358   canvas->Clear();
 359
 360   histoSignificance_Pi0->SetXTitle("p_{t} (GeV/c)");
 361   //  histoSignificance_Pi0->Divide(deltaPt);
 362   histoSignificance_Pi0->SetYTitle("Significance");
 363   histoSignificance_Pi0->DrawCopy("e1");
 364   histoSignificance_Pi0->Write();
 365   canvas->Update();
 366   canvas->Clear();
 367
 368   histoSBRatio_Pi0->Delete();
 369   histoNormYield_Pi0->Delete();
 370   histoRawYield_Pi0->Delete();
 371   histoMass_Pi0->Delete();
 372   histoFWHM_Pi0->Delete();
 373   histoSignificance_Pi0->Delete();
 374   histoSB_Pi0->Delete();
 375   deltaPt->Delete();
 376
 377   canvas->Delete();
 378
 379   outputFile->Write();
 380   outputFile->Close();
 381 }