PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.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 #include "AliAnalysisUtils.h"
  32
  33 #include "AliFourPion.h"
  34
  35 #define PI 3.1415927
  36 #define G_Coeff 0.006399 // 2*pi*alpha*M_pion
  37 #define kappa3 0.1 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
  38 #define kappa4 0.5 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
  39
  40
  41 // Author: Dhevan Gangadharan
  42
  43 ClassImp(AliFourPion)
  44
  45 //________________________________________________________________________
  46 AliFourPion::AliFourPion():
  47 AliAnalysisTaskSE(),
  48   fname(0),
  49   fAOD(0x0),
  50   fOutputList(0x0),
  51   fPIDResponse(0x0),
  52   fEC(0x0),
  53   fEvt(0x0),
  54   fTempStruct(0x0),
  55   fRandomNumber(0x0),
  56   fLEGO(kTRUE),
  57   fMCcase(kFALSE),
  58   fAODcase(kTRUE),
  59   fPbPbcase(kTRUE),
  60   fGenerateSignal(kFALSE),
  61   fGeneratorOnly(kFALSE),
  62   fPdensityPairCut(kTRUE),
  63   fTabulatePairs(kFALSE),
  64   fRMax(11),
  65   ffcSq(0.7),
  66   fFilterBit(7),
  67   fMaxChi2NDF(10),
  68   fMinTPCncls(0),
  69   fBfield(0),
  70   fMbin(0),
  71   fFSIindex(0),
  72   fEDbin(0),
  73   fMbins(fCentBins),
  74   fMultLimit(0),
  75   fCentBinLowLimit(0),
  76   fCentBinHighLimit(1),
  77   fEventCounter(0),
  78   fEventsToMix(0),
  79   fZvertexBins(0),
  80   fMultLimits(),
  81   fQcut(0),
  82   fQLowerCut(0),
  83   fNormQcutLow(0),
  84   fNormQcutHigh(0),
  85   fKupperBound(0),
  86   fQupperBoundQ2(0),
  87   fQupperBoundQ3(0),
  88   fQupperBoundQ4(0),
  89   fQbinsQ2(1),
  90   fQbinsQ3(1),
  91   fQbinsQ4(1),
  92   fQupperBoundWeights(0),
  93   fKstepT(),
  94   fKstepY(),
  95   fKmeanT(),
  96   fKmeanY(),
  97   fKmiddleT(),
  98   fKmiddleY(),
  99   fQstep(0),
 100   fQstepWeights(0),
 101   fQmean(),
 102   fDampStart(0),
 103   fDampStep(0),
 104   fTPCTOFboundry(0),
 105   fTOFboundry(0),
 106   fSigmaCutTPC(2.0),
 107   fSigmaCutTOF(2.0),
 108   fMinSepPairEta(0.03),
 109   fMinSepPairPhi(0.04),
 110   fShareQuality(0),
 111   fShareFraction(0),
 112   fTrueMassP(0),
 113   fTrueMassPi(0),
 114   fTrueMassK(0),
 115   fTrueMassKs(0),
 116   fTrueMassLam(0),
 117   fKtIndexL(0),
 118   fKtIndexH(0),
 119   fQoIndexL(0),
 120   fQoIndexH(0),
 121   fQsIndexL(0),
 122   fQsIndexH(0),
 123   fQlIndexL(0),
 124   fQlIndexH(0),
 125   fDummyB(0),
 126   fKT3transition(0.3),
 127   fKT4transition(0.3),
 128   fDefaultsCharSwitch(),
 129   fLowQPairSwitch_E0E0(),
 130   fLowQPairSwitch_E0E1(),
 131   fLowQPairSwitch_E0E2(),
 132   fLowQPairSwitch_E0E3(),
 133   fLowQPairSwitch_E1E1(),
 134   fLowQPairSwitch_E1E2(),
 135   fLowQPairSwitch_E1E3(),
 136   fLowQPairSwitch_E2E3(),
 137   fNormQPairSwitch_E0E0(),
 138   fNormQPairSwitch_E0E1(),
 139   fNormQPairSwitch_E0E2(),
 140   fNormQPairSwitch_E0E3(),
 141   fNormQPairSwitch_E1E1(),
 142   fNormQPairSwitch_E1E2(),
 143   fNormQPairSwitch_E1E3(),
 144   fNormQPairSwitch_E2E3(),
 145   fMomResC2(0x0),
 146   fWeightmuonCorrection(0x0)
 147 {
 148   // Default constructor
 149   for(Int_t mb=0; mb<fMbins; mb++){
 150     for(Int_t edB=0; edB<fEDbins; edB++){
 151       for(Int_t c1=0; c1<2; c1++){
 152         for(Int_t c2=0; c2<2; c2++){
 153           for(Int_t term=0; term<2; term++){
 154
 155             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
 156
 157             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
 158             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
 159             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
 160             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
 161             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
 162             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
 163
 164           }// term_2
 165
 166
 167           for(Int_t c3=0; c3<2; c3++){
 168             for(Int_t term=0; term<5; term++){
 169
 170               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
 171               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
 172               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
 173               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
 174
 175             }// term_3
 176
 177             for(Int_t c4=0; c4<2; c4++){
 178               for(Int_t term=0; term<13; term++){
 179
 180                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
 181                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
 182                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
 183                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
 184
 185               }// term_4
 186
 187             }// c4
 188           }//c3
 189         }//c2
 190       }//c1
 191       for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
 192         for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
 193           KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
 194           KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
 195         }
 196       }
 197
 198     }// ED
 199   }// Mbin
 200
 201   // Initialize FSI histograms
 202   for(Int_t i=0; i<12; i++){
 203     fFSIss[i]=0x0;
 204     fFSIos[i]=0x0;
 205   }
 206
 207
 208   // Initialize fNormWeight and fNormWeightErr to 0
 209   for(Int_t i=0; i<3; i++){// Kt iterator
 210     for(Int_t j=0; j<10; j++){// Mbin iterator
 211       fNormWeight[i][j]=0x0;
 212     }
 213   }
 214
 215
 216 }
 217 //________________________________________________________________________
 218 AliFourPion::AliFourPion(const Char_t *name)
 219 : AliAnalysisTaskSE(name),
 220   fname(name),
 221   fAOD(0x0),
 222   fOutputList(0x0),
 223   fPIDResponse(0x0),
 224   fEC(0x0),
 225   fEvt(0x0),
 226   fTempStruct(0x0),
 227   fRandomNumber(0x0),
 228   fLEGO(kTRUE),
 229   fMCcase(kFALSE),
 230   fAODcase(kTRUE),
 231   fPbPbcase(kTRUE),
 232   fGenerateSignal(kFALSE),
 233   fGeneratorOnly(kFALSE),
 234   fPdensityPairCut(kTRUE),
 235   fTabulatePairs(kFALSE),
 236   fRMax(11),
 237   ffcSq(0.7),
 238   fFilterBit(7),
 239   fMaxChi2NDF(10),
 240   fMinTPCncls(0),
 241   fBfield(0),
 242   fMbin(0),
 243   fFSIindex(0),
 244   fEDbin(0),
 245   fMbins(fCentBins),
 246   fMultLimit(0),
 247   fCentBinLowLimit(0),
 248   fCentBinHighLimit(1),
 249   fEventCounter(0),
 250   fEventsToMix(0),
 251   fZvertexBins(0),
 252   fMultLimits(),
 253   fQcut(0),
 254   fQLowerCut(0),
 255   fNormQcutLow(0),
 256   fNormQcutHigh(0),
 257   fKupperBound(0),
 258   fQupperBoundQ2(0),
 259   fQupperBoundQ3(0),
 260   fQupperBoundQ4(0),
 261   fQbinsQ2(1),
 262   fQbinsQ3(1),
 263   fQbinsQ4(1),
 264   fQupperBoundWeights(0),
 265   fKstepT(),
 266   fKstepY(),
 267   fKmeanT(),
 268   fKmeanY(),
 269   fKmiddleT(),
 270   fKmiddleY(),
 271   fQstep(0),
 272   fQstepWeights(0),
 273   fQmean(),
 274   fDampStart(0),
 275   fDampStep(0),
 276   fTPCTOFboundry(0),
 277   fTOFboundry(0),
 278   fSigmaCutTPC(2.0),
 279   fSigmaCutTOF(2.0),
 280   fMinSepPairEta(0.03),
 281   fMinSepPairPhi(0.04),
 282   fShareQuality(0),
 283   fShareFraction(0),
 284   fTrueMassP(0),
 285   fTrueMassPi(0),
 286   fTrueMassK(0),
 287   fTrueMassKs(0),
 288   fTrueMassLam(0),
 289   fKtIndexL(0),
 290   fKtIndexH(0),
 291   fQoIndexL(0),
 292   fQoIndexH(0),
 293   fQsIndexL(0),
 294   fQsIndexH(0),
 295   fQlIndexL(0),
 296   fQlIndexH(0),
 297   fDummyB(0),
 298   fKT3transition(0.3),
 299   fKT4transition(0.3),
 300   fDefaultsCharSwitch(),
 301   fLowQPairSwitch_E0E0(),
 302   fLowQPairSwitch_E0E1(),
 303   fLowQPairSwitch_E0E2(),
 304   fLowQPairSwitch_E0E3(),
 305   fLowQPairSwitch_E1E1(),
 306   fLowQPairSwitch_E1E2(),
 307   fLowQPairSwitch_E1E3(),
 308   fLowQPairSwitch_E2E3(),
 309   fNormQPairSwitch_E0E0(),
 310   fNormQPairSwitch_E0E1(),
 311   fNormQPairSwitch_E0E2(),
 312   fNormQPairSwitch_E0E3(),
 313   fNormQPairSwitch_E1E1(),
 314   fNormQPairSwitch_E1E2(),
 315   fNormQPairSwitch_E1E3(),
 316   fNormQPairSwitch_E2E3(),
 317   fMomResC2(0x0),
 318   fWeightmuonCorrection(0x0)
 319 {
 320   // Main constructor
 321   fAODcase=kTRUE;
 322   fPdensityPairCut = kTRUE;
 323
 324
 325   for(Int_t mb=0; mb<fMbins; mb++){
 326     for(Int_t edB=0; edB<fEDbins; edB++){
 327       for(Int_t c1=0; c1<2; c1++){
 328         for(Int_t c2=0; c2<2; c2++){
 329           for(Int_t term=0; term<2; term++){
 330
 331             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
 332
 333             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
 334             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
 335             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
 336             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
 337             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
 338             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
 339
 340           }// term_2
 341
 342           for(Int_t c3=0; c3<2; c3++){
 343             for(Int_t term=0; term<5; term++){
 344
 345               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
 346               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
 347               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
 348               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
 349
 350             }// term_3
 351
 352             for(Int_t c4=0; c4<2; c4++){
 353               for(Int_t term=0; term<13; term++){
 354
 355                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
 356                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
 357                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
 358                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
 359
 360               }// term_4
 361             }// c4
 362           }//c3
 363         }//c2
 364       }//c1
 365
 366       for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
 367         for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
 368           KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
 369           KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
 370         }
 371       }
 372
 373     }// ED
 374   }// Mbin
 375
 376   // Initialize FSI histograms
 377   for(Int_t i=0; i<12; i++){
 378     fFSIss[i]=0x0;
 379     fFSIos[i]=0x0;
 380   }
 381
 382   // Initialize fNormWeight and fNormWeightErr to 0
 383   for(Int_t i=0; i<3; i++){// Kt iterator
 384     for(Int_t j=0; j<10; j++){// Mbin iterator
 385       fNormWeight[i][j]=0x0;
 386     }
 387   }
 388
 389
 390   DefineOutput(1, TList::Class());
 391 }
 392 //________________________________________________________________________
 393 AliFourPion::AliFourPion(const AliFourPion &obj)
 394   : AliAnalysisTaskSE(obj.fname),
 395     fname(obj.fname),
 396     fAOD(obj.fAOD),
 397     //fESD(obj.fESD),
 398     fOutputList(obj.fOutputList),
 399     fPIDResponse(obj.fPIDResponse),
 400     fEC(obj.fEC),
 401     fEvt(obj.fEvt),
 402     fTempStruct(obj.fTempStruct),
 403     fRandomNumber(obj.fRandomNumber),
 404     fLEGO(obj.fLEGO),
 405     fMCcase(obj.fMCcase),
 406     fAODcase(obj.fAODcase),
 407     fPbPbcase(obj.fPbPbcase),
 408     fGenerateSignal(obj.fGenerateSignal),
 409     fGeneratorOnly(obj.fGeneratorOnly),
 410     fPdensityPairCut(obj.fPdensityPairCut),
 411     fTabulatePairs(obj.fTabulatePairs),
 412     fRMax(obj.fRMax),
 413     ffcSq(obj.ffcSq),
 414     fFilterBit(obj.fFilterBit),
 415     fMaxChi2NDF(obj.fMaxChi2NDF),
 416     fMinTPCncls(obj.fMinTPCncls),
 417     fBfield(obj.fBfield),
 418     fMbin(obj.fMbin),
 419     fFSIindex(obj.fFSIindex),
 420     fEDbin(obj.fEDbin),
 421     fMbins(obj.fMbins),
 422     fMultLimit(obj.fMultLimit),
 423     fCentBinLowLimit(obj.fCentBinLowLimit),
 424     fCentBinHighLimit(obj.fCentBinHighLimit),
 425     fEventCounter(obj.fEventCounter),
 426     fEventsToMix(obj.fEventsToMix),
 427     fZvertexBins(obj.fZvertexBins),
 428     fMultLimits(),
 429     fQcut(0),
 430     fQLowerCut(obj.fQLowerCut),
 431     fNormQcutLow(0),
 432     fNormQcutHigh(0),
 433     fKupperBound(obj.fKupperBound),
 434     fQupperBoundQ2(obj.fQupperBoundQ2),
 435     fQupperBoundQ3(obj.fQupperBoundQ3),
 436     fQupperBoundQ4(obj.fQupperBoundQ4),
 437     fQbinsQ2(obj.fQbinsQ2),
 438     fQbinsQ3(obj.fQbinsQ3),
 439     fQbinsQ4(obj.fQbinsQ4),
 440     fQupperBoundWeights(obj.fQupperBoundWeights),
 441     fKstepT(),
 442     fKstepY(),
 443     fKmeanT(),
 444     fKmeanY(),
 445     fKmiddleT(),
 446     fKmiddleY(),
 447     fQstep(obj.fQstep),
 448     fQstepWeights(obj.fQstepWeights),
 449     fQmean(),
 450     fDampStart(obj.fDampStart),
 451     fDampStep(obj.fDampStep),
 452     fTPCTOFboundry(obj.fTPCTOFboundry),
 453     fTOFboundry(obj.fTOFboundry),
 454     fSigmaCutTPC(obj.fSigmaCutTPC),
 455     fSigmaCutTOF(obj.fSigmaCutTOF),
 456     fMinSepPairEta(obj.fMinSepPairEta),
 457     fMinSepPairPhi(obj.fMinSepPairPhi),
 458     fShareQuality(obj.fShareQuality),
 459     fShareFraction(obj.fShareFraction),
 460     fTrueMassP(obj.fTrueMassP),
 461     fTrueMassPi(obj.fTrueMassPi),
 462     fTrueMassK(obj.fTrueMassK),
 463     fTrueMassKs(obj.fTrueMassKs),
 464     fTrueMassLam(obj.fTrueMassLam),
 465     fKtIndexL(obj.fKtIndexL),
 466     fKtIndexH(obj.fKtIndexH),
 467     fQoIndexL(obj.fQoIndexL),
 468     fQoIndexH(obj.fQoIndexH),
 469     fQsIndexL(obj.fQsIndexL),
 470     fQsIndexH(obj.fQsIndexH),
 471     fQlIndexL(obj.fQlIndexL),
 472     fQlIndexH(obj.fQlIndexH),
 473     fDummyB(obj.fDummyB),
 474     fKT3transition(obj.fKT3transition),
 475     fKT4transition(obj.fKT4transition),
 476     fDefaultsCharSwitch(),
 477     fLowQPairSwitch_E0E0(),
 478     fLowQPairSwitch_E0E1(),
 479     fLowQPairSwitch_E0E2(),
 480     fLowQPairSwitch_E0E3(),
 481     fLowQPairSwitch_E1E1(),
 482     fLowQPairSwitch_E1E2(),
 483     fLowQPairSwitch_E1E3(),
 484     fLowQPairSwitch_E2E3(),
 485     fNormQPairSwitch_E0E0(),
 486     fNormQPairSwitch_E0E1(),
 487     fNormQPairSwitch_E0E2(),
 488     fNormQPairSwitch_E0E3(),
 489     fNormQPairSwitch_E1E1(),
 490     fNormQPairSwitch_E1E2(),
 491     fNormQPairSwitch_E1E3(),
 492     fNormQPairSwitch_E2E3(),
 493     fMomResC2(obj.fMomResC2),
 494     fWeightmuonCorrection(obj.fWeightmuonCorrection)
 495 {
 496   // Copy Constructor
 497
 498   for(Int_t i=0; i<12; i++){
 499     fFSIss[i]=obj.fFSIss[i];
 500     fFSIos[i]=obj.fFSIos[i];
 501   }
 502
 503   // Initialize fNormWeight and fNormWeightErr to 0
 504   for(Int_t i=0; i<3; i++){// Kt iterator
 505     for(Int_t j=0; j<10; j++){// Mbin iterator
 506       fNormWeight[i][j]=0x0;
 507     }
 508   }
 509
 510
 511 }
 512 //________________________________________________________________________
 513 AliFourPion &AliFourPion::operator=(const AliFourPion &obj)
 514 {
 515   // Assignment operator
 516   if (this == &obj)
 517     return *this;
 518
 519   fname = obj.fname;
 520   fAOD = obj.fAOD;
 521   fOutputList = obj.fOutputList;
 522   fPIDResponse = obj.fPIDResponse;
 523   fEC = obj.fEC;
 524   fEvt = obj.fEvt;
 525   fTempStruct = obj.fTempStruct;
 526   fRandomNumber = obj.fRandomNumber;
 527   fLEGO = fLEGO;
 528   fMCcase = obj.fMCcase;
 529   fAODcase = obj.fAODcase;
 530   fPbPbcase = obj.fPbPbcase;
 531   fGenerateSignal = obj.fGenerateSignal;
 532   fGeneratorOnly = obj.fGeneratorOnly;
 533   fPdensityPairCut = obj.fPdensityPairCut;
 534   fTabulatePairs = obj.fTabulatePairs;
 535   fRMax = obj.fRMax;
 536   ffcSq = obj.ffcSq;
 537   fFilterBit = obj.fFilterBit;
 538   fMaxChi2NDF = obj.fMaxChi2NDF;
 539   fMinTPCncls = obj.fMinTPCncls;
 540   fBfield = obj.fBfield;
 541   fMbin = obj.fMbin;
 542   fFSIindex = obj.fFSIindex;
 543   fEDbin = obj.fEDbin;
 544   fMbins = obj.fMbins;
 545   fMultLimit = obj.fMultLimit;
 546   fCentBinLowLimit = obj.fCentBinLowLimit;
 547   fCentBinHighLimit = obj.fCentBinHighLimit;
 548   fEventCounter = obj.fEventCounter;
 549   fEventsToMix = obj.fEventsToMix;
 550   fZvertexBins = obj.fZvertexBins;
 551   fQLowerCut = obj.fQLowerCut;
 552   fKupperBound = obj.fKupperBound;
 553   fQupperBoundQ2 = obj.fQupperBoundQ2;
 554   fQupperBoundQ3 = obj.fQupperBoundQ3;
 555   fQupperBoundQ4 = obj.fQupperBoundQ4;
 556   fQbinsQ2 = obj.fQbinsQ2;
 557   fQbinsQ3 = obj.fQbinsQ3;
 558   fQbinsQ4 = obj.fQbinsQ4;
 559   fQupperBoundWeights = obj.fQupperBoundWeights;
 560   fQstep = obj.fQstep;
 561   fQstepWeights = obj.fQstepWeights;
 562   fDampStart = obj.fDampStart;
 563   fDampStep = obj.fDampStep;
 564   fTPCTOFboundry = obj.fTPCTOFboundry;
 565   fTOFboundry = obj.fTOFboundry;
 566   fSigmaCutTPC = obj.fSigmaCutTPC;
 567   fSigmaCutTOF = obj.fSigmaCutTOF;
 568   fMinSepPairEta = obj.fMinSepPairEta;
 569   fMinSepPairPhi = obj.fMinSepPairPhi;
 570   fShareQuality = obj.fShareQuality;
 571   fShareFraction = obj.fShareFraction;
 572   fTrueMassP = obj.fTrueMassP;
 573   fTrueMassPi = obj.fTrueMassPi;
 574   fTrueMassK = obj.fTrueMassK;
 575   fTrueMassKs = obj.fTrueMassKs;
 576   fTrueMassLam = obj.fTrueMassLam;
 577   fKtIndexL = obj.fKtIndexL;
 578   fKtIndexH = obj.fKtIndexH;
 579   fQoIndexL = obj.fQoIndexL;
 580   fQoIndexH = obj.fQoIndexH;
 581   fQsIndexL = obj.fQsIndexL;
 582   fQsIndexH = obj.fQsIndexH;
 583   fQlIndexL = obj.fQlIndexL;
 584   fQlIndexH = obj.fQlIndexH;
 585   fDummyB = obj.fDummyB;
 586   fKT3transition = obj.fKT3transition;
 587   fKT4transition = obj.fKT4transition;
 588   fMomResC2 = obj.fMomResC2;
 589   fWeightmuonCorrection = obj.fWeightmuonCorrection;
 590
 591   for(Int_t i=0; i<12; i++){
 592     fFSIss[i]=obj.fFSIss[i];
 593     fFSIos[i]=obj.fFSIos[i];
 594   }
 595   for(Int_t i=0; i<3; i++){// Kt iterator
 596     for(Int_t j=0; j<10; j++){// Mbin iterator
 597       fNormWeight[i][j]=obj.fNormWeight[i][j];
 598     }
 599   }
 600
 601   return (*this);
 602 }
 603 //________________________________________________________________________
 604 AliFourPion::~AliFourPion()
 605 {
 606   // Destructor
 607   if(fAOD) delete fAOD;
 608   //if(fESD) delete fESD;
 609   if(fOutputList) delete fOutputList;
 610   if(fPIDResponse) delete fPIDResponse;
 611   if(fEC) delete fEC;
 612   if(fEvt) delete fEvt;
 613   if(fTempStruct) delete [] fTempStruct;
 614   if(fRandomNumber) delete fRandomNumber;
 615   if(fMomResC2) delete fMomResC2;
 616   if(fWeightmuonCorrection) delete fWeightmuonCorrection;
 617
 618   for(Int_t j=0; j<kMultLimitPbPb; j++){
 619     if(fLowQPairSwitch_E0E0[j]) delete [] fLowQPairSwitch_E0E0[j];
 620     if(fLowQPairSwitch_E0E1[j]) delete [] fLowQPairSwitch_E0E1[j];
 621     if(fLowQPairSwitch_E0E2[j]) delete [] fLowQPairSwitch_E0E2[j];
 622     if(fLowQPairSwitch_E0E3[j]) delete [] fLowQPairSwitch_E0E3[j];
 623     if(fLowQPairSwitch_E1E1[j]) delete [] fLowQPairSwitch_E1E1[j];
 624     if(fLowQPairSwitch_E1E2[j]) delete [] fLowQPairSwitch_E1E2[j];
 625     if(fLowQPairSwitch_E1E3[j]) delete [] fLowQPairSwitch_E1E3[j];
 626     if(fLowQPairSwitch_E2E3[j]) delete [] fLowQPairSwitch_E2E3[j];
 627     //
 628     if(fNormQPairSwitch_E0E0[j]) delete [] fNormQPairSwitch_E0E0[j];
 629     if(fNormQPairSwitch_E0E1[j]) delete [] fNormQPairSwitch_E0E1[j];
 630     if(fNormQPairSwitch_E0E2[j]) delete [] fNormQPairSwitch_E0E2[j];
 631     if(fNormQPairSwitch_E0E3[j]) delete [] fNormQPairSwitch_E0E3[j];
 632     if(fNormQPairSwitch_E1E1[j]) delete [] fNormQPairSwitch_E1E1[j];
 633     if(fNormQPairSwitch_E1E2[j]) delete [] fNormQPairSwitch_E1E2[j];
 634     if(fNormQPairSwitch_E1E3[j]) delete [] fNormQPairSwitch_E1E3[j];
 635     if(fNormQPairSwitch_E2E3[j]) delete [] fNormQPairSwitch_E2E3[j];
 636   }
 637
 638   //
 639   for(Int_t mb=0; mb<fMbins; mb++){
 640     for(Int_t edB=0; edB<fEDbins; edB++){
 641       for(Int_t c1=0; c1<2; c1++){
 642         for(Int_t c2=0; c2<2; c2++){
 643           for(Int_t term=0; term<2; term++){
 644
 645             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2;
 646
 647             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal;
 648             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared;
 649             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL;
 650             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW;
 651             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL;
 652             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW;
 653             //
 654             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
 655             if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
 656           }// term_2
 657
 658           for(Int_t c3=0; c3<2; c3++){
 659             for(Int_t term=0; term<5; term++){
 660
 661               if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3;
 662               if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3;
 663               if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor;
 664               //
 665               if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm;
 666
 667             }// term_3
 668
 669             for(Int_t c4=0; c4<2; c4++){
 670               for(Int_t term=0; term<13; term++){
 671
 672                 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4;
 673                 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4;
 674                 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor;
 675                 //
 676                 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm;
 677               }// term_4
 678
 679             }//c4
 680           }//c3
 681         }//c2
 682       }//c1
 683       for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
 684         for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
 685           if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD;
 686           if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD;
 687         }
 688       }
 689
 690     }// ED
 691   }// Mbin
 692
 693
 694   for(Int_t i=0; i<12; i++){
 695     if(fFSIss[i]) delete fFSIss[i];
 696     if(fFSIos[i]) delete fFSIos[i];
 697   }
 698   for(Int_t i=0; i<3; i++){// Kt iterator
 699     for(Int_t j=0; j<10; j++){// Mbin iterator
 700       if(fNormWeight[i][j]) delete fNormWeight[i][j];
 701     }
 702   }
 703
 704 }
 705 //________________________________________________________________________
 706 void AliFourPion::ParInit()
 707 {
 708   cout<<"AliFourPion MyInit() call"<<endl;
 709   cout<<"lego:"<<fLEGO<<"  MCcase:"<<fMCcase<<"  PbPbcase:"<<fPbPbcase<<"  TabulatePairs:"<<fTabulatePairs<<"  GenSignal:"<<fGenerateSignal<<"  CentLow:"<<fCentBinLowLimit<<"  CentHigh:"<<fCentBinHighLimit<<"  RMax:"<<fRMax<<"  fc^2:"<<ffcSq<<"  FB:"<<fFilterBit<<"  MaxChi2/NDF:"<<fMaxChi2NDF<<"  MinTPCncls:"<<fMinTPCncls<<"  MinPairSepEta:"<<fMinSepPairEta<<"  MinPairSepPhi:"<<fMinSepPairPhi<<"  NsigTPC:"<<fSigmaCutTPC<<"  NsigTOF:"<<fSigmaCutTOF<<endl;
 710
 711   fRandomNumber = new TRandom3();
 712   fRandomNumber->SetSeed(0);
 713
 714   //
 715   fEventCounter=0;
 716   fEventsToMix=3;
 717   fZvertexBins=2;//2
 718
 719   fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
 720   fTOFboundry = 2.1;// TOF pid used below this momentum
 721
 722   ////////////////////////////////////////////////
 723   // PadRow Pair Cuts
 724   fShareQuality = .5;// max
 725   fShareFraction = .05;// max
 726   ////////////////////////////////////////////////
 727
 728
 729   fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
 730   fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
 731
 732
 733
 734   if(fPbPbcase) {// PbPb
 735     fMultLimit=kMultLimitPbPb;
 736     fMbins=fCentBins;
 737     fQcut=0.1;
 738     fNormQcutLow = 0.15;// 0.15
 739     fNormQcutHigh = 0.2;// 0.175
 740   }else {// pp
 741     fMultLimit=kMultLimitpp;
 742     fMbins=kMultBinspp;
 743     fQcut=0.6;
 744     fNormQcutLow = 1.0;
 745     fNormQcutHigh = 1.5;
 746   }
 747
 748   fQLowerCut = 0.005;// was 0.005
 749   fKupperBound = 1.0;
 750   //
 751   fKstepY[0] = 1.6;
 752   fKmeanY[0] = 0;// central y
 753   fKmiddleY[0] = 0;
 754
 755   // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
 756   if(fKbinsT==4){
 757     fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
 758     fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
 759     fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
 760   }
 761   // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
 762   if(fKbinsT==3){
 763     fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
 764     fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
 765     fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
 766   }
 767   // 2x1 (Kt: 0-0.35, 0.35-1.0)
 768   if(fKbinsT==2){
 769     fKstepT[0] = 0.35; fKstepT[1] = 0.65;
 770     fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
 771     fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
 772   }
 773
 774   //
 775   fQupperBoundWeights = 0.2;
 776   fQupperBoundQ2 = 2.0;
 777   fQupperBoundQ3 = 0.6;
 778   fQupperBoundQ4 = 0.6;
 779   fQbinsQ2 = fQupperBoundQ2/0.005;
 780   fQbinsQ3 = fQupperBoundQ3/0.005;
 781   fQbinsQ4 = fQupperBoundQ4/0.005;
 782   fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
 783   for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
 784   //
 785   fDampStart = 0.5;// was 0.3, then 0.5
 786   fDampStep = 0.02;
 787
 788   //
 789   fKT3transition = 0.3;
 790   fKT4transition = 0.3;
 791
 792   fEC = new AliFourPionEventCollection **[fZvertexBins];
 793   for(UShort_t i=0; i<fZvertexBins; i++){
 794
 795     fEC[i] = new AliFourPionEventCollection *[fMbins];
 796
 797     for(UShort_t j=0; j<fMbins; j++){
 798
 799       fEC[i][j] = new AliFourPionEventCollection(fEventsToMix+1, fMultLimit, kMCarrayLimit, fMCcase);
 800     }
 801   }
 802
 803   for(Int_t i=0; i<kMultLimitPbPb; i++) fDefaultsCharSwitch[i]='0';
 804   for(Int_t i=0; i<kMultLimitPbPb; i++) {
 805     fLowQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 806     fLowQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 807     fLowQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 808     fLowQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 809     fLowQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 810     fLowQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 811     fLowQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 812     fLowQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 813     //
 814     fNormQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 815     fNormQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 816     fNormQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 817     fNormQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 818     fNormQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 819     fNormQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 820     fNormQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 821     fNormQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
 822   }
 823
 824   fTempStruct = new AliFourPionTrackStruct[fMultLimit];
 825
 826
 827   fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
 828
 829
 830
 831   // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
 832   if(!fLEGO) {
 833     SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
 834     if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
 835     if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
 836     if(!fMCcase && !fTabulatePairs) SetMuonCorrections(fLEGO);// Read Muon corrections
 837   }
 838
 839   /////////////////////////////////////////////
 840   /////////////////////////////////////////////
 841
 842 }
 843 //________________________________________________________________________
 844 void AliFourPion::UserCreateOutputObjects()
 845 {
 846   // Create histograms
 847   // Called once
 848
 849   ParInit();// Initialize my settings
 850
 851
 852   fOutputList = new TList();
 853   fOutputList->SetOwner();
 854
 855   TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
 856   fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
 857   fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
 858   fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
 859   fOutputList->Add(fVertexDist);
 860
 861
 862   TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
 863   fOutputList->Add(fDCAxyDistPlus);
 864   TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
 865   fOutputList->Add(fDCAzDistPlus);
 866   TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
 867   fOutputList->Add(fDCAxyDistMinus);
 868   TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
 869   fOutputList->Add(fDCAzDistMinus);
 870
 871
 872   TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
 873   fOutputList->Add(fEvents1);
 874   TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
 875   fOutputList->Add(fEvents2);
 876
 877   TH1F *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
 878   fMultDist0->GetXaxis()->SetTitle("Multiplicity");
 879   fOutputList->Add(fMultDist0);
 880
 881   TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
 882   fMultDist1->GetXaxis()->SetTitle("Multiplicity");
 883   fOutputList->Add(fMultDist1);
 884
 885   TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
 886   fMultDist2->GetXaxis()->SetTitle("Multiplicity");
 887   fOutputList->Add(fMultDist2);
 888
 889   TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
 890   fMultDist3->GetXaxis()->SetTitle("Multiplicity");
 891   fOutputList->Add(fMultDist3);
 892
 893   TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
 894   fOutputList->Add(fPtEtaDist);
 895
 896   TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
 897   fOutputList->Add(fPhiPtDist);
 898
 899   TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
 900   fOutputList->Add(fTOFResponse);
 901   TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
 902   fOutputList->Add(fTPCResponse);
 903
 904   TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
 905   fOutputList->Add(fRejectedPairs);
 906   TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
 907   fOutputList->Add(fRejectedEvents);
 908
 909   TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
 910   if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
 911   TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
 912   if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
 913   TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
 914   if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
 915   TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
 916   if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
 917   TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
 918   if(fMCcase) fOutputList->Add(fPairsPadRowNum);
 919   TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
 920   if(fMCcase) fOutputList->Add(fPairsPadRowDen);
 921
 922
 923
 924   TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
 925   if(fMCcase) fOutputList->Add(fResonanceOSPairs);
 926   TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
 927   if(fMCcase) fOutputList->Add(fAllOSPairs);
 928
 929   TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 930   if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
 931   TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 932   if(fMCcase) fOutputList->Add(fAllSCPionPairs);
 933   TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 934   if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
 935   TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 936   if(fMCcase) fOutputList->Add(fAllMCPionPairs);
 937
 938   //
 939   TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
 940   if(fMCcase) fOutputList->Add(fMuonParents);
 941   TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
 942   if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
 943   TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
 944   if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
 945   TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
 946   if(fMCcase) fOutputList->Add(fPionCandidates);
 947   //
 948
 949
 950   TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
 951   fOutputList->Add(fAvgMult);
 952
 953   TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
 954   fOutputList->Add(fTrackChi2NDF);
 955   TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
 956   fOutputList->Add(fTrackTPCncls);
 957
 958
 959   TH1D *fTPNRejects3pion1 = new TH1D("fTPNRejects3pion1","",fQbinsQ3,0,fQupperBoundQ3);
 960   fOutputList->Add(fTPNRejects3pion1);
 961   TH1D *fTPNRejects3pion2 = new TH1D("fTPNRejects3pion2","",fQbinsQ3,0,fQupperBoundQ3);
 962   fOutputList->Add(fTPNRejects3pion2);
 963   TH1D *fTPNRejects4pion1 = new TH1D("fTPNRejects4pion1","",fQbinsQ4,0,fQupperBoundQ4);
 964   fOutputList->Add(fTPNRejects4pion1);
 965
 966   TH3D *fKT3DistTerm1 = new TH3D("fKT3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 967   TH3D *fKT3DistTerm5 = new TH3D("fKT3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 968   fOutputList->Add(fKT3DistTerm1);
 969   fOutputList->Add(fKT3DistTerm5);
 970   TH3D *fKT4DistTerm1 = new TH3D("fKT4DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 971   TH3D *fKT4DistTerm13 = new TH3D("fKT4DistTerm13","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
 972   fOutputList->Add(fKT4DistTerm1);
 973   fOutputList->Add(fKT4DistTerm13);
 974
 975
 976   TProfile2D *fKT3AvgpT = new TProfile2D("fKT3AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
 977   fOutputList->Add(fKT3AvgpT);
 978   TProfile2D *fKT4AvgpT = new TProfile2D("fKT4AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
 979   fOutputList->Add(fKT4AvgpT);
 980
 981
 982   TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
 983   TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
 984   TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
 985   TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
 986   TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
 987   TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
 988   TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
 989   TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
 990   TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
 991   TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
 992   TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
 993   TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
 994   fOutputList->Add(fMCWeight3DTerm1SC);
 995   fOutputList->Add(fMCWeight3DTerm1SCden);
 996   fOutputList->Add(fMCWeight3DTerm2SC);
 997   fOutputList->Add(fMCWeight3DTerm2SCden);
 998   fOutputList->Add(fMCWeight3DTerm1MC);
 999   fOutputList->Add(fMCWeight3DTerm1MCden);
1000   fOutputList->Add(fMCWeight3DTerm2MC);
1001   fOutputList->Add(fMCWeight3DTerm2MCden);
1002   fOutputList->Add(fMCWeight3DTerm3MC);
1003   fOutputList->Add(fMCWeight3DTerm3MCden);
1004   fOutputList->Add(fMCWeight3DTerm4MC);
1005   fOutputList->Add(fMCWeight3DTerm4MCden);
1006
1007
1008   if(fPdensityPairCut){
1009
1010     for(Int_t mb=0; mb<fMbins; mb++){
1011       if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
1012
1013       for(Int_t edB=0; edB<fEDbins; edB++){
1014         for(Int_t c1=0; c1<2; c1++){
1015           for(Int_t c2=0; c2<2; c2++){
1016             for(Int_t term=0; term<2; term++){
1017
1018               TString *nameEx2 = new TString("TwoParticle_Charge1_");
1019               *nameEx2 += c1;
1020               nameEx2->Append("_Charge2_");
1021               *nameEx2 += c2;
1022               nameEx2->Append("_M_");
1023               *nameEx2 += mb;
1024               nameEx2->Append("_ED_");
1025               *nameEx2 += edB;
1026               nameEx2->Append("_Term_");
1027               *nameEx2 += term+1;
1028
1029               if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
1030
1031
1032               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1033               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2);
1034               TString *nameEx2QW=new TString(nameEx2->Data());
1035               nameEx2QW->Append("_QW");
1036               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1037               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW);
1038               TString *nameAvgP=new TString(nameEx2->Data());
1039               nameAvgP->Append("_AvgP");
1040               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, fQbinsQ2,0,fQupperBoundQ2, 0.,1.0,"");
1041               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP);
1042
1043               TString *nameUnitMult=new TString(nameEx2->Data());
1044               nameUnitMult->Append("_UnitMult");
1045               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin = new TH2D(nameUnitMult->Data(),"Two Particle Distribution",21,0.5,21.5, fQbinsQ2,0,fQupperBoundQ2);
1046               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin);
1047
1048               if(fMCcase){
1049                 // Momentum resolution histos
1050                 TString *nameIdeal = new TString(nameEx2->Data());
1051                 nameIdeal->Append("_Ideal");
1052                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1053                 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal);
1054                 TString *nameSmeared = new TString(nameEx2->Data());
1055                 nameSmeared->Append("_Smeared");
1056                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1057                 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared);
1058                 //
1059                 // Muon correction histos
1060                 TString *nameMuonIdeal=new TString(nameEx2->Data());
1061                 nameMuonIdeal->Append("_MuonIdeal");
1062                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1063                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal);
1064                 TString *nameMuonSmeared=new TString(nameEx2->Data());
1065                 nameMuonSmeared->Append("_MuonSmeared");
1066                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1067                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared);
1068                 //
1069                 TString *nameMuonPionK2=new TString(nameEx2->Data());
1070                 nameMuonPionK2->Append("_MuonPionK2");
1071                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2 = new TH2D(nameMuonPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1072                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2);
1073                 //
1074                 TString *namePionPionK2=new TString(nameEx2->Data());
1075                 namePionPionK2->Append("_PionPionK2");
1076                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2 = new TH2D(namePionPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1077                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2);
1078                 //
1079                 //
1080                 TString *nameEx2MC=new TString(nameEx2->Data());
1081                 nameEx2MC->Append("_MCqinv");
1082                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1083                 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
1084                 TString *nameEx2MCQW=new TString(nameEx2->Data());
1085                 nameEx2MCQW->Append("_MCqinvQW");
1086                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1087                 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
1088                 //
1089                 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
1090                 nameEx2PIDpurityDen->Append("_PIDpurityDen");
1091                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1092                 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
1093                 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
1094                 nameEx2PIDpurityNum->Append("_PIDpurityNum");
1095                 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH3D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",16,0.5,16.5, 20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1096                 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
1097               }
1098               TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1099               nameEx2OSLB1->Append("_osl_b1");
1100               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1101               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL);
1102               nameEx2OSLB1->Append("_QW");
1103               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1104               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW);
1105               //
1106               TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1107               nameEx2OSLB2->Append("_osl_b2");
1108               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1109               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL);
1110               nameEx2OSLB2->Append("_QW");
1111               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1112               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW);
1113
1114             }// term_2
1115
1116
1117
1118             // skip 3-particle if Tabulate6DPairs is true
1119             if(fTabulatePairs) continue;
1120
1121             for(Int_t c3=0; c3<2; c3++){
1122               for(Int_t term=0; term<5; term++){
1123
1124                 TString *namePC3 = new TString("ThreeParticle_Charge1_");
1125                 *namePC3 += c1;
1126                 namePC3->Append("_Charge2_");
1127                 *namePC3 += c2;
1128                 namePC3->Append("_Charge3_");
1129                 *namePC3 += c3;
1130                 namePC3->Append("_M_");
1131                 *namePC3 += mb;
1132                 namePC3->Append("_ED_");
1133                 *namePC3 += edB;
1134                 namePC3->Append("_Term_");
1135                 *namePC3 += term+1;
1136
1137                 ///////////////////////////////////////
1138                 // skip degenerate histograms
1139                 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1140                 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1141                 /////////////////////////////////////////
1142
1143
1144                 TString *nameNorm=new TString(namePC3->Data());
1145                 nameNorm->Append("_Norm");
1146                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1147                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1148                 //
1149
1150                 TString *name1DQ=new TString(namePC3->Data());
1151                 name1DQ->Append("_1D");
1152                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH1D(name1DQ->Data(),"", fQbinsQ3,0,fQupperBoundQ3);
1153                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1154                 //
1155                 TString *nameKfactor=new TString(namePC3->Data());
1156                 nameKfactor->Append("_Kfactor");
1157                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1158                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor);
1159                 //
1160                 TString *nameMeanQinv=new TString(namePC3->Data());
1161                 nameMeanQinv->Append("_MeanQinv");
1162                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv = new TProfile(nameMeanQinv->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,.2, "");
1163                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv);
1164
1165                 if(fMCcase==kTRUE){
1166                   // Momentum resolution correction histos
1167                   TString *nameMomResIdeal=new TString(namePC3->Data());
1168                   nameMomResIdeal->Append("_Ideal");
1169                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1170                   if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1171                   TString *nameMomResSmeared=new TString(namePC3->Data());
1172                   nameMomResSmeared->Append("_Smeared");
1173                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1174                   if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared);
1175                   // Muon correction histos
1176                   TString *nameMuonIdeal=new TString(namePC3->Data());
1177                   nameMuonIdeal->Append("_MuonIdeal");
1178                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1179                   if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal);
1180                   TString *nameMuonSmeared=new TString(namePC3->Data());
1181                   nameMuonSmeared->Append("_MuonSmeared");
1182                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1183                   if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared);
1184                   //
1185                   TString *nameMuonPionK3=new TString(namePC3->Data());
1186                   nameMuonPionK3->Append("_MuonPionK3");
1187                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH3D(nameMuonPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1188                   if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3);
1189                   //
1190                   TString *namePionPionK3=new TString(namePC3->Data());
1191                   namePionPionK3->Append("_PionPionK3");
1192                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH3D(namePionPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1193                   if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3);
1194
1195                 }// MCcase
1196                 //
1197                 if(c1==c2 && c1==c3 && term==4 ){
1198                   TString *nameTwoPartNorm=new TString(namePC3->Data());
1199                   nameTwoPartNorm->Append("_TwoPartNorm");
1200                   Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1201                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm);
1202                 }// term=4
1203
1204               }// term_3
1205
1206               for(Int_t c4=0; c4<2; c4++){
1207                 for(Int_t term=0; term<13; term++){
1208
1209                   TString *namePC4 = new TString("FourParticle_Charge1_");
1210                   *namePC4 += c1;
1211                   namePC4->Append("_Charge2_");
1212                   *namePC4 += c2;
1213                   namePC4->Append("_Charge3_");
1214                   *namePC4 += c3;
1215                   namePC4->Append("_Charge4_");
1216                   *namePC4 += c4;
1217                   namePC4->Append("_M_");
1218                   *namePC4 += mb;
1219                   namePC4->Append("_ED_");
1220                   *namePC4 += edB;
1221                   namePC4->Append("_Term_");
1222                   *namePC4 += term+1;
1223
1224                   ///////////////////////////////////////
1225                   // skip degenerate histograms
1226                   if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
1227                   if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
1228                   if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
1229                   /////////////////////////////////////////
1230
1231                   TString *nameNorm=new TString(namePC4->Data());
1232                   nameNorm->Append("_Norm");
1233                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1234                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
1235                   //
1236                   TString *name1DQ=new TString(namePC4->Data());
1237                   name1DQ->Append("_1D");
1238                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
1239                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
1240                   //
1241                   TString *nameKfactor=new TString(namePC4->Data());
1242                   nameKfactor->Append("_Kfactor");
1243                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1244                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
1245                   //
1246                   if(c1==c2 && c1==c3 && c1==c4 && term==12 ){
1247                     TString *nameTwoPartNorm=new TString(namePC4->Data());
1248                     nameTwoPartNorm->Append("_TwoPartNorm");
1249                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1250                     fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
1251                   }
1252
1253                   if(fMCcase==kTRUE){
1254                     // Momentum resolution correction histos
1255                     TString *nameMomResIdeal=new TString(namePC4->Data());
1256                     nameMomResIdeal->Append("_Ideal");
1257                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1258                     if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
1259                     TString *nameMomResSmeared=new TString(namePC4->Data());
1260                     nameMomResSmeared->Append("_Smeared");
1261                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1262                     if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
1263                     // Muon correction histos
1264                     TString *nameMuonIdeal=new TString(namePC4->Data());
1265                     nameMuonIdeal->Append("_MuonIdeal");
1266                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1267                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
1268                     TString *nameMuonSmeared=new TString(namePC4->Data());
1269                     nameMuonSmeared->Append("_MuonSmeared");
1270                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1271                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared);
1272                     //
1273                     TString *nameMuonPionK4=new TString(namePC4->Data());
1274                     nameMuonPionK4->Append("_MuonPionK4");
1275                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH3D(nameMuonPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1276                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4);
1277                     //
1278                     TString *namePionPionK4=new TString(namePC4->Data());
1279                     namePionPionK4->Append("_PionPionK4");
1280                     Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH3D(namePionPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1281                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4);
1282
1283                   }// MCcase
1284
1285
1286                 }
1287               }
1288
1289             }//c3
1290           }//c2
1291         }//c1
1292       }// ED
1293     }// mbin
1294   }// Pdensity Method
1295
1296
1297   if(fTabulatePairs){
1298
1299     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1300       for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1301         for(Int_t mb=0; mb<fMbins; mb++){
1302           for(Int_t edB=0; edB<fEDbins; edB++){
1303
1304             TString *nameNum = new TString("TPN_num_Kt_");
1305             *nameNum += tKbin;
1306             nameNum->Append("_Ky_");
1307             *nameNum += yKbin;
1308             nameNum->Append("_M_");
1309             *nameNum += mb;
1310             nameNum->Append("_ED_");
1311             *nameNum += edB;
1312
1313             TString *nameDen = new TString("TPN_den_Kt_");
1314             *nameDen += tKbin;
1315             nameDen->Append("_Ky_");
1316             *nameDen += yKbin;
1317             nameDen->Append("_M_");
1318             *nameDen += mb;
1319             nameDen->Append("_ED_");
1320             *nameDen += edB;
1321
1322             if(edB==0){
1323               KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1324               fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
1325
1326               KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1327               fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
1328             }
1329
1330           }
1331         }
1332       }
1333     }
1334
1335   }
1336
1337
1338   TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1339   fOutputList->Add(fQsmearMean);
1340   TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1341   fOutputList->Add(fQsmearSq);
1342   TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
1343   fOutputList->Add(fQ2Res);
1344   TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
1345   fOutputList->Add(fQ3Res);
1346   TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
1347   fOutputList->Add(fQ4Res);
1348
1349   TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
1350   fOutputList->Add(DistQinv4pion);
1351   TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
1352   if(fMCcase) fOutputList->Add(DistQinvMC4pion);
1353
1354   TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
1355   fOutputList->Add(fAvgQ12VersusQ3);
1356   TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
1357   fOutputList->Add(fAvgQ13VersusQ3);
1358   TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
1359   fOutputList->Add(fAvgQ23VersusQ3);
1360
1361   TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5);
1362   fOutputList->Add(fDistPionParents4);
1363
1364   ////////////////////////////////////
1365   ///////////////////////////////////
1366
1367   PostData(1, fOutputList);
1368 }
1369
1370 //________________________________________________________________________
1371 void AliFourPion::UserExec(Option_t *)
1372 {
1373   // Main loop
1374   // Called for each event
1375   //cout<<"===========  Event # "<<fEventCounter+1<<"  ==========="<<endl;
1376   fEventCounter++;
1377
1378   if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1379
1380   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1381   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1382
1383
1384   // Trigger Cut
1385   if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1386   Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1387   if(!isSelected1 && !fMCcase) {return;}
1388   }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1389     Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1390     Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1391     if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1392   }else {return;}
1393
1394   ///////////////////////////////////////////////////////////
1395   const AliAODVertex *primaryVertexAOD;
1396   AliCentrality *centrality;// for AODs and ESDs
1397
1398
1399   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1400   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1401   fPIDResponse = inputHandler->GetPIDResponse();
1402
1403
1404   TClonesArray *mcArray = 0x0;
1405   if(fMCcase){
1406     if(fAODcase){
1407       mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1408       if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1409         cout<<"No MC particle branch found or Array too large!!"<<endl;
1410         return;
1411       }
1412     }
1413   }
1414
1415
1416   UInt_t status=0;
1417   Int_t positiveTracks=0, negativeTracks=0;
1418   Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1419
1420   Double_t vertex[3]={0};
1421   Int_t zbin=0;
1422   Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1423   /////////////////////////////////////////////////
1424
1425
1426   Float_t centralityPercentile=0;
1427   Float_t cStep=5.0, cStart=0;
1428
1429
1430   if(fAODcase){// AOD case
1431
1432     if(fPbPbcase){
1433       centrality = fAOD->GetCentrality();
1434       centralityPercentile = centrality->GetCentralityPercentile("V0M");
1435       if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1436       if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range.  Skipping Event"<<endl;*/ return;}
1437       cout<<"Centrality % = "<<centralityPercentile<<endl;
1438     }
1439
1440
1441     ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
1442
1443     // Pile-up rejection
1444     AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1445     if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1446     else AnaUtil->SetUseMVPlpSelection(kFALSE);
1447     Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1448     if(pileUpCase) return;
1449
1450     ////////////////////////////////
1451     // Vertexing
1452     ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1453     primaryVertexAOD = fAOD->GetPrimaryVertex();
1454     vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1455
1456     if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1457     ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1458
1459     if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1460
1461     ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1462
1463     fBfield = fAOD->GetMagneticField();
1464
1465     for(Int_t i=0; i<fZvertexBins; i++){
1466       if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1467         zbin=i;
1468         break;
1469       }
1470     }
1471
1472
1473
1474     /////////////////////////////
1475     // Create Shuffled index list
1476     Int_t randomIndex[fAOD->GetNumberOfTracks()];
1477     for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1478     Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1479     /////////////////////////////
1480
1481
1482     // Track loop
1483     for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1484       AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1485       if (!aodtrack) continue;
1486       if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1487
1488       status=aodtrack->GetStatus();
1489
1490       if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1491       ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1492       ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1493       if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
1494       if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1495
1496       if(fFilterBit != 7){
1497         Bool_t goodTrackOtherFB = kFALSE;
1498         for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1499           AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1500           if(!aodtrack2) continue;
1501           if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1502
1503           if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1504
1505         }
1506         if(!goodTrackOtherFB) continue;
1507       }
1508
1509
1510       if(aodtrack->Pt() < 0.16) continue;
1511       if(fabs(aodtrack->Eta()) > 0.8) continue;
1512
1513
1514       Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1515       if(!goodMomentum) continue;
1516       aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1517
1518
1519       Double_t dca2[2]={0};
1520       dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1521       dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1522       Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1523
1524       fTempStruct[myTracks].fStatus = status;
1525       fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1526       fTempStruct[myTracks].fId = aodtrack->GetID();
1527       //
1528       fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1529       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1530       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1531       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1532       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1533       fTempStruct[myTracks].fEta = aodtrack->Eta();
1534       fTempStruct[myTracks].fCharge = aodtrack->Charge();
1535       fTempStruct[myTracks].fDCAXY = dca2[0];
1536       fTempStruct[myTracks].fDCAZ = dca2[1];
1537       fTempStruct[myTracks].fDCA = dca3d;
1538       fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1539       fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1540
1541
1542
1543       if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1544       //if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1545       //if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1546
1547
1548       // PID section
1549       fTempStruct[myTracks].fElectron = kFALSE;
1550       fTempStruct[myTracks].fPion = kFALSE;
1551       fTempStruct[myTracks].fKaon = kFALSE;
1552       fTempStruct[myTracks].fProton = kFALSE;
1553
1554       Float_t nSigmaTPC[5];
1555       Float_t nSigmaTOF[5];
1556       nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1557       nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1558       fTempStruct[myTracks].fTOFhit = kFALSE;// default
1559       Float_t signalTPC=0, signalTOF=0;
1560       Double_t integratedTimesTOF[10]={0};
1561
1562
1563       Bool_t DoPIDWorkAround=kTRUE;
1564       //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1565       if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1566       if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1567         nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1568         nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1569         nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1570         nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1571         nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1572         //
1573         nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1574         nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1575         nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1576         nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1577         nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1578         signalTPC = aodtrack->GetTPCsignal();
1579         if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1580           fTempStruct[myTracks].fTOFhit = kTRUE;
1581           signalTOF = aodtrack->GetTOFsignal();
1582           aodtrack->GetIntegratedTimes(integratedTimesTOF);
1583         }else fTempStruct[myTracks].fTOFhit = kFALSE;
1584
1585       }else {// FilterBit 7 PID workaround
1586
1587         for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1588           AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1589           if (!aodTrack2) continue;
1590           if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1591
1592           UInt_t status2=aodTrack2->GetStatus();
1593
1594           nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1595           nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1596           nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1597           nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1598           nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1599           //
1600           nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1601           nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1602           nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1603           nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1604           nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1605           signalTPC = aodTrack2->GetTPCsignal();
1606
1607           if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1608             fTempStruct[myTracks].fTOFhit = kTRUE;
1609             signalTOF = aodTrack2->GetTOFsignal();
1610             aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1611           }else fTempStruct[myTracks].fTOFhit = kFALSE;
1612
1613         }// aodTrack2
1614       }// FilterBit 7 PID workaround
1615
1616
1617       ///////////////////
1618       ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1619       if(fTempStruct[myTracks].fTOFhit) {
1620         ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1621       }
1622       ///////////////////
1623
1624       // Use TOF if good hit and above threshold
1625       if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1626         if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1627         if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1628         if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1629         if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1630       }else {// TPC info instead
1631         if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1632         if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1633         if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1634         if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1635       }
1636
1637
1638       // Ensure there is only 1 candidate per track
1639       if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1640       if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1641       if(!fTempStruct[myTracks].fPion) continue;// only pions
1642       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1643       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1644       if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1645       //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1646       ////////////////////////
1647       //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1648
1649
1650
1651       if(fTempStruct[myTracks].fCharge==+1) {
1652         ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1653         ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1654       }else {
1655         ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1656         ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1657       }
1658
1659       ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1660       ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1661
1662
1663
1664       if(fTempStruct[myTracks].fPion) {// pions
1665         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1666         fTempStruct[myTracks].fKey = 1;
1667       }else if(fTempStruct[myTracks].fKaon){// kaons
1668         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1669         fTempStruct[myTracks].fKey = 10;
1670       }else{// protons
1671         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1672         fTempStruct[myTracks].fKey = 100;
1673       }
1674
1675
1676
1677       if(aodtrack->Charge() > 0) positiveTracks++;
1678       else negativeTracks++;
1679
1680       if(fTempStruct[myTracks].fPion) pionCount++;
1681       if(fTempStruct[myTracks].fKaon) kaonCount++;
1682       if(fTempStruct[myTracks].fProton) protonCount++;
1683
1684       myTracks++;
1685
1686       if(fMCcase){// muon mothers
1687         AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1688         if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1689           AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1690           if(parent->IsPhysicalPrimary()){
1691             ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1692           }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1693         }
1694         ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1695       }
1696     }
1697     //cout<<"kinkcount = "<<kinkcount<<"   pionkinks = "<<pionkinks<<"   primarypionkinks = "<<primarypionkinks<<endl;
1698   }else {// ESD tracks
1699     cout<<"ESDs not supported currently"<<endl;
1700     return;
1701   }
1702
1703   // Generator info only
1704   if(fMCcase && fGeneratorOnly){
1705     myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1706     for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1707       if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1708       if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1709
1710       AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1711       if(!mcParticle) continue;
1712       if(fabs(mcParticle->Eta())>0.8) continue;
1713       if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
1714       if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
1715       if(!mcParticle->IsPrimary()) continue;
1716       if(!mcParticle->IsPhysicalPrimary()) continue;
1717       if(abs(mcParticle->GetPdgCode())!=211) continue;
1718
1719       fTempStruct[myTracks].fP[0] = mcParticle->Px();
1720       fTempStruct[myTracks].fP[1] = mcParticle->Py();
1721       fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1722       fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1723
1724       fTempStruct[myTracks].fId = myTracks;// use my track counter
1725       fTempStruct[myTracks].fLabel = mctrackN;
1726       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1727       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1728       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1729       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1730       fTempStruct[myTracks].fEta = mcParticle->Eta();
1731       fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1732       fTempStruct[myTracks].fDCAXY = 0.;
1733       fTempStruct[myTracks].fDCAZ = 0.;
1734       fTempStruct[myTracks].fDCA = 0.;
1735       fTempStruct[myTracks].fPion = kTRUE;
1736       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1737       fTempStruct[myTracks].fKey = 1;
1738
1739       myTracks++;
1740       pionCount++;
1741     }
1742   }
1743
1744   if(myTracks >= 1) {
1745     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1746   }
1747
1748
1749   //cout<<"There are "<<myTracks<<"  myTracks"<<endl;
1750   //cout<<"pionCount = "<<pionCount<<"   kaonCount = "<<kaonCount<<"   protonCount = "<<protonCount<<endl;
1751
1752
1753   /////////////////////////////////////////
1754   // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1755   if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1756   /////////////////////////////////////////
1757
1758
1759   ////////////////////////////////
1760   ///////////////////////////////
1761   // Mbin determination
1762   //
1763   // Mbin set to Pion Count Only for pp!!!!!!!
1764   fMbin=-1;
1765   if(!fPbPbcase){
1766     for(Int_t i=0; i<kMultBinspp; i++){
1767       if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1768       // Mbin 0 has 1 pion
1769     }
1770   }else{
1771     for(Int_t i=0; i<fCentBins; i++){
1772       if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1773         fMbin=i;// 0 = most central
1774         break;
1775       }
1776     }
1777   }
1778
1779   if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1780
1781   ///////////////////
1782   // can only be called after fMbin has been set
1783   // Radius parameter only matters for Monte-Carlo data
1784   SetFSIindex(fRMax);
1785   ///////////////////
1786
1787   Int_t rBinForTPNMomRes = 10;
1788   if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1789   else if(fMbin==1) {rBinForTPNMomRes=9;}
1790   else if(fMbin<=3) {rBinForTPNMomRes=8;}
1791   else if(fMbin<=5) {rBinForTPNMomRes=7;}
1792   else {rBinForTPNMomRes=6;}
1793
1794   //////////////////////////////////////////////////
1795   fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1796   //////////////////////////////////////////////////
1797
1798
1799
1800   ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1801   ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1802
1803   ////////////////////////////////////
1804   // Add event to buffer if > 0 tracks
1805   if(myTracks > 0){
1806     fEC[zbin][fMbin]->FIFOShift();
1807     (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1808     (fEvt)->fNtracks = myTracks;
1809     (fEvt)->fFillStatus = 1;
1810     for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1811     if(fMCcase){
1812       (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1813       for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1814         AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1815         (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1816         (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1817         (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1818         (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1819         (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1820         (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1821       }
1822     }
1823   }
1824
1825
1826
1827   Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1828   Float_t qout=0, qside=0, qlong=0;
1829   Float_t kT12=0;
1830   Float_t q3=0, q3MC=0;
1831   Float_t q4=0, q4MC=0;
1832   Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1833   Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1834   Float_t pVect1[4]={0};
1835   Float_t pVect2[4]={0};
1836   Float_t pVect3[4]={0};
1837   Float_t pVect4[4]={0};
1838   Float_t pVect1MC[4]={0};
1839   Float_t pVect2MC[4]={0};
1840   Float_t pVect3MC[4]={0};
1841   Float_t pVect4MC[4]={0};
1842   Float_t Pparent1[4]={0};
1843   Float_t Pparent2[4]={0};
1844   Float_t Pparent3[4]={0};
1845   Float_t Pparent4[4]={0};
1846   Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1847   Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
1848   Float_t weight12CC[3]={0};
1849   Float_t weight13CC[3]={0};
1850   Float_t weight14CC[3]={0};
1851   Float_t weight23CC[3]={0};
1852   Float_t weight24CC[3]={0};
1853   Float_t weight34CC[3]={0};
1854   Float_t weightTotal=0;//, weightTotalErr=0;
1855   Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0;
1856   Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1857   Float_t parentQ3=0;
1858   Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1859   Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1860   Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
1861   Bool_t GoodTripletWeights=kFALSE;
1862   //
1863   AliAODMCParticle *mcParticle1=0x0;
1864   AliAODMCParticle *mcParticle2=0x0;
1865
1866
1867   ////////////////////
1868   //Int_t PairCount[7]={0};
1869   //Int_t NormPairCount[7]={0};
1870   Int_t KT3index=0, KT4index=0;
1871
1872   // reset to defaults
1873   for(Int_t i=0; i<kMultLimitPbPb; i++) {
1874     fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1875     fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1876     fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1877     fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1878     fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1879     fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1880     fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1881     fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1882     //
1883     fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1884     fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1885     fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1886     fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1887     fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1888     fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1889     fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1890     fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1891   }
1892
1893
1894   //////////////////////////////////////////
1895   // make low-q pair storage and normalization-pair storage
1896   //
1897   for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
1898     for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
1899       if(en1>1 && en1==en2) continue;
1900
1901       for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
1902         for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
1903
1904
1905           pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1906           pVect1[1]=(fEvt+en1)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1907           pVect1[2]=(fEvt+en1)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1908           pVect1[3]=(fEvt+en1)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1909           ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
1910           ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1911
1912           qinv12 = GetQinv(pVect1, pVect2);
1913           kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1914           SetFillBins2(ch1, ch2, bin1, bin2);
1915
1916           if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1917           if(ch1 == ch2 && !fGeneratorOnly){
1918             if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1919               if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1920               continue;
1921             }
1922           }
1923
1924           GetQosl(pVect1, pVect2, qout, qside, qlong);
1925           if( (en1+en2==0)) {
1926             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
1927             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
1928             // osl frame
1929             if((kT12 > 0.2) && (kT12 < 0.3)){
1930               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1931               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1932             }
1933             if((kT12 > 0.6) && (kT12 < 0.7)){
1934               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1935               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1936             }
1937             // unit mult bins
1938             if( (fEvt+en1)->fNtracks%100==0){
1939               Int_t kTindex=0;
1940               if(kT12>0.3) kTindex=1;
1941               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1942               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
1943             }
1944
1945           }
1946           if( (en1+en2==1)) {
1947             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
1948             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
1949             // osl frame
1950             if((kT12 > 0.2) && (kT12 < 0.3)){
1951               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1952               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1953             }
1954             if((kT12 > 0.6) && (kT12 < 0.7)){
1955               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1956               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1957             }
1958             // unit mult bins
1959             if( (fEvt+en1)->fNtracks%100==0){
1960               Int_t kTindex=0;
1961               if(kT12>0.3) kTindex=1;
1962               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1963               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
1964             }
1965           }
1966           //////////////////////////////////////////
1967           if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
1968             Float_t kY = 0;
1969             Int_t kTbin=-1, kYbin=-1;
1970
1971             for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}}
1972             for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
1973             if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1974             if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1975             if(fGenerateSignal && en2==0) {
1976               Int_t chGroup2[2]={ch1,ch2};
1977               Float_t WInput = MCWeight(chGroup2, fRMax, 0.7, qinv12, kT12);
1978               KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
1979             }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1980           }
1981
1982           //////////////////////////////////////////////////////////////////////////////
1983
1984           if(qinv12 <= fQcut) {
1985             if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
1986             if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
1987             if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
1988             if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
1989             if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
1990             if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
1991             if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
1992             if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
1993           }
1994           if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
1995             if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
1996             if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
1997             if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
1998             if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
1999             if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
2000             if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2001             if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2002             if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
2003           }
2004
2005         }
2006       }
2007     }
2008   }
2009
2010   //cout<<PairCount[0]<<"  "<<PairCount[1]<<"  "<<PairCount[2]<<"  "<<PairCount[3]<<"  "<<PairCount[4]<<"  "<<PairCount[5]<<"  "<<PairCount[6]<<endl;
2011   //cout<<NormPairCount[0]<<"  "<<NormPairCount[1]<<"  "<<NormPairCount[2]<<"  "<<NormPairCount[3]<<"  "<<NormPairCount[4]<<"  "<<NormPairCount[5]<<"  "<<NormPairCount[6]<<endl;
2012   ///////////////////////////////////////////////////
2013   // Do not use pairs from events with too many pairs
2014
2015   ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2016
2017   ///////////////////////////////////////////////////
2018
2019
2020   if(fTabulatePairs) return;
2021
2022
2023
2024   ////////////////////////////////////////////////////
2025   ////////////////////////////////////////////////////
2026   // Normalization counting of 3- and 4-particle terms
2027   for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2028     for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2029       if(en2==0 && en3>2) continue;// not needed config
2030       if(en2==1 && en3==en2) continue;// not needed config
2031       for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2032         if(en3==0 && en4>1) continue;// not needed config
2033         if(en3==1 && en4==3) continue;// not needed configs
2034         if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2035
2036         for (Int_t i=0; i<myTracks; i++) {// 1st particle
2037           pVect1[1]=(fEvt)->fTracks[i].fP[0];
2038           pVect1[2]=(fEvt)->fTracks[i].fP[1];
2039           pVect1[3]=(fEvt)->fTracks[i].fP[2];
2040           ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2041
2042           for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2043             if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2044             else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2045
2046             pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2047             pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2048             pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2049             ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2050
2051             for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2052               if(en3==0) {
2053                 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2054                 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2055               }else if(en3==1){
2056                 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2057                 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2058               }else{
2059                 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2060                 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2061               }
2062
2063               pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2064               pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2065               pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2066               ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2067               Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2068               SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2069
2070               Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2071               if(KT3<=fKT3transition) KT3index=0;
2072               else KT3index=1;
2073
2074               if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2075               if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2076               if(en2==0 && en3==1 && en4==2) {
2077                 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2078                 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2079                 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2080               }
2081
2082
2083               for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2084                 if(en4==0){
2085                   if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2086                   if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2087                   if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2088                 }else if(en4==1){
2089                   if(en3==0){
2090                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2091                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2092                     if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2093                   }else{
2094                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2095                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2096                     if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2097                   }
2098                 }else if(en4==2){
2099                   if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2100                   if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2101                   if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
2102                 }else{
2103                   if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2104                   if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2105                   if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
2106                 }
2107
2108                 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2109                 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2110                 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2111                 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2112                 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2113                 if(KT4<=fKT4transition) KT4index=0;
2114                 else KT4index=1;
2115
2116                 Bool_t FillTerms[13]={kFALSE};
2117                 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2118                 //
2119                 for(int ft=0; ft<13; ft++) {
2120                   if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.);
2121                 }
2122
2123
2124               }
2125             }
2126           }
2127         }
2128
2129       }
2130     }
2131   }
2132
2133
2134
2135
2136     ///////////////////////////////////////////////////////////////////////
2137     ///////////////////////////////////////////////////////////////////////
2138     ///////////////////////////////////////////////////////////////////////
2139     //
2140     //
2141     // Start the Main Correlation Analysis
2142     //
2143     //
2144     ///////////////////////////////////////////////////////////////////////
2145
2146
2147
2148     ////////////////////////////////////////////////////
2149     ////////////////////////////////////////////////////
2150     for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2151       for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2152         if(en2==0 && en3>2) continue;// not needed config
2153         if(en2==1 && en3==en2) continue;// not needed config
2154         for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2155           if(en3==0 && en4>1) continue;// not needed config
2156           if(en3==1 && en4==3) continue;// not needed configs
2157           if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2158
2159           Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
2160
2161           /////////////////////////////////////////////////////////////
2162           for (Int_t i=0; i<myTracks; i++) {// 1st particle
2163             pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2164             pVect1[1]=(fEvt)->fTracks[i].fP[0];
2165             pVect1[2]=(fEvt)->fTracks[i].fP[1];
2166             pVect1[3]=(fEvt)->fTracks[i].fP[2];
2167             ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2168
2169             /////////////////////////////////////////////////////////////
2170             for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2171               if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2172               else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2173
2174               pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2175               pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2176               pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2177               pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2178               ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2179               qinv12 = GetQinv(pVect1, pVect2);
2180               kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2181               SetFillBins2(ch1, ch2, bin1, bin2);
2182               Int_t kTindex=0;
2183               if(kT12<=0.3) kTindex=0;
2184               else kTindex=1;
2185
2186               FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2187
2188               // two particle terms filled during tabulation of low-q pairs
2189
2190
2191               if(fMCcase){
2192                 FilledMCpair12=kFALSE;
2193
2194                 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
2195                   for(Int_t rstep=0; rstep<10; rstep++){
2196                     Float_t coeff = (rstep)*0.2*(0.18/1.2);
2197                     Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2198                     if(phi1 > 2*PI) phi1 -= 2*PI;
2199                     if(phi1 < 0) phi1 += 2*PI;
2200                     Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2201                     if(phi2 > 2*PI) phi2 -= 2*PI;
2202                     if(phi2 < 0) phi2 += 2*PI;
2203                     Float_t deltaphi = phi1 - phi2;
2204                     if(deltaphi > PI) deltaphi -= PI;
2205                     if(deltaphi < -PI) deltaphi += PI;
2206
2207                     if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2208                     else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2209                   }
2210
2211                 }// pair selection
2212
2213                 // Check that label does not exceed stack size
2214                 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2215                   if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2216                   pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2217                   pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2218                   pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2219                   pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2220                   pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2221                   qinv12MC = GetQinv(pVect1MC, pVect2MC);
2222
2223                   if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2224                     ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
2225                     ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2226                     ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2227                   }
2228
2229                   // secondary contamination
2230                   if(ENsum==0){
2231                     mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2232                     mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2233                     if(!mcParticle1 || !mcParticle2) continue;
2234                     if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2235                       if(ch1==ch2) {
2236                         ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2237                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2238                           ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2239                         }
2240                       }else{
2241                         ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2242                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2243                           ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2244                         }
2245                       }
2246                     }
2247                   }
2248
2249                   if(ENsum==6){// all mixed events
2250                     Int_t chGroup2[2]={ch1,ch2};
2251
2252                     Float_t rForQW=5.0;
2253                     if(fFSIindex<=1) rForQW=10;
2254                     else if(fFSIindex==2) rForQW=9;
2255                     else if(fFSIindex==3) rForQW=8;
2256                     else if(fFSIindex==4) rForQW=7;
2257                     else if(fFSIindex==5) rForQW=6;
2258                     else if(fFSIindex==6) rForQW=5;
2259                     else if(fFSIindex==7) rForQW=4;
2260                     else if(fFSIindex==8) rForQW=3;
2261                     else rForQW=2;
2262
2263
2264                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
2265                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
2266                     // pion purity
2267                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2268                     Int_t SCNumber = 1;
2269                     Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2270                     if(PdgCodeSum==22) SCNumber=1;// e-e
2271                     else if(PdgCodeSum==24) SCNumber=2;// e-mu
2272                     else if(PdgCodeSum==222) SCNumber=3;// e-pi
2273                     else if(PdgCodeSum==332) SCNumber=4;// e-k
2274                     else if(PdgCodeSum==2223) SCNumber=5;// e-p
2275                     else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2276                     else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2277                     else if(PdgCodeSum==334) SCNumber=8;// mu-k
2278                     else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2279                     else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2280                     else if(PdgCodeSum==532) SCNumber=11;// pi-k
2281                     else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2282                     else if(PdgCodeSum==642) SCNumber=13;// k-k
2283                     else if(PdgCodeSum==2533) SCNumber=14;// k-p
2284                     else if(PdgCodeSum==4424) SCNumber=15;// p-p
2285                     else {SCNumber=16;}
2286
2287                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2288
2289                     ///////////////////////
2290                     // muon contamination
2291                     Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2292                     Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2293                     pionParent1=kFALSE; pionParent2=kFALSE;
2294                     FilledMCpair12=kTRUE;
2295                     //
2296                     if(abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode)==13 || abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13){// muon check
2297                       Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2298                       if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2299                         pionParent1=kTRUE;
2300                         Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2301                         Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2302                       }
2303                       //
2304                       if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2305                         Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2306                         if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2307                           pionParent2=kTRUE;
2308                           Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2309                           Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2310                         }
2311                       }
2312
2313                       parentQinv12 = GetQinv(Pparent1, Pparent2);
2314
2315                       if(pionParent1 || pionParent2){
2316                         if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2317                           Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2318                           Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2319                           for(Int_t term=1; term<=2; term++){
2320                             for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2321                               Float_t Rvalue = 5+Riter;
2322                               Float_t WInput = 1.0;
2323                               if(term==1) {
2324                                 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
2325                               }else{
2326                                 muonPionK12 = 1.0; pionPionK12=1.0;
2327                               }
2328
2329                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2330                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2331                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2332                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2333                             }// Riter
2334                           }// term loop
2335
2336                           if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2337                           else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2338                         }// parentQ check
2339                       }// pion parent check
2340                     }// muon check
2341
2342
2343                     // momentum resolution
2344                     for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2345                       Float_t Rvalue = 5+Riter;
2346                       Float_t WInput = MCWeight(chGroup2, Rvalue, 0.7, qinv12MC, 0.);
2347                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2348                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2349                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2350                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2351                     }
2352
2353                   }// ENsum check
2354                 }// MC array check
2355               }// MC case
2356
2357
2358
2359               /////////////////////////////////////////////////////////////
2360               for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2361                 if(en3==0) {
2362                   if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2363                   if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2364                 }else if(en3==1){
2365                   if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2366                   if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2367                 }else{
2368                   if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2369                   if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
2370                 }
2371
2372                 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2373                 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2374                 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2375                 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2376                 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2377                 qinv13 = GetQinv(pVect1, pVect3);
2378                 qinv23 = GetQinv(pVect2, pVect3);
2379                 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2380
2381                 FilledMCtriplet123 = kFALSE;
2382
2383                 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2384                 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2385
2386                 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2387                 if(KT3<=fKT3transition) KT3index=0;
2388                 else KT3index=1;
2389
2390                 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2391                 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
2392
2393                 if(ENsum==0) {
2394                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3);
2395                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23));
2396                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2397                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2398                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2399                 }
2400                 if(ENsum==6) {
2401                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2402                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2403                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2404                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2405                 }
2406                 if(ENsum==3){
2407                   if(fill2) {
2408                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
2409                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
2410                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2411                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2412                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2413                   }if(fill3) {
2414                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
2415                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
2416                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2417                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2418                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2419                   }if(fill4) {
2420                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
2421                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
2422                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2423                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2424                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2425                   }
2426                 }
2427
2428                 // r3 denominator
2429                 if(ENsum==6 && ch1==ch2 && ch1==ch3){
2430                   GoodTripletWeights = kFALSE;
2431                   //
2432                   GetWeight(pVect1, pVect2, weight12, weight12Err);
2433                   GetWeight(pVect1, pVect3, weight13, weight13Err);
2434                   GetWeight(pVect2, pVect3, weight23, weight23Err);
2435
2436                   if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2437                     if(fMbin==0 && bin1==0) {
2438                       ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2439                     }
2440                   }else{
2441
2442                     Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
2443                     Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
2444                     if(!fGenerateSignal && !fMCcase) {
2445                       Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2446                       Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2447                       Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
2448                       if(momBin12 >= 20) momBin12 = 19;
2449                       if(momBin13 >= 20) momBin13 = 19;
2450                       if(momBin23 >= 20) momBin23 = 19;
2451                       MomResCorr12 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin12);
2452                       MomResCorr13 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin13);
2453                       MomResCorr23 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin23);
2454                       MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2455                       MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2456                       MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
2457                     }
2458                     // no MRC, no Muon Correction
2459                     weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2460                     weight12CC[0] /= FSICorr12*ffcSq;
2461                     weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2462                     weight13CC[0] /= FSICorr13*ffcSq;
2463                     weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2464                     weight23CC[0] /= FSICorr23*ffcSq;
2465                     if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2466                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2467                     }
2468                     // no Muon Correction
2469                     weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2470                     weight12CC[1] /= FSICorr12*ffcSq;
2471                     weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2472                     weight13CC[1] /= FSICorr13*ffcSq;
2473                     weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2474                     weight23CC[1] /= FSICorr23*ffcSq;
2475                     if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2476                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2477                     }
2478                     // both Corrections
2479                     weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2480                     weight12CC[2] /= FSICorr12*ffcSq;
2481                     weight12CC[2] *= MuonCorr12;
2482                     weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2483                     weight13CC[2] /= FSICorr13*ffcSq;
2484                     weight13CC[2] *= MuonCorr13;
2485                     weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2486                     weight23CC[2] /= FSICorr23*ffcSq;
2487                     weight23CC[2] *= MuonCorr23;
2488                     if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
2489                       if(fMbin==0 && bin1==0) {
2490                         ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
2491                       }
2492                     }else{
2493                       GoodTripletWeights = kTRUE;
2494                       /////////////////////////////////////////////////////
2495                       weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2496                       /////////////////////////////////////////////////////
2497                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
2498                     }// 2nd r3 den check
2499
2500
2501                   }// 1st r3 den check
2502
2503                 }// r3 den
2504
2505
2506                 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2507                   ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
2508                   if(q3<0.06){
2509                     Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2510                     Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2511                     Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2512                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2513                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2514                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
2515                   }
2516                 }
2517                 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2518
2519
2520
2521
2522                 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2523                   if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2524
2525                     pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2526                     pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2527                     pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2528                     pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2529                     qinv13MC = GetQinv(pVect1MC, pVect3MC);
2530                     qinv23MC = GetQinv(pVect2MC, pVect3MC);
2531
2532                     q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2533                     if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2534
2535                     Int_t chGroup3[3]={ch1,ch2,ch3};
2536                     Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
2537                     //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
2538                     //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
2539                     //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
2540                     Float_t kTGroup3[3]={0};
2541
2542                     Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2543                     pionParent3=kFALSE;
2544
2545                     if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2546                       Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2547                       if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2548                         pionParent3=kTRUE;
2549                         Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2550                         Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2551                       }
2552                     }
2553
2554                     parentQinv13 = GetQinv(Pparent1, Pparent3);
2555                     parentQinv23 = GetQinv(Pparent2, Pparent3);
2556                     parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2557
2558                     if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2559                       FilledMCtriplet123=kTRUE;
2560                       if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2561
2562                         Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2563                         //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2564                         //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2565                         //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2566                         Float_t parentkTGroup3[3]={0};
2567
2568                         ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2569                         ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2570                         ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2571
2572                         for(Int_t term=1; term<=4; term++){
2573                           if(term==1) {}
2574                           else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2575                           else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2576                           else {if(!pionParent2 && !pionParent3) continue;}
2577                           for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2578                             Float_t Rvalue = 5+Riter;
2579                             Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2580                             Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2581                             Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
2582                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput);
2583                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput);
2584                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI);
2585                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI);
2586                             //
2587                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC);
2588                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3);
2589                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC);
2590                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3);
2591                           }// Riter
2592                         }// term loop
2593
2594                       }// pion parent check
2595                     }// parentQ check (muon correction)
2596
2597                     // 3-pion momentum resolution
2598                     for(Int_t term=1; term<=5; term++){
2599                       for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2600                         Float_t Rvalue = 5+Riter;
2601                         Float_t WInput = MCWeight3(term, Rvalue, 0.7, chGroup3, QinvMCGroup3, kTGroup3);
2602                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput);
2603                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput);
2604                       }
2605                     }
2606
2607                   }// 3rd particle label check
2608                 }// MCcase and ENsum==6
2609
2610
2611
2612
2613                 /////////////////////////////////////////////////////////////
2614                 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2615                   if(en4==0){
2616                     if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2617                     if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2618                     if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2619                   }else if(en4==1){
2620                     if(en3==0){
2621                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2622                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2623                       if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2624                     }else{
2625                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2626                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2627                       if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2628                     }
2629                   }else if(en4==2){
2630                     if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2631                     if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2632                     if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2633                   }else{
2634                     if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2635                     if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2636                     if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
2637                   }
2638
2639                   pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2640                   pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2641                   pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2642                   pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2643                   ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2644                   qinv14 = GetQinv(pVect1, pVect4);
2645                   qinv24 = GetQinv(pVect2, pVect4);
2646                   qinv34 = GetQinv(pVect3, pVect4);
2647                   q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
2648
2649                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2650                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13);
2651                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23);
2652                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2653                   }
2654
2655                   Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2656                   if(KT4<=fKT4transition) KT4index=0;
2657                   else KT4index=1;
2658
2659                   FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2660                   FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2661                   FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
2662
2663                   Bool_t FillTerms[13]={kFALSE};
2664                   SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2665                   //
2666                   for(int ft=0; ft<13; ft++) {
2667                     Float_t FSIfactor = 1.0;
2668                     if(ft==0) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2669                     else if(ft<=4) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2670                     else if(ft<=10) FSIfactor = 1/(FSICorr12);
2671                     else if(ft==11) FSIfactor = 1/(FSICorr12 * FSICorr34);
2672                     else FSIfactor = 1.0;
2673                     if(FillTerms[ft]) {
2674                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
2675                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
2676                     }
2677                   }
2678
2679                   /////////////////////////////////////////////////////////////
2680                   // r4{2}
2681                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6 && GoodTripletWeights){
2682                     GetWeight(pVect1, pVect4, weight14, weight14Err);
2683                     GetWeight(pVect2, pVect4, weight24, weight24Err);
2684                     GetWeight(pVect3, pVect4, weight34, weight34Err);
2685
2686                     Float_t MomResCorr14=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
2687                     Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
2688                     if(!fGenerateSignal && !fMCcase) {
2689                       Int_t momBin14 = fMomResC2->GetYaxis()->FindBin(qinv14);
2690                       Int_t momBin24 = fMomResC2->GetYaxis()->FindBin(qinv24);
2691                       Int_t momBin34 = fMomResC2->GetYaxis()->FindBin(qinv34);
2692                       if(momBin14 >= 20) momBin14 = 19;
2693                       if(momBin24 >= 20) momBin24 = 19;
2694                       if(momBin34 >= 20) momBin34 = 19;
2695                       MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14);
2696                       MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24);
2697                       MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34);
2698                       MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2699                       MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2700                       MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
2701                     }
2702
2703                     // no MRC, no Muon Correction
2704                     weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2705                     weight14CC[0] /= FSICorr14*ffcSq;
2706                     weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2707                     weight24CC[0] /= FSICorr24*ffcSq;
2708                     weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2709                     weight34CC[0] /= FSICorr34*ffcSq;
2710                     if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2711                       weightTotal  = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2712                       weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2713                       weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2714                       weightTotal /= 3.;
2715                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2716                     }
2717                     // no Muon Correction
2718                     weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2719                     weight14CC[1] /= FSICorr14*ffcSq;
2720                     weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2721                     weight24CC[1] /= FSICorr24*ffcSq;
2722                     weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2723                     weight34CC[1] /= FSICorr34*ffcSq;
2724                     if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2725                       weightTotal  = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]);
2726                       weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]);
2727                       weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]);
2728                       weightTotal /= 3.;
2729                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal);
2730                     }
2731                     // both corrections
2732                     weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2733                     weight14CC[2] /= FSICorr14*ffcSq;
2734                     weight14CC[2] *= MuonCorr14;
2735                     weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2736                     weight24CC[2] /= FSICorr24*ffcSq;
2737                     weight24CC[2] *= MuonCorr24;
2738                     weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2739                     weight34CC[2] /= FSICorr34*ffcSq;
2740                     weight34CC[2] *= MuonCorr34;
2741
2742                     if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity
2743                       if(fMbin==0 && bin1==0 && KT4index==0) {
2744                         ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2])));
2745                       }
2746                     }else{
2747                       /////////////////////////////////////////////////////
2748                       weightTotal  = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]);
2749                       weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]);
2750                       weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]);
2751                       weightTotal /= 3.;
2752                       /////////////////////////////////////////////////////
2753                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal);
2754                     }
2755                   }
2756                   /////////////////////////////////////////////////////////////
2757
2758                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==0){
2759                     ((TH3D*)fOutputList->FindObject("fKT4DistTerm1"))->Fill(fMbin+1, KT4, q4);
2760                     if(q4<0.085){
2761                       Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2762                       Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2763                       Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2764                       Float_t pt4=sqrt(pow(pVect4[1],2)+pow(pVect4[2],2));
2765                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt1);
2766                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt2);
2767                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt3);
2768                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt4);
2769                     }
2770                   }
2771                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT4DistTerm13"))->Fill(fMbin+1, KT4, q4);
2772
2773
2774                   // momenumtum resolution and muon corrections
2775                   if(fMCcase && ENsum==6 && FilledMCtriplet123){// for momentum resolution and muon correction
2776                     if((fEvt+en4)->fTracks[l].fLabel < (fEvt+en4)->fMCarraySize){
2777
2778                       pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2779                       pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
2780                       pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
2781                       pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
2782                       qinv14MC = GetQinv(pVect1MC, pVect4MC);
2783                       qinv24MC = GetQinv(pVect2MC, pVect4MC);
2784                       qinv34MC = GetQinv(pVect3MC, pVect4MC);
2785
2786                       q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) +  pow(qinv24MC,2) +  pow(qinv34MC,2));
2787                       if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC);
2788                       if(ch1==ch2 && ch1==ch3 && ch1==ch4) {
2789                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC);
2790                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC);
2791                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC);
2792                       }
2793                       Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
2794                       Float_t QinvMCGroup4[6]={qinv12MC, qinv13MC, qinv14MC, qinv23MC, qinv24MC, qinv34MC};
2795                       Float_t kTGroup4[6]={0};
2796
2797                       Pparent4[0]=pVect4MC[0]; Pparent4[1]=pVect4MC[1]; Pparent4[2]=pVect4MC[2]; Pparent4[3]=pVect4MC[3];
2798                       pionParent4=kFALSE;
2799                       if(abs((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPdgCode)==13){// muon check
2800                         Int_t MotherLabel4 = (fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fMotherLabel;
2801                         if(abs((fEvt+en4)->fMCtracks[MotherLabel4].fPdgCode)==211) {
2802                           pionParent4=kTRUE;
2803                           Pparent4[1] = (fEvt+en4)->fMCtracks[MotherLabel4].fPx; Pparent4[2] = (fEvt+en4)->fMCtracks[MotherLabel4].fPy; Pparent4[3] = (fEvt+en4)->fMCtracks[MotherLabel4].fPz;
2804                           Pparent4[0] = sqrt(pow(Pparent4[1],2)+pow(Pparent4[2],2)+pow(Pparent4[3],2)+pow(fTrueMassPi,2));
2805                         }
2806                       }
2807
2808                       parentQinv14 = GetQinv(Pparent1, Pparent4);
2809                       parentQinv24 = GetQinv(Pparent2, Pparent4);
2810                       parentQinv34 = GetQinv(Pparent3, Pparent4);
2811                       Float_t parentQ4 = sqrt(pow(parentQ3,2) + pow(parentQinv14,2) + pow(parentQinv24,2) + pow(parentQinv34,2));
2812
2813                       if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){
2814                         if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon
2815
2816                           if(pionParent1) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(1);
2817                           if(pionParent2) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(2);
2818                           if(pionParent3) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(3);
2819                           if(pionParent4) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(4);
2820                           Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34};
2821                           Float_t parentkTGroup4[6]={0};
2822
2823                           for(Int_t term=1; term<=12; term++){
2824                             if(term==1) {}
2825                             else if(term==2) {if(!pionParent1 && !pionParent2 && !pionParent3) continue;}
2826                             else if(term==3) {if(!pionParent1 && !pionParent2 && !pionParent4) continue;}
2827                             else if(term==4) {if(!pionParent1 && !pionParent3 && !pionParent4) continue;}
2828                             else if(term==5) {if(!pionParent2 && !pionParent3 && !pionParent4) continue;}
2829                             else if(term==6) {if(!pionParent1 && !pionParent2) continue;}
2830                             else if(term==7) {if(!pionParent1 && !pionParent3) continue;}
2831                             else if(term==8) {if(!pionParent1 && !pionParent4) continue;}
2832                             else if(term==9) {if(!pionParent2 && !pionParent3) continue;}
2833                             else if(term==10) {if(!pionParent2 && !pionParent4) continue;}
2834                             else if(term==11) {if(!pionParent3 && !pionParent4) continue;}
2835                             else {}
2836                             for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2837                               Float_t Rvalue = 5+Riter;
2838                               Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
2839                               Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
2840                               Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
2841                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(1, Rvalue, q4MC, WInput);
2842                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ4, WInput);
2843                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(1, Rvalue, q4MC, WInputFSI);
2844                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(1, Rvalue, parentQ4, WInputParentFSI);
2845                               //
2846                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(2, Rvalue, q4MC);
2847                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ4);
2848                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(2, Rvalue, q4MC);
2849                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(2, Rvalue, parentQ4);
2850                             }// Riter
2851                           }// term loop
2852
2853                         }// pion parent check
2854                       }// parentQ check (muon correction)
2855
2856                       // 4-pion momentum resolution
2857                       for(Int_t term=1; term<=13; term++){
2858                         for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2859                           Float_t Rvalue = 5+Riter;
2860                           Float_t WInput = MCWeight4(term, Rvalue, 0.7, chGroup4, QinvMCGroup4, kTGroup4);
2861                           Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fIdeal->Fill(Rvalue, q4MC, WInput);
2862                           Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fSmeared->Fill(Rvalue, q4, WInput);
2863                         }
2864                       }
2865
2866                     }// label check particle 4
2867                   }// MCcase
2868
2869                 }// 4th particle
2870               }// 3rd particle
2871             }// 2nd particle
2872           }// 1st particle
2873
2874         }// en4
2875       }// en3
2876     }// en2
2877
2878
2879
2880
2881
2882
2883
2884   // Post output data.
2885   PostData(1, fOutputList);
2886
2887 }
2888 //________________________________________________________________________
2889 void AliFourPion::Terminate(Option_t *)
2890 {
2891   // Called once at the end of the query
2892
2893   cout<<"Done"<<endl;
2894
2895 }
2896 //________________________________________________________________________
2897 Bool_t AliFourPion::AcceptPair(AliFourPionTrackStruct first, AliFourPionTrackStruct second)
2898 {
2899
2900   if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
2901
2902   // propagate through B field to r=1m
2903   Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
2904   if(phi1 > 2*PI) phi1 -= 2*PI;
2905   if(phi1 < 0) phi1 += 2*PI;
2906   Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
2907   if(phi2 > 2*PI) phi2 -= 2*PI;
2908   if(phi2 < 0) phi2 += 2*PI;
2909
2910   Float_t deltaphi = phi1 - phi2;
2911   if(deltaphi > PI) deltaphi -= 2*PI;
2912   if(deltaphi < -PI) deltaphi += 2*PI;
2913   deltaphi = fabs(deltaphi);
2914
2915   if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2916
2917
2918   // propagate through B field to r=1.6m
2919   phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
2920   if(phi1 > 2*PI) phi1 -= 2*PI;
2921   if(phi1 < 0) phi1 += 2*PI;
2922   phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
2923   if(phi2 > 2*PI) phi2 -= 2*PI;
2924   if(phi2 < 0) phi2 += 2*PI;
2925
2926   deltaphi = phi1 - phi2;
2927   if(deltaphi > PI) deltaphi -= 2*PI;
2928   if(deltaphi < -PI) deltaphi += 2*PI;
2929   deltaphi = fabs(deltaphi);
2930
2931   if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2932
2933
2934
2935   //
2936
2937   Int_t ncl1 = first.fClusterMap.GetNbits();
2938   Int_t ncl2 = second.fClusterMap.GetNbits();
2939   Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2940   Double_t shfrac = 0; Double_t qfactor = 0;
2941   for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2942     if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2943       if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2944         sumQ++;
2945         sumCls+=2;
2946         sumSha+=2;}
2947       else {sumQ--; sumCls+=2;}
2948     }
2949     else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2950       sumQ++;
2951       sumCls++;}
2952   }
2953   if (sumCls>0) {
2954     qfactor = sumQ*1.0/sumCls;
2955     shfrac = sumSha*1.0/sumCls;
2956   }
2957
2958   if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2959
2960
2961   return kTRUE;
2962
2963
2964 }
2965 //________________________________________________________________________
2966 Float_t AliFourPion::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2967 {
2968   Float_t arg = G_Coeff/qinv;
2969
2970   if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2971   else {return (exp(-arg)-1)/(-arg);}
2972
2973 }
2974 //________________________________________________________________________
2975 void AliFourPion::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2976 {
2977   Int_t j, k;
2978   Int_t a = i2 - i1;
2979   for (Int_t i = i1; i < i2+1; i++) {
2980     j = (Int_t) (gRandom->Rndm() * a);
2981     k = iarr[j];
2982     iarr[j] = iarr[i];
2983     iarr[i] = k;
2984   }
2985 }
2986
2987
2988 //________________________________________________________________________
2989 Float_t AliFourPion::GetQinv(Float_t track1[], Float_t track2[]){
2990
2991   Float_t qinv = sqrt( fabs(pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2)) );
2992   return qinv;
2993
2994 }
2995 //________________________________________________________________________
2996 void AliFourPion::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
2997
2998   Float_t p0 = track1[0] + track2[0];
2999   Float_t px = track1[1] + track2[1];
3000   Float_t py = track1[2] + track2[2];
3001   Float_t pz = track1[3] + track2[3];
3002
3003   Float_t mt = sqrt(p0*p0 - pz*pz);
3004   Float_t pt = sqrt(px*px + py*py);
3005
3006   Float_t v0 = track1[0] - track2[0];
3007   Float_t vx = track1[1] - track2[1];
3008   Float_t vy = track1[2] - track2[2];
3009   Float_t vz = track1[3] - track2[3];
3010
3011   qout = (px*vx + py*vy)/pt;
3012   qside = (px*vy - py*vx)/pt;
3013   qlong = (p0*vz - pz*v0)/mt;
3014 }
3015 //________________________________________________________________________
3016 void AliFourPion::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]){
3017
3018   if(legoCase){
3019     cout<<"LEGO call to SetWeightArrays"<<endl;
3020
3021     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3022       for(Int_t mb=0; mb<fCentBins; mb++){
3023         fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
3024         fNormWeight[tKbin][mb]->SetDirectory(0);
3025       }
3026     }
3027
3028   }else{
3029
3030     TFile *wFile = new TFile("WeightFile.root","READ");
3031     if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3032     else cout<<"Good Weight File Found!"<<endl;
3033
3034     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3035       for(Int_t mb=0; mb<fCentBins; mb++){
3036
3037         TString *name = new TString("Weight_Kt_");
3038         *name += tKbin;
3039         name->Append("_Ky_0");
3040         name->Append("_M_");
3041         *name += mb;
3042         name->Append("_ED_0");
3043
3044
3045         fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
3046         fNormWeight[tKbin][mb]->SetDirectory(0);
3047
3048
3049       }//mb
3050     }//kt
3051
3052     wFile->Close();
3053   }
3054
3055   cout<<"Done reading weight file"<<endl;
3056
3057 }
3058 //________________________________________________________________________
3059 void AliFourPion::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3060
3061   Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3062   //
3063   Float_t qOut=0,qSide=0,qLong=0;
3064   GetQosl(track1, track2, qOut, qSide, qLong);
3065   qOut = fabs(qOut);
3066   qSide = fabs(qSide);
3067   qLong = fabs(qLong);
3068   Float_t wd=0, xd=0, yd=0, zd=0;
3069   //Float_t qinvtemp=GetQinv(0,track1, track2);
3070   //
3071
3072   if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
3073   else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
3074   else {
3075     for(Int_t i=0; i<fKbinsT-1; i++){
3076       if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
3077     }
3078   }
3079   wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
3080   //
3081   /////////
3082   if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
3083   else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
3084   else {
3085     for(Int_t i=0; i<kQbinsWeights-1; i++){
3086       if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
3087     }
3088     xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
3089   }
3090   //
3091   if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
3092   else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
3093   else {
3094     for(Int_t i=0; i<kQbinsWeights-1; i++){
3095       if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
3096     }
3097     yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
3098   }
3099   //
3100   if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
3101   else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
3102   else {
3103     for(Int_t i=0; i<kQbinsWeights-1; i++){
3104       if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
3105     }
3106     zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
3107   }
3108   //
3109
3110
3111   //Float_t min = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3112   Float_t minErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3113   /*
3114   Float_t deltaW=0;
3115   // kt
3116   deltaW += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - min)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3117   // Qo
3118   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - min)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3119   // Qs
3120   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - min)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3121   // Ql
3122   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - min)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3123   //
3124   wgt = min + deltaW;
3125   */
3126
3127
3128   //
3129   // w interpolation (kt)
3130   Float_t c000 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*wd;
3131   Float_t c100 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*wd;
3132   Float_t c010 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*wd;
3133   Float_t c001 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*wd;
3134   Float_t c110 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*wd;
3135   Float_t c101 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*wd;
3136   Float_t c011 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*wd;
3137   Float_t c111 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*wd;
3138   // x interpolation (qOut)
3139   Float_t c00 = c000*(1-xd) + c100*xd;
3140   Float_t c10 = c010*(1-xd) + c110*xd;
3141   Float_t c01 = c001*(1-xd) + c101*xd;
3142   Float_t c11 = c011*(1-xd) + c111*xd;
3143   // y interpolation (qSide)
3144   Float_t c0 = c00*(1-yd) + c10*yd;
3145   Float_t c1 = c01*(1-yd) + c11*yd;
3146   // z interpolation (qLong)
3147   wgt = (c0*(1-zd) + c1*zd);
3148
3149   ////
3150
3151   // Denominator errors negligible compared to numerator so do not waste cpu time below.
3152   //Float_t deltaWErr=0;
3153   // Kt
3154   /*
3155   deltaWErr += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3156   // Qo
3157   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3158   // Qs
3159   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - minErr)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3160   // Ql
3161   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - minErr)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3162   */
3163   wgtErr = minErr;
3164
3165
3166 }
3167 //________________________________________________________________________
3168 Float_t AliFourPion::MCWeight(Int_t c[2], Float_t R, Float_t fcSq, Float_t qinv, Float_t k12){
3169
3170   Float_t radius = R/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3171   Float_t r12=radius*(1-k12/2.0);
3172   SetFSIindex(R);
3173   Float_t coulCorr12 = FSICorrelation(c[0], c[1], qinv);
3174   if(c[0]==c[1]){
3175     Float_t arg=qinv*r12;
3176     Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3177     EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3178     return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
3179   }else {
3180     return ((1-fcSq) + fcSq*coulCorr12);
3181   }
3182
3183 }
3184 //________________________________________________________________________
3185 Float_t AliFourPion::MCWeightOSL(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t qo, Float_t qs, Float_t ql){
3186
3187   Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3188   Float_t radiusSide = radiusOut;
3189   Float_t radiusLong = radiusOut;
3190   Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3191   Float_t coulCorr12 = FSICorrelation(charge1, charge2, qinv);
3192   if(charge1==charge2){
3193     return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
3194   }else {
3195     return ((1-myDamp) + myDamp*coulCorr12);
3196   }
3197
3198 }
3199
3200 //________________________________________________________________________
3201 Float_t AliFourPion::MCWeight3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3], Float_t kT[3]){
3202   // FSI + QS correlations
3203   if(term==5) return 1.0;
3204
3205   Float_t radius=R/0.19733;
3206   Float_t r12=radius*(1-kT[0]/2.0);
3207   Float_t r13=radius*(1-kT[1]/2.0);
3208   Float_t r23=radius*(1-kT[2]/2.0);
3209
3210   Float_t fc = sqrt(fcSq);
3211
3212   SetFSIindex(R);
3213   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3214   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3215   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3216
3217   if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3218     Float_t arg12=qinv[0]*r12;
3219     Float_t arg13=qinv[1]*r13;
3220     Float_t arg23=qinv[2]*r23;
3221     Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
3222     EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
3223     Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
3224     EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
3225     Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
3226     EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
3227     if(term==1){
3228       Float_t C3QS = 1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2) + exp(-pow(qinv[1]*r13,2))*pow(EW13,2) + exp(-pow(qinv[2]*r23,2))*pow(EW23,2);
3229       C3QS += 2*exp(-(pow(r12,2)*pow(qinv[0],2) + pow(r13,2)*pow(qinv[1],2) + pow(r23,2)*pow(qinv[2],2))/2.)*EW12*EW13*EW23;
3230       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3231       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12;
3232       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13;
3233       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23;
3234       C3 += pow(fc,3)*C3QS*Kfactor12*Kfactor13*Kfactor23;
3235       return C3;
3236     }else if(term==2){
3237       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12);
3238     }else if(term==3){
3239       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13);
3240     }else if(term==4){
3241       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23);
3242     }else return 1.0;
3243
3244   }else{// mixed charge case
3245     Float_t arg=qinv[0]*r12;
3246     Float_t KfactorSC = Kfactor12;
3247     Float_t KfactorMC1 = Kfactor13;
3248     Float_t KfactorMC2 = Kfactor23;
3249     if(c[0]==c[2]) {arg=qinv[1]*r13; KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3250     if(c[1]==c[2]) {arg=qinv[2]*r23; KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3251     Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3252     EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3253     if(term==1){
3254       Float_t C3QS = 1 + exp(-pow(arg,2))*pow(EW,2);
3255       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3256       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(arg,2))*pow(EW,2))*KfactorSC;
3257       C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3258       C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3259       C3 += pow(fc,3)*C3QS*KfactorSC*KfactorMC1*KfactorMC2;
3260       return C3;
3261     }else if(term==2){
3262       if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3263       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3264     }else if(term==3){
3265       return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3266     }else if(term==4){
3267       if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3268       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3269     }else return 1.0;
3270   }
3271
3272 }
3273 //________________________________________________________________________
3274 Float_t AliFourPion::MCWeightFSI3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3]){
3275   // FSI only (no QS correlations)
3276   if(term==5) return 1.0;
3277
3278   Float_t fc = sqrt(fcSq);
3279   SetFSIindex(R);
3280   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3281   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3282   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3283
3284   if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3285     if(term==1){
3286       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3287       C3 += pow(fc,2)*(1-fc)*Kfactor12;
3288       C3 += pow(fc,2)*(1-fc)*Kfactor13;
3289       C3 += pow(fc,2)*(1-fc)*Kfactor23;
3290       C3 += pow(fc,3)*Kfactor12*Kfactor13*Kfactor23;
3291       return C3;
3292     }else if(term==2){
3293       return ((1-fcSq) + fcSq*Kfactor12);
3294     }else if(term==3){
3295       return ((1-fcSq) + fcSq*Kfactor13);
3296     }else if(term==4){
3297       return ((1-fcSq) + fcSq*Kfactor23);
3298     }else return 1.0;
3299
3300   }else{// mixed charge case
3301     Float_t KfactorSC = Kfactor12;
3302     Float_t KfactorMC1 = Kfactor13;
3303     Float_t KfactorMC2 = Kfactor23;
3304     if(c[0]==c[2]) {KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3305     if(c[1]==c[2]) {KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3306     if(term==1){
3307       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3308       C3 += pow(fc,2)*(1-fc)*KfactorSC;
3309       C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3310       C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3311       C3 += pow(fc,3)*KfactorSC*KfactorMC1*KfactorMC2;
3312       return C3;
3313     }else if(term==2){
3314       if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*KfactorSC);
3315       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3316     }else if(term==3){
3317       return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3318     }else if(term==4){
3319       if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*KfactorSC);
3320       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3321     }else return 1.0;
3322   }
3323
3324 }
3325 //________________________________________________________________________
3326 Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6], Float_t kT[6]){
3327   if(term==13) return 1.0;
3328
3329   // Charge ordering:
3330   // ----, ---+, --++, -+++, ++++
3331   //
3332   // term ordering:
3333   // Term 1: 1-2-3-4  (all same-event)
3334   // Term 2: 1-2-3 4 (particle 4 from different event)
3335   // Term 3: 1-2-4 3 (particle 3 from different event)
3336   // Term 4: 1-3-4 2 (particle 2 from different event)
3337   // Term 5: 2-3-4 1 (particle 1 from different event)
3338   // Term 6: 1-2 3 4 (particle 1 and 2 from same event)
3339   // Term 7: 1-3 2 4
3340   // Term 8: 1-4 2 3
3341   // Term 9: 2-3 1 4
3342   // Term 10: 2-4 1 3
3343   // Term 11: 3-4 1 2
3344   // Term 12: 1 2 3 4 (all from different events)
3345
3346   Float_t radius = R/0.19733;
3347   Float_t r[6]={0};
3348   r[0]=radius*(1-kT[0]/2.0);
3349   r[1]=radius*(1-kT[1]/2.0);
3350   r[2]=radius*(1-kT[2]/2.0);
3351   r[3]=radius*(1-kT[3]/2.0);
3352   r[4]=radius*(1-kT[4]/2.0);
3353   r[5]=radius*(1-kT[5]/2.0);
3354
3355   Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3356
3357   Float_t fc = sqrt(fcSq);
3358   SetFSIindex(R);
3359   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3360   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3361   Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3362   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3363   Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3364   Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3365   Float_t arg12=qinv[0]*r[0];
3366   Float_t arg13=qinv[1]*r[1];
3367   Float_t arg14=qinv[2]*r[2];
3368   Float_t arg23=qinv[3]*r[3];
3369   Float_t arg24=qinv[4]*r[4];
3370   Float_t arg34=qinv[5]*r[5];
3371   // Exchange Amplitudes
3372   Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
3373   Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
3374   Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
3375   Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
3376   Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
3377   Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));
3378
3379   if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3380
3381     if(term==1){
3382       Float_t C4QS = 1 + pow(EA12,2) + pow(EA13,2) + pow(EA14,2) + pow(EA23,2) + pow(EA24,2) + pow(EA34,2);// baseline + single pairs
3383       C4QS += pow(EA12,2) * pow(EA34,2);// 2-pairs
3384       C4QS += pow(EA13,2) * pow(EA24,2);// 2-pairs
3385       C4QS += pow(EA14,2) * pow(EA23,2);// 2-pairs
3386       C4QS += 2*EA12*EA13*EA23 + 2*EA12*EA14*EA24 + 2*EA13*EA14*EA34 + 2*EA23*EA24*EA34;// 3-particle exhange
3387       C4QS += 3*EA12*EA23*EA34*EA14 + 3*EA12*EA13*EA34*EA24;// 4-particle exchange
3388       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3389       C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA12,2))*Kfactor12 + (1 + pow(EA13,2))*Kfactor13 + (1 + pow(EA14,2))*Kfactor14 );
3390       C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA23,2))*Kfactor23 + (1 + pow(EA24,2))*Kfactor24 + (1 + pow(EA34,2))*Kfactor34);
3391       C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA13,2) + pow(EA23,2) + 2*EA12*EA13*EA23) * Kfactor12*Kfactor13*Kfactor23;
3392       C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA14,2) + pow(EA24,2) + 2*EA12*EA14*EA24) * Kfactor12*Kfactor14*Kfactor24;
3393       C4 += pow(fc,3)*(1-fc)*(1 + pow(EA13,2) + pow(EA14,2) + pow(EA34,2) + 2*EA13*EA14*EA34) * Kfactor13*Kfactor14*Kfactor34;
3394       C4 += pow(fc,3)*(1-fc)*(1 + pow(EA23,2) + pow(EA24,2) + pow(EA34,2) + 2*EA23*EA24*EA34) * Kfactor23*Kfactor24*Kfactor34;
3395       C4 += pow(fc,4)*C4QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3396       return C4;
3397     }else if(term<=5){
3398       Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0;
3399       if(term==2)      {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3400       else if(term==3) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3401       else if(term==4) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3402       else             {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3403       Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2);
3404       C3QS += 2*EA1*EA2*EA3;
3405       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3406       C3 += pow(fc,2)*(1-fc)*( (1+pow(EA1,2))*Kpair1 + (1+pow(EA2,2))*Kpair2 + (1+pow(EA3,2))*Kpair3 );
3407       C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3408       return C3;
3409     }else if(term<=11){
3410       if(term==6)       return ((1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12);
3411       else if(term==7)  return ((1-fcSq) + fcSq*(1 + pow(EA13,2))*Kfactor13);
3412       else if(term==8)  return ((1-fcSq) + fcSq*(1 + pow(EA14,2))*Kfactor14);
3413       else if(term==9)  return ((1-fcSq) + fcSq*(1 + pow(EA23,2))*Kfactor23);
3414       else if(term==10) return ((1-fcSq) + fcSq*(1 + pow(EA24,2))*Kfactor24);
3415       else              return ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34);
3416     }else if(term==12){
3417       Float_t C22 = (1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12;
3418       C22 *= (1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34;
3419       return C22;
3420     }else return 1.0;
3421
3422   }else{// mixed charge case
3423     if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3424       Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3425       Int_t c_OddOneOut = 1;
3426       if(ChargeSum==3) c_OddOneOut = 0;
3427       //
3428       if(c[0]==c_OddOneOut) {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3429       else if(c[1]==c_OddOneOut) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3430       else if(c[2]==c_OddOneOut) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;   Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3431       else {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;   Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3432
3433       if(term==1){
3434         Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3435         Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3436         C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 + (1 + pow(EA3,2))*Kpair3 );
3437         C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3438         C4 += pow(fc,3)*(1-fc)*(1 + pow(EA1,2) + pow(EA2,2) +  pow(EA3,2) + 2*EA1*EA2*EA3) * Kpair1*Kpair2*Kpair3;
3439         C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair4*Kpair5 + (1+pow(EA2,2))*Kpair2*Kpair4*Kpair6 + (1+pow(EA3,2))*Kpair3*Kpair5*Kpair6);// doesn't matter which MC K's
3440         C4 += pow(fc,4)*C3QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3441         return C4;
3442       }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3443         Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3444         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3445         C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;
3446         C3 += pow(fc,2)*(1-fc)*(1+pow(EA2,2))*Kpair2;
3447         C3 += pow(fc,2)*(1-fc)*(1+pow(EA3,2))*Kpair3;
3448         C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3449         return C3;
3450       }else if(term<=5){// one SC pair, two MC pairs
3451         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3452         C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3453         C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3454         C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3455         C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair5;
3456         return C3;
3457       }else if(term==6 || term==7){
3458         if(ChargeSum==1) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3459         else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3460       }else if(term==8){
3461         return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3462       }else if(term==9){
3463         return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3464       }else if(term==10 || term==11){
3465         if(ChargeSum==3) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3466         else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3467       }else return 1.0;// for 12 and 13
3468     }else{// --++ configuration
3469       Float_t EA1=0, EA2=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3470       if(c[0]==c[1]) {EA1=EA12; EA2=EA34; Kpair1=Kfactor12; Kpair2=Kfactor34;    Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3471       else if(c[0]==c[2]) {EA1=EA13; EA2=EA24; Kpair1=Kfactor13; Kpair2=Kfactor24;    Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3472       else {EA1=EA14; EA2=EA23; Kpair1=Kfactor14; Kpair2=Kfactor23;    Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3473       //
3474       if(term==1){
3475         Float_t C2QS = 1 + pow(EA1,2)*pow(EA2,2);
3476         Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3477         C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 );
3478         C4 += pow(fc,2)*pow(1-fc,2)*( Kpair3 + Kpair4 + Kpair5 + Kpair6 );
3479         C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair3*Kpair4 + (1 + pow(EA2,2))*Kpair2*Kpair3*Kpair4);
3480         C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair5*Kpair6 + (1 + pow(EA2,2))*Kpair2*Kpair5*Kpair6);// doesn't matter which two MC K's used
3481         C4 += pow(fc,4)*C2QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3482         return C4;
3483       }else if(term<=5){
3484         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3485         C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3486         C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3487         C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3488         C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair6;
3489         return C3;
3490       }else if(term==6 || term==11){
3491         return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3492       }else if(term !=12 && term !=13){
3493         return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3494       }else return 1.0;
3495     }
3496   }
3497
3498 }
3499 //________________________________________________________________________
3500 Float_t AliFourPion::MCWeightFSI4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6]){
3501   if(term==13) return 1.0;
3502
3503   Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3504
3505   Float_t fc = sqrt(fcSq);
3506   SetFSIindex(R);
3507   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3508   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3509   Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3510   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3511   Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3512   Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3513
3514   if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3515
3516     if(term==1){
3517       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3518       C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor12 + Kfactor13 + Kfactor14 );
3519       C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor23 + Kfactor24 + Kfactor34 );
3520       C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor13*Kfactor23;
3521       C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor14*Kfactor24;
3522       C4 += pow(fc,3)*(1-fc) * Kfactor13*Kfactor14*Kfactor34;
3523       C4 += pow(fc,3)*(1-fc) * Kfactor23*Kfactor24*Kfactor34;
3524       C4 += pow(fc,4) * Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3525       return C4;
3526     }else if(term<=5){
3527       Float_t Kpair1=0, Kpair2=0, Kpair3=0;
3528       if(term==2) {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3529       else if(term==3) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3530       else if(term==4) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3531       else {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3532       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3533       C3 += pow(fc,2)*(1-fc)*( Kpair1 + Kpair2 + Kpair3 );
3534       C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3535       return C3;
3536     }else if(term<=11){
3537       if(term==6) return ((1-fcSq) + fcSq*Kfactor12);
3538       else if(term==7) return ((1-fcSq) + fcSq*Kfactor13);
3539       else if(term==8) return ((1-fcSq) + fcSq*Kfactor14);
3540       else if(term==9) return ((1-fcSq) + fcSq*Kfactor23);
3541       else if(term==10) return ((1-fcSq) + fcSq*Kfactor24);
3542       else return ((1-fcSq) + fcSq*Kfactor34);
3543     }else if(term==12){
3544       Float_t C22 = (1-fcSq) + fcSq*Kfactor12;
3545       C22 *= (1-fcSq) + fcSq*Kfactor34;
3546       return C22;
3547     }else return 1.0;
3548
3549   }else{// mixed charge case
3550     if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3551       Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3552       Int_t c_OddOneOut = 1;
3553       if(ChargeSum==3) c_OddOneOut = 0;
3554       //
3555       if(c[0]==c_OddOneOut) {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3556       else if(c[1]==c_OddOneOut) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3557       else if(c[2]==c_OddOneOut) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;   Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3558       else {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;   Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3559
3560       if(term==1){
3561         Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3562         C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 );
3563         C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3564         C4 += pow(fc,3)*(1-fc)*Kpair1*Kpair2*Kpair3;
3565         C4 += pow(fc,3)*(1-fc)*(Kpair1*Kpair4*Kpair5 + Kpair2*Kpair4*Kpair6 + Kpair3*Kpair5*Kpair6);// doesn't matter which two MC K's used
3566         C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3567         return C4;
3568       }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3569         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3570         C3 += pow(fc,2)*(1-fc)*Kpair1;
3571         C3 += pow(fc,2)*(1-fc)*Kpair2;
3572         C3 += pow(fc,2)*(1-fc)*Kpair3;
3573         C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3574         return C3;
3575       }else if(term<=5){// one SC pair, two MC pairs
3576         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3577         C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3578         C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3579         C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3580         C3 += pow(fc,3)*Kpair1*Kpair4*Kpair5;
3581         return C3;
3582       }else if(term==6 || term==7){
3583         if(ChargeSum==1) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3584         else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3585       }else if(term==8){
3586         return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3587       }else if(term==9){
3588         return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3589       }else if(term==10 || term==11){
3590         if(ChargeSum==3) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3591         else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3592       }else return 1.0;// 12 and 13
3593     }else{// --++ configuration
3594       Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3595       if(c[0]==c[1]) {Kpair1=Kfactor12; Kpair2=Kfactor34;    Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3596       else if(c[0]==c[2]) {Kpair1=Kfactor13; Kpair2=Kfactor24;    Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3597       else {Kpair1=Kfactor14; Kpair2=Kfactor23;    Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3598       //
3599       if(term==1){
3600         Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3601         C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 + Kpair4 + Kpair5 + Kpair6);
3602         C4 += pow(fc,3)*(1-fc)*( Kpair1*Kpair3*Kpair4 + Kpair2*Kpair3*Kpair4 + Kpair1*Kpair5*Kpair6 + Kpair2*Kpair5*Kpair6);
3603         C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3604         return C4;
3605       }else if(term<=5){
3606         Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3607         C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3608         C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3609         C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3610         C3 += pow(fc,3)*Kpair1*Kpair4*Kpair6;
3611         return C3;
3612       }else if(term==6 || term==11){
3613         return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3614       }else if(term !=12 && term !=13){
3615         return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3616       }else return 1.0;
3617     }
3618   }
3619
3620 }
3621 //________________________________________________________________________
3622 void AliFourPion::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
3623
3624
3625   if(legoCase){
3626     cout<<"LEGO call to SetMomResCorrections"<<endl;
3627     fMomResC2 = (TH2D*)temp2D->Clone();
3628     fMomResC2->SetDirectory(0);
3629   }else {
3630     TFile *momResFile = new TFile("MomResFile.root","READ");
3631     if(!momResFile->IsOpen()) {
3632       cout<<"No momentum resolution file found"<<endl;
3633       AliFatal("No momentum resolution file found.  Kill process.");
3634     }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3635
3636     TH2D *temp2D2 = (TH2D*)momResFile->Get("MRC_C2_SC");
3637     fMomResC2 = (TH2D*)temp2D2->Clone();
3638     fMomResC2->SetDirectory(0);
3639
3640     momResFile->Close();
3641   }
3642
3643
3644   for(Int_t bx=1; bx<=fMomResC2->GetNbinsX(); bx++){
3645     for(Int_t by=1; by<=fMomResC2->GetNbinsY(); by++){
3646       if(fMomResC2->GetBinContent(bx,by) > 1.5) fMomResC2->SetBinContent(bx,by, 1.5);// Maximum is ~1.02
3647       if(fMomResC2->GetBinContent(bx,by) < 0.95) fMomResC2->SetBinContent(bx,by, 0.95);// Minimum is ~0.98
3648     }
3649   }
3650
3651   cout<<"Done reading momentum resolution file"<<endl;
3652 }
3653 //________________________________________________________________________
3654 void AliFourPion::SetFSICorrelations(Bool_t legoCase, TH1D *tempss[12], TH1D *tempos[12]){
3655   // read in 2-particle and 3-particle FSI correlations = K2 & K3
3656   // 2-particle input histo from file is binned in qinv.  3-particle in qinv of each pair
3657   if(legoCase){
3658     cout<<"LEGO call to SetFSICorrelations"<<endl;
3659     for(Int_t MB=0; MB<12; MB++) {
3660       fFSIss[MB] = (TH1D*)tempss[MB]->Clone();
3661       fFSIos[MB] = (TH1D*)tempos[MB]->Clone();
3662       //
3663       fFSIss[MB]->SetDirectory(0);
3664       fFSIos[MB]->SetDirectory(0);
3665     }
3666   }else {
3667     cout<<"non LEGO call to SetFSICorrelations"<<endl;
3668     TFile *fsifile = new TFile("KFile.root","READ");
3669     if(!fsifile->IsOpen()) {
3670       cout<<"No FSI file found"<<endl;
3671       AliFatal("No FSI file found.  Kill process.");
3672     }else {cout<<"Good FSI File Found!"<<endl;}
3673
3674     TH1D *temphistoSS[12];
3675     TH1D *temphistoOS[12];
3676     for(Int_t MB=0; MB<12; MB++) {
3677       TString *nameK2SS = new TString("K2ss_");
3678       *nameK2SS += MB;
3679       temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
3680       //
3681       TString *nameK2OS = new TString("K2os_");
3682       *nameK2OS += MB;
3683       temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
3684       //
3685       fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
3686       fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
3687       fFSIss[MB]->SetDirectory(0);
3688       fFSIos[MB]->SetDirectory(0);
3689     }
3690     //
3691
3692     fsifile->Close();
3693   }
3694
3695   cout<<"Done reading FSI file"<<endl;
3696 }
3697 //________________________________________________________________________
3698 Float_t AliFourPion::FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
3699   // returns 2-particle Coulomb correlations = K2
3700   Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
3701   Int_t qbinH = qbinL+1;
3702   if(qbinL <= 0) return 1.0;
3703   if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) return 1.0;
3704
3705   Float_t slope=0;
3706   if(charge1==charge2){
3707     slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
3708     slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
3709     return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
3710   }else {
3711     slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
3712     slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
3713     return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
3714   }
3715 }
3716 //________________________________________________________________________
3717 void AliFourPion::SetFillBins2(Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3718   if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3719   else {b1=c1; b2=c2;}
3720 }
3721 //________________________________________________________________________
3722 void AliFourPion::SetFillBins3(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){
3723
3724   // "part" specifies which pair is from the same event.  Only relevant for terms 2-4
3725   Bool_t seSS=kFALSE;
3726   if(part==1) {// default case (irrelevant for term 1 and term 5)
3727     if(c1==c2) seSS=kTRUE;
3728   }
3729   if(part==2){
3730     if(c1==c3) seSS=kTRUE;
3731   }
3732
3733
3734   // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3735   if( (c1+c2+c3)==1) {
3736     b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3737     //
3738     if(seSS) fill2=kTRUE;
3739     else {fill3=kTRUE; fill4=kTRUE;}
3740     //
3741   }else if( (c1+c2+c3)==2) {
3742     b1=0; b2=1; b3=1;
3743     //
3744     if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3745     else fill4=kTRUE;
3746     //
3747   }else {
3748     b1=c1; b2=c2; b3=c3;
3749     fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3750   }
3751
3752 }
3753 //________________________________________________________________________
3754 void AliFourPion::SetFillBins4(Int_t c1, Int_t c2, Int_t c3, Int_t c4, Int_t &b1, Int_t &b2, Int_t &b3, Int_t &b4, Int_t ENsum, Bool_t fillTerm[13]){
3755
3756   // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3757   if( (c1+c2+c3+c4)==0 || (c1+c2+c3+c4)==4) {// all of the same charge: ---- or ++++
3758
3759     b1=c1; b2=c2; b3=c3; b4=c4;
3760     if(ENsum==0) fillTerm[0]=kTRUE;
3761     else if(ENsum==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
3762     else if(ENsum==2) {fillTerm[11]=kTRUE;}
3763     else if(ENsum==3) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
3764     else fillTerm[12]=kTRUE;
3765
3766   }else if( (c1+c2+c3+c4)==1) {// one positive charge: ---+
3767
3768     b1=0; b2=0; b3=0; b4=1;// Re-assign to merge degenerate histos
3769     if(ENsum==0) fillTerm[0]=kTRUE;
3770     else if(ENsum==1){
3771       if(c4==1) fillTerm[1]=kTRUE;
3772       else {fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
3773     }else if(ENsum==2){
3774       fillTerm[11]=kTRUE;
3775     }else if(ENsum==3){
3776       if(c3==1 || c4==1) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[8]=kTRUE;}
3777       else {fillTerm[7]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
3778     }else fillTerm[12]=kTRUE;
3779
3780   }else if( (c1+c2+c3+c4)==2) {// two positive charges: --++
3781
3782     b1=0; b2=0; b3=1; b4=1;// Re-assign to merge degenerate histos
3783     if(ENsum==0) fillTerm[0]=kTRUE;
3784     else if(ENsum==1){
3785       if(c4==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE;}
3786       else {fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
3787     }else if(ENsum==2){
3788       if( (c1+c2)==0) fillTerm[11]=kTRUE;
3789     }else if(ENsum==3){
3790       if( (c1+c2)==0) fillTerm[5]=kTRUE;
3791       else if( (c1+c2)==1) {fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE;}
3792       else fillTerm[10]=kTRUE;
3793     }else fillTerm[12]=kTRUE;
3794
3795   }else{// three positive charges
3796
3797     b1=0; b2=1; b3=1; b4=1;// Re-assign to merge degenerate histos
3798     if(ENsum==0) fillTerm[0]=kTRUE;
3799     else if(ENsum==1){
3800       if(c4==0) fillTerm[4]=kTRUE;
3801       else {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE;}
3802     }else if(ENsum==2){
3803       fillTerm[11]=kTRUE;
3804     }else if(ENsum==3){
3805       if(c3==0 || c4==0) {fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
3806       else {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE;}
3807     }else fillTerm[12]=kTRUE;
3808
3809   }
3810
3811 }
3812 //________________________________________________________________________
3813 void AliFourPion::SetFSIindex(Float_t R){
3814   if(!fMCcase){
3815     if(fPbPbcase){
3816       if(fMbin==0) fFSIindex = 0;//0-5%
3817       else if(fMbin==1) fFSIindex = 1;//5-10%
3818       else if(fMbin<=3) fFSIindex = 2;//10-20%
3819       else if(fMbin<=5) fFSIindex = 3;//20-30%
3820       else if(fMbin<=7) fFSIindex = 4;//30-40%
3821       else if(fMbin<=9) fFSIindex = 5;//40-50%
3822       else if(fMbin<=12) fFSIindex = 6;//40-50%
3823       else if(fMbin<=15) fFSIindex = 7;//40-50%
3824       else if(fMbin<=18) fFSIindex = 8;//40-50%
3825       else fFSIindex = 8;//90-100%
3826     }else fFSIindex = 9;// pp and pPb
3827   }else{// FSI binning for MC
3828     if(R>=10.) fFSIindex = 0;
3829     else if(R>=9.) fFSIindex = 1;
3830     else if(R>=8.) fFSIindex = 2;
3831     else if(R>=7.) fFSIindex = 3;
3832     else if(R>=6.) fFSIindex = 4;
3833     else if(R>=5.) fFSIindex = 5;
3834     else if(R>=4.) fFSIindex = 6;
3835     else if(R>=3.) fFSIindex = 7;
3836     else if(R>=2.) fFSIindex = 8;
3837     else fFSIindex = 9;
3838   }
3839 }
3840 //________________________________________________________________________
3841 void AliFourPion::SetMuonCorrections(Bool_t legoCase, TH2D *tempMuon){
3842   if(legoCase){
3843     cout<<"LEGO call to SetMuonCorrections"<<endl;
3844     fWeightmuonCorrection = (TH2D*)tempMuon->Clone();
3845     fWeightmuonCorrection->SetDirectory(0);
3846   }else {
3847     cout<<"non LEGO call to SetMuonCorrections"<<endl;
3848     TFile *MuonFile=new TFile("MuonCorrection.root","READ");
3849     if(!MuonFile->IsOpen()) {
3850       cout<<"No Muon file found"<<endl;
3851       AliFatal("No Muon file found.  Kill process.");
3852     }else {cout<<"Good Muon File Found!"<<endl;}
3853
3854     fWeightmuonCorrection = (TH2D*)MuonFile->Get("WeightmuonCorrection");
3855     fWeightmuonCorrection->SetDirectory(0);
3856     //
3857     MuonFile->Close();
3858   }
3859   cout<<"Done reading Muon file"<<endl;
3860 }