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
38 // Author: Dhevan Gangadharan
40 ClassImp(AliChaoticity)
42 //________________________________________________________________________
43 AliChaoticity::AliChaoticity():
58 fPdensityExplicitLoop(kFALSE),
59 fPdensityPairCut(kTRUE),
60 fTabulatePairs(kFALSE),
78 fQupperBoundWeights(0),
97 fMinSepTPCEntrancePhi(0),
98 fMinSepTPCEntranceEta(0),
127 fOtherPairLocation1(),
128 fOtherPairLocation2(),
133 // Default constructor
134 for(Int_t mb=0; mb<fMbins; mb++){
135 for(Int_t edB=0; edB<kEDbins; edB++){
136 for(Int_t c1=0; c1<2; c1++){
137 for(Int_t c2=0; c2<2; c2++){
138 for(Int_t sc=0; sc<kSCLimit2; sc++){
139 for(Int_t term=0; term<2; term++){
141 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
143 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
144 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
145 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
146 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
147 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
148 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
153 for(Int_t c3=0; c3<2; c3++){
154 for(Int_t sc=0; sc<kSCLimit3; sc++){
155 for(Int_t term=0; term<5; term++){
157 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
158 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
159 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
160 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
161 for(Int_t dt=0; dt<kDENtypes; dt++){
162 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
170 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
171 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
172 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
173 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
181 //________________________________________________________________________
182 AliChaoticity::AliChaoticity(const Char_t *name, Bool_t MCdecision, Bool_t Tabulatedecision, Bool_t PbPbdecision, Int_t lowCentBin, Int_t highCentBin, Bool_t lego)
183 : AliAnalysisTaskSE(name),
196 fPbPbcase(PbPbdecision),
197 fPdensityExplicitLoop(kFALSE),
198 fPdensityPairCut(kTRUE),
199 fTabulatePairs(Tabulatedecision),
205 fCentBinLowLimit(lowCentBin),
206 fCentBinHighLimit(highCentBin),
217 fQupperBoundWeights(0),
236 fMinSepTPCEntrancePhi(0),
237 fMinSepTPCEntranceEta(0),
266 fOtherPairLocation1(),
267 fOtherPairLocation2(),
276 fTabulatePairs=Tabulatedecision;
277 fPbPbcase=PbPbdecision;
278 fPdensityExplicitLoop = kFALSE;
279 fPdensityPairCut = kTRUE;
280 fCentBinLowLimit = lowCentBin;
281 fCentBinHighLimit = highCentBin;
283 for(Int_t mb=0; mb<fMbins; mb++){
284 for(Int_t edB=0; edB<kEDbins; edB++){
285 for(Int_t c1=0; c1<2; c1++){
286 for(Int_t c2=0; c2<2; c2++){
287 for(Int_t sc=0; sc<kSCLimit2; sc++){
288 for(Int_t term=0; term<2; term++){
290 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
292 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
293 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
294 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
295 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
296 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
297 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
302 for(Int_t c3=0; c3<2; c3++){
303 for(Int_t sc=0; sc<kSCLimit3; sc++){
304 for(Int_t term=0; term<5; term++){
306 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
307 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
308 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
309 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
310 for(Int_t dt=0; dt<kDENtypes; dt++){
311 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
319 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
320 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
321 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
322 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
330 DefineOutput(1, TList::Class());
332 //________________________________________________________________________
333 AliChaoticity::AliChaoticity(const AliChaoticity &obj)
334 : AliAnalysisTaskSE(obj.fname),
338 fOutputList(obj.fOutputList),
339 fPIDResponse(obj.fPIDResponse),
342 fTempStruct(obj.fTempStruct),
343 fRandomNumber(obj.fRandomNumber),
345 fMCcase(obj.fMCcase),
346 fAODcase(obj.fAODcase),
347 fPbPbcase(obj.fPbPbcase),
348 fPdensityExplicitLoop(obj.fPdensityExplicitLoop),
349 fPdensityPairCut(obj.fPdensityPairCut),
350 fTabulatePairs(obj.fTabulatePairs),
351 fBfield(obj.fBfield),
355 fMultLimit(obj.fMultLimit),
356 fCentBinLowLimit(obj.fCentBinLowLimit),
357 fCentBinHighLimit(obj.fCentBinHighLimit),
358 fEventCounter(obj.fEventCounter),
359 fEventsToMix(obj.fEventsToMix),
360 fZvertexBins(obj.fZvertexBins),
363 fQLowerCut(obj.fQLowerCut),
366 fKupperBound(obj.fKupperBound),
367 fQupperBound(obj.fQupperBound),
368 fQupperBoundWeights(obj.fQupperBoundWeights),
377 fDampStart(obj.fDampStart),
378 fDampStep(obj.fDampStep),
383 fTPCTOFboundry(obj.fTPCTOFboundry),
384 fTOFboundry(obj.fTOFboundry),
385 fSigmaCutTPC(obj.fSigmaCutTPC),
386 fSigmaCutTOF(obj.fSigmaCutTOF),
387 fMinSepTPCEntrancePhi(obj.fMinSepTPCEntrancePhi),
388 fMinSepTPCEntranceEta(obj.fMinSepTPCEntranceEta),
389 fShareQuality(obj.fShareQuality),
390 fShareFraction(obj.fShareFraction),
391 fTrueMassP(obj.fTrueMassP),
392 fTrueMassPi(obj.fTrueMassPi),
393 fTrueMassK(obj.fTrueMassK),
394 fTrueMassKs(obj.fTrueMassKs),
395 fTrueMassLam(obj.fTrueMassLam),
396 fKtbinL(obj.fKtbinL),
397 fKtbinH(obj.fKtbinH),
398 fKybinL(obj.fKybinL),
399 fKybinH(obj.fKybinH),
400 fQobinL(obj.fQobinL),
401 fQobinH(obj.fQobinH),
402 fQsbinL(obj.fQsbinL),
403 fQsbinH(obj.fQsbinH),
404 fQlbinL(obj.fQlbinL),
405 fQlbinH(obj.fQlbinH),
406 fDummyB(obj.fDummyB),
417 fOtherPairLocation1(),
418 fOtherPairLocation2(),
426 //________________________________________________________________________
427 AliChaoticity &AliChaoticity::operator=(const AliChaoticity &obj)
429 // Assignment operator
436 fOutputList = obj.fOutputList;
437 fPIDResponse = obj.fPIDResponse;
440 fTempStruct = obj.fTempStruct;
441 fRandomNumber = obj.fRandomNumber;
443 fMCcase = obj.fMCcase;
444 fAODcase = obj.fAODcase;
445 fPbPbcase = obj.fPbPbcase;
446 fPdensityExplicitLoop = obj.fPdensityExplicitLoop;
447 fPdensityPairCut = obj.fPdensityPairCut;
448 fTabulatePairs = obj.fTabulatePairs;
449 fBfield = obj.fBfield;
453 fMultLimit = obj.fMultLimit;
454 fCentBinLowLimit = obj.fCentBinLowLimit;
455 fCentBinHighLimit = obj.fCentBinHighLimit;
456 fEventCounter = obj.fEventCounter;
457 fEventsToMix = obj.fEventsToMix;
458 fZvertexBins = obj.fZvertexBins;
461 fQLowerCut = obj.fQLowerCut;
464 fKupperBound = obj.fKupperBound;
465 fQupperBound = obj.fQupperBound;
466 fQupperBoundWeights = obj.fQupperBoundWeights;
475 fDampStart = obj.fDampStart;
476 fDampStep = obj.fDampStep;
481 fTPCTOFboundry = obj.fTPCTOFboundry;
482 fTOFboundry = obj.fTOFboundry;
483 fSigmaCutTPC = obj.fSigmaCutTPC;
484 fSigmaCutTOF = obj.fSigmaCutTOF;
485 fMinSepTPCEntrancePhi = obj.fMinSepTPCEntrancePhi;
486 fMinSepTPCEntranceEta = obj.fMinSepTPCEntranceEta;
487 fShareQuality = obj.fShareQuality;
488 fShareFraction = obj.fShareFraction;
489 fTrueMassP = obj.fTrueMassP;
490 fTrueMassPi = obj.fTrueMassPi;
491 fTrueMassK = obj.fTrueMassK;
492 fTrueMassKs = obj.fTrueMassKs;
493 fTrueMassLam = obj.fTrueMassLam;
494 fKtbinL = obj.fKtbinL;
495 fKtbinH = obj.fKtbinH;
496 fKybinL = obj.fKybinL;
497 fKybinH = obj.fKybinH;
498 fQobinL = obj.fQobinL;
499 fQobinH = obj.fQobinH;
500 fQsbinL = obj.fQsbinL;
501 fQsbinH = obj.fQsbinH;
502 fQlbinL = obj.fQlbinL;
503 fQlbinH = obj.fQlbinH;
504 fDummyB = obj.fDummyB;
505 fNormWeight = obj.fNormWeight;
506 fNormWeightErr = obj.fNormWeightErr;
510 //________________________________________________________________________
511 AliChaoticity::~AliChaoticity()
514 if(fAOD) delete fAOD;
515 if(fESD) delete fESD;
516 if(fOutputList) delete fOutputList;
517 if(fPIDResponse) delete fPIDResponse;
519 if(fEvt) delete fEvt;
520 if(fTempStruct) delete [] fTempStruct;
521 if(fRandomNumber) delete fRandomNumber;
523 for(int i=0; i<fMultLimit; i++){
524 if(fPairLocationSE[i]) delete [] fPairLocationSE[i];
525 if(fPairLocationME[i]) delete [] fPairLocationME[i];
526 for(int j=0; j<2; j++){
527 if(fOtherPairLocation1[j][i]) delete [] fOtherPairLocation1[j][i];
528 if(fOtherPairLocation2[j][i]) delete [] fOtherPairLocation2[j][i];
530 for(int j=0; j<3; j++) if(fNormPairSwitch[j][i]) delete [] fNormPairSwitch[j][i];
531 for(int j=0; j<4; j++) if(fPairSplitCut[j][i]) delete [] fPairSplitCut[j][i];
533 for(int i=0; i<kPairLimit; i++) if(fTripletSkip1[i]) delete [] fTripletSkip1[i];
534 for(int i=0; i<2*kPairLimit; i++) if(fTripletSkip2[i]) delete [] fTripletSkip2[i];
535 for(int i=0; i<3; i++) if(fNormPairs[i]) delete [] fNormPairs[i];
537 for(Int_t mb=0; mb<fMbins; mb++){
538 for(Int_t edB=0; edB<kEDbins; edB++){
539 for(Int_t c1=0; c1<2; c1++){
540 for(Int_t c2=0; c2<2; c2++){
541 for(Int_t sc=0; sc<kSCLimit2; sc++){
542 for(Int_t term=0; term<2; term++){
544 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;
546 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;
547 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;
548 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;
549 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;
550 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;
551 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;
556 for(Int_t c3=0; c3<2; c3++){
557 for(Int_t sc=0; sc<kSCLimit3; sc++){
558 for(Int_t term=0; term<5; term++){
560 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;
561 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;
562 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;
563 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;
564 for(Int_t dt=0; dt<kDENtypes; dt++){
565 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;
573 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
574 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
575 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD;
576 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD;
585 //________________________________________________________________________
586 void AliChaoticity::ParInit()
588 cout<<"AliChaoticity MyInit() call"<<endl;
590 fRandomNumber = new TRandom3();
591 fRandomNumber->SetSeed(0);
596 if(fPdensityExplicitLoop) fEventsToMix=3;
597 else if(fPdensityPairCut && !fPdensityExplicitLoop) fEventsToMix=2;
601 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
602 fTOFboundry = 2.1;// TOF pid used below this momentum
603 fSigmaCutTPC = 2.0;// 2.0 (norm). 1.5 Sys Deviation
604 fSigmaCutTOF = 2.0;// 2.0 (norm). 1.5 Sys Deviation
606 ////////////////////////////////////////////////
608 fMinSepTPCEntrancePhi = 0.035;// 0.035(norm). 0.04 Sys Deviation
609 fMinSepTPCEntranceEta = 0.035;// 0.035(norm). 0.04 Sys Deviation
610 fShareQuality = .5;// max
611 fShareFraction = .05;// max
612 ////////////////////////////////////////////////
615 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
616 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
620 if(fPbPbcase) {// PbPb
621 fMultLimit=kMultLimitPbPb;
626 fNormQcutLow[0] = 0.15;//1.06 (test also at 0.15)
627 fNormQcutHigh[0] = 0.175;//1.1 (test also at 0.175)
628 fNormQcutLow[1] = 1.34;//1.34
629 fNormQcutHigh[1] = 1.4;//1.4
630 fNormQcutLow[2] = 1.1;//1.1
631 fNormQcutHigh[2] = 1.4;//1.4
634 fMultLimit=kMultLimitpp;
639 fNormQcutLow[0] = 1.0;
640 fNormQcutHigh[0] = 1.5;
641 fNormQcutLow[1] = 1.0;
642 fNormQcutHigh[1] = 1.5;
643 fNormQcutLow[2] = 1.0;
644 fNormQcutHigh[2] = 1.5;
647 fQLowerCut = 0.002;// was 0.002
650 // 4x1 (Kt: 0-0.2, 0.2-0.4, 0.4-0.6, 0.6-1.0)
651 //fKstepT[0] = 0.2; fKstepT[1] = 0.2; fKstepT[2] = 0.2; fKstepT[3] = 0.4;
652 //fKstepY[0] = 1.6;// -0.8 to +0.8
653 //fKmeanT[0] = 0.1; fKmeanT[1] = 0.3; fKmeanT[2] = 0.5; fKmeanT[3] = 0.8;
654 //fKmeanY[0] = 0;// central y
655 //fKmiddleT[0] = 0.1; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.5; fKmiddleT[3] = 0.8;
657 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
658 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
659 fKstepY[0] = 1.6;// -0.8 to +0.8
660 fKmeanT[0] = 0.241; fKmeanT[1] = 0.369; fKmeanT[2] = 0.573;
661 fKmeanY[0] = 0;// central y
662 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
664 // 2x1 (Kt: 0-0.4, 0.4-1.0)
665 //fKstepT[0] = 0.4; fKstepT[1] = 0.6;
666 //fKstepY[0] = 1.6;// -0.8 to +0.8
667 //fKmeanT[0] = 0.255; fKmeanT[1] = 0.505;
668 //fKmiddleT[0] = 0.2; fKmiddleT[1] = 0.7;
669 //fKmeanY[0] = 0;// central y
673 //fKstepY[0] = 1.6;// -0.8 to +0.8
674 //fKmeanT[0] = 0.306;
675 //fKmiddleT[0] = 0.5;
676 //fKmeanY[0] = 0;// central y
679 fQupperBoundWeights = 0.2;
681 fQstep = fQupperBoundWeights/Float_t(kQbinsWeights);
682 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstep;}
689 fEC = new AliChaoticityEventCollection **[fZvertexBins];
690 for(UShort_t i=0; i<fZvertexBins; i++){
692 fEC[i] = new AliChaoticityEventCollection *[fMbins];
694 for(UShort_t j=0; j<fMbins; j++){
696 fEC[i][j] = new AliChaoticityEventCollection(fEventsToMix+1, fMultLimit, kPairLimit, fMCcase);
701 for(Int_t i=0; i<fMultLimit; i++) fDefaultsCharMult[i]='0';
702 for(Int_t i=0; i<kPairLimit; i++) fDefaultsCharSE[i]='0';
703 for(Int_t i=0; i<2*kPairLimit; i++) fDefaultsCharME[i]='0';
704 for(Int_t i=0; i<fMultLimit; i++) fDefaultsInt[i]=-1;
705 for(Int_t i=0; i<fMultLimit; i++) fPairLocationSE[i] = new TArrayI(fMultLimit,fDefaultsInt);
706 for(Int_t i=0; i<fMultLimit; i++) fPairLocationME[i] = new TArrayI(fMultLimit,fDefaultsInt);
707 for(Int_t i=0; i<kPairLimit; i++) fTripletSkip1[i] = new TArrayC(fMultLimit,fDefaultsCharSE);
708 for(Int_t i=0; i<2*kPairLimit; i++) fTripletSkip2[i] = new TArrayC(fMultLimit,fDefaultsCharME);
711 // Normalization utilities
712 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
713 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
714 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
715 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
716 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
717 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
718 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
720 // Track Merging/Splitting utilities
721 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P11
722 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P12
723 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P13
724 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[3][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P23
727 fNormPairs[0] = new AliChaoticityNormPairStruct[kNormPairLimit];
728 fNormPairs[1] = new AliChaoticityNormPairStruct[kNormPairLimit];
731 fTempStruct = new AliChaoticityTrackStruct[fMultLimit];
734 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
737 // Initialize Weight Array
738 fNormWeight = new Float_t******[fMbins];// fMbin
739 fNormWeightErr = new Float_t******[fMbins];// fMbin
740 for(Int_t i=0; i<fMbins; i++){// Mbin iterator
741 fNormWeight[i] = new Float_t*****[kEDbins];// create ED bins
742 fNormWeightErr[i] = new Float_t*****[kEDbins];// create ED bins
744 for(Int_t j=0; j<kEDbins; j++){// ED iterator
745 fNormWeight[i][j] = new Float_t****[kKbinsT];// create Kt bins
746 fNormWeightErr[i][j] = new Float_t****[kKbinsT];// create Kt bins
748 for(Int_t k=0; k<kKbinsT; k++){// Kt iterator
749 fNormWeight[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
750 fNormWeightErr[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
752 for(Int_t l=0; l<kKbinsY; l++){// Ky iterator
753 fNormWeight[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
754 fNormWeightErr[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
756 for(Int_t m=0; m<kQbinsWeights; m++){// Qout iterator
757 fNormWeight[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
758 fNormWeightErr[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
760 for(Int_t n=0; n<kQbinsWeights; n++){// Qside iterator
761 fNormWeight[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
762 fNormWeightErr[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
772 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
773 if(!fLEGO && !fTabulatePairs) {
774 SetWeightArrays();// Set Weight Array
775 SetCoulCorrections();// Read in 2-particle coulomb corrections
776 if(!fMCcase) SetMomResCorrections();// Read Momentum resolution file
780 //________________________________________________________________________
781 void AliChaoticity::UserCreateOutputObjects()
786 ParInit();// Initialize my settings
789 fOutputList = new TList();
790 fOutputList->SetOwner();
792 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
793 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
794 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
795 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
796 fOutputList->Add(fVertexDist);
799 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
800 fOutputList->Add(fDCAxyDistPlus);
801 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
802 fOutputList->Add(fDCAzDistPlus);
803 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
804 fOutputList->Add(fDCAxyDistMinus);
805 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
806 fOutputList->Add(fDCAzDistMinus);
809 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
810 fOutputList->Add(fEvents1);
811 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
812 fOutputList->Add(fEvents2);
814 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
815 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
816 fOutputList->Add(fMultDist1);
818 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
819 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
820 fOutputList->Add(fMultDist2);
822 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
823 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
824 fOutputList->Add(fMultDist3);
826 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
827 fOutputList->Add(fPtEtaDist);
829 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
830 fOutputList->Add(fPhiPtDist);
832 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
833 fOutputList->Add(fTOFResponse);
834 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
835 fOutputList->Add(fTPCResponse);
837 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
838 fOutputList->Add(fRejectedPairs);
839 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
840 fOutputList->Add(fRejectedEvents);
842 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
843 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
844 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
845 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
847 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
848 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
849 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
850 if(fMCcase) fOutputList->Add(fAllOSPairs);
852 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
853 fOutputList->Add(fAvgMult);
855 TH3D *fTPNRejects = new TH3D("fTPNRejects","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
856 fOutputList->Add(fTPNRejects);
858 TH2D *fKt3Dist = new TH2D("fKt3Dist","",fMbins,.5,fMbins+.5, 10,0,1);
859 fOutputList->Add(fKt3Dist);
863 if(fPdensityExplicitLoop || fPdensityPairCut){
865 for(Int_t mb=0; mb<fMbins; mb++){
866 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
868 for(Int_t edB=0; edB<kEDbins; edB++){
869 for(Int_t c1=0; c1<2; c1++){
870 for(Int_t c2=0; c2<2; c2++){
871 for(Int_t sc=0; sc<kSCLimit2; sc++){
872 for(Int_t term=0; term<2; term++){
874 TString *nameEx2 = new TString("Explicit2_Charge1_");
876 nameEx2->Append("_Charge2_");
878 nameEx2->Append("_SC_");
880 nameEx2->Append("_M_");
882 nameEx2->Append("_ED_");
884 nameEx2->Append("_Term_");
887 if(sc==0 || sc==3 || sc==5){
888 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
891 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);
892 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
893 // Momentum resolution histos
894 if(fMCcase && sc==0){
895 TString *nameIdeal = new TString(nameEx2->Data());
896 nameIdeal->Append("_Ideal");
897 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",kDENtypes,-0.5,kDENtypes-0.5, kQbinsWeights,0,fQupperBoundWeights);
898 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
899 TString *nameSmeared = new TString(nameEx2->Data());
900 nameSmeared->Append("_Smeared");
901 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",kDENtypes,-0.5,kDENtypes-0.5, kQbinsWeights,0,fQupperBoundWeights);
902 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
906 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
907 nameEx2OSLB1->Append("_osl_b1");
908 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
909 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL);
910 nameEx2OSLB1->Append("_QW");
911 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
912 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW);
914 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
915 nameEx2OSLB2->Append("_osl_b2");
916 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
917 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL);
918 nameEx2OSLB2->Append("_QW");
919 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",50,0,.2, 50,0,.2, 50,0,.2);
920 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW);
927 // skip 3-particle if Tabulate6DPairs is true
928 if(fTabulatePairs) continue;
930 for(Int_t c3=0; c3<2; c3++){
931 for(Int_t sc=0; sc<kSCLimit3; sc++){
932 for(Int_t term=0; term<5; term++){
933 TString *nameEx3 = new TString("Explicit3_Charge1_");
935 nameEx3->Append("_Charge2_");
937 nameEx3->Append("_Charge3_");
939 nameEx3->Append("_SC_");
941 nameEx3->Append("_M_");
943 nameEx3->Append("_ED_");
945 nameEx3->Append("_Term_");
948 TString *namePC3 = new TString("PairCut3_Charge1_");
950 namePC3->Append("_Charge2_");
952 namePC3->Append("_Charge3_");
954 namePC3->Append("_SC_");
956 namePC3->Append("_M_");
958 namePC3->Append("_ED_");
960 namePC3->Append("_Term_");
963 ///////////////////////////////////////
964 // skip degenerate histograms
965 if(sc==0 || sc==6 || sc==9){// Identical species
966 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
967 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
969 if( (c1+c2)==1) {if(c1!=0) continue;}
970 }else {}// do nothing for pi-k-p case
972 /////////////////////////////////////////
974 if(fPdensityExplicitLoop){
975 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);
976 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3);
978 nameEx3->Append("_Norm");
979 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);
980 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3);
982 if(fPdensityPairCut){
983 TString *nameNorm=new TString(namePC3->Data());
984 nameNorm->Append("_Norm");
985 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
986 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
989 TString *name3DQ=new TString(namePC3->Data());
990 name3DQ->Append("_3D");
991 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);
992 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
995 if(c1==c2 && c1==c3 && term==4 && sc==0 && fMCcase==kFALSE){
996 for(Int_t dt=0; dt<kDENtypes; dt++){
997 TString *nameDenType=new TString("PairCut3_Charge1_");
999 nameDenType->Append("_Charge2_");
1001 nameDenType->Append("_Charge3_");
1003 nameDenType->Append("_SC_");
1005 nameDenType->Append("_M_");
1007 nameDenType->Append("_ED_");
1008 *nameDenType += edB;
1009 nameDenType->Append("_TPN_");
1012 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);
1013 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm);
1014 // neglect errors for TPN
1015 //nameDenType->Append("_Err");
1016 //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);
1017 //fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr);
1022 }// c and sc exclusion
1036 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
1037 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
1038 for(Int_t mb=0; mb<fMbins; mb++){
1039 for(Int_t edB=0; edB<kEDbins; edB++){
1041 TString *nameNum = new TString("TwoPart_num_Kt_");
1043 nameNum->Append("_Ky_");
1045 nameNum->Append("_M_");
1047 nameNum->Append("_ED_");
1050 TString *nameDen = new TString("TwoPart_den_Kt_");
1052 nameDen->Append("_Ky_");
1054 nameDen->Append("_M_");
1056 nameDen->Append("_ED_");
1060 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = new TH3I(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1061 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD);
1063 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = new TH3I(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1064 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD);
1073 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1074 fOutputList->Add(fQsmearMean);
1075 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1076 fOutputList->Add(fQsmearSq);
1077 TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
1078 fOutputList->Add(fQDist);
1081 ////////////////////////////////////
1082 ///////////////////////////////////
1084 PostData(1, fOutputList);
1087 //________________________________________________________________________
1088 void AliChaoticity::Exec(Option_t *)
1091 // Called for each event
1092 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1095 if(fAODcase) fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1096 else fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1098 if(fAODcase) {if (!fAOD) {Printf("ERROR: fAOD not available"); return;}}
1099 else {if (!fESD) {Printf("ERROR: fESD not available"); return;}}
1101 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1104 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1105 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1106 if(!isSelected1 && !fMCcase) {/*cout<<"Event Rejected"<<endl;*/ return;}
1107 }if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1108 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1109 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1110 if(!isSelected1 && !isSelected2 && !fMCcase) {/*cout<<"Event Rejected"<<endl;*/ return;}
1111 }else {/*cout<<"Event Rejected"<<endl;*/ return;}
1113 ///////////////////////////////////////////////////////////
1114 const AliAODVertex *primaryVertexAOD;
1115 AliCentrality *centrality;// for AODs and ESDs
1118 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1119 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1120 fPIDResponse = inputHandler->GetPIDResponse();
1123 TClonesArray *mcArray = 0x0;
1126 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1127 if(!mcArray || mcArray->GetEntriesFast() >= 110000){
1128 cout<<"No MC particle branch found or Array too large!!"<<endl;
1129 cout<<mcArray->GetEntriesFast()<<endl;
1137 Int_t positiveTracks=0, negativeTracks=0;
1138 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1140 Double_t vertex[3]={0};
1142 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1143 /////////////////////////////////////////////////
1146 Float_t centralityPercentile=0;
1147 Float_t cStep=5.0, cStart=0;
1150 if(fAODcase){// AOD case
1153 centrality = fAOD->GetCentrality();
1154 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1155 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1156 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {cout<<"Centrality out of Range. Skipping Event"<<endl; return;}
1157 cout<<"Centrality % = "<<centralityPercentile<<endl;
1163 ////////////////////////////////
1165 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1166 primaryVertexAOD = fAOD->GetPrimaryVertex();
1167 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1169 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1170 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1172 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1173 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1175 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1177 fBfield = fAOD->GetMagneticField();
1179 for(Int_t i=0; i<fZvertexBins; i++){
1180 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1188 /////////////////////////////
1189 // Create Shuffled index list
1190 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1191 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1192 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1193 /////////////////////////////
1196 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1197 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1198 if (!aodtrack) continue;
1199 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1201 status=aodtrack->GetStatus();
1203 if((Bool_t)(((1<<7) & aodtrack->GetFilterMap()) == 0)) continue;// AOD filterBit cut
1204 // 1<<5 is for "standard cuts with tight dca cut"
1205 // 1<<7 is for TPC only tracking
1206 if(aodtrack->Pt() < 0.16) continue;
1207 if(fabs(aodtrack->Eta()) > 0.8) continue;
1210 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1211 if(!goodMomentum) continue;
1212 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1217 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1218 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1219 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1221 fTempStruct[myTracks].fStatus = status;
1222 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1223 fTempStruct[myTracks].fId = aodtrack->GetID();
1224 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1225 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1226 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1227 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1228 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1229 fTempStruct[myTracks].fEta = aodtrack->Eta();
1230 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1231 fTempStruct[myTracks].fDCAXY = dca2[0];
1232 fTempStruct[myTracks].fDCAZ = dca2[1];
1233 fTempStruct[myTracks].fDCA = dca3d;
1234 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1235 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1239 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1240 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1241 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1243 if(fTempStruct[myTracks].fCharge==+1) {
1244 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1245 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1247 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1248 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1251 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1252 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1255 // nSimga PID workaround
1256 fTempStruct[myTracks].fElectron = kFALSE;
1257 fTempStruct[myTracks].fPion = kFALSE;
1258 fTempStruct[myTracks].fKaon = kFALSE;
1259 fTempStruct[myTracks].fProton = kFALSE;
1261 Float_t nSigmaTPC[5];
1262 Float_t nSigmaTOF[5];
1263 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1264 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1265 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1266 Float_t signalTPC=0, signalTOF=0;
1267 Double_t integratedTimesTOF[10]={0};
1268 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1269 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1270 if (!aodTrack2) continue;
1271 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1273 UInt_t status2=aodTrack2->GetStatus();
1275 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1276 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1277 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1278 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1279 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1281 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1282 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1283 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1284 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1285 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1286 signalTPC = aodTrack2->GetTPCsignal();
1288 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1289 fTempStruct[myTracks].fTOFhit = kTRUE;
1290 signalTOF = aodTrack2->GetTOFsignal();
1291 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1292 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1297 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1298 if(fTempStruct[myTracks].fTOFhit) {
1299 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1303 // Use TOF if good hit and above threshold
1304 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1305 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1306 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1307 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1308 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1309 }else {// TPC info instead
1310 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1311 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1312 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1313 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1317 //////////////////////////////////////
1318 // Bayesian PIDs for TPC only tracking
1319 //const Double_t* PID = aodtrack->PID();
1320 //fTempStruct[myTracks].fElectron = kFALSE;
1321 //fTempStruct[myTracks].Pion = kFALSE;
1322 //fTempStruct[myTracks].Kaon = kFALSE;
1323 //fTempStruct[myTracks].Proton = kFALSE;
1326 //if(PID[2] > 0.2) fTempStruct[myTracks].Pion = kTRUE;// pions: 0.2 --> 0.5
1329 //if(PID[3] <= 0.5) fTempStruct[myTracks].Kaon = kFALSE;// kaons
1332 //if(PID[4] <= 0.5) fTempStruct[myTracks].Proton = kFALSE;// protons
1333 //////////////////////////////////////
1336 // Ensure there is only 1 candidate per track
1337 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1338 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1339 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1340 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1341 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1342 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1343 ////////////////////////
1344 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1346 if(!fTempStruct[myTracks].fPion) continue;// only pions
1351 if(fTempStruct[myTracks].fPion) {// pions
1352 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1353 fTempStruct[myTracks].fKey = 1;
1354 }else if(fTempStruct[myTracks].fKaon){// kaons
1355 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1356 fTempStruct[myTracks].fKey = 10;
1358 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1359 fTempStruct[myTracks].fKey = 100;
1364 if(aodtrack->Charge() > 0) positiveTracks++;
1365 else negativeTracks++;
1367 if(fTempStruct[myTracks].fPion) pionCount++;
1368 if(fTempStruct[myTracks].fKaon) kaonCount++;
1369 if(fTempStruct[myTracks].fProton) protonCount++;
1374 }else {// ESD tracks
1375 cout<<"ESDs not supported currently"<<endl;
1381 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1385 cout<<"There are "<<myTracks<<" myTracks"<<endl;
1386 cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
1388 /////////////////////////////////////////
1389 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1390 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1391 /////////////////////////////////////////
1394 ////////////////////////////////
1395 ///////////////////////////////
1396 // Mbin determination
1398 // Mbin set to Pion Count Only for pp!!!!!!!
1401 for(Int_t i=0; i<kMultBinspp; i++){
1402 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1403 // Mbin 0 has 1 pion
1406 for(Int_t i=0; i<kCentBins; i++){
1407 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1408 fMbin=i;// 0 = most central
1414 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1416 //////////////////////////////////////////////////
1417 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1418 //////////////////////////////////////////////////
1422 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1423 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1425 ////////////////////////////////////
1426 // Add event to buffer if > 0 tracks
1428 fEC[zbin][fMbin]->FIFOShift();
1429 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1430 (fEvt)->fNtracks = myTracks;
1431 (fEvt)->fFillStatus = 1;
1432 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1434 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1435 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1436 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1437 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1438 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1439 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1440 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1447 Float_t qinv12=0, qinv13=0, qinv23=0;
1448 Float_t qout=0, qside=0, qlong=0;
1449 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1450 Float_t transK12=0, rapK12=0, transK3=0;
1451 Int_t transKbin=0, rapKbin=0;
1453 Int_t ch1=0, ch2=0, ch3=0;
1454 Short_t key1=0, key2=0, key3=0;
1455 Int_t bin1=0, bin2=0, bin3=0;
1456 Float_t pVect1[4]={0};
1457 Float_t pVect2[4]={0};
1458 Float_t pVect3[4]={0};
1459 Float_t pVect1MC[4]={0};
1460 Float_t pVect2MC[4]={0};
1461 Int_t index1=0, index2=0, index3=0;
1462 Float_t weight12=0, weight13=0, weight23=0;
1463 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1464 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1465 Float_t weightTotal=0;//, weightTotalErr=0;
1468 AliAODMCParticle *mcParticle1=0x0;
1469 AliAODMCParticle *mcParticle2=0x0;
1472 if(fPdensityPairCut){
1473 ////////////////////
1474 Int_t pairCountSE=0, pairCountME=0;
1475 Int_t normPairCount[2]={0};
1476 Int_t numOtherPairs1[2][fMultLimit];
1477 Int_t numOtherPairs2[2][fMultLimit];
1478 Bool_t exitCode=kFALSE;
1479 Int_t tempNormFillCount[2][2][2][10][5];
1482 // reset to defaults
1483 for(Int_t i=0; i<fMultLimit; i++) {
1484 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1485 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1487 // Normalization Utilities
1488 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1489 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1490 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1491 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1492 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1493 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1494 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1495 numOtherPairs1[0][i]=0;
1496 numOtherPairs1[1][i]=0;
1497 numOtherPairs2[0][i]=0;
1498 numOtherPairs2[1][i]=0;
1500 // Track Merging/Splitting Utilities
1501 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1502 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1503 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1504 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1507 // Reset the temp Normalization counters
1508 for(Int_t i=0; i<2; i++){// Charge1
1509 for(Int_t j=0; j<2; j++){// Charge2
1510 for(Int_t k=0; k<2; k++){// Charge3
1511 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1512 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1513 tempNormFillCount[i][j][k][l][m] = 0;
1521 ///////////////////////////////////////////////////////
1522 // Start the pairing process
1526 for (Int_t i=0; i<myTracks; i++) {
1531 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1533 key1 = (fEvt)->fTracks[i].fKey;
1534 key2 = (fEvt+en2)->fTracks[j].fKey;
1535 Short_t fillIndex2 = FillIndex2part(key1+key2);
1536 Short_t qCutBin = SetQcutBin(fillIndex2);
1537 Short_t normBin = SetNormBin(fillIndex2);
1538 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1539 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1540 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1541 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1544 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1545 GetQosl(pVect1, pVect2, qout, qside, qlong);
1546 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1548 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1550 ///////////////////////////////
1551 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1552 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1553 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1555 if(fMCcase && ch1==ch2 && fMbin==0){
1556 for(Int_t rstep=0; rstep<10; rstep++){
1557 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1558 // propagate through B field to r=1.2m
1559 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1560 if(phi1 > 2*PI) phi1 -= 2*PI;
1561 if(phi1 < 0) phi1 += 2*PI;
1562 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1563 if(phi2 > 2*PI) phi2 -= 2*PI;
1564 if(phi2 < 0) phi2 += 2*PI;
1565 Float_t deltaphi = phi1 - phi2;
1566 if(deltaphi > PI) deltaphi -= PI;
1567 if(deltaphi < -PI) deltaphi += PI;
1568 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1571 // Pair Splitting/Merging cut
1573 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1574 fPairSplitCut[0][i]->AddAt('1',j);
1575 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1581 if(fMCcase && fillIndex2==0){
1583 // Check that label does not exceed stack size
1584 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1585 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1586 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1587 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1588 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1589 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1590 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1591 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1592 if(qinv12<0.1 && ch1==ch2) {
1593 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1594 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1595 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
1599 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 10fm
1600 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
1601 Int_t denIndex = rIter*kNDampValues + myDampIt;
1603 Float_t WInput = MCWeight(ch1,ch2, rIter, myDampIt, qinv12MC);
1604 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
1605 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
1606 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
1607 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
1611 //HIJING resonance test
1613 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
1614 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
1615 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
1616 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
1617 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
1618 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
1625 //////////////////////////////////////////
1627 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
1630 if(bin1==bin2 && fillIndex2==0){
1631 if((transK12 > 0.2) && (transK12 < 0.3)){
1632 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1633 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1635 if((transK12 > 0.6) && (transK12 < 0.7)){
1636 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1637 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1640 //////////////////////////////////////////
1642 if(fillIndex2==0 && bin1==bin2){
1644 transKbin=-1; rapKbin=-1;
1645 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1646 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1647 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1648 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1649 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1655 // exit out of loop if there are too many pairs
1656 if(pairCountSE >= kPairLimit) {cout<<"Too many SE pairs"<<endl; exitCode=kTRUE; continue;}
1657 if(exitCode) continue;
1659 //////////////////////////
1661 if(qinv12 <= fQcut[qCutBin]) {
1662 ///////////////////////////
1664 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
1665 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
1666 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
1667 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
1668 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
1669 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
1670 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
1671 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
1672 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1673 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1674 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1675 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
1678 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1679 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1680 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1681 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1682 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1683 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
1684 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
1685 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1686 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1687 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1688 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1689 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1692 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
1694 fPairLocationSE[i]->AddAt(pairCountSE,j);
1701 /////////////////////////////////////////////////////////
1702 // Normalization Region
1704 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1706 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1707 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1708 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1710 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1711 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1712 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1715 //other past pairs with particle j
1716 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
1717 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
1718 if(locationOtherPair < 0) continue;// no pair there
1719 Int_t indexOther1 = i;
1720 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
1722 // Both possible orderings of other indexes
1723 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
1725 // 1 and 2 are from SE
1726 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
1727 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
1728 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1729 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
1731 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
1734 }// pastpair P11 loop
1737 fNormPairSwitch[en2][i]->AddAt('1',j);
1738 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1739 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1741 numOtherPairs1[en2][i]++;
1742 numOtherPairs2[en2][j]++;
1745 normPairCount[en2]++;
1746 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1755 //////////////////////////////////////////////
1758 for (Int_t i=0; i<myTracks; i++) {
1763 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1765 key1 = (fEvt)->fTracks[i].fKey;
1766 key2 = (fEvt+en2)->fTracks[j].fKey;
1767 Short_t fillIndex2 = FillIndex2part(key1+key2);
1768 Short_t qCutBin = SetQcutBin(fillIndex2);
1769 Short_t normBin = SetNormBin(fillIndex2);
1770 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1771 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1772 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1773 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1775 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1776 GetQosl(pVect1, pVect2, qout, qside, qlong);
1777 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1779 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
1781 ///////////////////////////////
1782 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1783 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1784 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1786 if(fMCcase && ch1==ch2 && fMbin==0){
1787 for(Int_t rstep=0; rstep<10; rstep++){
1788 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1789 // propagate through B field to r=1.2m
1790 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1791 if(phi1 > 2*PI) phi1 -= 2*PI;
1792 if(phi1 < 0) phi1 += 2*PI;
1793 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1794 if(phi2 > 2*PI) phi2 -= 2*PI;
1795 if(phi2 < 0) phi2 += 2*PI;
1796 Float_t deltaphi = phi1 - phi2;
1797 if(deltaphi > PI) deltaphi -= PI;
1798 if(deltaphi < -PI) deltaphi += PI;
1799 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1804 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1805 fPairSplitCut[1][i]->AddAt('1',j);
1810 //////////////////////////////////////////
1812 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
1814 if(bin1==bin2 && fillIndex2==0){
1815 if((transK12 > 0.2) && (transK12 < 0.3)){
1816 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1817 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1819 if((transK12 > 0.6) && (transK12 < 0.7)){
1820 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1821 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1824 //////////////////////////////////////////
1826 if(fillIndex2==0 && bin1==bin2){
1828 transKbin=-1; rapKbin=-1;
1829 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1830 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1831 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1832 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1833 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1839 if(pairCountME >= 2*kPairLimit) {cout<<"Too many ME pairs"<<endl; exitCode=kTRUE; continue;}
1840 if(exitCode) continue;
1842 if(qinv12 <= fQcut[qCutBin]) {
1843 ///////////////////////////
1846 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
1847 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
1848 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
1849 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
1850 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
1851 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
1852 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
1853 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
1854 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1855 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1856 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1857 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
1860 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1861 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1862 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1863 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1864 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1865 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
1866 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
1867 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1868 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1869 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1870 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1871 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1874 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
1876 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
1882 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1884 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1885 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1886 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1888 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1889 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1890 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1892 //other past pairs in P11 with particle i
1893 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
1894 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
1895 if(locationOtherPairP11 < 0) continue;// no pair there
1896 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
1898 //Check other past pairs in P12
1899 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
1901 // 1 and 3 are from SE
1902 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
1903 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
1904 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1905 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
1906 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
1909 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
1910 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
1911 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
1917 fNormPairSwitch[en2][i]->AddAt('1',j);
1918 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1919 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1921 numOtherPairs1[en2][i]++;
1922 numOtherPairs2[en2][j]++;
1924 normPairCount[en2]++;
1925 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
1934 ///////////////////////////////////////
1935 // P13 pairing (just for Norm counting of term5)
1936 for (Int_t i=0; i<myTracks; i++) {
1938 // exit out of loop if there are too many pairs
1939 // dont bother with this loop if exitCode is set.
1945 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1947 key1 = (fEvt)->fTracks[i].fKey;
1948 key2 = (fEvt+en2)->fTracks[j].fKey;
1949 Short_t fillIndex2 = FillIndex2part(key1+key2);
1950 Short_t normBin = SetNormBin(fillIndex2);
1951 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1952 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1953 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1954 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1956 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1958 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
1960 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1961 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1964 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1965 fPairSplitCut[2][i]->AddAt('1',j);
1970 /////////////////////////////////////////////////////////
1971 // Normalization Region
1973 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1975 fNormPairSwitch[en2][i]->AddAt('1',j);
1983 ///////////////////////////////////////
1984 // P23 pairing (just for Norm counting of term5)
1986 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
1988 // exit out of loop if there are too many pairs
1989 // dont bother with this loop if exitCode is set.
1995 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
1999 key1 = (fEvt+en1)->fTracks[i].fKey;
2000 key2 = (fEvt+en2)->fTracks[j].fKey;
2001 Short_t fillIndex2 = FillIndex2part(key1+key2);
2002 Short_t normBin = SetNormBin(fillIndex2);
2003 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2004 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2005 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2006 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2008 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2010 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2012 ///////////////////////////////
2013 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2014 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2017 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2018 fPairSplitCut[3][i]->AddAt('1',j);
2023 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2025 Int_t index1P23 = i;
2026 Int_t index2P23 = j;
2028 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2029 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2030 if(locationOtherPairP12 < 0) continue; // no pair there
2031 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2034 //Check other past pair status in P13
2035 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2037 // all from different event
2038 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2039 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2040 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2041 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2043 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2051 ///////////////////////////////////////////////////
2052 // Do not use pairs from events with too many pairs
2054 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2055 (fEvt)->fNpairsSE = 0;
2056 (fEvt)->fNpairsME = 0;
2057 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2058 return;// Skip event
2060 (fEvt)->fNpairsSE = pairCountSE;
2061 (fEvt)->fNpairsME = pairCountME;
2062 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2064 ///////////////////////////////////////////////////
2068 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
2069 cout<<"Start Main analysis"<<endl;
2071 ///////////////////////////////////////////////////////////////////////
2072 ///////////////////////////////////////////////////////////////////////
2073 ///////////////////////////////////////////////////////////////////////
2076 // Start the Correlation Analysis
2079 ///////////////////////////////////////////////////////////////////////
2082 /////////////////////////////////////////////////////////
2083 // Skip 3-particle part if Tabulate6DPairs is set to true
2084 if(fTabulatePairs) return;
2085 /////////////////////////////////////////////////////////
2087 // Set the Normalization counters
2088 for(Int_t termN=0; termN<5; termN++){
2091 if((fEvt)->fNtracks ==0) continue;
2093 if((fEvt)->fNtracks ==0) continue;
2094 if((fEvt+1)->fNtracks ==0) continue;
2096 if((fEvt)->fNtracks ==0) continue;
2097 if((fEvt+1)->fNtracks ==0) continue;
2098 if((fEvt+2)->fNtracks ==0) continue;
2101 for(Int_t sc=0; sc<kSCLimit3; sc++){
2103 for(Int_t c1=0; c1<2; c1++){
2104 for(Int_t c2=0; c2<2; c2++){
2105 for(Int_t c3=0; c3<2; c3++){
2107 if(sc==0 || sc==6 || sc==9){// Identical species
2108 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2109 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2111 if( (c1+c2)==1) {if(c1!=0) continue;}
2112 }else {}// do nothing for pi-k-p case
2113 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2122 /////////////////////////////////////////////
2123 // Calculate Pair-Cut Correlations
2124 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2127 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2128 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2131 for(Int_t p1=0; p1<nump1; p1++){
2134 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2135 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2136 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2137 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2138 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2139 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2140 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2141 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2142 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2143 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2144 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2145 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2146 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2148 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2151 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2152 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2153 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2154 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2155 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2156 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2157 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2158 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2159 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2160 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2161 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2162 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2163 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2165 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2171 for(Int_t en2=0; en2<3; en2++){
2172 //////////////////////////////////////
2174 Bool_t skipcase=kTRUE;
2175 Short_t config=-1, part=-1;
2176 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2177 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2178 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2179 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2181 if(skipcase) continue;
2186 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2190 // remove auto-correlations and duplicate triplets
2192 if( index1 == index3) continue;
2193 if( index2 == index3) continue;
2194 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2196 // skip the switched off triplets
2197 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2198 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2201 ///////////////////////////////
2202 // Turn off 1st possible degenerate triplet
2203 if(index1 < index3){// verify correct id ordering ( index1 < k )
2204 if(fPairLocationSE[index1]->At(index3) >= 0){
2205 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2207 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2208 }else {// or k < index1
2209 if(fPairLocationSE[index3]->At(index1) >= 0){
2210 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2212 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2214 // turn off 2nd possible degenerate triplet
2215 if(index2 < index3){// verify correct id ordering (index2 < k)
2216 if(fPairLocationSE[index2]->At(index3) >= 0){
2217 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2219 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2220 }else {// or k < index2
2221 if(fPairLocationSE[index3]->At(index2) >= 0){
2222 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2224 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2229 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2230 ///////////////////////////////
2231 // Turn off 1st possible degenerate triplet
2232 if(fPairLocationME[index1]->At(index3) >= 0){
2233 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2236 // turn off 2nd possible degenerate triplet
2237 if(fPairLocationME[index2]->At(index3) >= 0){
2238 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2241 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2242 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2243 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2244 }// end config 2 part 1
2246 if(config==2 && part==2){// P12T1
2247 if( index1 == index3) continue;
2249 // skip the switched off triplets
2250 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2251 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2254 // turn off another possible degenerate
2255 if(fPairLocationME[index3]->At(index2) >= 0){
2256 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2257 }// end config 2 part 2
2259 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2260 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2261 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2262 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2264 if(config==3){// P12T3
2265 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2266 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2267 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2272 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2273 key3 = (fEvt+en2)->fTracks[k].fKey;
2274 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2275 Short_t fillIndex13 = FillIndex2part(key1+key3);
2276 Short_t fillIndex23 = FillIndex2part(key2+key3);
2277 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2278 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2279 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2280 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2281 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2282 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
2283 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2284 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2286 if(qinv13 < fQLowerCut) continue;
2287 if(qinv23 < fQLowerCut) continue;
2288 if(qinv13 > 3.*fQcut[qCutBin13]) continue;
2289 if(qinv23 > 3.*fQcut[qCutBin23]) continue;
2291 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2292 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2293 Float_t firstQ=0, secondQ=0, thirdQ=0;
2296 if(config==1) {// 123
2297 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2299 if(fillIndex3 <= 2){
2300 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
2301 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ);
2302 ((TH2F*)fOutputList->FindObject("fKt3Dist"))->Fill(fMbin+1, transK3);
2305 }else if(config==2){// 12, 13, 23
2307 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2308 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2310 // loop over terms 2-4
2311 for(Int_t jj=0; jj<3; jj++){
2312 if(jj==0) {if(!fill2) continue;}//12
2313 if(jj==1) {if(!fill3) continue;}//13
2314 if(jj==2) {if(!fill4) continue;}//23
2316 if(fillIndex3 <= 2){
2317 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj+2, firstQ, secondQ, thirdQ);
2318 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj+1].fTerms3->Fill(firstQ, secondQ, thirdQ);
2322 }else {// config 3: All particles from different events
2323 //Simulate Filling of other event-mixed lowQ pairs
2325 Float_t enhancement=1.0;
2327 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2328 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2330 if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2331 if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2332 if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
2334 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2337 if(fillIndex3 <= 2){
2338 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
2339 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ, enhancement);
2343 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2344 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
2345 if(fMCcase) continue;// only calcualte TPN for real data
2347 GetWeight(pVect1, pVect2, weight12, weight12Err);
2348 GetWeight(pVect1, pVect3, weight13, weight13Err);
2349 GetWeight(pVect2, pVect3, weight23, weight23Err);
2350 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2354 // Coul corrections from Wave-functions
2355 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 5fm to 20fm
2356 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
2357 Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2358 Int_t denIndex = rIter*kNDampValues + myDampIt;
2361 Float_t coulCorr12 = CoulCorr(+1,+1, rIter, qinv12);
2362 Float_t coulCorr13 = CoulCorr(+1,+1, rIter, qinv13);
2363 Float_t coulCorr23 = CoulCorr(+1,+1, rIter, qinv23);
2364 weight12CC = ((weight12+1)*GetMomRes(denIndex, qinv12) - myDamp/coulCorr12 - (1-myDamp));
2365 weight12CC *= coulCorr12/myDamp;
2366 weight13CC = ((weight13+1)*GetMomRes(denIndex, qinv13) - myDamp/coulCorr13 - (1-myDamp));
2367 weight13CC *= coulCorr13/myDamp;
2368 weight23CC = ((weight23+1)*GetMomRes(denIndex, qinv23) - myDamp/coulCorr23 - (1-myDamp));
2369 weight23CC *= coulCorr23/myDamp;
2371 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {
2372 if(denIndex==71 && fMbin==0 && ch1==-1) {
2373 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2374 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2376 continue;// C2^QS can never be less than unity
2379 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
2380 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2383 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2385 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2386 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2387 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2388 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2389 weightTotalErr /= pow(2*weightTotal,2);
2391 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, enhancement*weightTotalErr);
2401 }// end 3rd particle
2410 }// end of PdensityPairs
2418 ////////////////////////////////////////////////////////
2419 // Pdensity Method with Explicit Loops
2420 if(fPdensityExplicitLoop){
2422 ////////////////////////////////////
2423 // 2nd, 3rd, and 4th order Correlations
2426 for (Int_t i=0; i<myTracks; i++) {
2427 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2428 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2429 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2430 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2431 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2432 key1 = (fEvt)->fTracks[i].fKey;
2435 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2437 if(en2==0) startbin2=i+1;
2440 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2441 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2442 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2443 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2444 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2445 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2446 key2 = (fEvt+en2)->fTracks[j].fKey;
2448 Short_t fillIndex12 = FillIndex2part(key1+key2);
2449 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2451 if(qinv12 < fQLowerCut) continue;
2454 // 2-particle part is filled always during pair creator
2457 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2459 if(en3==en2) startbin3=j+1;
2464 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2465 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2466 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2467 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2468 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2469 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2470 key3 = (fEvt+en3)->fTracks[k].fKey;
2472 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2473 Short_t fillIndex13 = FillIndex2part(key1+key3);
2474 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2475 Short_t fillIndex23 = FillIndex2part(key2+key3);
2476 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2479 if(qinv13 < fQLowerCut) continue;
2480 if(qinv23 < fQLowerCut) continue;
2483 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2485 Short_t normBin12 = SetNormBin(fillIndex12);
2486 Short_t normBin13 = SetNormBin(fillIndex13);
2487 Short_t normBin23 = SetNormBin(fillIndex23);
2490 if(en3==0 && en2==0) {// 123
2491 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2493 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2495 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2496 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2497 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
2501 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2504 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2505 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2509 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2510 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2511 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2512 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2517 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2518 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2519 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2520 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2525 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2526 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2527 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2528 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2533 }else if(en2==1 && en3==2){// all uncorrelated events
2534 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2536 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
2537 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2538 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2539 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
2542 Short_t qCutBin12 = SetQcutBin(fillIndex12);
2543 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2544 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2546 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
2549 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
2550 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2551 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2569 }// End of PdensityExplicit
2574 // Post output data.
2575 PostData(1, fOutputList);
2578 //________________________________________________________________________
2579 void AliChaoticity::Terminate(Option_t *)
2581 // Called once at the end of the query
2586 //________________________________________________________________________
2587 Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
2590 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
2592 // propagate through B field to r=1m
2593 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// mine. 0.1798 for D=1.2m
2594 if(phi1 > 2*PI) phi1 -= 2*PI;
2595 if(phi1 < 0) phi1 += 2*PI;
2596 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// mine. 0.1798 for D=1.2m
2597 if(phi2 > 2*PI) phi2 -= 2*PI;
2598 if(phi2 < 0) phi2 += 2*PI;
2600 Float_t deltaphi = phi1 - phi2;
2601 if(deltaphi > PI) deltaphi -= 2*PI;
2602 if(deltaphi < -PI) deltaphi += 2*PI;
2603 deltaphi = fabs(deltaphi);
2605 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Sep just after TPC entrance
2607 // propagate through B field to r=1.6m
2608 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.1798 for D=1.2m
2609 if(phi1 > 2*PI) phi1 -= 2*PI;
2610 if(phi1 < 0) phi1 += 2*PI;
2611 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.1798 for D=1.2m
2612 if(phi2 > 2*PI) phi2 -= 2*PI;
2613 if(phi2 < 0) phi2 += 2*PI;
2615 deltaphi = phi1 - phi2;
2616 if(deltaphi > PI) deltaphi -= 2*PI;
2617 if(deltaphi < -PI) deltaphi += 2*PI;
2618 deltaphi = fabs(deltaphi);
2620 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Sep just after TPC entrance
2625 Int_t ncl1 = first.fClusterMap.GetNbits();
2626 Int_t ncl2 = second.fClusterMap.GetNbits();
2627 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2628 Double_t shfrac = 0; Double_t qfactor = 0;
2629 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2630 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2631 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2635 else {sumQ--; sumCls+=2;}
2637 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2642 qfactor = sumQ*1.0/sumCls;
2643 shfrac = sumSha*1.0/sumCls;
2646 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2653 //________________________________________________________________________
2654 Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2656 Float_t arg = G_Coeff/qinv;
2658 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2659 else {return (exp(-arg)-1)/(-arg);}
2662 //________________________________________________________________________
2663 void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2667 for (Int_t i = i1; i < i2+1; i++) {
2668 j = (Int_t) (gRandom->Rndm() * a);
2674 //________________________________________________________________________
2675 Short_t AliChaoticity::FillIndex2part(Short_t key){
2677 if(key==2) return 0;// pi-pi
2678 else if(key==11) return 1;// pi-k
2679 else if(key==101) return 2;// pi-p
2680 else if(key==20) return 3;// k-k
2681 else if(key==110) return 4;// k-p
2682 else return 5;// p-p
2684 //________________________________________________________________________
2685 Short_t AliChaoticity::FillIndex3part(Short_t key){
2687 if(key==3) return 0;// pi-pi-pi
2688 else if(key==12) return 1;// pi-pi-k
2689 else if(key==21) return 2;// k-k-pi
2690 else if(key==102) return 3;// pi-pi-p
2691 else if(key==201) return 4;// p-p-pi
2692 else if(key==111) return 5;// pi-k-p
2693 else if(key==30) return 6;// k-k-k
2694 else if(key==120) return 7;// k-k-p
2695 else if(key==210) return 8;// p-p-k
2696 else return 9;// p-p-p
2699 //________________________________________________________________________
2700 Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
2701 if(fi <= 2) return 0;
2702 else if(fi==3) return 1;
2705 //________________________________________________________________________
2706 Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
2708 else if(fi <= 2) return 1;
2711 //________________________________________________________________________
2712 void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
2714 if(fi==0 || fi==3 || fi==5){// Identical species
2715 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
2716 else {b1=c1; b2=c2;}
2717 }else {// Mixed species
2718 if(key1 < key2) { b1=c1; b2=c2;}
2719 else {b1=c2; b2=c1;}
2723 //________________________________________________________________________
2724 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){
2727 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
2730 Short_t seKeySum=0;// only used for pi-k-p case
2731 if(part==1) {// default case (irrelevant for term 1 and term 5)
2732 if(c1==c2) seSS=kTRUE;
2733 if(key1==key2) seSK=kTRUE;
2734 seKeySum = key1+key2;
2737 if(c1==c3) seSS=kTRUE;
2738 if(key1==key3) seSK=kTRUE;
2739 seKeySum = key1+key3;
2743 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
2745 if(fi==0 || fi==6 || fi==9){// Identical species
2746 if( (c1+c2+c3)==1) {
2747 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
2749 if(seSS) fill2=kTRUE;
2750 else {fill3=kTRUE; fill4=kTRUE;}
2752 }else if( (c1+c2+c3)==2) {
2755 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
2759 b1=c1; b2=c2; b3=c3;
2760 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
2762 }else if(fi != 5){// all the rest except pi-k-p
2765 if( (c1+c2)==1) {b1=0; b2=1;}
2766 else {b1=c1; b2=c2;}
2767 }else if(key1==key3){
2769 if( (c1+c3)==1) {b1=0; b2=1;}
2770 else {b1=c1; b2=c3;}
2771 }else {// Key2==Key3
2773 if( (c2+c3)==1) {b1=0; b2=1;}
2774 else {b1=c2; b2=c3;}
2776 //////////////////////////////
2777 if(seSK) fill2=kTRUE;// Same keys from Same Event
2778 else {// Different keys from Same Event
2779 if( (c1+c2+c3)==1) {
2781 if(seSS) fill3=kTRUE;
2783 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
2784 }else if( (c1+c2+c3)==2) {
2786 if(seSS) fill4=kTRUE;
2788 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
2789 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
2791 /////////////////////////////
2792 }else {// pi-k-p (no charge ordering applies since all are unique)
2794 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
2795 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
2797 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
2798 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
2800 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
2801 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
2803 ////////////////////////////////////
2804 if(seKeySum==11) fill2=kTRUE;
2805 else if(seKeySum==101) fill3=kTRUE;
2807 ////////////////////////////////////
2811 //________________________________________________________________________
2812 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){
2814 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
2815 if(fi==0 || fi==6 || fi==9){// Identical species
2816 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
2817 if(term==1 || term==5){
2818 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
2819 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
2820 else {fQ=q23; sQ=q12; tQ=q13;}
2821 }else if(term==2 && part==1){
2822 fQ=q12; sQ=q13; tQ=q23;
2823 }else if(term==2 && part==2){
2824 fQ=q13; sQ=q12; tQ=q23;
2825 }else if(term==3 && part==1){
2827 if(c1==c3) {fQ=q13; tQ=q23;}
2828 else {fQ=q23; tQ=q13;}
2829 }else if(term==3 && part==2){
2831 if(c1==c2) {fQ=q12; tQ=q23;}
2832 else {fQ=q23; tQ=q12;}
2833 }else if(term==4 && part==1){
2835 if(c1==c3) {fQ=q13; sQ=q23;}
2836 else {fQ=q23; sQ=q13;}
2837 }else if(term==4 && part==2){
2839 if(c1==c2) {fQ=q12; sQ=q23;}
2840 else {fQ=q23; sQ=q12;}
2841 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
2842 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
2843 if(term==1 || term==5){
2844 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
2845 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
2846 else {tQ=q23; sQ=q12; fQ=q13;}
2847 }else if(term==2 && part==1){
2849 if(c1==c3) {tQ=q13; sQ=q23;}
2850 else {tQ=q23; sQ=q13;}
2851 }else if(term==2 && part==2){
2853 if(c1==c2) {tQ=q12; sQ=q23;}
2854 else {tQ=q23; sQ=q12;}
2855 }else if(term==3 && part==1){
2857 if(c1==c3) {tQ=q13; fQ=q23;}
2858 else {tQ=q23; fQ=q13;}
2859 }else if(term==3 && part==2){
2861 if(c1==c2) {tQ=q12; fQ=q23;}
2862 else {tQ=q23; fQ=q12;}
2863 }else if(term==4 && part==1){
2864 tQ=q12; sQ=q13; fQ=q23;
2865 }else if(term==4 && part==2){
2866 tQ=q13; sQ=q12; fQ=q23;
2867 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
2868 }else {// fQ=ss, sQ=ss, tQ=ss
2869 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
2870 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
2871 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
2872 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
2873 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
2874 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
2875 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
2877 }else if(fi != 5){// all the rest except pi-k-p
2881 // cases not explicity shown below are not possible
2882 if(term==1 || term==5) {sQ=q13; tQ=q23;}
2883 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
2884 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
2885 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
2886 else cout<<"problem!!!!!!!!!!!!!"<<endl;
2888 if(c1==c3) {sQ=q13; tQ=q23;}
2889 else {sQ=q23; tQ=q13;}
2891 if(c1==c3) {tQ=q13; sQ=q23;}
2892 else {tQ=q23; sQ=q13;}
2894 }else if(key1==key3){
2897 // cases not explicity shown below are not possible
2898 if(term==1 || term==5) {sQ=q12; tQ=q23;}
2899 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
2900 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
2901 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
2902 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
2904 if(c1==c2) {sQ=q12; tQ=q23;}
2905 else {sQ=q23; tQ=q12;}
2907 if(c1==c2) {tQ=q12; sQ=q23;}
2908 else {tQ=q23; sQ=q12;}
2910 }else {// key2==key3
2913 // cases not explicity shown below are not possible
2914 if(term==1 || term==5) {sQ=q12; tQ=q13;}
2915 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
2916 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
2917 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
2918 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
2919 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
2921 if(c1==c2) {sQ=q12; tQ=q13;}
2922 else {sQ=q13; tQ=q12;}
2924 if(c1==c2) {tQ=q12; sQ=q13;}
2925 else {tQ=q13; sQ=q12;}
2930 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
2931 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
2933 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
2934 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
2936 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
2937 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
2944 //________________________________________________________________________
2945 Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
2949 if(fi==0 || fi==3 || fi==5){// identical masses
2950 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));
2951 }else{// different masses
2952 Float_t px = track1[1] + track2[1];
2953 Float_t py = track1[2] + track2[2];
2954 Float_t pz = track1[3] + track2[3];
2955 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
2956 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
2957 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
2959 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
2960 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
2961 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
2962 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
2969 //________________________________________________________________________
2970 void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
2972 Float_t p0 = track1[0] + track2[0];
2973 Float_t px = track1[1] + track2[1];
2974 Float_t py = track1[2] + track2[2];
2975 Float_t pz = track1[3] + track2[3];
2977 Float_t mt = sqrt(p0*p0 - pz*pz);
2978 Float_t pt = sqrt(px*px + py*py);
2980 Float_t v0 = track1[0] - track2[0];
2981 Float_t vx = track1[1] - track2[1];
2982 Float_t vy = track1[2] - track2[2];
2983 Float_t vz = track1[3] - track2[3];
2985 qout = (px*vx + py*vy)/pt;
2986 qside = (px*vy - py*vx)/pt;
2987 qlong = (p0*vz - pz*v0)/mt;
2989 //________________________________________________________________________
2990 void AliChaoticity::SetWeightArrays(){
2992 for(Int_t mb=0; mb<fMbins; mb++){
2993 for(Int_t edB=0; edB<kEDbins; edB++){
2994 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
2995 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
2997 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
2998 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
2999 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3001 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3002 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3014 TFile *wFile = new TFile("WeightFile.root","READ");
3015 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3016 else cout<<"Good Weight File Found!"<<endl;
3018 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3019 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3020 for(Int_t mb=0; mb<fMbins; mb++){
3021 for(Int_t edB=0; edB<kEDbins; edB++){
3023 TString *name = new TString("Weight_Kt_");
3025 name->Append("_Ky_");
3027 name->Append("_M_");
3029 name->Append("_ED_");
3032 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3034 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3035 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3036 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3038 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3039 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3052 cout<<"Done Setting Weights"<<endl;
3055 //________________________________________________________________________
3056 void AliChaoticity::SetWeightArraysLEGO(TH3F *histos[kKbinsT][kCentBins]){
3058 for(Int_t mb=0; mb<fMbins; mb++){
3059 for(Int_t edB=0; edB<kEDbins; edB++){
3060 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3061 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3063 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3064 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3065 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3067 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3068 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
3078 //________________________________________________________________________
3079 void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3081 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3082 Float_t ky=0;// central rapidity
3084 Float_t qOut=0,qSide=0,qLong=0;
3085 GetQosl(track1, track2, qOut, qSide, qLong);
3087 qSide = fabs(qSide);
3088 qLong = fabs(qLong);
3091 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3092 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3094 for(Int_t i=0; i<kKbinsT-1; i++){
3095 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3099 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3100 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3102 for(Int_t i=0; i<kKbinsY-1; i++){
3103 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3108 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3109 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3111 for(Int_t i=0; i<kQbinsWeights-1; i++){
3112 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3116 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3117 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3119 for(Int_t i=0; i<kQbinsWeights-1; i++){
3120 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3124 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3125 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3127 for(Int_t i=0; i<kQbinsWeights-1; i++){
3128 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3134 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3135 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3140 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3142 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3144 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstep;
3146 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstep;
3148 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstep;
3156 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3157 Float_t deltaWErr=0;
3160 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3162 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3164 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstep;
3166 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstep;
3168 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstep;
3170 wgtErr = minErr + deltaWErr;
3175 //________________________________________________________________________
3176 Float_t AliChaoticity::CoulCorr(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
3178 if(rIndex >= kRVALUES) return 1.0;
3179 Int_t indexL = Int_t(fabs(qinv*1000. - 2.)/2.);// starts at 2MeV
3180 Int_t indexH = indexL+1;
3181 if(indexL >= kNlinesCoulFile-1) return 1.0;
3184 if(charge1==charge2){
3185 slope = (fCoulSS[rIndex][indexL] - fCoulSS[rIndex][indexH])/(fQCoul[indexL]-fQCoul[indexH]);
3186 return (slope*(qinv - fQCoul[indexL]) + fCoulSS[rIndex][indexL]);
3188 slope = (fCoulOS[rIndex][indexL] - fCoulOS[rIndex][indexH])/(fQCoul[indexL]-fQCoul[indexH]);
3189 return (slope*(qinv - fQCoul[indexL]) + fCoulOS[rIndex][indexL]);
3192 //________________________________________________________________________
3193 void AliChaoticity::SetCoulCorrections(){
3195 // Set default values
3196 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3198 for(Int_t j=0; j<kRVALUES; j++) {// radii columns
3199 fCoulSS[j][i]=-1;// should cause a crash in the sqrt() function if ifstream not successful below = (exit the code).
3200 fCoulOS[j][i]=-1;// should cause a crash in the sqrt() function if ifstream not successful below = (exit the code).
3204 ifstream mystream("cc2_l100_all.txt");
3205 for(Int_t j=0; j<kRVALUES; j++) {// radii columns (3-10fm in cc2_l100_all.txt)
3206 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3207 mystream >> fQCoul[i];// Q2
3209 mystream >> fCoulSS[j][i];
3210 mystream >> fCoulOS[j][i];
3215 //________________________________________________________________________
3216 void AliChaoticity::SetCoulCorrectionsLEGO(Float_t qPoints[kNlinesCoulFile], Float_t coulSS[kRVALUES][kNlinesCoulFile], Float_t coulOS[kRVALUES][kNlinesCoulFile]){
3218 // Set default values
3219 for(Int_t i=0; i<kNlinesCoulFile; i++) {
3220 fQCoul[i] = qPoints[i];
3221 for(Int_t j=0; j<kRVALUES; j++) {// radii columns
3222 fCoulSS[j][i] = coulSS[j][i];
3223 fCoulOS[j][i] = coulOS[j][i];
3228 //________________________________________________________________________
3229 Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3231 Float_t radius = Float_t(rIndex+3);
3232 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3233 Float_t coulCorr12 = CoulCorr(charge1, charge2, rIndex, qinv);
3234 if(charge1==charge2){
3235 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius/0.19733,2)))/coulCorr12);
3237 return ((1-myDamp) + myDamp/coulCorr12);
3241 //________________________________________________________________________
3242 Float_t AliChaoticity::MCWeightr3(Int_t term, Int_t rIndex, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
3243 if(term==5) return 1.0;
3244 Float_t radius = 3. + rIndex;//starts at 5fm
3245 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3246 Float_t fc = sqrt(myDamp);
3247 Float_t coulCorr12 = CoulCorr(+1,+1, rIndex, q12);
3248 Float_t coulCorr13 = CoulCorr(+1,+1, rIndex, q13);
3249 Float_t coulCorr23 = CoulCorr(+1,+1, rIndex, q23);
3252 Float_t c3QS = 1 + exp(-pow(q12*radius/0.19733,2)) + exp(-pow(q13*radius/0.19733,2)) + exp(-pow(q23*radius/0.19733,2));
3253 c3QS += 2*exp(-pow(radius/0.19733,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3254 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3255 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius/0.19733,2)))/coulCorr12;
3256 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius/0.19733,2)))/coulCorr13;
3257 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius/0.19733,2)))/coulCorr23;
3258 w123 += pow(fc,3)*c3QS/(coulCorr12*coulCorr13*coulCorr23);
3261 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius/0.19733,2)))/coulCorr12);
3263 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius/0.19733,2)))/coulCorr13);
3265 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius/0.19733,2)))/coulCorr23);
3269 //________________________________________________________________________
3270 Float_t AliChaoticity::MCWeightOSL(Float_t radius, Float_t myDamp, Float_t Q_out, Float_t Q_side, Float_t Q_long, Float_t qinv){
3271 Int_t rIndex = Int_t(radius+0.001-3);
3272 Float_t coulCorr12 = CoulCorr(+1,+1, rIndex, qinv);
3274 return ((1-myDamp) + myDamp*(1 + exp(-pow(radius/0.19733,2)*(Q_out*Q_out+Q_side*Q_side+Q_long*Q_long)))/coulCorr12);
3277 //________________________________________________________________________
3278 void AliChaoticity::SetMomResCorrections(){
3280 for(Int_t i=0; i<kDENtypes; i++){
3281 for(Int_t j=0; j<kQbinsWeights; j++){
3282 fMomResWeights[i][j]=1.0;
3286 TFile *momResFile = new TFile("MomResFile.root","READ");
3287 if(!momResFile->IsOpen()) {cout<<"No Momentum Resolution File!!!!!!!!!!"<<endl; return;}
3288 else cout<<"Good Momentum Resolution File Found!"<<endl;
3290 TH2D *fMomResWeights_pp = (TH2D*)momResFile->Get("MomResHisto_pp");
3292 for(Int_t i=0; i<fMomResWeights_pp->GetNbinsX(); i++){
3293 for(Int_t j=0; j<fMomResWeights_pp->GetNbinsY(); j++){
3295 fMomResWeights[i][j]=Float_t(fMomResWeights_pp->GetBinContent(i+1,j+1));
3300 momResFile->Close();
3303 //________________________________________________________________________
3304 void AliChaoticity::SetMomResCorrectionsLEGO(TH2D *histo){
3305 for(Int_t i=0; i<histo->GetNbinsX(); i++){
3306 for(Int_t j=0; j<histo->GetNbinsY(); j++){
3307 fMomResWeights[i][j]=Float_t(histo->GetBinContent(i+1,j+1));
3311 //________________________________________________________________________
3312 Float_t AliChaoticity::GetMomRes(Int_t denIndex, Float_t qinv){
3314 Int_t qbinL = Int_t(fabs(qinv*1000. - 2.5)/5.);// starts at 2.5MeV (center of bin)
3315 Int_t qbinH = qbinL+1;
3316 if(qbinL >=kQbinsWeights-1) return 1.0;
3317 Float_t slope = (fMomResWeights[denIndex][qbinL] - fMomResWeights[denIndex][qbinH])/(-0.005);
3318 return (slope*(qinv - (0.005*(qbinL+0.5))) + fMomResWeights[denIndex][qbinL]);
3321 //________________________________________________________________________