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(Int_t corrType=1) {
  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, 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   const Double_t yWin=2*0.5; // rapidity window
  92   const Double_t wbx=fK0sSi->GetBinWidth(3);
  93   const Int_t nbx=fLambdaFromXi->GetNbinsX();
  94   const Int_t nby=fLambdaFromXi->GetNbinsY();
  95   const Int_t nbz=fLambdaFromXi->GetNbinsZ();
  96
  97   fK0sMC->Scale(1/brch[0]);
  98   fLambdaMC->Scale(1/brch[1]);
  99   fLambdaBarMC->Scale(1/brch[2]);
 100
 101   const TH2 *in[]={
 102     fK0sSi,fK0sAs,fK0sMC,
 103     fLambdaSi,fLambdaAs,fLambdaMC,
 104     fLambdaBarSi,fLambdaBarAs,fLambdaBarMC
 105   };
 106   TH1 *fd[]={
 107     0,0,0,
 108     fLambdaFromXi,fXiSiP,fXiSiPMC,
 109     fLambdaBarFromXiBar,fXiBarSiP,fXiBarSiPMC
 110   };
 111
 112   for (Int_t i=0; i<iMax; i++) {
 113       const TH2 *ptRw=in[3*i+0], *ptAs=in[3*i+1], *ptMc=in[3*i+2];
 114       TH1 *ptAsPx=0, *ptMcPx=0, *pt=0;
 115
 116       TH2 *cr=(TH2*)ptRw->Clone();
 117       Correct(cr, ptAs, ptMc);
 118
 119       //++++ pT spectrum
 120       if (corrType==1) {
 121          pt    =ptRw->ProjectionX("_px",0,-1,"e");
 122          ptAsPx=ptAs->ProjectionX("_px",0,-1,"e");
 123          ptMcPx=ptMc->ProjectionX("_px",0,-1,"e");
 124          Correct(pt, ptAsPx, ptMcPx);
 125
 126          TString effName = "Efficiency for " + pnam[i] + " (" + name + ")";
 127          TH1 *eff = (TH1*)ptAsPx->Clone();
 128          eff->SetTitle(effName.Data());
 129          eff->Divide(ptAsPx,ptMcPx,1,1,"b");
 130          eff->SetName("eff");
 131          new TCanvas; eff->Draw();
 132       } else {
 133          pt=cr->ProjectionX("_px",0,-1,"e");
 134       }
 135       TString ptName = pnam[i] + " p_{T} spectrum";
 136       pt->SetTitle(ptName.Data());
 137       Normalise(yWin, wbx, nEvents, pt);
 138
 139       new TCanvas;
 140       pt->Draw();
 141
 142       if (i>0) {
 143       //++++ feeddown
 144          TH3 *fd3=(TH3*)fd[3*i + 0];
 145          TH1 *rl =(TH1*)fd[3*i + 1];
 146          TH1 *mc =(TH1*)fd[3*i + 2];
 147
 148          Double_t nLambdaFromXiMc=fd3->Integral();
 149          Double_t nXiMc=mc->Integral();
 150          Double_t nXiRl=rl->Integral();
 151          Double_t f=(nLambdaFromXiMc*nXiRl/nXiMc)/ptRw->Integral();
 152          Double_t ef=f*TMath::Sqrt(nXiMc)/nXiMc;
 153          cout<<endl<<"Global FD correction: "<<f<<"+/-"<<ef<<endl;
 154
 155          rl->Divide(mc);
 156
 157          for (Int_t k=1; k<=nbx; k++) {
 158              for (Int_t l=1; l<=nby; l++) {
 159                  for (Int_t m=1; m<=nbz; m++) {
 160                      Float_t c=fd3->GetBinContent(k,l,m);
 161                      c *= rl->GetBinContent(m);
 162                      fd3->SetBinContent(k,l,m,c);
 163                  }
 164              }
 165          }
 166
 167          TH2 *fd2=(TH2*)fd3->Project3D("yxe");
 168          Correct(fd2, ptAs, ptMc);
 169
 170          TH1 *fd1=0;
 171          if (corrType==1) {
 172             fd1=(TH1*)fd3->Project3D("x");
 173             Correct(fd1, ptAsPx, ptMcPx);
 174          } else {
 175             fd1=fd2->ProjectionX("_px",0,-1,"e");
 176          }
 177          Normalise(yWin, wbx, nEvents, fd1);
 178
 179          //new TCanvas();
 180          fd1->Draw("same");
 181
 182          cr->Add(fd2,-1);
 183       }
 184       continue;
 185
 186       //++++ c*tau
 187       TF2 *f2=new TF2("myexpo2",myExp2,0.,10.,0.,100.,1+1+1+nbx);
 188       f2->SetParameter(0, ctau[i]);
 189       f2->FixParameter(1, ctau[i]);
 190       f2->FixParameter(2, mass[i]);
 191       for (Int_t p=1+1+1; p<=nbx+1+1; p++)
 192           f2->SetParameter(p,fK0sSi->GetBinContent(p,1));
 193
 194       new TCanvas;
 195       TFitResultPtr r=cr->Fit(f2,"SQ");
 196
 197       Int_t status   = r;
 198       Double_t chi2  = r->Chi2()/r->Ndf();
 199       Double_t slope = r->Parameter(0);
 200       Double_t error = r->ParError(0);
 201       cout<<endl<<pnam[i]<<"  \tstatus: "<<status<<"   chi2/ndf: "<<chi2<<
 202         "\tc*tau: "<<slope<<"+/-"<<error<<endl<<endl;
 203   }
 204
 205   return;
 206 }
 207
 208 Double_t myExp2(Double_t *v, Double_t *p) {
 209    Double_t pt=v[0];
 210    Double_t lt=v[1];
 211
 212    Double_t ctau=p[1];
 213    Double_t mass=p[2];
 214    Double_t ct=mass*lt/pt;
 215    if ((lt < 2) || (ct > 2.5*ctau)) {
 216       TF1::RejectPoint();
 217       return 0;
 218    }
 219
 220    Int_t i=pt/0.1 + 1 + 1;
 221    return p[i]*TMath::Exp(-ct/p[0]);
 222 }
 223
 224 void Correct(TH1 *rw, const TH1 *as, const TH1 *mc) {
 225   TH1 *eff = (TH1*)as->Clone();
 226   eff->Divide(as,mc,1,1,"b");
 227   rw->Divide(eff);
 228   delete eff;
 229 }
 230
 231 void Normalise(Double_t yw, Double_t bw, Double_t ne, TH1 *pt) {
 232    pt->Scale(1/yw);    // rapidity window
 233    pt->Scale(1/bw);    // bin width
 234    pt->Scale(1/ne);    // number of events
 235 }