]>
Commit | Line | Data |
---|---|---|
390dc2ca | 1 | #if !defined(__CINT__) || defined(__MAKECINT__) |
2 | #include <TMath.h> | |
3 | #include <TROOT.h> | |
4 | #include <Riostream.h> | |
5 | #include <TCanvas.h> | |
6 | ||
7 | #include <TString.h> | |
8 | ||
9 | #include <TFile.h> | |
10 | #include <TList.h> | |
11 | #include <TH1F.h> | |
12 | #include <TH1D.h> | |
13 | #include <TF2.h> | |
14 | #include <TFitResult.h> | |
15 | #include <TFitResultPtr.h> | |
16 | #include <TH2F.h> | |
17 | #include <TH3F.h> | |
18 | #endif | |
19 | ||
20 | extern TROOT *gROOT; | |
21 | //extern TFile *gFile; | |
22 | ||
23 | ||
24 | void PostProcessCTau() { | |
25 | //----------------------------------------------------------------- | |
26 | // PostProcess the histograms produced by AliAnalysisTaskCTau* | |
27 | // Produce: | |
28 | // -- V0 particle spectra (not corrected for the feeddown) | |
29 | // -- Estimated (anti)Lambda feed-down spectra (from Xi's only) | |
30 | // -- V0 particle c*tau distributions (feeddown corrected) | |
31 | //----------------------------------------------------------------- | |
32 | ||
89808fd4 | 33 | //TString name("cTau_0090"); //centrality |
390dc2ca | 34 | //TString name("cTau_0005"); //centrality |
35 | //TString name("cTau_2040"); //centrality | |
89808fd4 | 36 | TString name("cTau_4060"); //centrality |
390dc2ca | 37 | //TString name("cTau_6080"); //centrality |
38 | //TString name("cTau_8090"); //centrality | |
39 | ||
40 | ||
41 | //+++++ Real data histograms | |
42 | TFile::Open("new/real/all/" + name + ".root"); | |
43 | TList *list=(TList *)gFile->Get(name); | |
44 | ||
45 | TH2F *fK0sSi=(TH2F*)list->FindObject("fK0sSi");fK0sSi->Sumw2(); | |
46 | TH2F *fLambdaSi=(TH2F*)list->FindObject("fLambdaSi");fLambdaSi->Sumw2(); | |
47 | TH1F *fXiSiP=(TH1F*)list->FindObject("fXiSiP"); fXiSiP->Sumw2(); | |
89808fd4 | 48 | TH2F *fLambdaBarSi= |
49 | (TH2F*)list->FindObject("fLambdaBarSi"); fLambdaBarSi->Sumw2(); | |
50 | TH1F *fXiBarSiP=(TH1F*)list->FindObject("fXiBarSiP"); fXiBarSiP->Sumw2(); | |
390dc2ca | 51 | |
52 | TH1F *fMult=(TH1F*)list->FindObject("fMult"); | |
53 | //+++++ | |
54 | ||
55 | const Double_t nEvents=fMult->Integral(); | |
56 | ||
57 | //+++++ MC histograms | |
58 | TFile::Open("new/mc_b_plus/all/" + name + "_mc.root"); | |
59 | TList *listmc=(TList *)gFile->Get(name+"_mc"); | |
60 | ||
61 | TH2F *fK0sMC=(TH2F*)listmc->FindObject("fK0sMC"); fK0sMC->Sumw2(); | |
62 | TH2F *fK0sAs=(TH2F*)listmc->FindObject("fK0sAs"); fK0sAs->Sumw2(); | |
63 | ||
64 | TH2F *fLambdaMC=(TH2F*)listmc->FindObject("fLambdaMC"); fLambdaMC->Sumw2(); | |
65 | TH2F *fLambdaAs=(TH2F*)listmc->FindObject("fLambdaAs"); fLambdaAs->Sumw2(); | |
66 | ||
89808fd4 | 67 | TH2F *fLambdaBarMC= |
68 | (TH2F*)listmc->FindObject("fLambdaBarMC"); fLambdaBarMC->Sumw2(); | |
69 | TH2F *fLambdaBarAs= | |
70 | (TH2F*)listmc->FindObject("fLambdaBarAs"); fLambdaBarAs->Sumw2(); | |
71 | ||
390dc2ca | 72 | TH1F *fXiSiPMC=(TH1F*)listmc->FindObject("fXiSiP");fXiSiPMC->Sumw2(); |
73 | TH3F *fLambdaFromXi=(TH3F*)listmc->FindObject("fLambdaFromXi"); | |
74 | fLambdaFromXi->Sumw2(); | |
89808fd4 | 75 | |
76 | TH1F *fXiBarSiPMC=(TH1F*)listmc->FindObject("fXiBarSiP");fXiBarSiPMC->Sumw2(); | |
77 | TH3F *fLambdaBarFromXiBar=(TH3F*)listmc->FindObject("fLambdaBarFromXiBar"); | |
78 | fLambdaBarFromXiBar->Sumw2(); | |
390dc2ca | 79 | //+++++ |
80 | ||
81 | ||
82 | Double_t myExp2(Double_t *v, Double_t *p); | |
83 | void Correct(TH1 *rw, const TH1 *as, const TH1 *mc); | |
8bd51e5f | 84 | void Normalise(Double_t, Double_t, Double_t, Double_t, TH1 *h); |
390dc2ca | 85 | |
89808fd4 | 86 | const Int_t iMax=3; // Number of V0 particles |
87 | const TString pnam[]={"K^{0}_{S}", "#Lambda", "#bar{#Lambda}"}; | |
88 | const Double_t brch[]={0.69, 0.64, 0.64}; | |
89 | const Double_t mass[]={0.4977, 1.115, 1.115}; | |
90 | const Double_t ctau[]={2.68, 7.89, 7.89}; | |
390dc2ca | 91 | |
92 | const TH2 *in[]={ | |
93 | fK0sSi,fK0sAs,fK0sMC, | |
89808fd4 | 94 | fLambdaSi,fLambdaAs,fLambdaMC, |
95 | fLambdaBarSi,fLambdaBarAs,fLambdaBarMC | |
390dc2ca | 96 | }; |
8bd51e5f | 97 | TH1 *fd[]={ |
98 | 0,0,0, | |
89808fd4 | 99 | fLambdaFromXi,fXiSiP,fXiSiPMC, |
100 | fLambdaBarFromXiBar,fXiBarSiP,fXiBarSiPMC | |
8bd51e5f | 101 | }; |
390dc2ca | 102 | Double_t wbx=fK0sSi->GetBinWidth(3); |
8bd51e5f | 103 | Int_t nbx=fLambdaFromXi->GetNbinsX(); |
104 | Int_t nby=fLambdaFromXi->GetNbinsY(); | |
105 | Int_t nbz=fLambdaFromXi->GetNbinsZ(); | |
106 | ||
390dc2ca | 107 | for (Int_t i=0; i<iMax; i++) { |
108 | TH2 *cr=(TH2*)in[i*3 + 0]->Clone(); | |
109 | Correct(cr, in[i*3 + 1], in[i*3 + 2]); | |
110 | ||
111 | //++++ pT spectrum | |
112 | TH1 *pt=cr->ProjectionX("_px",0,-1,"e"); | |
113 | TString ptName = pnam[i] + " p_{T} spectrum"; | |
114 | pt->SetTitle(ptName.Data()); | |
8bd51e5f | 115 | Normalise(brch[i], 2*0.5, wbx, nEvents, pt); |
116 | Double_t eipt=0., ipt=pt->IntegralAndError(1, nbx, eipt); | |
390dc2ca | 117 | |
390dc2ca | 118 | new TCanvas; |
119 | pt->Draw(); | |
120 | ||
8bd51e5f | 121 | if (i>0) { |
122 | //++++ feeddown | |
123 | TH3 *fd3=(TH3*)fd[3*i + 0]; | |
124 | TH1 *rl =(TH1*)fd[3*i + 1]; | |
125 | TH1 *mc =(TH1*)fd[3*i + 2]; | |
126 | rl->Divide(mc); | |
127 | ||
128 | for (Int_t k=1; k<=nbx; k++) { | |
129 | for (Int_t l=1; l<=nby; l++) { | |
130 | for (Int_t m=1; m<=nbz; m++) { | |
89808fd4 | 131 | Float_t c=fd3->GetBinContent(k,l,m); |
8bd51e5f | 132 | c *= rl->GetBinContent(m); |
133 | fd3->SetBinContent(k,l,m,c); | |
134 | } | |
135 | } | |
136 | } | |
137 | ||
138 | TH2 *fd2=(TH2*)fd3->Project3D("yxe"); | |
139 | TH1 *fd1=fd2->ProjectionX("_px",0,-1,"e"); | |
140 | Normalise(brch[i], 2*0.5, wbx, nEvents, fd1); | |
141 | Double_t eifd=0., ifd=fd1->IntegralAndError(1, nbx, eifd); | |
142 | ||
143 | Double_t f=ifd/ipt; | |
144 | Double_t ef=1/ipt*TMath::Sqrt(eifd*eifd + f*f*eipt*eipt); | |
145 | cout<<endl<<"Global FD correction: "<<f<<"+/-"<<ef<<endl; | |
146 | ||
147 | //new TCanvas(); | |
148 | fd1->Draw("same"); | |
149 | ||
150 | cr->Add(fd2,-1); | |
151 | } | |
89808fd4 | 152 | //continue; |
390dc2ca | 153 | |
154 | //++++ c*tau | |
155 | TF2 *f2=new TF2("myexpo2",myExp2,0.,10.,0.,100.,1+1+1+nbx); | |
156 | f2->SetParameter(0, ctau[i]); | |
157 | f2->FixParameter(1, ctau[i]); | |
158 | f2->FixParameter(2, mass[i]); | |
159 | for (Int_t p=1+1+1; p<=nbx+1+1; p++) | |
160 | f2->SetParameter(p,fK0sSi->GetBinContent(p,1)); | |
161 | ||
162 | new TCanvas; | |
163 | TFitResultPtr r=cr->Fit(f2,"SQ"); | |
164 | ||
165 | Int_t status = r; | |
166 | Double_t chi2 = r->Chi2()/r->Ndf(); | |
167 | Double_t slope = r->Parameter(0); | |
168 | Double_t error = r->ParError(0); | |
8bd51e5f | 169 | cout<<endl<<pnam[i]<<" \tstatus: "<<status<<" chi2/ndf: "<<chi2<< |
170 | "\tc*tau: "<<slope<<"+/-"<<error<<endl<<endl; | |
390dc2ca | 171 | } |
172 | ||
173 | return; | |
174 | } | |
175 | ||
176 | Double_t myExp2(Double_t *v, Double_t *p) { | |
177 | Double_t pt=v[0]; | |
178 | Double_t lt=v[1]; | |
179 | ||
180 | Double_t ctau=p[1]; | |
181 | Double_t mass=p[2]; | |
182 | Double_t ct=mass*lt/pt; | |
183 | if ((lt < 2) || (ct > 2.5*ctau)) { | |
184 | TF1::RejectPoint(); | |
185 | return 0; | |
186 | } | |
187 | ||
188 | Int_t i=pt/0.1 + 1 + 1; | |
189 | return p[i]*TMath::Exp(-ct/p[0]); | |
190 | } | |
191 | ||
192 | void Correct(TH1 *rw, const TH1 *as, const TH1 *mc) { | |
193 | TH1 *eff = (TH1*)as->Clone(); | |
194 | eff->Divide(as,mc,1,1,"b"); | |
195 | rw->Divide(eff); | |
196 | delete eff; | |
197 | } | |
198 | ||
8bd51e5f | 199 | void Normalise(Double_t br, Double_t yw, Double_t bw, Double_t ne, TH1 *pt) { |
200 | pt->Scale(1/br); // branching ratio | |
201 | pt->Scale(1/yw); // rapidity window | |
202 | pt->Scale(1/bw); // bin width | |
203 | pt->Scale(1/ne); // number of events | |
204 | } |