PWGCF/Correlations/DPhi/AliAnalysisTaskpxpx.cxx

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