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