PWGGA/PHOSTasks/PHOS_PbPb/macros/Flow/RP.C

   1 void RP(){
   2
   3   gStyle->SetOptFit(0);
   4   gStyle->SetOptStat(0);
   5
   6   TFile *f = new TFile("../flow11h.root");
   7
   8   if(!f->IsOpen()){
   9     printf("No such file\n");
  10     return 1;
  11   }
  12
  13   char key[55] ;
  14 //Get histograms
  15
  16   //Read resolution
  17   TH2F * hcos2AC = (TH2F*)f->Get("cos2AC");
  18   TH2F * hcos2V0AC = (TH2F*)f->Get("cos2V0AC");
  19   TH2F * hcos2V0ATPC = (TH2F*)f->Get("cos2V0ATPC");
  20   TH2F * hcos2V0CTPC = (TH2F*)f->Get("cos2V0CTPC");
  21
  22 //  TH1F * hresCh = (TH1F*)f2->Get("hEvPlResEta4");
  23
  24   TH1D * hres = 0 ;
  25   TH1D * hresF = 0 ;
  26
  27   TH1D* rpT2  = new TH1D("rpTPC", "rpTPC 2 sub", hcos2AC->GetNbinsY(),0,100.);
  28   TH1D* rpT3 = new TH1D("rpTPC2", "rpTPC 3 sub", hcos2AC->GetNbinsY(),0,100.);
  29   TH1D* rpA  = new TH1D("rpV0A", "rpV0A", hcos2AC->GetNbinsY(),0,100.);
  30   TH1D* rpC  = new TH1D("rpV0C", "rpV0C", hcos2AC->GetNbinsY(),0,100.);
  31
  32   Double_t resT2, resTA, resTC, resAC;
  33
  34   for(int i=0;i<hcos2AC->GetNbinsY();i++){
  35     hres = hcos2AC->ProjectionX("res",i,i) ;
  36     resT2 = GetRes(hres);
  37     if(resT2>0)rpT2->SetBinContent(i, resT2);
  38
  39     hresF = hcos2V0AC->ProjectionX("res2",i,i) ;
  40     resAC = GetCos(hresF);
  41     hres = hcos2V0ATPC->ProjectionX("res3",i,i) ;
  42     resTA = GetCos(hres);
  43     hres = hcos2V0CTPC->ProjectionX("res4",i,i) ;
  44     resTC = GetCos(hres);
  45
  46     if(resAC>0)rpT3->SetBinContent(i, TMath::Sqrt(resTA*resTC/resAC));
  47     if(resTC>0)rpA->SetBinContent(i, TMath::Sqrt(resTA*resAC/resTC));
  48     if(resTA>0)rpC->SetBinContent(i, TMath::Sqrt(resTC*resAC/resTA));
  49
  50 //cout<<"i="<<i<<", rpA="<<TMath::Sqrt(resTA*resAC/resTC)<< endl;
  51   }
  52
  53   rpT2->SetTitle("Reaction plane resolution with different ALICE detectors");
  54   rpT2->GetXaxis()->SetTitle("centrality percentage");
  55   rpT2->GetYaxis()->SetTitle("resolution");
  56
  57   rpT2->SetMarkerStyle(20);
  58   rpT2->SetMarkerColor(2);
  59   rpT2->SetLineColor(2);
  60   rpT2->Draw("p");
  61
  62   rpT3->SetMarkerStyle(21);
  63   rpT3->SetMarkerColor(2);
  64   rpT3->SetLineColor(2);
  65   rpT3->Draw("psame");
  66
  67   rpA->SetMarkerStyle(22);
  68   rpA->SetMarkerColor(kMagenta);
  69   rpA->SetLineColor(kMagenta);
  70   rpA->Draw("psame");
  71
  72   rpC->SetMarkerStyle(22);
  73   rpC->SetMarkerColor(kGreen);
  74   rpC->SetLineColor(kGreen);
  75   rpC->Draw("psame");
  76
  77 /*
  78   hresCh->SetMarkerStyle(21);
  79   hresCh->SetMarkerColor(kBlue);
  80   hresCh->SetLineColor(kBlue);
  81   hresCh->Draw("same");
  82 */
  83
  84   TLegend * l = new TLegend(0.12,0.15,0.6,0.35) ;
  85   l->SetFillColor(0) ;
  86   l->SetTextSize(0.03) ;
  87   l->AddEntry(rpT2,"RPRes with TPC 2 subs","p") ;
  88   l->AddEntry(rpT3,"RPRes with TPC 3 subs","p") ;
  89   l->AddEntry(rpA,"RPRes with V0A","p") ;
  90   l->AddEntry(rpC,"RPRes with V0C","p") ;
  91 //  l->AddEntry(hresCh,"RPRes TPC (Alex Dobrin)","p") ;
  92   l->Draw() ;
  93
  94 }
  95
  96 Double_t Ollitrault(Double_t chi){
  97   Double_t x = 0.25*chi*chi;
  98   Double_t resk1 = 0.626657 * chi * exp(-x) * (TMath::BesselI0((float)x) + TMath::BesselI1((float)x));
  99   Double_t resk2 = 0.626657 * chi * exp(-x) * (TMath::Sqrt(2./TMath::Pi()/x)*TMath::SinH(x) + TMath::Sqrt(2./TMath::Pi()/x)*(TMath::CosH(x) - TMath::SinH(x)/x));
 100
 101   return resk1;
 102 }
 103
 104 Double_t GetCos(TH1D* hres){
 105   Double_t resMean = hres->GetMean() ;
 106   cout<<resMean<<endl;
 107   if(resMean>0)return resMean ;
 108   else return 0. ;
 109 }
 110
 111 Double_t GetRes(TH1D* hres){
 112   Double_t resMean = hres->GetMean() ;
 113   Double_t resOld=0 ;
 114   if(resMean>0)resOld=1./TMath::Sqrt(resMean) ;
 115
 116   float rxn2=resMean;
 117
 118   double chi = 2;
 119   for(int j=0; j<15; j++){
 120     double resSub = sqrt(rxn2);
 121     Double_t resEve = Ollitrault(chi);
 122     chi= (resEve < resSub) ? chi+1.0*pow(0.5, j) : chi-1.0*pow(0.5, j);
 123   }
 124   return Ollitrault(TMath::Sqrt(2)*chi);
 125 }
 126