prod/acrcaf/qa_pp/AliAnalysisTaskQASym.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TCanvas.h"
   6 #include "TList.h"
   7 #include "TParticle.h"
   8 #include "TParticlePDG.h"
   9 #include "TProfile.h"
  10 #include "TNtuple.h"
  11 #include "TFile.h"
  12
  13 #include "AliAnalysisTask.h"
  14 #include "AliAnalysisManager.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliLog.h"
  18 #include "AliESDVertex.h"
  19 #include "AliESDInputHandler.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliMultiplicity.h"
  22
  23
  24 #include "AliAnalysisTaskQASym.h"
  25 #include "AliExternalTrackParam.h"
  26 #include "AliTrackReference.h"
  27
  28 #include "AliHeader.h"
  29 #include "AliGenEventHeader.h"
  30 #include "AliGenDPMjetEventHeader.h"
  31
  32 // Analysis Task for basic QA on the ESD
  33
  34 // Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing,
  35 //          Andreas Morsch, Eva Sicking
  36
  37 ClassImp(AliAnalysisTaskQASym)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskQASym::AliAnalysisTaskQASym(const char *name)
  41   : AliAnalysisTaskSE(name)
  42     ,fFieldOn(kTRUE)
  43     ,fHists(0)
  44     ,fHistRECpt(0)
  45     ,fEta(0)
  46     ,fEtaPhi(0)
  47     ,fEtaPt(0)
  48     ,fQPt(0)
  49     ,fDca(0)
  50     ,fqRec(0)
  51     ,fsigmaPt(0)
  52
  53     ,fRecPtPos(0)
  54     ,fRecPtNeg(0)
  55     ,fRecPhiPos(0)
  56     ,fRecPhiNeg(0)
  57     ,fRecEtaPos(0)
  58     ,fRecEtaNeg(0)
  59     ,fRecEtaPtPos(0)
  60     ,fRecEtaPtNeg(0)
  61     ,fRecDcaPos(0)
  62     ,fRecDcaNeg(0)
  63     ,fRecDcaNegInv(0)
  64     ,fRecDPos(0)
  65     ,fRecDNeg(0)
  66
  67     ,fRecQPtPosEta(0)
  68     ,fRecQPtNegEta(0)
  69     ,fRecPtPosEta(0)
  70     ,fRecPtNegEta(0)
  71     ,fRecPhiPosEta(0)
  72     ,fRecPhiNegEta(0)
  73     ,fRecDcaPosEta(0)
  74     ,fRecDcaNegEta(0)
  75     ,fRecDPosEta(0)
  76     ,fRecDNegEta(0)
  77
  78     ,fRecPtPosVz(0)
  79     ,fRecPtNegVz(0)
  80     ,fRecEtaPosVz(0)
  81     ,fRecEtaNegVz(0)
  82     ,fRecPhiPosVz(0)
  83     ,fRecPhiNegVz(0)
  84     ,fSignedDcaPosVz(0)
  85     ,fSignedDcaNegVz(0)
  86     ,fRecQPtPosEtaVz(0)
  87     ,fRecQPtNegEtaVz(0)
  88     ,fRecEtaPtPosVz(0)
  89     ,fRecEtaPtNegVz(0)
  90
  91
  92     ,fDeltaPhiAll(0)
  93     ,fDeltaPhiLeading(0)
  94     ,fDiffDcaD(0)
  95     ,fPhiRec(0)
  96     ,fThetaRec(0)
  97     ,fNumber(0)
  98     ,fVx(0)
  99     ,fVy(0)
 100     ,fVz(0)
 101     ,test(0)
 102
 103
 104     ,fCuts(0)
 105
 106
 107
 108 {
 109     // Constructor
 110   for(Int_t i = 0;i<18;++i){
 111     fRecPtTpcSector[i] = 0;
 112     fRecEtaTpcSector[i] = 0;
 113     fSignedDcaTpcSector[i] = 0;
 114     fRecQPtTpcSector[i] = 0;
 115     fRecEtaPtTpcSector[i] = 0;
 116   }
 117
 118   for(Int_t i = 0;i< 7;++i){
 119     fRecPtPosLadder[i] = 0;
 120     fRecPtNegLadder[i] = 0;
 121     fRecPhiPosLadder[i] = 0;
 122     fRecPhiNegLadder[i] = 0;
 123     fRecEtaPosLadder[i] = 0;
 124     fRecEtaNegLadder[i] = 0;
 125     fSignDcaPos[i] = 0;
 126     fSignDcaNeg[i] = 0;
 127     fSignDcaNegInv[i] = 0;
 128     fPtSigmaPos[i] =0;
 129     fPtSigmaNeg[i] =0;
 130     fqPtRec[i] =0;
 131     fDcaSigmaPos[i] =0;
 132     fDcaSigmaNeg[i] =0;
 133   }
 134
 135   DefineOutput(1,  TList::Class());
 136
 137
 138
 139 }
 140
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskQASym::UserCreateOutputObjects()
 144 {
 145   // Create histograms
 146   // Called once
 147
 148   Double_t range = 300.;
 149   Double_t pt = 20.;
 150
 151   fHists = new TList();
 152   //  test   = new TNtuple("test","test",
 153   //                      "pt:phi:theta:x:y:z:charge");
 154   fHistRECpt   = new TH1F("fHistRECpt",
 155                           " p_{T}",
 156                           100, 0., pt);
 157   fEta   = new TH1F("fEta",
 158                           " #eta",
 159                           200, -2., 2.);
 160   fEtaPhi   = new TH2F("fEtaPhi",
 161                           " #eta - #phi",
 162                           200, -2., 2., 128, 0., 2. * TMath::Pi());
 163
 164   fThetaRec   = new TH1F("fThetaRec",
 165                           " #theta",
 166                           180, 0., TMath::Pi());
 167   fPhiRec   = new TH1F("fPhiRec",
 168                           " #phi",
 169                           180, 0., 2*TMath::Pi());
 170   fNumber   = new TH1F("fNumber",
 171                        "number of tracks per event",
 172                        50, 0.5, 49.5);
 173   fVx   = new TH1F("fVx",
 174                    "X of vertex",
 175                    500, -500., 500.);
 176   fVy   = new TH1F("fVy",
 177                    "Y of vertex",
 178                    500, -500., 500.);
 179   fVz   = new TH1F("fVz",
 180                    "Z of vertex",
 181                    500, -500., 500.);
 182
 183   fEtaPt   = new TH1F("fEtaPt",
 184                           " #eta/p_{T} ",
 185                           100, -1., 1.);
 186
 187   fQPt   = new TH1F("fQPt",
 188                           " charge/p_{T} ",
 189                           100, -1., 1.);
 190
 191   fDca   = new TH1F("fDca",
 192                           " dca ",
 193                           200, -range, range);
 194
 195
 196   fqRec    = new TH1F("fqRec",
 197                           " charge all reconstructed particle",
 198                           21, -9.5, 10.5);
 199
 200   fsigmaPt    = new TH1F("fsigmaPt",
 201                           "Log_{10}(#sigma_{p_{T}})",
 202                           200, -2., 8.);
 203
 204
 205
 206
 207   //------------
 208   for(Int_t ITSlayer_case=0;ITSlayer_case<7;ITSlayer_case++){
 209
 210     fSignDcaPos[ITSlayer_case]   = new TH1F(Form("fSignDcaPos%d", ITSlayer_case),
 211                                    " Signed dca",
 212                                    200, -range, range);
 213     fSignDcaPos[ITSlayer_case]->GetXaxis()->SetTitle("dca");
 214     fSignDcaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
 215
 216
 217     fSignDcaNeg[ITSlayer_case]   = new TH1F(Form("fSignDcaNeg%d", ITSlayer_case),
 218                                    " Signed dcas",
 219                                    200, -range, range);
 220     fSignDcaNeg[ITSlayer_case]->GetXaxis()->SetTitle("dca");
 221     fSignDcaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
 222
 223     fSignDcaNegInv[ITSlayer_case]   = new TH1F(Form("fSignDcaNegInv%d", ITSlayer_case),
 224                                    " inverse Signed dca ",
 225                                    200, -range, range);
 226     fSignDcaNegInv[ITSlayer_case]->GetXaxis()->SetTitle("-dca");
 227     fSignDcaNegInv[ITSlayer_case]->GetYaxis()->SetTitle("");
 228
 229
 230
 231
 232     fPtSigmaPos[ITSlayer_case]   = new TH1F(Form("fPtSigmaPos%d", ITSlayer_case),
 233                                    " #sigma_{pT} ",
 234                                    208, -2., 8.);
 235     fPtSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
 236     fPtSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
 237
 238
 239     fPtSigmaNeg[ITSlayer_case]   = new TH1F(Form("fPtSigmaNeg%d",ITSlayer_case),
 240                                   " #sigma_{pT}",
 241                                    208, -2., 8.);
 242     fPtSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
 243     fPtSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
 244
 245
 246
 247
 248
 249     fqPtRec[ITSlayer_case]   = new TH1F(Form("fqPtRec%d",ITSlayer_case),
 250                                   "q/ p_{T}",
 251                                    200, -100., 100.);
 252     fqPtRec[ITSlayer_case]->GetXaxis()->SetTitle("q_{tr}/p_{T, tr} (GeV/c)");
 253     fqPtRec[ITSlayer_case]->GetYaxis()->SetTitle("");
 254
 255
 256
 257
 258
 259     fDcaSigmaPos[ITSlayer_case]   = new TH2F(Form("fDcaSigmaPos%d", ITSlayer_case),
 260                                    " p_{T} shift vs #sigma_{pT} ",
 261                                    200, -range, range,200, -4., 4. );
 262     fDcaSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA)");
 263     fDcaSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
 264
 265
 266     fDcaSigmaNeg[ITSlayer_case]   = new TH2F(Form("fDcaSigmaNeg%d", ITSlayer_case),
 267                                    " p_{T} shift vs #sigma_{pT} ",
 268                                    200, -range, range,200, -4., 4. );
 269     fDcaSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA");
 270     fDcaSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
 271
 272
 273
 274
 275
 276
 277     // YIELDs---------- positive and negative particles
 278
 279     fRecPtPos   = new TH1F("fRecPtPos",
 280                            " p_{T}",
 281                            100, 0.,pt);
 282     fRecPtPos->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 283     fRecPtNeg   = new TH1F("fRecPtNeg",
 284                            " p_{T} ",
 285                            100, 0., pt);
 286     fRecPtNeg->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 287
 288
 289     fRecPhiPos   = new TH1F("fRecPhiPos",
 290                             "#phi",
 291                             361, 0., 360.);
 292     fRecPhiPos->GetXaxis()->SetTitle("#phi (deg)");
 293
 294     fRecPhiNeg   = new TH1F("fRecPhiNeg",
 295                             "#phi ",
 296                             361, 0., 360.);
 297     fRecPhiNeg->GetXaxis()->SetTitle("#phi (deg)");
 298
 299     fRecEtaPos   = new TH1F("fRecEtaPos",
 300                             "#eta",
 301                             200, -2., 2.);
 302     fRecEtaPos->GetXaxis()->SetTitle("#eta");
 303
 304     fRecEtaNeg   = new TH1F("fRecEtaNeg",
 305                             "#eta",
 306                             200, -2., 2.);
 307     fRecEtaNeg->GetXaxis()->SetTitle("#eta");
 308
 309     fRecEtaPtPos   = new TH1F("fRecEtaPtPos",
 310                               "#eta/p_{T}",
 311                               200, -0.1, .1);
 312     fRecEtaPtPos->GetXaxis()->SetTitle("#eta/p_{T}");
 313
 314     fRecEtaPtNeg   = new TH1F("fRecEtaPtNeg",
 315                               "#eta/p_{T}",
 316                               200, -.1, .1);
 317     fRecEtaPtNeg->GetXaxis()->SetTitle("#eta/p_{T}");
 318
 319     fRecDcaPos   = new TH1F("fRecDcaPos",
 320                          " dca",
 321                            100, -range, range);
 322     fRecDcaPos->GetXaxis()->SetTitle("dca (cm)");
 323     fRecDcaNeg   = new TH1F("fRecDcaNeg",
 324                            " dca",
 325                            100, -range, range);
 326     fRecDcaNeg->GetXaxis()->SetTitle("dca (cm)");
 327
 328     fRecDcaNegInv   = new TH1F("fRecDcaNegInv",
 329                            " dca",
 330                            100, -range, range);
 331     fRecDcaNegInv->GetXaxis()->SetTitle("dca (cm)");
 332
 333
 334     fRecDPos   = new TH1F("fRecDPos",
 335                          " d",
 336                            100, -range, range);
 337     fRecDPos->GetXaxis()->SetTitle("d (cm)");
 338     fRecDNeg   = new TH1F("fRecDNeg",
 339                            "d",
 340                            100, -range, range);
 341     fRecDNeg->GetXaxis()->SetTitle("d (cm)");
 342
 343
 344     //  YIELDs ---------------- positive and negative eta
 345
 346
 347     fRecQPtPosEta   = new TH1F("fRecQPtPosEta",
 348                                "q/p_{T}",
 349                                200, -0.5, 0.5);
 350     fRecQPtPosEta->GetXaxis()->SetTitle("q/p_{T} ");
 351
 352     fRecQPtNegEta   = new TH1F("fRecQPtNegEta",
 353                                "q/p_{T}",
 354                                200, -0.5, 0.5);
 355     fRecQPtNegEta->GetXaxis()->SetTitle("q/p_{T}");
 356
 357     fRecPtPosEta   = new TH1F("fRecPtPosEta",
 358                               " p_{T} ",
 359                               100, 0., pt);
 360     fRecPtPosEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 361
 362     fRecPtNegEta   = new TH1F("fRecPtNegEta",
 363                               " p_{T}",
 364                               100, 0., pt);
 365     fRecPtNegEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 366
 367     fRecPhiPosEta   = new TH1F("fRecPhiPosEta",
 368                             "#phi",
 369                             361, 0., 360);
 370     fRecPhiPosEta->GetXaxis()->SetTitle("#phi (deg)");
 371
 372     fRecPhiNegEta   = new TH1F("fRecPhiNegEta",
 373                             "#phi ",
 374                             361, 0, 360);
 375     fRecPhiNegEta->GetXaxis()->SetTitle("#phi (deg)");
 376
 377     fRecDcaPosEta   = new TH1F("fRecDcaPosEta",
 378                          " dca ",
 379                            100, -range, range);
 380     fRecDcaPosEta->GetXaxis()->SetTitle("dca (cm)");
 381     fRecDcaNegEta   = new TH1F("fRecDcaNegEta",
 382                            " dca",
 383                            100, -range, range);
 384     fRecDcaNegEta->GetXaxis()->SetTitle("dca (cm)");
 385
 386     fRecDPosEta   = new TH1F("fRecDPosEta",
 387                          " d",
 388                            100, -range, range);
 389     fRecDPosEta->GetXaxis()->SetTitle("d (cm)");
 390     fRecDNegEta   = new TH1F("fRecDNegEta",
 391                            "d",
 392                            100, -5., 5.);
 393     fRecDNegEta->GetXaxis()->SetTitle("d (cm)");
 394
 395
 396
 397     //  YIELDs ---------------- for TPC sectors
 398     for(Int_t sector=0; sector<18;sector++){
 399
 400
 401       fRecPtTpcSector[sector]   = new TH1F(Form("fRecPtTpcSector%02d",sector),
 402                                            Form("p_{T} distribution: TPC sector %d",
 403                                                 sector),100, 0., pt);
 404       fRecPtTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 405
 406       fRecEtaTpcSector[sector]   = new TH1F(Form("fRecEtaTpcSector%02d",sector),
 407                                            Form("#eta distribution: TPC sector %d",
 408                                                 sector),200, -2., 2.);
 409       fRecEtaTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 410
 411
 412       fSignedDcaTpcSector[sector]   = new TH1F(Form("fSignedDcaTpcSector%02d",sector),
 413                                            Form("dca distribution: TPC sector %d",
 414                                                 sector),200, -range, range );
 415       fSignedDcaTpcSector[sector]->GetXaxis()->SetTitle("dca");
 416
 417       fRecQPtTpcSector[sector]   = new TH1F(Form("fRecQPtTpcSector%02d",sector),
 418                                            Form("Q/ p_{T} distribution: TPC sector %d",
 419                                                 sector),100, -1., 1.);
 420       fRecQPtTpcSector[sector]->GetXaxis()->SetTitle("Q/p_{T} (GeV/c)");
 421
 422       fRecEtaPtTpcSector[sector]   = new TH1F(Form("fRecEtaPtTpcSector%02d",sector),
 423                                            Form("#eta/ p_{T} distribution: TPC sector %d",
 424                                                 sector),100, -1., 1.);
 425       fRecEtaPtTpcSector[sector]->GetXaxis()->SetTitle("#eta/p_{T} (GeV/c)");
 426
 427     }
 428     // YIELDS ITS ladder
 429     for(Int_t i=0;i<7;i++){
 430       fRecPtPosLadder[i]   = new TH1F(Form("fRecPtPosLadder%d", i),
 431                              " p_{T} distribution",
 432                              100, 0., pt);
 433       fRecPtPosLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 434       fRecPtNegLadder[i]   = new TH1F(Form("fRecPtNegLadder%d",i),
 435                              " p_{T} distribution ",
 436                              100, 0., pt);
 437       fRecPtNegLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 438
 439
 440       fRecPhiPosLadder[i]   = new TH1F(Form("fRecPhiPosLadder%d",i),
 441                                  "#phi distribution: all pos eta",
 442                                  361, 0., 360);
 443       fRecPhiPosLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
 444
 445       fRecPhiNegLadder[i]   = new TH1F(Form("fRecPhiNegLadder%d", i),
 446                                  "#phi distribution: all neg eta",
 447                                  361, 0, 360);
 448       fRecPhiNegLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
 449
 450
 451
 452       fRecEtaPosLadder[i]   = new TH1F(Form("fRecEtaPosLadder%d",i),
 453                                        "#eta distribution",
 454                                  200, -2., 2.);
 455       fRecEtaPosLadder[i]->GetXaxis()->SetTitle("#eta)");
 456
 457       fRecEtaNegLadder[i]   = new TH1F(Form("fRecEtaNegLadder%d", i),
 458                                  "#eta distribution",
 459                                  200, -2., 2.);
 460       fRecEtaNegLadder[i]->GetXaxis()->SetTitle("#eta");
 461     }
 462
 463     Double_t vzmax = 30;
 464
 465     fRecPtPosVz = new TH2F("fRecPtPosVz",
 466                            "p_{T} distribution vs Vz()",
 467                            100, -1., 2., 200,-vzmax,vzmax);
 468     fRecPtPosVz->GetXaxis()->SetTitle("log_{10}(p_{T})");
 469
 470     fRecPtNegVz = new TH2F("fRecPtNegVz",
 471                            "p_{T} distribution vs Vz()",
 472                            100, -1., 2.,200,-vzmax,vzmax);
 473     fRecPtNegVz->GetXaxis()->SetTitle("Log_{10}(p_{T})");
 474
 475
 476     fRecEtaPosVz= new TH2F("fRecEtaPosVz",
 477                           "#eta distribution vs Vz()",
 478                           100, -2., 2., 200,-vzmax,vzmax);
 479     fRecEtaPosVz->GetXaxis()->SetTitle("#eta");
 480     fRecEtaNegVz = new TH2F("fRecEtaNegVz",
 481                            "#eta distribution vs Vz()",
 482                            100, -2., 2.,200,-vzmax,vzmax);
 483     fRecEtaNegVz->GetXaxis()->SetTitle("#eta");
 484
 485     fRecPhiPosVz= new TH2F("fRecPhiPosVz",
 486                           "#eta distribution vs Vz()",
 487                           361, 0., 360., 200,-vzmax,vzmax);
 488     fRecPhiPosVz->GetXaxis()->SetTitle("#phi (deg)");
 489     fRecPhiNegVz = new TH2F("fRecPhiNegVz",
 490                            "dca vs Vz()",
 491                            361, 0., 360.,200,-vzmax,vzmax);
 492     fRecPhiNegVz->GetXaxis()->SetTitle("#phi (deg)");
 493
 494     fSignedDcaPosVz= new TH2F("fSignedDcaPosVz",
 495                           "#eta distribution vs Vz()",
 496                           200, -range, range, 200,-vzmax,vzmax);
 497     fSignedDcaPosVz->GetXaxis()->SetTitle("dca (cm)");
 498     fSignedDcaNegVz = new TH2F("fSignedDcaNegVz",
 499                            "dca vs Vz()",
 500                            200, -range, range,200,-vzmax,vzmax);
 501     fSignedDcaNegVz->GetXaxis()->SetTitle("dca (cm)");
 502
 503     fRecQPtPosEtaVz= new TH2F("fRecQPtPosEtaVz",
 504                           " Q/p_{T} distribution vs Vz()",
 505                           100, -1., 1., 200,-vzmax,vzmax);
 506     fRecQPtPosEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
 507     fRecQPtNegEtaVz = new TH2F("fRecQPtNegEtaVz",
 508                            " Q/p_{T} distribution vs Vz()",
 509                            100, -1., 1.,200,-vzmax,vzmax);
 510     fRecQPtNegEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
 511
 512
 513     fRecEtaPtPosVz= new TH2F("fRecEtaPtPosVz",
 514                           " #eta/p_{T} distribution vs Vz()",
 515                           100, -1., 1., 200,-vzmax,vzmax);
 516     fRecEtaPtPosVz->GetXaxis()->SetTitle("#eta/p_{T");
 517     fRecEtaPtNegVz = new TH2F("fRecEtaPtNegVz",
 518                            " #eta/p_{T} distribution vs Vz()",
 519                            100, -1., 1.,200,-vzmax,vzmax);
 520     fRecEtaPtNegVz->GetXaxis()->SetTitle("#eta/p_{T}");
 521
 522     //new
 523     fDeltaPhiAll = new TH1F("fDeltaPhiAll",
 524                            " #Delta #phi",200,-360,360);
 525     fDeltaPhiAll->GetXaxis()->SetTitle("#Delta #phi");
 526
 527
 528     fDeltaPhiLeading = new TH2F("fDeltaPhiLeading",
 529                            " #Delta #phi",361,-360,360, 361,0, 360);
 530     fDeltaPhiLeading->GetXaxis()->SetTitle("#Delta #phi (deg.)");
 531     fDeltaPhiLeading->GetYaxis()->SetTitle("#phi_{leading particle} (deg.)");
 532
 533     fDiffDcaD    = new TH1F("fDiffDcaD",
 534                             "dca-d",
 535                             200, -5., 5.);
 536
 537   }
 538
 539   fHists->SetOwner();
 540
 541   fHists->Add(fHistRECpt);
 542   fHists->Add(fEta);
 543   fHists->Add(fEtaPhi);
 544   fHists->Add(fThetaRec);
 545   fHists->Add(fPhiRec);
 546   fHists->Add(fNumber);
 547   fHists->Add(fVx);
 548   fHists->Add(fVy);
 549   fHists->Add(fVz);
 550
 551   fHists->Add(fEtaPt);
 552   fHists->Add(fQPt);
 553   fHists->Add(fDca);
 554
 555   fHists->Add(fDeltaPhiAll);
 556   fHists->Add(fDeltaPhiLeading);
 557   fHists->Add(fDiffDcaD);
 558
 559   fHists->Add(fqRec);
 560   fHists->Add(fsigmaPt);
 561
 562   fHists->Add(fRecPtPos);
 563   fHists->Add(fRecPtNeg);
 564   fHists->Add(fRecPhiPos);
 565   fHists->Add(fRecPhiNeg);
 566   fHists->Add(fRecEtaPos);
 567   fHists->Add(fRecEtaNeg);
 568   fHists->Add(fRecEtaPtPos);
 569   fHists->Add(fRecEtaPtNeg);
 570   fHists->Add(fRecDcaPos);
 571   fHists->Add(fRecDcaNeg);
 572   fHists->Add(fRecDcaNegInv);
 573   fHists->Add(fRecDPos);
 574   fHists->Add(fRecDNeg);
 575
 576
 577   fHists->Add(fRecQPtPosEta);
 578   fHists->Add(fRecQPtNegEta);
 579   fHists->Add(fRecPtPosEta);
 580   fHists->Add(fRecPtNegEta);
 581   fHists->Add(fRecPhiPosEta);
 582   fHists->Add(fRecPhiNegEta);
 583   fHists->Add(fRecDcaPosEta);
 584   fHists->Add(fRecDcaNegEta);
 585   fHists->Add(fRecDPosEta);
 586   fHists->Add(fRecDNegEta);
 587
 588
 589   for(Int_t i=0;i<18;i++){
 590     fHists->Add(fRecPtTpcSector[i]);
 591     fHists->Add(fRecEtaTpcSector[i]);
 592     fHists->Add(fSignedDcaTpcSector[i]);
 593     fHists->Add(fRecQPtTpcSector[i]);
 594     fHists->Add(fRecEtaPtTpcSector[i]);
 595   }
 596
 597   for(Int_t i=0;i<7;i++){
 598     fHists->Add(fRecPtPosLadder[i]);
 599     fHists->Add(fRecPtNegLadder[i]);
 600     fHists->Add(fRecPhiPosLadder[i]);
 601     fHists->Add(fRecPhiNegLadder[i]);
 602     fHists->Add(fRecEtaPosLadder[i]);
 603     fHists->Add(fRecEtaNegLadder[i]);
 604   }
 605
 606   fHists->Add(fRecPtPosVz);
 607   fHists->Add(fRecPtNegVz);
 608   fHists->Add(fRecEtaPosVz);
 609   fHists->Add(fRecEtaNegVz);
 610   fHists->Add(fRecPhiPosVz);
 611   fHists->Add(fRecPhiNegVz);
 612   fHists->Add(fSignedDcaPosVz);
 613   fHists->Add(fSignedDcaNegVz);
 614   fHists->Add(fRecQPtPosEtaVz);
 615   fHists->Add(fRecQPtNegEtaVz);
 616   fHists->Add(fRecEtaPtPosVz);
 617   fHists->Add(fRecEtaPtNegVz);
 618
 619
 620   for(Int_t i=0;i<7;i++){
 621     fHists->Add(fSignDcaPos[i]);
 622     fHists->Add(fSignDcaNeg[i]);
 623     fHists->Add(fSignDcaNegInv[i]);
 624
 625     fHists->Add(fPtSigmaPos[i]);
 626     fHists->Add(fPtSigmaNeg[i]);
 627     fHists->Add(fqPtRec[i]);
 628
 629     fHists->Add(fDcaSigmaPos[i]);
 630     fHists->Add(fDcaSigmaNeg[i]);
 631
 632
 633   }
 634
 635
 636
 637   for (Int_t i=0; i<fHists->GetEntries(); ++i) {
 638     TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));
 639     if (h1){
 640       // Printf("%s ",h1->GetName());
 641       h1->Sumw2();
 642     }
 643   }
 644     // BKC
 645
 646
 647 }
 648
 649 //__________________________________________________________
 650
 651 void AliAnalysisTaskQASym::UserExec(Option_t *)
 652 {
 653   printf("I'm here \n");
 654   AliVEvent *event = InputEvent();
 655   if (!event) {
 656      Printf("ERROR: Could not retrieve event");
 657      return;
 658   }
 659
 660
 661   if(Entry()==0){
 662     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 663     if(esd){
 664       Printf("We are reading from ESD");
 665     }
 666
 667   }
 668
 669   Printf("There are %d tracks in this event", event->GetNumberOfTracks());
 670
 671
 672   Int_t   LeadingTrack  =   0;
 673   Float_t LeadingEnergy = -20.;
 674   Float_t LeadingPhi    =   0;//TMath::Pi();
 675
 676
 677   if(event->GetNumberOfTracks()!=0) fNumber->Fill(event->GetNumberOfTracks());
 678
 679   const AliVVertex* vertex = event->GetPrimaryVertex();
 680   Float_t vz = vertex->GetZ();
 681   if (TMath::Abs(vz) > 10.) return;
 682
 683   for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
 684
 685
 686     AliVParticle *track = event->GetTrack(iTrack);
 687     AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
 688     if (!track) {
 689       Printf("ERROR: Could not receive track %d", iTrack);
 690       continue;
 691     }
 692
 693     if (!fCuts->AcceptTrack(esdtrack)) continue;
 694     //const AliExternalTrackParam * tpcP = esdtrack->GetTPCInnerParam();
 695     const AliExternalTrackParam *  tpcP = esdtrack;
 696     if (!tpcP) continue;
 697
 698
 699     if (tpcP->Pt() > 50.) {
 700       printf("%5d %5d %13.3f %13.3f \n", event->GetPeriodNumber(), event->GetOrbitNumber(),
 701              tpcP->Pt(), esdtrack->Pt());
 702       //      AliFatal("Jet");
 703     }
 704
 705     if(tpcP->E()>LeadingEnergy){
 706       LeadingTrack=iTrack;
 707       LeadingEnergy=tpcP->E();
 708       LeadingPhi=tpcP->Phi();
 709     }
 710
 711
 712
 713
 714     //propagate to dca
 715     esdtrack->PropagateToDCA(event->GetPrimaryVertex(),
 716                              event->GetMagneticField(), 10000.);
 717
 718
 719     // if(tpcP->Pt()<2.)continue;
 720
 721
 722     fqRec->Fill(tpcP->Charge());
 723
 724
 725     Double_t sigmapt = tpcP->GetSigma1Pt2();
 726     sigmapt= sqrt(sigmapt);
 727     sigmapt= sigmapt *(tpcP->Pt()*tpcP->Pt());
 728
 729     if(sigmapt == 0.)continue;
 730     fsigmaPt->Fill(TMath::Log10(sigmapt));
 731
 732
 733
 734
 735
 736
 737
 738     // hits in ITS layer
 739     Int_t cas=-1;
 740     if(esdtrack->HasPointOnITSLayer(0))
 741       cas=0;
 742     else if(!esdtrack->HasPointOnITSLayer(0)
 743             &&  esdtrack->HasPointOnITSLayer(1))
 744       cas=1;
 745     else if(!esdtrack->HasPointOnITSLayer(0)
 746             && !esdtrack->HasPointOnITSLayer(1)
 747             &&  esdtrack->HasPointOnITSLayer(2))
 748       cas=2;
 749     else if(!esdtrack->HasPointOnITSLayer(0)
 750             && !esdtrack->HasPointOnITSLayer(1)
 751             && !esdtrack->HasPointOnITSLayer(2)
 752             &&  esdtrack->HasPointOnITSLayer(3))
 753       cas=3;
 754     else if(!esdtrack->HasPointOnITSLayer(0)
 755             && !esdtrack->HasPointOnITSLayer(1)
 756             && !esdtrack->HasPointOnITSLayer(2)
 757             && !esdtrack->HasPointOnITSLayer(3)
 758             &&  esdtrack->HasPointOnITSLayer(4))
 759       cas=4;
 760     else if(   !esdtrack->HasPointOnITSLayer(0)
 761             && !esdtrack->HasPointOnITSLayer(1)
 762             && !esdtrack->HasPointOnITSLayer(2)
 763             && !esdtrack->HasPointOnITSLayer(3)
 764             && !esdtrack->HasPointOnITSLayer(4)
 765             &&  esdtrack->HasPointOnITSLayer(5))
 766       cas=5;
 767     else
 768       cas=6;
 769
 770
 771
 772     //-------------------
 773
 774     sdca_tr = (tpcP->Py()*tpcP->Xv()
 775                - tpcP->Px()*tpcP->Yv())/tpcP->Pt();
 776
 777
 778
 779     fqPtRec[cas]->Fill(tpcP->Charge()/tpcP->Pt());
 780
 781
 782     //    test->Fill(tpcP->Pt(),tpcP->Phi(), tpcP->Theta(),
 783     //         tpcP->Xv(),tpcP->Yv(), tpcP->Zv(), tpcP->Charge());
 784     fHistRECpt->Fill(tpcP->Pt());
 785     fEta->Fill(tpcP->Eta());
 786     fThetaRec->Fill(tpcP->Theta());
 787     fPhiRec->Fill(tpcP->Phi());
 788     fVx->Fill(tpcP->Xv());
 789     fVy->Fill(tpcP->Yv());
 790     fVz->Fill(tpcP->Zv());
 791
 792
 793     fEtaPt->Fill(tpcP->Eta()/tpcP->Pt());
 794     fQPt->Fill(tpcP->Charge()/tpcP->Pt());
 795     fDca->Fill(sdca_tr);
 796
 797
 798
 799
 800     esdtrack->GetImpactParameters(xy,z);
 801     fDiffDcaD->Fill(sdca_tr+xy);
 802
 803
 804     //for positive particles
 805     if(tpcP->Charge()>0){
 806       fRecPtPos->Fill(tpcP->Pt());
 807       fRecPtPosLadder[cas]->Fill(tpcP->Pt());
 808       fRecPtPosVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
 809       fRecPhiPos->Fill(TMath::RadToDeg()*tpcP->Phi());
 810
 811
 812       fRecPhiPosLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
 813       fRecPhiPosVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
 814       fSignedDcaPosVz->Fill(sdca_tr,tpcP->Zv());
 815
 816       fRecEtaPos->Fill(tpcP->Eta());
 817       fRecEtaPosLadder[cas]->Fill(tpcP->Eta());
 818       fRecEtaPtPos->Fill(tpcP->Eta()/tpcP->Pt());
 819       fRecEtaPosVz->Fill(tpcP->Eta(),tpcP->Zv());
 820       fRecEtaPtPosVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
 821
 822       fRecDcaPos->Fill(sdca_tr);
 823       fRecDPos->Fill(xy);
 824       fSignDcaPos[cas]->Fill(sdca_tr);
 825
 826
 827       fDcaSigmaPos[cas]->Fill(sdca_tr, TMath::Log10(sigmapt));
 828
 829       fPtSigmaPos[cas]->Fill(TMath::Log10(sigmapt));
 830     }
 831     //and negative particles
 832     else {
 833       fRecPtNeg->Fill(tpcP->Pt());
 834       fRecPtNegLadder[cas]->Fill(tpcP->Pt());
 835       fRecPtNegVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
 836
 837       fRecPhiNeg->Fill(TMath::RadToDeg()*tpcP->Phi());
 838       fRecPhiNegLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
 839       fRecPhiNegVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
 840       fSignedDcaNegVz->Fill(sdca_tr,tpcP->Zv());
 841       fRecEtaPtNegVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
 842
 843       fRecEtaNeg->Fill(tpcP->Eta());
 844       fRecEtaNegLadder[cas]->Fill(tpcP->Eta());
 845       fRecEtaPtNeg->Fill(tpcP->Eta()/tpcP->Pt());
 846       fRecEtaNegVz->Fill(tpcP->Eta(),tpcP->Zv());
 847
 848       fRecDcaNeg->Fill(sdca_tr);
 849       fRecDcaNegInv->Fill(-sdca_tr);
 850       fRecDNeg->Fill(xy);
 851       fSignDcaNeg[cas]->Fill(sdca_tr);
 852       fSignDcaNegInv[cas]->Fill(-sdca_tr);
 853
 854
 855       fDcaSigmaNeg[cas]->Fill(sdca_tr,TMath::Log10(sigmapt));
 856
 857       fPtSigmaNeg[cas]->Fill(TMath::Log10(sigmapt));
 858     }
 859
 860
 861
 862     //all particles with positive eta
 863     if(tpcP->Eta()>0){
 864       fRecQPtPosEta->Fill(tpcP->Charge()/tpcP->Pt());
 865       fRecPtPosEta->Fill(tpcP->Pt());
 866       fRecPhiPosEta->Fill(TMath::RadToDeg()*tpcP->Phi());
 867       fRecQPtPosEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
 868       fRecDcaPosEta->Fill(sdca_tr);
 869       fRecDPosEta->Fill(xy);
 870     }
 871     //all particles with negative eta (and eta==0)
 872     else{
 873       fRecQPtNegEta->Fill(tpcP->Charge()/tpcP->Pt());
 874       fRecPtNegEta->Fill(tpcP->Pt());
 875       fRecPhiNegEta->Fill(TMath::RadToDeg()*tpcP->Phi());
 876       fRecQPtNegEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
 877       fRecDcaNegEta->Fill(sdca_tr);
 878       fRecDNegEta->Fill(xy);
 879
 880     }
 881
 882
 883     //spectren detected by TPC sectors
 884     //pt cut on 1 GeV/c ?!
 885     // if(tpcP->Pt()<1.) continue;
 886     fRecPtTpcSector[Int_t(tpcP->Phi()*
 887                           TMath::RadToDeg()/20)]->Fill(tpcP->Pt());
 888     fRecEtaTpcSector[Int_t(tpcP->Phi()*
 889                           TMath::RadToDeg()/20)]->Fill(tpcP->Eta());
 890     fSignedDcaTpcSector[Int_t(tpcP->Phi()*
 891                           TMath::RadToDeg()/20)]->Fill(sdca_tr);
 892     fRecQPtTpcSector[Int_t(tpcP->Phi()*
 893                           TMath::RadToDeg()/20)]->Fill(tpcP->Charge()/tpcP->Pt());
 894     fRecEtaPtTpcSector[Int_t(tpcP->Phi()*
 895                           TMath::RadToDeg()/20)]->Fill(tpcP->Eta()/tpcP->Pt());
 896
 897
 898
 899
 900
 901
 902   // another track loop
 903     for (Int_t iTrack2 = 0; iTrack2 < event->GetNumberOfTracks(); iTrack2++) {
 904
 905       if(LeadingTrack==iTrack2) continue;
 906
 907       AliVParticle *track2 = event->GetTrack(iTrack2);
 908       AliESDtrack* esdtrack2 =  dynamic_cast<AliESDtrack*>(track2);
 909       if (!track2) {
 910         Printf("ERROR: Could not receive track %d", iTrack);
 911         continue;
 912       }
 913       if (!fCuts->AcceptTrack(esdtrack2)) continue;
 914       //propagate to dca
 915       esdtrack2->PropagateToDCA(event->GetPrimaryVertex(),
 916                                event->GetMagneticField(), 10000.);
 917
 918       fDeltaPhiLeading->Fill((LeadingPhi-esdtrack2->Phi())*TMath::RadToDeg(),
 919                              LeadingPhi*TMath::RadToDeg() );
 920
 921
 922
 923     }//second track loop
 924   }//first track loop
 925
 926
 927
 928
 929   // Post output data.
 930   // PostData(1, fHistPt);
 931   PostData(1, fHists);
 932 }
 933
 934
 935
 936
 937
 938 //________________________________________________________________________
 939 void AliAnalysisTaskQASym::Terminate(Option_t *)
 940 {
 941
 942
 943 }
 944
 945
 946
 947
 948