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