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