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