PWGCF/Correlations/DPhi/AliDptDptInMC.cxx

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