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