]> 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
code cleanup, documentation, placement of 'using' statements
[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),
47 fESD(0x0),
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),
186 fESD(0x0),
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),
337 fESD(obj.fESD),
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;
435 fESD = obj.fESD;
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;
515 if(fESD) delete fESD;
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
1092 cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1093 fEventCounter++;
1094
1095 if(fAODcase) fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1096 else fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1097
1098 if(fAODcase) {if (!fAOD) {Printf("ERROR: fAOD not available"); return;}}
1099 else {if (!fESD) {Printf("ERROR: fESD not available"); return;}}
1100
1101 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1102
1103 // Trigger Cut
1104 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1105 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1106 if(!isSelected1 && !fMCcase) {cout<<"Event Rejected"<<endl; return;}
1107 }if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1108 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1109 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1110 if(!isSelected1 && !isSelected2 && !fMCcase) {cout<<"Event Rejected"<<endl; return;}
1111 }else {cout<<"Event Rejected"<<endl; return;}
1112
1113 ///////////////////////////////////////////////////////////
1114 const AliAODVertex *primaryVertexAOD;
1115 AliCentrality *centrality;// for AODs and ESDs
1116
1117
1118 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1119 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1120 fPIDResponse = inputHandler->GetPIDResponse();
1121
1122
1123 TClonesArray *mcArray = 0x0;
1124 if(fMCcase){
1125 if(fAODcase){
1126 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1127 if(!mcArray || mcArray->GetEntriesFast() >= 110000){
1128 cout<<"No MC particle branch found or Array too large!!"<<endl;
1129 cout<<mcArray->GetEntriesFast()<<endl;
1130 return;
1131 }
1132 }
1133 }
1134
1135
1136 UInt_t status=0;
1137 Int_t positiveTracks=0, negativeTracks=0;
1138 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1139
1140 Double_t vertex[3]={0};
1141 Int_t zbin=0;
1142 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1143 /////////////////////////////////////////////////
1144
1145
1146 Float_t centralityPercentile=0;
1147 Float_t cStep=5.0, cStart=0;
1148
1149
1150 if(fAODcase){// AOD case
1151
1152 if(fPbPbcase){
1153 centrality = fAOD->GetCentrality();
1154 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1155 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1156 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {cout<<"Centrality out of Range. Skipping Event"<<endl; return;}
1157 cout<<"Centrality % = "<<centralityPercentile<<endl;
1158 }
1159
1160
1161
1162
1163 ////////////////////////////////
1164 // Vertexing
1165 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1166 primaryVertexAOD = fAOD->GetPrimaryVertex();
1167 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1168
1169 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1170 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1171
1172 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1173 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1174
1175 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1176
1177 fBfield = fAOD->GetMagneticField();
1178
1179 for(Int_t i=0; i<fZvertexBins; i++){
1180 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1181 zbin=i;
1182 break;
1183 }
1184 }
1185
1186
1187
1188 /////////////////////////////
1189 // Create Shuffled index list
1190 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1191 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1192 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1193 /////////////////////////////
1194
1195 // Track loop
1196 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1197 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1198 if (!aodtrack) continue;
1199 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1200
1201 status=aodtrack->GetStatus();
1202
1203 if((Bool_t)(((1<<7) & aodtrack->GetFilterMap()) == 0)) continue;// AOD filterBit cut
1204 // 1<<5 is for "standard cuts with tight dca cut"
1205 // 1<<7 is for TPC only tracking
1206 if(aodtrack->Pt() < 0.16) continue;
1207 if(fabs(aodtrack->Eta()) > 0.8) continue;
1208
1209
1210 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1211 if(!goodMomentum) continue;
1212 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1213
1214 Float_t dca2[2];
1215 Float_t dca3d;
1216
1217 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1218 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1219 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1220
1221 fTempStruct[myTracks].fStatus = status;
1222 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1223 fTempStruct[myTracks].fId = aodtrack->GetID();
1224 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1225 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1226 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1227 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1228 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1229 fTempStruct[myTracks].fEta = aodtrack->Eta();
1230 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1231 fTempStruct[myTracks].fDCAXY = dca2[0];
1232 fTempStruct[myTracks].fDCAZ = dca2[1];
1233 fTempStruct[myTracks].fDCA = dca3d;
1234 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1235 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1236
1237
1238
1239 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1240 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1241 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1242
1243 if(fTempStruct[myTracks].fCharge==+1) {
1244 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1245 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1246 }else {
1247 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1248 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1249 }
1250
1251 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1252 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1253
1254
1255 // nSimga PID workaround
1256 fTempStruct[myTracks].fElectron = kFALSE;
1257 fTempStruct[myTracks].fPion = kFALSE;
1258 fTempStruct[myTracks].fKaon = kFALSE;
1259 fTempStruct[myTracks].fProton = kFALSE;
1260
1261 Float_t nSigmaTPC[5];
1262 Float_t nSigmaTOF[5];
1263 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1264 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1265 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1266 Float_t signalTPC=0, signalTOF=0;
1267 Double_t integratedTimesTOF[10]={0};
1268 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1269 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1270 if (!aodTrack2) continue;
1271 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1272
1273 UInt_t status2=aodTrack2->GetStatus();
1274
1275 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1276 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1277 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1278 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1279 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1280 //
1281 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1282 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1283 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1284 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1285 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1286 signalTPC = aodTrack2->GetTPCsignal();
1287
1288 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1289 fTempStruct[myTracks].fTOFhit = kTRUE;
1290 signalTOF = aodTrack2->GetTOFsignal();
1291 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1292 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1293
1294 }// aodTrack2
1295
1296 ///////////////////
1297 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1298 if(fTempStruct[myTracks].fTOFhit) {
1299 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1300 }
1301 ///////////////////
1302
1303 // Use TOF if good hit and above threshold
1304 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1305 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1306 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1307 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1308 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1309 }else {// TPC info instead
1310 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1311 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1312 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1313 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1314 }
1315
1316
1317 //////////////////////////////////////
1318 // Bayesian PIDs for TPC only tracking
1319 //const Double_t* PID = aodtrack->PID();
1320 //fTempStruct[myTracks].fElectron = kFALSE;
1321 //fTempStruct[myTracks].Pion = kFALSE;
1322 //fTempStruct[myTracks].Kaon = kFALSE;
1323 //fTempStruct[myTracks].Proton = kFALSE;
1324
1325 // Pions
1326 //if(PID[2] > 0.2) fTempStruct[myTracks].Pion = kTRUE;// pions: 0.2 --> 0.5
1327 //
1328 // Kaons
1329 //if(PID[3] <= 0.5) fTempStruct[myTracks].Kaon = kFALSE;// kaons
1330 //
1331 // Protons
1332 //if(PID[4] <= 0.5) fTempStruct[myTracks].Proton = kFALSE;// protons
1333 //////////////////////////////////////
1334
1335
1336 // Ensure there is only 1 candidate per track
1337 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1338 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1339 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1340 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1341 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1342 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1343 ////////////////////////
1344 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1345
1346 if(!fTempStruct[myTracks].fPion) continue;// only pions
1347
1348
1349
1350
1351 if(fTempStruct[myTracks].fPion) {// pions
1352 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1353 fTempStruct[myTracks].fKey = 1;
1354 }else if(fTempStruct[myTracks].fKaon){// kaons
1355 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1356 fTempStruct[myTracks].fKey = 10;
1357 }else{// protons
1358 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1359 fTempStruct[myTracks].fKey = 100;
1360 }
1361
1362
1363
1364 if(aodtrack->Charge() > 0) positiveTracks++;
1365 else negativeTracks++;
1366
1367 if(fTempStruct[myTracks].fPion) pionCount++;
1368 if(fTempStruct[myTracks].fKaon) kaonCount++;
1369 if(fTempStruct[myTracks].fProton) protonCount++;
1370
1371 myTracks++;
1372
1373 }
1374 }else {// ESD tracks
1375 cout<<"ESDs not supported currently"<<endl;
1376 return;
1377 }
1378
1379
1380 if(myTracks >= 1) {
1381 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1382 }
1383
1384
1385 cout<<"There are "<<myTracks<<" myTracks"<<endl;
1386 cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
1387
1388 /////////////////////////////////////////
1389 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1390 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1391 /////////////////////////////////////////
1392
1393
1394 ////////////////////////////////
1395 ///////////////////////////////
1396 // Mbin determination
1397 //
1398 // Mbin set to Pion Count Only for pp!!!!!!!
1399 fMbin=-1;
1400 if(!fPbPbcase){
1401 for(Int_t i=0; i<kMultBinspp; i++){
1402 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1403 // Mbin 0 has 1 pion
1404 }
1405 }else{
1406 for(Int_t i=0; i<kCentBins; i++){
1407 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1408 fMbin=i;// 0 = most central
1409 break;
1410 }
1411 }
1412 }
1413
1414 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1415
1416 //////////////////////////////////////////////////
1417 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1418 //////////////////////////////////////////////////
1419
1420
1421
1422 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1423 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1424
1425 ////////////////////////////////////
1426 // Add event to buffer if > 0 tracks
1427 if(myTracks > 0){
1428 fEC[zbin][fMbin]->FIFOShift();
1429 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1430 (fEvt)->fNtracks = myTracks;
1431 (fEvt)->fFillStatus = 1;
1432 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1433 if(fMCcase){
1434 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1435 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1436 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1437 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1438 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1439 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1440 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1441 }
1442 }
1443 }
1444
1445
1446
1447 Float_t qinv12=0, qinv13=0, qinv23=0;
1448 Float_t qout=0, qside=0, qlong=0;
1449 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1450 Float_t transK12=0, rapK12=0, transK3=0;
1451 Int_t transKbin=0, rapKbin=0;
1452 Float_t q3=0;
1453 Int_t ch1=0, ch2=0, ch3=0;
1454 Short_t key1=0, key2=0, key3=0;
1455 Int_t bin1=0, bin2=0, bin3=0;
1456 Float_t pVect1[4]={0};
1457 Float_t pVect2[4]={0};
1458 Float_t pVect3[4]={0};
1459 Float_t pVect1MC[4]={0};
1460 Float_t pVect2MC[4]={0};
1461 Int_t index1=0, index2=0, index3=0;
1462 Float_t weight12=0, weight13=0, weight23=0;
1463 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1464 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1465 Float_t weightTotal=0;//, weightTotalErr=0;
1466 //
1467 Float_t qinv12MC=0;
1468 AliAODMCParticle *mcParticle1=0x0;
1469 AliAODMCParticle *mcParticle2=0x0;
1470
1471
1472 if(fPdensityPairCut){
1473 ////////////////////
1474 Int_t pairCountSE=0, pairCountME=0;
1475 Int_t normPairCount[2]={0};
1476 Int_t numOtherPairs1[2][fMultLimit];
1477 Int_t numOtherPairs2[2][fMultLimit];
1478 Bool_t exitCode=kFALSE;
1479 Int_t tempNormFillCount[2][2][2][10][5];
1480
1481
1482 // reset to defaults
1483 for(Int_t i=0; i<fMultLimit; i++) {
1484 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1485 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1486
1487 // Normalization Utilities
1488 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1489 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1490 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1491 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1492 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1493 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1494 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1495 numOtherPairs1[0][i]=0;
1496 numOtherPairs1[1][i]=0;
1497 numOtherPairs2[0][i]=0;
1498 numOtherPairs2[1][i]=0;
1499
1500 // Track Merging/Splitting Utilities
1501 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1502 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1503 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1504 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1505 }
1506
1507 // Reset the temp Normalization counters
1508 for(Int_t i=0; i<2; i++){// Charge1
1509 for(Int_t j=0; j<2; j++){// Charge2
1510 for(Int_t k=0; k<2; k++){// Charge3
1511 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1512 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1513 tempNormFillCount[i][j][k][l][m] = 0;
1514 }
1515 }
1516 }
1517 }
1518 }
1519
1520
1521 ///////////////////////////////////////////////////////
1522 // Start the pairing process
1523 // P11 pairing
1524 // 1st Particle
1525
1526 for (Int_t i=0; i<myTracks; i++) {
1527
1528 Int_t en2=0;
1529
1530 // 2nd particle
1531 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1532
1533 key1 = (fEvt)->fTracks[i].fKey;
1534 key2 = (fEvt+en2)->fTracks[j].fKey;
1535 Short_t fillIndex2 = FillIndex2part(key1+key2);
1536 Short_t qCutBin = SetQcutBin(fillIndex2);
1537 Short_t normBin = SetNormBin(fillIndex2);
1538 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1539 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1540 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1541 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1542 //
1543
1544 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1545 GetQosl(pVect1, pVect2, qout, qside, qlong);
1546 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1547
1548 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1549
1550 ///////////////////////////////
1551 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1552 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1553 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1554
1555 if(fMCcase && ch1==ch2 && fMbin==0){
1556 for(Int_t rstep=0; rstep<10; rstep++){
1557 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1558 // propagate through B field to r=1.2m
1559 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1560 if(phi1 > 2*PI) phi1 -= 2*PI;
1561 if(phi1 < 0) phi1 += 2*PI;
1562 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1563 if(phi2 > 2*PI) phi2 -= 2*PI;
1564 if(phi2 < 0) phi2 += 2*PI;
1565 Float_t deltaphi = phi1 - phi2;
1566 if(deltaphi > PI) deltaphi -= PI;
1567 if(deltaphi < -PI) deltaphi += PI;
1568 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1569 }
1570 }
1571 // Pair Splitting/Merging cut
1572 if(ch1 == ch2){
1573 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1574 fPairSplitCut[0][i]->AddAt('1',j);
1575 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1576 continue;
1577 }
1578 }
1579
1580 // HIJING tests
1581 if(fMCcase && fillIndex2==0){
1582
1583 // Check that label does not exceed stack size
1584 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1585 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1586 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1587 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1588 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1589 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1590 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1591 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1592 if(qinv12<0.1 && ch1==ch2) {
1593 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1594 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1595 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
1596 }
1597
1598
1599 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 10fm
1600 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
1601 Int_t denIndex = rIter*kNDampValues + myDampIt;
1602
1603 Float_t WInput = MCWeight(ch1,ch2, rIter, myDampIt, qinv12MC);
1604 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
1605 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
1606 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
1607 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
1608 }
1609 }
1610
1611 //HIJING resonance test
1612 if(ch1 != ch2){
1613 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
1614 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
1615 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
1616 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
1617 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
1618 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
1619 }
1620 }
1621 }
1622 }// label check 2
1623 }// MC case
1624
1625 //////////////////////////////////////////
1626 // 2-particle term
1627 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
1628
1629 // osl frame
1630 if(bin1==bin2 && fillIndex2==0){
1631 if((transK12 > 0.2) && (transK12 < 0.3)){
1632 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1633 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1634 }
1635 if((transK12 > 0.6) && (transK12 < 0.7)){
1636 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1637 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1638 }
1639 }
1640 //////////////////////////////////////////
1641 if(fTabulatePairs){
1642 if(fillIndex2==0 && bin1==bin2){
1643 rapK12 = 0;
1644 transKbin=-1; rapKbin=-1;
1645 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1646 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1647 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1648 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1649 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1650 continue;
1651 }
1652 }
1653
1654
1655 // exit out of loop if there are too many pairs
1656 if(pairCountSE >= kPairLimit) {cout<<"Too many SE pairs"<<endl; exitCode=kTRUE; continue;}
1657 if(exitCode) continue;
1658
1659 //////////////////////////
1660 // Enforce the Qcut
1661 if(qinv12 <= fQcut[qCutBin]) {
1662 ///////////////////////////
1663 // particle 1
1664 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
1665 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
1666 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
1667 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
1668 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
1669 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
1670 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
1671 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
1672 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1673 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1674 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1675 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
1676 }
1677 // particle 2
1678 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1679 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1680 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1681 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1682 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1683 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
1684 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
1685 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1686 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1687 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1688 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1689 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1690 }
1691
1692 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
1693
1694 fPairLocationSE[i]->AddAt(pairCountSE,j);
1695
1696 pairCountSE++;
1697
1698 }
1699
1700
1701 /////////////////////////////////////////////////////////
1702 // Normalization Region
1703
1704 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1705 // particle 1
1706 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1707 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1708 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1709 // particle 2
1710 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1711 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1712 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1713
1714
1715 //other past pairs with particle j
1716 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
1717 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
1718 if(locationOtherPair < 0) continue;// no pair there
1719 Int_t indexOther1 = i;
1720 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
1721
1722 // Both possible orderings of other indexes
1723 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
1724
1725 // 1 and 2 are from SE
1726 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
1727 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
1728 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1729 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
1730
1731 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
1732 }
1733
1734 }// pastpair P11 loop
1735
1736
1737 fNormPairSwitch[en2][i]->AddAt('1',j);
1738 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1739 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1740
1741 numOtherPairs1[en2][i]++;
1742 numOtherPairs2[en2][j]++;
1743
1744
1745 normPairCount[en2]++;
1746 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1747
1748 }// Norm Region
1749
1750 }// j particle
1751 }// i particle
1752
1753
1754
1755 //////////////////////////////////////////////
1756 // P12 pairing
1757 // 1st Particle
1758 for (Int_t i=0; i<myTracks; i++) {
1759
1760 Int_t en2=1;
1761
1762 // 2nd particle
1763 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1764
1765 key1 = (fEvt)->fTracks[i].fKey;
1766 key2 = (fEvt+en2)->fTracks[j].fKey;
1767 Short_t fillIndex2 = FillIndex2part(key1+key2);
1768 Short_t qCutBin = SetQcutBin(fillIndex2);
1769 Short_t normBin = SetNormBin(fillIndex2);
1770 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1771 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1772 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1773 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1774
1775 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1776 GetQosl(pVect1, pVect2, qout, qside, qlong);
1777 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1778
1779 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
1780
1781 ///////////////////////////////
1782 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1783 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1784 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1785
1786 if(fMCcase && ch1==ch2 && fMbin==0){
1787 for(Int_t rstep=0; rstep<10; rstep++){
1788 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1789 // propagate through B field to r=1.2m
1790 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1791 if(phi1 > 2*PI) phi1 -= 2*PI;
1792 if(phi1 < 0) phi1 += 2*PI;
1793 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1794 if(phi2 > 2*PI) phi2 -= 2*PI;
1795 if(phi2 < 0) phi2 += 2*PI;
1796 Float_t deltaphi = phi1 - phi2;
1797 if(deltaphi > PI) deltaphi -= PI;
1798 if(deltaphi < -PI) deltaphi += PI;
1799 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1800 }
1801 }
1802
1803 if(ch1 == ch2){
1804 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1805 fPairSplitCut[1][i]->AddAt('1',j);
1806 continue;
1807 }
1808 }
1809
1810 //////////////////////////////////////////
1811 // 2-particle term
1812 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
1813 // osl frame
1814 if(bin1==bin2 && fillIndex2==0){
1815 if((transK12 > 0.2) && (transK12 < 0.3)){
1816 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1817 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1818 }
1819 if((transK12 > 0.6) && (transK12 < 0.7)){
1820 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1821 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1822 }
1823 }
1824 //////////////////////////////////////////
1825 if(fTabulatePairs){
1826 if(fillIndex2==0 && bin1==bin2){
1827 rapK12 = 0;
1828 transKbin=-1; rapKbin=-1;
1829 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1830 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1831 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1832 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1833 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1834 continue;
1835 }
1836 }
1837
1838
1839 if(pairCountME >= 2*kPairLimit) {cout<<"Too many ME pairs"<<endl; exitCode=kTRUE; continue;}
1840 if(exitCode) continue;
1841
1842 if(qinv12 <= fQcut[qCutBin]) {
1843 ///////////////////////////
1844
1845 // particle 1
1846 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
1847 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
1848 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
1849 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
1850 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
1851 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
1852 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
1853 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
1854 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1855 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1856 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1857 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
1858 }
1859 // particle 2
1860 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1861 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1862 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1863 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1864 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1865 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
1866 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
1867 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1868 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1869 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1870 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1871 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1872 }
1873
1874 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
1875
1876 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
1877
1878 pairCountME++;
1879
1880 }
1881
1882 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1883 // particle 1
1884 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1885 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1886 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1887 // particle 2
1888 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1889 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1890 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1891
1892 //other past pairs in P11 with particle i
1893 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
1894 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
1895 if(locationOtherPairP11 < 0) continue;// no pair there
1896 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
1897
1898 //Check other past pairs in P12
1899 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
1900
1901 // 1 and 3 are from SE
1902 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
1903 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
1904 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1905 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
1906 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
1907
1908
1909 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
1910 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
1911 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
1912
1913
1914 }// P11 loop
1915
1916
1917 fNormPairSwitch[en2][i]->AddAt('1',j);
1918 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1919 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1920
1921 numOtherPairs1[en2][i]++;
1922 numOtherPairs2[en2][j]++;
1923
1924 normPairCount[en2]++;
1925 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1926
1927 }// Norm Region
1928
1929
1930 }
1931 }
1932
1933
1934 ///////////////////////////////////////
1935 // P13 pairing (just for Norm counting of term5)
1936 for (Int_t i=0; i<myTracks; i++) {
1937
1938 // exit out of loop if there are too many pairs
1939 // dont bother with this loop if exitCode is set.
1940 if(exitCode) break;
1941
1942 // 2nd particle
1943 Int_t en2=2;
1944
1945 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1946
1947 key1 = (fEvt)->fTracks[i].fKey;
1948 key2 = (fEvt+en2)->fTracks[j].fKey;
1949 Short_t fillIndex2 = FillIndex2part(key1+key2);
1950 Short_t normBin = SetNormBin(fillIndex2);
1951 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1952 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1953 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1954 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1955
1956 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1957
1958 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
1959
1960 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1961 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1962
1963 if(ch1 == ch2){
1964 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1965 fPairSplitCut[2][i]->AddAt('1',j);
1966 continue;
1967 }
1968 }
1969
1970 /////////////////////////////////////////////////////////
1971 // Normalization Region
1972
1973 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1974
1975 fNormPairSwitch[en2][i]->AddAt('1',j);
1976
1977 }// Norm Region
1978 }
1979 }
1980
1981
1982
1983 ///////////////////////////////////////
1984 // P23 pairing (just for Norm counting of term5)
1985 Int_t en1=1;
1986 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
1987
1988 // exit out of loop if there are too many pairs
1989 // dont bother with this loop if exitCode is set.
1990 if(exitCode) break;
1991
1992 // 2nd event
1993 Int_t en2=2;
1994 // 2nd particle
1995 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1996
1997 if(exitCode) break;
1998
1999 key1 = (fEvt+en1)->fTracks[i].fKey;
2000 key2 = (fEvt+en2)->fTracks[j].fKey;
2001 Short_t fillIndex2 = FillIndex2part(key1+key2);
2002 Short_t normBin = SetNormBin(fillIndex2);
2003 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2004 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2005 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2006 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2007
2008 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2009
2010 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2011
2012 ///////////////////////////////
2013 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2014 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2015
2016 if(ch1 == ch2){
2017 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2018 fPairSplitCut[3][i]->AddAt('1',j);
2019 continue;
2020 }
2021 }
2022
2023 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2024
2025 Int_t index1P23 = i;
2026 Int_t index2P23 = j;
2027
2028 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2029 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2030 if(locationOtherPairP12 < 0) continue; // no pair there
2031 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2032
2033
2034 //Check other past pair status in P13
2035 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2036
2037 // all from different event
2038 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2039 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2040 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2041 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2042
2043 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2044 }
2045 }
2046 }
2047
2048
2049
2050
2051 ///////////////////////////////////////////////////
2052 // Do not use pairs from events with too many pairs
2053 if(exitCode) {
2054 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2055 (fEvt)->fNpairsSE = 0;
2056 (fEvt)->fNpairsME = 0;
2057 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2058 return;// Skip event
2059 }else{
2060 (fEvt)->fNpairsSE = pairCountSE;
2061 (fEvt)->fNpairsME = pairCountME;
2062 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2063 }
2064 ///////////////////////////////////////////////////
2065
2066
2067
2068 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
2069 cout<<"Start Main analysis"<<endl;
2070
2071 ///////////////////////////////////////////////////////////////////////
2072 ///////////////////////////////////////////////////////////////////////
2073 ///////////////////////////////////////////////////////////////////////
2074 //
2075 //
2076 // Start the Correlation Analysis
2077 //
2078 //
2079 ///////////////////////////////////////////////////////////////////////
2080
2081
2082 /////////////////////////////////////////////////////////
2083 // Skip 3-particle part if Tabulate6DPairs is set to true
2084 if(fTabulatePairs) return;
2085 /////////////////////////////////////////////////////////
2086
2087 // Set the Normalization counters
2088 for(Int_t termN=0; termN<5; termN++){
2089
2090 if(termN==0){
2091 if((fEvt)->fNtracks ==0) continue;
2092 }else if(termN<4){
2093 if((fEvt)->fNtracks ==0) continue;
2094 if((fEvt+1)->fNtracks ==0) continue;
2095 }else {
2096 if((fEvt)->fNtracks ==0) continue;
2097 if((fEvt+1)->fNtracks ==0) continue;
2098 if((fEvt+2)->fNtracks ==0) continue;
2099 }
2100
2101 for(Int_t sc=0; sc<kSCLimit3; sc++){
2102
2103 for(Int_t c1=0; c1<2; c1++){
2104 for(Int_t c2=0; c2<2; c2++){
2105 for(Int_t c3=0; c3<2; c3++){
2106
2107 if(sc==0 || sc==6 || sc==9){// Identical species
2108 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2109 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2110 }else if(sc!=5){
2111 if( (c1+c2)==1) {if(c1!=0) continue;}
2112 }else {}// do nothing for pi-k-p case
2113 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2114 }
2115 }
2116 }
2117 }
2118 }
2119
2120
2121
2122 /////////////////////////////////////////////
2123 // Calculate Pair-Cut Correlations
2124 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2125
2126 Int_t nump1=0;
2127 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2128 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2129
2130 // 1st pair
2131 for(Int_t p1=0; p1<nump1; p1++){
2132
2133 if(en1case==0){
2134 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2135 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2136 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2137 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2138 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2139 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2140 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2141 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2142 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2143 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2144 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2145 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2146 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2147
2148 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2149 }
2150 if(en1case==1){
2151 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2152 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2153 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2154 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2155 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2156 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2157 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2158 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2159 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2160 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2161 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2162 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2163 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2164
2165 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2166 }
2167
2168
2169
2170 // en2 buffer
2171 for(Int_t en2=0; en2<3; en2++){
2172 //////////////////////////////////////
2173
2174 Bool_t skipcase=kTRUE;
2175 Short_t config=-1, part=-1;
2176 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2177 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2178 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2179 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2180
2181 if(skipcase) continue;
2182
2183
2184 // 3-particle terms
2185 // 3rd particle
2186 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2187 index3 = k;
2188
2189
2190 // remove auto-correlations and duplicate triplets
2191 if(config==1){
2192 if( index1 == index3) continue;
2193 if( index2 == index3) continue;
2194 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2195
2196 // skip the switched off triplets
2197 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2198 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2199 continue;
2200 }
2201 ///////////////////////////////
2202 // Turn off 1st possible degenerate triplet
2203 if(index1 < index3){// verify correct id ordering ( index1 < k )
2204 if(fPairLocationSE[index1]->At(index3) >= 0){
2205 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2206 }
2207 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2208 }else {// or k < index1
2209 if(fPairLocationSE[index3]->At(index1) >= 0){
2210 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2211 }
2212 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2213 }
2214 // turn off 2nd possible degenerate triplet
2215 if(index2 < index3){// verify correct id ordering (index2 < k)
2216 if(fPairLocationSE[index2]->At(index3) >= 0){
2217 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2218 }
2219 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2220 }else {// or k < index2
2221 if(fPairLocationSE[index3]->At(index2) >= 0){
2222 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2223 }
2224 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2225 }
2226
2227 }// end config 1
2228
2229 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2230 ///////////////////////////////
2231 // Turn off 1st possible degenerate triplet
2232 if(fPairLocationME[index1]->At(index3) >= 0){
2233 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2234 }
2235
2236 // turn off 2nd possible degenerate triplet
2237 if(fPairLocationME[index2]->At(index3) >= 0){
2238 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2239 }
2240
2241 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2242 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2243 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2244 }// end config 2 part 1
2245
2246 if(config==2 && part==2){// P12T1
2247 if( index1 == index3) continue;
2248
2249 // skip the switched off triplets
2250 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2251 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2252 continue;
2253 }
2254 // turn off another possible degenerate
2255 if(fPairLocationME[index3]->At(index2) >= 0){
2256 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2257 }// end config 2 part 2
2258
2259 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2260 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2261 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2262 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2263 }
2264 if(config==3){// P12T3
2265 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2266 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2267 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2268 }// end config 3
2269
2270
2271
2272 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2273 key3 = (fEvt+en2)->fTracks[k].fKey;
2274 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2275 Short_t fillIndex13 = FillIndex2part(key1+key3);
2276 Short_t fillIndex23 = FillIndex2part(key2+key3);
2277 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2278 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2279 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2280 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2281 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2282 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
2283 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2284 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2285
2286 if(qinv13 < fQLowerCut) continue;
2287 if(qinv23 < fQLowerCut) continue;
2288 if(qinv13 > 3.*fQcut[qCutBin13]) continue;
2289 if(qinv23 > 3.*fQcut[qCutBin23]) continue;
2290
2291 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2292 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2293 Float_t firstQ=0, secondQ=0, thirdQ=0;
2294
2295
2296 if(config==1) {// 123
2297 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2298
2299 if(fillIndex3 <= 2){
2300 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
2301 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ);
2302 ((TH2F*)fOutputList->FindObject("fKt3Dist"))->Fill(fMbin+1, transK3);
2303 }
2304
2305 }else if(config==2){// 12, 13, 23
2306
2307 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2308 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2309
2310 // loop over terms 2-4
2311 for(Int_t jj=0; jj<3; jj++){
2312 if(jj==0) {if(!fill2) continue;}//12
2313 if(jj==1) {if(!fill3) continue;}//13
2314 if(jj==2) {if(!fill4) continue;}//23
2315
2316 if(fillIndex3 <= 2){
2317 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj+2, firstQ, secondQ, thirdQ);
2318 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj+1].fTerms3->Fill(firstQ, secondQ, thirdQ);
2319 }
2320 }
2321
2322 }else {// config 3: All particles from different events
2323 //Simulate Filling of other event-mixed lowQ pairs
2324
2325 Float_t enhancement=1.0;
2326 Int_t nUnderCut=0;
2327 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2328 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2329
2330 if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2331 if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2332 if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
2333
2334 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2335
2336
2337 if(fillIndex3 <= 2){
2338 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
2339 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ, enhancement);
2340 }
2341
2342
2343 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2344 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
2345 if(fMCcase) continue;// only calcualte TPN for real data
2346
2347 GetWeight(pVect1, pVect2, weight12, weight12Err);
2348 GetWeight(pVect1, pVect3, weight13, weight13Err);
2349 GetWeight(pVect2, pVect3, weight23, weight23Err);
2350 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2351
2352
2353
2354 // Coul corrections from Wave-functions
2355 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 5fm to 20fm
2356 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
2357 Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2358 Int_t denIndex = rIter*kNDampValues + myDampIt;
2359
2360
2361 Float_t coulCorr12 = CoulCorr(+1,+1, rIter, qinv12);
2362 Float_t coulCorr13 = CoulCorr(+1,+1, rIter, qinv13);
2363 Float_t coulCorr23 = CoulCorr(+1,+1, rIter, qinv23);
2364 weight12CC = ((weight12+1)*GetMomRes(denIndex, qinv12) - myDamp/coulCorr12 - (1-myDamp));
2365 weight12CC *= coulCorr12/myDamp;
2366 weight13CC = ((weight13+1)*GetMomRes(denIndex, qinv13) - myDamp/coulCorr13 - (1-myDamp));
2367 weight13CC *= coulCorr13/myDamp;
2368 weight23CC = ((weight23+1)*GetMomRes(denIndex, qinv23) - myDamp/coulCorr23 - (1-myDamp));
2369 weight23CC *= coulCorr23/myDamp;
2370
2371 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {
2372 if(denIndex==71 && fMbin==0 && ch1==-1) {
2373 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2374 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2375 }
2376 continue;// C2^QS can never be less than unity
2377 }
2378
2379 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
2380 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2381
2382
2383 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2384 /*
2385 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2386 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2387 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2388 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2389 weightTotalErr /= pow(2*weightTotal,2);
2390
2391 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, enhancement*weightTotalErr);
2392 }
2393 */
2394
2395 }// damp iter
2396 }// R iter
2397
2398
2399
2400 }
2401 }// end 3rd particle
2402 }// en2
2403
2404
2405 }
2406 }
2407
2408
2409 ///////////////////
2410 }// end of PdensityPairs
2411
2412
2413
2414
2415
2416
2417
2418 ////////////////////////////////////////////////////////
2419 // Pdensity Method with Explicit Loops
2420 if(fPdensityExplicitLoop){
2421
2422 ////////////////////////////////////
2423 // 2nd, 3rd, and 4th order Correlations
2424
2425 // First Particle
2426 for (Int_t i=0; i<myTracks; i++) {
2427 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2428 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2429 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2430 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2431 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2432 key1 = (fEvt)->fTracks[i].fKey;
2433
2434 // Second Event
2435 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2436 Int_t startbin2=0;
2437 if(en2==0) startbin2=i+1;
2438
2439 // Second Particle
2440 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2441 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2442 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2443 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2444 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2445 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2446 key2 = (fEvt+en2)->fTracks[j].fKey;
2447
2448 Short_t fillIndex12 = FillIndex2part(key1+key2);
2449 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2450
2451 if(qinv12 < fQLowerCut) continue;
2452
2453
2454 // 2-particle part is filled always during pair creator
2455
2456 // Third Event
2457 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2458 Int_t startbin3=0;
2459 if(en3==en2) startbin3=j+1;
2460 else startbin3=0;
2461
2462
2463 // Third Particle
2464 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2465 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2466 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2467 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2468 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2469 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2470 key3 = (fEvt+en3)->fTracks[k].fKey;
2471
2472 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2473 Short_t fillIndex13 = FillIndex2part(key1+key3);
2474 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2475 Short_t fillIndex23 = FillIndex2part(key2+key3);
2476 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2477
2478
2479 if(qinv13 < fQLowerCut) continue;
2480 if(qinv23 < fQLowerCut) continue;
2481
2482
2483 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2484
2485 Short_t normBin12 = SetNormBin(fillIndex12);
2486 Short_t normBin13 = SetNormBin(fillIndex13);
2487 Short_t normBin23 = SetNormBin(fillIndex23);
2488
2489
2490 if(en3==0 && en2==0) {// 123
2491 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2492
2493 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2494
2495 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2496 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2497 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
2498 }
2499 }
2500
2501 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2502 Float_t gFact=1;
2503
2504 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2505 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2506
2507
2508 if(fill2){
2509 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2510 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2511 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2512 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2513 }
2514 }
2515 }
2516 if(fill3){
2517 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2518 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2519 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2520 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2521 }
2522 }
2523 }
2524 if(fill4){
2525 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2526 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2527 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2528 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2529 }
2530 }
2531 }
2532
2533 }else if(en2==1 && en3==2){// all uncorrelated events
2534 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2535
2536 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
2537 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2538 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2539 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
2540 }
2541 }
2542 Short_t qCutBin12 = SetQcutBin(fillIndex12);
2543 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2544 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2545
2546 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
2547
2548 Int_t nUnderCut=0;
2549 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
2550 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2551 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2552
2553 }
2554
2555 }else {}
2556
2557
2558 }// 3rd particle
2559 }// 3rd event
2560
2561 }// 2nd particle
2562 }// 2nd event
2563
2564 }// 1st particle
2565
2566
2567
2568
2569 }// End of PdensityExplicit
2570
2571
2572
2573
2574 // Post output data.
2575 PostData(1, fOutputList);
2576
2577}
2578//________________________________________________________________________
2579void AliChaoticity::Terminate(Option_t *)
2580{
2581 // Called once at the end of the query
2582
2583 cout<<"Done"<<endl;
2584
2585}
2586//________________________________________________________________________
2587Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
2588{
2589
2590 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
2591
2592 // propagate through B field to r=1m
2593 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// mine. 0.1798 for D=1.2m
2594 if(phi1 > 2*PI) phi1 -= 2*PI;
2595 if(phi1 < 0) phi1 += 2*PI;
2596 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// mine. 0.1798 for D=1.2m
2597 if(phi2 > 2*PI) phi2 -= 2*PI;
2598 if(phi2 < 0) phi2 += 2*PI;
2599
2600 Float_t deltaphi = phi1 - phi2;
2601 if(deltaphi > PI) deltaphi -= 2*PI;
2602 if(deltaphi < -PI) deltaphi += 2*PI;
2603 deltaphi = fabs(deltaphi);
2604
2605 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Sep just after TPC entrance
2606
2607 // propagate through B field to r=1.6m
2608 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.1798 for D=1.2m
2609 if(phi1 > 2*PI) phi1 -= 2*PI;
2610 if(phi1 < 0) phi1 += 2*PI;
2611 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.1798 for D=1.2m
2612 if(phi2 > 2*PI) phi2 -= 2*PI;
2613 if(phi2 < 0) phi2 += 2*PI;
2614
2615 deltaphi = phi1 - phi2;
2616 if(deltaphi > PI) deltaphi -= 2*PI;
2617 if(deltaphi < -PI) deltaphi += 2*PI;
2618 deltaphi = fabs(deltaphi);
2619
2620 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Sep just after TPC entrance
2621
2622
2623 //
2624
2625 Int_t ncl1 = first.fClusterMap.GetNbits();
2626 Int_t ncl2 = second.fClusterMap.GetNbits();
2627 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2628 Double_t shfrac = 0; Double_t qfactor = 0;
2629 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2630 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2631 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2632 sumQ++;
2633 sumCls+=2;
2634 sumSha+=2;}
2635 else {sumQ--; sumCls+=2;}
2636 }
2637 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2638 sumQ++;
2639 sumCls++;}
2640 }
2641 if (sumCls>0) {
2642 qfactor = sumQ*1.0/sumCls;
2643 shfrac = sumSha*1.0/sumCls;
2644 }
2645
2646 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2647
2648
2649 return kTRUE;
2650
2651
2652}
2653//________________________________________________________________________
2654Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2655{
2656 Float_t arg = G_Coeff/qinv;
2657
2658 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2659 else {return (exp(-arg)-1)/(-arg);}
2660
2661}
2662//________________________________________________________________________
2663void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2664{
2665 Int_t j, k;
2666 Int_t a = i2 - i1;
2667 for (Int_t i = i1; i < i2+1; i++) {
2668 j = (Int_t) (gRandom->Rndm() * a);
2669 k = iarr[j];
2670 iarr[j] = iarr[i];
2671 iarr[i] = k;
2672 }
2673}
2674//________________________________________________________________________
2675Short_t AliChaoticity::FillIndex2part(Short_t key){
2676
2677 if(key==2) return 0;// pi-pi
2678 else if(key==11) return 1;// pi-k
2679 else if(key==101) return 2;// pi-p
2680 else if(key==20) return 3;// k-k
2681 else if(key==110) return 4;// k-p
2682 else return 5;// p-p
2683}
2684//________________________________________________________________________
2685Short_t AliChaoticity::FillIndex3part(Short_t key){
2686
2687 if(key==3) return 0;// pi-pi-pi
2688 else if(key==12) return 1;// pi-pi-k
2689 else if(key==21) return 2;// k-k-pi
2690 else if(key==102) return 3;// pi-pi-p
2691 else if(key==201) return 4;// p-p-pi
2692 else if(key==111) return 5;// pi-k-p
2693 else if(key==30) return 6;// k-k-k
2694 else if(key==120) return 7;// k-k-p
2695 else if(key==210) return 8;// p-p-k
2696 else return 9;// p-p-p
2697
2698}
2699//________________________________________________________________________
2700Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
2701 if(fi <= 2) return 0;
2702 else if(fi==3) return 1;
2703 else return 2;
2704}
2705//________________________________________________________________________
2706Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
2707 if(fi==0) return 0;
2708 else if(fi <= 2) return 1;
2709 else return 2;
2710}
2711//________________________________________________________________________
2712void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
2713
2714 if(fi==0 || fi==3 || fi==5){// Identical species
2715 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
2716 else {b1=c1; b2=c2;}
2717 }else {// Mixed species
2718 if(key1 < key2) { b1=c1; b2=c2;}
2719 else {b1=c2; b2=c1;}
2720 }
2721
2722}
2723//________________________________________________________________________
2724void 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){
2725
2726
2727 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
2728 Bool_t seSS=kFALSE;
2729 Bool_t seSK=kFALSE;
2730 Short_t seKeySum=0;// only used for pi-k-p case
2731 if(part==1) {// default case (irrelevant for term 1 and term 5)
2732 if(c1==c2) seSS=kTRUE;
2733 if(key1==key2) seSK=kTRUE;
2734 seKeySum = key1+key2;
2735 }
2736 if(part==2){
2737 if(c1==c3) seSS=kTRUE;
2738 if(key1==key3) seSK=kTRUE;
2739 seKeySum = key1+key3;
2740 }
2741
2742
2743 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
2744
2745 if(fi==0 || fi==6 || fi==9){// Identical species
2746 if( (c1+c2+c3)==1) {
2747 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
2748 //
2749 if(seSS) fill2=kTRUE;
2750 else {fill3=kTRUE; fill4=kTRUE;}
2751 //
2752 }else if( (c1+c2+c3)==2) {
2753 b1=0; b2=1; b3=1;
2754 //
2755 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
2756 else fill4=kTRUE;
2757 //
2758 }else {
2759 b1=c1; b2=c2; b3=c3;
2760 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
2761 }
2762 }else if(fi != 5){// all the rest except pi-k-p
2763 if(key1==key2){
2764 b3=c3;
2765 if( (c1+c2)==1) {b1=0; b2=1;}
2766 else {b1=c1; b2=c2;}
2767 }else if(key1==key3){
2768 b3=c2;
2769 if( (c1+c3)==1) {b1=0; b2=1;}
2770 else {b1=c1; b2=c3;}
2771 }else {// Key2==Key3
2772 b3=c1;
2773 if( (c2+c3)==1) {b1=0; b2=1;}
2774 else {b1=c2; b2=c3;}
2775 }
2776 //////////////////////////////
2777 if(seSK) fill2=kTRUE;// Same keys from Same Event
2778 else {// Different keys from Same Event
2779 if( (c1+c2+c3)==1) {
2780 if(b3==0) {
2781 if(seSS) fill3=kTRUE;
2782 else fill4=kTRUE;
2783 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
2784 }else if( (c1+c2+c3)==2) {
2785 if(b3==1) {
2786 if(seSS) fill4=kTRUE;
2787 else fill3=kTRUE;
2788 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
2789 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
2790 }
2791 /////////////////////////////
2792 }else {// pi-k-p (no charge ordering applies since all are unique)
2793 if(key1==1){
2794 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
2795 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
2796 }else if(key1==10){
2797 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
2798 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
2799 }else {// key1==100
2800 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
2801 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
2802 }
2803 ////////////////////////////////////
2804 if(seKeySum==11) fill2=kTRUE;
2805 else if(seKeySum==101) fill3=kTRUE;
2806 else fill4=kTRUE;
2807 ////////////////////////////////////
2808 }
2809
2810}
2811//________________________________________________________________________
2812void 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){
2813
2814 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
2815 if(fi==0 || fi==6 || fi==9){// Identical species
2816 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
2817 if(term==1 || term==5){
2818 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
2819 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
2820 else {fQ=q23; sQ=q12; tQ=q13;}
2821 }else if(term==2 && part==1){
2822 fQ=q12; sQ=q13; tQ=q23;
2823 }else if(term==2 && part==2){
2824 fQ=q13; sQ=q12; tQ=q23;
2825 }else if(term==3 && part==1){
2826 sQ=q12;
2827 if(c1==c3) {fQ=q13; tQ=q23;}
2828 else {fQ=q23; tQ=q13;}
2829 }else if(term==3 && part==2){
2830 sQ=q13;
2831 if(c1==c2) {fQ=q12; tQ=q23;}
2832 else {fQ=q23; tQ=q12;}
2833 }else if(term==4 && part==1){
2834 tQ=q12;
2835 if(c1==c3) {fQ=q13; sQ=q23;}
2836 else {fQ=q23; sQ=q13;}
2837 }else if(term==4 && part==2){
2838 tQ=q13;
2839 if(c1==c2) {fQ=q12; sQ=q23;}
2840 else {fQ=q23; sQ=q12;}
2841 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
2842 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
2843 if(term==1 || term==5){
2844 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
2845 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
2846 else {tQ=q23; sQ=q12; fQ=q13;}
2847 }else if(term==2 && part==1){
2848 fQ=q12;
2849 if(c1==c3) {tQ=q13; sQ=q23;}
2850 else {tQ=q23; sQ=q13;}
2851 }else if(term==2 && part==2){
2852 fQ=q13;
2853 if(c1==c2) {tQ=q12; sQ=q23;}
2854 else {tQ=q23; sQ=q12;}
2855 }else if(term==3 && part==1){
2856 sQ=q12;
2857 if(c1==c3) {tQ=q13; fQ=q23;}
2858 else {tQ=q23; fQ=q13;}
2859 }else if(term==3 && part==2){
2860 sQ=q13;
2861 if(c1==c2) {tQ=q12; fQ=q23;}
2862 else {tQ=q23; fQ=q12;}
2863 }else if(term==4 && part==1){
2864 tQ=q12; sQ=q13; fQ=q23;
2865 }else if(term==4 && part==2){
2866 tQ=q13; sQ=q12; fQ=q23;
2867 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
2868 }else {// fQ=ss, sQ=ss, tQ=ss
2869 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
2870 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
2871 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
2872 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
2873 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
2874 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
2875 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
2876 }
2877 }else if(fi != 5){// all the rest except pi-k-p
2878 if(key1==key2){
2879 fQ=q12;
2880 if(c1==c2){
2881 // cases not explicity shown below are not possible
2882 if(term==1 || term==5) {sQ=q13; tQ=q23;}
2883 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
2884 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
2885 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
2886 else cout<<"problem!!!!!!!!!!!!!"<<endl;
2887 }else if(c3==0){
2888 if(c1==c3) {sQ=q13; tQ=q23;}
2889 else {sQ=q23; tQ=q13;}
2890 }else {//c3==1
2891 if(c1==c3) {tQ=q13; sQ=q23;}
2892 else {tQ=q23; sQ=q13;}
2893 }
2894 }else if(key1==key3){
2895 fQ=q13;
2896 if(c1==c3){
2897 // cases not explicity shown below are not possible
2898 if(term==1 || term==5) {sQ=q12; tQ=q23;}
2899 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
2900 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
2901 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
2902 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
2903 }else if(c2==0){
2904 if(c1==c2) {sQ=q12; tQ=q23;}
2905 else {sQ=q23; tQ=q12;}
2906 }else {//c2==1
2907 if(c1==c2) {tQ=q12; sQ=q23;}
2908 else {tQ=q23; sQ=q12;}
2909 }
2910 }else {// key2==key3
2911 fQ=q23;
2912 if(c2==c3){
2913 // cases not explicity shown below are not possible
2914 if(term==1 || term==5) {sQ=q12; tQ=q13;}
2915 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
2916 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
2917 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
2918 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
2919 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
2920 }else if(c1==0){
2921 if(c1==c2) {sQ=q12; tQ=q13;}
2922 else {sQ=q13; tQ=q12;}
2923 }else {//c1==1
2924 if(c1==c2) {tQ=q12; sQ=q13;}
2925 else {tQ=q13; sQ=q12;}
2926 }
2927 }
2928 }else {// pi-k-p
2929 if(key1==1){
2930 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
2931 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
2932 }else if(key1==10){
2933 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
2934 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
2935 }else {// key1==100
2936 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
2937 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
2938 }
2939
2940 }
2941
2942
2943}
2944//________________________________________________________________________
2945Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
2946
2947 Float_t qinv=1.0;
2948
2949 if(fi==0 || fi==3 || fi==5){// identical masses
2950 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));
2951 }else{// different masses
2952 Float_t px = track1[1] + track2[1];
2953 Float_t py = track1[2] + track2[2];
2954 Float_t pz = track1[3] + track2[3];
2955 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
2956 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
2957 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
2958
2959 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
2960 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
2961 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
2962 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
2963 qinv = sqrt(qinv);
2964 }
2965
2966 return qinv;
2967
2968}
2969//________________________________________________________________________
2970void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
2971
2972 Float_t p0 = track1[0] + track2[0];
2973 Float_t px = track1[1] + track2[1];
2974 Float_t py = track1[2] + track2[2];
2975 Float_t pz = track1[3] + track2[3];
2976
2977 Float_t mt = sqrt(p0*p0 - pz*pz);
2978 Float_t pt = sqrt(px*px + py*py);
2979
2980 Float_t v0 = track1[0] - track2[0];
2981 Float_t vx = track1[1] - track2[1];
2982 Float_t vy = track1[2] - track2[2];
2983 Float_t vz = track1[3] - track2[3];
2984
2985 qout = (px*vx + py*vy)/pt;
2986 qside = (px*vy - py*vx)/pt;
2987 qlong = (p0*vz - pz*v0)/mt;
2988}
2989//________________________________________________________________________
2990void AliChaoticity::SetWeightArrays(){
2991
2992 for(Int_t mb=0; mb<fMbins; mb++){
2993 for(Int_t edB=0; edB<kEDbins; edB++){
2994 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
2995 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
2996 //
2997 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
2998 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
2999 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3000
3001 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3002 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3003 }
3004 }
3005 }
3006 }
3007 }
3008 //
3009 }
3010 }
3011
3012
3013
3014 TFile *wFile = new TFile("WeightFile.root","READ");
3015 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3016 else cout<<"Good Weight File Found!"<<endl;
3017
3018 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3019 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3020 for(Int_t mb=0; mb<fMbins; mb++){
3021 for(Int_t edB=0; edB<kEDbins; edB++){
3022
3023 TString *name = new TString("Weight_Kt_");
3024 *name += tKbin;
3025 name->Append("_Ky_");
3026 *name += yKbin;
3027 name->Append("_M_");
3028 *name += mb;
3029 name->Append("_ED_");
3030 *name += edB;
3031
3032 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3033
3034 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3035 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3036 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3037
3038 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3039 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3040 }
3041 }
3042 }
3043 delete fTempHisto;
3044 }//ED
3045 }//mb
3046 }//ky
3047 }//kt
3048
3049 wFile->Close();
3050
3051
3052 cout<<"Done Setting Weights"<<endl;
3053
3054}
3055//________________________________________________________________________
3056void AliChaoticity::SetWeightArraysLEGO(TH3F *histos[kKbinsT][kCentBins]){
3057
3058 for(Int_t mb=0; mb<fMbins; mb++){
3059 for(Int_t edB=0; edB<kEDbins; edB++){
3060 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3061 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3062 //
3063 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3064 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3065 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3066
3067 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3068 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
3069 }
3070 }
3071 }
3072 }
3073 }
3074 //
3075 }
3076 }
3077}
3078//________________________________________________________________________
3079void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3080
3081 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3082 Float_t ky=0;// central rapidity
3083 //
3084 Float_t qOut=0,qSide=0,qLong=0;
3085 GetQosl(track1, track2, qOut, qSide, qLong);
3086 qOut = fabs(qOut);
3087 qSide = fabs(qSide);
3088 qLong = fabs(qLong);
3089 //
3090
3091 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3092 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3093 else {
3094 for(Int_t i=0; i<kKbinsT-1; i++){
3095 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3096 }
3097 }
3098 //
3099 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3100 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3101 else {
3102 for(Int_t i=0; i<kKbinsY-1; i++){
3103 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3104 }
3105 }
3106 //
3107 /////////
3108 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3109 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3110 else {
3111 for(Int_t i=0; i<kQbinsWeights-1; i++){
3112 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3113 }
3114 }
3115 //
3116 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3117 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3118 else {
3119 for(Int_t i=0; i<kQbinsWeights-1; i++){
3120 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3121 }
3122 }
3123 //
3124 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3125 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3126 else {
3127 for(Int_t i=0; i<kQbinsWeights-1; i++){
3128 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3129 }
3130 }
3131 //
3132
3133
3134 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3135 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3136
3137
3138 Float_t deltaW=0;
3139 // kt
3140 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3141 // Ky
3142 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3143 // Qo
3144 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstep;
3145 // Qs
3146 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstep;
3147 // Ql
3148 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstep;
3149
3150 //
3151 wgt = min + deltaW;
3152
3153
3154 ////
3155
3156 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3157 Float_t deltaWErr=0;
3158 // Kt
3159 /*
3160 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3161 // Ky
3162 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3163 // Qo
3164 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstep;
3165 // Qs
3166 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstep;
3167 // Ql
3168 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstep;
3169 */
3170 wgtErr = minErr + deltaWErr;
3171
3172
3173
3174}
3175//________________________________________________________________________
3176Float_t AliChaoticity::CoulCorr(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
3177
3178 if(rIndex >= kRVALUES) return 1.0;
3179 Int_t indexL = Int_t(fabs(qinv*1000. - 2.)/2.);// starts at 2MeV
3180 Int_t indexH = indexL+1;
3181 if(indexL >= kNlinesCoulFile-1) return 1.0;
3182
3183 Float_t slope=0;
3184 if(charge1==charge2){
3185 slope = (fCoulSS[rIndex][indexL] - fCoulSS[rIndex][indexH])/(fQCoul[indexL]-fQCoul[indexH]);
3186 return (slope*(qinv - fQCoul[indexL]) + fCoulSS[rIndex][indexL]);
3187 }else {
3188 slope = (fCoulOS[rIndex][indexL] - fCoulOS[rIndex][indexH])/(fQCoul[indexL]-fQCoul[indexH]);
3189 return (slope*(qinv - fQCoul[indexL]) + fCoulOS[rIndex][indexL]);
3190 }
3191}
3192//________________________________________________________________________
3193void AliChaoticity::SetCoulCorrections(){
3194
3195 // Set default values
3196 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3197 fQCoul[i]=0;
3198 for(Int_t j=0; j<kRVALUES; j++) {// radii columns
3199 fCoulSS[j][i]=-1;// should cause a crash in the sqrt() function if ifstream not successful below = (exit the code).
3200 fCoulOS[j][i]=-1;// should cause a crash in the sqrt() function if ifstream not successful below = (exit the code).
3201 }
3202 }
3203
3204 ifstream mystream("cc2_l100_all.txt");
3205 for(Int_t j=0; j<kRVALUES; j++) {// radii columns (3-10fm in cc2_l100_all.txt)
3206 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3207 mystream >> fQCoul[i];// Q2
3208 //
3209 mystream >> fCoulSS[j][i];
3210 mystream >> fCoulOS[j][i];
3211 }
3212 }
3213 mystream.close();
3214}
3215//________________________________________________________________________
3216void AliChaoticity::SetCoulCorrectionsLEGO(Float_t qPoints[kNlinesCoulFile], Float_t coulSS[kRVALUES][kNlinesCoulFile], Float_t coulOS[kRVALUES][kNlinesCoulFile]){
3217
3218 // Set default values
3219 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3220 fQCoul[i] = qPoints[i];
3221 for(Int_t j=0; j<kRVALUES; j++) {// radii columns
3222 fCoulSS[j][i] = coulSS[j][i];
3223 fCoulOS[j][i] = coulOS[j][i];
3224 }
3225 }
3226
3227}
3228//________________________________________________________________________
3229Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3230
3231 Float_t radius = Float_t(rIndex+3);
3232 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3233 Float_t coulCorr12 = CoulCorr(charge1, charge2, rIndex, qinv);
3234 if(charge1==charge2){
3235 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius/0.19733,2)))/coulCorr12);
3236 }else {
3237 return ((1-myDamp) + myDamp/coulCorr12);
3238 }
3239
3240}
3241//________________________________________________________________________
3242Float_t AliChaoticity::MCWeightr3(Int_t term, Int_t rIndex, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
3243 if(term==5) return 1.0;
3244 Float_t radius = 3. + rIndex;//starts at 5fm
3245 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3246 Float_t fc = sqrt(myDamp);
3247 Float_t coulCorr12 = CoulCorr(+1,+1, rIndex, q12);
3248 Float_t coulCorr13 = CoulCorr(+1,+1, rIndex, q13);
3249 Float_t coulCorr23 = CoulCorr(+1,+1, rIndex, q23);
3250
3251 if(term==1){
3252 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));
3253 c3QS += 2*exp(-pow(radius/0.19733,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3254 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3255 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius/0.19733,2)))/coulCorr12;
3256 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius/0.19733,2)))/coulCorr13;
3257 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius/0.19733,2)))/coulCorr23;
3258 w123 += pow(fc,3)*c3QS/(coulCorr12*coulCorr13*coulCorr23);
3259 return w123;
3260 }else if(term==2){
3261 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius/0.19733,2)))/coulCorr12);
3262 }else if(term==3){
3263 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius/0.19733,2)))/coulCorr13);
3264 }else if(term==4){
3265 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius/0.19733,2)))/coulCorr23);
3266 }else return 1.0;
3267
3268}
3269//________________________________________________________________________
3270Float_t AliChaoticity::MCWeightOSL(Float_t radius, Float_t myDamp, Float_t Q_out, Float_t Q_side, Float_t Q_long, Float_t qinv){
3271 Int_t rIndex = Int_t(radius+0.001-3);
3272 Float_t coulCorr12 = CoulCorr(+1,+1, rIndex, qinv);
3273
3274 return ((1-myDamp) + myDamp*(1 + exp(-pow(radius/0.19733,2)*(Q_out*Q_out+Q_side*Q_side+Q_long*Q_long)))/coulCorr12);
3275
3276}
3277//________________________________________________________________________
3278void AliChaoticity::SetMomResCorrections(){
3279
3280 for(Int_t i=0; i<kDENtypes; i++){
3281 for(Int_t j=0; j<kQbinsWeights; j++){
3282 fMomResWeights[i][j]=1.0;
3283 }
3284 }
3285
3286 TFile *momResFile = new TFile("MomResFile.root","READ");
3287 if(!momResFile->IsOpen()) {cout<<"No Momentum Resolution File!!!!!!!!!!"<<endl; return;}
3288 else cout<<"Good Momentum Resolution File Found!"<<endl;
3289
3290 TH2D *fMomResWeights_pp = (TH2D*)momResFile->Get("MomResHisto_pp");
3291
3292 for(Int_t i=0; i<fMomResWeights_pp->GetNbinsX(); i++){
3293 for(Int_t j=0; j<fMomResWeights_pp->GetNbinsY(); j++){
3294
3295 fMomResWeights[i][j]=Float_t(fMomResWeights_pp->GetBinContent(i+1,j+1));
3296
3297 }
3298 }
3299
3300 momResFile->Close();
3301
3302}
3303//________________________________________________________________________
3304void AliChaoticity::SetMomResCorrectionsLEGO(TH2D *histo){
3305 for(Int_t i=0; i<histo->GetNbinsX(); i++){
3306 for(Int_t j=0; j<histo->GetNbinsY(); j++){
3307 fMomResWeights[i][j]=Float_t(histo->GetBinContent(i+1,j+1));
3308 }
3309 }
3310}
3311//________________________________________________________________________
3312Float_t AliChaoticity::GetMomRes(Int_t denIndex, Float_t qinv){
3313
3314 Int_t qbinL = Int_t(fabs(qinv*1000. - 2.5)/5.);// starts at 2.5MeV (center of bin)
3315 Int_t qbinH = qbinL+1;
3316 if(qbinL >=kQbinsWeights-1) return 1.0;
3317 Float_t slope = (fMomResWeights[denIndex][qbinL] - fMomResWeights[denIndex][qbinH])/(-0.005);
3318 return (slope*(qinv - (0.005*(qbinL+0.5))) + fMomResWeights[denIndex][qbinL]);
3319
3320}
3321//________________________________________________________________________