PWGCF/Correlations/DPhi/AliDptDptInMC.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
  16 /* $Id:$ */
  17
  18 #include "TChain.h"
  19 #include "TList.h"
  20 #include "TFile.h"
  21 #include "TTree.h"
  22 #include "TH1D.h"
  23 #include "TH2D.h"
  24 #include "TH3D.h"
  25 #include "THnSparse.h"
  26 #include "TCanvas.h"
  27
  28 #include <TROOT.h>
  29 #include <TChain.h>
  30 #include <TFile.h>
  31 #include <TList.h>
  32 #include <TMath.h>
  33 #include <TTree.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TH3F.h>
  37 #include <TProfile.h>
  38 #include <TH1D.h>
  39 #include <TH2D.h>
  40 #include <TH3D.h>
  41 #include <TRandom.h>
  42 #include "AliAnalysisManager.h"
  43
  44 #include "AliAODHandler.h"
  45 #include "AliAODInputHandler.h"
  46 #include "AliInputEventHandler.h"
  47 #include "AliLog.h"
  48 #include "AliESDEvent.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliMultiplicity.h"
  51 #include "AliCentrality.h"
  52 #include "AliDptDptInMC.h"
  53
  54 #include "AliPID.h"
  55 #include "AliPIDResponse.h"
  56
  57 #include "AliESDVertex.h"
  58 #include "AliESDEvent.h"
  59 #include "AliESDInputHandler.h"
  60 #include "AliAODEvent.h"
  61 #include "AliAODTrack.h"
  62 #include "AliAODInputHandler.h"
  63 #include "AliESD.h"
  64 #include "AliESDEvent.h"
  65 #include "AliAODEvent.h"
  66 #include "AliStack.h"
  67 #include "AliESDtrackCuts.h"
  68 #include "AliAODMCHeader.h"
  69 #include "AliAODMCParticle.h"
  70 #include "AliAODMCParticle.h"
  71 #include "AliAODMCHeader.h"
  72 #include "AliAODHeader.h"
  73
  74 #include "AliGenHijingEventHeader.h"
  75 #include "AliGenEventHeader.h"
  76 #include "AliPID.h"
  77 #include "AliAODPid.h"
  78 #include "AliPIDResponse.h"
  79 #include "AliAODpidUtil.h"
  80 #include "AliPIDCombined.h"
  81
  82 using std::cout;
  83 using std::endl;
  84
  85 ClassImp(AliDptDptInMC)
  86
  87 AliDptDptInMC::AliDptDptInMC()
  88 : AliAnalysisTaskSE(),
  89 fAODEvent(0),
  90 fESDEvent(0),             //! ESD Event
  91 fInputHandler(0),
  92 fPIDResponse(0x0),
  93 _outputHistoList(0),
  94 fArrayMC (0),
  95 _twoPi         ( 2.0 * 3.1415927),
  96 _eventCount    ( 0),
  97 _debugLevel    ( 0),
  98 _singlesOnly   ( 0),
  99 _useWeights    ( 0),
 100 _sameFilter    ( false),
 101 _rejectPileup  ( 1),
 102 _rejectPairConversion ( 0),
 103 _vertexZMin           ( -10),
 104 _vertexZMax           (  10),
 105 _vertexXYMin          ( -10),
 106 _vertexXYMax          (  10),
 107 _centralityMethod     (  4),
 108 _centralityMin        (  0.),
 109 _centralityMax        (  0.),
 110 _requestedCharge_1    (   1),
 111 _requestedCharge_2    (  -1),
 112 _dcaZMin              ( -3),
 113 _dcaZMax              (  3.),
 114 _dcaXYMin             ( -2.4),
 115 _dcaXYMax             (  2.4),
 116 _dedxMin              ( 0),
 117 _dedxMax              ( 100000),
 118 _nClusterMin          ( 80),
 119 _trackFilterBit       (0),
 120 fNSigmaCut            (3.),
 121 fAnalysisType ("MCAOD"),
 122 fExcludeResonancesInMC (kFALSE),
 123 fExcludeElectronsInMC (kFALSE),
 124 _field    ( 1.),
 125 _nTracks  ( 0 ),
 126 _mult0    ( 0 ),
 127 _mult1    ( 0 ),
 128 _mult2    ( 0 ),
 129 _mult3    ( 0 ),
 130 _mult4    ( 0 ),
 131 _mult4a    ( 0 ),
 132 _mult5    ( 0 ),
 133 _mult6    ( 0 ),
 134 arraySize ( 2000),
 135 _id_1(0),
 136 _charge_1(0),
 137 _iEtaPhi_1(0),
 138 _iPt_1(0),
 139 _pt_1(0),
 140 _px_1(0),
 141 _py_1(0),
 142 _pz_1(0),
 143 _correction_1(0),
 144 _dedx_1(0),
 145 _id_2(0),
 146 _charge_2(0),
 147 _iEtaPhi_2(0),
 148 _iPt_2(0),
 149 _pt_2(0),
 150 _px_2(0),
 151 _py_2(0),
 152 _pz_2(0),
 153 _correction_2(0),
 154 _dedx_2(0),
 155 _correctionWeight_1(0),
 156 _correctionWeight_2(0),
 157 _nBins_M0(500),       _min_M0(0),        _max_M0(10000),          _width_M0(20),
 158 _nBins_M1(500),       _min_M1(0),        _max_M1(10000),          _width_M1(20),
 159 _nBins_M2(500),       _min_M2(0),        _max_M2(10000),          _width_M2(20),
 160 _nBins_M3(500),       _min_M3(0),        _max_M3(10000),          _width_M3(20),
 161 _nBins_M4(100),       _min_M4(0),        _max_M4(1),              _width_M4(0.01),
 162 _nBins_M5(100),       _min_M5(0),        _max_M5(1),              _width_M5(0.01),
 163 _nBins_M6(100),       _min_M6(0),        _max_M6(1),              _width_M6(0.01),
 164 _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ(0.5),
 165
 166 _nBins_pt_1(18),      _min_pt_1(0.2),    _max_pt_1(2.0),          _width_pt_1(0.1),
 167 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
 168 //_nBins_eta_1(18),
 169 //_min_eta_1(-0.9),  _max_eta_1(0.9),           _width_eta_1(0.1),
 170   _nBins_eta_1(0),
 171   _min_eta_1(0),  _max_eta_1(0),           _width_eta_1(0.1),
 172
 173 _nBins_etaPhi_1(0),
 174 _nBins_etaPhiPt_1(0),
 175 _nBins_zEtaPhiPt_1(0),
 176
 177 _nBins_pt_2(18),     _min_pt_2(0.2),     _max_pt_2(2.0),          _width_pt_2(0.1),
 178 _nBins_phi_2(72),    _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
 179 //_nBins_eta_2(18),
 180 //_min_eta_2(-0.9),     _max_eta_2(0.9),           _width_eta_2(0.1),
 181   _nBins_eta_2(0),
 182   _min_eta_2(0),     _max_eta_2(0),           _width_eta_2(0.1),
 183
 184
 185
 186 _nBins_etaPhi_2(0),
 187 _nBins_etaPhiPt_2(0),
 188 _nBins_zEtaPhiPt_2(0),
 189 _nBins_etaPhi_12(0),
 190 __n1_1(0),
 191 __n1_2(0),
 192 __n2_12(0),
 193 __s1pt_1(0),
 194 __s1pt_2(0),
 195 __s2ptpt_12(0),
 196 __s2NPt_12(0),
 197 __s2PtN_12(0),
 198 __n1Nw_1(0),
 199 __n1Nw_2(0),
 200 __n2Nw_12(0),
 201 __s1ptNw_1(0),
 202 __s1ptNw_2(0),
 203 __s2ptptNw_12(0),
 204 __s2NPtNw_12(0),
 205 __s2PtNNw_12(0),
 206 __n1_1_vsPt(0),
 207 __n1_1_vsEtaPhi(0),
 208 __s1pt_1_vsEtaPhi(0),
 209 __n1_1_vsZEtaPhiPt(0),
 210 __n1_2_vsPt(0),
 211 __n1_2_vsEtaPhi(0),
 212 __s1pt_2_vsEtaPhi(0),
 213 __n1_2_vsZEtaPhiPt(0),
 214 __n2_12_vsPtPt(0),
 215 __n2_12_vsEtaPhi(0),
 216 __s2ptpt_12_vsEtaPhi(0),
 217 __s2PtN_12_vsEtaPhi(0),
 218 __s2NPt_12_vsEtaPhi(0),
 219 _weight_1      ( 0    ),
 220 _weight_2      ( 0    ),
 221 _eventAccounting ( 0),
 222 _m0 ( 0),
 223 _m1 ( 0),
 224 _m2 ( 0),
 225 _m3 ( 0),
 226 _m4 ( 0),
 227 _m5 ( 0),
 228 _m6 ( 0),
 229 _vertexZ ( 0),
 230
 231 _Ncluster1  ( 0),
 232 _Ncluster2  ( 0),
 233 _etadis ( 0),
 234 _phidis ( 0),
 235 _dcaz   ( 0),
 236 _dcaxy  ( 0),
 237 _etadis2 ( 0),
 238 _phidis2 ( 0),
 239 _dcaz2   ( 0),
 240 _dcaxy2  ( 0),
 241 _etadis3 ( 0),
 242 _phidis3 ( 0),
 243 _dcaz3   ( 0),
 244 _dcaxy3  ( 0),
 245 _n1_1_vsPt         ( 0),
 246 _n1_1_vsEtaVsPhi   ( 0),
 247 _s1pt_1_vsEtaVsPhi ( 0),
 248 _n1_1_vsZVsEtaVsPhiVsPt ( 0),
 249 _n1_1_vsM          ( 0),  // w/ weight
 250 _s1pt_1_vsM        ( 0),
 251 _n1Nw_1_vsM        ( 0), // w/o weight
 252 _s1ptNw_1_vsM      ( 0),
 253 _dedxVsP_1         ( 0),
 254 _corrDedxVsP_1     ( 0),
 255 _betaVsP_1         ( 0),
 256 _n1_2_vsPt         ( 0),
 257 _n1_2_vsEtaVsPhi   ( 0),
 258 _s1pt_2_vsEtaVsPhi ( 0),
 259 _n1_2_vsZVsEtaVsPhiVsPt ( 0),
 260 _n1_2_vsM          ( 0),
 261 _s1pt_2_vsM        ( 0),
 262 _n1Nw_2_vsM        ( 0),
 263 _s1ptNw_2_vsM      ( 0),
 264 _dedxVsP_2         ( 0),
 265 _corrDedxVsP_2     ( 0),
 266 _betaVsP_2         ( 0),
 267 _n2_12_vsEtaPhi    ( 0),
 268 _n2_12_vsPtVsPt    ( 0),
 269 _s2PtPt_12_vsEtaPhi( 0),
 270 _s2PtN_12_vsEtaPhi ( 0),
 271 _s2NPt_12_vsEtaPhi ( 0),
 272 _n2_12_vsM         ( 0),
 273 _s2PtPt_12_vsM     ( 0),
 274 _s2PtN_12_vsM      ( 0),
 275 _s2NPt_12_vsM      ( 0),
 276 _n2Nw_12_vsM       ( 0),
 277 _s2PtPtNw_12_vsM   ( 0),
 278 _s2PtNNw_12_vsM    ( 0),
 279 _s2NPtNw_12_vsM    ( 0),
 280 _invMass           ( 0),
 281 _invMassElec       ( 0),
 282 n1Name("NA"),
 283 n1NwName("NA"),
 284 n2Name("NA"),
 285 n2NwName("NA"),
 286 n3Name("NA"),
 287 n1n1Name("NA"),
 288 n1n1n1Name("NA"),
 289 n2n1Name("NA"),
 290 r1Name("NA"),
 291 r2Name("NA"),
 292 r3Name("NA"),
 293 r2r1Name("NA"),
 294 c2Name("NA"),
 295 c3Name("NA"),
 296 d3Name("NA"),
 297 p3Name("NA"),
 298 cName("NA"),
 299
 300 intR2Name("NA"),
 301 binCorrName("NA"),
 302 intBinCorrName("NA"),
 303
 304 countsName("NA"),
 305 part_1_Name("NA"),
 306 part_2_Name("NA"),
 307 part_3_Name("NA"),
 308 pair_12_Name("NA"),
 309 pair_13_Name("NA"),
 310 pair_23_Name("NA"),
 311 tripletName("NA"),
 312
 313 avg("NA"),
 314 avgName("NA"),
 315 sumName("NA"),
 316 s1ptName("NA"),
 317 s1ptNwName("NA"),
 318 s1DptName("NA"),
 319
 320 s2PtPtName("NA"),
 321 s2NPtName("NA"),
 322 s2PtNName("NA"),
 323 s2DptDptName("NA"),
 324
 325 s2PtPtNwName("NA"),
 326 s2NPtNwName("NA"),
 327 s2PtNNwName("NA"),
 328
 329 ptName("NA"),
 330 ptptName("NA"),
 331 pt1pt1Name("NA"),
 332 DptName("NA"),
 333 DptDptName("NA"),
 334 RDptDptName("NA"),
 335 nPtName("NA"),
 336 ptNName("NA"),
 337 seanName("NA"),
 338
 339 _title_counts("NA"),
 340
 341 _title_m0("NA"),
 342 _title_m1("NA"),
 343 _title_m2("NA"),
 344 _title_m3("NA"),
 345 _title_m4("NA"),
 346 _title_m5("NA"),
 347 _title_m6("NA"),
 348
 349 _title_eta_1("NA"),
 350 _title_phi_1("NA"),
 351 _title_pt_1("NA"),
 352 _title_etaPhi_1("NA"),
 353 _title_n_1("NA"),
 354 _title_SumPt_1("NA"),
 355 _title_AvgPt_1("NA"),
 356 _title_AvgN_1("NA"),
 357 _title_AvgSumPt_1("NA"),
 358
 359 _title_eta_2("NA"),
 360 _title_phi_2("NA"),
 361 _title_pt_2("NA"),
 362 _title_etaPhi_2("NA"),
 363 _title_n_2("NA"),
 364 _title_SumPt_2("NA"),
 365 _title_AvgPt_2("NA"),
 366 _title_AvgN_2("NA"),
 367 _title_AvgSumPt_2("NA"),
 368
 369 _title_etaPhi_12("NA"),
 370
 371 _title_AvgN2_12("NA"),
 372 _title_AvgSumPtPt_12("NA"),
 373 _title_AvgSumPtN_12("NA"),
 374 _title_AvgNSumPt_12("NA"),
 375
 376 vsZ("NA"),
 377 vsM("NA"),
 378 vsPt("NA"),
 379 vsPhi("NA"),
 380 vsEta("NA"),
 381 vsEtaPhi("NA"),
 382 vsPtVsPt("NA")
 383 {
 384   printf("Default constructor called \n");
 385
 386   printf("passed \n ");
 387
 388 }
 389
 390 AliDptDptInMC::AliDptDptInMC(const TString & name)
 391 : AliAnalysisTaskSE(name),
 392 fAODEvent(0),
 393 fESDEvent(0),
 394 fInputHandler(0),
 395 fPIDResponse(0),
 396 _outputHistoList(0),
 397 fArrayMC (0),
 398 _twoPi         ( 2.0 * 3.1415927),
 399 _eventCount    ( 0),
 400 _debugLevel    ( 0),
 401 _singlesOnly   ( 0),
 402 _useWeights    ( 0),
 403 _sameFilter    ( false),
 404 _rejectPileup  ( 1),
 405 _rejectPairConversion ( 0),
 406 _vertexZMin           ( -10.),
 407 _vertexZMax           (  10.),
 408 _vertexXYMin          ( -10.),
 409 _vertexXYMax          (  10.),
 410 _centralityMethod     (  4),
 411 _centralityMin        (  0.),
 412 _centralityMax        (  1.),
 413 _requestedCharge_1    (   1),
 414 _requestedCharge_2    (  -1),
 415 _dcaZMin              ( -3),
 416 _dcaZMax              (  3.),
 417 _dcaXYMin             ( -2.4),
 418 _dcaXYMax             (  2.4),
 419 _dedxMin              ( 0),
 420 _dedxMax              ( 100000),
 421 _nClusterMin          ( 80),
 422 _trackFilterBit       ( 0),
 423 fNSigmaCut            ( 3.),
 424
 425 fAnalysisType ("MCAOD"),
 426 fExcludeResonancesInMC (kFALSE),
 427 fExcludeElectronsInMC (kFALSE),
 428
 429 _field    ( 1.),
 430 _nTracks  ( 0 ),
 431 _mult0    ( 0 ),
 432 _mult1    ( 0 ),
 433 _mult2    ( 0 ),
 434 _mult3    ( 0 ),
 435 _mult4    ( 0 ),
 436 _mult4a    ( 0 ),
 437 _mult5    ( 0 ),
 438 _mult6    ( 0 ),
 439 arraySize ( 2000),
 440 _id_1(0),
 441 _charge_1(0),
 442 _iEtaPhi_1(0),
 443 _iPt_1(0),
 444 _pt_1(0),
 445 _px_1(0),
 446 _py_1(0),
 447 _pz_1(0),
 448 _correction_1(0),
 449 _dedx_1(0),
 450 _id_2(0),
 451 _charge_2(0),
 452 _iEtaPhi_2(0),
 453 _iPt_2(0),
 454 _pt_2(0),
 455 _px_2(0),
 456 _py_2(0),
 457 _pz_2(0),
 458 _correction_2(0),
 459 _dedx_2(0),
 460 _correctionWeight_1(0),
 461 _correctionWeight_2(0),
 462 _nBins_M0(500),       _min_M0(0),        _max_M0(10000),          _width_M0(20),
 463 _nBins_M1(500),       _min_M1(0),        _max_M1(10000),          _width_M1(20),
 464 _nBins_M2(500),       _min_M2(0),        _max_M2(10000),          _width_M2(20),
 465 _nBins_M3(500),       _min_M3(0),        _max_M3(10000),          _width_M3(20),
 466 _nBins_M4(100),       _min_M4(0),        _max_M4(1),              _width_M4(0.01),
 467 _nBins_M5(100),       _min_M5(0),        _max_M5(1),              _width_M5(0.01),
 468 _nBins_M6(100),       _min_M6(0),        _max_M6(1),              _width_M6(0.01),
 469 _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ(0.5),
 470
 471 _nBins_pt_1(18),      _min_pt_1(0.2),    _max_pt_1(2.0),          _width_pt_1(0.1),
 472 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
 473   //_nBins_eta_1(18),
 474   //_min_eta_1(-0.9),    _max_eta_1(0.9),           _width_eta_1(0.1),
 475   _nBins_eta_1(0),
 476   _min_eta_1(0),    _max_eta_1(0),           _width_eta_1(0.1),
 477
 478 _nBins_etaPhi_1(0),
 479 _nBins_etaPhiPt_1(0),
 480 _nBins_zEtaPhiPt_1(0),
 481
 482 _nBins_pt_2(18),     _min_pt_2(0.2),     _max_pt_2(2.0),          _width_pt_2(0.1),
 483 _nBins_phi_2(72),    _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
 484   //_nBins_eta_2(18),
 485   //_min_eta_2(-0.9),     _max_eta_2(0.9),           _width_eta_2(0.1),
 486   _nBins_eta_2(0),
 487   _min_eta_2(0),     _max_eta_2(0),           _width_eta_2(0.1),
 488
 489
 490 _nBins_etaPhi_2(0),
 491 _nBins_etaPhiPt_2(0),
 492 _nBins_zEtaPhiPt_2(0),
 493 _nBins_etaPhi_12(0),
 494 __n1_1(0),
 495 __n1_2(0),
 496 __n2_12(0),
 497 __s1pt_1(0),
 498 __s1pt_2(0),
 499 __s2ptpt_12(0),
 500 __s2NPt_12(0),
 501 __s2PtN_12(0),
 502 __n1Nw_1(0),
 503 __n1Nw_2(0),
 504 __n2Nw_12(0),
 505 __s1ptNw_1(0),
 506 __s1ptNw_2(0),
 507 __s2ptptNw_12(0),
 508 __s2NPtNw_12(0),
 509 __s2PtNNw_12(0),
 510 __n1_1_vsPt(0),
 511 __n1_1_vsEtaPhi(0),
 512 __s1pt_1_vsEtaPhi(0),
 513 __n1_1_vsZEtaPhiPt(0),
 514 __n1_2_vsPt(0),
 515 __n1_2_vsEtaPhi(0),
 516 __s1pt_2_vsEtaPhi(0),
 517 __n1_2_vsZEtaPhiPt(0),
 518 __n2_12_vsPtPt(0),
 519 __n2_12_vsEtaPhi(0),
 520 __s2ptpt_12_vsEtaPhi(0),
 521 __s2PtN_12_vsEtaPhi(0),
 522 __s2NPt_12_vsEtaPhi(0),
 523 _weight_1        ( 0    ),
 524 _weight_2        ( 0    ),
 525 _eventAccounting ( 0),
 526 _m0 ( 0),
 527 _m1 ( 0),
 528 _m2 ( 0),
 529 _m3 ( 0),
 530 _m4 ( 0),
 531 _m5 ( 0),
 532 _m6 ( 0),
 533 _vertexZ ( 0),
 534 _Ncluster1  ( 0),
 535 _Ncluster2  ( 0),
 536 _etadis ( 0),
 537 _phidis ( 0),
 538
 539 _dcaz ( 0),
 540 _dcaxy ( 0),
 541 _etadis2 ( 0),
 542 _phidis2 ( 0),
 543 _dcaz2   ( 0),
 544 _dcaxy2  ( 0),
 545 _etadis3 ( 0),
 546 _phidis3 ( 0),
 547 _dcaz3   ( 0),
 548 _dcaxy3  ( 0),
 549 _n1_1_vsPt         ( 0),
 550 _n1_1_vsEtaVsPhi   ( 0),
 551 _s1pt_1_vsEtaVsPhi ( 0),
 552 _n1_1_vsZVsEtaVsPhiVsPt ( 0),
 553 _n1_1_vsM          ( 0),  // w/ weight
 554 _s1pt_1_vsM        ( 0),
 555 _n1Nw_1_vsM        ( 0), // w/o weight
 556 _s1ptNw_1_vsM      ( 0),
 557 _dedxVsP_1         ( 0),
 558 _corrDedxVsP_1     ( 0),
 559 _betaVsP_1         ( 0),
 560 _n1_2_vsPt         ( 0),
 561 _n1_2_vsEtaVsPhi   ( 0),
 562 _s1pt_2_vsEtaVsPhi ( 0),
 563 _n1_2_vsZVsEtaVsPhiVsPt ( 0),
 564 _n1_2_vsM          ( 0),
 565 _s1pt_2_vsM        ( 0),
 566 _n1Nw_2_vsM        ( 0),
 567 _s1ptNw_2_vsM      ( 0),
 568 _dedxVsP_2         ( 0),
 569 _corrDedxVsP_2     ( 0),
 570 _betaVsP_2         ( 0),
 571 _n2_12_vsEtaPhi    ( 0),
 572 _n2_12_vsPtVsPt    ( 0),
 573 _s2PtPt_12_vsEtaPhi( 0),
 574 _s2PtN_12_vsEtaPhi ( 0),
 575 _s2NPt_12_vsEtaPhi ( 0),
 576 _n2_12_vsM         ( 0),
 577 _s2PtPt_12_vsM     ( 0),
 578 _s2PtN_12_vsM      ( 0),
 579 _s2NPt_12_vsM      ( 0),
 580 _n2Nw_12_vsM       ( 0),
 581 _s2PtPtNw_12_vsM   ( 0),
 582 _s2PtNNw_12_vsM    ( 0),
 583 _s2NPtNw_12_vsM    ( 0),
 584 _invMass           ( 0),
 585 _invMassElec       ( 0),
 586 n1Name("NA"),
 587 n1NwName("NA"),
 588 n2Name("NA"),
 589 n2NwName("NA"),
 590 n3Name("NA"),
 591 n1n1Name("NA"),
 592 n1n1n1Name("NA"),
 593 n2n1Name("NA"),
 594 r1Name("NA"),
 595 r2Name("NA"),
 596 r3Name("NA"),
 597 r2r1Name("NA"),
 598 c2Name("NA"),
 599 c3Name("NA"),
 600 d3Name("NA"),
 601 p3Name("NA"),
 602 cName("NA"),
 603
 604 intR2Name("NA"),
 605 binCorrName("NA"),
 606 intBinCorrName("NA"),
 607
 608 countsName("NA"),
 609 part_1_Name("NA"),
 610 part_2_Name("NA"),
 611 part_3_Name("NA"),
 612 pair_12_Name("NA"),
 613 pair_13_Name("NA"),
 614 pair_23_Name("NA"),
 615 tripletName("NA"),
 616
 617 avg("NA"),
 618 avgName("NA"),
 619 sumName("NA"),
 620 s1ptName("NA"),
 621 s1ptNwName("NA"),
 622 s1DptName("NA"),
 623
 624 s2PtPtName("NA"),
 625 s2NPtName("NA"),
 626 s2PtNName("NA"),
 627 s2DptDptName("NA"),
 628
 629 s2PtPtNwName("NA"),
 630 s2NPtNwName("NA"),
 631 s2PtNNwName("NA"),
 632
 633 ptName("NA"),
 634 ptptName("NA"),
 635 pt1pt1Name("NA"),
 636 DptName("NA"),
 637 DptDptName("NA"),
 638 RDptDptName("NA"),
 639 nPtName("NA"),
 640 ptNName("NA"),
 641 seanName("NA"),
 642
 643 _title_counts("NA"),
 644
 645 _title_m0("NA"),
 646 _title_m1("NA"),
 647 _title_m2("NA"),
 648 _title_m3("NA"),
 649 _title_m4("NA"),
 650 _title_m5("NA"),
 651 _title_m6("NA"),
 652
 653 _title_eta_1("NA"),
 654 _title_phi_1("NA"),
 655 _title_pt_1("NA"),
 656 _title_etaPhi_1("NA"),
 657 _title_n_1("NA"),
 658 _title_SumPt_1("NA"),
 659 _title_AvgPt_1("NA"),
 660 _title_AvgN_1("NA"),
 661 _title_AvgSumPt_1("NA"),
 662
 663 _title_eta_2("NA"),
 664 _title_phi_2("NA"),
 665 _title_pt_2("NA"),
 666 _title_etaPhi_2("NA"),
 667 _title_n_2("NA"),
 668 _title_SumPt_2("NA"),
 669 _title_AvgPt_2("NA"),
 670 _title_AvgN_2("NA"),
 671 _title_AvgSumPt_2("NA"),
 672
 673 _title_etaPhi_12("NA"),
 674
 675 _title_AvgN2_12("NA"),
 676 _title_AvgSumPtPt_12("NA"),
 677 _title_AvgSumPtN_12("NA"),
 678 _title_AvgNSumPt_12("NA"),
 679
 680 vsZ("NA"),
 681 vsM("NA"),
 682 vsPt("NA"),
 683 vsPhi("NA"),
 684 vsEta("NA"),
 685 vsEtaPhi("NA"),
 686 vsPtVsPt("NA")
 687 {
 688   printf("2nd constructor called ");
 689
 690   DefineOutput(0, TList::Class());
 691
 692   printf("passed  ");
 693
 694 }
 695
 696 AliDptDptInMC::~AliDptDptInMC()
 697 {
 698   /*
 699    delete _id_1;
 700    delete _charge_1;
 701    delete _iEtaPhi_1;
 702    delete _iPt_1;
 703    delete _pt_1;
 704    delete _px_1;
 705    delete _py_1;
 706    delete _pz_1;
 707    delete _correction_1;
 708    delete _dedx_1;
 709    delete __n1_1_vsPt;
 710    delete __n1_1_vsEtaPhi;
 711    delete __s1pt_1_vsEtaPhi;
 712    delete __n1_1_vsZEtaPhiPt;
 713    if (_correctionWeight_1) delete _correctionWeight_1;
 714
 715    if (!_sameFilter)
 716    {
 717    delete _id_2;
 718    delete _charge_2;
 719    delete _iEtaPhi_2;
 720    delete _iPt_2;
 721    delete _pt_2;
 722    delete _px_2;
 723    delete _py_2;
 724    delete _pz_2;
 725    delete _correction_2;
 726    delete _dedx_2;
 727    delete __n1_2_vsPt;
 728    delete __n1_2_vsEtaPhi;
 729    delete __s1pt_2_vsEtaPhi;
 730    delete __n1_2_vsZEtaPhiPt;
 731    if (_correctionWeight_2) delete _correctionWeight_2;
 732    }
 733
 734    if (!_singlesOnly)
 735    {
 736    delete __n2_12_vsPtPt;
 737    delete __n2_12_vsEtaPhi;
 738    delete __s2ptpt_12_vsEtaPhi;
 739    delete __s2PtN_12_vsEtaPhi;
 740    delete __s2NPt_12_vsEtaPhi;
 741    }
 742    */
 743 }
 744
 745 void AliDptDptInMC::UserCreateOutputObjects()
 746 {
 747   //cout<< "AliDptDptInMC::CreateOutputObjects() Starting " << endl;
 748   OpenFile(0);
 749   _outputHistoList = new TList();
 750   _outputHistoList->SetOwner();
 751
 752   //if (_useWeights) DefineInput(2, TList::Class());
 753
 754   //Setup the parameters of the histograms
 755   _nBins_M0 = 500; _min_M0   = 0.;    _max_M0    = 5000.;  _width_M0 = (_max_M0-_min_M0)/_nBins_M0;
 756   _nBins_M1 = 500; _min_M1   = 0.;    _max_M1    = 5000.;  _width_M1 = (_max_M1-_min_M1)/_nBins_M1;
 757   _nBins_M2 = 500; _min_M2   = 0.;    _max_M2    = 5000.;  _width_M2 = (_max_M2-_min_M2)/_nBins_M2;
 758   _nBins_M3 = 500; _min_M3   = 0.;    _max_M3    = 5000.;  _width_M3 = (_max_M3-_min_M3)/_nBins_M3;
 759   _nBins_M4 = 100; _min_M4   = 0.;    _max_M4    = 100.;   _width_M4 = (_max_M4-_min_M4)/_nBins_M4;
 760   _nBins_M5 = 100; _min_M5   = 0.;    _max_M5    = 100.;   _width_M5 = (_max_M5-_min_M5)/_nBins_M5;
 761   _nBins_M6 = 100; _min_M6   = 0.;    _max_M6    = 100.;   _width_M6 = (_max_M6-_min_M6)/_nBins_M6;
 762
 763   _min_vertexZ       = _vertexZMin;
 764   _max_vertexZ       = _vertexZMax;
 765   _width_vertexZ     = 0.5;
 766   _nBins_vertexZ     = int(0.5+ (_max_vertexZ - _min_vertexZ)/_width_vertexZ);
 767   _nBins_pt_1        = int(0.5+ (_max_pt_1 -_min_pt_1 )/_width_pt_1);
 768   _nBins_eta_1       = int(0.5+ (_max_eta_1-_min_eta_1)/_width_eta_1);
 769   _width_phi_1       = (_max_phi_1  - _min_phi_1)  /_nBins_phi_1;
 770   _nBins_etaPhi_1    = _nBins_phi_1    * _nBins_eta_1;
 771   _nBins_etaPhiPt_1  = _nBins_etaPhi_1 * _nBins_pt_1;
 772   _nBins_zEtaPhiPt_1 = _nBins_vertexZ  * _nBins_etaPhiPt_1;
 773
 774   _nBins_pt_2        =  int(0.5+ (_max_pt_2 -_min_pt_2 )/_width_pt_2);
 775   _nBins_eta_2       = int(0.5+ (_max_eta_2-_min_eta_2)/_width_eta_2);
 776   _width_phi_2       = (_max_phi_2  - _min_phi_2)  /_nBins_phi_2;
 777   _nBins_etaPhi_2    = _nBins_phi_2    * _nBins_eta_2;
 778   _nBins_etaPhiPt_2  = _nBins_etaPhi_2 * _nBins_pt_2;
 779   _nBins_zEtaPhiPt_2 = _nBins_vertexZ  * _nBins_etaPhiPt_2;
 780   _nBins_etaPhi_12   = _nBins_etaPhi_1 * _nBins_etaPhi_2;
 781
 782   //setup the work arrays
 783
 784   _id_1       = new int[arraySize];
 785   _charge_1   = new int[arraySize];
 786   //_iPhi_1 = new int[arraySize];
 787   //_iEta_1 = new int[arraySize];
 788   _iEtaPhi_1  = new int[arraySize];
 789   _iPt_1      = new int[arraySize];
 790   _pt_1       = new float[arraySize];
 791   _px_1       = new float[arraySize];
 792   _py_1       = new float[arraySize];
 793   _pz_1       = new float[arraySize];
 794   //_phi_1 = new float[arraySize];
 795   //_eta_1 = new float[arraySize];
 796   _correction_1 = new float[arraySize];
 797   //_dedx_1     = new float[arraySize];
 798
 799   __n1_1_vsPt              = getDoubleArray(_nBins_pt_1,        0.);
 800
 801   __n1_1_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_1,    0.);
 802   __s1pt_1_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_1,    0.);
 803   __n1_1_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_1,  0.);
 804
 805   //cout << "==========================================================================================" << endl;
 806   //cout << "=============== Booking for particle 1 done." << endl;
 807   //cout << "_requestedCharge_1: " << _requestedCharge_1 << endl;
 808   //cout << "_requestedCharge_2: " << _requestedCharge_2 << endl;
 809
 810   if (_requestedCharge_2!=_requestedCharge_1)
 811     {
 812       //cout << " creating arrays for particle 2 with size: " << arraySize << endl;
 813     _sameFilter = 0;
 814     //particle 2
 815     _id_2       = new int[arraySize];
 816     _charge_2   = new int[arraySize];
 817     //_iPhi_2   = new int[arraySize];
 818     //_iEta_2   = new int[arraySize];
 819     _iEtaPhi_2  = new int[arraySize];
 820     _iPt_2      = new int[arraySize];
 821     _pt_2       = new float[arraySize];
 822     _px_2       = new float[arraySize];
 823     _py_2       = new float[arraySize];
 824     _pz_2       = new float[arraySize];
 825     //_phi_2    = new float[arraySize];
 826     //_eta_2    = new float[arraySize];
 827     _correction_2 = new float[arraySize];
 828     //_dedx_2       = new float[arraySize];
 829
 830     __n1_2_vsPt              = getDoubleArray(_nBins_pt_2,        0.);
 831     __n1_2_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_2,    0.);
 832     __s1pt_2_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_2,    0.);
 833     __n1_2_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_2, 0.);
 834
 835     }
 836
 837   __n2_12_vsPtPt           = getDoubleArray(_nBins_pt_1*_nBins_pt_2,0.);
 838   __n2_12_vsEtaPhi         = getFloatArray(_nBins_etaPhi_12,       0.);
 839   __s2ptpt_12_vsEtaPhi     = getFloatArray(_nBins_etaPhi_12,       0.);
 840   __s2PtN_12_vsEtaPhi      = getFloatArray(_nBins_etaPhi_12,       0.);
 841   __s2NPt_12_vsEtaPhi      = getFloatArray(_nBins_etaPhi_12,       0.);
 842
 843   // Setup all the labels needed.
 844
 845   part_1_Name   = "_1";
 846   part_2_Name   = "_2";
 847   pair_12_Name  = "_12";
 848
 849   n1Name     = "n1";
 850   n2Name     = "n2";
 851   n1NwName   = "n1Nw";
 852   n2NwName   = "n2Nw";
 853   r1Name     = "r1";
 854   r2Name     = "r2";
 855   r3Name     = "r3";
 856   r2r1Name   = "r2r1";
 857   c2Name     = "c2";
 858   c3Name     = "c3";
 859   d3Name     = "d3";
 860   p3Name     = "p3";
 861   cName      = "sean";
 862
 863   intR2Name       = "intR2";
 864   binCorrName     = "binCorr";
 865   intBinCorrName  = "intBinCorr";
 866
 867   avgName      = "avg";
 868   sumName      = "sum";
 869   s1ptName     = "sumPt";
 870   s1ptNwName   = "sumPtNw";
 871   s1DptName    = "sumDpt";
 872   s2PtPtName   = "sumPtPt";
 873   s2PtPtNwName = "sumPtPtNw";
 874   s2DptDptName = "sumDptDpt";
 875   s2NPtName    = "sumNPt";
 876   s2NPtNwName  = "sumNPtNw";
 877   s2PtNName    = "sumPtN";
 878   s2NPtNwName  = "sumNPtNw";
 879   s2PtNNwName  = "sumPtNNw";
 880   ptName       = "avgPt";
 881   ptptName     = "avgPtPt";
 882   pt1pt1Name   = "avgPtavgPt";
 883   DptName      = "avgDpt";
 884   DptDptName   = "avgDptDpt";
 885   RDptDptName  = "relDptDpt"; // ratio of avgDptDpt by avgPt*avgPt
 886   nPtName      = "avgNpt";
 887   ptNName      = "avgPtN";
 888   seanName     = "seanC";
 889
 890   _title_counts = "yield";
 891
 892   _title_m0     = "M_{0}";
 893   _title_m1     = "M_{1}";
 894   _title_m2     = "M_{2}";
 895   _title_m3     = "M_{3}";
 896   _title_m4     = "V0Centrality";
 897   _title_m5     = "TrkCentrality";
 898   _title_m6     = "SpdCentrality";
 899
 900   _title_eta_1       = "#eta_{1}";
 901   _title_phi_1       = "#varphi_{1} (radian)";
 902   _title_etaPhi_1    = "#eta_{1}#times#varphi_{1}";
 903   _title_pt_1        = "p_{t,1} (GeV/c)";
 904   _title_n_1         = "n_{1}";
 905   _title_SumPt_1     = "#Sigma p_{t,1} (GeV/c)";
 906   _title_AvgPt_1     = "#LT p_{t,1} #GT (GeV/c)";
 907   _title_AvgN_1      = "#LT n_{1} #GT";
 908   _title_AvgSumPt_1  = "#LT #Sigma p_{t,1} #GT (GeV/c)";
 909
 910   _title_eta_2       = "#eta_{2}";
 911   _title_phi_2       = "#varphi_{2} (radian)";
 912   _title_etaPhi_2    = "#eta_{2}#times#varphi_{2}";
 913   _title_pt_2        = "p_{t,2} (GeV/c)";
 914   _title_n_2         = "n_{2}";
 915   _title_SumPt_2     = "#Sigma p_{t,1} (GeV/c)";
 916   _title_AvgPt_2     = "#LT p_{t,2} #GT (GeV/c)";
 917   _title_AvgN_2      = "#LT n_{2} #GT";
 918   _title_AvgSumPt_2  = "#LT #Sigma p_{t,2} #GT (GeV/c)";
 919
 920   _title_etaPhi_12   = "#eta_{1}#times#varphi_{1}#times#eta_{2}#times#varphi_{2}";
 921
 922   _title_AvgN2_12       = "#LT n_{2} #GT";;
 923   _title_AvgSumPtPt_12  = "#LT #Sigma p_{t,1}p_{t,2} #GT";;
 924   _title_AvgSumPtN_12   = "#LT #Sigma p_{t,1}N #GT";;
 925   _title_AvgNSumPt_12   = "#LT N#Sigma p_{t,2} #GT";;
 926
 927
 928   vsZ         = "_vsZ";
 929   vsM         = "_vsM";
 930   vsPt         = "_vsPt";
 931   vsPhi        = "_vsPhi";
 932   vsEta        = "_vsEta";
 933   vsEtaPhi     = "_vsEtaPhi";
 934   vsPtVsPt     = "_vsPtVsPt";
 935
 936
 937   if (_useWeights)
 938     {
 939     int iZ, iEtaPhi, iPt;
 940     int iZ1,iEtaPhi1,iPt1;
 941     int a, b;
 942     if (_weight_1)
 943       {
 944       _correctionWeight_1 = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
 945       a = _nBins_pt_1;
 946       b = _nBins_etaPhi_1*_nBins_pt_1;
 947       for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 948         {
 949         for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_1; iEtaPhi++, iEtaPhi1++)
 950           {
 951           for (iPt=0,iPt1=1; iPt<_nBins_pt_1; iPt++, iPt1++)
 952             {
 953             _correctionWeight_1[iZ*b+iEtaPhi*a+iPt] = _weight_1->GetBinContent(iZ1,iEtaPhi1,iPt1);
 954             }
 955           }
 956         }
 957       } // _weight_1
 958     else
 959       {
 960       AliError("AliDptDptInMC:: _weight_1 is a null pointer.");
 961       return;
 962       }
 963     if (!_sameFilter)
 964       {
 965       if (_weight_2)
 966         {
 967         _correctionWeight_2 = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
 968         a = _nBins_pt_2;
 969         b = _nBins_etaPhi_2*_nBins_pt_2;
 970         for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 971           {
 972           for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_2; iEtaPhi++, iEtaPhi1++)
 973             {
 974             for (iPt=0,iPt1=1; iPt<_nBins_pt_2; iPt++, iPt1++)
 975               {
 976               _correctionWeight_2[iZ*b+iEtaPhi*a+iPt] = _weight_2->GetBinContent(iZ1,iEtaPhi1,iPt1);
 977               }
 978             }
 979           }
 980         } // _weight_2
 981       else
 982         {
 983         AliError("AliDptDptInMC:: _weight_1 is a null pointer.");
 984         return;
 985         }
 986       }
 987     }
 988
 989   createHistograms();
 990   PostData(0,_outputHistoList);
 991
 992   //cout<< "AliDptDptInMC::CreateOutputObjects() DONE " << endl;
 993
 994 }
 995
 996 void  AliDptDptInMC::createHistograms()
 997 {
 998   AliInfo(" AliDptDptInMC::createHistoHistograms() Creating Event Histos");
 999   TString name;
1000
1001   name = "eventAccounting";
1002
1003   // bin index : what it is...
1004   //         0 :  number of event submitted
1005   //         1 :  number accepted by centrality cut
1006   //         2 :  number accepted by centrality cut and z cut
1007   //         3 :  total number of particles that satisfied filter 1
1008   //         4 :  total number of particles that satisfied filter 2
1009   _eventAccounting      = createHisto1D(name,name,10, -0.5, 9.5, "event Code", _title_counts);
1010
1011   name = "m0"; _m0      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m0, _title_counts);
1012   name = "m1"; _m1      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m1, _title_counts);
1013   name = "m2"; _m2      = createHisto1D(name,name,_nBins_M2, _min_M2, _max_M2, _title_m2, _title_counts);
1014   name = "m3"; _m3      = createHisto1D(name,name,_nBins_M3, _min_M3, _max_M3, _title_m3, _title_counts);
1015   name = "m4"; _m4      = createHisto1D(name,name,_nBins_M4, _min_M4, _max_M4, _title_m4, _title_counts);
1016   name = "m5"; _m5      = createHisto1D(name,name,_nBins_M5, _min_M5, _max_M5, _title_m5, _title_counts);
1017   name = "m6"; _m6      = createHisto1D(name,name,_nBins_M6, _min_M6, _max_M6, _title_m6, _title_counts);
1018   name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
1019
1020
1021   name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
1022   name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
1023   name = "DCAz";    _dcaz     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
1024   name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
1025
1026   name = "Eta2";     _etadis2   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
1027   name = "Phi2";     _phidis2   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
1028   name = "DCAz2";    _dcaz2     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
1029   name = "DCAxy2";   _dcaxy2    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
1030
1031   name = "Eta3";     _etadis3   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
1032   name = "Phi3";     _phidis3   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
1033   name = "DCAz3";    _dcaz3     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
1034   name = "DCAxy3";   _dcaxy3    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
1035
1036   //name = "Eta";     _etadis   = createHisto1F(name,name, 250, 0.0, 2.5, "#eta","counts"); //temporaryly
1037   //name = "Phi";     _phidis   = createHisto1F(name,name, 250, 0.0, 2.5, "#phi","counts");
1038   //name = "DCAz";    _dcaz     = createHisto1F(name,name, 250, 0.0, 2.5, "dcaZ","counts");
1039
1040   //name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
1041   //name = "Nclus2";   _Ncluster2    = createHisto1F(name,name, 200, 0, 200, "Ncluster2","counts");
1042
1043   if (_singlesOnly)
1044     {
1045     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);
1046
1047     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);
1048     //name = "dedxVsP_1";                          _dedxVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1049     //name = "corrDedxVsP_1";                      _corrDedxVsP_1          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1050     //name = "betaVsP_1";                          _betaVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
1051
1052     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);
1053     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);
1054     //name = "dedxVsP_2";                          _dedxVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1055     //name = "corrDedxVsP_2";                      _corrDedxVsP_2          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1056     //name = "betaVsP_2";                          _betaVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
1057
1058     }
1059   else
1060     {
1061     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);
1062     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);
1063     name = n1Name+part_1_Name+vsM;            _n1_1_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1064     name = s1ptName+part_1_Name+vsM;          _s1pt_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1065     name = n1NwName+part_1_Name+vsM;          _n1Nw_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1066     name = s1ptNwName+part_1_Name+vsM;        _s1ptNw_1_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1067
1068     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);
1069     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);
1070     name = n1Name+part_2_Name + vsM;          _n1_2_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_2);
1071     name = s1ptName+part_2_Name + vsM;        _s1pt_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_2);
1072     name = n1NwName+part_2_Name+vsM;          _n1Nw_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1073     name = s1ptNwName+part_2_Name+vsM;        _s1ptNw_2_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1074
1075     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);
1076     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);
1077     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);
1078     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);
1079     name = n2Name+pair_12_Name+vsPtVsPt;      _n2_12_vsPtVsPt       = createHisto2F(name,name, _nBins_pt_1, _min_pt_1, _max_pt_1, _nBins_pt_2, _min_pt_2, _max_pt_2, _title_pt_1, _title_pt_2, _title_AvgN2_12);
1080
1081     name = n2Name+pair_12_Name + vsM;         _n2_12_vsM            = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN2_12);
1082     name = s2PtPtName+pair_12_Name + vsM;     _s2PtPt_12_vsM        = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtPt_12);
1083     name = s2PtNName+pair_12_Name + vsM;      _s2PtN_12_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);
1084     name = s2NPtName+pair_12_Name + vsM;      _s2NPt_12_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);
1085
1086     name = n2NwName+pair_12_Name + vsM;       _n2Nw_12_vsM          = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN2_12);
1087     name = s2PtPtNwName+pair_12_Name + vsM;   _s2PtPtNw_12_vsM      = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtPt_12);
1088     name = s2PtNNwName+pair_12_Name + vsM;    _s2PtNNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);
1089     name = s2NPtNwName+pair_12_Name + vsM;    _s2NPtNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);
1090
1091     name = "mInv";     _invMass     = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts");
1092     name = "mInvElec"; _invMassElec = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts");
1093     }
1094
1095   AliInfo(" AliDptDptInMC::createHistoHistograms() All Done");
1096 }
1097 //-----------------------//
1098
1099 void  AliDptDptInMC::finalizeHistograms()
1100 {
1101
1102   AliInfo("AliDptDptInMC::finalizeHistograms() starting");
1103   AliInfo(Form("CorrelationAnalyzers::finalizeHistograms()   _eventCount : %d",int(_eventCount)));
1104   if (_singlesOnly)
1105     {
1106     if (_sameFilter)
1107       {
1108       fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1109       fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1110       fillHistoWithArray(_n1_2_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1111       fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1112       }
1113     else
1114       {
1115       fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1116       fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1117       fillHistoWithArray(_n1_2_vsPt,              __n1_2_vsPt,        _nBins_pt_2);
1118       fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_2_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_2, _nBins_pt_2);
1119       }
1120     }
1121   else
1122     {
1123     if (_sameFilter)
1124       {
1125       fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1126       fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1127       fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1128       fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1129       }
1130     else
1131       {
1132       fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1133       fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1134       fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_2_vsEtaPhi,    _nBins_eta_2,   _nBins_phi_2);
1135       fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_2_vsEtaPhi,  _nBins_eta_2,   _nBins_phi_2);
1136       }
1137     fillHistoWithArray(_n2_12_vsEtaPhi,     __n2_12_vsEtaPhi,     _nBins_etaPhi_12);
1138     fillHistoWithArray(_s2PtPt_12_vsEtaPhi, __s2ptpt_12_vsEtaPhi, _nBins_etaPhi_12);
1139     fillHistoWithArray(_s2PtN_12_vsEtaPhi,  __s2PtN_12_vsEtaPhi,  _nBins_etaPhi_12);
1140     fillHistoWithArray(_s2NPt_12_vsEtaPhi,  __s2NPt_12_vsEtaPhi,  _nBins_etaPhi_12);
1141     fillHistoWithArray(_n2_12_vsPtVsPt,     __n2_12_vsPtPt,       _nBins_pt_1,    _nBins_pt_2);
1142
1143     }
1144   AliInfo("AliDptDptInMC::finalizeHistograms()  Done ");
1145 }
1146 //--------------//
1147
1148
1149 void  AliDptDptInMC::UserExec(Option_t */*option*/)
1150 {
1151
1152
1153   int    k1,k2;
1154   int    iPhi, iEta, iEtaPhi, iPt, charge;
1155   float  q, phi, pt, eta, corr, corrPt, px, py, pz;
1156   int    ij;
1157   int    id_1, q_1, iEtaPhi_1, iPt_1;
1158   float  pt_1, px_1, py_1, pz_1, corr_1;
1159   int    id_2, q_2, iEtaPhi_2, iPt_2;
1160   float  pt_2, px_2, py_2, pz_2, corr_2;
1161   float  ptpt;
1162   int    iVertex, iVertexP1, iVertexP2;
1163   int    iZEtaPhiPt;
1164   float  massElecSq = 2.5e-7;
1165   const  AliAODVertex*  vertex;
1166   //int    nClus;
1167   bool   bitOK;
1168
1169   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
1170   if (!manager) {
1171     //cout<<"ERROR: Analysis manager not found."<<endl;
1172     return;
1173   }
1174   //coneect to the inputHandler------------
1175   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
1176   if (!inputHandler) {
1177     return;
1178   }
1179
1180   fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
1181   if (!fAODEvent)
1182     {
1183       return;
1184     }
1185   fPIDResponse =inputHandler->GetPIDResponse();
1186   if (!fPIDResponse){
1187     AliFatal("This Task needs the PID response attached to the inputHandler");
1188     return;
1189   }
1190   // count all events looked at here
1191   _eventCount++;
1192
1193   if (_eventAccounting)
1194     {
1195       _eventAccounting->Fill(0);// count all calls to this function
1196     }
1197   else
1198     {
1199       return;
1200     }
1201
1202   _eventAccounting->Fill(1);// count all calls to this function with a valid pointer
1203   //reset single particle counters
1204   k1 = k2 = 0;
1205   __n1_1 = __n1_2 = __s1pt_1 = __s1pt_2 = __n1Nw_1 = __n1Nw_2 = __s1ptNw_1 = __s1ptNw_2 = 0;
1206
1207   float v0Centr  = -999.;
1208   float v0ACentr  = -999.;
1209   float trkCentr = -999.;
1210   float spdCentr = -999.;
1211   //float vertexX  = -999;
1212   //float vertexY  = -999;
1213   float vertexZ  = -999;
1214   //float vertexXY = -999;
1215
1216   float centrality = -999;
1217   //Double_t nSigma =-999;
1218   if(fAODEvent)
1219     {
1220
1221       //Centrality
1222       //AliAODHeader* centralityObject = ((AliVAODHeader*)fAODEvent->GetHeader())->GetCentralityP();
1223       AliCentrality* centralityObject =  ((AliVAODHeader*)fAODEvent->GetHeader())->GetCentralityP();
1224       if (centralityObject)
1225         {
1226           v0Centr  = centralityObject->GetCentralityPercentile("V0M");
1227           v0ACentr  = centralityObject->GetCentralityPercentile("V0A");
1228           trkCentr = centralityObject->GetCentralityPercentile("TRK");
1229           spdCentr = centralityObject->GetCentralityPercentile("CL1");
1230
1231         }
1232
1233       _mult4    = v0Centr;
1234       _mult4a   = v0ACentr;
1235       _mult5    = trkCentr;
1236       _mult6    = spdCentr;
1237       _field    = fAODEvent->GetMagneticField();
1238       //_centralityMethod
1239       switch (_centralityMethod)
1240         {
1241         case 0: centrality = _mult0; break;
1242         case 1: centrality = _mult1; break;
1243         case 2: centrality = _mult2; break;
1244         case 3: centrality = _mult3; break;
1245         case 4: centrality = _mult4; break;
1246         case 5: centrality = _mult5; break;
1247         case 6: centrality = _mult6; break;
1248         case 7: centrality = _mult4a; break;
1249         }
1250
1251       if ( centrality < _centralityMin ||  centrality > _centralityMax ) return;
1252
1253       _eventAccounting->Fill(2);// count all events with right centrality
1254
1255       vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertex();
1256
1257       if(vertex->GetNContributors() > 0)
1258         {
1259           //vertexX = vertex->GetX();
1260           //vertexY = vertex->GetY();
1261           vertexZ = vertex->GetZ();
1262           //vertexXY = sqrt(vertexX*vertexX+vertexY*vertexY);
1263         }
1264       if (vertexZ  < _vertexZMin  || vertexZ  > _vertexZMax ) return;
1265
1266       iVertex = int((vertexZ-_min_vertexZ)/_width_vertexZ);
1267       iVertexP1 = iVertex*_nBins_etaPhiPt_1;
1268       iVertexP2 = iVertex*_nBins_etaPhiPt_2;
1269       if (iVertex<0 || iVertex>=_nBins_vertexZ)
1270         {
1271           AliError("AliTaskDptCorrMC::Exec(Option_t *option) iVertex<0 || iVertex>=_nBins_vertexZ ");
1272           return;
1273         }
1274       _eventAccounting->Fill(3);// count all calls to this function with a valid pointer
1275       //======================
1276       //***********************************
1277       //MC AOD Truth
1278       if(fAnalysisType == "MCAOD")
1279         {
1280           AliMCEvent* mcEvent = MCEvent();
1281           _nTracks = mcEvent->GetNumberOfTracks();
1282           _mult3    = _nTracks;
1283           //loop over tracks starts here
1284           for (int iTrack=0; iTrack< _nTracks; iTrack++)
1285             {
1286               AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTrack);
1287
1288               if (!aodTrack)
1289                 {
1290                   AliError(Form("AliTaskDptCorrMC::Exec(Option_t *option) No track ofr iTrack=%d", iTrack));
1291                   continue;
1292                 }
1293
1294               if(!aodTrack->IsPhysicalPrimary()) continue;
1295
1296               q      = aodTrack->Charge();
1297               charge = int(q);
1298               phi    = aodTrack->Phi();
1299               pt     = aodTrack->Pt();
1300               px     = aodTrack->Px();
1301               py     = aodTrack->Py();
1302               pz     = aodTrack->Pz();
1303               eta    = aodTrack->Eta();
1304               // Kinematics cuts from ESD track cuts
1305               if( pt < 0.2 || pt > 2.0)      continue;
1306               if( eta < -0.8 || eta > 0.8)  continue;
1307
1308               if(q == 0) continue;
1309               _dcaz->Fill(pt); //AllCh
1310
1311               if(TMath::Abs(aodTrack->GetPdgCode()) == 211)
1312                 {
1313                   _etadis->Fill(pt); //pion
1314                 }
1315               if(TMath::Abs(aodTrack->GetPdgCode()) == 321)
1316                 {
1317                   _phidis->Fill(pt); //kaon
1318                 }
1319
1320               //_etadis->Fill(eta);
1321               //_phidis->Fill(phi);
1322
1323               if(fExcludeResonancesInMC)
1324                 {
1325                   //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
1326                   Int_t gMotherIndex = aodTrack->GetMother();
1327                   if(gMotherIndex != -1) {
1328                     AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1329                     if(motherTrack) {
1330                       Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1331
1332                       if(pdgCodeOfMother == 311  ||
1333                          pdgCodeOfMother == -311 ||
1334                          pdgCodeOfMother == 310  ||
1335                          pdgCodeOfMother == 3122 ||
1336                          pdgCodeOfMother == -3122 ||
1337                          pdgCodeOfMother == 111 ||
1338                          pdgCodeOfMother == 22 ) continue;
1339                     }
1340                   }
1341                 }
1342               //Exclude electrons with PDG
1343
1344               if(fExcludeElectronsInMC) {
1345                 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1346               }
1347
1348               //Particle loop 1st particle
1349               if (q > 0)
1350                 {
1351
1352                   iPhi   = int( phi/_width_phi_1);
1353
1354                   if (iPhi<0 || iPhi>=_nBins_phi_1 )
1355                     {
1356                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1357                       return;
1358                     }
1359
1360                   iEta    = int((eta-_min_eta_1)/_width_eta_1);
1361                   if (iEta<0 || iEta>=_nBins_eta_1)
1362                     {
1363                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1364                       continue;
1365                     }
1366
1367                   iPt     = int((pt -_min_pt_1 )/_width_pt_1 );
1368                   if (iPt<0  || iPt >=_nBins_pt_1)
1369                     {
1370                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1371                       continue;
1372                     }
1373                   iEtaPhi = iEta*_nBins_phi_1+iPhi;
1374                   iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
1375
1376                   if (_correctionWeight_1)
1377                     corr = _correctionWeight_1[iZEtaPhiPt];
1378                   else
1379                     corr = 1;
1380                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1)
1381                     {
1382                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1");
1383                       continue;
1384                     }
1385
1386                   if (_singlesOnly)
1387                     {
1388
1389                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1390
1391                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1392
1393                     }
1394
1395                   else
1396                     {
1397                       corrPt                      = corr*pt;
1398                       _id_1[k1]                   = iTrack;
1399                       _charge_1[k1]               = charge;
1400                       _iEtaPhi_1[k1]              = iEtaPhi;
1401                       _iPt_1[k1]                  = iPt;
1402                       _pt_1[k1]                   = pt;
1403                       _px_1[k1]                   = px;
1404                       _py_1[k1]                   = py;
1405                       _pz_1[k1]                   = pz;
1406                       _correction_1[k1]           = corr;
1407                       __n1_1                      += corr;
1408                       __n1_1_vsEtaPhi[iEtaPhi]    += corr;
1409                       __s1pt_1                    += corrPt;
1410                       __s1pt_1_vsEtaPhi[iEtaPhi]  += corrPt;
1411                       __n1Nw_1                    += 1;
1412                       __s1ptNw_1                  += pt;
1413                       ++k1;
1414                       if (k1>=arraySize)
1415                         {
1416                           AliError(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
1417                           return;
1418                         }
1419                     }
1420                 }
1421
1422               if (q < 0)
1423                 {
1424                   iPhi   = int( phi/_width_phi_2);
1425
1426                   if (iPhi<0 || iPhi>=_nBins_phi_2 )
1427                     {
1428                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1429                       return;
1430                     }
1431
1432                   iEta    = int((eta-_min_eta_2)/_width_eta_2);
1433                   if (iEta<0 || iEta>=_nBins_eta_2)
1434                     {
1435                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1436                       continue;
1437                     }
1438
1439                   iPt     = int((pt -_min_pt_2 )/_width_pt_2 );
1440                   if (iPt<0  || iPt >=_nBins_pt_2)
1441                     {
1442                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1443                       continue;
1444                     }
1445
1446                   iEtaPhi = iEta*_nBins_phi_2+iPhi;
1447                   iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
1448                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
1449                     {
1450                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
1451                       continue;
1452                     }
1453
1454                   if (_correctionWeight_2)
1455                     corr = _correctionWeight_2[iZEtaPhiPt];
1456                   else
1457                     corr = 1;
1458                   //dpt = pt - (charge>0) ? _avgPt_vsEtaPhi_2p[iEtaPhi] : _avgPt_vsEtaPhi_2m[iEtaPhi];
1459
1460                   if (_singlesOnly)
1461                     {
1462                       __n1_2_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1463                       __n1_2_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1464                     }
1465                   else
1466                     {
1467                       corrPt                      = corr*pt;
1468                       _id_2[k2]                   = iTrack;         //cout << "step 1" << endl;
1469                       _charge_2[k2]               = charge;         //cout << "step 2" << endl;
1470                       _iEtaPhi_2[k2]              = iEtaPhi;        //cout << "step 3" << endl;
1471                       _iPt_2[k2]                  = iPt;            //cout << "step 4" << endl;
1472                       _pt_2[k2]                   = pt;             //cout << "step 5" << endl;
1473                       _px_2[k2]                   = px;             //cout << "step 6" << endl;
1474                       _py_2[k2]                   = py;             //cout << "step 7" << endl;
1475                       _pz_2[k2]                   = pz;             //cout << "step 8" << endl;
1476                       _correction_2[k2]           = corr;           //cout << "step 9" << endl;
1477                       __n1_2                      += corr;          //cout << "step 10" << endl;
1478                       __s1pt_2                    += corrPt;        //cout << "step 13" << endl;
1479                       __n1Nw_2                    += 1;
1480                       __n1_2_vsEtaPhi[iEtaPhi]    += corr;          //cout << "step 11" << endl;
1481                       __s1pt_2_vsEtaPhi[iEtaPhi]  += corrPt;        //cout << "step 14" << endl;
1482                       __s1ptNw_2                  += pt;
1483                       ++k2;
1484                       if (k2>=arraySize)
1485                         {
1486                           AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize));
1487                           return;
1488                         }
1489                     }
1490
1491                   //cout << "done with track" << endl;
1492                 } //iTrack
1493             } //data aod loop
1494           } //MC AOD loop ends
1495
1496       //***********************************************
1497       //MC AOD Reconstructed tracks
1498       if(fAnalysisType == "MCAODreco")
1499         {
1500
1501           fArrayMC = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1502           if (!fArrayMC) {
1503             AliError("No array of MC particles found !!!");
1504           }
1505
1506           AliMCEvent* mcEvent = MCEvent();
1507           if (!mcEvent) {
1508             AliError("ERROR: Could not retrieve MC event");
1509           }
1510
1511           TExMap *trackMap = new TExMap();//Mapping matrix----
1512           //1st loop track for Global tracks
1513           for(Int_t i = 0; i < fAODEvent->GetNumberOfTracks(); i++)
1514             {
1515               AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
1516               if(!aodTrack) {
1517                 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
1518                 continue;
1519               }
1520               Int_t gID = aodTrack->GetID();
1521               if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks
1522             }
1523
1524           AliAODTrack* newAodTrack;
1525
1526           for (int iTrack=0; iTrack< fAODEvent->GetNumberOfTracks(); iTrack++)
1527             {
1528
1529               AliAODTrack *t = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(iTrack));
1530
1531               if(!t) {
1532                 AliError(Form("ERROR: Could not retrieve AODtrack %d",iTrack));
1533                 continue;
1534               }
1535
1536               bitOK  = t->TestFilterBit(_trackFilterBit);
1537               if (!bitOK) continue;
1538               Int_t gID = t->GetID();
1539               newAodTrack = gID >= 0 ?t : dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(trackMap->GetValue(-1-gID)));
1540               if(!newAodTrack) {
1541                 AliFatal("Not a standard AOD?");
1542               }
1543
1544               q      = t->Charge();
1545               charge = int(q);
1546               phi    = t->Phi();
1547               pt     = t->Pt();
1548               px     = t->Px();
1549               py     = t->Py();
1550               pz     = t->Pz();
1551               eta    = t->Eta();
1552               //Float_t dcaXY = t->DCA();
1553               //Float_t dcaZ  = t->ZAtDCA();
1554
1555               // get the electron nsigma
1556               Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
1557               Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kPion));
1558               Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kKaon));
1559               Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kProton));
1560
1561               //Make the decision based on the n-sigma of electrons exclusively
1562               if(nSigma < fNSigmaCut
1563                  && nSigmaPions   > fNSigmaCut
1564                  && nSigmaKaons   > fNSigmaCut
1565                  && nSigmaProtons > fNSigmaCut ) continue;
1566
1567
1568               if(q == 0) continue;
1569               // Kinematics cuts
1570               if( pt < 0.2 || pt > 2.0)      continue;
1571               if( eta < _min_eta_1 || eta > _max_eta_1)  continue;
1572
1573
1574               /*
1575                 _dcaz->Fill(pt);
1576                 if (nSigmaPions < fNSigmaCut)
1577                 {
1578                 _etadis->Fill(pt);
1579                 }
1580                 if (nSigmaKaons < fNSigmaCut)
1581                 {
1582                 _phidis->Fill(pt);
1583                 }
1584               */ //particle ratio calculation
1585
1586               Double_t pos[3];
1587               newAodTrack->GetXYZ(pos);
1588
1589               //Double_t DCAX = pos[0] - vertexX;
1590               //Double_t DCAY = pos[1] - vertexY;
1591               Double_t DCAZ = pos[2] - vertexZ;
1592
1593               //Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
1594
1595               /*if (DCAZ     <  _dcaZMin ||
1596                 DCAZ     >  _dcaZMax ||
1597                 DCAXY    >  _dcaXYMax ) continue;
1598               */
1599
1600               Int_t label = TMath::Abs(t->GetLabel());
1601               AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1602               //------------Test of tracks QAs------------
1603               if (AODmcTrack->IsPhysicalPrimary())
1604                 {
1605                   _dcaz->Fill(DCAZ);
1606                   //_dcaxy->Fill(DCAXY);
1607                   //_etadis->Fill(eta);
1608                   _phidis->Fill(phi);
1609                 }
1610               if (AODmcTrack->IsSecondaryFromWeakDecay())
1611                 {
1612                   _dcaz2->Fill(DCAZ);
1613                   //_dcaxy2->Fill(DCAXY);
1614                   //_etadis2->Fill(eta);
1615                   _phidis2->Fill(phi);
1616                 }
1617
1618               if (AODmcTrack->IsSecondaryFromMaterial())
1619                 {
1620                   _dcaz3->Fill(DCAZ);
1621                   //_dcaxy3->Fill(DCAXY);
1622                   //_etadis3->Fill(eta);
1623                   _phidis3->Fill(phi);
1624                 }
1625               //---------------------------------
1626
1627               /* //W/Wo Secondaries
1628               //if (!AODmcTrack->IsPhysicalPrimary()) continue;
1629               //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
1630               if(fExcludeResonancesInMC)
1631               {
1632               Int_t label = TMath::Abs(t->GetLabel());
1633               AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1634
1635               Int_t gMotherIndex = AODmcTrack->GetMother();
1636               if(gMotherIndex != -1) {
1637               AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1638               if(motherTrack) {
1639               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1640
1641               if(pdgCodeOfMother == 311  ||
1642               pdgCodeOfMother == -311 ||
1643               pdgCodeOfMother == 310  ||
1644               pdgCodeOfMother == 3122 ||
1645               pdgCodeOfMother == -3122 ||
1646               pdgCodeOfMother == 111 ||
1647               pdgCodeOfMother == 22 ) continue;
1648               }
1649               }
1650               } */
1651
1652
1653               //Int_t label = TMath::Abs(t->GetLabel());
1654               //AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1655               if (AODmcTrack)
1656                 {
1657                   if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1658                 }
1659
1660               //Particle 1
1661               if (t->Charge() > 0)
1662                 {
1663
1664                   iPhi   = int( phi/_width_phi_1);
1665
1666                   if (iPhi<0 || iPhi>=_nBins_phi_1 )
1667                     {
1668                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1669                       return;
1670                     }
1671
1672                   iEta    = int((eta-_min_eta_1)/_width_eta_1);
1673                   if (iEta<0 || iEta>=_nBins_eta_1)
1674                     {
1675                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1676                       continue;
1677                     }
1678                   iPt     = int((pt -_min_pt_1 )/_width_pt_1 );
1679                   if (iPt<0  || iPt >=_nBins_pt_1)
1680                     {
1681                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1682                       continue;
1683                     }
1684                   iEtaPhi = iEta*_nBins_phi_1+iPhi;
1685                   iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
1686
1687                   if (_correctionWeight_1)
1688                     corr = _correctionWeight_1[iZEtaPhiPt];
1689                   else
1690                     corr = 1;
1691                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1)
1692                     {
1693                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1");
1694                       continue;
1695                     }
1696
1697
1698                   if (_singlesOnly)
1699                     {
1700
1701                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1702
1703                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1704
1705                     }
1706
1707                   else
1708                     {
1709                       corrPt                      = corr*pt;
1710                       _id_1[k1]                   = iTrack;
1711                       _charge_1[k1]               = charge;
1712                       _iEtaPhi_1[k1]              = iEtaPhi;
1713                       _iPt_1[k1]                  = iPt;
1714                       _pt_1[k1]                   = pt;
1715                       _px_1[k1]                   = px;
1716                       _py_1[k1]                   = py;
1717                       _pz_1[k1]                   = pz;
1718                       _correction_1[k1]           = corr;
1719                       __n1_1                      += corr;
1720                       __n1_1_vsEtaPhi[iEtaPhi]    += corr;
1721                       __s1pt_1                    += corrPt;
1722                       __s1pt_1_vsEtaPhi[iEtaPhi]  += corrPt;
1723                       __n1Nw_1                    += 1;
1724                       __s1ptNw_1                  += pt;
1725                       ++k1;
1726                       if (k1>=arraySize)
1727                         {
1728                           AliError(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
1729                           return;
1730                         }
1731                     }
1732                 }
1733
1734               if (t->Charge() < 0 )
1735                 {
1736                   //===================================
1737
1738                   iPhi   = int( phi/_width_phi_2);
1739
1740                   if (iPhi<0 || iPhi>=_nBins_phi_2 )
1741                     {
1742                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1743                       return;
1744                     }
1745
1746                   iEta    = int((eta-_min_eta_2)/_width_eta_2);
1747                   if (iEta<0 || iEta>=_nBins_eta_2)
1748                     {
1749                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1750                       continue;
1751                     }
1752                   iPt     = int((pt -_min_pt_2 )/_width_pt_2 );
1753                   if (iPt<0  || iPt >=_nBins_pt_2)
1754                     {
1755                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1756                       continue;
1757                     }
1758
1759                   iEtaPhi = iEta*_nBins_phi_2+iPhi;
1760                   iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
1761                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
1762                     {
1763                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
1764                       continue;
1765                     }
1766
1767
1768                   if (_correctionWeight_2)
1769                     corr = _correctionWeight_2[iZEtaPhiPt];
1770                   else
1771                     corr = 1;
1772                   //dpt = pt - (charge>0) ? _avgPt_vsEtaPhi_2p[iEtaPhi] : _avgPt_vsEtaPhi_2m[iEtaPhi];
1773
1774                   if (_singlesOnly)
1775                     {
1776                       __n1_2_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1777                       __n1_2_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1778                     }
1779
1780                   else
1781                     {
1782                       corrPt                      = corr*pt;
1783                       _id_2[k2]                   = iTrack;         //cout << "step 1" << endl;
1784                       _charge_2[k2]               = charge;         //cout << "step 2" << endl;
1785                       _iEtaPhi_2[k2]              = iEtaPhi;        //cout << "step 3" << endl;
1786                       _iPt_2[k2]                  = iPt;            //cout << "step 4" << endl;
1787                       _pt_2[k2]                   = pt;             //cout << "step 5" << endl;
1788                       _px_2[k2]                   = px;             //cout << "step 6" << endl;
1789                       _py_2[k2]                   = py;             //cout << "step 7" << endl;
1790                       _pz_2[k2]                   = pz;             //cout << "step 8" << endl;
1791                       _correction_2[k2]           = corr;           //cout << "step 9" << endl;
1792                       __n1_2                      += corr;          //cout << "step 10" << endl;
1793                       __s1pt_2                    += corrPt;        //cout << "step 13" << endl;
1794                       __n1Nw_2                    += 1;
1795                       __n1_2_vsEtaPhi[iEtaPhi]    += corr;          //cout << "step 11" << endl;
1796                       __s1pt_2_vsEtaPhi[iEtaPhi]  += corrPt;        //cout << "step 14" << endl;
1797                       __s1ptNw_2                  += pt;
1798                       ++k2;
1799
1800                       if (k2>=arraySize)
1801                         {
1802                           AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize));
1803                           return;
1804                         }
1805                     }
1806
1807                   //cout << "done with track" << endl;
1808                 } //iTrack
1809             } //data aod loop
1810         } //MC AOD loop ends
1811
1812       //************************************************
1813     } //AOD events
1814
1815
1816   //  _m0->Fill(_mult0);
1817   //_m1->Fill(_mult1);
1818   //_m2->Fill(_mult2);
1819   //_m3->Fill(_mult3);
1820   //_m4->Fill(_mult4);
1821   //_m5->Fill(_mult5);
1822   //_m6->Fill(_mult6);
1823   _vertexZ->Fill(vertexZ);
1824
1825   if (_singlesOnly)
1826     {
1827       // nothing to do here.
1828     }
1829   else
1830     {
1831       if (_sameFilter)
1832         {
1833           _n1_1_vsM->Fill(centrality,      __n1_1);
1834           _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
1835           _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
1836           _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
1837           _n1_2_vsM->Fill(centrality,      __n1_1);
1838           _s1pt_2_vsM->Fill(centrality,    __s1pt_1);
1839           _n1Nw_2_vsM->Fill(centrality,    __n1Nw_1);
1840           _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_1);
1841           // reset pair counters
1842           __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
1843           __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
1844           if (_field>0)
1845             {
1846               for (int i1=0; i1<k1; i1++)
1847                 {
1848                   ////cout << "         i1:" << i1 << endl;
1849                   id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1850                   q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1851                   iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1852                   iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1853                   corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1854                   pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1855
1856
1857                   //1 and 2
1858                   for (int i2=i1+1; i2<k1; i2++)
1859                     {
1860                       ////cout << "         i2:" << i2 << endl;
1861                       id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
1862                       if (id_1!=id_2)
1863                         {
1864                           q_2       = _charge_1[i2];     ////cout << "        q_1:" << q_1 << endl;
1865                           iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1866                           iPt_2     = _iPt_1[i2];        ////cout << "      iPt_1:" << iPt_1 << endl;
1867                           corr_2    = _correction_1[i2]; ////cout << "     corr_1:" << corr_1 << endl;
1868                           pt_2      = _pt_1[i2];         ////cout << "       pt_1:" << pt_1 << endl;
1869
1870                           corr      = corr_1*corr_2;
1871                           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))
1872                             {
1873                               ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
1874                             }
1875                           else // swap particles
1876                             {
1877                               ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
1878                             }
1879
1880                           __n2_12                  += corr;
1881                           __n2_12_vsEtaPhi[ij]     += corr;
1882                           ptpt                     = pt_1*pt_2;
1883                           __s2ptpt_12              += corr*ptpt;
1884                           __s2PtN_12               += corr*pt_1;
1885                           __s2NPt_12               += corr*pt_2;
1886                           __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
1887                           __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
1888                           __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
1889                           __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
1890
1891                           __n2Nw_12                  += 1;
1892                           __s2ptptNw_12              += ptpt;
1893                           __s2PtNNw_12               += pt_1;
1894                           __s2NPtNw_12               += pt_2;
1895
1896                         }
1897                     } //i2
1898                 } //i1
1899             }
1900
1901           else // field<0
1902             {
1903               for (int i1=0; i1<k1; i1++)
1904                 {
1905                   ////cout << "         i1:" << i1 << endl;
1906                   id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1907                   q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1908                   iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1909                   iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1910                   corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1911                   pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1912
1913                   //1 and 2
1914                   for (int i2=i1+1; i2<k1; i2++)
1915                     {
1916                       ////cout << "         i2:" << i2 << endl;
1917                       id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
1918                       if (id_1!=id_2)
1919                         {
1920                           q_2       = _charge_1[i2];     ////cout << "        q_2:" << q_2 << endl;
1921                           iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
1922                           iPt_2     = _iPt_1[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
1923                           corr_2    = _correction_1[i2]; ////cout << "     corr_2:" << corr_2 << endl;
1924                           pt_2      = _pt_1[i2];         ////cout << "       pt_2:" << pt_2 << endl;
1925                           corr      = corr_1*corr_2;
1926                           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))
1927                             {
1928                               ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
1929                             }
1930                           else // swap particles
1931                             {
1932                               ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
1933                             }
1934
1935                           __n2_12                  += corr;
1936                           __n2_12_vsEtaPhi[ij]     += corr;
1937                           ptpt                     = pt_1*pt_2;
1938                           __s2ptpt_12              += corr*ptpt;
1939                           __s2PtN_12               += corr*pt_1;
1940                           __s2NPt_12               += corr*pt_2;
1941                           __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
1942                           __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
1943                           __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
1944                           __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
1945
1946                           __n2Nw_12                  += 1;
1947                           __s2ptptNw_12              += ptpt;
1948                           __s2PtNNw_12               += pt_1;
1949                           __s2NPtNw_12               += pt_2;
1950
1951                         }
1952                     } //i2
1953                 } //i1
1954             }
1955         }
1956       else  // filter 1 and 2 are different -- must do all particle pairs...
1957         {
1958           _n1_1_vsM->Fill(centrality,      __n1_1);
1959           _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
1960           _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
1961           _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
1962           _n1_2_vsM->Fill(centrality,      __n1_2);
1963           _s1pt_2_vsM->Fill(centrality,    __s1pt_2);
1964           _n1Nw_2_vsM->Fill(centrality,    __n1Nw_2);
1965           _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_2);
1966           // reset pair counters
1967           __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
1968           __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
1969           for (int i1=0; i1<k1; i1++)
1970             {
1971               ////cout << "         i1:" << i1 << endl;
1972               id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1973               q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1974               iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1975               iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1976               corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1977               pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1978               px_1      = _px_1[i1];          ////cout << "      px_1:" << px_1 << endl;
1979               py_1      = _py_1[i1];          ////cout << "      py_1:" << py_1 << endl;
1980               pz_1      = _pz_1[i1];          ////cout << "      pz_1:" << pz_1 << endl;
1981               //1 and 2
1982               for (int i2=0; i2<k2; i2++)
1983                 {
1984                   ////cout << "         i2:" << i2 << endl;
1985                   id_2   = _id_2[i2];              ////cout << "       id_2:" << id_2 << endl;
1986                   if (id_1!=id_2)  // exclude auto correlation
1987                     {
1988                       q_2       = _charge_2[i2];     ////cout << "        q_2:" << q_2 << endl;
1989                       iEtaPhi_2 = _iEtaPhi_2[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
1990                       iPt_2     = _iPt_2[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
1991                       corr_2    = _correction_2[i2]; ////cout << "     corr_2:" << corr_2 << endl;
1992                       pt_2      = _pt_2[i2];         ////cout << "       pt_2:" << pt_2 << endl;
1993                       px_2      = _px_2[i2];          ////cout << "      px_2:" << px_2 << endl;
1994                       py_2      = _py_2[i2];          ////cout << "      py_2:" << py_2 << endl;
1995                       pz_2      = _pz_2[i2];          ////cout << "      pz_2:" << pz_2 << endl;
1996
1997
1998                       if (_rejectPairConversion)
1999                         {
2000                           float e1Sq = massElecSq + pt_1*pt_1 + pz_1*pz_1;
2001                           float e2Sq = massElecSq + pt_2*pt_2 + pz_2*pz_2;
2002                           float mInvSq = 2*(massElecSq + sqrt(e1Sq*e2Sq) - px_1*px_2 - py_1*py_2 - pz_1*pz_2 );
2003                           float mInv = sqrt(mInvSq);
2004                           _invMass->Fill(mInv);
2005                         }
2006
2007                       corr      = corr_1*corr_2;
2008                       ij        = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
2009                       __n2_12                  += corr;
2010                       __n2_12_vsEtaPhi[ij]     += corr;
2011                       ptpt                     = pt_1*pt_2;
2012                       __s2ptpt_12              += corr*ptpt;
2013                       __s2PtN_12               += corr*pt_1;
2014                       __s2NPt_12               += corr*pt_2;
2015                       __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
2016                       __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
2017                       __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
2018                       __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
2019                       __n2Nw_12                  += 1;
2020                       __s2ptptNw_12              += ptpt;
2021                       __s2PtNNw_12               += pt_1;
2022                       __s2NPtNw_12               += pt_2;
2023
2024                     }
2025                 } //i2
2026             } //i1
2027         }
2028
2029       _n2_12_vsM->Fill(centrality,     __n2_12);
2030       _s2PtPt_12_vsM->Fill(centrality, __s2ptpt_12);
2031       _s2PtN_12_vsM->Fill(centrality,  __s2NPt_12);
2032       _s2NPt_12_vsM->Fill(centrality,  __s2PtN_12);
2033
2034       _n2Nw_12_vsM->Fill(centrality,     __n2Nw_12);
2035       _s2PtPtNw_12_vsM->Fill(centrality, __s2ptptNw_12);
2036       _s2PtNNw_12_vsM->Fill(centrality,  __s2NPtNw_12);
2037       _s2NPtNw_12_vsM->Fill(centrality,  __s2PtNNw_12);
2038
2039     }
2040
2041
2042   AliInfo("AliTaskDptCorrMC::UserExec()   -----------------Event Done ");
2043   PostData(0,_outputHistoList);
2044
2045 } //End of UserExec
2046
2047
2048
2049
2050 void   AliDptDptInMC::FinishTaskOutput()
2051 {
2052   AliInfo("AliDptDptInMC::FinishTaskOutput() Starting.");
2053   Printf("= 0 ====================================================================");
2054   finalizeHistograms();
2055   AliInfo("= 1 ====================================================================");
2056   PostData(0,_outputHistoList);
2057   AliInfo("= 2 ====================================================================");
2058   AliInfo("AliDptDptInMC::FinishTaskOutput() Done.");
2059 }
2060
2061 void   AliDptDptInMC::Terminate(Option_t* /*option*/)
2062 {
2063   AliInfo("AliDptDptInMC::Terminate() Starting/Done.");
2064 }
2065
2066
2067 //Tools
2068 //===================================================================================================
2069 void  AliDptDptInMC::fillHistoWithArray(TH1 * h, float * array, int size)
2070 {
2071   int i, i1;
2072   float v1, ev1, v2, ev2, sum, esum;
2073   for (i=0, i1=1; i<size; ++i,++i1)
2074     {
2075     v1  = array[i]; ev1 = sqrt(v1);
2076     v2  = h->GetBinContent(i1);
2077     ev2 = h->GetBinError(i1);
2078     sum = v1 + v2;
2079     esum = sqrt(ev1*ev1+ev2*ev2);
2080     h->SetBinContent(i1,sum);
2081     h->SetBinError(i1,esum);
2082     }
2083 }
2084
2085 void  AliDptDptInMC::fillHistoWithArray(TH2 * h, float * array, int size1, int size2)
2086 {
2087   int i, i1;
2088   int j, j1;
2089   float v1, ev1, v2, ev2, sum, esum;
2090   for (i=0, i1=1; i<size1; ++i,++i1)
2091     {
2092     for (j=0, j1=1; j<size2; ++j,++j1)
2093       {
2094       v1  = array[i*size2+j]; ev1 = sqrt(v1);
2095       v2  = h->GetBinContent(i1,j1);
2096       ev2 = h->GetBinError(i1,j1);
2097       sum = v1 + v2;
2098       esum = sqrt(ev1*ev1+ev2*ev2);
2099       h->SetBinContent(i1,j1,sum);
2100       h->SetBinError(i1,j1,esum);
2101       }
2102     }
2103 }
2104
2105 void  AliDptDptInMC::fillHistoWithArray(TH3 * h, float * array, int size1, int size2, int size3)
2106 {
2107   int i, i1;
2108   int j, j1;
2109   int k, k1;
2110   float v1, ev1, v2, ev2, sum, esum;
2111   int size23 = size2*size3;
2112   for (i=0, i1=1; i<size1; ++i,++i1)
2113     {
2114     for (j=0, j1=1; j<size2; ++j,++j1)
2115       {
2116       for (k=0, k1=1; k<size3; ++k,++k1)
2117         {
2118         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
2119         v2  = h->GetBinContent(i1,j1,k1);
2120         ev2 = h->GetBinError(i1,j1,k1);
2121         sum = v1 + v2;
2122         esum = sqrt(ev1*ev1+ev2*ev2);
2123         h->SetBinContent(i1,j1,k1,sum);
2124         h->SetBinError(i1,j1,k1,esum);
2125         }
2126       }
2127     }
2128 }
2129
2130 void  AliDptDptInMC::fillHistoWithArray(TH1 * h, double * array, int size)
2131 {
2132   int i, i1;
2133   double v1, ev1, v2, ev2, sum, esum;
2134   for (i=0, i1=1; i<size; ++i,++i1)
2135     {
2136     v1  = array[i]; ev1 = sqrt(v1);
2137     v2  = h->GetBinContent(i1);
2138     ev2 = h->GetBinError(i1);
2139     sum = v1 + v2;
2140     esum = sqrt(ev1*ev1+ev2*ev2);
2141     h->SetBinContent(i1,sum);
2142     h->SetBinError(i1,esum);
2143     }
2144 }
2145
2146 void  AliDptDptInMC::fillHistoWithArray(TH2 * h, double * array, int size1, int size2)
2147 {
2148   int i, i1;
2149   int j, j1;
2150   double v1, ev1, v2, ev2, sum, esum;
2151   for (i=0, i1=1; i<size1; ++i,++i1)
2152     {
2153     for (j=0, j1=1; j<size2; ++j,++j1)
2154       {
2155       v1  = array[i*size2+j]; ev1 = sqrt(v1);
2156       v2  = h->GetBinContent(i1,j1);
2157       ev2 = h->GetBinError(i1,j1);
2158       sum = v1 + v2;
2159       esum = sqrt(ev1*ev1+ev2*ev2);
2160       h->SetBinContent(i1,j1,sum);
2161       h->SetBinError(i1,j1,esum);
2162       }
2163     }
2164 }
2165
2166 void  AliDptDptInMC::fillHistoWithArray(TH3 * h, double * array, int size1, int size2, int size3)
2167 {
2168   int i, i1;
2169   int j, j1;
2170   int k, k1;
2171   double v1, ev1, v2, ev2, sum, esum;
2172   int size23 = size2*size3;
2173   for (i=0, i1=1; i<size1; ++i,++i1)
2174     {
2175     for (j=0, j1=1; j<size2; ++j,++j1)
2176       {
2177       for (k=0, k1=1; k<size3; ++k,++k1)
2178         {
2179         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
2180         v2  = h->GetBinContent(i1,j1,k1);
2181         ev2 = h->GetBinError(i1,j1,k1);
2182         sum = v1 + v2;
2183         esum = sqrt(ev1*ev1+ev2*ev2);
2184         h->SetBinContent(i1,j1,k1,sum);
2185         h->SetBinError(i1,j1,k1,esum);
2186         }
2187       }
2188     }
2189 }
2190
2191 //________________________________________________________________________
2192 double *  AliDptDptInMC::getDoubleArray(int size, double v)
2193 {
2194   /// Allocate an array of type double with n values
2195   /// Initialize the array to the given value
2196   double * array = new double [size];
2197   for (int i=0;i<size;++i) array[i]=v;
2198   return array;
2199 }
2200
2201 //________________________________________________________________________
2202 float *  AliDptDptInMC::getFloatArray(int size, float v)
2203 {
2204   /// Allocate an array of type float with n values
2205   /// Initialize the array to the given value
2206   float * array = new float [size];
2207   for (int i=0;i<size;++i) array[i]=v;
2208   return array;
2209 }
2210
2211
2212 //________________________________________________________________________
2213 TH1D * AliDptDptInMC::createHisto1D(const TString &  name, const TString &  title,
2214                                                       int n, double xMin, double xMax,
2215                                                       const TString &  xTitle, const TString &  yTitle)
2216 {
2217   //CreateHisto new 1D historgram
2218   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
2219   TH1D * h = new TH1D(name,title,n,xMin,xMax);
2220   h->GetXaxis()->SetTitle(xTitle);
2221   h->GetYaxis()->SetTitle(yTitle);
2222   addToList(h);
2223   return h;
2224 }
2225
2226
2227 //________________________________________________________________________
2228 TH1D * AliDptDptInMC::createHisto1D(const TString &  name, const TString &  title,
2229                                                       int n, double * bins,
2230                                                       const TString &  xTitle, const TString &  yTitle)
2231 {
2232   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
2233   TH1D * h = new TH1D(name,title,n,bins);
2234   h->GetXaxis()->SetTitle(xTitle);
2235   h->GetYaxis()->SetTitle(yTitle);
2236   addToList(h);
2237   return h;
2238 }
2239
2240
2241 //________________________________________________________________________
2242 TH2D * AliDptDptInMC::createHisto2D(const TString &  name, const TString &  title,
2243                                                       int nx, double xMin, double xMax, int ny, double yMin, double yMax,
2244                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2245 {
2246   AliInfo(Form("createHisto 2D histo %s  nx: %d  xMin: %f10.4 xMax: %f10.4  ny: %d   yMin: %f10.4 yMax: %f10.4",name.Data(),nx,xMin,xMax,ny,yMin,yMax));
2247   TH2D * h = new TH2D(name,title,nx,xMin,xMax,ny,yMin,yMax);
2248   h->GetXaxis()->SetTitle(xTitle);
2249   h->GetYaxis()->SetTitle(yTitle);
2250   h->GetZaxis()->SetTitle(zTitle);
2251   addToList(h);
2252   return h;
2253 }
2254
2255 //________________________________________________________________________
2256 TH2D * AliDptDptInMC::createHisto2D(const TString &  name, const TString &  title,
2257                                                       int nx, double* xbins, int ny, double yMin, double yMax,
2258                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2259 {
2260   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
2261   TH2D * h;
2262   h = new TH2D(name,title,nx,xbins,ny,yMin,yMax);
2263   h->GetXaxis()->SetTitle(xTitle);
2264   h->GetYaxis()->SetTitle(yTitle);
2265   h->GetZaxis()->SetTitle(zTitle);
2266   addToList(h);
2267   return h;
2268 }
2269
2270 //// F /////
2271 //________________________________________________________________________
2272 TH1F * AliDptDptInMC::createHisto1F(const TString &  name, const TString &  title,
2273                                                         int n, double xMin, double xMax,
2274                                                         const TString &  xTitle, const TString &  yTitle)
2275 {
2276   //CreateHisto new 1D historgram
2277   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
2278   TH1F * h = new TH1F(name,title,n,xMin,xMax);
2279   h->GetXaxis()->SetTitle(xTitle);
2280   h->GetYaxis()->SetTitle(yTitle);
2281   addToList(h);
2282   return h;
2283 }
2284
2285
2286 //________________________________________________________________________
2287 TH1F * AliDptDptInMC::createHisto1F(const TString &  name, const TString &  title,
2288                                                         int n, double * bins,
2289                                                         const TString &  xTitle, const TString &  yTitle)
2290 {
2291   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
2292   TH1F * h = new TH1F(name,title,n,bins);
2293   h->GetXaxis()->SetTitle(xTitle);
2294   h->GetYaxis()->SetTitle(yTitle);
2295   addToList(h);
2296   return h;
2297 }
2298
2299
2300 //________________________________________________________________________
2301 TH2F * AliDptDptInMC::createHisto2F(const TString &  name, const TString &  title,
2302                                                         int nx, double xMin, double xMax, int ny, double yMin, double yMax,
2303                                                         const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2304 {
2305   AliInfo(Form("createHisto 2D histo %s  nx: %d  xMin: %f10.4 xMax: %f10.4  ny: %d   yMin: %f10.4 yMax: %f10.4",name.Data(),nx,xMin,xMax,ny,yMin,yMax));
2306   TH2F * h = new TH2F(name,title,nx,xMin,xMax,ny,yMin,yMax);
2307   h->GetXaxis()->SetTitle(xTitle);
2308   h->GetYaxis()->SetTitle(yTitle);
2309   h->GetZaxis()->SetTitle(zTitle);
2310   addToList(h);
2311   return h;
2312 }
2313
2314 //________________________________________________________________________
2315 TH2F * AliDptDptInMC::createHisto2F(const TString &  name, const TString &  title,
2316                                                         int nx, double* xbins, int ny, double yMin, double yMax,
2317                                                         const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2318 {
2319   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
2320   TH2F * h;
2321   h = new TH2F(name,title,nx,xbins,ny,yMin,yMax);
2322   h->GetXaxis()->SetTitle(xTitle);
2323   h->GetYaxis()->SetTitle(yTitle);
2324   h->GetZaxis()->SetTitle(zTitle);
2325   addToList(h);
2326   return h;
2327 }
2328
2329
2330 //________________________________________________________________________
2331 TH3F * AliDptDptInMC::createHisto3F(const TString &  name, const TString &  title,
2332                                                       int nx, double xMin, double xMax,
2333                                                       int ny, double yMin, double yMax,
2334                                                       int nz, double zMin, double zMax,
2335                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2336 {
2337   AliInfo(Form("createHisto 2D histo %s  nx: %d  xMin: %f10.4 xMax: %f10.4  ny: %d   yMin: %f10.4 yMax: %f10.4 nz: %d   zMin: %f10.4 zMax: %f10.4",name.Data(),nx,xMin,xMax,ny,yMin,yMax,nz,zMin,zMax));
2338   TH3F * h = new TH3F(name,title,nx,xMin,xMax,ny,yMin,yMax,nz,zMin,zMax);
2339   h->GetXaxis()->SetTitle(xTitle);
2340   h->GetYaxis()->SetTitle(yTitle);
2341   h->GetZaxis()->SetTitle(zTitle);
2342   addToList(h);
2343   return h;
2344 }
2345
2346
2347 //________________________________________________________________________
2348 TProfile * AliDptDptInMC::createProfile(const TString & name, const TString & description,
2349                                                             int nx,double xMin,double xMax,
2350                                                             const TString &  xTitle, const TString &  yTitle)
2351 {
2352   AliInfo(Form("createHisto 1D profile %s   nBins: %d  xMin: %f10.4 xMax: %f10.4",name.Data(),nx,xMin,xMax));
2353   TProfile * h = new TProfile(name,description,nx,xMin,xMax);
2354   h->GetXaxis()->SetTitle(xTitle);
2355   h->GetYaxis()->SetTitle(yTitle);
2356   addToList(h);
2357   return h;
2358 }
2359
2360 //________________________________________________________________________
2361 TProfile * AliDptDptInMC::createProfile(const TString &  name,const TString &  description,
2362                                                             int nx,  double* bins,
2363                                                             const TString &  xTitle, const TString &  yTitle)
2364 {
2365   AliInfo(Form("createHisto 1D profile %s  with %d non uniform bins",name.Data(),nx));
2366   TProfile * h = new TProfile(name,description,nx,bins);
2367   h->GetXaxis()->SetTitle(xTitle);
2368   h->GetYaxis()->SetTitle(yTitle);
2369   addToList(h);
2370   return h;
2371 }
2372
2373
2374 void   AliDptDptInMC::addToList(TH1 *h)
2375 {
2376   if (_outputHistoList)
2377     {
2378     _outputHistoList->Add(h);
2379     }
2380   else
2381     AliInfo("addToList(TH1 *h) _outputHistoList is null!!!!! Should abort ship");
2382
2383 }
2384
2385
2386