]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/macros/MakeMomResHisto_Merged.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / MakeMomResHisto_Merged.C
CommitLineData
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
24using namespace std;
25
26void MakeMomResHisto_Merged(){
27
28 int MBmerged = 0;
29 //
30
31 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11.root","READ");
32 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R10.root","READ");
33 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_lego_L0p68R9andR8.root","READ");
34 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R7.root","READ");
35 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R6.root","READ");
36 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_lego_L0p68R11_PID1p5_and_FB5.root","READ");
37 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_chi3p1.root","READ");
38 TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_highKt.root","READ");
39 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17d_fix_genSignal_Rinv11.root","READ");
40 //TFile *infile1 = new TFile("Results/PDC_HIJING_12a17a_myRun_L0p68R11_ncls.root","READ");
41 //TFile *infile1 = new TFile("MyOutput.root","READ");
42 //
43 //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_12a17e_noPadCut_lego.root","READ");
44 //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_newr3_lego.root","READ");
45 //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_GRS_lego.root","READ");
46 //TFile *infile2 = new TFile("Results/PDC_HIJING_12a17b_myRun_L0p4R8.root","READ");
47
48 //TDirectoryFile *tdir1 = (TDirectoryFile*)infile1->Get("PWGCF.outputChaoticityAnalysis.root");
49 //TList *list1=(TList*)tdir1->Get("ChaoticityOutput_2");
50 TList *list1=(TList*)infile1->Get("MyList");
51 infile1->Close();
52 //TDirectoryFile *tdir2 = (TDirectoryFile*)infile2->Get("PWGCF.outputChaoticityAnalysis.root");
53 //TList *list2=(TList*)tdir2->Get("ChaoticityOutput");
54 //TList *list2=(TList*)infile2->Get("MyList");
55 //infile2->Close();
56
57 TList *list;
58
59 const int Mbins=2;
60
61 TH2D *num_i_2[2];// 2-particle
62 TH2D *den_i_2[2];// 2-particle
63 TH2D *num_s_2[2];// 2-particle
64 TH2D *den_s_2[2];// 2-particle
65 //
66 TH2D *merged_num_i_2[2];
67 TH2D *merged_den_i_2[2];
68 TH2D *merged_num_s_2[2];
69 TH2D *merged_den_s_2[2];
70 //
71
72
73 TString *names[2][4];
74 // 3d
75 TString *names_terms[2][5][2];// charge-combo (ss or os), terms, ideal/smeared, Mbins
76 TH3D *merged_terms_i[2][5];// charge-combo (ss or os), terms, Mbins
77 TH3D *merged_terms_s[2][5];// charge-combo (ss or os), terms, Mbins
78 TH1D *merged_terms1D_i[2][5];
79 TH1D *merged_terms1D_s[2][5];
80 //
81 TString *names_termsK3[2][2];// charge-combo (ss or os), Sum/En, Mbins
82 TH3D *merged_terms_SumK3[2];// charge-combo (ss or os), Mbins
83 TH3D *merged_terms_EnK3[2];// charge-combo (ss or os), Mbins
84 //
85 TString *names_terms_FVP[5][2][2];// terms, ideal/smeared, Projection 1 or 2, Mbins
86 TH3D *merged_terms_FVP_i[5][2];// terms, Projection 1 or 2, Mbins
87 TH3D *merged_terms_FVP_s[5][2];// terms, Projection 1 or 2, Mbins
88 //
89 TString *names_terms_FVP_K[5][2][2];// terms, Sum/En, Projection 1 or 2, Mbins
90 TH3D *merged_terms_FVP_SumK[5][2];// terms, Projection 1 or 2, Mbins
91 TH3D *merged_terms_FVP_EnK[5][2];// terms, Projection 1 or 2, Mbins
92 //
93 TString *names_terms_FVPTPN[2][2];// ideal/smeared, Projection 1 or 2, Mbins
94 TH3D *merged_terms_FVPTPN_i[2];// Projection 1 or 2, Mbins
95 TH3D *merged_terms_FVPTPN_s[2];// Projection 1 or 2, Mbins
96
97
98 for(int mb=0; mb<Mbins; mb++){// Mbin loop
99 cout<<"Cent bin "<<mb<<endl;
100 for(int ch=0; ch<2; ch++){// + or - loop
101
102 int MixedCbin1=0, MixedCbin2=0, MixedCbin3=0;
103 // 3d histos
104 for(int term=0; term<5; term++){// term loop
105 names_terms[0][term][0]=new TString("PairCut3_Charge1_"); *names_terms[0][term][0] += ch;
106 names_terms[0][term][0]->Append("_Charge2_"); *names_terms[0][term][0] += ch;
107 names_terms[0][term][0]->Append("_Charge3_"); *names_terms[0][term][0] += ch;
108 names_terms[0][term][0]->Append("_SC_0_M_");
109 *names_terms[0][term][0] += mb;
110 names_terms[0][term][0]->Append("_ED_0_");
111 TString *TPNnameBase = new TString(names_terms[0][term][0]->Data());
112 names_terms[0][term][0]->Append("Term_"); *names_terms[0][term][0] += term+1;
113 names_terms[0][term][1]=new TString(names_terms[0][term][0]->Data());
114 if(term==0) {
115 names_termsK3[0][0]=new TString(names_terms[0][term][0]->Data());
116 names_termsK3[0][1]=new TString(names_terms[0][term][0]->Data());
117 names_termsK3[0][0]->Append("_SumK3");
118 names_termsK3[0][1]->Append("_EnK3");
119 }
120 //
121
122 names_terms_FVP[term][0][0] = new TString(names_terms[0][term][0]->Data());
123 names_terms_FVP[term][1][0] = new TString(names_terms[0][term][0]->Data());
124 names_terms_FVP[term][0][0]->Append("_4VectProd1Ideal");
125 names_terms_FVP[term][1][0]->Append("_4VectProd1Smeared");
126 //
127 names_terms_FVP[term][0][1] = new TString(names_terms[0][term][1]->Data());
128 names_terms_FVP[term][1][1] = new TString(names_terms[0][term][1]->Data());
129 names_terms_FVP[term][0][1]->Append("_4VectProd2Ideal");
130 names_terms_FVP[term][1][1]->Append("_4VectProd2Smeared");
131 //
132 names_terms_FVP_K[term][0][0] = new TString(names_terms[0][term][1]->Data());// Prod1 SumK
133 names_terms_FVP_K[term][1][0] = new TString(names_terms[0][term][1]->Data());// Prod1 EnK
134 names_terms_FVP_K[term][0][1] = new TString(names_terms[0][term][1]->Data());// Prod2 SumK
135 names_terms_FVP_K[term][1][1] = new TString(names_terms[0][term][1]->Data());// Prod2 EnK
136 if(term==0) {
137 names_terms_FVP_K[term][0][0]->Append("_4VectProd1SumK3");
138 names_terms_FVP_K[term][1][0]->Append("_4VectProd1EnK3");
139 names_terms_FVP_K[term][0][1]->Append("_4VectProd2SumK3");
140 names_terms_FVP_K[term][1][1]->Append("_4VectProd2EnK3");
141 names_terms_FVPTPN[0][0] = new TString(TPNnameBase->Data());
142 names_terms_FVPTPN[0][1] = new TString(TPNnameBase->Data());
143 names_terms_FVPTPN[1][0] = new TString(TPNnameBase->Data());
144 names_terms_FVPTPN[1][1] = new TString(TPNnameBase->Data());
145 names_terms_FVPTPN[0][0]->Append("TPN_0_4VectProd1Ideal");
146 names_terms_FVPTPN[0][1]->Append("TPN_0_4VectProd2Ideal");
147 names_terms_FVPTPN[1][0]->Append("TPN_0_4VectProd1Smeared");
148 names_terms_FVPTPN[1][1]->Append("TPN_0_4VectProd2Smeared");
149 }
150 if(term>0 && term<4){
151 names_terms_FVP_K[term][0][0]->Append("_4VectProd1SumK2");
152 names_terms_FVP_K[term][1][0]->Append("_4VectProd1EnK2");
153 names_terms_FVP_K[term][0][1]->Append("_4VectProd2SumK2");
154 names_terms_FVP_K[term][1][1]->Append("_4VectProd2EnK2");
155 }
156
157 names_terms[0][term][0]->Append("_Ideal");
158 names_terms[0][term][1]->Append("_Smeared");
159
160 //
161 if(ch==0) {MixedCbin1=0; MixedCbin2=0; MixedCbin3=1;}
162 else {MixedCbin1=0; MixedCbin2=1; MixedCbin3=1;}
163 names_terms[1][term][0]=new TString("PairCut3_Charge1_"); *names_terms[1][term][0] += MixedCbin1;
164 names_terms[1][term][0]->Append("_Charge2_"); *names_terms[1][term][0] += MixedCbin2;
165 names_terms[1][term][0]->Append("_Charge3_"); *names_terms[1][term][0] += MixedCbin3;
166 names_terms[1][term][0]->Append("_SC_0_M_");
167 *names_terms[1][term][0] += mb;
168 names_terms[1][term][0]->Append("_ED_0_Term_"); *names_terms[1][term][0] += term+1;
169 names_terms[1][term][1]=new TString(names_terms[1][term][0]->Data());
170 if(term==0) {
171 names_termsK3[1][0]=new TString(names_terms[1][term][0]->Data());
172 names_termsK3[1][1]=new TString(names_terms[1][term][0]->Data());
173 names_termsK3[1][0]->Append("_SumK3");
174 names_termsK3[1][1]->Append("_EnK3");
175 }
176 names_terms[1][term][0]->Append("_Ideal");
177 names_terms[1][term][1]->Append("_Smeared");
178 }// term loop
179
180
181 TString *base1 = new TString("Explicit2_Charge1_"); *base1 += ch;
182 base1->Append("_Charge2_"); *base1 += ch; base1->Append("_SC_0_M_");
183 names[0][0]=new TString(base1->Data());
184 *names[0][0] += mb;
185 names[0][0]->Append("_ED_0_Term_1_Ideal");
186 names[0][1]=new TString(base1->Data());
187 *names[0][1] += mb;
188 names[0][1]->Append("_ED_0_Term_2_Ideal");
189 names[0][2]=new TString(base1->Data());
190 *names[0][2] += mb;
191 names[0][2]->Append("_ED_0_Term_1_Smeared");
192 names[0][3]=new TString(base1->Data());
193 *names[0][3] += mb;
194 names[0][3]->Append("_ED_0_Term_2_Smeared");
195 // mixed charge 2-particle case is double counted here. It's ok (only the stat errors are wrong...not used anyway).
196 TString *base2 = new TString("Explicit2_Charge1_0_Charge2_1"); base2->Append("_SC_0_M_");
197 names[1][0]=new TString(base2->Data());
198 *names[1][0] += mb;
199 names[1][0]->Append("_ED_0_Term_1_Ideal");
200 names[1][1]=new TString(base2->Data());
201 *names[1][1] += mb;
202 names[1][1]->Append("_ED_0_Term_2_Ideal");
203 names[1][2]=new TString(base2->Data());
204 *names[1][2] += mb;
205 names[1][2]->Append("_ED_0_Term_1_Smeared");
206 names[1][3]=new TString(base2->Data());
207 *names[1][3] += mb;
208 names[1][3]->Append("_ED_0_Term_2_Smeared");
209 //
210 if(mb<2) list=(TList*)list1->Clone();
211 //else list=(TList*)list2->Clone();
212
213
214 for(int ChProd=0; ChProd<2; ChProd++){// charge product loop
215 //
216 num_i_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][0]->Data());
217 den_i_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][1]->Data());
218 num_s_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][2]->Data());
219 den_s_2[ChProd]=(TH2D*)list->FindObject(names[ChProd][3]->Data());
220 //
221
222 if(ch==0 && mb==0){
223 for(int term=0; term<5; term++){
224 merged_terms_i[ChProd][term] = (TH3D*)(list->FindObject(names_terms[ChProd][term][0]->Data()))->Clone();
225 merged_terms_s[ChProd][term] = (TH3D*)(list->FindObject(names_terms[ChProd][term][1]->Data()))->Clone();
226
227 if(term==0){
228 merged_terms_SumK3[ChProd] = (TH3D*)(list->FindObject(names_termsK3[ChProd][0]->Data()))->Clone();
229 merged_terms_EnK3[ChProd] = (TH3D*)(list->FindObject(names_termsK3[ChProd][1]->Data()))->Clone();
230 }
231 if(ChProd==0) {
232 //merged_terms_FVP_i[term][0] = (TH3D*)(list->FindObject(names_terms_FVP[term][0][0]->Data()))->Clone();// ideal, Prod 1
233 //merged_terms_FVP_s[term][0] = (TH3D*)(list->FindObject(names_terms_FVP[term][1][0]->Data()))->Clone();// smeared, Prod 1
234 //merged_terms_FVP_i[term][1] = (TH3D*)(list->FindObject(names_terms_FVP[term][0][1]->Data()))->Clone();// ideal, Prod 2
235 //merged_terms_FVP_s[term][1] = (TH3D*)(list->FindObject(names_terms_FVP[term][1][1]->Data()))->Clone();// smeared, Prod 2
236 if(term !=4){
237 //merged_terms_FVP_SumK[term][0] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][0][0]->Data()))->Clone();// Prod 1 Sum K
238 //merged_terms_FVP_EnK[term][0] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][1][0]->Data()))->Clone();// Prod 1 En K
239 //merged_terms_FVP_SumK[term][1] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][0][1]->Data()))->Clone();// Prod 2 Sum K
240 //merged_terms_FVP_EnK[term][1] = (TH3D*)(list->FindObject(names_terms_FVP_K[term][1][1]->Data()))->Clone();// Prod 2 En K
241 }
242 if(term==0){
243 //merged_terms_FVPTPN_i[0] = (TH3D*)(list->FindObject(names_terms_FVPTPN[0][0]->Data()))->Clone();// ideal, Prod 1
244 //merged_terms_FVPTPN_s[0] = (TH3D*)(list->FindObject(names_terms_FVPTPN[1][0]->Data()))->Clone();// smeared, Prod 1
245 //merged_terms_FVPTPN_i[1] = (TH3D*)(list->FindObject(names_terms_FVPTPN[0][1]->Data()))->Clone();// ideal, Prod 2
246 //merged_terms_FVPTPN_s[1] = (TH3D*)(list->FindObject(names_terms_FVPTPN[1][1]->Data()))->Clone();// smeared, Prod 2
247 }
248 }
249 }
250 }else{// ch==1 || mb!=0
251 for(int term=0; term<5; term++){
252 if(ChProd==0) {
253 merged_terms_i[0][term]->Add((TH3D*)(list->FindObject(names_terms[0][term][0]->Data())));
254 merged_terms_s[0][term]->Add((TH3D*)(list->FindObject(names_terms[0][term][1]->Data())));
255 if(term==0) {
256 merged_terms_SumK3[0]->Add((TH3D*)(list->FindObject(names_termsK3[0][0]->Data())));
257 merged_terms_EnK3[0]->Add((TH3D*)(list->FindObject(names_termsK3[0][1]->Data())));
258 /*merged_terms_FVPTPN_i[0]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[0][0]->Data())));
259 merged_terms_FVPTPN_i[1]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[0][1]->Data())));
260 merged_terms_FVPTPN_s[0]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[1][0]->Data())));
261 merged_terms_FVPTPN_s[1]->Add((TH3D*)(list->FindObject(names_terms_FVPTPN[1][1]->Data())));*/
262 }
263 /* merged_terms_FVP_i[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][0][0]->Data())));
264 merged_terms_FVP_s[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][1][0]->Data())));
265 merged_terms_FVP_i[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][0][1]->Data())));
266 merged_terms_FVP_s[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP[term][1][1]->Data())));*/
267
268 if(term<4){
269 /*merged_terms_FVP_SumK[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][0][0]->Data())));
270 merged_terms_FVP_EnK[term][0]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][1][0]->Data())));
271 merged_terms_FVP_SumK[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][0][1]->Data())));
272 merged_terms_FVP_EnK[term][1]->Add((TH3D*)(list->FindObject(names_terms_FVP_K[term][1][1]->Data())));*/
273 }
274 }else {// different placement of same-charge pair
275 int termOther=term;
276 if(ch==1){
277 if(term==0) termOther=0;
278 else if(term==1) termOther=3;
279 else if(term==2) termOther=2;
280 else if(term==3) termOther=1;
281 else termOther=4;
282 }
283 TH3D *temphisto_i=((TH3D*)(list->FindObject(names_terms[1][termOther][0]->Data())));
284 TH3D *temphisto_s=((TH3D*)(list->FindObject(names_terms[1][termOther][1]->Data())));
285 TH3D *temphisto_SumK3=((TH3D*)(list->FindObject(names_termsK3[1][0]->Data())));
286 TH3D *temphisto_EnK3=((TH3D*)(list->FindObject(names_termsK3[1][1]->Data())));
287 for(int bin1=1; bin1<merged_terms_i[1][term]->GetNbinsX(); bin1++){
288 for(int bin2=1; bin2<merged_terms_i[1][term]->GetNbinsY(); bin2++){
289 for(int bin3=1; bin3<merged_terms_i[1][term]->GetNbinsZ(); bin3++){
290
291 if(ch==0){
292 merged_terms_i[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_i[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_i->GetBinContent(bin1,bin2,bin3));
293 merged_terms_s[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_s[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_s->GetBinContent(bin1,bin2,bin3));
294 }else {
295 merged_terms_i[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_i[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_i->GetBinContent(bin3,bin2,bin1));
296 merged_terms_s[1][term]->SetBinContent(bin1,bin2,bin3, merged_terms_s[1][term]->GetBinContent(bin1,bin2,bin3)+temphisto_s->GetBinContent(bin3,bin2,bin1));
297 }
298
299 //
300 if(term==0){
301 //merged_terms_SumK3[1]->SetBinContent(bin1,bin2,bin3, merged_terms_SumK3[1]->GetBinContent(bin1,bin2,bin3)+temphisto_SumK3->GetBinContent(bin3,bin2,bin1));
302 //merged_terms_EnK3[1]->SetBinContent(bin1,bin2,bin3, merged_terms_EnK3[1]->GetBinContent(bin1,bin2,bin3)+temphisto_EnK3->GetBinContent(bin3,bin2,bin1));
303 }
304 }
305 }
306 }
307 }
308 }
309 }// else for 3-particle part
310
311 // 2-particles
312 if(ch==0 && mb==0) {
313 merged_num_i_2[ChProd] = (TH2D*)num_i_2[ChProd]->Clone();
314 merged_den_i_2[ChProd] = (TH2D*)den_i_2[ChProd]->Clone();
315 merged_num_s_2[ChProd] = (TH2D*)num_s_2[ChProd]->Clone();
316 merged_den_s_2[ChProd] = (TH2D*)den_s_2[ChProd]->Clone();
317 }else {
318 merged_num_i_2[ChProd]->Add(num_i_2[ChProd]);
319 merged_den_i_2[ChProd]->Add(den_i_2[ChProd]);
320 merged_num_s_2[ChProd]->Add(num_s_2[ChProd]);
321 merged_den_s_2[ChProd]->Add(den_s_2[ChProd]);
322
323 }
324 }// ch loop
325
326
327 }// mb loop
328
329 }// ChProd loop
330
331
332 cout<<"Start Writing to Output File"<<endl;
333 TFile *outfile = new TFile("MomResFile_temp.root","RECREATE");
334 TFile *outfileOffline = new TFile("MomResFile_Offline_temp.root","RECREATE");
335 int NbinsX = merged_num_i_2[0]->GetNbinsX();
336 int NbinsY = merged_num_i_2[0]->GetNbinsY();
337 float Xlow = merged_num_i_2[0]->GetXaxis()->GetBinLowEdge(1);
338 float Xhigh = merged_num_i_2[0]->GetXaxis()->GetBinUpEdge(NbinsX);
339 float Ylow = merged_num_i_2[0]->GetYaxis()->GetBinLowEdge(1);
340 float Yhigh = merged_num_i_2[0]->GetYaxis()->GetBinUpEdge(NbinsY);
341 TH2D *MomResHisto_pp = new TH2D("MomResHisto_pp","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
342 TH2D *MomResHisto_mp = new TH2D("MomResHisto_mp","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
343 //
344 TH2D *MomResHisto_pp_term1 = new TH2D("MomResHisto_pp_term1","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
345 TH2D *MomResHisto_pp_term2 = new TH2D("MomResHisto_pp_term2","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
346 TH2D *MomResHisto_mp_term1 = new TH2D("MomResHisto_mp_term1","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
347 TH2D *MomResHisto_mp_term2 = new TH2D("MomResHisto_mp_term2","",NbinsX,Xlow,Xhigh, NbinsY,Ylow,Yhigh);
348 //
349 // 3d part
350 int NbinsX_3d = merged_terms_i[0][0]->GetNbinsX();
351 float Xlow_3d = merged_terms_i[0][0]->GetXaxis()->GetBinLowEdge(1);
352 float Xhigh_3d = merged_terms_i[0][0]->GetXaxis()->GetBinUpEdge(NbinsX_3d);
353 TH3D *MomResHisto_3d_SC[5];
354 TH3D *MomResHisto_3d_MC[5];
355 TH1D *MomResHisto_1d_SC[5];
356 TH1D *MomResHisto_1d_MC[5];
357 TH3D *MomResHisto_3d_SC_K3;
358 TH3D *MomResHisto_3d_MC_K3;
359 // 4vectProd
360 /*const int NEdges_FVP_temp = merged_terms_FVP_i[0][0]->GetNbinsX() + 1;
361 const int NEdges_FVP = NEdges_FVP_temp;
362 double Edges[NEdges_FVP]={0};
363 for(int edg=0; edg<NEdges_FVP; edg++){
364 Edges[edg] = merged_terms_FVP_i[0][0]->GetXaxis()->GetBinLowEdge(edg+1);
365 //cout<<Edges[edg]<<endl;
366 }
367 TH3D *MomResHisto_3d_FVP1[5];
368 TH3D *MomResHisto_3d_FVP2[5];
369 TH3D *MomResHisto_3d_FVP1K[5];
370 TH3D *MomResHisto_3d_FVP2K[5];
371 TH3D *MomResHisto_TPN_FVP1;
372 TH3D *MomResHisto_TPN_FVP2;*/
373
374
375 for(int term=0; term<5; term++){
376 TString *nameSC = new TString("MomResHisto3_SC_term");
377 TString *nameMC = new TString("MomResHisto3_MC_term");
378 TString *nameSC1D = new TString("MomResHisto1D_SC_term");
379 TString *nameMC1D = new TString("MomResHisto1D_MC_term");
380 TString *nameFVP1 = new TString("MomResHisto3_FVP1_term");
381 TString *nameFVP2 = new TString("MomResHisto3_FVP2_term");
382 TString *nameFVP1K = new TString("MomResHisto3_FVP1K_term");
383 TString *nameFVP2K = new TString("MomResHisto3_FVP2K_term");
384 *nameSC += term+1;
385 *nameMC += term+1;
386 *nameSC1D += term+1;
387 *nameMC1D += term+1;
388 *nameFVP1 += term+1;
389 *nameFVP2 += term+1;
390 *nameFVP1K += term+1;
391 *nameFVP2K += term+1;
392 nameSC->Append("_M");
393 nameMC->Append("_M");
394 nameSC1D->Append("_M");
395 nameMC1D->Append("_M");
396 nameFVP1->Append("_M");
397 nameFVP2->Append("_M");
398 nameFVP1K->Append("_M");
399 nameFVP2K->Append("_M");
400 *nameSC += MBmerged;
401 *nameMC += MBmerged;
402 *nameSC1D += MBmerged;
403 *nameMC1D += MBmerged;
404 *nameFVP1 += MBmerged;
405 *nameFVP2 += MBmerged;
406 *nameFVP1K += MBmerged;
407 *nameFVP2K += MBmerged;
408 MomResHisto_3d_SC[term] = new TH3D(nameSC->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
409 MomResHisto_3d_MC[term] = new TH3D(nameMC->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
410 MomResHisto_1d_SC[term] = new TH1D(nameSC1D->Data(),"",20,0,0.2);
411 MomResHisto_1d_MC[term] = new TH1D(nameMC1D->Data(),"",20,0,0.2);
412 /*MomResHisto_3d_FVP1[term] = new TH3D(nameFVP1->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
413 MomResHisto_3d_FVP2[term] = new TH3D(nameFVP2->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
414 MomResHisto_3d_FVP1K[term] = new TH3D(nameFVP1K->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
415 MomResHisto_3d_FVP2K[term] = new TH3D(nameFVP2K->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);*/
416 }
417 //TString *nameTPN_FVP1 = new TString("MomResHisto3_TPN_FVP1_M");
418 //TString *nameTPN_FVP2 = new TString("MomResHisto3_TPN_FVP2_M");
419 TString *nameSC_K3 = new TString("AvgK3_SC_M");
420 TString *nameMC_K3 = new TString("AvgK3_MC_M");
421 *nameSC_K3 += MBmerged;
422 *nameMC_K3 += MBmerged;
423 //*nameTPN_FVP1 += MBmerged;
424 //*nameTPN_FVP2 += MBmerged;
425 MomResHisto_3d_SC_K3 = new TH3D(nameSC_K3->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
426 MomResHisto_3d_MC_K3 = new TH3D(nameMC_K3->Data(),"",NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d, NbinsX_3d,Xlow_3d,Xhigh_3d);
427 //MomResHisto_TPN_FVP1 = new TH3D(nameTPN_FVP1->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
428 //MomResHisto_TPN_FVP2 = new TH3D(nameTPN_FVP2->Data(),"",NEdges_FVP-1,Edges, NEdges_FVP-1,Edges, NEdges_FVP-1,Edges);
429
430
431
432 for(int term=0; term<5; term++){
433
434 // 1D Q3 projection
435 double Sum_i_SC[20]={0};
436 double Sum_i_MC[20]={0};
437 double Sum_s_SC[20]={0};
438 double Sum_s_MC[20]={0};
439
440 for(int i=1; i<NbinsX_3d; i++){
441 for(int j=1; j<NbinsX_3d; j++){
442 for(int l=1; l<NbinsX_3d; l++){
443 //if((i+j)==2 || (i+l)==2 || (j+l)==2) continue;
444 double q3 = sqrt(pow(merged_terms_i[0][term]->GetXaxis()->GetBinCenter(i+1),2)+pow(merged_terms_i[0][term]->GetYaxis()->GetBinCenter(j+1),2)+pow(merged_terms_i[0][term]->GetZaxis()->GetBinCenter(l+1),2));
445 int q3bin = MomResHisto_1d_SC[term]->GetXaxis()->FindBin(q3);
446 if(q3bin>20) continue;
447 Sum_i_SC[q3bin-1] += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1);
448 Sum_i_MC[q3bin-1] += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1);
449 Sum_s_SC[q3bin-1] += merged_terms_s[0][term]->GetBinContent(i+1,j+1,l+1);
450 Sum_s_MC[q3bin-1] += merged_terms_s[1][term]->GetBinContent(i+1,j+1,l+1);
451 }
452 }
453 }
454 for(int i=0; i<20; i++){
455
456 if(Sum_s_SC[i]>0) MomResHisto_1d_SC[term]->SetBinContent(i+1, Sum_i_SC[i]/Sum_s_SC[i]);
457 if(Sum_s_MC[i]>0) MomResHisto_1d_MC[term]->SetBinContent(i+1, Sum_i_MC[i]/Sum_s_MC[i]);
458 }
459
460 // full 3d method
461 merged_terms_i[0][term]->Divide(merged_terms_s[0][term]);
462 merged_terms_i[1][term]->Divide(merged_terms_s[1][term]);
463 //merged_terms_FVP_i[term][0]->Divide(merged_terms_FVP_s[term][0]);// Prod 1
464 //merged_terms_FVP_i[term][1]->Divide(merged_terms_FVP_s[term][1]);// Prod 2
465 if(term < 4){
466 //merged_terms_FVP_SumK[term][0]->Divide(merged_terms_FVP_EnK[term][0]);// Prod 1
467 //merged_terms_FVP_SumK[term][1]->Divide(merged_terms_FVP_EnK[term][1]);// Prod 2
468 }
469 }
470 //merged_terms_FVPTPN_i[0]->Divide(merged_terms_FVPTPN_s[0]);// Prod 1
471 //merged_terms_FVPTPN_i[1]->Divide(merged_terms_FVPTPN_s[1]);// Prod 2
472 merged_terms_SumK3[0]->Divide(merged_terms_EnK3[0]);// Avg SC K3 for Q12,Q13,Q23
473 merged_terms_SumK3[1]->Divide(merged_terms_EnK3[1]);// Avg MC K3 for Q12,Q13,Q23
474
475
476
477 for(int term=0; term<5; term++){
478 for(int i=0; i<NbinsX_3d; i++){
479 for(int j=0; j<NbinsX_3d; j++){
480 for(int l=0; l<NbinsX_3d; l++){
481
482 double MomResCorr=0; double terms=0;
483 if(term==0 || term==4){
484 if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
485 if(merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1); terms++;}
486 if(merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1); terms++;}
487 if(merged_terms_i[0][term]->GetBinContent(j+1,l+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,l+1,i+1); terms++;}
488 if(merged_terms_i[0][term]->GetBinContent(l+1,i+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,i+1,j+1); terms++;}
489 if(merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1); terms++;}
490 }else if(term==1){
491 if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
492 if(merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,l+1,j+1); terms++;}
493 }else if(term==2){
494 if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
495 if(merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(l+1,j+1,i+1); terms++;}
496 }else {
497 if(merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(i+1,j+1,l+1); terms++;}
498 if(merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1)>0) {MomResCorr += merged_terms_i[0][term]->GetBinContent(j+1,i+1,l+1); terms++;}
499 }
500
501 if(terms > 0) MomResCorr /= terms;
502 MomResHisto_3d_SC[term]->SetBinContent(i+1,j+1,l+1, MomResCorr);
503
504 // Mixed-Charge
505 MomResCorr=0; terms=0;
506 if(term==0 || term==1 || term==4){
507 if(merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1); terms++;}
508 if(merged_terms_i[1][term]->GetBinContent(i+1,l+1,j+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,l+1,j+1); terms++;}
509 }else{
510 if(merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1)>0) {MomResCorr += merged_terms_i[1][term]->GetBinContent(i+1,j+1,l+1); terms++;}
511 }
512
513 if(terms > 0) MomResCorr /= terms;
514 MomResHisto_3d_MC[term]->SetBinContent(i+1,j+1,l+1, MomResCorr);
515
516 if(term==0){
517 MomResHisto_3d_SC_K3->SetBinContent(i+1,j+1,l+1, merged_terms_SumK3[0]->GetBinContent(i+1,j+1,l+1));
518 MomResHisto_3d_MC_K3->SetBinContent(i+1,j+1,l+1, merged_terms_SumK3[1]->GetBinContent(i+1,j+1,l+1));
519 }
520 // condition edge effects
521 if(MomResHisto_3d_SC[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_SC[term]->SetBinContent(i+1,j+1,l+1, 1.0);
522 if(MomResHisto_3d_MC[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_MC[term]->SetBinContent(i+1,j+1,l+1, 1.0);
523 }
524 }
525 }
526 }// term
527
528 // FVP
529 /*for(int mb=0; mb<1; mb++){
530 for(int term=0; term<5; term++){
531 for(int i=0; i<NEdges_FVP-1; i++){
532 for(int j=0; j<NEdges_FVP-1; j++){
533 for(int l=0; l<NEdges_FVP-1; l++){
534 MomResHisto_3d_FVP1[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_i[term][0]->GetBinContent(i+1,j+1,l+1));
535 MomResHisto_3d_FVP2[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_i[term][1]->GetBinContent(i+1,j+1,l+1));
536 if(term < 4){
537 MomResHisto_3d_FVP1K[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_SumK[term][0]->GetBinContent(i+1,j+1,l+1));
538 MomResHisto_3d_FVP2K[term]->SetBinContent(i+1,j+1,l+1, merged_terms_FVP_SumK[term][1]->GetBinContent(i+1,j+1,l+1));
539 }
540 if(term==0){
541 MomResHisto_TPN_FVP1->SetBinContent(i+1,j+1,l+1, merged_terms_FVPTPN_i[0]->GetBinContent(i+1,j+1,l+1));
542 MomResHisto_TPN_FVP2->SetBinContent(i+1,j+1,l+1, merged_terms_FVPTPN_i[1]->GetBinContent(i+1,j+1,l+1));
543 }
544 // condition edge effects
545 if(MomResHisto_3d_FVP1[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_FVP1[term]->SetBinContent(i+1,j+1,l+1, 1.0);
546 if(MomResHisto_3d_FVP2[term]->GetBinContent(i+1,j+1,l+1) > 10) MomResHisto_3d_FVP2[term]->SetBinContent(i+1,j+1,l+1, 1.0);
547 }
548 }
549 }
550 }
551 }*/
552
553 // 2-particles
554 for(int ChProd=0; ChProd<2; ChProd++){
555 for(int i=0; i<NbinsX; i++){
556 for(int j=0; j<NbinsY; j++){
557
558 double weight = 1.0, weight_term1=1.0, weight_term2=1.0;
559 double weight_e = 0, weight_term1_e=0, weight_term2_e=0;
560 bool goodbin=kTRUE;
561 if(merged_num_i_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
562 if(merged_den_i_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
563 if(merged_num_s_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
564 if(merged_den_s_2[ChProd]->GetBinContent(i+1,j+1) == 0) goodbin=kFALSE;
565
566 if(goodbin==kTRUE){
567 double NUM_i=merged_num_i_2[ChProd]->GetBinContent(i+1,j+1);
568 double NUM_s=merged_num_s_2[ChProd]->GetBinContent(i+1,j+1);
569 double DEN_i=merged_den_i_2[ChProd]->GetBinContent(i+1,j+1);
570 double DEN_s=merged_den_s_2[ChProd]->GetBinContent(i+1,j+1);
571 weight = (NUM_i/DEN_i)/(NUM_s/DEN_s);
572 weight_e = pow((sqrt(NUM_s)*NUM_i/DEN_i)/(pow(NUM_s,2)/DEN_s),2);// approximate. does not take lambda into account!!
573 weight_e += pow((sqrt(DEN_s)*NUM_i/DEN_i)/(NUM_s),2);
574 weight_e = sqrt(weight_e);
575 //
576 weight_term1 = NUM_i/NUM_s;
577 weight_term1_e = pow(sqrt(NUM_i)/NUM_s,2);
578 weight_term1_e += pow(sqrt(NUM_s)*NUM_i/pow(NUM_s,2),2);
579 weight_term1_e = sqrt(weight_term1_e);
580 weight_term2 = DEN_i/DEN_s;
581 weight_term2_e = pow(sqrt(DEN_i)/DEN_s,2);
582 weight_term2_e += pow(sqrt(DEN_s)*DEN_i/pow(DEN_s,2),2);
583 weight_term2_e = sqrt(weight_term2_e);
584 }
585
586 if(ChProd==0) {
587 MomResHisto_pp->SetBinContent(i+1,j+1, weight);
588 MomResHisto_pp->SetBinError(i+1,j+1, weight_e);
589 MomResHisto_pp_term1->SetBinContent(i+1,j+1, weight_term1);
590 MomResHisto_pp_term2->SetBinContent(i+1,j+1, weight_term2);
591 MomResHisto_pp_term1->SetBinError(i+1,j+1, weight_term1_e);
592 MomResHisto_pp_term2->SetBinError(i+1,j+1, weight_term2_e);
593 }else {
594 MomResHisto_mp->SetBinContent(i+1,j+1, weight);
595 MomResHisto_mp->SetBinError(i+1,j+1, weight_e);
596 MomResHisto_mp_term1->SetBinContent(i+1,j+1, weight_term1);
597 MomResHisto_mp_term2->SetBinContent(i+1,j+1, weight_term2);
598 MomResHisto_mp_term1->SetBinError(i+1,j+1, weight_term1_e);
599 MomResHisto_mp_term2->SetBinError(i+1,j+1, weight_term2_e);
600 }
601 }
602 }
603 }
604
605
606
607 outfile->cd();
608 MomResHisto_pp->Write("MomResHisto_pp");
609 MomResHisto_mp->Write("MomResHisto_mp");
610 MomResHisto_pp_term1->Write("MomResHisto_pp_term1");
611 MomResHisto_pp_term2->Write("MomResHisto_pp_term2");
612 MomResHisto_mp_term1->Write("MomResHisto_mp_term1");
613 MomResHisto_mp_term2->Write("MomResHisto_mp_term2");
614 outfile->Close();
615
616 outfileOffline->cd();
617 if(MBmerged==0){
618 MomResHisto_pp->Write("MomResHisto_pp");
619 MomResHisto_mp->Write("MomResHisto_mp");
620 MomResHisto_pp_term1->Write("MomResHisto_pp_term1");
621 MomResHisto_pp_term2->Write("MomResHisto_pp_term2");
622 MomResHisto_mp_term1->Write("MomResHisto_mp_term1");
623 MomResHisto_mp_term2->Write("MomResHisto_mp_term2");
624 }
625 //
626
627 for(int term=0; term<5; term++){
628 MomResHisto_3d_SC[term]->Write();
629 MomResHisto_3d_MC[term]->Write();
630 MomResHisto_1d_SC[term]->Write();
631 MomResHisto_1d_MC[term]->Write();
632 //MomResHisto_3d_FVP1[term]->Write();
633 //MomResHisto_3d_FVP2[term]->Write();
634 if(term < 4){
635 //MomResHisto_3d_FVP1K[term]->Write();
636 //MomResHisto_3d_FVP2K[term]->Write();
637 }
638 }
639 MomResHisto_3d_SC_K3->Write();
640 MomResHisto_3d_MC_K3->Write();
641 //MomResHisto_TPN_FVP1->Write();
642 //MomResHisto_TPN_FVP2->Write();
643
644
645 outfileOffline->Close();
646
647}