prod/acrcaf/qa_pp/AliAnalysisTaskQASym.cxx

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