Fixed issue in PWGCF/Correlations/DPhi/AliAnalysisTaskDptDptQA.cxx
authorafestant <andrea.festanti@cern.ch>
Fri, 20 Jun 2014 17:59:52 +0000 (19:59 +0200)
committerhristov <Peter.Hristov@cern.ch>
Mon, 27 Oct 2014 12:42:42 +0000 (13:42 +0100)
PWGCF/Correlations/DPhi/AliAnalysisTaskDptDptQA.cxx

index 792c79e..784e453 100644 (file)
@@ -1,12 +1,19 @@
-#include "TChain.h"
-#include "TList.h"
-#include "TFile.h"
-#include "TTree.h"
-#include "TH1D.h"
-#include "TH2D.h"
-#include "TH3D.h"
-#include "THnSparse.h"
-#include "TCanvas.h"
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id:$ */
 
 #include <TROOT.h>
 #include <TChain.h>
@@ -26,6 +33,7 @@
 
 #include "AliAODHandler.h"
 #include "AliAODInputHandler.h"
+//#include "AliAODMCParticle.h"
 #include "AliInputEventHandler.h"
 #include "AliLog.h"
 #include "AliESDEvent.h"
 #include "AliCentrality.h"
 #include "AliAnalysisTaskDptDptQA.h"
 
-#include "AliPID.h"
-#include "AliPIDResponse.h"
-
-#include "AliESDVertex.h"
-#include "AliESDEvent.h"
-#include "AliESDInputHandler.h"
-#include "AliAODEvent.h"
-#include "AliAODTrack.h"
-#include "AliAODInputHandler.h"
-#include "AliESD.h"
-#include "AliESDEvent.h"
-#include "AliAODEvent.h"
-#include "AliStack.h"
-#include "AliESDtrackCuts.h"
-#include "AliAODMCHeader.h"
-
-
-#include "AliGenHijingEventHeader.h"
-#include "AliGenEventHeader.h"
-#include "AliPID.h"
-#include "AliAODPid.h"
-#include "AliPIDResponse.h"
-#include "AliAODpidUtil.h"
-#include "AliPIDCombined.h"
-
-
-
 ClassImp(AliAnalysisTaskDptDptQA)
 
 AliAnalysisTaskDptDptQA::AliAnalysisTaskDptDptQA()
@@ -68,7 +49,6 @@ AliAnalysisTaskDptDptQA::AliAnalysisTaskDptDptQA()
 fAODEvent(0), 
 fESDEvent(0),             //! ESD Event 
 fInputHandler(0),
-fPIDResponse(0x0),
 _outputHistoList(0),
 _twoPi         ( 2.0 * 3.1415927),
 _eventCount    ( 0), 
@@ -89,13 +69,12 @@ _requestedCharge_1    (   1),
 _requestedCharge_2    (  -1),
 _dcaZMin              ( -3), 
 _dcaZMax              (  3.), 
-_dcaXYMin             ( -2.4), 
-_dcaXYMax             (  2.4),
+_dcaXYMin             ( -3.), 
+_dcaXYMax             (  3.),
 _dedxMin              ( 0),
 _dedxMax              ( 100000),
-_nClusterMin          ( 80), 
-_trackFilterBit       (0),
-fNSigmaCut            (3.),
+_nClusterMin          ( 70), 
+_trackFilterBit       ( 128),
 _field    ( 1.),
 _nTracks  ( 0 ),
 _mult0    ( 0 ),
@@ -103,7 +82,6 @@ _mult1    ( 0 ),
 _mult2    ( 0 ),
 _mult3    ( 0 ),
 _mult4    ( 0 ),
-_mult4a    ( 0 ),
 _mult5    ( 0 ),
 _mult6    ( 0 ),
 arraySize ( 2000),
@@ -140,7 +118,7 @@ _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ
 
 _nBins_pt_1(18),      _min_pt_1(0.2),    _max_pt_1(2.0),          _width_pt_1(0.1),
 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
-_nBins_eta_1(0),      _min_eta_1(0),  _max_eta_1(0),           _width_eta_1(0.1),
+_nBins_eta_1(0),      _min_eta_1(0),    _max_eta_1(0),            _width_eta_1(0.1),
 
 _nBins_etaPhi_1(0), 
 _nBins_etaPhiPt_1(0),
@@ -148,7 +126,7 @@ _nBins_zEtaPhiPt_1(0),
 
 _nBins_pt_2(18),     _min_pt_2(0.2),     _max_pt_2(2.0),          _width_pt_2(0.1),
 _nBins_phi_2(72),    _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
-_nBins_eta_2(0),     _min_eta_2(0),     _max_eta_2(0),           _width_eta_2(0.1),
+_nBins_eta_2(0),    _min_eta_2(0),     _max_eta_2(0),           _width_eta_2(0.1),
 
 _nBins_etaPhi_2(0), 
 _nBins_etaPhiPt_2(0),
@@ -193,22 +171,18 @@ _m3 ( 0),
 _m4 ( 0),
 _m5 ( 0),
 _m6 ( 0),
-_vertexZ ( 0),
-  
-_Ncluster1  ( 0),
-_Ncluster2  ( 0),
+_vertexZ           ( 0),
 _etadis ( 0),
 _phidis ( 0),
 _dcaz   ( 0),
 _dcaxy  ( 0),
-_spectra  ( 0),
 _n1_1_vsPt         ( 0),         
 _n1_1_vsEtaVsPhi   ( 0),
 _s1pt_1_vsEtaVsPhi ( 0), 
 _n1_1_vsZVsEtaVsPhiVsPt ( 0),
-_n1_1_vsM          ( 0), 
+_n1_1_vsM          ( 0),  // w/ weight
 _s1pt_1_vsM        ( 0),
-_n1Nw_1_vsM        ( 0),
+_n1Nw_1_vsM        ( 0), // w/o weight
 _s1ptNw_1_vsM      ( 0),
 _dedxVsP_1         ( 0),
 _corrDedxVsP_1     ( 0),
@@ -238,6 +212,7 @@ _s2PtPtNw_12_vsM   ( 0),
 _s2PtNNw_12_vsM    ( 0),       
 _s2NPtNw_12_vsM    ( 0), 
 _invMass           ( 0),
+_invMassElec       ( 0),
 n1Name("NA"),
 n1NwName("NA"),
 n2Name("NA"),
@@ -351,7 +326,6 @@ AliAnalysisTaskDptDptQA::AliAnalysisTaskDptDptQA(const TString & name)
 fAODEvent(0), 
 fESDEvent(0),  
 fInputHandler(0),
-fPIDResponse(0),
 _outputHistoList(0),
 _twoPi         ( 2.0 * 3.1415927),
 _eventCount    ( 0), 
@@ -372,13 +346,12 @@ _requestedCharge_1    (   1),
 _requestedCharge_2    (  -1),
 _dcaZMin              ( -3), 
 _dcaZMax              (  3.), 
-_dcaXYMin             ( -2.4), 
-_dcaXYMax             (  2.4),
+_dcaXYMin             ( -3.), 
+_dcaXYMax             (  3.),
 _dedxMin              ( 0),
 _dedxMax              ( 100000),
-_nClusterMin          ( 80), 
+_nClusterMin          ( 70), 
 _trackFilterBit       ( 0),
-fNSigmaCut            ( 3.),
 _field    ( 1.),
 _nTracks  ( 0 ),
 _mult0    ( 0 ),
@@ -386,7 +359,6 @@ _mult1    ( 0 ),
 _mult2    ( 0 ),
 _mult3    ( 0 ),
 _mult4    ( 0 ),
-_mult4a    ( 0 ),
 _mult5    ( 0 ),
 _mult6    ( 0 ),
 arraySize ( 2000),
@@ -423,7 +395,7 @@ _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ
 
 _nBins_pt_1(18),      _min_pt_1(0.2),    _max_pt_1(2.0),          _width_pt_1(0.1),
 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
