[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / macros / Flow / Methods.h

char inname[256];

//  const Int_t nbin=11 ;
//  Double_t xa[14]={1.,1.5,2.,2.5,3.,4.,5.,6.,8.,10.,14.,20.} ;
  const Int_t nbin=12 ;
  Double_t xa[13]={0.6,1.,1.5,2.,3.,4.,5.,7.,9.,11.,13.,16.,20.} ;
//  const Int_t nbin=11 ;
//  Double_t xa[12]={0.6,1.,1.5,2.,2.5,3.,4.,5.,7.,10.,15.,20.} ;
    

Double_t Ollitrault(Double_t chi){
  Double_t x = 0.25*chi*chi;
  Double_t resk1 = 0.626657 * chi * exp(-x) * (TMath::BesselI0((float)x) + TMath::BesselI1((float)x));
  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));

  return resk1;
}

Double_t GetCos(TH2F* hcos2AC, Int_t cen){
  TH1D *hres = 0 ;
  TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  Double_t resMean = 0, weight=0. ;

  if(cen==21){
  TFile * f3 = new TFile("../flow11h_Apr14.root") ;
  hPHOS = (TH2F*)f3->Get("hCenPHOS");
  }

Int_t bin0, bin1;

if(cen==10){ bin0=1; bin1=3;}
if(cen==11){ bin0=5; bin1=9;}
if(cen==0){ bin0=1; bin1=2;}
if(cen==1){ bin0=2; bin1=3;}
if(cen==2){ bin0=3; bin1=5;}
if(cen==3){ bin0=5; bin1=7;}
if(cen==4){ bin0=7; bin1=9;}
if(cen==5){ bin0=9; bin1=11;}
if(cen==20){ bin0=3; bin1=11;} //10-50%
if(cen==21){ bin0=1; bin1=5;} //0-20%

  for(Int_t i=bin0; i<bin1; i++){
    TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ; 
    Double_t wi = hwi->GetMean() ;
    hres = hcos2AC->ProjectionX("res",i,i) ; 
    resMean += wi*hres->GetMean() ;
    weight += wi ;
    delete hwi ;
    delete hres ;
  }
  if(weight>0.){
    resMean/=weight ;
    return resMean ;
  }
  else{
    return 0. ;
  }
}
Double_t GetRes(TH2F* hcos2AC, Int_t cen){

  TH1D *hres = 0 ;
  Double_t resMean = 0, weight=0. ;

Int_t bin0, bin1;

if(cen==10){ bin0=1; bin1=3;} //0-10 (Central)
if(cen==11){ bin0=5; bin1=9;} //20-40 (SemiCentral)
if(cen==0){ bin0=1; bin1=2;}
if(cen==1){ bin0=2; bin1=3;}
if(cen==2){ bin0=3; bin1=5;}
if(cen==3){ bin0=5; bin1=7;}
if(cen==4){ bin0=7; bin1=9;}
if(cen==5){ bin0=9; bin1=11;}
if(cen==20){ bin0=3; bin1=11;} //10-50%
if(cen==21){ bin0=1; bin1=5;} //0-20%

  TH1D *hres = 0 ;
  TH2F * hPHOS = (TH2F*)gROOT->FindObjectAny("hCenPHOS");
  Double_t resMean = 0, weight=0. ;
  if(cen==21){
  TFile * f3 = new TFile("../flow11h_Apr14.root") ;
  hPHOS = (TH2F*)f3->Get("hCenPHOS");
  }

  for(Int_t i=bin0; i<bin1; i++){
    TH1D * hwi = hPHOS->ProjectionX("wi",5*i-4,5*i) ; 
    Double_t wi = hwi->GetMean() ;
    hres = hcos2AC->ProjectionX("res",i,i) ; 
    resMean += wi*hres->GetMean() ;
    weight += wi ;
    delete hwi ;
    delete hres ;
  }
  if(weight>0.){
    resMean/=weight ;
  }
  else{
    resMean = 0. ;
  }

    
  double chi = 2;
  double resSub = TMath::Sqrt(resMean);
  for(int j=0; j<15; j++){
    double resEve = Ollitrault(chi);
    chi= (resEve < resSub) ? chi+1.0*pow(0.5, j) : chi-1.0*pow(0.5, j);
  }
  return Ollitrault(TMath::Sqrt(2.)*chi);
}

