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