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