]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
variable normalization range of the event mixing histograms for correlation functions
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawCorrelationFunctionPsi.C
CommitLineData
437a453a 1const Int_t numberOfCentralityBins = 13;
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"};
a38fd7f3 3
4const Int_t gRebin = 1;
07d0a35c 5
2b4abbc5 6void drawCorrelationFunctionPsi(const char* filename = "AnalysisResultsPsi_train_06feb.root",
52daf7b2 7 Int_t gCentrality = 1,
5de9ad1a 8 Int_t gBit = -1,
9 const char* gCentralityEstimator = 0x0,
52daf7b2 10 Bool_t kShowShuffled = kFALSE,
11 Bool_t kShowMixed = kTRUE,
6acdbcb2 12 Double_t psiMin = -0.5,
5de9ad1a 13 Double_t psiMax = 3.5,
20006629 14 Double_t vertexZMin = -10.,
15 Double_t vertexZMax = 10.,
6acdbcb2 16 Double_t ptTriggerMin = -1.,
17 Double_t ptTriggerMax = -1.,
18 Double_t ptAssociatedMin = -1.,
742af4bd 19 Double_t ptAssociatedMax = -1.,
2b4abbc5 20 Bool_t normToTrig = kTRUE,
3aa5487c 21 Double_t normalizationRangePhi = TMath::Pi()/6.,
20006629 22 Bool_t kUseVzBinning = kTRUE,
742af4bd 23 Int_t rebinEta = 1,
32e94079 24 Int_t rebinPhi = 1,
dd46a0c3 25 TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
633996e7 26 Bool_t bToy = kFALSE,
27 TString listNameAdd = "")
32e94079 28{
a38fd7f3 29 //Macro that draws the correlation functions from the balance function
30 //analysis vs the reaction plane
31 //Author: Panos.Christakoglou@nikhef.nl
2b4abbc5 32 //gROOT->LoadMacro("~/SetPlotStyle.C");
33 //SetPlotStyle();
52daf7b2 34 gStyle->SetPalette(1,0);
35
a38fd7f3 36 //Load the PWG2ebye library
2b4abbc5 37 gSystem->Load("libCore.so");
38 gSystem->Load("libGeom.so");
39 gSystem->Load("libVMC.so");
40 gSystem->Load("libPhysics.so");
41 gSystem->Load("libTree.so");
42 gSystem->Load("libSTEERBase.so");
43 gSystem->Load("libESD.so");
44 gSystem->Load("libAOD.so");
45
a38fd7f3 46 gSystem->Load("libANALYSIS.so");
47 gSystem->Load("libANALYSISalice.so");
48 gSystem->Load("libEventMixing.so");
49 gSystem->Load("libCORRFW.so");
50 gSystem->Load("libPWGTools.so");
51 gSystem->Load("libPWGCFebye.so");
52
53 //Prepare the objects and return them
3ce45c12 54 TList *listQA = NULL;
633996e7 55 TList *list = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy,listNameAdd);
52daf7b2 56 TList *listShuffled = NULL;
633996e7 57 if(kShowShuffled) listShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy,listNameAdd);
52daf7b2 58 TList *listMixed = NULL;
633996e7 59 if(kShowMixed) listMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy,listNameAdd);
52daf7b2 60
94415ad1 61 //get the QA histograms (for convolution)
3ce45c12 62
94415ad1 63 //Open the file again
64 TFile *f = TFile::Open(filename,"UPDATE");
65 if((!f)||(!f->IsOpen())) {
66 Printf("The file %s is not found.",filename);
67 }
a9737921 68
94415ad1 69
70 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
71 if(!dir) {
72 Printf("The TDirectoryFile is not found.",filename);
3ce45c12 73 }
94415ad1 74
75 TString listQAName = "listQAPsi_";
76
77 listQAName += centralityArray[gCentrality-1];
78 if(gBit > -1) {
79 listQAName += "_Bit"; listQAName += gBit; }
80 if(gCentralityEstimator) {
81 listQAName += "_"; listQAName += gCentralityEstimator;}
82 listQAName += listNameAdd;
83
84
85 listQA = (TList*)dir->Get(Form("%s",listQAName.Data()));
86 if(!listQA) {
87 Printf("TList %s not found!!!",listQAName.Data());
88 }
89
90
a38fd7f3 91 if(!list) {
92 Printf("The TList object was not created");
93 return;
94 }
95 else
3ce45c12 96 draw(list,listShuffled,listMixed,listQA,
20006629 97 gCentralityEstimator,gCentrality,psiMin,psiMax,vertexZMin,vertexZMax,
3aa5487c 98 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig, normalizationRangePhi,
99 kUseVzBinning,rebinEta,rebinPhi,eventClass);
a38fd7f3 100}
101
102//______________________________________________________//
52daf7b2 103TList *GetListOfObjects(const char* filename,
104 Int_t gCentrality,
5de9ad1a 105 Int_t gBit,
106 const char *gCentralityEstimator,
dd46a0c3 107 Int_t kData = 1,
633996e7 108 Bool_t bToy = kFALSE,
109 TString listNameAdd = "") {
a38fd7f3 110 //Get the TList objects (QA, bf, bf shuffled)
a38fd7f3 111 TList *listBF = 0x0;
a38fd7f3 112
113 //Open the file
ff0805a7 114 TFile *f = TFile::Open(filename,"UPDATE");
a38fd7f3 115 if((!f)||(!f->IsOpen())) {
116 Printf("The file %s is not found. Aborting...",filename);
117 return listBF;
118 }
119 //f->ls();
120
6acdbcb2 121 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
a38fd7f3 122 if(!dir) {
123 Printf("The TDirectoryFile is not found. Aborting...",filename);
124 return listBF;
125 }
126 //dir->ls();
2b4abbc5 127
52daf7b2 128 TString listBFName;
129 if(kData == 0) {
130 //cout<<"no shuffling - no mixing"<<endl;
dd46a0c3 131 listBFName = "listBFPsi";
52daf7b2 132 }
133 else if(kData == 1) {
134 //cout<<"shuffling - no mixing"<<endl;
dd46a0c3 135 listBFName = "listBFPsiShuffled";
52daf7b2 136 }
137 else if(kData == 2) {
138 //cout<<"no shuffling - mixing"<<endl;
dd46a0c3 139 listBFName = "listBFPsiMixed";
140 }
141
142 // different list names in case of toy model
143 if(!bToy){
144 listBFName += "_";
145 listBFName += centralityArray[gCentrality-1];
146 if(gBit > -1) {
147 listBFName += "_Bit"; listBFName += gBit; }
148 if(gCentralityEstimator) {
149 listBFName += "_"; listBFName += gCentralityEstimator;}
150 }
151 else{
152 listBFName.ReplaceAll("Psi","");
52daf7b2 153 }
633996e7 154 listBFName += listNameAdd;
a38fd7f3 155
5365d1d7 156 // histograms were already retrieved (in first iteration)
157 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
158 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
a38fd7f3 159 }
a38fd7f3 160
5365d1d7 161 // histograms were not yet retrieved (this is the first iteration)
162 else{
163
164 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
633996e7 165 if(!listBF) dir->ls();
5365d1d7 166 cout<<"======================================================="<<endl;
32e94079 167 cout<<"List name (control): "<<listBFName.Data()<<endl;
17d10bdd 168 cout<<"List name: "<<listBF->GetName()<<endl;
5365d1d7 169 //listBF->ls();
a38fd7f3 170
5365d1d7 171 //Get the histograms
172 TString histoName;
17d10bdd 173 if(kData == 0)
174 histoName = "fHistP";
5365d1d7 175 else if(kData == 1)
17d10bdd 176 histoName = "fHistP_shuffle";
5365d1d7 177 else if(kData == 2)
17d10bdd 178 histoName = "fHistP";
179 if(gCentralityEstimator)
180 histoName += gCentralityEstimator;
5365d1d7 181 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
182 if(!fHistP) {
183 Printf("fHistP %s not found!!!",histoName.Data());
184 break;
185 }
186 fHistP->FillParent(); fHistP->DeleteContainers();
187
188 if(kData == 0)
17d10bdd 189 histoName = "fHistN";
5365d1d7 190 if(kData == 1)
17d10bdd 191 histoName = "fHistN_shuffle";
5365d1d7 192 if(kData == 2)
17d10bdd 193 histoName = "fHistN";
194 if(gCentralityEstimator)
195 histoName += gCentralityEstimator;
5365d1d7 196 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
197 if(!fHistN) {
198 Printf("fHistN %s not found!!!",histoName.Data());
199 break;
200 }
201 fHistN->FillParent(); fHistN->DeleteContainers();
202
203 if(kData == 0)
17d10bdd 204 histoName = "fHistPN";
5365d1d7 205 if(kData == 1)
17d10bdd 206 histoName = "fHistPN_shuffle";
5365d1d7 207 if(kData == 2)
17d10bdd 208 histoName = "fHistPN";
209 if(gCentralityEstimator)
210 histoName += gCentralityEstimator;
5365d1d7 211 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
212 if(!fHistPN) {
213 Printf("fHistPN %s not found!!!",histoName.Data());
214 break;
215 }
216 fHistPN->FillParent(); fHistPN->DeleteContainers();
217
218 if(kData == 0)
17d10bdd 219 histoName = "fHistNP";
5365d1d7 220 if(kData == 1)
17d10bdd 221 histoName = "fHistNP_shuffle";
5365d1d7 222 if(kData == 2)
17d10bdd 223 histoName = "fHistNP";
224 if(gCentralityEstimator)
225 histoName += gCentralityEstimator;
5365d1d7 226 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
227 if(!fHistNP) {
228 Printf("fHistNP %s not found!!!",histoName.Data());
229 break;
230 }
231 fHistNP->FillParent(); fHistNP->DeleteContainers();
232
233 if(kData == 0)
17d10bdd 234 histoName = "fHistPP";
5365d1d7 235 if(kData == 1)
17d10bdd 236 histoName = "fHistPP_shuffle";
5365d1d7 237 if(kData == 2)
17d10bdd 238 histoName = "fHistPP";
239 if(gCentralityEstimator)
240 histoName += gCentralityEstimator;
5365d1d7 241 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
242 if(!fHistPP) {
243 Printf("fHistPP %s not found!!!",histoName.Data());
244 break;
245 }
246 fHistPP->FillParent(); fHistPP->DeleteContainers();
247
248 if(kData == 0)
17d10bdd 249 histoName = "fHistNN";
5365d1d7 250 if(kData == 1)
17d10bdd 251 histoName = "fHistNN_shuffle";
5365d1d7 252 if(kData == 2)
17d10bdd 253 histoName = "fHistNN";
254 if(gCentralityEstimator)
255 histoName += gCentralityEstimator;
5365d1d7 256 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
257 if(!fHistNN) {
258 Printf("fHistNN %s not found!!!",histoName.Data());
259 break;
260 }
261 fHistNN->FillParent(); fHistNN->DeleteContainers();
17d10bdd 262
5365d1d7 263 dir->cd();
264 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
265
266 }// first iteration
a38fd7f3 267
ff0805a7 268 f->Close();
5365d1d7 269
a38fd7f3 270 return listBF;
271}
272
273//______________________________________________________//
20006629 274void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
275 TList *listQA,
17d10bdd 276 const char *gCentralityEstimator,
20006629 277 Int_t gCentrality, Double_t psiMin, Double_t psiMax,
278 Double_t vertexZMin,Double_t vertexZMax,
6acdbcb2 279 Double_t ptTriggerMin, Double_t ptTriggerMax,
742af4bd 280 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
3aa5487c 281 Bool_t normToTrig,
282 Double_t normalizationRangePhi,
20006629 283 Bool_t kUseVzBinning,
284 Int_t rebinEta, Int_t rebinPhi,TString eventClass) {
a38fd7f3 285 //Draws the correlation functions for every centrality bin
52daf7b2 286 //(+-), (-+), (++), (--)
a38fd7f3 287 AliTHn *hP = NULL;
288 AliTHn *hN = NULL;
289 AliTHn *hPN = NULL;
290 AliTHn *hNP = NULL;
291 AliTHn *hPP = NULL;
292 AliTHn *hNN = NULL;
293
17d10bdd 294 TString gHistPname = "fHistP";
295 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
296 TString gHistNname = "fHistN";
297 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
298 TString gHistPNname = "fHistPN";
299 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
300 TString gHistNPname = "fHistNP";
301 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
302 TString gHistPPname = "fHistPP";
303 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
304 TString gHistNNname = "fHistNN";
305 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
306
307 hP = (AliTHn*) list->FindObject(gHistPname.Data());
308 hN = (AliTHn*) list->FindObject(gHistNname.Data());
309 hPN = (AliTHn*) list->FindObject(gHistPNname.Data());
310 hNP = (AliTHn*) list->FindObject(gHistNPname.Data());
311 hPP = (AliTHn*) list->FindObject(gHistPPname.Data());
312 hNN = (AliTHn*) list->FindObject(gHistNNname.Data());
32e94079 313 hNN->Print();
314
a38fd7f3 315
316 //Create the AliBalancePsi object and fill it with the AliTHn objects
317 AliBalancePsi *b = new AliBalancePsi();
32e94079 318 b->SetEventClass(eventClass);
a38fd7f3 319 b->SetHistNp(hP);
320 b->SetHistNn(hN);
321 b->SetHistNpn(hPN);
322 b->SetHistNnp(hNP);
323 b->SetHistNpp(hPP);
324 b->SetHistNnn(hNN);
20006629 325 if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
52daf7b2 326
327 //balance function shuffling
328 AliTHn *hPShuffled = NULL;
329 AliTHn *hNShuffled = NULL;
330 AliTHn *hPNShuffled = NULL;
331 AliTHn *hNPShuffled = NULL;
332 AliTHn *hPPShuffled = NULL;
333 AliTHn *hNNShuffled = NULL;
334 if(listBFShuffled) {
335 //listBFShuffled->ls();
336
17d10bdd 337 gHistPname = "fHistP_shuffle";
338 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
339 gHistNname = "fHistN_shuffle";
340 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
341 gHistPNname = "fHistPN_shuffle";
342 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
343 gHistNPname = "fHistNP_shuffle";
344 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
345 gHistPPname = "fHistPP_shuffle";
346 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
347 gHistNNname = "fHistNN_shuffle";
348 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
349
350 hPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPname.Data());
52daf7b2 351 hPShuffled->SetName("gHistPShuffled");
17d10bdd 352 hNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNname.Data());
52daf7b2 353 hNShuffled->SetName("gHistNShuffled");
17d10bdd 354 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPNname.Data());
52daf7b2 355 hPNShuffled->SetName("gHistPNShuffled");
17d10bdd 356 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNPname.Data());
52daf7b2 357 hNPShuffled->SetName("gHistNPShuffled");
17d10bdd 358 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPPname.Data());
52daf7b2 359 hPPShuffled->SetName("gHistPPShuffled");
17d10bdd 360 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNNname.Data());
52daf7b2 361 hNNShuffled->SetName("gHistNNShuffled");
362
363 AliBalancePsi *bShuffled = new AliBalancePsi();
32e94079 364 bShuffled->SetEventClass(eventClass);
52daf7b2 365 bShuffled->SetHistNp(hPShuffled);
366 bShuffled->SetHistNn(hNShuffled);
367 bShuffled->SetHistNpn(hPNShuffled);
368 bShuffled->SetHistNnp(hNPShuffled);
369 bShuffled->SetHistNpp(hPPShuffled);
370 bShuffled->SetHistNnn(hNNShuffled);
20006629 371 if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
52daf7b2 372 }
373
374 //balance function mixing
375 AliTHn *hPMixed = NULL;
376 AliTHn *hNMixed = NULL;
377 AliTHn *hPNMixed = NULL;
378 AliTHn *hNPMixed = NULL;
379 AliTHn *hPPMixed = NULL;
380 AliTHn *hNNMixed = NULL;
381
382 if(listBFMixed) {
383 //listBFMixed->ls();
384
17d10bdd 385 gHistPname = "fHistP";
386 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
387 gHistNname = "fHistN";
388 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
389 gHistPNname = "fHistPN";
390 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
391 gHistNPname = "fHistNP";
392 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
393 gHistPPname = "fHistPP";
394 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
395 gHistNNname = "fHistNN";
396 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
397 hPMixed = (AliTHn*) listBFMixed->FindObject(gHistPname.Data());
52daf7b2 398 hPMixed->SetName("gHistPMixed");
17d10bdd 399 hNMixed = (AliTHn*) listBFMixed->FindObject(gHistNname.Data());
52daf7b2 400 hNMixed->SetName("gHistNMixed");
17d10bdd 401 hPNMixed = (AliTHn*) listBFMixed->FindObject(gHistPNname.Data());
52daf7b2 402 hPNMixed->SetName("gHistPNMixed");
17d10bdd 403 hNPMixed = (AliTHn*) listBFMixed->FindObject(gHistNPname.Data());
52daf7b2 404 hNPMixed->SetName("gHistNPMixed");
17d10bdd 405 hPPMixed = (AliTHn*) listBFMixed->FindObject(gHistPPname.Data());
52daf7b2 406 hPPMixed->SetName("gHistPPMixed");
17d10bdd 407 hNNMixed = (AliTHn*) listBFMixed->FindObject(gHistNNname.Data());
52daf7b2 408 hNNMixed->SetName("gHistNNMixed");
409
410 AliBalancePsi *bMixed = new AliBalancePsi();
32e94079 411 bMixed->SetEventClass(eventClass);
52daf7b2 412 bMixed->SetHistNp(hPMixed);
413 bMixed->SetHistNn(hNMixed);
414 bMixed->SetHistNpn(hPNMixed);
415 bMixed->SetHistNnp(hNPMixed);
416 bMixed->SetHistNpp(hPPMixed);
417 bMixed->SetHistNnn(hNNMixed);
20006629 418 if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
52daf7b2 419 }
420
3ce45c12 421
77b9ff18 422 TH2D *gHistPN[4];
423 TH2D *gHistNP[4];
424 TH2D *gHistPP[4];
425 TH2D *gHistNN[4];
a38fd7f3 426
77b9ff18 427 TCanvas *cPN[4];
428 TCanvas *cNP[4];
429 TCanvas *cPP[4];
430 TCanvas *cNN[4];
a38fd7f3 431 TString histoTitle, pngName;
3ce45c12 432
94415ad1 433 // need event statistics for per-trigger yield
434 TH2D* hVertexCentrality = NULL;
435 if( normToTrig && !listQA){
436 Printf("per-trigger yield option chosen, but QA list for event statistics not available");
437 return;
438 }
439 else{
440 hVertexCentrality = (TH2D*)listQA->FindObject("fHistVz");
441 if(!hVertexCentrality){
442 Printf("per-trigger yield option chosen, but QA histogram for event statistics not available");
443 return;
444 }
445 }
446
3ce45c12 447 // if no mixing then divide by convoluted histograms
448 if(!listBFMixed && listQA){
3ce45c12 449
e5c5c33f 450 if(!listQA->FindObject("fHistEtaPhiPos") || !listQA->FindObject("fHistEtaPhiNeg")){
3ce45c12 451
e5c5c33f 452 Printf("fHistEtaPhiPos or fHistEtaPhiNeg not found! --> no convolution");
453 listQA = NULL;
3ce45c12 454
e5c5c33f 455 }
456 else{
3ce45c12 457
e5c5c33f 458 TH2D* fHistPos = (TH2D*)((TH3D*)listQA->FindObject("fHistEtaPhiPos"))->Project3D("xy");
459 fHistPos->GetYaxis()->SetRangeUser(-0.79,0.79);
e5c5c33f 460
461 TH2D* fHistNeg = (TH2D*)((TH3D*)listQA->FindObject("fHistEtaPhiNeg"))->Project3D("xy");
462 fHistNeg->GetYaxis()->SetRangeUser(-0.79,0.79);
e5c5c33f 463
464 gHistPN[2] = convolute2D(fHistPos, fHistNeg, "hConvPN");
465 gHistPN[2]->Scale(1./gHistPN[2]->GetBinContent(gHistPN[2]->FindBin(0,0)));
466
467 gHistNP[2] = convolute2D(fHistNeg, fHistPos, "hConvNP");
468 gHistNP[2]->Scale(1./gHistNP[2]->GetBinContent(gHistNP[2]->FindBin(0,0)));
469
470 gHistNN[2] = convolute2D(fHistNeg, fHistNeg, "hConvNN");
471 gHistNN[2]->Scale(1./gHistNN[2]->GetBinContent(gHistNN[2]->FindBin(0,0)));
472
473 gHistPP[2] = convolute2D(fHistPos, fHistPos, "hConvPP");
474 gHistPP[2]->Scale(1./gHistPP[2]->GetBinContent(gHistPP[2]->FindBin(0,0)));
475 }
3ce45c12 476 }
6acdbcb2 477
52daf7b2 478 //(+-)
32e94079 479 if(eventClass == "Centrality"){
480 histoTitle = "(+-) | Centrality: ";
481 histoTitle += psiMin;
482 histoTitle += " - ";
483 histoTitle += psiMax;
484 histoTitle += " % ";
3b8a460b 485 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 486 }
487 else if(eventClass == "Multiplicity"){
488 histoTitle = "(+-) | Multiplicity: ";
489 histoTitle += psiMin;
490 histoTitle += " - ";
491 histoTitle += psiMax;
492 histoTitle += " tracks";
3b8a460b 493 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 494 }
495 else{ // "EventPlane" (default)
496 histoTitle = "(+-) | Centrality: ";
497 histoTitle += centralityArray[gCentrality-1];
498 histoTitle += "%";
499 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 500 histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
32e94079 501 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 502 histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
32e94079 503 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 504 histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
32e94079 505 else
3b8a460b 506 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 507 }
20006629 508 gHistPN[0] = b->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 509 if(rebinEta > 1 || rebinPhi > 1){
510 gHistPN[0]->Rebin2D(rebinEta,rebinPhi);
511 gHistPN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
512 }
52daf7b2 513 gHistPN[0]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 514 gHistPN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 515 gHistPN[0]->SetTitle(histoTitle.Data());
516 cPN[0] = new TCanvas("cPN0","",0,0,600,500);
517 cPN[0]->SetFillColor(10); cPN[0]->SetHighLightColor(10);
5de9ad1a 518 gHistPN[0]->DrawCopy("surf1fb");
519 gPad->SetTheta(30); // default is 30
520 //gPad->SetPhi(130); // default is 30
521 gPad->SetPhi(-60); // default is 30
522 gPad->Update();
6acdbcb2 523 pngName = "DeltaPhiDeltaEta.Centrality";
52daf7b2 524 pngName += centralityArray[gCentrality-1];
6acdbcb2 525 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
526 pngName += ".PositiveNegative.png";
5c9b4733 527 //cPN[0]->SaveAs(pngName.Data());
6acdbcb2 528
52daf7b2 529 if(listBFShuffled) {
34616eda 530
531 histoTitle.ReplaceAll("(+-)","(+-) shuffled");
52daf7b2 532
20006629 533 gHistPN[1] = bShuffled->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 534 if(rebinEta > 1 || rebinPhi > 1){
535 gHistPN[1]->Rebin2D(rebinEta,rebinPhi);
536 gHistPN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
537 }
52daf7b2 538 gHistPN[1]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 539 gHistPN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 540 gHistPN[1]->SetTitle(histoTitle.Data());
541 cPN[1] = new TCanvas("cPN1","",0,100,600,500);
542 cPN[1]->SetFillColor(10);
543 cPN[1]->SetHighLightColor(10);
5de9ad1a 544 gHistPN[1]->DrawCopy("surf1fb");
545 gPad->SetTheta(30); // default is 30
546 //gPad->SetPhi(130); // default is 30
547 gPad->SetPhi(-60); // default is 30
548 gPad->Update();
52daf7b2 549 pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
550 pngName += centralityArray[gCentrality-1];
551 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
552 pngName += ".PositiveNegative.png";
5c9b4733 553 //cPN[1]->SaveAs(pngName.Data());
52daf7b2 554 }
555
556 if(listBFMixed) {
34616eda 557
558 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) mixed");
559 else histoTitle.ReplaceAll("(+-)","(+-) mixed");
52daf7b2 560
742af4bd 561 // if normalization to trigger then do not divide Event mixing by number of trigger particles
20006629 562 gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 563 if(rebinEta > 1 || rebinPhi > 1){
564 gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
565 gHistPN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
566 }
742af4bd 567
568 // normalization to 1 at (0,0) --> Jan Fietes method
569 if(normToTrig){
bd7bf11d
MW		 570 Double_t mixedNorm = gHistPN[2]->Integral(gHistPN[2]->GetXaxis()->FindBin(0-10e-5),gHistPN[2]->GetXaxis()->FindBin(0+10e-5),gHistPN[2]->GetNbinsY()/2 + 1,gHistPN[2]->GetNbinsY());
571 mixedNorm /= 0.5 * gHistPN[2]->GetNbinsY()*(gHistPN[2]->GetXaxis()->FindBin(0.01) - gHistPN[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 572
573 // finite bin correction
574 Double_t binWidthEta = gHistPN[2]->GetXaxis()->GetBinWidth(gHistPN[2]->GetNbinsX());
575 Double_t maxEta = gHistPN[2]->GetXaxis()->GetBinUpEdge(gHistPN[2]->GetNbinsX());
576
577 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
578 //Printf("Finite bin correction: %f", finiteBinCorrection);
579 mixedNorm /= finiteBinCorrection;
580
742af4bd 581 gHistPN[2]->Scale(1./mixedNorm);
582 }
583
52daf7b2 584 gHistPN[2]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 585 gHistPN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 586 gHistPN[2]->SetTitle(histoTitle.Data());
587 cPN[2] = new TCanvas("cPN2","",0,200,600,500);
588 cPN[2]->SetFillColor(10);
589 cPN[2]->SetHighLightColor(10);
5de9ad1a 590 gHistPN[2]->DrawCopy("surf1fb");
591 gPad->SetTheta(30); // default is 30
592 //gPad->SetPhi(130); // default is 30
593 gPad->SetPhi(-60); // default is 30
594 gPad->Update();
52daf7b2 595 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
596 pngName += centralityArray[gCentrality-1];
597 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
598 pngName += ".PositiveNegative.png";
5c9b4733 599 //cPN[2]->SaveAs(pngName.Data());
77b9ff18 600
601 //Correlation function (+-)
3aa5487c 602 gHistPN[3] = b->GetCorrelationFunction("PN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed,normToTrig,normalizationRangePhi,hVertexCentrality);
77b9ff18 603 gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 604 if(normToTrig)
605 gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
606 else
607 gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
2b4abbc5 608 //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
77b9ff18 609 cPN[3] = new TCanvas("cPN3","",0,300,600,500);
610 cPN[3]->SetFillColor(10);
611 cPN[3]->SetHighLightColor(10);
612 gHistPN[3]->DrawCopy("surf1fb");
613 gPad->SetTheta(30); // default is 30
614 //gPad->SetPhi(130); // default is 30
615 gPad->SetPhi(-60); // default is 30
616 gPad->Update();
617 pngName = "CorrelationFunction.Centrality";
618 pngName += centralityArray[gCentrality-1];
619 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
620 pngName += ".PositiveNegative.png";
5c9b4733 621 //cPN[3]->SaveAs(pngName.Data());
52daf7b2 622 }
3ce45c12 623 // if no mixing then divide by convoluted histograms
624 else if(listQA){
34616eda 625
626 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) convoluted");
627 else histoTitle.ReplaceAll("(+-)","(+-) convoluted");
3ce45c12 628
faa03677 629 if(rebinEta > 1 || rebinPhi > 1){
630 gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
631 gHistPN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
632 }
633
3ce45c12 634 // normalization to 1 at (0,0) --> Jan Fietes method
635 if(normToTrig){
bd7bf11d
MW		 636 Double_t mixedNorm = gHistPN[2]->Integral(gHistPN[2]->GetXaxis()->FindBin(0-10e-5),gHistPN[2]->GetXaxis()->FindBin(0+10e-5),gHistPN[2]->GetNbinsY()/2 + 1,gHistPN[2]->GetNbinsY());
637 mixedNorm /= 0.5 * gHistPN[2]->GetNbinsY()*(gHistPN[2]->GetXaxis()->FindBin(0.01) - gHistPN[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 638
639 // finite bin correction
640 Double_t binWidthEta = gHistPN[2]->GetXaxis()->GetBinWidth(gHistPN[2]->GetNbinsX());
641 Double_t maxEta = gHistPN[2]->GetXaxis()->GetBinUpEdge(gHistPN[2]->GetNbinsX());
642
643 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
644 //Printf("Finite bin correction: %f", finiteBinCorrection);
645 mixedNorm /= finiteBinCorrection;
646
3ce45c12 647 gHistPN[2]->Scale(1./mixedNorm);
648 }
649
650 gHistPN[2]->GetYaxis()->SetTitleOffset(1.5);
651 gHistPN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
652 gHistPN[2]->SetTitle(histoTitle.Data());
653 cPN[2] = new TCanvas("cPN2","",0,200,600,500);
654 cPN[2]->SetFillColor(10);
655 cPN[2]->SetHighLightColor(10);
656 gHistPN[2]->DrawCopy("surf1fb");
657 gPad->SetTheta(30); // default is 30
658 //gPad->SetPhi(130); // default is 30
659 gPad->SetPhi(-60); // default is 30
660 gPad->Update();
661 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
662 pngName += centralityArray[gCentrality-1];
663 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
664 pngName += ".PositiveNegative.png";
665 //cPN[2]->SaveAs(pngName.Data());
666
667 //Correlation function (+-)
668 gHistPN[3] = dynamic_cast<TH2D *>(gHistPN[0]->Clone());
669 gHistPN[3]->Divide(gHistPN[2]);
670 gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 671 if(normToTrig)
672 gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
673 else
674 gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
2b4abbc5 675 //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
3ce45c12 676 cPN[3] = new TCanvas("cPN3","",0,300,600,500);
677 cPN[3]->SetFillColor(10);
678 cPN[3]->SetHighLightColor(10);
679 gHistPN[3]->DrawCopy("surf1fb");
680 gPad->SetTheta(30); // default is 30
681 //gPad->SetPhi(130); // default is 30
682 gPad->SetPhi(-60); // default is 30
683 gPad->Update();
684 pngName = "CorrelationFunction.Centrality";
685 pngName += centralityArray[gCentrality-1];
686 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
687 pngName += ".PositiveNegative.png";
688 //cPN[3]->SaveAs(pngName.Data());
689 }
52daf7b2 690
691 //(-+)
32e94079 692 if(eventClass == "Centrality"){
693 histoTitle = "(-+) | Centrality: ";
694 histoTitle += psiMin;
695 histoTitle += " - ";
696 histoTitle += psiMax;
697 histoTitle += " % ";
3b8a460b 698 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 699 }
700 else if(eventClass == "Multiplicity"){
701 histoTitle = "(-+) | Multiplicity: ";
702 histoTitle += psiMin;
703 histoTitle += " - ";
704 histoTitle += psiMax;
705 histoTitle += " tracks";
3b8a460b 706 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 707 }
708 else{ // "EventPlane" (default)
709 histoTitle = "(-+) | Centrality: ";
710 histoTitle += centralityArray[gCentrality-1];
711 histoTitle += "%";
712 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 713 histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
32e94079 714 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 715 histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
32e94079 716 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 717 histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
32e94079 718 else
3b8a460b 719 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 720 }
52daf7b2 721
20006629 722 gHistNP[0] = b->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 723 if(rebinEta > 1 || rebinPhi > 1){
724 gHistNP[0]->Rebin2D(rebinEta,rebinPhi);
725 gHistNP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
726 }
52daf7b2 727 gHistNP[0]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 728 gHistNP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 729 gHistNP[0]->SetTitle(histoTitle.Data());
730 cNP[0] = new TCanvas("cNP0","",100,0,600,500);
731 cNP[0]->SetFillColor(10);
732 cNP[0]->SetHighLightColor(10);
5de9ad1a 733 gHistNP[0]->DrawCopy("surf1fb");
734 gPad->SetTheta(30); // default is 30
735 //gPad->SetPhi(130); // default is 30
736 gPad->SetPhi(-60); // default is 30
737 gPad->Update();
6acdbcb2 738 pngName = "DeltaPhiDeltaEta.Centrality";
52daf7b2 739 pngName += centralityArray[gCentrality-1];
6acdbcb2 740 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
741 pngName += ".NegativePositive.png";
5c9b4733 742 //cNP[0]->SaveAs(pngName.Data());
52daf7b2 743
744 if(listBFShuffled) {
34616eda 745
746 histoTitle.ReplaceAll("(-+)","(-+) shuffled");
52daf7b2 747
20006629 748 gHistNP[1] = bShuffled->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 749 if(rebinEta > 1 || rebinPhi > 1){
750 gHistNP[1]->Rebin2D(rebinEta,rebinPhi);
751 gHistNP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
752 }
52daf7b2 753 gHistNP[1]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 754 gHistNP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 755 gHistNP[1]->SetTitle(histoTitle.Data());
756 cNP[1] = new TCanvas("cNP1","",100,100,600,500);
757 cNP[1]->SetFillColor(10);
758 cNP[1]->SetHighLightColor(10);
5de9ad1a 759 gHistNP[1]->DrawCopy("surf1fb");
760 gPad->SetTheta(30); // default is 30
761 //gPad->SetPhi(130); // default is 30
762 gPad->SetPhi(-60); // default is 30
763 gPad->Update();
52daf7b2 764 pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
765 pngName += centralityArray[gCentrality-1];
766 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
767 pngName += ".NegativePositive.png";
5c9b4733 768 //cNP[1]->SaveAs(pngName.Data());
52daf7b2 769 }
770
771 if(listBFMixed) {
34616eda 772
773 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) mixed");
774 else histoTitle.ReplaceAll("(-+)","(-+) mixed");
775
742af4bd 776 // if normalization to trigger then do not divide Event mixing by number of trigger particles
20006629 777 gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 778 if(rebinEta > 1 || rebinPhi > 1){
779 gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
780 gHistNP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
781 }
742af4bd 782 // normalization to 1 at (0,0) --> Jan Fietes method
783 if(normToTrig){
bd7bf11d
MW		 784 Double_t mixedNorm = gHistNP[2]->Integral(gHistNP[2]->GetXaxis()->FindBin(0-10e-5),gHistNP[2]->GetXaxis()->FindBin(0+10e-5),gHistNP[2]->GetNbinsY()/2 + 1,gHistNP[2]->GetNbinsY());
785 mixedNorm /= 0.5 * gHistNP[2]->GetNbinsY()*(gHistNP[2]->GetXaxis()->FindBin(0.01) - gHistNP[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 786
787 // finite bin correction
788 Double_t binWidthEta = gHistNP[2]->GetXaxis()->GetBinWidth(gHistNP[2]->GetNbinsX());
789 Double_t maxEta = gHistNP[2]->GetXaxis()->GetBinUpEdge(gHistNP[2]->GetNbinsX());
790
791 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
792 //Printf("Finite bin correction: %f", finiteBinCorrection);
793 mixedNorm /= finiteBinCorrection;
794
742af4bd 795 gHistNP[2]->Scale(1./mixedNorm);
796 }
797
52daf7b2 798 gHistNP[2]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 799 gHistNP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 800 gHistNP[2]->SetTitle(histoTitle.Data());
801 cNP[2] = new TCanvas("cNP2","",100,200,600,500);
802 cNP[2]->SetFillColor(10);
803 cNP[2]->SetHighLightColor(10);
5de9ad1a 804 gHistNP[2]->DrawCopy("surf1fb");
805 gPad->SetTheta(30); // default is 30
806 //gPad->SetPhi(130); // default is 30
807 gPad->SetPhi(-60); // default is 30
808 gPad->Update();
52daf7b2 809 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
810 pngName += centralityArray[gCentrality-1];
811 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
812 pngName += ".NegativePositive.png";
5c9b4733 813 //cNP[2]->SaveAs(pngName.Data());
77b9ff18 814
815 //Correlation function (-+)
3aa5487c 816 gHistNP[3] = b->GetCorrelationFunction("NP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed,normToTrig,normalizationRangePhi,hVertexCentrality);
77b9ff18 817 gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 818 if(normToTrig)
819 gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
820 else
821 gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
2b4abbc5 822 //gHistNP[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
77b9ff18 823 cNP[3] = new TCanvas("cNP3","",100,300,600,500);
824 cNP[3]->SetFillColor(10);
825 cNP[3]->SetHighLightColor(10);
826 gHistNP[3]->DrawCopy("surf1fb");
827 gPad->SetTheta(30); // default is 30
828 //gPad->SetPhi(130); // default is 30
829 gPad->SetPhi(-60); // default is 30
830 gPad->Update();
831 pngName = "CorrelationFunction.Centrality";
832 pngName += centralityArray[gCentrality-1];
833 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
834 pngName += ".NegativePositive.png";
5c9b4733 835 //cNP[3]->SaveAs(pngName.Data());
52daf7b2 836 }
3ce45c12 837 // if no mixing then divide by convoluted histograms
838 else if(listQA){
34616eda 839
840 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) convoluted");
841 else histoTitle.ReplaceAll("(-+)","(-+) convoluted");
3ce45c12 842
faa03677 843 if(rebinEta > 1 || rebinPhi > 1){
844 gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
845 gHistNP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
846 }
3ce45c12 847
848 // normalization to 1 at (0,0) --> Jan Fietes method
849 if(normToTrig){
bd7bf11d
MW		 850 Double_t mixedNorm = gHistNP[2]->Integral(gHistNP[2]->GetXaxis()->FindBin(0-10e-5),gHistNP[2]->GetXaxis()->FindBin(0+10e-5),gHistNP[2]->GetNbinsY()/2 + 1,gHistNP[2]->GetNbinsY());
851 mixedNorm /= 0.5 * gHistNP[2]->GetNbinsY()*(gHistNP[2]->GetXaxis()->FindBin(0.01) - gHistNP[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 852
853 // finite bin correction
854 Double_t binWidthEta = gHistNP[2]->GetXaxis()->GetBinWidth(gHistNP[2]->GetNbinsX());
855 Double_t maxEta = gHistNP[2]->GetXaxis()->GetBinUpEdge(gHistNP[2]->GetNbinsX());
856
857 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
858 //Printf("Finite bin correction: %f", finiteBinCorrection);
859 mixedNorm /= finiteBinCorrection;
860
3ce45c12 861 gHistNP[2]->Scale(1./mixedNorm);
862 }
863
864 gHistNP[2]->GetYaxis()->SetTitleOffset(1.5);
865 gHistNP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
866 gHistNP[2]->SetTitle(histoTitle.Data());
867 cNP[2] = new TCanvas("cNP2","",100,200,600,500);
868 cNP[2]->SetFillColor(10);
869 cNP[2]->SetHighLightColor(10);
870 gHistNP[2]->DrawCopy("surf1fb");
871 gPad->SetTheta(30); // default is 30
872 //gPad->SetPhi(130); // default is 30
873 gPad->SetPhi(-60); // default is 30
874 gPad->Update();
875 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
876 pngName += centralityArray[gCentrality-1];
877 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
878 pngName += ".NegativePositive.png";
879 //cNP[2]->SaveAs(pngName.Data());
880
881 //Correlation function (-+)
882 gHistNP[3] = dynamic_cast<TH2D *>(gHistNP[0]->Clone());
883 gHistNP[3]->Divide(gHistNP[2]);
884 gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 885 if(normToTrig)
886 gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
887 else
888 gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
2b4abbc5 889 //gHistNP[3]->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
3ce45c12 890 cNP[3] = new TCanvas("cNP3","",100,300,600,500);
891 cNP[3]->SetFillColor(10);
892 cNP[3]->SetHighLightColor(10);
893 gHistNP[3]->DrawCopy("surf1fb");
894 gPad->SetTheta(30); // default is 30
895 //gPad->SetPhi(130); // default is 30
896 gPad->SetPhi(-60); // default is 30
897 gPad->Update();
898 pngName = "CorrelationFunction.Centrality";
899 pngName += centralityArray[gCentrality-1];
900 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
901 pngName += ".NegativePositive.png";
902 //cNP[3]->SaveAs(pngName.Data());
903 }
904
905
52daf7b2 906 //(++)
32e94079 907 if(eventClass == "Centrality"){
908 histoTitle = "(++) | Centrality: ";
909 histoTitle += psiMin;
910 histoTitle += " - ";
911 histoTitle += psiMax;
912 histoTitle += " % ";
3b8a460b 913 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 914 }
915 else if(eventClass == "Multiplicity"){
916 histoTitle = "(++) | Multiplicity: ";
917 histoTitle += psiMin;
918 histoTitle += " - ";
919 histoTitle += psiMax;
920 histoTitle += " tracks";
3b8a460b 921 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 922 }
923 else{ // "EventPlane" (default)
924 histoTitle = "(++) | Centrality: ";
925 histoTitle += centralityArray[gCentrality-1];
926 histoTitle += "%";
927 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 928 histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
32e94079 929 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 930 histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
32e94079 931 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 932 histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
32e94079 933 else
3b8a460b 934 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 935 }
52daf7b2 936
20006629 937 gHistPP[0] = b->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 938 if(rebinEta > 1 || rebinPhi > 1){
939 gHistPP[0]->Rebin2D(rebinEta,rebinPhi);
940 gHistPP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
941 }
52daf7b2 942 gHistPP[0]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 943 gHistPP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 944 gHistPP[0]->SetTitle(histoTitle.Data());
945 cPP[0] = new TCanvas("cPP0","",200,0,600,500);
946 cPP[0]->SetFillColor(10);
947 cPP[0]->SetHighLightColor(10);
5de9ad1a 948 gHistPP[0]->DrawCopy("surf1fb");
949 gPad->SetTheta(30); // default is 30
950 //gPad->SetPhi(130); // default is 30
951 gPad->SetPhi(-60); // default is 30
952 gPad->Update();
6acdbcb2 953 pngName = "DeltaPhiDeltaEta.Centrality";
52daf7b2 954 pngName += centralityArray[gCentrality-1];
6acdbcb2 955 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
956 pngName += ".PositivePositive.png";
5c9b4733 957 //cPP[0]->SaveAs(pngName.Data());
6acdbcb2 958
52daf7b2 959 if(listBFShuffled) {
34616eda 960
961 histoTitle.ReplaceAll("(++)","(++) shuffled");
52daf7b2 962
20006629 963 gHistPP[1] = bShuffled->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 964 if(rebinEta > 1 || rebinPhi > 1){
965 gHistPP[1]->Rebin2D(rebinEta,rebinPhi);
966 gHistPP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
967 }
52daf7b2 968 gHistPP[1]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 969 gHistPP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 970 gHistPP[1]->SetTitle(histoTitle.Data());
971 cPP[1] = new TCanvas("cPP1","",200,100,600,500);
972 cPP[1]->SetFillColor(10);
973 cPP[1]->SetHighLightColor(10);
5de9ad1a 974 gHistPP[1]->DrawCopy("surf1fb");
975 gPad->SetTheta(30); // default is 30
976 //gPad->SetPhi(130); // default is 30
977 gPad->SetPhi(-60); // default is 30
978 gPad->Update();
52daf7b2 979 pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
980 pngName += centralityArray[gCentrality-1];
981 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
982 pngName += ".PositivePositive.png";
5c9b4733 983 //cPP[1]->SaveAs(pngName.Data());
52daf7b2 984 }
985
986 if(listBFMixed) {
34616eda 987
988 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) mixed");
989 else histoTitle.ReplaceAll("(++)","(++) mixed");
52daf7b2 990
742af4bd 991 // if normalization to trigger then do not divide Event mixing by number of trigger particles
20006629 992 gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 993 if(rebinEta > 1 || rebinPhi > 1){
994 gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
995 gHistPP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
996 }
742af4bd 997 // normalization to 1 at (0,0) --> Jan Fietes method
998 if(normToTrig){
bd7bf11d
MW		 999 Double_t mixedNorm = gHistPP[2]->Integral(gHistPP[2]->GetXaxis()->FindBin(0-10e-5),gHistPP[2]->GetXaxis()->FindBin(0+10e-5),gHistPP[2]->GetNbinsY()/2 + 1,gHistPP[2]->GetNbinsY());
1000 mixedNorm /= 0.5 * gHistPP[2]->GetNbinsY()*(gHistPP[2]->GetXaxis()->FindBin(0.01) - gHistPP[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 1001
1002 // finite bin correction
1003 Double_t binWidthEta = gHistPP[2]->GetXaxis()->GetBinWidth(gHistPP[2]->GetNbinsX());
1004 Double_t maxEta = gHistPP[2]->GetXaxis()->GetBinUpEdge(gHistPP[2]->GetNbinsX());
1005
1006 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1007 //Printf("Finite bin correction: %f", finiteBinCorrection);
1008 mixedNorm /= finiteBinCorrection;
1009
742af4bd 1010 gHistPP[2]->Scale(1./mixedNorm);
1011 }
1012
52daf7b2 1013 gHistPP[2]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 1014 gHistPP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 1015 gHistPP[2]->SetTitle(histoTitle.Data());
1016 cPP[2] = new TCanvas("cPP2","",200,200,600,500);
1017 cPP[2]->SetFillColor(10);
1018 cPP[2]->SetHighLightColor(10);
5de9ad1a 1019 gHistPP[2]->DrawCopy("surf1fb");
1020 gPad->SetTheta(30); // default is 30
1021 //gPad->SetPhi(130); // default is 30
1022 gPad->SetPhi(-60); // default is 30
1023 gPad->Update();
52daf7b2 1024 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1025 pngName += centralityArray[gCentrality-1];
1026 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1027 pngName += ".PositivePositive.png";
5c9b4733 1028 //cPP[2]->SaveAs(pngName.Data());
77b9ff18 1029
1030 //Correlation function (++)
3aa5487c 1031 gHistPP[3] = b->GetCorrelationFunction("PP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed,normToTrig,normalizationRangePhi,hVertexCentrality);
77b9ff18 1032 gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 1033 if(normToTrig)
1034 gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1035 else
1036 gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
2b4abbc5 1037 // gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
77b9ff18 1038 cPP[3] = new TCanvas("cPP3","",200,300,600,500);
1039 cPP[3]->SetFillColor(10);
1040 cPP[3]->SetHighLightColor(10);
1041 gHistPP[3]->DrawCopy("surf1fb");
1042 gPad->SetTheta(30); // default is 30
1043 //gPad->SetPhi(130); // default is 30
1044 gPad->SetPhi(-60); // default is 30
1045 gPad->Update();
1046 pngName = "CorrelationFunction.Centrality";
1047 pngName += centralityArray[gCentrality-1];
1048 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1049 pngName += ".PositivePositive.png";
5c9b4733 1050 //cPP[3]->SaveAs(pngName.Data());
52daf7b2 1051 }
3ce45c12 1052 // if no mixing then divide by convoluted histograms
1053 else if(listQA){
34616eda 1054
1055 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) convoluted");
1056 else histoTitle.ReplaceAll("(++)","(++) convoluted");
3ce45c12 1057
faa03677 1058 if(rebinEta > 1 || rebinPhi > 1){
1059 gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
1060 gHistPP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1061 }
3ce45c12 1062 // normalization to 1 at (0,0) --> Jan Fietes method
1063 if(normToTrig){
bd7bf11d
MW		 1064 Double_t mixedNorm = gHistPP[2]->Integral(gHistPP[2]->GetXaxis()->FindBin(0-10e-5),gHistPP[2]->GetXaxis()->FindBin(0+10e-5),gHistPP[2]->GetNbinsY()/2 + 1,gHistPP[2]->GetNbinsY());
1065 mixedNorm /= 0.5 * gHistPP[2]->GetNbinsY()*(gHistPP[2]->GetXaxis()->FindBin(0.01) - gHistPP[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 1066
1067 // finite bin correction
1068 Double_t binWidthEta = gHistPP[2]->GetXaxis()->GetBinWidth(gHistPP[2]->GetNbinsX());
1069 Double_t maxEta = gHistPP[2]->GetXaxis()->GetBinUpEdge(gHistPP[2]->GetNbinsX());
1070
1071 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1072 //Printf("Finite bin correction: %f", finiteBinCorrection);
1073 mixedNorm /= finiteBinCorrection;
1074
3ce45c12 1075 gHistPP[2]->Scale(1./mixedNorm);
1076 }
1077
1078 gHistPP[2]->GetYaxis()->SetTitleOffset(1.5);
1079 gHistPP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1080 gHistPP[2]->SetTitle(histoTitle.Data());
1081 cPP[2] = new TCanvas("cPP2","",200,200,600,500);
1082 cPP[2]->SetFillColor(10);
1083 cPP[2]->SetHighLightColor(10);
1084 gHistPP[2]->DrawCopy("surf1fb");
1085 gPad->SetTheta(30); // default is 30
1086 //gPad->SetPhi(130); // default is 30
1087 gPad->SetPhi(-60); // default is 30
1088 gPad->Update();
1089 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1090 pngName += centralityArray[gCentrality-1];
1091 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1092 pngName += ".PositivePositive.png";
1093 //cPP[2]->SaveAs(pngName.Data());
1094
1095 //Correlation function (++)
1096 gHistPP[3] = dynamic_cast<TH2D *>(gHistPP[0]->Clone());
1097 gHistPP[3]->Divide(gHistPP[2]);
1098 gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 1099 if(normToTrig)
1100 gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1101 else
1102 gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
2b4abbc5 1103 //gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
3ce45c12 1104 cPP[3] = new TCanvas("cPP3","",200,300,600,500);
1105 cPP[3]->SetFillColor(10);
1106 cPP[3]->SetHighLightColor(10);
1107 gHistPP[3]->DrawCopy("surf1fb");
1108 gPad->SetTheta(30); // default is 30
1109 //gPad->SetPhi(130); // default is 30
1110 gPad->SetPhi(-60); // default is 30
1111 gPad->Update();
1112 pngName = "CorrelationFunction.Centrality";
1113 pngName += centralityArray[gCentrality-1];
1114 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1115 pngName += ".PositivePositive.png";
1116 //cPP[3]->SaveAs(pngName.Data());
1117 }
52daf7b2 1118
1119 //(--)
32e94079 1120 if(eventClass == "Centrality"){
1121 histoTitle = "(--) | Centrality: ";
1122 histoTitle += psiMin;
1123 histoTitle += " - ";
1124 histoTitle += psiMax;
1125 histoTitle += " % ";
3b8a460b 1126 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 1127 }
1128 else if(eventClass == "Multiplicity"){
1129 histoTitle = "(--) | Multiplicity: ";
1130 histoTitle += psiMin;
1131 histoTitle += " - ";
1132 histoTitle += psiMax;
1133 histoTitle += " tracks";
3b8a460b 1134 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 1135 }
1136 else{ // "EventPlane" (default)
1137 histoTitle = "(--) | Centrality: ";
1138 histoTitle += centralityArray[gCentrality-1];
1139 histoTitle += "%";
1140 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 1141 histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
32e94079 1142 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 1143 histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
32e94079 1144 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 1145 histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
32e94079 1146 else
3b8a460b 1147 histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
32e94079 1148 }
52daf7b2 1149
20006629 1150 gHistNN[0] = b->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 1151 if(rebinEta > 1 || rebinPhi > 1){
1152 gHistNN[0]->Rebin2D(rebinEta,rebinPhi);
1153 gHistNN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1154 }
52daf7b2 1155 gHistNN[0]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 1156 gHistNN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 1157 gHistNN[0]->SetTitle(histoTitle.Data());
1158 cNN[0] = new TCanvas("cNN0","",300,0,600,500);
1159 cNN[0]->SetFillColor(10);
1160 cNN[0]->SetHighLightColor(10);
5de9ad1a 1161 gHistNN[0]->DrawCopy("surf1fb");
1162 gPad->SetTheta(30); // default is 30
1163 gPad->SetPhi(-60); // default is 30
1164 //gPad->SetPhi(-60); // default is 30
1165 gPad->Update();
6acdbcb2 1166 pngName = "DeltaPhiDeltaEta.Centrality";
52daf7b2 1167 pngName += centralityArray[gCentrality-1];
6acdbcb2 1168 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1169 pngName += ".NegativeNegative.png";
5c9b4733 1170 //cNN[0]->SaveAs(pngName.Data());
52daf7b2 1171
1172 if(listBFShuffled) {
34616eda 1173
1174 histoTitle.ReplaceAll("(--)","(--) shuffled");
52daf7b2 1175
20006629 1176 gHistNN[1] = bShuffled->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 1177 if(rebinEta > 1 || rebinPhi > 1){
1178 gHistNN[1]->Rebin2D(rebinEta,rebinPhi);
1179 gHistNN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1180 }
52daf7b2 1181 gHistNN[1]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 1182 gHistNN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 1183 gHistNN[1]->SetTitle(histoTitle.Data());
1184 cNN[1] = new TCanvas("cNN1","",300,100,600,500);
1185 cNN[1]->SetFillColor(10);
1186 cNN[1]->SetHighLightColor(10);
5de9ad1a 1187 gHistNN[1]->DrawCopy("surf1fb");
1188 gPad->SetTheta(30); // default is 30
1189 //gPad->SetPhi(130); // default is 30
1190 gPad->SetPhi(-60); // default is 30
1191 gPad->Update();
52daf7b2 1192 pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
1193 pngName += centralityArray[gCentrality-1];
1194 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1195 pngName += ".NegativeNegative.png";
5c9b4733 1196 //cNN[1]->SaveAs(pngName.Data());
52daf7b2 1197 }
1198
1199 if(listBFMixed) {
34616eda 1200
1201 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) mixed");
1202 else histoTitle.ReplaceAll("(--)","(--) mixed");
52daf7b2 1203
742af4bd 1204 // if normalization to trigger then do not divide Event mixing by number of trigger particles
20006629 1205 gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 1206 if(rebinEta > 1 || rebinPhi > 1){
1207 gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
1208 gHistNN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1209 }
742af4bd 1210 // normalization to 1 at (0,0) --> Jan Fietes method
1211 if(normToTrig){
bd7bf11d
MW		 1212 Double_t mixedNorm = gHistNN[2]->Integral(gHistNN[2]->GetXaxis()->FindBin(0-10e-5),gHistNN[2]->GetXaxis()->FindBin(0+10e-5),gHistNN[2]->GetNbinsY()/2 + 1,gHistNN[2]->GetNbinsY());
1213 mixedNorm /= 0.5 * gHistNN[2]->GetNbinsY()*(gHistNN[2]->GetXaxis()->FindBin(0.01) - gHistNN[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 1214
1215 // finite bin correction
1216 Double_t binWidthEta = gHistNN[2]->GetXaxis()->GetBinWidth(gHistNN[2]->GetNbinsX());
1217 Double_t maxEta = gHistNN[2]->GetXaxis()->GetBinUpEdge(gHistNN[2]->GetNbinsX());
1218
1219 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1220 //Printf("Finite bin correction: %f", finiteBinCorrection);
1221 mixedNorm /= finiteBinCorrection;
1222
742af4bd 1223 gHistNN[2]->Scale(1./mixedNorm);
1224 }
1225
52daf7b2 1226 gHistNN[2]->GetYaxis()->SetTitleOffset(1.5);
77b9ff18 1227 gHistNN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
52daf7b2 1228 gHistNN[2]->SetTitle(histoTitle.Data());
1229 cNN[2] = new TCanvas("cNN2","",300,200,600,500);
1230 cNN[2]->SetFillColor(10);
1231 cNN[2]->SetHighLightColor(10);
5de9ad1a 1232 gHistNN[2]->DrawCopy("surf1fb");
1233 gPad->SetTheta(30); // default is 30
1234 //gPad->SetPhi(130); // default is 30
1235 gPad->SetPhi(-60); // default is 30
1236 gPad->Update();
52daf7b2 1237 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1238 pngName += centralityArray[gCentrality-1];
1239 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1240 pngName += ".NegativeNegative.png";
5c9b4733 1241 //cNN[2]->SaveAs(pngName.Data());
77b9ff18 1242
1243 //Correlation function (--)
3aa5487c 1244 gHistNN[3] = b->GetCorrelationFunction("NN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed,normToTrig,normalizationRangePhi,hVertexCentrality);
77b9ff18 1245 gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 1246 if(normToTrig)
1247 gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1248 else
1249 gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
2b4abbc5 1250 // gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
77b9ff18 1251 cNN[3] = new TCanvas("cNN3","",300,300,600,500);
1252 cNN[3]->SetFillColor(10);
1253 cNN[3]->SetHighLightColor(10);
1254 gHistNN[3]->DrawCopy("surf1fb");
1255 gPad->SetTheta(30); // default is 30
1256 //gPad->SetPhi(130); // default is 30
1257 gPad->SetPhi(-60); // default is 30
1258 gPad->Update();
1259 pngName = "CorrelationFunction.Centrality";
1260 pngName += centralityArray[gCentrality-1];
1261 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1262 pngName += ".NegativeNegative.png";
5c9b4733 1263 //cNN[3]->SaveAs(pngName.Data());
52daf7b2 1264 }
3ce45c12 1265 // if no mixing then divide by convoluted histograms
1266 else if(listQA){
34616eda 1267
1268 if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) convoluted");
1269 else histoTitle.ReplaceAll("(--)","(--) convoluted");
3ce45c12 1270
faa03677 1271 if(rebinEta > 1 || rebinPhi > 1){
1272 gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
1273 gHistNN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1274 }
3ce45c12 1275 // normalization to 1 at (0,0) --> Jan Fietes method
1276 if(normToTrig){
bd7bf11d
MW		 1277 Double_t mixedNorm = gHistNN[2]->Integral(gHistNN[2]->GetXaxis()->FindBin(0-10e-5),gHistNN[2]->GetXaxis()->FindBin(0+10e-5),gHistNN[2]->GetNbinsY()/2 + 1,gHistNN[2]->GetNbinsY());
1278 mixedNorm /= 0.5 * gHistNN[2]->GetNbinsY()*(gHistNN[2]->GetXaxis()->FindBin(0.01) - gHistNN[2]->GetXaxis()->FindBin(-0.01) + 1);
5bb3cd60
MW		 1279
1280 // finite bin correction
1281 Double_t binWidthEta = gHistNN[2]->GetXaxis()->GetBinWidth(gHistNN[2]->GetNbinsX());
1282 Double_t maxEta = gHistNN[2]->GetXaxis()->GetBinUpEdge(gHistNN[2]->GetNbinsX());
1283
1284 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1285 //Printf("Finite bin correction: %f", finiteBinCorrection);
1286 mixedNorm /= finiteBinCorrection;
1287
3ce45c12 1288 gHistNN[2]->Scale(1./mixedNorm);
1289 }
1290
1291 gHistNN[2]->GetYaxis()->SetTitleOffset(1.5);
1292 gHistNN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1293 gHistNN[2]->SetTitle(histoTitle.Data());
1294 cNN[2] = new TCanvas("cNN2","",300,200,600,500);
1295 cNN[2]->SetFillColor(10);
1296 cNN[2]->SetHighLightColor(10);
1297 gHistNN[2]->DrawCopy("surf1fb");
1298 gPad->SetTheta(30); // default is 30
1299 //gPad->SetPhi(130); // default is 30
1300 gPad->SetPhi(-60); // default is 30
1301 gPad->Update();
1302 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1303 pngName += centralityArray[gCentrality-1];
1304 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1305 pngName += ".NegativeNegative.png";
1306 //cNN[2]->SaveAs(pngName.Data());
1307
1308 //Correlation function (--)
1309 gHistNN[3] = dynamic_cast<TH2D *>(gHistNN[0]->Clone());
1310 gHistNN[3]->Divide(gHistNN[2]);
1311 gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
41f2bc59 1312 if(normToTrig)
1313 gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1314 else
1315 gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
2b4abbc5 1316 //gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
3ce45c12 1317 cNN[3] = new TCanvas("cNN3","",300,300,600,500);
1318 cNN[3]->SetFillColor(10);
1319 cNN[3]->SetHighLightColor(10);
1320 gHistNN[3]->DrawCopy("surf1fb");
1321 gPad->SetTheta(30); // default is 30
1322 //gPad->SetPhi(130); // default is 30
1323 gPad->SetPhi(-60); // default is 30
1324 gPad->Update();
1325 pngName = "CorrelationFunction.Centrality";
1326 pngName += centralityArray[gCentrality-1];
1327 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1328 pngName += ".NegativeNegative.png";
1329 //cNN[3]->SaveAs(pngName.Data());
1330 }
5c9b4733 1331
1332 //Write to output file
32e94079 1333 TString newFileName = "correlationFunction.";
1334 if(eventClass == "Centrality"){
1335 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1336 newFileName += ".PsiAll.PttFrom";
1337 }
1338 else if(eventClass == "Multiplicity"){
1339 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1340 newFileName += ".PsiAll.PttFrom";
1341 }
1342 else{ // "EventPlane" (default)
1343 newFileName += "Centrality";
1344 newFileName += gCentrality; newFileName += ".Psi";
1345 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1346 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1347 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1348 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1349 else newFileName += "All.PttFrom";
1350 }
5c9b4733 1351 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1352 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1353 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1354 newFileName += Form("%.1f",ptAssociatedMax);
09b14834 1355
1356 newFileName += "_";
1357 newFileName += Form("%.1f",psiMin);
1358 newFileName += "-";
1359 newFileName += Form("%.1f",psiMax);
5c9b4733 1360 newFileName += ".root";
09b14834 1361
5c9b4733 1362 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1363 gHistPN[0]->SetName("gHistPNRaw"); gHistPN[0]->Write();
1364 gHistNP[0]->SetName("gHistNPRaw"); gHistNP[0]->Write();
1365 gHistPP[0]->SetName("gHistPPRaw"); gHistPP[0]->Write();
1366 gHistNN[0]->SetName("gHistNNRaw"); gHistNN[0]->Write();
1367 if(listBFShuffled) {
1368 gHistPN[1]->SetName("gHistPNShuffled"); gHistPN[1]->Write();
1369 gHistNP[1]->SetName("gHistNPShuffled"); gHistNP[1]->Write();
1370 gHistPP[1]->SetName("gHistPPShuffled"); gHistPP[1]->Write();
1371 gHistNN[1]->SetName("gHistNNShuffled"); gHistNN[1]->Write();
1372 }
3ce45c12 1373 if(listBFMixed || (!listBFMixed&&listQA)) {
5c9b4733 1374 gHistPN[2]->SetName("gHistPNMixed"); gHistPN[2]->Write();
1375 gHistNP[2]->SetName("gHistNPMixed"); gHistNP[2]->Write();
1376 gHistPP[2]->SetName("gHistPPMixed"); gHistPP[2]->Write();
1377 gHistNN[2]->SetName("gHistNNMixed"); gHistNN[2]->Write();
1378
1379 gHistPN[3]->SetName("gHistPNCorrelationFunctions"); gHistPN[3]->Write();
1380 gHistNP[3]->SetName("gHistNPCorrelationFunctions"); gHistNP[3]->Write();
1381 gHistPP[3]->SetName("gHistPPCorrelationFunctions"); gHistPP[3]->Write();
1382 gHistNN[3]->SetName("gHistNNCorrelationFunctions"); gHistNN[3]->Write();
1383 }
1384 newFile->Close();
07d0a35c 1385
1386 // some cleaning
1387 for(Int_t i = 0; i < 4; i++){
1388
1389 if(!listBFShuffled && i == 1) continue;
1390 if(!listBFMixed && (i == 2 || i == 3)) continue;
1391
1392 if(gHistPP[i]) delete gHistPP[i];
1393 if(gHistPN[i]) delete gHistPN[i];
1394 if(gHistNP[i]) delete gHistNP[i];
1395 if(gHistNN[i]) delete gHistNN[i];
1396
1397 if(cPN[i]) delete cPN[i];
1398 if(cNP[i]) delete cNP[i];
1399 if(cPP[i]) delete cPP[i];
1400 if(cNN[i]) delete cNN[i];
1401 }
1402
1403 delete hP;
1404 delete hN;
1405 delete hPP;
1406 delete hPN;
1407 delete hNP;
1408 delete hNN;
1409
1410 delete hPMixed;
1411 delete hNMixed;
1412 delete hPPMixed;
1413 delete hPNMixed;
1414 delete hNPMixed;
1415 delete hNNMixed;
1416
1417 delete hPShuffled;
1418 delete hNShuffled;
1419 delete hPPShuffled;
1420 delete hPNShuffled;
1421 delete hNPShuffled;
1422 delete hNNShuffled;
1423
5c9b4733 1424}
1425
1426//____________________________________________________________//
3b8a460b 1427void drawCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1428 const char* gCentralityEstimator = "V0M",
1429 Int_t gBit = 128,
1430 const char* gEventPlaneEstimator = "VZERO",
5c9b4733 1431 Int_t gCentrality = 1,
1432 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
20006629 1433 Double_t vertexZMin = -10.,
1434 Double_t vertexZMax = 10.,
5c9b4733 1435 Double_t ptTriggerMin = -1.,
1436 Double_t ptTriggerMax = -1.,
1437 Double_t ptAssociatedMin = -1.,
3b8a460b 1438 Double_t ptAssociatedMax = -1.,
1439 Bool_t kFit = kFALSE) {
5c9b4733 1440 //Macro that draws the charge dependent correlation functions
1441 //for each centrality bin for the different pT of trigger and
1442 //associated particles
1443 //Author: Panos.Christakoglou@nikhef.nl
1444 TGaxis::SetMaxDigits(3);
1445
1446 //Get the input file
09b14834 1447 /* TString filename = lhcPeriod;
3b8a460b 1448 filename += "/Centrality"; filename += gCentralityEstimator;
1449 filename += "_Bit"; filename += gBit;
1450 filename += "_"; filename += gEventPlaneEstimator;
1451 filename +="/PttFrom";
1452 filename += Form("%.1f",ptTriggerMin); filename += "To";
1453 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1454 filename += Form("%.1f",ptAssociatedMin); filename += "To";
09b14834 1455 filename += Form("%.1f",ptAssociatedMax); */
1456
1457 TString filename = "correlationFunction.Centrality";
5c9b4733 1458 filename += gCentrality; filename += ".Psi";
1459 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1460 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1461 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1462 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
1463 else filename += "All.PttFrom";
1464 filename += Form("%.1f",ptTriggerMin); filename += "To";
1465 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1466 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1467 filename += Form("%.1f",ptAssociatedMax);
09b14834 1468
1469 filename += "_";
1470 filename += Form("%.1f",psiMin);
1471 filename += "-";
1472 filename += Form("%.1f",psiMax);
5c9b4733 1473 filename += ".root";
1474
1475 //Open the file
1476 TFile *f = TFile::Open(filename.Data());
1477 if((!f)||(!f->IsOpen())) {
1478 Printf("The file %s is not found. Aborting...",filename);
3b8a460b 1479 return;
5c9b4733 1480 }
1481 //f->ls();
1482
1483 //Latex
1484 TString centralityLatex = "Centrality: ";
1485 centralityLatex += centralityArray[gCentrality-1];
1486 centralityLatex += "%";
1487
1488 TString psiLatex;
1489 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 1490 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
5c9b4733 1491 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 1492 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
5c9b4733 1493 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 1494 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
5c9b4733 1495 else
3b8a460b 1496 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
5c9b4733 1497
1498 TString pttLatex = Form("%.1f",ptTriggerMin);
1499 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1500 pttLatex += " GeV/c";
1501
1502 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1503 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1504 ptaLatex += " GeV/c";
1505
1506 TLatex *latexInfo1 = new TLatex();
1507 latexInfo1->SetNDC();
1508 latexInfo1->SetTextSize(0.045);
1509 latexInfo1->SetTextColor(1);
1510
1511 TString pngName;
1512
1513 //============================================================//
1514 //Get the +- correlation function
1515 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1516 gHistPN->SetStats(kFALSE);
2b4abbc5 1517 gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
3b8a460b 1518 gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
5c9b4733 1519 gHistPN->GetXaxis()->CenterTitle();
1520 gHistPN->GetXaxis()->SetTitleOffset(1.2);
1521 gHistPN->GetYaxis()->CenterTitle();
1522 gHistPN->GetYaxis()->SetTitleOffset(1.2);
2b4abbc5 1523 gHistPN->GetZaxis()->SetTitleOffset(1.5);
5c9b4733 1524 TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
1525 cPN->SetFillColor(10); cPN->SetHighLightColor(10);
1526 cPN->SetLeftMargin(0.15);
1527 gHistPN->DrawCopy("surf1fb");
2b4abbc5 1528
5c9b4733 1529 gPad->SetTheta(30); // default is 30
1530 gPad->SetPhi(-60); // default is 30
1531 gPad->Update();
1532
3b8a460b 1533 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1534 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1535 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1536 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
5c9b4733 1537
1538 pngName = "CorrelationFunction.Centrality";
1539 pngName += centralityArray[gCentrality-1];
1540 pngName += ".Psi";
1541 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1542 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1543 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1544 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1545 else pngName += "All.PttFrom";
1546 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1547 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1548 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1549 pngName += Form("%.1f",ptAssociatedMax);
1550 pngName += ".PositiveNegative.png";
1551 cPN->SaveAs(pngName.Data());
3b8a460b 1552 if(kFit)
20006629 1553 fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
3b8a460b 1554 ptTriggerMin,ptTriggerMax,
1555 ptAssociatedMin, ptAssociatedMax,gHistPN);
1556
5c9b4733 1557 //============================================================//
1558 //Get the -+ correlation function
1559 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1560 gHistNP->SetStats(kFALSE);
2b4abbc5 1561 gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
3b8a460b 1562 gHistNP->GetXaxis()->SetRangeUser(-1.4,1);
5c9b4733 1563 gHistNP->GetXaxis()->CenterTitle();
1564 gHistNP->GetXaxis()->SetTitleOffset(1.2);
1565 gHistNP->GetYaxis()->CenterTitle();
1566 gHistNP->GetYaxis()->SetTitleOffset(1.2);
2b4abbc5 1567 gHistNP->GetZaxis()->SetTitleOffset(1.5);
5c9b4733 1568 TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
1569 cNP->SetFillColor(10); cNP->SetHighLightColor(10);
1570 cNP->SetLeftMargin(0.15);
1571 gHistNP->DrawCopy("surf1fb");
2b4abbc5 1572
5c9b4733 1573 gPad->SetTheta(30); // default is 30
1574 gPad->SetPhi(-60); // default is 30
1575 gPad->Update();
1576
3b8a460b 1577 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1578 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1579 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1580 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
5c9b4733 1581
1582 pngName = "CorrelationFunction.Centrality";
1583 pngName += centralityArray[gCentrality-1];
1584 pngName += ".Psi";
1585 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1586 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1587 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1588 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1589 else pngName += "All.PttFrom";
1590 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1591 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1592 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1593 pngName += Form("%.1f",ptAssociatedMax);
1594 pngName += ".NegativePositive.png";
1595 cNP->SaveAs(pngName.Data());
1596
3b8a460b 1597 if(kFit)
20006629 1598 fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
3b8a460b 1599 ptTriggerMin,ptTriggerMax,
1600 ptAssociatedMin, ptAssociatedMax,gHistNP);
1601
5c9b4733 1602 //============================================================//
1603 //Get the ++ correlation function
1604 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1605 gHistPP->SetStats(kFALSE);
2b4abbc5 1606 gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
3b8a460b 1607 gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
5c9b4733 1608 gHistPP->GetXaxis()->CenterTitle();
1609 gHistPP->GetXaxis()->SetTitleOffset(1.2);
1610 gHistPP->GetYaxis()->CenterTitle();
1611 gHistPP->GetYaxis()->SetTitleOffset(1.2);
2b4abbc5 1612 gHistPP->GetZaxis()->SetTitleOffset(1.5);
5c9b4733 1613 TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
1614 cPP->SetFillColor(10); cPP->SetHighLightColor(10);
1615 cPP->SetLeftMargin(0.15);
1616 gHistPP->DrawCopy("surf1fb");
2b4abbc5 1617
5c9b4733 1618 gPad->SetTheta(30); // default is 30
1619 gPad->SetPhi(-60); // default is 30
1620 gPad->Update();
1621
3b8a460b 1622 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1623 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1624 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1625 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
5c9b4733 1626
1627 pngName = "CorrelationFunction.Centrality";
1628 pngName += centralityArray[gCentrality-1];
1629 pngName += ".Psi";
1630 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1631 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1632 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1633 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1634 else pngName += "All.PttFrom";
1635 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1636 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1637 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1638 pngName += Form("%.1f",ptAssociatedMax);
1639 pngName += ".PositivePositive.png";
1640 cPP->SaveAs(pngName.Data());
1641
3b8a460b 1642 if(kFit)
20006629 1643 fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
3b8a460b 1644 ptTriggerMin,ptTriggerMax,
1645 ptAssociatedMin, ptAssociatedMax,gHistPP);
1646
5c9b4733 1647 //============================================================//
1648 //Get the -- correlation function
1649 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1650 gHistNN->SetStats(kFALSE);
2b4abbc5 1651 gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
3b8a460b 1652 gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
5c9b4733 1653 gHistNN->GetXaxis()->CenterTitle();
1654 gHistNN->GetXaxis()->SetTitleOffset(1.2);
1655 gHistNN->GetYaxis()->CenterTitle();
1656 gHistNN->GetYaxis()->SetTitleOffset(1.2);
2b4abbc5 1657 gHistNN->GetZaxis()->SetTitleOffset(1.5);
5c9b4733 1658 TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
1659 cNN->SetFillColor(10); cNN->SetHighLightColor(10);
1660 cNN->SetLeftMargin(0.15);
1661 gHistNN->DrawCopy("surf1fb");
2b4abbc5 1662
5c9b4733 1663 gPad->SetTheta(30); // default is 30
1664 gPad->SetPhi(-60); // default is 30
1665 gPad->Update();
1666
3b8a460b 1667 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1668 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1669 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1670 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
5c9b4733 1671
1672 pngName = "CorrelationFunction.Centrality";
1673 pngName += centralityArray[gCentrality-1];
1674 pngName += ".Psi";
1675 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1676 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1677 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1678 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1679 else pngName += "All.PttFrom";
1680 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1681 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1682 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1683 pngName += Form("%.1f",ptAssociatedMax);
1684 pngName += ".NegativeNegative.png";
1685 cNN->SaveAs(pngName.Data());
24ca71a6 1686
3b8a460b 1687 if(kFit)
20006629 1688 fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
3b8a460b 1689 ptTriggerMin,ptTriggerMax,
1690 ptAssociatedMin, ptAssociatedMax,gHistNN);
a38fd7f3 1691}
1692
3b8a460b 1693//____________________________________________________________//
2340240c 1694void drawProjections(Bool_t kProjectInEta = kFALSE,
41f2bc59 1695 Int_t binMin = 1,
1696 Int_t binMax = 80,
3b8a460b 1697 Int_t gCentrality = 1,
1698 Double_t psiMin = -0.5,
1699 Double_t psiMax = 3.5,
1700 Double_t ptTriggerMin = -1.,
1701 Double_t ptTriggerMax = -1.,
1702 Double_t ptAssociatedMin = -1.,
41f2bc59 1703 Double_t ptAssociatedMax = -1.,
2340240c 1704 Bool_t kUseZYAM = kFALSE,
1705 TString eventClass = "Centrality") {
2b4abbc5 1706 //Macro that draws the charge dependent correlation functions PROJECTIONS
1707 //for each centrality bin for the different pT of trigger and
1708 //associated particles
1709 TGaxis::SetMaxDigits(3);
1710
1711 //Get the input file
09b14834 1712 /*TString filename = lhcPeriod;
3b8a460b 1713 filename += "/Centrality"; filename += gCentralityEstimator;
1714 filename += "_Bit"; filename += gBit;
1715 filename += "_"; filename += gEventPlaneEstimator;
1716 filename +="/PttFrom";
1717 filename += Form("%.1f",ptTriggerMin); filename += "To";
1718 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1719 filename += Form("%.1f",ptAssociatedMin); filename += "To";
09b14834 1720 filename += Form("%.1f",ptAssociatedMax); */
1721
2340240c 1722 TString filename = "correlationFunction.";
1723 if(eventClass == "Centrality"){
1724 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1725 filename += ".PsiAll.PttFrom";
1726 }
1727 else if(eventClass == "Multiplicity"){
1728 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1729 filename += ".PsiAll.PttFrom";
1730 }
1731 else{ // "EventPlane" (default)
1732 filename += "Centrality";
1733 filename += gCentrality; filename += ".Psi";
1734 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1735 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1736 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1737 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1738 else filename += "All.PttFrom";
1739 }
2b4abbc5 1740 filename += Form("%.1f",ptTriggerMin); filename += "To";
1741 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1742 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1743 filename += Form("%.1f",ptAssociatedMax);
2340240c 1744 //if(k2pMethod) filename += "_2pMethod";
09b14834 1745
1746 filename += "_";
1747 filename += Form("%.1f",psiMin);
1748 filename += "-";
1749 filename += Form("%.1f",psiMax);
2340240c 1750
2b4abbc5 1751 filename += ".root";
1752
1753 //Open the file
1754 TFile *f = TFile::Open(filename.Data());
1755 if((!f)||(!f->IsOpen())) {
1756 Printf("The file %s is not found. Aborting...",filename);
1757 return listBF;
1758 }
1759 //f->ls();
1760
1761 //Latex
1762 TString centralityLatex = "Centrality: ";
2340240c 1763 if(eventClass == "Centrality")
1764 centralityLatex += Form("%.0f - %.0f",psiMin,psiMax);
1765 else
1766 centralityLatex += centralityArray[gCentrality-1];
1767 centralityLatex += " %";
2b4abbc5 1768
1769 TString psiLatex;
1770 if((psiMin == -0.5)&&(psiMax == 0.5))
3b8a460b 1771 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
2b4abbc5 1772 else if((psiMin == 0.5)&&(psiMax == 1.5))
3b8a460b 1773 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
2b4abbc5 1774 else if((psiMin == 1.5)&&(psiMax == 2.5))
3b8a460b 1775 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
2b4abbc5 1776 else
3b8a460b 1777 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
2b4abbc5 1778
1779 TString pttLatex = Form("%.1f",ptTriggerMin);
1780 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1781 pttLatex += " GeV/c";
1782
1783 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1784 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1785 ptaLatex += " GeV/c";
1786
1787 TLatex *latexInfo1 = new TLatex();
1788 latexInfo1->SetNDC();
1789 latexInfo1->SetTextSize(0.045);
1790 latexInfo1->SetTextColor(1);
1791
1792 TString pngName;
1793
1794 //============================================================//
1795 //Get the +- correlation function
1796 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1797 gHistPN->SetStats(kFALSE);
1798 gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
3b8a460b 1799 gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
2b4abbc5 1800 gHistPN->GetXaxis()->CenterTitle();
1801 gHistPN->GetXaxis()->SetTitleOffset(1.2);
1802 gHistPN->GetYaxis()->CenterTitle();
1803 gHistPN->GetYaxis()->SetTitleOffset(1.2);
1804 gHistPN->GetZaxis()->SetTitleOffset(1.5);
1805 TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
1806 cPN->SetFillColor(10);
1807 cPN->SetHighLightColor(10);
1808 cPN->SetLeftMargin(0.15);
2b4abbc5 1809
1810 //================//
3b8a460b 1811 TH1D* gHistPNprojection = 0x0;
1812 Double_t sum = 0.0;
1813 Double_t gError = 0.0;
1814 Int_t nCounter = 0;
41f2bc59 1815 //projection in delta eta
3b8a460b 1816 if(kProjectInEta) {
1817 gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsX(),gHistPN->GetXaxis()->GetXmin(),gHistPN->GetXaxis()->GetXmax());
1818 for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
1819 sum = 0.; gError = 0.0; nCounter = 0;
41f2bc59 1820 for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
3b8a460b 1821 sum += gHistPN->GetBinContent(iBinX,iBinY);
1822 if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1823 Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
1824 gError += exy*exy;
1825 }
1826 if(nCounter != 0) {
1827 sum /= nCounter;
1828 gError = TMath::Sqrt(gError)/nCounter;
1829 }
1830 gHistPNprojection->SetBinContent(iBinX,sum);
1831 gHistPNprojection->SetBinError(iBinX,gError);
1832 }
1833 gHistPNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
41f2bc59 1834 if(kUseZYAM)
1835 gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
1836 else
1837 gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#eta)");
3b8a460b 1838 gHistPNprojection->GetXaxis()->SetTitle("#Delta#eta");
41f2bc59 1839 }//projection in delta eta
1840 //projection in delta phi
3b8a460b 1841 else {
1842 gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsY(),gHistPN->GetYaxis()->GetXmin(),gHistPN->GetYaxis()->GetXmax());
1843 for(Int_t iBinY = 1; iBinY <= gHistPN->GetNbinsY(); iBinY++) {
1844 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 1845 for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
1846 //for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
3b8a460b 1847 sum += gHistPN->GetBinContent(iBinX,iBinY);
1848 if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1849 Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
1850 gError += exy*exy;
1851 }
1852 if(nCounter != 0) {
1853 sum /= nCounter;
1854 gError = TMath::Sqrt(gError)/nCounter;
1855 }
1856 gHistPNprojection->SetBinContent(iBinY,sum);
1857 gHistPNprojection->SetBinError(iBinY,gError);
1858 }
41f2bc59 1859 if(kUseZYAM)
1860 gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
1861 else
1862 gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#varphi)");
3b8a460b 1863 gHistPNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2b4abbc5 1864 }
41f2bc59 1865
3b8a460b 1866 //ZYAM
41f2bc59 1867 if(kUseZYAM) {
1868 Double_t reference = gHistPNprojection->GetBinContent(gHistPNprojection->GetMinimumBin());
1869 for(Int_t iBinX = 1; iBinX <= gHistPNprojection->GetNbinsX(); iBinX++)
1870 gHistPNprojection->SetBinContent(iBinX,gHistPNprojection->GetBinContent(iBinX) - reference);
1871 }
3b8a460b 1872
2b4abbc5 1873 gHistPNprojection->GetYaxis()->SetTitleOffset(1.5);
3b8a460b 1874 gHistPNprojection->SetMarkerStyle(20);
1875 gHistPNprojection->SetStats(kFALSE);
1876 gHistPNprojection->DrawCopy("E");
2b4abbc5 1877 //=================//
1878
1879 gPad->SetTheta(30); // default is 30
1880 gPad->SetPhi(-60); // default is 30
1881 gPad->Update();
1882
3b8a460b 1883 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1884 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1885 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1886 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
2b4abbc5 1887
2340240c 1888 pngName = filename;
1889 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
1890 else pngName.ReplaceAll(".root","_DeltaPhiProjection");
2b4abbc5 1891 pngName += ".PositiveNegative.png";
1892 cPN->SaveAs(pngName.Data());
1893
1894 //============================================================//
1895 //Get the -+ correlation function
1896 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1897 gHistNP->SetStats(kFALSE);
1898 gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
3b8a460b 1899 gHistNP->GetXaxis()->SetRangeUser(-1.4,1.4);
2b4abbc5 1900 gHistNP->GetXaxis()->CenterTitle();
1901 gHistNP->GetXaxis()->SetTitleOffset(1.2);
1902 gHistNP->GetYaxis()->CenterTitle();
1903 gHistNP->GetYaxis()->SetTitleOffset(1.2);
1904 gHistNP->GetZaxis()->SetTitleOffset(1.5);
1905 TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
1906 cNP->SetFillColor(10); cNP->SetHighLightColor(10);
1907 cNP->SetLeftMargin(0.15);
2b4abbc5 1908
1909 //================//
3b8a460b 1910 TH1D* gHistNPprojection = 0x0;
1911 Double_t sum = 0.0;
1912 Double_t gError = 0.0;
1913 Int_t nCounter = 0;
1914 if(kProjectInEta) {
1915 gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsX(),gHistNP->GetXaxis()->GetXmin(),gHistNP->GetXaxis()->GetXmax());
1916 for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
1917 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 1918 for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
1919 //for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
3b8a460b 1920 sum += gHistNP->GetBinContent(iBinX,iBinY);
1921 if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1922 Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
1923 gError += exy*exy;
1924 }
1925 if(nCounter != 0) {
1926 sum /= nCounter;
1927 gError = TMath::Sqrt(gError)/nCounter;
1928 }
1929 gHistNPprojection->SetBinContent(iBinX,sum);
1930 gHistNPprojection->SetBinError(iBinX,gError);
1931 }
1932 gHistNPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
41f2bc59 1933 if(kUseZYAM)
1934 gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
1935 else
1936 gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#eta)");
3b8a460b 1937 gHistNPprojection->GetXaxis()->SetTitle("#Delta#eta");
2b4abbc5 1938 }
3b8a460b 1939 else {
1940 gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsY(),gHistNP->GetYaxis()->GetXmin(),gHistNP->GetYaxis()->GetXmax());
1941 for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
1942 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 1943 for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
1944 //for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
3b8a460b 1945 sum += gHistNP->GetBinContent(iBinX,iBinY);
1946 if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1947 Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
1948 gError += exy*exy;
1949 }
1950 if(nCounter != 0) {
1951 sum /= nCounter;
1952 gError = TMath::Sqrt(gError)/nCounter;
1953 }
1954 gHistNPprojection->SetBinContent(iBinY,sum);
1955 gHistNPprojection->SetBinError(iBinY,gError);
1956 }
41f2bc59 1957 if(kUseZYAM)
1958 gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
1959 else
1960 gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#varphi)");
3b8a460b 1961 gHistNPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2b4abbc5 1962 }
3b8a460b 1963 //ZYAM
41f2bc59 1964 if(kUseZYAM) {
1965 Double_t reference = gHistNPprojection->GetBinContent(gHistNPprojection->GetMinimumBin());
3b8a460b 1966 for(Int_t iBinX = 1; iBinX <= gHistNPprojection->GetNbinsX(); iBinX++)
1967 gHistNPprojection->SetBinContent(iBinX,gHistNPprojection->GetBinContent(iBinX) - reference);
41f2bc59 1968 }
3b8a460b 1969
2b4abbc5 1970 gHistNPprojection->GetYaxis()->SetTitleOffset(1.5);
3b8a460b 1971 gHistNPprojection->SetMarkerStyle(20);
1972 gHistNPprojection->SetStats(kFALSE);
1973 gHistNPprojection->DrawCopy("E");
2b4abbc5 1974 //================//
1975
1976 gPad->SetTheta(30); // default is 30
1977 gPad->SetPhi(-60); // default is 30
1978 gPad->Update();
1979
3b8a460b 1980 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1981 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1982 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1983 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
2b4abbc5 1984
2340240c 1985 pngName = filename;
1986 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
1987 else pngName.ReplaceAll(".root","_DeltaPhiProjection");
2b4abbc5 1988 pngName += ".NegativePositive.png";
1989 cNP->SaveAs(pngName.Data());
1990
1991 //============================================================//
1992 //Get the ++ correlation function
1993 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1994 gHistPP->SetStats(kFALSE);
1995 gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
3b8a460b 1996 gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
2b4abbc5 1997 gHistPP->GetXaxis()->CenterTitle();
1998 gHistPP->GetXaxis()->SetTitleOffset(1.2);
1999 gHistPP->GetYaxis()->CenterTitle();
2000 gHistPP->GetYaxis()->SetTitleOffset(1.2);
2001 gHistPP->GetZaxis()->SetTitleOffset(1.5);
2002 TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
2003 cPP->SetFillColor(10); cPP->SetHighLightColor(10);
2004 cPP->SetLeftMargin(0.15);
3b8a460b 2005
2b4abbc5 2006 //=================//
2007 TH1D* gHistPPprojection = 0x0;
3b8a460b 2008 Double_t sum = 0.0;
2009 Double_t gError = 0.0;
2010 Int_t nCounter = 0;
2011 if(kProjectInEta) {
2012 gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsX(),gHistPP->GetXaxis()->GetXmin(),gHistPP->GetXaxis()->GetXmax());
2013 for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
2014 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 2015 for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
2016 //for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
3b8a460b 2017 sum += gHistPP->GetBinContent(iBinX,iBinY);
2018 if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2019 Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
2020 gError += exy*exy;
2021 }
2022 if(nCounter != 0) {
2023 sum /= nCounter;
2024 gError = TMath::Sqrt(gError)/nCounter;
2025 }
2026 gHistPPprojection->SetBinContent(iBinX,sum);
2027 gHistPPprojection->SetBinError(iBinX,gError);
2028 }
2029 gHistPPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
41f2bc59 2030 if(kUseZYAM)
2031 gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
2032 else
2033 gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#eta)");
3b8a460b 2034 gHistPPprojection->GetXaxis()->SetTitle("#Delta#eta");
2b4abbc5 2035 }
3b8a460b 2036 else {
2037 gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsY(),gHistPP->GetYaxis()->GetXmin(),gHistPP->GetYaxis()->GetXmax());
2038 for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
2039 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 2040 for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
2041 //for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
3b8a460b 2042 sum += gHistPP->GetBinContent(iBinX,iBinY);
2043 if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2044 Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
2045 gError += exy*exy;
2046 }
2047 if(nCounter != 0) {
2048 sum /= nCounter;
2049 gError = TMath::Sqrt(gError)/nCounter;
2050 }
2051 gHistPPprojection->SetBinContent(iBinY,sum);
2052 gHistPPprojection->SetBinError(iBinY,gError);
2053 }
41f2bc59 2054 if(kUseZYAM)
2055 gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
2056 else
2057 gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#varphi)");
3b8a460b 2058 gHistPPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2b4abbc5 2059 }
3b8a460b 2060 //ZYAM
41f2bc59 2061 if(kUseZYAM) {
2062 Double_t reference = gHistPPprojection->GetBinContent(gHistPPprojection->GetMinimumBin());
2063 for(Int_t iBinX = 1; iBinX <= gHistPPprojection->GetNbinsX(); iBinX++)
2064 gHistPPprojection->SetBinContent(iBinX,gHistPPprojection->GetBinContent(iBinX) - reference);
2065 }
2066
2b4abbc5 2067 gHistPPprojection->GetYaxis()->SetTitleOffset(1.5);
3b8a460b 2068 gHistPPprojection->SetMarkerStyle(20);
2069 gHistPPprojection->SetStats(kFALSE);
2070 gHistPPprojection->DrawCopy("E");
2b4abbc5 2071 //================//
2072
2073 gPad->SetTheta(30); // default is 30
2074 gPad->SetPhi(-60); // default is 30
2075 gPad->Update();
2076
3b8a460b 2077 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
2078 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
2079 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
2080 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
2b4abbc5 2081
2340240c 2082 pngName = filename;
2083 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
2084 else pngName.ReplaceAll(".root","_DeltaPhiProjection");
2b4abbc5 2085 pngName += ".PositivePositive.png";
2086 cPP->SaveAs(pngName.Data());
2087
2088 //============================================================//
2089 //Get the -- correlation function
2090 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
2091 gHistNN->SetStats(kFALSE);
2092 gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
3b8a460b 2093 gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
2b4abbc5 2094 gHistNN->GetXaxis()->CenterTitle();
2095 gHistNN->GetXaxis()->SetTitleOffset(1.2);
2096 gHistNN->GetYaxis()->CenterTitle();
2097 gHistNN->GetYaxis()->SetTitleOffset(1.2);
2098 gHistNN->GetZaxis()->SetTitleOffset(1.5);
2099 TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
2100 cNN->SetFillColor(10); cNN->SetHighLightColor(10);
2101 cNN->SetLeftMargin(0.15);
2b4abbc5 2102
2103 //=================//
3b8a460b 2104 TH1D* gHistNNprojection = 0x0;
2105 Double_t sum = 0.0;
2106 Double_t gError = 0.0;
2107 Int_t nCounter = 0;
2108 if(kProjectInEta) {
2109 gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsX(),gHistNN->GetXaxis()->GetXmin(),gHistNN->GetXaxis()->GetXmax());
2110 for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
2111 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 2112 for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
2113 //for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
3b8a460b 2114 sum += gHistNN->GetBinContent(iBinX,iBinY);
2115 if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2116 Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
2117 gError += exy*exy;
2118 }
2119 if(nCounter != 0) {
2120 sum /= nCounter;
2121 gError = TMath::Sqrt(gError)/nCounter;
2122 }
2123 gHistNNprojection->SetBinContent(iBinX,sum);
2124 gHistNNprojection->SetBinError(iBinX,gError);
2125 }
2126 gHistNNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
41f2bc59 2127 if(kUseZYAM)
2128 gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
2129 else
2130 gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#eta)");
3b8a460b 2131 gHistNNprojection->GetXaxis()->SetTitle("#Delta#eta");
2b4abbc5 2132 }
3b8a460b 2133 else {
2134 gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsY(),gHistNN->GetYaxis()->GetXmin(),gHistNN->GetYaxis()->GetXmax());
2135 for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
2136 sum = 0.; gError = 0.0; nCounter = 0;
09b14834 2137 for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
2138 //for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
3b8a460b 2139 sum += gHistNN->GetBinContent(iBinX,iBinY);
2140 if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2141 Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
2142 gError += exy*exy;
2143 }
2144 if(nCounter != 0) {
2145 sum /= nCounter;
2146 gError = TMath::Sqrt(gError)/nCounter;
2147 }
2148 gHistNNprojection->SetBinContent(iBinY,sum);
2149 gHistNNprojection->SetBinError(iBinY,gError);
2150 }
41f2bc59 2151 if(kUseZYAM)
2152 gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
2153 else
2154 gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#varphi)");
3b8a460b 2155 gHistNNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2b4abbc5 2156 }
3b8a460b 2157 //ZYAM
41f2bc59 2158 if(kUseZYAM) {
2159 Double_t reference = gHistNNprojection->GetBinContent(gHistNNprojection->GetMinimumBin());
2160 for(Int_t iBinX = 1; iBinX <= gHistNNprojection->GetNbinsX(); iBinX++)
2161 gHistNNprojection->SetBinContent(iBinX,gHistNNprojection->GetBinContent(iBinX) - reference);
2162 }
2b4abbc5 2163
2b4abbc5 2164 gHistNNprojection->GetYaxis()->SetTitleOffset(1.5);
3b8a460b 2165 gHistNNprojection->SetMarkerStyle(20);
2166 gHistNNprojection->SetStats(kFALSE);
2167 gHistNNprojection->DrawCopy("E");
2b4abbc5 2168 //=================//
2169
2170 gPad->SetTheta(30); // default is 30
2171 gPad->SetPhi(-60); // default is 30
2172 gPad->Update();
2173
3b8a460b 2174 latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
2175 latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
2176 latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
2177 latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
2b4abbc5 2178
2340240c 2179 pngName = filename;
2180 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
2181 else pngName.ReplaceAll(".root","_DeltaPhiProjection");
2b4abbc5 2182 pngName += ".NegativeNegative.png";
2183 cNN->SaveAs(pngName.Data());
2184
17c69b4e 2185 TString outFileName = filename;
2186 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
2187 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
2188 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
3b8a460b 2189 gHistNPprojection->Write();
2190 gHistPNprojection->Write();
2191 gHistNNprojection->Write();
2192 gHistPPprojection->Write();
2b4abbc5 2193 fProjection->Close();
2b4abbc5 2194}
742af4bd 2195
5c9b4733 2196//____________________________________________________________//
2197void fitCorrelationFunctions(Int_t gCentrality = 1,
2198 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
20006629 2199 Double_t vertexZMin = -10.,Double_t vertexZMax = 10.,
5c9b4733 2200 Double_t ptTriggerMin = -1.,
2201 Double_t ptTriggerMax = -1.,
2202 Double_t ptAssociatedMin = -1.,
2203 Double_t ptAssociatedMax = -1.,
2204 TH2D *gHist) {
07d0a35c 2205
2206 cout<<"FITTING FUNCTION"<<endl;
2207
5c9b4733 2208 //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
2209 //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
2210 //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
2211 //wing structures: [11]*TMath::Power(x,2)
2212 //flow contribution (v1 up to v4): 2.*([12]*TMath::Cos(y) + [13]*TMath::Cos(2.*y) + [14]*TMath::Cos(3.*y) + [15]*TMath::Cos(4.*y))
24ca71a6 2213 TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))+[5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[17])/[9]),2)),[10]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[17])/[9]),2)),[10]))+[11]*TMath::Power(x,2)+2.*[12]*([13]*TMath::Cos(y) + [14]*TMath::Cos(2.*y) + [15]*TMath::Cos(3.*y) + [16]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
5c9b4733 2214 gFitFunction->SetName("gFitFunction");
742af4bd 2215
2216
5c9b4733 2217 //Normalization
742af4bd 2218 gFitFunction->SetParName(0,"N1");
5c9b4733 2219 //near side peak
742af4bd 2220 gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
2221 gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
2222 gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
2223 gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
5c9b4733 2224 //away side ridge
742af4bd 2225 gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
2226 gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
2227 gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
24ca71a6 2228 //longitudinal ridge
742af4bd 2229 gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
2230 gFitFunction->FixParameter(8,0.0);
2231 gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
2232 gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
5c9b4733 2233 //wing structures
742af4bd 2234 gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
5c9b4733 2235 //flow contribution
742af4bd 2236 gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);gFitFunction->SetParLimits(12,0.0,10);
2237 gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);gFitFunction->SetParLimits(13,0.0,10);
2238 gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);gFitFunction->SetParLimits(14,0.0,10);
2239 gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);gFitFunction->SetParLimits(15,0.0,10);
2240 gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);gFitFunction->SetParLimits(16,0.0,10);
2241 gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
2242
2243 // flow parameters
2244 Double_t fNV = 0.;
2245 Double_t fV1 = 0.;
2246 Double_t fV2 = 0.;
2247 Double_t fV3 = 0.;
2248 Double_t fV4 = 0.;
2249
2250 //Fitting the correlation function (first the edges to extract flow)
2251 gHist->Fit("gFitFunction","nm","",1.0,1.6);
2252
2253 fNV += gFitFunction->GetParameter(12);
2254 fV1 += gFitFunction->GetParameter(13);
2255 fV2 += gFitFunction->GetParameter(14);
2256 fV3 += gFitFunction->GetParameter(15);
2257 fV4 += gFitFunction->GetParameter(16);
2258
2259 gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
2260
2261 fNV += gFitFunction->GetParameter(12);
2262 fV1 += gFitFunction->GetParameter(13);
2263 fV2 += gFitFunction->GetParameter(14);
2264 fV3 += gFitFunction->GetParameter(15);
2265 fV4 += gFitFunction->GetParameter(16);
2266
2267 // Now fit the whole with fixed flow
2268 gFitFunction->FixParameter(12,fNV/2.);
2269 gFitFunction->FixParameter(13,fV1/2.);
2270 gFitFunction->FixParameter(14,fV2/2.);
2271 gFitFunction->FixParameter(15,fV3/2.);
2272 gFitFunction->FixParameter(16,fV4/2.);
2273
2274 gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
2275 gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
2276 gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
2277 gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
2278 //gFitFunction->ReleaseParameter(4);gFitFunction->SetParameter(4,1.0);gFitFunction->SetParLimits(4,0.0,2.0);
2279 gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,1.0);//gFitFunction->SetParLimits(5,0.0,10);
2280 gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.5);//gFitFunction->SetParLimits(6,0.0,10);
2281 //gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(7,1.0);gFitFunction->SetParLimits(7,0.0,2.0);
478c95cf 2282 //gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(8,0.05);
2283 //gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(9,0.6);gFitFunction->SetParLimits(9,0.1,10.0);
742af4bd 2284 //gFitFunction->ReleaseParameter(10);gFitFunction->SetParameter(10,1.0);gFitFunction->SetParLimits(10,0.0,2.0);
2285 gFitFunction->ReleaseParameter(17);gFitFunction->SetParameter(17,0.7);gFitFunction->SetParLimits(17,0.0,0.9);
2286
5c9b4733 2287 gHist->Fit("gFitFunction","nm");
2288
742af4bd 2289
5c9b4733 2290 //Cloning the histogram
2291 TString histoName = gHist->GetName(); histoName += "Fit";
2292 TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
2293 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
2294 gHistResidual->SetName("gHistResidual");
2295 gHistResidual->Sumw2();
2296
2297 for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
2298 for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
2299 gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
2300 }
2301 }
2302 gHistResidual->Add(gHistFit,-1);
2303
2304 //Write to output file
2305 TString newFileName = "correlationFunctionFit";
2306 if(histoName.Contains("PN")) newFileName += "PN";
2307 else if(histoName.Contains("NP")) newFileName += "NP";
2308 else if(histoName.Contains("PP")) newFileName += "PP";
2309 else if(histoName.Contains("NN")) newFileName += "NN";
2310 newFileName += ".Centrality";
2311 newFileName += gCentrality; newFileName += ".Psi";
2312 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
2313 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
2314 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
2315 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
2316 else newFileName += "All.PttFrom";
2317 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
2318 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
2319 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
2320 newFileName += Form("%.1f",ptAssociatedMax);
09b14834 2321
2322 newFileName += "_";
2323 newFileName += Form("%.1f",psiMin);
2324 newFileName += "-";
2325 newFileName += Form("%.1f",psiMax);
2326
5c9b4733 2327 newFileName += ".root";
2328 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
5ea3b761 2329 gHist->Write();
5c9b4733 2330 gHistFit->Write();
2331 gHistResidual->Write();
2332 gFitFunction->Write();
2333 newFile->Close();
2334
2335
2336}
24ca71a6 2337
3b8a460b 2338//____________________________________________________________//
2339TH2D *convolute2D(TH2D* h1, TH2D* h2, TString hname) {
3ce45c12 2340 //
2341 // this function does the convolution of 2 eta or phi "efficiencies" in a deta or dphi "efficiency"
2342 // and returns a new histogram which is normalized to the number of bin combinations
2343
2344 // histogram definition
2345 TH2D *hConv = NULL;
1b5c6300 2346 hConv = new TH2D(hname.Data(),hname.Data(), h1->GetNbinsY(),-2,2,h1->GetNbinsX(),-TMath::Pi()/2.,3*TMath::Pi()/2.);
3ce45c12 2347
2348 Double_t x1 = 0.;
2349 Double_t x2 = 0.;
2350 Double_t y1 = 0.;
2351 Double_t y2 = 0.;
2352 Double_t z1 = 0.;
2353 Double_t z2 = 0.;
2354 Double_t dx = 0.;
2355 Double_t dy = 0.;
2356 Double_t dz = 0.;
2357
2358
2359 // convolution
2360 for(Int_t i = 0; i < h1->GetNbinsX(); i ++){
2361 cout<<"CONVOLUTING BIN = "<<i<<" of "<<h1->GetNbinsX()<<endl;
2362 for(Int_t j = 0; j < h1->GetNbinsY(); j ++){
2363 for(Int_t k = 0; k < h2->GetNbinsX(); k ++){
2364 for(Int_t l = 0; l < h2->GetNbinsY(); l ++){
2365
2366 x1 = (Double_t)h1->GetXaxis()->GetBinCenter(i+1);
2367 y1 = (Double_t)h1->GetYaxis()->GetBinCenter(j+1);
2368 x2 = (Double_t)h2->GetXaxis()->GetBinCenter(k+1);
2369 y2 = (Double_t)h2->GetYaxis()->GetBinCenter(l+1);
2370
2371 z1 = (Double_t)h1->GetBinContent(i+1,j+1);
2372 z2 = (Double_t)h2->GetBinContent(k+1,l+1);
2373
1b5c6300 2374 // need the gymnastics to keep the same binning
2375 dx = x1 - x2 - (h1->GetXaxis()->GetBinWidth(1)/2.);
2376 if(gRandom->Gaus() > 0) dy = y1 - y2 + (h1->GetYaxis()->GetBinWidth(1)/2.);
2377 else dy = y1 - y2 - (h1->GetYaxis()->GetBinWidth(1)/2.);
3ce45c12 2378
2379 if(dx>3./2.*TMath::Pi()) dx = dx - 2.*TMath::Pi();
2380 if(dx<-1./2.*TMath::Pi()) dx = 2*TMath::Pi() + dx;
2381
2382 dz = z1*z2;
2383
2384 hConv->Fill(dy,dx,dz);
2385
2386 }
2387 }
2388 }
2389 }
2390
2391 return hConv;
2392
2393}