]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/drawBalanceFunction2DPsi.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 / EBYE / macros / drawBalanceFunction2DPsi.C
CommitLineData
095c0aaf 1const Int_t numberOfCentralityBins = 14;
2TString centralityArray[numberOfCentralityBins] = {"0-80","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3","92-8500","0-10000"};
e4391c02 3
62416a3a 4// reduced ranges for moments determination
5// (Average over all event classes of 3 sigma of Gauss + Pol0 Fit)
6Double_t rangePbPbEtaMedMom = 3. * 0.233;// PbPb, Delta Eta, 2<pT,assoc<3<pT,trig<4
7Double_t rangePbPbEtaHighMom = 3. * 0.096;// PbPb, Delta Eta, 3<pT,assoc<8<pT,trig<15
8Double_t rangePbPbPhiMedMom = 3. * 0.273;// PbPb, Delta Phi, 2<pT,assoc<3<pT,trig<4
9Double_t rangePbPbPhiHighMom = 3. * 0.097;// PbPb, Delta Eta, 3<pT,assoc<8<pT,trig<15
6acdbcb2 10
11const Int_t gRebin = 1;
12void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
eb63b883 13 Int_t gCentrality = 1,
5de9ad1a 14 Int_t gBit = -1,
15 const char* gCentralityEstimator = 0x0,
52daf7b2 16 Bool_t kShowShuffled = kFALSE,
17 Bool_t kShowMixed = kTRUE,
eb63b883 18 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
34616eda 19 Double_t vertexZMin = -10.,
20 Double_t vertexZMax = 10.,
6acdbcb2 21 Double_t ptTriggerMin = -1.,
22 Double_t ptTriggerMax = -1.,
23 Double_t ptAssociatedMin = -1.,
f0c5040c 24 Double_t ptAssociatedMax = -1.,
832b21d8 25 Bool_t kUseVzBinning = kTRUE,
26 Bool_t k2pMethod = kTRUE,
3de51a46 27 TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
28 Bool_t bToy = kFALSE)
d67eae53 29{
6acdbcb2 30 //Macro that draws the BF distributions for each centrality bin
31 //for reaction plane dependent analysis
32 //Author: Panos.Christakoglou@nikhef.nl
33 //Load the PWG2ebye library
4070f709 34 gSystem->Load("libANALYSIS");
35 gSystem->Load("libANALYSISalice");
36 gSystem->Load("libEventMixing");
37 gSystem->Load("libCORRFW");
38 gSystem->Load("libPWGTools");
39 gSystem->Load("libPWGCFebye");
6acdbcb2 40
9a0495b0 41 //gROOT->LoadMacro("~/SetPlotStyle.C");
42 //SetPlotStyle();
8795e993 43 gStyle->SetPalette(1,0);
15dd45a0 44
6acdbcb2 45 //Prepare the objects and return them
3de51a46 46 TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy);
52daf7b2 47 TList *listBFShuffled = NULL;
3de51a46 48 if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy);
52daf7b2 49 TList *listBFMixed = NULL;
3de51a46 50 if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy);
6acdbcb2 51 if(!listBF) {
52 Printf("The TList object was not created");
53 return;
54 }
55 else
27634c13 56 draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
34616eda 57 psiMin,psiMax,vertexZMin,vertexZMax,
f0c5040c 58 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
832b21d8 59 kUseVzBinning,k2pMethod,eventClass);
6acdbcb2 60}
61
62//______________________________________________________//
5365d1d7 63TList *GetListOfObjects(const char* filename,
64 Int_t gCentrality,
65 Int_t gBit,
66 const char *gCentralityEstimator,
3de51a46 67 Int_t kData = 1,
68 Bool_t bToy = kFALSE) {
6acdbcb2 69 //Get the TList objects (QA, bf, bf shuffled)
6acdbcb2 70 TList *listBF = 0x0;
6acdbcb2 71
72 //Open the file
5365d1d7 73 TFile *f = TFile::Open(filename,"UPDATE");
6acdbcb2 74 if((!f)||(!f->IsOpen())) {
75 Printf("The file %s is not found. Aborting...",filename);
76 return listBF;
77 }
78 //f->ls();
79
80 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
81 if(!dir) {
82 Printf("The TDirectoryFile is not found. Aborting...",filename);
83 return listBF;
84 }
85 //dir->ls();
86
87 TString listBFName;
eb63b883 88 if(kData == 0) {
89 //cout<<"no shuffling - no mixing"<<endl;
3de51a46 90 listBFName = "listBFPsi";
eb63b883 91 }
92 else if(kData == 1) {
93 //cout<<"shuffling - no mixing"<<endl;
3de51a46 94 listBFName = "listBFPsiShuffled";
eb63b883 95 }
96 else if(kData == 2) {
97 //cout<<"no shuffling - mixing"<<endl;
3de51a46 98 listBFName = "listBFPsiMixed";
99 }
100
101 // different list names in case of toy model
102 if(!bToy){
dd46a0c3 103 listBFName += "_";
3de51a46 104 listBFName += centralityArray[gCentrality-1];
105 if(gBit > -1) {
106 listBFName += "_Bit"; listBFName += gBit; }
107 if(gCentralityEstimator) {
108 listBFName += "_"; listBFName += gCentralityEstimator;}
109 }
110 else{
111 listBFName.ReplaceAll("Psi","");
eb63b883 112 }
6acdbcb2 113
5365d1d7 114 // histograms were already retrieved (in first iteration)
115 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
116 //cout<<"second iteration"<<endl;
117 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
6acdbcb2 118 }
6acdbcb2 119
5365d1d7 120 // histograms were not yet retrieved (this is the first iteration)
121 else{
122
123 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
124 cout<<"======================================================="<<endl;
5365d1d7 125 cout<<"List name: "<<listBF->GetName()<<endl;
126 //listBF->ls();
6acdbcb2 127
5365d1d7 128 //Get the histograms
129 TString histoName;
130 if(kData == 0)
27634c13 131 histoName = "fHistP";
5365d1d7 132 else if(kData == 1)
27634c13 133 histoName = "fHistP_shuffle";
5365d1d7 134 else if(kData == 2)
27634c13 135 histoName = "fHistP";
136 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 137 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
138 if(!fHistP) {
139 Printf("fHistP %s not found!!!",histoName.Data());
140 break;
141 }
142 fHistP->FillParent(); fHistP->DeleteContainers();
143
144 if(kData == 0)
27634c13 145 histoName = "fHistN";
5365d1d7 146 if(kData == 1)
27634c13 147 histoName = "fHistN_shuffle";
5365d1d7 148 if(kData == 2)
27634c13 149 histoName = "fHistN";
150 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 151 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
152 if(!fHistN) {
153 Printf("fHistN %s not found!!!",histoName.Data());
154 break;
155 }
156 fHistN->FillParent(); fHistN->DeleteContainers();
157
158 if(kData == 0)
27634c13 159 histoName = "fHistPN";
5365d1d7 160 if(kData == 1)
27634c13 161 histoName = "fHistPN_shuffle";
5365d1d7 162 if(kData == 2)
27634c13 163 histoName = "fHistPN";
164 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 165 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
166 if(!fHistPN) {
167 Printf("fHistPN %s not found!!!",histoName.Data());
168 break;
169 }
170 fHistPN->FillParent(); fHistPN->DeleteContainers();
171
172 if(kData == 0)
27634c13 173 histoName = "fHistNP";
5365d1d7 174 if(kData == 1)
27634c13 175 histoName = "fHistNP_shuffle";
5365d1d7 176 if(kData == 2)
27634c13 177 histoName = "fHistNP";
178 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 179 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
180 if(!fHistNP) {
181 Printf("fHistNP %s not found!!!",histoName.Data());
182 break;
183 }
184 fHistNP->FillParent(); fHistNP->DeleteContainers();
185
186 if(kData == 0)
27634c13 187 histoName = "fHistPP";
5365d1d7 188 if(kData == 1)
27634c13 189 histoName = "fHistPP_shuffle";
5365d1d7 190 if(kData == 2)
27634c13 191 histoName = "fHistPP";
192 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 193 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
194 if(!fHistPP) {
195 Printf("fHistPP %s not found!!!",histoName.Data());
196 break;
197 }
198 fHistPP->FillParent(); fHistPP->DeleteContainers();
199
200 if(kData == 0)
27634c13 201 histoName = "fHistNN";
5365d1d7 202 if(kData == 1)
27634c13 203 histoName = "fHistNN_shuffle";
5365d1d7 204 if(kData == 2)
27634c13 205 histoName = "fHistNN";
206 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 207 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
208 if(!fHistNN) {
209 Printf("fHistNN %s not found!!!",histoName.Data());
210 break;
211 }
212 fHistNN->FillParent(); fHistNN->DeleteContainers();
213
214 dir->cd();
215 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
216
217 }// first iteration
6acdbcb2 218
5365d1d7 219 f->Close();
6acdbcb2 220
221 return listBF;
222}
223
224//______________________________________________________//
eb63b883 225void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
27634c13 226 Int_t gCentrality, const char* gCentralityEstimator,
227 Double_t psiMin, Double_t psiMax,
34616eda 228 Double_t vertexZMin,
229 Double_t vertexZMax,
6acdbcb2 230 Double_t ptTriggerMin, Double_t ptTriggerMax,
f0c5040c 231 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
832b21d8 232 Bool_t kUseVzBinning=kFALSE,
d67eae53 233 Bool_t k2pMethod = kFALSE, TString eventClass) {
6acdbcb2 234 //balance function
235 AliTHn *hP = NULL;
236 AliTHn *hN = NULL;
237 AliTHn *hPN = NULL;
238 AliTHn *hNP = NULL;
239 AliTHn *hPP = NULL;
240 AliTHn *hNN = NULL;
241 //listBF->ls();
242 //Printf("=================");
27634c13 243 TString histoName = "fHistP";
244 if(gCentralityEstimator) histoName += gCentralityEstimator;
245 hP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 246 hP->SetName("gHistP");
27634c13 247 histoName = "fHistN";
248 if(gCentralityEstimator) histoName += gCentralityEstimator;
249 hN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 250 hN->SetName("gHistN");
27634c13 251 histoName = "fHistPN";
252 if(gCentralityEstimator) histoName += gCentralityEstimator;
253 hPN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 254 hPN->SetName("gHistPN");
27634c13 255 histoName = "fHistNP";
256 if(gCentralityEstimator) histoName += gCentralityEstimator;
257 hNP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 258 hNP->SetName("gHistNP");
27634c13 259 histoName = "fHistPP";
260 if(gCentralityEstimator) histoName += gCentralityEstimator;
261 hPP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 262 hPP->SetName("gHistPP");
27634c13 263 histoName = "fHistNN";
264 if(gCentralityEstimator) histoName += gCentralityEstimator;
265 hNN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 266 hNN->SetName("gHistNN");
6acdbcb2 267
268 AliBalancePsi *b = new AliBalancePsi();
d67eae53 269 b->SetEventClass(eventClass);
6acdbcb2 270 b->SetHistNp(hP);
271 b->SetHistNn(hN);
272 b->SetHistNpn(hPN);
273 b->SetHistNnp(hNP);
274 b->SetHistNpp(hPP);
275 b->SetHistNnn(hNN);
832b21d8 276 if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
277
6acdbcb2 278
279 //balance function shuffling
280 AliTHn *hPShuffled = NULL;
281 AliTHn *hNShuffled = NULL;
282 AliTHn *hPNShuffled = NULL;
283 AliTHn *hNPShuffled = NULL;
284 AliTHn *hPPShuffled = NULL;
285 AliTHn *hNNShuffled = NULL;
52daf7b2 286 if(listBFShuffled) {
287 //listBFShuffled->ls();
abf0c200 288 histoName = "fHistP_shuffle";
27634c13 289 if(gCentralityEstimator) histoName += gCentralityEstimator;
290 hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 291 hPShuffled->SetName("gHistPShuffled");
abf0c200 292 histoName = "fHistN_shuffle";
27634c13 293 if(gCentralityEstimator) histoName += gCentralityEstimator;
294 hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 295 hNShuffled->SetName("gHistNShuffled");
abf0c200 296 histoName = "fHistPN_shuffle";
27634c13 297 if(gCentralityEstimator) histoName += gCentralityEstimator;
298 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 299 hPNShuffled->SetName("gHistPNShuffled");
abf0c200 300 histoName = "fHistNP_shuffle";
27634c13 301 if(gCentralityEstimator) histoName += gCentralityEstimator;
302 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 303 hNPShuffled->SetName("gHistNPShuffled");
abf0c200 304 histoName = "fHistPP_shuffle";
27634c13 305 if(gCentralityEstimator) histoName += gCentralityEstimator;
306 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 307 hPPShuffled->SetName("gHistPPShuffled");
abf0c200 308 histoName = "fHistNN_shuffle";
27634c13 309 if(gCentralityEstimator) histoName += gCentralityEstimator;
310 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 311 hNNShuffled->SetName("gHistNNShuffled");
312
313 AliBalancePsi *bShuffled = new AliBalancePsi();
d67eae53 314 bShuffled->SetEventClass(eventClass);
52daf7b2 315 bShuffled->SetHistNp(hPShuffled);
316 bShuffled->SetHistNn(hNShuffled);
317 bShuffled->SetHistNpn(hPNShuffled);
318 bShuffled->SetHistNnp(hNPShuffled);
319 bShuffled->SetHistNpp(hPPShuffled);
320 bShuffled->SetHistNnn(hNNShuffled);
832b21d8 321 if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
322
52daf7b2 323 }
6acdbcb2 324
eb63b883 325 //balance function mixing
326 AliTHn *hPMixed = NULL;
327 AliTHn *hNMixed = NULL;
328 AliTHn *hPNMixed = NULL;
329 AliTHn *hNPMixed = NULL;
330 AliTHn *hPPMixed = NULL;
331 AliTHn *hNNMixed = NULL;
eb63b883 332
52daf7b2 333 if(listBFMixed) {
334 //listBFMixed->ls();
27634c13 335 histoName = "fHistP";
336 if(gCentralityEstimator) histoName += gCentralityEstimator;
337 hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 338 hPMixed->SetName("gHistPMixed");
27634c13 339 histoName = "fHistN";
340 if(gCentralityEstimator) histoName += gCentralityEstimator;
341 hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 342 hNMixed->SetName("gHistNMixed");
27634c13 343 histoName = "fHistPN";
344 if(gCentralityEstimator) histoName += gCentralityEstimator;
345 hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 346 hPNMixed->SetName("gHistPNMixed");
27634c13 347 histoName = "fHistNP";
348 if(gCentralityEstimator) histoName += gCentralityEstimator;
349 hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
350 histoName = "fHistNP";
351 if(gCentralityEstimator) histoName += gCentralityEstimator;
52daf7b2 352 hNPMixed->SetName("gHistNPMixed");
27634c13 353 histoName = "fHistPP";
354 if(gCentralityEstimator) histoName += gCentralityEstimator;
355 hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 356 hPPMixed->SetName("gHistPPMixed");
27634c13 357 histoName = "fHistNN";
358 if(gCentralityEstimator) histoName += gCentralityEstimator;
359 hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 360 hNNMixed->SetName("gHistNNMixed");
361
362 AliBalancePsi *bMixed = new AliBalancePsi();
d67eae53 363 bMixed->SetEventClass(eventClass);
52daf7b2 364 bMixed->SetHistNp(hPMixed);
365 bMixed->SetHistNn(hNMixed);
366 bMixed->SetHistNpn(hPNMixed);
367 bMixed->SetHistNnp(hNPMixed);
368 bMixed->SetHistNpp(hPPMixed);
369 bMixed->SetHistNnn(hNNMixed);
832b21d8 370 if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
371
52daf7b2 372 }
eb63b883 373
6acdbcb2 374 TH2D *gHistBalanceFunction;
eb63b883 375 TH2D *gHistBalanceFunctionSubtracted;
6acdbcb2 376 TH2D *gHistBalanceFunctionShuffled;
eb63b883 377 TH2D *gHistBalanceFunctionMixed;
6acdbcb2 378 TString histoTitle, pngName;
6acdbcb2 379
d67eae53 380 if(eventClass == "Centrality"){
381 histoTitle = "Centrality: ";
382 histoTitle += psiMin;
383 histoTitle += " - ";
384 histoTitle += psiMax;
385 histoTitle += " % ";
52daf7b2 386 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
d67eae53 387 }
388 else if(eventClass == "Multiplicity"){
389 histoTitle = "Multiplicity: ";
390 histoTitle += psiMin;
391 histoTitle += " - ";
392 histoTitle += psiMax;
393 histoTitle += " tracks";
394 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
395 }
396 else{ // "EventPlane" (default)
397 histoTitle = "Centrality: ";
398 histoTitle += centralityArray[gCentrality-1];
399 histoTitle += "%";
400 if((psiMin == -0.5)&&(psiMax == 0.5))
401 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
402 else if((psiMin == 0.5)&&(psiMax == 1.5))
403 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
404 else if((psiMin == 1.5)&&(psiMax == 2.5))
405 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
406 else
407 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
408 }
409
f0c5040c 410 if(k2pMethod)
411 if(bMixed)
34616eda 412 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 413 else{
414 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
415 return;
416 }
417 else
34616eda 418 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 419 gHistBalanceFunction->SetTitle(histoTitle.Data());
420 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
eb63b883 421 gHistBalanceFunction->SetName("gHistBalanceFunction");
422
52daf7b2 423 if(listBFShuffled) {
f0c5040c 424
425 if(k2pMethod)
426 if(bMixed)
34616eda 427 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 428 else{
429 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
430 return;
431 }
432 else
34616eda 433 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 434 gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
435 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
436 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
437 }
eb63b883 438
52daf7b2 439 if(listBFMixed) {
f0c5040c 440 if(k2pMethod)
441 if(bMixed)
34616eda 442 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 443 else{
444 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
445 return;
446 }
447 else
34616eda 448 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 449 gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
450 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
451 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
452
453 gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
454 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
455 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
456 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
457 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
458 }
eb63b883 459
460 //Draw the results
461 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
6acdbcb2 462 c1->SetFillColor(10);
463 c1->SetHighLightColor(10);
464 c1->SetLeftMargin(0.15);
eb63b883 465 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 466 gPad->SetTheta(30); // default is 30
467 gPad->SetPhi(-60); // default is 30
468 gPad->Update();
eb63b883 469 TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
470 c1a->SetFillColor(10);
471 c1a->SetHighLightColor(10);
472 c1a->SetLeftMargin(0.15);
473 gHistBalanceFunction->DrawCopy("colz");
474
52daf7b2 475 if(listBFShuffled) {
476 TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
477 c2->SetFillColor(10);
478 c2->SetHighLightColor(10);
479 c2->SetLeftMargin(0.15);
480 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 481 gPad->SetTheta(30); // default is 30
482 gPad->SetPhi(-60); // default is 30
483 gPad->Update();
52daf7b2 484 TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
485 c2a->SetFillColor(10);
486 c2a->SetHighLightColor(10);
487 c2a->SetLeftMargin(0.15);
488 gHistBalanceFunctionShuffled->DrawCopy("colz");
489 }
eb63b883 490
52daf7b2 491 if(listBFMixed) {
492 TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
493 c3->SetFillColor(10);
494 c3->SetHighLightColor(10);
495 c3->SetLeftMargin(0.15);
496 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 497 gPad->SetTheta(30); // default is 30
498 gPad->SetPhi(-60); // default is 30
499 gPad->Update();
52daf7b2 500 TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
501 c3a->SetFillColor(10);
502 c3a->SetHighLightColor(10);
503 c3a->SetLeftMargin(0.15);
504 gHistBalanceFunctionMixed->DrawCopy("colz");
eb63b883 505
52daf7b2 506 TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
507 c4->SetFillColor(10);
508 c4->SetHighLightColor(10);
509 c4->SetLeftMargin(0.15);
510 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 511 gPad->SetTheta(30); // default is 30
512 gPad->SetPhi(-60); // default is 30
513 gPad->Update();
52daf7b2 514 TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
515 c4a->SetFillColor(10);
516 c4a->SetHighLightColor(10);
517 c4a->SetLeftMargin(0.15);
518 gHistBalanceFunctionSubtracted->DrawCopy("colz");
742af4bd 519
bd36d661 520 fitbalanceFunction(gCentrality, psiMin , psiMax,
742af4bd 521 ptTriggerMin, ptTriggerMax,
522 ptAssociatedMin, ptAssociatedMax,
d67eae53 523 gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
52daf7b2 524 }
eb63b883 525
d67eae53 526 TString newFileName = "balanceFunction2D.";
527 if(eventClass == "Centrality"){
528 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
529 newFileName += ".PsiAll.PttFrom";
530 }
531 else if(eventClass == "Multiplicity"){
532 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
533 newFileName += ".PsiAll.PttFrom";
534 }
535 else{ // "EventPlane" (default)
536 newFileName += "Centrality";
537 newFileName += gCentrality; newFileName += ".Psi";
538 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
539 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
540 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
541 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
542 else newFileName += "All.PttFrom";
543 }
648f1a5a 544 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
545 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
546 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
547 newFileName += Form("%.1f",ptAssociatedMax);
f0c5040c 548 if(k2pMethod) newFileName += "_2pMethod";
a9f20288 549
a2da950e 550 newFileName += "_";
551 newFileName += Form("%.1f",psiMin);
552 newFileName += "-";
553 newFileName += Form("%.1f",psiMax);
554 newFileName += ".root";
eb63b883 555
556 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
557 fOutput->cd();
52daf7b2 558 /*hP->Write(); hN->Write();
559 hPN->Write(); hNP->Write();
560 hPP->Write(); hNN->Write();
561 hPShuffled->Write(); hNShuffled->Write();
562 hPNShuffled->Write(); hNPShuffled->Write();
563 hPPShuffled->Write(); hNNShuffled->Write();
564 hPMixed->Write(); hNMixed->Write();
565 hPNMixed->Write(); hNPMixed->Write();
566 hPPMixed->Write(); hNNMixed->Write();*/
eb63b883 567 gHistBalanceFunction->Write();
52daf7b2 568 if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
569 if(listBFMixed) {
570 gHistBalanceFunctionMixed->Write();
571 gHistBalanceFunctionSubtracted->Write();
572 }
eb63b883 573 fOutput->Close();
6acdbcb2 574}
575
742af4bd 576//____________________________________________________________//
577void fitbalanceFunction(Int_t gCentrality = 1,
578 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
579 Double_t ptTriggerMin = -1.,
580 Double_t ptTriggerMax = -1.,
581 Double_t ptAssociatedMin = -1.,
582 Double_t ptAssociatedMax = -1.,
d67eae53 583 TH2D *gHist,
584 Bool_t k2pMethod = kFALSE,
585 TString eventClass="EventPlane") {
48f521eb 586 //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
587 //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
742af4bd 588 cout<<"FITTING FUNCTION"<<endl;
589
48f521eb 590 TF2 *gFitFunction = new TF2("gFitFunction","[0] + [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2)) + [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))",-1.2,1.2,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
742af4bd 591 gFitFunction->SetName("gFitFunction");
592
742af4bd 593 //Normalization
594 gFitFunction->SetParName(0,"N1");
48f521eb 595 gFitFunction->SetParameter(0,1.0);
596
597 //2D balance function
598 gFitFunction->SetParName(1,"N_{BF}");
599 gFitFunction->SetParameter(1,1.0);
600 gFitFunction->SetParLimits(1, 0., 100.);
601 gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
602 gFitFunction->SetParameter(2,0.6);
603 gFitFunction->SetParLimits(2, 0., 1.);
604 gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
605 gFitFunction->SetParameter(3,0.0);
606 gFitFunction->SetParLimits(3, -0.2, 0.2);
607 gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
608 gFitFunction->SetParameter(4,0.6);
609 gFitFunction->SetParLimits(4, 0., 1.);
610 gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
611 gFitFunction->SetParameter(5,0.0);
612 gFitFunction->SetParLimits(5, -0.2, 0.2);
613
614 //short range structure
615 gFitFunction->SetParName(6,"N_{SR}");
616 gFitFunction->SetParameter(6,5.0);
617 gFitFunction->SetParLimits(6, 0., 100.);
618 gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
619 gFitFunction->SetParameter(7,0.01);
620 gFitFunction->SetParLimits(7, 0.0, 0.1);
621 gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
622 gFitFunction->SetParameter(8,0.0);
623 gFitFunction->SetParLimits(8, -0.01, 0.01);
624 gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
625 gFitFunction->SetParameter(9,0.01);
626 gFitFunction->SetParLimits(9, 0.0, 0.1);
627 gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
628 gFitFunction->SetParameter(10,0.0);
629 gFitFunction->SetParLimits(10, -0.01, 0.01);
742af4bd 630
631
632 //Cloning the histogram
742af4bd 633 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
634 gHistResidual->SetName("gHistResidual");
635 gHistResidual->Sumw2();
636
48f521eb 637 //Fitting the 2D bf
638 for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
639 gHist->Fit("gFitFunction","nm");
640 for(Int_t iParam = 0; iParam < 11; iParam++)
641 gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
742af4bd 642 }
48f521eb 643 cout<<"======================================================"<<endl;
644 cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
645 cout<<"======================================================"<<endl;
646
647 //Getting the residual
648 gHistResidual->Add(gFitFunction,-1);
742af4bd 649
650 //Write to output file
d67eae53 651 TString newFileName = "balanceFunctionFit2D.";
652 if(eventClass == "Centrality"){
653 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
654 newFileName += ".PsiAll.PttFrom";
655 }
656 else if(eventClass == "Multiplicity"){
657 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
658 newFileName += ".PsiAll.PttFrom";
659 }
660 else{ // "EventPlane" (default)
661 newFileName += "Centrality";
662 newFileName += gCentrality; newFileName += ".Psi";
663 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
664 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
665 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
666 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
667 else newFileName += "All.PttFrom";
668 }
742af4bd 669 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
670 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
671 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
672 newFileName += Form("%.1f",ptAssociatedMax);
d67eae53 673 if(k2pMethod) newFileName += "_2pMethod";
a9f20288 674
675 newFileName += "_";
676 newFileName += Form("%.1f",psiMin);
677 newFileName += "-";
678 newFileName += Form("%.1f",psiMax);
742af4bd 679 newFileName += ".root";
a9f20288 680
742af4bd 681 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
682 gHist->Write();
742af4bd 683 gHistResidual->Write();
684 gFitFunction->Write();
685 newFile->Close();
742af4bd 686}
687
eb63b883 688//____________________________________________________________//
648f1a5a 689void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
bd36d661 690 const char* gCentralityEstimator = "V0M",
691 Int_t gBit = 128,
692 const char* gEventPlaneEstimator = "VZERO",
648f1a5a 693 Int_t gCentrality = 1,
694 Bool_t kShowShuffled = kFALSE,
695 Bool_t kShowMixed = kFALSE,
eb63b883 696 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
697 Double_t ptTriggerMin = -1.,
698 Double_t ptTriggerMax = -1.,
699 Double_t ptAssociatedMin = -1.,
bd36d661 700 Double_t ptAssociatedMax = -1.,
701 Bool_t k2pMethod = kTRUE) {
eb63b883 702 //Macro that draws the BF distributions for each centrality bin
703 //for reaction plane dependent analysis
704 //Author: Panos.Christakoglou@nikhef.nl
db7174c0 705 TGaxis::SetMaxDigits(3);
eb63b883 706
707 //Get the input file
bd36d661 708 TString filename = lhcPeriod;
709 filename += "/Centrality"; filename += gCentralityEstimator;
710 filename += "_Bit"; filename += gBit;
711 filename += "_"; filename += gEventPlaneEstimator;
712 filename +="/PttFrom";
648f1a5a 713 filename += Form("%.1f",ptTriggerMin); filename += "To";
714 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
715 filename += Form("%.1f",ptAssociatedMin); filename += "To";
716 filename += Form("%.1f",ptAssociatedMax);
717 filename += "/balanceFunction2D.Centrality";
eb63b883 718 filename += gCentrality; filename += ".Psi";
719 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
720 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
721 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
648f1a5a 722 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
723 else filename += "All.PttFrom";
724 filename += Form("%.1f",ptTriggerMin); filename += "To";
db7174c0 725 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
648f1a5a 726 filename += Form("%.1f",ptAssociatedMin); filename += "To";
bd36d661 727 filename += Form("%.1f",ptAssociatedMax);
728 if(k2pMethod) filename += "_2pMethod";
a9f20288 729
730 filename += "_";
731 filename += Form("%.1f",psiMin);
732 filename += "-";
733 filename += Form("%.1f",psiMax);
bd36d661 734 filename += ".root";
eb63b883 735
736 //Open the file
737 TFile *f = TFile::Open(filename.Data());
738 if((!f)||(!f->IsOpen())) {
739 Printf("The file %s is not found. Aborting...",filename);
740 return listBF;
6acdbcb2 741 }
eb63b883 742 //f->ls();
743
744 //Raw balance function
745 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
746 gHistBalanceFunction->SetStats(kFALSE);
747 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
748 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
749 gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
750 gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
751 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
752 gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
753 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
27634c13 754 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
15dd45a0 755 gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
eb63b883 756
757 //Shuffled balance function
648f1a5a 758 if(kShowShuffled) {
759 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
760 gHistBalanceFunctionShuffled->SetStats(kFALSE);
761 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
762 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
763 gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
764 gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
765 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
766 gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
767 gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
27634c13 768 gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 769 gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
770 }
eb63b883 771
772 //Mixed balance function
648f1a5a 773 if(kShowMixed) {
774 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
775 gHistBalanceFunctionMixed->SetStats(kFALSE);
776 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
777 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
778 gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
779 gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
780 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
781 gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
782 gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
27634c13 783 gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 784 gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
785 }
eb63b883 786
787 //Subtracted balance function
648f1a5a 788 if(kShowMixed) {
789 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
790 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
791 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
792 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
793 gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
794 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
795 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
796 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
797 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
27634c13 798 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 799 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
800 }
801
eb63b883 802 TString pngName;
803
804 TString centralityLatex = "Centrality: ";
805 centralityLatex += centralityArray[gCentrality-1];
806 centralityLatex += "%";
807
808 TString psiLatex;
809 if((psiMin == -0.5)&&(psiMax == 0.5))
810 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
811 else if((psiMin == 0.5)&&(psiMax == 1.5))
812 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
813 else if((psiMin == 1.5)&&(psiMax == 2.5))
814 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
815 else
816 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
817
818 TString pttLatex = Form("%.1f",ptTriggerMin);
819 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
820 pttLatex += " GeV/c";
821
822 TString ptaLatex = Form("%.1f",ptAssociatedMin);
823 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
824 ptaLatex += " GeV/c";
825
826 TLatex *latexInfo1 = new TLatex();
827 latexInfo1->SetNDC();
15dd45a0 828 latexInfo1->SetTextSize(0.045);
eb63b883 829 latexInfo1->SetTextColor(1);
830
831 //Draw the results
832 TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
833 c1->SetFillColor(10); c1->SetHighLightColor(10);
834 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
bd36d661 835 gHistBalanceFunction->SetTitle("");
eb63b883 836 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
837 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
53ff9638 838 //gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
eb63b883 839 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
27634c13 840 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
eb63b883 841 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 842 gPad->SetTheta(30); // default is 30
843 gPad->SetPhi(-60); // default is 30
844 gPad->Update();
eb63b883 845
15dd45a0 846 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
847 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
848 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
849 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
eb63b883 850
bd36d661 851 TString pngName = "BalanceFunction2D.";
852 pngName += "Centrality";
853 pngName += gCentrality;
854 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
855 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
856 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
857 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
858 else pngName += "All.PttFrom";
859 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
860 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
861 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
862 pngName += Form("%.1f",ptAssociatedMax);
863 if(k2pMethod) pngName += "_2pMethod";
864 pngName += ".png";
865
866 c1->SaveAs(pngName.Data());
867
648f1a5a 868 if(kShowShuffled) {
869 TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
870 c2->SetFillColor(10); c2->SetHighLightColor(10);
871 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
872 gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
873 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
874 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
875 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
876 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 877 gPad->SetTheta(30); // default is 30
878 gPad->SetPhi(-60); // default is 30
879 gPad->Update();
880
648f1a5a 881 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
882 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
883 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
884 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
885 }
eb63b883 886
648f1a5a 887 if(kShowMixed) {
888 TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
889 c3->SetFillColor(10); c3->SetHighLightColor(10);
890 c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
891 gHistBalanceFunctionMixed->SetTitle("Mixed events");
892 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
893 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
894 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
895 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 896 gPad->SetTheta(30); // default is 30
897 gPad->SetPhi(-60); // default is 30
898 gPad->Update();
899
648f1a5a 900 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
901 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
902 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
903 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
904 }
eb63b883 905
648f1a5a 906 if(kShowMixed) {
907 TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
908 c4->SetFillColor(10); c4->SetHighLightColor(10);
909 c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
910 gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
911 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
912 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
913 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
914 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 915 gPad->SetTheta(30); // default is 30
916 gPad->SetPhi(-60); // default is 30
917 gPad->Update();
918
648f1a5a 919 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
920 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
921 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
922 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
923 }
6acdbcb2 924}
17c69b4e 925
926//____________________________________________________________//
79dc837f 927void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
928 Bool_t kProjectInEta = kFALSE,
17c69b4e 929 Int_t binMin = 1,
930 Int_t binMax = 80,
931 Int_t gCentrality = 1,
932 Double_t psiMin = -0.5,
933 Double_t psiMax = 3.5,
17c69b4e 934 Double_t ptTriggerMin = -1.,
935 Double_t ptTriggerMax = -1.,
936 Double_t ptAssociatedMin = -1.,
937 Double_t ptAssociatedMax = -1.,
17c69b4e 938 Bool_t k2pMethod = kTRUE,
939 TString eventClass = "Centrality",
1822002d 940 Bool_t bRootMoments = kFALSE,
40a8b898 941 Bool_t kUseZYAM = kFALSE,
7827f01e 942 Bool_t bReduceRangeForMoments = kFALSE,
943 Bool_t bFWHM = kFALSE) {
438dab2c 944 //Macro that draws the balance functions PROJECTIONS
17c69b4e 945 //for each centrality bin for the different pT of trigger and
946 //associated particles
947 TGaxis::SetMaxDigits(3);
948
949 //first we need some libraries
a9f20288 950 gSystem->Load("libTree");
951 gSystem->Load("libGeom");
952 gSystem->Load("libVMC");
953 gSystem->Load("libXMLIO");
954 gSystem->Load("libPhysics");
955
956 gSystem->Load("libSTEERBase");
957 gSystem->Load("libESD");
958 gSystem->Load("libAOD");
959
4070f709 960 gSystem->Load("libANALYSIS");
961 gSystem->Load("libANALYSISalice");
962 gSystem->Load("libEventMixing");
963 gSystem->Load("libCORRFW");
964 gSystem->Load("libPWGTools");
965 gSystem->Load("libPWGCFebye");
17c69b4e 966
438dab2c 967 gStyle->SetOptStat(0);
968
17c69b4e 969 //Get the input file
970 TString filename = "balanceFunction2D.";
971 if(eventClass == "Centrality"){
972 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
973 filename += ".PsiAll.PttFrom";
974 }
975 else if(eventClass == "Multiplicity"){
976 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
977 filename += ".PsiAll.PttFrom";
978 }
979 else{ // "EventPlane" (default)
980 filename += "Centrality";
981 filename += gCentrality; filename += ".Psi";
982 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
983 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
984 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
985 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
986 else filename += "All.PttFrom";
987 }
988 filename += Form("%.1f",ptTriggerMin); filename += "To";
989 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
990 filename += Form("%.1f",ptAssociatedMin); filename += "To";
991 filename += Form("%.1f",ptAssociatedMax);
992 if(k2pMethod) filename += "_2pMethod";
a9f20288 993
a2da950e 994 filename += "_";
995 filename += Form("%.1f",psiMin);
996 filename += "-";
997 filename += Form("%.1f",psiMax);
17c69b4e 998 filename += ".root";
999
1000 //Open the file
79dc837f 1001 TFile *f = 0x0;
1002 if(!gHistBalanceFunction2D) {
1003 TFile::Open(filename.Data());
1004 if((!f)||(!f->IsOpen())) {
1005 Printf("The file %s is not found. Aborting...",filename);
1006 return listBF;
1007 }
1008 //f->ls();
17c69b4e 1009 }
79dc837f 1010
17c69b4e 1011 //Latex
1012 TString centralityLatex = "Centrality: ";
1013 centralityLatex += centralityArray[gCentrality-1];
1014 centralityLatex += "%";
1015
1016 TString psiLatex;
1017 if((psiMin == -0.5)&&(psiMax == 0.5))
1018 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1019 else if((psiMin == 0.5)&&(psiMax == 1.5))
1020 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1021 else if((psiMin == 1.5)&&(psiMax == 2.5))
1022 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1023 else
1024 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1025
1026 TString pttLatex = Form("%.1f",ptTriggerMin);
1027 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1028 pttLatex += " GeV/c";
1029
1030 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1031 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1032 ptaLatex += " GeV/c";
1033
1034 TLatex *latexInfo1 = new TLatex();
1035 latexInfo1->SetNDC();
1036 latexInfo1->SetTextSize(0.045);
1037 latexInfo1->SetTextColor(1);
1038
17c69b4e 1039
1040 //============================================================//
1041 //Get subtracted and mixed balance function
79dc837f 1042 TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1043 TH2D *gHistBalanceFunctionMixed2D = 0x0;
1044 if(!gHistBalanceFunction2D) {
1045 gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1046 gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1047 }
1048 else {
1049 gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1050 gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1051 }
17c69b4e 1052
1053 TH1D *gHistBalanceFunctionSubtracted = NULL;
1054 TH1D *gHistBalanceFunctionMixed = NULL;
1055
a9f20288 1056 TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1057 TH1D *gHistBalanceFunctionMixed_scale = NULL;
1058
17c69b4e 1059 if(kProjectInEta){
86044267 1060 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
a9f20288 1061 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
86044267 1062 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
a9f20288 1063 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
17c69b4e 1064 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1065 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1066 }
1067 else{
86044267 1068 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
a9f20288 1069 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
86044267 1070 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
a9f20288 1071 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
17c69b4e 1072 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1073 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1074 }
1075
1076 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1077 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
17c69b4e 1078
1079 gHistBalanceFunctionMixed->SetMarkerStyle(25);
17c69b4e 1080
1081 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1082 c1->SetFillColor(10);
1083 c1->SetHighLightColor(10);
1084 c1->SetLeftMargin(0.15);
1085 gHistBalanceFunctionSubtracted->DrawCopy("E");
a9f20288 1086 gHistBalanceFunctionMixed->DrawCopy("ESAME");
17c69b4e 1087
1088 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1089 legend->SetTextSize(0.045);
1090 legend->SetTextFont(42);
1091 legend->SetBorderSize(0);
1092 legend->SetFillStyle(0);
1093 legend->SetFillColor(10);
1094 legend->SetMargin(0.25);
1095 legend->SetShadowColor(10);
1096 legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1097 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1098 legend->Draw();
1099
72ed0eda 1100 Double_t sigmaGaus = 0.;
1101 Double_t sigmaGausError = 0.;
17c69b4e 1102
17c69b4e 1103 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1104
1105 if(bRootMoments){
732acb5b 1106
1107 // need to restrict to near side in case of dphi
4331a520 1108 if(!kProjectInEta && !bReduceRangeForMoments) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
732acb5b 1109
53ff9638 1110 if(bReduceRangeForMoments){
40a8b898 1111
62416a3a 1112 Double_t rangeReduced = 0.;
1113
1114 if(kProjectInEta){
1115 if(ptTriggerMax>10){
1116 //for 3<pT,assoc<8<pT,trig<15
1117 rangeReduced = rangePbPbEtaHighMom;
1118 }
1119 else{
1120 // default (for 2<pT,assoc<3<pT,trig<4)
1121 rangeReduced = rangePbPbEtaMedMom;
1122 }
1123 }
1124 else{
1125 if(ptTriggerMax>10){
1126 rangeReduced = rangePbPbPhiHighMom;
1127 }
1128 else{
1129 rangeReduced = rangePbPbPhiMedMom;
1130 }
1131 }
40a8b898
MW		 1132
1133 cout<<"Use reduced range = "<<rangeReduced<<endl;
4331a520 1134 Int_t binLow = gHistBalanceFunctionSubtracted->GetXaxis()->FindBin(-rangeReduced);
1135 Int_t binHigh = gHistBalanceFunctionSubtracted->GetXaxis()->FindBin(rangeReduced);
62416a3a 1136 gHistBalanceFunctionSubtracted->GetXaxis()->SetRange(binLow,binHigh);
40a8b898
MW		 1137 }
1138
17c69b4e 1139 meanLatex = "#mu = ";
1140 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1141 meanLatex += " #pm ";
1142 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1143
1144 rmsLatex = "#sigma = ";
1145 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1146 rmsLatex += " #pm ";
1147 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1148
1149 skewnessLatex = "S = ";
1150 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1151 skewnessLatex += " #pm ";
1152 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1153
1154 kurtosisLatex = "K = ";
1155 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1156 kurtosisLatex += " #pm ";
1157 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
a2da950e 1158
732acb5b 1159 Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1160 Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1161 Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1162 Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
438dab2c 1163
732acb5b 1164
438dab2c 1165 // store in txt files
438dab2c 1166 TString meanFileName = filename;
79dc837f 1167 if(kProjectInEta)
23a0fc15 1168 meanFileName= "deltaEtaProjection_Mean.txt";
1169 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1170 else
23a0fc15 1171 meanFileName = "deltaPhiProjection_Mean.txt";
1172 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1173 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1174 fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1175 fileMean.close();
1176
1177 TString rmsFileName = filename;
79dc837f 1178 if(kProjectInEta)
23a0fc15 1179 rmsFileName = "deltaEtaProjection_Rms.txt";
1180 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1181 else
732acb5b 1182 rmsFileName = "deltaPhiProjection_Rms.txt";
23a0fc15 1183 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1184 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1185 fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
72ed0eda 1186 //fileRms << " " << sigmaGaus << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl; // just for testing
438dab2c 1187 fileRms.close();
1188
1189 TString skewnessFileName = filename;
79dc837f 1190 if(kProjectInEta)
23a0fc15 1191 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1192 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1193 else
23a0fc15 1194 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1195 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1196 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1197 fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1198 fileSkewness.close();
1199
1200 TString kurtosisFileName = filename;
79dc837f 1201 if(kProjectInEta)
23a0fc15 1202 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1203 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1204 else
23a0fc15 1205 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1206 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1207 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1208 fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1209 fileKurtosis.close();
732acb5b 1210
1211 if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
40a8b898 1212 else gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-1.6,1.6);
17c69b4e 1213 }
1214 // calculate the moments by hand
1215 else{
1216
1217 Double_t meanAnalytical, meanAnalyticalError;
1218 Double_t sigmaAnalytical, sigmaAnalyticalError;
1219 Double_t skewnessAnalytical, skewnessAnalyticalError;
1220 Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1822002d 1221
17c69b4e 1222 Int_t gDeltaEtaPhi = 2;
1223 if(kProjectInEta) gDeltaEtaPhi = 1;
1224
1225 AliBalancePsi *bHelper = new AliBalancePsi;
1822002d 1226 bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
17c69b4e 1227
1228 meanLatex = "#mu = ";
1229 meanLatex += Form("%.3f",meanAnalytical);
1230 meanLatex += " #pm ";
1231 meanLatex += Form("%.3f",meanAnalyticalError);
1232
1233 rmsLatex = "#sigma = ";
1234 rmsLatex += Form("%.3f",sigmaAnalytical);
1235 rmsLatex += " #pm ";
1236 rmsLatex += Form("%.3f",sigmaAnalyticalError);
1237
1238 skewnessLatex = "S = ";
1239 skewnessLatex += Form("%.3f",skewnessAnalytical);
1240 skewnessLatex += " #pm ";
1241 skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1242
1243 kurtosisLatex = "K = ";
1244 kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1245 kurtosisLatex += " #pm ";
1246 kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1247 Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1248 Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1249 Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1250 Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
438dab2c 1251
1252 // store in txt files
1253 TString meanFileName = filename;
79dc837f 1254 if(kProjectInEta)
23a0fc15 1255 meanFileName = "deltaEtaProjection_Mean.txt";
1256 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1257 else
23a0fc15 1258 meanFileName = "deltaPhiProjection_Mean.txt";
1259 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1260 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1261 fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1262 fileMean.close();
1263
1264 TString rmsFileName = filename;
79dc837f 1265 if(kProjectInEta)
23a0fc15 1266 rmsFileName = "deltaEtaProjection_Rms.txt";
1267 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1268 else
23a0fc15 1269 rmsFileName = "deltaPhiProjection_Rms.txt";
1270//rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1271 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1272 fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1273 fileRms.close();
1274
1275 TString skewnessFileName = filename;
79dc837f 1276 if(kProjectInEta)
23a0fc15 1277 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1278 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1279 else
23a0fc15 1280 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1281 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1282 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1283 fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1284 fileSkewness.close();
1285
1286 TString kurtosisFileName = filename;
79dc837f 1287 if(kProjectInEta)
23a0fc15 1288 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1289 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1290 else
23a0fc15 1291 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1292 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1293 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1294 fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1295 fileKurtosis.close();
17c69b4e 1296 }
1297
3de51a46 1298 // Weighted mean as calculated for 1D analysis
1299 Double_t weightedMean, weightedMeanError;
40a8b898
MW		 1300 if(bReduceRangeForMoments){
1301 GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,rangeReduced);
1302 }
1303 else{
1304 GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,-1.);
1305 }
3de51a46 1306 Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1307
1308 // store in txt files
1309 TString weightedMeanFileName = filename;
1310 if(kProjectInEta)
1311 weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1312 //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1313 else
1314 weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1315 //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1316 ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1317 fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1318 fileWeightedMean.close();
438dab2c 1319
17c69b4e 1320 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1321 c2->SetFillColor(10);
1322 c2->SetHighLightColor(10);
1323 c2->SetLeftMargin(0.15);
1324 gHistBalanceFunctionSubtracted->DrawCopy("E");
7827f01e 1325
1326 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1327 Double_t fwhm_spline = 0.;
1328 Double_t fwhmError = 0.;
1329 if (bFWHM){
1330 AliBalancePsi *bHelper_1 = new AliBalancePsi;
1331 bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError);
1332 Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError);
1333
1334 //store and in txt files FWHM
1335 TString sigmaFileName = filename;
1336 if(kProjectInEta) {
1337 sigmaFileName = "deltaEtaProjection_FWHM.txt";
1338 }
1339 else{
1340 sigmaFileName = "deltaPhiProjection_FWHM.txt";
1341 }
1342 ofstream fileSigma(sigmaFileName.Data(),ios::app);
1343 fileSigma << " " << fwhm_spline << " " <<fwhmError<<endl;
1344 fileSigma.close();
1345
1346 gHistBalanceFunctionSubtracted->SetLineColor(2);
1347 gHistBalanceFunctionSubtracted->SetMarkerColor(2);
1348 gHistBalanceFunctionSubtracted->Draw("SAME");
1349 }
1350 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
438dab2c 1351
17c69b4e 1352 TLatex *latex = new TLatex();
1353 latex->SetNDC();
1354 latex->SetTextSize(0.035);
1355 latex->SetTextColor(1);
1356 latex->DrawLatex(0.64,0.85,meanLatex.Data());
1357 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1358 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1359 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1360
438dab2c 1361 TString pngName = filename;
1362 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1363 else pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1364
1365 c2->SaveAs(pngName.Data());
1366
17c69b4e 1367 TString outFileName = filename;
1368 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1369 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1370 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1371 gHistBalanceFunctionSubtracted->Write();
1372 gHistBalanceFunctionMixed->Write();
1373 fProjection->Close();
1374}
79dc837f 1375
1376//____________________________________________________________//
1377void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1378 Int_t gTrainNumber = 64,
1379 const char* gCentralityEstimator = "V0M",
1380 Int_t gBit = 128,
1381 const char* gEventPlaneEstimator = "VZERO",
1382 Int_t gCentrality = 1,
1383 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1384 Double_t vertexZMin = -10.,
1385 Double_t vertexZMax = 10.,
1386 Double_t ptTriggerMin = -1.,
1387 Double_t ptTriggerMax = -1.,
1388 Double_t ptAssociatedMin = -1.,
86044267 1389 Double_t ptAssociatedMax = -1.,
732acb5b 1390 TString eventClass = "Multiplicity",
a91351d9 1391 Bool_t bRootMoments = kFALSE,
40a8b898 1392 Bool_t bFullPhiForEtaProjection = kFALSE,
7827f01e 1393 Bool_t bReduceRangeForMoments = kFALSE,
1394 Bool_t bFWHM = kFALSE) {
79dc837f 1395 //Macro that draws the BF distributions for each centrality bin
1396 //for reaction plane dependent analysis
1397 //Author: Panos.Christakoglou@nikhef.nl
1398 TGaxis::SetMaxDigits(3);
1399 gStyle->SetPalette(55,0);
1400
1401 //Get the input file
1402 TString filename = lhcPeriod;
86044267 1403 if(lhcPeriod != ""){
1404 //filename += "/Train"; filename += gTrainNumber;
1405 filename +="/PttFrom";
1406 filename += Form("%.1f",ptTriggerMin); filename += "To";
1407 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1408 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1409 filename += Form("%.1f",ptAssociatedMax);
1410 filename += "/correlationFunction.";
1411 }
1412 else{
1413 filename += "correlationFunction.";
1414 }
1415 if(eventClass == "Centrality"){
1416 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1417 filename += ".PsiAll.PttFrom";
1418 }
1419 else if(eventClass == "Multiplicity"){
1420 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1421 filename += ".PsiAll.PttFrom";
1422 }
1423 else{ // "EventPlane" (default)
1424 filename += "Centrality";
1425 filename += gCentrality; filename += ".Psi";
1426 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1427 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1428 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1429 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1430 else filename += "All.PttFrom";
1431 }
79dc837f 1432 filename += Form("%.1f",ptTriggerMin); filename += "To";
1433 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1434 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1435 filename += Form("%.1f",ptAssociatedMax);
1436 filename += "_";
1437 filename += Form("%.1f",psiMin);
1438 filename += "-";
1439 filename += Form("%.1f",psiMax);
1440 filename += ".root";
1441
1442 //Open the file
1443 TFile *f = TFile::Open(filename.Data());
1444 if((!f)||(!f->IsOpen())) {
1445 Printf("The file %s is not found. Aborting...",filename);
1446 return listBF;
1447 }
1448 //f->ls();
1449
1450 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1451 if(!gHistPN) return;
1452 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1453 if(!gHistNP) return;
1454 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1455 if(!gHistPP) return;
1456 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1457 if(!gHistNN) return;
1458
53ff9638 1459 // in order to get unzoomed (in older versions used smaller user ranger)
1460 Int_t binMinX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmin()+0.00001);
1461 Int_t binMaxX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmax()-0.00001);
1462 Int_t binMinY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmin()+0.00001);
1463 Int_t binMaxY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmax()-0.00001);
1464 gHistPN->GetXaxis()->SetRange(binMinX,binMaxX); gHistPN->GetYaxis()->SetRange(binMinY,binMaxY);
1465 gHistNP->GetXaxis()->SetRange(binMinX,binMaxX); gHistNP->GetYaxis()->SetRange(binMinY,binMaxY);
1466 gHistPP->GetXaxis()->SetRange(binMinX,binMaxX); gHistPP->GetYaxis()->SetRange(binMinY,binMaxY);
1467 gHistNN->GetXaxis()->SetRange(binMinX,binMaxX); gHistNN->GetYaxis()->SetRange(binMinY,binMaxY);
1468
79dc837f 1469 gHistPN->Sumw2();
1470 gHistPP->Sumw2();
1471 gHistPN->Add(gHistPP,-1);
1472 gHistNP->Sumw2();
1473 gHistNN->Sumw2();
1474 gHistNP->Add(gHistNN,-1);
1475 gHistPN->Add(gHistNP);
3aa5487c 1476 //gHistPN->Scale(0.5);//not needed anymore, since pT,assoc < pT,trig in all pT bins
79dc837f 1477 TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1478 gHistBalanceFunction2D->SetStats(kFALSE);
1479 gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1480 gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1481 gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1482
1483 //Draw the results
1484 TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1485 c0->SetFillColor(10); c0->SetHighLightColor(10);
1486 c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1487 gHistBalanceFunction2D->SetTitle("");
1488 gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1489 gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1490 gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1491 gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
53ff9638 1492 //gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
79dc837f 1493 gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1494 gHistBalanceFunction2D->DrawCopy("lego2");
1495 gPad->SetTheta(30); // default is 30
1496 gPad->SetPhi(-60); // default is 30
1497 gPad->Update();
1498
1499 TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1500
1501 TString pttLatex = Form("%.1f",ptTriggerMin);
1502 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1503 pttLatex += " GeV/c";
1504
1505 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1506 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1507 ptaLatex += " GeV/c";
1508
1509 TLatex *latexInfo1 = new TLatex();
1510 latexInfo1->SetNDC();
1511 latexInfo1->SetTextSize(0.045);
1512 latexInfo1->SetTextColor(1);
1513 latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1514 latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1515 latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1516
1517 TString pngName = "BalanceFunction2D.";
86044267 1518 pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
79dc837f 1519 pngName += ".PttFrom";
1520 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1521 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1522 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1523 pngName += Form("%.1f",ptAssociatedMax);
1524 pngName += ".png";
1525
1526 c0->SaveAs(pngName.Data());
1527
a91351d9
MW		 1528 // do the full range for the projection on eta (for cross checking with published results)
1529 if(bFullPhiForEtaProjection){
1530
1531 drawProjections(gHistBalanceFunction2D,
1532 kTRUE,
1533 1,72,
1534 gCentrality,
1535 psiMin,psiMax,
1536 ptTriggerMin,ptTriggerMax,
1537 ptAssociatedMin,ptAssociatedMax,
1538 kTRUE,
1539 eventClass.Data(),
40a8b898
MW		 1540 bRootMoments,
1541 kFALSE,
7827f01e 1542 bReduceRangeForMoments,
1543 bFWHM);
a91351d9
MW		 1544 }
1545 else{
7827f01e 1546 drawProjections(gHistBalanceFunction2D,
a91351d9
MW		 1547 kTRUE,
1548 1,36,
1549 gCentrality,
1550 psiMin,psiMax,
1551 ptTriggerMin,ptTriggerMax,
1552 ptAssociatedMin,ptAssociatedMax,
1553 kTRUE,
1554 eventClass.Data(),
7827f01e 1555 bRootMoments,
1556 kFALSE,
1557 bReduceRangeForMoments,
1558 bFWHM);
a91351d9 1559 }
79dc837f 1560
1561 drawProjections(gHistBalanceFunction2D,
86044267 1562 kFALSE,
1563 1,80,
1564 gCentrality,
1565 psiMin,psiMax,
1566 ptTriggerMin,ptTriggerMax,
1567 ptAssociatedMin,ptAssociatedMax,
1568 kTRUE,
1569 eventClass.Data(),
40a8b898
MW		 1570 bRootMoments,
1571 kFALSE,
7827f01e 1572 bReduceRangeForMoments,
1573 bFWHM);
3de51a46 1574
1575 TString outFileName = filename;
1576 outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1577 gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1578 TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1579 gHistBalanceFunction2D->Write();
1580 fOut->Close();
1581
1582}
1583
6f14379e 1584
1585//____________________________________________________________//
1586void mergeBFPsi2D(TString momDirectory = "./",
1587 TString directory1 = "",
1588 TString directory2 = "",
1589 TString directory3 = "",
1590 TString directory4 = "",
1591 Int_t gCentrality = 1,
1592 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1593 TString eventClass = "Centrality",
1594 Double_t ptTriggerMin = -1.,
1595 Double_t ptTriggerMax = -1.,
1596 Double_t ptAssociatedMin = -1.,
40a8b898 1597 Double_t ptAssociatedMax = -1.,
864839b8 1598 Bool_t bReduceRangeForMoments = kFALSE,
1599 Bool_t bFWHM = kFALSE
6f14379e 1600) {
1601 //Macro that draws the BF distributions for each centrality bin
1602 //for reaction plane dependent analysis
1603 //Author: Panos.Christakoglou@nikhef.nl
1604 TGaxis::SetMaxDigits(3);
1605 gStyle->SetPalette(55,0);
1606
1607 const Int_t nMaxDirectories = 4; // maximum number of directories to merge (set to 4 for now)
1608 TString sDirectory[nMaxDirectories];
1609 Int_t nDirectories = nMaxDirectories;
1610
1611 TString sFileName[nMaxDirectories];
1612 TFile *fBF[nMaxDirectories];
1613 TH2D *hBF[nMaxDirectories];
1614 Double_t entries[nMaxDirectories];
1615
1616 TFile *fOut;
1617 TH2D *hBFOut;
1618 Double_t entriesOut = 0.;
1619
1620 // find out how many directories we have to merge
1621 if(directory1==""){
1622 nDirectories=0;
1623 }
1624 else if(directory2==""){
1625 nDirectories=1;
1626 sDirectory[0]=directory1;
1627 }
1628 else if(directory3==""){
1629 nDirectories=2;
1630 sDirectory[0]=directory1;
1631 sDirectory[1]=directory2;
1632 }
1633 else if(directory4==""){
1634 nDirectories=3;
1635 sDirectory[0]=directory1;
1636 sDirectory[1]=directory2;
1637 sDirectory[2]=directory3;
1638 }
1639 else{
1640 nDirectories=nMaxDirectories;
1641 sDirectory[0]=directory1;
1642 sDirectory[1]=directory2;
1643 sDirectory[2]=directory3;
1644 sDirectory[3]=directory4;
1645 }
1646
1647 for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1648
1649 //Get the input file
1650 sFileName[iDir] = sDirectory[iDir];
1651 sFileName[iDir] += "/Centrality"; sFileName[iDir] += gCentrality;
1652 sFileName[iDir] += "/";
1653sFileName[iDir] += momDirectory;
1654 sFileName[iDir] += "/balanceFunction2D.";
1655
1656 if(eventClass == "Centrality"){
1657 sFileName[iDir] += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1658 sFileName[iDir] += ".PsiAll.PttFrom";
1659 }
1660 else if(eventClass == "Multiplicity"){
1661 sFileName[iDir] += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1662 sFileName[iDir] += ".PsiAll.PttFrom";
1663 }
1664 else{ // "EventPlane" (default)
1665 sFileName[iDir] += "Centrality";
1666 sFileName[iDir] += gCentrality; sFileName[iDir] += ".Psi";
1667 if((psiMin == -0.5)&&(psiMax == 0.5)) sFileName[iDir] += "InPlane.Ptt";
1668 else if((psiMin == 0.5)&&(psiMax == 1.5)) sFileName[iDir] += "Intermediate.Ptt";
1669 else if((psiMin == 1.5)&&(psiMax == 2.5)) sFileName[iDir] += "OutOfPlane.Ptt";
1670 else if((psiMin == 2.5)&&(psiMax == 3.5)) sFileName[iDir] += "Rest.PttFrom";
1671 else sFileName[iDir] += "All.PttFrom";
1672 }
1673 sFileName[iDir] += Form("%.1f",ptTriggerMin); sFileName[iDir] += "To";
1674 sFileName[iDir] += Form("%.1f",ptTriggerMax); sFileName[iDir] += "PtaFrom";
1675 sFileName[iDir] += Form("%.1f",ptAssociatedMin); sFileName[iDir] += "To";
1676 sFileName[iDir] += Form("%.1f",ptAssociatedMax);
1677 sFileName[iDir] += "_";
1678 sFileName[iDir] += Form("%.1f",psiMin);
1679 sFileName[iDir] += "-";
1680 sFileName[iDir] += Form("%.1f",psiMax);
1681 sFileName[iDir] += ".root";
1682
1683 //Open the file
1684 fBF[iDir] = TFile::Open(sFileName[iDir].Data());
1685 if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
6952fc3f 1686 Printf("The file %s is not found. Not used...",sFileName[iDir].Data());
6f14379e 1687 continue;
1688 }
1689 //fBF[iDir]->ls();
1690
1691 hBF[iDir] = (TH2D*)fBF[iDir]->Get("gHistBalanceFunctionSubtracted");
1692 if(!hBF[iDir]) continue;
6952fc3f 1693 entries[iDir] = (Double_t) hBF[iDir]->GetEffectiveEntries(); // weighted histograms
6f14379e 1694 entriesOut += entries[iDir];
1695 cout<<" BF histogram "<<iDir<<" has "<<entries[iDir] <<" entries."<<endl;
6952fc3f 1696 }
6f14379e 1697
6952fc3f 1698 // second loop for scaling and adding (for average)
1699 for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1700
1701 if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
1702 Printf("The file %s is not found. Not used...",sFileName[iDir].Data());
1703 continue;
1704 }
1705
1706 hBF[iDir]->Scale(entries[iDir]/entriesOut);
1707 if(!hBFOut) hBFOut = (TH2D*)hBF[iDir]->Clone("gHistBalanceFunctionSubtractedOut");
1708 else hBFOut->Add(hBF[iDir]);
6f14379e 1709 }
1710
40a8b898 1711 drawProjections(hBFOut,
7827f01e 1712 kTRUE,
1713 1,36,
1714 gCentrality,
1715 psiMin,psiMax,
1716 ptTriggerMin,ptTriggerMax,
1717 ptAssociatedMin,ptAssociatedMax,
1718 kTRUE,
1719 eventClass,
1720 kTRUE,
1721 kFALSE,
1722 bReduceRangeForMoments,
1723 bFWHM);
6f14379e 1724
1725 drawProjections(hBFOut,
1726 kFALSE,
1727 1,80,
1728 gCentrality,
1729 psiMin,psiMax,
1730 ptTriggerMin,ptTriggerMax,
1731 ptAssociatedMin,ptAssociatedMax,
1732 kTRUE,
1733 eventClass.Data(),
40a8b898
MW		 1734 kTRUE,
1735 kFALSE,
7827f01e 1736 bReduceRangeForMoments,
1737 bFWHM);
6f14379e 1738
1739 // output
1740 TString outFileName = "balanceFunction2D.";
1741
1742 if(eventClass == "Centrality"){
1743 outFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1744 outFileName += ".PsiAll.PttFrom";
1745 }
1746 else if(eventClass == "Multiplicity"){
1747 outFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1748 outFileName += ".PsiAll.PttFrom";
1749 }
1750 else{ // "EventPlane" (default)
1751 outFileName += "Centrality";
1752 outFileName += gCentrality; outFileName += ".Psi";
1753 if((psiMin == -0.5)&&(psiMax == 0.5)) outFileName += "InPlane.Ptt";
1754 else if((psiMin == 0.5)&&(psiMax == 1.5)) outFileName += "Intermediate.Ptt";
1755 else if((psiMin == 1.5)&&(psiMax == 2.5)) outFileName += "OutOfPlane.Ptt";
1756 else if((psiMin == 2.5)&&(psiMax == 3.5)) outFileName += "Rest.PttFrom";
1757 else outFileName += "All.PttFrom";
1758 }
1759 outFileName += Form("%.1f",ptTriggerMin); outFileName += "To";
1760 outFileName += Form("%.1f",ptTriggerMax); outFileName += "PtaFrom";
1761 outFileName += Form("%.1f",ptAssociatedMin); outFileName += "To";
1762 outFileName += Form("%.1f",ptAssociatedMax);
1763 outFileName += "_";
1764 outFileName += Form("%.1f",psiMin);
1765 outFileName += "-";
1766 outFileName += Form("%.1f",psiMax);
1767 outFileName += ".root";
1768
1769 hBFOut->SetName("gHistBalanceFunctionSubtracted");
1770 fBFOut = TFile::Open(outFileName.Data(),"recreate");
1771 hBFOut->Write();
1772 fBFOut->Close();
1773
1774}
1775
3de51a46 1776//____________________________________________________________________//
40a8b898 1777void GetWeightedMean1D(TH1D *gHistBalance, Bool_t kProjectInEta = kTRUE, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME,Double_t rangeReduced = -1.) {
3de51a46 1778
1779 //Prints the calculated width of the BF and its error
1780 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1781 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1782 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1783 Double_t deltaBalP2 = 0.0, integral = 0.0;
1784 Double_t deltaErrorNew = 0.0;
1785
1786 //Retrieve this variables from Histogram
1787 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1788 if(fStopBin > -1) fNumberOfBins = fStopBin;
1789 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
1790 Double_t currentBinCenter = 0.;
1791
1792 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1793
1794 // in order to recover the |abs| in the 1D analysis
1795 currentBinCenter = gHistBalance->GetBinCenter(i);
1796 if(kProjectInEta){
1797 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1798 }
1799 else{
1800 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1801 if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1802 }
1803
40a8b898
MW		 1804 if(rangeReduced > 0 && currentBinCenter > rangeReduced )
1805 continue;
1806
3de51a46 1807 gSumXi += currentBinCenter;
1808 gSumBi += gHistBalance->GetBinContent(i);
1809 gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1810 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1811 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1812 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
1813 gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1814
1815 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1816 integral += fP2Step*gHistBalance->GetBinContent(i);
1817 }
1818 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1819 deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1820
1821 Double_t integralError = TMath::Sqrt(deltaBalP2);
1822
1823 Double_t delta = gSumBiXi / gSumBi;
1824 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1825
1826 WM = delta;
1827 WME = deltaError;
79dc837f 1828}
3de51a46 1829