PWGLF/STRANGENESS/LambdaK0PbPb/SpectraV0CutVariations.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 TFile *fAss;
  21 TFile *fGen;
  22 TFile *fRaw;
  23
  24 void Generated();
  25 void Corrections(Float_t cMin, Float_t cMax, TString centr);
  26 void RawYields(Float_t cMin, Float_t cMax, TString centr);
  27 void Spectra(TString centr);
  28
  29 void SpectraV0CutVariations(Float_t cMin, Float_t cMax, TString centr) {
  30
  31   fRaw=TFile::Open((centr+"/AliV0CutVariations.root").Data());
  32   fAss=TFile::Open((centr+"/AliV0CutVariationsMC.root").Data());
  33   fGen=TFile::Open("Generated.root");
  34   if (!fGen) {
  35      Generated();
  36      fGen=TFile::Open("Generated.root");
  37   }
  38
  39   Corrections(cMin,cMax,centr);
  40   RawYields(cMin,cMax,centr);
  41   Spectra(centr);
  42
  43   fAss->Close();
  44   fGen->Close();
  45   fRaw->Close();
  46 }
  47
  48 TH1 *GetEfficiency(Float_t, Float_t, const Char_t *, const Char_t *);
  49 TH1 *GetFeedDown(Float_t, Float_t, TString);
  50
  51 void Corrections(Float_t cmin, Float_t cmax, TString centr) {
  52
  53   TString name;
  54
  55   TH1 *effK0s=
  56   GetEfficiency(cmin, cmax, "fK0sAs","f3dHistPrimRawPtVsYVsMultK0Short");
  57   name="eff_K0s_";
  58   name+=centr;
  59   effK0s->SetName(name.Data());
  60
  61   TH1 *effLambda=
  62   GetEfficiency(cmin, cmax, "fLambdaAs","f3dHistPrimRawPtVsYVsMultLambda");
  63   name="eff_Lambda_";
  64   name+=centr;
  65   effLambda->SetName(name.Data());
  66     TH1 *fdLambda=GetFeedDown(cmin, cmax, centr);
  67     name="fd_Lambda_";
  68     name+=centr;
  69     fdLambda->SetName(name.Data());
  70
  71
  72   TH1 *effLambdaBar=
  73   GetEfficiency(cmin,cmax,"fLambdaBarAs","f3dHistPrimRawPtVsYVsMultAntiLambda");
  74   name="eff_LambdaBar_";
  75   name+=centr;
  76   effLambdaBar->SetName(name.Data());
  77
  78   TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
  79     effK0s->Write();
  80     effLambda->Write();    fdLambda->Write();
  81     effLambdaBar->Write();
  82   f->Close();
  83 }
  84
  85 TH1 *GetEfficiency(Float_t cMin, Float_t cMax,
  86                    const Char_t *chis, const Char_t *znam) {
  87   Double_t xbins[]={
  88    0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,
  89    1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
  90    2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
  91    4.5,5.0,5.5,6.5,8.0,10.0,12.0
  92   };
  93   Int_t nb=sizeof(xbins)/sizeof(Double_t);
  94   Int_t nb1=nb-1;
  95
  96   // Numerator
  97   fAss->cd();
  98   TH2F *f2d=(TH2F*)gDirectory->Get(chis); f2d->Sumw2();
  99   TH1D *hAs=f2d->ProjectionX("hAs",0,-1,"e");
 100
 101   // Denominator
 102   fGen->cd();
 103   TH3F *f3d = (TH3F*)gDirectory->Get(znam);
 104   f3d->Sumw2();
 105
 106   TH1D *fpt = f3d->ProjectionX("fpt",
 107                   f3d->GetYaxis()->FindBin(-0.5+1e-2),
 108                   f3d->GetYaxis()->FindBin(+0.5-1e-2),
 109                   f3d->GetZaxis()->FindBin(cMin+1e-2),
 110                   f3d->GetZaxis()->FindBin(cMax-1e-2)
 111               );
 112    TH1 *hMc=fpt->Rebin(nb1,"hMc",xbins);
 113
 114   //Efficiency
 115   TH1 *eff = (TH1*)hAs->Clone();
 116   eff->Divide(hAs,hMc,1,1,"b");
 117
 118   return eff;
 119 }
 120
 121 Bool_t GetBinContentError(const TH1 *h, Double_t x, Double_t &c, Double_t &e) {
 122    Int_t i1=h->GetXaxis()->FindFixBin(x);
 123    if (i1 <=0) return kFALSE;
 124    Int_t nb=h->GetNbinsX();
 125    if (i1 >nb) return kFALSE;
 126
 127    Double_t x1=h->GetBinCenter(i1);
 128
 129    if (TMath::Abs(x1-x) < 1e-13) {
 130      c=h->GetBinContent(i1);
 131      e=h->GetBinError(i1);
 132      return kTRUE;
 133    }
 134
 135    Int_t i2 = (x1 < x) ? i1+1 : i1-1;
 136    if (i2 <=0) i2=i1+1;
 137    if (i2 >nb) i2=i1-1;
 138    Double_t x2=h->GetBinCenter(i2);
 139
 140    Double_t c1=h->GetBinContent(i1);
 141    Double_t c2=h->GetBinContent(i2);
 142    c = c1 + (x - x1)*(c2 - c1)/(x2 - x1);
 143    e = 0.5*(h->GetBinError(i1) + h->GetBinError(i2));
 144
 145    return kTRUE;
 146 }
 147
 148 TH1 *GetFeedDown(Float_t cMin, Float_t cMax, TString cent) {
 149   // Get the FD matrix
 150   fAss->cd();
 151   TH3F *f3d = (TH3F*)gDirectory->Get("fLambdaFromXi");
 152   f3d->Sumw2();
 153   TH2D *fdMat = (TH2D*)f3d->Project3D("zxe");
 154   const Int_t nbx=fdMat->GetNbinsX();
 155   const Int_t nby=fdMat->GetNbinsY();
 156   TAxis *xiPtAxis=fdMat->GetYaxis();
 157
 158   // Re-normalise the FD matrix with MC Xi spectrum
 159   fGen->cd();
 160   f3d = (TH3F*)gDirectory->Get("f3dHistGenPtVsYVsMultXiMinus");
 161   f3d->Sumw2();
 162
 163   TH1D *hMcXi = f3d->ProjectionX("fpt",
 164                   f3d->GetYaxis()->FindBin(-0.5+1e-2),
 165                   f3d->GetYaxis()->FindBin(+0.5-1e-2),
 166                   f3d->GetZaxis()->FindBin(cMin+1e-2),
 167                   f3d->GetZaxis()->FindBin(cMax-1e-2)
 168               );
 169
 170   TH2D *h2McXi=new TH2D(*fdMat);
 171   h2McXi->Reset();
 172   h2McXi->Sumw2();
 173
 174   for (Int_t i=1; i<=nbx; i++) {
 175       for (Int_t j=1; j<=nby; j++) {
 176           Double_t pt=xiPtAxis->GetBinCenter(j);
 177           Double_t c=0.,e=0.;
 178           if (!GetBinContentError(hMcXi,pt,c,e)) continue;
 179           h2McXi->SetBinContent(i,j,c);
 180           h2McXi->SetBinError(i,j,e);
 181       }
 182   }
 183   fdMat->Divide(h2McXi);
 184   delete h2McXi;
 185
 186
 187   // Multiply the re-normalised matrix with the real Xi spectrum
 188   TString fileName("systcorrectedpt_cent");
 189      //*** "Dictionary" ***
 190      if (cent=="0005") cent="0010";
 191      if (cent=="6080") cent="6090";
 192      if (cent=="8090") cent="6090";
 193   TFile *f=TFile::Open((fileName+cent+".root").Data());
 194      TH1F *hReXi=(TH1F*)gDirectory->Get("correctedpt_0");
 195
 196      TH2D *h2ReXi=new TH2D(*fdMat);
 197      h2ReXi->Reset();
 198      h2ReXi->Sumw2();
 199
 200      for (Int_t i=1; i<=nbx; i++) {
 201          for (Int_t j=1; j<=nby; j++) {
 202             Double_t pt=xiPtAxis->GetBinCenter(j);
 203             Double_t c=0.,e=0.;
 204             if (!GetBinContentError(hReXi,pt,c,e)) continue;
 205             h2ReXi->SetBinContent(i,j,c);
 206             h2ReXi->SetBinError(i,j,e);
 207          }
 208      }
 209      fdMat->Multiply(h2ReXi);
 210      delete h2ReXi;
 211      f->Close();
 212
 213   TH1 *fd=fdMat->ProjectionX("_px",0,-1,"e");
 214   fd->Scale(0.1,"width");
 215   return fd;
 216 }
 217
 218
 219 void RawYields(Float_t cMin, Float_t cMax, TString centr) {
 220   TString name;
 221   TH2F *f2d=0;
 222
 223   //+++ Number of events for normalisation
 224   TFile *file=TFile::Open("LHC10h_pass2/Merged.root");
 225     TList *v0list=(TList *)gFile->Get("PWGLFExtractV0_PP/clistV0");
 226     TH1F *fMult=(TH1F*)v0list->FindObject("fHistMultiplicity");
 227     Int_t i1=fMult->GetXaxis()->FindBin(cMin+1e-2);
 228     Int_t i2=fMult->GetXaxis()->FindBin(cMax-1e-2);
 229     Float_t nEvents=fMult->Integral(i1,i2);
 230   file->Close();
 231
 232   fRaw->cd();
 233
 234   name="raw_K0s_";
 235   name+=centr;
 236   f2d=(TH2F*)gDirectory->Get("fK0sSi"); f2d->Sumw2();
 237   TH1D *rawK0s=f2d->ProjectionX(name,0,-1,"e");
 238   rawK0s->Scale(1/nEvents,"width");
 239
 240   name="raw_Lambda_";
 241   name+=centr;
 242   f2d=(TH2F*)gDirectory->Get("fLambdaSi"); f2d->Sumw2();
 243   TH1D *rawLambda=f2d->ProjectionX(name,0,-1,"e");
 244   rawLambda->Scale(1/nEvents,"width");
 245
 246   name="raw_LambdaBar_";
 247   name+=centr;
 248   f2d=(TH2F*)gDirectory->Get("fLambdaBarSi"); f2d->Sumw2();
 249   TH1D *rawLambdaBar=f2d->ProjectionX(name,0,-1,"e");
 250   rawLambdaBar->Scale(1/nEvents,"width");
 251
 252   TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
 253     rawK0s->Write();
 254     rawLambda->Write();
 255     rawLambdaBar->Write();
 256   f->Close();
 257 }
 258
 259 Double_t fd(Double_t x) {
 260   //Effective FD correction
 261   return 0.1619 + 0.05295*x - 0.01749*x*x + 0.001425*x*x*x - 3.446e-05*x*x*x*x;
 262 }
 263
 264 void Spectra(TString centr) {
 265   TString name;
 266   TH1 *eff=0;
 267   TH1D *raw=0;
 268   TH1D *spe=0;
 269
 270   TFile *f=TFile::Open("SpectraV0CutVariations.root","update");
 271     name="eff_K0s_";
 272     name+=centr;
 273     eff = (TH1*)gDirectory->Get(name.Data());
 274     name="raw_K0s_";
 275     name+=centr;
 276     raw = (TH1D*)gDirectory->Get(name.Data());
 277     spe = new TH1D(*raw);
 278     spe->Divide(eff);
 279     name="K0s_";
 280     name+=centr;
 281     spe->SetName(name.Data());
 282     spe->Write();
 283
 284     name="eff_Lambda_";
 285     name+=centr;
 286     eff = (TH1*)gDirectory->Get(name.Data());
 287     name="raw_Lambda_";
 288     name+=centr;
 289     raw = (TH1D*)gDirectory->Get(name.Data());
 290     spe = new TH1D(*raw);
 291     for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
 292         Double_t pt=spe->GetBinCenter(i);
 293         Double_t c=spe->GetBinContent(i);
 294         c -= (c*fd(pt));
 295         spe->SetBinContent(i,c);
 296     }
 297     spe->Divide(eff);
 298     name="Lambda_";
 299     name+=centr;
 300     spe->SetName(name.Data());
 301     spe->Write();
 302
 303     name="eff_LambdaBar_";
 304     name+=centr;
 305     eff = (TH1*)gDirectory->Get(name.Data());
 306     name="raw_LambdaBar_";
 307     name+=centr;
 308     raw = (TH1D*)gDirectory->Get(name.Data());
 309     spe = new TH1D(*raw);
 310     for (Int_t i=1; i<=spe->GetNbinsX(); i++) {
 311         Double_t pt=spe->GetBinCenter(i);
 312         Double_t c=spe->GetBinContent(i);
 313         c -= (c*fd(pt));
 314         spe->SetBinContent(i,c);
 315     }
 316     spe->Divide(eff);
 317     name="LambdaBar_";
 318     name+=centr;
 319     spe->SetName(name.Data());
 320     spe->Write();
 321   f->Close();
 322 }
 323
 324 void Generated() {
 325   TList *v0listMC=0;
 326   TH3F *h3=0;
 327   //
 328   TFile::Open("LHC11a10b_plus/Merged.root");
 329   v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
 330
 331   TH3F *k0s =
 332     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
 333   TH3F *k0s_nonInj =
 334     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
 335
 336   TH3F *lam =
 337     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
 338   TH3F *lam_nonInj =
 339     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
 340
 341   TH3F *alam =
 342     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
 343   TH3F *alam_nonInj =
 344     (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
 345
 346   TH3F *xiMinus =
 347     (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
 348   TH3F *xiPlus =
 349     (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
 350
 351
 352   //
 353   TFile::Open("LHC11a10b_bis/Merged.root");
 354   v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
 355
 356   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
 357   k0s->Add(h3);
 358   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjK0Short");
 359   k0s_nonInj->Add(h3);
 360
 361   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
 362   lam->Add(h3);
 363   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjLambda");
 364   lam_nonInj->Add(h3);
 365
 366   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
 367   alam->Add(h3);
 368   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultNonInjAntiLambda");
 369   alam_nonInj->Add(h3);
 370
 371   h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
 372   xiMinus->Add(h3);
 373   h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
 374   xiPlus->Add(h3);
 375
 376
 377   //
 378   TFile::Open("LHC11a10a_bis/Merged.root");
 379   v0listMC=(TList *)gFile->Get("PWGLFExtractPerformanceV0_PP_MC/clistV0MC");
 380
 381   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
 382   k0s_nonInj->Add(h3);
 383   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
 384   lam_nonInj->Add(h3);
 385   h3 = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
 386   alam_nonInj->Add(h3);
 387   h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiMinus");
 388   xiMinus->Add(h3);
 389   h3 = (TH3F*)v0listMC->FindObject("f3dHistGenPtVsYVsMultXiPlus");
 390   xiPlus->Add(h3);
 391
 392   TFile *f=TFile::Open("Generated.root","new");
 393   k0s->Write(); k0s_nonInj->Write();
 394   lam->Write(); lam_nonInj->Write();
 395   alam->Write(); alam_nonInj->Write();
 396   xiMinus->Write();
 397   xiPlus->Write();
 398   f->Close();
 399 }
 400