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