PWGCF/Correlations/DPhi/AliAnalysisTask3PCorrelations.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 <TROOT.h>
  19 #include <TChain.h>
  20 #include <TFile.h>
  21 #include <TList.h>
  22 #include <TMath.h>
  23 #include <TTree.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TH3F.h>
  27 #include <TProfile.h>
  28 #include <TH1D.h>
  29 #include <TH2D.h>
  30 #include <TH3D.h>
  31 #include <TRandom.h>
  32 #include "AliAnalysisManager.h"
  33
  34 #include "AliAODHandler.h"
  35 #include "AliAODInputHandler.h"
  36 //#include "AliAODMCParticle.h"
  37 #include "AliInputEventHandler.h"
  38 #include "AliLog.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliMultiplicity.h"
  42 #include "AliCentrality.h"
  43 #include "AliAnalysisTask3PCorrelations.h"
  44
  45 using std::cout;
  46 using std::endl;
  47
  48 ClassImp(AliAnalysisTask3PCorrelations)
  49
  50 AliAnalysisTask3PCorrelations::AliAnalysisTask3PCorrelations()
  51 : AliAnalysisTaskSE(),
  52 fAODEvent(0),
  53 fESDEvent(0),             //! ESD Event
  54 fInputHandler(0),
  55 _outputHistoList(0),
  56 _twoPi         ( 2.0 * 3.1415927),
  57 _eventCount    ( 0),
  58 _debugLevel    ( 0),
  59 _singlesOnly   ( 0),
  60 _useWeights    ( 0),
  61 _rejectPileup  ( 1),
  62 _rejectPairConversion ( 0),
  63 _vertexZMin           ( -10),
  64 _vertexZMax           (  10),
  65 _vertexXYMin          ( -10),
  66 _vertexXYMax          (  10),
  67 _centralityMethod     (  4),
  68 _centralityMin        (  0.),
  69 _centralityMax        (  0.),
  70 _dcaZMin              ( -3),
  71 _dcaZMax              (  3.),
  72 _dcaXYMin             ( -3.),
  73 _dcaXYMax             (  3.),
  74 _dedxMin              ( 0),
  75 _dedxMax              ( 100000),
  76 _nClusterMin          ( 70),
  77 _trackFilterBit       ( 128),
  78 _field    ( 1.),
  79 _nTracks  ( 0 ),
  80 _mult0    ( 0 ),
  81 _mult1    ( 0 ),
  82 _mult2    ( 0 ),
  83 _mult3    ( 0 ),
  84 _mult4    ( 0 ),
  85 _mult5    ( 0 ),
  86 _mult6    ( 0 ),
  87 arraySize ( 2000),
  88 _id_1(0),
  89 _charge_1(0),
  90 _iPhi_1(0),
  91 _pt_1(0),
  92 _px_1(0),
  93 _py_1(0),
  94 _pz_1(0),
  95 _correction_1(0),
  96 _dedx_1(0),
  97 _id_2(0),
  98 _charge_2(0),
  99 _iPhi_2(0),
 100 _pt_2(0),
 101 _px_2(0),
 102 _py_2(0),
 103 _pz_2(0),
 104 _correction_2(0),
 105 _dedx_2(0),
 106 _id_3(0),
 107 _charge_3(0),
 108 _iPhi_3(0),
 109 _pt_3(0),
 110 _px_3(0),
 111 _py_3(0),
 112 _pz_3(0),
 113 _correction_3(0),
 114 _dedx_3(0),
 115
 116 _correctionWeight_1P(0),
 117 _correctionWeight_1M(0),
 118 _correctionWeight_2P(0),
 119 _correctionWeight_2M(0),
 120 _correctionWeight_3P(0),
 121 _correctionWeight_3M(0),
 122 _nBins_M0(500),       _min_M0(0),        _max_M0(10000),          _width_M0(20),
 123 _nBins_M1(500),       _min_M1(0),        _max_M1(10000),          _width_M1(20),
 124 _nBins_M2(500),       _min_M2(0),        _max_M2(10000),          _width_M2(20),
 125 _nBins_M3(500),       _min_M3(0),        _max_M3(10000),          _width_M3(20),
 126 _nBins_M4(100),       _min_M4(0),        _max_M4(1),              _width_M4(0.01),
 127 _nBins_M5(100),       _min_M5(0),        _max_M5(1),              _width_M5(0.01),
 128 _nBins_M6(100),       _min_M6(0),        _max_M6(1),              _width_M6(0.01),
 129 _nBins_vsM(0),        _min_vsM(0),       _max_vsM(0),             _width_vsM(0),
 130 _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ(0.5),
 131
 132 _nBins_pt_1(8),       _min_pt_1(2.),    _max_pt_1(10.0),          _width_pt_1(1.0),
 133 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
 134 _nBins_eta_1(1),      _min_eta_1(-0.25), _max_eta_1(0.25),        _width_eta_1(0.5),
 135 _nBins_etaPhi_1(0),
 136 _nBins_etaPhiPt_1(0),
 137 _nBins_zEtaPhiPt_1(0),
 138
 139 _nBins_pt_2(5),       _min_pt_2(0.5),     _max_pt_2(1.0),          _width_pt_2(0.1),
 140 _nBins_phi_2(72),     _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
 141 _nBins_eta_2(1),      _min_eta_2(-1),     _max_eta_2(-0.7),        _width_eta_2(0.3),
 142 _nBins_etaPhi_2(0),
 143 _nBins_etaPhiPt_2(0),
 144 _nBins_zEtaPhiPt_2(0),
 145
 146 _nBins_pt_3(5),       _min_pt_3(0.5),     _max_pt_3(1.0),          _width_pt_3(0.1),
 147 _nBins_phi_3(72),     _min_phi_3(0),      _max_phi_3(2.*3.1415927),_width_phi_3(2.*3.1415927/72),
 148 _nBins_eta_3(1),      _min_eta_3(0.7),    _max_eta_3(1.),          _width_eta_3(0.3),
 149 _nBins_etaPhi_3(0),
 150 _nBins_etaPhiPt_3(0),
 151 _nBins_zEtaPhiPt_3(0),
 152
 153 _nBins_phi_12(0),
 154 _nBins_phi_13(0),
 155 _nBins_phi_23(0),
 156 _nBins_phi_123(0),
 157
 158
 159 _nBins_etaPhi_12(0),
 160 _nBins_etaPhi_13(0),
 161 _nBins_etaPhi_23(0),
 162 _nBins_etaPhi_123(0),
 163
 164 __n1_1(0),
 165 __n1_2(0),
 166 __n1_3(0),
 167 __n2_12(0),
 168 __n2_13(0),
 169 __n2_23(0),
 170 __n3_123(0),
 171 __n1Nw_1(0),
 172 __n1Nw_2(0),
 173 __n1Nw_3(0),
 174 __n2Nw_12(0),
 175 __n2Nw_13(0),
 176 __n2Nw_23(0),
 177 __n3Nw_123(0),
 178
 179 __n1_1_vsPt(0),
 180 __n1_1_vsPhi(0),
 181 __n1_1P_vsZEtaPhiPt(0),
 182 __n1_1M_vsZEtaPhiPt(0),
 183
 184 __n1_2_vsPt(0),
 185 __n1_2_vsPhi(0),
 186 __n1_2P_vsZEtaPhiPt(0),
 187 __n1_2M_vsZEtaPhiPt(0),
 188
 189 __n1_3_vsPt(0),
 190 __n1_3_vsPhi(0),
 191 __n1_3P_vsZEtaPhiPt(0),
 192 __n1_3M_vsZEtaPhiPt(0),
 193
 194 __n2_12_vsPhi(0),
 195 __n2_13_vsPhi(0),
 196 __n2_23_vsPhi(0),
 197 __n3_123_vsPhi(0),
 198
 199 _weight_1P      ( 0    ),
 200 _weight_1M      ( 0    ),
 201 _weight_2P      ( 0    ),
 202 _weight_2M      ( 0    ),
 203 _weight_3P      ( 0    ),
 204 _weight_3M      ( 0    ),
 205 _eventAccounting ( 0),
 206 _m0 ( 0),
 207 _m1 ( 0),
 208 _m2 ( 0),
 209 _m3 ( 0),
 210 _m4 ( 0),
 211 _m5 ( 0),
 212 _m6 ( 0),
 213 _vertexZ        ( 0),
 214
 215 _n1_1_vsPhi     ( 0),
 216 _n1_1_vsM          ( 0), // w/ weight
 217 _n1Nw_1_vsM        ( 0), // w/o weight
 218
 219 _n1_1_vsPt      ( 0),
 220 _n1_1P_vsZVsEtaVsPhiVsPt ( 0),
 221 _n1_1M_vsZVsEtaVsPhiVsPt ( 0),
 222 _dedxVsP_1         ( 0),
 223 _corrDedxVsP_1     ( 0),
 224 _betaVsP_1         ( 0),
 225
 226 _n1_2_vsPhi        ( 0),
 227 _n1_2_vsM          ( 0),
 228 _n1Nw_2_vsM        ( 0),
 229 _n1_2_vsPt         ( 0),
 230 _n1_2P_vsZVsEtaVsPhiVsPt ( 0),
 231 _n1_2M_vsZVsEtaVsPhiVsPt ( 0),
 232 _dedxVsP_2         ( 0),
 233 _corrDedxVsP_2     ( 0),
 234 _betaVsP_2         ( 0),
 235
 236 _n1_3_vsPhi        ( 0),
 237 _n1_3_vsM          ( 0),
 238 _n1Nw_3_vsM        ( 0),
 239 _n1_3_vsPt         ( 0),
 240 _n1_3P_vsZVsEtaVsPhiVsPt ( 0),
 241 _n1_3M_vsZVsEtaVsPhiVsPt ( 0),
 242 _dedxVsP_3         ( 0),
 243 _corrDedxVsP_3     ( 0),
 244 _betaVsP_3         ( 0),
 245
 246 _n2_12_vsPhi    ( 0),
 247 _n2_13_vsPhi    ( 0),
 248 _n2_23_vsPhi    ( 0),
 249 _n3_123_vsPhi   ( 0),
 250
 251 _n2_12_vsM         ( 0),
 252 _n2_13_vsM         ( 0),
 253 _n2_23_vsM         ( 0),
 254 _n3_123_vsM        ( 0),
 255
 256 _n2Nw_12_vsM       ( 0),
 257 _n2Nw_13_vsM       ( 0),
 258 _n2Nw_23_vsM       ( 0),
 259 _n3Nw_123_vsM      ( 0),
 260
 261 n1Name("NA"),
 262 n2Name("NA"),
 263 n3Name("NA"),
 264 n1NwName("NA"),
 265 n2NwName("NA"),
 266 n3NwName("NA"),
 267 n1n1Name("NA"),
 268 n1n1n1Name("NA"),
 269 n2n1Name("NA"),
 270 r1Name("NA"),
 271 r2Name("NA"),
 272 r3Name("NA"),
 273 r2r1Name("NA"),
 274 c2Name("NA"),
 275 c3Name("NA"),
 276 d3Name("NA"),
 277 p3Name("NA"),
 278
 279 intR2Name("NA"),
 280 binCorrName("NA"),
 281 intBinCorrName("NA"),
 282
 283 countsName("NA"),
 284 part_1_Name("NA"),
 285 part_1P_Name("NA"),
 286 part_1M_Name("NA"),
 287 part_2_Name("NA"),
 288 part_2P_Name("NA"),
 289 part_2M_Name("NA"),
 290 part_3_Name("NA"),
 291 part_3P_Name("NA"),
 292 part_3M_Name("NA"),
 293 pair_12_Name("NA"),
 294 pair_13_Name("NA"),
 295 pair_23_Name("NA"),
 296 triplet_Name("NA"),
 297
 298 _title_counts("NA"),
 299 _title_m0("NA"),
 300 _title_m1("NA"),
 301 _title_m2("NA"),
 302 _title_m3("NA"),
 303 _title_m4("NA"),
 304 _title_m5("NA"),
 305 _title_m6("NA"),
 306
 307 _title_phi_1("NA"),
 308 _title_phi_2("NA"),
 309 _title_phi_3("NA"),
 310 _title_phi_12("NA"),
 311 _title_phi_13("NA"),
 312 _title_phi_23("NA"),
 313 _title_phi_123("NA"),
 314
 315 vsZ("NA"),
 316 vsM("NA"),
 317 vsPt("NA"),
 318 vsPhi("NA"),
 319 vsEta("NA"),
 320 vsEtaPhi("NA")
 321 {
 322   printf("Default constructor called \n");
 323
 324   printf("passed \n ");
 325
 326 }
 327
 328 AliAnalysisTask3PCorrelations::AliAnalysisTask3PCorrelations(const TString & name)
 329 : AliAnalysisTaskSE(name),
 330 fAODEvent(0),
 331 fESDEvent(0),             //! ESD Event
 332 fInputHandler(0),
 333 _outputHistoList(0),
 334 _twoPi         ( 2.0 * 3.1415927),
 335 _eventCount    ( 0),
 336 _debugLevel    ( 0),
 337 _singlesOnly   ( 0),
 338 _useWeights    ( 0),
 339 _rejectPileup  ( 1),
 340 _rejectPairConversion ( 0),
 341 _vertexZMin           ( -10),
 342 _vertexZMax           (  10),
 343 _vertexXYMin          ( -10),
 344 _vertexXYMax          (  10),
 345 _centralityMethod     (  4),
 346 _centralityMin        (  0.),
 347 _centralityMax        (  0.),
 348 _dcaZMin              ( -3),
 349 _dcaZMax              (  3.),
 350 _dcaXYMin             ( -3.),
 351 _dcaXYMax             (  3.),
 352 _dedxMin              ( 0),
 353 _dedxMax              ( 100000),
 354 _nClusterMin          ( 70),
 355 _trackFilterBit       ( 128),
 356 _field    ( 1.),
 357 _nTracks  ( 0 ),
 358 _mult0    ( 0 ),
 359 _mult1    ( 0 ),
 360 _mult2    ( 0 ),
 361 _mult3    ( 0 ),
 362 _mult4    ( 0 ),
 363 _mult5    ( 0 ),
 364 _mult6    ( 0 ),
 365 arraySize ( 2000),
 366 _id_1(0),
 367 _charge_1(0),
 368 _iPhi_1(0),
 369 _pt_1(0),
 370 _px_1(0),
 371 _py_1(0),
 372 _pz_1(0),
 373 _correction_1(0),
 374 _dedx_1(0),
 375 _id_2(0),
 376 _charge_2(0),
 377 _iPhi_2(0),
 378 _pt_2(0),
 379 _px_2(0),
 380 _py_2(0),
 381 _pz_2(0),
 382 _correction_2(0),
 383 _dedx_2(0),
 384 _id_3(0),
 385 _charge_3(0),
 386 _iPhi_3(0),
 387 _pt_3(0),
 388 _px_3(0),
 389 _py_3(0),
 390 _pz_3(0),
 391 _correction_3(0),
 392 _dedx_3(0),
 393
 394 _correctionWeight_1P(0),
 395 _correctionWeight_1M(0),
 396 _correctionWeight_2P(0),
 397 _correctionWeight_2M(0),
 398 _correctionWeight_3P(0),
 399 _correctionWeight_3M(0),
 400 _nBins_M0(500),       _min_M0(0),        _max_M0(10000),          _width_M0(20),
 401 _nBins_M1(500),       _min_M1(0),        _max_M1(10000),          _width_M1(20),
 402 _nBins_M2(500),       _min_M2(0),        _max_M2(10000),          _width_M2(20),
 403 _nBins_M3(500),       _min_M3(0),        _max_M3(10000),          _width_M3(20),
 404 _nBins_M4(100),       _min_M4(0),        _max_M4(1),              _width_M4(0.01),
 405 _nBins_M5(100),       _min_M5(0),        _max_M5(1),              _width_M5(0.01),
 406 _nBins_M6(100),       _min_M6(0),        _max_M6(1),              _width_M6(0.01),
 407 _nBins_vsM(0),        _min_vsM(0),       _max_vsM(0),             _width_vsM(0),
 408
 409 _nBins_vertexZ(40),   _min_vertexZ(-10), _max_vertexZ(10),        _width_vertexZ(0.5),
 410
 411 _nBins_pt_1(8),       _min_pt_1(2.),    _max_pt_1(10.0),          _width_pt_1(1.0),
 412 _nBins_phi_1(72),     _min_phi_1(0),     _max_phi_1(2.*3.1415927),_width_phi_1(2.*3.1415927/72.),
 413 _nBins_eta_1(1),      _min_eta_1(-0.25), _max_eta_1(0.25),        _width_eta_1(0.5),
 414 _nBins_etaPhi_1(0),
 415 _nBins_etaPhiPt_1(0),
 416 _nBins_zEtaPhiPt_1(0),
 417
 418 _nBins_pt_2(5),       _min_pt_2(0.5),     _max_pt_2(1.0),          _width_pt_2(0.1),
 419 _nBins_phi_2(72),     _min_phi_2(0),      _max_phi_2(2.*3.1415927),_width_phi_2(2.*3.1415927/72),
 420 _nBins_eta_2(1),      _min_eta_2(-1),     _max_eta_2(-0.7),        _width_eta_2(0.3),
 421 _nBins_etaPhi_2(0),
 422 _nBins_etaPhiPt_2(0),
 423 _nBins_zEtaPhiPt_2(0),
 424
 425 _nBins_pt_3(5),       _min_pt_3(0.5),     _max_pt_3(1.0),          _width_pt_3(0.1),
 426 _nBins_phi_3(72),     _min_phi_3(0),      _max_phi_3(2.*3.1415927),_width_phi_3(2.*3.1415927/72),
 427 _nBins_eta_3(1),      _min_eta_3(0.7),    _max_eta_3(1.),          _width_eta_3(0.3),
 428 _nBins_etaPhi_3(0),
 429 _nBins_etaPhiPt_3(0),
 430 _nBins_zEtaPhiPt_3(0),
 431
 432 _nBins_phi_12(0),
 433 _nBins_phi_13(0),
 434 _nBins_phi_23(0),
 435 _nBins_phi_123(0),
 436
 437 _nBins_etaPhi_12(0),
 438 _nBins_etaPhi_13(0),
 439 _nBins_etaPhi_23(0),
 440 _nBins_etaPhi_123(0),
 441
 442 __n1_1(0),
 443 __n1_2(0),
 444 __n1_3(0),
 445 __n2_12(0),
 446 __n2_13(0),
 447 __n2_23(0),
 448 __n3_123(0),
 449 __n1Nw_1(0),
 450 __n1Nw_2(0),
 451 __n1Nw_3(0),
 452 __n2Nw_12(0),
 453 __n2Nw_13(0),
 454 __n2Nw_23(0),
 455 __n3Nw_123(0),
 456
 457 __n1_1_vsPt(0),
 458 __n1_1_vsPhi(0),
 459 __n1_1P_vsZEtaPhiPt(0),
 460 __n1_1M_vsZEtaPhiPt(0),
 461
 462 __n1_2_vsPt(0),
 463 __n1_2_vsPhi(0),
 464 __n1_2P_vsZEtaPhiPt(0),
 465 __n1_2M_vsZEtaPhiPt(0),
 466
 467 __n1_3_vsPt(0),
 468 __n1_3_vsPhi(0),
 469 __n1_3P_vsZEtaPhiPt(0),
 470 __n1_3M_vsZEtaPhiPt(0),
 471
 472 __n2_12_vsPhi(0),
 473 __n2_13_vsPhi(0),
 474 __n2_23_vsPhi(0),
 475 __n3_123_vsPhi(0),
 476
 477 _weight_1P      ( 0    ),
 478 _weight_1M      ( 0    ),
 479 _weight_2P      ( 0    ),
 480 _weight_2M      ( 0    ),
 481 _weight_3P      ( 0    ),
 482 _weight_3M      ( 0    ),
 483 _eventAccounting ( 0),
 484 _m0 ( 0),
 485 _m1 ( 0),
 486 _m2 ( 0),
 487 _m3 ( 0),
 488 _m4 ( 0),
 489 _m5 ( 0),
 490 _m6 ( 0),
 491 _vertexZ        ( 0),
 492
 493 _n1_1_vsPhi     ( 0),
 494 _n1_1_vsM          ( 0), // w/ weight
 495 _n1Nw_1_vsM        ( 0), // w/o weight
 496
 497 _n1_1_vsPt      ( 0),
 498 _n1_1P_vsZVsEtaVsPhiVsPt ( 0),
 499 _n1_1M_vsZVsEtaVsPhiVsPt ( 0),
 500 _dedxVsP_1         ( 0),
 501 _corrDedxVsP_1     ( 0),
 502 _betaVsP_1         ( 0),
 503
 504 _n1_2_vsPhi        ( 0),
 505 _n1_2_vsM          ( 0),
 506 _n1Nw_2_vsM        ( 0),
 507 _n1_2_vsPt         ( 0),
 508 _n1_2P_vsZVsEtaVsPhiVsPt ( 0),
 509 _n1_2M_vsZVsEtaVsPhiVsPt ( 0),
 510 _dedxVsP_2         ( 0),
 511 _corrDedxVsP_2     ( 0),
 512 _betaVsP_2         ( 0),
 513
 514 _n1_3_vsPhi        ( 0),
 515 _n1_3_vsM          ( 0),
 516 _n1Nw_3_vsM        ( 0),
 517 _n1_3_vsPt         ( 0),
 518 _n1_3P_vsZVsEtaVsPhiVsPt ( 0),
 519 _n1_3M_vsZVsEtaVsPhiVsPt ( 0),
 520 _dedxVsP_3         ( 0),
 521 _corrDedxVsP_3     ( 0),
 522 _betaVsP_3         ( 0),
 523
 524 _n2_12_vsPhi    ( 0),
 525 _n2_13_vsPhi    ( 0),
 526 _n2_23_vsPhi    ( 0),
 527 _n3_123_vsPhi   ( 0),
 528
 529 _n2_12_vsM         ( 0),
 530 _n2_13_vsM         ( 0),
 531 _n2_23_vsM         ( 0),
 532 _n3_123_vsM        ( 0),
 533
 534 _n2Nw_12_vsM       ( 0),
 535 _n2Nw_13_vsM       ( 0),
 536 _n2Nw_23_vsM       ( 0),
 537 _n3Nw_123_vsM      ( 0),
 538
 539 n1Name("NA"),
 540 n2Name("NA"),
 541 n3Name("NA"),
 542 n1NwName("NA"),
 543 n2NwName("NA"),
 544 n3NwName("NA"),
 545 n1n1Name("NA"),
 546 n1n1n1Name("NA"),
 547 n2n1Name("NA"),
 548 r1Name("NA"),
 549 r2Name("NA"),
 550 r3Name("NA"),
 551 r2r1Name("NA"),
 552 c2Name("NA"),
 553 c3Name("NA"),
 554 d3Name("NA"),
 555 p3Name("NA"),
 556
 557 intR2Name("NA"),
 558 binCorrName("NA"),
 559 intBinCorrName("NA"),
 560
 561 countsName("NA"),
 562 part_1_Name("NA"),
 563 part_1P_Name("NA"),
 564 part_1M_Name("NA"),
 565 part_2_Name("NA"),
 566 part_2P_Name("NA"),
 567 part_2M_Name("NA"),
 568 part_3_Name("NA"),
 569 part_3P_Name("NA"),
 570 part_3M_Name("NA"),
 571 pair_12_Name("NA"),
 572 pair_13_Name("NA"),
 573 pair_23_Name("NA"),
 574 triplet_Name("NA"),
 575
 576 _title_counts("NA"),
 577 _title_m0("NA"),
 578 _title_m1("NA"),
 579 _title_m2("NA"),
 580 _title_m3("NA"),
 581 _title_m4("NA"),
 582 _title_m5("NA"),
 583 _title_m6("NA"),
 584
 585 _title_phi_1("NA"),
 586 _title_phi_2("NA"),
 587 _title_phi_3("NA"),
 588
 589 _title_phi_12("NA"),
 590 _title_phi_13("NA"),
 591 _title_phi_23("NA"),
 592 _title_phi_123("NA"),
 593
 594 vsZ("NA"),
 595 vsM("NA"),
 596 vsPt("NA"),
 597 vsPhi("NA"),
 598 vsEta("NA"),
 599 vsEtaPhi("NA")
 600 {
 601   printf("2nd constructor called ");
 602
 603   DefineOutput(0, TList::Class());
 604
 605   printf("passed  ");
 606
 607 }
 608
 609 AliAnalysisTask3PCorrelations::~AliAnalysisTask3PCorrelations()
 610 {
 611   // no ops --- at least for now....
 612 }
 613
 614 void AliAnalysisTask3PCorrelations::UserCreateOutputObjects()
 615 {
 616   cout<< "AliAnalysisTask3PCorrelations::CreateOutputObjects() Starting " << endl;
 617   OpenFile(0);
 618   _outputHistoList = new TList();
 619   _outputHistoList->SetOwner();
 620
 621   //if (_useWeights) DefineInput(2, TList::Class());
 622
 623   //Setup the parameters of the histograms
 624   _nBins_M0 = 500; _min_M0   = 0.;    _max_M0    = 5000.;  _width_M0 = (_max_M0-_min_M0)/_nBins_M0;
 625   _nBins_M1 = 500; _min_M1   = 0.;    _max_M1    = 5000.;  _width_M1 = (_max_M1-_min_M1)/_nBins_M1;
 626   _nBins_M2 = 500; _min_M2   = 0.;    _max_M2    = 5000.;  _width_M2 = (_max_M2-_min_M2)/_nBins_M2;
 627   _nBins_M3 = 500; _min_M3   = 0.;    _max_M3    = 5000.;  _width_M3 = (_max_M3-_min_M3)/_nBins_M3;
 628   _nBins_M4 = 100; _min_M4   = 0.;    _max_M4    = 100.;   _width_M4 = (_max_M4-_min_M4)/_nBins_M4;
 629   _nBins_M5 = 100; _min_M5   = 0.;    _max_M5    = 100.;   _width_M5 = (_max_M5-_min_M5)/_nBins_M5;
 630   _nBins_M6 = 100; _min_M6   = 0.;    _max_M6    = 100.;   _width_M6 = (_max_M6-_min_M6)/_nBins_M6;
 631
 632   _min_vertexZ       = _vertexZMin;
 633   _max_vertexZ       = _vertexZMax;
 634   _width_vertexZ     = 0.5;
 635   _nBins_vertexZ     = int(0.5+ (_max_vertexZ - _min_vertexZ)/_width_vertexZ);
 636
 637   _nBins_pt_1        = int(0.5+ (_max_pt_1 -_min_pt_1 )/_width_pt_1);
 638   _nBins_eta_1       = int(0.5+ (_max_eta_1-_min_eta_1)/_width_eta_1);
 639   _width_phi_1       = (_max_phi_1  - _min_phi_1)  /_nBins_phi_1;
 640   _nBins_etaPhi_1    = _nBins_phi_1    * _nBins_eta_1;
 641   _nBins_zEtaPhiPt_1 = _nBins_vertexZ  * _nBins_phi_1 * _nBins_eta_1 * _nBins_pt_1;
 642
 643   _nBins_pt_2        =  int(0.5+ (_max_pt_2 -_min_pt_2 )/_width_pt_2);
 644   _nBins_eta_2       = int(0.5+ (_max_eta_2-_min_eta_2)/_width_eta_2);
 645   _width_phi_2       = (_max_phi_2  - _min_phi_2)  /_nBins_phi_2;
 646   _nBins_etaPhi_2    = _nBins_phi_2    * _nBins_eta_2;
 647   _nBins_zEtaPhiPt_2 = _nBins_vertexZ   * _nBins_phi_2 * _nBins_eta_2 * _nBins_pt_2;
 648
 649   _nBins_pt_3        =  int(0.5+ (_max_pt_3 -_min_pt_3 )/_width_pt_3);
 650   _nBins_eta_3       = int(0.5+ (_max_eta_3-_min_eta_3)/_width_eta_3);
 651   _width_phi_3       = (_max_phi_3  - _min_phi_3)  /_nBins_phi_3;
 652   _nBins_etaPhi_3    = _nBins_phi_3    * _nBins_eta_3;
 653   _nBins_zEtaPhiPt_3 = _nBins_vertexZ   * _nBins_phi_3 * _nBins_eta_3 * _nBins_pt_3;
 654
 655
 656   _nBins_etaPhi_12   = _nBins_etaPhi_1 * _nBins_etaPhi_2;
 657   _nBins_etaPhi_13   = _nBins_etaPhi_1 * _nBins_etaPhi_3;
 658   _nBins_etaPhi_23   = _nBins_etaPhi_2 * _nBins_etaPhi_3;
 659   _nBins_etaPhi_123  = _nBins_etaPhi_1 * _nBins_etaPhi_2 * _nBins_etaPhi_3;
 660
 661   _nBins_phi_12   = _nBins_phi_1 * _nBins_phi_2;
 662   _nBins_phi_13   = _nBins_phi_1 * _nBins_phi_3;
 663   _nBins_phi_23   = _nBins_phi_2 * _nBins_phi_3;
 664   _nBins_phi_123  = _nBins_phi_1 * _nBins_phi_2 * _nBins_phi_3;
 665
 666
 667   //setup the work arrays
 668   if (_singlesOnly)
 669     {
 670     __n1_1_vsPt         = getDoubleArray(_nBins_pt_1,         0.);
 671     __n1_1P_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_1,  0.);
 672     __n1_1M_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_1,  0.);
 673     __n1_2_vsPt         = getDoubleArray(_nBins_pt_2,         0.);
 674     __n1_2P_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_2,  0.);
 675     __n1_2M_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_2,  0.);
 676     __n1_3_vsPt         = getDoubleArray(_nBins_pt_3,         0.);
 677     __n1_3P_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_3,  0.);
 678     __n1_3M_vsZEtaPhiPt = getFloatArray (_nBins_zEtaPhiPt_3,  0.);
 679     }
 680   else
 681     {
 682     _id_1          = new int[arraySize];
 683     //_charge_1      = new int[arraySize];
 684     _iPhi_1        = new int[arraySize];
 685     //_pt_1          = new float[arraySize];
 686     //_px_1          = new float[arraySize];
 687     //_py_1          = new float[arraySize];
 688     //_pz_1          = new float[arraySize];
 689     _correction_1  = new float[arraySize];
 690     //_dedx_1        = new float[arraySize];
 691     __n1_1_vsPhi   = getDoubleArray(_nBins_phi_1,    0.);
 692
 693     _id_2          = new int[arraySize];
 694     //_charge_2      = new int[arraySize];
 695     _iPhi_2        = new int[arraySize];
 696     //_pt_2          = new float[arraySize];
 697     //_px_2          = new float[arraySize];
 698     //_py_2          = new float[arraySize];
 699     //_pz_2          = new float[arraySize];
 700     _correction_2  = new float[arraySize];
 701     //_dedx_2        = new float[arraySize];
 702     __n1_2_vsPhi   = getDoubleArray(_nBins_phi_2,       0.);
 703
 704     _id_3          = new int[arraySize];
 705     //_charge_3      = new int[arraySize];
 706     _iPhi_3        = new int[arraySize];
 707     //_pt_3          = new float[arraySize];
 708     //_px_3          = new float[arraySize];
 709     //_py_3          = new float[arraySize];
 710     //_pz_3          = new float[arraySize];
 711     _correction_3  = new float[arraySize];
 712     //_dedx_3        = new float[arraySize];
 713     __n1_3_vsPhi        = getDoubleArray(_nBins_phi_3,       0.);
 714
 715     __n2_12_vsPhi    = getFloatArray(_nBins_phi_12,    0.);
 716     __n2_13_vsPhi    = getFloatArray(_nBins_phi_13,    0.);
 717     __n2_23_vsPhi    = getFloatArray(_nBins_phi_23,    0.);
 718     __n3_123_vsPhi   = getFloatArray(_nBins_phi_123,   0.);
 719     }
 720
 721
 722   // Setup all the labels needed.
 723   part_1_Name   = "_1";
 724   part_1P_Name  = "_1P";
 725   part_1M_Name  = "_1M";
 726   part_2_Name   = "_2";
 727   part_2P_Name  = "_2P";
 728   part_2M_Name  = "_2M";
 729   part_3_Name   = "_3";
 730   part_3P_Name  = "_3P";
 731   part_3M_Name  = "_3M";
 732   pair_12_Name  = "_12";
 733   pair_13_Name  = "_13";
 734   pair_23_Name  = "_23";
 735   triplet_Name   = "_123";
 736
 737   n1Name     = "n1";
 738   n2Name     = "n2";
 739   n3Name     = "n3";
 740   n1NwName   = "n1Nw";
 741   n2NwName   = "n2Nw";
 742   n3NwName   = "n3Nw";
 743
 744   r1Name     = "r1";
 745   r2Name     = "r2";
 746   r3Name     = "r3";
 747   r2r1Name   = "r2r1";
 748   c2Name     = "c2";
 749   c3Name     = "c3";
 750   d3Name     = "d3";
 751   p3Name     = "p3";
 752
 753   intR2Name       = "intR2";
 754   binCorrName     = "binCorr";
 755   intBinCorrName  = "intBinCorr";
 756
 757
 758   _title_counts = "yield";
 759   _title_m0     = "M_{0}";
 760   _title_m1     = "M_{1}";
 761   _title_m2     = "M_{2}";
 762   _title_m3     = "M_{3}";
 763   _title_m4     = "V0Centrality";
 764   _title_m5     = "TrkCentrality";
 765   _title_m6     = "SpdCentrality";
 766
 767   _title_phi_1       = "#varphi_{1} (radian)";
 768   _title_phi_2       = "#varphi_{2} (radian)";
 769   _title_phi_3       = "#varphi_{3} (radian)";
 770   _title_phi_12   = "#eta_{1}#times#varphi_{1}#times#eta_{2}#times#varphi_{2}";
 771   _title_phi_13   = "#eta_{1}#times#varphi_{1}#times#eta_{3}#times#varphi_{3}";
 772   _title_phi_23   = "#eta_{2}#times#varphi_{2}#times#eta_{3}#times#varphi_{3}";
 773
 774   vsZ         = "_vsZ";
 775   vsM         = "_vsM";
 776   vsPt         = "_vsPt";
 777   vsPhi        = "_vsPhi";
 778   vsEta        = "_vsEta";
 779   vsEtaPhi     = "vsPhi";
 780
 781
 782   if (_useWeights)
 783     {
 784     int iZ, iPhi, iPt;
 785     int iZ1,iPhi1,iPt1;
 786     int a, b;
 787     if (_weight_1P)
 788       {
 789       _correctionWeight_1P = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
 790       _correctionWeight_1M = new float[_nBins_vertexZ*_nBins_etaPhi_1*_nBins_pt_1];
 791       a = _nBins_pt_1;
 792       b = _nBins_etaPhi_1*_nBins_pt_1;
 793       for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 794         {
 795         for (iPhi=0,iPhi1=1; iPhi<_nBins_etaPhi_1; iPhi++, iPhi1++)
 796           {
 797           for (iPt=0,iPt1=1; iPt<_nBins_pt_1; iPt++, iPt1++)
 798             {
 799             _correctionWeight_1P[iZ*b+iPhi*a+iPt] = _weight_1P->GetBinContent(iZ1,iPhi1,iPt1);
 800             _correctionWeight_1M[iZ*b+iPhi*a+iPt] = _weight_1M->GetBinContent(iZ1,iPhi1,iPt1);
 801             }
 802           }
 803         }
 804       } // _weight_1
 805     else
 806       {
 807       AliError("AliAnalysisTask3PCorrelations:: _weight_1 is a null pointer.");
 808       return;
 809       }
 810     if (_weight_2P)
 811       {
 812       _correctionWeight_2P = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
 813       _correctionWeight_2M = new float[_nBins_vertexZ*_nBins_etaPhi_2*_nBins_pt_2];
 814       a = _nBins_pt_2;
 815       b = _nBins_etaPhi_2*_nBins_pt_2;
 816       for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 817         {
 818         for (iPhi=0,iPhi1=1; iPhi<_nBins_etaPhi_2; iPhi++, iPhi1++)
 819           {
 820           for (iPt=0,iPt1=1; iPt<_nBins_pt_2; iPt++, iPt1++)
 821             {
 822             _correctionWeight_2P[iZ*b+iPhi*a+iPt] = _weight_2P->GetBinContent(iZ1,iPhi1,iPt1);
 823             _correctionWeight_2M[iZ*b+iPhi*a+iPt] = _weight_2M->GetBinContent(iZ1,iPhi1,iPt1);
 824             }
 825           }
 826         }
 827       } // _weight_2
 828     else
 829       {
 830       AliError("AliAnalysisTask3PCorrelations:: _weight_2 is a null pointer.");
 831       return;
 832       }
 833     if (_weight_3P)
 834       {
 835       _correctionWeight_3P = new float[_nBins_vertexZ*_nBins_etaPhi_3*_nBins_pt_3];
 836       _correctionWeight_3M = new float[_nBins_vertexZ*_nBins_etaPhi_3*_nBins_pt_3];
 837       a = _nBins_pt_3;
 838       b = _nBins_etaPhi_3*_nBins_pt_3;
 839       for (iZ=0,iZ1=1; iZ<_nBins_vertexZ; iZ++, iZ1++)
 840         {
 841         for (iPhi=0,iPhi1=1; iPhi<_nBins_etaPhi_2; iPhi++, iPhi1++)
 842           {
 843           for (iPt=0,iPt1=1; iPt<_nBins_pt_3; iPt++, iPt1++)
 844             {
 845             _correctionWeight_3P[iZ*b+iPhi*a+iPt] = _weight_3P->GetBinContent(iZ1,iPhi1,iPt1);
 846             _correctionWeight_3M[iZ*b+iPhi*a+iPt] = _weight_3M->GetBinContent(iZ1,iPhi1,iPt1);
 847             }
 848           }
 849         }
 850       } // _weight_2
 851     else
 852       {
 853       AliError("AliAnalysisTask3PCorrelations:: _weight_2 is a null pointer.");
 854       return;
 855       }
 856     }
 857
 858   createHistograms();
 859   PostData(0,_outputHistoList);
 860
 861   cout<< "AliAnalysisTask3PCorrelations::CreateOutputObjects() DONE " << endl;
 862
 863 }
 864
 865 void  AliAnalysisTask3PCorrelations::createHistograms()
 866 {
 867   AliInfo(" AliAnalysisTask3PCorrelations::createHistoHistograms() Creating Event Histos");
 868   TString name;
 869
 870   // bin index : what it is...
 871   //         0 :  number of event submitted
 872   //         1 :  number accepted by centrality cut
 873   //         2 :  number accepted by centrality cut and z cut
 874   //         3 :  total number of particles that satisfied filter 1
 875   //         4 :  total number of particles that satisfied filter 2
 876   //         5 :  total number of particles that satisfied filter 3
 877   name = "eventAccounting";
 878   _eventAccounting      = createHisto1D(name,name,10, -0.5, 9.5, "event Code", _title_counts);
 879
 880   name = "m0"; _m0      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m0, _title_counts);
 881   name = "m1"; _m1      = createHisto1D(name,name,_nBins_M1, _min_M1, _max_M1, _title_m1, _title_counts);
 882   name = "m2"; _m2      = createHisto1D(name,name,_nBins_M2, _min_M2, _max_M2, _title_m2, _title_counts);
 883   name = "m3"; _m3      = createHisto1D(name,name,_nBins_M3, _min_M3, _max_M3, _title_m3, _title_counts);
 884   name = "m4"; _m4      = createHisto1D(name,name,_nBins_M4, _min_M4, _max_M4, _title_m4, _title_counts);
 885   name = "m5"; _m5      = createHisto1D(name,name,_nBins_M5, _min_M5, _max_M5, _title_m5, _title_counts);
 886   name = "m6"; _m6      = createHisto1D(name,name,_nBins_M6, _min_M6, _max_M6, _title_m6, _title_counts);
 887   name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
 888
 889   TString title_vsM;
 890   //_centralityMethod
 891   switch (_centralityMethod)
 892   {
 893     case 0: _nBins_vsM = _nBins_M0; _min_vsM = _min_M0; _max_vsM = _max_M0; _width_vsM = _width_M0; title_vsM = _title_m0; break;
 894     case 1: _nBins_vsM = _nBins_M1; _min_vsM = _min_M1; _max_vsM = _max_M1; _width_vsM = _width_M1; title_vsM = _title_m1; break;
 895     case 2: _nBins_vsM = _nBins_M2; _min_vsM = _min_M2; _max_vsM = _max_M2; _width_vsM = _width_M2; title_vsM = _title_m2; break;
 896     case 3: _nBins_vsM = _nBins_M3; _min_vsM = _min_M3; _max_vsM = _max_M3; _width_vsM = _width_M3; title_vsM = _title_m3; break;
 897     case 4: _nBins_vsM = _nBins_M4; _min_vsM = _min_M4; _max_vsM = _max_M4; _width_vsM = _width_M4; title_vsM = _title_m4; break;
 898     case 5: _nBins_vsM = _nBins_M5; _min_vsM = _min_M5; _max_vsM = _max_M5; _width_vsM = _width_M5; title_vsM = _title_m5; break;
 899     case 6: _nBins_vsM = _nBins_M6; _min_vsM = _min_M6; _max_vsM = _max_M6; _width_vsM = _width_M6; title_vsM = _title_m6; break;
 900   }
 901
 902   AliInfo(" AliAnalysisTask3PCorrelations::createHistoHistograms() Creating Part 1 Histos");
 903
 904   if (_singlesOnly)
 905     {
 906     name = n1Name+part_1_Name+vsPt;               _n1_1_vsPt               = createHisto1F(name,name, _nBins_pt_1,  _min_pt_1,  _max_pt_1,   "pt1", "counts");
 907     name = "dedxVsP_1";                           _dedxVsP_1               = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
 908     name = "betaVsP_1";                           _betaVsP_1               = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
 909     name = n1Name+part_1P_Name+vsZ+vsEtaPhi+vsPt; _n1_1P_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_phi_1,  "pt1");
 910     name = n1Name+part_1M_Name+vsZ+vsEtaPhi+vsPt; _n1_1M_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_phi_1,  "pt1");
 911
 912     name = n1Name+part_2_Name+vsPt;               _n1_2_vsPt               = createHisto1F(name,name, _nBins_pt_2,  _min_pt_2,  _max_pt_2,   "pt2",  "counts");
 913     name = "dedxVsP_2";                           _dedxVsP_2               = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
 914     name = "betaVsP_2";                           _betaVsP_2               = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
 915     name = n1Name+part_2P_Name+vsZ+vsEtaPhi+vsPt; _n1_2P_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_phi_2,  "pt2");
 916     name = n1Name+part_2M_Name+vsZ+vsEtaPhi+vsPt; _n1_2M_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_phi_2,  "pt2");
 917
 918     name = n1Name+part_3_Name+vsPt;               _n1_3_vsPt               = createHisto1F(name,name, _nBins_pt_3,  _min_pt_3,  _max_pt_3,   "pt3",  "counts");
 919     name = "dedxVsP_3";                           _dedxVsP_3               = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
 920     name = "betaVsP_3";                           _betaVsP_3               = createHisto2F(name,name,400,-2.,2.,120,0.5,1.1,"p (GeV/c)", "beta", "counts");
 921     name = n1Name+part_3P_Name+vsZ+vsEtaPhi+vsPt; _n1_3P_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_3, 0., double(_nBins_etaPhi_3), _nBins_pt_3, _min_pt_3, _max_pt_3, "zVertex", _title_phi_3,  "pt3");
 922     name = n1Name+part_3M_Name+vsZ+vsEtaPhi+vsPt; _n1_3M_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_3, 0., double(_nBins_etaPhi_3), _nBins_pt_3, _min_pt_3, _max_pt_3, "zVertex", _title_phi_3,  "pt3");
 923     }
 924   else
 925     {
 926     name = n1Name+part_1_Name+vsPhi;  _n1_1_vsPhi   = createHisto1F(name,name, _nBins_phi_1, _min_phi_1, _max_phi_1, _title_phi_1,  "counts");
 927     name = n1Name+part_2_Name+vsPhi;  _n1_2_vsPhi   = createHisto1F(name,name, _nBins_phi_2, _min_phi_2, _max_phi_2, _title_phi_2,  "counts");
 928     name = n1Name+part_3_Name+vsPhi;  _n1_3_vsPhi   = createHisto1F(name,name, _nBins_phi_3, _min_phi_3, _max_phi_3, _title_phi_3,  "counts");
 929
 930     name = n1Name+part_1_Name + vsM;  _n1_1_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{1}>");
 931     name = n1Name+part_2_Name + vsM;  _n1_2_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{2}>");
 932     name = n1Name+part_3_Name + vsM;  _n1_3_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{3}>");
 933
 934     name = n1NwName+part_1_Name+vsM;  _n1Nw_1_vsM   = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{1}>");
 935     name = n1NwName+part_2_Name+vsM;  _n1Nw_2_vsM   = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{2}>");
 936     name = n1NwName+part_3_Name+vsM;  _n1Nw_3_vsM   = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{3}>");
 937
 938     name = n2Name+pair_12_Name+vsPhi; _n2_12_vsPhi  = createHisto2F(name,name, _nBins_phi_1, _min_phi_1, _max_phi_1, _nBins_phi_2, _min_phi_2, _max_phi_2, _title_phi_1, _title_phi_2, "<N_{12}>");
 939     name = n2Name+pair_13_Name+vsPhi; _n2_13_vsPhi  = createHisto2F(name,name, _nBins_phi_1, _min_phi_1, _max_phi_1, _nBins_phi_3, _min_phi_3, _max_phi_3, _title_phi_2, _title_phi_3, "<N_{13}>");
 940     name = n2Name+pair_23_Name+vsPhi; _n2_23_vsPhi  = createHisto2F(name,name, _nBins_phi_2, _min_phi_2, _max_phi_2, _nBins_phi_3, _min_phi_3, _max_phi_3, _title_phi_1, _title_phi_3, "<N_{23}>");
 941     name = n3Name+triplet_Name+vsPhi; _n3_123_vsPhi = createHisto3F(name,name, _nBins_phi_1, _min_phi_1, _max_phi_1, _nBins_phi_2, _min_phi_2, _max_phi_2, _nBins_phi_3, _min_phi_3, _max_phi_3, _title_phi_1, _title_phi_2, _title_phi_3);
 942
 943     name = n2Name+pair_12_Name+vsM;  _n2_12_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{12}>");
 944     name = n2Name+pair_13_Name+vsM;  _n2_13_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{13}>");
 945     name = n2Name+pair_23_Name+vsM;  _n2_23_vsM     = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{23}>");
 946     name = n3Name+pair_23_Name+vsM;  _n3_123_vsM    = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{123}>");
 947     name = n2NwName+pair_12_Name+vsM; _n2Nw_12_vsM  = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{12}>");
 948     name = n2NwName+pair_13_Name+vsM; _n2Nw_13_vsM  = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{13}>");
 949     name = n2NwName+pair_23_Name+vsM; _n2Nw_23_vsM  = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{23}>");
 950     name = n3NwName+triplet_Name+vsM; _n3Nw_123_vsM = createProfile(name,name, _nBins_vsM,   _min_vsM,   _max_vsM,  title_vsM,      "<N_{123}>");
 951
 952     }
 953
 954   AliInfo(" AliAnalysisTask3PCorrelations::createHistoHistograms() All Done");
 955 }
 956 //-----------------------//
 957
 958 void  AliAnalysisTask3PCorrelations::finalizeHistograms()
 959 {
 960
 961   AliInfo("AliAnalysisTask3PCorrelations::finalizeHistograms() starting");
 962   AliInfo(Form("CorrelationAnalyzers::finalizeHistograms()   _eventCount : %d",int(_eventCount)));
 963
 964   if (_singlesOnly)
 965     {
 966     fillHistoWithArray(_n1_1_vsPt,               __n1_1_vsPt,         _nBins_pt_1);
 967     fillHistoWithArray(_n1_1P_vsZVsEtaVsPhiVsPt, __n1_1P_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
 968     fillHistoWithArray(_n1_1M_vsZVsEtaVsPhiVsPt, __n1_1M_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_1, _nBins_pt_1);
 969
 970     fillHistoWithArray(_n1_2_vsPt,               __n1_2_vsPt,         _nBins_pt_2);
 971     fillHistoWithArray(_n1_2P_vsZVsEtaVsPhiVsPt, __n1_2P_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_2, _nBins_pt_2);
 972     fillHistoWithArray(_n1_2M_vsZVsEtaVsPhiVsPt, __n1_2M_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_2, _nBins_pt_2);
 973
 974     fillHistoWithArray(_n1_3_vsPt,               __n1_3_vsPt,         _nBins_pt_3);
 975     fillHistoWithArray(_n1_3P_vsZVsEtaVsPhiVsPt, __n1_3P_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_3, _nBins_pt_3);
 976     fillHistoWithArray(_n1_3M_vsZVsEtaVsPhiVsPt, __n1_3M_vsZEtaPhiPt, _nBins_vertexZ, _nBins_etaPhi_3, _nBins_pt_3);
 977
 978     }
 979   else
 980     {
 981     fillHistoWithArray(_n1_1_vsPhi,    __n1_1_vsPhi,   _nBins_phi_1);
 982     fillHistoWithArray(_n1_2_vsPhi,    __n1_2_vsPhi,   _nBins_phi_2);
 983     fillHistoWithArray(_n1_3_vsPhi,    __n1_3_vsPhi,   _nBins_phi_3);
 984     fillHistoWithArray(_n2_12_vsPhi,   __n2_12_vsPhi,  _nBins_phi_1, _nBins_phi_2);
 985     fillHistoWithArray(_n2_13_vsPhi,   __n2_13_vsPhi,  _nBins_phi_1, _nBins_phi_3);
 986     fillHistoWithArray(_n2_23_vsPhi,   __n2_23_vsPhi,  _nBins_phi_2, _nBins_phi_3);
 987     fillHistoWithArray(_n3_123_vsPhi,  __n3_123_vsPhi, _nBins_phi_1, _nBins_phi_2, _nBins_phi_3);
 988     }
 989   AliInfo("AliAnalysisTask3PCorrelations::finalizeHistograms()  Done ");
 990 }
 991 //--------------//
 992
 993
 994 void  AliAnalysisTask3PCorrelations::UserExec(Option_t */*option*/)
 995 {
 996   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - Starting!!!!" << endl;
 997
 998   int    k1=0,k2=0,k3=0;
 999   int    iPhi=0, iEta=0, iEtaPhi=0, iPt=0, charge=0;
1000   float  q=0, p=0, phi=0, pt=0, eta=0, corr=0, dedx=0, px=0, py=0, pz=0;
1001   int    ij_12=0, ij_13=0, ij_23=0, ijk_123=0;
1002   int    id_1=0, id_2=0, id_3=0, iPhi_1=0, iPhi_2=0, iPhi_3=0;
1003   float  corr_1=0, corr_2=0, corr_3=0;
1004   float  corr_12=0, corr_13=0, corr_23=0, corr_123=0;
1005   int    iVertex=0, iVertexP1=0, iVertexP2=0, iVertexP3=0;
1006   int    iZEtaPhiPt=0;
1007   //float  massElecSq = 2.5e-7=0;
1008   double b[2];
1009   double bCov[3];
1010   const  AliAODVertex*  vertex=0;
1011   int    nClus=0;
1012   bool   bitOK=0;
1013
1014   // count all events looked at here
1015   _eventCount++;
1016
1017   if (_eventAccounting)
1018     {
1019     _eventAccounting->Fill(0);// count all calls to this function
1020     }
1021   else
1022     {
1023     cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - !_eventAccounting" << endl;
1024     cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - Skip this task" << endl;
1025     return;
1026     }
1027
1028   //Get pointer to current event
1029   //fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
1030   fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
1031
1032   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 2" << endl;
1033
1034   if(!fAODEvent)
1035     {
1036     AliError("AliAnalysisTask3PCorrelations::Exec(Option_t *option) !fAODEvent");
1037     return;
1038     }
1039
1040
1041   float v0Centr  = -999.;
1042   float trkCentr = -999.;
1043   float spdCentr = -999.;
1044   float vertexX  = -999;
1045   float vertexY  = -999;
1046   float vertexZ  = -999;
1047   float vertexXY = -999;
1048   float dcaZ     = -999;
1049   float dcaXY    = -999;
1050   float centrality = -999;
1051   k1 = k2 = k3 = 0;
1052   __n1_1 = __n1_2 = __n1_3 = __n1Nw_1 = __n1Nw_2 = __n1Nw_3 = 0;
1053
1054   _eventAccounting->Fill(1);// count all calls to this function with a valid pointer
1055
1056   //Centrality
1057   AliCentrality* centralityObject =  fAODEvent->GetHeader()->GetCentralityP();
1058   if (centralityObject)
1059     {
1060     // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 6" << endl;
1061
1062     v0Centr  = centralityObject->GetCentralityPercentile("V0M");
1063     trkCentr = centralityObject->GetCentralityPercentile("TRK");
1064     spdCentr = centralityObject->GetCentralityPercentile("CL1");
1065     // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 7" << endl;
1066
1067     }
1068   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 8" << endl;
1069
1070   _nTracks  = fAODEvent->GetNTracks();
1071   _mult0    = 0;
1072   _mult1    = 0;
1073   _mult2    = 0;
1074   _mult3    = _nTracks;
1075   _mult4    = v0Centr;
1076   _mult5    = trkCentr;
1077   _mult6    = spdCentr;
1078   _field    = fAODEvent->GetMagneticField();
1079
1080   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - _field:" << _field << endl;
1081
1082   //_centralityMethod
1083   switch (_centralityMethod)
1084   {
1085     case 0: centrality = _mult0; break;
1086     case 1: centrality = _mult1; break;
1087     case 2: centrality = _mult2; break;
1088     case 3: centrality = _mult3; break;
1089     case 4: centrality = _mult4; break;
1090     case 5: centrality = _mult5; break;
1091     case 6: centrality = _mult6; break;
1092   }
1093
1094   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 10" << endl;
1095
1096   //filter on centrality
1097   if ( centrality < _centralityMin ||
1098       centrality > _centralityMax)
1099     {
1100     // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 11" << endl;
1101
1102     // we dont analyze this event here...
1103     return;
1104     }
1105   _eventAccounting->Fill(2);// count all events with right centrality
1106                             // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 12" << endl;
1107
1108   // filter on z and xy vertex
1109   vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertexSPD();
1110   if (!vertex || vertex->GetNContributors()<1)
1111     {
1112     vertexZ  = -999;
1113     vertexXY = -999;
1114     // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - No valid vertex object or poor vertex" << endl;
1115     }
1116   else
1117     {
1118       vertexX = vertex->GetX();
1119       vertexY = vertex->GetY();
1120       vertexZ = vertex->GetZ();
1121       vertexXY = sqrt(vertexX*vertexX+vertexY*vertexY);
1122     }
1123   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 13" << endl;
1124
1125   // cout << "vertexZ : " << vertexZ  << endl;
1126   // cout << "vertexXY: " << vertexXY << endl;
1127   // cout << "_vertexZMin: " << _vertexZMin << endl;
1128   // cout << "_vertexZMax: " << _vertexZMax << endl;
1129   // cout << "_vertexXYMin: " << _vertexXYMin << endl;
1130   // cout << "_vertexXYMax: " << _vertexXYMax << endl;
1131
1132   if (!vertex ||
1133       vertexZ  < _vertexZMin  ||
1134       vertexZ  > _vertexZMax  ||
1135       vertexXY < _vertexXYMin ||
1136       vertexXY > _vertexXYMax)
1137     return;
1138
1139
1140   // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 14" << endl;
1141
1142   iVertex = int((vertexZ-_min_vertexZ)/_width_vertexZ);
1143   iVertexP1 = iVertex*_nBins_etaPhiPt_1;
1144   iVertexP2 = iVertex*_nBins_etaPhiPt_2;
1145   iVertexP3 = iVertex*_nBins_etaPhiPt_3;
1146   if (iVertex<0 || iVertex>=_nBins_vertexZ)
1147     {
1148     AliError("AliAnalysisTask3PCorrelations::Exec(Option_t *option) iVertex<0 || iVertex>=_nBins_vertexZ ");
1149     return;
1150     }
1151   _eventAccounting->Fill(3);// count all calls to this function with a valid pointer
1152                             // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 15" << endl;
1153
1154
1155   // loop over all particles for singles
1156   //debug() << "_nTracks:"<< _nTracks << endl;
1157   //reset single particle counters
1158   for (int iTrack=0; iTrack< _nTracks; iTrack++)
1159     {
1160     // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - 16 iTrack:" << iTrack << endl;
1161
1162     AliAODTrack * t = (AliAODTrack *) fAODEvent->GetTrack(iTrack);
1163     if (!t)
1164       {
1165       AliError(Form("AliAnalysisTask3PCorrelations::Exec(Option_t *option) No track ofr iTrack=%d", iTrack));
1166       continue;
1167       }
1168
1169     q      = t->Charge();
1170     charge = int(q);
1171     phi    = t->Phi();
1172     p      = t->P();
1173     pt     = t->Pt();
1174     px     = t->Px();
1175     py     = t->Py();
1176     pz     = t->Pz();
1177     eta    = t->Eta();
1178     dedx   = t->GetTPCsignal();
1179     nClus  = t->GetTPCNcls();
1180     bitOK  = t->TestFilterBit(_trackFilterBit);
1181     if (!bitOK ||
1182         dedx  < _dedxMin   ||
1183         dedx  > _dedxMax   ||
1184         nClus < _nClusterMin) continue;
1185
1186     // cout << "_trackFilterBit:" << _trackFilterBit << " Track returns:" << bitOK << endl;
1187     // cout << "   pt:" << pt   << " _min_pt_1:" << _min_pt_1 << " _max_pt_1:" << _max_pt_1<< endl;
1188     // cout << "  phi:" << phi  << endl;
1189     // cout << "  eta:" << eta  << " _min_eta_1:" << _min_eta_1 << " _max_eta_1:" << _max_eta_1<< endl;
1190     // cout << " dedx:" << dedx << " _dedxMin:" << _dedxMin << " _dedxMax:" << _dedxMax << endl;
1191     // cout << "nclus:" << nClus<< " _nClusterMin:" << _nClusterMin << endl;
1192     if (pt       >=  _min_pt_1  &&
1193         pt       <   _max_pt_1  &&
1194         eta      >=  _min_eta_1 &&
1195         eta      <   _max_eta_1)
1196       {
1197       // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - check vertex  for 1:" << endl;
1198       // Get the dca information
1199       if (t->PropagateToDCA(vertex, _field, 100., b, bCov) )
1200         {
1201         dcaXY = b[0];
1202         dcaZ  = b[1];
1203         }
1204       else
1205         {
1206         dcaXY = -999999;
1207         dcaZ  = -999999;
1208         }
1209       // cout << "1 dcaZ:" << dcaZ << " _dcaZMin:" << _dcaZMin << " _dcaZMax:" << _dcaZMax << endl;
1210       // cout << "1 dcaXY:" << dcaXY << " _dcaXYMin:" << _dcaXYMin << " _dcaXYMax:" << _dcaXYMax << endl;
1211
1212       // skip track if DCA too large
1213       if (dcaZ     >=  _dcaZMin &&
1214           dcaZ     <   _dcaZMax &&
1215           dcaXY    >=  _dcaXYMin &&
1216           dcaXY    <   _dcaXYMax)
1217         continue; //track does not have a valid DCA
1218                   // cout << "keep track:" << endl;
1219
1220       iPhi   = int( phi/_width_phi_1);
1221       // cout << " AliAnalysisTask3PCorrelations::analyze(Event * event) -1- iTrack:" << iTrack<< endl<< "pt:" << pt << " phi:" <<  phi << " eta:" << eta << endl;
1222       if (iPhi<0 || iPhi>=_nBins_phi_1 )
1223         {
1224         AliWarning("AliAnalysisTask3PCorrelations::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1225         return;
1226         }
1227
1228       iEta    = int((eta-_min_eta_1)/_width_eta_1);
1229       if (iEta<0 || iEta>=_nBins_eta_1)
1230         {
1231         AliWarning(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1232         continue;
1233         }
1234       iPt     = int((pt -_min_pt_1 )/_width_pt_1 );
1235       if (iPt<0  || iPt >=_nBins_pt_1)
1236         {
1237         AliWarning(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1238         continue;
1239         }
1240       iEtaPhi    = iEta*_nBins_phi_1+iPhi;
1241       iZEtaPhiPt = iVertexP1 + iEtaPhi*_nBins_pt_1 + iPt;
1242
1243       if (charge>0 && _correctionWeight_1P)
1244         corr = _correctionWeight_1P[iZEtaPhiPt];
1245       else if (charge<0 && _correctionWeight_1M)
1246         corr = _correctionWeight_1M[iZEtaPhiPt];
1247       else
1248         corr = 1;
1249       // cout << "all good; process track:" << endl;
1250       if (_singlesOnly)
1251         {
1252         __n1_1_vsPt[iPt]                += corr;
1253         if (charge>0)
1254           __n1_1P_vsZEtaPhiPt[iZEtaPhiPt] += corr;
1255         else if (charge>0)
1256           __n1_1M_vsZEtaPhiPt[iZEtaPhiPt] += corr;
1257         }
1258       else
1259         {
1260         _id_1[k1]            = iTrack;
1261         //_charge_1[k1]        = charge;
1262         _iPhi_1[k1]          = iPhi;
1263         _correction_1[k1]    = corr;
1264         __n1_1               += corr;
1265         __n1_1_vsPhi[iPhi]   += corr;
1266         __n1Nw_1             += 1;
1267         ++k1;
1268         if (k1>=arraySize)
1269           {
1270           AliError(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) k1 >=arraySize; arraySize: %d",arraySize));
1271           return;
1272           }
1273         }
1274       }
1275
1276     else if (pt  >= _min_pt_2  &&
1277              pt  <  _max_pt_2  &&
1278              eta >= _min_eta_2 &&
1279              eta <  _max_eta_2)
1280       {
1281       // Get the dca information
1282       // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - check vertex for 2:" << endl;
1283       if (t->PropagateToDCA(vertex, _field, 100., b, bCov) )
1284         {
1285         dcaXY = b[0];
1286         dcaZ  = b[1];
1287         }
1288       else
1289         {
1290         dcaXY = -999999;
1291         dcaZ  = -999999;
1292         }
1293       // cout << "2 dcaZ:" << dcaZ << " _dcaZMin:" << _dcaZMin << " _dcaZMax:" << _dcaZMax << endl;
1294       // cout << "2 dcaXY:" << dcaXY << " _dcaXYMin:" << _dcaXYMin << " _dcaXYMax:" << _dcaXYMax << endl;
1295       // skip track if DCA too large
1296       if (dcaZ     >=  _dcaZMin &&
1297           dcaZ     <   _dcaZMax &&
1298           dcaXY    >=  _dcaXYMin &&
1299           dcaXY    <   _dcaXYMax)
1300         continue; //track does not have a valid DCA
1301
1302       iPhi   = int( phi/_width_phi_2);
1303       //cout << " AliAnalysisTask3PCorrelations::analyze(Event * event) -1- iTrack:" << iTrack  << endl
1304       //<< "pt:" << pt << " phi:" <<  phi << " eta:" << eta << endl;
1305       if (iPhi<0 || iPhi>=_nBins_phi_2 )
1306         {
1307         AliWarning("AliAnalysisTask3PCorrelations::analyze() iPhi<0 || iPhi>=_nBins_phi_1");
1308         return;
1309         }
1310
1311       iEta    = int((eta-_min_eta_2)/_width_eta_2);
1312       if (iEta<0 || iEta>=_nBins_eta_2)
1313         {
1314         AliWarning(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1315         continue;
1316         }
1317       iPt     = int((pt -_min_pt_2 )/_width_pt_2 );
1318       if (iPt<0  || iPt >=_nBins_pt_2)
1319         {
1320         AliWarning(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) Mismatched iPt: %d",iPt));
1321         continue;
1322         }
1323
1324       iEtaPhi    = iEta*_nBins_phi_2+iPhi;
1325       iZEtaPhiPt = iVertexP2 + iEtaPhi*_nBins_pt_2 + iPt;
1326       if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2)
1327         {
1328         AliWarning("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_2");
1329         continue;
1330         }
1331
1332       if (charge>0 && _correctionWeight_2P)
1333         corr = _correctionWeight_2P[iZEtaPhiPt];
1334       else if (charge<0 && _correctionWeight_2M)
1335         corr = _correctionWeight_2M[iZEtaPhiPt];
1336       else
1337         corr = 1;
1338
1339       if (_singlesOnly)
1340         {
1341         //_dedxVsP_2->Fill(p*q,dedx);
1342         __n1_2_vsPt[iPt]               += corr;          // cout << "step 15" << endl;
1343         if (charge>0)
1344           __n1_2P_vsZEtaPhiPt[iZEtaPhiPt] += corr;
1345         else if (charge>0)
1346           __n1_2M_vsZEtaPhiPt[iZEtaPhiPt] += corr;          }
1347       else
1348         {
1349         _id_2[k2]           = iTrack;         // cout << "step 1" << endl;
1350                                               //_charge_2[k2]       = charge;         // cout << "step 2" << endl;
1351         _iPhi_2[k2]         = iPhi;        // cout << "step 3" << endl;
1352         _correction_2[k2]   = corr;           // cout << "step 9" << endl;
1353         __n1_2              += corr;          // cout << "step 10" << endl;
1354         __n1Nw_2            += 1;
1355         __n1_2_vsPhi[iPhi]  += corr;          // cout << "step 11" << endl;
1356         ++k2;
1357         if (k2>=arraySize)
1358           {
1359           AliWarning(Form("-W- k2 >=arraySize; arraySize: %d",arraySize));
1360           return;
1361           }
1362         }
1363       }
1364
1365     else if (pt  >= _min_pt_3  &&
1366              pt  <  _max_pt_3  &&
1367              eta >= _min_eta_3 &&
1368              eta <  _max_eta_3)
1369       {
1370       // Get the dca information
1371       // cout << "AliAnalysisTask3PCorrelations::UserExec(Option_t *option) - check vertex for 2:" << endl;
1372       if (t->PropagateToDCA(vertex, _field, 100., b, bCov) )
1373         {
1374         dcaXY = b[0];
1375         dcaZ  = b[1];
1376         }
1377       else
1378         {
1379         dcaXY = -999999;
1380         dcaZ  = -999999;
1381         }
1382       // skip track if DCA too large
1383       if (dcaZ     >=  _dcaZMin &&
1384           dcaZ     <   _dcaZMax &&
1385           dcaXY    >=  _dcaXYMin &&
1386           dcaXY    <   _dcaXYMax)
1387         continue; //track does not have a valid DCA
1388
1389       iPhi   = int( phi/_width_phi_3);
1390       if (iPhi<0 || iPhi>=_nBins_phi_3 )
1391         {
1392         AliWarning("AliAnalysisTask3PCorrelations::analyze() iPhi<0 || iPhi>=_nBins_phi_3");
1393         return;
1394         }
1395
1396       iEta    = int((eta-_min_eta_3)/_width_eta_3);
1397       if (iEta<0 || iEta>=_nBins_eta_3)
1398         {
1399         AliWarning(Form("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) Mismatched iEta: %d", iEta));
1400         continue;
1401         }
1402       iEtaPhi    = iEta*_nBins_phi_3+iPhi;
1403       iZEtaPhiPt = iVertexP3 + iEtaPhi*_nBins_pt_3 + iPt;
1404       if (iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_3)
1405         {
1406         AliWarning("AliAnalysisTask3PCorrelations::analyze(AliceEvent * event) iZEtaPhiPt<0 || iZEtaPhiPt>=_nBins_zEtaPhiPt_3");
1407         continue;
1408         }
1409       if (charge>0 && _correctionWeight_3P)
1410         corr = _correctionWeight_3P[iZEtaPhiPt];
1411       else if (charge<0 && _correctionWeight_3M)
1412         corr = _correctionWeight_3M[iZEtaPhiPt];
1413       else
1414         corr = 1;
1415
1416       if (_singlesOnly)
1417         {
1418         __n1_3_vsPt[iPt]               += corr;          // cout << "step 15" << endl;
1419         if (charge>0)
1420           __n1_3P_vsZEtaPhiPt[iZEtaPhiPt] += corr;
1421         else if (charge>0)
1422           __n1_3M_vsZEtaPhiPt[iZEtaPhiPt] += corr;          }
1423       else
1424         {
1425         _id_3[k3]            = iTrack;         // cout << "step 1" << endl;
1426                                                //_charge_3[k3]        = charge;         // cout << "step 2" << endl;
1427         _iPhi_3[k3]          = iPhi;        // cout << "step 3" << endl;
1428         _correction_3[k3]    = corr;           // cout << "step 9" << endl;
1429         __n1_3               += corr;          // cout << "step 10" << endl;
1430         __n1Nw_3             += 1;
1431         __n1_3_vsPhi[iPhi]   += corr;          // cout << "step 11" << endl;
1432         ++k3;
1433         if (k3>=arraySize)
1434           {
1435           AliWarning(Form("-W- k3 >=arraySize; arraySize: %d",arraySize));
1436           return;
1437           }
1438         }
1439       }
1440
1441     // cout << "done with track" << endl;
1442     } //iTrack
1443
1444   // cout << "Filling histograms now" << endl;
1445   _m0->Fill(_mult0);
1446   _m1->Fill(_mult1);
1447   _m2->Fill(_mult2);
1448   _m3->Fill(_mult3);
1449   _m4->Fill(_mult4);
1450   _m5->Fill(_mult5);
1451   _m6->Fill(_mult6);
1452   _vertexZ->Fill(vertexZ);
1453
1454
1455   //if singels only selected, do not fill pair histograms.
1456   if (_singlesOnly)
1457     {
1458
1459     }
1460   else
1461     {
1462     // reset pair counters
1463     _n1_1_vsM->Fill(    centrality, __n1_1);
1464     _n1Nw_1_vsM->Fill(  centrality, __n1Nw_1);
1465     _n1_2_vsM->Fill(    centrality, __n1_2);
1466     _n1Nw_2_vsM->Fill(  centrality, __n1Nw_2);
1467     _n1_3_vsM->Fill(    centrality, __n1_3);
1468     _n1Nw_3_vsM->Fill(  centrality, __n1Nw_3);
1469
1470     __n2_12 = __n2Nw_12 = __n2_13 = __n2Nw_13 = __n2_23 = __n2Nw_23 = __n3_123 = __n3Nw_123 = 0;
1471
1472     for (int i1=0; i1<k1; i1++)
1473       {
1474       // cout << "         i1:" << i1 << endl;
1475       id_1    = _id_1[i1];
1476       iPhi_1  = _iPhi_1[i1];
1477       corr_1  = _correction_1[i1];
1478       //1 and 2
1479       for (int i2=0;i2<k2;i2++)
1480         {
1481           id_2 = _id_2[i2];
1482           if (id_1!=id_2)
1483             {
1484             iPhi_2  = _iPhi_2[i2];
1485             corr_2  = _correction_2[i2];
1486             corr_12 = corr_1*corr_2;
1487             ij_12 = iPhi_1*_nBins_phi_1 + iPhi_2;
1488             __n2_12_vsPhi[ij_12] += corr_12;
1489             __n2_12              += corr_12;
1490             __n2Nw_12            += 1;
1491             for (int i3=0;i3<k3;i3++)
1492               {
1493                 id_3 = _id_3[i3];
1494                 if (id_1!=id_3 && id_2!=id_3)
1495                   {
1496                   iPhi_3   = _iPhi_3[i3];
1497                   corr_3   = _correction_3[i3];
1498                   corr_123 = corr_12*corr_3;
1499                   ijk_123  = ij_12 *_nBins_phi_2 + iPhi_3;
1500                   __n3_123_vsPhi[ijk_123] += corr_123;
1501                   __n3_123                += corr_123;
1502                   __n3Nw_123              += 1;
1503                   }
1504               } //i3
1505             }
1506         } //i2
1507           //1 and 3
1508       for (int i3=0;i3<k3;i3++)
1509         {
1510           id_3 = _id_3[i3];
1511           if (id_1!=id_3)
1512             {
1513             iPhi_3  = _iPhi_3[i3];
1514             corr_3  = _correction_3[i3];
1515             corr_13 = corr_1*corr_3;
1516             ij_13   = iPhi_1*_nBins_phi_1 + iPhi_3;
1517             __n2_13_vsPhi[ij_13] += corr_13;
1518             __n2_13              += corr_13;
1519             __n2Nw_13            += 1;
1520             }
1521         } //i3
1522       } //i1
1523
1524     //2 & 3
1525     for (int i2=0; i2<k2; i2++)
1526       {
1527       id_2    = _id_2[i2];
1528       iPhi_2  = _iPhi_2[i2];
1529       corr_2  = _correction_2[i2];
1530       for (int i3=0;i3<k3;i3++)
1531         {
1532           id_3 = _id_3[i3];
1533           if (id_2!=id_3)
1534             {
1535             iPhi_3  = _iPhi_3[i3];
1536             corr_3  = _correction_3[i3];
1537             corr_23 = corr_2*corr_3;
1538             ij_23   = iPhi_2*_nBins_phi_2 + iPhi_3;
1539             __n2_23_vsPhi[ij_23] += corr_23;
1540             __n2_23              += corr_23;
1541             __n2Nw_23            += 1;
1542             }
1543         } //i2
1544       }
1545     _n2_12_vsM->Fill(   centrality, __n2_12);
1546     _n2_13_vsM->Fill(   centrality, __n2_13);
1547     _n2_23_vsM->Fill(   centrality, __n2_23);
1548     _n3_123_vsM->Fill(  centrality, __n3_123);
1549     _n2Nw_12_vsM->Fill( centrality, __n2Nw_12);
1550     _n2Nw_13_vsM->Fill( centrality, __n2Nw_13);
1551     _n2Nw_23_vsM->Fill( centrality, __n2Nw_23);
1552     _n3Nw_123_vsM->Fill(centrality, __n3_123);
1553     }
1554
1555   // cout << "Event Done " << endl;
1556   PostData(0,_outputHistoList);
1557
1558 }
1559
1560 void   AliAnalysisTask3PCorrelations::FinishTaskOutput()
1561 {
1562   cout << "AliAnalysisTask3PCorrelations::FinishTaskOutput() Starting." << endl;
1563   finalizeHistograms();
1564   PostData(0,_outputHistoList);
1565   cout << "AliAnalysisTask3PCorrelations::FinishTaskOutput() Done." << endl;
1566 }
1567
1568 void   AliAnalysisTask3PCorrelations::Terminate(Option_t* /*option*/)
1569 {
1570   // no ops
1571 }
1572
1573
1574 //Tools
1575 //===================================================================================================
1576 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH1 * h, float * array, int size)
1577 {
1578   int i, i1;
1579   float v1, ev1, v2, ev2, sum, esum;
1580   for (i=0, i1=1; i<size; ++i,++i1)
1581     {
1582     v1  = array[i]; ev1 = sqrt(v1);
1583     v2  = h->GetBinContent(i1);
1584     ev2 = h->GetBinError(i1);
1585     sum = v1 + v2;
1586     esum = sqrt(ev1*ev1+ev2*ev2);
1587     h->SetBinContent(i1,sum);
1588     h->SetBinError(i1,esum);
1589     }
1590 }
1591
1592 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH2 * h, float * array, int size1, int size2)
1593 {
1594   int i, i1;
1595   int j, j1;
1596   float v1, ev1, v2, ev2, sum, esum;
1597   for (i=0, i1=1; i<size1; ++i,++i1)
1598     {
1599     for (j=0, j1=1; j<size2; ++j,++j1)
1600       {
1601       v1  = array[i*size2+j]; ev1 = sqrt(v1);
1602       v2  = h->GetBinContent(i1,j1);
1603       ev2 = h->GetBinError(i1,j1);
1604       sum = v1 + v2;
1605       esum = sqrt(ev1*ev1+ev2*ev2);
1606       h->SetBinContent(i1,j1,sum);
1607       h->SetBinError(i1,j1,esum);
1608       }
1609     }
1610 }
1611
1612 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH3 * h, float * array, int size1, int size2, int size3)
1613 {
1614   int i, i1;
1615   int j, j1;
1616   int k, k1;
1617   float v1, ev1, v2, ev2, sum, esum;
1618   int size23 = size2*size3;
1619   for (i=0, i1=1; i<size1; ++i,++i1)
1620     {
1621     for (j=0, j1=1; j<size2; ++j,++j1)
1622       {
1623       for (k=0, k1=1; k<size3; ++k,++k1)
1624         {
1625         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
1626         v2  = h->GetBinContent(i1,j1,k1);
1627         ev2 = h->GetBinError(i1,j1,k1);
1628         sum = v1 + v2;
1629         esum = sqrt(ev1*ev1+ev2*ev2);
1630         h->SetBinContent(i1,j1,k1,sum);
1631         h->SetBinError(i1,j1,k1,esum);
1632         }
1633       }
1634     }
1635 }
1636
1637 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH1 * h, double * array, int size)
1638 {
1639   int i, i1;
1640   double v1, ev1, v2, ev2, sum, esum;
1641   for (i=0, i1=1; i<size; ++i,++i1)
1642     {
1643     v1  = array[i]; ev1 = sqrt(v1);
1644     v2  = h->GetBinContent(i1);
1645     ev2 = h->GetBinError(i1);
1646     sum = v1 + v2;
1647     esum = sqrt(ev1*ev1+ev2*ev2);
1648     h->SetBinContent(i1,sum);
1649     h->SetBinError(i1,esum);
1650     }
1651 }
1652
1653 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH2 * h, double * array, int size1, int size2)
1654 {
1655   int i, i1;
1656   int j, j1;
1657   double v1, ev1, v2, ev2, sum, esum;
1658   for (i=0, i1=1; i<size1; ++i,++i1)
1659     {
1660     for (j=0, j1=1; j<size2; ++j,++j1)
1661       {
1662       v1  = array[i*size2+j]; ev1 = sqrt(v1);
1663       v2  = h->GetBinContent(i1,j1);
1664       ev2 = h->GetBinError(i1,j1);
1665       sum = v1 + v2;
1666       esum = sqrt(ev1*ev1+ev2*ev2);
1667       h->SetBinContent(i1,j1,sum);
1668       h->SetBinError(i1,j1,esum);
1669       }
1670     }
1671 }
1672
1673 void  AliAnalysisTask3PCorrelations::fillHistoWithArray(TH3 * h, double * array, int size1, int size2, int size3)
1674 {
1675   int i, i1;
1676   int j, j1;
1677   int k, k1;
1678   double v1, ev1, v2, ev2, sum, esum;
1679   int size23 = size2*size3;
1680   for (i=0, i1=1; i<size1; ++i,++i1)
1681     {
1682     for (j=0, j1=1; j<size2; ++j,++j1)
1683       {
1684       for (k=0, k1=1; k<size3; ++k,++k1)
1685         {
1686         v1  = array[i*size23+j*size3+k]; ev1 = sqrt(v1);
1687         v2  = h->GetBinContent(i1,j1,k1);
1688         ev2 = h->GetBinError(i1,j1,k1);
1689         sum = v1 + v2;
1690         esum = sqrt(ev1*ev1+ev2*ev2);
1691         h->SetBinContent(i1,j1,k1,sum);
1692         h->SetBinError(i1,j1,k1,esum);
1693         }
1694       }
1695     }
1696 }
1697
1698 //________________________________________________________________________
1699 double *  AliAnalysisTask3PCorrelations::getDoubleArray(int size, double v)
1700 {
1701   /// Allocate an array of type double with n values
1702   /// Initialize the array to the given value
1703   double * array = new double [size];
1704   for (int i=0;i<size;++i) array[i]=v;
1705   return array;
1706 }
1707
1708 //________________________________________________________________________
1709 float *  AliAnalysisTask3PCorrelations::getFloatArray(int size, float v)
1710 {
1711   /// Allocate an array of type float with n values
1712   /// Initialize the array to the given value
1713   float * array = new float [size];
1714   for (int i=0;i<size;++i) array[i]=v;
1715   return array;
1716 }
1717
1718
1719 //________________________________________________________________________
1720 TH1D * AliAnalysisTask3PCorrelations::createHisto1D(const TString &  name, const TString &  title,
1721                                                     int n, double xMin, double xMax,
1722                                                     const TString &  xTitle, const TString &  yTitle)
1723 {
1724   //CreateHisto new 1D historgram
1725   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
1726   TH1D * h = new TH1D(name,title,n,xMin,xMax);
1727   h->GetXaxis()->SetTitle(xTitle);
1728   h->GetYaxis()->SetTitle(yTitle);
1729   addToList(h);
1730   return h;
1731 }
1732
1733
1734 //________________________________________________________________________
1735 TH1D * AliAnalysisTask3PCorrelations::createHisto1D(const TString &  name, const TString &  title,
1736                                                     int n, double * bins,
1737                                                     const TString &  xTitle, const TString &  yTitle)
1738 {
1739   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
1740   TH1D * h = new TH1D(name,title,n,bins);
1741   h->GetXaxis()->SetTitle(xTitle);
1742   h->GetYaxis()->SetTitle(yTitle);
1743   addToList(h);
1744   return h;
1745 }
1746
1747
1748 //________________________________________________________________________
1749 TH2D * AliAnalysisTask3PCorrelations::createHisto2D(const TString &  name, const TString &  title,
1750                                                     int nx, double xMin, double xMax, int ny, double yMin, double yMax,
1751                                                     const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
1752 {
1753   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));
1754   TH2D * h = new TH2D(name,title,nx,xMin,xMax,ny,yMin,yMax);
1755   h->GetXaxis()->SetTitle(xTitle);
1756   h->GetYaxis()->SetTitle(yTitle);
1757   h->GetZaxis()->SetTitle(zTitle);
1758   addToList(h);
1759   return h;
1760 }
1761
1762 //________________________________________________________________________
1763 TH2D * AliAnalysisTask3PCorrelations::createHisto2D(const TString &  name, const TString &  title,
1764                                                     int nx, double* xbins, int ny, double yMin, double yMax,
1765                                                     const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
1766 {
1767   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
1768   TH2D * h;
1769   h = new TH2D(name,title,nx,xbins,ny,yMin,yMax);
1770   h->GetXaxis()->SetTitle(xTitle);
1771   h->GetYaxis()->SetTitle(yTitle);
1772   h->GetZaxis()->SetTitle(zTitle);
1773   addToList(h);
1774   return h;
1775 }
1776
1777 //// F /////
1778 //________________________________________________________________________
1779 TH1F * AliAnalysisTask3PCorrelations::createHisto1F(const TString &  name, const TString &  title,
1780                                                     int n, double xMin, double xMax,
1781                                                     const TString &  xTitle, const TString &  yTitle)
1782 {
1783   //CreateHisto new 1D historgram
1784   AliInfo(Form("createHisto 1D histo %s   nBins: %d  xMin: %f   xMax: %f",name.Data(),n,xMin,xMax));
1785   TH1F * h = new TH1F(name,title,n,xMin,xMax);
1786   h->GetXaxis()->SetTitle(xTitle);
1787   h->GetYaxis()->SetTitle(yTitle);
1788   addToList(h);
1789   return h;
1790 }
1791
1792
1793 //________________________________________________________________________
1794 TH1F * AliAnalysisTask3PCorrelations::createHisto1F(const TString &  name, const TString &  title,
1795                                                     int n, double * bins,
1796                                                     const TString &  xTitle, const TString &  yTitle)
1797 {
1798   AliInfo(Form("createHisto 1D histo %s   with %d non uniform nBins",name.Data(),n));
1799   TH1F * h = new TH1F(name,title,n,bins);
1800   h->GetXaxis()->SetTitle(xTitle);
1801   h->GetYaxis()->SetTitle(yTitle);
1802   addToList(h);
1803   return h;
1804 }
1805
1806
1807 //________________________________________________________________________
1808 TH2F * AliAnalysisTask3PCorrelations::createHisto2F(const TString &  name, const TString &  title,
1809                                                     int nx, double xMin, double xMax, int ny, double yMin, double yMax,
1810                                                     const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
1811 {
1812   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));
1813   TH2F * h = new TH2F(name,title,nx,xMin,xMax,ny,yMin,yMax);
1814   h->GetXaxis()->SetTitle(xTitle);
1815   h->GetYaxis()->SetTitle(yTitle);
1816   h->GetZaxis()->SetTitle(zTitle);
1817   addToList(h);
1818   return h;
1819 }
1820
1821 //________________________________________________________________________
1822 TH2F * AliAnalysisTask3PCorrelations::createHisto2F(const TString &  name, const TString &  title,
1823                                                     int nx, double* xbins, int ny, double yMin, double yMax,
1824                                                     const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
1825 {
1826   AliInfo(Form("createHisto 2D histo %s   with %d non uniform nBins",name.Data(),nx));
1827   TH2F * h;
1828   h = new TH2F(name,title,nx,xbins,ny,yMin,yMax);
1829   h->GetXaxis()->SetTitle(xTitle);
1830   h->GetYaxis()->SetTitle(yTitle);
1831   h->GetZaxis()->SetTitle(zTitle);
1832   addToList(h);
1833   return h;
1834 }
1835
1836
1837 //________________________________________________________________________
1838 TH3F * AliAnalysisTask3PCorrelations::createHisto3F(const TString &  name, const TString &  title,
1839                                                     int nx, double xMin, double xMax,
1840                                                     int ny, double yMin, double yMax,
1841                                                     int nz, double zMin, double zMax,
1842                                                     const TString &  xTitle, const TString &  yTitle, const TString &  zTitle)
1843 {
1844   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));
1845   TH3F * h = new TH3F(name,title,nx,xMin,xMax,ny,yMin,yMax,nz,zMin,zMax);
1846   h->GetXaxis()->SetTitle(xTitle);
1847   h->GetYaxis()->SetTitle(yTitle);
1848   h->GetZaxis()->SetTitle(zTitle);
1849   addToList(h);
1850   return h;
1851 }
1852
1853 //________________________________________________________________________
1854 TProfile * AliAnalysisTask3PCorrelations::createProfile(const TString & name, const TString & description,
1855                                                         int nx,double xMin,double xMax,
1856                                                         const TString &  xTitle, const TString &  yTitle)
1857 {
1858   AliInfo(Form("createHisto 1D profile %s   nBins: %d  xMin: %f10.4 xMax: %f10.4",name.Data(),nx,xMin,xMax));
1859   TProfile * h = new TProfile(name,description,nx,xMin,xMax);
1860   h->GetXaxis()->SetTitle(xTitle);
1861   h->GetYaxis()->SetTitle(yTitle);
1862   addToList(h);
1863   return h;
1864 }
1865
1866 //________________________________________________________________________
1867 TProfile * AliAnalysisTask3PCorrelations::createProfile(const TString &  name,const TString &  description,
1868                                                         int nx,  double* bins,
1869                                                         const TString &  xTitle, const TString &  yTitle)
1870 {
1871   AliInfo(Form("createHisto 1D profile %s  with %d non uniform bins",name.Data(),nx));
1872   TProfile * h = new TProfile(name,description,nx,bins);
1873   h->GetXaxis()->SetTitle(xTitle);
1874   h->GetYaxis()->SetTitle(yTitle);
1875   addToList(h);
1876   return h;
1877 }
1878
1879
1880 void   AliAnalysisTask3PCorrelations::addToList(TH1 *h)
1881 {
1882   if (_outputHistoList)
1883     {
1884     _outputHistoList->Add(h);
1885     }
1886   else
1887     cout << "addToList(TH1 *h) _outputHistoList is null!!!!! Shoudl abort ship" << endl;
1888
1889 }
1890
1891
1892