Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / final_spectra / makeSpectra.C
1 #include <TFile.h>
2 #include <TH1.h>
3 #include <TF1.h>
4 #include <TCanvas.h>
5 #include <TList.h>
6 #include <TLegend.h>
7 #include <TStyle.h>
8 #include <TMath.h>
9 #include <TSystem.h>
10 #include <TLatex.h>
11 #include <TASImage.h>
12 #include <TLine.h>
13 #include <TBox.h>
14
15 #include "my_tools.C"
16
17
18 /*
19   Info:
20   Multiplies the fraction with the charged spectrum. 
21   Applies a rapidity correction in case of K and p.
22
23   To run:
24
25   gSystem->AddIncludePath("-I../macros")
26   gSystem->AddIncludePath("-I../lib")
27   gROOT->SetMacroPath(".:../macros")
28   .L my_tools.C+
29   .L makeSpectra.C+
30
31   MakeSpectra("../corrected_fraction/fraction_7tev_b_all.root", "NOPID_pp_7TeV.root", "final_spectra_7tev_b.root");
32
33
34  */
35
36 Double_t rap_correction(Double_t* x, Double_t* par);
37
38 void MakeSpectra(const Char_t* fileName,
39                  const Char_t* fileNameCharged,
40                  const Char_t* outFileName,
41                  const Char_t* endName="PP")
42 {
43   TFile* fileData = FindFileFresh(fileName);
44
45   if(!fileData)
46     return;
47
48   TFile* fileDataCharged = FindFileFresh(fileNameCharged);
49   
50   if(!fileDataCharged)
51     return;
52   
53   const Int_t nPid = 3;
54   const Char_t* pidNames[nPid] = {"Pion", "Kaon", "Proton"}; 
55   const Char_t* titleNames[nPid] = {"#pi^{+} + #pi^{-}", "K^{+} + K^{-}", "p + #bar{p}"}; 
56   const Double_t mass[nPid] = {0.140, 0.494, 0.938};
57
58   TH1D* hPtSpectrum[nPid] = {0, 0, 0};
59   TH1D* hPtSpectrumSyst[nPid] = {0, 0, 0};
60
61   TH1D* hPtCharged = (TH1D*)fileDataCharged->Get(Form("hDnDptCharged_%s", endName));
62   TH1D* hPtChargedSyst = (TH1D*)fileDataCharged->Get(Form("hDnDptChargedSyst_%s", endName));
63
64   TF1* fRap = new TF1("fRap", rap_correction, 0.0, 50.0, 2);
65
66   for(Int_t pid = 0; pid < nPid; pid++) {
67     
68     // Pions
69     TH1D* hFraction     = (TH1D*)fileData->Get(Form("h%sFraction_%s", pidNames[pid], endName));
70     TH1D* hFractionSyst = (TH1D*)fileData->Get(Form("h%sFractionSyst_%s", pidNames[pid], endName));
71     
72     hPtSpectrum[pid] = (TH1D*)hFraction->Clone(Form("h%sSpectrum_%s", pidNames[pid], endName));
73     hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid]));
74     hPtSpectrum[pid]->Multiply(hPtCharged);
75
76     hPtSpectrumSyst[pid] = (TH1D*)hFractionSyst->Clone(Form("h%sSpectrumSyst_%s", pidNames[pid], endName));
77     hPtSpectrumSyst[pid]->SetMarkerStyle(1);
78     hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid]));
79     hPtSpectrumSyst[pid]->Multiply(hPtChargedSyst);
80     hPtSpectrumSyst[pid]->SetFillStyle(1001);
81     hPtSpectrumSyst[pid]->SetFillColor(kGray);
82
83     //
84     // Rapidity correction
85     //
86     fRap->SetParameters(0.8, mass[pid]);
87     hPtSpectrum[pid]->Divide(fRap);
88     hPtSpectrumSyst[pid]->Divide(fRap);
89   }
90
91   TFile* fileOut = new TFile(outFileName, "RECREATE");
92   for(Int_t pid = 0; pid < nPid; pid++) {
93     hPtSpectrumSyst[pid]->Write();
94     hPtSpectrum[pid]->Write();
95   }
96   fileOut->Close();
97 }
98
99 //______________________________________________________________________________
100 Double_t rap_correction(Double_t* x, Double_t* par)
101 {
102   Double_t pt = x[0];  
103
104   Double_t eta  = par[0];
105   Double_t mass = par[1];
106
107   const Double_t mt = TMath::Sqrt(pt*pt + mass*mass);
108   
109   const Double_t rap = TMath::ASinH(pt/mt*TMath::SinH(eta));
110   
111   return rap/eta;
112 }