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

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

#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 MakeWeightFile()
{
  
  gStyle->SetOptStat(0);
  gStyle->SetOptDate(0);
  

  //TFile *InputFile = new TFile("Results/RawWeightFile_11h.root","READ");
  //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11.root","READ");
  //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11_Smeared.root","READ");
  //TFile *InputFile = new TFile("Results/RawWeightFile_11h_PID1p5.root","READ");
  //TFile *InputFile = new TFile("Results/RawWeightFile_FB5.root","READ");
  //TFile *InputFile = new TFile("Results/RawWeightFile_11h_Chi3p1.root","READ");
  //TFile *InputFile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// standard weights
  //TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_new.root","READ");
  TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_2sigma_3sigmaTTC.root","READ");// _1(2sigmaTTC), _2(3sigmaTTC)

  //TFile *InputFile = new TFile("MyOutput.root","READ");
  TDirectoryFile *tdir = (TDirectoryFile*)InputFile->Get("PWGCF.outputChaoticityAnalysis.root");
  TList *MyList=(TList*)tdir->Get("ChaoticityOutput_2");
  //TList *MyList=(TList*)InputFile->Get("MyList");
  InputFile->Close();

  const int KtBins=4;
  const int KyBins=1;
  const int EDBins=1;
  const int MBins=10;//10, make sure MB assignment below is correct!!!!!!!!!!!!!!!
  
  TFile *OutFile = new TFile("WeightFile_temp.root","RECREATE");
  TH3F *WeightHistos[KtBins][MBins];
  
  for(int ktB=1; ktB<=KtBins; ktB++){
    for(int MB=1; MB<=MBins; MB++){
      TString *InNameNum = new TString("TwoPart_num_Kt_");
      *InNameNum += ktB-1;
      InNameNum->Append("_Ky_0_M_");
      if(MB<=10) *InNameNum += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
      else *InNameNum += 1;
      InNameNum->Append("_ED_0");
      //
      TString *InNameDen = new TString("TwoPart_den_Kt_");
      *InNameDen += ktB-1;
      InNameDen->Append("_Ky_0_M_");
      if(MB<=10) *InNameDen += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
      else *InNameDen += 1;
      InNameDen->Append("_ED_0");
      //
      TH3D *tempNum = (TH3D*)MyList->FindObject(InNameNum->Data());
      TH3D *tempDen = (TH3D*)MyList->FindObject(InNameDen->Data());
      //
      int NormBinStart = tempNum->GetXaxis()->FindBin(0.135);//0.135
      int NormBinEnd = tempNum->GetXaxis()->FindBin(0.2);//0.2
      double Norm = tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
      Norm /= tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
      cout<<"Normalization = "<<Norm<<endl;
      cout<<tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<"  "<<tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<endl;
      //
      TString *OutNameWeight = new TString("Weight_Kt_");
      *OutNameWeight += ktB-1;
      OutNameWeight->Append("_Ky_0_M_");
      *OutNameWeight += MB-1;
      OutNameWeight->Append("_ED_0");
      
      int Nbins=tempNum->GetNbinsX();
      double QLimit = tempNum->GetXaxis()->GetBinUpEdge(Nbins);
      WeightHistos[ktB-1][MB-1] = new TH3F(OutNameWeight->Data(),"r3 Weights", Nbins,0,QLimit, Nbins,0,QLimit, Nbins,0,QLimit);
      WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitle("out");
      WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitle("side");
      WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitle("long");
      WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitleOffset(1.8);
      WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitleOffset(1.8);
      WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitleOffset(1.8);
      for(int outB=1; outB<=Nbins; outB++){
	for(int sideB=1; sideB<=Nbins; sideB++){
	  for(int longB=1; longB<=Nbins; longB++){
	    double weight=1, weight_e=0;
	    if(tempDen->GetBinContent(outB,sideB,longB) > 0 && tempNum->GetBinContent(outB,sideB,longB) > 0) {
	      //if(outB==sideB && outB==longB) cout<<outB<<"  "<<tempNum->GetBinContent(outB,sideB,longB)<<"  "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
	      weight = double(tempNum->GetBinContent(outB,sideB,longB))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm;
	      weight_e = pow(sqrt(double(tempNum->GetBinContent(outB,sideB,longB)))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm,2);
	      weight_e += pow(sqrt(double(tempDen->GetBinContent(outB,sideB,longB)))*double(tempNum->GetBinContent(outB,sideB,longB))/pow(double(tempDen->GetBinContent(outB,sideB,longB)),2) / Norm,2);
	      weight_e = sqrt(weight_e);
	      //if(weight < 0.8 && weight > 0.1) cout<<outB<<"  "<<sideB<<"  "<<longB<<"  "<<tempNum->GetBinContent(outB,sideB,longB)<<"  "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
	    }
	    if(weight==0){
	      WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, 0);
	      WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, 0);
	    }else {
	      WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, weight-1.0);// difference from unity
	      WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, weight_e);
	    }

	  }
	}
      }
      
      WeightHistos[ktB-1][MB-1]->Write();

    }// MB
  }// ktB

  int BOI1=2;
  int BOI2=2;
  TH3D *histoOI = WeightHistos[0][0]->Clone();
  histoOI->SetDirectory(0);
  //OutFile->Close();
  
  TH1D *pro = WeightHistos[0][0]->ProjectionY("pro",BOI1,BOI1,BOI2,BOI2);
  pro->Draw();

  /*TFile *projections=new TFile("projections.root","RECREATE");
  histoOI->GetXaxis()->SetRange(BOI1,BOI1);
  TH2D *proYZ = (TH2D*)histoOI->Project3D("yz");
  histoOI->GetXaxis()->SetRange(1,40);
  histoOI->GetYaxis()->SetRange(BOI1,BOI1);
  TH2D *proXZ = (TH2D*)histoOI->Project3D("xz");
  histoOI->GetYaxis()->SetRange(1,40);
  histoOI->GetZaxis()->SetRange(BOI1,BOI1);
  TH2D *proXY = (TH2D*)histoOI->Project3D("xy");
  proYZ->GetXaxis()->SetRangeUser(0,0.1);
  proYZ->GetYaxis()->SetRangeUser(0,0.1);
  proXZ->GetXaxis()->SetRangeUser(0,0.1);
  proXZ->GetYaxis()->SetRangeUser(0,0.1);
  proXY->GetXaxis()->SetRangeUser(0,0.1);
  proXY->GetYaxis()->SetRangeUser(0,0.1);
  proXY->Draw("lego");
  
  proYZ->Write();
  proXZ->Write();
  proXY->Write();
  */

}
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
	7	#include "TVector2.h"
	8	#include "TFile.h"
	9	#include "TString.h"
	10	#include "TF1.h"
	11	#include "TH1.h"
	12	#include "TH2.h"
	13	#include "TH3.h"
	14	#include "TProfile.h"
	15	#include "TProfile2D.h"
	16	#include "TMath.h"
	17	#include "TText.h"
	18	#include "TRandom3.h"
	19	#include "TArray.h"
	20	#include "TLegend.h"
	21	#include "TStyle.h"
	22	#include "TMinuit.h"
	23
	24	using namespace std;
	25
	26
	27	void MakeWeightFile()
	28	{
	29
	30	gStyle->SetOptStat(0);
	31	gStyle->SetOptDate(0);
	32
	33
	34	//TFile *InputFile = new TFile("Results/RawWeightFile_11h.root","READ");
	35	//TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11.root","READ");
	36	//TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11_Smeared.root","READ");
	37	//TFile *InputFile = new TFile("Results/RawWeightFile_11h_PID1p5.root","READ");
	38	//TFile *InputFile = new TFile("Results/RawWeightFile_FB5.root","READ");
	39	//TFile *InputFile = new TFile("Results/RawWeightFile_11h_Chi3p1.root","READ");
	40	//TFile *InputFile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// standard weights
	41	//TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_new.root","READ");
	42	TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_2sigma_3sigmaTTC.root","READ");// _1(2sigmaTTC), _2(3sigmaTTC)
	43
	44	//TFile *InputFile = new TFile("MyOutput.root","READ");
	45	TDirectoryFile tdir = (TDirectoryFile)InputFile->Get("PWGCF.outputChaoticityAnalysis.root");
	46	TList MyList=(TList)tdir->Get("ChaoticityOutput_2");
	47	//TList MyList=(TList)InputFile->Get("MyList");
	48	InputFile->Close();
	49
	50	const int KtBins=4;
	51	const int KyBins=1;
	52	const int EDBins=1;
	53	const int MBins=10;//10, make sure MB assignment below is correct!!!!!!!!!!!!!!!
	54
	55	TFile *OutFile = new TFile("WeightFile_temp.root","RECREATE");
	56	TH3F *WeightHistos[KtBins][MBins];
	57
	58	for(int ktB=1; ktB<=KtBins; ktB++){
	59	for(int MB=1; MB<=MBins; MB++){
	60	TString *InNameNum = new TString("TwoPart_num_Kt_");
	61	*InNameNum += ktB-1;
	62	InNameNum->Append("_Ky_0_M_");
	63	if(MB<=10) *InNameNum += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
	64	else *InNameNum += 1;
65	InNameNum->Append("_ED_0");
66	//
67	TString *InNameDen = new TString("TwoPart_den_Kt_");
68	*InNameDen += ktB-1;
69	InNameDen->Append("_Ky_0_M_");
70	if(MB<=10) *InNameDen += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
71	else *InNameDen += 1;
72	InNameDen->Append("_ED_0");
73	//
74	TH3D tempNum = (TH3D)MyList->FindObject(InNameNum->Data());
75	TH3D tempDen = (TH3D)MyList->FindObject(InNameDen->Data());
76	//
77	int NormBinStart = tempNum->GetXaxis()->FindBin(0.135);//0.135
78	int NormBinEnd = tempNum->GetXaxis()->FindBin(0.2);//0.2
79	double Norm = tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
80	Norm /= tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
81	cout<<"Normalization = "<<Norm<<endl;
82	cout<<tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<" "<<tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<endl;
83	//
84	TString *OutNameWeight = new TString("Weight_Kt_");
85	*OutNameWeight += ktB-1;
86	OutNameWeight->Append("_Ky_0_M_");
87	*OutNameWeight += MB-1;
88	OutNameWeight->Append("_ED_0");
89
90	int Nbins=tempNum->GetNbinsX();
91	double QLimit = tempNum->GetXaxis()->GetBinUpEdge(Nbins);
92	WeightHistos[ktB-1][MB-1] = new TH3F(OutNameWeight->Data(),"r3 Weights", Nbins,0,QLimit, Nbins,0,QLimit, Nbins,0,QLimit);
93	WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitle("out");
94	WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitle("side");
95	WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitle("long");
96	WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitleOffset(1.8);
97	WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitleOffset(1.8);
98	WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitleOffset(1.8);
99	for(int outB=1; outB<=Nbins; outB++){
100	for(int sideB=1; sideB<=Nbins; sideB++){
101	for(int longB=1; longB<=Nbins; longB++){
102	double weight=1, weight_e=0;
103	if(tempDen->GetBinContent(outB,sideB,longB) > 0 && tempNum->GetBinContent(outB,sideB,longB) > 0) {
104	//if(outB==sideB && outB==longB) cout<<outB<<" "<<tempNum->GetBinContent(outB,sideB,longB)<<" "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
105	weight = double(tempNum->GetBinContent(outB,sideB,longB))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm;
106	weight_e = pow(sqrt(double(tempNum->GetBinContent(outB,sideB,longB)))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm,2);
107	weight_e += pow(sqrt(double(tempDen->GetBinContent(outB,sideB,longB)))*double(tempNum->GetBinContent(outB,sideB,longB))/pow(double(tempDen->GetBinContent(outB,sideB,longB)),2) / Norm,2);
108	weight_e = sqrt(weight_e);
109	//if(weight < 0.8 && weight > 0.1) cout<<outB<<" "<<sideB<<" "<<longB<<" "<<tempNum->GetBinContent(outB,sideB,longB)<<" "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
110	}
111	if(weight==0){
112	WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, 0);
113	WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, 0);
114	}else {
115	WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, weight-1.0);// difference from unity
116	WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, weight_e);
117	}
118
119	}
120	}
121	}
122
123	WeightHistos[ktB-1][MB-1]->Write();
124
125	}// MB
126	}// ktB
127
128	int BOI1=2;
129	int BOI2=2;
130	TH3D *histoOI = WeightHistos[0][0]->Clone();
131	histoOI->SetDirectory(0);
132	//OutFile->Close();
133
134	TH1D *pro = WeightHistos[0][0]->ProjectionY("pro",BOI1,BOI1,BOI2,BOI2);
135	pro->Draw();
136
137	/TFile projections=new TFile("projections.root","RECREATE");
138	histoOI->GetXaxis()->SetRange(BOI1,BOI1);
139	TH2D proYZ = (TH2D)histoOI->Project3D("yz");
140	histoOI->GetXaxis()->SetRange(1,40);
141	histoOI->GetYaxis()->SetRange(BOI1,BOI1);
142	TH2D proXZ = (TH2D)histoOI->Project3D("xz");
143	histoOI->GetYaxis()->SetRange(1,40);
144	histoOI->GetZaxis()->SetRange(BOI1,BOI1);
145	TH2D proXY = (TH2D)histoOI->Project3D("xy");
146	proYZ->GetXaxis()->SetRangeUser(0,0.1);
147	proYZ->GetYaxis()->SetRangeUser(0,0.1);
148	proXZ->GetXaxis()->SetRangeUser(0,0.1);
149	proXZ->GetYaxis()->SetRangeUser(0,0.1);
150	proXY->GetXaxis()->SetRangeUser(0,0.1);
151	proXY->GetYaxis()->SetRangeUser(0,0.1);
152	proXY->Draw("lego");
153
154	proYZ->Write();
155	proXZ->Write();
156	proXY->Write();
157	*/
158
159	}