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   __n1_1_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_1,    0.);
 785   __s1pt_1_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_1,    0.);
 786   __n1_1_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_1,  0.);
 787
 788   //cout << "==========================================================================================" << endl;
 789   //cout << "=============== Booking for particle 1 done." << endl;
 790   //cout << "_requestedCharge_1: " << _requestedCharge_1 << endl;
 791   //cout << "_requestedCharge_2: " << _requestedCharge_2 << endl;
 792
 793   if (_requestedCharge_2!=_requestedCharge_1)
 794     {
 795       //cout << " creating arrays for particle 2 with size: " << arraySize << endl;
 796     _sameFilter = 0;
 797     //particle 2
 798     _id_2       = new int[arraySize];
 799     _charge_2   = new int[arraySize];
 800     //_iPhi_2   = new int[arraySize];
 801     //_iEta_2   = new int[arraySize];
 802     _iEtaPhi_2  = new int[arraySize];
 803     _iPt_2      = new int[arraySize];
 804     _pt_2       = new float[arraySize];
 805     _px_2       = new float[arraySize];
 806     _py_2       = new float[arraySize];
 807     _pz_2       = new float[arraySize];
 808     //_phi_2    = new float[arraySize];
 809     //_eta_2    = new float[arraySize];
 810     _correction_2 = new float[arraySize];
 811     //_dedx_2       = new float[arraySize];
 812
 813     __n1_2_vsPt              = getDoubleArray(_nBins_pt_2,        0.);
 814     __n1_2_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_2,    0.);
 815     __s1pt_2_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_2,    0.);
 816     __n1_2_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_2, 0.);
 817
 818     }
 819
 820   __n2_12_vsPtPt           = getDoubleArray(_nBins_pt_1*_nBins_pt_2,0.);
 821   __n2_12_vsEtaPhi         = getFloatArray(_nBins_etaPhi_12,       0.);
 822   __s2ptpt_12_vsEtaPhi     = getFloatArray(_nBins_etaPhi_12,       0.);
 823   __s2PtN_12_vsEtaPhi      = getFloatArray(_nBins_etaPhi_12,       0.);
 824   __s2NPt_12_vsEtaPhi      = getFloatArray(_nBins_etaPhi_12,       0.);
 825
 826   // Setup all the labels needed.
 827
 828   part_1_Name   = "_1";
 829   part_2_Name   = "_2";
 830   pair_12_Name  = "_12";
 831
 832   n1Name     = "n1";
 833   n2Name     = "n2";
 834   n1NwName   = "n1Nw";
 835   n2NwName   = "n2Nw";
 836   r1Name     = "r1";
 837   r2Name     = "r2";
 838   r3Name     = "r3";
 839   r2r1Name   = "r2r1";
 840   c2Name     = "c2";
 841   c3Name     = "c3";
 842   d3Name     = "d3";
 843   p3Name     = "p3";
 844   cName      = "sean";
 845
 846   intR2Name       = "intR2";
 847   binCorrName     = "binCorr";
 848   intBinCorrName  = "intBinCorr";
 849
 850   avgName      = "avg";
 851   sumName      = "sum";
 852   s1ptName     = "sumPt";
 853   s1ptNwName   = "sumPtNw";
 854   s1DptName    = "sumDpt";
 855   s2PtPtName   = "sumPtPt";
 856   s2PtPtNwName = "sumPtPtNw";
 857   s2DptDptName = "sumDptDpt";
 858   s2NPtName    = "sumNPt";
 859   s2NPtNwName  = "sumNPtNw";
 860   s2PtNName    = "sumPtN";
 861   s2NPtNwName  = "sumNPtNw";
 862   s2PtNNwName  = "sumPtNNw";
 863   ptName       = "avgPt";
 864   ptptName     = "avgPtPt";
 865   pt1pt1Name   = "avgPtavgPt";
 866   DptName      = "avgDpt";
 867   DptDptName   = "avgDptDpt";
 868   RDptDptName  = "relDptDpt"; // ratio of avgDptDpt by avgPt*avgPt
 869   nPtName      = "avgNpt";
 870   ptNName      = "avgPtN";
 871   seanName     = "seanC";
 872
 873   _title_counts = "yield";
 874
 875   _title_m0     = "M_{0}";
 876   _title_m1     = "M_{1}";
 877   _title_m2     = "M_{2}";
 878   _title_m3     = "M_{3}";
 879   _title_m4     = "V0Centrality";
 880   _title_m5     = "TrkCentrality";
 881   _title_m6     = "SpdCentrality";
 882
 883   _title_eta_1       = "#eta_{1}";
 884   _title_phi_1       = "#varphi_{1} (radian)";
 885   _title_etaPhi_1    = "#eta_{1}#times#varphi_{1}";
 886   _title_pt_1        = "p_{t,1} (GeV/c)";
 887   _title_n_1         = "n_{1}";
 888   _title_SumPt_1     = "#Sigma p_{t,1} (GeV/c)";
 889   _title_AvgPt_1     = "#LT p_{t,1} #GT (GeV/c)";
 890   _title_AvgN_1      = "#LT n_{1} #GT";
 891   _title_AvgSumPt_1  = "#LT #Sigma p_{t,1} #GT (GeV/c)";
 892
 893   _title_eta_2       = "#eta_{2}";
 894   _title_phi_2       = "#varphi_{2} (radian)";
 895   _title_etaPhi_2    = "#eta_{2}#times#varphi_{2}";
 896   _title_pt_2        = "p_{t,2} (GeV/c)";
 897   _title_n_2         = "n_{2}";
 898   _title_SumPt_2     = "#Sigma p_{t,1} (GeV/c)";
 899   _title_AvgPt_2     = "#LT p_{t,2} #GT (GeV/c)";
 900   _title_AvgN_2      = "#LT n_{2} #GT";
 901   _title_AvgSumPt_2  = "#LT #Sigma p_{t,2} #GT (GeV/c)";
 902
 903   _title_etaPhi_12   = "#eta_{1}#times#varphi_{1}#times#eta_{2}#times#varphi_{2}";
 904
 905   _title_AvgN2_12       = "#LT n_{2} #GT";;
 906   _title_AvgSumPtPt_12  = "#LT #Sigma p_{t,1}p_{t,2} #GT";;
 907   _title_AvgSumPtN_12   = "#LT #Sigma p_{t,1}N #GT";;
 908   _title_AvgNSumPt_12   = "#LT N#Sigma p_{t,2} #GT";;
 909
 910
 911   vsZ         = "_vsZ";
 912   vsM         = "_vsM";
 913   vsPt         = "_vsPt";
 914   vsPhi        = "_vsPhi";
 915   vsEta        = "_vsEta";
 916   vsEtaPhi     = "_vsEtaPhi";
 917   vsPtVsPt     = "_vsPtVsPt";
 918
 919
 920   if (_useWeights)
 921     {
 922     int iZ, iEtaPhi, iPt;
 923     int iZ1,iEtaPhi1,iPt1;
 924     int a, b;
 925     if (_weight_1)
 926       {
 927       _correctionWeight_1 = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
 928       a = _nBins_pt_1;
 929       b = _nBins_etaPhi_1*_nBins_pt_1;
 930       for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 931         {
 932         for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_1; iEtaPhi++, iEtaPhi1++)
 933           {
 934           for (iPt=0,iPt1=1; iPt<_nBins_pt_1; iPt++, iPt1++)
 935             {
 936             _correctionWeight_1[iZ*b+iEtaPhi*a+iPt] = _weight_1->GetBinContent(iZ1,iEtaPhi1,iPt1);
 937             }
 938           }
 939         }
 940       } // _weight_1
 941     else
 942       {
 943       AliError("AliDptDptInMC:: _weight_1 is a null pointer.");
 944       return;
 945       }
 946     if (!_sameFilter)
 947       {
 948       if (_weight_2)
 949         {
 950         _correctionWeight_2 = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
 951         a = _nBins_pt_2;
 952         b = _nBins_etaPhi_2*_nBins_pt_2;
 953         for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 954           {
 955           for (iEtaPhi=0,iEtaPhi1=1; iEtaPhi<_nBins_etaPhi_2; iEtaPhi++, iEtaPhi1++)
 956             {
 957             for (iPt=0,iPt1=1; iPt<_nBins_pt_2; iPt++, iPt1++)
 958               {
 959               _correctionWeight_2[iZ*b+iEtaPhi*a+iPt] = _weight_2->GetBinContent(iZ1,iEtaPhi1,iPt1);
 960               }
 961             }
 962           }
 963         } // _weight_2
 964       else
 965         {
 966         AliError("AliDptDptInMC:: _weight_1 is a null pointer.");
 967         return;
 968         }
 969       }
 970     }
 971
 972   createHistograms();
 973   PostData(0,_outputHistoList);
 974
 975   //cout<< "AliDptDptInMC::CreateOutputObjects() DONE " << endl;
 976
 977 }
 978
 979 void  AliDptDptInMC::createHistograms()
 980 {
 981   AliInfo(" AliDptDptInMC::createHistoHistograms() Creating Event Histos");
 982   TString name;
 983
 984   name = "eventAccounting";
 985
 986   // bin index : what it is...
 987   //         0 :  number of event submitted
 988   //         1 :  number accepted by centrality cut
 989   //         2 :  number accepted by centrality cut and z cut
 990   //         3 :  total number of particles that satisfied filter 1
 991   //         4 :  total number of particles that satisfied filter 2
 992   _eventAccounting      = createHisto1D(name,name,10, -0.5, 9.5, "event Code", _title_counts);
 993
 994   name = "m0"; _m0      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m0, _title_counts);
 995   name = "m1"; _m1      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m1, _title_counts);
 996   name = "m2"; _m2      = createHisto1D(name,name,_nBins_M2, _min_M2, _max_M2, _title_m2, _title_counts);
 997   name = "m3"; _m3      = createHisto1D(name,name,_nBins_M3, _min_M3, _max_M3, _title_m3, _title_counts);
 998   name = "m4"; _m4      = createHisto1D(name,name,_nBins_M4, _min_M4, _max_M4, _title_m4, _title_counts);
 999   name = "m5"; _m5      = createHisto1D(name,name,_nBins_M5, _min_M5, _max_M5, _title_m5, _title_counts);
