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