-_nBins_eta_1(0),      _min_eta_1(0),    _max_eta_1(0),           _width_eta_1(0.1),
+_nBins_eta_1(0),     _min_eta_1(0),      _max_eta_1(0),           _width_eta_1(0.1),
 
 _nBins_etaPhi_1(0), 
 _nBins_etaPhiPt_1(0),
@@ -431,7 +403,7 @@ _nBins_zEtaPhiPt_1(0),
 
 _nBins_pt_2(18),     _min_pt_2(0.2),     _max_pt_2(2.0),          _width_pt_2(0.1),
 _nBins_phi_2(72),    _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
-_nBins_eta_2(0),    _min_eta_2(0),     _max_eta_2(0),           _width_eta_2(0.1),
+_nBins_eta_2(0),     _min_eta_2(0),      _max_eta_2(0),           _width_eta_2(0.1),
 
 _nBins_etaPhi_2(0), 
 _nBins_etaPhiPt_2(0),
@@ -476,22 +448,18 @@ _m3 ( 0),
 _m4 ( 0),
 _m5 ( 0),
 _m6 ( 0),
-_vertexZ ( 0),
-_Ncluster1  ( 0),
-_Ncluster2  ( 0),
+_vertexZ           ( 0),
 _etadis ( 0),
 _phidis ( 0),
-
-_dcaz ( 0),
-_dcaxy ( 0),
-_spectra  ( 0),
+_dcaz   ( 0),
+_dcaxy  ( 0),
 _n1_1_vsPt         ( 0),         
 _n1_1_vsEtaVsPhi   ( 0),
 _s1pt_1_vsEtaVsPhi ( 0), 
 _n1_1_vsZVsEtaVsPhiVsPt ( 0),
-_n1_1_vsM          ( 0), 
+_n1_1_vsM          ( 0),  // w/ weight
 _s1pt_1_vsM        ( 0),
-_n1Nw_1_vsM        ( 0), 
+_n1Nw_1_vsM        ( 0), // w/o weight
 _s1ptNw_1_vsM      ( 0),
 _dedxVsP_1         ( 0),
 _corrDedxVsP_1     ( 0),
@@ -521,6 +489,7 @@ _s2PtPtNw_12_vsM   ( 0),
 _s2PtNNw_12_vsM    ( 0),       
 _s2NPtNw_12_vsM    ( 0), 
 _invMass           ( 0),
+_invMassElec       ( 0),
 n1Name("NA"),
 n1NwName("NA"),
 n2Name("NA"),
@@ -633,15 +602,17 @@ vsPtVsPt("NA")
 
 AliAnalysisTaskDptDptQA::~AliAnalysisTaskDptDptQA()
 {
-  
+
 }
 
 void AliAnalysisTaskDptDptQA::UserCreateOutputObjects()
 {
+  cout<< "AliAnalysisTaskDptDptQA::CreateOutputObjects() Starting " << endl;
   OpenFile(0);
   _outputHistoList = new TList();
   _outputHistoList->SetOwner();
   
+  //Setup the parameters of the histograms
   _nBins_M0 = 500; _min_M0   = 0.;    _max_M0    = 5000.;  _width_M0 = (_max_M0-_min_M0)/_nBins_M0;
   _nBins_M1 = 500; _min_M1   = 0.;    _max_M1    = 5000.;  _width_M1 = (_max_M1-_min_M1)/_nBins_M1;
   _nBins_M2 = 500; _min_M2   = 0.;    _max_M2    = 5000.;  _width_M2 = (_max_M2-_min_M2)/_nBins_M2;
@@ -668,37 +639,54 @@ void AliAnalysisTaskDptDptQA::UserCreateOutputObjects()
   _nBins_etaPhiPt_2  = _nBins_etaPhi_2 * _nBins_pt_2;
   _nBins_zEtaPhiPt_2 = _nBins_vertexZ  * _nBins_etaPhiPt_2;
   _nBins_etaPhi_12   = _nBins_etaPhi_1 * _nBins_etaPhi_2;
-    
+  
+  //setup the work arrays
+  
   _id_1       = new int[arraySize];   
   _charge_1   = new int[arraySize]; 
+  //_iPhi_1 = new int[arraySize]; 
+  //_iEta_1 = new int[arraySize]; 
   _iEtaPhi_1  = new int[arraySize]; 
   _iPt_1      = new int[arraySize];  
   _pt_1       = new float[arraySize];    
   _px_1       = new float[arraySize];   
   _py_1       = new float[arraySize];   
   _pz_1       = new float[arraySize];   
+  //_phi_1 = new float[arraySize];  
+  //_eta_1 = new float[arraySize];  
   _correction_1 = new float[arraySize];
-    
+  _dedx_1     = new float[arraySize];
+  
   __n1_1_vsPt              = getDoubleArray(_nBins_pt_1,        0.);
   __n1_1_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_1,    0.);
   __s1pt_1_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_1,    0.);
   __n1_1_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_1,  0.);
   
-    
+  cout << "==========================================================================================" << endl;
+  cout << "=============== Booking for particle 1 done." << endl;
+  cout << "_requestedCharge_1: " << _requestedCharge_1 << endl;
+  cout << "_requestedCharge_2: " << _requestedCharge_2 << endl;
+  
   if (_requestedCharge_2!=_requestedCharge_1)
     {
-      _sameFilter = 0;
+    cout << " creating arrays for particle 2 with size: " << arraySize << endl;
+    _sameFilter = 0;
     //particle 2
     _id_2       = new int[arraySize];   
     _charge_2   = new int[arraySize]; 
+    //_iPhi_2   = new int[arraySize]; 
+    //_iEta_2   = new int[arraySize]; 
     _iEtaPhi_2  = new int[arraySize]; 
     _iPt_2      = new int[arraySize];  
     _pt_2       = new float[arraySize];   
     _px_2       = new float[arraySize];   
     _py_2       = new float[arraySize];   
     _pz_2       = new float[arraySize];   
+    //_phi_2    = new float[arraySize];  
+    //_eta_2    = new float[arraySize];  
     _correction_2 = new float[arraySize];
-        
+    _dedx_2       = new float[arraySize];
+    
     __n1_2_vsPt              = getDoubleArray(_nBins_pt_2,        0.);
     __n1_2_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_2,    0.);
     __s1pt_2_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_2,    0.);
@@ -808,60 +796,63 @@ void AliAnalysisTaskDptDptQA::UserCreateOutputObjects()
   
   if (_useWeights)
     {
-    int iZ, iEtaPhi, iPt;
-    int iZ1,iEtaPhi1,iPt1;
-    int a, b;
-    if (_weight_1)
-      {
-      _correctionWeight_1 = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
-      a = _nBins_pt_1;
-      b = _nBins_etaPhi_1*_nBins_pt_1;
-      for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
-        {
-        for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_1; iEtaPhi++, iEtaPhi1++)
-          {
-          for (iPt=0,iPt1=1; iPt<_nBins_pt_1; iPt++, iPt1++)
-            {
-            _correctionWeight_1[iZ*b+iEtaPhi*a+iPt] = _weight_1->GetBinContent(iZ1,iEtaPhi1,iPt1);
-            }      
-          }
-        }
-      } // _weight_1
-    else
-      {
-      AliError("AliAnalysisTaskDptDptQA:: _weight_1 is a null pointer.");
-      return;
-      }
-    if (!_sameFilter) 
-      {
-      if (_weight_2)
-        {
-        _correctionWeight_2 = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
-        a = _nBins_pt_2;
-        b = _nBins_etaPhi_2*_nBins_pt_2;
-        for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
-          {
-          for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_2; iEtaPhi++, iEtaPhi1++)
-            {
-            for (iPt=0,iPt1=1; iPt<_nBins_pt_2; iPt++, iPt1++)
-              {
-              _correctionWeight_2[iZ*b+iEtaPhi*a+iPt] = _weight_2->GetBinContent(iZ1,iEtaPhi1,iPt1);
-              }      
-            }
-          }
-        } // _weight_2
+      int iZ, iEtaPhi, iPt;
+      int iZ1,iEtaPhi1,iPt1;
+      int a, b;
+      if (_weight_1)
+       {
+
+         cout<<"Prabhat Sanu weeding is here *************"<<"   "<<endl;
+
+         _correctionWeight_1 = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
+         a = _nBins_pt_1;
+         b = _nBins_etaPhi_1*_nBins_pt_1;
+         for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
+           {
+             for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_1; iEtaPhi++, iEtaPhi1++)
+               {
+                 for (iPt=0,iPt1=1; iPt<_nBins_pt_1; iPt++, iPt1++)
+                   {
+                     _correctionWeight_1[iZ*b+iEtaPhi*a+iPt] = _weight_1->GetBinContent(iZ1,iEtaPhi1,iPt1);
+                   }      
+               }
+           }
+       } // _weight_1
       else
