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