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