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