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