PWGGA/GammaConv/macros/Extract_IntegratedPi0Yield.C

   1 // provided by Gamma Conversion Group, PWGGA, Kathrin Koch, kkoch@physi.uni-heidelberg.de
   2
   3 #include <fstream>
   4 #include <Riostream.h>
   5 //#include <PlottingGammaHistograms.h>
   6
   7 extern TRandom *gRandom;
   8 extern TBenchmark *gBenchmark;
   9 extern TSystem *gSystem;
  10
  11 void Extract_IntegratedPi0Yield(const char* cutSelection="", const char *inputRootFile = "AnalysisResults",const char *path = "./",const char* outputDir="./Output/",const char* outputDatFileBase= "RB-data",const char* suffix=""){
  12
  13   Bool_t writePlotFiles=kFALSE; // if one wants to write f.ex .gif files from the plots
  14   if(suffix!=""){
  15     writePlotFiles=kTRUE;
  16   }
  17
  18   gROOT->Reset();
  19   gROOT->SetStyle("Plain");
  20
  21   // rebinning
  22   Int_t RebinMass = 2;
  23
  24   Double_t BGFit_range[2] = {0.17,0.3};
  25
  26   Double_t Pi0Mass = 0.135;
  27   Double_t Pi0Mass_range[2] = {0.12,0.15};
  28   Double_t Pi0Width = 0.003;
  29   Double_t Pi0Width_range[2] = {0.001,0.007};
  30   Double_t Pi0Slope = 0.007;
  31   Double_t Pi0Slope_range[2] = {0.001,0.016};
  32
  33   Double_t Pi0FitRange[2] = {0.1,0.16};
  34
  35   Bool_t MC = kFALSE;
  36
  37
  38   // --------------------------------- end of self definitions --------------------------------
  39
  40   TString filename = Form("%s%s.root",path,inputRootFile);
  41   TFile f(filename.Data());
  42
  43   TCanvas* c_mass = new TCanvas("c_mass","",200,10,600,600);  // gives the page size
  44   c_mass->SetFillColor(0);
  45
  46
  47   TDirectory *pwg4dir =(TDirectory*)f.Get(Form("PWGGA_GammaConversion_%s",cutSelection));
  48   TList *fHistosGammaConversion = (TList*)pwg4dir->Get(Form("histogramsAliGammaConversion_%s",cutSelection));
  49   TList *fESDContainer = (TList*)fHistosGammaConversion->FindObject("ESD histograms");
  50   TList *fMCContainer = (TList*)fHistosGammaConversion->FindObject("MC histograms");
  51   TList *fBackgroundContainer = (TList*)fHistosGammaConversion->FindObject("Background histograms");
  52   TH1F *ESD_Mother_Mass= (TH1F*)fESDContainer->FindObject("ESD_Mother_InvMass");
  53
  54   TH1F *ESD_Background_Mass=(TH1F*)fBackgroundContainer->FindObject("ESD_Background_InvMass");
  55
  56   TH1F * ESD_NumberOfContributorsVtx=(TH1F*)fESDContainer->FindObject("ESD_NumberOfContributorsVtx");
  57   ESD_NumberOfContributorsVtx->SetAxisRange(1.,100.);
  58
  59   Float_t nGoodEvents = ESD_NumberOfContributorsVtx->Integral();
  60
  61   ESD_Mother_Mass->Sumw2();
  62   ESD_Background_Mass->Sumw2();
  63
  64   // for Pi0 reco efficiency
  65   Double_t nPi0_MC = 0;
  66   Double_t nPi0_MC_error = 0;
  67   if(MC){
  68     TH1D *MC_Pi0_Pt=fMCContainer->FindObject("MC_Pi0_Pt");
  69     nPi0_MC = MC_Pi0_Pt->GetEntries();
  70     nPi0_MC_error = sqrt(nPi0_MC);
  71   }
  72
  73   // get counts for reco and background within pi0 range
  74   TAxis *xaxis_reco = ESD_Mother_Mass->GetXaxis();
  75   Int_t r_1 = xaxis_reco->FindBin(BGFit_range[0]);
  76   Int_t r_2 = xaxis_reco->FindBin(BGFit_range[1]);
  77   Double_t r =  ESD_Mother_Mass->Integral(r_1,r_2); // Integral(75,125)
  78   TAxis *xaxis_back =ESD_Background_Mass->GetXaxis();
  79   Int_t b_1 = xaxis_back->FindBin(BGFit_range[0]);
  80   Int_t b_2 = xaxis_back->FindBin(BGFit_range[1]);
  81   Double_t b =  ESD_Background_Mass->Integral(b_1,b_2);
  82   Double_t norm = 1;
  83   if(b != 0) norm = r/b;
  84   ESD_Background_Mass->Sumw2();
  85   ESD_Background_Mass->Scale(norm);
  86
  87   ESD_Mother_Mass->Rebin(RebinMass);
  88   ESD_Background_Mass->Rebin(RebinMass);
  89
  90   ESD_Mother_Mass->SetTitle(Form("%s - cut: %s",inputRootFile,cutSelection));
  91   ESD_Mother_Mass->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
  92   ESD_Mother_Mass->GetXaxis()->SetRangeUser(0.0,0.3);
  93   ESD_Mother_Mass->SetLineColor(2);
  94   ESD_Mother_Mass->Draw("");
  95
  96   ESD_Background_Mass->GetXaxis()->SetRangeUser(0.0,0.3);
  97   ESD_Background_Mass->SetLineColor(4);
  98   ESD_Background_Mass->Draw("same");
  99
 100   if(writePlotFiles){
 101     c_mass->Print(Form("%sSpectra-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
 102   }
 103   c_mass->Clear();
 104
 105   TAxis *xaxis_Reco = ESD_Mother_Mass->GetXaxis();
 106   Double_t ilowmassPi0_bin = xaxis_Reco->FindBin(0.1);
 107   Double_t ihighmassPi0_bin = xaxis_Reco->FindBin(0.16);
 108   Double_t Reco_error;
 109   Double_t Reco = ESD_Mother_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Reco_error);
 110   //        cout << "Reco:   " << Reco << "      error:  " << Reco_error << endl;
 111
 112   TAxis *xaxis_Back = ESD_Background_Mass->GetXaxis();
 113   Double_t ilowmassPi0_bin = xaxis_Back->FindBin(0.1);
 114   Double_t ihighmassPi0_bin = xaxis_Back->FindBin(0.16);
 115   Double_t Back_error;
 116   Double_t Back = ESD_Background_Mass->IntegralAndError(ilowmassPi0_bin,ihighmassPi0_bin,Back_error);
 117   //        cout << "Back:   " << Back << "      error:  " << Back_error << endl;
 118
 119
 120   TH1D* Signal = (TH1D*)ESD_Mother_Mass->Clone();
 121   Signal->Sumw2();
 122   Signal->Add(ESD_Background_Mass,-1.);
 123   Signal->Draw();
 124
 125   Double_t Pi0Amplitude = Signal->GetMaximum();
 126   Double_t Pi0Amplitude_Min = Pi0Amplitude*80/100;
 127   Double_t Pi0Amplitude_Max = Pi0Amplitude*100/100;
 128
 129   TF1 *fReco = 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]);
 130
 131   // fit the peak around pi0 mass and draw
 132   fReco->SetParameter(0,Pi0Amplitude);
 133   fReco->SetParameter(1,Pi0Mass);
 134   fReco->SetParameter(2,Pi0Width);
 135   fReco->SetParameter(3,Pi0Slope);
 136
 137   fReco->SetParLimits(0,Pi0Amplitude_Min,Pi0Amplitude_Max);
 138   fReco->SetParLimits(1,Pi0Mass_range[0],Pi0Mass_range[1]);
 139   fReco->SetParLimits(2,Pi0Width_range[0],Pi0Width_range[1]);
 140   fReco->SetParLimits(3,Pi0Slope_range[0],Pi0Slope_range[1]);
 141
 142   Signal->Fit(fReco,"RME");
 143   fReco->SetLineColor(3);
 144   fReco->SetLineWidth(0.7);
 145   fReco->SetLineStyle(4);
 146   fReco->Draw("same");
 147
 148   Double_t parameter[4];
 149   fReco->GetParameters(parameter);
 150
 151   Double_t Mass = parameter[1];
 152   Double_t Mass_error = fReco->GetParError(1);
 153
 154   // calculation of FWHM + error
 155   Int_t Pi0Amplitude_fit = fReco->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 156   Int_t Pi0Amplitude_bin = fReco->GetXaxis()->FindBin(Pi0Amplitude_fit);
 157   Double_t MassFWHM1 = fReco->GetX(Pi0Amplitude_fit*0.5,0.1,0.135);
 158   Double_t MassFWHM2 = fReco->GetX(Pi0Amplitude_fit*0.5,0.135,0.16);
 159   Double_t FWHM = MassFWHM2 - MassFWHM1;
 160
 161   // + error
 162   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]);
 163   fReco_EPlus->SetParameter(0,parameter[0] + fReco->GetParError(0));
 164   fReco_EPlus->SetParameter(1,parameter[1] + fReco->GetParError(1));
 165   fReco_EPlus->SetParameter(2,parameter[2] + fReco->GetParError(2));
 166   fReco_EPlus->SetParameter(3,parameter[3] + fReco->GetParError(3));
 167
 168   Int_t Pi0Amplitude_EPlus_fit = fReco_EPlus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 169   Int_t Pi0Amplitude_EPlus_bin = fReco_EPlus->GetXaxis()->FindBin(Pi0Amplitude_EPlus_fit);
 170   Double_t MassFWHM1_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.1,0.135);
 171   Double_t MassFWHM2_EPlus = fReco_EPlus->GetX(Pi0Amplitude_EPlus_fit*0.5,0.135,0.16);
 172   Double_t FWHM_EPlus = MassFWHM2_EPlus - MassFWHM1_EPlus;
 173
 174   // - error
 175   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]);
 176   fReco_EMinus->SetParameter(0,parameter[0] - fReco->GetParError(0));
 177   fReco_EMinus->SetParameter(1,parameter[1] - fReco->GetParError(1));
 178   fReco_EMinus->SetParameter(2,parameter[2] - fReco->GetParError(2));
 179   fReco_EMinus->SetParameter(3,parameter[3] - fReco->GetParError(3));
 180
 181   Int_t Pi0Amplitude_EMinus_fit = fReco_EMinus->GetMaximum(Pi0FitRange[0],Pi0FitRange[1]);
 182   Int_t Pi0Amplitude_EMinus_bin = fReco_EMinus->GetXaxis()->FindBin(Pi0Amplitude_EMinus_fit);
 183   Double_t MassFWHM1_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.1,0.135);
 184   Double_t MassFWHM2_EMinus = fReco_EMinus->GetX(Pi0Amplitude_EMinus_fit*0.5,0.135,0.16);
 185   Double_t FWHM_EMinus = MassFWHM2_EMinus - MassFWHM1_EMinus;
 186
 187   Double_t Error1 = TMath::Abs(FWHM-FWHM_EPlus);
 188   Double_t Error2 = TMath::Abs(FWHM-FWHM_EMinus);
 189
 190   Double_t FWHM_error;
 191   if(Error1>=Error2) FWHM_error = Error1;
 192   if(Error1 < Error2) FWHM_error = Error2;
 193
 194   fstream file;
 195   file.open(Form("%s%s-%s.dat",outputDir,outputDatFileBase,inputRootFile), ios::out|ios::app);
 196   file << inputRootFile << "    " << cutSelection << "    " << nGoodEvents<< "  " << Reco << "  " << Reco_error  << "   " << Back << "  " << Back_error << "    " << Mass << "  " << Mass_error << "    " << FWHM << "  " << FWHM_error << "    " << nPi0_MC<< "    " << nPi0_MC_error << endl;
 197   file.close();
 198
 199   if(writePlotFiles){
 200     c_mass->Print(Form("%sInvMass-%s-%s.%s",outputDir,inputRootFile,cutSelection,suffix));
 201   }
 202   c_mass->Update();
 203
 204   delete fReco;
 205   delete fReco_EPlus;
 206   delete fReco_EMinus;
 207   delete c_mass;
 208 }