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