1000   name = "m6"; _m6      = createHisto1D(name,name,_nBins_M6, _min_M6, _max_M6, _title_m6, _title_counts);
1001   name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
1002
1003
1004   name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
1005   name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
1006   name = "DCAz";    _dcaz     = createHisto1F(name,name, 500, -5.0, 5.0, "dcaZ","counts");
1007   name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 500, -5.0, 5.0, "dcaXY","counts");
1008
1009   /*name = "dedxVsP_1"; _dedxVsP_1  = createHisto2D(name,name,1000,-10.,10.,1000,0.,1000.,"p (GeV/c)", "dedx", "counts");
1010   name = "dedxVsP_2"; _dedxVsP_2  = createHisto2D(name,name,1000,-10.,10.,1000,0.,1000.,"p (GeV/c)", "dedx", "counts");
1011   name = "corrDedxVsP_1"; _corrDedxVsP_1 = createHisto2D(name,name,1000,-10.,10.,1000,0.,500,"p (GeV/c)", "dedx", "counts");
1012   name = "corrDedxVsP_2"; _corrDedxVsP_2 = createHisto2D(name,name,1000,-10.,10.,1000,0.,500,"p (GeV/c)", "dedx", "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     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);
1021     //name = "dedxVsP_1";                          _dedxVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1022     //name = "corrDedxVsP_1";                      _corrDedxVsP_1          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1023     //name = "betaVsP_1";                          _betaVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
1024
1025     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);
1026     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);
1027     //name = "dedxVsP_2";                          _dedxVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1028     //name = "corrDedxVsP_2";                      _corrDedxVsP_2          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
1029     //name = "betaVsP_2";                          _betaVsP_2              = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
1030
1031     }
1032   else
1033     {
1034     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);
1035     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);
1036     name = n1Name+part_1_Name+vsM;            _n1_1_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1037     name = s1ptName+part_1_Name+vsM;          _s1pt_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1038     name = n1NwName+part_1_Name+vsM;          _n1Nw_1_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1039     name = s1ptNwName+part_1_Name+vsM;        _s1ptNw_1_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1040
1041     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);
1042     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);
1043     name = n1Name+part_2_Name + vsM;          _n1_2_vsM             = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_2);
1044     name = s1ptName+part_2_Name + vsM;        _s1pt_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_2);
1045     name = n1NwName+part_2_Name+vsM;          _n1Nw_2_vsM           = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN_1);
1046     name = s1ptNwName+part_2_Name+vsM;        _s1ptNw_2_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPt_1);
1047
1048     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);
1049     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);
1050     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);
1051     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);
1052     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);
1053
1054     name = n2Name+pair_12_Name + vsM;         _n2_12_vsM            = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN2_12);
1055     name = s2PtPtName+pair_12_Name + vsM;     _s2PtPt_12_vsM        = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtPt_12);
1056     name = s2PtNName+pair_12_Name + vsM;      _s2PtN_12_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);
1057     name = s2NPtName+pair_12_Name + vsM;      _s2NPt_12_vsM         = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);
1058
1059     name = n2NwName+pair_12_Name + vsM;       _n2Nw_12_vsM          = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgN2_12);
1060     name = s2PtPtNwName+pair_12_Name + vsM;   _s2PtPtNw_12_vsM      = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtPt_12);
1061     name = s2PtNNwName+pair_12_Name + vsM;    _s2PtNNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);
1062     name = s2NPtNwName+pair_12_Name + vsM;    _s2NPtNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);
1063
1064     name = "mInv";     _invMass     = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts");
1065     name = "mInvElec"; _invMassElec = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts");
1066     }
1067
1068   AliInfo(" AliDptDptInMC::createHistoHistograms() All Done");
1069 }
1070 //-----------------------//
1071
1072 void  AliDptDptInMC::finalizeHistograms()
1073 {
1074
1075   AliInfo("AliDptDptInMC::finalizeHistograms() starting");
1076   AliInfo(Form("CorrelationAnalyzers::finalizeHistograms()   _eventCount : %d",int(_eventCount)));
1077   if (_singlesOnly)
1078     {
1079     if (_sameFilter)
1080       {
1081       fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1082       fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1083       fillHistoWithArray(_n1_2_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1084       fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1085       }
1086     else
1087       {
1088       fillHistoWithArray(_n1_1_vsPt,              __n1_1_vsPt,        _nBins_pt_1);
1089       fillHistoWithArray(_n1_1_vsZVsEtaVsPhiVsPt, __n1_1_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
1090       fillHistoWithArray(_n1_2_vsPt,              __n1_2_vsPt,        _nBins_pt_2);
1091       fillHistoWithArray(_n1_2_vsZVsEtaVsPhiVsPt, __n1_2_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_2, _nBins_pt_2);
1092       }
1093     }
1094   else
1095     {
1096     if (_sameFilter)
1097       {
1098       fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1099       fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1100       fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1101       fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1102       }
1103     else
1104       {
1105       fillHistoWithArray(_n1_1_vsEtaVsPhi,        __n1_1_vsEtaPhi,    _nBins_eta_1,   _nBins_phi_1);
1106       fillHistoWithArray(_s1pt_1_vsEtaVsPhi,      __s1pt_1_vsEtaPhi,  _nBins_eta_1,   _nBins_phi_1);
1107       fillHistoWithArray(_n1_2_vsEtaVsPhi,        __n1_2_vsEtaPhi,    _nBins_eta_2,   _nBins_phi_2);
1108       fillHistoWithArray(_s1pt_2_vsEtaVsPhi,      __s1pt_2_vsEtaPhi,  _nBins_eta_2,   _nBins_phi_2);
1109       }
1110     fillHistoWithArray(_n2_12_vsEtaPhi,     __n2_12_vsEtaPhi,     _nBins_etaPhi_12);
1111     fillHistoWithArray(_s2PtPt_12_vsEtaPhi, __s2ptpt_12_vsEtaPhi, _nBins_etaPhi_12);
1112     fillHistoWithArray(_s2PtN_12_vsEtaPhi,  __s2PtN_12_vsEtaPhi,  _nBins_etaPhi_12);
1113     fillHistoWithArray(_s2NPt_12_vsEtaPhi,  __s2NPt_12_vsEtaPhi,  _nBins_etaPhi_12);
1114     fillHistoWithArray(_n2_12_vsPtVsPt,     __n2_12_vsPtPt,       _nBins_pt_1,    _nBins_pt_2);
1115
1116     }
1117   AliInfo("AliDptDptInMC::finalizeHistograms()  Done ");
1118 }
1119 //--------------//
1120
1121
1122 void  AliDptDptInMC::UserExec(Option_t */*option*/)
1123 {
1124
1125
1126   int    k1,k2;
1127   int    iPhi, iEta, iEtaPhi, iPt, charge;
1128   float  q, phi, pt, eta, corr, corrPt, px, py, pz;
1129   int    ij;
1130   int    id_1, q_1, iEtaPhi_1, iPt_1;
1131   float  pt_1, px_1, py_1, pz_1, corr_1;
1132   int    id_2, q_2, iEtaPhi_2, iPt_2;
1133   float  pt_2, px_2, py_2, pz_2, corr_2;
1134   float  ptpt;
1135   int    iVertex, iVertexP1, iVertexP2;
1136   int    iZEtaPhiPt;
1137   float  massElecSq = 2.5e-7;
1138   const  AliAODVertex*  vertex;
1139   //int    nClus;
1140   bool   bitOK;
1141
1142   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
1143   if (!manager) {
1144     //cout<<"ERROR: Analysis manager not found."<<endl;
1145     return;
1146   }
1147   //coneect to the inputHandler------------
1148   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
1149   if (!inputHandler) {
1150     return;
1151   }
1152
1153   fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
1154   if (!fAODEvent)
1155     {
1156       return;
1157     }
1158   fPIDResponse =inputHandler->GetPIDResponse();
1159   if (!fPIDResponse){
1160     AliFatal("This Task needs the PID response attached to the inputHandler");
1161     return;
1162   }
1163   // count all events looked at here
1164   _eventCount++;
1165
1166   if (_eventAccounting)
1167     {
1168       _eventAccounting->Fill(0);// count all calls to this function
1169     }
1170   else
1171     {
1172       return;
1173     }
1174
1175   _eventAccounting->Fill(1);// count all calls to this function with a valid pointer
1176   //reset single particle counters
1177   k1 = k2 = 0;
1178   __n1_1 = __n1_2 = __s1pt_1 = __s1pt_2 = __n1Nw_1 = __n1Nw_2 = __s1ptNw_1 = __s1ptNw_2 = 0;
1179
1180   float v0Centr  = -999.;
1181   float v0ACentr  = -999.;
1182   float trkCentr = -999.;
1183   float spdCentr = -999.;
1184   float vertexX  = -999;
1185   float vertexY  = -999;
1186   float vertexZ  = -999;
1187   float vertexXY = -999;
1188
1189   float centrality = -999;
1190   //Double_t nSigma =-999;
1191   if(fAODEvent)
1192     {
1193
1194       //Centrality
1195       //AliAODHeader* centralityObject = fAODEvent->GetHeader()->GetCentralityP();
1196       AliCentrality* centralityObject =  fAODEvent->GetHeader()->GetCentralityP();
1197       if (centralityObject)
1198         {
1199           v0Centr  = centralityObject->GetCentralityPercentile("V0M");
1200           v0ACentr  = centralityObject->GetCentralityPercentile("V0A");
1201           trkCentr = centralityObject->GetCentralityPercentile("TRK");
1202           spdCentr = centralityObject->GetCentralityPercentile("CL1");
1203
1204         }
1205
1206       _mult4    = v0Centr;
1207       _mult4a   = v0ACentr;
1208       _mult5    = trkCentr;
1209       _mult6    = spdCentr;
1210       _field    = fAODEvent->GetMagneticField();
1211       //_centralityMethod
1212       switch (_centralityMethod)
1213         {
1214         case 0: centrality = _mult0; break;
1215         case 1: centrality = _mult1; break;
1216         case 2: centrality = _mult2; break;
1217         case 3: centrality = _mult3; break;
1218         case 4: centrality = _mult4; break;
1219         case 5: centrality = _mult5; break;
1220         case 6: centrality = _mult6; break;
1221         case 7: centrality = _mult4a; break;
1222         }
1223
1224       if ( centrality < _centralityMin ||  centrality > _centralityMax ) return;
1225
1226       _eventAccounting->Fill(2);// count all events with right centrality
1227
1228       vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertex();
1229
1230       if(vertex->GetNContributors() > 0)
1231         {
1232           vertexX = vertex->GetX();
1233           vertexY = vertex->GetY();
1234           vertexZ = vertex->GetZ();
1235           vertexXY = sqrt(vertexX*vertexX+vertexY*vertexY);
1236         }
1237       if (vertexZ  < _vertexZMin  || vertexZ  > _vertexZMax ) return;
1238
1239       iVertex = int((vertexZ-_min_vertexZ)/_width_vertexZ);
1240       iVertexP1 = iVertex*_nBins_etaPhiPt_1;
1241       iVertexP2 = iVertex*_nBins_etaPhiPt_2;
1242       if (iVertex<0 || iVertex>=_nBins_vertexZ)
1243         {
1244           AliError("AliTaskDptCorrMC::Exec(Option_t *option) iVertex<0 || iVertex>=_nBins_vertexZ ");
1245           return;
1246         }
1247       _eventAccounting->Fill(3);// count all calls to this function with a valid pointer
1248       //======================
1249       //***********************************
1250       //MC AOD Truth
1251       if(fAnalysisType == "MCAOD")
1252         {
1253           AliMCEvent* mcEvent = MCEvent();
1254           _nTracks = mcEvent->GetNumberOfTracks();
1255           _mult3    = _nTracks;
1256           //loop over tracks starts here
1257           for (int iTrack=0; iTrack< _nTracks; iTrack++)
1258             {
1259               AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTrack);
1260
1261               if (!aodTrack)
1262                 {
1263                   AliError(Form("AliTaskDptCorrMC::Exec(Option_t *option) No track ofr iTrack=%d", iTrack));
1264                   continue;
1265                 }
1266
1267               //if(!aodTrack->IsPhysicalPrimary()) continue;
1268
1269               q      = aodTrack->Charge();
1270               charge = int(q);
1271               phi    = aodTrack->Phi();
1272               pt     = aodTrack->Pt();
1273               px     = aodTrack->Px();
1274               py     = aodTrack->Py();
1275               pz     = aodTrack->Pz();
1276               eta    = aodTrack->Eta();
1277               // Kinematics cuts from ESD track cuts
1278               if( pt < 0.2 || pt > 2.0)      continue;
1279               if( eta < -0.8 || eta > 0.8)  continue;
1280
1281               _etadis->Fill(eta);
1282               _phidis->Fill(phi);
1283               // Remove neutral tracks
1284                if(q == 0) continue;
1285
1286               if(fExcludeResonancesInMC)
1287                 {
1288                   //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
1289                   Int_t gMotherIndex = aodTrack->GetMother();
1290                   if(gMotherIndex != -1) {
1291                     AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1292                     if(motherTrack) {
1293                       Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1294
1295                       if(pdgCodeOfMother == 311  ||
1296                          pdgCodeOfMother == -311 ||
1297                          pdgCodeOfMother == 310  ||
1298                          pdgCodeOfMother == 3122 ||
1299                          pdgCodeOfMother == -3122 ||
1300                          pdgCodeOfMother == 111 ||
1301                          pdgCodeOfMother == 22 ) continue;
1302                     }
1303                   }
1304                 }
1305               //Exclude electrons with PDG
1306
1307               if(fExcludeElectronsInMC) {
1308                 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1309               }
1310
1311               //Particle loop 1st particle
1312               if (q > 0)
1313                 {
1314
1315                   iPhi   = int( phi/_width_phi_1);
1316
1317                   if (iPhi<0 || iPhi>=_nBins_phi_1 )
1318                     {
1319                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1320                       return;
1321                     }
1322
1323                   iEta    = int((eta-_min_eta_1)/_width_eta_1);
1324                   if (iEta<0 || iEta>=_nBins_eta_1)
1325                     {
1326                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1327                       continue;
1328                     }
1329
1330                   iPt     = int((pt -_min_pt_1 )/_width_pt_1 );
1331                   if (iPt<0  || iPt >=_nBins_pt_1)
1332                     {
1333                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1334                       continue;
1335                     }
1336                   iEtaPhi = iEta*_nBins_phi_1+iPhi;
1337                   iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
1338
1339                   if (_correctionWeight_1)
1340                     corr = _correctionWeight_1[iZEtaPhiPt];
1341                   else
1342                     corr = 1;
1343                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1)
1344                     {
1345                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1");
1346                       continue;
1347                     }
1348
1349                   if (_singlesOnly)
1350                     {
1351
1352                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1353                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1354
1355                     }
1356
1357                   else
1358                     {
1359                       corrPt                      = corr*pt;
1360                       _id_1[k1]                   = iTrack;
1361                       _charge_1[k1]               = charge;
1362                       _iEtaPhi_1[k1]              = iEtaPhi;
1363                       _iPt_1[k1]                  = iPt;
1364                       _pt_1[k1]                   = pt;
1365                       _px_1[k1]                   = px;
1366                       _py_1[k1]                   = py;
1367                       _pz_1[k1]                   = pz;
1368                       _correction_1[k1]           = corr;
1369                       __n1_1                      += corr;
1370                       __n1_1_vsEtaPhi[iEtaPhi]    += corr;
1371                       __s1pt_1                    += corrPt;
1372                       __s1pt_1_vsEtaPhi[iEtaPhi]  += corrPt;
1373                       __n1Nw_1                    += 1;
1374                       __s1ptNw_1                  += pt;
1375                       ++k1;
1376                       if (k1>=arraySize)
1377                         {
1378                           AliError(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
1379                           return;
1380                         }
1381                     }
1382                 }
1383
1384               if (q < 0)
1385                 {
1386                   iPhi   = int( phi/_width_phi_2);
1387
1388                   if (iPhi<0 || iPhi>=_nBins_phi_2 )
1389                     {
1390                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1391                       return;
1392                     }
1393
1394                   iEta    = int((eta-_min_eta_2)/_width_eta_2);
1395                   if (iEta<0 || iEta>=_nBins_eta_2)
1396                     {
1397                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1398                       continue;
1399                     }
1400
1401                   iPt     = int((pt -_min_pt_2 )/_width_pt_2 );
1402                   if (iPt<0  || iPt >=_nBins_pt_2)
1403                     {
1404                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1405                       continue;
1406                     }
1407
1408                   iEtaPhi = iEta*_nBins_phi_2+iPhi;
1409                   iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
1410                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
1411                     {
1412                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
1413                       continue;
1414                     }
1415
1416                   if (_correctionWeight_2)
1417                     corr = _correctionWeight_2[iZEtaPhiPt];
1418                   else
1419                     corr = 1;
1420                   //dpt = pt - (charge>0) ? _avgPt_vsEtaPhi_2p[iEtaPhi] : _avgPt_vsEtaPhi_2m[iEtaPhi];
1421
1422                   if (_singlesOnly)
1423                     {
1424                       __n1_2_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1425                       __n1_2_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1426                     }
1427                   else
1428                     {
1429                       corrPt                      = corr*pt;
1430                       _id_2[k2]                   = iTrack;         //cout << "step 1" << endl;
1431                       _charge_2[k2]               = charge;         //cout << "step 2" << endl;
1432                       _iEtaPhi_2[k2]              = iEtaPhi;        //cout << "step 3" << endl;
1433                       _iPt_2[k2]                  = iPt;            //cout << "step 4" << endl;
1434                       _pt_2[k2]                   = pt;             //cout << "step 5" << endl;
1435                       _px_2[k2]                   = px;             //cout << "step 6" << endl;
1436                       _py_2[k2]                   = py;             //cout << "step 7" << endl;
1437                       _pz_2[k2]                   = pz;             //cout << "step 8" << endl;
1438                       _correction_2[k2]           = corr;           //cout << "step 9" << endl;
1439                       __n1_2                      += corr;          //cout << "step 10" << endl;
1440                       __s1pt_2                    += corrPt;        //cout << "step 13" << endl;
1441                       __n1Nw_2                    += 1;
1442                       __n1_2_vsEtaPhi[iEtaPhi]    += corr;          //cout << "step 11" << endl;
1443                       __s1pt_2_vsEtaPhi[iEtaPhi]  += corrPt;        //cout << "step 14" << endl;
1444                       __s1ptNw_2                  += pt;
1445                       ++k2;
1446                       if (k2>=arraySize)
1447                         {
1448                           AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize));
1449                           return;
1450                         }
1451                     }
1452
1453                   //cout << "done with track" << endl;
1454                 } //iTrack
1455             } //data aod loop
1456           } //MC AOD loop ends
1457
1458       //***********************************************
1459       //MC AOD Reconstructed tracks
1460       if(fAnalysisType == "MCAODreco")
1461         {
1462
1463           fArrayMC = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1464           if (!fArrayMC) {
1465             AliError("No array of MC particles found !!!");
1466           }
1467
1468           AliMCEvent* mcEvent = MCEvent();
1469           if (!mcEvent) {
1470             AliError("ERROR: Could not retrieve MC event");
1471           }
1472
1473           for (int iTrack=0; iTrack< fAODEvent->GetNumberOfTracks(); iTrack++)
1474             {
1475
1476               AliAODTrack *t = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(iTrack));
1477
1478               if(!t) {
1479                 AliError(Form("ERROR: Could not retrieve AODtrack %d",iTrack));
1480                 continue;
1481               }
1482
1483               bitOK  = t->TestFilterBit(_trackFilterBit);
1484               if (!bitOK) continue;
1485
1486               q      = t->Charge();
1487               charge = int(q);
1488               phi    = t->Phi();
1489               pt     = t->Pt();
1490               px     = t->Px();
1491               py     = t->Py();
1492               pz     = t->Pz();
1493               eta    = t->Eta();
1494
1495               Double_t nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(t,(AliPID::EParticleType)AliPID::kElectron));
1496               Double_t nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(t,(AliPID::EParticleType)AliPID::kPion));
1497               Double_t nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(t,(AliPID::EParticleType)AliPID::kKaon));
1498               Double_t nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(t,(AliPID::EParticleType)AliPID::kProton));
1499
1500               if(nsigmaelectron  < fNSigmaCut &&
1501                  nsigmapion      > fNSigmaCut &&
1502                  nsigmakaon      > fNSigmaCut &&
1503                  nsigmaproton    > fNSigmaCut ) continue;
1504
1505               Float_t dcaXY = t->DCA();
1506               Float_t dcaZ  = t->ZAtDCA();
1507
1508               // Kinematics cuts
1509               if( pt < 0.2 || pt > 2.0)  continue;
1510               if( eta < -0.8 || eta > 0.8)  continue;
1511
1512               Int_t label = TMath::Abs(t->GetLabel());
1513               AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1514
1515               //W/Wo Secondaries
1516               //if (!AODmcTrack->IsPhysicalPrimary()) continue;
1517
1518
1519                if(fExcludeResonancesInMC)
1520                  {
1521                    //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
1522                    Int_t gMotherIndex = AODmcTrack->GetMother();
1523                    if(gMotherIndex != -1) {
1524                      AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1525                      if(motherTrack) {
1526                        Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1527
1528                        if(pdgCodeOfMother == 311  ||
1529                           pdgCodeOfMother == -311 ||
1530                           pdgCodeOfMother == 310  ||
1531                           pdgCodeOfMother == 3122 ||
1532                           pdgCodeOfMother == -3122 ||
1533                           pdgCodeOfMother == 111 ||
1534                           pdgCodeOfMother == 22 ) continue;
1535                      }
1536                    }
1537                  }
1538
1539                if (AODmcTrack)
1540                  {
1541                    if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1542                  }
1543
1544                _etadis->Fill(eta);
1545                _phidis->Fill(phi);
1546                _dcaz->Fill(dcaZ);
1547                _dcaxy->Fill(dcaXY);
1548
1549               //Particle 1
1550               if (t->Charge() > 0)
1551                 {
1552
1553                   iPhi   = int( phi/_width_phi_1);
1554
1555                   if (iPhi<0 || iPhi>=_nBins_phi_1 )
1556                     {
1557                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1558                       return;
1559                     }
1560
1561                   iEta    = int((eta-_min_eta_1)/_width_eta_1);
1562                   if (iEta<0 || iEta>=_nBins_eta_1)
1563                     {
1564                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1565                       continue;
1566                     }
1567                   iPt     = int((pt -_min_pt_1 )/_width_pt_1 );
1568                   if (iPt<0  || iPt >=_nBins_pt_1)
1569                     {
1570                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1571                       continue;
1572                     }
1573                   iEtaPhi = iEta*_nBins_phi_1+iPhi;
1574                   iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
1575
1576                   if (_correctionWeight_1)
1577                     corr = _correctionWeight_1[iZEtaPhiPt];
1578                   else
1579                     corr = 1;
1580                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1)
1581                     {
1582                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_1");
1583                       continue;
1584                     }
1585
1586
1587                   if (_singlesOnly)
1588                     {
1589
1590                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1591                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1592
1593                     }
1594
1595                   else
1596                     {
1597                       corrPt                      = corr*pt;
1598                       _id_1[k1]                   = iTrack;
1599                       _charge_1[k1]               = charge;
1600                       _iEtaPhi_1[k1]              = iEtaPhi;
1601                       _iPt_1[k1]                  = iPt;
1602                       _pt_1[k1]                   = pt;
1603                       _px_1[k1]                   = px;
1604                       _py_1[k1]                   = py;
1605                       _pz_1[k1]                   = pz;
1606                       _correction_1[k1]           = corr;
1607                       __n1_1                      += corr;
1608                       __n1_1_vsEtaPhi[iEtaPhi]    += corr;
1609                       __s1pt_1                    += corrPt;
1610                       __s1pt_1_vsEtaPhi[iEtaPhi]  += corrPt;
1611                       __n1Nw_1                    += 1;
1612                       __s1ptNw_1                  += pt;
1613                       ++k1;
1614                       if (k1>=arraySize)
1615                         {
1616                           AliError(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
1617                           return;
1618                         }
1619                     }
1620                 }
1621
1622               if (t->Charge() < 0 )
1623                 {
1624                   //===================================
1625
1626                   iPhi   = int( phi/_width_phi_2);
1627
1628                   if (iPhi<0 || iPhi>=_nBins_phi_2 )
1629                     {
1630                       AliWarning("AliTaskDptCorrMC::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1631                       return;
1632                     }
1633
1634                   iEta    = int((eta-_min_eta_2)/_width_eta_2);
1635                   if (iEta<0 || iEta>=_nBins_eta_2)
1636                     {
1637                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1638                       continue;
1639                     }
1640                   iPt     = int((pt -_min_pt_2 )/_width_pt_2 );
1641                   if (iPt<0  || iPt >=_nBins_pt_2)
1642                     {
1643                       AliWarning(Form("AliTaskDptCorrMC::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1644                       continue;
1645                     }
1646
1647                   iEtaPhi = iEta*_nBins_phi_2+iPhi;
1648                   iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
1649                   if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
1650                     {
1651                       AliWarning("AliTaskDptCorrMC::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
1652                       continue;
1653                     }
1654
1655
1656                   if (_correctionWeight_2)
1657                     corr = _correctionWeight_2[iZEtaPhiPt];
1658                   else
1659                     corr = 1;
1660                   //dpt = pt - (charge>0) ? _avgPt_vsEtaPhi_2p[iEtaPhi] : _avgPt_vsEtaPhi_2m[iEtaPhi];
1661
1662                   if (_singlesOnly)
1663                     {
1664                       __n1_2_vsPt[iPt]               += corr;          //cout << "step 15" << endl;
1665                       __n1_2_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;
1666                     }
1667
1668                   else
1669                     {
1670                       corrPt                      = corr*pt;
1671                       _id_2[k2]                   = iTrack;         //cout << "step 1" << endl;
1672                       _charge_2[k2]               = charge;         //cout << "step 2" << endl;
1673                       _iEtaPhi_2[k2]              = iEtaPhi;        //cout << "step 3" << endl;
1674                       _iPt_2[k2]                  = iPt;            //cout << "step 4" << endl;
1675                       _pt_2[k2]                   = pt;             //cout << "step 5" << endl;
1676                       _px_2[k2]                   = px;             //cout << "step 6" << endl;
1677                       _py_2[k2]                   = py;             //cout << "step 7" << endl;
1678                       _pz_2[k2]                   = pz;             //cout << "step 8" << endl;
1679                       _correction_2[k2]           = corr;           //cout << "step 9" << endl;
1680                       __n1_2                      += corr;          //cout << "step 10" << endl;
1681                       __s1pt_2                    += corrPt;        //cout << "step 13" << endl;
1682                       __n1Nw_2                    += 1;
1683                       __n1_2_vsEtaPhi[iEtaPhi]    += corr;          //cout << "step 11" << endl;
1684                       __s1pt_2_vsEtaPhi[iEtaPhi]  += corrPt;        //cout << "step 14" << endl;
1685                       __s1ptNw_2                  += pt;
1686                       ++k2;
1687
1688                       if (k2>=arraySize)
1689                         {
1690                           AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize));
1691                           return;
1692                         }
1693                     }
1694
1695                   //cout << "done with track" << endl;
1696                 } //iTrack
1697             } //data aod loop
1698         } //MC AOD loop ends
1699
1700       //************************************************
1701     } //AOD events
1702
1703
1704   _m0->Fill(_mult0);
1705   _m1->Fill(_mult1);
1706   _m2->Fill(_mult2);
1707   _m3->Fill(_mult3);
1708   _m4->Fill(_mult4);
1709   _m5->Fill(_mult5);
1710   _m6->Fill(_mult6);
1711   _vertexZ->Fill(vertexZ);
1712
1713   if (_singlesOnly)
1714     {
1715       // nothing to do here.
1716     }
1717   else
1718     {
1719       if (_sameFilter)
1720         {
1721           _n1_1_vsM->Fill(centrality,      __n1_1);
1722           _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
1723           _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
1724           _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
1725           _n1_2_vsM->Fill(centrality,      __n1_1);
1726           _s1pt_2_vsM->Fill(centrality,    __s1pt_1);
1727           _n1Nw_2_vsM->Fill(centrality,    __n1Nw_1);
1728           _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_1);
1729           // reset pair counters
1730           __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
1731           __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
1732           if (_field>0)
1733             {
1734               for (int i1=0; i1<k1; i1++)
1735                 {
1736                   ////cout << "         i1:" << i1 << endl;
1737                   id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1738                   q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1739                   iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1740                   iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1741                   corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1742                   pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1743
1744
1745                   //1 and 2
1746                   for (int i2=i1+1; i2<k1; i2++)
1747                     {
1748                       ////cout << "         i2:" << i2 << endl;
1749                       id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
1750                       if (id_1!=id_2)
1751                         {
1752                           q_2       = _charge_1[i2];     ////cout << "        q_1:" << q_1 << endl;
1753                           iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1754                           iPt_2     = _iPt_1[i2];        ////cout << "      iPt_1:" << iPt_1 << endl;
1755                           corr_2    = _correction_1[i2]; ////cout << "     corr_1:" << corr_1 << endl;
1756                           pt_2      = _pt_1[i2];         ////cout << "       pt_1:" << pt_1 << endl;
1757
1758                           corr      = corr_1*corr_2;
1759                           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))
1760                             {
1761                               ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
1762                             }
1763                           else // swap particles
1764                             {
1765                               ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
1766                             }
1767
1768                           __n2_12                  += corr;
1769                           __n2_12_vsEtaPhi[ij]     += corr;
1770                           ptpt                     = pt_1*pt_2;
1771                           __s2ptpt_12              += corr*ptpt;
1772                           __s2PtN_12               += corr*pt_1;
1773                           __s2NPt_12               += corr*pt_2;
1774                           __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
1775                           __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
1776                           __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
1777                           __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
1778
1779                           __n2Nw_12                  += 1;
1780                           __s2ptptNw_12              += ptpt;
1781                           __s2PtNNw_12               += pt_1;
1782                           __s2NPtNw_12               += pt_2;
1783
1784                         }
1785                     } //i2
1786                 } //i1
1787             }
1788
1789           else // field<0
1790             {
1791               for (int i1=0; i1<k1; i1++)
1792                 {
1793                   ////cout << "         i1:" << i1 << endl;
1794                   id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1795                   q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1796                   iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1797                   iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1798                   corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1799                   pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1800
1801                   //1 and 2
1802                   for (int i2=i1+1; i2<k1; i2++)
1803                     {
1804                       ////cout << "         i2:" << i2 << endl;
1805                       id_2      = _id_1[i2];              ////cout << "       id_2:" << id_2 << endl;
1806                       if (id_1!=id_2)
1807                         {
1808                           q_2       = _charge_1[i2];     ////cout << "        q_2:" << q_2 << endl;
1809                           iEtaPhi_2 = _iEtaPhi_1[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
1810                           iPt_2     = _iPt_1[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
1811                           corr_2    = _correction_1[i2]; ////cout << "     corr_2:" << corr_2 << endl;
1812                           pt_2      = _pt_1[i2];         ////cout << "       pt_2:" << pt_2 << endl;
1813                           corr      = corr_1*corr_2;
1814                           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))
1815                             {
1816                               ij = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
1817                             }
1818                           else // swap particles
1819                             {
1820                               ij = iEtaPhi_2*_nBins_etaPhi_1 + iEtaPhi_1;   ////cout << " ij:" << ij<< endl;
1821                             }
1822
1823                           __n2_12                  += corr;
1824                           __n2_12_vsEtaPhi[ij]     += corr;
1825                           ptpt                     = pt_1*pt_2;
1826                           __s2ptpt_12              += corr*ptpt;
1827                           __s2PtN_12               += corr*pt_1;
1828                           __s2NPt_12               += corr*pt_2;
1829                           __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
1830                           __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
1831                           __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
1832                           __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
1833
1834                           __n2Nw_12                  += 1;
1835                           __s2ptptNw_12              += ptpt;
1836                           __s2PtNNw_12               += pt_1;
1837                           __s2NPtNw_12               += pt_2;
1838
1839                         }
1840                     } //i2
1841                 } //i1
1842             }
1843         }
1844       else  // filter 1 and 2 are different -- must do all particle pairs...
1845         {
1846           _n1_1_vsM->Fill(centrality,      __n1_1);
1847           _s1pt_1_vsM->Fill(centrality,    __s1pt_1);
1848           _n1Nw_1_vsM->Fill(centrality,    __n1Nw_1);
1849           _s1ptNw_1_vsM->Fill(centrality,  __s1ptNw_1);
1850           _n1_2_vsM->Fill(centrality,      __n1_2);
1851           _s1pt_2_vsM->Fill(centrality,    __s1pt_2);
1852           _n1Nw_2_vsM->Fill(centrality,    __n1Nw_2);
1853           _s1ptNw_2_vsM->Fill(centrality,  __s1ptNw_2);
1854           // reset pair counters
1855           __n2_12   = __s2ptpt_12   = __s2NPt_12    = __s2PtN_12    = 0;
1856           __n2Nw_12 = __s2ptptNw_12 = __s2NPtNw_12  = __s2PtNNw_12  = 0;
1857           for (int i1=0; i1<k1; i1++)
1858             {
1859               ////cout << "         i1:" << i1 << endl;
1860               id_1      = _id_1[i1];           ////cout << "       id_1:" << id_1 << endl;
1861               q_1       = _charge_1[i1];       ////cout << "        q_1:" << q_1 << endl;
1862               iEtaPhi_1 = _iEtaPhi_1[i1];      ////cout << "  iEtaPhi_1:" << iEtaPhi_1 << endl;
1863               iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
1864               corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
1865               pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
1866               px_1      = _px_1[i1];          ////cout << "      px_1:" << px_1 << endl;
1867               py_1      = _py_1[i1];          ////cout << "      py_1:" << py_1 << endl;
1868               pz_1      = _pz_1[i1];          ////cout << "      pz_1:" << pz_1 << endl;
1869               //1 and 2
1870               for (int i2=0; i2<k2; i2++)
1871                 {
1872                   ////cout << "         i2:" << i2 << endl;
1873                   id_2   = _id_2[i2];              ////cout << "       id_2:" << id_2 << endl;
1874                   if (id_1!=id_2)  // exclude auto correlation
1875                     {
1876                       q_2       = _charge_2[i2];     ////cout << "        q_2:" << q_2 << endl;
1877                       iEtaPhi_2 = _iEtaPhi_2[i2];    ////cout << "  iEtaPhi_2:" << iEtaPhi_2 << endl;
1878                       iPt_2     = _iPt_2[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
1879                       corr_2    = _correction_2[i2]; ////cout << "     corr_2:" << corr_2 << endl;
1880                       pt_2      = _pt_2[i2];         ////cout << "       pt_2:" << pt_2 << endl;
1881                       px_2      = _px_2[i2];          ////cout << "      px_2:" << px_2 << endl;
1882                       py_2      = _py_2[i2];          ////cout << "      py_2:" << py_2 << endl;
1883                       pz_2      = _pz_2[i2];          ////cout << "      pz_2:" << pz_2 << endl;
1884
1885
1886                       if (_rejectPairConversion)
1887                         {
1888                           float e1Sq = massElecSq + pt_1*pt_1 + pz_1*pz_1;
1889                           float e2Sq = massElecSq + pt_2*pt_2 + pz_2*pz_2;
1890                           float mInvSq = 2*(massElecSq + sqrt(e1Sq*e2Sq) - px_1*px_2 - py_1*py_2 - pz_1*pz_2 );
1891                           float mInv = sqrt(mInvSq);
1892                           _invMass->Fill(mInv);
1893                         }
1894
1895                       corr      = corr_1*corr_2;
1896                       ij        = iEtaPhi_1*_nBins_etaPhi_1 + iEtaPhi_2;   ////cout << " ij:" << ij<< endl;
1897                       __n2_12                  += corr;
1898                       __n2_12_vsEtaPhi[ij]     += corr;
1899                       ptpt                     = pt_1*pt_2;
1900                       __s2ptpt_12              += corr*ptpt;
1901                       __s2PtN_12               += corr*pt_1;
1902                       __s2NPt_12               += corr*pt_2;
1903                       __s2ptpt_12_vsEtaPhi[ij] += corr*ptpt;
1904                       __s2PtN_12_vsEtaPhi[ij]  += corr*pt_1;
1905                       __s2NPt_12_vsEtaPhi[ij]  += corr*pt_2;
1906                       __n2_12_vsPtPt[iPt_1*_nBins_pt_2 + iPt_2] += corr;
1907                       __n2Nw_12                  += 1;
1908                       __s2ptptNw_12              += ptpt;
1909                       __s2PtNNw_12               += pt_1;
1910                       __s2NPtNw_12               += pt_2;
1911
1912                     }
1913                 } //i2
1914             } //i1
1915         }
1916
1917       _n2_12_vsM->Fill(centrality,     __n2_12);
1918       _s2PtPt_12_vsM->Fill(centrality, __s2ptpt_12);
1919       _s2PtN_12_vsM->Fill(centrality,  __s2NPt_12);
1920       _s2NPt_12_vsM->Fill(centrality,  __s2PtN_12);
1921
1922       _n2Nw_12_vsM->Fill(centrality,     __n2Nw_12);
1923       _s2PtPtNw_12_vsM->Fill(centrality, __s2ptptNw_12);
1924       _s2PtNNw_12_vsM->Fill(centrality,  __s2NPtNw_12);
1925       _s2NPtNw_12_vsM->Fill(centrality,  __s2PtNNw_12);
1926
1927     }
1928
1929
1930   AliInfo("AliTaskDptCorrMC::UserExec()   -----------------Event Done ");
1931   PostData(0,_outputHistoList);
1932
1933 } //End of UserExec
1934
1935
1936
1937
1938 void   AliDptDptInMC::FinishTaskOutput()
1939 {
1940   AliInfo("AliDptDptInMC::FinishTaskOutput() Starting.");
1941   Printf("= 0 ====================================================================");
1942   finalizeHistograms();
1943   AliInfo("= 1 ====================================================================");
1944   PostData(0,_outputHistoList);
1945   AliInfo("= 2 ====================================================================");
1946   AliInfo("AliDptDptInMC::FinishTaskOutput() Done.");
1947 }
1948
1949 void   AliDptDptInMC::Terminate(Option_t* /*option*/)
1950 {
1951   AliInfo("AliDptDptInMC::Terminate() Starting/Done.");
1952 }
1953
1954
1955 //Tools
1956 //===================================================================================================
1957 void  AliDptDptInMC::fillHistoWithArray(TH1 * h, float * array, int size)
1958 {
1959   int i, i1;
1960   float v1, ev1, v2, ev2, sum, esum;
1961   for (i=0, i1=1; i<size; ++i,++i1)
1962     {
1963     v1  = array[i]; ev1 = sqrt(v1);
1964     v2  = h->GetBinContent(i1);
1965     ev2 = h->GetBinError(i1);
1966     sum = v1 + v2;
1967     esum = sqrt(ev1*ev1+ev2*ev2);
1968     h->SetBinContent(i1,sum);
1969     h->SetBinError(i1,esum);
1970     }
1971 }
1972
1973 void  AliDptDptInMC::fillHistoWithArray(TH2 * h, float * array, int size1, int size2)
1974 {
1975   int i, i1;
1976   int j, j1;
1977   float v1, ev1, v2, ev2, sum, esum;
1978   for (i=0, i1=1; i<size1; ++i,++i1)
1979     {
1980     for (j=0, j1=1; j<size2; ++j,++j1)
1981       {
1982       v1  = array[i*size2+j]; ev1 = sqrt(v1);
1983       v2  = h->GetBinContent(i1,j1);
1984       ev2 = h->GetBinError(i1,j1);
1985       sum = v1 + v2;
1986       esum = sqrt(ev1*ev1+ev2*ev2);
1987       h->SetBinContent(i1,j1,sum);
1988       h->SetBinError(i1,j1,esum);
1989       }
1990     }
1991 }
1992
1993 void  AliDptDptInMC::fillHistoWithArray(TH3 * h, float * array, int size1, int size2, int size3)
1994 {
1995   int i, i1;
1996   int j, j1;
1997   int k, k1;
1998   float v1, ev1, v2, ev2, sum, esum;
1999   int size23 = size2*size3;
2000   for (i=0, i1=1; i<size1; ++i,++i1)
2001     {
2002     for (j=0, j1=1; j<size2; ++j,++j1)
2003       {
2004       for (k=0, k1=1; k<size3; ++k,++k1)
2005         {
2006         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
2007         v2  = h->GetBinContent(i1,j1,k1);
2008         ev2 = h->GetBinError(i1,j1,k1);
2009         sum = v1 + v2;
2010         esum = sqrt(ev1*ev1+ev2*ev2);
2011         h->SetBinContent(i1,j1,k1,sum);
2012         h->SetBinError(i1,j1,k1,esum);
2013         }
2014       }
2015     }
2016 }
2017
2018 void  AliDptDptInMC::fillHistoWithArray(TH1 * h, double * array, int size)
2019 {
2020   int i, i1;
2021   double v1, ev1, v2, ev2, sum, esum;
2022   for (i=0, i1=1; i<size; ++i,++i1)
2023     {
2024     v1  = array[i]; ev1 = sqrt(v1);
2025     v2  = h->GetBinContent(i1);
2026     ev2 = h->GetBinError(i1);
2027     sum = v1 + v2;
2028     esum = sqrt(ev1*ev1+ev2*ev2);
2029     h->SetBinContent(i1,sum);
2030     h->SetBinError(i1,esum);
2031     }
2032 }
2033
2034 void  AliDptDptInMC::fillHistoWithArray(TH2 * h, double * array, int size1, int size2)
2035 {
2036   int i, i1;
2037   int j, j1;
2038   double v1, ev1, v2, ev2, sum, esum;
2039   for (i=0, i1=1; i<size1; ++i,++i1)
2040     {
2041     for (j=0, j1=1; j<size2; ++j,++j1)
2042       {
2043       v1  = array[i*size2+j]; ev1 = sqrt(v1);
2044       v2  = h->GetBinContent(i1,j1);
2045       ev2 = h->GetBinError(i1,j1);
2046       sum = v1 + v2;
2047       esum = sqrt(ev1*ev1+ev2*ev2);
2048       h->SetBinContent(i1,j1,sum);
2049       h->SetBinError(i1,j1,esum);
2050       }
2051     }
2052 }
2053
2054 void  AliDptDptInMC::fillHistoWithArray(TH3 * h, double * array, int size1, int size2, int size3)
2055 {
2056   int i, i1;
2057   int j, j1;
2058   int k, k1;
2059   double v1, ev1, v2, ev2, sum, esum;
2060   int size23 = size2*size3;
2061   for (i=0, i1=1; i<size1; ++i,++i1)
2062     {
2063     for (j=0, j1=1; j<size2; ++j,++j1)
2064       {
2065       for (k=0, k1=1; k<size3; ++k,++k1)
2066         {
2067         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
2068         v2  = h->GetBinContent(i1,j1,k1);
2069         ev2 = h->GetBinError(i1,j1,k1);
2070         sum = v1 + v2;
2071         esum = sqrt(ev1*ev1+ev2*ev2);
2072         h->SetBinContent(i1,j1,k1,sum);
2073         h->SetBinError(i1,j1,k1,esum);
2074         }
2075       }
2076     }
2077 }
2078
2079 //________________________________________________________________________
2080 double *  AliDptDptInMC::getDoubleArray(int size, double v)
2081 {
2082   /// Allocate an array of type double with n values
2083   /// Initialize the array to the given value
2084   double * array = new double [size];
2085   for (int i=0;i<size;++i) array[i]=v;
2086   return array;
2087 }
2088
2089 //________________________________________________________________________
2090 float *  AliDptDptInMC::getFloatArray(int size, float v)
2091 {
2092   /// Allocate an array of type float with n values
2093   /// Initialize the array to the given value
2094   float * array = new float [size];
2095   for (int i=0;i<size;++i) array[i]=v;
2096   return array;
2097 }
2098
2099
2100 //________________________________________________________________________
2101 TH1D * AliDptDptInMC::createHisto1D(const TString &  name, const TString &  title,
2102                                                       int n, double xMin, double xMax,
2103                                                       const TString &  xTitle, const TString &  yTitle)
2104 {
2105   //CreateHisto new 1D historgram
2106   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
2107   TH1D * h = new TH1D(name,title,n,xMin,xMax);
2108   h->GetXaxis()->SetTitle(xTitle);
2109   h->GetYaxis()->SetTitle(yTitle);
2110   addToList(h);
2111   return h;
2112 }
2113
2114
2115 //________________________________________________________________________
2116 TH1D * AliDptDptInMC::createHisto1D(const TString &  name, const TString &  title,
2117                                                       int n, double * bins,
2118                                                       const TString &  xTitle, const TString &  yTitle)
2119 {
2120   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
2121   TH1D * h = new TH1D(name,title,n,bins);
2122   h->GetXaxis()->SetTitle(xTitle);
2123   h->GetYaxis()->SetTitle(yTitle);
2124   addToList(h);
2125   return h;
2126 }
2127
2128
2129 //________________________________________________________________________
2130 TH2D * AliDptDptInMC::createHisto2D(const TString &  name, const TString &  title,
2131                                                       int nx, double xMin, double xMax, int ny, double yMin, double yMax,
2132                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2133 {
2134   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));
2135   TH2D * h = new TH2D(name,title,nx,xMin,xMax,ny,yMin,yMax);
2136   h->GetXaxis()->SetTitle(xTitle);
2137   h->GetYaxis()->SetTitle(yTitle);
2138   h->GetZaxis()->SetTitle(zTitle);
2139   addToList(h);
2140   return h;
2141 }
2142
2143 //________________________________________________________________________
2144 TH2D * AliDptDptInMC::createHisto2D(const TString &  name, const TString &  title,
2145                                                       int nx, double* xbins, int ny, double yMin, double yMax,
2146                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2147 {
2148   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
2149   TH2D * h;
2150   h = new TH2D(name,title,nx,xbins,ny,yMin,yMax);
2151   h->GetXaxis()->SetTitle(xTitle);
2152   h->GetYaxis()->SetTitle(yTitle);
2153   h->GetZaxis()->SetTitle(zTitle);
2154   addToList(h);
2155   return h;
2156 }
2157
2158 //// F /////
2159 //________________________________________________________________________
2160 TH1F * AliDptDptInMC::createHisto1F(const TString &  name, const TString &  title,
2161                                                         int n, double xMin, double xMax,
2162                                                         const TString &  xTitle, const TString &  yTitle)
2163 {
2164   //CreateHisto new 1D historgram
2165   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
2166   TH1F * h = new TH1F(name,title,n,xMin,xMax);
2167   h->GetXaxis()->SetTitle(xTitle);
2168   h->GetYaxis()->SetTitle(yTitle);
2169   addToList(h);
2170   return h;
2171 }
2172
2173
2174 //________________________________________________________________________
2175 TH1F * AliDptDptInMC::createHisto1F(const TString &  name, const TString &  title,
2176                                                         int n, double * bins,
2177                                                         const TString &  xTitle, const TString &  yTitle)
2178 {
2179   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
2180   TH1F * h = new TH1F(name,title,n,bins);
2181   h->GetXaxis()->SetTitle(xTitle);
2182   h->GetYaxis()->SetTitle(yTitle);
2183   addToList(h);
2184   return h;
2185 }
2186
2187
2188 //________________________________________________________________________
2189 TH2F * AliDptDptInMC::createHisto2F(const TString &  name, const TString &  title,
2190                                                         int nx, double xMin, double xMax, int ny, double yMin, double yMax,
2191                                                         const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2192 {
2193   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));
2194   TH2F * h = new TH2F(name,title,nx,xMin,xMax,ny,yMin,yMax);
2195   h->GetXaxis()->SetTitle(xTitle);
2196   h->GetYaxis()->SetTitle(yTitle);
2197   h->GetZaxis()->SetTitle(zTitle);
2198   addToList(h);
2199   return h;
2200 }
2201
2202 //________________________________________________________________________
2203 TH2F * AliDptDptInMC::createHisto2F(const TString &  name, const TString &  title,
2204                                                         int nx, double* xbins, int ny, double yMin, double yMax,
2205                                                         const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2206 {
2207   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
2208   TH2F * h;
2209   h = new TH2F(name,title,nx,xbins,ny,yMin,yMax);
2210   h->GetXaxis()->SetTitle(xTitle);
2211   h->GetYaxis()->SetTitle(yTitle);
2212   h->GetZaxis()->SetTitle(zTitle);
2213   addToList(h);
2214   return h;
2215 }
2216
2217
2218 //________________________________________________________________________
2219 TH3F * AliDptDptInMC::createHisto3F(const TString &  name, const TString &  title,
2220                                                       int nx, double xMin, double xMax,
2221                                                       int ny, double yMin, double yMax,
2222                                                       int nz, double zMin, double zMax,
2223                                                       const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
2224 {
2225   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));
2226   TH3F * h = new TH3F(name,title,nx,xMin,xMax,ny,yMin,yMax,nz,zMin,zMax);
2227   h->GetXaxis()->SetTitle(xTitle);
2228   h->GetYaxis()->SetTitle(yTitle);
2229   h->GetZaxis()->SetTitle(zTitle);
2230   addToList(h);
2231   return h;
2232 }
2233
2234
2235 //________________________________________________________________________
2236 TProfile * AliDptDptInMC::createProfile(const TString & name, const TString & description,
2237                                                             int nx,double xMin,double xMax,
2238                                                             const TString &  xTitle, const TString &  yTitle)
2239 {
2240   AliInfo(Form("createHisto 1D profile %s   nBins: %d  xMin: %f10.4 xMax: %f10.4",name.Data(),nx,xMin,xMax));
2241   TProfile * h = new TProfile(name,description,nx,xMin,xMax);
2242   h->GetXaxis()->SetTitle(xTitle);
2243   h->GetYaxis()->SetTitle(yTitle);
2244   addToList(h);
2245   return h;
2246 }
2247
2248 //________________________________________________________________________
2249 TProfile * AliDptDptInMC::createProfile(const TString &  name,const TString &  description,
2250                                                             int nx,  double* bins,
2251                                                             const TString &  xTitle, const TString &  yTitle)
2252 {
2253   AliInfo(Form("createHisto 1D profile %s  with %d non uniform bins",name.Data(),nx));
2254   TProfile * h = new TProfile(name,description,nx,bins);
2255   h->GetXaxis()->SetTitle(xTitle);
2256   h->GetYaxis()->SetTitle(yTitle);
2257   addToList(h);
2258   return h;
2259 }
2260
2261
2262 void   AliDptDptInMC::addToList(TH1 *h)
2263 {
2264   if (_outputHistoList)
2265     {
2266     _outputHistoList->Add(h);
2267     }
2268   else
2269     AliInfo("addToList(TH1 *h) _outputHistoList is null!!!!! Should abort ship");
2270
2271 }
2272
2273
2274