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