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