7 #include "TObjString.h"
15 #include "TProfile2D.h"
20 #include "AliAnalysisTask.h"
21 #include "AliAnalysisManager.h"
24 #include "AliESDEvent.h"
25 #include "AliESDInputHandler.h"
26 #include "AliESDtrackCuts.h"
28 #include "AliAODEvent.h"
29 #include "AliAODInputHandler.h"
30 #include "AliAODMCParticle.h"
31 #include "AliAnalysisUtils.h"
33 #include "AliFourPion.h"
36 #define G_Coeff 0.006399 // 2*pi*alpha*M_pion
37 #define FmToGeV 0.19733 // conversion of Fm to GeV
38 #define kappa3 0.15 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
39 #define kappa4 0.32 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
40 #define kappa3Fit 0.1 // kappa3 for c4QS fit
41 #define kappa4Fit 0.5 // kappa4 for c4QS fit
43 // Author: Dhevan Gangadharan
47 //________________________________________________________________________
48 AliFourPion::AliFourPion():
62 fGenerateSignal(kFALSE),
63 fGeneratorOnly(kFALSE),
64 fTabulatePairs(kFALSE),
65 fLinearInterpolation(kTRUE),
66 fMixedChargeCut(kFALSE),
98 fQupperBoundWeights(0),
114 fMinSepPairEta(0.03),
115 fMinSepPairPhi(0.04),
136 fDefaultsCharSwitch(),
137 fLowQPairSwitch_E0E0(),
138 fLowQPairSwitch_E0E1(),
139 fLowQPairSwitch_E0E2(),
140 fLowQPairSwitch_E0E3(),
141 fLowQPairSwitch_E1E1(),
142 fLowQPairSwitch_E1E2(),
143 fLowQPairSwitch_E1E3(),
144 fLowQPairSwitch_E2E3(),
145 fNormQPairSwitch_E0E0(),
146 fNormQPairSwitch_E0E1(),
147 fNormQPairSwitch_E0E2(),
148 fNormQPairSwitch_E0E3(),
149 fNormQPairSwitch_E1E1(),
150 fNormQPairSwitch_E1E2(),
151 fNormQPairSwitch_E1E3(),
152 fNormQPairSwitch_E2E3(),
155 fWeightmuonCorrection(0x0)
157 // Default constructor
158 for(Int_t mb=0; mb<fMbins; mb++){
159 for(Int_t edB=0; edB<fEDbins; edB++){
160 for(Int_t c1=0; c1<2; c1++){
161 for(Int_t c2=0; c2<2; c2++){
162 for(Int_t term=0; term<2; term++){
164 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
166 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
167 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
168 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
169 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
170 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
171 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
176 for(Int_t c3=0; c3<2; c3++){
177 for(Int_t term=0; term<5; term++){
179 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
180 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
181 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
182 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
186 for(Int_t c4=0; c4<2; c4++){
187 for(Int_t term=0; term<13; term++){
189 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
190 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
191 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
192 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
200 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
201 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
202 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
203 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
210 // Initialize FSI histograms
211 for(Int_t i=0; i<12; i++){
217 // Initialize fNormWeight and fNormWeightErr to 0
218 for(Int_t i=0; i<3; i++){// Kt iterator
219 for(Int_t j=0; j<10; j++){// Mbin iterator
220 fNormWeight[i][j]=0x0;
226 //________________________________________________________________________
227 AliFourPion::AliFourPion(const Char_t *name)
228 : AliAnalysisTaskSE(name),
241 fGenerateSignal(kFALSE),
242 fGeneratorOnly(kFALSE),
243 fTabulatePairs(kFALSE),
244 fLinearInterpolation(kTRUE),
245 fMixedChargeCut(kFALSE),
259 fCentBinHighLimit(1),
277 fQupperBoundWeights(0),
293 fMinSepPairEta(0.03),
294 fMinSepPairPhi(0.04),
315 fDefaultsCharSwitch(),
316 fLowQPairSwitch_E0E0(),
317 fLowQPairSwitch_E0E1(),
318 fLowQPairSwitch_E0E2(),
319 fLowQPairSwitch_E0E3(),
320 fLowQPairSwitch_E1E1(),
321 fLowQPairSwitch_E1E2(),
322 fLowQPairSwitch_E1E3(),
323 fLowQPairSwitch_E2E3(),
324 fNormQPairSwitch_E0E0(),
325 fNormQPairSwitch_E0E1(),
326 fNormQPairSwitch_E0E2(),
327 fNormQPairSwitch_E0E3(),
328 fNormQPairSwitch_E1E1(),
329 fNormQPairSwitch_E1E2(),
330 fNormQPairSwitch_E1E3(),
331 fNormQPairSwitch_E2E3(),
334 fWeightmuonCorrection(0x0)
341 for(Int_t mb=0; mb<fMbins; mb++){
342 for(Int_t edB=0; edB<fEDbins; edB++){
343 for(Int_t c1=0; c1<2; c1++){
344 for(Int_t c2=0; c2<2; c2++){
345 for(Int_t term=0; term<2; term++){
347 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
349 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
350 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
351 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
352 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
353 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
354 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
358 for(Int_t c3=0; c3<2; c3++){
359 for(Int_t term=0; term<5; term++){
361 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
362 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
363 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
364 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
368 for(Int_t c4=0; c4<2; c4++){
369 for(Int_t term=0; term<13; term++){
371 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
372 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
373 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
374 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
382 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
383 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
384 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
385 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
392 // Initialize FSI histograms
393 for(Int_t i=0; i<12; i++){
398 // Initialize fNormWeight and fNormWeightErr to 0
399 for(Int_t i=0; i<3; i++){// Kt iterator
400 for(Int_t j=0; j<10; j++){// Mbin iterator
401 fNormWeight[i][j]=0x0;
406 DefineOutput(1, TList::Class());
408 //________________________________________________________________________
409 AliFourPion::AliFourPion(const AliFourPion &obj)
410 : AliAnalysisTaskSE(obj.fname),
414 fOutputList(obj.fOutputList),
415 fPIDResponse(obj.fPIDResponse),
418 fTempStruct(obj.fTempStruct),
419 fRandomNumber(obj.fRandomNumber),
421 fMCcase(obj.fMCcase),
422 fAODcase(obj.fAODcase),
423 fPbPbcase(obj.fPbPbcase),
424 fGenerateSignal(obj.fGenerateSignal),
425 fGeneratorOnly(obj.fGeneratorOnly),
426 fTabulatePairs(obj.fTabulatePairs),
427 fLinearInterpolation(obj.fLinearInterpolation),
428 fMixedChargeCut(obj.fMixedChargeCut),
431 ffcSqMRC(obj.ffcSqMRC),
432 fFilterBit(obj.fFilterBit),
433 fMaxChi2NDF(obj.fMaxChi2NDF),
434 fMinTPCncls(obj.fMinTPCncls),
435 fBfield(obj.fBfield),
437 fFSIindex(obj.fFSIindex),
440 fMultLimit(obj.fMultLimit),
441 fCentBinLowLimit(obj.fCentBinLowLimit),
442 fCentBinHighLimit(obj.fCentBinHighLimit),
443 fEventCounter(obj.fEventCounter),
444 fEventsToMix(obj.fEventsToMix),
445 fZvertexBins(obj.fZvertexBins),
450 fQLowerCut(obj.fQLowerCut),
453 fKupperBound(obj.fKupperBound),
454 fQupperBoundQ2(obj.fQupperBoundQ2),
455 fQupperBoundQ3(obj.fQupperBoundQ3),
456 fQupperBoundQ4(obj.fQupperBoundQ4),
457 fQbinsQ2(obj.fQbinsQ2),
458 fQbinsQ3(obj.fQbinsQ3),
459 fQbinsQ4(obj.fQbinsQ4),
460 fQupperBoundWeights(obj.fQupperBoundWeights),
468 fQstepWeights(obj.fQstepWeights),
470 fDampStart(obj.fDampStart),
471 fDampStep(obj.fDampStep),
472 fTPCTOFboundry(obj.fTPCTOFboundry),
473 fTOFboundry(obj.fTOFboundry),
474 fSigmaCutTPC(obj.fSigmaCutTPC),
475 fSigmaCutTOF(obj.fSigmaCutTOF),
476 fMinSepPairEta(obj.fMinSepPairEta),
477 fMinSepPairPhi(obj.fMinSepPairPhi),
478 fShareQuality(obj.fShareQuality),
479 fShareFraction(obj.fShareFraction),
480 fTrueMassP(obj.fTrueMassP),
481 fTrueMassPi(obj.fTrueMassPi),
482 fTrueMassK(obj.fTrueMassK),
483 fTrueMassKs(obj.fTrueMassKs),
484 fTrueMassLam(obj.fTrueMassLam),
485 fKtIndexL(obj.fKtIndexL),
486 fKtIndexH(obj.fKtIndexH),
487 fQoIndexL(obj.fQoIndexL),
488 fQoIndexH(obj.fQoIndexH),
489 fQsIndexL(obj.fQsIndexL),
490 fQsIndexH(obj.fQsIndexH),
491 fQlIndexL(obj.fQlIndexL),
492 fQlIndexH(obj.fQlIndexH),
493 fDummyB(obj.fDummyB),
494 fKT3transition(obj.fKT3transition),
495 fKT4transition(obj.fKT4transition),
498 fDefaultsCharSwitch(),
499 fLowQPairSwitch_E0E0(),
500 fLowQPairSwitch_E0E1(),
501 fLowQPairSwitch_E0E2(),
502 fLowQPairSwitch_E0E3(),
503 fLowQPairSwitch_E1E1(),
504 fLowQPairSwitch_E1E2(),
505 fLowQPairSwitch_E1E3(),
506 fLowQPairSwitch_E2E3(),
507 fNormQPairSwitch_E0E0(),
508 fNormQPairSwitch_E0E1(),
509 fNormQPairSwitch_E0E2(),
510 fNormQPairSwitch_E0E3(),
511 fNormQPairSwitch_E1E1(),
512 fNormQPairSwitch_E1E2(),
513 fNormQPairSwitch_E1E3(),
514 fNormQPairSwitch_E2E3(),
515 fMomResC2SC(obj.fMomResC2SC),
516 fMomResC2MC(obj.fMomResC2MC),
517 fWeightmuonCorrection(obj.fWeightmuonCorrection)
521 for(Int_t i=0; i<12; i++){
522 fFSIss[i]=obj.fFSIss[i];
523 fFSIos[i]=obj.fFSIos[i];
526 // Initialize fNormWeight and fNormWeightErr to 0
527 for(Int_t i=0; i<3; i++){// Kt iterator
528 for(Int_t j=0; j<10; j++){// Mbin iterator
529 fNormWeight[i][j]=0x0;
535 //________________________________________________________________________
536 AliFourPion &AliFourPion::operator=(const AliFourPion &obj)
538 // Assignment operator
544 fOutputList = obj.fOutputList;
545 fPIDResponse = obj.fPIDResponse;
548 fTempStruct = obj.fTempStruct;
549 fRandomNumber = obj.fRandomNumber;
551 fMCcase = obj.fMCcase;
552 fAODcase = obj.fAODcase;
553 fPbPbcase = obj.fPbPbcase;
554 fGenerateSignal = obj.fGenerateSignal;
555 fGeneratorOnly = obj.fGeneratorOnly;
556 fTabulatePairs = obj.fTabulatePairs;
557 fLinearInterpolation = obj.fLinearInterpolation;
558 fMixedChargeCut = obj.fMixedChargeCut;
561 ffcSqMRC = obj.ffcSqMRC;
562 fFilterBit = obj.fFilterBit;
563 fMaxChi2NDF = obj.fMaxChi2NDF;
564 fMinTPCncls = obj.fMinTPCncls;
565 fBfield = obj.fBfield;
567 fFSIindex = obj.fFSIindex;
570 fMultLimit = obj.fMultLimit;
571 fCentBinLowLimit = obj.fCentBinLowLimit;
572 fCentBinHighLimit = obj.fCentBinHighLimit;
573 fEventCounter = obj.fEventCounter;
574 fEventsToMix = obj.fEventsToMix;
575 fZvertexBins = obj.fZvertexBins;
579 fQLowerCut = obj.fQLowerCut;
580 fKupperBound = obj.fKupperBound;
581 fQupperBoundQ2 = obj.fQupperBoundQ2;
582 fQupperBoundQ3 = obj.fQupperBoundQ3;
583 fQupperBoundQ4 = obj.fQupperBoundQ4;
584 fQbinsQ2 = obj.fQbinsQ2;
585 fQbinsQ3 = obj.fQbinsQ3;
586 fQbinsQ4 = obj.fQbinsQ4;
587 fQupperBoundWeights = obj.fQupperBoundWeights;
589 fQstepWeights = obj.fQstepWeights;
590 fDampStart = obj.fDampStart;
591 fDampStep = obj.fDampStep;
592 fTPCTOFboundry = obj.fTPCTOFboundry;
593 fTOFboundry = obj.fTOFboundry;
594 fSigmaCutTPC = obj.fSigmaCutTPC;
595 fSigmaCutTOF = obj.fSigmaCutTOF;
596 fMinSepPairEta = obj.fMinSepPairEta;
597 fMinSepPairPhi = obj.fMinSepPairPhi;
598 fShareQuality = obj.fShareQuality;
599 fShareFraction = obj.fShareFraction;
600 fTrueMassP = obj.fTrueMassP;
601 fTrueMassPi = obj.fTrueMassPi;
602 fTrueMassK = obj.fTrueMassK;
603 fTrueMassKs = obj.fTrueMassKs;
604 fTrueMassLam = obj.fTrueMassLam;
605 fKtIndexL = obj.fKtIndexL;
606 fKtIndexH = obj.fKtIndexH;
607 fQoIndexL = obj.fQoIndexL;
608 fQoIndexH = obj.fQoIndexH;
609 fQsIndexL = obj.fQsIndexL;
610 fQsIndexH = obj.fQsIndexH;
611 fQlIndexL = obj.fQlIndexL;
612 fQlIndexH = obj.fQlIndexH;
613 fDummyB = obj.fDummyB;
614 fKT3transition = obj.fKT3transition;
615 fKT4transition = obj.fKT4transition;
616 fMomResC2SC = obj.fMomResC2SC;
617 fMomResC2MC = obj.fMomResC2MC;
618 fWeightmuonCorrection = obj.fWeightmuonCorrection;
620 for(Int_t i=0; i<12; i++){
621 fFSIss[i]=obj.fFSIss[i];
622 fFSIos[i]=obj.fFSIos[i];
624 for(Int_t i=0; i<3; i++){// Kt iterator
625 for(Int_t j=0; j<10; j++){// Mbin iterator
626 fNormWeight[i][j]=obj.fNormWeight[i][j];
632 //________________________________________________________________________
633 AliFourPion::~AliFourPion()
636 if(fAOD) delete fAOD;
637 //if(fESD) delete fESD;
638 if(fOutputList) delete fOutputList;
639 if(fPIDResponse) delete fPIDResponse;
641 if(fEvt) delete fEvt;
642 if(fTempStruct) delete [] fTempStruct;
643 if(fRandomNumber) delete fRandomNumber;
644 if(fMomResC2SC) delete fMomResC2SC;
645 if(fMomResC2MC) delete fMomResC2MC;
646 if(fWeightmuonCorrection) delete fWeightmuonCorrection;
648 for(Int_t j=0; j<kMultLimitPbPb; j++){
649 if(fLowQPairSwitch_E0E0[j]) delete [] fLowQPairSwitch_E0E0[j];
650 if(fLowQPairSwitch_E0E1[j]) delete [] fLowQPairSwitch_E0E1[j];
651 if(fLowQPairSwitch_E0E2[j]) delete [] fLowQPairSwitch_E0E2[j];
652 if(fLowQPairSwitch_E0E3[j]) delete [] fLowQPairSwitch_E0E3[j];
653 if(fLowQPairSwitch_E1E1[j]) delete [] fLowQPairSwitch_E1E1[j];
654 if(fLowQPairSwitch_E1E2[j]) delete [] fLowQPairSwitch_E1E2[j];
655 if(fLowQPairSwitch_E1E3[j]) delete [] fLowQPairSwitch_E1E3[j];
656 if(fLowQPairSwitch_E2E3[j]) delete [] fLowQPairSwitch_E2E3[j];
658 if(fNormQPairSwitch_E0E0[j]) delete [] fNormQPairSwitch_E0E0[j];
659 if(fNormQPairSwitch_E0E1[j]) delete [] fNormQPairSwitch_E0E1[j];
660 if(fNormQPairSwitch_E0E2[j]) delete [] fNormQPairSwitch_E0E2[j];
661 if(fNormQPairSwitch_E0E3[j]) delete [] fNormQPairSwitch_E0E3[j];
662 if(fNormQPairSwitch_E1E1[j]) delete [] fNormQPairSwitch_E1E1[j];
663 if(fNormQPairSwitch_E1E2[j]) delete [] fNormQPairSwitch_E1E2[j];
664 if(fNormQPairSwitch_E1E3[j]) delete [] fNormQPairSwitch_E1E3[j];
665 if(fNormQPairSwitch_E2E3[j]) delete [] fNormQPairSwitch_E2E3[j];
669 for(Int_t mb=0; mb<fMbins; mb++){
670 for(Int_t edB=0; edB<fEDbins; edB++){
671 for(Int_t c1=0; c1<2; c1++){
672 for(Int_t c2=0; c2<2; c2++){
673 for(Int_t term=0; term<2; term++){
675 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2;
677 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal;
678 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared;
679 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL;
680 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW;
681 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL;
682 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW;
684 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
685 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
688 for(Int_t c3=0; c3<2; c3++){
689 for(Int_t term=0; term<5; term++){
691 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3;
692 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3;
693 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor;
694 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted;
696 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm;
700 for(Int_t c4=0; c4<2; c4++){
701 for(Int_t term=0; term<13; term++){
703 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4;
704 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4;
705 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor;
706 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted;
708 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm;
715 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
716 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
717 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD;
718 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD;
726 for(Int_t i=0; i<12; i++){
727 if(fFSIss[i]) delete fFSIss[i];
728 if(fFSIos[i]) delete fFSIos[i];
730 for(Int_t i=0; i<3; i++){// Kt iterator
731 for(Int_t j=0; j<10; j++){// Mbin iterator
732 if(fNormWeight[i][j]) delete fNormWeight[i][j];
737 //________________________________________________________________________
738 void AliFourPion::ParInit()
740 cout<<"AliFourPion MyInit() call"<<endl;
741 cout<<"lego:"<<fLEGO<<" MCcase:"<<fMCcase<<" PbPbcase:"<<fPbPbcase<<" TabulatePairs:"<<fTabulatePairs<<" GenSignal:"<<fGenerateSignal<<" CentLow:"<<fCentBinLowLimit<<" CentHigh:"<<fCentBinHighLimit<<" RMax:"<<fRMax<<" fc^2:"<<ffcSq<<" FB:"<<fFilterBit<<" MaxChi2/NDF:"<<fMaxChi2NDF<<" MinTPCncls:"<<fMinTPCncls<<" MinPairSepEta:"<<fMinSepPairEta<<" MinPairSepPhi:"<<fMinSepPairPhi<<" NsigTPC:"<<fSigmaCutTPC<<" NsigTOF:"<<fSigmaCutTOF<<endl;
743 fRandomNumber = new TRandom3();
744 fRandomNumber->SetSeed(0);
751 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
752 fTOFboundry = 2.1;// TOF pid used below this momentum
754 ////////////////////////////////////////////////
756 fShareQuality = .5;// max
757 fShareFraction = .05;// max
758 ////////////////////////////////////////////////
761 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
762 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
766 if(fPbPbcase) {// PbPb
767 fMultLimit=kMultLimitPbPb;
770 fNormQcutLow = 0.15;// 0.15
771 fNormQcutHigh = 0.2;// 0.175
773 fMultLimit=kMultLimitpp;
780 fQLowerCut = 0.005;// was 0.005
784 fKmeanY[0] = 0;// central y
787 // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
789 fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
790 fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
791 fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
793 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
795 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
796 fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
797 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
799 // 2x1 (Kt: 0-0.35, 0.35-1.0)
801 fKstepT[0] = 0.35; fKstepT[1] = 0.65;
802 fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
803 fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
807 fQupperBoundWeights = 0.2;
808 fQupperBoundQ2 = 2.0;
809 fQupperBoundQ3 = 0.6;
810 fQupperBoundQ4 = 0.6;
811 fQbinsQ2 = fQupperBoundQ2/0.005;
812 fQbinsQ3 = fQupperBoundQ3/0.005;
813 fQbinsQ4 = fQupperBoundQ4/0.005;
814 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
815 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
817 fDampStart = 0.5;// was 0.3, then 0.5
822 fEC = new AliFourPionEventCollection **[fZvertexBins];
823 for(UShort_t i=0; i<fZvertexBins; i++){
825 fEC[i] = new AliFourPionEventCollection *[fMbins];
827 for(UShort_t j=0; j<fMbins; j++){
829 fEC[i][j] = new AliFourPionEventCollection(fEventsToMix+1, fMultLimit, kMCarrayLimit, fMCcase);
833 for(Int_t i=0; i<kMultLimitPbPb; i++) fDefaultsCharSwitch[i]='0';
834 for(Int_t i=0; i<kMultLimitPbPb; i++) {
835 fLowQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
836 fLowQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
837 fLowQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
838 fLowQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
839 fLowQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
840 fLowQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
841 fLowQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
842 fLowQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
844 fNormQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
845 fNormQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
846 fNormQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
847 fNormQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
848 fNormQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
849 fNormQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
850 fNormQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
851 fNormQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
854 fTempStruct = new AliFourPionTrackStruct[fMultLimit];
857 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
861 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
863 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
864 if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
865 if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
866 if(!fMCcase && !fTabulatePairs) SetMuonCorrections(fLEGO);// Read Muon corrections
869 /////////////////////////////////////////////
870 /////////////////////////////////////////////
873 //________________________________________________________________________
874 void AliFourPion::UserCreateOutputObjects()
879 ParInit();// Initialize my settings
882 fOutputList = new TList();
883 fOutputList->SetOwner();
885 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
886 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
887 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
888 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
889 fOutputList->Add(fVertexDist);
892 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
893 fOutputList->Add(fDCAxyDistPlus);
894 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
895 fOutputList->Add(fDCAzDistPlus);
896 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
897 fOutputList->Add(fDCAxyDistMinus);
898 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
899 fOutputList->Add(fDCAzDistMinus);
902 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
903 fOutputList->Add(fEvents1);
904 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
905 fOutputList->Add(fEvents2);
907 TH1F *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
908 fMultDist0->GetXaxis()->SetTitle("Multiplicity");
909 fOutputList->Add(fMultDist0);
911 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
912 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
913 fOutputList->Add(fMultDist1);
915 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
916 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
917 fOutputList->Add(fMultDist2);
919 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
920 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
921 fOutputList->Add(fMultDist3);
923 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
924 fOutputList->Add(fPtEtaDist);
926 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
927 fOutputList->Add(fPhiPtDist);
929 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
930 fOutputList->Add(fTOFResponse);
931 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
932 fOutputList->Add(fTPCResponse);
934 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",400,0,2);
935 fOutputList->Add(fRejectedPairs);
936 TH1F *fRejectedPairsWeighting = new TH1F("fAcceptedPairsWeighting","",400,0,2);
937 fOutputList->Add(fRejectedPairsWeighting);
938 TH1F *fTotalPairsWeighting = new TH1F("fTotalPairsWeighting","",400,0,2);
939 fOutputList->Add(fTotalPairsWeighting);
941 TH1F *fRejectedPairsMC = new TH1F("fRejectedPairsMC","",400,0,2);
942 fOutputList->Add(fRejectedPairsMC);
943 TH1F *fRejectedPairsWeightingMC = new TH1F("fAcceptedPairsWeightingMC","",400,0,2);
944 fOutputList->Add(fRejectedPairsWeightingMC);
945 TH1F *fTotalPairsWeightingMC = new TH1F("fTotalPairsWeightingMC","",400,0,2);
946 fOutputList->Add(fTotalPairsWeightingMC);
948 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
949 fOutputList->Add(fRejectedEvents);
951 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
952 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
953 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
954 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
955 TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
956 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
957 TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
958 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
959 TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
960 if(fMCcase) fOutputList->Add(fPairsPadRowNum);
961 TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
962 if(fMCcase) fOutputList->Add(fPairsPadRowDen);
966 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
967 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
968 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
969 if(fMCcase) fOutputList->Add(fAllOSPairs);
971 TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
972 if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
973 TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
974 if(fMCcase) fOutputList->Add(fAllSCPionPairs);
975 TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
976 if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
977 TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
978 if(fMCcase) fOutputList->Add(fAllMCPionPairs);
981 TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
982 if(fMCcase) fOutputList->Add(fMuonParents);
983 TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
984 if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
985 TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
986 if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
987 TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
988 if(fMCcase) fOutputList->Add(fPionCandidates);
992 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
993 fOutputList->Add(fAvgMult);
995 TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
996 fOutputList->Add(fTrackChi2NDF);
997 TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
998 fOutputList->Add(fTrackTPCncls);
1001 TH1D *fTPNRejects3pion1 = new TH1D("fTPNRejects3pion1","",fQbinsQ3,0,fQupperBoundQ3);
1002 fOutputList->Add(fTPNRejects3pion1);
1003 TH1D *fTPNRejects3pion2 = new TH1D("fTPNRejects3pion2","",fQbinsQ3,0,fQupperBoundQ3);
1004 fOutputList->Add(fTPNRejects3pion2);
1005 TH1D *fTPNRejects4pion1 = new TH1D("fTPNRejects4pion1","",fQbinsQ4,0,fQupperBoundQ4);
1006 fOutputList->Add(fTPNRejects4pion1);
1008 TH3D *fKT3DistTerm1 = new TH3D("fKT3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1009 TH3D *fKT3DistTerm5 = new TH3D("fKT3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1010 fOutputList->Add(fKT3DistTerm1);
1011 fOutputList->Add(fKT3DistTerm5);
1012 TH3D *fKT4DistTerm1 = new TH3D("fKT4DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1013 TH3D *fKT4DistTerm13 = new TH3D("fKT4DistTerm13","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1014 fOutputList->Add(fKT4DistTerm1);
1015 fOutputList->Add(fKT4DistTerm13);
1018 TProfile2D *fKT3AvgpT = new TProfile2D("fKT3AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
1019 fOutputList->Add(fKT3AvgpT);
1020 TProfile2D *fKT4AvgpT = new TProfile2D("fKT4AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
1021 fOutputList->Add(fKT4AvgpT);
1022 TH3D* fQ3AvgpT = new TH3D("fQ3AvgpT","", 2,-0.5,1.5, fQbinsQ3,0,fQupperBoundQ3, 180,0.1,1.0);
1023 fOutputList->Add(fQ3AvgpT);
1024 TH3D* fQ4AvgpT = new TH3D("fQ4AvgpT","", 2,-0.5,1.5, fQbinsQ4,0,fQupperBoundQ4, 180,0.1,1.0);
1025 fOutputList->Add(fQ4AvgpT);
1028 TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
1029 TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
1030 TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
1031 TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
1032 TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
1033 TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
1034 TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
1035 TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
1036 TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
1037 TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
1038 TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
1039 TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
1040 fOutputList->Add(fMCWeight3DTerm1SC);
1041 fOutputList->Add(fMCWeight3DTerm1SCden);
1042 fOutputList->Add(fMCWeight3DTerm2SC);
1043 fOutputList->Add(fMCWeight3DTerm2SCden);
1044 fOutputList->Add(fMCWeight3DTerm1MC);
1045 fOutputList->Add(fMCWeight3DTerm1MCden);
1046 fOutputList->Add(fMCWeight3DTerm2MC);
1047 fOutputList->Add(fMCWeight3DTerm2MCden);
1048 fOutputList->Add(fMCWeight3DTerm3MC);
1049 fOutputList->Add(fMCWeight3DTerm3MCden);
1050 fOutputList->Add(fMCWeight3DTerm4MC);
1051 fOutputList->Add(fMCWeight3DTerm4MCden);
1055 for(Int_t mb=0; mb<fMbins; mb++){
1056 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
1058 for(Int_t edB=0; edB<fEDbins; edB++){
1059 for(Int_t c1=0; c1<2; c1++){
1060 for(Int_t c2=0; c2<2; c2++){
1061 for(Int_t term=0; term<2; term++){
1063 TString *nameEx2 = new TString("TwoParticle_Charge1_");
1065 nameEx2->Append("_Charge2_");
1067 nameEx2->Append("_M_");
1069 nameEx2->Append("_ED_");
1071 nameEx2->Append("_Term_");
1074 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
1077 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1078 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2);
1079 TString *nameEx2QW=new TString(nameEx2->Data());
1080 nameEx2QW->Append("_QW");
1081 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1082 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW);
1083 TString *nameAvgP=new TString(nameEx2->Data());
1084 nameAvgP->Append("_AvgP");
1085 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, fQbinsQ2,0,fQupperBoundQ2, 0.,1.0,"");
1086 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP);
1088 TString *nameUnitMult=new TString(nameEx2->Data());
1089 nameUnitMult->Append("_UnitMult");
1090 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin = new TH2D(nameUnitMult->Data(),"Two Particle Distribution",21,0.5,21.5, fQbinsQ2,0,fQupperBoundQ2);
1091 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin);
1094 // Momentum resolution histos
1095 TString *nameIdeal = new TString(nameEx2->Data());
1096 nameIdeal->Append("_Ideal");
1097 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1098 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal);
1099 TString *nameSmeared = new TString(nameEx2->Data());
1100 nameSmeared->Append("_Smeared");
1101 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1102 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared);
1104 // Muon correction histos
1105 TString *nameMuonIdeal=new TString(nameEx2->Data());
1106 nameMuonIdeal->Append("_MuonIdeal");
1107 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1108 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal);
1109 TString *nameMuonSmeared=new TString(nameEx2->Data());
1110 nameMuonSmeared->Append("_MuonSmeared");
1111 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1112 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared);
1114 TString *nameMuonPionK2=new TString(nameEx2->Data());
1115 nameMuonPionK2->Append("_MuonPionK2");
1116 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2 = new TH2D(nameMuonPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1117 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2);
1119 TString *namePionPionK2=new TString(nameEx2->Data());
1120 namePionPionK2->Append("_PionPionK2");
1121 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2 = new TH2D(namePionPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1122 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2);
1125 TString *nameEx2MC=new TString(nameEx2->Data());
1126 nameEx2MC->Append("_MCqinv");
1127 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1128 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
1129 TString *nameEx2MCQW=new TString(nameEx2->Data());
1130 nameEx2MCQW->Append("_MCqinvQW");
1131 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1132 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
1134 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
1135 nameEx2PIDpurityDen->Append("_PIDpurityDen");
1136 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1137 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
1138 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
1139 nameEx2PIDpurityNum->Append("_PIDpurityNum");
1140 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH3D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",16,0.5,16.5, 20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1141 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
1143 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1144 nameEx2OSLB1->Append("_osl_b1");
1145 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1146 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL);
1147 nameEx2OSLB1->Append("_QW");
1148 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1149 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW);
1151 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1152 nameEx2OSLB2->Append("_osl_b2");
1153 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1154 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL);
1155 nameEx2OSLB2->Append("_QW");
1156 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1157 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW);
1163 // skip 3-particle if Tabulate6DPairs is true
1164 if(fTabulatePairs) continue;
1166 for(Int_t c3=0; c3<2; c3++){
1167 for(Int_t term=0; term<5; term++){
1169 TString *namePC3 = new TString("ThreeParticle_Charge1_");
1171 namePC3->Append("_Charge2_");
1173 namePC3->Append("_Charge3_");
1175 namePC3->Append("_M_");
1177 namePC3->Append("_ED_");
1179 namePC3->Append("_Term_");
1182 ///////////////////////////////////////
1183 // skip degenerate histograms
1184 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1185 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1186 /////////////////////////////////////////
1189 TString *nameNorm=new TString(namePC3->Data());
1190 nameNorm->Append("_Norm");
1191 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1192 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1195 TString *name1DQ=new TString(namePC3->Data());
1196 name1DQ->Append("_1D");
1197 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH1D(name1DQ->Data(),"", fQbinsQ3,0,fQupperBoundQ3);
1198 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1200 TString *nameKfactor=new TString(namePC3->Data());
1201 nameKfactor->Append("_Kfactor");
1202 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1203 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor);
1205 TString *nameKfactorW=new TString(namePC3->Data());
1206 nameKfactorW->Append("_KfactorWeighted");
1207 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted = new TProfile(nameKfactorW->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1208 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted);
1210 TString *nameMeanQinv=new TString(namePC3->Data());
1211 nameMeanQinv->Append("_MeanQinv");
1212 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv = new TProfile(nameMeanQinv->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,.2, "");
1213 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv);
1216 // Momentum resolution correction histos
1217 TString *nameMomResIdeal=new TString(namePC3->Data());
1218 nameMomResIdeal->Append("_Ideal");
1219 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1220 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1221 TString *nameMomResSmeared=new TString(namePC3->Data());
1222 nameMomResSmeared->Append("_Smeared");
1223 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1224 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared);
1225 // Muon correction histos
1226 TString *nameMuonIdeal=new TString(namePC3->Data());
1227 nameMuonIdeal->Append("_MuonIdeal");
1228 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1229 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal);
1230 TString *nameMuonSmeared=new TString(namePC3->Data());
1231 nameMuonSmeared->Append("_MuonSmeared");
1232 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1233 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared);
1235 TString *nameMuonPionK3=new TString(namePC3->Data());
1236 nameMuonPionK3->Append("_MuonPionK3");
1237 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH3D(nameMuonPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1238 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3);
1240 TString *namePionPionK3=new TString(namePC3->Data());
1241 namePionPionK3->Append("_PionPionK3");
1242 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH3D(namePionPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1243 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3);
1247 if(c1==c2 && c1==c3 && term==4 ){
1248 TString *nameTwoPartNorm=new TString(namePC3->Data());
1249 nameTwoPartNorm->Append("_TwoPartNorm");
1250 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1251 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm);
1253 TString *nameTwoPartNegNorm=new TString(namePC3->Data());
1254 nameTwoPartNegNorm->Append("_TwoPartNegNorm");
1255 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNegNorm = new TH2D(nameTwoPartNegNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1256 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNegNorm);
1258 TString *nameTwoPartNormErr=new TString(namePC3->Data());
1259 nameTwoPartNormErr->Append("_TwoPartNormErr");
1260 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNormErr = new TH2D(nameTwoPartNormErr->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1261 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNormErr);
1266 for(Int_t c4=0; c4<2; c4++){
1267 for(Int_t term=0; term<13; term++){
1269 TString *namePC4 = new TString("FourParticle_Charge1_");
1271 namePC4->Append("_Charge2_");
1273 namePC4->Append("_Charge3_");
1275 namePC4->Append("_Charge4_");
1277 namePC4->Append("_M_");
1279 namePC4->Append("_ED_");
1281 namePC4->Append("_Term_");
1284 ///////////////////////////////////////
1285 // skip degenerate histograms
1286 if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
1287 if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
1288 if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
1289 /////////////////////////////////////////
1291 TString *nameNorm=new TString(namePC4->Data());
1292 nameNorm->Append("_Norm");
1293 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1294 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
1296 TString *name1DQ=new TString(namePC4->Data());
1297 name1DQ->Append("_1D");
1298 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
1299 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
1301 TString *nameKfactor=new TString(namePC4->Data());
1302 nameKfactor->Append("_Kfactor");
1303 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1304 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
1306 TString *nameKfactorW=new TString(namePC4->Data());
1307 nameKfactorW->Append("_KfactorWeighted");
1308 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted = new TProfile(nameKfactorW->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1309 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted);
1311 if(c1==c2 && c1==c3 && c1==c4 && term==12 ){
1312 TString *nameTwoPartNorm=new TString(namePC4->Data());
1313 nameTwoPartNorm->Append("_TwoPartNorm");
1314 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1315 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
1317 TString *nameTwoPartNegNorm=new TString(namePC4->Data());
1318 nameTwoPartNegNorm->Append("_TwoPartNegNorm");
1319 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNegNorm = new TH2D(nameTwoPartNegNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1320 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNegNorm);
1322 TString *nameTwoPartNormErr=new TString(namePC4->Data());
1323 nameTwoPartNormErr->Append("_TwoPartNormErr");
1324 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNormErr = new TH2D(nameTwoPartNormErr->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1325 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNormErr);
1329 // Momentum resolution correction histos
1330 TString *nameMomResIdeal=new TString(namePC4->Data());
1331 nameMomResIdeal->Append("_Ideal");
1332 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1333 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
1334 TString *nameMomResSmeared=new TString(namePC4->Data());
1335 nameMomResSmeared->Append("_Smeared");
1336 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1337 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
1338 // Muon correction histos
1339 TString *nameMuonIdeal=new TString(namePC4->Data());
1340 nameMuonIdeal->Append("_MuonIdeal");
1341 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1342 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
1343 TString *nameMuonSmeared=new TString(namePC4->Data());
1344 nameMuonSmeared->Append("_MuonSmeared");
1345 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1346 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared);
1348 TString *nameMuonPionK4=new TString(namePC4->Data());
1349 nameMuonPionK4->Append("_MuonPionK4");
1350 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH3D(nameMuonPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1351 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4);
1353 TString *namePionPionK4=new TString(namePC4->Data());
1354 namePionPionK4->Append("_PionPionK4");
1355 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH3D(namePionPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1356 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4);
1374 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1375 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1376 for(Int_t mb=0; mb<fMbins; mb++){
1377 for(Int_t edB=0; edB<fEDbins; edB++){
1379 TString *nameNum = new TString("TPN_num_Kt_");
1381 nameNum->Append("_Ky_");
1383 nameNum->Append("_M_");
1385 nameNum->Append("_ED_");
1388 TString *nameDen = new TString("TPN_den_Kt_");
1390 nameDen->Append("_Ky_");
1392 nameDen->Append("_M_");
1394 nameDen->Append("_ED_");
1398 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1399 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
1401 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1402 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
1413 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1414 fOutputList->Add(fQsmearMean);
1415 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1416 fOutputList->Add(fQsmearSq);
1417 TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
1418 fOutputList->Add(fQ2Res);
1419 TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
1420 fOutputList->Add(fQ3Res);
1421 TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
1422 fOutputList->Add(fQ4Res);
1424 TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
1425 fOutputList->Add(DistQinv4pion);
1426 TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
1427 if(fMCcase) fOutputList->Add(DistQinvMC4pion);
1429 TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
1430 fOutputList->Add(fAvgQ12VersusQ3);
1431 TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
1432 fOutputList->Add(fAvgQ13VersusQ3);
1433 TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
1434 fOutputList->Add(fAvgQ23VersusQ3);
1436 TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5);
1437 fOutputList->Add(fDistPionParents4);
1439 TH2D *fDistTPCNclsFindable = new TH2D("fDistTPCNclsFindable","", 100,0,0.5, 201,-0.5,200.5);
1440 fDistTPCNclsFindable->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsFindable->GetYaxis()->SetTitle("Ncls Findable");
1441 fOutputList->Add(fDistTPCNclsFindable);
1442 TProfile *fProfileTPCNclsFindable = new TProfile("fProfileTPCNclsFindable","",100,0,0.5, 0,200, "");
1443 fProfileTPCNclsFindable->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsFindable->GetYaxis()->SetTitle("<Ncls Findable>");
1444 fOutputList->Add(fProfileTPCNclsFindable);
1446 TH2D *fDistTPCNclsCrossed = new TH2D("fDistTPCNclsCrossed","",100,0,0.5, 201,-0.5,200.5);
1447 fDistTPCNclsCrossed->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsCrossed->GetYaxis()->SetTitle("Ncls Crossed");
1448 fOutputList->Add(fDistTPCNclsCrossed);
1449 TProfile *fProfileTPCNclsCrossed = new TProfile("fProfileTPCNclsCrossed","",100,0,0.5, 0,200, "");
1450 fProfileTPCNclsCrossed->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsCrossed->GetYaxis()->SetTitle("<Ncls Crossed>");
1451 fOutputList->Add(fProfileTPCNclsCrossed);
1453 TH2D *fDistTPCNclsFindableRatio = new TH2D("fDistTPCNclsFindableRatio","",100,0,0.5, 100,0,1);
1454 fDistTPCNclsFindableRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsFindableRatio->GetYaxis()->SetTitle("Ncls / Ncls Findable");
1455 fOutputList->Add(fDistTPCNclsFindableRatio);
1456 TProfile *fProfileTPCNclsFindableRatio = new TProfile("fProfileTPCNclsFindableRatio","",100,0,0.5, 0,1, "");
1457 fProfileTPCNclsFindableRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsFindableRatio->GetYaxis()->SetTitle("<Ncls / Ncls Findable>");
1458 fOutputList->Add(fProfileTPCNclsFindableRatio);
1460 TH2D *fDistTPCNclsCrossedRatio = new TH2D("fDistTPCNclsCrossedRatio","",100,0,0.5, 100,0,1);
1461 fDistTPCNclsCrossedRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsCrossedRatio->GetYaxis()->SetTitle("Ncls / Ncls Crossed");
1462 fOutputList->Add(fDistTPCNclsCrossedRatio);
1463 TProfile *fProfileTPCNclsCrossedRatio = new TProfile("fProfileTPCNclsCrossedRatio","",100,0,0.5, 0,1, "");
1464 fProfileTPCNclsCrossedRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsCrossedRatio->GetYaxis()->SetTitle("<Ncls / Ncls Crossed>");
1465 fOutputList->Add(fProfileTPCNclsCrossedRatio);
1467 TH2D *fc4QSFitNum = new TH2D("fc4QSFitNum","",7,0.5,7.5, fQbinsQ4,0,fQupperBoundQ4);
1468 fOutputList->Add(fc4QSFitNum);
1469 TH2D *fc4QSFitDen = new TH2D("fc4QSFitDen","",7,0.5,7.5, fQbinsQ4,0,fQupperBoundQ4);
1470 fOutputList->Add(fc4QSFitDen);
1472 ////////////////////////////////////
1473 ///////////////////////////////////
1475 PostData(1, fOutputList);
1478 //________________________________________________________________________
1479 void AliFourPion::UserExec(Option_t *)
1482 // Called for each event
1483 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1486 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1488 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1489 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1493 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1494 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1495 if(!isSelected1 && !fMCcase) {return;}
1496 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1497 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1498 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1499 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1502 ///////////////////////////////////////////////////////////
1503 const AliAODVertex *primaryVertexAOD;
1504 AliCentrality *centrality;// for AODs and ESDs
1507 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1508 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1509 fPIDResponse = inputHandler->GetPIDResponse();
1512 TClonesArray *mcArray = 0x0;
1515 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1516 if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1517 cout<<"No MC particle branch found or Array too large!!"<<endl;
1525 Int_t positiveTracks=0, negativeTracks=0;
1526 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1528 Double_t vertex[3]={0};
1530 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1531 /////////////////////////////////////////////////
1534 Float_t centralityPercentile=0;
1535 Float_t cStep=5.0, cStart=0;
1538 if(fAODcase){// AOD case
1541 centrality = fAOD->GetCentrality();
1542 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1543 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1544 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
1545 cout<<"Centrality % = "<<centralityPercentile<<endl;
1549 ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
1551 // Pile-up rejection
1552 AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1553 if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1554 else AnaUtil->SetUseMVPlpSelection(kFALSE);
1555 Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1556 if(pileUpCase) return;
1558 ////////////////////////////////
1560 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1561 primaryVertexAOD = fAOD->GetPrimaryVertex();
1562 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1564 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1565 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1567 if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1569 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1571 fBfield = fAOD->GetMagneticField();
1573 for(Int_t i=0; i<fZvertexBins; i++){
1574 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1582 /////////////////////////////
1583 // Create Shuffled index list
1584 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1585 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1586 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1587 /////////////////////////////
1591 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1592 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1593 if (!aodtrack) continue;
1594 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1596 status=aodtrack->GetStatus();
1598 if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1599 ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1600 ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1601 if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
1602 if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1604 if(fFilterBit != 7){
1605 Bool_t goodTrackOtherFB = kFALSE;
1606 for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1607 AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1608 if(!aodtrack2) continue;
1609 if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1611 if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1614 if(!goodTrackOtherFB) continue;
1618 if(aodtrack->Pt() < fMinPt) continue;
1619 if(fabs(aodtrack->Eta()) > 0.8) continue;
1622 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1623 if(!goodMomentum) continue;
1624 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1627 Double_t dca2[2]={0};
1628 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1629 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1630 Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1632 fTempStruct[myTracks].fStatus = status;
1633 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1634 fTempStruct[myTracks].fId = aodtrack->GetID();
1636 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1637 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1638 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1639 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1640 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1641 fTempStruct[myTracks].fEta = aodtrack->Eta();
1642 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1643 fTempStruct[myTracks].fDCAXY = dca2[0];
1644 fTempStruct[myTracks].fDCAZ = dca2[1];
1645 fTempStruct[myTracks].fDCA = dca3d;
1646 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1647 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1651 if(fTempStruct[myTracks].fMom > fMaxPt) continue;// upper P bound
1652 //if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1653 //if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1657 fTempStruct[myTracks].fElectron = kFALSE;
1658 fTempStruct[myTracks].fPion = kFALSE;
1659 fTempStruct[myTracks].fKaon = kFALSE;
1660 fTempStruct[myTracks].fProton = kFALSE;
1662 Float_t nSigmaTPC[5];
1663 Float_t nSigmaTOF[5];
1664 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1665 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1666 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1667 Float_t signalTPC=0, signalTOF=0;
1668 Double_t integratedTimesTOF[10]={0};
1671 Bool_t DoPIDWorkAround=kTRUE;
1672 //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1673 if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1674 if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1675 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1676 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1677 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1678 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1679 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1681 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1682 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1683 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1684 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1685 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1686 signalTPC = aodtrack->GetTPCsignal();
1687 if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1688 fTempStruct[myTracks].fTOFhit = kTRUE;
1689 signalTOF = aodtrack->GetTOFsignal();
1690 aodtrack->GetIntegratedTimes(integratedTimesTOF);
1691 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1693 }else {// FilterBit 7 PID workaround
1695 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1696 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1697 if (!aodTrack2) continue;
1698 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1700 UInt_t status2=aodTrack2->GetStatus();
1702 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1703 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1704 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1705 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1706 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1708 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1709 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1710 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1711 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1712 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1713 signalTPC = aodTrack2->GetTPCsignal();
1715 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1716 fTempStruct[myTracks].fTOFhit = kTRUE;
1717 signalTOF = aodTrack2->GetTOFsignal();
1718 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1719 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1721 //if(aodTrack2->Pt()<0.2) cout<<aodTrack2->GetTPCNclsF()<<" "<<aodTrack2->GetTPCNCrossedRows()<<" "<<aodTrack2->GetTPCNcls()<<" "<<aodTrack2->GetTPCFoundFraction()<<endl;
1725 }// FilterBit 7 PID workaround
1729 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1730 if(fTempStruct[myTracks].fTOFhit) {
1731 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1735 // Use TOF if good hit and above threshold
1736 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1737 if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1738 if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1739 if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1740 if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1741 }else {// TPC info instead
1742 if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1743 if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1744 if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1745 if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1749 // Ensure there is only 1 candidate per track
1750 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1751 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1752 if(!fTempStruct[myTracks].fPion) continue;// only pions
1753 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1754 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1755 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1756 //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1757 ////////////////////////
1758 //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1762 if(fTempStruct[myTracks].fCharge==+1) {
1763 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1764 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1766 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1767 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1770 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1771 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1773 ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1774 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1776 ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1777 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1779 if(aodtrack->GetTPCNclsF() > 0){
1780 ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1781 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1784 ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1785 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1788 if(fTempStruct[myTracks].fPion) {// pions
1789 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1790 fTempStruct[myTracks].fKey = 1;
1791 }else if(fTempStruct[myTracks].fKaon){// kaons
1792 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1793 fTempStruct[myTracks].fKey = 10;
1795 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1796 fTempStruct[myTracks].fKey = 100;
1801 if(aodtrack->Charge() > 0) positiveTracks++;
1802 else negativeTracks++;
1804 if(fTempStruct[myTracks].fPion) pionCount++;
1805 if(fTempStruct[myTracks].fKaon) kaonCount++;
1806 if(fTempStruct[myTracks].fProton) protonCount++;
1810 if(fMCcase){// muon mothers
1811 AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1812 if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1813 AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1814 if(parent->IsPhysicalPrimary()){
1815 ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1816 }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1818 ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1821 //cout<<"kinkcount = "<<kinkcount<<" pionkinks = "<<pionkinks<<" primarypionkinks = "<<primarypionkinks<<endl;
1822 }else {// ESD tracks
1823 cout<<"ESDs not supported currently"<<endl;
1827 // Generator info only
1828 if(fMCcase && fGeneratorOnly){
1829 myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1830 for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1831 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1832 if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1834 AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1835 if(!mcParticle) continue;
1836 if(fabs(mcParticle->Eta())>0.8) continue;
1837 if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
1838 if(mcParticle->Pt() < fMinPt || mcParticle->Pt() > fMaxPt) continue;
1839 if(!mcParticle->IsPrimary()) continue;
1840 if(!mcParticle->IsPhysicalPrimary()) continue;
1841 if(abs(mcParticle->GetPdgCode())!=211) continue;
1843 fTempStruct[myTracks].fP[0] = mcParticle->Px();
1844 fTempStruct[myTracks].fP[1] = mcParticle->Py();
1845 fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1846 fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1848 fTempStruct[myTracks].fId = myTracks;// use my track counter
1849 fTempStruct[myTracks].fLabel = mctrackN;
1850 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1851 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1852 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1853 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1854 fTempStruct[myTracks].fEta = mcParticle->Eta();
1855 fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1856 fTempStruct[myTracks].fDCAXY = 0.;
1857 fTempStruct[myTracks].fDCAZ = 0.;
1858 fTempStruct[myTracks].fDCA = 0.;
1859 fTempStruct[myTracks].fPion = kTRUE;
1860 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1861 fTempStruct[myTracks].fKey = 1;
1869 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1873 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1874 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
1877 /////////////////////////////////////////
1878 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1879 if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1880 /////////////////////////////////////////
1883 ////////////////////////////////
1884 ///////////////////////////////
1885 // Mbin determination
1887 // Mbin set to Pion Count Only for pp!!!!!!!
1890 for(Int_t i=0; i<kMultBinspp; i++){
1891 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1892 // Mbin 0 has 1 pion
1895 for(Int_t i=0; i<fCentBins; i++){
1896 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1897 fMbin=i;// 0 = most central
1903 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1906 // can only be called after fMbin has been set
1907 // Radius parameter only matters for Monte-Carlo data
1911 Int_t rBinForTPNMomRes = 10;
1912 if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1913 else if(fMbin==1) {rBinForTPNMomRes=9;}
1914 else if(fMbin<=3) {rBinForTPNMomRes=8;}
1915 else if(fMbin<=5) {rBinForTPNMomRes=7;}
1916 else {rBinForTPNMomRes=6;}
1918 //////////////////////////////////////////////////
1919 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1920 //////////////////////////////////////////////////
1924 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1925 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1927 ////////////////////////////////////
1928 // Add event to buffer if > 0 tracks
1930 fEC[zbin][fMbin]->FIFOShift();
1931 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1932 (fEvt)->fNtracks = myTracks;
1933 (fEvt)->fFillStatus = 1;
1934 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1936 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1937 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1938 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1939 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1940 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1941 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1942 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1943 (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1944 (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1951 Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1952 Float_t qout=0, qside=0, qlong=0;
1954 Float_t q3=0, q3MC=0;
1955 Float_t q4=0, q4MC=0;
1956 Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1957 Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1958 Float_t pVect1[4]={0};
1959 Float_t pVect2[4]={0};
1960 Float_t pVect3[4]={0};
1961 Float_t pVect4[4]={0};
1962 Float_t pVect1MC[4]={0};
1963 Float_t pVect2MC[4]={0};
1964 Float_t pVect3MC[4]={0};
1965 Float_t pVect4MC[4]={0};
1966 Float_t Pparent1[4]={0};
1967 Float_t Pparent2[4]={0};
1968 Float_t Pparent3[4]={0};
1969 Float_t Pparent4[4]={0};
1970 Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1971 Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
1972 Float_t weight12CC[3]={0};
1973 Float_t weight13CC[3]={0};
1974 Float_t weight14CC[3]={0};
1975 Float_t weight23CC[3]={0};
1976 Float_t weight24CC[3]={0};
1977 Float_t weight34CC[3]={0};
1978 //Float_t weight12CC_e=0, weight13CC_e=0, weight14CC_e=0, weight23CC_e=0, weight24CC_e=0, weight34CC_e=0;
1979 Float_t weightTotal=0;//, weightTotalErr=0;
1980 Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0;
1981 Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1983 Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1984 Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1985 Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
1986 Bool_t Positive1stTripletWeights=kTRUE, Positive2ndTripletWeights=kTRUE;
1987 Float_t T12=0, T13=0, T14=0, T23=0, T24=0, T34=0;
1988 Int_t momBin12=1, momBin13=1, momBin14=1, momBin23=1, momBin24=1, momBin34=1;
1989 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr14=1.0, MomResCorr23=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
1991 AliAODMCParticle *mcParticle1=0x0;
1992 AliAODMCParticle *mcParticle2=0x0;
1995 ////////////////////
1996 //Int_t PairCount[7]={0};
1997 //Int_t NormPairCount[7]={0};
1998 Int_t KT3index=0, KT4index=0;
2000 // reset to defaults
2001 for(Int_t i=0; i<kMultLimitPbPb; i++) {
2002 fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2003 fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2004 fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2005 fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2006 fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2007 fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2008 fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2009 fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2011 fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2012 fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2013 fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2014 fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2015 fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2016 fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2017 fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2018 fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2022 //////////////////////////////////////////
2023 // make low-q pair storage and normalization-pair storage
2025 for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
2026 for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
2027 if(en1>1 && en1==en2) continue;
2029 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
2030 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2033 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2034 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2035 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2036 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2037 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2038 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2040 qinv12 = GetQinv(pVect1, pVect2);
2041 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2042 SetFillBins2(ch1, ch2, bin1, bin2);
2044 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
2045 if(ch1 == ch2 && !fGeneratorOnly){
2046 Int_t tempChGroup[2]={0,0};
2047 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2048 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2049 if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
2052 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2054 if(fMixedChargeCut && ch1 != ch2 && !fGeneratorOnly && !fMCcase){// remove +- low-q pairs to keep balance between ++ and +- contributions to multi-particle Q3,Q4 projections
2055 Int_t tempChGroup[2]={0,1};
2056 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2057 if(!AcceptPairPM((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2058 if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairsMC"))->Fill(qinv12);
2061 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2064 GetQosl(pVect1, pVect2, qout, qside, qlong);
2066 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
2067 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
2069 if((kT12 > 0.2) && (kT12 < 0.3)){
2070 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2071 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2073 if((kT12 > 0.6) && (kT12 < 0.7)){
2074 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2075 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2078 if( (fEvt+en1)->fNtracks%100==0){
2080 if(kT12>0.3) kTindex=1;
2081 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2082 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
2087 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
2088 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
2090 if((kT12 > 0.2) && (kT12 < 0.3)){
2091 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2092 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2094 if((kT12 > 0.6) && (kT12 < 0.7)){
2095 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2096 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2099 if( (fEvt+en1)->fNtracks%100==0){
2101 if(kT12>0.3) kTindex=1;
2102 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2103 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
2106 //////////////////////////////////////////
2107 if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
2109 Int_t kTbin=-1, kYbin=-1;
2111 for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}}
2112 for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
2113 if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2114 if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2115 if(fGenerateSignal && en2==0) {
2116 Int_t chGroup2[2]={ch1,ch2};
2117 Float_t WInput = MCWeight(chGroup2, fRMax, ffcSqMRC, qinv12, kT12);
2118 KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
2119 }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2122 //////////////////////////////////////////////////////////////////////////////
2124 if(qinv12 <= fQcut) {
2125 if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
2126 if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
2127 if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
2128 if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
2129 if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
2130 if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
2131 if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
2132 if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
2134 if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
2135 if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
2136 if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
2137 if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
2138 if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
2139 if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
2140 if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2141 if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2142 if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
2150 //cout<<PairCount[0]<<" "<<PairCount[1]<<" "<<PairCount[2]<<" "<<PairCount[3]<<" "<<PairCount[4]<<" "<<PairCount[5]<<" "<<PairCount[6]<<endl;
2151 //cout<<NormPairCount[0]<<" "<<NormPairCount[1]<<" "<<NormPairCount[2]<<" "<<NormPairCount[3]<<" "<<NormPairCount[4]<<" "<<NormPairCount[5]<<" "<<NormPairCount[6]<<endl;
2152 ///////////////////////////////////////////////////
2153 // Do not use pairs from events with too many pairs
2155 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2157 ///////////////////////////////////////////////////
2160 if(fTabulatePairs) return;
2162 /*TF1 *SCpairWeight = new TF1("SCpairWeight","[0] + [1]*x + [2]*exp(-[3]*x)",0,0.2);// same-charge pair weight for monte-carlo data without two-track cuts.
2163 SCpairWeight->FixParameter(0, 0.959);
2164 SCpairWeight->FixParameter(1, 0.278);
2165 SCpairWeight->FixParameter(2, -1.759);
2166 SCpairWeight->FixParameter(3, 115.107);*/
2168 ////////////////////////////////////////////////////
2169 ////////////////////////////////////////////////////
2170 // Normalization counting of 3- and 4-particle terms
2171 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2172 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2173 if(en2==0 && en3>2) continue;// not needed config
2174 if(en2==1 && en3==en2) continue;// not needed config
2175 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2176 if(en3==0 && en4>1) continue;// not needed config
2177 if(en3==1 && en4==3) continue;// not needed configs
2178 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2180 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2181 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2182 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2183 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2184 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2186 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2187 if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2188 else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2190 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2191 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2192 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2193 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2195 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2197 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2198 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2200 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2201 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2203 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2204 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2207 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2208 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2209 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2210 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2211 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2212 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2214 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2215 if(KT3<=fKT3transition) KT3index=0;
2218 if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2219 if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2220 if(en2==0 && en3==1 && en4==2) {
2221 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2222 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2223 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2227 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2229 if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2230 if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2231 if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2234 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2235 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2236 if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2238 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2239 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2240 if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2243 if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2244 if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2245 if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
2247 if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2248 if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2249 if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
2252 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2253 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2254 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2255 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2256 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2257 if(KT4<=fKT4transition) KT4index=0;
2260 Bool_t FillTerms[13]={kFALSE};
2261 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2263 for(int ft=0; ft<13; ft++) {
2264 if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.);
2280 ///////////////////////////////////////////////////////////////////////
2281 ///////////////////////////////////////////////////////////////////////
2282 ///////////////////////////////////////////////////////////////////////
2285 // Start the Main Correlation Analysis
2288 ///////////////////////////////////////////////////////////////////////
2292 ////////////////////////////////////////////////////
2293 ////////////////////////////////////////////////////
2294 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2295 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2296 if(en2==0 && en3>2) continue;// not needed config
2297 if(en2==1 && en3==en2) continue;// not needed config
2298 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2299 if(en3==0 && en4>1) continue;// not needed config
2300 if(en3==1 && en4==3) continue;// not needed configs
2301 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2303 Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
2305 /////////////////////////////////////////////////////////////
2306 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2307 pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2308 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2309 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2310 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2311 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2313 /////////////////////////////////////////////////////////////
2314 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2315 if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2316 else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2318 pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2319 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2320 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2321 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2322 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2323 qinv12 = GetQinv(pVect1, pVect2);
2324 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2325 SetFillBins2(ch1, ch2, bin1, bin2);
2327 if(kT12<=0.3) kTindex=0;
2330 FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2332 // two particle terms filled during tabulation of low-q pairs
2336 FilledMCpair12=kFALSE;
2338 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
2339 for(Int_t rstep=0; rstep<10; rstep++){
2340 Float_t coeff = (rstep)*0.2*(0.18/1.2);
2341 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2342 if(phi1 > 2*PI) phi1 -= 2*PI;
2343 if(phi1 < 0) phi1 += 2*PI;
2344 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2345 if(phi2 > 2*PI) phi2 -= 2*PI;
2346 if(phi2 < 0) phi2 += 2*PI;
2347 Float_t deltaphi = phi1 - phi2;
2348 if(deltaphi > PI) deltaphi -= PI;
2349 if(deltaphi < -PI) deltaphi += PI;
2351 if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2352 else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2357 // Check that label does not exceed stack size
2358 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2359 if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2360 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2361 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2362 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2363 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2364 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2365 qinv12MC = GetQinv(pVect1MC, pVect2MC);
2367 if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2368 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
2369 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2370 ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2373 // secondary contamination
2375 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2376 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2377 if(!mcParticle1 || !mcParticle2) continue;
2378 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2380 ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2381 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2382 ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2385 ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2386 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2387 ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2393 if(ENsum==6){// all mixed events
2394 Int_t chGroup2[2]={ch1,ch2};
2397 if(fFSIindex<=1) rForQW=10;
2398 else if(fFSIindex==2) rForQW=9;
2399 else if(fFSIindex==3) rForQW=8;
2400 else if(fFSIindex==4) rForQW=7;
2401 else if(fFSIindex==5) rForQW=6;
2402 else if(fFSIindex==6) rForQW=5;
2403 else if(fFSIindex==7) rForQW=4;
2404 else if(fFSIindex==8) rForQW=3;
2408 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
2409 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
2411 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2413 Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2414 if(PdgCodeSum==22) SCNumber=1;// e-e
2415 else if(PdgCodeSum==24) SCNumber=2;// e-mu
2416 else if(PdgCodeSum==222) SCNumber=3;// e-pi
2417 else if(PdgCodeSum==332) SCNumber=4;// e-k
2418 else if(PdgCodeSum==2223) SCNumber=5;// e-p
2419 else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2420 else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2421 else if(PdgCodeSum==334) SCNumber=8;// mu-k
2422 else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2423 else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2424 else if(PdgCodeSum==532) SCNumber=11;// pi-k
2425 else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2426 else if(PdgCodeSum==642) SCNumber=13;// k-k
2427 else if(PdgCodeSum==2533) SCNumber=14;// k-p
2428 else if(PdgCodeSum==4424) SCNumber=15;// p-p
2431 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2433 ///////////////////////
2434 // muon contamination
2435 Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2436 Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2437 pionParent1=kFALSE; pionParent2=kFALSE;
2438 FilledMCpair12=kTRUE;
2440 if(abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode)==13 || abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13){// muon check
2441 Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2442 if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2444 Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2445 Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2448 if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2449 Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2450 if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2452 Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2453 Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2457 parentQinv12 = GetQinv(Pparent1, Pparent2);
2459 if(pionParent1 || pionParent2){
2460 if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2461 Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2462 Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2463 for(Int_t term=1; term<=2; term++){
2464 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2465 Float_t Rvalue = 5+Riter;
2466 Float_t WInput = 1.0;
2468 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
2470 muonPionK12 = 1.0; pionPionK12=1.0;
2473 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2474 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2475 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2476 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2480 if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2481 else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2483 }// pion parent check
2487 // momentum resolution
2488 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2489 Float_t Rvalue = 5+Riter;
2490 Float_t WInput = MCWeight(chGroup2, Rvalue, ffcSqMRC, qinv12MC, 0.);
2491 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2492 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2493 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2494 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2503 /////////////////////////////////////////////////////////////
2504 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2506 if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2507 if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2509 if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2510 if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2512 if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2513 if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
2516 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2517 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2518 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2519 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2520 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2521 qinv13 = GetQinv(pVect1, pVect3);
2522 qinv23 = GetQinv(pVect2, pVect3);
2523 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2525 FilledMCtriplet123 = kFALSE;
2527 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2528 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2530 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2531 if(KT3<=fKT3transition) KT3index=0;
2534 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2535 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
2536 if(!fGenerateSignal && !fMCcase) {
2537 momBin12 = fMomResC2SC->GetYaxis()->FindBin(qinv12);
2538 momBin13 = fMomResC2SC->GetYaxis()->FindBin(qinv13);
2539 momBin23 = fMomResC2SC->GetYaxis()->FindBin(qinv23);
2540 if(momBin12 >= 20) momBin12 = 19;
2541 if(momBin13 >= 20) momBin13 = 19;
2542 if(momBin23 >= 20) momBin23 = 19;
2544 if(ch1==ch2) MomResCorr12 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin12);
2545 else MomResCorr12 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin12);
2546 if(ch1==ch3) MomResCorr13 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin13);
2547 else MomResCorr13 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin13);
2548 if(ch2==ch3) MomResCorr23 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin23);
2549 else MomResCorr23 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin23);
2553 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3);
2554 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23));
2555 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactorWeighted->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23), MomResCorr12*MomResCorr13*MomResCorr23);
2556 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2557 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2558 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2561 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2562 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2563 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2564 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2568 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
2569 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
2570 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2571 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2572 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2573 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2575 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
2576 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
2577 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2578 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2579 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2580 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2582 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
2583 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
2584 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2585 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2586 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2587 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2592 if(ENsum==6 && ch1==ch2 && ch1==ch3){
2593 Positive1stTripletWeights = kTRUE;
2595 GetWeight(pVect1, pVect2, weight12, weight12Err);
2596 GetWeight(pVect1, pVect3, weight13, weight13Err);
2597 GetWeight(pVect2, pVect3, weight23, weight23Err);
2599 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2600 if(fMbin==0 && bin1==0) {
2601 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2605 Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
2606 if(!fGenerateSignal && !fMCcase) {
2607 MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2608 MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2609 MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
2612 // no MRC, no Muon Correction
2613 weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2614 weight12CC[0] /= FSICorr12*ffcSq;
2615 weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2616 weight13CC[0] /= FSICorr13*ffcSq;
2617 weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2618 weight23CC[0] /= FSICorr23*ffcSq;
2619 if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2620 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2622 // no Muon Correction
2623 weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2624 weight12CC[1] /= FSICorr12*ffcSq;
2625 weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2626 weight13CC[1] /= FSICorr13*ffcSq;
2627 weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2628 weight23CC[1] /= FSICorr23*ffcSq;
2629 if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2630 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2633 weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2634 weight12CC[2] /= FSICorr12*ffcSq;
2635 weight12CC[2] *= MuonCorr12;
2636 weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2637 weight13CC[2] /= FSICorr13*ffcSq;
2638 weight13CC[2] *= MuonCorr13;
2639 weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2640 weight23CC[2] /= FSICorr23*ffcSq;
2641 weight23CC[2] *= MuonCorr23;
2643 if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
2644 if(fMbin==0 && bin1==0) {
2645 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
2647 if(weight12CC[2] < 0) weight12CC[2]=0;
2648 if(weight13CC[2] < 0) weight13CC[2]=0;
2649 if(weight23CC[2] < 0) weight23CC[2]=0;
2650 Positive1stTripletWeights = kFALSE;
2652 /////////////////////////////////////////////////////
2653 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2654 /////////////////////////////////////////////////////
2655 if(Positive1stTripletWeights){
2656 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
2657 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(4, q3, 1);
2659 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(4, q3, 1);
2662 // Full Weight reconstruction
2664 for(Int_t RcohIndex=0; RcohIndex<2; RcohIndex++){// Rcoh=0, then Rcoh=Rch
2665 for(Int_t GIndex=0; GIndex<50; GIndex++){
2666 Int_t FillBin = 5 + RcohIndex*50 + GIndex;
2667 Float_t G = 0.02*GIndex;
2669 T12 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight12CC[2])) / (2*pow(1-G,2));
2670 T13 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight13CC[2])) / (2*pow(1-G,2));
2671 T23 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight23CC[2])) / (2*pow(1-G,2));
2672 weightTotal = 2*G*(1-G)*(T12 + T13 + T23) + pow(1-G,2)*(T12*T12 + T13*T13 + T23*T23);
2673 weightTotal += 2*G*pow(1-G,2)*(T12*T13 + T12*T23 + T13*T23) + 2*pow(1-G,3)*T12*T13*T23;
2675 T12 = sqrt(weight12CC[2] / (1-G*G));
2676 T13 = sqrt(weight13CC[2] / (1-G*G));
2677 T23 = sqrt(weight23CC[2] / (1-G*G));
2678 weightTotal = (1-G*G)*(T12*T12 + T13*T13 + T23*T23);
2679 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3)) * T12*T13*T23;
2681 if(Positive1stTripletWeights){
2682 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(FillBin, q3, weightTotal);
2684 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(FillBin, q3, weightTotal);
2689 /*weight12CC_e = weight12Err*MomResCorr12 / FSICorr12 / ffcSq * MuonCorr12;
2690 weight13CC_e = weight13Err*MomResCorr13 / FSICorr13 / ffcSq * MuonCorr13;
2691 weight23CC_e = weight23Err*MomResCorr23 / FSICorr23 / ffcSq * MuonCorr23;
2692 if(weight12CC[2]*weight13CC[2]*weight23CC[2] > 0){
2693 weightTotalErr = pow(2 * sqrt(3) * weight12CC_e*weight13CC[2]*weight23CC[2] / sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]),2);
2695 weightTotalErr += pow(weight12CC_e,2) + pow(weight13CC_e,2) + pow(weight23CC_e,2);
2696 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNormErr->Fill(4, q3, weightTotalErr);*/
2698 }// 1st r3 den check
2703 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2704 ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
2706 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2707 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2708 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2709 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2710 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2711 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
2713 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt1);
2714 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt2);
2715 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt3);
2720 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2725 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2726 if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2728 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2729 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2730 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2731 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2732 qinv13MC = GetQinv(pVect1MC, pVect3MC);
2733 qinv23MC = GetQinv(pVect2MC, pVect3MC);
2735 q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2736 if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2738 Float_t TripletWeightTTC=1.0;// same-charge weights to mimic two-track depletion of same-charge pairs
2739 //if(ch1==ch2 && qinv12>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv12);
2740 //if(ch1==ch3 && qinv13>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv13);
2741 //if(ch2==ch3 && qinv23>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv23);
2743 Int_t chGroup3[3]={ch1,ch2,ch3};
2744 Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
2745 //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
2746 //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
2747 //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
2748 Float_t kTGroup3[3]={0};
2750 Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2753 if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2754 Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2755 if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2757 Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2758 Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2762 parentQinv13 = GetQinv(Pparent1, Pparent3);
2763 parentQinv23 = GetQinv(Pparent2, Pparent3);
2764 parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2766 if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2767 FilledMCtriplet123=kTRUE;
2768 if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2770 Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2771 //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2772 //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2773 //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2774 Float_t parentkTGroup3[3]={0};
2776 ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2777 ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2778 ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2780 for(Int_t term=1; term<=4; term++){
2782 else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2783 else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2784 else {if(!pionParent2 && !pionParent3) continue;}
2785 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2786 Float_t Rvalue = 5+Riter;
2787 Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2788 Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2789 Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
2790 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput*TripletWeightTTC);
2791 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput*TripletWeightTTC);
2792 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI*TripletWeightTTC);
2793 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI*TripletWeightTTC);
2795 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2796 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
2797 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2798 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
2802 }// pion parent check
2803 }// parentQ check (muon correction)
2805 // 3-pion momentum resolution
2806 for(Int_t term=1; term<=5; term++){
2807 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2808 Float_t Rvalue = 5+Riter;
2809 Float_t WInput = MCWeight3(term, Rvalue, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2810 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput*TripletWeightTTC);
2811 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput*TripletWeightTTC);
2815 }// 3rd particle label check
2816 }// MCcase and ENsum==6
2821 /////////////////////////////////////////////////////////////
2822 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2824 if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2825 if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2826 if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2829 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2830 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2831 if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2833 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2834 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2835 if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2838 if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2839 if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2840 if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2842 if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2843 if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2844 if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
2847 pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2848 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2849 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2850 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2851 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2852 qinv14 = GetQinv(pVect1, pVect4);
2853 qinv24 = GetQinv(pVect2, pVect4);
2854 qinv34 = GetQinv(pVect3, pVect4);
2855 q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
2857 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2858 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13);
2859 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23);
2860 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2863 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2864 if(KT4<=fKT4transition) KT4index=0;
2867 FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2868 FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2869 FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
2871 if(!fGenerateSignal && !fMCcase) {
2872 momBin14 = fMomResC2SC->GetYaxis()->FindBin(qinv14);
2873 momBin24 = fMomResC2SC->GetYaxis()->FindBin(qinv24);
2874 momBin34 = fMomResC2SC->GetYaxis()->FindBin(qinv34);
2875 if(momBin14 >= 20) momBin14 = 19;
2876 if(momBin24 >= 20) momBin24 = 19;
2877 if(momBin34 >= 20) momBin34 = 19;
2879 if(ch1==ch4) MomResCorr14 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin14);
2880 else MomResCorr14 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin14);
2881 if(ch2==ch4) MomResCorr24 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin24);
2882 else MomResCorr24 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin24);
2883 if(ch3==ch4) MomResCorr34 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin34);
2884 else MomResCorr34 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin34);
2887 Bool_t FillTerms[13]={kFALSE};
2888 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2890 for(int ft=0; ft<13; ft++) {
2891 Float_t FSIfactor = 1.0;
2892 Float_t MomResWeight = 1.0;
2894 FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2895 MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr14 * MomResCorr23 * MomResCorr24 * MomResCorr34;
2897 FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2898 MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr23;
2900 FSIfactor = 1/(FSICorr12);
2901 MomResWeight = MomResCorr12;
2903 FSIfactor = 1/(FSICorr12 * FSICorr34);
2904 MomResWeight = MomResCorr12 * MomResCorr34;
2905 }else {FSIfactor = 1.0; MomResWeight = 1.0;}
2907 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
2908 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
2909 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactorWeighted->Fill(q4, FSIfactor, MomResWeight);
2913 /////////////////////////////////////////////////////////////
2915 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2916 Positive2ndTripletWeights=kTRUE;
2918 GetWeight(pVect1, pVect4, weight14, weight14Err);
2919 GetWeight(pVect2, pVect4, weight24, weight24Err);
2920 GetWeight(pVect3, pVect4, weight34, weight34Err);
2922 Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
2923 if(!fGenerateSignal && !fMCcase) {
2924 MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2925 MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2926 MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
2929 // no MRC, no Muon Correction
2930 weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2931 weight14CC[0] /= FSICorr14*ffcSq;
2932 weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2933 weight24CC[0] /= FSICorr24*ffcSq;
2934 weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2935 weight34CC[0] /= FSICorr34*ffcSq;
2936 if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2937 weightTotal = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2938 weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2939 weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2941 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2943 // no Muon Correction
2944 weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2945 weight14CC[1] /= FSICorr14*ffcSq;
2946 weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2947 weight24CC[1] /= FSICorr24*ffcSq;
2948 weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2949 weight34CC[1] /= FSICorr34*ffcSq;
2950 if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2951 weightTotal = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]);
2952 weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]);
2953 weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]);
2955 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal);
2958 weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2959 weight14CC[2] /= FSICorr14*ffcSq;
2960 weight14CC[2] *= MuonCorr14;
2961 weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2962 weight24CC[2] /= FSICorr24*ffcSq;
2963 weight24CC[2] *= MuonCorr24;
2964 weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2965 weight34CC[2] /= FSICorr34*ffcSq;
2966 weight34CC[2] *= MuonCorr34;
2968 if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity
2969 if(fMbin==0 && bin1==0 && KT4index==0) {
2970 ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2])));
2972 if(weight14CC[2] < 0) weight14CC[2]=0;
2973 if(weight24CC[2] < 0) weight24CC[2]=0;
2974 if(weight34CC[2] < 0) weight34CC[2]=0;
2975 Positive2ndTripletWeights=kFALSE;
2977 /////////////////////////////////////////////////////
2978 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]);
2979 weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]);
2980 weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]);
2982 if(Positive1stTripletWeights && Positive2ndTripletWeights){
2983 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal);
2984 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(4, q4, 1);
2986 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNegNorm->Fill(4, q4, 1);
2988 // Full Weight reconstruction
2989 for(Int_t RcohIndex=0; RcohIndex<2; RcohIndex++){// Rcoh=0, then Rcoh=Rch
2990 for(Int_t GIndex=0; GIndex<50; GIndex++){
2991 Int_t FillBin = 5 + RcohIndex*50 + GIndex;
2992 Float_t G = 0.02*GIndex;
2993 if(RcohIndex==0){// Rcoh=0
2994 Float_t a = pow(1-G,2);
2995 Float_t b = 2*G*(1-G);
2996 T12 = (-b + sqrt(pow(b,2) + 4*a*weight12CC[2])) / (2*a);
2997 T13 = (-b + sqrt(pow(b,2) + 4*a*weight13CC[2])) / (2*a);
2998 T14 = (-b + sqrt(pow(b,2) + 4*a*weight14CC[2])) / (2*a);
2999 T23 = (-b + sqrt(pow(b,2) + 4*a*weight23CC[2])) / (2*a);
3000 T24 = (-b + sqrt(pow(b,2) + 4*a*weight24CC[2])) / (2*a);
3001 T34 = (-b + sqrt(pow(b,2) + 4*a*weight34CC[2])) / (2*a);
3002 weightTotal = 2*G*(1-G)*(T12 + T13 + T14 + T23 + T24 + T34) + pow(1-G,2)*(T12*T12 + T13*T13 + T14*T14 + T23*T23 + T24*T24 + T34*T34);// 2-pion
3003 weightTotal += 2*G*pow(1-G,3)*(T12*T34*T34 + T12*T12*T34 + T13*T24*T24 + T13*T13*T24 + T14*T23*T23 + T14*T14*T23);// 2-pair
3004 weightTotal += pow(1-G,4)*(pow(T12,2)*pow(T34,2) + pow(T13,2)*pow(T24,2) + pow(T14,2)*pow(T23,2));// 2-pair fully chaotic
3005 weightTotal += 2*G*pow(1-G,2)*(T12*T13 + T12*T23 + T13*T23 + T12*T14 + T12*T24 + T14*T24);// 3-pion
3006 weightTotal += 2*G*pow(1-G,2)*(T13*T14 + T13*T34 + T14*T34 + T23*T24 + T23*T34 + T24*T34);// 3-pion
3007 weightTotal += 2*pow(1-G,3)*(T12*T13*T23 + T12*T14*T24 + T13*T14*T34 + T23*T24*T34);// 3-pion fully chaotic
3008 weightTotal += 2*G*pow(1-G,3)*(T12*T14*T34 + T12*T14*T23 + T12*T23*T34 + T14*T23*T34);// 4-pion
3009 weightTotal += 2*G*pow(1-G,3)*(T12*T13*T34 + T12*T34*T24 + T12*T24*T13 + T13*T24*T34);// 4-pion
3010 weightTotal += 2*G*pow(1-G,3)*(T14*T13*T23 + T14*T13*T24 + T13*T23*T24 + T14*T24*T23);// 4-pion
3011 weightTotal += 2*pow(1-G,4)*(T12*T13*T24*T34 + T12*T14*T23*T34 + T13*T14*T23*T24);// 4-pion fully chaotic
3013 T12 = sqrt(weight12CC[2] / (1-G*G));
3014 T13 = sqrt(weight13CC[2] / (1-G*G));
3015 T14 = sqrt(weight14CC[2] / (1-G*G));
3016 T23 = sqrt(weight23CC[2] / (1-G*G));
3017 T24 = sqrt(weight24CC[2] / (1-G*G));
3018 T34 = sqrt(weight34CC[2] / (1-G*G));
3019 weightTotal = (1-G*G)*(T12*T12 + T13*T13 + T14*T14 + T23*T23 + T24*T24 + T34*T34);// 2-pion
3020 weightTotal += (4*G*pow(1-G,3)+pow(1-G,4))*(pow(T12,2)*pow(T34,2) + pow(T13,2)*pow(T24,2) + pow(T14,2)*pow(T23,2));// 2-pair
3021 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3))*(T12*T13*T23);// 3-pion
3022 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3))*(T12*T14*T24);// 3-pion
3023 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3))*(T13*T14*T34);// 3-pion
3024 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3))*(T23*T24*T34);// 3-pion
3025 weightTotal += (8*G*pow(1-G,3) + 2*pow(1-G,4))*(T12*T13*T24*T34);// 4-pion
3026 weightTotal += (8*G*pow(1-G,3) + 2*pow(1-G,4))*(T12*T14*T23*T34);// 4-pion
3027 weightTotal += (8*G*pow(1-G,3) + 2*pow(1-G,4))*(T13*T14*T23*T24);// 4-pion
3029 if(Positive1stTripletWeights && Positive2ndTripletWeights){
3030 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(FillBin, q4, weightTotal);
3032 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNegNorm->Fill(FillBin, q4, weightTotal);
3037 /*weight14CC_e = weight14Err*MomResCorr14 / FSICorr14 / ffcSq * MuonCorr14;
3038 weight24CC_e = weight24Err*MomResCorr24 / FSICorr24 / ffcSq * MuonCorr24;
3039 weight34CC_e = weight34Err*MomResCorr34 / FSICorr34 / ffcSq * MuonCorr34;
3040 if(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2] > 0){
3041 weightTotalErr = pow( 6 * 2 * weight12CC_e*weight13CC[2]*weight24CC[2]*weight34CC[2] / sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]),2);
3043 if(weight12CC[2]*weight13CC[2]*weight23CC[2] > 0){
3044 weightTotalErr += pow( 8 * sqrt(3) * weight12CC_e*weight13CC[2]*weight23CC[2] / sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]),2);
3046 weightTotalErr += 2*(pow(weight12CC_e*weight34CC[2],2) + pow(weight13CC_e*weight24CC[2],2) + pow(weight14CC_e*weight23CC[2],2));
3047 weightTotalErr += pow(weight12CC_e,2) + pow(weight13CC_e,2) + pow(weight14CC_e,2) + pow(weight23CC_e,2) + pow(weight24CC_e,2) + pow(weight34CC_e,2);
3048 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNormErr->Fill(4, q4, weightTotalErr);
3050 if(fMbin==0 && KT4index==0){
3051 for(Int_t Rindex=0; Rindex<7; Rindex++){
3052 Float_t R = (6. + Rindex)/FmToGeV;
3053 Float_t arg12=qinv12*R;
3054 Float_t arg13=qinv13*R;
3055 Float_t arg14=qinv14*R;
3056 Float_t arg23=qinv23*R;
3057 Float_t arg24=qinv24*R;
3058 Float_t arg34=qinv34*R;
3059 // Exchange Amplitudes
3060 Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
3061 Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
3062 Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
3063 Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
3064 Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
3065 Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3Fit/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4Fit/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));
3067 Float_t TotalCorrelation = 1 + 2*(EA12*EA13*EA24*EA34 + EA12*EA14*EA23*EA34 + EA13*EA14*EA23*EA24);
3068 ((TH2D*)fOutputList->FindObject("fc4QSFitNum"))->Fill(Rindex+1, q4, TotalCorrelation);
3069 ((TH2D*)fOutputList->FindObject("fc4QSFitDen"))->Fill(Rindex+1, q4);
3073 /////////////////////////////////////////////////////////////
3075 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==0){
3076 ((TH3D*)fOutputList->FindObject("fKT4DistTerm1"))->Fill(fMbin+1, KT4, q4);
3078 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
3079 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
3080 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
3081 Float_t pt4=sqrt(pow(pVect4[1],2)+pow(pVect4[2],2));
3082 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt1);
3083 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt2);
3084 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt3);
3085 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt4);
3087 ((TH3D*)fOutputList->FindObject("fQ4AvgpT"))->Fill(KT4index, q4, pt1);
3088 ((TH3D*)fOutputList->FindObject("fQ4AvgpT"))->Fill(KT4index, q4, pt2);
3089 ((TH3D*)fOutputList->FindObject("fQ4AvgpT"))->Fill(KT4index, q4, pt3);
3090 ((TH3D*)fOutputList->FindObject("fQ4AvgpT"))->Fill(KT4index, q4, pt4);
3094 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT4DistTerm13"))->Fill(fMbin+1, KT4, q4);
3097 // momenumtum resolution and muon corrections
3098 if(fMCcase && ENsum==6 && FilledMCtriplet123){// for momentum resolution and muon correction
3099 if((fEvt+en4)->fTracks[l].fLabel < (fEvt+en4)->fMCarraySize){
3101 pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
3102 pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
3103 pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
3104 pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
3105 qinv14MC = GetQinv(pVect1MC, pVect4MC);
3106 qinv24MC = GetQinv(pVect2MC, pVect4MC);
3107 qinv34MC = GetQinv(pVect3MC, pVect4MC);
3109 q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) + pow(qinv24MC,2) + pow(qinv34MC,2));
3110 if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC);
3111 if(ch1==ch2 && ch1==ch3 && ch1==ch4) {
3112 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC);
3113 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC);
3114 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC);
3117 Float_t QuadWeightTTC=1.0;// same-charge weights to mimic two-track depletion of same-charge pairs
3118 //if(ch1==ch2 && qinv12>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv12);
3119 //if(ch1==ch3 && qinv13>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv13);
3120 //if(ch1==ch4 && qinv14>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv14);
3121 //if(ch2==ch3 && qinv23>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv23);
3122 //if(ch2==ch4 && qinv24>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv24);
3123 //if(ch3==ch4 && qinv34>0.006) QuadWeightTTC *= SCpairWeight->Eval(qinv34);
3126 Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
3127 Float_t QinvMCGroup4[6]={qinv12MC, qinv13MC, qinv14MC, qinv23MC, qinv24MC, qinv34MC};
3128 Float_t kTGroup4[6]={0};
3130 Pparent4[0]=pVect4MC[0]; Pparent4[1]=pVect4MC[1]; Pparent4[2]=pVect4MC[2]; Pparent4[3]=pVect4MC[3];
3132 if(abs((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPdgCode)==13){// muon check
3133 Int_t MotherLabel4 = (fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fMotherLabel;
3134 if(abs((fEvt+en4)->fMCtracks[MotherLabel4].fPdgCode)==211) {
3136 Pparent4[1] = (fEvt+en4)->fMCtracks[MotherLabel4].fPx; Pparent4[2] = (fEvt+en4)->fMCtracks[MotherLabel4].fPy; Pparent4[3] = (fEvt+en4)->fMCtracks[MotherLabel4].fPz;
3137 Pparent4[0] = sqrt(pow(Pparent4[1],2)+pow(Pparent4[2],2)+pow(Pparent4[3],2)+pow(fTrueMassPi,2));
3141 parentQinv14 = GetQinv(Pparent1, Pparent4);
3142 parentQinv24 = GetQinv(Pparent2, Pparent4);
3143 parentQinv34 = GetQinv(Pparent3, Pparent4);
3144 Float_t parentQ4 = sqrt(pow(parentQ3,2) + pow(parentQinv14,2) + pow(parentQinv24,2) + pow(parentQinv34,2));
3146 if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){
3147 if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon
3149 if(pionParent1) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(1);
3150 if(pionParent2) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(2);
3151 if(pionParent3) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(3);
3152 if(pionParent4) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(4);
3153 Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34};
3154 Float_t parentkTGroup4[6]={0};
3156 for(Int_t term=1; term<=12; term++){
3158 else if(term==2) {if(!pionParent1 && !pionParent2 && !pionParent3) continue;}
3159 else if(term==3) {if(!pionParent1 && !pionParent2 && !pionParent4) continue;}
3160 else if(term==4) {if(!pionParent1 && !pionParent3 && !pionParent4) continue;}
3161 else if(term==5) {if(!pionParent2 && !pionParent3 && !pionParent4) continue;}
3162 else if(term==6) {if(!pionParent1 && !pionParent2) continue;}
3163 else if(term==7) {if(!pionParent1 && !pionParent3) continue;}
3164 else if(term==8) {if(!pionParent1 && !pionParent4) continue;}
3165 else if(term==9) {if(!pionParent2 && !pionParent3) continue;}
3166 else if(term==10) {if(!pionParent2 && !pionParent4) continue;}
3167 else if(term==11) {if(!pionParent3 && !pionParent4) continue;}
3169 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
3170 Float_t Rvalue = 5+Riter;
3171 Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
3172 Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
3173 Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
3174 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(1, Rvalue, q4MC, WInput*QuadWeightTTC);
3175 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ4, WInput*QuadWeightTTC);
3176 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(1, Rvalue, q4MC, WInputFSI*QuadWeightTTC);
3177 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(1, Rvalue, parentQ4, WInputParentFSI*QuadWeightTTC);
3179 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(2, Rvalue, q4MC, QuadWeightTTC);
3180 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ4, QuadWeightTTC);
3181 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(2, Rvalue, q4MC, QuadWeightTTC);
3182 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(2, Rvalue, parentQ4, QuadWeightTTC);
3186 }// pion parent check
3187 }// parentQ check (muon correction)
3189 // 4-pion momentum resolution
3190 for(Int_t term=1; term<=13; term++){
3191 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
3192 Float_t Rvalue = 5+Riter;
3193 Float_t WInput = MCWeight4(term, Rvalue, ffcSqMRC, chGroup4, QinvMCGroup4, kTGroup4);
3194 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fIdeal->Fill(Rvalue, q4MC, WInput*QuadWeightTTC);
3195 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fSmeared->Fill(Rvalue, q4, WInput*QuadWeightTTC);
3199 }// label check particle 4
3217 // Post output data.
3218 PostData(1, fOutputList);
3221 //________________________________________________________________________
3222 void AliFourPion::Terminate(Option_t *)
3224 // Called once at the end of the query
3229 //________________________________________________________________________
3230 Bool_t AliFourPion::AcceptPair(AliFourPionTrackStruct first, AliFourPionTrackStruct second)
3233 if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
3235 // propagate through B field to r=1m
3236 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
3237 if(phi1 > 2*PI) phi1 -= 2*PI;
3238 if(phi1 < 0) phi1 += 2*PI;
3239 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
3240 if(phi2 > 2*PI) phi2 -= 2*PI;
3241 if(phi2 < 0) phi2 += 2*PI;
3243 Float_t deltaphi = phi1 - phi2;
3244 if(deltaphi > PI) deltaphi -= 2*PI;
3245 if(deltaphi < -PI) deltaphi += 2*PI;
3246 deltaphi = fabs(deltaphi);
3248 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3251 // propagate through B field to r=1.6m
3252 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
3253 if(phi1 > 2*PI) phi1 -= 2*PI;
3254 if(phi1 < 0) phi1 += 2*PI;
3255 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
3256 if(phi2 > 2*PI) phi2 -= 2*PI;
3257 if(phi2 < 0) phi2 += 2*PI;
3259 deltaphi = phi1 - phi2;
3260 if(deltaphi > PI) deltaphi -= 2*PI;
3261 if(deltaphi < -PI) deltaphi += 2*PI;
3262 deltaphi = fabs(deltaphi);
3264 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3270 /* Int_t ncl1 = first.fClusterMap.GetNbits();
3271 Int_t ncl2 = second.fClusterMap.GetNbits();
3272 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
3273 Double_t shfrac = 0; Double_t qfactor = 0;
3274 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
3275 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
3276 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
3280 else {sumQ--; sumCls+=2;}
3282 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
3287 qfactor = sumQ*1.0/sumCls;
3288 shfrac = sumSha*1.0/sumCls;
3291 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
3298 //________________________________________________________________________
3299 Bool_t AliFourPion::AcceptPairPM(AliFourPionTrackStruct first, AliFourPionTrackStruct second)
3300 {// optional pair cuts for +- pairs
3302 if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
3304 // propagate through B field to r=1m
3305 Float_t phi1 = first.fPhi - asin(1.*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
3306 if(phi1 > 2*PI) phi1 -= 2*PI;
3307 if(phi1 < 0) phi1 += 2*PI;
3308 Float_t phi2 = second.fPhi - asin(1.*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
3309 if(phi2 > 2*PI) phi2 -= 2*PI;
3310 if(phi2 < 0) phi2 += 2*PI;
3312 Float_t deltaphi = phi1 - phi2;
3313 if(deltaphi > PI) deltaphi -= 2*PI;
3314 if(deltaphi < -PI) deltaphi += 2*PI;
3315 deltaphi = fabs(deltaphi);
3317 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3320 // propagate through B field to r=1.6m
3321 phi1 = first.fPhi - asin(1.*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
3322 if(phi1 > 2*PI) phi1 -= 2*PI;
3323 if(phi1 < 0) phi1 += 2*PI;
3324 phi2 = second.fPhi - asin(1.*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
3325 if(phi2 > 2*PI) phi2 -= 2*PI;
3326 if(phi2 < 0) phi2 += 2*PI;
3328 deltaphi = phi1 - phi2;
3329 if(deltaphi > PI) deltaphi -= 2*PI;
3330 if(deltaphi < -PI) deltaphi += 2*PI;
3331 deltaphi = fabs(deltaphi);
3333 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3338 //________________________________________________________________________
3339 Float_t AliFourPion::Gamov(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
3341 Float_t arg = G_Coeff/qinv;
3343 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
3344 else {return (exp(-arg)-1)/(-arg);}
3347 //________________________________________________________________________
3348 void AliFourPion::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
3352 for (Int_t i = i1; i < i2+1; i++) {
3353 j = (Int_t) (gRandom->Rndm() * a);
3361 //________________________________________________________________________
3362 Float_t AliFourPion::GetQinv(Float_t track1[], Float_t track2[]){
3364 Float_t qinv = sqrt( fabs(pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2)) );
3368 //________________________________________________________________________
3369 void AliFourPion::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3371 Float_t p0 = track1[0] + track2[0];
3372 Float_t px = track1[1] + track2[1];
3373 Float_t py = track1[2] + track2[2];
3374 Float_t pz = track1[3] + track2[3];
3376 Float_t mt = sqrt(p0*p0 - pz*pz);
3377 Float_t pt = sqrt(px*px + py*py);
3379 Float_t v0 = track1[0] - track2[0];
3380 Float_t vx = track1[1] - track2[1];
3381 Float_t vy = track1[2] - track2[2];
3382 Float_t vz = track1[3] - track2[3];
3384 qout = (px*vx + py*vy)/pt;
3385 qside = (px*vy - py*vx)/pt;
3386 qlong = (p0*vz - pz*v0)/mt;
3388 //________________________________________________________________________
3389 void AliFourPion::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]){
3392 cout<<"LEGO call to SetWeightArrays"<<endl;
3394 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3395 for(Int_t mb=0; mb<fCentBins; mb++){
3396 fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
3397 fNormWeight[tKbin][mb]->SetDirectory(0);
3403 TFile *wFile = new TFile("WeightFile.root","READ");
3404 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3405 else cout<<"Good Weight File Found!"<<endl;
3407 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3408 for(Int_t mb=0; mb<fCentBins; mb++){
3410 TString *name = new TString("Weight_Kt_");
3412 name->Append("_Ky_0");
3413 name->Append("_M_");
3415 name->Append("_ED_0");
3418 fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
3419 fNormWeight[tKbin][mb]->SetDirectory(0);
3428 cout<<"Done reading weight file"<<endl;
3431 //________________________________________________________________________
3432 void AliFourPion::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3434 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3436 Float_t qOut=0,qSide=0,qLong=0;
3437 GetQosl(track1, track2, qOut, qSide, qLong);
3439 qSide = fabs(qSide);
3440 qLong = fabs(qLong);
3441 Float_t wd=0, xd=0, yd=0, zd=0;
3442 //Float_t qinvtemp=GetQinv(0,track1, track2);
3445 if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
3446 else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
3448 for(Int_t i=0; i<fKbinsT-1; i++){
3449 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
3452 wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
3453 if(fMaxPt<=0.251) {fKtIndexL=0; fKtIndexH=0; wd=0;}
3454 if(fMinPt>0.249 && fKtIndexL==0) {fKtIndexL=1; wd=0;}
3456 if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
3457 else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
3459 for(Int_t i=0; i<kQbinsWeights-1; i++){
3460 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
3462 xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
3465 if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
3466 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
3468 for(Int_t i=0; i<kQbinsWeights-1; i++){
3469 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
3471 yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
3474 if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
3475 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
3477 for(Int_t i=0; i<kQbinsWeights-1; i++){
3478 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
3480 zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
3483 if(fLinearInterpolation){// Linear Interpolation of osl
3484 // w interpolation (kt)
3485 Float_t c000 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*wd;
3486 Float_t c100 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*wd;
3487 Float_t c010 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*wd;
3488 Float_t c001 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*wd;
3489 Float_t c110 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*wd;
3490 Float_t c101 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*wd;
3491 Float_t c011 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*wd;
3492 Float_t c111 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*wd;
3493 // x interpolation (qOut)
3494 Float_t c00 = c000*(1-xd) + c100*xd;
3495 Float_t c10 = c010*(1-xd) + c110*xd;
3496 Float_t c01 = c001*(1-xd) + c101*xd;
3497 Float_t c11 = c011*(1-xd) + c111*xd;
3498 // y interpolation (qSide)
3499 Float_t c0 = c00*(1-yd) + c10*yd;
3500 Float_t c1 = c01*(1-yd) + c11*yd;
3501 // z interpolation (qLong)
3502 wgt = (c0*(1-zd) + c1*zd);
3503 }else{// cubic interpolation of osl
3505 for(Int_t x=0; x<4; x++){
3506 for(Int_t y=0; y<4; y++){
3507 for(Int_t z=0; z<4; z++){
3508 Int_t binO = fQoIndexL + x;
3509 Int_t binS = fQsIndexL + y;
3510 Int_t binL = fQlIndexL + z;
3511 if(binO<=0) binO = 1;
3512 if(binS<=0) binS = 1;
3513 if(binL<=0) binL = 1;
3514 if(binO>kQbinsWeights) binO = kQbinsWeights;
3515 if(binS>kQbinsWeights) binS = kQbinsWeights;
3516 if(binL>kQbinsWeights) binL = kQbinsWeights;
3517 farrP1[x][y][z] = fNormWeight[fKtIndexL][fMbin]->GetBinContent(binO,binS,binL);
3518 farrP2[x][y][z] = fNormWeight[fKtIndexH][fMbin]->GetBinContent(binO,binS,binL);
3522 Float_t coord[3]={xd, yd, zd};
3523 Float_t c0 = nCubicInterpolate(3, (Float_t*) farrP1, coord);
3524 Float_t c1 = nCubicInterpolate(3, (Float_t*) farrP2, coord);
3526 wgt = c0*(1-wd) + c1*wd;
3530 // simplified stat error
3531 Float_t avgErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3532 avgErr += fNormWeight[fKtIndexH][fMbin]->GetBinError(fQoIndexL+1,fQsIndexL+1,fQlIndexL+1);
3539 //________________________________________________________________________
3540 Float_t AliFourPion::MCWeight(Int_t c[2], Float_t R, Float_t fcSq, Float_t qinv, Float_t k12){
3542 Float_t radius = R/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3543 Float_t r12=radius*(1-k12/2.0);
3545 Float_t coulCorr12 = FSICorrelation(c[0], c[1], qinv);
3547 Float_t arg=qinv*r12;
3548 Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3549 EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3550 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
3552 return ((1-fcSq) + fcSq*coulCorr12);
3556 //________________________________________________________________________
3557 Float_t AliFourPion::MCWeightOSL(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t qo, Float_t qs, Float_t ql){
3559 Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3560 Float_t radiusSide = radiusOut;
3561 Float_t radiusLong = radiusOut;
3562 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3563 Float_t coulCorr12 = FSICorrelation(charge1, charge2, qinv);
3564 if(charge1==charge2){
3565 return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
3567 return ((1-myDamp) + myDamp*coulCorr12);
3572 //________________________________________________________________________
3573 Float_t AliFourPion::MCWeight3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3], Float_t kT[3]){
3574 // FSI + QS correlations
3575 if(term==5) return 1.0;
3577 Float_t radius=R/0.19733;
3578 Float_t r12=radius*(1-kT[0]/2.0);
3579 Float_t r13=radius*(1-kT[1]/2.0);
3580 Float_t r23=radius*(1-kT[2]/2.0);
3582 Float_t fc = sqrt(fcSq);
3585 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3586 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3587 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3589 if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3590 Float_t arg12=qinv[0]*r12;
3591 Float_t arg13=qinv[1]*r13;
3592 Float_t arg23=qinv[2]*r23;
3593 Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
3594 EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
3595 Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
3596 EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
3597 Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
3598 EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
3600 Float_t C3QS = 1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2) + exp(-pow(qinv[1]*r13,2))*pow(EW13,2) + exp(-pow(qinv[2]*r23,2))*pow(EW23,2);
3601 C3QS += 2*exp(-(pow(r12,2)*pow(qinv[0],2) + pow(r13,2)*pow(qinv[1],2) + pow(r23,2)*pow(qinv[2],2))/2.)*EW12*EW13*EW23;
3602 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3603 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12;
3604 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13;
3605 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23;
3606 C3 += pow(fc,3)*C3QS*Kfactor12*Kfactor13*Kfactor23;
3609 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12);
3611 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13);
3613 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23);
3616 }else{// mixed charge case
3617 Float_t arg=qinv[0]*r12;
3618 Float_t KfactorSC = Kfactor12;
3619 Float_t KfactorMC1 = Kfactor13;
3620 Float_t KfactorMC2 = Kfactor23;
3621 if(c[0]==c[2]) {arg=qinv[1]*r13; KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3622 if(c[1]==c[2]) {arg=qinv[2]*r23; KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3623 Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3624 EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3626 Float_t C3QS = 1 + exp(-pow(arg,2))*pow(EW,2);
3627 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3628 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(arg,2))*pow(EW,2))*KfactorSC;
3629 C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3630 C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3631 C3 += pow(fc,3)*C3QS*KfactorSC*KfactorMC1*KfactorMC2;
3634 if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3635 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3637 return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3639 if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3640 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3645 //________________________________________________________________________
3646 Float_t AliFourPion::MCWeightFSI3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3]){
3647 // FSI only (no QS correlations)
3648 if(term==5) return 1.0;
3650 Float_t fc = sqrt(fcSq);
3652 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3653 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3654 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3656 if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3658 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3659 C3 += pow(fc,2)*(1-fc)*Kfactor12;
3660 C3 += pow(fc,2)*(1-fc)*Kfactor13;
3661 C3 += pow(fc,2)*(1-fc)*Kfactor23;
3662 C3 += pow(fc,3)*Kfactor12*Kfactor13*Kfactor23;
3665 return ((1-fcSq) + fcSq*Kfactor12);
3667 return ((1-fcSq) + fcSq*Kfactor13);
3669 return ((1-fcSq) + fcSq*Kfactor23);
3672 }else{// mixed charge case
3673 Float_t KfactorSC = Kfactor12;
3674 Float_t KfactorMC1 = Kfactor13;
3675 Float_t KfactorMC2 = Kfactor23;
3676 if(c[0]==c[2]) {KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3677 if(c[1]==c[2]) {KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3679 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3680 C3 += pow(fc,2)*(1-fc)*KfactorSC;
3681 C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3682 C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3683 C3 += pow(fc,3)*KfactorSC*KfactorMC1*KfactorMC2;
3686 if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*KfactorSC);
3687 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3689 return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3691 if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*KfactorSC);
3692 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3697 //________________________________________________________________________
3698 Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6], Float_t kT[6]){
3699 if(term==13) return 1.0;
3702 // ----, ---+, --++, -+++, ++++
3705 // Term 1: 1-2-3-4 (all same-event)
3706 // Term 2: 1-2-3 4 (particle 4 from different event)
3707 // Term 3: 1-2-4 3 (particle 3 from different event)
3708 // Term 4: 1-3-4 2 (particle 2 from different event)
3709 // Term 5: 2-3-4 1 (particle 1 from different event)
3710 // Term 6: 1-2 3 4 (particle 1 and 2 from same event)
3716 // Term 12: 1 2 3 4 (all from different events)
3718 Float_t radius = R/0.19733;
3720 r[0]=radius*(1-kT[0]/2.0);
3721 r[1]=radius*(1-kT[1]/2.0);
3722 r[2]=radius*(1-kT[2]/2.0);
3723 r[3]=radius*(1-kT[3]/2.0);
3724 r[4]=radius*(1-kT[4]/2.0);
3725 r[5]=radius*(1-kT[5]/2.0);
3727 Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3729 Float_t fc = sqrt(fcSq);
3731 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3732 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3733 Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3734 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3735 Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3736 Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3737 Float_t arg12=qinv[0]*r[0];
3738 Float_t arg13=qinv[1]*r[1];
3739 Float_t arg14=qinv[2]*r[2];
3740 Float_t arg23=qinv[3]*r[3];
3741 Float_t arg24=qinv[4]*r[4];
3742 Float_t arg34=qinv[5]*r[5];
3743 // Exchange Amplitudes
3744 Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
3745 Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
3746 Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
3747 Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
3748 Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
3749 Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));
3751 if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3754 Float_t C4QS = 1 + pow(EA12,2) + pow(EA13,2) + pow(EA14,2) + pow(EA23,2) + pow(EA24,2) + pow(EA34,2);// baseline + single pairs
3755 C4QS += pow(EA12,2) * pow(EA34,2);// 2-pairs
3756 C4QS += pow(EA13,2) * pow(EA24,2);// 2-pairs
3757 C4QS += pow(EA14,2) * pow(EA23,2);// 2-pairs
3758 C4QS += 2*EA12*EA13*EA23 + 2*EA12*EA14*EA24 + 2*EA13*EA14*EA34 + 2*EA23*EA24*EA34;// 3-particle exhange
3759 C4QS += 3*EA12*EA23*EA34*EA14 + 3*EA12*EA13*EA34*EA24;// 4-particle exchange
3760 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3761 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA12,2))*Kfactor12 + (1 + pow(EA13,2))*Kfactor13 + (1 + pow(EA14,2))*Kfactor14 );
3762 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA23,2))*Kfactor23 + (1 + pow(EA24,2))*Kfactor24 + (1 + pow(EA34,2))*Kfactor34);
3763 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA13,2) + pow(EA23,2) + 2*EA12*EA13*EA23) * Kfactor12*Kfactor13*Kfactor23;
3764 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA14,2) + pow(EA24,2) + 2*EA12*EA14*EA24) * Kfactor12*Kfactor14*Kfactor24;
3765 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA13,2) + pow(EA14,2) + pow(EA34,2) + 2*EA13*EA14*EA34) * Kfactor13*Kfactor14*Kfactor34;
3766 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA23,2) + pow(EA24,2) + pow(EA34,2) + 2*EA23*EA24*EA34) * Kfactor23*Kfactor24*Kfactor34;
3767 C4 += pow(fc,4)*C4QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3770 Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0;
3771 if(term==2) {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3772 else if(term==3) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3773 else if(term==4) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3774 else {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3775 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2);
3776 C3QS += 2*EA1*EA2*EA3;
3777 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3778 C3 += pow(fc,2)*(1-fc)*( (1+pow(EA1,2))*Kpair1 + (1+pow(EA2,2))*Kpair2 + (1+pow(EA3,2))*Kpair3 );
3779 C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3782 if(term==6) return ((1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12);
3783 else if(term==7) return ((1-fcSq) + fcSq*(1 + pow(EA13,2))*Kfactor13);
3784 else if(term==8) return ((1-fcSq) + fcSq*(1 + pow(EA14,2))*Kfactor14);
3785 else if(term==9) return ((1-fcSq) + fcSq*(1 + pow(EA23,2))*Kfactor23);
3786 else if(term==10) return ((1-fcSq) + fcSq*(1 + pow(EA24,2))*Kfactor24);
3787 else return ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34);
3789 Float_t C22 = (1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12;
3790 C22 *= (1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34;
3794 }else{// mixed charge case
3795 if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3796 Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3797 Int_t c_OddOneOut = 1;
3798 if(ChargeSum==3) c_OddOneOut = 0;
3800 if(c[0]==c_OddOneOut) {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3801 else if(c[1]==c_OddOneOut) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3802 else if(c[2]==c_OddOneOut) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24; Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3803 else {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23; Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3806 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3807 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3808 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 + (1 + pow(EA3,2))*Kpair3 );
3809 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3810 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3) * Kpair1*Kpair2*Kpair3;
3811 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair4*Kpair5 + (1+pow(EA2,2))*Kpair2*Kpair4*Kpair6 + (1+pow(EA3,2))*Kpair3*Kpair5*Kpair6);// doesn't matter which MC K's
3812 C4 += pow(fc,4)*C3QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3814 }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3815 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3816 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3817 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;
3818 C3 += pow(fc,2)*(1-fc)*(1+pow(EA2,2))*Kpair2;
3819 C3 += pow(fc,2)*(1-fc)*(1+pow(EA3,2))*Kpair3;
3820 C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3822 }else if(term<=5){// one SC pair, two MC pairs
3823 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3824 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3825 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3826 C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3827 C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair5;
3829 }else if(term==6 || term==7){
3830 if(ChargeSum==1) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3831 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3833 return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3835 return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3836 }else if(term==10 || term==11){
3837 if(ChargeSum==3) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3838 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3839 }else return 1.0;// for 12 and 13
3840 }else{// --++ configuration
3841 Float_t EA1=0, EA2=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3842 if(c[0]==c[1]) {EA1=EA12; EA2=EA34; Kpair1=Kfactor12; Kpair2=Kfactor34; Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3843 else if(c[0]==c[2]) {EA1=EA13; EA2=EA24; Kpair1=Kfactor13; Kpair2=Kfactor24; Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3844 else {EA1=EA14; EA2=EA23; Kpair1=Kfactor14; Kpair2=Kfactor23; Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3847 Float_t C2QS = 1 + pow(EA1,2)*pow(EA2,2);
3848 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3849 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 );
3850 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair3 + Kpair4 + Kpair5 + Kpair6 );
3851 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair3*Kpair4 + (1 + pow(EA2,2))*Kpair2*Kpair3*Kpair4);
3852 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair5*Kpair6 + (1 + pow(EA2,2))*Kpair2*Kpair5*Kpair6);// doesn't matter which two MC K's used
3853 C4 += pow(fc,4)*C2QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3856 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3857 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3858 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3859 C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3860 C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair6;
3862 }else if(term==6 || term==11){
3863 return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3864 }else if(term!=12 && term !=13){
3865 return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3867 Float_t C22 = (1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1;
3868 C22 *= (1-fcSq) + fcSq*(1 + pow(EA2,2))*Kpair2;
3875 //________________________________________________________________________
3876 Float_t AliFourPion::MCWeightFSI4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6]){
3877 if(term==13) return 1.0;
3879 Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3881 Float_t fc = sqrt(fcSq);
3883 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3884 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3885 Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3886 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3887 Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3888 Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3890 if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3893 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3894 C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor12 + Kfactor13 + Kfactor14 );
3895 C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor23 + Kfactor24 + Kfactor34 );
3896 C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor13*Kfactor23;
3897 C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor14*Kfactor24;
3898 C4 += pow(fc,3)*(1-fc) * Kfactor13*Kfactor14*Kfactor34;
3899 C4 += pow(fc,3)*(1-fc) * Kfactor23*Kfactor24*Kfactor34;
3900 C4 += pow(fc,4) * Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3903 Float_t Kpair1=0, Kpair2=0, Kpair3=0;
3904 if(term==2) {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3905 else if(term==3) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3906 else if(term==4) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3907 else {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3908 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3909 C3 += pow(fc,2)*(1-fc)*( Kpair1 + Kpair2 + Kpair3 );
3910 C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3913 if(term==6) return ((1-fcSq) + fcSq*Kfactor12);
3914 else if(term==7) return ((1-fcSq) + fcSq*Kfactor13);
3915 else if(term==8) return ((1-fcSq) + fcSq*Kfactor14);
3916 else if(term==9) return ((1-fcSq) + fcSq*Kfactor23);
3917 else if(term==10) return ((1-fcSq) + fcSq*Kfactor24);
3918 else return ((1-fcSq) + fcSq*Kfactor34);
3920 Float_t C22 = (1-fcSq) + fcSq*Kfactor12;
3921 C22 *= (1-fcSq) + fcSq*Kfactor34;
3925 }else{// mixed charge case
3926 if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3927 Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3928 Int_t c_OddOneOut = 1;
3929 if(ChargeSum==3) c_OddOneOut = 0;
3931 if(c[0]==c_OddOneOut) {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3932 else if(c[1]==c_OddOneOut) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3933 else if(c[2]==c_OddOneOut) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24; Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3934 else {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23; Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3937 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3938 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 );
3939 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3940 C4 += pow(fc,3)*(1-fc)*Kpair1*Kpair2*Kpair3;
3941 C4 += pow(fc,3)*(1-fc)*(Kpair1*Kpair4*Kpair5 + Kpair2*Kpair4*Kpair6 + Kpair3*Kpair5*Kpair6);// doesn't matter which two MC K's used
3942 C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3944 }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3945 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3946 C3 += pow(fc,2)*(1-fc)*Kpair1;
3947 C3 += pow(fc,2)*(1-fc)*Kpair2;
3948 C3 += pow(fc,2)*(1-fc)*Kpair3;
3949 C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3951 }else if(term<=5){// one SC pair, two MC pairs
3952 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3953 C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3954 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3955 C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3956 C3 += pow(fc,3)*Kpair1*Kpair4*Kpair5;
3958 }else if(term==6 || term==7){
3959 if(ChargeSum==1) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3960 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3962 return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3964 return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3965 }else if(term==10 || term==11){
3966 if(ChargeSum==3) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3967 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3968 }else return 1.0;// 12 and 13
3969 }else{// --++ configuration
3970 Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3971 if(c[0]==c[1]) {Kpair1=Kfactor12; Kpair2=Kfactor34; Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3972 else if(c[0]==c[2]) {Kpair1=Kfactor13; Kpair2=Kfactor24; Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3973 else {Kpair1=Kfactor14; Kpair2=Kfactor23; Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3976 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3977 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 + Kpair4 + Kpair5 + Kpair6);
3978 C4 += pow(fc,3)*(1-fc)*( Kpair1*Kpair3*Kpair4 + Kpair2*Kpair3*Kpair4 + Kpair1*Kpair5*Kpair6 + Kpair2*Kpair5*Kpair6);
3979 C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3982 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3983 C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3984 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3985 C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3986 C3 += pow(fc,3)*Kpair1*Kpair4*Kpair6;
3988 }else if(term==6 || term==11){
3989 return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3990 }else if(term !=12 && term !=13){
3991 return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3993 Float_t C22 = (1-fcSq) + fcSq*Kpair1;
3994 C22 *= (1-fcSq) + fcSq*Kpair2;
4001 //________________________________________________________________________
4002 void AliFourPion::SetMomResCorrections(Bool_t legoCase, TH2D *temp2DSC, TH2D *temp2DMC){
4006 cout<<"LEGO call to SetMomResCorrections"<<endl;
4007 fMomResC2SC = (TH2D*)temp2DSC->Clone();
4008 fMomResC2SC->SetDirectory(0);
4009 fMomResC2MC = (TH2D*)temp2DMC->Clone();
4010 fMomResC2MC->SetDirectory(0);
4012 TFile *momResFile = new TFile("MomResFile.root","READ");
4013 if(!momResFile->IsOpen()) {
4014 cout<<"No momentum resolution file found"<<endl;
4015 AliFatal("No momentum resolution file found. Kill process.");
4016 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
4018 TH2D *temp2DSC2 = (TH2D*)momResFile->Get("MRC_C2_SC");
4019 fMomResC2SC = (TH2D*)temp2DSC2->Clone();
4020 fMomResC2SC->SetDirectory(0);
4022 TH2D *temp2DMC2 = (TH2D*)momResFile->Get("MRC_C2_MC");
4023 fMomResC2MC = (TH2D*)temp2DMC2->Clone();
4024 fMomResC2MC->SetDirectory(0);
4026 momResFile->Close();
4030 for(Int_t bx=1; bx<=fMomResC2SC->GetNbinsX(); bx++){
4031 for(Int_t by=1; by<=fMomResC2SC->GetNbinsY(); by++){
4032 if(fMomResC2SC->GetBinContent(bx,by) > 1.5) fMomResC2SC->SetBinContent(bx,by, 1.0);// Maximum is ~1.02
4033 if(fMomResC2SC->GetBinContent(bx,by) < 0.8) fMomResC2SC->SetBinContent(bx,by, 1.0);// Minimum is ~0.8
4034 if(fMomResC2MC->GetBinContent(bx,by) > 1.5) fMomResC2MC->SetBinContent(bx,by, 1.0);// Maximum is ~1.02
4035 if(fMomResC2MC->GetBinContent(bx,by) < 0.8) fMomResC2MC->SetBinContent(bx,by, 1.0);// Minimum is ~0.8
4039 cout<<"Done reading momentum resolution file"<<endl;
4041 //________________________________________________________________________
4042 void AliFourPion::SetFSICorrelations(Bool_t legoCase, TH1D *tempss[12], TH1D *tempos[12]){
4043 // read in 2-particle and 3-particle FSI correlations = K2 & K3
4044 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
4046 cout<<"LEGO call to SetFSICorrelations"<<endl;
4047 for(Int_t MB=0; MB<12; MB++) {
4048 fFSIss[MB] = (TH1D*)tempss[MB]->Clone();
4049 fFSIos[MB] = (TH1D*)tempos[MB]->Clone();
4051 fFSIss[MB]->SetDirectory(0);
4052 fFSIos[MB]->SetDirectory(0);
4055 cout<<"non LEGO call to SetFSICorrelations"<<endl;
4056 TFile *fsifile = new TFile("KFile.root","READ");
4057 if(!fsifile->IsOpen()) {
4058 cout<<"No FSI file found"<<endl;
4059 AliFatal("No FSI file found. Kill process.");
4060 }else {cout<<"Good FSI File Found!"<<endl;}
4062 TH1D *temphistoSS[12];
4063 TH1D *temphistoOS[12];
4064 for(Int_t MB=0; MB<12; MB++) {
4065 TString *nameK2SS = new TString("K2ss_");
4067 temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
4069 TString *nameK2OS = new TString("K2os_");
4071 temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
4073 fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
4074 fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
4075 fFSIss[MB]->SetDirectory(0);
4076 fFSIos[MB]->SetDirectory(0);
4083 cout<<"Done reading FSI file"<<endl;
4085 //________________________________________________________________________
4086 Float_t AliFourPion::FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
4087 // returns 2-particle Coulomb correlations = K2
4088 Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
4089 Int_t qbinH = qbinL+1;
4090 if(qbinL <= 0) return 1.0;
4091 if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) {
4092 if(charge1!=charge2) {
4093 Float_t ScaleFac = (fFSIos[fFSIindex]->GetBinContent(fFSIos[fFSIindex]->GetNbinsX()-1) - 1);
4094 ScaleFac /= (Gamov(charge1, charge2, fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(fFSIos[fFSIindex]->GetNbinsX()-1)) - 1);
4095 return ( (Gamov(charge1, charge2, qinv)-1)*ScaleFac + 1);
4097 Float_t ScaleFac = (fFSIss[fFSIindex]->GetBinContent(fFSIss[fFSIindex]->GetNbinsX()-1) - 1);
4098 ScaleFac /= (Gamov(charge1, charge2, fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(fFSIss[fFSIindex]->GetNbinsX()-1)) - 1);
4099 return ( (Gamov(charge1, charge2, qinv)-1)*ScaleFac + 1);
4104 if(charge1==charge2){
4105 slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
4106 slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
4107 return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
4109 slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
4110 slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
4111 return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
4114 //________________________________________________________________________
4115 void AliFourPion::SetFillBins2(Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
4116 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
4117 else {b1=c1; b2=c2;}
4119 //________________________________________________________________________
4120 void AliFourPion::SetFillBins3(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){
4122 // "part" specifies which pair is from the same event. Only relevant for terms 2-4
4124 if(part==1) {// default case (irrelevant for term 1 and term 5)
4125 if(c1==c2) seSS=kTRUE;
4128 if(c1==c3) seSS=kTRUE;
4132 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
4133 if( (c1+c2+c3)==1) {
4134 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
4136 if(seSS) fill2=kTRUE;
4137 else {fill3=kTRUE; fill4=kTRUE;}
4139 }else if( (c1+c2+c3)==2) {
4142 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
4146 b1=c1; b2=c2; b3=c3;
4147 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
4151 //________________________________________________________________________
4152 void AliFourPion::SetFillBins4(Int_t c1, Int_t c2, Int_t c3, Int_t c4, Int_t &b1, Int_t &b2, Int_t &b3, Int_t &b4, Int_t ENsum, Bool_t fillTerm[13]){
4154 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
4155 if( (c1+c2+c3+c4)==0 || (c1+c2+c3+c4)==4) {// all of the same charge: ---- or ++++
4157 b1=c1; b2=c2; b3=c3; b4=c4;
4158 if(ENsum==0) fillTerm[0]=kTRUE;
4159 else if(ENsum==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
4160 else if(ENsum==2) {fillTerm[11]=kTRUE;}
4161 else if(ENsum==3) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
4162 else fillTerm[12]=kTRUE;
4164 }else if( (c1+c2+c3+c4)==1) {// one positive charge: ---+
4166 b1=0; b2=0; b3=0; b4=1;// Re-assign to merge degenerate histos
4167 if(ENsum==0) fillTerm[0]=kTRUE;
4169 if(c4==1) fillTerm[1]=kTRUE;
4170 else {fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
4174 if(c3==1 || c4==1) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[8]=kTRUE;}
4175 else {fillTerm[7]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
4176 }else fillTerm[12]=kTRUE;
4178 }else if( (c1+c2+c3+c4)==2) {// two positive charges: --++
4180 b1=0; b2=0; b3=1; b4=1;// Re-assign to merge degenerate histos
4181 if(ENsum==0) fillTerm[0]=kTRUE;
4183 if(c4==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE;}
4184 else {fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
4186 if( (c1+c2)==0) fillTerm[11]=kTRUE;
4188 if( (c1+c2)==0) fillTerm[5]=kTRUE;
4189 else if( (c1+c2)==1) {fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE;}
4190 else fillTerm[10]=kTRUE;
4191 }else fillTerm[12]=kTRUE;
4193 }else{// three positive charges
4195 b1=0; b2=1; b3=1; b4=1;// Re-assign to merge degenerate histos
4196 if(ENsum==0) fillTerm[0]=kTRUE;
4198 if(c4==0) fillTerm[4]=kTRUE;
4199 else {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE;}
4203 if(c3==0 || c4==0) {fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
4204 else {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE;}
4205 }else fillTerm[12]=kTRUE;
4210 //________________________________________________________________________
4211 void AliFourPion::SetFSIindex(Float_t R){
4214 if(fMbin==0) fFSIindex = 0;//0-5%
4215 else if(fMbin==1) fFSIindex = 1;//5-10%
4216 else if(fMbin<=3) fFSIindex = 2;//10-20%
4217 else if(fMbin<=5) fFSIindex = 3;//20-30%
4218 else if(fMbin<=7) fFSIindex = 4;//30-40%
4219 else if(fMbin<=9) fFSIindex = 5;//40-50%
4220 else if(fMbin<=12) fFSIindex = 6;//40-50%
4221 else if(fMbin<=15) fFSIindex = 7;//40-50%
4222 else if(fMbin<=18) fFSIindex = 8;//40-50%
4223 else fFSIindex = 8;//90-100%
4224 }else fFSIindex = 9;// pp and pPb
4225 }else{// FSI binning for MC
4226 if(R>=10.) fFSIindex = 0;
4227 else if(R>=9.) fFSIindex = 1;
4228 else if(R>=8.) fFSIindex = 2;
4229 else if(R>=7.) fFSIindex = 3;
4230 else if(R>=6.) fFSIindex = 4;
4231 else if(R>=5.) fFSIindex = 5;
4232 else if(R>=4.) fFSIindex = 6;
4233 else if(R>=3.) fFSIindex = 7;
4234 else if(R>=2.) fFSIindex = 8;
4238 //________________________________________________________________________
4239 void AliFourPion::SetMuonCorrections(Bool_t legoCase, TH2D *tempMuon){
4241 cout<<"LEGO call to SetMuonCorrections"<<endl;
4242 fWeightmuonCorrection = (TH2D*)tempMuon->Clone();
4243 fWeightmuonCorrection->SetDirectory(0);
4245 cout<<"non LEGO call to SetMuonCorrections"<<endl;
4246 TFile *MuonFile=new TFile("MuonCorrection.root","READ");
4247 if(!MuonFile->IsOpen()) {
4248 cout<<"No Muon file found"<<endl;
4249 AliFatal("No Muon file found. Kill process.");
4250 }else {cout<<"Good Muon File Found!"<<endl;}
4252 fWeightmuonCorrection = (TH2D*)MuonFile->Get("WeightmuonCorrection");
4253 fWeightmuonCorrection->SetDirectory(0);
4257 cout<<"Done reading Muon file"<<endl;
4259 //________________________________________________________________________
4260 Float_t AliFourPion::cubicInterpolate (Float_t p[4], Float_t x) {
4261 return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));// Paulinternet
4263 //________________________________________________________________________
4264 Float_t AliFourPion::nCubicInterpolate (Int_t n, Float_t* p, Float_t coordinates[]) {
4267 return cubicInterpolate(p, *coordinates);
4271 Int_t skip = 1 << (n - 1) * 2;
4272 arr[0] = nCubicInterpolate(n - 1, p, coordinates + 1);
4273 arr[1] = nCubicInterpolate(n - 1, p + skip, coordinates + 1);
4274 arr[2] = nCubicInterpolate(n - 1, p + 2*skip, coordinates + 1);
4275 arr[3] = nCubicInterpolate(n - 1, p + 3*skip, coordinates + 1);
4276 return cubicInterpolate(arr, *coordinates);