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