PWGGA/PHOSTasks/PHOS_PbPb/macros/Flow/Methods.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]={0.6,1.,1.5,2.,3.,4.,5.,7.,9.,11.,13.,16.,20.} ;
   7 //  const Int_t nbin=11 ;
   8 //  Double_t xa[12]={0.6,1.,1.5,2.,2.5,3.,4.,5.,7.,10.,15.,20.} ;
   9
  10
  11 Double_t Ollitrault(Double_t chi){
  12   Double_t x = 0.25*chi*chi;
  13   Double_t resk1 = 0.626657 * chi * exp(-x) * (TMath::BesselI0((float)x) + TMath::BesselI1((float)x));
  14   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));
  15
  16   return resk1;
  17 }
  18
  19 Double_t GetCos(TH2F* hcos2AC, Int_t cen){
  20   TH1D *hres = 0 ;
  21   TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  22   Double_t resMean = 0, weight=0. ;
  23
  24   if(cen==21){
  25   TFile * f3 = new TFile("../flow11h_Apr14.root") ;
  26   hPHOS = (TH2F*)f3->Get("hCenPHOS");
  27   }
  28
  29 Int_t bin0, bin1;
  30
  31 if(cen==10){ bin0=1; bin1=3;}
  32 if(cen==11){ bin0=5; bin1=9;}
  33 if(cen==0){ bin0=1; bin1=2;}
  34 if(cen==1){ bin0=2; bin1=3;}
  35 if(cen==2){ bin0=3; bin1=5;}
  36 if(cen==3){ bin0=5; bin1=7;}
  37 if(cen==4){ bin0=7; bin1=9;}
  38 if(cen==5){ bin0=9; bin1=11;}
  39 if(cen==20){ bin0=3; bin1=11;} //10-50%
  40 if(cen==21){ bin0=1; bin1=5;} //0-20%
  41
  42   for(Int_t i=bin0; i<bin1; i++){
  43     TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ;
  44     Double_t wi = hwi->GetMean() ;
  45     hres = hcos2AC->ProjectionX("res",i,i) ;
  46     resMean += wi*hres->GetMean() ;
  47     weight += wi ;
  48     delete hwi ;
  49     delete hres ;
  50   }
  51   if(weight>0.){
  52     resMean/=weight ;
  53     return resMean ;
  54   }
  55   else{
  56     return 0. ;
  57   }
  58 }
  59 Double_t GetRes(TH2F* hcos2AC, Int_t cen){
  60
  61   TH1D *hres = 0 ;
  62   Double_t resMean = 0, weight=0. ;
  63
  64 Int_t bin0, bin1;
  65
  66 if(cen==10){ bin0=1; bin1=3;} //0-10 (Central)
  67 if(cen==11){ bin0=5; bin1=9;} //20-40 (SemiCentral)
  68 if(cen==0){ bin0=1; bin1=2;}
  69 if(cen==1){ bin0=2; bin1=3;}
  70 if(cen==2){ bin0=3; bin1=5;}
  71 if(cen==3){ bin0=5; bin1=7;}
  72 if(cen==4){ bin0=7; bin1=9;}
  73 if(cen==5){ bin0=9; bin1=11;}
  74 if(cen==20){ bin0=3; bin1=11;} //10-50%
  75 if(cen==21){ bin0=1; bin1=5;} //0-20%
  76
  77   TH1D *hres = 0 ;
  78   TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  79   Double_t resMean = 0, weight=0. ;
  80   if(cen==21){
  81   TFile * f3 = new TFile("../flow11h_Apr14.root") ;
  82   hPHOS = (TH2F*)f3->Get("hCenPHOS");
  83   }
  84
  85   for(Int_t i=bin0; i<bin1; i++){
  86     TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ;
  87     Double_t wi = hwi->GetMean() ;
  88     hres = hcos2AC->ProjectionX("res",i,i) ;
  89     resMean += wi*hres->GetMean() ;
  90     weight += wi ;
  91     delete hwi ;
  92     delete hres ;
  93   }
  94   if(weight>0.){
  95     resMean/=weight ;
  96   }
  97   else{
  98     resMean = 0. ;
  99   }
 100
 101
 102   double chi = 2;
 103   double resSub = TMath::Sqrt(resMean);
 104   for(int j=0; j<15; j++){
 105     double resEve = Ollitrault(chi);
 106     chi= (resEve < resSub) ? chi+1.0*pow(0.5, j) : chi-1.0*pow(0.5, j);
 107   }
 108   return Ollitrault(TMath::Sqrt(2.)*chi);
 109 }
 110
 111 TH3F* GetRealMixed(Int_t cen, char w[255], char d[255], char kind[255], char kind2[255], char PID[255], Int_t Real=1){
 112   if(cen==10||cen==0||cen==1||cen==21) sprintf(inname,"../flow11hCentral_Apr14.root");
 113   else sprintf(inname,"../flow11hSemiCentral_Apr14.root");
 114   char inname2[255];
 115   if(cen==21) sprintf(inname2,"../flow11hSemiCentral_Apr14.root");
 116
 117   TFile * f = new TFile(inname) ;
 118   if(cen==21) TFile * f2 = new TFile(inname2) ;
 119
 120   char key[125] ;
 121   TH3F * h3DR =0 ;TH3F * h3DM=0;
 122
 123   if(cen==10){//centrality 0-10
 124     sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
 125     h3DR = (TH3F*)f->Get(key) ;
 126     sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
 127     h3DM= (TH3F*)f->Get(key) ;
 128
 129     sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
 130     TH3F *h3DRa = (TH3F*)f->Get(key) ;
 131     sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
 132     TH3F *h3DMa= (TH3F*)f->Get(key) ;
 133
 134     sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
 135     h3DRb = (TH3F*)f->Get(key) ;
 136     sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
 137     h3DMb= (TH3F*)f->Get(key) ;
 138     sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
 139     h3DRc = (TH3F*)f->Get(key) ;
 140     sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
 141     h3DMc= (TH3F*)f->Get(key) ;
 142
 143     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 144     Double_t e1 = hCen->Integral(6,8)/3.;
 145     Double_t e2 = hCen->Integral(9,9);
 146     Double_t e3 = hCen->Integral(10,10);
 147
 148 cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;
 149
 150     h3DR->Add(h3DRa); delete h3DRa;
 151     h3DM->Add(h3DMa); delete h3DMa;
 152     h3DR->Add(h3DRb,e1/e2); delete h3DRb;
 153     h3DM->Add(h3DMb,e1/e2); delete h3DMb;
 154     h3DR->Add(h3DRc,e1/e3); delete h3DRc;
 155     h3DM->Add(h3DMc,e1/e3); delete h3DMc;
 156   }
 157   if(cen==21){//centrality 0-20
 158     sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
 159     h3DR = (TH3F*)f->Get(key) ;
 160     sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
 161     h3DM= (TH3F*)f->Get(key) ;
 162
 163     sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
 164     TH3F *h3DRa = (TH3F*)f->Get(key) ;
 165     sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
 166     TH3F *h3DMa= (TH3F*)f->Get(key) ;
 167
 168     sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
 169     h3DRb = (TH3F*)f->Get(key) ;
 170     sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
 171     h3DMb= (TH3F*)f->Get(key) ;
 172     sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
 173     h3DRc = (TH3F*)f->Get(key) ;
 174     sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
 175     h3DMc= (TH3F*)f->Get(key) ;
 176
 177     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 178     Double_t e1 = hCen->Integral(6,8)/3.;
 179     Double_t e2 = hCen->Integral(9,9);
 180     Double_t e3 = hCen->Integral(10,10);
 181
 182 cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;
 183
 184     h3DR->Add(h3DRa); delete h3DRa;
 185     h3DM->Add(h3DMa); delete h3DMa;
 186     h3DR->Add(h3DRb,e1/e2); delete h3DRb;
 187     h3DM->Add(h3DMb,e1/e2); delete h3DMb;
 188     h3DR->Add(h3DRc,e1/e3); delete h3DRc;
 189     h3DM->Add(h3DMc,e1/e3); delete h3DMc;
 190
 191 //now 10-20
 192     sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
 193     TH3F* h3DR1 = (TH3F*)f2->Get(key) ;
 194     sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
 195     TH3F* h3DM1= (TH3F*)f2->Get(key) ;
 196
 197     sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
 198     TH3F* h3DRa = (TH3F*)f2->Get(key) ;
 199     sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
 200     TH3F* h3DMa= (TH3F*)f2->Get(key) ;
 201     sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
 202     TH3F* h3DRb = (TH3F*)f2->Get(key) ;
 203     sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
 204     TH3F* h3DMb= (TH3F*)f2->Get(key) ;
 205     sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
 206     TH3F* h3DRc = (TH3F*)f2->Get(key) ;
 207     sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
 208     TH3F* h3DMc= (TH3F*)f2->Get(key) ;
 209     sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
 210     TH3F* h3DRd = (TH3F*)f2->Get(key) ;
 211     sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
 212     TH3F* h3DMd= (TH3F*)f2->Get(key) ;
 213
 214     TH1D* hCen2 = ((TH2F*)f2->Get("hCentrality"))->ProjectionX();
 215     Double_t e1 = hCen2->Integral(11,11);
 216     Double_t e2 = hCen2->Integral(12,12);
 217     Double_t e3 = hCen2->Integral(13,13);
 218     Double_t e4 = hCen2->Integral(14,15)/2.;
 219     Double_t e5 = hCen2->Integral(16,20)/5.;
 220
 221 cout<<"w_11="<<e5/e1<<", w_12="<<e5/e2<<", w_13="<<e5/e3<<", w_14-15="<<e5/e4<<endl;
 222
 223     h3DR->Add(h3DR1,e5/e1); delete h3DR1;
 224     h3DM->Add(h3DM1,e5/e1); delete h3DM1;
 225     h3DR->Add(h3DRa,e5/e1); delete h3DRa;
 226     h3DM->Add(h3DMa,e5/e1); delete h3DMa;
 227     h3DR->Add(h3DRb,e5/e2); delete h3DRb;
 228     h3DM->Add(h3DMb,e5/e2); delete h3DMb;
 229     h3DR->Add(h3DRc,e5/e3); delete h3DRc;
 230     h3DM->Add(h3DMc,e5/e3); delete h3DMc;
 231     h3DR->Add(h3DRd,e5/e4); delete h3DRd;
 232   }
 233   if(cen==11){//centrality 20-40
 234     sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
 235     h3DR = (TH3F*)f->Get(key) ;
 236     sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
 237     h3DM= (TH3F*)f->Get(key) ;
 238
 239     sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
 240     TH3F *h3DRa = (TH3F*)f->Get(key) ;
 241     sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
 242     TH3F *h3DMa= (TH3F*)f->Get(key) ;
 243
 244     h3DR->Add(h3DRa); delete h3DRa;
 245     h3DM->Add(h3DMa); delete h3DMa;
 246   }
 247   if(cen==20){//centrality 10-50
 248     sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
 249     h3DR = (TH3F*)f->Get(key) ;
 250     sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
 251     h3DM= (TH3F*)f->Get(key) ;
 252
 253     sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
 254     TH3F *h3DRa = (TH3F*)f->Get(key) ;
 255     sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
 256     TH3F *h3DMa= (TH3F*)f->Get(key) ;
 257
 258     sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
 259     TH3F *h3DRb = (TH3F*)f->Get(key) ;
 260     sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
 261     TH3F *h3DMb= (TH3F*)f->Get(key) ;
 262
 263     sprintf(key,"h%sMassPt%s%s%s_cen7",w,d,kind,PID) ;
 264     TH3F *h3DRc = (TH3F*)f->Get(key) ;
 265     sprintf(key,"hMiMassPt%s%s%s_cen7",d,kind2,PID) ;
 266     TH3F *h3DMc= (TH3F*)f->Get(key) ;
 267
 268     h3DR->Add(h3DRa); delete h3DRa;
 269     h3DM->Add(h3DMa); delete h3DMa;
 270     h3DR->Add(h3DRb); delete h3DRb;
 271     h3DM->Add(h3DMb); delete h3DMb;
 272     h3DR->Add(h3DRc); delete h3DRc;
 273     h3DM->Add(h3DMc); delete h3DMc;
 274   }
 275 //small bins
 276   if(cen==0){//centrality 0-5
 277     sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
 278     h3DR = (TH3F*)f->Get(key) ;
 279     sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
 280     h3DM= (TH3F*)f->Get(key) ;
 281   }
 282   if(cen==1){//centrality 5-10
 283     sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
 284     h3DR = (TH3F*)f->Get(key) ;
 285     sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
 286     h3DM= (TH3F*)f->Get(key) ;
 287     sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
 288     h3DRa = (TH3F*)f->Get(key) ;
 289     sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
 290     h3DMa= (TH3F*)f->Get(key) ;
 291     sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
 292     h3DRb = (TH3F*)f->Get(key) ;
 293     sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
 294     h3DMb= (TH3F*)f->Get(key) ;
 295
 296     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 297     Double_t e1 = hCen->Integral(6,8)/3.;
 298     Double_t e2 = hCen->Integral(9,9);
 299     Double_t e3 = hCen->Integral(10,10);
 300
 301 cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;
 302
 303     h3DR->Add(h3DRa,e1/e2); delete h3DRa;
 304     h3DM->Add(h3DMa,e1/e2); delete h3DMa;
 305     h3DR->Add(h3DRb,e1/e3); delete h3DRb;
 306     h3DM->Add(h3DMb,e1/e3); delete h3DMb;
 307
 308 /*
 309     h3DR->Add(h3DRa); delete h3DRa;
 310     h3DM->Add(h3DMa); delete h3DMa;
 311     h3DR->Add(h3DRb); delete h3DRb;
 312     h3DM->Add(h3DMb); delete h3DMb;
 313 */
 314   }
 315   if(cen==2){//centrality 10-20
 316     sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
 317     h3DR = (TH3F*)f->Get(key) ;
 318     sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
 319     h3DM= (TH3F*)f->Get(key) ;
 320
 321     sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
 322     TH3F* h3DRa = (TH3F*)f->Get(key) ;
 323     sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
 324     TH3F* h3DMa= (TH3F*)f->Get(key) ;
 325     sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
 326     TH3F* h3DRb = (TH3F*)f->Get(key) ;
 327     sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
 328     TH3F* h3DMb= (TH3F*)f->Get(key) ;
 329     sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
 330     TH3F* h3DRc = (TH3F*)f->Get(key) ;
 331     sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
 332     TH3F* h3DMc= (TH3F*)f->Get(key) ;
 333     sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
 334     TH3F* h3DRd = (TH3F*)f->Get(key) ;
 335     sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
 336     TH3F* h3DMd= (TH3F*)f->Get(key) ;
 337
 338     TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
 339     Double_t e1 = hCen->Integral(11,11);
 340     Double_t e2 = hCen->Integral(12,12);
 341     Double_t e3 = hCen->Integral(13,13);
 342     Double_t e4 = hCen->Integral(14,15)/2.;
 343     Double_t e5 = hCen->Integral(16,20)/5.;
 344
 345 cout<<"w_11="<<e5/e1<<", w_12="<<e5/e2<<", w_13="<<e5/e3<<", w_14-15="<<e5/e4<<endl;
 346
 347     h3DR->Add(h3DRa,e5/e1); delete h3DRa;
 348     h3DM->Add(h3DMa,e5/e1); delete h3DMa;
 349     h3DR->Add(h3DRb,e5/e2); delete h3DRb;
 350     h3DM->Add(h3DMb,e5/e2); delete h3DMb;
 351     h3DR->Add(h3DRc,e5/e3); delete h3DRc;
 352     h3DM->Add(h3DMc,e5/e3); delete h3DMc;
 353     h3DR->Add(h3DRd,e5/e4); delete h3DRd;
 354     h3DM->Add(h3DMd,e5/e4); delete h3DMd;
 355   }
 356   if(cen==3){//centrality 20-30
 357     sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
 358     h3DR = (TH3F*)f->Get(key) ;
 359     sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
 360     h3DM= (TH3F*)f->Get(key) ;
 361   }
 362   if(cen==4){//centrality 30-40
 363     sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
 364     h3DR = (TH3F*)f->Get(key) ;
 365     sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
 366     h3DM= (TH3F*)f->Get(key) ;
 367   }
 368   if(cen==5){//centrality 40-50
 369     sprintf(key,"h%sMassPt%s%s%s_cen7",w,d,kind,PID) ;
 370     h3DR = (TH3F*)f->Get(key) ;
 371     sprintf(key,"hMiMassPt%s%s%s_cen7",d,kind2,PID) ;
 372     h3DM= (TH3F*)f->Get(key) ;
 373   }
 374
 375
 376   if(Real)return h3DR;
 377   else return h3DM;
 378 }
 379
 380 TH1F * GetCen(){
 381 cout<<inname<<endl;
 382   TFile * f = new TFile(inname) ;
 383   return (TH1F*)f->Get("hCentrality") ;
 384 }
 385
 386 TH2F * GetCenPHOS(){
 387   TFile * f = new TFile(inname) ;
 388   return (TH2F*)f->Get("hCenPHOS") ;
 389 }
 390
 391