]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/Chaoticity/AliChaoticity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
forgot init new data member histogram in ctor
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / Chaoticity / AliChaoticity.cxx
CommitLineData
cd12341d 1#include <iostream>
2#include <math.h>
3#include "TChain.h"
4#include "TFile.h"
5#include "TKey.h"
6#include "TObject.h"
7#include "TObjString.h"
8#include "TList.h"
9#include "TTree.h"
10#include "TH1F.h"
11#include "TH1D.h"
12#include "TH2D.h"
13#include "TH3D.h"
14#include "TProfile.h"
15#include "TProfile2D.h"
16#include "TCanvas.h"
17#include "TRandom3.h"
18#include "TF1.h"
19
20#include "AliAnalysisTask.h"
21#include "AliAnalysisManager.h"
22
23
24#include "AliESDEvent.h"
25#include "AliESDInputHandler.h"
26#include "AliESDtrackCuts.h"
27
28#include "AliAODEvent.h"
29#include "AliAODInputHandler.h"
30#include "AliAODMCParticle.h"
31
32#include "AliChaoticity.h"
33
34#define PI 3.1415927
35#define G_Coeff 0.006399 // 2*pi*alpha*M_pion
36
37
38// Author: Dhevan Gangadharan
39
40ClassImp(AliChaoticity)
41
42//________________________________________________________________________
43AliChaoticity::AliChaoticity():
44AliAnalysisTaskSE(),
45 fname(0),
46 fAOD(0x0),
b6e5ec54 47//fESD(0x0),
cd12341d 48 fOutputList(0x0),
49 fPIDResponse(0x0),
50 fEC(0x0),
51 fEvt(0x0),
52 fTempStruct(0x0),
53 fRandomNumber(0x0),
54 fLEGO(kTRUE),
55 fMCcase(kFALSE),
56 fAODcase(kTRUE),
57 fPbPbcase(kTRUE),
58 fPdensityExplicitLoop(kFALSE),
59 fPdensityPairCut(kTRUE),
60 fTabulatePairs(kFALSE),
61 fBfield(0),
62 fMbin(0),
63 fEDbin(0),
64 fMbins(0),
65 fMultLimit(0),
66 fCentBinLowLimit(0),
67 fCentBinHighLimit(1),
68 fEventCounter(0),
69 fEventsToMix(0),
70 fZvertexBins(0),
71 fMultLimits(),
72 fQcut(),
73 fQLowerCut(0),
74 fNormQcutLow(),
75 fNormQcutHigh(),
76 fKupperBound(0),
77 fQupperBound(0),
78 fQupperBoundWeights(0),
79 fKstepT(),
80 fKstepY(),
81 fKmeanT(),
82 fKmeanY(),
83 fKmiddleT(),
84 fKmiddleY(),
85 fQstep(0),
5e3e77d6 86 fQstepWeights(0),
cd12341d 87 fQmean(),
88 fDampStart(0),
89 fDampStep(0),
5e3e77d6 90 fMomResC2(),
91 fMomRes3DTerm1(),
92 fMomRes3DTerm2(),
93 fMomRes3DTerm3(),
94 fMomRes3DTerm4(),
95 fMomRes3DTerm5(),
cd12341d 96 fTPCTOFboundry(0),
97 fTOFboundry(0),
98 fSigmaCutTPC(0),
99 fSigmaCutTOF(0),
100 fMinSepTPCEntrancePhi(0),
101 fMinSepTPCEntranceEta(0),
102 fShareQuality(0),
103 fShareFraction(0),
104 fTrueMassP(0),
105 fTrueMassPi(0),
106 fTrueMassK(0),
107 fTrueMassKs(0),
108 fTrueMassLam(0),
109 fKtbinL(0),
110 fKtbinH(0),
111 fKybinL(0),
112 fKybinH(0),
113 fQobinL(0),
114 fQobinH(0),
115 fQsbinL(0),
116 fQsbinH(0),
117 fQlbinL(0),
118 fQlbinH(0),
119 fDummyB(0),
120 fNormWeight(0x0),
121 fNormWeightErr(0x0),
122 fDefaultsCharMult(),
123 fDefaultsCharSE(),
124 fDefaultsCharME(),
125 fDefaultsInt(),
126 fPairLocationSE(),
127 fPairLocationME(),
128 fTripletSkip1(),
129 fTripletSkip2(),
130 fOtherPairLocation1(),
131 fOtherPairLocation2(),
132 fNormPairSwitch(),
133 fPairSplitCut(),
5e3e77d6 134 fNormPairs(),
135 fFSI2SS(),
136 fFSI2OS(),
137 fFSIOmega0SS(),
138 fFSIOmega0OS()
cd12341d 139{
140 // Default constructor
141 for(Int_t mb=0; mb<fMbins; mb++){
142 for(Int_t edB=0; edB<kEDbins; edB++){
143 for(Int_t c1=0; c1<2; c1++){
144 for(Int_t c2=0; c2<2; c2++){
145 for(Int_t sc=0; sc<kSCLimit2; sc++){
146 for(Int_t term=0; term<2; term++){
147
148 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
149
150 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
151 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
152 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
153 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
154 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
155 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
156
157 }// term_2
158 }// SC_2
159
160 for(Int_t c3=0; c3<2; c3++){
161 for(Int_t sc=0; sc<kSCLimit3; sc++){
162 for(Int_t term=0; term<5; term++){
163
164 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
165 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
166 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
167 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
168 for(Int_t dt=0; dt<kDENtypes; dt++){
169 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
170 }//dt
171
172 }// term_3
173 }// SC_3
174 }//c3
175 }//c2
176 }//c1
177 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
178 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
179 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
180 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
181 }
182 }
183
184 }// ED
185 }// Mbin
186
187}
188//________________________________________________________________________
189AliChaoticity::AliChaoticity(const Char_t *name, Bool_t MCdecision, Bool_t Tabulatedecision, Bool_t PbPbdecision, Int_t lowCentBin, Int_t highCentBin, Bool_t lego)
190: AliAnalysisTaskSE(name),
191 fname(name),
192 fAOD(0x0),
b6e5ec54 193 //fESD(0x0),
cd12341d 194 fOutputList(0x0),
195 fPIDResponse(0x0),
196 fEC(0x0),
197 fEvt(0x0),
198 fTempStruct(0x0),
199 fRandomNumber(0x0),
200 fLEGO(lego),
201 fMCcase(MCdecision),
202 fAODcase(kTRUE),
203 fPbPbcase(PbPbdecision),
204 fPdensityExplicitLoop(kFALSE),
205 fPdensityPairCut(kTRUE),
206 fTabulatePairs(Tabulatedecision),
207 fBfield(0),
208 fMbin(0),
209 fEDbin(0),
210 fMbins(0),
211 fMultLimit(0),
212 fCentBinLowLimit(lowCentBin),
213 fCentBinHighLimit(highCentBin),
214 fEventCounter(0),
215 fEventsToMix(0),
216 fZvertexBins(0),
217 fMultLimits(),
218 fQcut(),
219 fQLowerCut(0),
220 fNormQcutLow(),
221 fNormQcutHigh(),
222 fKupperBound(0),
223 fQupperBound(0),
224 fQupperBoundWeights(0),
225 fKstepT(),
226 fKstepY(),
227 fKmeanT(),
228 fKmeanY(),
229 fKmiddleT(),
230 fKmiddleY(),
231 fQstep(0),
5e3e77d6 232 fQstepWeights(0),
cd12341d 233 fQmean(),
234 fDampStart(0),
235 fDampStep(0),
5e3e77d6 236 fMomResC2(),
237 fMomRes3DTerm1(),
238 fMomRes3DTerm2(),
239 fMomRes3DTerm3(),
240 fMomRes3DTerm4(),
241 fMomRes3DTerm5(),
cd12341d 242 fTPCTOFboundry(0),
243 fTOFboundry(0),
244 fSigmaCutTPC(0),
245 fSigmaCutTOF(0),
246 fMinSepTPCEntrancePhi(0),
247 fMinSepTPCEntranceEta(0),
248 fShareQuality(0),
249 fShareFraction(0),
250 fTrueMassP(0),
251 fTrueMassPi(0),
252 fTrueMassK(0),
253 fTrueMassKs(0),
254 fTrueMassLam(0),
255 fKtbinL(0),
256 fKtbinH(0),
257 fKybinL(0),
258 fKybinH(0),
259 fQobinL(0),
260 fQobinH(0),
261 fQsbinL(0),
262 fQsbinH(0),
263 fQlbinL(0),
264 fQlbinH(0),
265 fDummyB(0),
266 fNormWeight(0x0),
267 fNormWeightErr(0x0),
268 fDefaultsCharMult(),
269 fDefaultsCharSE(),
270 fDefaultsCharME(),
271 fDefaultsInt(),
272 fPairLocationSE(),
273 fPairLocationME(),
274 fTripletSkip1(),
275 fTripletSkip2(),
276 fOtherPairLocation1(),
277 fOtherPairLocation2(),
278 fNormPairSwitch(),
279 fPairSplitCut(),
5e3e77d6 280 fNormPairs(),
281 fFSI2SS(),
282 fFSI2OS(),
283 fFSIOmega0SS(),
284 fFSIOmega0OS()
cd12341d 285{
286 // Main constructor
287 fLEGO = lego;
288 fAODcase=kTRUE;
289 fMCcase=MCdecision;
290 fTabulatePairs=Tabulatedecision;
291 fPbPbcase=PbPbdecision;
292 fPdensityExplicitLoop = kFALSE;
293 fPdensityPairCut = kTRUE;
294 fCentBinLowLimit = lowCentBin;
295 fCentBinHighLimit = highCentBin;
296
297 for(Int_t mb=0; mb<fMbins; mb++){
298 for(Int_t edB=0; edB<kEDbins; edB++){
299 for(Int_t c1=0; c1<2; c1++){
300 for(Int_t c2=0; c2<2; c2++){
301 for(Int_t sc=0; sc<kSCLimit2; sc++){
302 for(Int_t term=0; term<2; term++){
303
304 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
305
306 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
307 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
308 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
309 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
310 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
311 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
312
313 }// term_2
314 }// SC_2
315
316 for(Int_t c3=0; c3<2; c3++){
317 for(Int_t sc=0; sc<kSCLimit3; sc++){
318 for(Int_t term=0; term<5; term++){
319
320 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
321 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
322 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
323 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
324 for(Int_t dt=0; dt<kDENtypes; dt++){
325 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
326 }//dt
327
328 }// term_3
329 }// SC_3
330 }//c3
331 }//c2
332 }//c1
333 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
334 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
335 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
336 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
337 }
338 }
339
340 }// ED
341 }// Mbin
342
343
344 DefineOutput(1, TList::Class());
345}
346//________________________________________________________________________
347AliChaoticity::AliChaoticity(const AliChaoticity &obj)
348 : AliAnalysisTaskSE(obj.fname),
349 fname(obj.fname),
350 fAOD(obj.fAOD),
b6e5ec54 351 //fESD(obj.fESD),
cd12341d 352 fOutputList(obj.fOutputList),
353 fPIDResponse(obj.fPIDResponse),
354 fEC(obj.fEC),
355 fEvt(obj.fEvt),
356 fTempStruct(obj.fTempStruct),
357 fRandomNumber(obj.fRandomNumber),
358 fLEGO(obj.fLEGO),
359 fMCcase(obj.fMCcase),
360 fAODcase(obj.fAODcase),
361 fPbPbcase(obj.fPbPbcase),
362 fPdensityExplicitLoop(obj.fPdensityExplicitLoop),
363 fPdensityPairCut(obj.fPdensityPairCut),
364 fTabulatePairs(obj.fTabulatePairs),
365 fBfield(obj.fBfield),
366 fMbin(obj.fMbin),
367 fEDbin(obj.fEDbin),
368 fMbins(obj.fMbins),
369 fMultLimit(obj.fMultLimit),
370 fCentBinLowLimit(obj.fCentBinLowLimit),
371 fCentBinHighLimit(obj.fCentBinHighLimit),
372 fEventCounter(obj.fEventCounter),
373 fEventsToMix(obj.fEventsToMix),
374 fZvertexBins(obj.fZvertexBins),
375 fMultLimits(),
376 fQcut(),
377 fQLowerCut(obj.fQLowerCut),
378 fNormQcutLow(),
379 fNormQcutHigh(),
380 fKupperBound(obj.fKupperBound),
381 fQupperBound(obj.fQupperBound),
382 fQupperBoundWeights(obj.fQupperBoundWeights),
383 fKstepT(),
384 fKstepY(),
385 fKmeanT(),
386 fKmeanY(),
387 fKmiddleT(),
388 fKmiddleY(),
389 fQstep(obj.fQstep),
5e3e77d6 390 fQstepWeights(obj.fQstepWeights),
cd12341d 391 fQmean(),
392 fDampStart(obj.fDampStart),
393 fDampStep(obj.fDampStep),
5e3e77d6 394 fMomResC2(),
395 fMomRes3DTerm1(),
396 fMomRes3DTerm2(),
397 fMomRes3DTerm3(),
398 fMomRes3DTerm4(),
399 fMomRes3DTerm5(),
cd12341d 400 fTPCTOFboundry(obj.fTPCTOFboundry),
401 fTOFboundry(obj.fTOFboundry),
402 fSigmaCutTPC(obj.fSigmaCutTPC),
403 fSigmaCutTOF(obj.fSigmaCutTOF),
404 fMinSepTPCEntrancePhi(obj.fMinSepTPCEntrancePhi),
405 fMinSepTPCEntranceEta(obj.fMinSepTPCEntranceEta),
406 fShareQuality(obj.fShareQuality),
407 fShareFraction(obj.fShareFraction),
408 fTrueMassP(obj.fTrueMassP),
409 fTrueMassPi(obj.fTrueMassPi),
410 fTrueMassK(obj.fTrueMassK),
411 fTrueMassKs(obj.fTrueMassKs),
412 fTrueMassLam(obj.fTrueMassLam),
413 fKtbinL(obj.fKtbinL),
414 fKtbinH(obj.fKtbinH),
415 fKybinL(obj.fKybinL),
416 fKybinH(obj.fKybinH),
417 fQobinL(obj.fQobinL),
418 fQobinH(obj.fQobinH),
419 fQsbinL(obj.fQsbinL),
420 fQsbinH(obj.fQsbinH),
421 fQlbinL(obj.fQlbinL),
422 fQlbinH(obj.fQlbinH),
423 fDummyB(obj.fDummyB),
424 fNormWeight(),
425 fNormWeightErr(),
426 fDefaultsCharMult(),
427 fDefaultsCharSE(),
428 fDefaultsCharME(),
429 fDefaultsInt(),
430 fPairLocationSE(),
431 fPairLocationME(),
432 fTripletSkip1(),
433 fTripletSkip2(),
434 fOtherPairLocation1(),
435 fOtherPairLocation2(),
436 fNormPairSwitch(),
437 fPairSplitCut(),
5e3e77d6 438 fNormPairs(),
439 fFSI2SS(),
440 fFSI2OS(),
441 fFSIOmega0SS(),
442 fFSIOmega0OS()
cd12341d 443{
444 // Copy constructor
cd12341d 445}
446//________________________________________________________________________
447AliChaoticity &AliChaoticity::operator=(const AliChaoticity &obj)
448{
449 // Assignment operator
450 if (this == &obj)
451 return *this;
452
453 fname = obj.fname;
454 fAOD = obj.fAOD;
b6e5ec54 455 //fESD = obj.fESD;
cd12341d 456 fOutputList = obj.fOutputList;
457 fPIDResponse = obj.fPIDResponse;
458 fEC = obj.fEC;
459 fEvt = obj.fEvt;
460 fTempStruct = obj.fTempStruct;
461 fRandomNumber = obj.fRandomNumber;
462 fLEGO = fLEGO;
463 fMCcase = obj.fMCcase;
464 fAODcase = obj.fAODcase;
465 fPbPbcase = obj.fPbPbcase;
466 fPdensityExplicitLoop = obj.fPdensityExplicitLoop;
467 fPdensityPairCut = obj.fPdensityPairCut;
468 fTabulatePairs = obj.fTabulatePairs;
469 fBfield = obj.fBfield;
470 fMbin = obj.fMbin;
471 fEDbin = obj.fEDbin;
472 fMbins = obj.fMbins;
473 fMultLimit = obj.fMultLimit;
474 fCentBinLowLimit = obj.fCentBinLowLimit;
475 fCentBinHighLimit = obj.fCentBinHighLimit;
476 fEventCounter = obj.fEventCounter;
477 fEventsToMix = obj.fEventsToMix;
478 fZvertexBins = obj.fZvertexBins;
479 //fMultLimits = ;
480 //fQcut = ;
481 fQLowerCut = obj.fQLowerCut;
482 //fNormQcutLow = ;
483 //fNormQcutHigh = ;
484 fKupperBound = obj.fKupperBound;
485 fQupperBound = obj.fQupperBound;
486 fQupperBoundWeights = obj.fQupperBoundWeights;
487 //fKstepT = ;
488 //fKstepY = ;
489 //fKmeanT = ;
490 //fKmeanY = ;
491 //fKmiddleT = ;
492 //fKmiddleY = ;
493 fQstep = obj.fQstep;
5e3e77d6 494 fQstepWeights = obj.fQstepWeights;
cd12341d 495 //fQmean = ;
496 fDampStart = obj.fDampStart;
497 fDampStep = obj.fDampStep;
cd12341d 498 fTPCTOFboundry = obj.fTPCTOFboundry;
499 fTOFboundry = obj.fTOFboundry;
500 fSigmaCutTPC = obj.fSigmaCutTPC;
501 fSigmaCutTOF = obj.fSigmaCutTOF;
502 fMinSepTPCEntrancePhi = obj.fMinSepTPCEntrancePhi;
503 fMinSepTPCEntranceEta = obj.fMinSepTPCEntranceEta;
504 fShareQuality = obj.fShareQuality;
505 fShareFraction = obj.fShareFraction;
506 fTrueMassP = obj.fTrueMassP;
507 fTrueMassPi = obj.fTrueMassPi;
508 fTrueMassK = obj.fTrueMassK;
509 fTrueMassKs = obj.fTrueMassKs;
510 fTrueMassLam = obj.fTrueMassLam;
511 fKtbinL = obj.fKtbinL;
512 fKtbinH = obj.fKtbinH;
513 fKybinL = obj.fKybinL;
514 fKybinH = obj.fKybinH;
515 fQobinL = obj.fQobinL;
516 fQobinH = obj.fQobinH;
517 fQsbinL = obj.fQsbinL;
518 fQsbinH = obj.fQsbinH;
519 fQlbinL = obj.fQlbinL;
520 fQlbinH = obj.fQlbinH;
521 fDummyB = obj.fDummyB;
522 fNormWeight = obj.fNormWeight;
523 fNormWeightErr = obj.fNormWeightErr;
524
525 return (*this);
526}
527//________________________________________________________________________
528AliChaoticity::~AliChaoticity()
529{
530 // Destructor
531 if(fAOD) delete fAOD;
b6e5ec54 532 //if(fESD) delete fESD;
cd12341d 533 if(fOutputList) delete fOutputList;
534 if(fPIDResponse) delete fPIDResponse;
535 if(fEC) delete fEC;
536 if(fEvt) delete fEvt;
68b2d0b9 537 if(fTempStruct) delete [] fTempStruct;
cd12341d 538 if(fRandomNumber) delete fRandomNumber;
5e3e77d6 539 if(fFSI2SS) delete fFSI2SS;
540 if(fFSI2OS) delete fFSI2OS;
541 if(fFSIOmega0SS) delete fFSIOmega0SS;
542 if(fFSIOmega0OS) delete fFSIOmega0OS;
543 if(fMomResC2) delete fMomResC2;
544 if(fMomRes3DTerm1) delete fMomRes3DTerm1;
545 if(fMomRes3DTerm2) delete fMomRes3DTerm2;
546 if(fMomRes3DTerm3) delete fMomRes3DTerm3;
547 if(fMomRes3DTerm4) delete fMomRes3DTerm4;
548 if(fMomRes3DTerm5) delete fMomRes3DTerm5;
549
68b2d0b9 550 for(int i=0; i<fMultLimit; i++){
551 if(fPairLocationSE[i]) delete [] fPairLocationSE[i];
552 if(fPairLocationME[i]) delete [] fPairLocationME[i];
553 for(int j=0; j<2; j++){
554 if(fOtherPairLocation1[j][i]) delete [] fOtherPairLocation1[j][i];
555 if(fOtherPairLocation2[j][i]) delete [] fOtherPairLocation2[j][i];
556 }
557 for(int j=0; j<3; j++) if(fNormPairSwitch[j][i]) delete [] fNormPairSwitch[j][i];
558 for(int j=0; j<4; j++) if(fPairSplitCut[j][i]) delete [] fPairSplitCut[j][i];
559 }
560 for(int i=0; i<kPairLimit; i++) if(fTripletSkip1[i]) delete [] fTripletSkip1[i];
561 for(int i=0; i<2*kPairLimit; i++) if(fTripletSkip2[i]) delete [] fTripletSkip2[i];
562 for(int i=0; i<3; i++) if(fNormPairs[i]) delete [] fNormPairs[i];
cd12341d 563 //
564 for(Int_t mb=0; mb<fMbins; mb++){
565 for(Int_t edB=0; edB<kEDbins; edB++){
566 for(Int_t c1=0; c1<2; c1++){
567 for(Int_t c2=0; c2<2; c2++){
568 for(Int_t sc=0; sc<kSCLimit2; sc++){
569 for(Int_t term=0; term<2; term++){
570
571 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2;
572
573 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal;
574 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared;
575 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL;
576 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW;
577 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL;
578 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW;
579
580 }// term_2
581 }// SC_2
582
583 for(Int_t c3=0; c3<2; c3++){
584 for(Int_t sc=0; sc<kSCLimit3; sc++){
585 for(Int_t term=0; term<5; term++){
586
587 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3;
588 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3;
589 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3;
590 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3;
591 for(Int_t dt=0; dt<kDENtypes; dt++){
592 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm;
593 }//dt
594
595 }// term_3
596 }// SC_3
597 }//c3
598 }//c2
599 }//c1
600 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
601 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
602 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD;
603 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD;
604 }
605 }
606
607 }// ED
608 }// Mbin
609
610
611}
612//________________________________________________________________________
613void AliChaoticity::ParInit()
614{
615 cout<<"AliChaoticity MyInit() call"<<endl;
616
617 fRandomNumber = new TRandom3();
618 fRandomNumber->SetSeed(0);
cd12341d 619
620 //
621 fEventCounter=0;
622 if(fPdensityExplicitLoop) fEventsToMix=3;
623 else if(fPdensityPairCut && !fPdensityExplicitLoop) fEventsToMix=2;
624 else fEventsToMix=0;
625 fZvertexBins=2;//2
626
627 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
628 fTOFboundry = 2.1;// TOF pid used below this momentum
629 fSigmaCutTPC = 2.0;// 2.0 (norm). 1.5 Sys Deviation
630 fSigmaCutTOF = 2.0;// 2.0 (norm). 1.5 Sys Deviation
631
632 ////////////////////////////////////////////////
633 // Proton Pair Cuts
634 fMinSepTPCEntrancePhi = 0.035;// 0.035(norm). 0.04 Sys Deviation
635 fMinSepTPCEntranceEta = 0.035;// 0.035(norm). 0.04 Sys Deviation
636 fShareQuality = .5;// max
637 fShareFraction = .05;// max
638 ////////////////////////////////////////////////
5e3e77d6 639
cd12341d 640
641 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
642 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
643
644
645
646 if(fPbPbcase) {// PbPb
647 fMultLimit=kMultLimitPbPb;
648 fMbins=kCentBins;
649 fQcut[0]=0.1;
650 fQcut[1]=0.1;
651 fQcut[2]=0.6;
652 fNormQcutLow[0] = 0.15;//1.06 (test also at 0.15)
653 fNormQcutHigh[0] = 0.175;//1.1 (test also at 0.175)
654 fNormQcutLow[1] = 1.34;//1.34
655 fNormQcutHigh[1] = 1.4;//1.4
656 fNormQcutLow[2] = 1.1;//1.1
657 fNormQcutHigh[2] = 1.4;//1.4
658 }
659 else {// pp
660 fMultLimit=kMultLimitpp;
661 fMbins=kMultBinspp;
662 fQcut[0]=0.6;
663 fQcut[1]=0.6;
664 fQcut[2]=0.6;
665 fNormQcutLow[0] = 1.0;
666 fNormQcutHigh[0] = 1.5;
667 fNormQcutLow[1] = 1.0;
668 fNormQcutHigh[1] = 1.5;
669 fNormQcutLow[2] = 1.0;
670 fNormQcutHigh[2] = 1.5;
671 }
672
5e3e77d6 673 fQLowerCut = 0.005;// was 0.002
cd12341d 674 fKupperBound = 1.0;
675 //
676 // 4x1 (Kt: 0-0.2, 0.2-0.4, 0.4-0.6, 0.6-1.0)
677 //fKstepT[0] = 0.2; fKstepT[1] = 0.2; fKstepT[2] = 0.2; fKstepT[3] = 0.4;
678 //fKstepY[0] = 1.6;// -0.8 to +0.8
679 //fKmeanT[0] = 0.1; fKmeanT[1] = 0.3; fKmeanT[2] = 0.5; fKmeanT[3] = 0.8;
680 //fKmeanY[0] = 0;// central y
681 //fKmiddleT[0] = 0.1; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.5; fKmiddleT[3] = 0.8;
682 //fKmiddleY[0] = 0;
683 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
684 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
685 fKstepY[0] = 1.6;// -0.8 to +0.8
686 fKmeanT[0] = 0.241; fKmeanT[1] = 0.369; fKmeanT[2] = 0.573;
687 fKmeanY[0] = 0;// central y
688 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
689 fKmiddleY[0] = 0;
690 // 2x1 (Kt: 0-0.4, 0.4-1.0)
691 //fKstepT[0] = 0.4; fKstepT[1] = 0.6;
692 //fKstepY[0] = 1.6;// -0.8 to +0.8
693 //fKmeanT[0] = 0.255; fKmeanT[1] = 0.505;
694 //fKmiddleT[0] = 0.2; fKmiddleT[1] = 0.7;
695 //fKmeanY[0] = 0;// central y
696 //fKmiddleY[0] = 0;
697 // 1x1 (Kt: 0-1.0)
698 //fKstepT[0] = 1.0;
699 //fKstepY[0] = 1.6;// -0.8 to +0.8
700 //fKmeanT[0] = 0.306;
701 //fKmiddleT[0] = 0.5;
702 //fKmeanY[0] = 0;// central y
703 //fKmiddleY[0] = 0;
704 //
705 fQupperBoundWeights = 0.2;
706 fQupperBound = 0.1;
5e3e77d6 707 fQstep = fQupperBound/Float_t(kQbins);
708 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
709 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
cd12341d 710 //
711 fDampStart = 0.3;
712 fDampStep = 0.02;
713
714
715
716 fEC = new AliChaoticityEventCollection **[fZvertexBins];
717 for(UShort_t i=0; i<fZvertexBins; i++){
718
719 fEC[i] = new AliChaoticityEventCollection *[fMbins];
720
721 for(UShort_t j=0; j<fMbins; j++){
722
723 fEC[i][j] = new AliChaoticityEventCollection(fEventsToMix+1, fMultLimit, kPairLimit, fMCcase);
724 }
725 }
726
727
728 for(Int_t i=0; i<fMultLimit; i++) fDefaultsCharMult[i]='0';
729 for(Int_t i=0; i<kPairLimit; i++) fDefaultsCharSE[i]='0';
730 for(Int_t i=0; i<2*kPairLimit; i++) fDefaultsCharME[i]='0';
731 for(Int_t i=0; i<fMultLimit; i++) fDefaultsInt[i]=-1;
732 for(Int_t i=0; i<fMultLimit; i++) fPairLocationSE[i] = new TArrayI(fMultLimit,fDefaultsInt);
733 for(Int_t i=0; i<fMultLimit; i++) fPairLocationME[i] = new TArrayI(fMultLimit,fDefaultsInt);
734 for(Int_t i=0; i<kPairLimit; i++) fTripletSkip1[i] = new TArrayC(fMultLimit,fDefaultsCharSE);
735 for(Int_t i=0; i<2*kPairLimit; i++) fTripletSkip2[i] = new TArrayC(fMultLimit,fDefaultsCharME);
736
737
738 // Normalization utilities
739 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
740 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
741 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
742 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
743 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
744 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
745 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
746
747 // Track Merging/Splitting utilities
748 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P11
749 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P12
750 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P13
751 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[3][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P23
752
753
754 fNormPairs[0] = new AliChaoticityNormPairStruct[kNormPairLimit];
755 fNormPairs[1] = new AliChaoticityNormPairStruct[kNormPairLimit];
756
757
758 fTempStruct = new AliChaoticityTrackStruct[fMultLimit];
759
760
761 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
762
763
764 // Initialize Weight Array
765 fNormWeight = new Float_t******[fMbins];// fMbin
766 fNormWeightErr = new Float_t******[fMbins];// fMbin
767 for(Int_t i=0; i<fMbins; i++){// Mbin iterator
768 fNormWeight[i] = new Float_t*****[kEDbins];// create ED bins
769 fNormWeightErr[i] = new Float_t*****[kEDbins];// create ED bins
770
771 for(Int_t j=0; j<kEDbins; j++){// ED iterator
772 fNormWeight[i][j] = new Float_t****[kKbinsT];// create Kt bins
773 fNormWeightErr[i][j] = new Float_t****[kKbinsT];// create Kt bins
774
775 for(Int_t k=0; k<kKbinsT; k++){// Kt iterator
776 fNormWeight[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
777 fNormWeightErr[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
778
779 for(Int_t l=0; l<kKbinsY; l++){// Ky iterator
780 fNormWeight[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
781 fNormWeightErr[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
782
783 for(Int_t m=0; m<kQbinsWeights; m++){// Qout iterator
784 fNormWeight[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
785 fNormWeightErr[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
786
787 for(Int_t n=0; n<kQbinsWeights; n++){// Qside iterator
788 fNormWeight[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
789 fNormWeightErr[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
790 }
791 }
792 }
793 }
794 }
795 }
796
797
798
799 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
800 if(!fLEGO && !fTabulatePairs) {
5e3e77d6 801 SetWeightArrays(fLEGO);// Set Weight Array
802 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
803 if(!fMCcase) SetMomResCorrections(fLEGO);// Read Momentum resolution file
cd12341d 804 }
805
806}
807//________________________________________________________________________
808void AliChaoticity::UserCreateOutputObjects()
809{
810 // Create histograms
811 // Called once
812
813 ParInit();// Initialize my settings
5e3e77d6 814
cd12341d 815
816 fOutputList = new TList();
817 fOutputList->SetOwner();
818
819 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
820 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
821 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
822 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
823 fOutputList->Add(fVertexDist);
824
825
826 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
827 fOutputList->Add(fDCAxyDistPlus);
828 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
829 fOutputList->Add(fDCAzDistPlus);
830 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
831 fOutputList->Add(fDCAxyDistMinus);
832 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
833 fOutputList->Add(fDCAzDistMinus);
834
835
836 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
837 fOutputList->Add(fEvents1);
838 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
839 fOutputList->Add(fEvents2);
840
841 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
842 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
843 fOutputList->Add(fMultDist1);
844
845 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
846 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
847 fOutputList->Add(fMultDist2);
848
849 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
850 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
851 fOutputList->Add(fMultDist3);
852
853 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
854 fOutputList->Add(fPtEtaDist);
855
856 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
857 fOutputList->Add(fPhiPtDist);
858
859 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
860 fOutputList->Add(fTOFResponse);
861 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
862 fOutputList->Add(fTPCResponse);
863
864 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
865 fOutputList->Add(fRejectedPairs);
866 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
867 fOutputList->Add(fRejectedEvents);
868
869 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
870 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
871 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
872 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
873
874 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
875 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
876 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
877 if(fMCcase) fOutputList->Add(fAllOSPairs);
878
879 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
880 fOutputList->Add(fAvgMult);
881
882 TH3D *fTPNRejects = new TH3D("fTPNRejects","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
883 fOutputList->Add(fTPNRejects);
884
885 TH2D *fKt3Dist = new TH2D("fKt3Dist","",fMbins,.5,fMbins+.5, 10,0,1);
886 fOutputList->Add(fKt3Dist);
887
888
889
890 if(fPdensityExplicitLoop || fPdensityPairCut){
891
892 for(Int_t mb=0; mb<fMbins; mb++){
893 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
894
895 for(Int_t edB=0; edB<kEDbins; edB++){
896 for(Int_t c1=0; c1<2; c1++){
897 for(Int_t c2=0; c2<2; c2++){
898 for(Int_t sc=0; sc<kSCLimit2; sc++){
899 for(Int_t term=0; term<2; term++){
900
901 TString *nameEx2 = new TString("Explicit2_Charge1_");
902 *nameEx2 += c1;
903 nameEx2->Append("_Charge2_");
904 *nameEx2 += c2;
905 nameEx2->Append("_SC_");
906 *nameEx2 += sc;
907 nameEx2->Append("_M_");
908 *nameEx2 += mb;
909 nameEx2->Append("_ED_");
910 *nameEx2 += edB;
911 nameEx2->Append("_Term_");
912 *nameEx2 += term+1;
913
914 if(sc==0 || sc==3 || sc==5){
915 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
916 }
917
918 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
919 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
5e3e77d6 920 TString *nameEx2QW=new TString(nameEx2->Data());
921 nameEx2QW->Append("_QW");
922 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
923 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW);
cd12341d 924 // Momentum resolution histos
925 if(fMCcase && sc==0){
926 TString *nameIdeal = new TString(nameEx2->Data());
927 nameIdeal->Append("_Ideal");
928 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",kDENtypes,-0.5,kDENtypes-0.5, kQbinsWeights,0,fQupperBoundWeights);
929 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
930 TString *nameSmeared = new TString(nameEx2->Data());
931 nameSmeared->Append("_Smeared");
932 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",kDENtypes,-0.5,kDENtypes-0.5, kQbinsWeights,0,fQupperBoundWeights);
933 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
934 }
935 if(c1==c2 && sc==0){
936
937 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
938 nameEx2OSLB1->Append("_osl_b1");
939 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
940 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL);
941 nameEx2OSLB1->Append("_QW");
942 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
943 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW);
944 //
945 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
946 nameEx2OSLB2->Append("_osl_b2");
947 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
948 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL);
949 nameEx2OSLB2->Append("_QW");
950 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
951 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW);
952
953 }
954
955 }// term_2
956 }// SC_2
957
958 // skip 3-particle if Tabulate6DPairs is true
959 if(fTabulatePairs) continue;
960
961 for(Int_t c3=0; c3<2; c3++){
962 for(Int_t sc=0; sc<kSCLimit3; sc++){
963 for(Int_t term=0; term<5; term++){
964 TString *nameEx3 = new TString("Explicit3_Charge1_");
965 *nameEx3 += c1;
966 nameEx3->Append("_Charge2_");
967 *nameEx3 += c2;
968 nameEx3->Append("_Charge3_");
969 *nameEx3 += c3;
970 nameEx3->Append("_SC_");
971 *nameEx3 += sc;
972 nameEx3->Append("_M_");
973 *nameEx3 += mb;
974 nameEx3->Append("_ED_");
975 *nameEx3 += edB;
976 nameEx3->Append("_Term_");
977 *nameEx3 += term+1;
978
979 TString *namePC3 = new TString("PairCut3_Charge1_");
980 *namePC3 += c1;
981 namePC3->Append("_Charge2_");
982 *namePC3 += c2;
983 namePC3->Append("_Charge3_");
984 *namePC3 += c3;
985 namePC3->Append("_SC_");
986 *namePC3 += sc;
987 namePC3->Append("_M_");
988 *namePC3 += mb;
989 namePC3->Append("_ED_");
990 *namePC3 += edB;
991 namePC3->Append("_Term_");
992 *namePC3 += term+1;
993
994 ///////////////////////////////////////
995 // skip degenerate histograms
996 if(sc==0 || sc==6 || sc==9){// Identical species
997 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
998 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
999 }else if(sc!=5){
1000 if( (c1+c2)==1) {if(c1!=0) continue;}
1001 }else {}// do nothing for pi-k-p case
1002
1003 /////////////////////////////////////////
1004
1005 if(fPdensityExplicitLoop){
1006 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = new TH1D(nameEx3->Data(),"Three Particle Distribution",200,0,2);
1007 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3);
1008 //
1009 nameEx3->Append("_Norm");
1010 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = new TH1D(nameEx3->Data(),"Explicit_3 Norm",1,-0.5,0.5);
1011 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3);
1012 }
1013 if(fPdensityPairCut){
1014 TString *nameNorm=new TString(namePC3->Data());
1015 nameNorm->Append("_Norm");
1016 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1017 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1018 //
1019 if(sc<=2){
1020 TString *name3DQ=new TString(namePC3->Data());
1021 name3DQ->Append("_3D");
1022 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH3D(name3DQ->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1023 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1024 //
5e3e77d6 1025 if(c1==c2 && c1==c3 && sc==0 && fMCcase==kFALSE){
1026 TString *name4vectCC3=new TString(namePC3->Data());
1027 name4vectCC3->Append("_4VectProdCC3");
1028 // use 3.75e6 MeV^4 as the resolution on QprodSum
1029 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC3 = new TH2D(name4vectCC3->Data(),"",200,0,0.0003, 200,0,0.0003);
1030 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC3);
1031 if(term!=0){
1032 TString *name4vectCC2=new TString(namePC3->Data());
1033 name4vectCC2->Append("_4VectProdCC2");
1034 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC2 = new TH2D(name4vectCC2->Data(),"",200,0,0.0003, 200,0,0.0003);
1035 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC2);
1036 }
1037 if(term==4){
1038 TString *name4vectCC0=new TString(namePC3->Data());
1039 name4vectCC0->Append("_4VectProdCC0");
1040 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC0 = new TH2D(name4vectCC0->Data(),"",200,0,0.0003, 200,0,0.0003);
1041 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProdTermsCC0);
1042 }
1043 }
1044 if(sc==0 && fMCcase==kTRUE){
1045 TString *name3DMomResIdeal=new TString(namePC3->Data());
1046 name3DMomResIdeal->Append("_Ideal");
1047 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH3D(name3DMomResIdeal->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1048 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1049 TString *name3DMomResSmeared=new TString(namePC3->Data());
1050 name3DMomResSmeared->Append("_Smeared");
1051 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH3D(name3DMomResSmeared->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1052 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared);
1053 }
1054
1055 //
cd12341d 1056 if(c1==c2 && c1==c3 && term==4 && sc==0 && fMCcase==kFALSE){
1057 for(Int_t dt=0; dt<kDENtypes; dt++){
1058 TString *nameDenType=new TString("PairCut3_Charge1_");
1059 *nameDenType += c1;
1060 nameDenType->Append("_Charge2_");
1061 *nameDenType += c2;
1062 nameDenType->Append("_Charge3_");
1063 *nameDenType += c3;
1064 nameDenType->Append("_SC_");
1065 *nameDenType += sc;
1066 nameDenType->Append("_M_");
1067 *nameDenType += mb;
1068 nameDenType->Append("_ED_");
1069 *nameDenType += edB;
1070 nameDenType->Append("_TPN_");
1071 *nameDenType += dt;
1072
1073 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1074 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm);
1075 // neglect errors for TPN
1076 //nameDenType->Append("_Err");
1077 //Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1078 //fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr);
5e3e77d6 1079 //
1080 TString *name4vectTPN=new TString(nameDenType->Data());
1081 name4vectTPN->Append("_4VectProd");
1082 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProdTwoPartNorm = new TH2D(name4vectTPN->Data(),"",200,0,0.0003, 200,0,0.0003);
1083 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProdTwoPartNorm);
1084
cd12341d 1085 }
1086
1087 }// term=4
1088 }// c and sc exclusion
1089 }// PdensityPairCut
1090 }// term_3
1091 }// SC_3
1092 }//c3
1093 }//c2
1094 }//c1
1095 }// ED
1096 }// mbin
1097 }// Pdensity Method
1098
1099
1100 if(fTabulatePairs){
1101
1102 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
1103 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
1104 for(Int_t mb=0; mb<fMbins; mb++){
1105 for(Int_t edB=0; edB<kEDbins; edB++){
1106
1107 TString *nameNum = new TString("TwoPart_num_Kt_");
1108 *nameNum += tKbin;
1109 nameNum->Append("_Ky_");
1110 *nameNum += yKbin;
1111 nameNum->Append("_M_");
1112 *nameNum += mb;
1113 nameNum->Append("_ED_");
1114 *nameNum += edB;
1115
1116 TString *nameDen = new TString("TwoPart_den_Kt_");
1117 *nameDen += tKbin;
1118 nameDen->Append("_Ky_");
1119 *nameDen += yKbin;
1120 nameDen->Append("_M_");
1121 *nameDen += mb;
1122 nameDen->Append("_ED_");
1123 *nameDen += edB;
1124
1125
1126 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = new TH3I(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1127 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD);
1128
1129 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = new TH3I(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1130 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD);
1131 }
1132 }
1133 }
1134 }
1135
1136 }
1137
1138
1139 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1140 fOutputList->Add(fQsmearMean);
1141 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1142 fOutputList->Add(fQsmearSq);
1143 TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
1144 fOutputList->Add(fQDist);
1145
1146
1147 ////////////////////////////////////
1148 ///////////////////////////////////
1149
1150 PostData(1, fOutputList);
1151}
1152
1153//________________________________________________________________________
1154void AliChaoticity::Exec(Option_t *)
1155{
1156 // Main loop
1157 // Called for each event
1ccd6f0d 1158 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
cd12341d 1159 fEventCounter++;
1160
b6e5ec54 1161 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
cd12341d 1162
b6e5ec54 1163 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1164 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
cd12341d 1165
cd12341d 1166
1167 // Trigger Cut
1168 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1169 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1ccd6f0d 1170 if(!isSelected1 && !fMCcase) {/*cout<<"Event Rejected"<<endl;*/ return;}
5e3e77d6 1171 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
cd12341d 1172 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1173 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1ccd6f0d 1174 if(!isSelected1 && !isSelected2 && !fMCcase) {/*cout<<"Event Rejected"<<endl;*/ return;}
1175 }else {/*cout<<"Event Rejected"<<endl;*/ return;}
cd12341d 1176
1177 ///////////////////////////////////////////////////////////
1178 const AliAODVertex *primaryVertexAOD;
1179 AliCentrality *centrality;// for AODs and ESDs
1180
1181
1182 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1183 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1184 fPIDResponse = inputHandler->GetPIDResponse();
1185
1186
1187 TClonesArray *mcArray = 0x0;
1188 if(fMCcase){
1189 if(fAODcase){
1190 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1191 if(!mcArray || mcArray->GetEntriesFast() >= 110000){
1192 cout<<"No MC particle branch found or Array too large!!"<<endl;
1193 cout<<mcArray->GetEntriesFast()<<endl;
1194 return;
1195 }
1196 }
1197 }
1198
1199
1200 UInt_t status=0;
1201 Int_t positiveTracks=0, negativeTracks=0;
1202 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1203
1204 Double_t vertex[3]={0};
1205 Int_t zbin=0;
1206 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1207 /////////////////////////////////////////////////
1208
1209
1210 Float_t centralityPercentile=0;
1211 Float_t cStep=5.0, cStart=0;
1212
1213
1214 if(fAODcase){// AOD case
1215
1216 if(fPbPbcase){
1217 centrality = fAOD->GetCentrality();
1218 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1219 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
5e3e77d6 1220 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
cd12341d 1221 cout<<"Centrality % = "<<centralityPercentile<<endl;
1222 }
1223
1224
1225
1226
1227 ////////////////////////////////
1228 // Vertexing
1229 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1230 primaryVertexAOD = fAOD->GetPrimaryVertex();
1231 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1232
1233 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1234 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1235
1236 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1237 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1238
1239 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1240
1241 fBfield = fAOD->GetMagneticField();
1242
1243 for(Int_t i=0; i<fZvertexBins; i++){
1244 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1245 zbin=i;
1246 break;
1247 }
1248 }
1249
1250
1251
1252 /////////////////////////////
1253 // Create Shuffled index list
1254 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1255 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1256 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1257 /////////////////////////////
1258
1259 // Track loop
1260 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1261 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1262 if (!aodtrack) continue;
1263 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1264
1265 status=aodtrack->GetStatus();
1266
1267 if((Bool_t)(((1<<7) & aodtrack->GetFilterMap()) == 0)) continue;// AOD filterBit cut
1268 // 1<<5 is for "standard cuts with tight dca cut"
1269 // 1<<7 is for TPC only tracking
1270 if(aodtrack->Pt() < 0.16) continue;
1271 if(fabs(aodtrack->Eta()) > 0.8) continue;
1272
1273
1274 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1275 if(!goodMomentum) continue;
1276 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1277
1278 Float_t dca2[2];
1279 Float_t dca3d;
1280
1281 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1282 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1283 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1284
1285 fTempStruct[myTracks].fStatus = status;
1286 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1287 fTempStruct[myTracks].fId = aodtrack->GetID();
1288 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1289 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1290 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1291 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1292 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1293 fTempStruct[myTracks].fEta = aodtrack->Eta();
1294 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1295 fTempStruct[myTracks].fDCAXY = dca2[0];
1296 fTempStruct[myTracks].fDCAZ = dca2[1];
1297 fTempStruct[myTracks].fDCA = dca3d;
1298 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1299 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1300
1301
1302
1303 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1304 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1305 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1306
1307 if(fTempStruct[myTracks].fCharge==+1) {
1308 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1309 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1310 }else {
1311 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1312 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1313 }
1314
1315 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1316 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1317
1318
1319 // nSimga PID workaround
1320 fTempStruct[myTracks].fElectron = kFALSE;
1321 fTempStruct[myTracks].fPion = kFALSE;
1322 fTempStruct[myTracks].fKaon = kFALSE;
1323 fTempStruct[myTracks].fProton = kFALSE;
1324
1325 Float_t nSigmaTPC[5];
1326 Float_t nSigmaTOF[5];
1327 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1328 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1329 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1330 Float_t signalTPC=0, signalTOF=0;
1331 Double_t integratedTimesTOF[10]={0};
1332 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1333 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1334 if (!aodTrack2) continue;
1335 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1336
1337 UInt_t status2=aodTrack2->GetStatus();
1338
1339 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1340 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1341 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1342 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1343 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1344 //
1345 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1346 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1347 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1348 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1349 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1350 signalTPC = aodTrack2->GetTPCsignal();
1351
1352 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1353 fTempStruct[myTracks].fTOFhit = kTRUE;
1354 signalTOF = aodTrack2->GetTOFsignal();
1355 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1356 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1357
1358 }// aodTrack2
5e3e77d6 1359
cd12341d 1360 ///////////////////
1361 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1362 if(fTempStruct[myTracks].fTOFhit) {
1363 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1364 }
1365 ///////////////////
1366
1367 // Use TOF if good hit and above threshold
1368 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1369 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1370 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1371 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1372 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1373 }else {// TPC info instead
1374 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1375 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1376 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1377 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1378 }
1379
1380
1381 //////////////////////////////////////
1382 // Bayesian PIDs for TPC only tracking
1383 //const Double_t* PID = aodtrack->PID();
1384 //fTempStruct[myTracks].fElectron = kFALSE;
1385 //fTempStruct[myTracks].Pion = kFALSE;
1386 //fTempStruct[myTracks].Kaon = kFALSE;
1387 //fTempStruct[myTracks].Proton = kFALSE;
1388
1389 // Pions
1390 //if(PID[2] > 0.2) fTempStruct[myTracks].Pion = kTRUE;// pions: 0.2 --> 0.5
1391 //
1392 // Kaons
1393 //if(PID[3] <= 0.5) fTempStruct[myTracks].Kaon = kFALSE;// kaons
1394 //
1395 // Protons
1396 //if(PID[4] <= 0.5) fTempStruct[myTracks].Proton = kFALSE;// protons
1397 //////////////////////////////////////
1398
1399
1400 // Ensure there is only 1 candidate per track
1401 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1402 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1403 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1404 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1405 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1406 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1407 ////////////////////////
1408 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1409
1410 if(!fTempStruct[myTracks].fPion) continue;// only pions
1411
1412
1413
1414
1415 if(fTempStruct[myTracks].fPion) {// pions
1416 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1417 fTempStruct[myTracks].fKey = 1;
1418 }else if(fTempStruct[myTracks].fKaon){// kaons
1419 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1420 fTempStruct[myTracks].fKey = 10;
1421 }else{// protons
1422 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1423 fTempStruct[myTracks].fKey = 100;
1424 }
1425
1426
1427
1428 if(aodtrack->Charge() > 0) positiveTracks++;
1429 else negativeTracks++;
1430
1431 if(fTempStruct[myTracks].fPion) pionCount++;
1432 if(fTempStruct[myTracks].fKaon) kaonCount++;
1433 if(fTempStruct[myTracks].fProton) protonCount++;
1434
1435 myTracks++;
1436
1437 }
1438 }else {// ESD tracks
1439 cout<<"ESDs not supported currently"<<endl;
1440 return;
1441 }
1442
1443
1444 if(myTracks >= 1) {
1445 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1446 }
1447
1448
b6e5ec54 1449 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1450 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
cd12341d 1451
1452 /////////////////////////////////////////
1453 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1454 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1455 /////////////////////////////////////////
1456
1457
1458 ////////////////////////////////
1459 ///////////////////////////////
1460 // Mbin determination
1461 //
1462 // Mbin set to Pion Count Only for pp!!!!!!!
1463 fMbin=-1;
1464 if(!fPbPbcase){
1465 for(Int_t i=0; i<kMultBinspp; i++){
1466 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1467 // Mbin 0 has 1 pion
1468 }
1469 }else{
1470 for(Int_t i=0; i<kCentBins; i++){
1471 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1472 fMbin=i;// 0 = most central
1473 break;
1474 }
1475 }
1476 }
1477
1478 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1479
1480 //////////////////////////////////////////////////
1481 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1482 //////////////////////////////////////////////////
1483
1484
1485
1486 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1487 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1488
1489 ////////////////////////////////////
1490 // Add event to buffer if > 0 tracks
1491 if(myTracks > 0){
1492 fEC[zbin][fMbin]->FIFOShift();
1493 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1494 (fEvt)->fNtracks = myTracks;
1495 (fEvt)->fFillStatus = 1;
1496 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1497 if(fMCcase){
1498 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1499 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1500 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1501 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1502 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1503 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1504 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1505 }
1506 }
1507 }
1508
1509
1510
1511 Float_t qinv12=0, qinv13=0, qinv23=0;
5e3e77d6 1512 Int_t qinv12Bin=0, qinv13Bin=0, qinv23Bin=0;
cd12341d 1513 Float_t qout=0, qside=0, qlong=0;
1514 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1515 Float_t transK12=0, rapK12=0, transK3=0;
1516 Int_t transKbin=0, rapKbin=0;
1517 Float_t q3=0;
1518 Int_t ch1=0, ch2=0, ch3=0;
1519 Short_t key1=0, key2=0, key3=0;
1520 Int_t bin1=0, bin2=0, bin3=0;
1521 Float_t pVect1[4]={0};
1522 Float_t pVect2[4]={0};
1523 Float_t pVect3[4]={0};
1524 Float_t pVect1MC[4]={0};
1525 Float_t pVect2MC[4]={0};
5e3e77d6 1526 Float_t pVect3MC[4]={0};
cd12341d 1527 Int_t index1=0, index2=0, index3=0;
1528 Float_t weight12=0, weight13=0, weight23=0;
1529 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1530 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1531 Float_t weightTotal=0;//, weightTotalErr=0;
1532 //
5e3e77d6 1533 Float_t qinv12MC=0, qinv13MC=0, qinv23MC=0;
cd12341d 1534 AliAODMCParticle *mcParticle1=0x0;
1535 AliAODMCParticle *mcParticle2=0x0;
1536
1537
1538 if(fPdensityPairCut){
1539 ////////////////////
1540 Int_t pairCountSE=0, pairCountME=0;
1541 Int_t normPairCount[2]={0};
1542 Int_t numOtherPairs1[2][fMultLimit];
1543 Int_t numOtherPairs2[2][fMultLimit];
1544 Bool_t exitCode=kFALSE;
1545 Int_t tempNormFillCount[2][2][2][10][5];
1546
1547
1548 // reset to defaults
1549 for(Int_t i=0; i<fMultLimit; i++) {
1550 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1551 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1552
1553 // Normalization Utilities
1554 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1555 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1556 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1557 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1558 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1559 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1560 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1561 numOtherPairs1[0][i]=0;
1562 numOtherPairs1[1][i]=0;
1563 numOtherPairs2[0][i]=0;
1564 numOtherPairs2[1][i]=0;
1565
1566 // Track Merging/Splitting Utilities
1567 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1568 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1569 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1570 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1571 }
1572
1573 // Reset the temp Normalization counters
1574 for(Int_t i=0; i<2; i++){// Charge1
1575 for(Int_t j=0; j<2; j++){// Charge2
1576 for(Int_t k=0; k<2; k++){// Charge3
1577 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1578 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1579 tempNormFillCount[i][j][k][l][m] = 0;
1580 }
1581 }
1582 }
1583 }
1584 }
1585
1586
1587 ///////////////////////////////////////////////////////
1588 // Start the pairing process
1589 // P11 pairing
1590 // 1st Particle
1591
1592 for (Int_t i=0; i<myTracks; i++) {
1593
1594 Int_t en2=0;
1595
1596 // 2nd particle
1597 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1598
1599 key1 = (fEvt)->fTracks[i].fKey;
1600 key2 = (fEvt+en2)->fTracks[j].fKey;
1601 Short_t fillIndex2 = FillIndex2part(key1+key2);
1602 Short_t qCutBin = SetQcutBin(fillIndex2);
1603 Short_t normBin = SetNormBin(fillIndex2);
1604 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1605 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1606 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1607 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1608 //
1609
1610 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1611 GetQosl(pVect1, pVect2, qout, qside, qlong);
1612 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1613
1614 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1615
1616 ///////////////////////////////
1617 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1618 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1619 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1620
1621 if(fMCcase && ch1==ch2 && fMbin==0){
1622 for(Int_t rstep=0; rstep<10; rstep++){
1623 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1624 // propagate through B field to r=1.2m
1625 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1626 if(phi1 > 2*PI) phi1 -= 2*PI;
1627 if(phi1 < 0) phi1 += 2*PI;
1628 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1629 if(phi2 > 2*PI) phi2 -= 2*PI;
1630 if(phi2 < 0) phi2 += 2*PI;
1631 Float_t deltaphi = phi1 - phi2;
1632 if(deltaphi > PI) deltaphi -= PI;
1633 if(deltaphi < -PI) deltaphi += PI;
1634 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1635 }
1636 }
5e3e77d6 1637
cd12341d 1638 // Pair Splitting/Merging cut
1639 if(ch1 == ch2){
1640 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1641 fPairSplitCut[0][i]->AddAt('1',j);
1642 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1643 continue;
1644 }
1645 }
5e3e77d6 1646
cd12341d 1647 // HIJING tests
1648 if(fMCcase && fillIndex2==0){
1649
1650 // Check that label does not exceed stack size
1651 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1652 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1653 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1654 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1655 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1656 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1657 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1658 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1659 if(qinv12<0.1 && ch1==ch2) {
1660 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1661 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1662 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
1663 }
1664
1665
5e3e77d6 1666 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm + 1 Therminator setting
cd12341d 1667 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
1668 Int_t denIndex = rIter*kNDampValues + myDampIt;
cd12341d 1669 Float_t WInput = MCWeight(ch1,ch2, rIter, myDampIt, qinv12MC);
1670 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
1671 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
1672 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
1673 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
1674 }
1675 }
1676
1677 //HIJING resonance test
1678 if(ch1 != ch2){
1679 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
1680 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
1681 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
1682 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
1683 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
1684 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
1685 }
1686 }
1687 }
1688 }// label check 2
1689 }// MC case
1690
1691 //////////////////////////////////////////
1692 // 2-particle term
1693 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 1694 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
1695
cd12341d 1696 // osl frame
1697 if(bin1==bin2 && fillIndex2==0){
1698 if((transK12 > 0.2) && (transK12 < 0.3)){
1699 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1700 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1701 }
1702 if((transK12 > 0.6) && (transK12 < 0.7)){
1703 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1704 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1705 }
1706 }
5e3e77d6 1707
cd12341d 1708 //////////////////////////////////////////
1709 if(fTabulatePairs){
1710 if(fillIndex2==0 && bin1==bin2){
1711 rapK12 = 0;
1712 transKbin=-1; rapKbin=-1;
1713 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1714 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1715 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1716 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1717 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1718 continue;
1719 }
1720 }
1721
1722
1723 // exit out of loop if there are too many pairs
1724 if(pairCountSE >= kPairLimit) {cout<<"Too many SE pairs"<<endl; exitCode=kTRUE; continue;}
1725 if(exitCode) continue;
5e3e77d6 1726
cd12341d 1727 //////////////////////////
1728 // Enforce the Qcut
1729 if(qinv12 <= fQcut[qCutBin]) {
5e3e77d6 1730
cd12341d 1731 ///////////////////////////
1732 // particle 1
1733 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
1734 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
1735 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
1736 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
1737 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
1738 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
1739 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
1740 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
1741 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1742 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1743 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1744 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
5e3e77d6 1745 }
cd12341d 1746 // particle 2
1747 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1748 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1749 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1750 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1751 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1752 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
1753 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
1754 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1755 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1756 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1757 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1758 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1759 }
5e3e77d6 1760
cd12341d 1761 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
1762
1763 fPairLocationSE[i]->AddAt(pairCountSE,j);
1764
1765 pairCountSE++;
1766
1767 }
1768
1769
1770 /////////////////////////////////////////////////////////
1771 // Normalization Region
1772
1773 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1774 // particle 1
1775 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1776 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1777 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1778 // particle 2
1779 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1780 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1781 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1782
1783
1784 //other past pairs with particle j
1785 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
1786 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
1787 if(locationOtherPair < 0) continue;// no pair there
1788 Int_t indexOther1 = i;
1789 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
1790
1791 // Both possible orderings of other indexes
1792 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
1793
1794 // 1 and 2 are from SE
1795 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
1796 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
1797 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1798 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
1799
1800 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
1801 }
1802
1803 }// pastpair P11 loop
1804
1805
1806 fNormPairSwitch[en2][i]->AddAt('1',j);
1807 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1808 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1809
1810 numOtherPairs1[en2][i]++;
1811 numOtherPairs2[en2][j]++;
1812
1813
1814 normPairCount[en2]++;
1815 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1816
1817 }// Norm Region
1818
1819 }// j particle
1820 }// i particle
1821
1822
1823
1824 //////////////////////////////////////////////
1825 // P12 pairing
1826 // 1st Particle
1827 for (Int_t i=0; i<myTracks; i++) {
1828
1829 Int_t en2=1;
1830
1831 // 2nd particle
1832 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1833
1834 key1 = (fEvt)->fTracks[i].fKey;
1835 key2 = (fEvt+en2)->fTracks[j].fKey;
1836 Short_t fillIndex2 = FillIndex2part(key1+key2);
1837 Short_t qCutBin = SetQcutBin(fillIndex2);
1838 Short_t normBin = SetNormBin(fillIndex2);
1839 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1840 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1841 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1842 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1843
1844 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1845 GetQosl(pVect1, pVect2, qout, qside, qlong);
1846 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1847
1848 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
1849
1850 ///////////////////////////////
1851 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1852 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1853 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1854
1855 if(fMCcase && ch1==ch2 && fMbin==0){
1856 for(Int_t rstep=0; rstep<10; rstep++){
1857 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1858 // propagate through B field to r=1.2m
1859 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1860 if(phi1 > 2*PI) phi1 -= 2*PI;
1861 if(phi1 < 0) phi1 += 2*PI;
1862 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1863 if(phi2 > 2*PI) phi2 -= 2*PI;
1864 if(phi2 < 0) phi2 += 2*PI;
1865 Float_t deltaphi = phi1 - phi2;
1866 if(deltaphi > PI) deltaphi -= PI;
1867 if(deltaphi < -PI) deltaphi += PI;
1868 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1869 }
1870 }
1871
1872 if(ch1 == ch2){
1873 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1874 fPairSplitCut[1][i]->AddAt('1',j);
1875 continue;
1876 }
1877 }
1878
1879 //////////////////////////////////////////
1880 // 2-particle term
1881 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 1882 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
1883
cd12341d 1884 // osl frame
1885 if(bin1==bin2 && fillIndex2==0){
1886 if((transK12 > 0.2) && (transK12 < 0.3)){
1887 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1888 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1889 }
1890 if((transK12 > 0.6) && (transK12 < 0.7)){
1891 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1892 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1893 }
1894 }
1895 //////////////////////////////////////////
1896 if(fTabulatePairs){
1897 if(fillIndex2==0 && bin1==bin2){
1898 rapK12 = 0;
1899 transKbin=-1; rapKbin=-1;
1900 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1901 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1902 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1903 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1904 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1905 continue;
1906 }
1907 }
1908
1909
1910 if(pairCountME >= 2*kPairLimit) {cout<<"Too many ME pairs"<<endl; exitCode=kTRUE; continue;}
1911 if(exitCode) continue;
1912
1913 if(qinv12 <= fQcut[qCutBin]) {
1914 ///////////////////////////
1915
1916 // particle 1
1917 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
1918 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
1919 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
1920 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
1921 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
1922 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
1923 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
1924 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
1925 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1926 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1927 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1928 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
1929 }
1930 // particle 2
1931 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1932 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1933 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1934 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1935 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1936 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
1937 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
1938 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1939 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1940 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1941 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1942 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1943 }
1944
1945 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
1946
1947 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
1948
1949 pairCountME++;
1950
1951 }
1952
1953 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1954 // particle 1
1955 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1956 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1957 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1958 // particle 2
1959 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1960 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1961 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1962
1963 //other past pairs in P11 with particle i
1964 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
1965 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
1966 if(locationOtherPairP11 < 0) continue;// no pair there
1967 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
1968
1969 //Check other past pairs in P12
1970 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
1971
1972 // 1 and 3 are from SE
1973 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
1974 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
1975 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1976 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
1977 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
1978
1979
1980 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
1981 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
1982 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
1983
1984
1985 }// P11 loop
1986
1987
1988 fNormPairSwitch[en2][i]->AddAt('1',j);
1989 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1990 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1991
1992 numOtherPairs1[en2][i]++;
1993 numOtherPairs2[en2][j]++;
1994
1995 normPairCount[en2]++;
1996 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1997
1998 }// Norm Region
1999
2000
2001 }
2002 }
2003
2004
2005 ///////////////////////////////////////
2006 // P13 pairing (just for Norm counting of term5)
2007 for (Int_t i=0; i<myTracks; i++) {
2008
2009 // exit out of loop if there are too many pairs
2010 // dont bother with this loop if exitCode is set.
2011 if(exitCode) break;
2012
2013 // 2nd particle
2014 Int_t en2=2;
2015
2016 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2017
2018 key1 = (fEvt)->fTracks[i].fKey;
2019 key2 = (fEvt+en2)->fTracks[j].fKey;
2020 Short_t fillIndex2 = FillIndex2part(key1+key2);
2021 Short_t normBin = SetNormBin(fillIndex2);
2022 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2023 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2024 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2025 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2026
2027 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2028
2029 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2030
2031 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2032 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2033
2034 if(ch1 == ch2){
2035 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2036 fPairSplitCut[2][i]->AddAt('1',j);
2037 continue;
2038 }
2039 }
2040
2041 /////////////////////////////////////////////////////////
2042 // Normalization Region
2043
2044 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2045
2046 fNormPairSwitch[en2][i]->AddAt('1',j);
2047
2048 }// Norm Region
2049 }
2050 }
2051
2052
2053
2054 ///////////////////////////////////////
2055 // P23 pairing (just for Norm counting of term5)
2056 Int_t en1=1;
2057 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
2058
2059 // exit out of loop if there are too many pairs
2060 // dont bother with this loop if exitCode is set.
2061 if(exitCode) break;
2062
2063 // 2nd event
2064 Int_t en2=2;
2065 // 2nd particle
2066 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2067
2068 if(exitCode) break;
2069
2070 key1 = (fEvt+en1)->fTracks[i].fKey;
2071 key2 = (fEvt+en2)->fTracks[j].fKey;
2072 Short_t fillIndex2 = FillIndex2part(key1+key2);
2073 Short_t normBin = SetNormBin(fillIndex2);
2074 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2075 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2076 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2077 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2078
2079 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2080
2081 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2082
2083 ///////////////////////////////
2084 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2085 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2086
2087 if(ch1 == ch2){
2088 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2089 fPairSplitCut[3][i]->AddAt('1',j);
2090 continue;
2091 }
2092 }
2093
2094 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2095
2096 Int_t index1P23 = i;
2097 Int_t index2P23 = j;
2098
2099 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2100 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2101 if(locationOtherPairP12 < 0) continue; // no pair there
2102 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2103
2104
2105 //Check other past pair status in P13
2106 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2107
2108 // all from different event
2109 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2110 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2111 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2112 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2113
2114 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2115 }
2116 }
2117 }
2118
2119
2120
2121
2122 ///////////////////////////////////////////////////
2123 // Do not use pairs from events with too many pairs
2124 if(exitCode) {
2125 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2126 (fEvt)->fNpairsSE = 0;
2127 (fEvt)->fNpairsME = 0;
2128 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2129 return;// Skip event
2130 }else{
2131 (fEvt)->fNpairsSE = pairCountSE;
2132 (fEvt)->fNpairsME = pairCountME;
2133 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2134 }
2135 ///////////////////////////////////////////////////
2136
2137
2138
2139 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
b6e5ec54 2140 //cout<<"Start Main analysis"<<endl;
cd12341d 2141
2142 ///////////////////////////////////////////////////////////////////////
2143 ///////////////////////////////////////////////////////////////////////
2144 ///////////////////////////////////////////////////////////////////////
2145 //
2146 //
5e3e77d6 2147 // Start the Main Correlation Analysis
cd12341d 2148 //
2149 //
2150 ///////////////////////////////////////////////////////////////////////
2151
2152
2153 /////////////////////////////////////////////////////////
2154 // Skip 3-particle part if Tabulate6DPairs is set to true
2155 if(fTabulatePairs) return;
2156 /////////////////////////////////////////////////////////
2157
2158 // Set the Normalization counters
2159 for(Int_t termN=0; termN<5; termN++){
2160
2161 if(termN==0){
2162 if((fEvt)->fNtracks ==0) continue;
2163 }else if(termN<4){
2164 if((fEvt)->fNtracks ==0) continue;
2165 if((fEvt+1)->fNtracks ==0) continue;
2166 }else {
2167 if((fEvt)->fNtracks ==0) continue;
2168 if((fEvt+1)->fNtracks ==0) continue;
2169 if((fEvt+2)->fNtracks ==0) continue;
2170 }
2171
2172 for(Int_t sc=0; sc<kSCLimit3; sc++){
2173
2174 for(Int_t c1=0; c1<2; c1++){
2175 for(Int_t c2=0; c2<2; c2++){
2176 for(Int_t c3=0; c3<2; c3++){
2177
2178 if(sc==0 || sc==6 || sc==9){// Identical species
2179 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2180 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2181 }else if(sc!=5){
2182 if( (c1+c2)==1) {if(c1!=0) continue;}
2183 }else {}// do nothing for pi-k-p case
2184 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2185 }
2186 }
2187 }
2188 }
2189 }
2190
2191
2192
2193 /////////////////////////////////////////////
2194 // Calculate Pair-Cut Correlations
2195 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2196
2197 Int_t nump1=0;
2198 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2199 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2200
2201 // 1st pair
2202 for(Int_t p1=0; p1<nump1; p1++){
2203
2204 if(en1case==0){
2205 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2206 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2207 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2208 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2209 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2210 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2211 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2212 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2213 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2214 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2215 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2216 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2217 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2218
2219 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2220 }
2221 if(en1case==1){
2222 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2223 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2224 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2225 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2226 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2227 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2228 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2229 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2230 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2231 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2232 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2233 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2234 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2235
2236 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2237 }
2238
cd12341d 2239
2240 // en2 buffer
2241 for(Int_t en2=0; en2<3; en2++){
2242 //////////////////////////////////////
2243
2244 Bool_t skipcase=kTRUE;
2245 Short_t config=-1, part=-1;
2246 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2247 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2248 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2249 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2250
2251 if(skipcase) continue;
2252
2253
2254 // 3-particle terms
2255 // 3rd particle
2256 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2257 index3 = k;
2258
2259
2260 // remove auto-correlations and duplicate triplets
2261 if(config==1){
2262 if( index1 == index3) continue;
2263 if( index2 == index3) continue;
2264 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2265
2266 // skip the switched off triplets
2267 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2268 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2269 continue;
2270 }
2271 ///////////////////////////////
2272 // Turn off 1st possible degenerate triplet
2273 if(index1 < index3){// verify correct id ordering ( index1 < k )
2274 if(fPairLocationSE[index1]->At(index3) >= 0){
2275 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2276 }
2277 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2278 }else {// or k < index1
2279 if(fPairLocationSE[index3]->At(index1) >= 0){
2280 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2281 }
2282 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2283 }
2284 // turn off 2nd possible degenerate triplet
2285 if(index2 < index3){// verify correct id ordering (index2 < k)
2286 if(fPairLocationSE[index2]->At(index3) >= 0){
2287 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2288 }
2289 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2290 }else {// or k < index2
2291 if(fPairLocationSE[index3]->At(index2) >= 0){
2292 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2293 }
2294 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2295 }
2296
2297 }// end config 1
2298
2299 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2300 ///////////////////////////////
2301 // Turn off 1st possible degenerate triplet
2302 if(fPairLocationME[index1]->At(index3) >= 0){
2303 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2304 }
2305
2306 // turn off 2nd possible degenerate triplet
2307 if(fPairLocationME[index2]->At(index3) >= 0){
2308 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2309 }
2310
2311 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2312 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2313 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2314 }// end config 2 part 1
2315
2316 if(config==2 && part==2){// P12T1
2317 if( index1 == index3) continue;
2318
2319 // skip the switched off triplets
2320 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2321 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2322 continue;
2323 }
2324 // turn off another possible degenerate
2325 if(fPairLocationME[index3]->At(index2) >= 0){
2326 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2327 }// end config 2 part 2
2328
2329 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2330 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2331 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2332 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2333 }
2334 if(config==3){// P12T3
2335 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2336 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2337 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2338 }// end config 3
2339
2340
5e3e77d6 2341
cd12341d 2342 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2343 key3 = (fEvt+en2)->fTracks[k].fKey;
2344 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2345 Short_t fillIndex13 = FillIndex2part(key1+key3);
2346 Short_t fillIndex23 = FillIndex2part(key2+key3);
2347 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2348 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2349 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2350 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2351 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2352 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
5e3e77d6 2353 if(fMCcase){
2354 pVect3MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPx;
2355 pVect3MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPy;
2356 pVect3MC[3] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPz;
2357 pVect3MC[0] = sqrt(pow(pVect3MC[1],2)+pow(pVect3MC[2],2)+pow(pVect3MC[3],2)+pow(fTrueMassPi,2));
2358 qinv12MC = GetQinv(0, pVect1MC, pVect2MC);
2359 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
2360 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
2361 }
cd12341d 2362 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2363 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2364
5e3e77d6 2365
cd12341d 2366 if(qinv13 < fQLowerCut) continue;
2367 if(qinv23 < fQLowerCut) continue;
5e3e77d6 2368 if(qinv13 > fQcut[qCutBin13]) continue;// cut value was 3x higher before
2369 if(qinv23 > fQcut[qCutBin23]) continue;// cut value was 3x higher before
2370 // if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
2371
cd12341d 2372 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2373 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2374 Float_t firstQ=0, secondQ=0, thirdQ=0;
5e3e77d6 2375 //
2376 if(!fMCcase){
2377 qinv12Bin = fMomRes3DTerm1->GetXaxis()->FindBin(qinv12);
2378 qinv13Bin = fMomRes3DTerm1->GetYaxis()->FindBin(qinv13);
2379 qinv23Bin = fMomRes3DTerm1->GetZaxis()->FindBin(qinv23);
2380 }
2381 //4-vector product sums
2382 Double_t Qsum = (pVect1[0]-pVect2[0])*(pVect2[1]-pVect3[1]) - (pVect1[1]-pVect2[1])*(pVect2[0]-pVect3[0]);
2383 Qsum += (pVect1[0]-pVect2[0])*(pVect2[2]-pVect3[2]) - (pVect1[2]-pVect2[2])*(pVect2[0]-pVect3[0]);
2384 Qsum += (pVect1[0]-pVect2[0])*(pVect2[3]-pVect3[3]) - (pVect1[3]-pVect2[3])*(pVect2[0]-pVect3[0]);
2385 Qsum += (pVect1[1]-pVect2[1])*(pVect2[2]-pVect3[2]) - (pVect1[2]-pVect2[2])*(pVect2[1]-pVect3[1]);
2386 Qsum += (pVect1[1]-pVect2[1])*(pVect2[3]-pVect3[3]) - (pVect1[3]-pVect2[3])*(pVect2[1]-pVect3[1]);
2387 Qsum += (pVect1[2]-pVect2[2])*(pVect2[3]-pVect3[3]) - (pVect1[3]-pVect2[3])*(pVect2[2]-pVect3[2]);
2388 Qsum *= Qsum;
2389 // 4-vector inner product test
2390 Double_t QsumIn = (pVect1[0]-pVect2[0])*(pVect2[0]-pVect3[0]) - (pVect1[1]-pVect2[1])*(pVect2[1]-pVect3[1]);
2391 QsumIn += -(pVect1[2]-pVect2[2])*(pVect2[2]-pVect3[2]) - (pVect1[3]-pVect2[3])*(pVect2[3]-pVect3[3]);
2392 QsumIn *= QsumIn;
2393
2394 //
cd12341d 2395 if(config==1) {// 123
2396 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2397
2398 if(fillIndex3 <= 2){
2399 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
2400 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ);
2401 ((TH2F*)fOutputList->FindObject("fKt3Dist"))->Fill(fMbin+1, transK3);
5e3e77d6 2402 //
2403 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
2404 Float_t MomRes3 = fMomRes3DTerm1->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2405 Double_t K3 = FSICorrelationOmega0(kTRUE, firstQ, secondQ, thirdQ);// K3
2406 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProdTermsCC3->Fill(Qsum, QsumIn, MomRes3/K3);
2407 }
2408 //////////////////////////////////////
2409 // Momentum resolution calculation
2410 if(fillIndex3==0 && fMCcase){
2411 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2412 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 0, 1, firstQMC, secondQMC, thirdQMC);
2413 if(ch1==ch2 && ch1==ch3){// same charge
2414 Float_t WInput = MCWeight3D(kTRUE, 1, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2415 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2416 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2417 }else {// mixed charge
2418 Float_t WInput=1.0;
2419 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 1, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2420 else WInput = MCWeight3D(kFALSE, 1, kMCfixedRbin, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2421 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2422 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2423 }
2424 }// fMCcase
2425
cd12341d 2426 }
2427
2428 }else if(config==2){// 12, 13, 23
5e3e77d6 2429
cd12341d 2430 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2431 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2432
2433 // loop over terms 2-4
5e3e77d6 2434 for(Int_t jj=2; jj<5; jj++){
2435 if(jj==2) {if(!fill2) continue;}//12
2436 if(jj==3) {if(!fill3) continue;}//13
2437 if(jj==4) {if(!fill4) continue;}//23
cd12341d 2438
2439 if(fillIndex3 <= 2){
5e3e77d6 2440 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
2441 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ);
2442 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
2443 Float_t InteractingQ=0;
2444 if(part==1) {InteractingQ=qinv12;}// 12 from SE
2445 else {InteractingQ=qinv13;}// 13 from SE
2446 Float_t MomRes3 = fMomRes3DTerm1->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2447 Double_t K3 = FSICorrelationOmega0(kTRUE, qinv12, qinv13, qinv23);// K3
2448 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProdTermsCC3->Fill(Qsum, QsumIn, MomRes3/K3);
2449 //
2450 if(jj==2) MomRes3 = fMomRes3DTerm2->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2451 else if(jj==3) MomRes3 = fMomRes3DTerm3->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2452 else MomRes3 = fMomRes3DTerm4->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2453 Float_t K2 = FSICorrelation2(+1,+1, kRVALUES-1, InteractingQ);// K2 from Therminator source
2454 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProdTermsCC2->Fill(Qsum, QsumIn, MomRes3/K2);
2455 }
2456 //////////////////////////////////////
2457 // Momentum resolution calculation
2458 if(fillIndex3==0 && fMCcase){
2459 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2460 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
2461 if(ch1==ch2 && ch1==ch3){// same charge
2462 Float_t WInput = MCWeight3D(kTRUE, jj, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2463 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2464 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2465 }else {// mixed charge
2466 Float_t WInput=1.0;
2467 if(bin1==bin2) WInput = MCWeight3D(kFALSE, jj, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2468 else WInput = MCWeight3D(kFALSE, 6-jj, kMCfixedRbin, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2469 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2470 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2471 }
2472 }// fMCcase
2473
cd12341d 2474 }
2475 }
5e3e77d6 2476
cd12341d 2477 }else {// config 3: All particles from different events
2478 //Simulate Filling of other event-mixed lowQ pairs
2479
2480 Float_t enhancement=1.0;
2481 Int_t nUnderCut=0;
2482 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2483 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2484
2485 if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2486 if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2487 if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
2488
2489 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2490
2491
2492 if(fillIndex3 <= 2){
2493 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
2494 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ, enhancement);
5e3e77d6 2495 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
2496 Float_t MomRes3 = fMomRes3DTerm1->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2497 Double_t K3 = FSICorrelationOmega0(kTRUE, firstQ, secondQ, thirdQ);// K3
2498 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProdTermsCC3->Fill(Qsum, QsumIn, MomRes3/K3);
2499 //
2500 MomRes3 = fMomRes3DTerm2->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);// Term2 used but equivalent to 3 and 4
2501 Float_t K2 = FSICorrelation2(+1,+1, kRVALUES-1, firstQ);// firstQ used but any Q can be used here since all are on equal footing.
2502 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProdTermsCC2->Fill(Qsum, QsumIn, MomRes3/K2);
2503 //
2504 MomRes3 = fMomRes3DTerm5->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2505 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProdTermsCC0->Fill(Qsum, QsumIn, MomRes3);// untouched version
2506 }
2507 //////////////////////////////////////
2508 // Momentum resolution calculation
2509 if(fillIndex3==0 && fMCcase){
2510 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2511 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, 5, firstQMC, secondQMC, thirdQMC);
2512 if(ch1==ch2 && ch1==ch3){// same charge
2513 Float_t WInput = MCWeight3D(kTRUE, 5, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2514 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2515 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2516 }else {// mixed charge
2517 Float_t WInput=1.0;
2518 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 5, kMCfixedRbin, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2519 else WInput = MCWeight3D(kFALSE, 5, kMCfixedRbin, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss in this case. 1st Q argument is ss
2520 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2521 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2522 }
2523 }// fMCcase
2524
cd12341d 2525 }
5e3e77d6 2526
cd12341d 2527 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2528 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
5e3e77d6 2529
cd12341d 2530 if(fMCcase) continue;// only calcualte TPN for real data
2531
2532 GetWeight(pVect1, pVect2, weight12, weight12Err);
2533 GetWeight(pVect1, pVect3, weight13, weight13Err);
2534 GetWeight(pVect2, pVect3, weight23, weight23Err);
2535 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2536
5e3e77d6 2537
2538 // Coul correlations from Wave-functions
2539 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm
2540 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){// 0.3 to 0.6
cd12341d 2541 Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2542 Int_t denIndex = rIter*kNDampValues + myDampIt;
5e3e77d6 2543 Int_t momResIndex = denIndex;
2544
2545 Float_t coulCorr12 = FSICorrelation2(+1,+1, rIter, qinv12);
2546 Float_t coulCorr13 = FSICorrelation2(+1,+1, rIter, qinv13);
2547 Float_t coulCorr23 = FSICorrelation2(+1,+1, rIter, qinv23);
2548 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2549 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2550 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
cd12341d 2551
5e3e77d6 2552 if(momBin12 >= kQbins) momBin12 = kQbins-1;
2553 if(momBin13 >= kQbins) momBin13 = kQbins-1;
2554 if(momBin23 >= kQbins) momBin23 = kQbins-1;
cd12341d 2555
5e3e77d6 2556 weight12CC = ((weight12+1)*fMomResC2->GetBinContent(momResIndex+1, momBin12) - myDamp*coulCorr12 - (1-myDamp));
2557 weight12CC /= coulCorr12*myDamp;
2558 weight13CC = ((weight13+1)*fMomResC2->GetBinContent(momResIndex+1, momBin13) - myDamp*coulCorr13 - (1-myDamp));
2559 weight13CC /= coulCorr13*myDamp;
2560 weight23CC = ((weight23+1)*fMomResC2->GetBinContent(momResIndex+1, momBin23) - myDamp*coulCorr23 - (1-myDamp));
2561 weight23CC /= coulCorr23*myDamp;
cd12341d 2562
5e3e77d6 2563 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {// testing purposes only
cd12341d 2564 if(denIndex==71 && fMbin==0 && ch1==-1) {
2565 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2566 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2567 }
2568 continue;// C2^QS can never be less than unity
2569 }
2570
2571 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
2572 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
5e3e77d6 2573 weightTotal *= fMomRes3DTerm5->GetBinContent(qinv12Bin, qinv13Bin, qinv23Bin);
2574 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProdTwoPartNorm->Fill(Qsum, QsumIn, enhancement*weightTotal);
cd12341d 2575
2576
2577 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2578 /*
2579 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2580 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2581 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2582 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2583 weightTotalErr /= pow(2*weightTotal,2);
2584
2585 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, enhancement*weightTotalErr);
2586 }
2587 */
2588
2589 }// damp iter
2590 }// R iter
2591
2592
2593
2594 }
2595 }// end 3rd particle
2596 }// en2
2597
2598
2599 }
2600 }
2601
2602
2603 ///////////////////
2604 }// end of PdensityPairs
2605
2606
2607
2608
2609
2610
2611
2612 ////////////////////////////////////////////////////////
2613 // Pdensity Method with Explicit Loops
2614 if(fPdensityExplicitLoop){
2615
2616 ////////////////////////////////////
2617 // 2nd, 3rd, and 4th order Correlations
2618
2619 // First Particle
2620 for (Int_t i=0; i<myTracks; i++) {
2621 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2622 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2623 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2624 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2625 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2626 key1 = (fEvt)->fTracks[i].fKey;
2627
2628 // Second Event
2629 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2630 Int_t startbin2=0;
2631 if(en2==0) startbin2=i+1;
2632
2633 // Second Particle
2634 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2635 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2636 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2637 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2638 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2639 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2640 key2 = (fEvt+en2)->fTracks[j].fKey;
2641
2642 Short_t fillIndex12 = FillIndex2part(key1+key2);
2643 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2644
2645 if(qinv12 < fQLowerCut) continue;
2646
2647
2648 // 2-particle part is filled always during pair creator
2649
2650 // Third Event
2651 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2652 Int_t startbin3=0;
2653 if(en3==en2) startbin3=j+1;
2654 else startbin3=0;
2655
2656
2657 // Third Particle
2658 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2659 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2660 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2661 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2662 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2663 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2664 key3 = (fEvt+en3)->fTracks[k].fKey;
2665
2666 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2667 Short_t fillIndex13 = FillIndex2part(key1+key3);
2668 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2669 Short_t fillIndex23 = FillIndex2part(key2+key3);
2670 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2671
2672
2673 if(qinv13 < fQLowerCut) continue;
2674 if(qinv23 < fQLowerCut) continue;
2675
2676
2677 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2678
2679 Short_t normBin12 = SetNormBin(fillIndex12);
2680 Short_t normBin13 = SetNormBin(fillIndex13);
2681 Short_t normBin23 = SetNormBin(fillIndex23);
2682
2683
2684 if(en3==0 && en2==0) {// 123
2685 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2686
2687 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2688
2689 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2690 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2691 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
2692 }
2693 }
2694
2695 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2696 Float_t gFact=1;
2697
2698 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2699 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2700
2701
2702 if(fill2){
2703 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2704 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2705 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2706 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2707 }
2708 }
2709 }
2710 if(fill3){
2711 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2712 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2713 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2714 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2715 }
2716 }
2717 }
2718 if(fill4){
2719 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2720 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2721 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2722 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2723 }
2724 }
2725 }
2726
2727 }else if(en2==1 && en3==2){// all uncorrelated events
2728 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2729
2730 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
2731 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2732 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2733 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
2734 }
2735 }
2736 Short_t qCutBin12 = SetQcutBin(fillIndex12);
2737 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2738 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2739
2740 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
2741
2742 Int_t nUnderCut=0;
2743 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
2744 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2745 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2746
2747 }
2748
2749 }else {}
2750
2751
2752 }// 3rd particle
2753 }// 3rd event
2754
2755 }// 2nd particle
2756 }// 2nd event
2757
2758 }// 1st particle
2759
2760
2761
2762
2763 }// End of PdensityExplicit
2764
2765
2766
2767
2768 // Post output data.
2769 PostData(1, fOutputList);
2770
2771}
2772//________________________________________________________________________
2773void AliChaoticity::Terminate(Option_t *)
2774{
2775 // Called once at the end of the query
2776
2777 cout<<"Done"<<endl;
2778
2779}
2780//________________________________________________________________________
2781Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
2782{
2783
2784 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
2785
2786 // propagate through B field to r=1m
5e3e77d6 2787 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
cd12341d 2788 if(phi1 > 2*PI) phi1 -= 2*PI;
2789 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 2790 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
cd12341d 2791 if(phi2 > 2*PI) phi2 -= 2*PI;
2792 if(phi2 < 0) phi2 += 2*PI;
2793
2794 Float_t deltaphi = phi1 - phi2;
2795 if(deltaphi > PI) deltaphi -= 2*PI;
2796 if(deltaphi < -PI) deltaphi += 2*PI;
2797 deltaphi = fabs(deltaphi);
2798
5e3e77d6 2799 //cout<<deltaphi<<" "<<fMinSepTPCEntrancePhi<<" "<<fMinSepTPCEntranceEta<<endl;
2800 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 2801
2802 // propagate through B field to r=1.6m
5e3e77d6 2803 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
cd12341d 2804 if(phi1 > 2*PI) phi1 -= 2*PI;
2805 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 2806 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
cd12341d 2807 if(phi2 > 2*PI) phi2 -= 2*PI;
2808 if(phi2 < 0) phi2 += 2*PI;
2809
2810 deltaphi = phi1 - phi2;
2811 if(deltaphi > PI) deltaphi -= 2*PI;
2812 if(deltaphi < -PI) deltaphi += 2*PI;
2813 deltaphi = fabs(deltaphi);
2814
5e3e77d6 2815 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 2816
2817
2818 //
2819
2820 Int_t ncl1 = first.fClusterMap.GetNbits();
2821 Int_t ncl2 = second.fClusterMap.GetNbits();
2822 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2823 Double_t shfrac = 0; Double_t qfactor = 0;
2824 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2825 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2826 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2827 sumQ++;
2828 sumCls+=2;
2829 sumSha+=2;}
2830 else {sumQ--; sumCls+=2;}
2831 }
2832 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2833 sumQ++;
2834 sumCls++;}
2835 }
2836 if (sumCls>0) {
2837 qfactor = sumQ*1.0/sumCls;
2838 shfrac = sumSha*1.0/sumCls;
2839 }
2840
2841 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2842
2843
2844 return kTRUE;
5e3e77d6 2845
cd12341d 2846
2847}
2848//________________________________________________________________________
2849Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2850{
2851 Float_t arg = G_Coeff/qinv;
2852
2853 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2854 else {return (exp(-arg)-1)/(-arg);}
2855
2856}
2857//________________________________________________________________________
2858void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2859{
2860 Int_t j, k;
2861 Int_t a = i2 - i1;
2862 for (Int_t i = i1; i < i2+1; i++) {
2863 j = (Int_t) (gRandom->Rndm() * a);
2864 k = iarr[j];
2865 iarr[j] = iarr[i];
2866 iarr[i] = k;
2867 }
2868}
2869//________________________________________________________________________
2870Short_t AliChaoticity::FillIndex2part(Short_t key){
2871
2872 if(key==2) return 0;// pi-pi
2873 else if(key==11) return 1;// pi-k
2874 else if(key==101) return 2;// pi-p
2875 else if(key==20) return 3;// k-k
2876 else if(key==110) return 4;// k-p
2877 else return 5;// p-p
2878}
2879//________________________________________________________________________
2880Short_t AliChaoticity::FillIndex3part(Short_t key){
2881
2882 if(key==3) return 0;// pi-pi-pi
2883 else if(key==12) return 1;// pi-pi-k
2884 else if(key==21) return 2;// k-k-pi
2885 else if(key==102) return 3;// pi-pi-p
2886 else if(key==201) return 4;// p-p-pi
2887 else if(key==111) return 5;// pi-k-p
2888 else if(key==30) return 6;// k-k-k
2889 else if(key==120) return 7;// k-k-p
2890 else if(key==210) return 8;// p-p-k
2891 else return 9;// p-p-p
2892
2893}
2894//________________________________________________________________________
2895Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
2896 if(fi <= 2) return 0;
2897 else if(fi==3) return 1;
2898 else return 2;
2899}
2900//________________________________________________________________________
2901Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
2902 if(fi==0) return 0;
2903 else if(fi <= 2) return 1;
2904 else return 2;
2905}
2906//________________________________________________________________________
2907void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
2908
2909 if(fi==0 || fi==3 || fi==5){// Identical species
2910 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
2911 else {b1=c1; b2=c2;}
2912 }else {// Mixed species
2913 if(key1 < key2) { b1=c1; b2=c2;}
2914 else {b1=c2; b2=c1;}
2915 }
2916
2917}
2918//________________________________________________________________________
2919void AliChaoticity::SetFillBins3(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Short_t part, Int_t &b1, Int_t &b2, Int_t &b3, Bool_t &fill2, Bool_t &fill3, Bool_t &fill4){
2920
2921
2922 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
2923 Bool_t seSS=kFALSE;
2924 Bool_t seSK=kFALSE;
2925 Short_t seKeySum=0;// only used for pi-k-p case
2926 if(part==1) {// default case (irrelevant for term 1 and term 5)
2927 if(c1==c2) seSS=kTRUE;
2928 if(key1==key2) seSK=kTRUE;
2929 seKeySum = key1+key2;
2930 }
2931 if(part==2){
2932 if(c1==c3) seSS=kTRUE;
2933 if(key1==key3) seSK=kTRUE;
2934 seKeySum = key1+key3;
2935 }
2936
2937
2938 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
2939
2940 if(fi==0 || fi==6 || fi==9){// Identical species
2941 if( (c1+c2+c3)==1) {
2942 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
2943 //
2944 if(seSS) fill2=kTRUE;
2945 else {fill3=kTRUE; fill4=kTRUE;}
2946 //
2947 }else if( (c1+c2+c3)==2) {
2948 b1=0; b2=1; b3=1;
2949 //
2950 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
2951 else fill4=kTRUE;
2952 //
2953 }else {
2954 b1=c1; b2=c2; b3=c3;
2955 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
2956 }
2957 }else if(fi != 5){// all the rest except pi-k-p
2958 if(key1==key2){
2959 b3=c3;
2960 if( (c1+c2)==1) {b1=0; b2=1;}
2961 else {b1=c1; b2=c2;}
2962 }else if(key1==key3){
2963 b3=c2;
2964 if( (c1+c3)==1) {b1=0; b2=1;}
2965 else {b1=c1; b2=c3;}
2966 }else {// Key2==Key3
2967 b3=c1;
2968 if( (c2+c3)==1) {b1=0; b2=1;}
2969 else {b1=c2; b2=c3;}
2970 }
2971 //////////////////////////////
2972 if(seSK) fill2=kTRUE;// Same keys from Same Event
2973 else {// Different keys from Same Event
2974 if( (c1+c2+c3)==1) {
2975 if(b3==0) {
2976 if(seSS) fill3=kTRUE;
2977 else fill4=kTRUE;
2978 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
2979 }else if( (c1+c2+c3)==2) {
2980 if(b3==1) {
2981 if(seSS) fill4=kTRUE;
2982 else fill3=kTRUE;
2983 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
2984 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
2985 }
2986 /////////////////////////////
2987 }else {// pi-k-p (no charge ordering applies since all are unique)
2988 if(key1==1){
2989 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
2990 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
2991 }else if(key1==10){
2992 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
2993 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
2994 }else {// key1==100
2995 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
2996 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
2997 }
2998 ////////////////////////////////////
2999 if(seKeySum==11) fill2=kTRUE;
3000 else if(seKeySum==101) fill3=kTRUE;
3001 else fill4=kTRUE;
3002 ////////////////////////////////////
3003 }
3004
3005}
3006//________________________________________________________________________
3007void AliChaoticity::ArrangeQs(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Float_t q12, Float_t q13, Float_t q23, Short_t part, Short_t term, Float_t &fQ, Float_t &sQ, Float_t &tQ){
3008
3009 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
3010 if(fi==0 || fi==6 || fi==9){// Identical species
3011 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
3012 if(term==1 || term==5){
3013 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
3014 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
3015 else {fQ=q23; sQ=q12; tQ=q13;}
3016 }else if(term==2 && part==1){
3017 fQ=q12; sQ=q13; tQ=q23;
3018 }else if(term==2 && part==2){
3019 fQ=q13; sQ=q12; tQ=q23;
3020 }else if(term==3 && part==1){
3021 sQ=q12;
3022 if(c1==c3) {fQ=q13; tQ=q23;}
3023 else {fQ=q23; tQ=q13;}
3024 }else if(term==3 && part==2){
3025 sQ=q13;
3026 if(c1==c2) {fQ=q12; tQ=q23;}
3027 else {fQ=q23; tQ=q12;}
3028 }else if(term==4 && part==1){
3029 tQ=q12;
3030 if(c1==c3) {fQ=q13; sQ=q23;}
3031 else {fQ=q23; sQ=q13;}
3032 }else if(term==4 && part==2){
3033 tQ=q13;
3034 if(c1==c2) {fQ=q12; sQ=q23;}
3035 else {fQ=q23; sQ=q12;}
3036 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3037 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
3038 if(term==1 || term==5){
3039 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
3040 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
3041 else {tQ=q23; sQ=q12; fQ=q13;}
3042 }else if(term==2 && part==1){
3043 fQ=q12;
3044 if(c1==c3) {tQ=q13; sQ=q23;}
3045 else {tQ=q23; sQ=q13;}
3046 }else if(term==2 && part==2){
3047 fQ=q13;
3048 if(c1==c2) {tQ=q12; sQ=q23;}
3049 else {tQ=q23; sQ=q12;}
3050 }else if(term==3 && part==1){
3051 sQ=q12;
3052 if(c1==c3) {tQ=q13; fQ=q23;}
3053 else {tQ=q23; fQ=q13;}
3054 }else if(term==3 && part==2){
3055 sQ=q13;
3056 if(c1==c2) {tQ=q12; fQ=q23;}
3057 else {tQ=q23; fQ=q12;}
3058 }else if(term==4 && part==1){
3059 tQ=q12; sQ=q13; fQ=q23;
3060 }else if(term==4 && part==2){
3061 tQ=q13; sQ=q12; fQ=q23;
3062 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3063 }else {// fQ=ss, sQ=ss, tQ=ss
3064 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
3065 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
3066 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
3067 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
3068 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
3069 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
3070 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
3071 }
3072 }else if(fi != 5){// all the rest except pi-k-p
3073 if(key1==key2){
3074 fQ=q12;
3075 if(c1==c2){
3076 // cases not explicity shown below are not possible
3077 if(term==1 || term==5) {sQ=q13; tQ=q23;}
3078 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
3079 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
3080 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
3081 else cout<<"problem!!!!!!!!!!!!!"<<endl;
3082 }else if(c3==0){
3083 if(c1==c3) {sQ=q13; tQ=q23;}
3084 else {sQ=q23; tQ=q13;}
3085 }else {//c3==1
3086 if(c1==c3) {tQ=q13; sQ=q23;}
3087 else {tQ=q23; sQ=q13;}
3088 }
3089 }else if(key1==key3){
3090 fQ=q13;
3091 if(c1==c3){
3092 // cases not explicity shown below are not possible
3093 if(term==1 || term==5) {sQ=q12; tQ=q23;}
3094 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
3095 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
3096 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
3097 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3098 }else if(c2==0){
3099 if(c1==c2) {sQ=q12; tQ=q23;}
3100 else {sQ=q23; tQ=q12;}
3101 }else {//c2==1
3102 if(c1==c2) {tQ=q12; sQ=q23;}
3103 else {tQ=q23; sQ=q12;}
3104 }
3105 }else {// key2==key3
3106 fQ=q23;
3107 if(c2==c3){
3108 // cases not explicity shown below are not possible
3109 if(term==1 || term==5) {sQ=q12; tQ=q13;}
3110 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
3111 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
3112 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
3113 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
3114 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3115 }else if(c1==0){
3116 if(c1==c2) {sQ=q12; tQ=q13;}
3117 else {sQ=q13; tQ=q12;}
3118 }else {//c1==1
3119 if(c1==c2) {tQ=q12; sQ=q13;}
3120 else {tQ=q13; sQ=q12;}
3121 }
3122 }
3123 }else {// pi-k-p
3124 if(key1==1){
3125 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
3126 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
3127 }else if(key1==10){
3128 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
3129 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
3130 }else {// key1==100
3131 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
3132 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
3133 }
3134
3135 }
3136
3137
3138}
3139//________________________________________________________________________
3140Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
3141
3142 Float_t qinv=1.0;
3143
3144 if(fi==0 || fi==3 || fi==5){// identical masses
3145 qinv = sqrt( pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2));
3146 }else{// different masses
3147 Float_t px = track1[1] + track2[1];
3148 Float_t py = track1[2] + track2[2];
3149 Float_t pz = track1[3] + track2[3];
3150 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
3151 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
3152 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
3153
3154 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
3155 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
3156 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
3157 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
3158 qinv = sqrt(qinv);
3159 }
3160
3161 return qinv;
3162
3163}
3164//________________________________________________________________________
3165void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3166
3167 Float_t p0 = track1[0] + track2[0];
3168 Float_t px = track1[1] + track2[1];
3169 Float_t py = track1[2] + track2[2];
3170 Float_t pz = track1[3] + track2[3];
3171
3172 Float_t mt = sqrt(p0*p0 - pz*pz);
3173 Float_t pt = sqrt(px*px + py*py);
3174
3175 Float_t v0 = track1[0] - track2[0];
3176 Float_t vx = track1[1] - track2[1];
3177 Float_t vy = track1[2] - track2[2];
3178 Float_t vz = track1[3] - track2[3];
3179
3180 qout = (px*vx + py*vy)/pt;
3181 qside = (px*vy - py*vx)/pt;
3182 qlong = (p0*vz - pz*v0)/mt;
3183}
3184//________________________________________________________________________
5e3e77d6 3185void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[kKbinsT][kCentBins]){
cd12341d 3186
5e3e77d6 3187 if(legoCase){
3188 for(Int_t mb=0; mb<fMbins; mb++){
3189 for(Int_t edB=0; edB<kEDbins; edB++){
3190 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3191 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3192 //
3193 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3194 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3195 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3196
3197 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3198 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
3199 }
3200 }
3201 }
3202 }
3203 }
3204 //
3205 }
3206 }
3207 return;
3208 }// legoCase
3209
cd12341d 3210 for(Int_t mb=0; mb<fMbins; mb++){
3211 for(Int_t edB=0; edB<kEDbins; edB++){
3212 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3213 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3214 //
3215 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3216 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3217 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3218
3219 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3220 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3221 }
3222 }
3223 }
3224 }
3225 }
3226 //
3227 }
3228 }
3229
3230
5e3e77d6 3231
cd12341d 3232 TFile *wFile = new TFile("WeightFile.root","READ");
3233 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3234 else cout<<"Good Weight File Found!"<<endl;
5e3e77d6 3235
cd12341d 3236 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3237 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3238 for(Int_t mb=0; mb<fMbins; mb++){
3239 for(Int_t edB=0; edB<kEDbins; edB++){
3240
3241 TString *name = new TString("Weight_Kt_");
3242 *name += tKbin;
3243 name->Append("_Ky_");
3244 *name += yKbin;
3245 name->Append("_M_");
3246 *name += mb;
3247 name->Append("_ED_");
3248 *name += edB;
3249
3250 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3251
3252 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3253 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3254 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3255
3256 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3257 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3258 }
3259 }
3260 }
3261 delete fTempHisto;
3262 }//ED
3263 }//mb
3264 }//ky
3265 }//kt
3266
3267 wFile->Close();
5e3e77d6 3268
cd12341d 3269
5e3e77d6 3270 cout<<"Done reading weight file"<<endl;
cd12341d 3271
3272}
3273//________________________________________________________________________
cd12341d 3274void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3275
3276 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3277 Float_t ky=0;// central rapidity
3278 //
3279 Float_t qOut=0,qSide=0,qLong=0;
3280 GetQosl(track1, track2, qOut, qSide, qLong);
3281 qOut = fabs(qOut);
3282 qSide = fabs(qSide);
3283 qLong = fabs(qLong);
3284 //
3285
3286 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3287 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3288 else {
3289 for(Int_t i=0; i<kKbinsT-1; i++){
3290 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3291 }
3292 }
3293 //
3294 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3295 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3296 else {
3297 for(Int_t i=0; i<kKbinsY-1; i++){
3298 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3299 }
3300 }
3301 //
3302 /////////
3303 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3304 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3305 else {
3306 for(Int_t i=0; i<kQbinsWeights-1; i++){
3307 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3308 }
3309 }
3310 //
3311 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3312 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3313 else {
3314 for(Int_t i=0; i<kQbinsWeights-1; i++){
3315 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3316 }
3317 }
3318 //
3319 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3320 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3321 else {
3322 for(Int_t i=0; i<kQbinsWeights-1; i++){
3323 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3324 }
3325 }
3326 //
3327
3328
3329 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3330 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3331
3332
3333 Float_t deltaW=0;
3334 // kt
3335 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3336 // Ky
3337 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3338 // Qo
5e3e77d6 3339 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3340 // Qs
5e3e77d6 3341 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3342 // Ql
5e3e77d6 3343 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3344
3345 //
3346 wgt = min + deltaW;
3347
3348
3349 ////
3350
3351 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3352 Float_t deltaWErr=0;
3353 // Kt
3354 /*
3355 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3356 // Ky
3357 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3358 // Qo
5e3e77d6 3359 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3360 // Qs
5e3e77d6 3361 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3362 // Ql
5e3e77d6 3363 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3364 */
3365 wgtErr = minErr + deltaWErr;
3366
3367
3368
cd12341d 3369}
3370//________________________________________________________________________
3371Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3372
5e3e77d6 3373 Float_t radius = Float_t(rIndex+3.)/0.19733;// convert to GeV
3374 if(rIndex==kRVALUES-1) radius = 6.0/0.19733;// Therminator default radius
cd12341d 3375 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
5e3e77d6 3376 Float_t coulCorr12 = FSICorrelation2(charge1, charge2, rIndex, qinv);
cd12341d 3377 if(charge1==charge2){
5e3e77d6 3378 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius,2)))*coulCorr12);
cd12341d 3379 }else {
5e3e77d6 3380 return ((1-myDamp) + myDamp*coulCorr12);
cd12341d 3381 }
3382
3383}
3384//________________________________________________________________________
5e3e77d6 3385Float_t AliChaoticity::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t rIndex, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
cd12341d 3386 if(term==5) return 1.0;
5e3e77d6 3387
3388 Float_t radius = (3. + rIndex)/0.19733;//starts at 5fm. convert to GeV
cd12341d 3389 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3390 Float_t fc = sqrt(myDamp);
5e3e77d6 3391 if(SameCharge){// all three of the same charge
3392 Float_t coulCorr12 = FSICorrelation2(+1,+1, rIndex, q12);// K2
3393 Float_t coulCorr13 = FSICorrelation2(+1,+1, rIndex, q13);// K2
3394 Float_t coulCorr23 = FSICorrelation2(+1,+1, rIndex, q23);// K2
cd12341d 3395
5e3e77d6 3396 if(term==1){
3397 Float_t c3QS = 1 + exp(-pow(q12*radius,2)) + exp(-pow(q13*radius,2)) + exp(-pow(q23*radius,2));
3398 c3QS += 2*exp(-pow(radius,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3399 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3400 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3401 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius,2)))*coulCorr13;
3402 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius,2)))*coulCorr23;
3403 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kTRUE, q12, q13, q23);
3404 return w123;
3405 }else if(term==2){
3406 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3407 }else if(term==3){
3408 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius,2)))*coulCorr13);
3409 }else if(term==4){
3410 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius,2)))*coulCorr23);
3411 }else return 1.0;
3412
3413 }else{// mixed charge case pair 12 always treated as ss
3414 Float_t coulCorr12 = FSICorrelation2(+1,+1, rIndex, q12);// K2 ss
3415 Float_t coulCorr13 = FSICorrelation2(+1,-1, rIndex, q13);// K2 os
3416 Float_t coulCorr23 = FSICorrelation2(+1,-1, rIndex, q23);// K2 os
3417 if(term==1){
3418 Float_t c3QS = 1 + exp(-pow(q12*radius,2));
3419 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3420 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3421 w123 += pow(fc,2)*(1-fc)*coulCorr13;
3422 w123 += pow(fc,2)*(1-fc)*coulCorr23;
3423 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kFALSE, q12, q13, q23);
3424 return w123;
3425 }else if(term==2){
3426 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3427 }else if(term==3){
3428 return ((1-myDamp) + myDamp*coulCorr13);
3429 }else if(term==4){
3430 return ((1-myDamp) + myDamp*coulCorr23);
3431 }else return 1.0;
3432 }
3433
cd12341d 3434}
3435//________________________________________________________________________
3436Float_t AliChaoticity::MCWeightOSL(Float_t radius, Float_t myDamp, Float_t Q_out, Float_t Q_side, Float_t Q_long, Float_t qinv){
3437 Int_t rIndex = Int_t(radius+0.001-3);
5e3e77d6 3438 Float_t coulCorr12 = FSICorrelation2(+1,+1, rIndex, qinv);
cd12341d 3439
5e3e77d6 3440 return ((1-myDamp) + myDamp*(1 + exp(-pow(radius/0.19733,2)*(Q_out*Q_out+Q_side*Q_side+Q_long*Q_long)))*coulCorr12);
cd12341d 3441
3442}
3443//________________________________________________________________________
5e3e77d6 3444void AliChaoticity::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D, TH3D *temp3D[5]){
3445
3446
3447 if(legoCase){
3448 fMomResC2 = (TH2D*)temp2D->Clone();
3449 fMomRes3DTerm1=(TH3D*)temp3D[0]->Clone();
3450 fMomRes3DTerm2=(TH3D*)temp3D[1]->Clone();
3451 fMomRes3DTerm3=(TH3D*)temp3D[2]->Clone();
3452 fMomRes3DTerm4=(TH3D*)temp3D[3]->Clone();
3453 fMomRes3DTerm5=(TH3D*)temp3D[4]->Clone();
3454 //
3455 fMomResC2->SetDirectory(0);
3456 fMomRes3DTerm1->SetDirectory(0);
3457 fMomRes3DTerm2->SetDirectory(0);
3458 fMomRes3DTerm3->SetDirectory(0);
3459 fMomRes3DTerm4->SetDirectory(0);
3460 fMomRes3DTerm5->SetDirectory(0);
3461
3462 }else {
3463 TFile *momResFile = new TFile("MomResFile.root","READ");
3464 if(!momResFile->IsOpen()) {
3465 cout<<"No momentum resolution file found"<<endl;
3466 AliFatal("No momentum resolution file found. Kill process.");
3467 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3468
3469 TH2D *temp2D2 = (TH2D*)momResFile->Get("MomResHisto_pp");
3470 TH3D *temp3Dterm1 = (TH3D*)momResFile->Get("MomResHisto_3d_ppp_term1");
3471 TH3D *temp3Dterm2 = (TH3D*)momResFile->Get("MomResHisto_3d_ppp_term2");
3472 TH3D *temp3Dterm3 = (TH3D*)momResFile->Get("MomResHisto_3d_ppp_term3");
3473 TH3D *temp3Dterm4 = (TH3D*)momResFile->Get("MomResHisto_3d_ppp_term4");
3474 TH3D *temp3Dterm5 = (TH3D*)momResFile->Get("MomResHisto_3d_ppp_term5");
3475
3476 fMomResC2 = (TH2D*)temp2D2->Clone();
3477 fMomRes3DTerm1=(TH3D*)temp3Dterm1->Clone();
3478 fMomRes3DTerm2=(TH3D*)temp3Dterm2->Clone();
3479 fMomRes3DTerm3=(TH3D*)temp3Dterm3->Clone();
3480 fMomRes3DTerm4=(TH3D*)temp3Dterm4->Clone();
3481 fMomRes3DTerm5=(TH3D*)temp3Dterm5->Clone();
3482 //
3483 fMomResC2->SetDirectory(0);
3484 fMomRes3DTerm1->SetDirectory(0);
3485 fMomRes3DTerm2->SetDirectory(0);
3486 fMomRes3DTerm3->SetDirectory(0);
3487 fMomRes3DTerm4->SetDirectory(0);
3488 fMomRes3DTerm5->SetDirectory(0);
3489
3490 momResFile->Close();
3491 }
cd12341d 3492
5e3e77d6 3493 cout<<"Done reading momentum resolution file"<<endl;
3494}
3495//________________________________________________________________________
3496void AliChaoticity::SetFSICorrelations(Bool_t legoCase, TH2D *temp2D[2], TH3D *temp3D[2]){
3497 // read in 2-particle and 3-particle FSI correlations = K2 & K3
3498 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
3499 if(legoCase){
3500 fFSI2SS = (TH2D*)temp2D[0]->Clone();
3501 fFSI2OS = (TH2D*)temp2D[1]->Clone();
3502 fFSIOmega0SS = (TH3D*)temp3D[0]->Clone();
3503 fFSIOmega0OS = (TH3D*)temp3D[1]->Clone();
3504 //
3505 fFSI2SS->SetDirectory(0);
3506 fFSI2OS->SetDirectory(0);
3507 fFSIOmega0SS->SetDirectory(0);
3508 fFSIOmega0OS->SetDirectory(0);
3509 }else {
3510 TFile *fsifile = new TFile("KFile.root","READ");
3511 if(!fsifile->IsOpen()) {
3512 cout<<"No FSI file found"<<endl;
3513 AliFatal("No FSI file found. Kill process.");
3514 }else {cout<<"Good FSI File Found!"<<endl;}
3515
3516 TH2D *temphisto2SS = (TH2D*)fsifile->Get("K2ss");
3517 TH2D *temphisto2OS = (TH2D*)fsifile->Get("K2os");
3518 TH3D *temphisto3SS = (TH3D*)fsifile->Get("K3ss");
3519 TH3D *temphisto3OS = (TH3D*)fsifile->Get("K3os");
3520
3521 fFSI2SS = (TH2D*)temphisto2SS->Clone();
3522 fFSI2OS = (TH2D*)temphisto2OS->Clone();
3523 fFSIOmega0SS = (TH3D*)temphisto3SS->Clone();
3524 fFSIOmega0OS = (TH3D*)temphisto3OS->Clone();
3525 //
3526 fFSI2SS->SetDirectory(0);
3527 fFSI2SS->SetDirectory(0);
3528 fFSIOmega0SS->SetDirectory(0);
3529 fFSIOmega0OS->SetDirectory(0);
3530
3531 fsifile->Close();
cd12341d 3532 }
3533
5e3e77d6 3534 // condition FSI histogram for edge effects
3535 for(Int_t ii=1; ii<=fFSIOmega0SS->GetNbinsX(); ii++){
3536 for(Int_t jj=1; jj<=fFSIOmega0SS->GetNbinsY(); jj++){
3537 for(Int_t kk=1; kk<=fFSIOmega0SS->GetNbinsZ(); kk++){
3538
3539 if(fFSIOmega0SS->GetBinContent(ii,jj,kk) <=0){
3540 Double_t Q12 = fFSIOmega0SS->GetXaxis()->GetBinCenter(ii);
3541 Double_t Q23 = fFSIOmega0SS->GetYaxis()->GetBinCenter(jj);
3542 Double_t Q13 = fFSIOmega0SS->GetZaxis()->GetBinCenter(kk);
3543 //
3544 Int_t Q12bin=ii;
3545 Int_t Q23bin=jj;
3546 Int_t Q13bin=kk;
3547 Int_t AC=0;//Adjust Counter
3548 Int_t AClimit=10;// maximum bin shift
3549 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin++; AC++;}}
3550 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin--; AC++;}}
3551 //
3552 if(Q13 < sqrt(pow(Q12,2)+pow(Q23,2) - 2*Q12*Q23)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin++; AC++;}}
3553 if(Q13 > sqrt(pow(Q12,2)+pow(Q23,2) + 2*Q12*Q23)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin--; AC++;}}
3554 //
3555 if(Q23 < sqrt(pow(Q12,2)+pow(Q13,2) - 2*Q12*Q13)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin++; AC++;}}
3556 if(Q23 > sqrt(pow(Q12,2)+pow(Q13,2) + 2*Q12*Q13)) {while(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin--; AC++;}}
3557
3558 // Save cpu time by setting empty cell contents (edge effects) to nearest non-zero cell (these cells are not used very often anyway.)
3559 if(AC==AClimit) {
3560 fFSIOmega0SS->SetBinContent(ii,jj,kk, 1.0);
3561 fFSIOmega0OS->SetBinContent(ii,jj,kk, 1.0);
3562 }else {
3563 fFSIOmega0SS->SetBinContent(ii,jj,kk, fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin));
3564 fFSIOmega0OS->SetBinContent(ii,jj,kk, fFSIOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin));
3565 }
3566 }
3567
3568 }
cd12341d 3569 }
3570 }
3571
5e3e77d6 3572 cout<<"Done reading FSI file"<<endl;
cd12341d 3573}
3574//________________________________________________________________________
5e3e77d6 3575Float_t AliChaoticity::FSICorrelation2(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
3576 // returns 2-particle Coulomb correlations = K2
3577 if(rIndex >= kRVALUES) return 1.0;
3578 Int_t qbinL = fFSI2SS->GetYaxis()->FindBin(qinv-fFSI2SS->GetYaxis()->GetBinWidth(1)/2.);
3579 Int_t qbinH = qbinL+1;
3580 if(qbinL <= 0) return 1.0;
3581 if(qbinH > fFSI2SS->GetNbinsY()) return 1.0;
3582
3583 Float_t slope=0;
3584 if(charge1==charge2){
3585 slope = fFSI2SS->GetBinContent(rIndex+1, qbinL) - fFSI2SS->GetBinContent(rIndex+1, qbinH);
3586 slope /= fFSI2SS->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS->GetYaxis()->GetBinCenter(qbinH);
3587 return (slope*(qinv - fFSI2SS->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS->GetBinContent(rIndex+1, qbinL));
3588 }else {
3589 slope = fFSI2OS->GetBinContent(rIndex+1, qbinL) - fFSI2OS->GetBinContent(rIndex+1, qbinH);
3590 slope /= fFSI2OS->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS->GetYaxis()->GetBinCenter(qbinH);
3591 return (slope*(qinv - fFSI2OS->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS->GetBinContent(rIndex+1, qbinL));
cd12341d 3592 }
3593}
3594//________________________________________________________________________
5e3e77d6 3595Double_t AliChaoticity::FSICorrelationOmega0(Bool_t SameCharge, Double_t Q12, Double_t Q13, Double_t Q23){
3596 // remember: Omega0 histogram is in the following order (Q12, Q23, Q13)
3597 // returns 3d 3-particle Coulomb Correlation = K3
3598 Int_t Q12bin = fFSIOmega0SS->GetXaxis()->FindBin(Q12);
3599 Int_t Q13bin = fFSIOmega0SS->GetZaxis()->FindBin(Q13);
3600 Int_t Q23bin = fFSIOmega0SS->GetYaxis()->FindBin(Q23);
3601
3602 if(SameCharge){
3603 if(fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3604 else return fFSIOmega0SS->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
3605 }else{// mixed charge
3606 if(fFSIOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3607 else return fFSIOmega0OS->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
3608 }
cd12341d 3609}
3610//________________________________________________________________________