-        {
-        AliError("AliAnalysisTaskDptDptQA:: _weight_1 is a null pointer.");
-        return;
-        }
-      }
+       {
+         AliError("AliAnalysisTaskDptDptQA:: _weight_1 is a null pointer.");
+         return;
+       }
+      if (!_sameFilter) 
+       {
+         if (_weight_2)
+           {
+             _correctionWeight_2 = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
+             a = _nBins_pt_2;
+             b = _nBins_etaPhi_2*_nBins_pt_2;
+             for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
+               {
+                 for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_2; iEtaPhi++, iEtaPhi1++)
+                   {
+                     for (iPt=0,iPt1=1; iPt<_nBins_pt_2; iPt++, iPt1++)
+                       {
+                         _correctionWeight_2[iZ*b+iEtaPhi*a+iPt] = _weight_2->GetBinContent(iZ1,iEtaPhi1,iPt1);
+                       }      
+                   }
+               }
+           } // _weight_2
+         else
+           {
+             AliError("AliAnalysisTaskDptDptQA:: _weight_1 is a null pointer.");
+             return;
+           }
+       }
     }
   
   createHistograms();
   PostData(0,_outputHistoList);
   
-  //cout<< "AliAnalysisTaskDptDptQA::CreateOutputObjects() DONE " << endl;
+  cout<< "AliAnalysisTaskDptDptQA::CreateOutputObjects() DONE " << endl;
   
 }
 
@@ -872,7 +863,13 @@ void  AliAnalysisTaskDptDptQA::createHistograms()
   
   name = "eventAccounting";
   
-   _eventAccounting      = createHisto1D(name,name,10, -0.5, 9.5, "event Code", _title_counts);
+  // bin index : what it is...
+  //         0 :  number of event submitted
+  //         1 :  number accepted by centrality cut
+  //         2 :  number accepted by centrality cut and z cut
+  //         3 :  total number of particles that satisfied filter 1
+  //         4 :  total number of particles that satisfied filter 2
+  _eventAccounting      = createHisto1D(name,name,10, -0.5, 9.5, "event Code", _title_counts);
   
   name = "m0"; _m0      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m0, _title_counts);
   name = "m1"; _m1      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m1, _title_counts);
@@ -881,44 +878,42 @@ void  AliAnalysisTaskDptDptQA::createHistograms()
   name = "m4"; _m4      = createHisto1D(name,name,_nBins_M4, _min_M4, _max_M4, _title_m4, _title_counts);
   name = "m5"; _m5      = createHisto1D(name,name,_nBins_M5, _min_M5, _max_M5, _title_m5, _title_counts);
   name = "m6"; _m6      = createHisto1D(name,name,_nBins_M6, _min_M6, _max_M6, _title_m6, _title_counts);
-  name = "zV"; _vertexZ = createHisto1D(name,name,100, -10, 10, "z-Vertex (cm)", _title_counts);
+  name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
   
-
   name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
   name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
-  name = "DCAz";    _dcaz     = createHisto1F(name,name, 500, -5.0, 5.0, "dcaZ","counts");
-  name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 500, -5.0, 5.0, "dcaXY","counts");
-  name = "pTspectra";   _spectra    = createHisto1F(name,name, 300, 0.1, 2.5, "pT","counts");
-
-  //name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
-  //name = "Nclus2";   _Ncluster2    = createHisto1F(name,name, 200, 0, 200, "Ncluster2","counts");
   
   if (_singlesOnly)
     {
-    name = n1Name+part_1_Name+vsPt;              _n1_1_vsPt              = createHisto1F(name,name, _nBins_pt_1,  _min_pt_1,  _max_pt_1,   _title_pt_1,  _title_AvgN_1);
-    name = n1Name+part_1_Name+vsZ+vsEtaPhi+vsPt; _n1_1_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_1, 0., double(_nBins_etaPhi_1), _nBins_pt_1, _min_pt_1, _max_pt_1, "zVertex", _title_etaPhi_1,  _title_pt_1);
-
-    name = n1Name+part_2_Name+vsPt;              _n1_2_vsPt              = createHisto1F(name,name, _nBins_pt_2,  _min_pt_2,  _max_pt_2,   _title_pt_2,  _title_AvgN_2);
-    name = n1Name+part_2_Name+vsZ+vsEtaPhi+vsPt; _n1_2_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_2, 0., double(_nBins_etaPhi_2), _nBins_pt_2, _min_pt_2, _max_pt_2, "zVertex", _title_etaPhi_2,  _title_pt_2);
-
+      name = n1Name+part_1_Name+vsPt;              _n1_1_vsPt              = createHisto1F(name,name, _nBins_pt_1,  _min_pt_1,  _max_pt_1,   _title_pt_1,  _title_AvgN_1);
+      name = n1Name+part_1_Name+vsZ+vsEtaPhi+vsPt; _n1_1_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_1, 0., double(_nBins_etaPhi_1), _nBins_pt_1, _min_pt_1, _max_pt_1, "zVertex", _title_etaPhi_1,  _title_pt_1);
+      //name = "dedxVsP_1";                          _dedxVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
+      //name = "corrDedxVsP_1";                      _corrDedxVsP_1          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
+      //name = "betaVsP_1";                          _betaVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
+      
+      name = n1Name+part_2_Name+vsPt;              _n1_2_vsPt              = createHisto1F(name,name, _nBins_pt_2,  _min_pt_2,  _max_pt_2,   _title_pt_2,  _title_AvgN_2);
+      name = n1Name+part_2_Name+vsZ+vsEtaPhi+vsPt; _n1_2_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_2, 0., double(_nBins_etaPhi_2), _nBins_pt_2, _min_pt_2, _max_pt_2, "zVertex", _title_etaPhi_2,  _title_pt_2);
+      //name = "dedxVsP_2";                          _dedxVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
+      //name = "corrDedxVsP_2";                      _corrDedxVsP_2          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
+      //name = "betaVsP_2";                          _betaVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
     }
   else
     {
-    name = n1Name+part_1_Name+vsEtaPhi;       _n1_1_vsEtaVsPhi      = createHisto2F(name,name, _nBins_eta_1, _min_eta_1, _max_eta_1,  _nBins_phi_1, _min_phi_1, _max_phi_1,  _title_eta_1,  _title_phi_1,  _title_AvgN_1);
-    name = s1ptName+part_1_Name+vsEtaPhi;     _s1pt_1_vsEtaVsPhi    = createHisto2F(name,name, _nBins_eta_1, _min_eta_1, _max_eta_1,  _nBins_phi_1, _min_phi_1, _max_phi_1,  _title_eta_1,  _title_phi_1,  _title_AvgSumPt_1);
-    name = n1Name+part_1_Name+vsM;            _n1_1_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
-    name = s1ptName+part_1_Name+vsM;          _s1pt_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
-    name = n1NwName+part_1_Name+vsM;          _n1Nw_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
-    name = s1ptNwName+part_1_Name+vsM;        _s1ptNw_1_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
-
-    name = n1Name+part_2_Name+vsEtaPhi;       _n1_2_vsEtaVsPhi      = createHisto2F(name,name, _nBins_eta_2, _min_eta_2, _max_eta_2,  _nBins_phi_2, _min_phi_2, _max_phi_2,  _title_eta_2,  _title_phi_2,  _title_AvgN_2);
-    name = s1ptName+part_2_Name+vsEtaPhi;     _s1pt_2_vsEtaVsPhi    = createHisto2F(name,name, _nBins_eta_2, _min_eta_2, _max_eta_2,  _nBins_phi_2, _min_phi_2, _max_phi_2,  _title_eta_2,  _title_phi_2,  _title_AvgSumPt_2);
-    name = n1Name+part_2_Name + vsM;          _n1_2_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_2);
-    name = s1ptName+part_2_Name + vsM;        _s1pt_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_2);
-    name = n1NwName+part_2_Name+vsM;          _n1Nw_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
-    name = s1ptNwName+part_2_Name+vsM;        _s1ptNw_2_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
-
-    name = n2Name+pair_12_Name+vsEtaPhi;      _n2_12_vsEtaPhi       = createHisto1F(name,name, _nBins_etaPhi_12, 0.,        double(_nBins_etaPhi_12), _title_etaPhi_12, _title_AvgN2_12);        
+      name = n1Name+part_1_Name+vsEtaPhi;       _n1_1_vsEtaVsPhi      = createHisto2F(name,name, _nBins_eta_1, _min_eta_1, _max_eta_1,  _nBins_phi_1, _min_phi_1, _max_phi_1,  _title_eta_1,  _title_phi_1,  _title_AvgN_1);
+      name = s1ptName+part_1_Name+vsEtaPhi;     _s1pt_1_vsEtaVsPhi    = createHisto2F(name,name, _nBins_eta_1, _min_eta_1, _max_eta_1,  _nBins_phi_1, _min_phi_1, _max_phi_1,  _title_eta_1,  _title_phi_1,  _title_AvgSumPt_1);
+      name = n1Name+part_1_Name+vsM;            _n1_1_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
+      name = s1ptName+part_1_Name+vsM;          _s1pt_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
+      name = n1NwName+part_1_Name+vsM;          _n1Nw_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
+      name = s1ptNwName+part_1_Name+vsM;        _s1ptNw_1_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
+      
+      name = n1Name+part_2_Name+vsEtaPhi;       _n1_2_vsEtaVsPhi      = createHisto2F(name,name, _nBins_eta_2, _min_eta_2, _max_eta_2,  _nBins_phi_2, _min_phi_2, _max_phi_2,  _title_eta_2,  _title_phi_2,  _title_AvgN_2);
+      name = s1ptName+part_2_Name+vsEtaPhi;     _s1pt_2_vsEtaVsPhi    = createHisto2F(name,name, _nBins_eta_2, _min_eta_2, _max_eta_2,  _nBins_phi_2, _min_phi_2, _max_phi_2,  _title_eta_2,  _title_phi_2,  _title_AvgSumPt_2);
+      name = n1Name+part_2_Name + vsM;          _n1_2_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_2);
+      name = s1ptName+part_2_Name + vsM;        _s1pt_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_2);
+      name = n1NwName+part_2_Name+vsM;          _n1Nw_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
+      name = s1ptNwName+part_2_Name+vsM;        _s1ptNw_2_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
+      
+      name = n2Name+pair_12_Name+vsEtaPhi;      _n2_12_vsEtaPhi       = createHisto1F(name,name, _nBins_etaPhi_12, 0.,        double(_nBins_etaPhi_12), _title_etaPhi_12, _title_AvgN2_12);        
     name = s2PtPtName+pair_12_Name + vsEtaPhi;_s2PtPt_12_vsEtaPhi   = createHisto1F(name,name, _nBins_etaPhi_12, 0.,        double(_nBins_etaPhi_12), _title_etaPhi_12,  _title_AvgSumPtPt_12);    
     name = s2PtNName+pair_12_Name + vsEtaPhi; _s2PtN_12_vsEtaPhi    = createHisto1F(name,name, _nBins_etaPhi_12, 0.,        double(_nBins_etaPhi_12), _title_etaPhi_12,  _title_AvgSumPtN_12);    
     name = s2NPtName+pair_12_Name + vsEtaPhi; _s2NPt_12_vsEtaPhi    = createHisto1F(name,name, _nBins_etaPhi_12, 0.,        double(_nBins_etaPhi_12), _title_etaPhi_12,  _title_AvgNSumPt_12);    
