Fix: Correcting the feeddown for the efficiency
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / PostProcessCTau.C
CommitLineData
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
20extern TROOT *gROOT;
21//extern TFile *gFile;
22
23
24void 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
4a407824 33 TString name("cTau_0090"); //centrality
390dc2ca 34 //TString name("cTau_0005"); //centrality
35 //TString name("cTau_2040"); //centrality
4a407824 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");
4a407824 139 Correct(fd2, in[i*3 + 1], in[i*3 + 2]);
140
8bd51e5f 141 TH1 *fd1=fd2->ProjectionX("_px",0,-1,"e");
142 Normalise(brch[i], 2*0.5, wbx, nEvents, fd1);
143 Double_t eifd=0., ifd=fd1->IntegralAndError(1, nbx, eifd);
144
145 Double_t f=ifd/ipt;
146 Double_t ef=1/ipt*TMath::Sqrt(eifd*eifd + f*f*eipt*eipt);
147 cout<<endl<<"Global FD correction: "<<f<<"+/-"<<ef<<endl;
148
149 //new TCanvas();
150 fd1->Draw("same");
151
152 cr->Add(fd2,-1);
153 }
89808fd4 154 //continue;
390dc2ca 155
156 //++++ c*tau
157 TF2 *f2=new TF2("myexpo2",myExp2,0.,10.,0.,100.,1+1+1+nbx);
158 f2->SetParameter(0, ctau[i]);
159 f2->FixParameter(1, ctau[i]);
160 f2->FixParameter(2, mass[i]);
161 for (Int_t p=1+1+1; p<=nbx+1+1; p++)
162 f2->SetParameter(p,fK0sSi->GetBinContent(p,1));
163
164 new TCanvas;
165 TFitResultPtr r=cr->Fit(f2,"SQ");
166
167 Int_t status = r;
168 Double_t chi2 = r->Chi2()/r->Ndf();
169 Double_t slope = r->Parameter(0);
170 Double_t error = r->ParError(0);
8bd51e5f 171 cout<<endl<<pnam[i]<<" \tstatus: "<<status<<" chi2/ndf: "<<chi2<<
172 "\tc*tau: "<<slope<<"+/-"<<error<<endl<<endl;
390dc2ca 173 }
174
175 return;
176}
177
178Double_t myExp2(Double_t *v, Double_t *p) {
179 Double_t pt=v[0];
180 Double_t lt=v[1];
181
182 Double_t ctau=p[1];
183 Double_t mass=p[2];
184 Double_t ct=mass*lt/pt;
185 if ((lt < 2) || (ct > 2.5*ctau)) {
186 TF1::RejectPoint();
187 return 0;
188 }
189
190 Int_t i=pt/0.1 + 1 + 1;
191 return p[i]*TMath::Exp(-ct/p[0]);
192}
193
194void Correct(TH1 *rw, const TH1 *as, const TH1 *mc) {
195 TH1 *eff = (TH1*)as->Clone();
196 eff->Divide(as,mc,1,1,"b");
197 rw->Divide(eff);
198 delete eff;
199}
200
8bd51e5f 201void Normalise(Double_t br, Double_t yw, Double_t bw, Double_t ne, TH1 *pt) {
202 pt->Scale(1/br); // branching ratio
203 pt->Scale(1/yw); // rapidity window
204 pt->Scale(1/bw); // bin width
205 pt->Scale(1/ne); // number of events
206}