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