TH3F* GetRealMixed(Int_t cen, char w[255], char d[255], char kind[255], char kind2[255], char PID[255], Int_t Real=1){
  if(cen==10||cen==0||cen==1||cen==21) sprintf(inname,"../flow11hCentral_Apr14.root");
  else sprintf(inname,"../flow11hSemiCentral_Apr14.root");
  char inname2[255];
  if(cen==21) sprintf(inname2,"../flow11hSemiCentral_Apr14.root");

  TFile * f = new TFile(inname) ;
  if(cen==21) TFile * f2 = new TFile(inname2) ;

  char key[125] ;
  TH3F * h3DR =0 ;TH3F * h3DM=0;

  if(cen==10){//centrality 0-10 
    sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
    TH3F *h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
    TH3F *h3DMa= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
    h3DRb = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
    h3DMb= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
    h3DRc = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
    h3DMc= (TH3F*)f->Get(key) ;

    TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
    Double_t e1 = hCen->Integral(6,8)/3.;
    Double_t e2 = hCen->Integral(9,9);
    Double_t e3 = hCen->Integral(10,10);

cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;

    h3DR->Add(h3DRa); delete h3DRa;
    h3DM->Add(h3DMa); delete h3DMa;
    h3DR->Add(h3DRb,e1/e2); delete h3DRb;
    h3DM->Add(h3DMb,e1/e2); delete h3DMb;
    h3DR->Add(h3DRc,e1/e3); delete h3DRc;
    h3DM->Add(h3DMc,e1/e3); delete h3DMc;
  }
  if(cen==21){//centrality 0-20 
    sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
    TH3F *h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
    TH3F *h3DMa= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
    h3DRb = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
    h3DMb= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
    h3DRc = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
    h3DMc= (TH3F*)f->Get(key) ;

    TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
    Double_t e1 = hCen->Integral(6,8)/3.;
    Double_t e2 = hCen->Integral(9,9);
    Double_t e3 = hCen->Integral(10,10);

cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;

    h3DR->Add(h3DRa); delete h3DRa;
    h3DM->Add(h3DMa); delete h3DMa;
    h3DR->Add(h3DRb,e1/e2); delete h3DRb;
    h3DM->Add(h3DMb,e1/e2); delete h3DMb;
    h3DR->Add(h3DRc,e1/e3); delete h3DRc;
    h3DM->Add(h3DMc,e1/e3); delete h3DMc;

//now 10-20
    sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
    TH3F* h3DR1 = (TH3F*)f2->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
    TH3F* h3DM1= (TH3F*)f2->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
    TH3F* h3DRa = (TH3F*)f2->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
    TH3F* h3DMa= (TH3F*)f2->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
    TH3F* h3DRb = (TH3F*)f2->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
    TH3F* h3DMb= (TH3F*)f2->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
    TH3F* h3DRc = (TH3F*)f2->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
    TH3F* h3DMc= (TH3F*)f2->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
    TH3F* h3DRd = (TH3F*)f2->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
    TH3F* h3DMd= (TH3F*)f2->Get(key) ;

    TH1D* hCen2 = ((TH2F*)f2->Get("hCentrality"))->ProjectionX();
    Double_t e1 = hCen2->Integral(11,11);
    Double_t e2 = hCen2->Integral(12,12);
    Double_t e3 = hCen2->Integral(13,13);
    Double_t e4 = hCen2->Integral(14,15)/2.;
    Double_t e5 = hCen2->Integral(16,20)/5.;

cout<<"w_11="<<e5/e1<<", w_12="<<e5/e2<<", w_13="<<e5/e3<<", w_14-15="<<e5/e4<<endl;

    h3DR->Add(h3DR1,e5/e1); delete h3DR1;
    h3DM->Add(h3DM1,e5/e1); delete h3DM1;
    h3DR->Add(h3DRa,e5/e1); delete h3DRa;
    h3DM->Add(h3DMa,e5/e1); delete h3DMa;
    h3DR->Add(h3DRb,e5/e2); delete h3DRb;
    h3DM->Add(h3DMb,e5/e2); delete h3DMb;
    h3DR->Add(h3DRc,e5/e3); delete h3DRc;
    h3DM->Add(h3DMc,e5/e3); delete h3DMc;
    h3DR->Add(h3DRd,e5/e4); delete h3DRd;
  }
  if(cen==11){//centrality 20-40 
    sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
    
    sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
    TH3F *h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
    TH3F *h3DMa= (TH3F*)f->Get(key) ;

    h3DR->Add(h3DRa); delete h3DRa;
    h3DM->Add(h3DMa); delete h3DMa;
  }
  if(cen==20){//centrality 10-50
    sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
    TH3F *h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
    TH3F *h3DMa= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
    TH3F *h3DRb = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
    TH3F *h3DMb= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen7",w,d,kind,PID) ;
    TH3F *h3DRc = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen7",d,kind2,PID) ;
    TH3F *h3DMc= (TH3F*)f->Get(key) ;

    h3DR->Add(h3DRa); delete h3DRa;
    h3DM->Add(h3DMa); delete h3DMa;
    h3DR->Add(h3DRb); delete h3DRb;
    h3DM->Add(h3DMb); delete h3DMb;
    h3DR->Add(h3DRc); delete h3DRc;
    h3DM->Add(h3DMc); delete h3DMc;
  }
//small bins
  if(cen==0){//centrality 0-5 
    sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
  }
  if(cen==1){//centrality 5-10   
    sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
    h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
    h3DMa= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
    h3DRb = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
    h3DMb= (TH3F*)f->Get(key) ;

    TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
    Double_t e1 = hCen->Integral(6,8)/3.;
    Double_t e2 = hCen->Integral(9,9);
    Double_t e3 = hCen->Integral(10,10);

cout<<"w2="<<e1/e2<<", w3="<<e1/e3<<endl;

    h3DR->Add(h3DRa,e1/e2); delete h3DRa;
    h3DM->Add(h3DMa,e1/e2); delete h3DMa;
    h3DR->Add(h3DRb,e1/e3); delete h3DRb;
    h3DM->Add(h3DMb,e1/e3); delete h3DMb;

/*
    h3DR->Add(h3DRa); delete h3DRa;
    h3DM->Add(h3DMa); delete h3DMa;
    h3DR->Add(h3DRb); delete h3DRb;
    h3DM->Add(h3DMb); delete h3DMb;
*/
  }
  if(cen==2){//centrality 10-20   
    sprintf(key,"h%sMassPt%s%s%s_cen4",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen4",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;

    sprintf(key,"h%sMassPt%s%s%s_cen0",w,d,kind,PID) ;
    TH3F* h3DRa = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen0",d,kind2,PID) ;
    TH3F* h3DMa= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen1",w,d,kind,PID) ;
    TH3F* h3DRb = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen1",d,kind2,PID) ;
    TH3F* h3DMb= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen2",w,d,kind,PID) ;
    TH3F* h3DRc = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen2",d,kind2,PID) ;
    TH3F* h3DMc= (TH3F*)f->Get(key) ;
    sprintf(key,"h%sMassPt%s%s%s_cen3",w,d,kind,PID) ;
    TH3F* h3DRd = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen3",d,kind2,PID) ;
    TH3F* h3DMd= (TH3F*)f->Get(key) ;

    TH1D* hCen = ((TH2F*)f->Get("hCentrality"))->ProjectionX();
    Double_t e1 = hCen->Integral(11,11);
    Double_t e2 = hCen->Integral(12,12);
    Double_t e3 = hCen->Integral(13,13);
    Double_t e4 = hCen->Integral(14,15)/2.;
    Double_t e5 = hCen->Integral(16,20)/5.;

cout<<"w_11="<<e5/e1<<", w_12="<<e5/e2<<", w_13="<<e5/e3<<", w_14-15="<<e5/e4<<endl;

    h3DR->Add(h3DRa,e5/e1); delete h3DRa;
    h3DM->Add(h3DMa,e5/e1); delete h3DMa;
    h3DR->Add(h3DRb,e5/e2); delete h3DRb;
    h3DM->Add(h3DMb,e5/e2); delete h3DMb;
    h3DR->Add(h3DRc,e5/e3); delete h3DRc;
    h3DM->Add(h3DMc,e5/e3); delete h3DMc;
    h3DR->Add(h3DRd,e5/e4); delete h3DRd;
    h3DM->Add(h3DMd,e5/e4); delete h3DMd;
  }
  if(cen==3){//centrality 20-30   
    sprintf(key,"h%sMassPt%s%s%s_cen5",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen5",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
  }
  if(cen==4){//centrality 30-40   
    sprintf(key,"h%sMassPt%s%s%s_cen6",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen6",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
  }
  if(cen==5){//centrality 40-50   
    sprintf(key,"h%sMassPt%s%s%s_cen7",w,d,kind,PID) ;
    h3DR = (TH3F*)f->Get(key) ;
    sprintf(key,"hMiMassPt%s%s%s_cen7",d,kind2,PID) ;
    h3DM= (TH3F*)f->Get(key) ;
  }


  if(Real)return h3DR;
  else return h3DM;
}

TH1F * GetCen(){
cout<<inname<<endl;
  TFile * f = new TFile(inname) ;
  return (TH1F*)f->Get("hCentrality") ;
}

TH2F * GetCenPHOS(){
  TFile * f = new TFile(inname) ;
  return (TH2F*)f->Get("hCenPHOS") ;
}
Commit	Line	Data
e751e1d9	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.25chichi;
	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",5i-4,5i) ;
	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",5i-4,5i) ;
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.0pow(0.5, j) : chi-1.0pow(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