PWG4/JCORRAN/AliAnalysisTaskDiHadron.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //2- and 3-particle trigger particle correlation analysis
  16 //Author: Jason Glyndwr Ulery, ulery@uni-frankfurt.de
  17 //version: 3.4,  last revised: 2010/08/15
  18
  19 #include "Riostream.h"
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH1F.h"
  23 #include "TH2F.h"
  24 #include "TH3F.h"
  25 #include "TFormula.h"
  26 #include "TF1.h"
  27 #include "TF2.h"
  28 #include "TF3.h"
  29 #include "TVector3.h"
  30 #include "TMath.h"
  31
  32
  33 #include "AliAnalysisTask.h"
  34 #include "AliAnalysisManager.h"
  35
  36 #include "AliESDEvent.h"
  37 #include "AliESDInputHandler.h"
  38 #include "AliMCEvent.h"
  39 #include "AliMCEventHandler.h"
  40 #include "AliMCParticle.h"
  41 #include "AliStack.h"
  42 #include "AliESDVertex.h"
  43 #include "AliMultiplicity.h"
  44 #include "TParticle.h"
  45
  46
  47 //#include "AliHeader.h"
  48 //#include "AliGenEventHeader.h"
  49
  50
  51
  52 #include "AliAnalysisTaskDiHadron.h"
  53
  54
  55 ClassImp(AliAnalysisTaskDiHadron)
  56
  57 //----------------------------------------
  58 AliAnalysisTaskDiHadron::AliAnalysisTaskDiHadron(const char *name):
  59 AliAnalysisTask(name,""), fESD(0), fMC(0), fOutput(0),fMinClustersTPC(0),fMinClusterRatio(0),fMaxTPCchi2(0),fMinClustersITS(0),fEtaCut(0),fTrigEtaCut(0),fNearPhiCut(0),fXECut(0),fMaxDCA(0),fMaxDCAXY(0),fMaxDCAZ(0),fDCA2D(0),fTPCRefit(0),fITSRefit(0),fSPDCut(0),fMinPtAssoc(0),fMaxPtAssoc(0),fVzCut(0),fEfficiencyCorr(0),fDEBUG(0),fnBinPhi(0),fnBinEta(0),fnBinPhiEtaPhi(0),fnBinPhiEtaEta(0),fnBinPhi3(0),fnBinEta3(0),fPi(3.1415926535898),fdPhiMin(0),fdPhiMax(0),fNTPtBins(0),fNMix(0),fNCentBins(0),fNAPtBins(0),fNAPt3Bins(0),fNVertexBins(0),fNXEBins(0),fNIDs(0),fEffFitPt(0),fNFitLowParam(0),fNFitHighParam(0),fMCHistos(0),fFitLow(NULL),fFitHigh(NULL),fFitLowParam(NULL),fFitHighParam(NULL),fPtTrigArray(NULL),fPtAssocArray(NULL),fPtAssoc3Array1(NULL),fPtAssoc3Array2(NULL),fCentArrayMin(NULL),fCentArrayMax(NULL),fXEArray(NULL),fTrigIDArray(NULL),ftPhi(NULL),ftEta(NULL),ftPt(NULL),ftCharge(NULL),ftEff(NULL),ftPtAssoc3(NULL),ftNPtAssoc3(NULL)
  60
  61   {
  62
  63   //IO Slots
  64   DefineInput(0, TChain::Class());
  65   DefineOutput(0,TList::Class());
  66
  67
  68   for(int c=0;c<fNCentBins;c++){
  69     for(int v=0;v<fNVertexBins;v++){
  70       for(int jmc=0;jmc<2;jmc++){
  71         fMixPointer[c][v][jmc]=-1;
  72         fMixEnd[c][v][jmc]=-1;
  73         for(int ievts=0;ievts<fNMix;ievts++){
  74           fMPt[ievts][c][v][jmc]=NULL;
  75           fMPhi[ievts][c][v][jmc]=NULL;
  76           fMEta[ievts][c][v][jmc]=NULL;
  77           for(int dd=0;dd<10;dd++)fMPtAssoc3[ievts][c][v][jmc][dd]=NULL;
  78           fMNPtAssoc3[ievts][c][v][jmc]=NULL;
  79           fMixTrack[ievts][c][v][jmc]=0;
  80         }
  81       }
  82     }
  83   }
  84
  85   }
  86 //--------------------------------------
  87 void AliAnalysisTaskDiHadron::SetCuts(Int_t MinClustersTPC,  Float_t MinClusterRatio, Float_t MaxTPCchi2, Int_t MinClustersITS, Float_t EtaCut, Float_t TrigEtaCut, Float_t NearPhiCut, Float_t XECut, Float_t MaxDCA, Float_t MaxDCAXY, Float_t MaxDCAZ, Int_t DCA2D, Int_t TPCRefit, Int_t ITSRefit, Int_t SPDCut, Float_t MinPtAssoc, Float_t MaxPtAssoc, Float_t VzCut, Int_t NIDs, const char * TrigIDArray){
  88 //Sets the varibles for track and event cuts
  89   fMinClustersTPC=MinClustersTPC;
  90   fMinClusterRatio=MinClusterRatio;
  91   fMaxTPCchi2=MaxTPCchi2;
  92   fMinClustersITS=MinClustersITS;
  93   fEtaCut=EtaCut;
  94   fTrigEtaCut=TrigEtaCut;
  95   fNearPhiCut=NearPhiCut;
  96   fXECut=XECut;
  97   fMaxDCA=MaxDCA;
  98   fMaxDCAXY=MaxDCAXY;
  99   fMaxDCAZ=MaxDCAZ;
 100   fDCA2D=DCA2D;
 101   fTPCRefit=TPCRefit;
 102   fITSRefit=ITSRefit;
 103   fSPDCut=SPDCut;
 104   fMinPtAssoc=MinPtAssoc;
 105   fMaxPtAssoc=MaxPtAssoc;
 106   fVzCut=VzCut;
 107   fNIDs=NIDs;
 108   fTrigIDArray=(char*)TrigIDArray;
 109 }
 110 //--------------------------------------------------------
 111 void AliAnalysisTaskDiHadron::SetOptions(Int_t EfficiencyCorr, Int_t ffDEBUG,  Int_t MCHistos){
 112 //Sets some options
 113   fEfficiencyCorr=EfficiencyCorr;
 114   fDEBUG=ffDEBUG;
 115   fMCHistos=MCHistos;
 116 }
 117 //------------------------------------------------------
 118 void AliAnalysisTaskDiHadron::SetBins(Int_t nBinPhi, Int_t nBinEta, Int_t nBinPhiEtaPhi, Int_t nBinPhiEtaEta, Int_t nBinPhi3, Int_t nBinEta3,Float_t dPhiMin, Float_t dPhiMax, Int_t NTPtBins, Int_t NMixBins, Int_t NCentBins,Int_t NAPtBins, Int_t NAPt3Bins, Int_t NVertexBins, Int_t NXEBins,Float_t *PtTrigArray, Float_t *PtAssocArray,Float_t *PtAssoc3Array1, Float_t *PtAssoc3Array2, Int_t *CentArrayMin, Int_t *CentArrayMax, Float_t *XEArray){
 119 //sets up the histogram binning
 120   fnBinPhi=nBinPhi;
 121   fnBinEta=nBinEta;
 122   fnBinPhiEtaPhi=nBinPhiEtaPhi;
 123   fnBinPhiEtaEta=nBinPhiEtaEta;
 124   fnBinPhi3=nBinPhi3;
 125   fnBinEta3=nBinEta3;
 126   fdPhiMin=dPhiMin;
 127   fdPhiMax=dPhiMax;
 128   fNTPtBins=NTPtBins;
 129   fNMix=NMixBins;
 130   fNCentBins=NCentBins;
 131   fNAPtBins=NAPtBins;
 132   fNAPt3Bins=NAPt3Bins;
 133   fNVertexBins=NVertexBins;
 134   fNXEBins=NXEBins;
 135   fPtTrigArray=new Float_t [fNTPtBins];
 136   for(int i=0;i<fNTPtBins;i++)fPtTrigArray[i]=PtTrigArray[i];
 137   fPtAssocArray=new Float_t [fNAPtBins];
 138   for(int i=0;i<fNAPtBins;i++)fPtAssocArray[i]=PtAssocArray[i];
 139   fPtAssoc3Array1=new Float_t [fNAPt3Bins];
 140   for(int i=0;i<fNAPt3Bins;i++)fPtAssoc3Array1[i]=PtAssoc3Array1[i];
 141   fPtAssoc3Array2=new Float_t [fNAPt3Bins];
 142   for(int i=0;i<fNAPt3Bins;i++)fPtAssoc3Array2[i]=PtAssoc3Array2[i];
 143   fCentArrayMin=new Int_t [fNCentBins];
 144   for(int i=0;i<NCentBins;i++)fCentArrayMin[i]=CentArrayMin[i];
 145   fCentArrayMax=new Int_t [fNCentBins];
 146   for(int i=0;i<NCentBins;i++)fCentArrayMax[i]=CentArrayMax[i];
 147   fXEArray=new Float_t [fNXEBins];
 148   for(int i=0;i<fNXEBins;i++)fXEArray[i]=XEArray[i];
 149  for(int i=0;i<=fNVertexBins;i++)fVertexArray[i]=(2.*i/fNVertexBins-1)*fVzCut;
 150 }
 151 //-------------------------------------------------------
 152 void AliAnalysisTaskDiHadron::SetEfficiencies(Float_t EffFitPt, const TF1 *FitLow, const TF1 *FitHigh, Int_t NFitLowParam, Int_t NFitHighParam, Float_t *FitLowParam, Float_t *FitHighParam){
 153 //Sets up the efficiency corrections
 154   fEffFitPt=EffFitPt;
 155   fFitLow=(TF1*)FitLow;
 156   fFitHigh=(TF1*)FitHigh;
 157   fNFitLowParam=NFitLowParam;
 158   fNFitHighParam=NFitHighParam;
 159   fFitLowParam=new Float_t [fNFitLowParam*fNCentBins];
 160   for(int i=0;i<fNFitLowParam*fNCentBins;i++)fFitLowParam[i]=FitLowParam[i];
 161   fFitHighParam=new Float_t [fNFitHighParam*fNCentBins];
 162   for(int i=0;i<fNFitHighParam*fNCentBins;i++)fFitHighParam[i]=FitHighParam[i];
 163 }
 164
 165 //-----------------------------------------------------------
 166 void AliAnalysisTaskDiHadron::ConnectInputData(Option_t *){
 167   //Connect to ESD
 168   if(fDEBUG)Printf("Connecting");
 169    TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 170    if (!tree&&fDEBUG) {Printf("ERROR: Could not read chain from input slot 0");}
 171   else {
 172     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 173
 174     if (!esdH){
 175       if(fDEBUG)Printf("ERROR: Could not get ESDInputHandler");
 176     }
 177     else{
 178       fESD = esdH->GetEvent();
 179     }
 180
 181     //MC Data handler (so one can calcualte eff)
 182     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 183     if(mcH)fMC=mcH->MCEvent();
 184   }
 185   if(fDEBUG)Printf("Connected");
 186 }
 187
 188 //---------------------------------------------------------
 189 void AliAnalysisTaskDiHadron::CreateOutputObjects(){
 190   //Creates the histograms and list
 191   if(fDEBUG)Printf("Output");
 192   fOutput=new TList();
 193   fOutput->SetName(GetName());
 194   const Int_t buffersize = 256;
 195   char histname[buffersize];
 196   char histtitle[buffersize];
 197   int nptbins=fNAPtBins;
 198   int lptbins=0;
 199   const char *cmc1[2]={"","_MC"};
 200   const char *cmc2[2]={""," MC"};
 201   const char *sign1[3]={"","_LS","_ULS"};
 202   const char *sign2[3]={""," Like-Sign"," Unlike-Sign"};
 203   const char *sign31[4]={"","_LS","_ULT","_ULA"};
 204   const char *sign32[4]={""," Like-Sign"," Trigger-Diff"," Assoc-Diff"};
 205   Float_t etaEdge=fEtaCut+fTrigEtaCut;
 206   Float_t phiArray[fnBinPhi+1];
 207   Float_t etaArray[fnBinEta+1];
 208   Float_t phiEtaArrayPhi[fnBinPhiEtaPhi+1];
 209   Float_t phiEtaArrayEta[fnBinPhiEtaEta+1];
 210   for(int iphi=0;iphi<=fnBinPhi;iphi++){
 211     phiArray[iphi]=fdPhiMin+iphi*2*fPi/fnBinPhi;
 212   }
 213   for(int ieta=0;ieta<=fnBinEta;ieta++){
 214     etaArray[ieta]=-etaEdge+ieta*2*etaEdge/fnBinEta;
 215   }
 216   for(int iphi=0;iphi<=fnBinPhiEtaPhi;iphi++){
 217     phiEtaArrayPhi[iphi]=fdPhiMin+iphi*2*fPi/fnBinPhiEtaPhi;
 218   }
 219   for(int ieta=0;ieta<=fnBinPhiEtaEta;ieta++){
 220     phiEtaArrayEta[ieta]=-etaEdge+ieta*2*etaEdge/fnBinPhiEtaEta;
 221   }
 222   for(int imc=0;imc<=1;imc++){//MC loop
 223     if(imc==1&&!fMCHistos) continue;
 224     //Create the histograms
 225     snprintf(histname,buffersize,"fHistMult%s",cmc1[imc]);
 226     snprintf(histtitle,buffersize,"Multiplicity%s",cmc2[imc]);
 227     fHistMult[imc]=new TH1F(histname,histtitle,2000,-0.5,1999.5);
 228     fHistMult[imc]->Sumw2();
 229     fHistMult[imc]->GetXaxis()->SetTitle("Number of tracks");
 230     fHistMult[imc]->GetYaxis()->SetTitle("Counts");
 231     fOutput->Add(fHistMult[imc]);
 232
 233     for(int imult=0;imult<fNCentBins;imult++){//loop for multiplicity bins
 234
 235       //Histograms that are independent of the trigger
 236       snprintf(histname,buffersize,"fHistPt_C%d%s",imult,cmc1[imc]);
 237       snprintf(histtitle,buffersize,"P_{T} Distribution of Tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 238       fHistPt[imult][imc]=new TH1F(histname,histtitle,nptbins,fPtAssocArray);
 239       fHistPt[imult][imc]->Sumw2();
 240       fHistPt[imult][imc]->GetXaxis()->SetTitle("p_{T}");
 241       fHistPt[imult][imc]->GetYaxis()->SetTitle("Counts");
 242       fOutput->Add(fHistPt[imult][imc]);
 243
 244       //Histograms that are independent of the trigger
 245       snprintf(histname,buffersize,"fHistPtEff_C%d%s",imult,cmc1[imc]);
 246       snprintf(histtitle,buffersize,"P_{T} Distribution of Tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 247       fHistPtEff[imult][imc]=new TH1F(histname,histtitle,1000,0,100);
 248       fHistPtEff[imult][imc]->Sumw2();
 249       fHistPtEff[imult][imc]->GetXaxis()->SetTitle("p_{T}");
 250       fHistPtEff[imult][imc]->GetYaxis()->SetTitle("Counts");
 251       fOutput->Add(fHistPtEff[imult][imc]);
 252
 253       snprintf(histname,buffersize,"fHistPhi_C%d%s",imult,cmc1[imc]);
 254       snprintf(histtitle,buffersize,"#phi Distribution of Tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 255       fHistPhi[imult][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,nptbins,fPtAssocArray);
 256       fHistPhi[imult][imc]->Sumw2();
 257       fHistPhi[imult][imc]->GetXaxis()->SetTitle("#phi");
 258       fHistPhi[imult][imc]->GetYaxis()->SetTitle("P_{T}");
 259       fOutput->Add(fHistPhi[imult][imc]);
 260
 261       snprintf(histname,buffersize,"fHistPhiPt_C%d%s",imult,cmc1[imc]);
 262       snprintf(histtitle,buffersize,"P_{T} weighted #phi Distribution of tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 263       fHistPhiPt[imult][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,nptbins,fPtAssocArray);
 264       fHistPhiPt[imult][imc]->Sumw2();
 265       fHistPhiPt[imult][imc]->GetXaxis()->SetTitle("#phi");
 266       fHistPhiPt[imult][imc]->GetYaxis()->SetTitle("P_{T}");
 267       fOutput->Add(fHistPhiPt[imult][imc]);
 268
 269       snprintf(histname,buffersize,"fHistEta_C%d%s",imult,cmc1[imc]);
 270       snprintf(histtitle,buffersize,"#eta Distribution of Tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 271       fHistEta[imult][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,nptbins,fPtAssocArray);
 272       fHistEta[imult][imc]->Sumw2();
 273       fHistEta[imult][imc]->GetXaxis()->SetTitle("#eta");
 274       fHistEta[imult][imc]->GetYaxis()->SetTitle("P_{T}");
 275       fOutput->Add(fHistEta[imult][imc]);
 276
 277       snprintf(histname,buffersize,"fHistEtaPt_C%d%s",imult,cmc1[imc]);
 278       snprintf(histtitle,buffersize,"P_{T} weighted #eta Distribution of tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 279       fHistEtaPt[imult][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,nptbins,fPtAssocArray);
 280       fHistEtaPt[imult][imc]->Sumw2();
 281       fHistEtaPt[imult][imc]->GetXaxis()->SetTitle("#eta");
 282       fHistEtaPt[imult][imc]->GetYaxis()->SetTitle("P_{T}");
 283       fOutput->Add(fHistEtaPt[imult][imc]);
 284
 285       snprintf(histname,buffersize,"fHistNEvents_C%d%s",imult,cmc1[imc]);
 286       snprintf(histtitle,buffersize,"Number of Events and Number Passing Cuts %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 287       fHistNEvents[imult][imc]=new TH1F(histname,histtitle,2,-0.5,1.5);
 288       fHistNEvents[imult][imc]->Sumw2();
 289       fHistNEvents[imult][imc]->GetXaxis()->SetTitle("Events,Passing Cuts");
 290       fHistNEvents[imult][imc]->GetYaxis()->SetTitle("Number of Events");
 291       fOutput->Add(fHistNEvents[imult][imc]);
 292
 293       snprintf(histname,buffersize,"fHistNTrigger_C%d%s",imult,cmc1[imc]);
 294       snprintf(histtitle,buffersize,"Number of Triggers %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 295       fHistNTrigger[imult][imc]=new TH1F(histname,histtitle,fNTPtBins,-0.5,fNTPtBins-0.5);
 296       fHistNTrigger[imult][imc]->Sumw2();
 297       fHistNTrigger[imult][imc]->GetXaxis()->SetTitle("Trigger Number");
 298       fHistNTrigger[imult][imc]->GetYaxis()->SetTitle("Number of Triggers");
 299       fOutput->Add(fHistNTrigger[imult][imc]);
 300
 301       snprintf(histname,buffersize,"fHistNTriggerPt_C%d%s",imult,cmc1[imc]);
 302       snprintf(histtitle,buffersize,"P_{T} Weighted Number of Triggers %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 303       fHistNTriggerPt[imult][imc]=new TH1F(histname,histtitle,fNTPtBins,-0.5,fNTPtBins-0.5);
 304       fHistNTriggerPt[imult][imc]->Sumw2();
 305       fHistNTriggerPt[imult][imc]->GetXaxis()->SetTitle("Trigger Number");
 306       fHistNTriggerPt[imult][imc]->GetYaxis()->SetTitle("Number of Triggers");
 307       fOutput->Add(fHistNTriggerPt[imult][imc]);
 308
 309       snprintf(histname,buffersize,"fHistNMix_C%d%s",imult,cmc1[imc]);
 310       snprintf(histtitle,buffersize,"Number of Mixed Events %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 311       fHistNMix[imult][imc]=new TH1F(histname,histtitle,fNTPtBins,-0.5,fNTPtBins-0.5);
 312       fHistNMix[imult][imc]->Sumw2();
 313       fHistNMix[imult][imc]->GetXaxis()->SetTitle("Trigger Number");
 314       fHistNMix[imult][imc]->GetYaxis()->SetTitle("Number of Mixed Events");
 315       fOutput->Add(fHistNMix[imult][imc]);
 316
 317       snprintf(histname,buffersize,"fHistPhiEta_C%d%s",imult,cmc1[imc]);
 318       snprintf(histtitle,buffersize,"#phi-#eta distribution of tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 319       fHistPhiEta[imult][imc]=new TH3F(histname, histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,nptbins,fPtAssocArray);
 320       fHistPhiEta[imult][imc]->Sumw2();
 321       fHistPhiEta[imult][imc]->GetXaxis()->SetTitle("#phi");
 322       fHistPhiEta[imult][imc]->GetYaxis()->SetTitle("#eta");
 323       fHistPhiEta[imult][imc]->GetZaxis()->SetTitle("p_{T}");
 324       fOutput->Add(fHistPhiEta[imult][imc]);
 325
 326       snprintf(histname,buffersize,"fHistPhiEtaPt_C%d%s",imult,cmc1[imc]);
 327       snprintf(histtitle,buffersize,"Pt Weighted #phi-#eta distribution of tracks %dMult%d%s",fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 328       fHistPhiEtaPt[imult][imc]=new TH3F(histname, histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,nptbins,fPtAssocArray);
 329       fHistPhiEtaPt[imult][imc]->Sumw2();
 330       fHistPhiEtaPt[imult][imc]->GetXaxis()->SetTitle("#phi");
 331       fHistPhiEtaPt[imult][imc]->GetYaxis()->SetTitle("#eta");
 332       fHistPhiEtaPt[imult][imc]->GetZaxis()->SetTitle("p_{T}");
 333       fOutput->Add(fHistPhiEtaPt[imult][imc]);
 334
 335       //if(fDEBUG)Printf("OutPut2");
 336       //Histograms with a trigger dependence
 337
 338       for(int i=0;i<fNTPtBins;i++){
 339         for(int j=1;j<fNAPtBins;j++){
 340           if(fabs(fPtTrigArray[i]-fPtAssocArray[j])<1E-5)lptbins=j;
 341         }
 342         //if(fDEBUG)Printf("Loop: %d Pt %3.2f",i,fPtTrigArray[i]/fPtBinWidth);
 343
 344         //Ones with no centrality binning
 345         if(imult==0){
 346           snprintf(histname,buffersize,"fHistMultTrig_P%d%s",i,cmc1[imc]);
 347           snprintf(histtitle,buffersize,"Distrubition of number of tracks in triggered events with %3.1f<p_{T}^{Trig}<%3.1f%s",fPtTrigArray[i],fPtTrigArray[i+1],cmc2[imc]);
 348           fHistMultTrig[i][imc]=new TH1F(histname,histtitle,2000,0,2000);
 349           fHistMultTrig[i][imc]->Sumw2();
 350           fHistMultTrig[i][imc]->GetXaxis()->SetTitle("Number of Tracks");
 351           fHistMultTrig[i][imc]->GetYaxis()->SetTitle("Counts");
 352           fOutput->Add(fHistMultTrig[i][imc]);
 353         }
 354         snprintf(histname,buffersize,"fHistPtTrig_P%d_C%d%s",i,imult,cmc1[imc]);
 355         snprintf(histtitle,buffersize,"P_{T} distribution in triggered events with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 356         fHistPtTrig[i][imult][imc]=new TH1F(histname,histtitle,nptbins,fPtAssocArray);
 357         fHistPtTrig[i][imult][imc]->Sumw2();
 358         fHistPtTrig[i][imult][imc]->GetXaxis()->SetTitle("p_{T}");
 359         fHistPtTrig[i][imult][imc]->GetYaxis()->SetTitle("Counts");
 360         fOutput->Add(fHistPtTrig[i][imult][imc]);
 361
 362         snprintf(histname,buffersize,"fHistPhiTrig_P%d_C%d%s",i,imult,cmc1[imc]);
 363         snprintf(histtitle,buffersize,"Phi Distribution of triggered events with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 364         fHistPhiTrig[i][imult][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 365         fHistPhiTrig[i][imult][imc]->Sumw2();
 366         fHistPhiTrig[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 367         fHistPhiTrig[i][imult][imc]->GetYaxis()->SetTitle("p_{T}");
 368         fOutput->Add(fHistPhiTrig[i][imult][imc]);
 369
 370         snprintf(histname,buffersize,"fHistPhiTrigPt_P%d_C%d%s",i,imult,cmc1[imc]);
 371         snprintf(histtitle,buffersize,"P_{T} Weighted Phi Distribution of triggered events with %3.1f<p_{T}^{Trig}<%3.1f%s",fPtTrigArray[i],fPtTrigArray[i+1],cmc2[imc]);
 372         fHistPhiTrigPt[i][imult][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 373         fHistPhiTrigPt[i][imult][imc]->Sumw2();
 374         fHistPhiTrigPt[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 375         fHistPhiTrigPt[i][imult][imc]->GetYaxis()->SetTitle("p_{T}");
 376         fOutput->Add(fHistPhiTrigPt[i][imult][imc]);
 377
 378         snprintf(histname,buffersize,"fHistEtaTrig_P%d_C%d%s",i,imult,cmc1[imc]);
 379         snprintf(histtitle,buffersize,"Eta Distribution of triggered events with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 380         fHistEtaTrig[i][imult][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 381         fHistEtaTrig[i][imult][imc]->Sumw2();
 382         fHistEtaTrig[i][imult][imc]->GetXaxis()->SetTitle("#eta");
 383         fHistEtaTrig[i][imult][imc]->GetYaxis()->SetTitle("p_{T}");
 384         fOutput->Add(fHistEtaTrig[i][imult][imc]);
 385
 386         snprintf(histname,buffersize,"fHistEtaTrigPt_P%d_C%d%s",i,imult,cmc1[imc]);
 387         snprintf(histtitle,buffersize,"P_{T} Weighted Eta Distribution of triggered events with %3.1f<p_{T}^{Trig}<%3.1f%s",fPtTrigArray[i],fPtTrigArray[i+1],cmc2[imc]);
 388         fHistEtaTrigPt[i][imult][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 389         fHistEtaTrigPt[i][imult][imc]->Sumw2();
 390         fHistEtaTrigPt[i][imult][imc]->GetXaxis()->SetTitle("#eta");
 391         fHistEtaTrigPt[i][imult][imc]->GetYaxis()->SetTitle("p_{T}");
 392         fOutput->Add(fHistEtaTrigPt[i][imult][imc]);
 393
 394         snprintf(histname,buffersize,"fHistPhiEtaTrig_P%d_C%d%s",i,imult,cmc1[imc]);
 395         snprintf(histtitle,buffersize,"#phi-#eta distribution in triggered events %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 396         fHistPhiEtaTrig[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 397         fHistPhiEtaTrig[i][imult][imc]->Sumw2();
 398         fHistPhiEtaTrig[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 399         fHistPhiEtaTrig[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 400         fHistPhiEtaTrig[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 401         fOutput->Add(fHistPhiEtaTrig[i][imult][imc]);
 402
 403         snprintf(histname,buffersize,"fHistXEN_P%d_C%d%s",i,imult,cmc1[imc]);
 404         snprintf(histtitle,buffersize,"Near-Side X_{E} distribution for %3.1f<p_{T}^{Lead}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 405         fHistXEN[i][imult][imc]=new TH1F(histname,histtitle,fNXEBins,fXEArray);
 406         fHistXEN[i][imult][imc]->Sumw2();
 407         fHistXEN[i][imult][imc]->GetXaxis()->SetTitle("X_{E}");
 408         fOutput->Add(fHistXEN[i][imult][imc]);
 409
 410         snprintf(histname,buffersize,"fHistXENMixed_P%d_C%d%s",i,imult,cmc1[imc]);
 411         snprintf(histtitle,buffersize,"Mixed Near-Side X_{E} distribution for %3.1f<p_{T}^{Lead}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 412         fHistXENMix[i][imult][imc]=new TH1F(histname,histtitle,fNXEBins,fXEArray);
 413         fHistXENMix[i][imult][imc]->Sumw2();
 414         fHistXENMix[i][imult][imc]->GetXaxis()->SetTitle("X_{E}");
 415         fOutput->Add(fHistXENMix[i][imult][imc]);
 416
 417         snprintf(histname,buffersize,"fHistXEA_P%d_C%d%s",i,imult,cmc1[imc]);
 418         snprintf(histtitle,buffersize,"Away-Side X_{E} distribution for %3.1f<p_{T}^{Lead}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 419         fHistXEA[i][imult][imc]=new TH1F(histname,histtitle,fNXEBins,fXEArray);
 420         fHistXEA[i][imult][imc]->Sumw2();
 421         fHistXEA[i][imult][imc]->GetXaxis()->SetTitle("X_{E}");
 422         fOutput->Add(fHistXEA[i][imult][imc]);
 423
 424         snprintf(histname,buffersize,"fHistXEAMixed_P%d_C%d%s",i,imult,cmc1[imc]);
 425         snprintf(histtitle,buffersize,"Mixed Away-Side X_{E} distribution for %3.1f<p_{T}^{Lead}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 426         fHistXEAMix[i][imult][imc]=new TH1F(histname,histtitle,fNXEBins,fXEArray);
 427         fHistXEAMix[i][imult][imc]->Sumw2();
 428         fHistXEAMix[i][imult][imc]->GetXaxis()->SetTitle("X_{E}");
 429         fOutput->Add(fHistXEAMix[i][imult][imc]);
 430
 431         //signloop
 432         for(int isign=0;isign<3;isign++){
 433           snprintf(histname,buffersize,"fHistDeltaPhi_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 434           snprintf(histtitle,buffersize,"#Delta#phi Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 435           fHistDeltaPhi[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 436           fHistDeltaPhi[i][imult][isign][imc]->Sumw2();
 437           fHistDeltaPhi[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#phi");
 438           fHistDeltaPhi[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 439           fOutput->Add(fHistDeltaPhi[i][imult][isign][imc]);
 440
 441           snprintf(histname,buffersize,"fHistDeltaPhiPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 442           snprintf(histtitle,buffersize,"P_{T} Weighted #Delta#phi Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 443           fHistDeltaPhiPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 444           fHistDeltaPhiPt[i][imult][isign][imc]->Sumw2();
 445           fHistDeltaPhiPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#phi");
 446           fHistDeltaPhiPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 447           fOutput->Add(fHistDeltaPhiPt[i][imult][isign][imc]);
 448
 449           snprintf(histname,buffersize,"fHistDeltaPhiMix_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 450           snprintf(histtitle,buffersize,"#Delta#phi Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 451           fHistDeltaPhiMix[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 452           fHistDeltaPhiMix[i][imult][isign][imc]->Sumw2();
 453           fHistDeltaPhiMix[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#phi");
 454           fHistDeltaPhiMix[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 455           fOutput->Add(fHistDeltaPhiMix[i][imult][isign][imc]);
 456
 457           snprintf(histname,buffersize,"fHistDeltaPhiMixPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 458           snprintf(histtitle,buffersize,"P_{T} Weighted #Delta#phi Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 459           fHistDeltaPhiMixPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinPhi,phiArray,lptbins,fPtAssocArray);
 460           fHistDeltaPhiMixPt[i][imult][isign][imc]->Sumw2();
 461           fHistDeltaPhiMixPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#phi");
 462           fHistDeltaPhiMixPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 463           fOutput->Add(fHistDeltaPhiMixPt[i][imult][isign][imc]);
 464
 465           //etaNear
 466           snprintf(histname,buffersize,"fHistDeltaEtaN_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 467           snprintf(histtitle,buffersize,"Near-Side #Delta#eta Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 468           fHistDeltaEtaN[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 469           fHistDeltaEtaN[i][imult][isign][imc]->Sumw2();
 470           fHistDeltaEtaN[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 471           fHistDeltaEtaN[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 472           fOutput->Add(fHistDeltaEtaN[i][imult][isign][imc]);
 473
 474           snprintf(histname,buffersize,"fHistDeltaEtaNPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 475           snprintf(histtitle,buffersize,"Near-Side P_{T} Weighted #Delta#eta Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 476           fHistDeltaEtaNPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 477           fHistDeltaEtaNPt[i][imult][isign][imc]->Sumw2();
 478           fHistDeltaEtaNPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 479           fHistDeltaEtaNPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 480           fOutput->Add(fHistDeltaEtaNPt[i][imult][isign][imc]);
 481
 482           snprintf(histname,buffersize,"fHistDeltaEtaNMix_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 483           snprintf(histtitle,buffersize,"Near-Side #Delta#eta Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 484           fHistDeltaEtaNMix[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 485           fHistDeltaEtaNMix[i][imult][isign][imc]->Sumw2();
 486           fHistDeltaEtaNMix[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 487           fHistDeltaEtaNMix[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 488           fOutput->Add(fHistDeltaEtaNMix[i][imult][isign][imc]);
 489
 490           snprintf(histname,buffersize,"fHistDeltaEtaNMixPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 491           snprintf(histtitle,buffersize,"Near-Side P_{T} Weighted #Delta#eta Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 492           fHistDeltaEtaNMixPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 493           fHistDeltaEtaNMixPt[i][imult][isign][imc]->Sumw2();
 494           fHistDeltaEtaNMixPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 495           fHistDeltaEtaNMixPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 496           fOutput->Add(fHistDeltaEtaNMixPt[i][imult][isign][imc]);
 497
 498           //Away Eta
 499           snprintf(histname,buffersize,"fHistDeltaEtaA_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 500           snprintf(histtitle,buffersize,"Away-Side #Delta#eta Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 501           fHistDeltaEtaA[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 502           fHistDeltaEtaA[i][imult][isign][imc]->Sumw2();
 503           fHistDeltaEtaA[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 504           fHistDeltaEtaA[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 505           fOutput->Add(fHistDeltaEtaA[i][imult][isign][imc]);
 506
 507           snprintf(histname,buffersize,"fHistDeltaEtaAPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 508           snprintf(histtitle,buffersize,"Away-Side P_{T} Weighted #Delta#eta Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 509           fHistDeltaEtaAPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 510           fHistDeltaEtaAPt[i][imult][isign][imc]->Sumw2();
 511           fHistDeltaEtaAPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 512           fHistDeltaEtaAPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 513           fOutput->Add(fHistDeltaEtaAPt[i][imult][isign][imc]);
 514
 515           snprintf(histname,buffersize,"fHistDeltaEtaAMix_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 516           snprintf(histtitle,buffersize,"Away-Side #Delta#eta Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 517           fHistDeltaEtaAMix[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 518           fHistDeltaEtaAMix[i][imult][isign][imc]->Sumw2();
 519           fHistDeltaEtaAMix[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 520           fHistDeltaEtaAMix[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 521           fOutput->Add(fHistDeltaEtaAMix[i][imult][isign][imc]);
 522
 523           snprintf(histname,buffersize,"fHistDeltaEtaAMixPt_P%d_C%d%s%s",i,imult,sign1[isign],cmc1[imc]);
 524           snprintf(histtitle,buffersize,"Away-Side P_{T} Weighted #Delta#eta Mixed Event Distribution with %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],sign2[isign],cmc2[imc]);
 525           fHistDeltaEtaAMixPt[i][imult][isign][imc]=new TH2F(histname,histtitle,fnBinEta,etaArray,lptbins,fPtAssocArray);
 526           fHistDeltaEtaAMixPt[i][imult][isign][imc]->Sumw2();
 527           fHistDeltaEtaAMixPt[i][imult][isign][imc]->GetXaxis()->SetTitle("#Delta#eta");
 528           fHistDeltaEtaAMixPt[i][imult][isign][imc]->GetYaxis()->SetTitle("p_{T}");
 529           fOutput->Add(fHistDeltaEtaAMixPt[i][imult][isign][imc]);
 530
 531
 532       //====
 533         }//end isignloop
 534       snprintf(histname,buffersize,"fHistDeltaPhiEta_P%d_C%d%s",i,imult,cmc1[imc]);
 535       snprintf(histtitle,buffersize,"#Delta#phi-#Delta#eta %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 536       fHistDeltaPhiEta[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 537       fHistDeltaPhiEta[i][imult][imc]->Sumw2();
 538       fHistDeltaPhiEta[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 539       fHistDeltaPhiEta[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 540       fHistDeltaPhiEta[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 541       fOutput->Add(fHistDeltaPhiEta[i][imult][imc]);
 542
 543       snprintf(histname,buffersize,"fHistDeltaPhiEtaMix_P%d_C%d%s",i,imult,cmc1[imc]);
 544       snprintf(histtitle,buffersize,"#Delta#phi-#Delta#eta from Mixed Events %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 545       fHistDeltaPhiEtaMix[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 546       fHistDeltaPhiEtaMix[i][imult][imc]->Sumw2();
 547       fHistDeltaPhiEtaMix[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 548       fHistDeltaPhiEtaMix[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 549       fHistDeltaPhiEtaMix[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 550       fOutput->Add(fHistDeltaPhiEtaMix[i][imult][imc]);
 551
 552       snprintf(histname,buffersize,"fHistPhiEtaTrigPt_P%d_C%d%s",i,imult,cmc1[imc]);
 553       snprintf(histtitle,buffersize,"P_{T}-Weighted #phi-#eta distribution in triggered events %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 554       fHistPhiEtaTrigPt[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 555       fHistPhiEtaTrigPt[i][imult][imc]->Sumw2();
 556       fHistPhiEtaTrigPt[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 557       fHistPhiEtaTrigPt[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 558       fHistPhiEtaTrigPt[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 559       fOutput->Add(fHistPhiEtaTrigPt[i][imult][imc]);
 560
 561       snprintf(histname,buffersize,"fHistDeltaPhiEtaPt_P%d_C%d%s",i,imult,cmc1[imc]);
 562       snprintf(histtitle,buffersize,"P_{T}-Weighted #Delta#phi-#Delta#eta %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 563       fHistDeltaPhiEtaPt[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 564       fHistDeltaPhiEtaPt[i][imult][imc]->Sumw2();
 565       fHistDeltaPhiEtaPt[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 566       fHistDeltaPhiEtaPt[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 567       fHistDeltaPhiEtaPt[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 568       fOutput->Add(fHistDeltaPhiEtaPt[i][imult][imc]);
 569
 570       snprintf(histname,buffersize,"fHistDeltaPhiEtaMixPt_P%d_C%d%s",i,imult,cmc1[imc]);
 571       snprintf(histtitle,buffersize,"P_{T}-Weighted #Delta#phi-#Delta#eta from Mixed Events %3.1f<p_{T}^{Trig}<%3.1f %dMult%d%s",fPtTrigArray[i],fPtTrigArray[i+1],fCentArrayMin[imult],fCentArrayMax[imult],cmc2[imc]);
 572       fHistDeltaPhiEtaMixPt[i][imult][imc]=new TH3F(histname,histtitle,fnBinPhiEtaPhi,phiEtaArrayPhi,fnBinPhiEtaEta,phiEtaArrayEta,lptbins,fPtAssocArray);
 573       fHistDeltaPhiEtaMixPt[i][imult][imc]->Sumw2();
 574       fHistDeltaPhiEtaMixPt[i][imult][imc]->GetXaxis()->SetTitle("#phi");
 575       fHistDeltaPhiEtaMixPt[i][imult][imc]->GetYaxis()->SetTitle("#eta");
 576       fHistDeltaPhiEtaMixPt[i][imult][imc]->GetZaxis()->SetTitle("p_{T}");
 577       fOutput->Add(fHistDeltaPhiEtaMixPt[i][imult][imc]);
 578
 579       //Three-Particle Histograms
 580       for(int ipt=0;ipt<fNAPt3Bins;ipt++){
 581         for(int iSign=0;iSign<4;iSign++){
 582           snprintf(histname,buffersize,"fHistDeltaPhiPhi_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 583           snprintf(histtitle,buffersize,"Raw #Delta#phi-#Delta#phi %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 584           fHistDeltaPhiPhi[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinPhi3,fdPhiMin,fdPhiMax,fnBinPhi3,fdPhiMin,fdPhiMax);
 585           fHistDeltaPhiPhi[i][ipt][imult][iSign][imc]->Sumw2();
 586           fHistDeltaPhiPhi[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#phi_{1}");
 587           fHistDeltaPhiPhi[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#phi_{2}");
 588           fOutput->Add(fHistDeltaPhiPhi[i][ipt][imult][iSign][imc]);
 589
 590           snprintf(histname,buffersize,"fHistDeltaPhiPhiMix_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 591           snprintf(histtitle,buffersize,"Mixed #Delta#phi-#Delta#phi %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 592           fHistDeltaPhiPhiMix[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinPhi3,fdPhiMin,fdPhiMax,fnBinPhi3,fdPhiMin,fdPhiMax);
 593           fHistDeltaPhiPhiMix[i][ipt][imult][iSign][imc]->Sumw2();
 594           fHistDeltaPhiPhiMix[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#phi_{1}");
 595           fHistDeltaPhiPhiMix[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#phi_{2}");
 596           fOutput->Add(fHistDeltaPhiPhiMix[i][ipt][imult][iSign][imc]);
 597
 598           snprintf(histname,buffersize,"fHistDeltaPhiPhiSS_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 599           snprintf(histtitle,buffersize,"Soft-Soft #Delta#phi-#Delta#phi %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 600           fHistDeltaPhiPhiSS[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinPhi3,fdPhiMin,fdPhiMax,fnBinPhi3,fdPhiMin,fdPhiMax);
 601           fHistDeltaPhiPhiSS[i][ipt][imult][iSign][imc]->Sumw2();
 602           fHistDeltaPhiPhiSS[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#phi_{1}");
 603           fHistDeltaPhiPhiSS[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#phi_{2}");
 604           fOutput->Add(fHistDeltaPhiPhiSS[i][ipt][imult][iSign][imc]);
 605
 606           snprintf(histname,buffersize,"fHistDeltaEtaEta_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 607           snprintf(histtitle,buffersize,"Raw #Delta#eta-#Delta#eta %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 608           fHistDeltaEtaEta[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinEta3,-etaEdge,etaEdge,fnBinEta3,-etaEdge,etaEdge);
 609           fHistDeltaEtaEta[i][ipt][imult][iSign][imc]->Sumw2();
 610           fHistDeltaEtaEta[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#eta_{1}");
 611           fHistDeltaEtaEta[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#eta_{2}");
 612           fOutput->Add(fHistDeltaEtaEta[i][ipt][imult][iSign][imc]);
 613
 614     snprintf(histname,buffersize,"fHistDeltaEtaEtaMix_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 615           snprintf(histtitle,buffersize,"Mixed #Delta#eta-#Delta#eta %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 616           fHistDeltaEtaEtaMix[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinEta3,-etaEdge,etaEdge,fnBinEta3,-etaEdge,etaEdge);
 617           fHistDeltaEtaEtaMix[i][ipt][imult][iSign][imc]->Sumw2();
 618           fHistDeltaEtaEtaMix[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#eta_{1}");
 619           fHistDeltaEtaEtaMix[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#eta_{2}");
 620           fOutput->Add(fHistDeltaEtaEtaMix[i][ipt][imult][iSign][imc]);
 621
 622     snprintf(histname,buffersize,"fHistDeltaEtaEtaSS_P%dp%d_C%d%s%s",i,ipt,imult,cmc1[imc],sign31[iSign]);
 623           snprintf(histtitle,buffersize,"Soft-Soft #Delta#eta-#Delta#eta %3.1f<p_{T}^{Trig}<%3.1f %3.2f<p_{T}^{Assoc}<%3.2f %dMult%d%s%s",fPtTrigArray[i],fPtTrigArray[i+1],fPtAssoc3Array1[ipt],fPtAssoc3Array2[ipt],fCentArrayMin[imult],fCentArrayMax[imult],sign32[iSign],cmc2[imc]);
 624           fHistDeltaEtaEtaSS[i][ipt][imult][iSign][imc]=new TH2F(histname,histtitle,fnBinEta3,-etaEdge,etaEdge,fnBinEta3,-etaEdge,etaEdge);
 625           fHistDeltaEtaEtaSS[i][ipt][imult][iSign][imc]->Sumw2();
 626           fHistDeltaEtaEtaSS[i][ipt][imult][iSign][imc]->GetXaxis()->SetTitle("#Delta#eta_{1}");
 627           fHistDeltaEtaEtaSS[i][ipt][imult][iSign][imc]->GetYaxis()->SetTitle("#Delta#eta_{2}");
 628           fOutput->Add(fHistDeltaEtaEtaSS[i][ipt][imult][iSign][imc]);
 629
 630         }//iSign
 631       }//associated pt (ipt)
 632       }//pt loop (i)
 633     }//centrality loop (imult)
 634   }//imc
 635   if(fDEBUG)Printf("OutPut Created");
 636 }//CreateOutputObjects
 637
 638 Int_t AliAnalysisTaskDiHadron::CheckVertex(const AliESDEvent *rESD){
 639   //checks whether the vertex passes cuts
 640   Int_t rGood=-1;
 641   Float_t vtx[3];
 642   vtx[0]=rESD->GetPrimaryVertex()->GetX();
 643   vtx[1]=rESD->GetPrimaryVertex()->GetY();
 644   vtx[2]=rESD->GetPrimaryVertex()->GetZ();
 645   if((vtx[0]*vtx[0]+vtx[1]*vtx[1])<9) rGood=0; //vertex out of beam pipe
 646   if(fabs(vtx[2])<fVzCut)rGood=0;//Vertex Z cut
 647   if(fDEBUG)Printf("vtxZ %f",vtx[2]);
 648   for(int i=0;i<fNVertexBins;i++){
 649     if(vtx[2]>fVertexArray[i]&&vtx[2]<=fVertexArray[i+1]&&rGood==0)rGood=i;
 650   }
 651   return rGood;
 652 }
 653
 654 Int_t AliAnalysisTaskDiHadron::CheckTrigger(const AliESDEvent *rESD){
 655   //checks whether the trigger passes cuts
 656   Int_t rGood=0;
 657   TString trigID=rESD->GetFiredTriggerClasses();
 658   int count=0;
 659   char trigID2[50];
 660   int stop=0;//in as a safety
 661
 662   for(int i=0;i<fNIDs;i++){
 663     if(stop==1)continue;
 664     for(int j=0;j<50;j++){
 665       if(fTrigIDArray[count]==',')trigID2[j]='\0';
 666       else if(fTrigIDArray[count]=='\0'){trigID2[j]='\0';stop=1;}
 667       else trigID2[j]=fTrigIDArray[count];
 668       count++;
 669       if(trigID2[j]=='\0') break;
 670       }
 671       if(trigID.Contains(trigID2)) rGood=1;
 672   }
 673     return rGood;
 674 }
 675
 676 Int_t AliAnalysisTaskDiHadron::TrackCuts(const AliESDEvent *rESD, Float_t *rPt, Float_t *rEta, Float_t *rPhi, Short_t *rCharge, Float_t *rEff, Int_t **rPtAssoc3, Int_t *rNPtAssoc3, Int_t *rGoodTracks){
 677     //fills arrays with all of the tracks passing cuts
 678   rGoodTracks[0]=0;
 679   Int_t lead=0;
 680   Float_t leadPt=0;
 681   Int_t rTrack=fESD->GetNumberOfTracks();
 682   Float_t sPt, sEta, sPhi, sChi, sb[2], sbCov[3];
 683   Int_t sNcls, sNclsF, sITScls;
 684   Short_t sCharge;
 685   for(int iTrack=0;iTrack<rTrack;iTrack++){
 686     AliESDtrack *eSDtrack=rESD->GetTrack(iTrack);
 687     const AliExternalTrackParam *conTrack = eSDtrack->GetConstrainedParam();
 688     if(!conTrack)continue;
 689     sPt=conTrack->Pt();
 690     //if(fDEBUG)Printf("Pt%f",rPt);
 691     sEta=conTrack->Eta();
 692     sPhi=conTrack->Phi();
 693     sCharge=conTrack->Charge();
 694     if(sPhi<fdPhiMin)sPhi+=2*fPi;
 695     if(sPhi>fdPhiMax)sPhi-=2*fPi;
 696     if(sPt<fMinPtAssoc||sPt>fMaxPtAssoc)continue;//set Pt range
 697     if(fabs(sEta)>fEtaCut)continue;//set Eta Range
 698     if(!sCharge)continue;
 699     sNcls=eSDtrack->GetTPCNcls();
 700     //if(fDEBUG)Printf("NCLS%d",sNcls);
 701     if(sNcls<fMinClustersTPC)continue;
 702     sNclsF=eSDtrack->GetTPCNclsF();
 703     if((1.0*sNcls/sNclsF)<fMinClusterRatio)continue;//Clusters fit/ Possible
 704     sChi=(eSDtrack->GetTPCchi2())/sNcls;
 705     if(sChi>fMaxTPCchi2)continue;
 706     sITScls=eSDtrack->GetNcls(0);
 707     if(sITScls<fMinClustersITS)continue;
 708     eSDtrack->GetImpactParameters(sb,sbCov);
 709     if(!fDCA2D&&(sb[0]*sb[0]+sb[1]*sb[1])>(fMaxDCA*fMaxDCA))continue;//DCA cut
 710     if(fDCA2D==1&&(sb[0]*sb[0]/fMaxDCAXY/fMaxDCAXY+sb[1]*sb[1]/fMaxDCAZ/fMaxDCAZ)>1)continue;
 711     if(fDCA2D==2&&(0.35+0.42*std::pow(double(sPt),-0.9))<(sb[0]*sb[0]))continue;
 712     if(eSDtrack->GetKinkIndex(0)>0)continue;//removes kinked tracks
 713     if(!eSDtrack->GetStatus()&&AliESDtrack::kTPCrefit&&fTPCRefit)continue;//refit in TPC
 714     if((fITSRefit==1||(fITSRefit==2&&sPt>5))&&!eSDtrack->GetStatus()&&AliESDtrack::kITSrefit)continue;//refit of its tracks either for none,all, or >5 GeV/c
 715     if(fSPDCut&&!eSDtrack->HasPointOnITSLayer(0)&&!eSDtrack->HasPointOnITSLayer(1))continue;
 716     rPt[rGoodTracks[0]]=sPt;
 717     rEta[rGoodTracks[0]]=sEta;
 718     rPhi[rGoodTracks[0]]=sPhi;
 719     rCharge[rGoodTracks[0]]=sCharge;
 720     if(fEfficiencyCorr){
 721     if(sPt<fEffFitPt)rEff[rGoodTracks[0]]=1./fFitLow->Eval(sPt);
 722     else rEff[rGoodTracks[0]]=1./fFitHigh->Eval(sPt);
 723     }
 724     else rEff[rGoodTracks[0]]=1;
 725     if(sPt>leadPt)lead=rGoodTracks[0];
 726     //rPtAssoc3[rGoodTracks[0]]=new Int_t [10];
 727     rNPtAssoc3[rGoodTracks[0]]=0;
 728     for(int apt3=0;apt3<fNAPt3Bins;apt3++){
 729       if(sPt<fPtAssoc3Array2[apt3]&&sPt>=fPtAssoc3Array1[apt3]){
 730         rPtAssoc3[rGoodTracks[0]][rNPtAssoc3[rGoodTracks[0]]]=apt3;
 731         rNPtAssoc3[rGoodTracks[0]]++;
 732       }
 733     }
 734
 735     rGoodTracks[0]++;
 736
 737   }
 738   return lead;
 739 }
 740
 741 Int_t AliAnalysisTaskDiHadron::TrackCutsMC(AliMCEvent *rMC, Float_t *rPt, Float_t *rEta, Float_t *rPhi, Short_t *rCharge, Float_t *rEff, Int_t **rPtAssoc3, Int_t *rNPtAssoc3, Int_t *rGoodTracks){
 742 //Fills Arrays of MC particles
 743   rGoodTracks[1]=0;
 744   AliStack *rStack=rMC->Stack();
 745   Int_t rTrack=rStack->GetNtrack();
 746   Float_t sPt, sEta, sPhi;
 747   Short_t sCharge;
 748   Int_t lead=0;
 749   Float_t leadPt=0;
 750   for(int iTrack=0;iTrack<rTrack;iTrack++){
 751     TParticle *rParticle=rStack->Particle(iTrack);
 752     sPt=rParticle->Pt();
 753     //if(fDEBUG)Printf("MCPt%f",rPt);
 754     sEta=rParticle->Eta();
 755     sPhi=rParticle->Phi();
 756     sCharge=rMC->GetTrack(iTrack)->Charge();
 757     if(sPhi<fdPhiMin)sPhi+=2*fPi;
 758     if(sPhi>fdPhiMax)sPhi-=2*fPi;
 759     if(sPt<fMinPtAssoc||sPt>fMaxPtAssoc)continue;//set Pt range
 760     if(fabs(sEta)>fEtaCut)continue;//set Eta Range
 761     if(!rStack->IsPhysicalPrimary(iTrack))continue;//primary particles only
 762     if(!sCharge)continue;//only charged particles kept
 763     rPt[rGoodTracks[1]]=sPt;
 764     rEta[rGoodTracks[1]]=sEta;
 765     rPhi[rGoodTracks[1]]=sPhi;
 766     rCharge[rGoodTracks[1]]=sCharge;
 767     rEff[rGoodTracks[1]]=1;
 768     if(sPt>leadPt)lead=rGoodTracks[1];
 769     rNPtAssoc3[rGoodTracks[1]]=0;
 770     for(int apt3=0;apt3<fNAPt3Bins;apt3++){
 771       if(sPt<fPtAssoc3Array2[apt3]&&sPt>=fPtAssoc3Array1[apt3]){
 772         rPtAssoc3[rGoodTracks[1]][rNPtAssoc3[rGoodTracks[1]]]=apt3;
 773         rNPtAssoc3[rGoodTracks[1]]++;
 774       }
 775     }
 776     rGoodTracks[1]++;
 777   }
 778   return lead;
 779 }
 780 //------------------------------------------------------------
 781 void AliAnalysisTaskDiHadron::Exec(Option_t *)
 782 {
 783   //Main executable
 784   if(fDEBUG)Printf("Exec");
 785
 786   const int nTPtBins=fNTPtBins;
 787   const int nCentBins=fNCentBins;
 788   for(int ievent=0;ievent<=1;ievent++){
 789
 790   if(!fESD&&ievent==0){
 791     if(fDEBUG)Printf("Error: fESD not found");
 792     break;
 793   }
 794   if(!fMC&&ievent==1){
 795     break;
 796   }
 797   if(ievent==1&&!fMCHistos)break;//break out if MC event and we don't have fill of those set
 798   //Secondary check
 799   if(ievent==0){
 800     if(fESD->GetNumberOfTracks()<=0){
 801       if(fDEBUG)Printf("Error: no tracks");
 802       break;
 803     }
 804   }
 805   //The previous check doesn't seem to work as a fMC is bad not NULL
 806   if(ievent==1){
 807     if(fMC->GetNumberOfTracks()<=0){
 808       if(fDEBUG)Printf("<=0 MCTracks");
 809       break;
 810     }
 811   }
 812
 813   //Check for Trigger only on real data
 814   if(!fMC){
 815     if(!CheckTrigger(fESD)) break;
 816   }
 817   //I'll only cut on the reconstructed vertex since these are the events that will be used
 818   int vertexBin;
 819   vertexBin=CheckVertex(fESD);
 820   //else vertexBin=CheckVertex(fMC);
 821   if(vertexBin<0)break;
 822
 823   Int_t nGoodTracks[2]={0,0}, nTriggers[nTPtBins][nCentBins][2];
 824   Int_t nTrack;
 825   if(!ievent)nTrack=fESD->GetNumberOfTracks();
 826   else nTrack=fMC->Stack()->GetNtrack();
 827
 828   Float_t tdPhi, tdEta, tXE;
 829   Float_t tdPhi2, tdEta2;
 830   ftPhi=new Float_t [nTrack];
 831   ftEta=new Float_t [nTrack];
 832   ftPt=new Float_t [nTrack];
 833   ftCharge=new Short_t [nTrack];
 834   ftEff=new Float_t [nTrack];
 835   ftPtAssoc3=new Int_t *[nTrack];
 836   for(int i=0;i<nTrack;i++){
 837     ftPtAssoc3[i]=new Int_t [10];
 838   }
 839   ftNPtAssoc3=new Int_t [nTrack];
 840   Short_t sign;
 841
 842   //will need to do something exta for the effieiency once it comes from embedding
 843   for(int i=0;i<fNTPtBins;i++){
 844     for(int c=0;c<fNCentBins;c++){
 845     nTriggers[i][c][ievent]=0;
 846     }
 847   }
 848   Int_t tMult=fESD->GetMultiplicity()->GetNumberOfTracklets();//I think this is the correct multiplicity to use
 849   if(fDEBUG)Printf("Mult%d",tMult);
 850
 851   //Decide what multiplicy bins are filled with this event, note set to max of 4 as I didn't think more then 2 overlapping bins likely at one time, easliy changed
 852   Int_t multArray[4]={0,0,0,0};
 853   Int_t maxArray=0;
 854   for(int imult=0;imult<fNCentBins;imult++){
 855     if(tMult>=fCentArrayMin[imult]&&tMult<fCentArrayMax[imult]){
 856       multArray[maxArray]=imult;
 857       maxArray++;
 858     }
 859   }
 860   if(fDEBUG)Printf("maxArray%d",maxArray);
 861   //Set Efficiency for the centrality bin (lowest bin used in array if multiple overlap)
 862   for(int ipar=0;ipar<fNFitLowParam;ipar++){
 863     fFitLow->SetParameter(ipar,fFitLowParam[multArray[0]*fNCentBins+ipar]);
 864   }
 865    for(int ipar=0;ipar<fNFitHighParam;ipar++){
 866     fFitHigh->SetParameter(ipar,fFitHighParam[multArray[0]*fNCentBins+ipar]);
 867   }
 868   fHistMult[ievent]->Fill(tMult);
 869   for(int c=0;c<maxArray;c++){fHistNEvents[multArray[c]][ievent]->Fill(0);}//count the number of events used
 870   Int_t leadPart;
 871
 872   //returns arrays filled up to nGoodTracks with tracks passing cuts
 873   if(!ievent)leadPart=TrackCuts(fESD,ftPt,ftEta,ftPhi,ftCharge,ftEff,ftPtAssoc3,ftNPtAssoc3,nGoodTracks);
 874   else leadPart=TrackCutsMC(fMC,ftPt,ftEta,ftPhi,ftCharge,ftEff,ftPtAssoc3,ftNPtAssoc3,nGoodTracks);
 875   int nearEta=0,NearXE=0;
 876   int nearEta2=0;
 877   if(fDEBUG)Printf("Track Loop");
 878   for(int iTrack=0;iTrack<nGoodTracks[ievent];iTrack++){
 879     if(fDEBUG)Printf("Track%d Pt%f",iTrack,ftPt[iTrack]);
 880     //if(ftPhi[iTrack]<fdPhiMin)ftPhi[iTrack]+=2*fPi;
 881     //if(ftPhi[iTrack]>fdPhiMax)ftPhi[iTrack]-=2*fPi;
 882     for(int c=0;c<maxArray;c++){
 883       fHistPt[multArray[c]][ievent]->Fill(ftPt[iTrack],ftEff[iTrack]);
 884     fHistPtEff[multArray[c]][ievent]->Fill(ftPt[iTrack]);
 885     fHistPhi[multArray[c]][ievent]->Fill(ftPhi[iTrack],ftPt[iTrack],ftEff[iTrack]);
 886     fHistPhiPt[multArray[c]][ievent]->Fill(ftPhi[iTrack],ftPt[iTrack],ftPt[iTrack]*ftEff[iTrack]);
 887     fHistEta[multArray[c]][ievent]->Fill(ftEta[iTrack],ftPt[iTrack],ftEff[iTrack]);
 888     fHistEtaPt[multArray[c]][ievent]->Fill(ftEta[iTrack],ftPt[iTrack],ftPt[iTrack]*ftEff[iTrack]);
 889     fHistPhiEta[multArray[c]][ievent]->Fill(ftPhi[iTrack],ftEta[iTrack],ftPt[iTrack],ftEff[iTrack]);
 890     fHistPhiEtaPt[multArray[c]][ievent]->Fill(ftPhi[iTrack],ftEta[iTrack],ftPt[iTrack],ftPt[iTrack]*ftEff[iTrack]);
 891     }
 892     for(int i=0;i<fNTPtBins;i++){
 893       if(ftPt[iTrack]>fPtTrigArray[i]&&ftPt[iTrack]<fPtTrigArray[i+1]&&fabs(ftEta[iTrack])<fTrigEtaCut){
 894         if(fDEBUG)Printf("In %fpt%f",fPtTrigArray[i],fPtTrigArray[i+1]);
 895         fHistMultTrig[i][ievent]->Fill(tMult);
 896         for(int c=0;c<maxArray;c++){
 897           nTriggers[i][multArray[c]][ievent]++;
 898           fHistNTrigger[multArray[c]][ievent]->Fill(i);
 899           fHistNTriggerPt[multArray[c]][ievent]->Fill(i,ftPt[iTrack]);
 900         }
 901         if(fDEBUG)Printf("Assiciated Particle Loop");
 902         for(int iTrack2=0;iTrack2<nGoodTracks[ievent];iTrack2++){
 903           if(iTrack==iTrack2) continue;
 904           if(ftPt[iTrack2]>ftPt[iTrack])continue;
 905           tdPhi=ftPhi[iTrack]-ftPhi[iTrack2];
 906           if(tdPhi<-fPi)tdPhi+=2*fPi;
 907           if(tdPhi>fPi)tdPhi-=2*fPi;
 908           if(fabs(tdPhi)<fNearPhiCut)nearEta=1;
 909           else nearEta=0;
 910           if(fabs(tdPhi)<fXECut)NearXE=1;
 911           else NearXE=0;
 912           if(fabs(tdPhi)<(fPi/2))tdEta=ftEta[iTrack]-ftEta[iTrack2];
 913           else tdEta=ftEta[iTrack]+ftEta[iTrack2];
 914           if(tdPhi<fdPhiMin)tdPhi+=2*fPi;
 915           if(tdPhi>fdPhiMax)tdPhi-=2*fPi;
 916           if((ftCharge[iTrack]<0&&ftCharge[iTrack2]<0)||(ftCharge[iTrack]>0&&ftCharge[iTrack2]>0))sign=1;
 917           else sign=2;
 918           if(fDEBUG) Printf("dPhi %f  dEta %f",tdPhi,tdEta);
 919           for(int c=0;c<maxArray;c++){//loop over multiplicity bins
 920             fHistPtTrig[i][multArray[c]][ievent]->Fill(ftPt[iTrack2],ftEff[iTrack2]);
 921             fHistPhiTrig[i][multArray[c]][ievent]->Fill(ftPhi[iTrack2],ftPt[iTrack2],ftEff[iTrack2]);
 922             fHistPhiTrigPt[i][multArray[c]][ievent]->Fill(ftPhi[iTrack2],ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 923             fHistEtaTrig[i][multArray[c]][ievent]->Fill(ftEta[iTrack2],ftPt[iTrack2],ftEff[iTrack2]);
 924             fHistEtaTrigPt[i][multArray[c]][ievent]->Fill(ftEta[iTrack2],ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 925
 926             fHistPhiEtaTrig[i][multArray[c]][ievent]->Fill(ftPhi[iTrack2],ftEta[iTrack2],ftPt[iTrack2],ftEff[iTrack2]);
 927             fHistPhiEtaTrigPt[i][multArray[c]][ievent]->Fill(ftPhi[iTrack2],ftEta[iTrack2],ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 928             fHistDeltaPhi[i][multArray[c]][0][ievent]->Fill(tdPhi,ftPt[iTrack2],ftEff[iTrack2]);
 929             fHistDeltaPhiPt[i][multArray[c]][0][ievent]->Fill(tdPhi,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 930             fHistDeltaPhi[i][multArray[c]][sign][ievent]->Fill(tdPhi,ftPt[iTrack2],ftEff[iTrack2]);
 931             fHistDeltaPhiPt[i][multArray[c]][sign][ievent]->Fill(tdPhi,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 932
 933             if(nearEta){
 934             fHistDeltaEtaN[i][multArray[c]][0][ievent]->Fill(tdEta,ftPt[iTrack2],ftEff[iTrack2]);
 935             fHistDeltaEtaNPt[i][multArray[c]][0][ievent]->Fill(tdEta,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 936             fHistDeltaEtaN[i][multArray[c]][sign][ievent]->Fill(tdEta,ftPt[iTrack2],ftEff[iTrack2]);
 937             fHistDeltaEtaNPt[i][multArray[c]][sign][ievent]->Fill(tdEta,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 938             }
 939             else{
 940             fHistDeltaEtaA[i][multArray[c]][0][ievent]->Fill(tdEta,ftPt[iTrack2],ftEff[iTrack2]);
 941             fHistDeltaEtaAPt[i][multArray[c]][0][ievent]->Fill(tdEta,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 942             fHistDeltaEtaA[i][multArray[c]][sign][ievent]->Fill(tdEta,ftPt[iTrack2],ftEff[iTrack2]);
 943             fHistDeltaEtaAPt[i][multArray[c]][sign][ievent]->Fill(tdEta,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 944             }
 945             fHistDeltaPhiEta[i][multArray[c]][ievent]->Fill(tdPhi,tdEta,ftPt[iTrack2],ftEff[iTrack2]);
 946             fHistDeltaPhiEtaPt[i][multArray[c]][ievent]->Fill(tdPhi,tdEta,ftPt[iTrack2],ftPt[iTrack2]*ftEff[iTrack2]);
 947
 948             //only fill these if trigger particle is the leading particle
 949             if(iTrack==leadPart){
 950               if(NearXE){
 951                 tXE=ftPt[iTrack2]*cos(tdPhi)/ftPt[iTrack];
 952                 fHistXEN[i][multArray[c]][ievent]->Fill(tXE,ftEff[iTrack2]);
 953               }
 954               else{
 955                 tXE=ftPt[iTrack2]*cos(tdPhi+fPi)/ftPt[iTrack];
 956                 fHistXEA[i][multArray[c]][ievent]->Fill(tXE,ftEff[iTrack2]);
 957               }
 958             }
 959
 960           }//Centrality loop (c)
 961
 962           //3-particle Correlations
 963           for(int iTrack3=0;iTrack3<nGoodTracks[ievent];iTrack3++){
 964             if(iTrack2==iTrack3)continue;
 965             if(ftPt[iTrack3]>ftPt[iTrack])continue;
 966             tdPhi2=ftPhi[iTrack]-ftPhi[iTrack3];
 967             if(tdPhi2<-fPi)tdPhi2+=2*fPi;
 968             if(tdPhi2>fPi)tdPhi2-=2*fPi;
 969             if(fabs(tdPhi2)<fNearPhiCut&&nearEta==1)nearEta2=1;
 970             else nearEta2=0;
 971             //if(fabs(tdPhi)<fXECut)NearXE=1;
 972             //else NearXE=0;
 973             if(fabs(tdPhi2)<(fPi/2))tdEta2=ftEta[iTrack]-ftEta[iTrack3];
 974             else tdEta2=ftEta[iTrack]+ftEta[iTrack3];
 975             if(tdPhi2<fdPhiMin)tdPhi2+=2*fPi;
 976             if(tdPhi2>fdPhiMax)tdPhi2-=2*fPi;
 977             // if((ftCharge[iTrack]<0&&ftCharge[iTrack2]<0)||(ftCharge[iTrack]>0&&ftCharge[iTrack2]>0))sign=1;
 978             if((ftCharge[iTrack]<0&&ftCharge[iTrack2]<0&&ftCharge[iTrack3]<0)||(ftCharge[iTrack]>0&&ftCharge[iTrack2]>0&&ftCharge[iTrack3]>0))sign=1;
 979             else if((ftCharge[iTrack3]<0&&ftCharge[iTrack2]<0)||(ftCharge[iTrack3]>0&&ftCharge[iTrack2]>0))sign=2;
 980             else sign=3;
 981             for(int e=0;e<ftNPtAssoc3[iTrack2];e++){//check associated pT bin
 982               for(int f=0;f<ftNPtAssoc3[iTrack3];f++){
 983                 if(ftPtAssoc3[iTrack2][e]==ftPtAssoc3[iTrack3][f]){
 984                   for(int c=0;c<maxArray;c++){//loop over multiplicity bins
 985                     fHistDeltaPhiPhi[i][ftPtAssoc3[iTrack2][e]][multArray[c]][0][ievent]->Fill(tdPhi,tdPhi2,ftEff[iTrack2]*ftEff[iTrack3]);
 986                     fHistDeltaPhiPhi[i][ftPtAssoc3[iTrack2][e]][multArray[c]][sign][ievent]->Fill(tdPhi2,tdPhi,ftEff[iTrack2]*ftEff[iTrack3]);
 987
 988
 989                     if(nearEta2){
 990                       fHistDeltaEtaEta[i][ftPtAssoc3[iTrack2][e]][multArray[c]][0][ievent]->Fill(tdEta,tdEta2,ftEff[iTrack2]*ftEff[iTrack3]);
 991                       fHistDeltaEtaEta[i][ftPtAssoc3[iTrack2][e]][multArray[c]][sign][ievent]->Fill(tdEta,tdEta2,ftEff[iTrack2]*ftEff[iTrack3]);
 992                     }
 993                   }//multiplicity loop (c)
 994                 }
 995               }
 996             }//track checking loops
 997           }//iTrack3
 998         }//iTrack2 (associated track loop)
 999
1000         if(fDEBUG)Printf("Mixed Event Loop");
1001         for(int c=0;c<maxArray;c++){
1002           int d=multArray[c];//Centrality bin we are in
1003           if(fMixEnd[d][vertexBin][ievent]>=0){//check if there are any mixed events for this bin
1004             for(int imix=0;imix<=fMixEnd[d][vertexBin][ievent];imix++){//loop over the stored mixed events
1005               fHistNMix[d][ievent]->Fill(i);
1006               for(int iTrack2=0;iTrack2<fMixTrack[imix][d][vertexBin][ievent];iTrack2++){
1007                 if(ftPt[iTrack]<fMPt[imix][d][vertexBin][ievent][iTrack2])continue;
1008                 tdPhi=ftPhi[iTrack]-fMPhi[imix][d][vertexBin][ievent][iTrack2];
1009                 if(tdPhi<-fPi)tdPhi+=2*fPi;
1010                 if(tdPhi>fPi)tdPhi-=2*fPi;
1011                 if(fabs(tdPhi)<fNearPhiCut)nearEta=1;
1012                 else nearEta=0;
1013                 if(fabs(tdPhi)<fXECut)NearXE=1;
1014                 else NearXE=0;
1015                 if(fabs(tdPhi)<(fPi/2))tdEta=ftEta[iTrack]-fMEta[imix][d][vertexBin][ievent][iTrack2];
1016                 else tdEta=ftEta[iTrack]+fMEta[imix][d][vertexBin][ievent][iTrack2];
1017                 if(tdPhi<fdPhiMin)tdPhi+=2*fPi;
1018                 if(tdPhi>fdPhiMax)tdPhi-=2*fPi;
1019                 if((ftCharge[iTrack]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0)||(ftCharge[iTrack]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0))sign=1;
1020                 else sign=2;
1021
1022                 fHistDeltaPhiMix[i][d][0][ievent]->Fill(tdPhi,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1023                 fHistDeltaPhiMixPt[i][d][0][ievent]->Fill(tdPhi,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1024                 fHistDeltaPhiMix[i][d][sign][ievent]->Fill(tdPhi,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1025                 fHistDeltaPhiMixPt[i][d][sign][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1026                 if(nearEta){
1027                   fHistDeltaEtaNMix[i][d][0][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1028                 fHistDeltaEtaNMixPt[i][d][0][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1029                 fHistDeltaEtaNMix[i][d][sign][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1030                 fHistDeltaEtaNMixPt[i][d][sign][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1031                 }
1032                 else{
1033                   fHistDeltaEtaAMix[i][d][0][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1034                 fHistDeltaEtaAMixPt[i][d][0][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1035                 fHistDeltaEtaAMix[i][d][sign][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1036                 fHistDeltaEtaAMixPt[i][d][sign][ievent]->Fill(tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1037                 }
1038
1039                 fHistDeltaPhiEtaMix[i][d][ievent]->Fill(tdPhi,tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMEff[imix][d][vertexBin][ievent][iTrack2]);
1040                 fHistDeltaPhiEtaMixPt[i][d][ievent]->Fill(tdPhi,tdEta,fMPt[imix][d][vertexBin][ievent][iTrack2],fMPt[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack2]);
1041
1042                 if(iTrack==leadPart){
1043                   if(NearXE){
1044                     tXE=fMPt[imix][d][vertexBin][ievent][iTrack2]*cos(tdPhi)/ftPt[iTrack];
1045                     fHistXENMix[i][d][ievent]->Fill(tXE,fMEff[imix][d][vertexBin][ievent][iTrack2]);
1046                   }
1047                   else{
1048                     tXE=fMPt[imix][d][vertexBin][ievent][iTrack2]*cos(tdPhi+fPi)/ftPt[iTrack];
1049                     fHistXEAMix[i][multArray[c]][ievent]->Fill(tXE,fMEff[imix][d][vertexBin][ievent][iTrack2]);
1050                   }
1051                 }
1052                 //3-particle correlation soft-soft term (both associated from the same event)
1053                 for(int iTrack3=0;iTrack3<fMixTrack[imix][d][vertexBin][ievent];iTrack3++){
1054                   if(iTrack3==iTrack2)continue;
1055                   if(ftPt[iTrack]<fMPt[imix][d][vertexBin][ievent][iTrack3])continue;
1056                 tdPhi2=ftPhi[iTrack]-fMPhi[imix][d][vertexBin][ievent][iTrack3];
1057                 if(tdPhi2<-fPi)tdPhi2+=2*fPi;
1058                 if(tdPhi2>fPi)tdPhi2-=2*fPi;
1059                 if(fabs(tdPhi2)<fNearPhiCut&&nearEta)nearEta2=1;
1060                 else nearEta2=0;
1061                 if(fabs(tdPhi2)<(fPi/2))tdEta2=ftEta[iTrack]-fMEta[imix][d][vertexBin][ievent][iTrack3];
1062                 else tdEta2=ftEta[iTrack]+fMEta[imix][d][vertexBin][ievent][iTrack3];
1063                 if(tdPhi2<fdPhiMin)tdPhi2+=2*fPi;
1064                 if(tdPhi2>fdPhiMax)tdPhi2-=2*fPi;
1065                 //if((ftCharge[iTrack]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0)||(ftCharge[iTrack]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0))sign=1;
1066                 //else sign=2;
1067                 if((ftCharge[iTrack]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack3]<0)||(ftCharge[iTrack]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack3]>0))sign=1;
1068                 else if((fMCharge[imix][d][vertexBin][ievent][iTrack3]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0)||(fMCharge[imix][d][vertexBin][ievent][iTrack3]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0))sign=2;
1069                 else sign=3;
1070                 for(int e=0;e<fMNPtAssoc3[imix][d][vertexBin][ievent][iTrack2];e++){//check associated pT bin
1071                   for(int f=0;f<fMNPtAssoc3[imix][d][vertexBin][ievent][iTrack3];f++){
1072                     if(fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]==fMPtAssoc3[imix][d][vertexBin][ievent][f][iTrack3]){
1073                       fHistDeltaPhiPhiSS[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdPhi,tdPhi2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack3]);
1074                       fHistDeltaPhiPhiSS[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdPhi,tdPhi2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack3]);
1075
1076                       if(nearEta2){
1077                         fHistDeltaEtaEtaSS[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdEta,tdEta2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack3]);
1078                         fHistDeltaEtaEtaSS[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdEta,tdEta2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix][d][vertexBin][ievent][iTrack3]);
1079                       }//near-side
1080                     }
1081                   }
1082                 }//associated pt bin
1083                 }//iTrack3
1084
1085                 //3-particle mixed event (associated from different events)
1086                 //for(int imix2=0;imix2<=fMixEnd[d][vertexBin][ievent];imix2++){//loop over the stored mixed events
1087                 //if(imix2==imix)continue;
1088                   int imix2=imix+1;
1089                   if(imix2>=fMixEnd[d][vertexBin][ievent])imix2=0;
1090                   if(imix2==imix)continue;//will kill it when there is only 1 mixed event (remember to scale by 1 less then the number of mixed events in others (number of mixed-mixed is 2*(mixed-1)))
1091                   for(int iTrack3=0;iTrack3<fMixTrack[imix2][d][vertexBin][ievent];iTrack3++){
1092                     if(ftPt[iTrack]<fMPt[imix2][d][vertexBin][ievent][iTrack3])continue;
1093                     tdPhi2=ftPhi[iTrack]-fMPhi[imix2][d][vertexBin][ievent][iTrack3];
1094                     if(tdPhi2<-fPi)tdPhi2+=2*fPi;
1095                     if(tdPhi2>fPi)tdPhi2-=2*fPi;
1096                     if(fabs(tdPhi2)<fNearPhiCut&&nearEta)nearEta2=1;
1097                     else nearEta2=0;
1098                     if(fabs(tdPhi2)<(fPi/2))tdEta2=ftEta[iTrack]-fMEta[imix2][d][vertexBin][ievent][iTrack3];
1099                     else tdEta2=ftEta[iTrack]+fMEta[imix2][d][vertexBin][ievent][iTrack3];
1100                     if(tdPhi2<fdPhiMin)tdPhi2+=2*fPi;
1101                     if(tdPhi2>fdPhiMax)tdPhi2-=2*fPi;
1102                     //if((ftCharge[iTrack]<0&&fMCharge[imix2][d][vertexBin][ievent][iTrack2]<0)||(ftCharge[iTrack]>0&&fMCharge[imix2][d][vertexBin][ievent][iTrack2]>0))sign=1;
1103                     //else sign=2;
1104                     if((ftCharge[iTrack]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0&&fMCharge[imix2][d][vertexBin][ievent][iTrack3]<0)||(ftCharge[iTrack]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0&&fMCharge[imix2][d][vertexBin][ievent][iTrack3]>0))sign=1;
1105                     else if((fMCharge[imix2][d][vertexBin][ievent][iTrack3]<0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]<0)||(fMCharge[imix2][d][vertexBin][ievent][iTrack3]>0&&fMCharge[imix][d][vertexBin][ievent][iTrack2]>0))sign=2;
1106                     else sign=3;
1107                     for(int e=0;e<fMNPtAssoc3[imix][d][vertexBin][ievent][iTrack2];e++){//check associated pT bin
1108                       for(int f=0;f<fMNPtAssoc3[imix2][d][vertexBin][ievent][iTrack3];f++){
1109                         if(fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]==fMPtAssoc3[imix2][d][vertexBin][ievent][f][iTrack3]){
1110                           fHistDeltaPhiPhiMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdPhi,tdPhi2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1111                           fHistDeltaPhiPhiMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdPhi,tdPhi2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1112                           fHistDeltaPhiPhiMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdPhi2,tdPhi,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);//free factor of 2 in statistics
1113                           fHistDeltaPhiPhiMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdPhi2,tdPhi,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1114                       if(nearEta2){
1115                         fHistDeltaEtaEtaMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdEta,tdEta2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1116                         fHistDeltaEtaEtaMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdEta,tdEta2,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1117                         fHistDeltaEtaEtaMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][0][ievent]->Fill(tdEta2,tdEta,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1118                         fHistDeltaEtaEtaMix[i][fMPtAssoc3[imix][d][vertexBin][ievent][e][iTrack2]][d][sign][ievent]->Fill(tdEta2,tdEta,fMEff[imix][d][vertexBin][ievent][iTrack2]*fMEff[imix2][d][vertexBin][ievent][iTrack3]);
1119                       }//near-side
1120                         }
1121                       }
1122                     }//associated pt bin
1123                   }//iTrack3
1124
1125               }//iTrack2
1126               }//imix
1127             }//fMixEnd
1128           }//Centrality bins (c)
1129         }//pt trig cuts
1130       }//i Pt Trig
1131     }//itrack
1132
1133   //now store this event for mixing (using these dynamic arrays instead of the other fix to save memory)
1134     if(fDEBUG)Printf("Store Event For Mixing");
1135   for(int c=0;c<maxArray;c++){//loops over centrality bins
1136     int d=multArray[c];//too many nested arrays looked confusing d=which centrality bin
1137     if(fMixEnd[d][vertexBin][ievent]<(fNMix-1))fMixEnd[d][vertexBin][ievent]++;
1138     if(fMixPointer[d][vertexBin][ievent]<(fNMix-1))fMixPointer[d][vertexBin][ievent]++;
1139     else fMixPointer[d][vertexBin][ievent]=0;
1140     int e=fMixPointer[d][vertexBin][ievent];//nested arrays (e is event number in pool)
1141     delete [] fMPt[e][d][vertexBin][ievent];
1142     delete [] fMPhi[e][d][vertexBin][ievent];
1143     delete [] fMEta[e][d][vertexBin][ievent];
1144     delete [] fMCharge[e][d][vertexBin][ievent];
1145     delete [] fMEff[e][d][vertexBin][ievent];
1146     delete [] fMNPtAssoc3[e][d][vertexBin][ievent];
1147     for(int jj=0;jj<10;jj++){
1148     delete [] fMPtAssoc3[e][d][vertexBin][ievent][jj];
1149     }
1150     fMPt[e][d][vertexBin][ievent]=new Float_t [nGoodTracks[ievent]];
1151     fMPhi[e][d][vertexBin][ievent]=new Float_t [nGoodTracks[ievent]];
1152     fMEta[e][d][vertexBin][ievent]=new Float_t [nGoodTracks[ievent]];
1153     fMCharge[e][d][vertexBin][ievent]=new Short_t [nGoodTracks[ievent]];
1154     fMEff[e][d][vertexBin][ievent]=new Float_t [nGoodTracks[ievent]];
1155     fMixTrack[e][d][vertexBin][ievent]=nGoodTracks[ievent];
1156     fMNPtAssoc3[e][d][vertexBin][ievent]=new Short_t [nGoodTracks[ievent]];
1157     for(int jj=0;jj<10;jj++){
1158     fMPtAssoc3[e][d][vertexBin][ievent][jj]=new Short_t [nGoodTracks[ievent]];
1159     }
1160
1161     for(int iTrack=0;iTrack<nGoodTracks[ievent];iTrack++){
1162       fMPt[e][d][vertexBin][ievent][iTrack]=ftPt[iTrack];
1163       fMPhi[e][d][vertexBin][ievent][iTrack]=ftPhi[iTrack];
1164       fMEta[e][d][vertexBin][ievent][iTrack]=ftEta[iTrack];
1165       fMCharge[e][d][vertexBin][ievent][iTrack]=ftCharge[iTrack];
1166       fMEff[e][d][vertexBin][ievent][iTrack]=ftEff[iTrack];
1167       fMNPtAssoc3[e][d][vertexBin][ievent][iTrack]=ftNPtAssoc3[iTrack];
1168       // fMPtAssoc3[e][d][vertexBin][ievent][iTrack]=new Int_t [ftNPtAssoc3[iTrack]];
1169       for(int jj=0;jj<ftNPtAssoc3[iTrack];jj++){
1170         //if(fDEBUG) Printf("%d",ftPtAssoc3[iTrack][jj]);
1171         fMPtAssoc3[e][d][vertexBin][ievent][jj][iTrack]=ftPtAssoc3[iTrack][jj];
1172       }
1173
1174     }//iTracks
1175   }//Centrality (c)
1176   }//ievent
1177   //track=0;
1178   //track2=0;
1179
1180   PostData(0, fOutput);
1181   //get rid of these arrays from memory
1182   delete [] ftPhi;
1183   delete [] ftEta;
1184   delete [] ftPt;
1185   delete [] ftCharge;
1186   delete [] ftEff;
1187   delete [] ftNPtAssoc3;
1188   delete [] ftPtAssoc3;
1189   ftPhi=NULL;
1190   ftEta=NULL;
1191   ftPt=NULL;
1192   ftCharge=NULL;
1193   ftEff=NULL;
1194   ftNPtAssoc3=NULL;
1195   ftPtAssoc3=NULL;
1196
1197 }//Exec
1198
1199 //---------------------------------------------------
1200 void AliAnalysisTaskDiHadron::Terminate(Option_t *){
1201   //Terminates the code, frees up memory
1202   for(int ii=0;ii<fNMix;ii++){
1203     for(int cc=0;cc<fNCentBins;cc++){
1204       for(int vtx=0;vtx<fNVertexBins;vtx++){
1205         for(int jj=0;jj<2;jj++){
1206           delete [] fMPt[ii][cc][vtx][jj];
1207           delete [] fMPhi[ii][cc][vtx][jj];
1208           delete [] fMEta[ii][cc][vtx][jj];
1209           delete [] fMCharge[ii][cc][vtx][jj];
1210           delete [] fMEff[ii][cc][vtx][jj];
1211           delete [] fMNPtAssoc3[ii][cc][vtx][jj];
1212           for(int qq=0;qq<10;qq++){
1213             delete [] fMPtAssoc3[ii][cc][vtx][jj][qq];
1214             fMPtAssoc3[ii][cc][vtx][jj][qq]=NULL;
1215           }
1216           fMPt[ii][cc][vtx][jj]=NULL;
1217           fMPhi[ii][cc][vtx][jj]=NULL;
1218           fMEta[ii][cc][vtx][jj]=NULL;
1219           fMCharge[ii][cc][vtx][jj]=NULL;
1220           fMEff[ii][cc][vtx][jj]=NULL;
1221           fMNPtAssoc3[ii][cc][vtx][jj]=NULL;
1222
1223         }
1224       }
1225     }
1226   }
1227   delete [] fFitLowParam;
1228   delete [] fFitHighParam;
1229   delete [] fPtTrigArray;
1230   delete [] fPtAssocArray;
1231   delete [] fPtAssoc3Array1;
1232   delete [] fPtAssoc3Array2;
1233   delete [] fCentArrayMin;
1234   delete [] fCentArrayMax;
1235   delete [] fXEArray;
1236   fFitLowParam=NULL;
1237   fFitHighParam=NULL;
1238   fPtTrigArray=NULL;
1239   fPtAssocArray=NULL;
1240   fPtAssoc3Array1=NULL;
1241   fPtAssoc3Array2=NULL;
1242   fCentArrayMin=NULL;
1243   fCentArrayMax=NULL;
1244   Printf("Terminate AliAnalysisTaskDiHadron");
1245 }