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"
32 #include "AliChaoticity.h"
35 #define G_Coeff 0.006399 // 2*pi*alpha*M_pion
36 #define kappa3 0.24 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
37 #define kappa4 0.16 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
40 // Author: Dhevan Gangadharan
42 ClassImp(AliChaoticity)
44 //________________________________________________________________________
45 AliChaoticity::AliChaoticity():
59 fGenerateSignal(kFALSE),
60 fPdensityExplicitLoop(kFALSE),
61 fPdensityPairCut(kTRUE),
62 fTabulatePairs(kFALSE),
64 fFixedLambdaBinMomRes(9),
65 fFixedLambdaBinr3(10),
87 fQupperBoundWeights(0),
103 fMinSepPairEta(0.03),
104 fMinSepPairPhi(0.04),
129 fOtherPairLocation1(),
130 fOtherPairLocation2(),
138 // Default constructor
139 for(Int_t mb=0; mb<fMbins; mb++){
140 for(Int_t edB=0; edB<fEDbins; edB++){
141 for(Int_t c1=0; c1<2; c1++){
142 for(Int_t c2=0; c2<2; c2++){
143 for(Int_t sc=0; sc<kSCLimit2; sc++){
144 for(Int_t term=0; term<2; term++){
146 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
148 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
149 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
150 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
151 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
152 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
153 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
158 for(Int_t c3=0; c3<2; c3++){
159 for(Int_t sc=0; sc<kSCLimit3; sc++){
160 for(Int_t term=0; term<5; term++){
162 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
163 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
164 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
165 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
166 for(Int_t dt=0; dt<kDENtypes; dt++){
167 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
168 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = 0x0;
169 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = 0x0;
170 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = 0x0;
171 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = 0x0;
172 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = 0x0;
173 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = 0x0;
181 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
182 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
183 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
184 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
191 // Initialize 3-pion FSI histograms
192 for(Int_t i=0; i<6; i++){
198 // Initialize fNormWeight and fNormWeightErr to 0
199 for(Int_t i=0; i<3; i++){// Kt iterator
200 for(Int_t j=0; j<10; j++){// Mbin iterator
201 fNormWeight[i][j]=0x0;
207 //________________________________________________________________________
208 AliChaoticity::AliChaoticity(const Char_t *name)
209 : AliAnalysisTaskSE(name),
222 fGenerateSignal(kFALSE),
223 fPdensityExplicitLoop(kFALSE),
224 fPdensityPairCut(kTRUE),
225 fTabulatePairs(kFALSE),
227 fFixedLambdaBinMomRes(9),
228 fFixedLambdaBinr3(10),
239 fCentBinHighLimit(1),
250 fQupperBoundWeights(0),
266 fMinSepPairEta(0.03),
267 fMinSepPairPhi(0.04),
292 fOtherPairLocation1(),
293 fOtherPairLocation2(),
303 fPdensityExplicitLoop = kFALSE;
304 fPdensityPairCut = kTRUE;
307 for(Int_t mb=0; mb<fMbins; mb++){
308 for(Int_t edB=0; edB<fEDbins; edB++){
309 for(Int_t c1=0; c1<2; c1++){
310 for(Int_t c2=0; c2<2; c2++){
311 for(Int_t sc=0; sc<kSCLimit2; sc++){
312 for(Int_t term=0; term<2; term++){
314 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
316 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
317 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
318 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
319 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
320 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
321 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
326 for(Int_t c3=0; c3<2; c3++){
327 for(Int_t sc=0; sc<kSCLimit3; sc++){
328 for(Int_t term=0; term<5; term++){
330 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
331 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
332 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
333 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
334 for(Int_t dt=0; dt<kDENtypes; dt++){
335 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
336 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = 0x0;
337 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = 0x0;
338 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = 0x0;
339 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = 0x0;
340 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = 0x0;
341 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = 0x0;
349 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
350 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
351 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
352 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
359 // Initialize 3-pion FSI histograms
360 for(Int_t i=0; i<6; i++){
365 // Initialize fNormWeight and fNormWeightErr to 0
366 for(Int_t i=0; i<3; i++){// Kt iterator
367 for(Int_t j=0; j<10; j++){// Mbin iterator
368 fNormWeight[i][j]=0x0;
373 DefineOutput(1, TList::Class());
375 //________________________________________________________________________
376 AliChaoticity::AliChaoticity(const AliChaoticity &obj)
377 : AliAnalysisTaskSE(obj.fname),
381 fOutputList(obj.fOutputList),
382 fPIDResponse(obj.fPIDResponse),
385 fTempStruct(obj.fTempStruct),
386 fRandomNumber(obj.fRandomNumber),
388 fMCcase(obj.fMCcase),
389 fAODcase(obj.fAODcase),
390 fPbPbcase(obj.fPbPbcase),
391 fGenerateSignal(obj.fGenerateSignal),
392 fPdensityExplicitLoop(obj.fPdensityExplicitLoop),
393 fPdensityPairCut(obj.fPdensityPairCut),
394 fTabulatePairs(obj.fTabulatePairs),
396 fFixedLambdaBinMomRes(obj.fFixedLambdaBinMomRes),
397 fFixedLambdaBinr3(obj.fFixedLambdaBinr3),
398 fFilterBit(obj.fFilterBit),
399 fMaxChi2NDF(obj.fMaxChi2NDF),
400 fMinTPCncls(obj.fMinTPCncls),
401 fBfield(obj.fBfield),
403 fFSIbin(obj.fFSIbin),
406 fMultLimit(obj.fMultLimit),
407 fCentBinLowLimit(obj.fCentBinLowLimit),
408 fCentBinHighLimit(obj.fCentBinHighLimit),
409 fEventCounter(obj.fEventCounter),
410 fEventsToMix(obj.fEventsToMix),
411 fZvertexBins(obj.fZvertexBins),
414 fQLowerCut(obj.fQLowerCut),
417 fKupperBound(obj.fKupperBound),
418 fQupperBound(obj.fQupperBound),
419 fQupperBoundWeights(obj.fQupperBoundWeights),
427 fQstepWeights(obj.fQstepWeights),
429 fDampStart(obj.fDampStart),
430 fDampStep(obj.fDampStep),
431 fTPCTOFboundry(obj.fTPCTOFboundry),
432 fTOFboundry(obj.fTOFboundry),
433 fSigmaCutTPC(obj.fSigmaCutTPC),
434 fSigmaCutTOF(obj.fSigmaCutTOF),
435 fMinSepPairEta(obj.fMinSepPairEta),
436 fMinSepPairPhi(obj.fMinSepPairPhi),
437 fShareQuality(obj.fShareQuality),
438 fShareFraction(obj.fShareFraction),
439 fTrueMassP(obj.fTrueMassP),
440 fTrueMassPi(obj.fTrueMassPi),
441 fTrueMassK(obj.fTrueMassK),
442 fTrueMassKs(obj.fTrueMassKs),
443 fTrueMassLam(obj.fTrueMassLam),
444 fKtIndexL(obj.fKtIndexL),
445 fKtIndexH(obj.fKtIndexH),
446 fQoIndexL(obj.fQoIndexL),
447 fQoIndexH(obj.fQoIndexH),
448 fQsIndexL(obj.fQsIndexL),
449 fQsIndexH(obj.fQsIndexH),
450 fQlIndexL(obj.fQlIndexL),
451 fQlIndexH(obj.fQlIndexH),
452 fDummyB(obj.fDummyB),
461 fOtherPairLocation1(),
462 fOtherPairLocation2(),
466 fMomResC2(obj.fMomResC2),
467 fFSI2SS(obj.fFSI2SS),
472 for(Int_t i=0; i<6; i++){
473 fFSIOmega0SS[i]=obj.fFSIOmega0SS[i];
474 fFSIOmega0OS[i]=obj.fFSIOmega0OS[i];
477 // Initialize fNormWeight and fNormWeightErr to 0
478 for(Int_t i=0; i<3; i++){// Kt iterator
479 for(Int_t j=0; j<10; j++){// Mbin iterator
480 fNormWeight[i][j]=0x0;
486 //________________________________________________________________________
487 AliChaoticity &AliChaoticity::operator=(const AliChaoticity &obj)
489 // Assignment operator
495 fOutputList = obj.fOutputList;
496 fPIDResponse = obj.fPIDResponse;
499 fTempStruct = obj.fTempStruct;
500 fRandomNumber = obj.fRandomNumber;
502 fMCcase = obj.fMCcase;
503 fAODcase = obj.fAODcase;
504 fPbPbcase = obj.fPbPbcase;
505 fGenerateSignal = obj.fGenerateSignal;
506 fPdensityExplicitLoop = obj.fPdensityExplicitLoop;
507 fPdensityPairCut = obj.fPdensityPairCut;
508 fTabulatePairs = obj.fTabulatePairs;
510 fFixedLambdaBinMomRes = obj.fFixedLambdaBinMomRes;
511 fFixedLambdaBinr3 = obj.fFixedLambdaBinr3;
512 fFilterBit = obj.fFilterBit;
513 fMaxChi2NDF = obj.fMaxChi2NDF;
514 fMinTPCncls = obj.fMinTPCncls;
515 fBfield = obj.fBfield;
517 fFSIbin = obj.fFSIbin;
520 fMultLimit = obj.fMultLimit;
521 fCentBinLowLimit = obj.fCentBinLowLimit;
522 fCentBinHighLimit = obj.fCentBinHighLimit;
523 fEventCounter = obj.fEventCounter;
524 fEventsToMix = obj.fEventsToMix;
525 fZvertexBins = obj.fZvertexBins;
526 fQLowerCut = obj.fQLowerCut;
527 fKupperBound = obj.fKupperBound;
528 fQupperBound = obj.fQupperBound;
529 fQupperBoundWeights = obj.fQupperBoundWeights;
531 fQstepWeights = obj.fQstepWeights;
532 fDampStart = obj.fDampStart;
533 fDampStep = obj.fDampStep;
534 fTPCTOFboundry = obj.fTPCTOFboundry;
535 fTOFboundry = obj.fTOFboundry;
536 fSigmaCutTPC = obj.fSigmaCutTPC;
537 fSigmaCutTOF = obj.fSigmaCutTOF;
538 fMinSepPairEta = obj.fMinSepPairEta;
539 fMinSepPairPhi = obj.fMinSepPairPhi;
540 fShareQuality = obj.fShareQuality;
541 fShareFraction = obj.fShareFraction;
542 fTrueMassP = obj.fTrueMassP;
543 fTrueMassPi = obj.fTrueMassPi;
544 fTrueMassK = obj.fTrueMassK;
545 fTrueMassKs = obj.fTrueMassKs;
546 fTrueMassLam = obj.fTrueMassLam;
547 fKtIndexL = obj.fKtIndexL;
548 fKtIndexH = obj.fKtIndexH;
549 fQoIndexL = obj.fQoIndexL;
550 fQoIndexH = obj.fQoIndexH;
551 fQsIndexL = obj.fQsIndexL;
552 fQsIndexH = obj.fQsIndexH;
553 fQlIndexL = obj.fQlIndexL;
554 fQlIndexH = obj.fQlIndexH;
555 fDummyB = obj.fDummyB;
556 fMomResC2 = obj.fMomResC2;
557 fFSI2SS = obj.fFSI2SS;
558 fFSI2OS = obj.fFSI2OS;
560 for(Int_t i=0; i<6; i++){
561 fFSIOmega0SS[i]=obj.fFSIOmega0SS[i];
562 fFSIOmega0OS[i]=obj.fFSIOmega0OS[i];
564 for(Int_t i=0; i<3; i++){// Kt iterator
565 for(Int_t j=0; j<10; j++){// Mbin iterator
566 fNormWeight[i][j]=obj.fNormWeight[i][j];
572 //________________________________________________________________________
573 AliChaoticity::~AliChaoticity()
576 if(fAOD) delete fAOD;
577 //if(fESD) delete fESD;
578 if(fOutputList) delete fOutputList;
579 if(fPIDResponse) delete fPIDResponse;
581 if(fEvt) delete fEvt;
582 if(fTempStruct) delete [] fTempStruct;
583 if(fRandomNumber) delete fRandomNumber;
584 if(fMomResC2) delete fMomResC2;
585 if(fFSI2SS) delete fFSI2SS;
586 if(fFSI2OS) delete fFSI2OS;
588 for(Int_t i=0; i<fMultLimit; i++){
589 if(fPairLocationSE[i]) delete [] fPairLocationSE[i];
590 if(fPairLocationME[i]) delete [] fPairLocationME[i];
591 for(Int_t j=0; j<2; j++){
592 if(fOtherPairLocation1[j][i]) delete [] fOtherPairLocation1[j][i];
593 if(fOtherPairLocation2[j][i]) delete [] fOtherPairLocation2[j][i];
595 for(Int_t j=0; j<3; j++) if(fNormPairSwitch[j][i]) delete [] fNormPairSwitch[j][i];
596 for(Int_t j=0; j<4; j++) if(fPairSplitCut[j][i]) delete [] fPairSplitCut[j][i];
598 for(Int_t i=0; i<kPairLimit; i++) if(fTripletSkip1[i]) delete [] fTripletSkip1[i];
599 for(Int_t i=0; i<2*kPairLimit; i++) if(fTripletSkip2[i]) delete [] fTripletSkip2[i];
600 for(Int_t i=0; i<3; i++) if(fNormPairs[i]) delete [] fNormPairs[i];
602 for(Int_t mb=0; mb<fMbins; mb++){
603 for(Int_t edB=0; edB<fEDbins; edB++){
604 for(Int_t c1=0; c1<2; c1++){
605 for(Int_t c2=0; c2<2; c2++){
606 for(Int_t sc=0; sc<kSCLimit2; sc++){
607 for(Int_t term=0; term<2; term++){
609 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2;
611 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal;
612 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared;
613 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL;
614 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW;
615 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL;
616 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW;
618 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
619 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
623 for(Int_t c3=0; c3<2; c3++){
624 for(Int_t sc=0; sc<kSCLimit3; sc++){
625 for(Int_t term=0; term<5; term++){
627 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3;
628 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3;
629 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3;
630 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3;
631 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms;
632 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms;
633 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal;
634 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal;
635 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared;
636 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared;
637 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Q3W) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Q3W;
638 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Q3W) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Q3W;
640 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3;
641 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3;
642 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3;
643 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3;
645 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2;
646 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2;
647 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2;
648 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2;
651 for(Int_t dt=0; dt<kDENtypes; dt++){
652 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm;
653 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm;
654 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm;
655 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal;
656 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal;
657 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared;
658 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared;
667 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
668 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
669 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD;
670 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD;
678 for(Int_t i=0; i<6; i++){
679 if(fFSIOmega0SS[i]) delete fFSIOmega0SS[i];
680 if(fFSIOmega0OS[i]) delete fFSIOmega0OS[i];
682 for(Int_t i=0; i<3; i++){// Kt iterator
683 for(Int_t j=0; j<10; j++){// Mbin iterator
684 if(fNormWeight[i][j]) delete fNormWeight[i][j];
689 //________________________________________________________________________
690 void AliChaoticity::ParInit()
692 cout<<"AliChaoticity MyInit() call"<<endl;
693 cout<<"lego:"<<fLEGO<<" MCcase:"<<fMCcase<<" PbPbcase:"<<fPbPbcase<<" TabulatePairs:"<<fTabulatePairs<<" GenSignal:"<<fGenerateSignal<<" CentLow:"<<fCentBinLowLimit<<" CentHigh:"<<fCentBinHighLimit<<" RMax:"<<fRMax<<" LambdaBinMomRes:"<<fFixedLambdaBinMomRes<<" LambdaBinr3:"<<fFixedLambdaBinr3<<" FB:"<<fFilterBit<<" MaxChi2/NDF:"<<fMaxChi2NDF<<" MinTPCncls:"<<fMinTPCncls<<" MinPairSepEta:"<<fMinSepPairEta<<" MinPairSepPhi:"<<fMinSepPairPhi<<" NsigTPC:"<<fSigmaCutTPC<<" NsigTOF:"<<fSigmaCutTOF<<endl;
695 fRandomNumber = new TRandom3();
696 fRandomNumber->SetSeed(0);
700 if(fPdensityExplicitLoop) fEventsToMix=3;
701 else if(fPdensityPairCut && !fPdensityExplicitLoop) fEventsToMix=2;
705 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
706 fTOFboundry = 2.1;// TOF pid used below this momentum
708 ////////////////////////////////////////////////
710 fShareQuality = .5;// max
711 fShareFraction = .05;// max
712 ////////////////////////////////////////////////
715 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
716 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
720 if(fPbPbcase) {// PbPb
721 fMultLimit=kMultLimitPbPb;
723 fQcut[0]=0.1;//pi-pi, pi-k, pi-p
725 fQcut[2]=0.6;//the rest
726 fNormQcutLow[0] = 0.15;// 0.15 or 1.06
727 fNormQcutHigh[0] = 0.175;// 0.175 or 1.1
728 fNormQcutLow[1] = 1.34;//1.34
729 fNormQcutHigh[1] = 1.4;//1.4
730 fNormQcutLow[2] = 1.1;//1.1
731 fNormQcutHigh[2] = 1.4;//1.4
734 fMultLimit=kMultLimitpp;
739 fNormQcutLow[0] = 1.0;
740 fNormQcutHigh[0] = 1.5;
741 fNormQcutLow[1] = 1.0;
742 fNormQcutHigh[1] = 1.5;
743 fNormQcutLow[2] = 1.0;
744 fNormQcutHigh[2] = 1.5;
747 fQLowerCut = 0.005;// was 0.005
751 fKmeanY[0] = 0;// central y
754 // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
756 fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
757 fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
758 fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
760 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
762 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
763 fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
764 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
766 // 2x1 (Kt: 0-0.35, 0.35-1.0)
768 fKstepT[0] = 0.35; fKstepT[1] = 0.65;
769 fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
770 fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
774 fQupperBoundWeights = 0.2;
776 fQstep = fQupperBound/Float_t(kQbins);
777 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
778 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
780 fDampStart = 0.5;// was 0.3
786 fEC = new AliChaoticityEventCollection **[fZvertexBins];
787 for(UShort_t i=0; i<fZvertexBins; i++){
789 fEC[i] = new AliChaoticityEventCollection *[fMbins];
791 for(UShort_t j=0; j<fMbins; j++){
793 fEC[i][j] = new AliChaoticityEventCollection(fEventsToMix+1, fMultLimit, kPairLimit, kMCarrayLimit, fMCcase);
798 for(Int_t i=0; i<fMultLimit; i++) fDefaultsCharMult[i]='0';
799 for(Int_t i=0; i<kPairLimit; i++) fDefaultsCharSE[i]='0';
800 for(Int_t i=0; i<2*kPairLimit; i++) fDefaultsCharME[i]='0';
801 for(Int_t i=0; i<fMultLimit; i++) fDefaultsInt[i]=-1;
802 for(Int_t i=0; i<fMultLimit; i++) fPairLocationSE[i] = new TArrayI(fMultLimit,fDefaultsInt);
803 for(Int_t i=0; i<fMultLimit; i++) fPairLocationME[i] = new TArrayI(fMultLimit,fDefaultsInt);
804 for(Int_t i=0; i<kPairLimit; i++) fTripletSkip1[i] = new TArrayC(fMultLimit,fDefaultsCharSE);
805 for(Int_t i=0; i<2*kPairLimit; i++) fTripletSkip2[i] = new TArrayC(fMultLimit,fDefaultsCharME);
808 // Normalization utilities
809 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
810 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
811 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
812 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
813 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
814 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
815 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
817 // Track Merging/Splitting utilities
818 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P11
819 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P12
820 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P13
821 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[3][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P23
824 fNormPairs[0] = new AliChaoticityNormPairStruct[kNormPairLimit];
825 fNormPairs[1] = new AliChaoticityNormPairStruct[kNormPairLimit];
828 fTempStruct = new AliChaoticityTrackStruct[fMultLimit];
831 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
835 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
837 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
838 if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
839 if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
842 /////////////////////////////////////////////
843 /////////////////////////////////////////////
846 //________________________________________________________________________
847 void AliChaoticity::UserCreateOutputObjects()
852 ParInit();// Initialize my settings
855 fOutputList = new TList();
856 fOutputList->SetOwner();
858 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
859 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
860 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
861 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
862 fOutputList->Add(fVertexDist);
865 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
866 fOutputList->Add(fDCAxyDistPlus);
867 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
868 fOutputList->Add(fDCAzDistPlus);
869 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
870 fOutputList->Add(fDCAxyDistMinus);
871 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
872 fOutputList->Add(fDCAzDistMinus);
875 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
876 fOutputList->Add(fEvents1);
877 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
878 fOutputList->Add(fEvents2);
880 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
881 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
882 fOutputList->Add(fMultDist1);
884 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
885 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
886 fOutputList->Add(fMultDist2);
888 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
889 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
890 fOutputList->Add(fMultDist3);
892 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
893 fOutputList->Add(fPtEtaDist);
895 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
896 fOutputList->Add(fPhiPtDist);
898 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
899 fOutputList->Add(fTOFResponse);
900 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
901 fOutputList->Add(fTPCResponse);
903 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
904 fOutputList->Add(fRejectedPairs);
905 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
906 fOutputList->Add(fRejectedEvents);
908 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
909 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
910 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
911 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
912 TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
913 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
914 TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
915 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
916 TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
917 if(fMCcase) fOutputList->Add(fPairsPadRowNum);
918 TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
919 if(fMCcase) fOutputList->Add(fPairsPadRowDen);
923 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
924 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
925 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
926 if(fMCcase) fOutputList->Add(fAllOSPairs);
928 TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
929 if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
930 TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
931 if(fMCcase) fOutputList->Add(fAllSCPionPairs);
932 TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
933 if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
934 TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
935 if(fMCcase) fOutputList->Add(fAllMCPionPairs);
937 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
938 fOutputList->Add(fAvgMult);
940 TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
941 fOutputList->Add(fTrackChi2NDF);
942 TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
943 fOutputList->Add(fTrackTPCncls);
946 TH3D *fTPNRejects1 = new TH3D("fTPNRejects1","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
947 fOutputList->Add(fTPNRejects1);
948 TH3D *fTPNRejects2 = new TH3D("fTPNRejects2","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
949 fOutputList->Add(fTPNRejects2);
950 TH3D *fTPNRejects3 = new TH3D("fTPNRejects3","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
951 fOutputList->Add(fTPNRejects3);
952 TH3D *fTPNRejects4 = new TH3D("fTPNRejects4","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
953 fOutputList->Add(fTPNRejects4);
954 TH3D *fTPNRejects5 = new TH3D("fTPNRejects5","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
955 fOutputList->Add(fTPNRejects5);
958 TH3D *fKt3DistTerm1 = new TH3D("fKt3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
959 TH3D *fKt3DistTerm5 = new TH3D("fKt3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
960 fOutputList->Add(fKt3DistTerm1);
961 fOutputList->Add(fKt3DistTerm5);
963 TProfile2D *fKtTripletAvg = new TProfile2D("fKtTripletAvg","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
964 fOutputList->Add(fKtTripletAvg);
966 TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
967 TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
968 TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
969 TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
970 TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
971 TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
972 TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
973 TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
974 TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
975 TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
976 TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
977 TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
978 fOutputList->Add(fMCWeight3DTerm1SC);
979 fOutputList->Add(fMCWeight3DTerm1SCden);
980 fOutputList->Add(fMCWeight3DTerm2SC);
981 fOutputList->Add(fMCWeight3DTerm2SCden);
982 fOutputList->Add(fMCWeight3DTerm1MC);
983 fOutputList->Add(fMCWeight3DTerm1MCden);
984 fOutputList->Add(fMCWeight3DTerm2MC);
985 fOutputList->Add(fMCWeight3DTerm2MCden);
986 fOutputList->Add(fMCWeight3DTerm3MC);
987 fOutputList->Add(fMCWeight3DTerm3MCden);
988 fOutputList->Add(fMCWeight3DTerm4MC);
989 fOutputList->Add(fMCWeight3DTerm4MCden);
992 if(fPdensityExplicitLoop || fPdensityPairCut){
994 for(Int_t mb=0; mb<fMbins; mb++){
995 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
997 for(Int_t edB=0; edB<fEDbins; edB++){
998 for(Int_t c1=0; c1<2; c1++){
999 for(Int_t c2=0; c2<2; c2++){
1000 for(Int_t sc=0; sc<kSCLimit2; sc++){
1001 for(Int_t term=0; term<2; term++){
1003 TString *nameEx2 = new TString("Explicit2_Charge1_");
1005 nameEx2->Append("_Charge2_");
1007 nameEx2->Append("_SC_");
1009 nameEx2->Append("_M_");
1011 nameEx2->Append("_ED_");
1013 nameEx2->Append("_Term_");
1016 if(sc==0 || sc==3 || sc==5){
1017 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
1020 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
1021 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
1022 TString *nameEx2QW=new TString(nameEx2->Data());
1023 nameEx2QW->Append("_QW");
1024 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
1025 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW);
1026 TString *nameAvgP=new TString(nameEx2->Data());
1027 nameAvgP->Append("_AvgP");
1028 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, 400,0,2, 0.,1.0,"");
1029 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fAvgP);
1031 // Momentum resolution histos
1032 if(fMCcase && sc==0){
1033 TString *nameIdeal = new TString(nameEx2->Data());
1034 nameIdeal->Append("_Ideal");
1035 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",fRVALUES*kNDampValues,-0.5,fRVALUES*kNDampValues-0.5, kQbinsWeights,0,fQupperBoundWeights);
1036 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
1037 TString *nameSmeared = new TString(nameEx2->Data());
1038 nameSmeared->Append("_Smeared");
1039 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",fRVALUES*kNDampValues,-0.5,fRVALUES*kNDampValues-0.5, kQbinsWeights,0,fQupperBoundWeights);
1040 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
1042 TString *nameEx2MC=new TString(nameEx2->Data());
1043 nameEx2MC->Append("_MCqinv");
1044 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"",400,0,2);
1045 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
1046 TString *nameEx2MCQW=new TString(nameEx2->Data());
1047 nameEx2MCQW->Append("_MCqinvQW");
1048 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"",400,0,2);
1049 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
1051 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
1052 nameEx2PIDpurityDen->Append("_PIDpurityDen");
1053 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
1054 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
1055 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
1056 nameEx2PIDpurityNum->Append("_PIDpurityNum");
1057 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH2D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
1058 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
1062 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1063 nameEx2OSLB1->Append("_osl_b1");
1064 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1065 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL);
1066 nameEx2OSLB1->Append("_QW");
1067 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1068 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW);
1070 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1071 nameEx2OSLB2->Append("_osl_b2");
1072 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1073 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL);
1074 nameEx2OSLB2->Append("_QW");
1075 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1076 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW);
1083 // skip 3-particle if Tabulate6DPairs is true
1084 if(fTabulatePairs) continue;
1086 for(Int_t c3=0; c3<2; c3++){
1087 for(Int_t sc=0; sc<kSCLimit3; sc++){
1088 for(Int_t term=0; term<5; term++){
1089 TString *nameEx3 = new TString("Explicit3_Charge1_");
1091 nameEx3->Append("_Charge2_");
1093 nameEx3->Append("_Charge3_");
1095 nameEx3->Append("_SC_");
1097 nameEx3->Append("_M_");
1099 nameEx3->Append("_ED_");
1101 nameEx3->Append("_Term_");
1104 TString *namePC3 = new TString("PairCut3_Charge1_");
1106 namePC3->Append("_Charge2_");
1108 namePC3->Append("_Charge3_");
1110 namePC3->Append("_SC_");
1112 namePC3->Append("_M_");
1114 namePC3->Append("_ED_");
1116 namePC3->Append("_Term_");
1119 ///////////////////////////////////////
1120 // skip degenerate histograms
1121 if(sc==0 || sc==6 || sc==9){// Identical species
1122 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1123 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1125 if( (c1+c2)==1) {if(c1!=0) continue;}
1126 }else {}// do nothing for pi-k-p case
1128 /////////////////////////////////////////
1132 if(fPdensityExplicitLoop){
1133 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = new TH1D(nameEx3->Data(),"Three Particle Distribution",200,0,2);
1134 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3);
1136 nameEx3->Append("_Norm");
1137 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = new TH1D(nameEx3->Data(),"Explicit_3 Norm",1,-0.5,0.5);
1138 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3);
1140 if(fPdensityPairCut){
1141 TString *nameNorm=new TString(namePC3->Data());
1142 nameNorm->Append("_Norm");
1143 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1144 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1147 TString *name3DQ=new TString(namePC3->Data());
1148 name3DQ->Append("_3D");
1149 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH3D(name3DQ->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1150 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1153 const int NEdgesPos=16;
1154 double lowEdges4vectPos[NEdgesPos]={0};
1155 lowEdges4vectPos[0]=0.0;
1156 lowEdges4vectPos[1]=0.00005;// best resolution at low Q^2
1157 for(int edge=2; edge<NEdgesPos; edge++){
1158 lowEdges4vectPos[edge] = lowEdges4vectPos[edge-1] + lowEdges4vectPos[1]*(edge);
1160 const int NEdges=2*NEdgesPos-1;
1161 double lowEdges4vect[NEdges]={0};
1162 for(int edge=0; edge<NEdges; edge++){
1163 if(edge<NEdgesPos-1) lowEdges4vect[edge] = -lowEdges4vectPos[NEdgesPos-1-edge];
1164 else if(edge==NEdgesPos-1) lowEdges4vect[edge] = 0;
1165 else lowEdges4vect[edge] = lowEdges4vectPos[edge-NEdgesPos+1];
1166 //if(c1==c2 && c1==c3) cout<<lowEdges4vect[edge]<<endl;
1169 if(c1==c2 && c1==c3 && sc==0 && fMCcase==kFALSE){
1170 TString *name4vect1=new TString(namePC3->Data());
1171 TString *name4vect2=new TString(namePC3->Data());
1172 name4vect1->Append("_4VectProd1");
1173 name4vect2->Append("_4VectProd2");
1174 // use 3.75e6 MeV^4 as the resolution on QprodSum
1175 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms = new TH3D(name4vect1->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1176 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms);
1177 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms = new TH3D(name4vect2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1178 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms);
1180 if(sc==0 && fMCcase==kTRUE){
1181 TString *name3DMomResIdeal=new TString(namePC3->Data());
1182 name3DMomResIdeal->Append("_Ideal");
1183 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH3D(name3DMomResIdeal->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1184 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1185 TString *name3DMomResSmeared=new TString(namePC3->Data());
1186 name3DMomResSmeared->Append("_Smeared");
1187 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH3D(name3DMomResSmeared->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1188 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared);
1190 TString *name3DMomResQW12=new TString(namePC3->Data());
1191 name3DMomResQW12->Append("_QW12");
1192 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW12 = new TH3D(name3DMomResQW12->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1193 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW12);
1194 TString *name3DMomResQW13=new TString(namePC3->Data());
1195 name3DMomResQW13->Append("_QW13");
1196 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW13 = new TH3D(name3DMomResQW13->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1197 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW13);
1200 TString *name3DSumK3=new TString(namePC3->Data());
1201 name3DSumK3->Append("_SumK3");
1202 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSumK3 = new TH3D(name3DSumK3->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1203 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSumK3);
1204 TString *name3DEnK3=new TString(namePC3->Data());
1205 name3DEnK3->Append("_EnK3");
1206 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fEnK3 = new TH3D(name3DEnK3->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1207 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fEnK3);
1210 /*if(c1==c2 && c1==c3){
1211 TString *name4vect1Ideal=new TString(namePC3->Data());
1212 TString *name4vect1Smeared=new TString(namePC3->Data());
1213 TString *name4vect2Ideal=new TString(namePC3->Data());
1214 TString *name4vect2Smeared=new TString(namePC3->Data());
1215 TString *name4vect1Q3W=new TString(namePC3->Data());
1216 TString *name4vect2Q3W=new TString(namePC3->Data());
1217 name4vect1Ideal->Append("_4VectProd1Ideal");
1218 name4vect1Smeared->Append("_4VectProd1Smeared");
1219 name4vect2Ideal->Append("_4VectProd2Ideal");
1220 name4vect2Smeared->Append("_4VectProd2Smeared");
1221 name4vect1Q3W->Append("_4VectProd1Q3W");
1222 name4vect2Q3W->Append("_4VectProd2Q3W");
1224 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal = new TH3D(name4vect1Ideal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1225 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal);
1226 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared = new TH3D(name4vect1Smeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1227 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared);
1229 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal = new TH3D(name4vect2Ideal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1230 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal);
1231 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared = new TH3D(name4vect2Smeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1232 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared);
1234 if(term==0){// average Q3 in each FVP cell
1235 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Q3W = new TH3D(name4vect1Q3W->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1236 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Q3W);
1237 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Q3W = new TH3D(name4vect2Q3W->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1238 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Q3W);
1242 TString *name4vect1SumK3=new TString(namePC3->Data());
1243 TString *name4vect2SumK3=new TString(namePC3->Data());
1244 TString *name4vect1EnK3=new TString(namePC3->Data());
1245 TString *name4vect2EnK3=new TString(namePC3->Data());
1246 name4vect1SumK3->Append("_4VectProd1SumK3");
1247 name4vect2SumK3->Append("_4VectProd2SumK3");
1248 name4vect1EnK3->Append("_4VectProd1EnK3");
1249 name4vect2EnK3->Append("_4VectProd2EnK3");
1250 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3 = new TH3D(name4vect1SumK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1251 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3);
1252 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3 = new TH3D(name4vect2SumK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1253 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3);
1254 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3 = new TH3D(name4vect1EnK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1255 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3);
1256 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3 = new TH3D(name4vect2EnK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1257 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3);
1259 if(term > 0 && term < 4){
1260 TString *name4vect1SumK2=new TString(namePC3->Data());
1261 TString *name4vect2SumK2=new TString(namePC3->Data());
1262 TString *name4vect1EnK2=new TString(namePC3->Data());
1263 TString *name4vect2EnK2=new TString(namePC3->Data());
1264 name4vect1SumK2->Append("_4VectProd1SumK2");
1265 name4vect2SumK2->Append("_4VectProd2SumK2");
1266 name4vect1EnK2->Append("_4VectProd1EnK2");
1267 name4vect2EnK2->Append("_4VectProd2EnK2");
1268 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2 = new TH3D(name4vect1SumK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1269 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2);
1270 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2 = new TH3D(name4vect2SumK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1271 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2);
1272 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2 = new TH3D(name4vect1EnK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1273 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2);
1274 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2 = new TH3D(name4vect2EnK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1275 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2);
1280 if(c1==c2 && c1==c3 && term==4 && sc==0){
1281 for(Int_t dt=0; dt<kDENtypes; dt++){
1282 TString *nameDenType=new TString("PairCut3_Charge1_");
1284 nameDenType->Append("_Charge2_");
1286 nameDenType->Append("_Charge3_");
1288 nameDenType->Append("_SC_");
1290 nameDenType->Append("_M_");
1292 nameDenType->Append("_ED_");
1293 *nameDenType += edB;
1294 nameDenType->Append("_TPN_");
1297 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1298 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm);
1299 // neglect errors for TPN
1300 //nameDenType->Append("_Err");
1301 //Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1302 //fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr);
1304 /*TString *name4vect1TPN=new TString(nameDenType->Data());
1305 TString *name4vect2TPN=new TString(nameDenType->Data());
1306 name4vect1TPN->Append("_4VectProd1");
1307 name4vect2TPN->Append("_4VectProd2");
1308 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = new TH3D(name4vect1TPN->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1309 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm);
1310 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = new TH3D(name4vect2TPN->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1311 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm);
1314 TString *name4vect1TPNIdeal=new TString(nameDenType->Data());
1315 TString *name4vect2TPNIdeal=new TString(nameDenType->Data());
1316 TString *name4vect1TPNSmeared=new TString(nameDenType->Data());
1317 TString *name4vect2TPNSmeared=new TString(nameDenType->Data());
1318 name4vect1TPNIdeal->Append("_4VectProd1Ideal");
1319 name4vect2TPNIdeal->Append("_4VectProd2Ideal");
1320 name4vect1TPNSmeared->Append("_4VectProd1Smeared");
1321 name4vect2TPNSmeared->Append("_4VectProd2Smeared");
1322 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = new TH3D(name4vect1TPNIdeal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1323 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal);
1324 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = new TH3D(name4vect2TPNIdeal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1325 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal);
1326 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = new TH3D(name4vect1TPNSmeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1327 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared);
1328 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = new TH3D(name4vect2TPNSmeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1329 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared);
1335 }// c and sc exclusion
1349 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1350 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1351 for(Int_t mb=0; mb<fMbins; mb++){
1352 for(Int_t edB=0; edB<fEDbins; edB++){
1354 TString *nameNum = new TString("TwoPart_num_Kt_");
1356 nameNum->Append("_Ky_");
1358 nameNum->Append("_M_");
1360 nameNum->Append("_ED_");
1363 TString *nameDen = new TString("TwoPart_den_Kt_");
1365 nameDen->Append("_Ky_");
1367 nameDen->Append("_M_");
1369 nameDen->Append("_ED_");
1373 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1374 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD);
1376 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1377 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD);
1386 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1387 fOutputList->Add(fQsmearMean);
1388 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1389 fOutputList->Add(fQsmearSq);
1390 TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
1391 fOutputList->Add(fQDist);
1395 ////////////////////////////////////
1396 ///////////////////////////////////
1398 PostData(1, fOutputList);
1401 //________________________________________________________________________
1402 void AliChaoticity::Exec(Option_t *)
1405 // Called for each event
1406 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1409 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1411 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1412 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1416 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1417 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1418 if(!isSelected1 && !fMCcase) {return;}
1419 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1420 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1421 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1422 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1425 ///////////////////////////////////////////////////////////
1426 const AliAODVertex *primaryVertexAOD;
1427 AliCentrality *centrality;// for AODs and ESDs
1430 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1431 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1432 fPIDResponse = inputHandler->GetPIDResponse();
1435 TClonesArray *mcArray = 0x0;
1438 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1439 if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1440 cout<<"No MC particle branch found or Array too large!!"<<endl;
1448 Int_t positiveTracks=0, negativeTracks=0;
1449 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1451 Double_t vertex[3]={0};
1453 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1454 /////////////////////////////////////////////////
1457 Float_t centralityPercentile=0;
1458 Float_t cStep=5.0, cStart=0;
1461 if(fAODcase){// AOD case
1464 centrality = fAOD->GetCentrality();
1465 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1466 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1467 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
1468 cout<<"Centrality % = "<<centralityPercentile<<endl;
1474 ////////////////////////////////
1476 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1477 primaryVertexAOD = fAOD->GetPrimaryVertex();
1478 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1480 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1481 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1483 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1484 if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1486 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1488 fBfield = fAOD->GetMagneticField();
1490 for(Int_t i=0; i<fZvertexBins; i++){
1491 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1499 /////////////////////////////
1500 // Create Shuffled index list
1501 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1502 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1503 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1504 /////////////////////////////
1508 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1509 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1510 if (!aodtrack) continue;
1511 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1513 status=aodtrack->GetStatus();
1515 if(!aodtrack->TestFilterBit(BIT(7))) continue;// AOD filterBit cut
1518 if(fFilterBit != 7){
1519 Bool_t goodTrackOtherFB = kFALSE;
1520 for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1521 AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1522 if(!aodtrack2) continue;
1523 if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1525 if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1528 if(!goodTrackOtherFB) continue;
1532 if(aodtrack->Pt() < 0.16) continue;
1533 if(fabs(aodtrack->Eta()) > 0.8) continue;
1536 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1537 if(!goodMomentum) continue;
1538 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1541 Double_t dca2[2]={0};
1542 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1543 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1544 Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1546 fTempStruct[myTracks].fStatus = status;
1547 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1548 fTempStruct[myTracks].fId = aodtrack->GetID();
1549 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1550 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1551 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1552 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1553 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1554 fTempStruct[myTracks].fEta = aodtrack->Eta();
1555 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1556 fTempStruct[myTracks].fDCAXY = dca2[0];
1557 fTempStruct[myTracks].fDCAZ = dca2[1];
1558 fTempStruct[myTracks].fDCA = dca3d;
1559 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1560 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1564 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1565 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1566 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1570 fTempStruct[myTracks].fElectron = kFALSE;
1571 fTempStruct[myTracks].fPion = kFALSE;
1572 fTempStruct[myTracks].fKaon = kFALSE;
1573 fTempStruct[myTracks].fProton = kFALSE;
1575 Float_t nSigmaTPC[5];
1576 Float_t nSigmaTOF[5];
1577 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1578 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1579 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1580 Float_t signalTPC=0, signalTOF=0;
1581 Double_t integratedTimesTOF[10]={0};
1583 /*if(fFilterBit != 7 ) {
1584 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron));
1585 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon));
1586 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
1587 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
1588 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
1590 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron));
1591 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon));
1592 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion));
1593 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon));
1594 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton));
1595 signalTPC = aodtrack->GetTPCsignal();
1596 if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1597 fTempStruct[myTracks].fTOFhit = kTRUE;
1598 signalTOF = aodtrack->GetTOFsignal();
1599 aodtrack->GetIntegratedTimes(integratedTimesTOF);
1600 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1602 }else {// FilterBit 7 PID workaround
1604 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1605 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1606 if (!aodTrack2) continue;
1607 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1609 UInt_t status2=aodTrack2->GetStatus();
1611 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1612 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1613 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1614 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1615 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1617 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1618 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1619 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1620 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1621 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1622 signalTPC = aodTrack2->GetTPCsignal();
1624 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1625 fTempStruct[myTracks].fTOFhit = kTRUE;
1626 signalTOF = aodTrack2->GetTOFsignal();
1627 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1628 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1631 //}// FilterBit 7 PID workaround
1635 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1636 if(fTempStruct[myTracks].fTOFhit) {
1637 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1641 // Use TOF if good hit and above threshold
1642 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1643 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1644 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1645 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1646 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1647 }else {// TPC info instead
1648 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1649 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1650 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1651 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1655 // Ensure there is only 1 candidate per track
1656 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1657 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1658 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1659 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1660 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1661 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1662 ////////////////////////
1663 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1665 if(!fTempStruct[myTracks].fPion) continue;// only pions
1670 if(fTempStruct[myTracks].fPion) {// pions
1671 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1672 fTempStruct[myTracks].fKey = 1;
1673 }else if(fTempStruct[myTracks].fKaon){// kaons
1674 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1675 fTempStruct[myTracks].fKey = 10;
1677 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1678 fTempStruct[myTracks].fKey = 100;
1682 ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1683 ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1684 if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1685 if(aodtrack->GetTPCncls() < fMinTPCncls) continue;
1687 if(fTempStruct[myTracks].fCharge==+1) {
1688 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1689 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1691 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1692 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1695 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1696 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1699 if(aodtrack->Charge() > 0) positiveTracks++;
1700 else negativeTracks++;
1702 if(fTempStruct[myTracks].fPion) pionCount++;
1703 if(fTempStruct[myTracks].fKaon) kaonCount++;
1704 if(fTempStruct[myTracks].fProton) protonCount++;
1709 }else {// ESD tracks
1710 cout<<"ESDs not supported currently"<<endl;
1716 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1720 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1721 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
1723 /////////////////////////////////////////
1724 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1725 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1726 /////////////////////////////////////////
1729 ////////////////////////////////
1730 ///////////////////////////////
1731 // Mbin determination
1733 // Mbin set to Pion Count Only for pp!!!!!!!
1736 for(Int_t i=0; i<kMultBinspp; i++){
1737 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1738 // Mbin 0 has 1 pion
1741 for(Int_t i=0; i<fCentBins; i++){
1742 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1743 fMbin=i;// 0 = most central
1749 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1752 if(fMbin==0) fFSIbin = 0;//0-5%
1753 else if(fMbin==1) fFSIbin = 1;//5-10%
1754 else if(fMbin<=3) fFSIbin = 2;//10-20%
1755 else if(fMbin<=5) fFSIbin = 3;//20-30%
1756 else if(fMbin<=7) fFSIbin = 4;//30-40%
1757 else fFSIbin = 5;//40-50%
1759 Int_t rIndexForTPNMomRes = fRMax-6;
1760 if(fMbin==0) {rIndexForTPNMomRes=fRMax-6;}// 10 fm with EW (fRMax should be 11 for normal running)
1761 else if(fMbin==1) {rIndexForTPNMomRes=fRMax-7;}
1762 else if(fMbin<=3) {rIndexForTPNMomRes=fRMax-8;}
1763 else if(fMbin<=5) {rIndexForTPNMomRes=fRMax-9;}
1764 else {rIndexForTPNMomRes=fRMax-10;}
1766 //////////////////////////////////////////////////
1767 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1768 //////////////////////////////////////////////////
1772 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1773 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1775 ////////////////////////////////////
1776 // Add event to buffer if > 0 tracks
1778 fEC[zbin][fMbin]->FIFOShift();
1779 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1780 (fEvt)->fNtracks = myTracks;
1781 (fEvt)->fFillStatus = 1;
1782 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1784 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1785 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1786 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1787 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1788 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1789 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1790 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1797 Float_t qinv12=0, qinv13=0, qinv23=0;
1798 Float_t qout=0, qside=0, qlong=0;
1799 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1800 Float_t firstQ=0, secondQ=0, thirdQ=0;
1801 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
1802 Float_t firstKt=0, secondKt=0, thirdKt=0;
1803 Float_t transK12=0, rapK12=0, transK3=0;
1804 Int_t transKbin=0, rapKbin=0;
1805 Float_t q3=0, q3MC=0;
1806 Int_t ch1=0, ch2=0, ch3=0;
1807 Short_t key1=0, key2=0, key3=0;
1808 Int_t bin1=0, bin2=0, bin3=0;
1809 Float_t pVect1[4]={0};
1810 Float_t pVect2[4]={0};
1811 Float_t pVect3[4]={0};
1812 Float_t pVect1MC[4]={0};
1813 Float_t pVect2MC[4]={0};
1814 Float_t pVect3MC[4]={0};
1815 Int_t index1=0, index2=0, index3=0;
1816 Float_t weight12=0, weight13=0, weight23=0;
1817 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1818 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1819 Float_t weightTotal=0;//, weightTotalErr=0;
1820 Float_t qinv12MC=0, qinv13MC=0, qinv23MC=0;
1821 //Float_t Qsum1v1=0, Qsum2=0, Qsum3v1=0, Qsum1v2=0, Qsum3v2=0;
1822 //Float_t Qsum1v1MC=0, Qsum2MC=0, Qsum3v1MC=0, Qsum1v2MC=0, Qsum3v2MC=0;
1824 AliAODMCParticle *mcParticle1=0x0;
1825 AliAODMCParticle *mcParticle2=0x0;
1828 if(fPdensityPairCut){
1829 ////////////////////
1830 Int_t pairCountSE=0, pairCountME=0;
1831 Int_t normPairCount[2]={0};
1832 Int_t numOtherPairs1[2][fMultLimit];
1833 Int_t numOtherPairs2[2][fMultLimit];
1834 Bool_t exitCode=kFALSE;
1835 Int_t tempNormFillCount[2][2][2][10][5];
1838 // reset to defaults
1839 for(Int_t i=0; i<fMultLimit; i++) {
1840 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1841 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1843 // Normalization Utilities
1844 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1845 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1846 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1847 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1848 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1849 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1850 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1851 numOtherPairs1[0][i]=0;
1852 numOtherPairs1[1][i]=0;
1853 numOtherPairs2[0][i]=0;
1854 numOtherPairs2[1][i]=0;
1856 // Track Merging/Splitting Utilities
1857 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1858 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1859 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1860 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1863 // Reset the temp Normalization counters
1864 for(Int_t i=0; i<2; i++){// Charge1
1865 for(Int_t j=0; j<2; j++){// Charge2
1866 for(Int_t k=0; k<2; k++){// Charge3
1867 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1868 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1869 tempNormFillCount[i][j][k][l][m] = 0;
1877 ///////////////////////////////////////////////////////
1878 // Start the pairing process
1882 for (Int_t i=0; i<myTracks; i++) {
1887 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1889 key1 = (fEvt)->fTracks[i].fKey;
1890 key2 = (fEvt+en2)->fTracks[j].fKey;
1891 Short_t fillIndex2 = FillIndex2part(key1+key2);
1892 Short_t qCutBin = SetQcutBin(fillIndex2);
1893 Short_t normBin = SetNormBin(fillIndex2);
1894 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1895 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1896 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1897 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1901 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1902 GetQosl(pVect1, pVect2, qout, qside, qlong);
1903 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1908 ///////////////////////////////
1909 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1910 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1911 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1913 if(fMCcase && ch1==ch2 && fMbin==0 && qinv12<0.2){
1914 //////////////////////////
1915 // pad-row method testing
1916 Float_t coeff = (5)*0.2*(0.18/1.2);// 5 to evaluate at 1.0m in TPC
1917 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1918 if(phi1 > 2*PI) phi1 -= 2*PI;
1919 if(phi1 < 0) phi1 += 2*PI;
1920 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1921 if(phi2 > 2*PI) phi2 -= 2*PI;
1922 if(phi2 < 0) phi2 += 2*PI;
1923 Float_t deltaphi = phi1 - phi2;
1924 if(deltaphi > PI) deltaphi -= PI;
1925 if(deltaphi < -PI) deltaphi += PI;
1927 Int_t ncl1 = (fEvt)->fTracks[i].fClusterMap.GetNbits();
1928 Int_t ncl2 = (fEvt+en2)->fTracks[j].fClusterMap.GetNbits();
1929 Float_t sumCls = 0; Float_t sumSha = 0; Float_t sumQ = 0;
1930 Double_t shfrac = 0; //Double_t qfactor = 0;
1931 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
1932 if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Both clusters
1933 if ((fEvt)->fTracks[i].fSharedMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fSharedMap.TestBitNumber(imap)) { // Shared
1937 else {sumQ--; sumCls+=2;}
1939 else if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) || (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Non shared
1944 //qfactor = sumQ*1.0/sumCls;
1945 shfrac = sumSha*1.0/sumCls;
1947 if(fabs(deltaphi)<0.07 && fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
1948 ((TH3D*)fOutputList->FindObject("fPairsPadRowNum"))->Fill(transK12, shfrac, qinv12);
1951 for(Int_t rstep=0; rstep<10; rstep++){
1952 coeff = (rstep)*0.2*(0.18/1.2);
1953 phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1954 if(phi1 > 2*PI) phi1 -= 2*PI;
1955 if(phi1 < 0) phi1 += 2*PI;
1956 phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1957 if(phi2 > 2*PI) phi2 -= 2*PI;
1958 if(phi2 < 0) phi2 += 2*PI;
1959 deltaphi = phi1 - phi2;
1960 if(deltaphi > PI) deltaphi -= PI;
1961 if(deltaphi < -PI) deltaphi += PI;
1963 if(fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
1964 ((TH3F*)fOutputList->FindObject("fPairsShareFracDPhiNum"))->Fill(rstep, shfrac, deltaphi);
1966 //if(shfrac < 0.05){
1967 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1972 }// MCcase and pair selection
1974 // Pair Splitting/Merging cut
1975 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1977 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1978 fPairSplitCut[0][i]->AddAt('1',j);
1979 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1985 if(fMCcase && fillIndex2==0){
1987 // Check that label does not exceed stack size
1988 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1989 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1990 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1991 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1992 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1993 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1994 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1995 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1996 if(qinv12<0.1 && ch1==ch2) {
1997 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1998 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1999 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
2002 //if(transK12 <= 0.35) fEDbin=0;
2005 /*for(Int_t rIter=0; rIter<fRVALUES; rIter++){// 3fm to 8fm + 1 Therminator setting
2006 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
2007 Int_t denIndex = rIter*kNDampValues + myDampIt;
2008 Float_t WInput = MCWeight(ch1,ch2, rIter+kRmin, myDampIt, qinv12MC);
2009 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
2010 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
2011 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
2012 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
2017 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2018 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2020 //HIJING resonance test
2022 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
2023 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
2024 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
2025 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
2029 // secondary contamination
2030 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2032 ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
2033 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2034 ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
2037 ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
2038 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2039 ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
2044 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinv->Fill(qinv12MC, MCWeight(ch1,ch2,10,10,qinv12MC, 0.));// was 4,5
2045 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(ch1,ch2,10,10,qinv12MC, 0.));// was 4,5
2047 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityDen->Fill(transK12, qinv12);
2048 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
2049 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(transK12, qinv12);
2055 //////////////////////////////////////////
2057 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
2058 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
2059 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt)->fTracks[i].fMom);
2060 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt+en2)->fTracks[j].fMom);
2064 if((transK12 > 0.2) && (transK12 < 0.3)){
2065 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2066 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2068 if((transK12 > 0.6) && (transK12 < 0.7)){
2069 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2070 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2074 //////////////////////////////////////////
2076 if(fillIndex2==0 && bin1==bin2){
2078 transKbin=-1; rapKbin=-1;
2080 for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
2081 for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
2082 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2083 if((transKbin>=fKbinsT) || (rapKbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2084 Float_t WInput = 1.0;
2085 if(fGenerateSignal) {
2086 WInput = MCWeight(ch1,ch2, fRMax, fFixedLambdaBinr3, qinv12, transK12);
2087 //WInput = MCWeight(ch1,ch2, fRMax, fFixedLambdaBinr3, qinv12MC);
2088 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
2089 }else KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2096 // exit out of loop if there are too many pairs
2097 if(pairCountSE >= kPairLimit) {exitCode=kTRUE; continue;}// Too many SE pairs
2098 if(exitCode) continue;
2100 //////////////////////////
2102 if(qinv12 <= fQcut[qCutBin]) {
2104 ///////////////////////////
2106 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
2107 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
2108 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
2109 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
2110 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
2111 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
2112 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
2113 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
2114 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
2115 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
2116 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
2117 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
2120 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
2121 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
2122 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
2123 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
2124 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2125 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
2126 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
2127 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
2128 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
2129 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2130 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2131 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2134 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
2136 fPairLocationSE[i]->AddAt(pairCountSE,j);
2143 /////////////////////////////////////////////////////////
2144 // Normalization Region
2146 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2148 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
2149 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
2150 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
2152 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2153 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
2154 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
2157 //other past pairs with particle j
2158 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
2159 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
2160 if(locationOtherPair < 0) continue;// no pair there
2161 Int_t indexOther1 = i;
2162 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
2164 // Both possible orderings of other indexes
2165 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
2167 // 1 and 2 are from SE
2168 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
2169 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
2170 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2171 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2173 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
2176 }// pastpair P11 loop
2179 fNormPairSwitch[en2][i]->AddAt('1',j);
2180 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
2181 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
2183 numOtherPairs1[en2][i]++;
2184 numOtherPairs2[en2][j]++;
2187 normPairCount[en2]++;
2188 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2197 //////////////////////////////////////////////
2200 for (Int_t i=0; i<myTracks; i++) {
2205 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2207 key1 = (fEvt)->fTracks[i].fKey;
2208 key2 = (fEvt+en2)->fTracks[j].fKey;
2209 Short_t fillIndex2 = FillIndex2part(key1+key2);
2210 Short_t qCutBin = SetQcutBin(fillIndex2);
2211 Short_t normBin = SetNormBin(fillIndex2);
2212 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2213 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2214 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2215 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2217 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2218 GetQosl(pVect1, pVect2, qout, qside, qlong);
2219 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2220 //if(transK12 <= 0.35) fEDbin=0;
2225 ///////////////////////////////
2226 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2227 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2228 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
2231 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2232 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2233 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2234 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2235 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2236 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2237 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
2240 for(Int_t rIter=0; rIter<fRVALUES; rIter++){
2241 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
2242 Int_t denIndex = rIter*kNDampValues + myDampIt;
2243 Float_t WInput = MCWeight(ch1,ch2, rIter+kRmin, myDampIt, qinv12MC, 0.);
2244 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
2245 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
2246 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
2247 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
2252 /////////////////////////////////////////////////////
2253 if(!fTabulatePairs && qinv12 <= fQcut[qCutBin]) {// 3-particle MRC
2256 Short_t fillIndex3 = 0;
2257 key1=1; key2=1; key3=1;
2260 for (Int_t k=0; k<(fEvt+en3)->fNtracks; k++) {
2261 if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2262 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2263 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2264 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2265 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2266 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2267 qinv13 = GetQinv(0, pVect1, pVect3);
2268 qinv23 = GetQinv(0, pVect2, pVect3);
2270 if(qinv13 > fQcut[qCutBin] || qinv23 > fQcut[qCutBin]) continue;
2273 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2274 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2275 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2276 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2277 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
2278 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
2281 q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2282 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2286 // The below call to SetFillBins3 will work for all 3-particle terms since all are for fully mixed events. part is set to 1, but only matters for terms 2-4.
2287 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2288 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2291 for(Int_t jj=1; jj<=5; jj++){// term loop
2293 if(jj==2) {if(!fill2) continue;}//12
2294 if(jj==3) {if(!fill3) continue;}//13
2295 if(jj==4) {if(!fill4) continue;}//23
2299 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 1, jj, firstQ, secondQ, thirdQ);
2300 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 1, jj, firstQMC, secondQMC, thirdQMC);
2302 if(ch1==ch2 && ch1==ch3){// same charge
2303 WInput = MCWeight3D(kTRUE, jj, fFixedLambdaBinMomRes, firstQMC, secondQMC, thirdQMC, 0., 0., 0.);
2305 K3 = FSICorrelationTherm2(+1,+1, firstQMC)*FSICorrelationTherm2(+1,+1, secondQMC)*FSICorrelationTherm2(+1,+1, thirdQMC);// GRS
2306 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1SC"))->Fill(q3MC, WInput);
2307 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1SCden"))->Fill(q3MC);
2309 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2SC"))->Fill(q3MC, WInput);
2310 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2SCden"))->Fill(q3MC);
2312 }else {// mixed charge
2314 WInput = MCWeight3D(kFALSE, jj, fFixedLambdaBinMomRes, firstQMC, secondQMC, thirdQMC, 0., 0., 0.);
2315 if(jj==1) K3 = FSICorrelationTherm2(+1,+1, firstQMC)*FSICorrelationTherm2(+1,-1, secondQMC)*FSICorrelationTherm2(+1,-1, thirdQMC);// GRS
2317 if(jj==1 || jj==5) WInput = MCWeight3D(kFALSE, jj, fFixedLambdaBinMomRes, thirdQMC, secondQMC, firstQMC, 0., 0., 0.);// thirdQMC is ss
2318 else WInput = MCWeight3D(kFALSE, 6-jj, fFixedLambdaBinMomRes, thirdQMC, secondQMC, firstQMC, 0., 0., 0.);
2320 if(jj==1) K3 = FSICorrelationTherm2(+1,+1, thirdQMC)*FSICorrelationTherm2(+1,-1, secondQMC)*FSICorrelationTherm2(+1,-1, firstQMC);// GRS
2323 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1MC"))->Fill(q3MC, WInput);
2324 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1MCden"))->Fill(q3MC);
2328 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MC"))->Fill(q3MC, WInput);
2329 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MCden"))->Fill(q3MC);
2331 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MC"))->Fill(q3MC, WInput);
2332 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MCden"))->Fill(q3MC);
2334 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MC"))->Fill(q3MC, WInput);
2335 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MCden"))->Fill(q3MC);
2339 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MC"))->Fill(q3MC, WInput);
2340 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MCden"))->Fill(q3MC);
2342 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MC"))->Fill(q3MC, WInput);
2343 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MCden"))->Fill(q3MC);
2345 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MC"))->Fill(q3MC, WInput);
2346 ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MCden"))->Fill(q3MC);
2353 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2354 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2355 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fQW12->Fill(firstQMC, secondQMC, thirdQMC, WInput*firstQMC);
2356 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fQW13->Fill(firstQMC, secondQMC, thirdQMC, WInput*secondQMC);
2358 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSumK3->Fill(firstQMC, secondQMC, thirdQMC, WInput/K3);
2359 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fEnK3->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2362 /*if(ch1==ch2 && ch1==ch3){
2364 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2365 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2367 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2368 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2370 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2371 FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2373 FourVectProdTerms(pVect3, pVect1, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2374 FourVectProdTerms(pVect3MC, pVect1MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2376 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2377 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2380 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2381 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2382 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2383 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2385 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1Q3W->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput*q3);
2386 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2Q3W->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput*q3);
2389 if(qinv12MC > fQLowerCut && qinv13MC > fQLowerCut && qinv23MC > fQLowerCut){
2390 // does not really matter if MC or real data triplets are used to average 1/K3...but better to use umsmeared values
2392 WInput = MCWeight3D(kTRUE, 1, 25, firstQMC, secondQMC, thirdQMC);// pure 3-pion (lambda=1)
2393 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSumK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K3);
2394 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSumK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K3);
2395 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsEnK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2396 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsEnK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2398 Float_t InteractingQ=qinv12MC;
2399 Double_t K2 = FSICorrelationTherm2(+1,+1, InteractingQ);// K2 from Therminator source
2400 WInput = MCWeight3D(kTRUE, jj, 25, firstQMC, secondQMC, thirdQMC);// pure 2-pion (lambda=1)
2401 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSumK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K2);
2402 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSumK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K2);
2403 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsEnK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2404 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsEnK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2411 }// MCarray check, 3rd particle
2414 }// TabulatePairs Check
2416 }// MCarray check, 1st and 2nd particle
2418 // reset key's and fill bins (they were altered for 3 particle MRC calculation)
2419 key1 = (fEvt)->fTracks[i].fKey;
2420 key2 = (fEvt+en2)->fTracks[j].fKey;
2421 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
2424 if(ch1==ch2 && fMbin==0 && qinv12<0.2){
2425 //////////////////////////
2426 // pad-row method testing
2427 Float_t coeff = (5)*0.2*(0.18/1.2);// 5 to evaluate at 1.0m in TPC
2428 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2429 if(phi1 > 2*PI) phi1 -= 2*PI;
2430 if(phi1 < 0) phi1 += 2*PI;
2431 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2432 if(phi2 > 2*PI) phi2 -= 2*PI;
2433 if(phi2 < 0) phi2 += 2*PI;
2434 Float_t deltaphi = phi1 - phi2;
2435 if(deltaphi > PI) deltaphi -= PI;
2436 if(deltaphi < -PI) deltaphi += PI;
2438 Int_t ncl1 = (fEvt)->fTracks[i].fClusterMap.GetNbits();
2439 Int_t ncl2 = (fEvt+en2)->fTracks[j].fClusterMap.GetNbits();
2440 Float_t sumCls = 0; Float_t sumSha = 0; Float_t sumQ = 0;
2441 Double_t shfrac = 0; //Double_t qfactor = 0;
2442 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2443 if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Both clusters
2444 if ((fEvt)->fTracks[i].fSharedMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fSharedMap.TestBitNumber(imap)) { // Shared
2448 else {sumQ--; sumCls+=2;}
2450 else if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) || (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Non shared
2455 //qfactor = sumQ*1.0/sumCls;
2456 shfrac = sumSha*1.0/sumCls;
2458 if(fabs(deltaphi)<0.07 && fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
2459 ((TH3D*)fOutputList->FindObject("fPairsPadRowDen"))->Fill(transK12, shfrac, qinv12);
2462 for(Int_t rstep=0; rstep<10; rstep++){
2463 coeff = (rstep)*0.2*(0.18/1.2);
2464 // propagate through B field to r=1.2m
2465 phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2466 if(phi1 > 2*PI) phi1 -= 2*PI;
2467 if(phi1 < 0) phi1 += 2*PI;
2468 phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2469 if(phi2 > 2*PI) phi2 -= 2*PI;
2470 if(phi2 < 0) phi2 += 2*PI;
2471 deltaphi = phi1 - phi2;
2472 if(deltaphi > PI) deltaphi -= PI;
2473 if(deltaphi < -PI) deltaphi += PI;
2475 if(fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
2476 ((TH3F*)fOutputList->FindObject("fPairsShareFracDPhiDen"))->Fill(rstep, shfrac, deltaphi);
2478 //if(shfrac < 0.05){
2479 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2486 }// desired pair selection
2494 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2496 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2497 fPairSplitCut[1][i]->AddAt('1',j);
2502 //////////////////////////////////////////
2504 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
2505 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
2509 if((transK12 > 0.2) && (transK12 < 0.3)){
2510 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2511 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2513 if((transK12 > 0.6) && (transK12 < 0.7)){
2514 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2515 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2518 //////////////////////////////////////////
2520 if(fillIndex2==0 && bin1==bin2){
2522 transKbin=-1; rapKbin=-1;
2524 for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
2525 for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
2526 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2527 if((transKbin>=fKbinsT) || (rapKbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2529 if(fGenerateSignal) KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2530 else KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2537 if(pairCountME >= 2*kPairLimit) {exitCode=kTRUE; continue;}// Too many SE pairs
2538 if(exitCode) continue;
2540 if(qinv12 <= fQcut[qCutBin]) {
2541 ///////////////////////////
2544 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
2545 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
2546 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
2547 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
2548 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
2549 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
2550 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
2551 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
2552 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
2553 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
2554 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
2555 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
2558 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
2559 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
2560 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
2561 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
2562 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2563 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
2564 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
2565 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
2566 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
2567 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2568 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2569 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2572 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
2574 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
2580 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2582 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
2583 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
2584 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
2586 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2587 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
2588 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
2590 //other past pairs in P11 with particle i
2591 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
2592 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
2593 if(locationOtherPairP11 < 0) continue;// no pair there
2594 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
2596 //Check other past pairs in P12
2597 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
2599 // 1 and 3 are from SE
2600 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
2601 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
2602 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2603 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2604 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
2607 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
2608 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
2609 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
2615 fNormPairSwitch[en2][i]->AddAt('1',j);
2616 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
2617 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
2619 numOtherPairs1[en2][i]++;
2620 numOtherPairs2[en2][j]++;
2622 normPairCount[en2]++;
2623 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2632 ///////////////////////////////////////
2633 // P13 pairing (just for Norm counting of term5)
2634 for (Int_t i=0; i<myTracks; i++) {
2636 // exit out of loop if there are too many pairs
2637 // dont bother with this loop if exitCode is set.
2643 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2645 key1 = (fEvt)->fTracks[i].fKey;
2646 key2 = (fEvt+en2)->fTracks[j].fKey;
2647 Short_t fillIndex2 = FillIndex2part(key1+key2);
2648 Short_t normBin = SetNormBin(fillIndex2);
2649 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2650 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2651 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2652 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2654 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2656 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2658 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2659 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2662 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2663 fPairSplitCut[2][i]->AddAt('1',j);
2668 /////////////////////////////////////////////////////////
2669 // Normalization Region
2671 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2673 fNormPairSwitch[en2][i]->AddAt('1',j);
2681 ///////////////////////////////////////
2682 // P23 pairing (just for Norm counting of term5)
2684 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
2686 // exit out of loop if there are too many pairs
2687 // dont bother with this loop if exitCode is set.
2693 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2697 key1 = (fEvt+en1)->fTracks[i].fKey;
2698 key2 = (fEvt+en2)->fTracks[j].fKey;
2699 Short_t fillIndex2 = FillIndex2part(key1+key2);
2700 Short_t normBin = SetNormBin(fillIndex2);
2701 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2702 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2703 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2704 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2706 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2708 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2710 ///////////////////////////////
2711 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2712 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2715 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2716 fPairSplitCut[3][i]->AddAt('1',j);
2721 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2723 Int_t index1P23 = i;
2724 Int_t index2P23 = j;
2726 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2727 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2728 if(locationOtherPairP12 < 0) continue; // no pair there
2729 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2732 //Check other past pair status in P13
2733 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2735 // all from different event
2736 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2737 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2738 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2739 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2741 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2749 ///////////////////////////////////////////////////
2750 // Do not use pairs from events with too many pairs
2752 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2753 (fEvt)->fNpairsSE = 0;
2754 (fEvt)->fNpairsME = 0;
2755 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2756 return;// Skip event
2758 (fEvt)->fNpairsSE = pairCountSE;
2759 (fEvt)->fNpairsME = pairCountME;
2760 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2762 ///////////////////////////////////////////////////
2765 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
2766 //cout<<"Start Main analysis"<<endl;
2768 ///////////////////////////////////////////////////////////////////////
2769 ///////////////////////////////////////////////////////////////////////
2770 ///////////////////////////////////////////////////////////////////////
2773 // Start the Main Correlation Analysis
2776 ///////////////////////////////////////////////////////////////////////
2780 /////////////////////////////////////////////////////////
2781 // Skip 3-particle part if Tabulate6DPairs is set to true
2782 if(fTabulatePairs) return;
2783 /////////////////////////////////////////////////////////
2785 // Set the Normalization counters
2786 for(Int_t termN=0; termN<5; termN++){
2789 if((fEvt)->fNtracks ==0) continue;
2791 if((fEvt)->fNtracks ==0) continue;
2792 if((fEvt+1)->fNtracks ==0) continue;
2794 if((fEvt)->fNtracks ==0) continue;
2795 if((fEvt+1)->fNtracks ==0) continue;
2796 if((fEvt+2)->fNtracks ==0) continue;
2799 for(Int_t sc=0; sc<kSCLimit3; sc++){
2801 for(Int_t c1=0; c1<2; c1++){
2802 for(Int_t c2=0; c2<2; c2++){
2803 for(Int_t c3=0; c3<2; c3++){
2805 if(sc==0 || sc==6 || sc==9){// Identical species
2806 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2807 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2809 if( (c1+c2)==1) {if(c1!=0) continue;}
2810 }else {}// do nothing for pi-k-p case
2811 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2820 /////////////////////////////////////////////
2821 // Calculate Pair-Cut Correlations
2822 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2825 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2826 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2829 for(Int_t p1=0; p1<nump1; p1++){
2832 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2833 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2834 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2835 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2836 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2837 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2838 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2839 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2840 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2842 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2843 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2844 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2845 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2846 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2849 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2850 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2851 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2852 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2853 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2854 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2855 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2856 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2857 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2859 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2860 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2861 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2862 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2863 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2868 for(Int_t en2=0; en2<3; en2++){
2869 //////////////////////////////////////
2871 Bool_t skipcase=kTRUE;
2872 Short_t config=-1, part=-1;
2873 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2874 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2875 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2876 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2878 if(skipcase) continue;
2883 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2887 // remove auto-correlations and duplicate triplets
2889 if( index1 == index3) continue;
2890 if( index2 == index3) continue;
2891 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2893 // skip the switched off triplets
2894 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2895 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2898 ///////////////////////////////
2899 // Turn off 1st possible degenerate triplet
2900 if(index1 < index3){// verify correct id ordering ( index1 < k )
2901 if(fPairLocationSE[index1]->At(index3) >= 0){
2902 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2904 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2905 }else {// or k < index1
2906 if(fPairLocationSE[index3]->At(index1) >= 0){
2907 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2909 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2911 // turn off 2nd possible degenerate triplet
2912 if(index2 < index3){// verify correct id ordering (index2 < k)
2913 if(fPairLocationSE[index2]->At(index3) >= 0){
2914 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2916 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2917 }else {// or k < index2
2918 if(fPairLocationSE[index3]->At(index2) >= 0){
2919 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2921 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2926 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2927 ///////////////////////////////
2928 // Turn off 1st possible degenerate triplet
2929 if(fPairLocationME[index1]->At(index3) >= 0){
2930 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2933 // turn off 2nd possible degenerate triplet
2934 if(fPairLocationME[index2]->At(index3) >= 0){
2935 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2938 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2939 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2940 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2941 }// end config 2 part 1
2943 if(config==2 && part==2){// P12T1
2944 if( index1 == index3) continue;
2946 // skip the switched off triplets
2947 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2948 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2951 // turn off another possible degenerate
2952 if(fPairLocationME[index3]->At(index2) >= 0){
2953 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2954 }// end config 2 part 2
2956 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2957 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2958 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2959 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2961 if(config==3){// P12T3
2962 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2963 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2964 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2969 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2970 key3 = (fEvt+en2)->fTracks[k].fKey;
2971 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2972 Short_t fillIndex13 = FillIndex2part(key1+key3);
2973 Short_t fillIndex23 = FillIndex2part(key2+key3);
2974 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2975 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2976 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2977 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2978 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2979 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
2980 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2981 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2983 if(qinv13 < fQLowerCut) continue;
2984 if(qinv23 < fQLowerCut) continue;
2985 if(qinv13 > fQcut[qCutBin13]) continue;
2986 if(qinv23 > fQcut[qCutBin23]) continue;
2988 Float_t Kt12=sqrt( pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2989 Float_t Kt13=sqrt( pow(pVect1[1]+pVect3[1],2) + pow(pVect1[2]+pVect3[2],2))/2.;
2990 Float_t Kt23=sqrt( pow(pVect2[1]+pVect3[1],2) + pow(pVect2[2]+pVect3[2],2))/2.;
2993 pVect3MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPx;
2994 pVect3MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPy;
2995 pVect3MC[3] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPz;
2996 pVect3MC[0] = sqrt(pow(pVect3MC[1],2)+pow(pVect3MC[2],2)+pow(pVect3MC[3],2)+pow(fTrueMassPi,2));
2997 qinv12MC = GetQinv(0, pVect1MC, pVect2MC);
2998 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
2999 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
3004 // if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
3006 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
3007 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
3009 if(transK3<0.3) fEDbin=0;
3012 firstQ=0; secondQ=0; thirdQ=0;
3018 if(config==1) {// 123
3019 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
3021 if(fillIndex3 <= 2){
3022 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
3023 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, Kt12, Kt13, Kt23, 0, 1, firstKt, secondKt, thirdKt);
3024 if(fillIndex3==0 && fMCcase) ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 0, 1, firstQMC, secondQMC, thirdQMC);
3025 Float_t WInput = 1.0;
3026 if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, 1, fFixedLambdaBinr3, firstQ, secondQ, thirdQ, firstKt, secondKt, thirdKt);
3027 //if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, 1, fFixedLambdaBinr3, firstQMC, secondQMC, thirdQMC);
3030 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
3034 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
3035 ((TH3D*)fOutputList->FindObject("fKt3DistTerm1"))->Fill(fMbin+1, transK3, q3);
3037 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
3038 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
3039 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
3040 ((TProfile2D*)fOutputList->FindObject("fKtTripletAvg"))->Fill(fMbin+1, fEDbin, pt1);
3041 ((TProfile2D*)fOutputList->FindObject("fKtTripletAvg"))->Fill(fMbin+1, fEDbin, pt2);
3042 ((TProfile2D*)fOutputList->FindObject("fKtTripletAvg"))->Fill(fMbin+1, fEDbin, pt3);
3044 /*FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3045 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
3046 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
3052 }else if(config==2){// 12, 13, 23
3054 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
3055 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
3057 // loop over terms 2-4
3058 for(Int_t jj=2; jj<5; jj++){
3059 if(jj==2) {if(!fill2) continue;}//12
3060 if(jj==3) {if(!fill3) continue;}//13
3061 if(jj==4) {if(!fill4) continue;}//23
3063 if(fillIndex3 <= 2){
3064 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
3065 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, Kt12, Kt13, Kt23, part, jj, firstKt, secondKt, thirdKt);
3066 if(fillIndex3==0 && fMCcase) ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
3067 Float_t WInput = 1.0;
3068 if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, jj, fFixedLambdaBinr3, firstQ, secondQ, thirdQ, firstKt, secondKt, thirdKt);
3069 //if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, jj, fFixedLambdaBinr3, firstQMC, secondQMC, thirdQMC);
3071 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
3073 /*if(fillIndex3==0 && ch1==ch2 && ch1==ch3){
3074 if(part==1){// P11T2
3076 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3078 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3080 FourVectProdTerms(pVect3, pVect1, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3084 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3086 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3088 FourVectProdTerms(pVect2, pVect1, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3092 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
3093 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
3100 }else {// config 3: All particles from different events
3102 // "enhancement" differs from 1.0 only when Qinv goes over fQcut
3103 //Float_t enhancement=1.0;
3104 //Int_t nUnderCut=0;
3105 //if(qinv13<fQcut[qCutBin13]) nUnderCut++;
3106 //if(qinv23<fQcut[qCutBin23]) nUnderCut++;
3107 //if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
3108 //if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
3109 //if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
3111 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
3113 if(ch1==ch2 && ch1==ch3 && fillIndex3==0) {
3114 //FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
3115 if(!fMCcase) ((TH3D*)fOutputList->FindObject("fKt3DistTerm5"))->Fill(fMbin+1, transK3, q3);
3118 if(fillIndex3 <= 2){
3119 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
3120 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, Kt12, Kt13, Kt23, part, 5, firstKt, secondKt, thirdKt);
3121 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ);
3122 /*if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
3123 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
3124 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
3128 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
3129 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
3132 GetWeight(pVect1, pVect2, weight12, weight12Err);
3133 GetWeight(pVect1, pVect3, weight13, weight13Err);
3134 GetWeight(pVect2, pVect3, weight23, weight23Err);
3136 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {
3137 if(fMbin==0 && bin1==0) {
3138 ((TH3F*)fOutputList->FindObject("fTPNRejects1"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
3140 continue;// weight should never be larger than 1
3144 Float_t myDamp = fDampStart + (fDampStep)*fFixedLambdaBinr3;// lambdabin=0.52 for v1 draft, 0.7 is more realistic
3146 Int_t momResIndex = rIndexForTPNMomRes*kNDampValues + fFixedLambdaBinMomRes;// lambdabin=0.52 for v1 draft, 0.4 is more realistic
3148 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, qinv12);
3149 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, qinv13);
3150 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, qinv23);
3151 if(coulCorr12 < 0.1 || coulCorr13 < 0.1 || coulCorr23 < 0.1) {// Safety check
3152 if(fMbin==0 && bin1==0) {
3153 ((TH3F*)fOutputList->FindObject("fTPNRejects2"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
3157 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
3158 if(!fGenerateSignal && !fMCcase) {
3159 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
3160 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
3161 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
3162 if(momBin12 >= kQbins) momBin12 = kQbins-1;
3163 if(momBin13 >= kQbins) momBin13 = kQbins-1;
3164 if(momBin23 >= kQbins) momBin23 = kQbins-1;
3165 MomResCorr12 = fMomResC2->GetBinContent(momResIndex+1, momBin12);
3166 MomResCorr13 = fMomResC2->GetBinContent(momResIndex+1, momBin13);
3167 MomResCorr23 = fMomResC2->GetBinContent(momResIndex+1, momBin23);
3168 if(MomResCorr12 > 1.2 || MomResCorr13 > 1.2 || MomResCorr23 > 1.2) {// Safety check
3169 if(fMbin==0 && bin1==0) {
3170 ((TH3F*)fOutputList->FindObject("fTPNRejects3"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
3175 weight12CC = ((weight12+1)*MomResCorr12 - myDamp*coulCorr12 - (1-myDamp));
3176 weight12CC /= coulCorr12*myDamp;
3177 weight13CC = ((weight13+1)*MomResCorr13 - myDamp*coulCorr13 - (1-myDamp));
3178 weight13CC /= coulCorr13*myDamp;
3179 weight23CC = ((weight23+1)*MomResCorr23 - myDamp*coulCorr23 - (1-myDamp));
3180 weight23CC /= coulCorr23*myDamp;
3182 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {
3183 if(fMbin==0 && bin1==0) {
3184 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
3185 ((TH3F*)fOutputList->FindObject("fTPNRejects4"))->Fill(qinv12, qinv13, qinv23, weightTotal);
3187 continue;// C2^QS can never be less than unity
3190 /////////////////////////////////////////////////////
3191 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
3192 /////////////////////////////////////////////////////
3194 if(weightTotal > 1.5) {
3195 if(fMbin==0 && bin1==0) {
3196 ((TH3F*)fOutputList->FindObject("fTPNRejects5"))->Fill(qinv12, qinv13, qinv23, weightTotal);
3198 continue;// C2^QS never be greater than 1.0 in theory. Can be slightly larger than 1.0 with fluctuations
3203 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, weightTotal);
3205 //Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNorm->Fill(Qsum1v1, Qsum2, Qsum3v1, weightTotal);
3206 //Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNorm->Fill(Qsum1v2, Qsum2, Qsum3v2, weightTotal);
3209 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
3211 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
3212 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
3213 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
3214 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
3215 weightTotalErr /= pow(2*weightTotal,2);
3217 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, weightTotalErr);
3224 }// end 3rd particle
3232 }// end of PdensityPairs
3240 ////////////////////////////////////////////////////////
3241 // Pdensity Method with Explicit Loops
3242 if(fPdensityExplicitLoop){
3244 ////////////////////////////////////
3245 // 2nd, 3rd, and 4th order Correlations
3248 for (Int_t i=0; i<myTracks; i++) {
3249 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
3250 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
3251 pVect1[1] = (fEvt)->fTracks[i].fP[0];
3252 pVect1[2] = (fEvt)->fTracks[i].fP[1];
3253 pVect1[3] = (fEvt)->fTracks[i].fP[2];
3254 key1 = (fEvt)->fTracks[i].fKey;
3257 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
3259 if(en2==0) startbin2=i+1;
3262 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
3263 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
3264 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
3265 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
3266 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
3267 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
3268 key2 = (fEvt+en2)->fTracks[j].fKey;
3270 Short_t fillIndex12 = FillIndex2part(key1+key2);
3271 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
3273 if(qinv12 < fQLowerCut) continue;
3276 // 2-particle part is filled always during pair creator
3279 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
3281 if(en3==en2) startbin3=j+1;
3286 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
3287 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
3288 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
3289 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
3290 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
3291 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
3292 key3 = (fEvt+en3)->fTracks[k].fKey;
3294 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
3295 Short_t fillIndex13 = FillIndex2part(key1+key3);
3296 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
3297 Short_t fillIndex23 = FillIndex2part(key2+key3);
3298 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
3301 if(qinv13 < fQLowerCut) continue;
3302 if(qinv23 < fQLowerCut) continue;
3305 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
3307 Short_t normBin12 = SetNormBin(fillIndex12);
3308 Short_t normBin13 = SetNormBin(fillIndex13);
3309 Short_t normBin23 = SetNormBin(fillIndex23);
3312 if(en3==0 && en2==0) {// 123
3313 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
3315 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
3317 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3318 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3319 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
3323 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
3326 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
3327 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
3331 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
3332 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3333 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3334 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
3339 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
3340 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3341 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3342 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
3347 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
3348 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3349 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3350 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
3355 }else if(en2==1 && en3==2){// all uncorrelated events
3356 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
3358 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
3359 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3360 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3361 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
3364 Short_t qCutBin12 = SetQcutBin(fillIndex12);
3365 Short_t qCutBin13 = SetQcutBin(fillIndex13);
3366 Short_t qCutBin23 = SetQcutBin(fillIndex23);
3368 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
3371 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
3372 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
3373 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
3391 }// End of PdensityExplicit
3396 // Post output data.
3397 PostData(1, fOutputList);
3400 //________________________________________________________________________
3401 void AliChaoticity::Terminate(Option_t *)
3403 // Called once at the end of the query
3408 //________________________________________________________________________
3409 Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
3412 if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
3414 // propagate through B field to r=1m
3415 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
3416 if(phi1 > 2*PI) phi1 -= 2*PI;
3417 if(phi1 < 0) phi1 += 2*PI;
3418 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
3419 if(phi2 > 2*PI) phi2 -= 2*PI;
3420 if(phi2 < 0) phi2 += 2*PI;
3422 Float_t deltaphi = phi1 - phi2;
3423 if(deltaphi > PI) deltaphi -= 2*PI;
3424 if(deltaphi < -PI) deltaphi += 2*PI;
3425 deltaphi = fabs(deltaphi);
3427 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3430 // propagate through B field to r=1.6m
3431 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
3432 if(phi1 > 2*PI) phi1 -= 2*PI;
3433 if(phi1 < 0) phi1 += 2*PI;
3434 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
3435 if(phi2 > 2*PI) phi2 -= 2*PI;
3436 if(phi2 < 0) phi2 += 2*PI;
3438 deltaphi = phi1 - phi2;
3439 if(deltaphi > PI) deltaphi -= 2*PI;
3440 if(deltaphi < -PI) deltaphi += 2*PI;
3441 deltaphi = fabs(deltaphi);
3443 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
3449 Int_t ncl1 = first.fClusterMap.GetNbits();
3450 Int_t ncl2 = second.fClusterMap.GetNbits();
3451 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
3452 Double_t shfrac = 0; Double_t qfactor = 0;
3453 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
3454 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
3455 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
3459 else {sumQ--; sumCls+=2;}
3461 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
3466 qfactor = sumQ*1.0/sumCls;
3467 shfrac = sumSha*1.0/sumCls;
3470 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
3477 //________________________________________________________________________
3478 Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
3480 Float_t arg = G_Coeff/qinv;
3482 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
3483 else {return (exp(-arg)-1)/(-arg);}
3486 //________________________________________________________________________
3487 void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
3491 for (Int_t i = i1; i < i2+1; i++) {
3492 j = (Int_t) (gRandom->Rndm() * a);
3498 //________________________________________________________________________
3499 Short_t AliChaoticity::FillIndex2part(Short_t key){
3501 if(key==2) return 0;// pi-pi
3502 else if(key==11) return 1;// pi-k
3503 else if(key==101) return 2;// pi-p
3504 else if(key==20) return 3;// k-k
3505 else if(key==110) return 4;// k-p
3506 else return 5;// p-p
3508 //________________________________________________________________________
3509 Short_t AliChaoticity::FillIndex3part(Short_t key){
3511 if(key==3) return 0;// pi-pi-pi
3512 else if(key==12) return 1;// pi-pi-k
3513 else if(key==21) return 2;// k-k-pi
3514 else if(key==102) return 3;// pi-pi-p
3515 else if(key==201) return 4;// p-p-pi
3516 else if(key==111) return 5;// pi-k-p
3517 else if(key==30) return 6;// k-k-k
3518 else if(key==120) return 7;// k-k-p
3519 else if(key==210) return 8;// p-p-k
3520 else return 9;// p-p-p
3523 //________________________________________________________________________
3524 Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
3525 if(fi <= 2) return 0;
3526 else if(fi==3) return 1;
3529 //________________________________________________________________________
3530 Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
3532 else if(fi <= 2) return 1;
3535 //________________________________________________________________________
3536 void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3538 if(fi==0 || fi==3 || fi==5){// Identical species
3539 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3540 else {b1=c1; b2=c2;}
3541 }else {// Mixed species
3542 if(key1 < key2) { b1=c1; b2=c2;}
3543 else {b1=c2; b2=c1;}
3547 //________________________________________________________________________
3548 void AliChaoticity::SetFillBins3(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Short_t part, Int_t &b1, Int_t &b2, Int_t &b3, Bool_t &fill2, Bool_t &fill3, Bool_t &fill4){
3551 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
3552 // part only matters for terms 2-4
3555 Short_t seKeySum=0;// only used for pi-k-p case
3556 if(part==1) {// default case (irrelevant for term 1 and term 5)
3557 if(c1==c2) seSS=kTRUE;
3558 if(key1==key2) seSK=kTRUE;
3559 seKeySum = key1+key2;
3562 if(c1==c3) seSS=kTRUE;
3563 if(key1==key3) seSK=kTRUE;
3564 seKeySum = key1+key3;
3568 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3570 if(fi==0 || fi==6 || fi==9){// Identical species
3571 if( (c1+c2+c3)==1) {
3572 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3574 if(seSS) fill2=kTRUE;
3575 else {fill3=kTRUE; fill4=kTRUE;}
3577 }else if( (c1+c2+c3)==2) {
3580 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3584 b1=c1; b2=c2; b3=c3;
3585 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3587 }else if(fi != 5){// all the rest except pi-k-p
3590 if( (c1+c2)==1) {b1=0; b2=1;}
3591 else {b1=c1; b2=c2;}
3592 }else if(key1==key3){
3594 if( (c1+c3)==1) {b1=0; b2=1;}
3595 else {b1=c1; b2=c3;}
3596 }else {// Key2==Key3
3598 if( (c2+c3)==1) {b1=0; b2=1;}
3599 else {b1=c2; b2=c3;}
3601 //////////////////////////////
3602 if(seSK) fill2=kTRUE;// Same keys from Same Event
3603 else {// Different keys from Same Event
3604 if( (c1+c2+c3)==1) {
3606 if(seSS) fill3=kTRUE;
3608 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
3609 }else if( (c1+c2+c3)==2) {
3611 if(seSS) fill4=kTRUE;
3613 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
3614 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
3616 /////////////////////////////
3617 }else {// pi-k-p (no charge ordering applies since all are unique)
3619 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
3620 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
3622 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
3623 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
3625 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
3626 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
3628 ////////////////////////////////////
3629 if(seKeySum==11) fill2=kTRUE;
3630 else if(seKeySum==101) fill3=kTRUE;
3632 ////////////////////////////////////
3636 //________________________________________________________________________
3637 void AliChaoticity::ArrangeQs(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Float_t q12, Float_t q13, Float_t q23, Short_t part, Short_t term, Float_t &fQ, Float_t &sQ, Float_t &tQ){
3639 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
3640 if(fi==0 || fi==6 || fi==9){// Identical species
3641 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
3642 if(term==1 || term==5){
3643 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
3644 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
3645 else {fQ=q23; sQ=q12; tQ=q13;}
3646 }else if(term==2 && part==1){
3647 fQ=q12; sQ=q13; tQ=q23;
3648 }else if(term==2 && part==2){
3649 fQ=q13; sQ=q12; tQ=q23;
3650 }else if(term==3 && part==1){
3652 if(c1==c3) {fQ=q13; tQ=q23;}
3653 else {fQ=q23; tQ=q13;}
3654 }else if(term==3 && part==2){
3656 if(c1==c2) {fQ=q12; tQ=q23;}
3657 else {fQ=q23; tQ=q12;}
3658 }else if(term==4 && part==1){
3660 if(c1==c3) {fQ=q13; sQ=q23;}
3661 else {fQ=q23; sQ=q13;}
3662 }else if(term==4 && part==2){
3664 if(c1==c2) {fQ=q12; sQ=q23;}
3665 else {fQ=q23; sQ=q12;}
3666 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3667 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
3668 if(term==1 || term==5){
3669 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
3670 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
3671 else {tQ=q23; sQ=q12; fQ=q13;}
3672 }else if(term==2 && part==1){
3674 if(c1==c3) {tQ=q13; sQ=q23;}
3675 else {tQ=q23; sQ=q13;}
3676 }else if(term==2 && part==2){
3678 if(c1==c2) {tQ=q12; sQ=q23;}
3679 else {tQ=q23; sQ=q12;}
3680 }else if(term==3 && part==1){
3682 if(c1==c3) {tQ=q13; fQ=q23;}
3683 else {tQ=q23; fQ=q13;}
3684 }else if(term==3 && part==2){
3686 if(c1==c2) {tQ=q12; fQ=q23;}
3687 else {tQ=q23; fQ=q12;}
3688 }else if(term==4 && part==1){
3689 tQ=q12; sQ=q13; fQ=q23;
3690 }else if(term==4 && part==2){
3691 tQ=q13; sQ=q12; fQ=q23;
3692 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3693 }else {// fQ=ss, sQ=ss, tQ=ss
3694 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
3695 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
3696 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
3697 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
3698 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
3699 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
3700 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
3702 }else if(fi != 5){// all the rest except pi-k-p
3706 // cases not explicity shown below are not possible
3707 if(term==1 || term==5) {sQ=q13; tQ=q23;}
3708 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
3709 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
3710 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
3711 else cout<<"problem!!!!!!!!!!!!!"<<endl;
3713 if(c1==c3) {sQ=q13; tQ=q23;}
3714 else {sQ=q23; tQ=q13;}
3716 if(c1==c3) {tQ=q13; sQ=q23;}
3717 else {tQ=q23; sQ=q13;}
3719 }else if(key1==key3){
3722 // cases not explicity shown below are not possible
3723 if(term==1 || term==5) {sQ=q12; tQ=q23;}
3724 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
3725 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
3726 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
3727 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3729 if(c1==c2) {sQ=q12; tQ=q23;}
3730 else {sQ=q23; tQ=q12;}
3732 if(c1==c2) {tQ=q12; sQ=q23;}
3733 else {tQ=q23; sQ=q12;}
3735 }else {// key2==key3
3738 // cases not explicity shown below are not possible
3739 if(term==1 || term==5) {sQ=q12; tQ=q13;}
3740 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
3741 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
3742 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
3743 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
3744 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3746 if(c1==c2) {sQ=q12; tQ=q13;}
3747 else {sQ=q13; tQ=q12;}
3749 if(c1==c2) {tQ=q12; sQ=q13;}
3750 else {tQ=q13; sQ=q12;}
3755 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
3756 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
3758 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
3759 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
3761 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
3762 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
3769 //________________________________________________________________________
3770 Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
3774 if(fi==0 || fi==3 || fi==5){// identical masses
3775 qinv = sqrt( pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2));
3776 }else{// different masses
3777 Float_t px = track1[1] + track2[1];
3778 Float_t py = track1[2] + track2[2];
3779 Float_t pz = track1[3] + track2[3];
3780 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
3781 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
3782 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
3784 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
3785 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
3786 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
3787 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
3794 //________________________________________________________________________
3795 void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3797 Float_t p0 = track1[0] + track2[0];
3798 Float_t px = track1[1] + track2[1];
3799 Float_t py = track1[2] + track2[2];
3800 Float_t pz = track1[3] + track2[3];
3802 Float_t mt = sqrt(p0*p0 - pz*pz);
3803 Float_t pt = sqrt(px*px + py*py);
3805 Float_t v0 = track1[0] - track2[0];
3806 Float_t vx = track1[1] - track2[1];
3807 Float_t vy = track1[2] - track2[2];
3808 Float_t vz = track1[3] - track2[3];
3810 qout = (px*vx + py*vy)/pt;
3811 qside = (px*vy - py*vx)/pt;
3812 qlong = (p0*vz - pz*v0)/mt;
3814 //________________________________________________________________________
3815 void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliChaoticity::fKbinsT][AliChaoticity::fCentBins]){
3818 cout<<"LEGO call to SetWeightArrays"<<endl;
3820 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3821 for(Int_t mb=0; mb<fCentBins; mb++){
3822 fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
3823 fNormWeight[tKbin][mb]->SetDirectory(0);
3829 TFile *wFile = new TFile("WeightFile.root","READ");
3830 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3831 else cout<<"Good Weight File Found!"<<endl;
3833 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3834 for(Int_t mb=0; mb<fCentBins; mb++){
3836 TString *name = new TString("Weight_Kt_");
3838 name->Append("_Ky_0");
3839 name->Append("_M_");
3841 name->Append("_ED_0");
3844 fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
3845 fNormWeight[tKbin][mb]->SetDirectory(0);
3854 cout<<"Done reading weight file"<<endl;
3857 //________________________________________________________________________
3858 void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3860 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3862 Float_t qOut=0,qSide=0,qLong=0;
3863 GetQosl(track1, track2, qOut, qSide, qLong);
3865 qSide = fabs(qSide);
3866 qLong = fabs(qLong);
3867 Float_t wd=0, xd=0, yd=0, zd=0;
3868 //Float_t qinv_temp=GetQinv(0,track1, track2);
3871 if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
3872 else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
3874 for(Int_t i=0; i<fKbinsT-1; i++){
3875 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
3878 wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
3881 if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
3882 else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
3884 for(Int_t i=0; i<kQbinsWeights-1; i++){
3885 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
3887 xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
3890 if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
3891 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
3893 for(Int_t i=0; i<kQbinsWeights-1; i++){
3894 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
3896 yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
3899 if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
3900 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
3902 for(Int_t i=0; i<kQbinsWeights-1; i++){
3903 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
3905 zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
3910 //Float_t min = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3911 Float_t minErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3915 deltaW += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - min)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3917 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - min)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3919 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - min)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3921 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - min)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3928 // w interpolation (kt)
3929 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;
3930 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;
3931 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;
3932 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;
3933 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;
3934 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;
3935 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;
3936 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;
3937 // x interpolation (qOut)
3938 Float_t c00 = c000*(1-xd) + c100*xd;
3939 Float_t c10 = c010*(1-xd) + c110*xd;
3940 Float_t c01 = c001*(1-xd) + c101*xd;
3941 Float_t c11 = c011*(1-xd) + c111*xd;
3942 // y interpolation (qSide)
3943 Float_t c0 = c00*(1-yd) + c10*yd;
3944 Float_t c1 = c01*(1-yd) + c11*yd;
3945 // z interpolation (qLong)
3946 wgt = (c0*(1-zd) + c1*zd);
3950 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3951 //Float_t deltaWErr=0;
3954 deltaWErr += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3956 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3958 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - minErr)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3960 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - minErr)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3966 //________________________________________________________________________
3967 Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t k12){
3969 Float_t radius = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3970 Float_t r12=radius*(1-k12/2.0);
3972 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3973 Float_t coulCorr12 = FSICorrelationTherm2(charge1, charge2, qinv);
3974 if(charge1==charge2){
3975 Float_t arg=qinv*r12;
3976 Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3977 EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3978 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
3980 return ((1-myDamp) + myDamp*coulCorr12);
3984 //________________________________________________________________________
3985 Float_t AliChaoticity::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){
3987 Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3988 Float_t radiusSide = radiusOut;
3989 Float_t radiusLong = radiusOut;
3990 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3991 Float_t coulCorr12 = FSICorrelationTherm2(charge1, charge2, qinv);
3992 if(charge1==charge2){
3993 return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
3995 return ((1-myDamp) + myDamp*coulCorr12);
4000 //________________________________________________________________________
4001 Float_t AliChaoticity::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23, Float_t k12, Float_t k13, Float_t k23){
4002 if(term==5) return 1.0;
4004 Float_t radius=fRMax;
4006 Float_t r12=radius*(1-k12/2.0);
4007 Float_t r13=radius*(1-k13/2.0);
4008 Float_t r23=radius*(1-k23/2.0);
4010 //else if(fMbin<=3) {radius = radius-1;}
4011 //else if(fMbin<=5) {radius = radius-2;}
4012 //else {radius = radius-3;}
4015 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
4016 Float_t fc = sqrt(myDamp);
4018 if(SameCharge){// all three of the same charge
4019 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2
4020 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, q13);// K2
4021 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, q23);// K2
4022 Float_t arg12=q12*r12;
4023 Float_t arg13=q13*r13;
4024 Float_t arg23=q23*r23;
4025 Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
4026 EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
4027 Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
4028 EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
4029 Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
4030 EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
4032 Float_t c3QS = 1 + exp(-pow(q12*r12,2))*pow(EW12,2) + exp(-pow(q13*r13,2))*pow(EW13,2) + exp(-pow(q23*r23,2))*pow(EW23,2);
4033 c3QS += 2*exp(-(pow(r12,2)*pow(q12,2) + pow(r13,2)*pow(q13,2) + pow(r23,2)*pow(q23,2))/2.)*EW12*EW13*EW23;
4034 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
4035 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*r12,2))*pow(EW12,2))*coulCorr12;
4036 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*r13,2))*pow(EW13,2))*coulCorr13;
4037 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*r23,2))*pow(EW23,2))*coulCorr23;
4038 w123 += pow(fc,3)*c3QS*coulCorr12*coulCorr13*coulCorr23;
4041 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*r12,2))*pow(EW12,2))*coulCorr12);
4043 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*r13,2))*pow(EW13,2))*coulCorr13);
4045 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*r23,2))*pow(EW23,2))*coulCorr23);
4048 }else{// mixed charge case pair 12 always treated as ss
4049 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2 ss
4050 Float_t coulCorr13 = FSICorrelationTherm2(+1,-1, q13);// K2 os
4051 Float_t coulCorr23 = FSICorrelationTherm2(+1,-1, q23);// K2 os
4052 Float_t arg12=q12*r12;
4053 Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
4054 EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
4056 Float_t c3QS = 1 + exp(-pow(q12*r12,2))*pow(EW12,2);
4057 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
4058 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*r12,2))*pow(EW12,2))*coulCorr12;
4059 w123 += pow(fc,2)*(1-fc)*coulCorr13;
4060 w123 += pow(fc,2)*(1-fc)*coulCorr23;
4061 w123 += pow(fc,3)*c3QS*coulCorr12*coulCorr13*coulCorr23;
4064 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*r12,2))*pow(EW12,2))*coulCorr12);
4066 return ((1-myDamp) + myDamp*coulCorr13);
4068 return ((1-myDamp) + myDamp*coulCorr23);
4073 //________________________________________________________________________
4074 void AliChaoticity::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
4078 cout<<"LEGO call to SetMomResCorrections"<<endl;
4079 fMomResC2 = (TH2D*)temp2D->Clone();
4080 fMomResC2->SetDirectory(0);
4082 TFile *momResFile = new TFile("MomResFile.root","READ");
4083 if(!momResFile->IsOpen()) {
4084 cout<<"No momentum resolution file found"<<endl;
4085 AliFatal("No momentum resolution file found. Kill process.");
4086 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
4088 TH2D *temp2D2 = (TH2D*)momResFile->Get("MomResHisto_pp");
4089 fMomResC2 = (TH2D*)temp2D2->Clone();
4090 fMomResC2->SetDirectory(0);
4092 momResFile->Close();
4095 // fMomResC2->GetBinContent(1,5) should be ~1.007
4096 if(fMomResC2->GetBinContent(1,5) > 1.2) AliFatal("AliChaoticity: SetMomResCorrections Problem");// Additional Safety check
4097 if(fMomResC2->GetBinContent(1,5) < 0.95) AliFatal("AliChaoticity: SetMomResCorrections Problem");// Additional Safety check
4099 for(Int_t bx=1; bx<=fMomResC2->GetNbinsX(); bx++){
4100 for(Int_t by=1; by<=fMomResC2->GetNbinsX(); by++){
4101 if(fMomResC2->GetBinContent(bx,by) > 1.5) fMomResC2->SetBinContent(bx,by, 1.5);// Maximum is ~1.02
4102 if(fMomResC2->GetBinContent(bx,by) < 0.95) fMomResC2->SetBinContent(bx,by, 0.95);// Minimum is ~0.98
4106 cout<<"Done reading momentum resolution file"<<endl;
4108 //________________________________________________________________________
4109 void AliChaoticity::SetFSICorrelations(Bool_t legoCase, TH2D *temp2DTherm[2], TH3D *temp3Dos[6], TH3D *temp3Dss[6]){
4110 // read in 2-particle and 3-particle FSI correlations = K2 & K3
4111 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
4113 cout<<"LEGO call to SetFSICorrelations"<<endl;
4114 fFSI2SS = (TH2D*)temp2DTherm[0]->Clone();
4115 fFSI2OS = (TH2D*)temp2DTherm[1]->Clone();
4117 fFSI2SS->SetDirectory(0);
4118 fFSI2OS->SetDirectory(0);
4120 for(Int_t CB=0; CB<6; CB++) {
4121 fFSIOmega0OS[CB] = (TH3D*)temp3Dos[CB]->Clone();
4122 fFSIOmega0SS[CB] = (TH3D*)temp3Dss[CB]->Clone();
4124 fFSIOmega0OS[CB]->SetDirectory(0);
4125 fFSIOmega0SS[CB]->SetDirectory(0);
4128 cout<<"non LEGO call to SetFSICorrelations"<<endl;
4129 TFile *fsifile = new TFile("KFile.root","READ");
4130 if(!fsifile->IsOpen()) {
4131 cout<<"No FSI file found"<<endl;
4132 AliFatal("No FSI file found. Kill process.");
4133 }else {cout<<"Good FSI File Found!"<<endl;}
4135 TH2D *temphisto2ThermSS = (TH2D*)fsifile->Get("K2ssT");
4136 TH2D *temphisto2ThermOS = (TH2D*)fsifile->Get("K2osT");
4137 TH3D *temphisto3OS[6];
4138 TH3D *temphisto3SS[6];
4139 for(Int_t CB=0; CB<6; CB++) {
4140 TString *nameK3SS = new TString("K3ss_");
4142 temphisto3SS[CB] = (TH3D*)fsifile->Get(nameK3SS->Data());
4144 TString *nameK3OS = new TString("K3os_");
4146 temphisto3OS[CB] = (TH3D*)fsifile->Get(nameK3OS->Data());
4149 fFSI2SS = (TH2D*)temphisto2ThermSS->Clone();
4150 fFSI2OS = (TH2D*)temphisto2ThermOS->Clone();
4151 fFSI2SS->SetDirectory(0);
4152 fFSI2OS->SetDirectory(0);
4154 for(Int_t CB=0; CB<6; CB++) {
4155 fFSIOmega0SS[CB] = (TH3D*)temphisto3SS[CB]->Clone();
4156 fFSIOmega0OS[CB] = (TH3D*)temphisto3OS[CB]->Clone();
4157 fFSIOmega0SS[CB]->SetDirectory(0);
4158 fFSIOmega0OS[CB]->SetDirectory(0);
4165 // fFSI2SS->GetBinContent(1,2) should be ~0.32
4166 if(fFSI2SS->GetBinContent(1,2) > 1.0) AliFatal("AliChaoticity: SetFSICorrelations Problem");// Additional Safety check
4167 if(fFSI2SS->GetBinContent(1,2) < 0.1) AliFatal("AliChaoticity: SetFSICorrelations Problem");// Additional Safety check
4169 for(Int_t ii=1; ii<=fFSI2SS->GetNbinsX(); ii++){
4170 for(Int_t jj=1; jj<=fFSI2SS->GetNbinsY(); jj++){
4171 if(fFSI2SS->GetBinContent(ii,jj) > 1.0) fFSI2SS->SetBinContent(ii,jj, 1.0);
4172 if(fFSI2OS->GetBinContent(ii,jj) > 10.0) fFSI2OS->SetBinContent(ii,jj, 10.0);
4174 if(fFSI2SS->GetBinContent(ii,jj) < 0.05) fFSI2SS->SetBinContent(ii,jj, 0.05);
4175 if(fFSI2OS->GetBinContent(ii,jj) < 0.9) fFSI2OS->SetBinContent(ii,jj, 0.9);
4179 cout<<"Done reading FSI file"<<endl;
4181 //________________________________________________________________________
4182 Float_t AliChaoticity::FSICorrelationTherm2(Int_t charge1, Int_t charge2, Float_t qinv){
4183 // returns 2-particle Coulomb correlations = K2
4184 Int_t qbinL = fFSI2SS->GetYaxis()->FindBin(qinv-fFSI2SS->GetYaxis()->GetBinWidth(1)/2.);
4185 Int_t qbinH = qbinL+1;
4186 if(qbinL <= 0) return 1.0;
4187 if(qbinH > fFSI2SS->GetNbinsY()) return 1.0;
4190 if(charge1==charge2){
4191 slope = fFSI2SS->GetBinContent(fFSIbin+1, qbinL) - fFSI2SS->GetBinContent(fFSIbin+1, qbinH);
4192 slope /= fFSI2SS->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS->GetYaxis()->GetBinCenter(qbinH);
4193 return (slope*(qinv - fFSI2SS->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS->GetBinContent(fFSIbin+1, qbinL));
4195 slope = fFSI2OS->GetBinContent(fFSIbin+1, qbinL) - fFSI2OS->GetBinContent(fFSIbin+1, qbinH);
4196 slope /= fFSI2OS->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS->GetYaxis()->GetBinCenter(qbinH);
4197 return (slope*(qinv - fFSI2OS->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS->GetBinContent(fFSIbin+1, qbinL));
4200 //________________________________________________________________________
4201 Double_t AliChaoticity::FSICorrelationOmega0(Bool_t SameCharge, Double_t Q12, Double_t Q13, Double_t Q23){
4202 // returns 3d 3-particle Coulomb Correlation = K3
4203 Int_t Q12bin = fFSIOmega0SS[fFSIbin]->GetXaxis()->FindBin(Q12);
4204 Int_t Q13bin = fFSIOmega0SS[fFSIbin]->GetZaxis()->FindBin(Q13);
4205 Int_t Q23bin = fFSIOmega0SS[fFSIbin]->GetYaxis()->FindBin(Q23);
4206 Int_t index12L = int(fabs(Q12 - fFSI2SS->GetYaxis()->GetBinWidth(1)/2.)/(fFSI2SS->GetYaxis()->GetBinWidth(1)));
4207 Int_t index12H = index12L+1;
4208 Int_t index13L = int(fabs(Q13 - fFSI2SS->GetYaxis()->GetBinWidth(1)/2.)/(fFSI2SS->GetYaxis()->GetBinWidth(1)));
4209 Int_t index13H = index13L+1;
4210 Int_t index23L = int(fabs(Q23 - fFSI2SS->GetYaxis()->GetBinWidth(1)/2.)/(fFSI2SS->GetYaxis()->GetBinWidth(1)));
4211 Int_t index23H = index23L+1;
4214 if(fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
4215 Double_t base = fFSIOmega0SS[fFSIbin]->GetBinContent(index12L+1, index23L+1, index13L+1);
4216 Double_t InterPolated = 0;
4217 Double_t slope12 = fFSIOmega0SS[fFSIbin]->GetBinContent(index12H+1, index23L+1, index13L+1);
4219 slope12 /= fFSIOmega0SS[fFSIbin]->GetXaxis()->GetBinWidth(1);
4220 InterPolated += slope12*fabs(Q12 - fFSIOmega0SS[fFSIbin]->GetXaxis()->GetBinCenter(index12L+1));
4221 Double_t slope23 = fFSIOmega0SS[fFSIbin]->GetBinContent(index12L+1, index23H+1, index13L+1);
4223 slope23 /= fFSIOmega0SS[fFSIbin]->GetYaxis()->GetBinWidth(1);
4224 InterPolated += slope23*fabs(Q23 - fFSIOmega0SS[fFSIbin]->GetYaxis()->GetBinCenter(index23L+1));
4225 Double_t slope13 = fFSIOmega0SS[fFSIbin]->GetBinContent(index12L+1, index23L+1, index13H+1);
4227 slope13 /= fFSIOmega0SS[fFSIbin]->GetZaxis()->GetBinWidth(1);
4228 InterPolated += slope13*fabs(Q13 - fFSIOmega0SS[fFSIbin]->GetZaxis()->GetBinCenter(index13L+1));
4229 if( (base+InterPolated) <= 0) return 1.0;
4230 return (base+InterPolated);
4232 }else{// mixed charge. Q12 is always designated as the same-charge pair
4233 if(fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
4234 Double_t base = fFSIOmega0OS[fFSIbin]->GetBinContent(index12L+1, index23H+1, index13H+1);
4235 Double_t InterPolated = 0;
4236 Double_t slope12 = fFSIOmega0OS[fFSIbin]->GetBinContent(index12H+1, index23H+1, index13H+1);
4238 slope12 /= fFSIOmega0OS[fFSIbin]->GetXaxis()->GetBinWidth(1);
4239 InterPolated += slope12*fabs(Q12 - fFSIOmega0OS[fFSIbin]->GetXaxis()->GetBinCenter(index12L+1));
4240 Double_t slope23 = fFSIOmega0OS[fFSIbin]->GetBinContent(index12L+1, index23L+1, index13H+1);
4242 slope23 /= fFSIOmega0OS[fFSIbin]->GetYaxis()->GetBinWidth(1);
4243 InterPolated += slope23*fabs(Q23 - fFSIOmega0OS[fFSIbin]->GetYaxis()->GetBinCenter(index23L+1));
4244 Double_t slope13 = fFSIOmega0OS[fFSIbin]->GetBinContent(index12L+1, index23H+1, index13L+1);
4246 slope13 /= fFSIOmega0OS[fFSIbin]->GetZaxis()->GetBinWidth(1);
4247 InterPolated += slope13*fabs(Q13 - fFSIOmega0OS[fFSIbin]->GetZaxis()->GetBinCenter(index13L+1));
4248 if( (base+InterPolated) <= 0) return 1.0;
4249 return (base+InterPolated);
4253 //________________________________________________________________________
4254 void AliChaoticity::FourVectProdTerms(Float_t pV1[], Float_t pV2[], Float_t pV3[], Float_t& QS1v1, Float_t& QS2, Float_t& QS3v1, Float_t& QS1v2, Float_t& QS3v2){
4255 QS1v1 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
4256 QS1v1 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
4257 QS1v1 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
4258 QS2 = (pV1[1]-pV2[1])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[1]-pV3[1]);
4259 QS3v1 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
4260 QS3v1 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
4262 QS1v2 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
4263 QS1v2 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
4264 QS3v2 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
4265 QS3v2 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
4266 QS3v2 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
4268 //________________________________________________________________________