PWGLF/QATasks/AliAnalysisTaskQAMultistrange.cxx

   1 /**************************************************************************
   2  *  Authors : Domenico Colella                                            *
   3  *                                                                        *
   4  *                                                                        *
   5  * Derived from the:                                                      *
   6  *  - Original AliAnalysisTaskCheckCascade (A. Maire, G. Hippolyte)       *
   7  *  - Adapted to PbPb analysis (M. Nicassio)                              *
   8  *  - Adapted to work on all collisidng systems (pp, PbPb and pPb),       *
   9  *    ESD/AOD and experimental/MC data (D. Colella)                       *
  10  *                                                                        *
  11  * Permission to use, copy, modify and distribute this software and its   *
  12  * documentation strictly for non-commercial purposes is hereby granted   *
  13  * without fee, provided that the above copyright notice appears in all   *
  14  * copies and that both the copyright notice and this permission notice   *
  15  * appear in the supporting documentation. The authors make no claims     *
  16  * about the suitability of this software for any purpose. It is          *
  17  * provided "as is" without express or implied warranty.                  *
  18  **************************************************************************/
  19
  20 class TTree;
  21 class TParticle;
  22 class TVector3;
  23
  24 class AliESDVertex;
  25 class AliAODVertex;
  26 class AliESDv0;
  27 class AliAODv0;
  28
  29
  30 #include <Riostream.h>
  31 #include "THnSparse.h"
  32 #include "TVector3.h"
  33 #include "TMath.h"
  34
  35 #include "AliLog.h"
  36 #include "AliCentrality.h"
  37 #include "AliHeader.h"   //for MC
  38 #include "AliMCEvent.h"  //for MC
  39 #include "AliStack.h"    //for MC
  40 #include "AliESDEvent.h"
  41 #include "AliAODEvent.h"
  42 #include "AliESDtrackCuts.h"
  43 #include "AliPIDResponse.h"
  44
  45 #include "AliInputEventHandler.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliESDInputHandler.h"
  48 #include "AliAODInputHandler.h"
  49 #include "AliCFContainer.h"
  50 #include "AliMultiplicity.h"
  51
  52 #include "AliGenEventHeader.h"         //for MC
  53 #include "AliGenCocktailEventHeader.h" //for MC
  54 #include "AliGenHijingEventHeader.h"   //for MC
  55 #include "AliAODMCParticle.h"          //for MC
  56
  57 #include "AliESDcascade.h"
  58 #include "AliAODcascade.h"
  59 #include "AliESDtrack.h"
  60 #include "AliAODTrack.h"
  61 #include "AliESDpid.h"
  62
  63
  64 #include "AliAnalysisTaskQAMultistrange.h"
  65
  66 ClassImp(AliAnalysisTaskQAMultistrange)
  67
  68
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskQAMultistrange::AliAnalysisTaskQAMultistrange()
  72   : AliAnalysisTaskSE           (),
  73     fisMC                       (kFALSE),
  74     fAnalysisType               ("ESD"),
  75     fCollidingSystem            ("PbPb"),
  76     fPIDResponse                (0),
  77     fkSDDSelectionOn            (kTRUE),
  78     fkQualityCutZprimVtxPos     (kTRUE),
  79     fkQualityCutNoTPConlyPrimVtx(kTRUE),
  80     fkQualityCutTPCrefit        (kTRUE),
  81     fkQualityCutnTPCcls         (kTRUE),
  82     fkQualityCutPileup          (kTRUE),
  83     fwithSDD                    (kTRUE),
  84     fMinnTPCcls                 (0),
  85     fCentrLowLim                (0),
  86     fCentrUpLim                 (0),
  87     fCentrEstimator             (0),
  88     fkUseCleaning               (0),
  89     fVtxRange                   (0),
  90     fMinPtCutOnDaughterTracks   (0),
  91     fEtaCutOnDaughterTracks     (0),
  92
  93     fCFContCascadeCuts(0),
  94     fCFContCascadeMCgen(0)
  95
  96 {
  97   // Dummy Constructor
  98 }
  99
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskQAMultistrange::AliAnalysisTaskQAMultistrange(const char *name)
 103   : AliAnalysisTaskSE           (name),
 104     fisMC                       (kFALSE),
 105     fAnalysisType               ("ESD"),
 106     fCollidingSystem            ("PbPb"),
 107     fPIDResponse                (0),
 108     fkSDDSelectionOn            (kTRUE),
 109     fkQualityCutZprimVtxPos     (kTRUE),
 110     fkQualityCutNoTPConlyPrimVtx(kTRUE),
 111     fkQualityCutTPCrefit        (kTRUE),
 112     fkQualityCutnTPCcls         (kTRUE),
 113     fkQualityCutPileup          (kTRUE),
 114     fwithSDD                    (kTRUE),
 115     fMinnTPCcls                 (0),
 116     fCentrLowLim                (0),
 117     fCentrUpLim                 (0),
 118     fCentrEstimator             (0),
 119     fkUseCleaning               (0),
 120     fVtxRange                   (0),
 121     fMinPtCutOnDaughterTracks   (0),
 122     fEtaCutOnDaughterTracks     (0),
 123
 124     fCFContCascadeCuts(0),
 125     fCFContCascadeMCgen(0)
 126
 127 {
 128   // Constructor
 129   // Output slot #0 writes into a TList container (Cascade)
 130   DefineOutput(1, AliCFContainer::Class());
 131   DefineOutput(2, AliCFContainer::Class());
 132
 133   AliLog::SetClassDebugLevel("AliAnalysisTaskQAMultistrange",1);
 134 }
 135
 136
 137 AliAnalysisTaskQAMultistrange::~AliAnalysisTaskQAMultistrange() {
 138   //
 139   // Destructor
 140   //
 141   // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
 142   // They will be deleted when fListCascade is deleted by the TSelector dtor
 143   // Because of TList::SetOwner() ...
 144   if (fCFContCascadeCuts && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())  { delete fCFContCascadeCuts;     fCFContCascadeCuts = 0x0;  }
 145   if (fCFContCascadeMCgen && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fCFContCascadeMCgen;     fCFContCascadeMCgen = 0x0;  }
 146 }
 147
 148
 149
 150 //_____________________________________________________________
 151 void AliAnalysisTaskQAMultistrange::UserCreateOutputObjects() {
 152   // Create histograms
 153   // Called once
 154
 155   //----------------------------------------------------------------------------
 156   // Option for AliLog: to suppress the extensive info prompted by a run with MC
 157   //----------------------------------------------------------------------------
 158   if (fisMC) AliLog::SetGlobalLogLevel(AliLog::kError);
 159
 160   //-----------------------------------------------
 161   // Particle Identification Setup (new PID object)
 162   //-----------------------------------------------
 163   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 164   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 165   fPIDResponse = inputHandler->GetPIDResponse();
 166
 167   //---------------------------------------------------
 168   // Define the container for the topological variables
 169   //---------------------------------------------------
 170   if(! fCFContCascadeCuts) {
 171       // Container meant to store all the relevant distributions corresponding to the cut variables.
 172       // NB: overflow/underflow of variables on which we want to cut later should be 0!!!
 173       const Int_t  lNbSteps      =  4 ;
 174       const Int_t  lNbVariables  =  20;
 175       //Array for the number of bins in each dimension :
 176       Int_t lNbBinsPerVar[lNbVariables] = {0};
 177       lNbBinsPerVar[0]  = 25;     //DcaCascDaughters             :  [0.0,2.4,3.0]       -> Rec.Cut = 2.0;
 178       lNbBinsPerVar[1]  = 25;     //DcaBachToPrimVertex          :  [0.0,0.24,100.0]    -> Rec.Cut = 0.01;
 179       lNbBinsPerVar[2]  = 30;     //CascCosineOfPointingAngle    :  [0.97,1.0]          -> Rec.Cut = 0.98;
 180       lNbBinsPerVar[3]  = 40;     //CascRadius                   :  [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
 181       lNbBinsPerVar[4]  = 30;     //InvMassLambdaAsCascDghter    :  [1.1,1.3]           -> Rec.Cut = 0.008;
 182       lNbBinsPerVar[5]  = 20;     //DcaV0Daughters               :  [0.0,2.0]           -> Rec.Cut = 1.5;
 183       lNbBinsPerVar[6]  = 201;    //V0CosineOfPointingAngleToPV  :  [0.89,1.0]          -> Rec.Cut = 0.9;
 184       lNbBinsPerVar[7]  = 40;     //V0Radius                     :  [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
 185       lNbBinsPerVar[8]  = 40;     //DcaV0ToPrimVertex            :  [0.0,0.39,110.0]    -> Rec.Cut = 0.01;
 186       lNbBinsPerVar[9]  = 25;     //DcaPosToPrimVertex           :  [0.0,0.24,100.0]    -> Rec.Cut = 0.05;
 187       lNbBinsPerVar[10] = 25;     //DcaNegToPrimVertex           :  [0.0,0.24,100.0]    -> Rec.Cut = 0.05
 188       lNbBinsPerVar[11] = 150;    //InvMassXi                    :   2-MeV/c2 bins
 189       lNbBinsPerVar[12] = 120;    //InvMassOmega                 :   2-MeV/c2 bins
 190       lNbBinsPerVar[13] = 100;    //XiTransvMom                  :  [0.0,10.0]
 191       lNbBinsPerVar[14] = 110;    //Y(Xi)                        :   0.02 in rapidity units
 192       lNbBinsPerVar[15] = 110;    //Y(Omega)                     :   0.02 in rapidity units
 193       lNbBinsPerVar[16] = 112;    //Proper lenght of cascade
 194       lNbBinsPerVar[17] = 112;    //Proper lenght of V0
 195       lNbBinsPerVar[18] = 201;    //V0CosineOfPointingAngleToXiV
 196       lNbBinsPerVar[19] = 11;     //Centrality
 197       //define the container
 198       fCFContCascadeCuts = new AliCFContainer("fCFContCascadeCuts","Container for Cascade cuts", lNbSteps, lNbVariables, lNbBinsPerVar );
 199       //Setting the bin limits
 200       //0 -  DcaXiDaughters
 201       Double_t *lBinLim0  = new Double_t[ lNbBinsPerVar[0] + 1 ];
 202          for(Int_t i=0; i< lNbBinsPerVar[0]; i++) lBinLim0[i] = (Double_t)0.0 + (2.4 - 0.0)/(lNbBinsPerVar[0] - 1) * (Double_t)i;
 203          lBinLim0[ lNbBinsPerVar[0] ] = 3.0;
 204       fCFContCascadeCuts -> SetBinLimits(0, lBinLim0);
 205       delete [] lBinLim0;
 206       //1 - DcaToPrimVertexXi
 207       Double_t *lBinLim1  = new Double_t[ lNbBinsPerVar[1] + 1 ];
 208          for(Int_t i=0; i<lNbBinsPerVar[1]; i++) lBinLim1[i] = (Double_t)0.0 + (0.24  - 0.0)/(lNbBinsPerVar[1] - 1) * (Double_t)i;
 209          lBinLim1[ lNbBinsPerVar[1] ] = 100.0;
 210       fCFContCascadeCuts -> SetBinLimits(1, lBinLim1);
 211       delete [] lBinLim1;
 212       //2 - CascCosineOfPointingAngle
 213       fCFContCascadeCuts->SetBinLimits(2, 0.97, 1.);
 214       //3 - CascRadius
 215       Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
 216          for(Int_t i=0; i< lNbBinsPerVar[3]; i++)   lBinLim3[i]  = (Double_t)0.0   + (3.9  - 0.0 )/(lNbBinsPerVar[3] - 1)  * (Double_t)i ;
 217          lBinLim3[ lNbBinsPerVar[3] ] = 1000.0;
 218       fCFContCascadeCuts -> SetBinLimits(3,  lBinLim3 );
 219       delete [] lBinLim3;
 220       //4 - InvMassLambdaAsCascDghter
 221       fCFContCascadeCuts->SetBinLimits(4, 1.1, 1.13);
 222       //5 - DcaV0Daughters
 223       fCFContCascadeCuts -> SetBinLimits(5, 0., 2.);
 224       //6 - V0CosineOfPointingAngleToPV
 225       fCFContCascadeCuts -> SetBinLimits(6, 0.8, 1.001);
 226       //7 - V0Radius
 227       Double_t *lBinLim7 = new Double_t[ lNbBinsPerVar[7] + 1];
 228          for(Int_t i=0; i< lNbBinsPerVar[7];i++) lBinLim7[i] = (Double_t)0.0 + (3.9 - 0.0)/(lNbBinsPerVar[7] - 1) * (Double_t)i;
 229          lBinLim7[ lNbBinsPerVar[7] ] = 1000.0;
 230       fCFContCascadeCuts -> SetBinLimits(7, lBinLim7);
 231       delete [] lBinLim7;
 232       //8 - DcaV0ToPrimVertex
 233       Double_t *lBinLim8  = new Double_t[ lNbBinsPerVar[8]+1 ];
 234          for(Int_t i=0; i< lNbBinsPerVar[8];i++)   lBinLim8[i]  = (Double_t)0.0   + (0.39  - 0.0 )/(lNbBinsPerVar[8]-1)  * (Double_t)i ;
 235          lBinLim8[ lNbBinsPerVar[8]  ] = 100.0;
 236       fCFContCascadeCuts -> SetBinLimits(8,  lBinLim8 );
 237       delete [] lBinLim8;
 238       //9 - DcaPosToPrimVertex
 239       Double_t *lBinLim9  = new Double_t[ lNbBinsPerVar[9]+1 ];
 240          for(Int_t i=0; i< lNbBinsPerVar[9];i++)   lBinLim9[i]  = (Double_t)0.0   + (0.24  - 0.0 )/(lNbBinsPerVar[9]-1)  * (Double_t)i ;
 241          lBinLim9[ lNbBinsPerVar[9]  ] = 100.0;
 242       fCFContCascadeCuts -> SetBinLimits(9,  lBinLim9 );
 243       delete [] lBinLim9;
 244       //10 - DcaNegToPrimVertex
 245       Double_t *lBinLim10  = new Double_t[ lNbBinsPerVar[10]+1 ];
 246          for(Int_t i=0; i< lNbBinsPerVar[10];i++)   lBinLim10[i]  = (Double_t)0.0   + (0.24  - 0.0 )/(lNbBinsPerVar[10]-1)  * (Double_t)i ;
 247          lBinLim10[ lNbBinsPerVar[10]  ] = 100.0;
 248       fCFContCascadeCuts -> SetBinLimits(10,  lBinLim10 );
 249       delete [] lBinLim10;
 250       //11 - InvMassXi
 251       fCFContCascadeCuts->SetBinLimits(11, 1.25, 1.40);
 252       //12 - InvMassOmega
 253       fCFContCascadeCuts->SetBinLimits(12, 1.62, 1.74);
 254       //13 - XiTransvMom
 255       fCFContCascadeCuts->SetBinLimits(13, 0.0, 10.0);
 256       //14 - Y(Xi)
 257       fCFContCascadeCuts->SetBinLimits(14, -1.1, 1.1);
 258       //15 - Y(Omega)
 259       fCFContCascadeCuts->SetBinLimits(15, -1.1, 1.1);
 260       //16 - Proper time of cascade
 261       Double_t *lBinLim16  = new Double_t[ lNbBinsPerVar[16]+1 ];
 262          for(Int_t i=0; i< lNbBinsPerVar[16];i++) lBinLim16[i] = (Double_t) -1. + (110. + 1.0 ) / (lNbBinsPerVar[16] - 1) * (Double_t) i;
 263          lBinLim16[ lNbBinsPerVar[16] ] = 2000.0;
 264       fCFContCascadeCuts->SetBinLimits(16, lBinLim16);
 265       //17 - Proper time of V0
 266       fCFContCascadeCuts->SetBinLimits(17, lBinLim16);
 267       //18 - V0CosineOfPointingAngleToXiV
 268       fCFContCascadeCuts -> SetBinLimits(18, 0.8, 1.001);
 269       //19 - Centrality
 270       Double_t *lBinLim19  = new Double_t[ lNbBinsPerVar[19]+1 ];
 271          for(Int_t i=3; i< lNbBinsPerVar[19]+1;i++)   lBinLim19[i]  = (Double_t)(i-1)*10.;
 272          lBinLim19[0] = 0.0;
 273          lBinLim19[1] = 5.0;
 274          lBinLim19[2] = 10.0;
 275       fCFContCascadeCuts->SetBinLimits(19,  lBinLim19 );
 276       delete [] lBinLim19;
 277       // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
 278       fCFContCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates");
 279       fCFContCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates");
 280       fCFContCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates");
 281       fCFContCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates");
 282       // Setting the variable title, per axis
 283       fCFContCascadeCuts->SetVarTitle(0,  "Dca(cascade daughters) (cm)");
 284       fCFContCascadeCuts->SetVarTitle(1,  "ImpactParamToPV(bachelor) (cm)");
 285       fCFContCascadeCuts->SetVarTitle(2,  "cos(cascade PA)");
 286       fCFContCascadeCuts->SetVarTitle(3,  "R_{2d}(cascade decay) (cm)");
 287       fCFContCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(as casc dghter) (GeV/c^{2})");
 288       fCFContCascadeCuts->SetVarTitle(5,  "Dca(V0 daughters) in Xi (cm)");
 289       fCFContCascadeCuts->SetVarTitle(6,  "cos(V0 PA) in cascade to PV");
 290       fCFContCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
 291       fCFContCascadeCuts->SetVarTitle(8,  "ImpactParamToPV(V0) (cm)");
 292       fCFContCascadeCuts->SetVarTitle(9,  "ImpactParamToPV(Pos) (cm)");
 293       fCFContCascadeCuts->SetVarTitle(10, "ImpactParamToPV(Neg) (cm)");
 294       fCFContCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
 295       fCFContCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
 296       fCFContCascadeCuts->SetVarTitle(13, "pt(cascade) (GeV/c)");
 297       fCFContCascadeCuts->SetVarTitle(14, "Y(Xi)");
 298       fCFContCascadeCuts->SetVarTitle(15, "Y(Omega)");
 299       fCFContCascadeCuts->SetVarTitle(16, "mL/p (cascade) (cm)");
 300       fCFContCascadeCuts->SetVarTitle(17, "mL/p (V0) (cm)");
 301       fCFContCascadeCuts->SetVarTitle(18, "cos(V0 PA) in cascade to XiV");
 302       fCFContCascadeCuts->SetVarTitle(19, "Centrality");
 303   }
 304
 305
 306   //-----------------------------------------------
 307   // Define the Container for the MC generated info
 308   //-----------------------------------------------
 309   if(! fCFContCascadeMCgen) {
 310       // Container meant to store general distributions for MC generated particles.
 311       // NB: overflow/underflow of variables on which we want to cut later should be 0!!!
 312       const Int_t  lNbStepsMC      =  4 ;
 313       const Int_t  lNbVariablesMC  =  7;
 314       //Array for the number of bins in each dimension :
 315       Int_t lNbBinsPerVarMC[lNbVariablesMC] = {0};
 316       lNbBinsPerVarMC[0] = 100;    //Total momentum        : [0.0,10.0]
 317       lNbBinsPerVarMC[1] = 100;    //Transverse momentum   : [0.0,10.0]
 318       lNbBinsPerVarMC[2] = 110;    //Y                     : [-1.1,1.1]
 319       lNbBinsPerVarMC[3] = 200;    //eta                   : [-10, 10]
 320       lNbBinsPerVarMC[4] = 200;    //theta                 : [-10, 190]
 321       lNbBinsPerVarMC[5] = 360;    //Phi                   : [0., 360.]
 322       lNbBinsPerVarMC[6] = 11;     //Centrality
 323       //define the container
 324       fCFContCascadeMCgen = new AliCFContainer("fCFContCascadeMCgen","Container for MC gen cascade ", lNbStepsMC, lNbVariablesMC, lNbBinsPerVarMC );
 325       //Setting the bin limits
 326       //0 - Total Momentum
 327       fCFContCascadeMCgen->SetBinLimits(0, 0.0, 10.0);
 328       //1 - Transverse Momentum
 329       fCFContCascadeMCgen->SetBinLimits(1, 0.0, 10.0);
 330       //2 - Y
 331       fCFContCascadeMCgen->SetBinLimits(2, -1.1, 1.1);
 332       //3 - Eta
 333       fCFContCascadeMCgen->SetBinLimits(3, -10, 10);
 334       //4 - Theta
 335       fCFContCascadeMCgen->SetBinLimits(4, -10, 190);
 336       //5 - Phi
 337       fCFContCascadeMCgen->SetBinLimits(5, 0.0, 360.0);
 338       //6 - Centrality
 339       Double_t *lBinLim6MC  = new Double_t[ lNbBinsPerVarMC[6]+1 ];
 340          for(Int_t i=3; i< lNbBinsPerVarMC[6]+1;i++)   lBinLim6MC[i]  = (Double_t)(i-1)*10.;
 341          lBinLim6MC[0] = 0.0;
 342          lBinLim6MC[1] = 5.0;
 343          lBinLim6MC[2] = 10.0;
 344       fCFContCascadeMCgen->SetBinLimits(6,  lBinLim6MC);
 345       delete [] lBinLim6MC;
 346       // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
 347       fCFContCascadeMCgen->SetStepTitle(0, "#Xi^{-} candidates");
 348       fCFContCascadeMCgen->SetStepTitle(1, "#bar{#Xi}^{+} candidates");
 349       fCFContCascadeMCgen->SetStepTitle(2, "#Omega^{-} candidates");
 350       fCFContCascadeMCgen->SetStepTitle(3, "#bar{#Omega}^{+} candidates");
 351       // Setting the variable title, per axis
 352       fCFContCascadeMCgen->SetVarTitle(0,  "MC gen p_tot (GeV/c)");
 353       fCFContCascadeMCgen->SetVarTitle(1,  "MC gen p_T (GeV/c)");
 354       fCFContCascadeMCgen->SetVarTitle(2,  "MC gen Rapidity");
 355       fCFContCascadeMCgen->SetVarTitle(3,  "MC gen Pseudo-rapidity");
 356       fCFContCascadeMCgen->SetVarTitle(4,  "MC gen Theta");
 357       fCFContCascadeMCgen->SetVarTitle(5,  "MC gen Phi");
 358       fCFContCascadeMCgen->SetVarTitle(6,  "MC gen Centrality");
 359   }
 360
 361
 362   PostData(1, fCFContCascadeCuts);
 363   PostData(2, fCFContCascadeMCgen);
 364
 365 }// end UserCreateOutputObjects
 366
 367
 368 //________________________________________________________________________
 369 void AliAnalysisTaskQAMultistrange::UserExec(Option_t *) {
 370
 371
 372   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 373   // Main loop (called for each event)
 374   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 375
 376   //------------------------
 377   // Define ESD/AOD handlers
 378   //------------------------
 379   AliESDEvent  *lESDevent = 0x0;
 380   AliAODEvent  *lAODevent = 0x0;
 381   AliMCEvent   *lMCevent  = 0x0; //for MC
 382   AliStack     *lMCstack  = 0x0; //for MC
 383   TClonesArray *arrayMC   = 0;   //for MC
 384
 385   //++++++++++++++++
 386   // Starting checks
 387   //++++++++++++++++
 388
 389   //----------------------
 390   // Check the PIDresponse
 391   //----------------------
 392   if(!fPIDResponse) {
 393        AliError("Cannot get pid response");
 394        return;
 395   }
 396
 397   //---------------------
 398   // Check the InputEvent
 399   //---------------------
 400   if (fAnalysisType == "ESD") {
 401       lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 402       if (!lESDevent) {
 403           AliWarning("ERROR: lESDevent not available \n");
 404           return;
 405       }
 406       if (fisMC) {
 407           lMCevent = MCEvent();
 408           if (!lMCevent) {
 409               Printf("ERROR: Could not retrieve MC event \n");
 410               cout << "Name of the file with pb :" << CurrentFileName() << endl;
 411               return;
 412           }
 413           lMCstack = lMCevent->Stack();
 414           if (!lMCstack) {
 415               Printf("ERROR: Could not retrieve MC stack \n");
 416               cout << "Name of the file with pb :" << CurrentFileName() << endl;
 417               return;
 418           }
 419       }
 420   } else if (fAnalysisType == "AOD") {
 421       lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 422       if (!lAODevent) {
 423           AliWarning("ERROR: lAODevent not available \n");
 424           return;
 425       }
 426       if (fisMC) {
 427           arrayMC = (TClonesArray*) lAODevent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 428           if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
 429       }
 430   } else {
 431     Printf("Analysis type (ESD or AOD) not specified \n");
 432     return;
 433   }
 434
 435
 436   //+++++++++++++++++
 437   // Event Selections
 438   //+++++++++++++++++
 439
 440   //----------------------------------------------------------
 441   // Centrality/Multiplicity selection for PbPb/pPb collisions
 442   //----------------------------------------------------------
 443   // -- Centrality
 444   AliCentrality* centrality = 0;
 445   if      (fAnalysisType == "ESD" && (fCollidingSystem == "PbPb" || fCollidingSystem == "pPb")) centrality = lESDevent->GetCentrality();
 446   else if (fAnalysisType == "AOD" && (fCollidingSystem == "PbPb" || fCollidingSystem == "pPb")) centrality = lAODevent->GetCentrality();
 447   Float_t lcentrality = 0.;
 448   if (fCollidingSystem == "PbPb" || fCollidingSystem == "pPb") {
 449        if (fkUseCleaning) lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
 450        else {
 451            lcentrality = centrality->GetCentralityPercentileUnchecked(fCentrEstimator.Data());
 452            if (centrality->GetQuality()>1) {
 453                PostData(1, fCFContCascadeCuts);
 454                return;
 455            }
 456        }
 457        //if (lcentrality<fCentrLowLim||lcentrality>=fCentrUpLim) {
 458        //    PostData(1, fCFContCascadeCuts);
 459        //    return;
 460        //}
 461   } else if (fCollidingSystem == "pp") lcentrality = 0.;
 462
 463   //---------------------
 464   // SDD status selection
 465   //---------------------
 466   if (fkSDDSelectionOn && fAnalysisType == "ESD") {
 467        TString trcl = lESDevent->GetFiredTriggerClasses();
 468        if      (fwithSDD) { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
 469        else if (!fwithSDD){ if((trcl.Contains("ALLNOTRD")))  { PostData(1, fCFContCascadeCuts); return; } }
 470   } else if (fkSDDSelectionOn && fAnalysisType == "AOD") {
 471        TString trcl = lAODevent->GetFiredTriggerClasses();
 472        if      (fwithSDD)  { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
 473        else if (!fwithSDD) { if((trcl.Contains("ALLNOTRD")))  { PostData(1, fCFContCascadeCuts); return; } }
 474   }
 475
 476   //------------------
 477   // Physics selection
 478   //------------------
 479   UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 480   Bool_t isSelected = 0;
 481   if      (fCollidingSystem == "pp" || fCollidingSystem == "PbPb") isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
 482   else if (fCollidingSystem == "pPb")                              isSelected = (maskIsSelected & AliVEvent::kINT7) == AliVEvent::kINT7;
 483   if (!isSelected){ PostData(1, fCFContCascadeCuts); return; }
 484
 485   //------------------------------
 486   // Well-established PV selection
 487   //------------------------------
 488   Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 489   Double_t lMagneticField        = -10.;
 490   if (fCollidingSystem == "PbPb" || fCollidingSystem == "pp") {
 491       if (fAnalysisType == "ESD") {
 492           const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
 493           const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
 494           if (fkQualityCutNoTPConlyPrimVtx) {
 495               const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
 496               if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtx->GetStatus() ){
 497                   AliWarning(" No SPD prim. vertex nor prim. Tracking vertex ... return !");
 498                   PostData(1, fCFContCascadeCuts);
 499                   return;
 500               }
 501           }
 502           lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
 503       } else if (fAnalysisType == "AOD") {
 504           const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex();
 505           if (!lPrimaryBestAODVtx){
 506               AliWarning("No prim. vertex in AOD... return!");
 507               PostData(1, fCFContCascadeCuts);
 508               return;
 509           }
 510           lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
 511       }
 512   } else if (fCollidingSystem == "pPb") {
 513       if (fAnalysisType == "ESD") {
 514           Bool_t fHasVertex = kFALSE;
 515           const AliESDVertex *vertex = lESDevent->GetPrimaryVertexTracks();
 516           if (vertex->GetNContributors() < 1) {
 517               vertex = lESDevent->GetPrimaryVertexSPD();
 518               if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
 519               else fHasVertex = kTRUE;
 520               TString vtxTyp = vertex->GetTitle();
 521               Double_t cov[6]={0};
 522               vertex->GetCovarianceMatrix(cov);
 523               Double_t zRes = TMath::Sqrt(cov[5]);
 524               if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fHasVertex = kFALSE;
 525           }
 526           else fHasVertex = kTRUE;
 527           if (fHasVertex == kFALSE) { //Is First event in chunk rejection: Still present!
 528               AliWarning("Pb / | PV does not satisfy selection criteria!");
 529               PostData(1, fCFContCascadeCuts);
 530               return;
 531           }
 532           vertex->GetXYZ( lBestPrimaryVtxPos );
 533       } else if (fAnalysisType == "AOD") {
 534           Bool_t fHasVertex = kFALSE;
 535           const AliAODVertex *vertex = lAODevent->GetPrimaryVertex();
 536           if (vertex->GetNContributors() < 1) {
 537               vertex = lAODevent->GetPrimaryVertexSPD();
 538               if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
 539               else fHasVertex = kTRUE;
 540               TString vtxTyp = vertex->GetTitle();
 541               Double_t cov[6]={0};
 542               vertex->GetCovarianceMatrix(cov);
 543               Double_t zRes = TMath::Sqrt(cov[5]);
 544               if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fHasVertex = kFALSE;
 545           }
 546           else fHasVertex = kTRUE;
 547           if (fHasVertex == kFALSE) { //Is First event in chunk rejection: Still present!
 548               AliWarning("Pb / | PV does not satisfy selection criteria!");
 549               PostData(1, fCFContCascadeCuts);
 550               return;
 551           }
 552           vertex->GetXYZ( lBestPrimaryVtxPos );
 553       }
 554   }
 555   // -- Magnetic field
 556   if      (fAnalysisType == "ESD") lMagneticField = lESDevent->GetMagneticField();
 557   else if (fAnalysisType == "AOD") lMagneticField = lAODevent->GetMagneticField();
 558
 559   //----------------------------
 560   // Vertex Z position selection
 561   //----------------------------
 562   if (fkQualityCutZprimVtxPos) {
 563       if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
 564           PostData(1, fCFContCascadeCuts);
 565           return;
 566       }
 567   }
 568
 569   //-----------------------
 570   // Pilup selection for pp
 571   //-----------------------
 572   if (fCollidingSystem == "pp") {
 573       if (fAnalysisType == "ESD") {
 574           if (fkQualityCutPileup) { if(lESDevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
 575       } else if (fAnalysisType == "AOD") {
 576           if (fkQualityCutPileup) { if(lAODevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
 577       }
 578   }
 579
 580
 581   ////////////////////////////
 582   // MC GENERATED CASCADE PART
 583   ////////////////////////////
 584
 585   //%%%%%%%%%%%%%%%%%
 586   // Gen cascade loop
 587   if (fisMC) {
 588
 589       Int_t lNbMCPrimary = 0;
 590       if      (fAnalysisType == "ESD") lNbMCPrimary = lMCstack->GetNprimary();    //lMCstack->GetNprimary(); or lMCstack->GetNtrack();
 591       else if (fAnalysisType == "AOD") lNbMCPrimary = arrayMC->GetEntries();
 592
 593       for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) {
 594
 595            Double_t partP      = 0.;
 596            Double_t partPt     = 0.;
 597            Double_t partEta    = 0.;
 598            Double_t partTheta  = 0.;
 599            Double_t partPhi    = 0.;
 600            Double_t partRap    = 0.;
 601            Double_t partEnergy = 0.; //for Rapidity
 602            Double_t partPz     = 0.; //for Rapidity
 603            Int_t    PDGcode    = 0;
 604
 605            if ( fAnalysisType == "ESD" ) {
 606                TParticle* lCurrentParticlePrimary = 0x0;
 607                lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
 608                if (!lCurrentParticlePrimary) {
 609                    Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
 610                    continue;
 611                }
 612                if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
 613                TParticle* cascMC = 0x0;
 614                cascMC = lCurrentParticlePrimary;
 615                if (!cascMC) {
 616                    Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
 617                    continue;
 618                }
 619                partP      = cascMC->P();
 620                partPt     = cascMC->Pt();
 621                partEta    = cascMC->Eta();
 622                partTheta  = cascMC->Theta()*180.0/TMath::Pi();
 623                partPhi    = cascMC->Phi()*180.0/TMath::Pi();
 624                partEnergy = cascMC->Energy();
 625                partPz     = cascMC->Pz();
 626                PDGcode    = lCurrentParticlePrimary->GetPdgCode();
 627            } else if ( fAnalysisType == "AOD" ) {
 628                AliAODMCParticle *lCurrentParticleaod = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
 629                if (!lCurrentParticleaod) {
 630                    Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
 631                    continue;
 632                }
 633                if (!lCurrentParticleaod->IsPhysicalPrimary()) continue;
 634                partP      = lCurrentParticleaod->P();
 635                partPt     = lCurrentParticleaod->Pt();
 636                partEta    = lCurrentParticleaod->Eta();
 637                partTheta  = lCurrentParticleaod->Theta()*180.0/TMath::Pi();
 638                partPhi    = lCurrentParticleaod->Phi()*180.0/TMath::Pi();
 639                partEnergy = lCurrentParticleaod->E();
 640                partPz     = lCurrentParticleaod->Pz();
 641                PDGcode    = lCurrentParticleaod->GetPdgCode();
 642            }
 643            partRap = 0.5*TMath::Log((partEnergy + partPz) / (partEnergy - partPz + 1.e-13));
 644
 645            // Fill the AliCFContainer
 646            Double_t lContainerCutVarsMC[7] = {0.0};
 647            lContainerCutVarsMC[0]  = partP;
 648            lContainerCutVarsMC[1]  = partPt;
 649            lContainerCutVarsMC[2]  = partRap;
 650            lContainerCutVarsMC[3]  = partEta;
 651            lContainerCutVarsMC[4]  = partTheta;
 652            lContainerCutVarsMC[5]  = partPhi;
 653            lContainerCutVarsMC[6]  = lcentrality;
 654            if (PDGcode == 3312)  fCFContCascadeMCgen->Fill(lContainerCutVarsMC,0); // for Xi-
 655            if (PDGcode == -3312) fCFContCascadeMCgen->Fill(lContainerCutVarsMC,1); // for Xi+
 656            if (PDGcode == 3334)  fCFContCascadeMCgen->Fill(lContainerCutVarsMC,2); // for Omega-
 657            if (PDGcode == -3334) fCFContCascadeMCgen->Fill(lContainerCutVarsMC,3); // for Omega+
 658
 659       }
 660   }
 661
 662
 663   //////////////////////////////
 664   // CASCADE RECONSTRUCTION PART
 665   //////////////////////////////
 666
 667   //%%%%%%%%%%%%%
 668   // Cascade loop
 669   Int_t ncascades = 0;
 670   if      (fAnalysisType == "ESD") ncascades = lESDevent->GetNumberOfCascades();
 671   else if (fAnalysisType == "AOD") ncascades = lAODevent->GetNumberOfCascades();
 672
 673   for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop (ESD or AOD)
 674
 675     // -------------------------------------
 676     // - Initialisation of the local variables that will be needed for ESD/AOD
 677     // -- Container variables (1st round)
 678     Double_t lDcaXiDaughters              = -1. ;                   //[Container]
 679     Double_t lXiCosineOfPointingAngle     = -1. ;                   //[Container]
 680     Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };             //Useful to define other variables: radius fid. vol., ctau, etc. for cascade
 681     Double_t lXiRadius                    = -1000. ;                //[Container]
 682     UShort_t lPosTPCClusters              = -1;                     //To check the quality of the tracks. For ESD only ...
 683     UShort_t lNegTPCClusters              = -1;                     //To check the quality of the tracks. For ESD only ...
 684     UShort_t lBachTPCClusters             = -1;                     //To check the quality of the tracks. For ESD only ...
 685     Double_t lInvMassLambdaAsCascDghter   = 0.;                     //[Container]
 686     Double_t lDcaV0DaughtersXi            = -1.;                    //[Container]
 687     Double_t lDcaBachToPrimVertexXi       = -1.;                    //[Container]
 688     Double_t lDcaV0ToPrimVertexXi         = -1.;                    //[Container]
 689     Double_t lDcaPosToPrimVertexXi        = -1.;                    //[Container]
 690     Double_t lDcaNegToPrimVertexXi        = -1.;                    //[Container]
 691     Double_t lV0CosineOfPointingAngle     = -1.;                    //[Container]
 692     Double_t lV0toXiCosineOfPointingAngle = -1.;                    //[Container]
 693     Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. };             //Useful to define other variables: radius fid. vol., ctau, etc. for VO
 694     Double_t lV0RadiusXi                  = -1000.0;                //[Container]
 695     Double_t lV0quality                   = 0.;                     //  ??
 696     Double_t lInvMassXiMinus              = 0.;                     //[Container]
 697     Double_t lInvMassXiPlus               = 0.;                     //[Container]
 698     Double_t lInvMassOmegaMinus           = 0.;                     //[Container]
 699     Double_t lInvMassOmegaPlus            = 0.;                     //[Container]
 700     // -- PID treatment
 701     Bool_t   lIsBachelorKaonForTPC = kFALSE;
 702     Bool_t   lIsBachelorPionForTPC = kFALSE;
 703     Bool_t   lIsNegPionForTPC      = kFALSE;
 704     Bool_t   lIsPosPionForTPC      = kFALSE;
 705     Bool_t   lIsNegProtonForTPC    = kFALSE;
 706     Bool_t   lIsPosProtonForTPC    = kFALSE;
 707     // -- More container variables and quality checks
 708     Double_t lXiMomX           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
 709     Double_t lXiMomY           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
 710     Double_t lXiMomZ           = 0.;                               //Useful to define other variables: lXiTransvMom, lXiTotMom
 711     Double_t lXiTransvMom      = 0.;                               //[Container]
 712     Double_t lXiTotMom         = 0.;                               //Useful to define other variables: cTau
 713     Double_t lV0PMomX          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
 714     Double_t lV0PMomY          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
 715     Double_t lV0PMomZ          = 0.;                               //Useful to define other variables: lV0TotMom, lpTrackTransvMom
 716     Double_t lV0NMomX          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
 717     Double_t lV0NMomY          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
 718     Double_t lV0NMomZ          = 0.;                               //Useful to define other variables: lV0TotMom, lnTrackTransvMom
 719     Double_t lV0TotMom         = 0.;                               //Useful to define other variables: lctauV0
 720     Double_t lBachMomX         = 0.;                               //Useful to define other variables: lBachTransvMom
 721     Double_t lBachMomY         = 0.;                               //Useful to define other variables: lBachTransvMom
 722     Double_t lBachMomZ         = 0.;                               //Useful to define other variables: lBachTransvMom
 723     Double_t lBachTransvMom    = 0.;                               //Selection on the min bachelor pT
 724     Double_t lpTrackTransvMom  = 0.;                               //Selection on the min bachelor pT
 725     Double_t lnTrackTransvMom  = 0.;                               //Selection on the min bachelor pT
 726     Short_t  lChargeXi         = -2;                               //Useful to select the particles based on the charge
 727     Double_t lRapXi            = -20.0;                            //[Container]
 728     Double_t lRapOmega         = -20.0;                            //[Container]
 729     Float_t  etaBach           = 0.;                               //Selection on the eta range
 730     Float_t  etaPos            = 0.;                               //Selection on the eta range
 731     Float_t  etaNeg            = 0.;                               //Selection on the eta range
 732     // --  variables for the AliCFContainer dedicated to cascade cut optmisiation: ESD and AOD
 733     if (fAnalysisType == "ESD") {
 734
 735           // -------------------------------------
 736           // - Load the cascades from the handler
 737           AliESDcascade *xi = lESDevent->GetCascade(iXi);
 738           if (!xi) continue;
 739
 740           // ---------------------------------------------------------------------------
 741           // - Assigning the necessary variables for specific AliESDcascade data members
 742           lV0quality = 0.;
 743           xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
 744           lDcaXiDaughters          = xi->GetDcaXiDaughters();
 745           lXiCosineOfPointingAngle = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
 746                                        // Take care : the best available vertex should be used (like in AliCascadeVertexer)
 747           xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
 748           lXiRadius             = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
 749
 750           // -------------------------------------------------------------------------------------------------------------------------------
 751           // - Around the tracks : Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
 752           UInt_t lIdxPosXi      = (UInt_t) TMath::Abs( xi->GetPindex() );
 753           UInt_t lIdxNegXi      = (UInt_t) TMath::Abs( xi->GetNindex() );
 754           UInt_t lBachIdx       = (UInt_t) TMath::Abs( xi->GetBindex() );
 755                                     // Care track label can be negative in MC production (linked with the track quality)
 756                                     // However = normally, not the case for track index ...
 757           // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
 758           if (lBachIdx == lIdxNegXi) continue;
 759           if (lBachIdx == lIdxPosXi) continue;
 760           // - Get the track for the daughters
 761           AliESDtrack *pTrackXi         = lESDevent->GetTrack( lIdxPosXi );
 762           AliESDtrack *nTrackXi         = lESDevent->GetTrack( lIdxNegXi );
 763           AliESDtrack *bachTrackXi      = lESDevent->GetTrack( lBachIdx );
 764           if (!pTrackXi || !nTrackXi || !bachTrackXi )  continue;
 765           // - Get the TPCnumber of cluster for the daughters
 766           lPosTPCClusters   = pTrackXi->GetTPCNcls();
 767           lNegTPCClusters   = nTrackXi->GetTPCNcls();
 768           lBachTPCClusters  = bachTrackXi->GetTPCNcls();
 769
 770           // ------------------------------------
 771           // - Rejection of a poor quality tracks
 772           if (fkQualityCutTPCrefit) {
 773                 // 1 - Poor quality related to TPCrefit
 774                 ULong_t pStatus    = pTrackXi->GetStatus();
 775                 ULong_t nStatus    = nTrackXi->GetStatus();
 776                 ULong_t bachStatus = bachTrackXi->GetStatus();
 777                 if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
 778                 if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
 779                 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning(" Bach.   track has no TPCrefit ... continue!"); continue; }
 780           }
 781           if (fkQualityCutnTPCcls) {
 782                 // 2 - Poor quality related to TPC clusters
 783                 if (lPosTPCClusters  < fMinnTPCcls) { AliWarning(" V0 Pos. track has less than minn TPC clusters ... continue!"); continue; }
 784                 if (lNegTPCClusters  < fMinnTPCcls) { AliWarning(" V0 Neg. track has less than minn TPC clusters ... continue!"); continue; }
 785                 if (lBachTPCClusters < fMinnTPCcls) { AliWarning(" Bach.   track has less than minn TPC clusters ... continue!"); continue; }
 786           }
 787
 788           // ------------------------------
 789           etaPos  = pTrackXi->Eta();
 790           etaNeg  = nTrackXi->Eta();
 791           etaBach = bachTrackXi->Eta();
 792           lInvMassLambdaAsCascDghter = xi->GetEffMass();
 793           lDcaV0DaughtersXi          = xi->GetDcaV0Daughters();
 794           lV0CosineOfPointingAngle   = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
 795           lDcaV0ToPrimVertexXi       = xi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
 796           lDcaBachToPrimVertexXi     = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1],     lMagneticField) );
 797           xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
 798           lV0RadiusXi                = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
 799           lDcaPosToPrimVertexXi      = TMath::Abs( pTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) );
 800           lDcaNegToPrimVertexXi      = TMath::Abs( nTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) );
 801
 802            //----------------------------------------------------------------------------------------------------
 803            // - Around effective masses. Change mass hypotheses to cover all the possibilities:  Xi-/+, Omega -/+
 804            if (bachTrackXi->Charge() < 0) {
 805                 // Calculate the effective mass of the Xi- candidate. pdg code 3312 = Xi-
 806                 lV0quality = 0.;
 807                 xi->ChangeMassHypothesis(lV0quality , 3312);
 808                 lInvMassXiMinus = xi->GetEffMassXi();
 809                 // Calculate the effective mass of the Xi- candidate. pdg code 3334 = Omega-
 810                 lV0quality = 0.;
 811                 xi->ChangeMassHypothesis(lV0quality , 3334);
 812                 lInvMassOmegaMinus = xi->GetEffMassXi();
 813                 // Back to default hyp.
 814                 lV0quality = 0.;
 815                 xi->ChangeMassHypothesis(lV0quality , 3312);
 816            }// end if negative bachelor
 817            if ( bachTrackXi->Charge() >  0 ) {
 818                 // Calculate the effective mass of the Xi+ candidate. pdg code -3312 = Xi+
 819                 lV0quality = 0.;
 820                 xi->ChangeMassHypothesis(lV0quality , -3312);
 821                 lInvMassXiPlus = xi->GetEffMassXi();
 822                 // Calculate the effective mass of the Xi+ candidate. pdg code -3334  = Omega+
 823                 lV0quality = 0.;
 824                 xi->ChangeMassHypothesis(lV0quality , -3334);
 825                 lInvMassOmegaPlus = xi->GetEffMassXi();
 826                 // Back to "default" hyp.
 827                 lV0quality = 0.;
 828                 xi->ChangeMassHypothesis(lV0quality , -3312);
 829            }// end if positive bachelor
 830
 831            // ----------------------------------------------
 832            // - TPC PID : 3-sigma bands on Bethe-Bloch curve
 833            // Bachelor
 834            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
 835            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
 836            // Negative V0 daughter
 837            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
 838            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
 839            // Positive V0 daughter
 840            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
 841            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
 842
 843            // ------------------------------
 844            // - Miscellaneous pieces of info that may help regarding data quality assessment.
 845            xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
 846            lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
 847            lXiTotMom    = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
 848            xi->GetNPxPyPz(lV0NMomX,lV0NMomY,lV0NMomZ);
 849            xi->GetPPxPyPz(lV0PMomX,lV0PMomY,lV0PMomZ);
 850            lV0TotMom = TMath::Sqrt(TMath::Power(lV0NMomX+lV0PMomX,2)+TMath::Power(lV0NMomY+lV0PMomY,2)+TMath::Power(lV0NMomZ+lV0PMomZ,2));
 851            xi->GetBPxPyPz( lBachMomX, lBachMomY, lBachMomZ );
 852            lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
 853            lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX + lV0NMomY*lV0NMomY );
 854            lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX + lV0PMomY*lV0PMomY );
 855            lChargeXi = xi->Charge();
 856            lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
 857            lRapXi    = xi->RapXi();
 858            lRapOmega = xi->RapOmega();
 859
 860     } else if (fAnalysisType == "AOD") {
 861
 862            // -------------------------------------
 863            // - Load the cascades from the handler
 864            const AliAODcascade *xi = lAODevent->GetCascade(iXi);
 865            if (!xi) continue;
 866
 867            //----------------------------------------------------------------------------
 868            // - Assigning the necessary variables for specific AliESDcascade data members
 869            lDcaXiDaughters              = xi->DcaXiDaughters();
 870            lXiCosineOfPointingAngle     = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0],
 871                                                           lBestPrimaryVtxPos[1],
 872                                                           lBestPrimaryVtxPos[2] );
 873            lPosXi[0] = xi->DecayVertexXiX();
 874            lPosXi[1] = xi->DecayVertexXiY();
 875            lPosXi[2] = xi->DecayVertexXiZ();
 876            lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
 877
 878            //-------------------------------------------------------------------------------------------------------------------------------
 879            // - Around the tracks: Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
 880            AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
 881            AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
 882            AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
 883            if (!pTrackXi || !nTrackXi || !bachTrackXi ) continue;
 884            UInt_t lIdxPosXi  = (UInt_t) TMath::Abs( pTrackXi->GetID() );
 885            UInt_t lIdxNegXi  = (UInt_t) TMath::Abs( nTrackXi->GetID() );
 886            UInt_t lBachIdx   = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
 887                                 // Care track label can be negative in MC production (linked with the track quality)
 888                                 // However = normally, not the case for track index ...
 889
 890            // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
 891            if (lBachIdx == lIdxNegXi) continue;
 892            if (lBachIdx == lIdxPosXi) continue;
 893            // - Get the TPCnumber of cluster for the daughters
 894            lPosTPCClusters   = pTrackXi->GetTPCNcls();
 895            lNegTPCClusters   = nTrackXi->GetTPCNcls();
 896            lBachTPCClusters  = bachTrackXi->GetTPCNcls();
 897
 898            // ------------------------------------
 899            // - Rejection of a poor quality tracks
 900            if (fkQualityCutTPCrefit) {
 901                 // - Poor quality related to TPCrefit
 902                 if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit)))    { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
 903                 if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit)))    { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
 904                 if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning(" Bach.   track has no TPCrefit ... continue!"); continue; }
 905            }
 906            if (fkQualityCutnTPCcls) {
 907                 // - Poor quality related to TPC clusters
 908                 if (lPosTPCClusters  < fMinnTPCcls) continue;
 909                 if (lNegTPCClusters  < fMinnTPCcls) continue;
 910                 if (lBachTPCClusters < fMinnTPCcls) continue;
 911            }
 912
 913            // ------------------------------------------------------------------------------------------------------------------------------
 914            // - Around the tracks: Bach + V0 (AOD). Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
 915            etaPos  = pTrackXi->Eta();
 916            etaNeg  = nTrackXi->Eta();
 917            etaBach = bachTrackXi->Eta();
 918            lChargeXi                    = xi->ChargeXi();
 919            if ( lChargeXi < 0)  lInvMassLambdaAsCascDghter      = xi->MassLambda();
 920            else                 lInvMassLambdaAsCascDghter      = xi->MassAntiLambda();
 921            lDcaV0DaughtersXi            = xi->DcaV0Daughters();
 922            lDcaV0ToPrimVertexXi                 = xi->DcaV0ToPrimVertex();
 923            lDcaBachToPrimVertexXi               = xi->DcaBachToPrimVertex();
 924            lPosV0Xi[0] = xi->DecayVertexV0X();
 925            lPosV0Xi[1] = xi->DecayVertexV0Y();
 926            lPosV0Xi[2] = xi->DecayVertexV0Z();
 927            lV0RadiusXi  = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
 928            lV0CosineOfPointingAngle     = xi->CosPointingAngle( lBestPrimaryVtxPos );
 929            lDcaPosToPrimVertexXi        = xi->DcaPosToPrimVertex();
 930            lDcaNegToPrimVertexXi        = xi->DcaNegToPrimVertex();
 931
 932            // ---------------------------------------------------------------------------------------------------
 933            // - Around effective masses. Change mass hypotheses to cover all the possibilities:  Xi-/+, Omega -/+
 934            if ( lChargeXi < 0 )         lInvMassXiMinus         = xi->MassXi();
 935            if ( lChargeXi > 0 )         lInvMassXiPlus          = xi->MassXi();
 936            if ( lChargeXi < 0 )         lInvMassOmegaMinus      = xi->MassOmega();
 937            if ( lChargeXi > 0 )         lInvMassOmegaPlus       = xi->MassOmega();
 938
 939            // ----------------------------------------------
 940            // - TPC PID : 3-sigma bands on Bethe-Bloch curve
 941            // Bachelor
 942            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
 943            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
 944            // Negative V0 daughter
 945            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
 946            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
 947            // Positive V0 daughter
 948            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
 949            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
 950
 951            //---------------------------------
 952            // - Extra info for QA (AOD)
 953            // Miscellaneous pieces of info that may help regarding data quality assessment.
 954            // Cascade transverse and total momentum
 955            lXiMomX = xi->MomXiX();
 956            lXiMomY = xi->MomXiY();
 957            lXiMomZ = xi->MomXiZ();
 958            lXiTransvMom         = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
 959            lXiTotMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
 960            Double_t lV0MomX = xi->MomV0X();
 961            Double_t lV0MomY = xi->MomV0Y();
 962            Double_t lV0MomZ = xi->MomV0Z();
 963            lV0TotMom = TMath::Sqrt(TMath::Power(lV0MomX,2)+TMath::Power(lV0MomY,2)+TMath::Power(lV0MomZ,2));
 964            lBachMomX = xi->MomBachX();
 965            lBachMomY = xi->MomBachY();
 966            lBachMomZ = xi->MomBachZ();
 967            lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
 968            lV0NMomX = xi->MomNegX();
 969            lV0NMomY = xi->MomNegY();
 970            lV0PMomX = xi->MomPosX();
 971            lV0PMomY = xi->MomPosY();
 972            lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX   + lV0NMomY*lV0NMomY );
 973            lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX   + lV0PMomY*lV0PMomY );
 974            lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
 975            lRapXi    = xi->RapXi();
 976            lRapOmega = xi->RapOmega();
 977
 978     }// end of AOD treatment
 979
 980     //---------------------------------------
 981     // Cut on pt of the three daughter tracks
 982     if (lBachTransvMom<fMinPtCutOnDaughterTracks) continue;
 983     if (lpTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
 984     if (lnTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
 985
 986     //---------------------------------------------------
 987     // Cut on pseudorapidity of the three daughter tracks
 988     if (TMath::Abs(etaBach)>fEtaCutOnDaughterTracks) continue;
 989     if (TMath::Abs(etaPos)>fEtaCutOnDaughterTracks) continue;
 990     if (TMath::Abs(etaNeg)>fEtaCutOnDaughterTracks) continue;
 991
 992     //----------------------------------
 993     // Calculate proper time for cascade
 994     Double_t cascadeMass = 0.;
 995     if ( ( (lChargeXi<0) && lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) ||
 996          ( (lChargeXi>0) && lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC )  ) cascadeMass = 1.321;
 997     if ( ( (lChargeXi<0) && lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC ) ||
 998          ( (lChargeXi>0) && lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )  ) cascadeMass = 1.672;
 999     Double_t lctau =  TMath::Sqrt(TMath::Power((lPosXi[0]-lBestPrimaryVtxPos[0]),2)+TMath::Power((lPosXi[1]-lBestPrimaryVtxPos[1]),2)+TMath::Power(( lPosXi[2]-lBestPrimaryVtxPos[2]),2));
1000     if (lXiTotMom!=0)         lctau = lctau*cascadeMass/lXiTotMom;
1001     else lctau = -1.;
1002     // Calculate proper time for Lambda (reconstructed)
1003     Float_t lambdaMass = 1.115683; // PDG mass
1004     Float_t distV0Xi =  TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2)+TMath::Power((lPosV0Xi[2]-lPosXi[2]),2));
1005     Float_t lctauV0 = -1.;
1006     if (lV0TotMom!=0) lctauV0 = distV0Xi*lambdaMass/lV0TotMom;
1007
1008
1009     // Fill the AliCFContainer (optimisation of topological selections)
1010     Double_t lContainerCutVars[21] = {0.0};
1011     lContainerCutVars[0]  = lDcaXiDaughters;
1012     lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
1013     lContainerCutVars[2]  = lXiCosineOfPointingAngle;
1014     lContainerCutVars[3]  = lXiRadius;
1015     lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
1016     lContainerCutVars[5]  = lDcaV0DaughtersXi;
1017     lContainerCutVars[6]  = lV0CosineOfPointingAngle;
1018     lContainerCutVars[7]  = lV0RadiusXi;
1019     lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
1020     lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
1021     lContainerCutVars[10] = lDcaNegToPrimVertexXi;
1022     lContainerCutVars[13] = lXiTransvMom;
1023     lContainerCutVars[16] = lctau;
1024     lContainerCutVars[17] = lctauV0;
1025     lContainerCutVars[18] = lV0toXiCosineOfPointingAngle;
1026     lContainerCutVars[19] = lcentrality;
1027     if ( lChargeXi < 0 ) {
1028          lContainerCutVars[11] = lInvMassXiMinus;
1029          lContainerCutVars[12] = lInvMassOmegaMinus;
1030          lContainerCutVars[14] = lRapXi;
1031          lContainerCutVars[15] = -1.;
1032          if (lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
1033          lContainerCutVars[11] = lInvMassXiMinus;
1034          lContainerCutVars[12] = lInvMassOmegaMinus;
1035          lContainerCutVars[14] = -1.;
1036          lContainerCutVars[15] = lRapOmega;
1037          if (lIsBachelorKaonForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
1038     } else {
1039          lContainerCutVars[11] = lInvMassXiPlus;
1040          lContainerCutVars[12] = lInvMassOmegaPlus;
1041          lContainerCutVars[14] = lRapXi;
1042          lContainerCutVars[15] = -1.;
1043          if (lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
1044          lContainerCutVars[11] = lInvMassXiPlus;
1045          lContainerCutVars[12] = lInvMassOmegaPlus;
1046          lContainerCutVars[14] = -1.;
1047          lContainerCutVars[15] = lRapOmega;
1048          if (lIsBachelorKaonForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
1049     }
1050
1051
1052   }// end of the Cascade loop (ESD or AOD)
1053
1054
1055   // Post output data.
1056   PostData(1, fCFContCascadeCuts);
1057   PostData(2, fCFContCascadeMCgen);
1058 }
1059
1060 //________________________________________________________________________
1061 void AliAnalysisTaskQAMultistrange::Terminate(Option_t *) {
1062
1063 }
1064