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