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