PWGCF/FEMTOSCOPY/Chaoticity/AliChaoticity.cxx

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