@@ -934,8 +929,8 @@ void  AliAnalysisTaskDptDptQA::createHistograms()
     name = s2PtNNwName+pair_12_Name + vsM;    _s2PtNNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);       
     name = s2NPtNwName+pair_12_Name + vsM;    _s2NPtNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);     
     
-    name = "mInv";     _invMass     = createHisto2F(name,name, 10, 0.0, 2.0, 200, 0.4, 0.6, "pT","InvMass","ptvsmass"); 
-    
+    name = "mInv";     _invMass     = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts"); 
+    name = "mInvElec"; _invMassElec = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts"); 
     }
   
   AliInfo(" AliAnalysisTaskDptDptQA::createHistoHistograms() All Done"); 
@@ -947,7 +942,46 @@ void  AliAnalysisTaskDptDptQA::finalizeHistograms()
   
   AliInfo("AliAnalysisTaskDptDptQA::finalizeHistograms() starting");
   AliInfo(Form("CorrelationAnalyzers::finalizeHistograms()   _eventCount : %d",int(_eventCount)));
-  
+  if (_singlesOnly)
+    {
+    if (_sameFilter)
+      {
+      fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
+      fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
+      fillHistoWithArray(_n1_2_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
+      fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
+      }
+    else
+      {
+      fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
+      fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
+      fillHistoWithArray(_n1_2_vsPt,              __n1_2_vsPt,        _nBins_pt_2);
+      fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_2_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_2, _nBins_pt_2);
+      }
+    }
+  else
+    {
+    if (_sameFilter)
+      {
+      fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
+      fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
+      fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
+      fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
+      }
+    else
+      {
+      fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
+      fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
+      fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_2_vsEtaPhi,    _nBins_eta_2,   _nBins_phi_2);
+      fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_2_vsEtaPhi,  _nBins_eta_2,   _nBins_phi_2);
+      }
+    fillHistoWithArray(_n2_12_vsEtaPhi,     __n2_12_vsEtaPhi,     _nBins_etaPhi_12);
+    fillHistoWithArray(_s2PtPt_12_vsEtaPhi, __s2ptpt_12_vsEtaPhi, _nBins_etaPhi_12);
+    fillHistoWithArray(_s2PtN_12_vsEtaPhi,  __s2PtN_12_vsEtaPhi,  _nBins_etaPhi_12);
+    fillHistoWithArray(_s2NPt_12_vsEtaPhi,  __s2NPt_12_vsEtaPhi,  _nBins_etaPhi_12);
+    fillHistoWithArray(_n2_12_vsPtVsPt,     __n2_12_vsPtPt,       _nBins_pt_1,    _nBins_pt_2);
+
+    }  
   AliInfo("AliAnalysisTaskDptDptQA::finalizeHistograms()  Done ");
 }
 //--------------//
