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