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