PWGLF/RESONANCES/extra/AliRsnAnalysisPhi7TeV.cxx

   1 //
   2 // Implementation file for implementation of data analysis aft 900 GeV
   3 //
   4 // Author: A. Pulvirenti
   5 //
   6
   7 #include "Riostream.h"
   8 #include <iomanip>
   9
  10 #include "TH1.h"
  11 #include "TH3.h"
  12 #include "TFile.h"
  13 #include "TTree.h"
  14 #include "TParticle.h"
  15 #include "TRandom.h"
  16 #include "TLorentzVector.h"
  17 #include "TDatabasePDG.h"
  18
  19 #include "AliLog.h"
  20 #include "AliESDpid.h"
  21 #include "AliESDEvent.h"
  22 #include "AliESDVertex.h"
  23 #include "AliESDtrack.h"
  24 #include "AliESDtrackCuts.h"
  25 #include "AliStack.h"
  26 #include "AliMCEvent.h"
  27 #include "AliTOFT0maker.h"
  28 #include "AliTOFcalib.h"
  29 #include "AliCDBManager.h"
  30 #include "AliITSPIDResponse.h"
  31 #include "AliCFContainer.h"
  32 #include "AliVEventHandler.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliMultiInputEventHandler.h"
  35 #include "AliAnalysisManager.h"
  36
  37 #include "AliRsnEvent.h"
  38 #include "AliRsnTarget.h"
  39 #include "AliRsnDaughter.h"
  40
  41 #include "AliRsnAnalysisPhi7TeV.h"
  42
  43 //__________________________________________________________________________________________________
  44 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name, Bool_t isMC) :
  45    AliAnalysisTaskSE(name),
  46    fIsMC(isMC),
  47    fAddTPC(kTRUE),
  48    fAddITS(kFALSE),
  49    fCheckITS(kTRUE),
  50    fCheckTPC(kTRUE),
  51    fCheckTOF(kTRUE),
  52    fStep(1000),
  53    fPhiMass(1.019455),
  54    fKaonMass(0.49677),
  55    fMaxVz(10.0),
  56    fITSNSigma(3.0),
  57    fTPCMomentumThreshold(0.350),
  58    fTOFNSigma(3.0),
  59    fESDtrackCutsTPC(0),
  60    fESDtrackCutsITS(0),
  61    fESDpid(0),
  62    fOutList(0x0),
  63    fCFunlike(0),
  64    fCFlikePP(0),
  65    fCFlikeMM(0),
  66    fCFtrues(0),
  67    fCFkaons(0),
  68    fHEvents(0x0)
  69 {
  70 //
  71 // Constructor
  72 //
  73
  74    fTPCNSigma[0] = 5.0;
  75    fTPCNSigma[1] = 3.0;
  76
  77    fVertexX[0] = fVertexX[1] = 0;
  78    fVertexY[0] = fVertexZ[1] = 0;
  79    fVertexZ[0] = fVertexZ[1] = 0;
  80
  81    DefineOutput(1, TList::Class());
  82 }
  83
  84 //__________________________________________________________________________________________________
  85 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) :
  86    AliAnalysisTaskSE(copy),
  87    fIsMC(copy.fIsMC),
  88    fAddTPC(copy.fAddTPC),
  89    fAddITS(copy.fAddITS),
  90    fCheckITS(copy.fCheckITS),
  91    fCheckTPC(copy.fCheckTPC),
  92    fCheckTOF(copy.fCheckTOF),
  93    fStep(copy.fStep),
  94    fPhiMass(copy.fPhiMass),
  95    fKaonMass(copy.fKaonMass),
  96    fMaxVz(copy.fMaxVz),
  97    fITSNSigma(copy.fITSNSigma),
  98    fTPCMomentumThreshold(copy.fTPCMomentumThreshold),
  99    fTOFNSigma(copy.fTOFNSigma),
 100    fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
 101    fESDtrackCutsITS(copy.fESDtrackCutsITS),
 102    fESDpid(copy.fESDpid),
 103    fOutList(0x0),
 104    fCFunlike(0),
 105    fCFlikePP(0),
 106    fCFlikeMM(0),
 107    fCFtrues(0),
 108    fCFkaons(0),
 109    fHEvents(0x0)
 110 {
 111 //
 112 // Copy constructor
 113 //
 114
 115    fTPCNSigma[0] = copy.fTPCNSigma[0];
 116    fTPCNSigma[1] = copy.fTPCNSigma[1];
 117
 118    fVertexX[0] = fVertexX[1] = 0;
 119    fVertexY[0] = fVertexZ[1] = 0;
 120    fVertexZ[0] = fVertexZ[1] = 0;
 121 }
 122
 123 //__________________________________________________________________________________________________
 124 AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7TeV& copy)
 125 {
 126 //
 127 // Assignment operator
 128 //
 129   if (this == &copy)
 130     return *this;
 131    fIsMC = copy.fIsMC;
 132    fAddTPC = copy.fAddTPC;
 133    fAddITS = copy.fAddITS;
 134    fCheckITS = copy.fCheckITS;
 135    fCheckTPC = copy.fCheckTPC;
 136    fCheckTOF = copy.fCheckTOF;
 137    fStep = copy.fStep;
 138    fPhiMass = copy.fPhiMass;
 139    fKaonMass = copy.fKaonMass;
 140    fMaxVz = copy.fMaxVz;
 141    fITSNSigma = copy.fITSNSigma;
 142    fTPCNSigma[0] = copy.fTPCNSigma[0];
 143    fTPCNSigma[1] = copy.fTPCNSigma[1];
 144    fTPCMomentumThreshold = copy.fTPCMomentumThreshold;
 145    fTOFNSigma = copy.fTOFNSigma;
 146    fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
 147    fESDtrackCutsITS = copy.fESDtrackCutsITS;
 148    fESDpid = copy.fESDpid;
 149
 150    return (*this);
 151 }
 152
 153 //__________________________________________________________________________________________________
 154 AliRsnAnalysisPhi7TeV::~AliRsnAnalysisPhi7TeV()
 155 {
 156 //
 157 // Destructor
 158 //
 159 }
 160
 161 //__________________________________________________________________________________________________
 162 void AliRsnAnalysisPhi7TeV::UserCreateOutputObjects()
 163 {
 164 //
 165 // Create the output data container
 166 //
 167
 168    // create output list
 169    OpenFile(1);
 170    fOutList = new TList;
 171
 172    // event related histograms
 173    fHEvents    = new TH1I("hEvents", "Event details", kVertexTypes, 0, kVertexTypes);
 174    fVertexX[0] = new TH1F("hVertexTRKX", "X position of primary vertex (tracks)", 200,  -2,  2);
 175    fVertexY[0] = new TH1F("hVertexTRKY", "Y position of primary vertex (tracks)", 200,  -2,  2);
 176    fVertexZ[0] = new TH1F("hVertexTRKZ", "Z position of primary vertex (tracks)", 400, -20, 20);
 177    fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)"   , 200,  -2,  2);
 178    fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)"   , 200,  -2,  2);
 179    fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)"   , 400, -20, 20);
 180
 181    fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks");
 182    fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD");
 183    fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks");
 184    fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD");
 185    fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex");
 186    fHEvents->GetXaxis()->SetBinLabel(6, "Empty event");
 187
 188    // define axes:
 189    // pairs [0] = inv. mass
 190    //       [1] = inv. mass resolution
 191    //       [2] = pt
 192    //       [3] = rec. multiplicity
 193    //       [4] = MC multiplicity (if available)
 194    //
 195    // kaons [0] = pt
 196    //       [1] = pTPC
 197    //       [2] = DCAr
 198    //       [3] = DCAz
 199    //       [4] = signalITS
 200    //       [5] = signalTPC
 201    //       [6] = betaTOF
 202    //       [7] = rec. multiplicity
 203    //       [8] = MC multiplicity (if available)
 204    Double_t minIM     =  0.9, maxIM    =   1.4;
 205    Double_t minIMRes  = -5.0, maxIMRes =   5.0;
 206    Double_t minMom    =  0.0, maxMom   =   5.0;
 207    Double_t minDCAr   =  0.0, maxDCAr  =   1.0;
 208    Double_t minDCAz   =  0.0, maxDCAz  =   5.0;
 209    Double_t minITS    =  0.0, maxITS   = 800.0; // raw signal
 210    Double_t minTPC    =  0.0, maxTPC   = 800.0; // raw signal
 211    Double_t minTOF    =  0.0, maxTOF   =   1.0; // beta = length / TOF time / c
 212    Double_t mult[]    =  {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 30., 35., 40., 50., 60., 70., 80., 90., 100., 120., 140., 160., 180., 500.};
 213    Int_t    nIM       =  500;
 214    Int_t    nIMRes    =  100;
 215    Int_t    nMom      =   50;
 216    Int_t    nDCAr     =  100;
 217    Int_t    nDCAz     =  100;
 218    Int_t    nITS      =  800;
 219    Int_t    nTPC      =  800;
 220    Int_t    nTOF      =  100;
 221    Int_t    nMult     =  sizeof(mult) / sizeof(mult[0]) - 1;
 222    Int_t    nBins1[9] =  {nMom, nMom, nDCAr, nDCAz, nITS, nTPC, nTOF, nMult, nMult};
 223    Int_t    nBins2[5] =  {nIM , nIMRes, nMom, nMult, nMult};
 224
 225    // containers for analysis pairs have 2 steps (quality, quality+PID)
 226    // and those for efficiency have 3 (add full MC)
 227    // REMINDER: to use variable bins, call container->SetBinLimits(iaxis, Double_t *array);
 228    fCFunlike = new AliCFContainer("CFUnlike", "", 2, 5, nBins2);
 229    fCFlikePP = new AliCFContainer("CFLikePP", "", 2, 5, nBins2);
 230    fCFlikeMM = new AliCFContainer("CFLikeMM", "", 2, 5, nBins2);
 231    fCFtrues  = new AliCFContainer("CFtrues" , "", 3, 5, nBins2);
 232    fCFkaons  = new AliCFContainer("CFkaons" , "", 3, 9, nBins1);
 233
 234    // create a unique temp array to all pair containers
 235    AliCFContainer *CFpair[4];
 236    CFpair[0] = fCFunlike;
 237    CFpair[1] = fCFlikePP;
 238    CFpair[2] = fCFlikeMM;
 239    CFpair[3] = fCFtrues;
 240
 241    // define axes:
 242    // pairs [0] = inv. mass
 243    //       [1] = inv. mass resolution
 244    //       [2] = pt
 245    //       [3] = rec. multiplicity
 246    //       [4] = MC multiplicity (if available)
 247    //
 248    // kaons [0] = pt
 249    //       [1] = pTPC
 250    //       [2] = DCAr
 251    //       [3] = DCAz
 252    //       [4] = signalITS
 253    //       [5] = signalTPC
 254    //       [6] = betaTOF
 255    //       [7] = rec. multiplicity
 256    //       [8] = MC multiplicity (if available)
 257    Int_t i;
 258    for (i = 0; i < 4; i++) {
 259       CFpair[i]->SetBinLimits(0, minIM   , maxIM   );
 260       CFpair[i]->SetBinLimits(1, minIMRes, maxIMRes);
 261       CFpair[i]->SetBinLimits(2, minMom  , maxMom  );
 262       CFpair[i]->SetBinLimits(3, mult);
 263       CFpair[i]->SetBinLimits(4, mult);
 264    }
 265    fCFkaons->SetBinLimits(0, minMom , maxMom );
 266    fCFkaons->SetBinLimits(1, minMom , maxMom );
 267    fCFkaons->SetBinLimits(2, minDCAr, maxDCAr);
 268    fCFkaons->SetBinLimits(3, minDCAz, maxDCAz);
 269    fCFkaons->SetBinLimits(4, minITS , maxITS );
 270    fCFkaons->SetBinLimits(5, minTPC , maxTPC );
 271    fCFkaons->SetBinLimits(6, minTOF , maxTOF );
 272    fCFkaons->SetBinLimits(7, mult);
 273    fCFkaons->SetBinLimits(8, mult);
 274
 275    // add all outputs to the list
 276    fOutList->Add(fHEvents);
 277    fOutList->Add(fVertexX[0]);
 278    fOutList->Add(fVertexY[0]);
 279    fOutList->Add(fVertexZ[0]);
 280    fOutList->Add(fVertexX[1]);
 281    fOutList->Add(fVertexY[1]);
 282    fOutList->Add(fVertexZ[1]);
 283    fOutList->Add(fCFunlike);
 284    fOutList->Add(fCFlikePP);
 285    fOutList->Add(fCFlikeMM);
 286    fOutList->Add(fCFtrues);
 287    fOutList->Add(fCFkaons);
 288
 289    // update histogram container
 290    PostData(1, fOutList);
 291 }
 292
 293 //__________________________________________________________________________________________________
 294 void AliRsnAnalysisPhi7TeV::UserExec(Option_t *)
 295 {
 296 //
 297 // Main execution function.
 298 // Fills the fHEvents data member with the following legenda:
 299 // 0 -- event OK, prim vertex with tracks
 300 // 1 -- event OK, prim vertex with SPD
 301 // 2 -- event OK but vz large
 302 // 3 -- event bad
 303 //
 304
 305    static Int_t evNum = 0;
 306    evNum++;
 307    if (evNum % fStep == 0) AliInfo(Form("Processing event #%7d", evNum));
 308
 309    // retrieve ESD event and related stack (if available)
 310    AliESDEvent *esd = static_cast<AliESDEvent*>(fInputEvent);
 311    AliMCEvent  *mc  = static_cast<AliMCEvent*>(fMCEvent);
 312
 313    // if the AliESDpid object is not assigned, do this now
 314    // find an AliESDpid from input handler
 315    if (!fESDpid) {
 316       AliAnalysisManager *mgr   = AliAnalysisManager::GetAnalysisManager();
 317       AliVEventHandler   *evh   = mgr->GetInputEventHandler();
 318       AliESDInputHandler *esdin = 0x0;
 319       if (evh->IsA() == AliMultiInputEventHandler::Class()) {
 320          AliMultiInputEventHandler *multh = static_cast<AliMultiInputEventHandler*>(evh);
 321          AliInputEventHandler      *input = multh->GetFirstInputEventHandler();
 322          if (input->IsA() == AliESDInputHandler::Class()) {
 323             esdin = static_cast<AliESDInputHandler*>(input);
 324          }
 325       }
 326       else if (evh->IsA() == AliESDInputHandler::Class()) {
 327          esdin = static_cast<AliESDInputHandler*>(evh);
 328       }
 329       if (!esdin) {
 330          AliFatal("Required ESD input handler");
 331          return;
 332       }
 333       fESDpid = esdin->GetESDpid();
 334    }
 335
 336    // check the event
 337    EVertexType eval = EventEval(esd);
 338    fHEvents->Fill((Int_t)eval);
 339
 340    // if the event is good for analysis, process it
 341    if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex) {
 342       ProcessESD(esd, mc);
 343       if (mc) ProcessMC (esd, mc);
 344    }
 345
 346    // update histogram container
 347    PostData(1, fOutList);
 348 }
 349
 350 //__________________________________________________________________________________________________
 351 void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
 352 {
 353 //
 354 // Terminate
 355 //
 356 }
 357
 358 //__________________________________________________________________________________________________
 359 AliRsnAnalysisPhi7TeV::EVertexType AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
 360 {
 361 //
 362 // Checks if the event is good for analysis.
 363 // Returns one of the flag values defined in the header
 364 //
 365
 366    // reject empty events
 367    Int_t ntracks = esd->GetNumberOfTracks();
 368    if (!ntracks) return kEmptyEvent;
 369
 370    // get the best primary vertex:
 371    // first try the one with tracks
 372    const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 373    const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 374    Double_t            vzTrk = 1000000.0;
 375    Double_t            vzSPD = 1000000.0;
 376    Int_t               ncTrk = -1;
 377    Int_t               ncSPD = -1;
 378    if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
 379    if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
 380    if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
 381    if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
 382    if (vTrk && ncTrk > 0) {
 383       // fill the histograms
 384       fVertexX[0]->Fill(vTrk->GetXv());
 385       fVertexY[0]->Fill(vTrk->GetYv());
 386       fVertexZ[0]->Fill(vTrk->GetZv());
 387
 388       // check VZ position
 389       if (vzTrk <= fMaxVz)
 390          return kGoodTracksPrimaryVertex;
 391       else
 392          return kFarTracksPrimaryVertex;
 393    } else if (vSPD && ncSPD > 0) {
 394       // fill the histograms
 395       fVertexX[1]->Fill(vSPD->GetXv());
 396       fVertexY[1]->Fill(vSPD->GetYv());
 397       fVertexZ[1]->Fill(vSPD->GetZv());
 398
 399       // check VZ position
 400       if (vzSPD <= fMaxVz)
 401          return kGoodSPDPrimaryVertex;
 402       else
 403          return kFarSPDPrimaryVertex;
 404    } else
 405       return kNoGoodPrimaryVertex;
 406 }
 407
 408 //__________________________________________________________________________________________________
 409 void AliRsnAnalysisPhi7TeV::ProcessESD
 410 (AliESDEvent *esd, AliMCEvent *mc)
 411 {
 412 //
 413 // This function works with the ESD object
 414 // and fills all containers
 415 //
 416
 417    // some utility arrays, declared static
 418    // to avoid creting/destroying at each event
 419    static TArrayC quality; // will be 1 or 0 depending if cut is passed or not
 420    static TArrayC pid;     // will be 1 or 0 depending if cut is passed or not
 421
 422    // get stack, if possible
 423    AliStack *stack = 0x0;
 424    if (mc) stack = mc->Stack();
 425
 426    // get number of tracks, which is used for setting arrays
 427    Int_t itrack, ntracks = esd->GetNumberOfTracks();
 428    quality.Set(ntracks);
 429    pid    .Set(ntracks);
 430
 431    // loop on tracks and sets the flags
 432    for (itrack = 0; itrack < ntracks; itrack++) {
 433
 434       // get track
 435       AliESDtrack *track = esd->GetTrack(itrack);
 436
 437       // check cuts
 438       quality[itrack] = (Char_t)OkQuality(track);
 439       pid    [itrack] = (Char_t)OkPID(track);
 440    }
 441
 442    // now, we can fill the histograms
 443    Int_t i1, i2, label;
 444    Double_t var[5];
 445    TLorentzVector p1[2], p2[2], sum[2];
 446    AliESDtrack *tr1, *tr2;
 447    TParticle   *pr1, *pr2;
 448    ETrackType type1, type2;
 449
 450    // multiplicity is defined for whole event
 451    var[3] = AliESDtrackCuts::GetReferenceMultiplicity(esd, kTRUE);
 452    var[4] = (mc ? mc->GetNumberOfTracks() : 0.0);
 453
 454    for (i1 = 0; i1 < ntracks; i1++) {
 455
 456       // check only quality
 457       if (!quality[i1]) continue;
 458
 459       // skip tracks which are not of a good type
 460       tr1 = esd->GetTrack(i1);
 461       type1 = TrackType(tr1);
 462       if (type1 >= kTrackTypes) continue;
 463
 464       // skip neutral
 465       tr1 = esd->GetTrack(i1);
 466       if (tr1->Charge() == 0) continue;
 467
 468       // try to retrieve the MC
 469       pr1 = 0x0;
 470       label = TMath::Abs(tr1->GetLabel());
 471       if (stack) {
 472          if (label >= 0 && label < stack->GetNtrack()) {
 473             pr1 = stack->Particle(label);
 474          }
 475       }
 476
 477       // initialize 4-momenta
 478       p1[0].SetXYZM(tr1->Px(), tr1->Py(), tr1->Pz(), fKaonMass);
 479       if (pr1) p1[1].SetXYZM(pr1->Px(), pr1->Py(), pr1->Pz(), fKaonMass);
 480
 481       for (i2 = i1+1; i2 < ntracks; i2++) {
 482
 483          // check only quality
 484          if (!quality[i2]) continue;
 485
 486          // skip tracks which are not of a good type
 487          tr2 = esd->GetTrack(i2);
 488          type2 = TrackType(tr2);
 489          if (type2 >= kTrackTypes) continue;
 490
 491          // skip neutrals
 492          tr2 = esd->GetTrack(i2);
 493          if (tr2->Charge() == 0) continue;
 494
 495          // try to retrieve MC
 496          pr2 = 0x0;
 497          label = TMath::Abs(tr2->GetLabel());
 498          if (stack) {
 499             if (label >= 0 && label < stack->GetNtrack()) {
 500                pr2 = stack->Particle(label);
 501             }
 502          }
 503
 504          // initialize 4-momenta
 505          p2[0].SetXYZM(tr2->Px(), tr2->Py(), tr2->Pz(), fKaonMass);
 506          if (pr2) p1[1].SetXYZM(pr2->Px(), pr2->Py(), pr2->Pz(), fKaonMass);
 507
 508          // sum 4-momenta
 509          for (Int_t i = 0; i < 2; i++) sum[i] = p1[i] + p2[i];
 510
 511          // compute values
 512          // [0] = inv. mass
 513          // [1] = inv. mass resolution
 514          // [2] = pt
 515          // [3] = rec. multiplicity
 516          // [4] = MC multiplicity (if available)
 517          var[0] = sum[0].M();
 518          if (mc) var[1] = (sum[1].M() - sum[0].M()) / sum[1].M();
 519          var[2] = sum[0].Perp();
 520
 521          // fill containers
 522          AliCFContainer *cf = 0x0;
 523          if (tr1->Charge() != tr2->Charge())
 524             cf = fCFunlike;
 525          else if (tr1->Charge() > 0)
 526             cf = fCFlikePP;
 527          else
 528             cf = fCFlikeMM;
 529
 530          // fill steps: 0 = quality onl1, 1 = quality + PID
 531          cf->Fill(var, 0);
 532          if (pid[i1] && pid[i2]) cf->Fill(var, 1);
 533       }
 534    }
 535 }
 536
 537 //__________________________________________________________________________________________________
 538 void AliRsnAnalysisPhi7TeV::ProcessMC(AliESDEvent *esd, AliMCEvent *mc)
 539 {
 540 //
 541 // Searches true pairs and fill efficiency for kaons
 542 //
 543
 544    // get stack, if possible
 545    AliStack *stack = 0x0;
 546    if (mc) stack = mc->Stack(); else return;
 547
 548    // count tracks
 549    Int_t ipart, itrack, npart = stack->GetNprimary();
 550    AliESDtrack *matched[2] = {0x0, 0x0};
 551    Double_t var1[10], var2[5];
 552    Float_t  b[2], bcov[3];
 553
 554    var1[7] = var2[3] = AliESDtrackCuts::GetReferenceMultiplicity(esd, kTRUE);
 555    var1[8] = var2[4] = npart;
 556
 557    // loop on kaons
 558    for (ipart = 0; ipart < npart; ipart++) {
 559       TParticle *part = stack->Particle(ipart);
 560       if (!stack->IsPhysicalPrimary(ipart)) continue;
 561       if (TMath::Abs(part->GetPdgCode()) != 321) continue;
 562
 563       // search best reconstructed track
 564       // which matches that particle ID
 565       // and checks it against cuts
 566       matched[0] = 0x0;
 567       itrack = BestMatchedTrack(esd, ipart);
 568       if (itrack >= 0) matched[0] = esd->GetTrack(itrack);
 569
 570       // compute all variables fo single-tracks
 571       // [0] = pt
 572       // [1] = pTPC
 573       // [2] = DCAr
 574       // [3] = DCAz
 575       // [4] = signalITS
 576       // [5] = signalTPC
 577       // [6] = betaTOF
 578       // [7] = rec. multiplicity
 579       // [8] = MC multiplicity (if available)
 580       var1[0] = part->Pt();
 581       var1[1] = -1.0;
 582       for (Int_t i = 0; i < 6; i++) var1[i] = -1.0;
 583       if (matched[0]) {
 584          matched[0]->GetImpactParameters(b, bcov);
 585          if (IsTPC(matched[0])) var1[1] = matched[0]->GetInnerParam()->P();
 586          var1[2] = TMath::Abs((Double_t)b[0]);
 587          var1[3] = TMath::Abs((Double_t)b[1]);
 588          var1[4] = matched[0]->GetITSsignal();
 589          var1[5] = matched[0]->GetTPCsignal();
 590          if (MatchTOF(matched[0])) {
 591             var1[6]  = matched[0]->GetIntegratedLength();
 592             var1[6] /= matched[0]->GetTOFsignal();
 593             var1[6] /= 0.299792458E8;
 594          }
 595       }
 596
 597       // fill container
 598       fCFkaons->Fill(var1, 0);
 599
 600       // replace MC momentum with rec
 601       if (matched[0]) {
 602          var1[0] = matched[0]->Pt();
 603          if (OkQuality(matched[0])) {
 604             fCFkaons->Fill(var1, 1);
 605             if (OkPID(matched[0])) fCFkaons->Fill(var1, 2);
 606          }
 607       }
 608    }
 609
 610    // loop on phi's
 611    Int_t id1, id2;
 612    TParticle *d1, *d2;
 613    TLorentzVector p1[2], p2[2], sum[2];
 614    for (ipart = 0; ipart < npart; ipart++) {
 615       TParticle *part = stack->Particle(ipart);
 616
 617       // check that particle is a phi
 618       // and that has decayed into kaons
 619       if (TMath::Abs(part->GetPdgCode()) != 333) continue;
 620       id1 = part->GetFirstDaughter();
 621       id2 = part->GetLastDaughter();
 622       d1 = d2 = 0x0;
 623       if (id1 >= 0 && id1 < npart) d1 = stack->Particle(id1);
 624       if (id2 >= 0 && id2 < npart) d2 = stack->Particle(id2);
 625       if (!d1 || !d2) continue;
 626       if (TMath::Abs(d1->GetPdgCode()) != 321) continue;
 627       if (TMath::Abs(d2->GetPdgCode()) != 321) continue;
 628
 629       // get matched tracks
 630       matched[0] = matched[1] = 0x0;
 631       itrack = BestMatchedTrack(esd, id1);
 632       if (itrack >= 0) matched[0] = esd->GetTrack(itrack);
 633       itrack = BestMatchedTrack(esd, id2);
 634       if (itrack >= 0) matched[1] = esd->GetTrack(itrack);
 635
 636       // compute values
 637       // [0] = inv. mass
 638       // [1] = inv. mass resolution
 639       // [2] = pt
 640       // [3] = rec. multiplicity
 641       // [4] = MC multiplicity (if available)
 642       p1[0].SetXYZM(d1->Px(), d1->Py(), d1->Pz(), fKaonMass);
 643       p2[0].SetXYZM(d2->Px(), d2->Py(), d2->Pz(), fKaonMass);
 644       p1[1].SetXYZM(0.0, 0.0, 0.0, 0.0);
 645       p2[1].SetXYZM(0.0, 0.0, 0.0, 0.0);
 646       if (matched[0]) p1[1].SetXYZM(matched[0]->Px(), matched[0]->Py(), matched[0]->Pz(), fKaonMass);
 647       if (matched[1]) p2[1].SetXYZM(matched[1]->Px(), matched[1]->Py(), matched[1]->Pz(), fKaonMass);
 648       for (Int_t i = 0; i < 2; i++) sum[i] = p1[i] + p2[i];
 649       var2[0] = sum[0].M();
 650       var2[1] = (sum[0].M() - sum[1].M()) / sum[0].M();
 651       var2[2] = sum[0].Perp();
 652
 653       // fill container
 654       fCFtrues->Fill(var2, 0);
 655
 656       // replace MC momenta with rec
 657       if (matched[0] && matched[1]) {
 658          var2[0] = sum[1].M();
 659          var2[2] = sum[1].Perp();
 660          if (OkQuality(matched[0]) && OkQuality(matched[1])) {
 661             fCFtrues->Fill(var2, 1);
 662             if (OkPID(matched[0]) && OkPID(matched[1])) fCFtrues->Fill(var2, 2);
 663          }
 664       }
 665    }
 666 }
 667
 668 //______________________________________________________________________________
 669 Bool_t AliRsnAnalysisPhi7TeV::OkPIDITS(AliESDtrack *track, AliPID::EParticleType pid)
 670 {
 671 //
 672 // Check ITS particle identification with 3sigma cut
 673 //
 674
 675    // reject not ITS standalone tracks
 676    if (!IsITS(track)) return kFALSE;
 677
 678    // count PID layers and reject if they are too few
 679    Int_t   k, nITSpidLayers = 0;
 680    UChar_t itsCluMap = track->GetITSClusterMap();
 681    for (k = 2; k < 6; k++) if (itsCluMap & (1 << k)) ++nITSpidLayers;
 682    if (nITSpidLayers < 3) {
 683       AliDebug(AliLog::kDebug + 2, "Rejecting track with too few ITS pid layers");
 684       return kFALSE;
 685    }
 686
 687    // check the track type (ITS+TPC or ITS standalone)
 688    // and reject it if it is of none of the allowed types
 689    Bool_t isSA = kFALSE;
 690    if (IsTPC(track)) isSA = kFALSE;
 691    else if (IsITS(track)) isSA = kTRUE;
 692    else {
 693       AliWarning("Track is neither ITS+TPC nor ITS standalone");
 694       return kFALSE;
 695    }
 696
 697    // create the PID response object and compute nsigma
 698    AliITSPIDResponse &itsrsp = fESDpid->GetITSResponse();
 699    Double_t mom    = track->P();
 700    Double_t nSigma = itsrsp.GetNumberOfSigmas(mom, track->GetITSsignal(), pid, nITSpidLayers, isSA);
 701
 702    // evaluate the cut
 703    Bool_t ok = (TMath::Abs(nSigma) <= fITSNSigma);
 704
 705    // debug message
 706    AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", nSigma, fITSNSigma, (ok ? "passed" : "failed")));
 707
 708    // outcome
 709    return ok;
 710 }
 711
 712 //______________________________________________________________________________
 713 Bool_t AliRsnAnalysisPhi7TeV::OkPIDTPC(AliESDtrack *track, AliPID::EParticleType pid)
 714 {
 715 //
 716 // Check TPC particle identification with {3|5}sigmacut,
 717 // depending on the track total momentum.
 718 //
 719
 720    // reject not TPC tracks
 721    if (!IsTPC(track)) return kFALSE;
 722    if (!track->GetInnerParam()) return kFALSE;
 723
 724    // setup TPC PID response
 725    AliTPCPIDResponse &tpcrsp = fESDpid->GetTPCResponse();
 726
 727    // get momentum and number of sigmas and choose the reference band
 728    Double_t mom       = track->GetInnerParam()->P();
 729    Double_t nSigma    = tpcrsp.GetNumberOfSigmas(mom, track->GetTPCsignal(), track->GetTPCsignalN(), pid);
 730    Double_t maxNSigma = fTPCNSigma[0];
 731    if (mom > fTPCMomentumThreshold) maxNSigma = fTPCNSigma[1];
 732
 733    // evaluate the cut
 734    Bool_t ok = (TMath::Abs(nSigma) <= maxNSigma);
 735
 736    // debug message
 737    AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", nSigma, maxNSigma, (ok ? "passed" : "failed")));
 738
 739    // outcome
 740    return ok;
 741 }
 742
 743 //______________________________________________________________________________
 744 Bool_t AliRsnAnalysisPhi7TeV::OkPIDTOF(AliESDtrack *track, AliPID::EParticleType pid)
 745 {
 746 //
 747 // Check TOF particle identification if matched there.
 748 //
 749
 750    // reject not TOF-matched tracks
 751    if (!MatchTOF(track)) return kFALSE;
 752
 753    // setup TOF PID response
 754    AliTOFPIDResponse &tofrsp = fESDpid->GetTOFResponse();
 755
 756    // get info for computation
 757    Double_t momentum = track->P();
 758    Double_t time     = track->GetTOFsignal();
 759    Double_t timeint[AliPID::kSPECIES];
 760    tofrsp.GetStartTime(momentum);
 761    track->GetIntegratedTimes(timeint);
 762
 763    // check the cut
 764    Double_t timeDiff = time - timeint[(Int_t)pid];
 765    Double_t sigmaRef = tofrsp.GetExpectedSigma(momentum, timeint[(Int_t)pid], AliPID::ParticleMass(pid));
 766    Double_t nSigma   = timeDiff / sigmaRef;
 767
 768    // evaluate the cut
 769    Bool_t ok = (TMath::Abs(nSigma) <= fTOFNSigma);
 770
 771    // debug message
 772    AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- nsigma = %f -- cut %s", nSigma, fTOFNSigma, (ok ? "passed" : "failed")));
 773
 774    // outcome
 775    return ok;
 776 }
 777
 778 //______________________________________________________________________________
 779 Int_t AliRsnAnalysisPhi7TeV::MatchedTrack(AliESDEvent *esd, Int_t label, Int_t &npassed, Int_t start)
 780 {
 781 //
 782 // Starting from 'start' searches a track in ESD which has label 'label',
 783 // and if it is found, records how many cuts it passes:
 784 // 0 = none
 785 // 1 = only quality
 786 // 2 = quality and PID
 787 // Returns the ID of that track (-1 if not found)
 788 //
 789
 790    Int_t it, ntracks = esd->GetNumberOfTracks();
 791    AliESDtrack *track = 0x0;
 792
 793    for (it = start; it < ntracks; it++) {
 794       track = esd->GetTrack(it);
 795       if (TMath::Abs(track->GetLabel() != label)) continue;
 796
 797       npassed = 0;
 798       if (OkQuality(track)) npassed++;
 799       if (OkPID    (track)) npassed++;
 800       return it;
 801    }
 802
 803    return -1;
 804 }
 805
 806 //______________________________________________________________________________
 807 Int_t AliRsnAnalysisPhi7TeV::BestMatchedTrack(AliESDEvent *esd, Int_t label)
 808 {
 809 //
 810 // Searches the best track which matches the specified label
 811 // where 'best' means the one which passes most cuts
 812 //
 813
 814    // first attempt:
 815    // if it fails, there are no matched tracks
 816    Int_t ncuts = 0;
 817    Int_t itrack = MatchedTrack(esd, label, ncuts);
 818    /*if (itrack < 0) return -1;
 819
 820    // if it succeeds, use it as a start for a loop
 821    Int_t bestCuts = ncuts;
 822    Int_t start    = itrack + 1;
 823    for (;;) {
 824       it = MatchedTrack(esd, label, ncuts, start);
 825       if (it < 0) break;
 826
 827       if (ncuts > bestCuts) {
 828          bestCuts = ncuts;
 829          itrack = it;
 830          start = it + 1;
 831       }
 832    }*/
 833
 834    return itrack;
 835 }