[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / MakeKFileTherm.C

#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <Riostream.h>
#include <complex>

#include "TVector2.h"
#include "TFile.h"
#include "TString.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TMath.h"
#include "TText.h"
#include "TRandom3.h"
#include "TArray.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TMinuit.h"

using namespace std;

void MakeKFileTherm(){
  // 2-particles

  const int kbVALUES=6;// 6 b values (2,3,5,7,8,9)

  TFile *probeFile=new TFile("Therm_FSI_b2.root","READ");// all files have same binning.
  TH2D *probeHisto = (TH2D*)probeFile->Get("K2_ss");// kt x qinv
  int binsY= probeHisto->GetNbinsY();
  double lowY = probeHisto->GetYaxis()->GetBinLowEdge(1);
  double highY = probeHisto->GetYaxis()->GetBinUpEdge(binsY);
  probeFile->Close();

  TFile *OutFile=new TFile("KFile_temp.root","RECREATE");
  TH2D *K2ssT = new TH2D("K2ssT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
  TH2D *K2osT = new TH2D("K2osT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
  TH3D *K2ssT_kt = new TH3D("K2ssT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
  TH3D *K2osT_kt = new TH3D("K2osT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
  K2ssT_kt->GetXaxis()->SetTitle("b bin"); K2osT_kt->GetXaxis()->SetTitle("b bin");
  K2ssT_kt->GetYaxis()->SetTitle("kt bin"); K2osT_kt->GetYaxis()->SetTitle("kt bin");
  K2ssT_kt->GetZaxis()->SetTitle("qinv"); K2osT_kt->GetZaxis()->SetTitle("qinv");

  for(int r=0; r<kbVALUES; r++){ 
    
    TString *nameT=new TString("Therm_FSI_b");
        
     
    if(r==0) {*nameT += 2;}
    if(r==1) {*nameT += 3;}
    if(r==2) {*nameT += 5;}
    if(r==3) {*nameT += 7;}
    if(r==4) {*nameT += 8;}
    if(r==5) {*nameT += 9;}
    
    nameT->Append(".root");
        
   

    //
    TFile *file_T=new TFile(nameT->Data(),"READ");
    TH2D *Num2_ss = (TH2D*)file_T->Get("K2_ss");
    TH2D *Den2_ss = (TH2D*)file_T->Get("PlaneWF_ss");
    TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY();// kt integrated
    TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY();// kt integrated
    TH2D *Num2_os = (TH2D*)file_T->Get("K2_os");
    TH2D *Den2_os = (TH2D*)file_T->Get("PlaneWF_os");
    TH1D *Num_os = (TH1D*)Num2_os->ProjectionY();// kt integrated
    TH1D *Den_os = (TH1D*)Den2_os->ProjectionY();// kt integrated
    Num_ss->Divide(Den_ss);
    Num_os->Divide(Den_os);
    for(int i=1; i<=binsY; i++){
      K2ssT->SetBinContent(r+1, i, Num_ss->GetBinContent(i));
      K2osT->SetBinContent(r+1, i, Num_os->GetBinContent(i));
    }
   
    // kt differential
    for(int ktbin=1; ktbin<=6; ktbin++){
      TString *name=new TString("pro_");
      *name += ktbin;
      *name += 1;
      TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
      *name += 2;
      TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
      *name += 3;
      TH1D *Num_os = (TH1D*)Num2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
      *name += 4;
      TH1D *Den_os = (TH1D*)Den2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
      Num_ss->Divide(Den_ss);
      Num_os->Divide(Den_os);
       for(int i=1; i<=binsY; i++){
	 K2ssT_kt->SetBinContent(r+1, ktbin, i, Num_ss->GetBinContent(i));
	 K2osT_kt->SetBinContent(r+1, ktbin, i, Num_os->GetBinContent(i));
       }
    }  
    
    // 3-particles

    TH3D *Num3_ss=(TH3D*)file_T->Get("K3ss_3D");
    TH3D *Den3_ss=(TH3D*)file_T->Get("PlaneWF3ss_3D");
    Num3_ss->Divide(Den3_ss);
    TString *OutNameSS = new TString("K3ss_");
    *OutNameSS += r;
    OutFile->cd();
    TH3D *K3ss = (TH3D*)Num3_ss->Clone();
    TH3D *temp3ss = (TH3D*)Num3_ss->Clone();
    for(int i=1; i<=K3ss->GetNbinsX(); i++){
      for(int j=1; j<=K3ss->GetNbinsY(); j++){
	for(int k=1; k<=K3ss->GetNbinsZ(); k++){

	  // GRS
	  //double GRS = K2ssT->GetBinContent(r+1, i) * K2ssT->GetBinContent(r+1, j) * K2ssT->GetBinContent(r+1, k);
	  //K3ss->SetBinContent(i,j,k, GRS);

	  // Omega0
	  if(temp3ss->GetBinContent(i,j,k) > 1.0) K3ss->SetBinContent(i,j,k, 1.0);
	  else{
	    double mean=0;
	    double terms=0;
	    if(temp3ss->GetBinContent(i,j,k) < 1.0) {mean += temp3ss->GetBinContent(i,j,k); terms++;}
	    if(temp3ss->GetBinContent(i,k,j) < 1.0) {mean += temp3ss->GetBinContent(i,k,j); terms++;}
	    if(temp3ss->GetBinContent(j,i,k) < 1.0) {mean += temp3ss->GetBinContent(j,i,k); terms++;}
	    if(temp3ss->GetBinContent(j,k,i) < 1.0) {mean += temp3ss->GetBinContent(j,k,i); terms++;}
	    if(temp3ss->GetBinContent(k,i,j) < 1.0) {mean += temp3ss->GetBinContent(k,i,j); terms++;}
	    if(temp3ss->GetBinContent(k,j,i) < 1.0) {mean += temp3ss->GetBinContent(k,j,i); terms++;}
	    
	    if(terms > 0) {mean /= terms; K3ss->SetBinContent(i,j,k, mean);}
	    else K3ss->SetBinContent(i,j,k, 0);
	    }

	}
      }
    }
    K3ss->Write(OutNameSS->Data());
    //
    TH3D *Num3_os=(TH3D*)file_T->Get("K3os_3D");
    TH3D *Den3_os=(TH3D*)file_T->Get("PlaneWF3os_3D");
    Num3_os->Divide(Den3_os);
    //
    TString *OutNameOS = new TString("K3os_");
    *OutNameOS += r;
    OutFile->cd();
    TH3D *K3os = (TH3D*)Num3_os->Clone();
    TH3D *temp3os = (TH3D*)Num3_os->Clone();
    for(int i=1; i<=K3os->GetNbinsX(); i++){
      for(int j=1; j<=K3os->GetNbinsY(); j++){
	for(int k=1; k<=K3os->GetNbinsZ(); k++){

	  // GRS
	  //double GRS = K2ssT->GetBinContent(r+1, i) * K2osT->GetBinContent(r+1, j) * K2osT->GetBinContent(r+1, k);
	  //K3os->SetBinContent(i,j,k, GRS);
	  
	  // Omega0
	  if(temp3os->GetBinContent(i,j,k) > 3.0) K3os->SetBinContent(i,j,k, 1.0);
	  else {
	    double mean=0;
	    double terms=0;
	    if(temp3os->GetBinContent(i,j,k) < 3.0) {mean += temp3os->GetBinContent(i,j,k); terms++;}
	    if(temp3os->GetBinContent(i,k,j) < 3.0) {mean += temp3os->GetBinContent(i,k,j); terms++;}
	    
	    if(terms > 0) {mean /= terms; K3os->SetBinContent(i,j,k, mean);}
	    else K3os->SetBinContent(i,j,k, 0);
	    }
	  
	}
      }
    }
    
    K3os->Write(OutNameOS->Data());
    
    
    file_T->Close();
    
  }// r loop
  
  //
  OutFile->cd();
  //
  K2ssT->Write("K2ssT");
  K2osT->Write("K2osT");
  K2ssT_kt->Write("K2ssT_kt");
  K2osT_kt->Write("K2osT_kt");
  //
  
  OutFile->Close();

}
Commit	Line	Data
9129699a	1	#include <math.h>
	2	#include <time.h>
	3	#include <stdio.h>
	4	#include <stdlib.h>
	5	#include <Riostream.h>
	6	#include <complex>
	7
	8	#include "TVector2.h"
	9	#include "TFile.h"
	10	#include "TString.h"
	11	#include "TF1.h"
	12	#include "TH1.h"
	13	#include "TH2.h"
	14	#include "TH3.h"
	15	#include "TProfile.h"
	16	#include "TProfile2D.h"
	17	#include "TMath.h"
	18	#include "TText.h"
	19	#include "TRandom3.h"
	20	#include "TArray.h"
	21	#include "TLegend.h"
	22	#include "TStyle.h"
	23	#include "TMinuit.h"
	24
	25	using namespace std;
	26
	27	void MakeKFileTherm(){
	28	// 2-particles
	29
	30	const int kbVALUES=6;// 6 b values (2,3,5,7,8,9)
	31
	32	TFile *probeFile=new TFile("Therm_FSI_b2.root","READ");// all files have same binning.
	33	TH2D probeHisto = (TH2D)probeFile->Get("K2_ss");// kt x qinv
	34	int binsY= probeHisto->GetNbinsY();
	35	double lowY = probeHisto->GetYaxis()->GetBinLowEdge(1);
	36	double highY = probeHisto->GetYaxis()->GetBinUpEdge(binsY);
	37	probeFile->Close();
	38
	39	TFile *OutFile=new TFile("KFile_temp.root","RECREATE");
	40	TH2D *K2ssT = new TH2D("K2ssT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
	41	TH2D *K2osT = new TH2D("K2osT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
	42	TH3D *K2ssT_kt = new TH3D("K2ssT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
	43	TH3D *K2osT_kt = new TH3D("K2osT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
	44	K2ssT_kt->GetXaxis()->SetTitle("b bin"); K2osT_kt->GetXaxis()->SetTitle("b bin");
	45	K2ssT_kt->GetYaxis()->SetTitle("kt bin"); K2osT_kt->GetYaxis()->SetTitle("kt bin");
	46	K2ssT_kt->GetZaxis()->SetTitle("qinv"); K2osT_kt->GetZaxis()->SetTitle("qinv");
	47
	48	for(int r=0; r<kbVALUES; r++){
	49
	50	TString *nameT=new TString("Therm_FSI_b");
	51
	52
	53	if(r==0) {*nameT += 2;}
	54	if(r==1) {*nameT += 3;}
	55	if(r==2) {*nameT += 5;}
	56	if(r==3) {*nameT += 7;}
	57	if(r==4) {*nameT += 8;}
	58	if(r==5) {*nameT += 9;}
	59
	60	nameT->Append(".root");
	61
	62
	63
	64	//
65	TFile *file_T=new TFile(nameT->Data(),"READ");
66	TH2D Num2_ss = (TH2D)file_T->Get("K2_ss");
67	TH2D Den2_ss = (TH2D)file_T->Get("PlaneWF_ss");
68	TH1D Num_ss = (TH1D)Num2_ss->ProjectionY();// kt integrated
69	TH1D Den_ss = (TH1D)Den2_ss->ProjectionY();// kt integrated
70	TH2D Num2_os = (TH2D)file_T->Get("K2_os");
71	TH2D Den2_os = (TH2D)file_T->Get("PlaneWF_os");
72	TH1D Num_os = (TH1D)Num2_os->ProjectionY();// kt integrated
73	TH1D Den_os = (TH1D)Den2_os->ProjectionY();// kt integrated
74	Num_ss->Divide(Den_ss);
75	Num_os->Divide(Den_os);
76	for(int i=1; i<=binsY; i++){
77	K2ssT->SetBinContent(r+1, i, Num_ss->GetBinContent(i));
78	K2osT->SetBinContent(r+1, i, Num_os->GetBinContent(i));
79	}
80
81	// kt differential
82	for(int ktbin=1; ktbin<=6; ktbin++){
83	TString *name=new TString("pro_");
84	*name += ktbin;
85	*name += 1;
86	TH1D Num_ss = (TH1D)Num2_ss->ProjectionY(name->Data(), ktbin2+3, ktbin2+4);
87	*name += 2;
88	TH1D Den_ss = (TH1D)Den2_ss->ProjectionY(name->Data(), ktbin2+3, ktbin2+4);
89	*name += 3;
90	TH1D Num_os = (TH1D)Num2_os->ProjectionY(name->Data(), ktbin2+3, ktbin2+4);
91	*name += 4;
92	TH1D Den_os = (TH1D)Den2_os->ProjectionY(name->Data(), ktbin2+3, ktbin2+4);
93	Num_ss->Divide(Den_ss);
94	Num_os->Divide(Den_os);
95	for(int i=1; i<=binsY; i++){
96	K2ssT_kt->SetBinContent(r+1, ktbin, i, Num_ss->GetBinContent(i));
97	K2osT_kt->SetBinContent(r+1, ktbin, i, Num_os->GetBinContent(i));
98	}
99	}
100
101	// 3-particles
102
103	TH3D Num3_ss=(TH3D)file_T->Get("K3ss_3D");
104	TH3D Den3_ss=(TH3D)file_T->Get("PlaneWF3ss_3D");
105	Num3_ss->Divide(Den3_ss);
106	TString *OutNameSS = new TString("K3ss_");
107	*OutNameSS += r;
108	OutFile->cd();
109	TH3D K3ss = (TH3D)Num3_ss->Clone();
110	TH3D temp3ss = (TH3D)Num3_ss->Clone();
111	for(int i=1; i<=K3ss->GetNbinsX(); i++){
112	for(int j=1; j<=K3ss->GetNbinsY(); j++){
113	for(int k=1; k<=K3ss->GetNbinsZ(); k++){
114
115	// GRS
116	//double GRS = K2ssT->GetBinContent(r+1, i) * K2ssT->GetBinContent(r+1, j) * K2ssT->GetBinContent(r+1, k);
117	//K3ss->SetBinContent(i,j,k, GRS);
118
119	// Omega0
120	if(temp3ss->GetBinContent(i,j,k) > 1.0) K3ss->SetBinContent(i,j,k, 1.0);
121	else{
122	double mean=0;
123	double terms=0;
124	if(temp3ss->GetBinContent(i,j,k) < 1.0) {mean += temp3ss->GetBinContent(i,j,k); terms++;}
125	if(temp3ss->GetBinContent(i,k,j) < 1.0) {mean += temp3ss->GetBinContent(i,k,j); terms++;}
126	if(temp3ss->GetBinContent(j,i,k) < 1.0) {mean += temp3ss->GetBinContent(j,i,k); terms++;}
127	if(temp3ss->GetBinContent(j,k,i) < 1.0) {mean += temp3ss->GetBinContent(j,k,i); terms++;}
128	if(temp3ss->GetBinContent(k,i,j) < 1.0) {mean += temp3ss->GetBinContent(k,i,j); terms++;}
129	if(temp3ss->GetBinContent(k,j,i) < 1.0) {mean += temp3ss->GetBinContent(k,j,i); terms++;}
130
131	if(terms > 0) {mean /= terms; K3ss->SetBinContent(i,j,k, mean);}
132	else K3ss->SetBinContent(i,j,k, 0);
133	}
134
135	}
136	}
137	}
138	K3ss->Write(OutNameSS->Data());
139	//
140	TH3D Num3_os=(TH3D)file_T->Get("K3os_3D");
141	TH3D Den3_os=(TH3D)file_T->Get("PlaneWF3os_3D");
142	Num3_os->Divide(Den3_os);
143	//
144	TString *OutNameOS = new TString("K3os_");
145	*OutNameOS += r;
146	OutFile->cd();
147	TH3D K3os = (TH3D)Num3_os->Clone();
148	TH3D temp3os = (TH3D)Num3_os->Clone();
149	for(int i=1; i<=K3os->GetNbinsX(); i++){
150	for(int j=1; j<=K3os->GetNbinsY(); j++){
151	for(int k=1; k<=K3os->GetNbinsZ(); k++){
152
153	// GRS
154	//double GRS = K2ssT->GetBinContent(r+1, i) * K2osT->GetBinContent(r+1, j) * K2osT->GetBinContent(r+1, k);
155	//K3os->SetBinContent(i,j,k, GRS);
156
157	// Omega0
158	if(temp3os->GetBinContent(i,j,k) > 3.0) K3os->SetBinContent(i,j,k, 1.0);
159	else {
160	double mean=0;
161	double terms=0;
162	if(temp3os->GetBinContent(i,j,k) < 3.0) {mean += temp3os->GetBinContent(i,j,k); terms++;}
163	if(temp3os->GetBinContent(i,k,j) < 3.0) {mean += temp3os->GetBinContent(i,k,j); terms++;}
164
165	if(terms > 0) {mean /= terms; K3os->SetBinContent(i,j,k, mean);}
166	else K3os->SetBinContent(i,j,k, 0);
167	}
168
169	}
170	}
171	}
172
173	K3os->Write(OutNameOS->Data());
174
175
176	file_T->Close();
177
178	}// r loop
179
180	//
181	OutFile->cd();
182	//
183	K2ssT->Write("K2ssT");
184	K2osT->Write("K2osT");
185	K2ssT_kt->Write("K2ssT_kt");
186	K2osT_kt->Write("K2osT_kt");
187	//
188
189	OutFile->Close();
190
191	}