PWGGA/PHOSTasks/PHOS_PbPb/macros/Flow/MethodsP.h

   1 char inname[256];
   2
   3 //  const Int_t nbin=11 ;
   4 //  Double_t xa[14]={1.,1.5,2.,2.5,3.,4.,5.,6.,8.,10.,14.,20.} ;
   5 //  const Int_t nbin=12 ;
   6 //  Double_t xa[13]={1.,1.5,2.,2.5,3.,4.,5.,7.,9.,11.,13.,16.,20.} ;
   7
   8 //const Int_t nbin=16 ;
   9 //Double_t xa[17]={0.6,1,1.5,2,2.5,3,3.5,4,5,6,7,8,10,12,14,16,20};
  10
  11 //from PCM
  12 const Int_t nbin=29 ;
  13 Double_t xa[30]={0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,3.0,3.5,4.0,5,6,7,8,10,12,14,16,20};
  14
  15 //  const Int_t nbin=12 ;
  16 //  Double_t xa[13]={0.6,1.,1.5,2.,2.5,3.,4.,5.,7.,10.,15.,20.} ;
  17
  18
  19 Double_t Ollitrault(Double_t chi){
  20   Double_t x = 0.25*chi*chi;
  21   Double_t resk1 = 0.626657 * chi * exp(-x) * (TMath::BesselI0((float)x) + TMath::BesselI1((float)x));
  22   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));
  23
  24   return resk1;
  25 }
  26
  27 Double_t GetCos(TH2F* hcos2AC, Int_t cen){
  28   TH1D *hres = 0 ;
  29   TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  30   Double_t resMean = 0, weight=0. ;
  31
  32 Int_t bin0, bin1;
  33
  34 if(cen==10){ bin0=1; bin1=3;}
  35 if(cen==11){ bin0=5; bin1=9;}
  36 if(cen==0){ bin0=1; bin1=2;}
  37 if(cen==1){ bin0=2; bin1=3;}
  38 if(cen==2){ bin0=3; bin1=5;}
  39 if(cen==3){ bin0=5; bin1=7;}
  40 if(cen==4){ bin0=7; bin1=9;}
  41 if(cen==5){ bin0=9; bin1=11;}
  42
  43   for(Int_t i=bin0; i<bin1; i++){
  44     TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ;
  45     Double_t wi = hwi->GetMean() ;
  46     hres = hcos2AC->ProjectionX("res",i,i) ;
  47     resMean += wi*hres->GetMean() ;
  48     weight += wi ;
  49     delete hwi ;
  50     delete hres ;
  51   }
  52   if(weight>0.){
  53     resMean/=weight ;
  54     return resMean ;
  55   }
  56   else{
  57     return 0. ;
  58   }
  59 }
  60 Double_t GetRes(TH2F* hcos2AC, Int_t cen){
  61
  62   TH1D *hres = 0 ;
  63   Double_t resMean = 0, weight=0. ;
  64
  65 Int_t bin0, bin1;
  66
  67 if(cen==10){ bin0=1; bin1=3;} //0-10 (Central)
  68 if(cen==11){ bin0=5; bin1=9;} //20-40 (SemiCentral)
  69 if(cen==0){ bin0=1; bin1=2;}
  70 if(cen==1){ bin0=2; bin1=3;}
  71 if(cen==2){ bin0=3; bin1=5;}
  72 if(cen==3){ bin0=5; bin1=7;}
  73 if(cen==4){ bin0=7; bin1=9;}
  74 if(cen==5){ bin0=9; bin1=11;}
  75
  76   TH1D *hres = 0 ;
  77   TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  78   Double_t resMean = 0, weight=0. ;
  79
  80   for(Int_t i=bin0; i<bin1; i++){
  81     TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ;
  82     Double_t wi = hwi->GetMean() ;
  83     hres = hcos2AC->ProjectionX("res",i,i) ;
  84     resMean += wi*hres->GetMean() ;
  85     weight += wi ;
  86     delete hwi ;
  87     delete hres ;
  88   }
  89   if(weight>0.){
  90     resMean/=weight ;
  91   }
  92   else{
  93     resMean = 0. ;
  94   }
  95
  96
  97   double chi = 2;
  98   double resSub = TMath::Sqrt(resMean);
  99   for(int j=0; j<15; j++){
 100     double resEve = Ollitrault(chi);
 101     chi= (resEve < resSub) ? chi+1.0*pow(0.5, j) : chi-1.0*pow(0.5, j);
 102   }
 103   return Ollitrault(TMath::Sqrt(2.)*chi);
 104 }
 105
 106 TH2F* GetRealMixed(Int_t cen, char w[255], char d[255], char kind[255], char kind2[255], char PID[255], Int_t Real=1){
 107   if(cen==10||cen==0||cen==1) sprintf(inname,"../flow11hCentral_Apr14.root");
 108   else sprintf(inname,"../flow11hSemiCentral_Apr14.root");
 109   TFile * f = new TFile(inname) ;
 110   char key[125] ;
 111   TH2F * h2DR =0 ;TH2F * h2DM=0;
 112
 113   if(cen==10){//centrality 0-10
 114     sprintf(key,"h%sPhotPhi%s%s%s_cen0",w,d,kind,PID) ;
 115 //cout<<key<<endl;
 116     h2DR = (TH2F*)f->Get(key) ;
 117
 118     sprintf(key,"h%sPhotPhi%s%s%s_cen1",w,d,kind,PID) ;
 119     TH2F *h2DRa = (TH2F*)f->Get(key) ;
 120
 121     sprintf(key,"h%sPhotPhi%s%s%s_cen2",w,d,kind,PID) ;
 122     h2DRb = (TH2F*)f->Get(key) ;
 123     sprintf(key,"h%sPhotPhi%s%s%s_cen3",w,d,kind,PID) ;
 124     h2DRc = (TH2F*)f->Get(key) ;
 125
 126     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 127     Double_t e1 = hCen->Integral(6,8)/3.;
 128     Double_t e2 = hCen->Integral(9,9);
 129     Double_t e3 = hCen->Integral(10,10);
 130
 131 cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;
 132
 133 //    h2DR->Add(h2DRa); delete h2DRa;
 134     h2DRa->Add(h2DRb,e1/e2); delete h2DRb;
 135     h2DRa->Add(h2DRc,e1/e3); delete h2DRc;
 136     h2DR->Add(h2DRa); delete h2DRa;
 137
 138
 139 //    h2DR->Add(h2DRa); delete h2DRa;
 140 //    h2DR->Add(h2DRb); delete h2DRb;
 141 //    h2DR->Add(h2DRc); delete h2DRc;
 142
 143   }
 144   if(cen==11){//centrality 20-40
 145     sprintf(key,"h%sPhotPhi%s%s%s_cen5",w,d,kind,PID) ;
 146     h2DR = (TH2F*)f->Get(key) ;
 147
 148     sprintf(key,"h%sPhotPhi%s%s%s_cen6",w,d,kind,PID) ;
 149     TH2F *h2DRa = (TH2F*)f->Get(key) ;
 150
 151     h2DR->Add(h2DRa); delete h2DRa;
 152   }
 153   if(cen==0){//centrality 0-5
 154     sprintf(key,"h%sPhotPhi%s%s%s_cen0",w,d,kind,PID) ;
 155     h2DR = (TH2F*)f->Get(key) ;
 156   }
 157   if(cen==1){//centrality 5-10
 158     sprintf(key,"h%sPhotPhi%s%s%s_cen1",w,d,kind,PID) ;
 159     h2DR = (TH2F*)f->Get(key) ;
 160     sprintf(key,"h%sPhotPhi%s%s%s_cen2",w,d,kind,PID) ;
 161     TH2F* h2DRa = (TH2F*)f->Get(key) ;
 162     sprintf(key,"h%sPhotPhi%s%s%s_cen3",w,d,kind,PID) ;
 163     TH2F* h2DRb = (TH2F*)f->Get(key) ;
 164
 165     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 166     Double_t e1 = hCen->Integral(6,8)/3.;
 167     Double_t e2 = hCen->Integral(9,9);
 168     Double_t e3 = hCen->Integral(10,10);
 169
 170 cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;
 171
 172     h2DR->Add(h2DRa,e1/e2); delete h2DRa;
 173     h2DR->Add(h2DRb,e1/e3); delete h2DRb;
 174
 175 //    h2DR->Add(h2DRa); delete h2DRa;
 176 //    h2DR->Add(h2DRb); delete h2DRb;
 177   }
 178   if(cen==2){//centrality 10-20
 179     sprintf(key,"h%sPhotPhi%s%s%s_cen4",w,d,kind,PID) ;
 180     h2DR = (TH2F*)f->Get(key) ;
 181
 182     sprintf(key,"h%sPhotPhi%s%s%s_cen0",w,d,kind,PID) ;
 183     TH2F* h2DRa = (TH2F*)f->Get(key) ;
 184     sprintf(key,"h%sPhotPhi%s%s%s_cen1",w,d,kind,PID) ;
 185     TH2F* h2DRb = (TH2F*)f->Get(key) ;
 186     sprintf(key,"h%sPhotPhi%s%s%s_cen2",w,d,kind,PID) ;
 187     TH2F* h2DRc = (TH2F*)f->Get(key) ;
 188     sprintf(key,"h%sPhotPhi%s%s%s_cen3",w,d,kind,PID) ;
 189     TH2F* h2DRd = (TH2F*)f->Get(key) ;
 190
 191
 192     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 193     Double_t e1 = hCen->Integral(11,11);
 194     Double_t e2 = hCen->Integral(12,12);
 195     Double_t e3 = hCen->Integral(13,13);
 196     Double_t e4 = hCen->Integral(14,15)/2.;
 197     Double_t e5 = hCen->Integral(16,20)/5.;
 198
 199 cout<<"w_11="<<e5/e1<<", w_12="<<e5/e2<<", w_13="<<e5/e3<<", w_14-15="<<e5/e4<<endl;
 200
 201     h2DR->Add(h2DRa,e5/e1); delete h2DRa;
 202     h2DR->Add(h2DRb,e5/e2); delete h2DRb;
 203     h2DR->Add(h2DRc,e5/e3); delete h2DRc;
 204     h2DR->Add(h2DRd,e5/e4); delete h2DRd;
 205
 206   }
 207   if(cen==3){//centrality 20-30
 208     sprintf(key,"h%sPhotPhi%s%s%s_cen5",w,d,kind,PID) ;
 209     h2DR = (TH2F*)f->Get(key) ;
 210   }
 211   if(cen==4){//centrality 30-40
 212     sprintf(key,"h%sPhotPhi%s%s%s_cen6",w,d,kind,PID) ;
 213     h2DR = (TH2F*)f->Get(key) ;
 214   }
 215   if(cen==5){//centrality 40-50
 216     sprintf(key,"h%sPhotPhi%s%s%s_cen7",w,d,kind,PID) ;
 217     h2DR = (TH2F*)f->Get(key) ;
 218   }
 219
 220   if(Real)return h2DR;
 221   else return h2DM;
 222 }
 223
 224 TH1F * GetCen(){
 225 cout<<inname<<endl;
 226   TFile * f = new TFile(inname) ;
 227   return (TH1F*)f->Get("hCentrality") ;
 228 }
 229
 230 TH2F * GetCenPHOS(){
 231   TFile * f = new TFile(inname) ;
 232   return (TH2F*)f->Get("hCenPHOS") ;
 233 }
 234
 235