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