@@ -957,55 +991,47 @@ void  AliAnalysisTaskDptDptQA::UserExec(Option_t */*option*/)
 {
   
   int    k1,k2;
-  int  charge;
-  float  q, phi, pt, eta,  px, py, pz;
-  //int    ij;
-  //int    id_1, q_1, iPt_1;
-  //float  pt_1, px_1, py_1, pz_1, corr_1;
-  //int    id_2, q_2,  iPt_2;
-  //float  pt_2, px_2, py_2, pz_2, corr_2;
-  //float  ptpt;
+  int    iPhi, iEta, iEtaPhi, iPt, charge;
+  float  q, p, phi, pt, eta, corr, corrPt, dedx, px, py, pz;
+  int    ij;
+  int    id_1, q_1, iEtaPhi_1, iPt_1;
+  float  pt_1, px_1, py_1, pz_1, corr_1, dedx_1;
+  int    id_2, q_2, iEtaPhi_2, iPt_2;
+  float  pt_2, px_2, py_2, pz_2, corr_2, dedx_2;
+  float  ptpt;
   int    iVertex, iVertexP1, iVertexP2;
-  //int    iZEtaPhiPt;
-  //float  massElecSq = 1.94797849000000016e-02;
-  //double b[2];
-  //double bCov[3];
+  int    iZEtaPhiPt;
+  float  massElecSq = 2.5e-7;
   const  AliAODVertex* vertex;
+  int    nClus;
   bool   bitOK;
-  
-  AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
-  if (!manager) {
-    return;
-  }
-  AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
-  if (!inputHandler) {
-    return;
-  }
-  
-  fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
-  //AliAODEvent* fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
-  
-  if (!fAODEvent)
-    {
-      return;
-    }
-  fPIDResponse =inputHandler->GetPIDResponse();
-  if (!fPIDResponse){
-    AliFatal("This Task needs the PID response attached to the inputHandler");
-    return;
-  }
-  
+    
   // count all events looked at here
   _eventCount++;
-  
+
+  // A. Adare - Fix from above to avoid coverity complaint about null dereference 
   if (_eventAccounting)
     {
       _eventAccounting->Fill(0);// count all calls to this function
     }
   else 
     {
-      
-      return;
+    cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - !_eventAccounting" << endl;
+    cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - Skip this task" << endl;
+    return;
+    }
+    
+  //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 1" << endl;
+  
+  //Get pointer to current event
+  //fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
+  fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
+  
+  //  if(!fAODEvent && !fESDEvent)
+  if(!fAODEvent)
+    {
+    AliError("AliAnalysisTaskDptDptQA::Exec(Option_t *option) !fAODEvent");
+    return;
     }
   
   _eventAccounting->Fill(1);// count all calls to this function with a valid pointer
@@ -1014,244 +1040,560 @@ void  AliAnalysisTaskDptDptQA::UserExec(Option_t */*option*/)
   __n1_1 = __n1_2 = __s1pt_1 = __s1pt_2 = __n1Nw_1 = __n1Nw_2 = __s1ptNw_1 = __s1ptNw_2 = 0;
   
   float v0Centr  = -999.;
-  float v0ACentr  = -999.;
   float trkCentr = -999.;
   float spdCentr = -999.;
-  
   float vertexX  = -999;
   float vertexY  = -999;
   float vertexZ  = -999;
-  //float vertexXY = -999;
-  // float dcaZ     = -999;
-  //float dcaXY    = -999;
+  float vertexXY = -999;
   float centrality = -999;
   
-  float Kaon1Charge[50000];
-  float Kaon1Px[50000];
-  float Kaon1Py[50000];
-  float Kaon1Pz[50000];
-  float Kaon2Charge[50000];
-  float Kaon2Px[50000];
-  float Kaon2Py[50000];
-  float Kaon2Pz[50000];
-
-  int nKaon1s = 0;
-  int nKaon2s = 0;
-
-  //k1Mass = 0.139570;
-  //kMass = 0.139570;
-  
-
   if(fAODEvent)
     {
-      //Centrality
-      AliCentrality* centralityObject =  fAODEvent->GetHeader()->GetCentralityP();
-      if (centralityObject)
-       {
-         //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 6" << endl;
-         
-         v0Centr  = centralityObject->GetCentralityPercentile("V0M");
-         v0ACentr  = centralityObject->GetCentralityPercentile("V0A");
-         trkCentr = centralityObject->GetCentralityPercentile("TRK"); 
-         spdCentr = centralityObject->GetCentralityPercentile("CL1");
-         
-       }
+    //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 5" << endl;
+    
+    //Centrality
+    AliCentrality* centralityObject =  fAODEvent->GetHeader()->GetCentralityP();
+    if (centralityObject)
+      {
+      //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 6" << endl;
       
-      _nTracks  =fAODEvent->GetNumberOfTracks();//NEW Test
-      _mult3    = _nTracks; 
-      _mult4    = v0Centr;
-      _mult4a    = v0ACentr;
-      _mult5    = trkCentr;
-      _mult6    = spdCentr;
-      _field    = fAODEvent->GetMagneticField(); 
+      v0Centr  = centralityObject->GetCentralityPercentile("V0M");
+      trkCentr = centralityObject->GetCentralityPercentile("TRK"); 
+      spdCentr = centralityObject->GetCentralityPercentile("CL1");
+      //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 7" << endl;
       
-      //_centralityMethod
-      switch (_centralityMethod)
-       {
-       case 0: centrality = _mult0; break;
-       case 1: centrality = _mult1; break;
-       case 2: centrality = _mult2; break;
-       case 3: centrality = _mult3; break;
-       case 4: centrality = _mult4; break;
-       case 5: centrality = _mult5; break;
-       case 6: centrality = _mult6; break;
-       case 7: centrality = _mult4a; break;
-       }
+      }
+    //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 8" << endl;
+    
+    _nTracks  = fAODEvent->GetNumberOfTracks();
+    //_mult0    = 0;
+    //_mult1    = 0;
+    //_mult2    = 0;
+    _mult3    = _nTracks; 
+    _mult4    = v0Centr;
+    _mult5    = trkCentr;
+    _mult6    = spdCentr;
+    _field    = fAODEvent->GetMagneticField(); 
+    
+    //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - _field:" << _field << endl;
+    
+    //_centralityMethod
+    switch (_centralityMethod)
+      {
+        case 0: centrality = _mult0; break;
+        case 1: centrality = _mult1; break;
+        case 2: centrality = _mult2; break;
+        case 3: centrality = _mult3; break;
+        case 4: centrality = _mult4; break;
+        case 5: centrality = _mult5; break;
+        case 6: centrality = _mult6; break;
+      }
+    
+    //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 10" << endl;
+    
+    //filter on centrality
+    // require the v0 and trk centralities to agree within 5%.
+    
+    //if ((Float_t(GetTPCMult(fAOD)) > (-40.3+1.22*GetGlobalMult(fAOD))) && (Float_t(GetTPCMult(fAOD)) < (32.1+1.59*GetGlobalMult(fAOD)))) {//3 sigma
+
+    
+    if ( centrality < _centralityMin ||  
+         centrality > _centralityMax ||
+         fabs(v0Centr-trkCentr)>5.0)
+      {
+      //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 11" << endl;
       
+      // we dont analyze this event here... 
+      return;
+      }
+    _eventAccounting->Fill(2);// count all events with right centrality
+    //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 12" << endl;
+    
+    // filter on z and xy vertex
+    vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertexSPD();
+    if (!vertex || vertex->GetNContributors()<1)
+      {
+      vertexZ  = -999;
+      vertexXY = -999;
+      cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - No valid vertex object or poor vertex" << endl;
+      }
+    else
+      { 
+        vertexX = vertex->GetX();
+        vertexY = vertex->GetY();
+        vertexZ = vertex->GetZ();
+        vertexXY = sqrt(vertexX*vertexX+vertexY*vertexY);
+      }
+    if (!vertex ||
+        vertexZ  < _vertexZMin  || 
+        vertexZ  > _vertexZMax  ||
+        vertexXY < _vertexXYMin || 
+        vertexXY > _vertexXYMax)
+      return;
+    
+    iVertex = int((vertexZ-_min_vertexZ)/_width_vertexZ);
+    iVertexP1 = iVertex*_nBins_etaPhiPt_1;
+    iVertexP2 = iVertex*_nBins_etaPhiPt_2;
+    if (iVertex<0 || iVertex>=_nBins_vertexZ)
+      {
+      AliError("AliAnalysisTaskDptDptQA::Exec(Option_t *option) iVertex<0 || iVertex>=_nBins_vertexZ ");
+      return;
+      }
+    _eventAccounting->Fill(3);// count all calls to this function with a valid pointer
+    
+    // loop over all particles for singles
+    //debug() << "_nTracks:"<< _nTracks << endl;
+    for (int iTrack=0; iTrack< _nTracks; iTrack++)
+      {
+      //cout << "AliAnalysisTaskDptDptQA::UserExec(Option_t *option) - 16 iTrack:" << iTrack << endl;
       
-      if ( centrality < _centralityMin ||  centrality > _centralityMax )
-       {
-         return;
-       }
-      _eventAccounting->Fill(2);// count all events with right centrality
+      AliAODTrack * t = (AliAODTrack *) fAODEvent->GetTrack(iTrack);
+      if (!t) 
+        {
+        AliError(Form("AliAnalysisTaskDptDptQA::Exec(Option_t *option) No track ofr iTrack=%d", iTrack));
+        continue;
+        }
       
-      // filter on z and xy vertex
-      vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertex();
-      // Double_t V[2];
-      //vertex->GetXYZ(V);      
-
-      if(vertex)
-       {
-         Double32_t fCov[6];
-         vertex->GetCovarianceMatrix(fCov);
-         if(vertex->GetNContributors() > 0)
-           {
-             if(fCov[5] != 0)
-               {
-                 vertexX = vertex->GetX();
-                 vertexY = vertex->GetY();
-                 vertexZ = vertex->GetZ();
-                 
-                 if(TMath::Abs(vertexZ) > 10)
-                   {
-                     return;
-                   } // Z-Vertex Cut  
-               }
-           }
-       }
+      q      = t->Charge();
+      charge = int(q);
+      phi    = t->Phi();
+      p      = t->P();
+      pt     = t->Pt();
+      px     = t->Px();
+      py     = t->Py();
+      pz     = t->Pz();
+      eta    = t->Eta();
+      dedx   = t->GetTPCsignal();
+      nClus  = t->GetTPCNcls();
+      bitOK  = t->TestFilterBit(_trackFilterBit);
       
-      _vertexZ->Fill(vertexZ);
+      //cout << "_trackFilterBit:" << _trackFilterBit << " Track returns:" << bitOK << endl;
+      //cout << "    q:" << q    << " _requestedCharge_1:" << _requestedCharge_1 << endl;
+      //cout << "   pt:" << pt   << " _min_pt_1:" << _min_pt_1 << " _max_pt_1:" << _max_pt_1<< endl;
+      //cout << "  phi:" << phi  << endl;
+      //cout << "  eta:" << eta  << " _min_eta_1:" << _min_eta_1 << " _max_eta_1:" << _max_eta_1<< endl;
+      //cout << " dedx:" << dedx << " _dedxMin:" << _dedxMin << " _dedxMax:" << _dedxMax << endl;
+      //cout << "nclus:" << nClus<< " _nClusterMin:" << _nClusterMin << endl;
       
-      iVertex = int((vertexZ-_min_vertexZ)/_width_vertexZ);
-      iVertexP1 = iVertex*_nBins_etaPhiPt_1;
-      iVertexP2 = iVertex*_nBins_etaPhiPt_2;
-      if (iVertex<0 || iVertex>=_nBins_vertexZ)
-       {
-         AliError("AliAnalysisTaskDptDptQA::Exec(Option_t *option) iVertex<0 || iVertex>=_nBins_vertexZ ");
-         return;
-       }
-      _eventAccounting->Fill(3);// count all calls to this function with a valid pointer
-      //====================== 
+      //common cuts first
       
-      //*********************************************************
-      TExMap *trackMap = new TExMap();//Mapping matrix----                                            
-      //1st loop track                                                                                
-      for(Int_t i = 0; i < _nTracks; i++)
-       {
-         AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
-         if(!aodTrack) {
-           AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
-           continue;
-         }
-         Int_t gID = aodTrack->GetID();
-         if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                       
-       }
-           
-      AliAODTrack* newAodTrack;
+      if (!bitOK || dedx<_dedxMin || dedx>_dedxMax || nClus<_nClusterMin) continue;
+      // Kinematics cuts used 
+                                                                                                         
+      if( pt < _min_pt_1 || pt > _max_pt_1) continue;
+      if( eta < _min_eta_1 || eta > _max_eta_1) continue;
+
+      _etadis->Fill(eta);                                                                          
+      _phidis->Fill(phi);   
       
-      //Track Loop starts here
-      for (int iTrack=0; iTrack< _nTracks; iTrack++)
-       {
-         AliAODTrack* t = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(iTrack));
-         if (!t) {
-           AliError(Form("Could not receive track %d", iTrack));
-           continue;
-         }
+      //Particle 1 
+      if (_requestedCharge_1 == charge &&
+         eta      >=  0.2 &&
+         eta      <   0.3)
+        {
+         iPhi   = int( phi/_width_phi_1);
          
-         bitOK  = t->TestFilterBit(_trackFilterBit);
-         if (!bitOK) continue; //128bit or 272bit
-                 
-         Int_t gID = t->GetID();
-         newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
+         if (iPhi<0 || iPhi>=_nBins_phi_1 ) 
+           {
+             AliWarning("AliAnalysisTaskDptDptQA::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
+             return;
+           }
          
-         q      = t->Charge();
-         charge = int(q);
-         phi    = t->Phi();
-         pt     = t->Pt();
-         px     = t->Px();
-         py     = t->Py();
-         pz     = t->Pz();
-         eta    = t->Eta();
-         //dcaXY = t->DCA(); 
-         //dcaZ  = t->ZAtDCA();  
-
-         TVector3 mom(px, py, pz);
-
-         //Double_t nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
-         Double_t nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kPion));
-         //Double_t nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kKaon));
-         //Double_t nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kProton));
+         iEta    = int((eta-_min_eta_1)/_width_eta_1); 
+         if (iEta<0 || iEta>=_nBins_eta_1) 
+           {
+             AliWarning(Form("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
+             continue;
+           }
+         iPt     = int((pt -_min_pt_1 )/_width_pt_1 ); 
+         if (iPt<0  || iPt >=_nBins_pt_1)
+           {
+             AliWarning(Form("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
+             continue;
+           }
+         iEtaPhi = iEta*_nBins_phi_1+iPhi;
+         iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
          
-         if(charge == 0) continue;
+         if (_correctionWeight_1)
+           corr = _correctionWeight_1[iZEtaPhiPt];
+         else
+           corr = 1;
+         if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1)
+           {
+             AliWarning("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1");
+             continue;
+           }
+         //cout << "all good; process track:" << endl;
          
-         // Kinematics cuts used                                                                                        
-         if( pt < _min_pt_1 || pt > _max_pt_1) continue;
-         if( eta < _min_eta_1 || eta > _max_eta_1) continue;
-
-         Double_t pos[3];
-         newAodTrack->GetXYZ(pos);
-
-         Double_t DCAX = pos[0] - vertexX;
-         Double_t DCAY = pos[1] - vertexY;
-         Double_t DCAZ = pos[2] - vertexZ;
-               
-         Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+         if (_singlesOnly)
+           {
+             //_betaVsP_1->Fill(p*q,_trackFilter_1->getBeta() );
+             //_dedxVsP_1->Fill(p*q,dedx);
+             //_corrDedxVsP_1->Fill(p*q,_trackFilter_1->getCorrectedDedx() );
+             __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
+             __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
+             
+           }
+         else
+           {
+             corrPt                      = corr*pt;
+             _id_1[k1]                   = iTrack;     
+             _charge_1[k1]               = charge;
+             _iEtaPhi_1[k1]              = iEtaPhi; 
+             _iPt_1[k1]                  = iPt;   
+             _pt_1[k1]                   = pt;   
+             _px_1[k1]                   = px;   
+             _py_1[k1]                   = py;   
+             _pz_1[k1]                   = pz;   
+             _correction_1[k1]           = corr; 
+             __n1_1                      += corr;
+             __n1_1_vsEtaPhi[iEtaPhi]    += corr; 
+             __s1pt_1                    += corrPt;
+             __s1pt_1_vsEtaPhi[iEtaPhi]  += corrPt;
+             __n1Nw_1                    += 1;
+             __s1ptNw_1                  += pt;
+             ++k1;
+             if (k1>=arraySize)
+               {
+                 AliError(Form("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
+                 return;
+               }
+           }
+        }
+      
+      if (!_sameFilter && _requestedCharge_2 == charge &&
+          eta      >=  -0.8 &&
+         eta      <   0.8)        
+       {
+         iPhi   = int( phi/_width_phi_2);
+         if (iPhi<0 || iPhi>=_nBins_phi_2 ) 
+           {
+             AliWarning("AliAnalysisTaskDptDptQA::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
+             return;
+           }
          
-         if (DCAZ     <  _dcaZMin || 
-             DCAZ     >  _dcaZMax ||
-             DCAXY    >  _dcaXYMax ) continue; 
+         iEta    = int((eta-_min_eta_2)/_width_eta_2);
+         if (iEta<0 || iEta>=_nBins_eta_2) 
+           {
+             AliWarning(Form("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
+             continue;
+           }
+         iPt     = int((pt -_min_pt_2 )/_width_pt_2 ); 
+         if (iPt<0  || iPt >=_nBins_pt_2)
+           {
+             AliWarning(Form("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
+             continue;
+           }
          
-         //==== QA ===========================
-         _dcaz->Fill(DCAZ);
-         _dcaxy->Fill(DCAXY);
-         _etadis->Fill(eta);
-         _phidis->Fill(phi);
-         _spectra->Fill(pt);
-         //===================================
-
-         if(TMath::Abs(nsigmapion) < fNSigmaCut && charge == -1) {
-           Kaon1Charge[nKaon1s] = q;
-           Kaon1Px[nKaon1s] = mom.X();
-           Kaon1Py[nKaon1s] = mom.Y();
-           Kaon1Pz[nKaon1s] = mom.Z();
-           nKaon1s++;
-
-         }
-         if(TMath::Abs(nsigmapion) < fNSigmaCut && charge == 1) {
-           Kaon2Charge[nKaon2s] = q;
-           Kaon2Px[nKaon1s] = mom.X();
-           Kaon2Py[nKaon1s] = mom.Y();
-           Kaon2Pz[nKaon1s] = mom.Z();
-           nKaon2s++;
-
-         }
-       
-         TLorentzVector pion1(0,0,0,0);
-         TLorentzVector pion2(0,0,0,0);
-         TLorentzVector Kaon(0,0,0,0);
-
-         for(int i=0; i<nKaon1s; i++)
+         iEtaPhi = iEta*_nBins_phi_2+iPhi;
+         iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
+         if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
            {
-             pion1.SetXYZM(Kaon1Px[i],Kaon1Py[i],Kaon1Pz[i],0.139570);
-             for(int j=0; j<nKaon2s; j++)
+             AliWarning("AliAnalysisTaskDptDptQA::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
+             continue;
+           }
+         
+         //cout << " iEtaPhi:" << iEtaPhi << "  _nBins_etaPhi_1: "<< _nBins_etaPhi_1<< "  _nBins_etaPhi_2: "<< _nBins_etaPhi_2<< endl;
+         //if (_useEbyECorrections)  corr = (charge>0) ? _correction_2p[iEtaPhi] : _correction_2m[iEtaPhi];      else  
+         //cout << "_correctionWeight_2:" << _correctionWeight_2 << endl;
+         if (_correctionWeight_2)
+           corr = _correctionWeight_2[iZEtaPhiPt];
+         else
+           corr = 1;
+         //dpt = pt - (charge>0) ? _avgPt_vsEtaPhi_2p[iEtaPhi] : _avgPt_vsEtaPhi_2m[iEtaPhi];       
+         
+         if (_singlesOnly)
+           {
+             //_dedxVsP_2->Fill(p*q,dedx);
+             __n1_2_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
+             __n1_2_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
+           }
+         else
+           {
+             corrPt                      = corr*pt;
+             _id_2[k2]                   = iTrack;         //cout << "step 1" << endl;
+             _charge_2[k2]               = charge;         //cout << "step 2" << endl;
+             _iEtaPhi_2[k2]              = iEtaPhi;        //cout << "step 3" << endl;
+             _iPt_2[k2]                  = iPt;            //cout << "step 4" << endl;
+             _pt_2[k2]                   = pt;             //cout << "step 5" << endl;
+             _px_2[k2]                   = px;             //cout << "step 6" << endl;
+             _py_2[k2]                   = py;             //cout << "step 7" << endl;
+             _pz_2[k2]                   = pz;             //cout << "step 8" << endl;
+             _correction_2[k2]           = corr;           //cout << "step 9" << endl;
+             __n1_2                      += corr;          //cout << "step 10" << endl;
+             __s1pt_2                    += corrPt;        //cout << "step 13" << endl;
+             __n1Nw_2                    += 1;
+             __n1_2_vsEtaPhi[iEtaPhi]    += corr;          //cout << "step 11" << endl;
+             __s1pt_2_vsEtaPhi[iEtaPhi]  += corrPt;        //cout << "step 14" << endl;
+             __s1ptNw_2                  += pt;
+             ++k2;
+             if (k2>=arraySize)
                {
-                 pion2.SetXYZM(Kaon2Px[j],Kaon2Py[j],Kaon2Pz[j],0.139570);
-                 Kaon = pion1 + pion2;
-                 float rapidity = Kaon.Rapidity();
-                 float pairpt = Kaon.Pt();
-                 if( rapidity < -0.50 || rapidity > 0.50 ) continue;
-                 _invMass->Fill(pairpt, Kaon.M());
+                 AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize)); 
+                 return;
                }
            }
+         
+         //cout << "done with track" << endl;
+       } //iTrack
+      } //aod 
+    }
+  
+  
+  //cout << "Filling histograms now" << endl;
+  _m0->Fill(_mult0);
+  _m1->Fill(_mult1);
+  _m2->Fill(_mult2);
+  _m3->Fill(_mult3);
+  _m4->Fill(_mult4);
+  _m5->Fill(_mult5);
+  _m6->Fill(_mult6);
+  _vertexZ->Fill(vertexZ);
+  
+  if (_singlesOnly)
+    {
+      // nothing to do here.
+    }
+  else
+    {
+      if (_sameFilter)
+       {
+         _n1_1_vsM->Fill(centrality,      __n1_1);
+         _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
+         _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
+         _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
+         _n1_2_vsM->Fill(centrality,      __n1_1);
+         _s1pt_2_vsM->Fill(centrality,    __s1pt_1);
+         _n1Nw_2_vsM->Fill(centrality,    __n1Nw_1);
+         _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_1);
+         // reset pair counters
+         __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
+         __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
+         if (_field>0)
+           {
+             for (int i1=0; i1<k1; i1++)
+               {
+                 ////cout << "         i1:" << i1 << endl;
+                 id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
+                 q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
+                 iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
+                 iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
+                 corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
+                 pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
+                 dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
+                 //1 and 2
+                 for (int i2=i1+1; i2<k1; i2++)
+                   {        
+                     ////cout << "         i2:" << i2 << endl;
+                     id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
+                     if (id_1!=id_2)
+                       {
+                         q_2       = _charge_1[i2];     ////cout << "        q_1:" << q_1 << endl;
+                         iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
+                         iPt_2     = _iPt_1[i2];        ////cout << "      iPt_1:" << iPt_1 << endl;
+                         corr_2    = _correction_1[i2]; ////cout << "     corr_1:" << corr_1 << endl;
+                         pt_2      = _pt_1[i2];         ////cout << "       pt_1:" << pt_1 << endl;
+                         dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
+                         corr      = corr_1*corr_2;
+                         if (q_2>q_1 || (q_1>0 && q_2>0 && pt_2<=pt_1) || (q_1<0 && q_2<0 && pt_2>=pt_1))
+                           {
+                             ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
+                           }
+                         else // swap particles
+                           {
+                             ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
+                           }
+                         
+                         __n2_12                  += corr;
+                         __n2_12_vsEtaPhi[ij]     += corr;
+                         ptpt                     = pt_1*pt_2;
+                         __s2ptpt_12              += corr*ptpt;
+                         __s2PtN_12               += corr*pt_1;
+                         __s2NPt_12               += corr*pt_2;
+                         __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
+                         __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
+                         __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
+                         __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
+                         
+                         __n2Nw_12                  += 1;
+                         __s2ptptNw_12              += ptpt;
+                         __s2PtNNw_12               += pt_1;
+                         __s2NPtNw_12               += pt_2;
+                         
+                       }
+                   } //i2       
+               } //i1       
+           }
+         else // field<0
+           {
+             for (int i1=0; i1<k1; i1++)
+               {
+                 ////cout << "         i1:" << i1 << endl;
+                 id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
+                 q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
+                 iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
+                 iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
+                 corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
+                 pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
+                 dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
+                 //1 and 2
+                 for (int i2=i1+1; i2<k1; i2++)
+                   {        
+                     ////cout << "         i2:" << i2 << endl;
+                     id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
+                     if (id_1!=id_2)
+                       {
+                         q_2       = _charge_1[i2];     ////cout << "        q_2:" << q_2 << endl;
+                         iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
+                         iPt_2     = _iPt_1[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
+                         corr_2    = _correction_1[i2]; ////cout << "     corr_2:" << corr_2 << endl;
+                         pt_2      = _pt_1[i2];         ////cout << "       pt_2:" << pt_2 << endl;
+                         dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
+                         corr      = corr_1*corr_2;
+                         if ( q_2<q_1 || (q_1>0 && q_2>0 && pt_2>=pt_1) || (q_1<0 && q_2<0 && pt_2<=pt_1))
+                           {
+                             ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
+                           }
+                         else // swap particles
+                           {
+                             ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
+                           }
+                         
+                         __n2_12                  += corr;
+                         __n2_12_vsEtaPhi[ij]     += corr;
+                         ptpt                     = pt_1*pt_2;
+                         __s2ptpt_12              += corr*ptpt;
+                         __s2PtN_12               += corr*pt_1;
+                         __s2NPt_12               += corr*pt_2;
+                         __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
+                         __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
+                         __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
+                         __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
+                         
+                         __n2Nw_12                  += 1;
+                         __s2ptptNw_12              += ptpt;
+                         __s2PtNNw_12               += pt_1;
+                         __s2NPtNw_12               += pt_2;
+                         
+                       }
+                   } //i2       
+               } //i1  
+           }
+       }
+      else  // filter 1 and 2 are different -- must do all particle pairs...
+       {
+         _n1_1_vsM->Fill(centrality,      __n1_1);
+         _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
+         _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
+         _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
+         _n1_2_vsM->Fill(centrality,      __n1_2);
+         _s1pt_2_vsM->Fill(centrality,    __s1pt_2);
+         _n1Nw_2_vsM->Fill(centrality,    __n1Nw_2);
+         _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_2);
+         // reset pair counters
+         __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
+         __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
+         for (int i1=0; i1<k1; i1++)
+           {
+             ////cout << "         i1:" << i1 << endl;
+             id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
+             q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
+             iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
+             iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
+             corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
+             pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
+             px_1      = _px_1[i1];          ////cout << "      px_1:" << px_1 << endl;
+             py_1      = _py_1[i1];          ////cout << "      py_1:" << py_1 << endl;
+             pz_1      = _pz_1[i1];          ////cout << "      pz_1:" << pz_1 << endl;
+             dedx_1    = _dedx_1[i1];        ////cout << "     dedx_1:" << dedx_1 << endl;
+             
+             //1 and 2
+             for (int i2=0; i2<k2; i2++)
+               {        
+                 ////cout << "         i2:" << i2 << endl;
+                 id_2   = _id_2[i2];              ////cout << "       id_2:" << id_2 << endl;
+                 if (id_1!=id_2)  // exclude auto correlation
+                   {
+                     q_2       = _charge_2[i2];     ////cout << "        q_2:" << q_2 << endl;
+                     iEtaPhi_2 = _iEtaPhi_2[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
+                     iPt_2     = _iPt_2[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
+                     corr_2    = _correction_2[i2]; ////cout << "     corr_2:" << corr_2 << endl;
+                     pt_2      = _pt_2[i2];         ////cout << "       pt_2:" << pt_2 << endl;
+                     px_2      = _px_2[i2];          ////cout << "      px_2:" << px_2 << endl;
+                     py_2      = _py_2[i2];          ////cout << "      py_2:" << py_2 << endl;
+                     pz_2      = _pz_2[i2];          ////cout << "      pz_2:" << pz_2 << endl;
+                     dedx_2    = _dedx_2[i2];        ////cout << "     dedx_2:" << dedx_2 << endl;
+                     if (_rejectPairConversion)
+                       {
+                         float e1Sq = massElecSq + pt_1*pt_1 + pz_1*pz_1;
+                         float e2Sq = massElecSq + pt_2*pt_2 + pz_2*pz_2;
+                         float mInvSq = 2*(massElecSq + sqrt(e1Sq*e2Sq) - px_1*px_2 - py_1*py_2 - pz_1*pz_2 );
+                         float mInv = sqrt(mInvSq);
+                         _invMass->Fill(mInv);
+                         if (mInv<0.5)
+                           {
+                             if (dedx_1>75. && dedx_2>75.)
+                               {
+                                 _invMassElec->Fill(mInv);
+                                 if (mInv<0.05) continue;
+                               }
+                           }
+                       }
+                     corr      = corr_1*corr_2;
+                     ij        = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
+                     __n2_12                  += corr;
+                     __n2_12_vsEtaPhi[ij]     += corr;
+                     ptpt                     = pt_1*pt_2;
+                     __s2ptpt_12              += corr*ptpt;
+                     __s2PtN_12               += corr*pt_1;
+                     __s2NPt_12               += corr*pt_2;
+                     __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
+                     __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
+                     __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
+                     __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;         
+                     __n2Nw_12                  += 1;
+                     __s2ptptNw_12              += ptpt;
+                     __s2PtNNw_12               += pt_1;
+                     __s2NPtNw_12               += pt_2;
+                     
+                   }
+               } //i2       
+           } //i1         
+       }
+      
+      _n2_12_vsM->Fill(centrality,     __n2_12);
+      _s2PtPt_12_vsM->Fill(centrality, __s2ptpt_12); 
+      _s2PtN_12_vsM->Fill(centrality,  __s2NPt_12);
+      _s2NPt_12_vsM->Fill(centrality,  __s2PtN_12);
+      
+      _n2Nw_12_vsM->Fill(centrality,     __n2Nw_12);
+      _s2PtPtNw_12_vsM->Fill(centrality, __s2ptptNw_12); 
+      _s2PtNNw_12_vsM->Fill(centrality,  __s2NPtNw_12);
+      _s2NPtNw_12_vsM->Fill(centrality,  __s2PtNNw_12);
+      
+    }
+  
+  
+  cout << "AliAnalysisTaskDptDptQA::UserExec()   -----------------Event Done " << endl;
+  PostData(0,_outputHistoList);
   
-         } //iTrack
-    } //aod 
 }
 
 void   AliAnalysisTaskDptDptQA::FinishTaskOutput()
 {
-  AliInfo("AliAnalysisTaskDptDptQA::FinishTaskOutput() Starting.");
-  Printf("= 0 ====================================================================");
+  cout << "AliAnalysisTaskDptDptQA::FinishTaskOutput() Starting." << endl;
+  cout << "= 0 ====================================================================" << endl;
   finalizeHistograms();
-  AliInfo("= 1 ====================================================================");
+  cout << "= 1 ====================================================================" << endl;
   PostData(0,_outputHistoList);
-  AliInfo("= 2 ====================================================================");
-  AliInfo("AliAnalysisTaskDptDptQA::FinishTaskOutput() Done.");
+  cout << "= 2 ====================================================================" << endl;
+  cout << "AliAnalysisTaskDptDptQA::FinishTaskOutput() Done." << endl;
 }
 
 void   AliAnalysisTaskDptDptQA::Terminate(Option_t* /*option*/)
 {
-  AliInfo("AliAnalysisTaskDptDptQA::Terminate() Starting/Done.");
+  cout << "AliAnalysisTaskDptDptQA::Terminate() Starting/Done." << endl;
 }
 
 
@@ -1569,7 +1911,7 @@ void   AliAnalysisTaskDptDptQA::addToList(TH1 *h)
     _outputHistoList->Add(h);
     }
   else
-    AliInfo("addToList(TH1 *h) _outputHistoList is null!!!!! Should abort ship");
+    cout << "addToList(TH1 *h) _outputHistoList is null!!!!! Shoudl abort ship" << endl;
 
 }