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