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