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