Fix: Correcting the feeddown for the efficiency
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / PostProcessCTau.C
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
33   TString name("cTau_0090"); //centrality
34   //TString name("cTau_0005"); //centrality
35   //TString name("cTau_2040"); //centrality
36   //TString name("cTau_4060"); //centrality
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();
48   TH2F *fLambdaBarSi=
49         (TH2F*)list->FindObject("fLambdaBarSi"); fLambdaBarSi->Sumw2();
50   TH1F *fXiBarSiP=(TH1F*)list->FindObject("fXiBarSiP"); fXiBarSiP->Sumw2();
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
67   TH2F *fLambdaBarMC=
68         (TH2F*)listmc->FindObject("fLambdaBarMC"); fLambdaBarMC->Sumw2();
69   TH2F *fLambdaBarAs=
70         (TH2F*)listmc->FindObject("fLambdaBarAs"); fLambdaBarAs->Sumw2();
71
72   TH1F *fXiSiPMC=(TH1F*)listmc->FindObject("fXiSiP");fXiSiPMC->Sumw2();
73   TH3F *fLambdaFromXi=(TH3F*)listmc->FindObject("fLambdaFromXi");
74   fLambdaFromXi->Sumw2();
75
76   TH1F *fXiBarSiPMC=(TH1F*)listmc->FindObject("fXiBarSiP");fXiBarSiPMC->Sumw2();
77   TH3F *fLambdaBarFromXiBar=(TH3F*)listmc->FindObject("fLambdaBarFromXiBar");
78   fLambdaBarFromXiBar->Sumw2();
79   //+++++
80
81
82   Double_t myExp2(Double_t *v, Double_t *p);
83   void Correct(TH1 *rw, const TH1 *as, const TH1 *mc);
84   void Normalise(Double_t, Double_t, Double_t, Double_t, TH1 *h);
85
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};
91
92   const TH2 *in[]={
93     fK0sSi,fK0sAs,fK0sMC,
94     fLambdaSi,fLambdaAs,fLambdaMC,
95     fLambdaBarSi,fLambdaBarAs,fLambdaBarMC
96   };
97   TH1 *fd[]={
98     0,0,0,
99     fLambdaFromXi,fXiSiP,fXiSiPMC,
100     fLambdaBarFromXiBar,fXiBarSiP,fXiBarSiPMC
101   };
102   Double_t wbx=fK0sSi->GetBinWidth(3);
103   Int_t nbx=fLambdaFromXi->GetNbinsX();
104   Int_t nby=fLambdaFromXi->GetNbinsY();
105   Int_t nbz=fLambdaFromXi->GetNbinsZ();
106  
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());
115       Normalise(brch[i], 2*0.5, wbx, nEvents, pt);      
116       Double_t eipt=0., ipt=pt->IntegralAndError(1, nbx, eipt);
117
118       new TCanvas; 
119       pt->Draw(); 
120
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++) {
131                      Float_t c=fd3->GetBinContent(k,l,m);
132                      c *= rl->GetBinContent(m);
133                      fd3->SetBinContent(k,l,m,c);
134                  }
135              }
136          }
137
138          TH2 *fd2=(TH2*)fd3->Project3D("yxe");
139          Correct(fd2, in[i*3 + 1], in[i*3 + 2]);
140
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       } 
154       //continue;
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);  
171       cout<<endl<<pnam[i]<<"  \tstatus: "<<status<<"   chi2/ndf: "<<chi2<<
172         "\tc*tau: "<<slope<<"+/-"<<error<<endl<<endl;
173   }
174
175   return;
176
177
178 Double_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
194 void 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
201 void 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 }