PWGLF/RESONANCES/extra/AliAnalysisTaskSigma1385.cxx

   1 /**************************************************************************
   2  * Authors : Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it)             *
   3  *      Enrico Fragiacomo (enrico.fragiacomo@ts.infn.it)             *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //-----------------------------------------------------------------
  16 //                 AliAnalysisTaskSigma1385 class
  17 //-----------------------------------------------------------------
  18
  19 class TTree;
  20 class TParticle;
  21 class TVector3;
  22
  23 #include "AliAnalysisManager.h"
  24 #include <AliMCEventHandler.h>
  25 #include <AliMCEvent.h>
  26 #include <AliStack.h>
  27
  28 class AliESDVertex;
  29 class AliESDv0;
  30 class AliAODv0;
  31
  32 #include <iostream>
  33
  34 #include "TList.h"
  35 #include "TH1.h"
  36 #include "TNtuple.h"
  37 #include "TGraph.h"
  38 #include "TCanvas.h"
  39 #include "TMath.h"
  40 #include "TChain.h"
  41 #include "AliLog.h"
  42 #include "AliESDEvent.h"
  43 #include "AliAODEvent.h"
  44 #include "AliCascadeVertexer.h"
  45 #include "AliESDcascade.h"
  46 #include "AliAODcascade.h"
  47 #include "AliAnalysisTaskSigma1385.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliESDpid.h"
  50 #include "AliTOFT0maker.h"
  51 #include "AliTOFcalib.h"
  52 #include "AliCDBManager.h"
  53 #include "AliESDtrackCuts.h"
  54
  55 ClassImp(AliAnalysisTaskSigma1385)
  56
  57 //________________________________________________________________________
  58 AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385()
  59    : AliAnalysisTaskSE(),
  60
  61      //-------------------------------- For PID
  62
  63      fisMC(0),
  64      fIsMC(fisMC),
  65      fCheckITS(kTRUE),
  66      fCheckTPC(kTRUE),
  67      fCheckTOF(kTRUE),
  68      fUseGlobal(kTRUE),
  69      fUseITSSA(kTRUE),
  70      fMaxITSband(3.0),
  71      fTPCpLimit(0.35),
  72      fTPCpar(),
  73      fMinTPCband(3.0),
  74      fMaxTPCband(5.0),
  75      fESDpid(0x0),
  76      fTOFmaker(0x0),
  77      fTOFcalib(0x0),
  78      fTOFcalibrateESD(!fisMC),
  79      fTOFcorrectTExp(kTRUE),
  80      fTOFuseT0(kTRUE),
  81      fTOFtuneMC(fisMC),
  82      fTOFresolution(100.0),
  83      fMinTOF(-2.5),
  84      fMaxTOF(3.5),
  85      fLastRun(-1),
  86      fOkTrack(),
  87      fAnalysisType("ESD"), fCollidingSystems(0), fDataType("REAL"), fListHistCascade(0),
  88      fHistEventMultiplicity(0), fHistEventMultiplicityRAVS(0),
  89      fNtuple1(0), fNtuple2(0), fNtuple3(0), fNtuple4(0)
  90
  91
  92
  93      //--------------------------------
  94
  95 {
  96    // Dummy Constructor
  97    Int_t i;
  98    for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
  99 }
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385(const char *name)
 103    : AliAnalysisTaskSE(name),
 104
 105
 106      //-------------------------------- For PID
 107
 108      fisMC(0),
 109      fIsMC(fisMC),
 110      fCheckITS(kTRUE),
 111      fCheckTPC(kTRUE),
 112      fCheckTOF(kTRUE),
 113      fUseGlobal(kTRUE),
 114      fUseITSSA(kTRUE),
 115      fMaxITSband(3.0),
 116      fTPCpLimit(0.35),
 117      fTPCpar(),
 118      fMinTPCband(3.0),
 119      fMaxTPCband(5.0),
 120      fESDpid(0x0),
 121      fTOFmaker(0x0),
 122      fTOFcalib(0x0),
 123      fTOFcalibrateESD(!fisMC),
 124      fTOFcorrectTExp(kTRUE),
 125      fTOFuseT0(kTRUE),
 126      fTOFtuneMC(fisMC),
 127      fTOFresolution(100.0),
 128      fMinTOF(-2.5),
 129      fMaxTOF(3.5),
 130      fLastRun(-1),
 131      fOkTrack(),
 132      fAnalysisType("ESD"), fCollidingSystems(0), fDataType("REAL"), fListHistCascade(0),
 133      fHistEventMultiplicity(0), fHistEventMultiplicityRAVS(0),
 134      fNtuple1(0), fNtuple2(0), fNtuple3(0), fNtuple4(0)
 135
 136      //--------------------------------
 137 {
 138
 139    // Output slot #0 writes into a TList container (Cascade)
 140    DefineOutput(1, TList::Class());
 141
 142    Int_t i;
 143    for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
 144 }
 145
 146 //________________________________________________________________________
 147 void AliAnalysisTaskSigma1385::UserCreateOutputObjects()
 148 {
 149    fListHistCascade = new TList();
 150
 151    if (! fHistEventMultiplicity) {
 152       fHistEventMultiplicity   = new TH1F("fHistEventMultiplicity" , "Nb of Events" , 4, -1.0, 3.0);
 153       fListHistCascade->Add(fHistEventMultiplicity);
 154    }
 155
 156    if (! fHistEventMultiplicityRAVS) {
 157       fHistEventMultiplicityRAVS   = new TH1F("fHistEventMultiplicityRAVS" , "Nb of Events Rejected After Vertex selection" , 4, -1.0, 3.0);
 158       fListHistCascade->Add(fHistEventMultiplicityRAVS);
 159    }
 160
 161    if (! fNtuple1) {
 162       fNtuple1 = new TNtuple("fNtuple1", "Ntuple1", "TrkNmb");
 163       fNtuple1->SetDirectory(0);
 164       fListHistCascade->Add(fNtuple1);
 165    }
 166
 167    if (! fNtuple2) {
 168       fNtuple2 = new TNtuple("fNtuple2", "Ntuple2", "s:dcal:lCosPoinAn:lDaugDCA:lambdap:lambdapt:lambdamass");
 169       fNtuple2->SetDirectory(0);
 170       fListHistCascade->Add(fNtuple2);
 171    }
 172
 173    if (! fNtuple3) {
 174       fNtuple3 = new TNtuple("fNtuple3", "Ntuple3", "c:dcapi:ppi:ptpi:bachphi:bachtheta:okPiTPC:okPiTOF");
 175       fNtuple3->SetDirectory(0);
 176       fListHistCascade->Add(fNtuple3);
 177    }
 178
 179    if (! fNtuple4) {
 180       fNtuple4 = new TNtuple("fNtuple4", "Ntuple4", "dca:mc:phi:theta:eta:y:pt:p:opang:invmass");
 181       fListHistCascade->Add(fNtuple4);
 182    }
 183
 184
 185 }// end UserCreateOutputObjects
 186
 187
 188 //________________________________________________________________________
 189 void AliAnalysisTaskSigma1385::UserExec(Option_t *)
 190 {
 191
 192    // Main loop
 193    // Called for each event
 194
 195
 196    Info("AliAnalysisTaskSigma1385", "Starting UserExec");
 197
 198    AliMCEventHandler* eventHandler;
 199    AliMCEvent* mcEvent = 0;
 200
 201    if (fDataType == "SIM") {
 202
 203       eventHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 204       if (!eventHandler) {
 205          Printf("ERROR: Could not retrieve MC event handler");
 206          return;
 207       }
 208
 209       mcEvent = eventHandler->MCEvent();
 210       if (!mcEvent) {
 211          Printf("ERROR: Could not retrieve MC event");
 212          return;
 213       }
 214
 215    }
 216
 217    AliStack* stack = 0;
 218    if (fDataType == "SIM") {stack = mcEvent->Stack(); fIsMC = 1; fisMC = 1;}
 219
 220    AliESDEvent *lESDevent = 0x0;
 221    AliAODEvent *lAODevent = 0x0;
 222
 223
 224    // Connect to the InputEvent
 225    Int_t ncascades = -1;
 226    if (fAnalysisType == "ESD") {
 227       lESDevent = dynamic_cast<AliESDEvent*>(InputEvent());
 228       if (lESDevent)
 229          ncascades = lESDevent->GetNumberOfCascades();
 230       else {
 231          Printf("ERROR: lESDevent not available \n");
 232          return;
 233       }
 234    } else if (fAnalysisType == "AOD") {
 235       lAODevent = dynamic_cast<AliAODEvent*>(InputEvent());
 236       if (lAODevent)
 237          ncascades = lAODevent->GetNumberOfCascades();
 238       else {
 239          Printf("ERROR: lAODevent not available \n");
 240          return;
 241       }
 242    } else {
 243       Printf("Input type undefined \n");
 244       return;
 245    }
 246
 247
 248    //------------------------------for PID
 249
 250    SetCheckITS(kTRUE);
 251    SetCheckTPC(kTRUE);
 252    SetCheckTOF(kTRUE);
 253
 254 // ----> set TPC range for PID and calibration
 255    SetTPCrange(5.0, 3.0);
 256    SetTPCpLimit(0.35);
 257
 258    // ----> set ITS range for PID
 259    SetITSband(4.0);
 260
 261    // ----> set TPC calibration
 262    if (fDataType == "SIM") SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
 263    else       SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
 264
 265    // ----> set the TOF calibration depending on type of input (sim/data)
 266    SetTOFcorrectTExp(kTRUE);
 267    SetTOFuseT0(kTRUE);
 268    SetTOFresolution(100.0);
 269    if (fDataType == "SIM") {
 270       SetTOFcalibrateESD(kFALSE);
 271       SetTOFtuneMC(kTRUE);
 272    } else {
 273
 274       SetTOFcalibrateESD(kTRUE);
 275       SetTOFtuneMC(kFALSE);
 276    }
 277
 278
 279    if (!fESDpid) {
 280       fESDpid = new AliESDpid;
 281       fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0], fTPCpar[1], fTPCpar[2], fTPCpar[3], fTPCpar[4]);
 282    }
 283
 284    // initialize DB to current run
 285    if (!lESDevent) {
 286       Printf("ERROR: lESDevent not available \n");
 287       return;
 288    }
 289    Int_t run = lESDevent->GetRunNumber();
 290    if (run != fLastRun) {
 291       //cout << "Run = " << run << " -- LAST = " << fLastRun << endl;
 292       fLastRun = run;
 293
 294       // setup TOF maker & calibration
 295       if (!fTOFcalib) fTOFcalib = new AliTOFcalib;
 296       fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 297       if (!fTOFmaker) fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 298       fTOFmaker->SetTimeResolution(fTOFresolution);
 299
 300       AliCDBManager *cdb = AliCDBManager::Instance();
 301       cdb->ClearCache(); // suggestion by Annalisa
 302       cdb->Clear();      // suggestion by Annalisa
 303       cdb->SetDefaultStorage("raw://");
 304       cdb->SetRun(run);
 305       fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 306       fTOFcalib->Init();
 307    }
 308
 309    // if required, calibrate the TOF t0 maker with current event
 310    if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(lESDevent);
 311    if (fTOFtuneMC) fTOFmaker->TuneForMC(lESDevent);
 312    if (fTOFuseT0) {
 313       fTOFmaker->ComputeT0TOF(lESDevent);
 314       fTOFmaker->ApplyT0TOF(lESDevent);
 315       fESDpid->MakePID(lESDevent, kFALSE, 0.);
 316    }
 317
 318
 319    //--------------------------------------------------------
 320
 321
 322    fHistEventMultiplicity->Fill(1);
 323
 324    //Some Quantities to characterize the event
 325
 326    Double_t b = lESDevent->GetMagneticField();
 327    Int_t trackNumber = lESDevent->GetNumberOfTracks();
 328
 329
 330    //---------------------------
 331    // new part from AliCascadeVertexer (use vertexer for reconstructing Sigma(1385)
 332    //
 333    //const AliESDVertex *vtxT3D=lESDevent->GetPrimaryVertex();
 334    const AliESDVertex *vtxT3D = lESDevent->GetPrimaryVertexTracks();
 335    if (vtxT3D->GetNContributors() < 1) {
 336       // SPD vertex
 337       vtxT3D = lESDevent->GetPrimaryVertexSPD();
 338       if (vtxT3D->GetNContributors() < 1) {
 339
 340
 341          fHistEventMultiplicityRAVS->Fill(1);
 342          return;
 343       }
 344
 345
 346    }
 347
 348    Double_t xPrimaryVertex = vtxT3D->GetX();
 349    Double_t yPrimaryVertex = vtxT3D->GetY();
 350    Double_t zPrimaryVertex = vtxT3D->GetZ();
 351
 352    if (zPrimaryVertex > 10 || zPrimaryVertex < -10)  return;
 353
 354    //-------------------------------------------------------
 355    // New Part about tracks global selection criteria
 356    //-------------------------------------------------------
 357
 358    AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
 359
 360    esdTrackCuts->SetAcceptKinkDaughters(0); // 0 = kFalse
 361    esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 362    esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 363    esdTrackCuts->SetMinNClustersTPC(70);
 364
 365
 366    //-------------------------------------------------------
 367    // loops over V0s
 368    //stores relevant V0s in an array
 369    Int_t nV0 = (Int_t)lESDevent->GetNumberOfV0s();
 370
 371    TObjArray vtcs(nV0);
 372    for (Int_t i = 0; i < nV0; i++) {
 373       AliESDv0 *v = lESDevent->GetV0(i);
 374       if (v->GetOnFlyStatus()) continue; // if kTRUE, then this V0 is recontructed
 375
 376       vtcs.AddLast(v);
 377    }
 378    nV0 = vtcs.GetEntriesFast();
 379
 380    //-------------------------------------------------------
 381    // loops over bachelor tracks
 382    // stores relevant tracks in another array
 383
 384    Int_t nentr = (Int_t)lESDevent->GetNumberOfTracks();
 385    TArrayI trk(nentr); Int_t ntr = 0;
 386
 387    for (Int_t i = 0; i < nentr; i++) {
 388       AliESDtrack *esdtr = lESDevent->GetTrack(i);
 389
 390       if (!(esdtr->GetStatus() & AliESDtrack::kITSrefit)) continue;
 391       if (!(esdtr->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 392
 393
 394       if (!(esdTrackCuts->AcceptTrack(esdtr))) continue;
 395
 396       trk[ntr++] = i;
 397    }
 398
 399
 400
 401    //-----------------------------------------------------
 402    // nested loops over V0s and bachelors
 403    Float_t massLambda = 1.11568;
 404    Int_t ncasc = 0;
 405    AliCascadeVertexer *cascvert = new AliCascadeVertexer();
 406    for (Int_t i = 0; i < nV0; i++) { //loop on V0s
 407
 408       // Lambda
 409       AliESDv0 *v = (AliESDv0*)vtcs.UncheckedAt(i);
 410       Int_t strness = 0;
 411       Double_t lambdaMass = 0;
 412       Float_t lambdaP = 0;
 413       Float_t lambdaPt = 0;
 414       Float_t lambdaDCA = 0; // DCA between Lambda and Primary Vertex
 415       Double_t v0cospointangle = 0;
 416       Float_t v0daughtersDCA = 0;
 417
 418
 419       // Lambda quality cuts
 420       UInt_t lIdxPosXi    = (UInt_t) TMath::Abs(v->GetPindex());
 421       UInt_t lIdxNegXi    = (UInt_t) TMath::Abs(v->GetNindex());
 422
 423       AliESDtrack *pTrackXi     = lESDevent->GetTrack(lIdxPosXi);
 424       AliESDtrack *nTrackXi     = lESDevent->GetTrack(lIdxNegXi);
 425
 426       // Filter like-sign V0
 427       if ( (TMath::Abs(pTrackXi->GetSign()) - TMath::Abs(nTrackXi->GetSign()) ) < 0.1) continue;
 428
 429       // WARNING: the following selections cannot be done for AOD yet...
 430
 431
 432       v->ChangeMassHypothesis(kLambda0); // the v0 must be Lambda
 433
 434
 435       if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
 436
 437          if (!(pTrackXi->GetStatus() & AliESDtrack::kTPCrefit))  continue;
 438          if (!(nTrackXi->GetStatus() & AliESDtrack::kTPCrefit))  continue;
 439
 440
 441          if ((pTrackXi->GetTPCNcls()) < 70) continue;
 442          if ((nTrackXi->GetTPCNcls()) < 70) continue;
 443
 444          strness = 1; lambdaMass = v->GetEffMass(); lambdaP = v->P(); lambdaPt = v->Pt(); lambdaDCA = v->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex); v0cospointangle = v->GetV0CosineOfPointingAngle(); v0daughtersDCA = v->GetDcaV0Daughters();
 445
 446
 447
 448       }
 449
 450
 451       v->ChangeMassHypothesis(kLambda0Bar); // the v0 must be Anti Lambda
 452
 453
 454
 455       if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
 456
 457          if (!(pTrackXi->GetStatus() & AliESDtrack::kTPCrefit))  continue;
 458          if (!(nTrackXi->GetStatus() & AliESDtrack::kTPCrefit))  continue;
 459
 460
 461          if ((pTrackXi->GetTPCNcls()) < 70) continue;
 462          if ((nTrackXi->GetTPCNcls()) < 70) continue;
 463
 464          Int_t temp = strness + 1; strness = -1 * temp; lambdaMass = v->GetEffMass(); lambdaP = v->P(); lambdaPt = v->Pt(); lambdaDCA = v->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex); v0cospointangle = v->GetV0CosineOfPointingAngle(); v0daughtersDCA = v->GetDcaV0Daughters();
 465
 466
 467       }
 468
 469       if (strness == 0) continue;
 470
 471
 472       for (Int_t j = 0; j < ntr; j++) { //loop on tracks
 473
 474          // Pion bachelor
 475          Int_t bidx = trk[j];
 476          Int_t bachTPCcls = 0;
 477
 478
 479          if ((bidx == v->GetIndex(0)) || (bidx == v->GetIndex(1))) continue; // bach and V0's daughter's must be different!
 480
 481          AliESDtrack *btrk = lESDevent->GetTrack(bidx);
 482
 483          if (!(btrk->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 484
 485          if ((bachTPCcls = btrk->GetTPCNcls()) < 70) continue;
 486
 487          //Bool_t *IsOkTrack = IsSelected(btrk); //for Alberto's PID
 488
 489          //Bool_t *okTrack = new Bool_t[5];
 490
 491          //for (Int_t k = 0; k < 5; k++) okTrack[k] = IsOkTrack[k];
 492
 493          IsSelected(btrk);
 494
 495
 496          Int_t bachCharge = btrk->Charge();
 497          Float_t pionDCA = TMath::Abs(btrk->GetD(xPrimaryVertex, yPrimaryVertex, b)); // DCA between Bachelor and Primary Vertex
 498
 499          Double_t bachphi = btrk->Phi();
 500          Double_t bachtheta = btrk->Theta();
 501          Double_t bachmass = 0.13957;
 502          Double_t bachp = btrk->GetP();
 503          Double_t bachpt = btrk->Pt();
 504
 505
 506          // Distance between Lambda and Pion bachelor
 507          AliESDv0 v0(*v), *pv0 = &v0;
 508          AliExternalTrackParam bt(*btrk), *pbt = &bt;
 509          Double_t dca = cascvert->PropagateToDCA(pv0, pbt, b); // distance bwn V0 and bach (in cm)
 510          if (dca > 10.0) continue;  // Note: was 0.1! Enlarged->further filter in second pass analysis
 511
 512
 513          AliESDcascade cascade(*pv0, *pbt, bidx); //constucts a sigma1385 candidate
 514          AliESDcascade *xi = &cascade;
 515
 516
 517          UInt_t lBachIdx   = (UInt_t) TMath::Abs(xi->GetBindex());
 518          AliESDtrack *bachTrackXi   = lESDevent->GetTrack(lBachIdx);
 519
 520
 521
 522          Short_t mcTrue=0;
 523          if (fDataType == "SIM") {
 524
 525             Int_t pLabel    = TMath::Abs(pTrackXi->GetLabel());
 526             Int_t nLabel    = TMath::Abs(nTrackXi->GetLabel());
 527             Int_t bachLabel = TMath::Abs(bachTrackXi->GetLabel());
 528
 529             Int_t motherpLabel = stack->Particle(pLabel)->GetFirstMother();
 530             Int_t mothernLabel = stack->Particle(nLabel)->GetFirstMother();
 531             Int_t motherbachLabel = stack->Particle(bachLabel)->GetFirstMother();
 532
 533
 534             //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 535             //cout<< "plabel " << pLabel << " nlabel " << nLabel << " mother p " << motherpLabel << " mother n " << mothernLabel << " bachlabel " << bachLabel << " mother bach " << motherbachLabel << endl;
 536             //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 537
 538             if (motherpLabel > -1 && mothernLabel > -1) {
 539                TParticle *plambda = stack->Particle(motherpLabel);
 540                Int_t grandmother = plambda->GetFirstMother();
 541
 542                motherpLabel = TMath::Abs(motherpLabel);
 543                mothernLabel = TMath::Abs(mothernLabel);
 544                //motherbachLabel = TMath::Abs(motherbachLabel);
 545
 546                //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 547                //cout<< "plabel " << pLabel << " nlabel " << nLabel << " mother p " << motherpLabel << " mother n " << mothernLabel << " mother lambda " << grandmother << " bachlabel " << bachLabel << " mother bach " << motherbachLabel << endl;
 548                //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 549
 550                if (motherbachLabel > -1) {
 551
 552
 553                   //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 554                   //cout<< "mother lambda " << grandmother << " mother bach " << motherbachLabel << " PDG Code "<< stack->Particle(grandmother)->GetPdgCode() << endl;
 555                   //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 556
 557                   if ((motherpLabel == mothernLabel) && (grandmother == motherbachLabel) &&
 558                       ((TMath::Abs(stack->Particle(grandmother)->GetPdgCode()) == 3114) ||
 559                        (TMath::Abs(stack->Particle(grandmother)->GetPdgCode()) == 3224))) mcTrue = 1;
 560
 561                   if( (motherpLabel == mothernLabel) && (grandmother == motherbachLabel) &&
 562                          (TMath::Abs(stack->Particle(grandmother)->GetPdgCode())==3312) ) mcTrue = 2;
 563
 564                }
 565
 566             }
 567
 568          }
 569
 570          Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
 571          Double_t lPMom[3] = {0., 0., 0.,};
 572          Double_t lNMom[3] = {0., 0., 0.,};
 573          Float_t lLambdaMomX       = 0., lLambdaMomY  = 0., lLambdaMomZ   = 0.;
 574          xi->GetBPxPyPz(lBachMomX,  lBachMomY,  lBachMomZ);
 575          xi->GetPPxPyPz(lPMom[0],  lPMom[1],  lPMom[2]);
 576          xi->GetNPxPyPz(lNMom[0],  lNMom[1],  lNMom[2]);
 577          lLambdaMomX = lPMom[0] + lNMom[0];
 578          lLambdaMomY = lPMom[1] + lNMom[1];
 579          lLambdaMomZ = lPMom[2] + lNMom[2];
 580
 581          Float_t  lRapXi    = xi->RapXi();
 582          Float_t  lEta      = xi->Eta();
 583          Float_t  lTheta    = xi->Theta();
 584          Float_t  lPhi      = xi->Phi();
 585          Float_t  lPt       = xi->Pt();
 586          Float_t   lP     = xi->P();
 587
 588          // Support variables for invariant mass calculation
 589          TLorentzVector lambda, pion, sigma;
 590          Double_t eLambda = TMath::Sqrt(lLambdaMomX * lLambdaMomX + lLambdaMomY * lLambdaMomY + lLambdaMomZ * lLambdaMomZ + lambdaMass * lambdaMass);
 591          Double_t ePion = TMath::Sqrt(lBachMomX * lBachMomX + lBachMomY * lBachMomY + lBachMomZ * lBachMomZ + bachmass * bachmass) ;
 592
 593          lambda.SetPxPyPzE(lLambdaMomX, lLambdaMomY, lLambdaMomZ, eLambda);
 594          pion.SetPxPyPzE(lBachMomX, lBachMomY, lBachMomZ, ePion);
 595
 596          sigma = lambda + pion;
 597
 598          Double_t openingangle, invmass = 0;
 599
 600          openingangle = TMath::ACos((lBachMomX * lLambdaMomX + lBachMomY * lLambdaMomY + lBachMomZ * lLambdaMomZ) / ((TMath::Sqrt(lBachMomX * lBachMomX + lBachMomY * lBachMomY + lBachMomZ * lBachMomZ)) * (TMath::Sqrt(lLambdaMomX * lLambdaMomX + lLambdaMomY * lLambdaMomY +
 601                                     lLambdaMomZ * lLambdaMomZ))));
 602
 603          invmass = sigma.M();
 604
 605
 606          fNtuple1->Fill(trackNumber);
 607          fNtuple2->Fill((1.*strness), lambdaDCA, v0cospointangle, v0daughtersDCA , lambdaP, lambdaPt, lambdaMass);
 608          //fNtuple3->Fill((1.*bachCharge), pionDCA, bachp, bachpt, bachphi, bachtheta, okTrack[1], okTrack[2]);
 609          fNtuple3->Fill((1.*bachCharge), pionDCA, bachp, bachpt, bachphi, bachtheta, fOkTrack[1], fOkTrack[2]);
 610          fNtuple4->Fill(dca, (1.*mcTrue), lPhi, lTheta, lEta, lRapXi, lPt, lP, openingangle, invmass);
 611
 612          ncasc++;
 613
 614          //delete IsOkTrack;
 615          //delete okTrack;
 616
 617
 618       } // end loop tracks
 619    } // end loop V0s
 620
 621    Info("AliAnalysisTaskSigma1385", "Number of reconstructed Sigma(1385): %d", ncasc);
 622
 623
 624    // Post output data.
 625    PostData(1, fListHistCascade);
 626
 627 }
 628
 629 //________________________________________________________________________
 630
 631 //Bool_t *AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
 632 void AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
 633 {
 634 //
 635 //
 636    //Bool_t okTrack[5];
 637
 638    //for (Int_t i = 0; i < 5; i++) okTrack[i] = kFALSE;
 639    for (Int_t i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
 640
 641
 642    AliITSPIDResponse itsrsp(fIsMC);
 643
 644
 645    ULong_t  status;
 646    Int_t    k, nITS;
 647    Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofRel;
 648    Bool_t   okTOF = kFALSE;
 649    Bool_t   okTPC = kFALSE;
 650    Bool_t   okITS = kFALSE;
 651    Bool_t   isTPC = kFALSE;
 652    Bool_t   isITSSA = kFALSE;
 653    Bool_t   isTOF = kFALSE;
 654    UChar_t  itsCluMap;
 655
 656    // get commonly used variables
 657    status  = (ULong_t)track->GetStatus();
 658    mom     = track->P();
 659    isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 660    isITSSA = ((!isTPC) && (track->IsOn(AliESDtrack::kITSrefit)) && (!track->IsOn(AliESDtrack::kITSpureSA))); //(!isTPC && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 661    isTOF   = (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) /* && mom > TMath::Max(b1, b2)*/);
 662
 663
 664    // check if the track type matches what is required
 665    if (!isTPC && !isITSSA) {
 666       AliDebug(AliLog::kDebug + 2, "Track is not either a TPC track or a ITS standalone. Rejected");
 667       //return okTrack;
 668       return;
 669
 670    } else if (isTPC && !fUseGlobal) {
 671       AliDebug(AliLog::kDebug + 2, "Global tracks not used. Rejected");
 672       //return okTrack;
 673       return;
 674
 675    } else if (isITSSA && !fUseITSSA) {
 676       AliDebug(AliLog::kDebug + 2, "ITS standalone not used. Rejected");
 677       //return okTrack;
 678       return;
 679    }
 680
 681    // does a preliminary check on TOF values, if necessary
 682    // then, if the reference time or TOF signal are meaningless
 683    // even if the 'isTOF' flag is true, switch it to false
 684    if (isTOF) {
 685       track->GetIntegratedTimes(times);
 686       tofTime  = (Double_t)track->GetTOFsignal();
 687       tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kPion], AliPID::ParticleMass(AliPID::kPion));
 688       tofRef   = times[AliPID::kPion];
 689       if (tofRef <= 0.0 && tofSigma <= 0.0) isTOF = kFALSE;
 690    }
 691
 692
 693
 694    // check TPC dE/dx
 695    if (isTPC) { // this branch is entered by all global tracks
 696       // check TPC dE/dx:
 697       if (fCheckTPC) {
 698          tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kPion));
 699          if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
 700          okTPC = (tpcNSigma <= tpcMaxNSigma);
 701          AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", tpcNSigma, tpcMaxNSigma, (okTPC ? "passed" : "failed")));
 702       } else {
 703          // if TPC is not checked, it is as if all tracks do pass the cut
 704          okTPC = kTRUE;
 705          AliDebug(AliLog::kDebug + 2, "TPC not checked, track accepted");
 706       }
 707
 708       // check TOF (only if flags are OK)
 709       if (fCheckTOF) {
 710          if (isTOF) {
 711             // TOF can be checked only when track is matched there
 712             track->GetIntegratedTimes(times);
 713             tofTime  = (Double_t)track->GetTOFsignal();
 714             tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kPion], AliPID::ParticleMass(AliPID::kPion));
 715             tofRef   = times[AliPID::kPion];
 716
 717             tofRel   = (tofTime - tofRef) / tofSigma;
 718             okTOF    = ((tofRel >= fMinTOF) && (tofRel <= fMaxTOF));
 719             AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- range = %f %f -- cut %s", tofRel, fMinTOF, fMaxTOF, (okTOF ? "passed" : "failed")));
 720          } else {
 721             // if TOF is not matched, the answer depends on TPC:
 722             // - if TPC is required, track is checked only there and TOF simply ignored
 723             // - if TPC is not required, track is rejected when TOF does not match it, if TOF check is required
 724             if (fCheckTPC) okTOF = kTRUE; else okTOF = kFALSE;
 725          }
 726       } else {
 727          okTOF = kTRUE;
 728       }
 729
 730    } else if (isITSSA) { // this branch is entered by all ITS standalone tracks
 731       // check dE/dx only if this is required, otherwise ITS standalone are just added but not checked for PID
 732       if (fCheckITS) {
 733          itsSignal = track->GetITSsignal();
 734          itsCluMap = track->GetITSClusterMap();
 735          nITS      = 0;
 736          for (k = 2; k < 6; k++) if (itsCluMap & (1 << k)) ++nITS;
 737          //if (nITS < 3) return kFALSE;
 738          //itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kPion, nITS, kTRUE);
 739          //okITS = ((TMath::Abs(itsNSigma)) <= fMaxITSband);
 740          if (nITS < 3)
 741             okITS = kFALSE;
 742          else {
 743             itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kPion, nITS, kTRUE);
 744             okITS = ((TMath::Abs(itsNSigma)) <= fMaxITSband);
 745             AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", itsNSigma, fMaxITSband, (okITS ? "passed" : "failed")));
 746          }
 747       } else {
 748          okITS = kTRUE;
 749       }
 750    } else {
 751       // if we are here, the track is surely bad
 752       okITS = kFALSE;
 753       okTPC = kFALSE;
 754       okTOF = kFALSE;
 755    }
 756
 757
 758    //okTrack[0] = okITS;
 759    //okTrack[1] = okTPC;
 760    //okTrack[2] = okTOF;
 761    //okTrack[3] = isTPC;
 762    //okTrack[4] = isITSSA;
 763
 764    fOkTrack[0] = okITS;
 765    fOkTrack[1] = okTPC;
 766    fOkTrack[2] = okTOF;
 767    fOkTrack[3] = isTPC;
 768    fOkTrack[4] = isITSSA;
 769
 770    //cout<<"##########################################"<<endl;
 771    //cout<<"isITSSA "<<isITSSA<< " isTPC "<<isTPC<<endl;
 772    //cout<<"##########################################"<<endl;
 773
 774
 775    //return okTrack;
 776 }
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788 //________________________________________________________________________
 789 void AliAnalysisTaskSigma1385::Terminate(Option_t *)
 790 {
 791    // Draw result to the screen
 792    // Called once at the end of the query
 793 }
 794
 795
 796