]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/RESONANCES/extra/AliXiStar.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / extra / AliXiStar.cxx
CommitLineData
3311ad64 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////////////////////////////////////
17//
18// This class is used to reconstruct the neutral Xi(1530) resonance.
19// This class essentially combines charged Xi candidates from the Xi Vertexer
20// with primary charged pions.
21//
22// authors: Dhevan Gangadharan (dhevan.raja.gangadharan@cern.ch)
23//
24////////////////////////////////////////////////////////////////////////////////
25
26
27
28#include <iostream>
29#include <math.h>
30#include "TChain.h"
31#include "TFile.h"
32#include "TKey.h"
33#include "TObject.h"
34#include "TObjString.h"
35#include "TList.h"
36#include "TTree.h"
37#include "TH1F.h"
38#include "TH1D.h"
39#include "TH2D.h"
40#include "TH3D.h"
41#include "TProfile.h"
42#include "TCanvas.h"
43
44#include "AliAnalysisTask.h"
45#include "AliAnalysisManager.h"
46
47
48#include "AliESDEvent.h"
49#include "AliESDInputHandler.h"
50#include "AliESDtrackCuts.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54
55#include "AliAODEvent.h"
56#include "AliAODInputHandler.h"
57#include "AliAODMCParticle.h"
58#include "AliAODcascade.h"
59#include "AliESDcascade.h"
60#include "AliV0vertexer.h"
61#include "AliCascadeVertexer.h"
62
63#include "AliXiStar.h"
64
65#define PI 3.1415927
66
67
68// Author: Dhevan Gangadharan
69
70ClassImp(AliXiStar)
71
72//________________________________________________________________________
73AliXiStar::AliXiStar():
74AliAnalysisTaskSE(),
75 fname(0),
76 fAOD(0x0),
77 fESD(0x0),
78 fOutputList(0x0),
79 fTrackCut(0x0),
80 fPIDResponse(0x0),
81 fEC(0x0),
82 fEvt(0x0),
83 fTempStruct(0x0),
84 fZvertexBins(0),
85 fEventsToMix(0),
86 fMultBins(0),
3311ad64 87 fMCcase(0),
88 fAODcase(0),
89 fEventCounter(0),
3311ad64 90 fMaxDecayLength(0),
3311ad64 91 fMassWindow(0),
3311ad64 92 fTrueMassPr(0),
93 fTrueMassPi(0),
94 fTrueMassK(0),
95 fTrueMassLam(0),
96 fTrueMassXi(0),
97 fESDTrack4(0x0),
98 fXiTrack(0x0),
0a21e6f4 99 fCutList(0)
3311ad64 100{
0a21e6f4 101 // Default Constructor
102 for (Int_t i=0; i<21; i++){
103 fCovMatrix[i]=-99999.;
104 if (i<12) fMultLimits[i] = 0;
105 }
106 for (Int_t i=0; i<kNCuts; i++){
107 fDecayParameters[i]=0;
108 for (Int_t j=0; j<kNCutVariations; j++){
109 fCutValues[j][i]=0;
110 }
111 }
112 //
113 for (Int_t cv=0; cv<kNCutVariations; cv++){
114 CutVar[cv].fXi=0x0;
115 CutVar[cv].fXibar=0x0;
116 CutVar[cv].fXiMinusPiPlus=0x0;
117 CutVar[cv].fXiMinusPiMinus=0x0;
118 CutVar[cv].fXiPlusPiPlus=0x0;
119 CutVar[cv].fXiPlusPiMinus=0x0;
120 //
121 CutVar[cv].fXiMinusPiPlusbkg=0x0;
122 CutVar[cv].fXiMinusPiMinusbkg=0x0;
123 CutVar[cv].fXiPlusPiPlusbkg=0x0;
124 CutVar[cv].fXiPlusPiMinusbkg=0x0;
125 //
126 CutVar[cv].fMCrecXi=0x0;
127 CutVar[cv].fMCrecXibar=0x0;
128 CutVar[cv].fMCrecXiMinusPiPlus=0x0;
129 CutVar[cv].fMCrecXiPlusPiMinus=0x0;
130 }
131
3311ad64 132}
133//________________________________________________________________________
134AliXiStar::AliXiStar(const char *name, Bool_t AODdecision, Bool_t MCdecision, Int_t CutListOption)
135 : AliAnalysisTaskSE(name),
136 fname(name),
137 fAOD(0x0),
138 fESD(0x0),
139 fOutputList(0x0),
140 fTrackCut(0x0),
141 fPIDResponse(0x0),
142 fEC(0x0),
143 fEvt(0x0),
144 fTempStruct(0x0),
145 fZvertexBins(0),
146 fEventsToMix(0),
147 fMultBins(0),
3311ad64 148 fMCcase(MCdecision),
149 fAODcase(AODdecision),
150 fEventCounter(0),
3311ad64 151 fMaxDecayLength(0),
3311ad64 152 fMassWindow(0),
3311ad64 153 fTrueMassPr(0),
154 fTrueMassPi(0),
155 fTrueMassK(0),
156 fTrueMassLam(0),
157 fTrueMassXi(0),
158 fESDTrack4(0x0),
159 fXiTrack(0x0),
160 fCutList(CutListOption)
fe45e21b 161
3311ad64 162{
163 // Main Constructor
92650a3e 164 for (Int_t i=0; i<21; i++){
165 fCovMatrix[i]=-99999.;
166 if (i<12) fMultLimits[i] = 0;
167 }
fe45e21b 168 for (Int_t i=0; i<kNCuts; i++){
169 fDecayParameters[i]=0;
170 for (Int_t j=0; j<kNCutVariations; j++){
171 fCutValues[j][i]=0;
172 }
173 }
0a21e6f4 174 //
175 for (Int_t cv=0; cv<kNCutVariations; cv++){
176 CutVar[cv].fXi=0x0;
177 CutVar[cv].fXibar=0x0;
178 CutVar[cv].fXiMinusPiPlus=0x0;
179 CutVar[cv].fXiMinusPiMinus=0x0;
180 CutVar[cv].fXiPlusPiPlus=0x0;
181 CutVar[cv].fXiPlusPiMinus=0x0;
182 //
183 CutVar[cv].fXiMinusPiPlusbkg=0x0;
184 CutVar[cv].fXiMinusPiMinusbkg=0x0;
185 CutVar[cv].fXiPlusPiPlusbkg=0x0;
186 CutVar[cv].fXiPlusPiMinusbkg=0x0;
187 //
188 CutVar[cv].fMCrecXi=0x0;
189 CutVar[cv].fMCrecXibar=0x0;
190 CutVar[cv].fMCrecXiMinusPiPlus=0x0;
191 CutVar[cv].fMCrecXiPlusPiMinus=0x0;
192 }
193
fe45e21b 194
3311ad64 195 // Define output slots here
196 // Output slot #1
197 DefineOutput(1, TList::Class());
198
199}
200//________________________________________________________________________
201AliXiStar::AliXiStar(const AliXiStar &obj)
202 : AliAnalysisTaskSE(obj.fname),
203 fname(obj.fname),
204 fAOD(obj.fAOD),
205 fESD(obj.fESD),
206 fOutputList(obj.fOutputList),
207 fTrackCut(obj.fTrackCut),
208 fPIDResponse(obj.fPIDResponse),
209 fEC(obj.fEC),
210 fEvt(obj.fEvt),
211 fTempStruct(obj.fTempStruct),
212 fZvertexBins(obj.fZvertexBins),
213 fEventsToMix(obj.fEventsToMix),
214 fMultBins(obj.fMultBins),
215 fMultLimits(),
216 fMCcase(obj.fMCcase),
217 fAODcase(obj.fAODcase),
218 fEventCounter(obj.fEventCounter),
3311ad64 219 fMaxDecayLength(obj.fMaxDecayLength),
3311ad64 220 fMassWindow(obj.fMassWindow),
3311ad64 221 fTrueMassPr(obj.fTrueMassPr),
222 fTrueMassPi(obj.fTrueMassPi),
223 fTrueMassK(obj.fTrueMassK),
224 fTrueMassLam(obj.fTrueMassLam),
225 fTrueMassXi(obj.fTrueMassXi),
226 fESDTrack4(obj.fESDTrack4),
227 fXiTrack(obj.fXiTrack),
228 fCutList(obj.fCutList)
229
230{
231 // Copy constructor
92650a3e 232 for (Int_t i=0; i<21; i++){
233 fCovMatrix[i]=obj.fCovMatrix[i];
234 if (i<12) fMultLimits[i]=obj.fMultLimits[i];
fe45e21b 235 }
236 for (Int_t i=0; i<kNCuts; i++){
237 fDecayParameters[i]=obj.fDecayParameters[i];
238 for (Int_t j=0; j<kNCutVariations; j++){
239 fCutValues[j][i]=obj.fCutValues[j][i];
240 }
241 }
242
3311ad64 243}
244//________________________________________________________________________
245AliXiStar &AliXiStar::operator=(const AliXiStar &obj)
246{
247 // Assignment operator
248 if (this == &obj)
249 return *this;
92650a3e 250
251 fname = obj.fname;
3311ad64 252 fAOD = obj.fAOD;
253 fESD = obj.fESD;
254 fOutputList = obj.fOutputList;
255 fTrackCut = obj.fTrackCut;
256 fPIDResponse = obj.fPIDResponse;
257 fEC = obj.fEC;
258 fEvt = obj.fEvt;
259 fTempStruct = obj.fTempStruct;
260 fZvertexBins = obj.fZvertexBins;
261 fEventsToMix = obj.fEventsToMix;
262 fMultBins = obj.fMultBins;
92650a3e 263 for (Int_t i=0; i<12; i++){
264 fMultLimits[i]=obj.fMultLimits[i];
265 }
3311ad64 266 fMCcase = obj.fMCcase;
267 fAODcase = obj.fAODcase;
268 fEventCounter = obj.fEventCounter;
3311ad64 269 fMaxDecayLength = obj.fMaxDecayLength;
3311ad64 270 fMassWindow = obj.fMassWindow;
92650a3e 271 for (Int_t i=0; i<21; i++){
272 fCovMatrix[i]=obj.fCovMatrix[i];
273 }
3311ad64 274 fTrueMassPr = obj.fTrueMassPr;
275 fTrueMassPi = obj.fTrueMassPi;
276 fTrueMassK = obj.fTrueMassK;
277 fTrueMassLam = obj.fTrueMassLam;
278 fTrueMassXi = obj.fTrueMassXi;
279 fESDTrack4 = obj.fESDTrack4;
280 fXiTrack = obj.fXiTrack;
281 fCutList = obj.fCutList;
fe45e21b 282
283 for (Int_t i=0; i<kNCuts; i++){
284 fDecayParameters[i]=obj.fDecayParameters[i];
285 for (Int_t j=0; j<kNCutVariations; j++){
286 fCutValues[j][i]=obj.fCutValues[j][i];
287 }
288 }
289
290
3311ad64 291 return (*this);
292}
293//________________________________________________________________________
294AliXiStar::~AliXiStar()
295{
296 // Destructor
3311ad64 297 if(fAOD) delete fAOD;
298 if(fESD) delete fESD;
299 if(fOutputList) delete fOutputList;
300 if(fTrackCut) delete fTrackCut;
301 if(fPIDResponse) delete fPIDResponse;
302 if(fEC) delete fEC;
303 if(fEvt) delete fEvt;
304 if(fTempStruct) delete fTempStruct;
305 if(fESDTrack4) delete fESDTrack4;
306 if(fXiTrack) delete fXiTrack;
fe45e21b 307
308 for (Int_t cv=0; cv<kNCutVariations; cv++){
309 if(CutVar[cv].fXi) delete CutVar[cv].fXi;
310 if(CutVar[cv].fXibar) delete CutVar[cv].fXibar;
311 if(CutVar[cv].fXiMinusPiPlus) delete CutVar[cv].fXiMinusPiPlus;
312 if(CutVar[cv].fXiMinusPiMinus) delete CutVar[cv].fXiMinusPiMinus;
313 if(CutVar[cv].fXiPlusPiPlus) delete CutVar[cv].fXiPlusPiPlus;
314 if(CutVar[cv].fXiPlusPiMinus) delete CutVar[cv].fXiPlusPiMinus;
315 //
316 if(CutVar[cv].fXiMinusPiPlusbkg) delete CutVar[cv].fXiMinusPiPlusbkg;
317 if(CutVar[cv].fXiMinusPiMinusbkg) delete CutVar[cv].fXiMinusPiMinusbkg;
318 if(CutVar[cv].fXiPlusPiPlusbkg) delete CutVar[cv].fXiPlusPiPlusbkg;
319 if(CutVar[cv].fXiPlusPiMinusbkg) delete CutVar[cv].fXiPlusPiMinusbkg;
320 //
321 if(CutVar[cv].fMCrecXi) delete CutVar[cv].fMCrecXi;
322 if(CutVar[cv].fMCrecXibar) delete CutVar[cv].fMCrecXibar;
323 if(CutVar[cv].fMCrecXiMinusPiPlus) delete CutVar[cv].fMCrecXiMinusPiPlus;
324 if(CutVar[cv].fMCrecXiPlusPiMinus) delete CutVar[cv].fMCrecXiPlusPiMinus;
325 }
326
3311ad64 327}
328//________________________________________________________________________
329void AliXiStar::XiStarInit()
330{
331 //
332 //Inits cuts and analysis settings
333 //
334
335 fEventCounter=0;// event counter initialization
336 cout<<"AliXiStar XiStarInit() call"<<endl;
337
338
339 ///////////////////////////////////////////////
340 // Track Cuts for ESD analysis
341 fTrackCut = new AliESDtrackCuts();
342 fTrackCut->SetPtRange(.15,1000);
343 fTrackCut->SetAcceptKinkDaughters(kFALSE);
fe45e21b 344 //fTrackCut->SetMinNClustersTPC(70);
3311ad64 345 fTrackCut->SetRequireTPCRefit(kTRUE);
346 ////////////////////////////////////////////////
347
348 fZvertexBins = 20;
349 fMultBins = 11;// This must also be set in AliXiStar.h
350 if(fMCcase) fEventsToMix = 0;
351 else fEventsToMix = 40;
352
353 // multiplicity edges for event mixing bins
354 fMultLimits[0]=0, fMultLimits[1]=5, fMultLimits[2]=10, fMultLimits[3]=15, fMultLimits[4]=20, fMultLimits[5]=25;
355 fMultLimits[6]=30, fMultLimits[7]=35, fMultLimits[8]=40, fMultLimits[9]=45, fMultLimits[10]=50, fMultLimits[11]=150;
356
357
358 fEC = new AliXiStarEventCollection **[fZvertexBins];
359 for(unsigned short i=0; i<fZvertexBins; i++){
360
361 fEC[i] = new AliXiStarEventCollection *[fMultBins];
362
363 for(unsigned short j=0; j<fMultBins; j++){
364
365 fEC[i][j] = new AliXiStarEventCollection(fEventsToMix+1);
366 }
367 }
368
369
370 fTempStruct = new AliXiStarTrackStruct[kNbinsM];
371
372 fESDTrack4 = new AliESDtrack();
373 fXiTrack = new AliESDtrack();
3311ad64 374
3311ad64 375
fe45e21b 376 fMaxDecayLength = 100.;
377 fMassWindow = 0.006;
378
379 /////////////////////////////////////////////////////////////////////////////////////
380 ///////////////////////
381 // DecayParameters Key (number represents array index)
382 // NclustersTPC: 0=proton, 1=pion first, 2=pion second, 3=pion third
383 // DCAVtx: 4=proton, 5=pion first, 6=pion second, 7=lambda, 8=pion third
384 // 9 = DCA proton-pion
385 // 10 = DCA Lambda-pion
386 // 11 = Rxy Lambda
387 // 12 = Rxy Xi
388 // 13 = Cos PA Lambda
389 // 14 = Cos PA Xi
3311ad64 390
fe45e21b 391 // Set Standard Reconstruction cut values
392 fCutValues[0][0] = 70; fCutValues[0][1] = 70; fCutValues[0][2] = 70; fCutValues[0][3] = 70;
393 fCutValues[0][4] = 0.04; fCutValues[0][5] = 0.04; fCutValues[0][6] = 0.05; fCutValues[0][7] = 0.07; fCutValues[0][8] = 2.0;
394 fCutValues[0][9] = 1.6;
395 fCutValues[0][10] = 1.6;
396 fCutValues[0][11] = 1.4;
397 fCutValues[0][12] = 0.8;
398 fCutValues[0][13] = 0.97;
399 fCutValues[0][14] = 0.97;
400 for(int cv=1; cv<kNCutVariations; cv++){
401 for(int ct=0; ct<kNCuts; ct++){
402 fCutValues[cv][ct] = fCutValues[0][ct];
403 }
3311ad64 404 }
fe45e21b 405 // Set Systematic Variations
406 fCutValues[1][0] = 80; fCutValues[1][1] = 80; fCutValues[1][2] = 80; fCutValues[1][3] = 80;// 80
407 fCutValues[2][4] = 0.104;// 0.104
408 fCutValues[3][5] = 0.104;// 0.104
409 fCutValues[4][6] = 0.08;// 0.08
410 fCutValues[5][7] = 0.1;// 0.1
411 fCutValues[6][8] = 1.0;// 1.0
412 fCutValues[7][9] = 0.94;// 0.94
413 fCutValues[8][10] = 1.41;// 1.41
414 fCutValues[9][11] = 4.39;// 4.39
415 fCutValues[10][12] = 0.95;// 0.95
416 fCutValues[11][13] = 0.99;// 0.99
417 fCutValues[12][14] = 0.985;// 0.085
418
419
420
3311ad64 421
3311ad64 422
423 // PDG mass values
424 fTrueMassPr=.93827, fTrueMassPi=.13957, fTrueMassK=.493677, fTrueMassLam=1.11568, fTrueMassXi=1.32171;
425
426 // The following CovMatrix is set so that PropogateToDCA() ignores track errors. Only used to propagate Xi to third pion for XiStar reconstruction
427 for(Int_t i=0; i<21; i++) fCovMatrix[i]=0;
428 fCovMatrix[0]=1, fCovMatrix[2]=1, fCovMatrix[5]=1, fCovMatrix[9]=1, fCovMatrix[14]=1, fCovMatrix[20]=1;
429
430
431}
432//________________________________________________________________________
433void AliXiStar::UserCreateOutputObjects()
434{
435
436 XiStarInit();// Initialize settings
437
438 // Create histograms
439 fOutputList = new TList();
440 fOutputList->SetOwner();
441
442 TH3F *fVertexDist1 = new TH3F("fVertexDist1","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
443 fVertexDist1->GetXaxis()->SetTitle("X Vertex (cm)");
444 fVertexDist1->GetYaxis()->SetTitle("Y Vertex (cm)");
445 fVertexDist1->GetZaxis()->SetTitle("Z Vertex (cm)");
446 fOutputList->Add(fVertexDist1);
447
448 TH3F *fVertexDist3 = new TH3F("fVertexDist3","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
449 fVertexDist3->GetXaxis()->SetTitle("X Vertex (cm)");
450 fVertexDist3->GetYaxis()->SetTitle("Y Vertex (cm)");
451 fVertexDist3->GetZaxis()->SetTitle("Z Vertex (cm)");
452 fOutputList->Add(fVertexDist3);
453
454 TH2F *fDCADist = new TH2F("fDCADist","DCA distribution",kNbinsM,-.5,kNbinsM-.5, 100,0,10);
455 fOutputList->Add(fDCADist);
456
457
458 TH3F *fMultDist3d = new TH3F("fMultDist3d","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5, kNbinsM,-.5,kNbinsM-.5, kNbinsM,-.5,kNbinsM-.5);
459 fMultDist3d->GetXaxis()->SetTitle("Multiplicity");
460 fMultDist3d->GetYaxis()->SetTitle("Positive Multiplicity");
461 fMultDist3d->GetZaxis()->SetTitle("Negative Multiplicity");
462 fMultDist3d->SetMarkerStyle(kFullCircle);
463 fOutputList->Add(fMultDist3d);
464
465
466 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
467 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
468 fOutputList->Add(fMultDist1);
469
470 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
471 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
472 fOutputList->Add(fMultDist2);
473
474 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
475 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
476 fOutputList->Add(fMultDist3);
477
478 TH1F *fMultDist4 = new TH1F("fMultDist4","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
479 fMultDist4->GetXaxis()->SetTitle("Multiplicity");
480 fOutputList->Add(fMultDist4);
481
482 TH1F *fMultDist5 = new TH1F("fMultDist5","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
483 fMultDist5->GetXaxis()->SetTitle("Multiplicity");
484 fOutputList->Add(fMultDist5);
485
486
487 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","PtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
488 fOutputList->Add(fPtEtaDist);
489
490 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","PhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
491 fOutputList->Add(fPhiPtDist);
492
fe45e21b 493
494 for(Int_t cv=0; cv<kNCutVariations; cv++){
495
496 if(cv==0){
497 TString *nameXi=new TString("fXi_");
498 TString *nameXibar=new TString("fXibar_");
499 *nameXi += cv;
500 *nameXibar += cv;
501 CutVar[cv].fXi = new TH3F(nameXi->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
502 fOutputList->Add(CutVar[cv].fXi);
503 CutVar[cv].fXibar = new TH3F(nameXibar->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
504 fOutputList->Add(CutVar[cv].fXibar);
505 //
506 TString *nameMCrecXi = new TString("fMCrecXi_");
507 TString *nameMCrecXibar = new TString("fMCrecXi_");
508 *nameMCrecXi += cv;
509 *nameMCrecXibar += cv;
510 CutVar[cv].fMCrecXi = new TH3F(nameMCrecXi->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
511 CutVar[cv].fMCrecXibar = new TH3F(nameMCrecXibar->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
512 fOutputList->Add(CutVar[cv].fMCrecXi);
513 fOutputList->Add(CutVar[cv].fMCrecXibar);
514 }
515 //
516 TString *nameXiMinusPiPlus = new TString("fXiMinusPiPlus_");
517 TString *nameXiMinusPiMinus = new TString("fXiMinusPiMinus_");
518 TString *nameXiPlusPiPlus = new TString("fXiPlusPiPlus_");
519 TString *nameXiPlusPiMinus = new TString("fXiPlusPiMinus_");
520 TString *nameXiMinusPiPlusbkg = new TString("fXiMinusPiPlusbkg_");
521 TString *nameXiMinusPiMinusbkg = new TString("fXiMinusPiMinusbkg_");
522 TString *nameXiPlusPiPlusbkg = new TString("fXiPlusPiPlusbkg_");
523 TString *nameXiPlusPiMinusbkg = new TString("fXiPlusPiMinusbkg_");
524 *nameXiMinusPiPlus += cv;
525 *nameXiMinusPiMinus += cv;
526 *nameXiPlusPiPlus += cv;
527 *nameXiPlusPiMinus += cv;
528 *nameXiMinusPiPlusbkg += cv;
529 *nameXiMinusPiMinusbkg += cv;
530 *nameXiPlusPiPlusbkg += cv;
531 *nameXiPlusPiMinusbkg += cv;
532 CutVar[cv].fXiMinusPiPlus = new TH3F(nameXiMinusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
533 CutVar[cv].fXiMinusPiMinus = new TH3F(nameXiMinusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
534 CutVar[cv].fXiPlusPiPlus = new TH3F(nameXiPlusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
535 CutVar[cv].fXiPlusPiMinus = new TH3F(nameXiPlusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
536 CutVar[cv].fXiMinusPiPlusbkg = new TH3F(nameXiMinusPiPlusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
537 CutVar[cv].fXiMinusPiMinusbkg = new TH3F(nameXiMinusPiMinusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
538 CutVar[cv].fXiPlusPiPlusbkg = new TH3F(nameXiPlusPiPlusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
539 CutVar[cv].fXiPlusPiMinusbkg = new TH3F(nameXiPlusPiMinusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
540
541 fOutputList->Add(CutVar[cv].fXiMinusPiPlus);
542 fOutputList->Add(CutVar[cv].fXiMinusPiMinus);
543 fOutputList->Add(CutVar[cv].fXiPlusPiPlus);
544 fOutputList->Add(CutVar[cv].fXiPlusPiMinus);
545 fOutputList->Add(CutVar[cv].fXiMinusPiPlusbkg);
546 fOutputList->Add(CutVar[cv].fXiMinusPiMinusbkg);
547 fOutputList->Add(CutVar[cv].fXiPlusPiPlusbkg);
548 fOutputList->Add(CutVar[cv].fXiPlusPiMinusbkg);
549 //
3311ad64 550
fe45e21b 551
552 TString *nameMCrecXiMinusPiPlus = new TString("fMCrecXiMinusPiPlus_");
553 TString *nameMCrecXiPlusPiMinus = new TString("fMCrecXiPlusPiMinus_");
554 *nameMCrecXiMinusPiPlus += cv;
555 *nameMCrecXiPlusPiMinus += cv;
556 CutVar[cv].fMCrecXiMinusPiPlus = new TH3F(nameMCrecXiMinusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
557 CutVar[cv].fMCrecXiPlusPiMinus = new TH3F(nameMCrecXiPlusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
558 fOutputList->Add(CutVar[cv].fMCrecXiMinusPiPlus);
559 fOutputList->Add(CutVar[cv].fMCrecXiPlusPiMinus);
560 //
561
562 /*
563 CutVar[cv].fMCrecXiStarxiy = new TH2F("fMCrecXiStarxiy","y distribution",80,-2,2, 80,-2,2);
564 CutVar[cv].fMCrecXiStarpiony = new TH2F("fMCrecXiStarpiony","y distribution",80,-2,2, 80,-2,2);
565 fOutputList->Add(CutVar[cv].fMCrecXiStarxiy);
566 fOutputList->Add(CutVar[cv].fMCrecXiStarpiony);
567 CutVar[cv].fMCrecXilambday = new TH2F("fMCrecXilambday","y distribution",80,-2,2, 80,-2,2);
568 CutVar[cv].fMCrecXipiony = new TH2F("fMCrecXipiony","y distribution",80,-2,2, 80,-2,2);
569 fOutputList->Add(CutVar[cv].fMCrecXilambday);
570 fOutputList->Add(CutVar[cv].fMCrecXipiony);
571 CutVar[cv].fMCrecLamprotony = new TH2F("fMCrecLamprotony","y distribution",80,-2,2, 80,-2,2);
572 CutVar[cv].fMCrecLampiony = new TH2F("fMCrecLampiony","y distribution",80,-2,2, 80,-2,2);
573 fOutputList->Add(CutVar[cv].fMCrecLamprotony);
574 fOutputList->Add(CutVar[cv].fMCrecLampiony);
575 */
576 }
3311ad64 577
fe45e21b 578
3311ad64 579
580
fe45e21b 581 //////////////////////
582 // MC input histos
583 TH3F *fMCinputXiStar = new TH3F("fMCinputXiStar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
584 TH3F *fMCinputXiStarbar = new TH3F("fMCinputXiStarbar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
3311ad64 585 fOutputList->Add(fMCinputXiStar);
586 fOutputList->Add(fMCinputXiStarbar);
587
588 TH3F *fMCinputXi = new TH3F("fMCinputXi","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
589 TH3F *fMCinputXibar = new TH3F("fMCinputXibar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
590 fOutputList->Add(fMCinputXi);
591 fOutputList->Add(fMCinputXibar);
592
593 //
594
fe45e21b 595 TH3F *fMCinputTotalXiStar1 = new TH3F("fMCinputTotalXiStar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
596 TH3F *fMCinputTotalXiStarbar1 = new TH3F("fMCinputTotalXiStarbar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
3311ad64 597 fOutputList->Add(fMCinputTotalXiStar1);
598 fOutputList->Add(fMCinputTotalXiStarbar1);
599
600 TH3F *fMCinputTotalXi1 = new TH3F("fMCinputTotalXi1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
601 TH3F *fMCinputTotalXibar1 = new TH3F("fMCinputTotalXibar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
602 fOutputList->Add(fMCinputTotalXi1);
603 fOutputList->Add(fMCinputTotalXibar1);
604
605 //
606
fe45e21b 607 TH3F *fMCinputTotalXiStar3 = new TH3F("fMCinputTotalXiStar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
608 TH3F *fMCinputTotalXiStarbar3 = new TH3F("fMCinputTotalXiStarbar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
3311ad64 609 fOutputList->Add(fMCinputTotalXiStar3);
610 fOutputList->Add(fMCinputTotalXiStarbar3);
611
612 TH3F *fMCinputTotalXi3 = new TH3F("fMCinputTotalXi3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
613 TH3F *fMCinputTotalXibar3 = new TH3F("fMCinputTotalXibar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
614 fOutputList->Add(fMCinputTotalXi3);
615 fOutputList->Add(fMCinputTotalXibar3);
616
617 //
3311ad64 618
619 TH2F *fMCinputXiStarxiy = new TH2F("fMCinputXiStarxiy","y distribution",80,-2,2, 80,-2,2);
620 TH2F *fMCinputXiStarpiony = new TH2F("fMCinputXiStarpiony","y distribution",80,-2,2, 80,-2,2);
621 fOutputList->Add(fMCinputXiStarxiy);
622 fOutputList->Add(fMCinputXiStarpiony);
623 TH2F *fMCinputXilambday = new TH2F("fMCinputXilambday","y distribution",80,-2,2, 80,-2,2);
624 TH2F *fMCinputXipiony = new TH2F("fMCinputXipiony","y distribution",80,-2,2, 80,-2,2);
625 fOutputList->Add(fMCinputXilambday);
626 fOutputList->Add(fMCinputXipiony);
627 TH2F *fMCinputLamprotony = new TH2F("fMCinputLamprotony","y distribution",80,-2,2, 80,-2,2);
628 TH2F *fMCinputLampiony = new TH2F("fMCinputLampiony","y distribution",80,-2,2, 80,-2,2);
629 fOutputList->Add(fMCinputLamprotony);
630 fOutputList->Add(fMCinputLampiony);
fe45e21b 631
3311ad64 632
633 ///////////////////////////////////
634
635 PostData(1, fOutputList);
636}
637
638//________________________________________________________________________
639void AliXiStar::Exec(Option_t *)
640{
641 // Main loop
642 // Called for each event
643 cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
644 fEventCounter++;
645
646 if(fAODcase) {cout<<"AODs not fully supported! Exiting event."<<endl; return;}
647
648 if(fAODcase) fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
649 else fESD = dynamic_cast<AliESDEvent*> (InputEvent());
650
651 if(fAODcase) {if (!fAOD) {Printf("ERROR: fAOD not available"); return;}}
652 else {if (!fESD) {Printf("ERROR: fESD not available"); return;}}
653
654
655 // ESD Trigger Cut
656 if(!fAODcase){
657 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected())) {cout<<"Event Rejected"<<endl; return;}
658 }
659
660 ///////////////////////////////////////////////////////////
661 const AliAODVertex *PrimaryVertexAOD;
662 const AliESDVertex *PrimaryVertexESD;
663
664 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
665 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
666 fPIDResponse = inputHandler->GetPIDResponse();
667
668
669 TClonesArray *mcArray = 0x0;
670 //AliAODMCParticle *mcXi;
671 //AliAODMCParticle *mcXiStarD2;
672 //AliAODMCParticle *mcXiStar;
673 AliMCEvent *mcEvent = 0x0;
674 AliStack *mcstack = 0x0;
675 TParticle *MCLamD1esd = 0x0;
676 TParticle *MCLamD2esd = 0x0;
677 TParticle *MCLamesd = 0x0;
678 TParticle *MCXiD2esd = 0x0;
679 TParticle *MCXiesd = 0x0;
680 TParticle *MCXiStarD2esd = 0x0;
681 TParticle *MCXiStaresd = 0x0;
682
683 Double_t px1,py1,pz1, px2,py2,pz2;
fe45e21b 684 Double_t p1sq,p2sq,e1,e2,angle;
3311ad64 685 Double_t dca3d;
686 Float_t dca2[2];
fe45e21b 687 Double_t xiVtx[3];//, xiStarVtx[3];
3311ad64 688 Double_t xiP[3], xiStarP[3];
689 Double_t xiStarMom;
690 Double_t xiMass, xiStarMass;
691 Double_t xiPt, xiStarPt;
692 Double_t xiY, xiStarY;
693 Double_t xiCharge;
fe45e21b 694 Double_t decayLengthXY;
3311ad64 695 Double_t pDaughter1[3];
696 Double_t pDaughter2[3];
697 Double_t xDaughter1[3];
698 Double_t xDaughter2[3];
699 //
700 Double_t bField=0;
701 UInt_t status=0;
702 Int_t positiveTracks=0, negativeTracks=0;
703 Int_t myTracks=0;
704 //
705 Double_t primaryVtx[3]={0};
706 Int_t mBin=0;
707 Int_t zBin=0;
708 Double_t zStep=2*10/Double_t(fZvertexBins), zStart=-10.;
709 //
710 Bool_t mcXiFilled=kFALSE;// So that mctracks are never used uninitialized
711
712 if(fMCcase){
713 if(fAODcase){
714 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
715 if(!mcArray){
716 cout<<"No MC particle branch found"<<endl;
717 return;
718 }
719 }else {
720 mcEvent = MCEvent();
721 if (!mcEvent) {Printf("ERROR: Could not retrieve MC event"); return;}
722
723 mcstack = mcEvent->Stack();
724 if (!mcstack) {Printf("ERROR: Could not retrieve the stack"); return;}
725 }
726 }
727
728
729 /////////////////////////////////////////////////
730
731
732
733 if(fAODcase){// AOD case
734 ////////////////////////////////
735 // Vertexing
736 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
737 PrimaryVertexAOD = fAOD->GetPrimaryVertex();
738 if(!PrimaryVertexAOD) return;
739
740 if(fMCcase){
741 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
742 AliAODMCParticle *mcInputTrack = (AliAODMCParticle*)mcArray->At(it);
743 if (!mcInputTrack) {
744 Error("UserExec", "Could not receive track %d", it);
745 continue;
746 }
747
748 if(!mcInputTrack->IsPhysicalPrimary()) continue;
749
750 // Xi
751 if(mcInputTrack->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
752 if(mcInputTrack->GetPdgCode() == -kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
753
754 // XiStar
755 if(mcInputTrack->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
756 if(mcInputTrack->GetPdgCode() == -kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
757 }
758 }
759
760
761
762 primaryVtx[0]=PrimaryVertexAOD->GetX(); primaryVtx[1]=PrimaryVertexAOD->GetY(); primaryVtx[2]=PrimaryVertexAOD->GetZ();
763 ((TH3F*)fOutputList->FindObject("fVertexDist1"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
764 if(fabs(primaryVtx[2]) > 10.) return; // Z-Vertex Cut
765 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
766
767
768 if(fAOD->IsPileupFromSPD()) return; // Reject Pile-up events
769
770 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fAOD->GetNumberOfTracks());
771 ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
772
773 if(PrimaryVertexAOD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fAOD->GetNumberOfTracks());
774
775
776
777 Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
778 // fNtracks Cut
779 if(fAOD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
780 bField = fAOD->GetMagneticField();
781
782 // Track loop
783 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
f15c1f69 784 AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(i));
785 if(!aodtrack) AliFatal("Not a standard AOD");
3311ad64 786 if (!aodtrack) continue;
787
788 status=aodtrack->GetStatus();
789
790
791 if( (status&AliESDtrack::kTPCrefit)==0) continue;// Require tpcrefit
792 if( aodtrack->GetTPCNcls() < 70) continue;// TPC Ncluster cut
793 if(aodtrack->Pt() < 0.15) continue;
794 AliAODVertex *VtxType=aodtrack->GetProdVertex();
795 if((VtxType->GetType())==AliAODVertex::kKink) continue;// Reject Kinks
796
797 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
798 if(!goodMomentum) continue;
799 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
800
801
802 aodtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
803
804 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
805 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
806 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
807
808 ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fAOD->GetNumberOfTracks(), dca3d);
809 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
810 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
811
812
813
814
815 fTempStruct[myTracks].fStatus = status;
816 fTempStruct[myTracks].fFilterMap = aodtrack->GetFilterMap();
817 fTempStruct[myTracks].fID = aodtrack->GetID();
818 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
819 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
820 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
821 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
822 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
823 fTempStruct[myTracks].fEta = aodtrack->Eta();
824 fTempStruct[myTracks].fCharge = aodtrack->Charge();
825 fTempStruct[myTracks].fDCAXY = dca2[0];
826 fTempStruct[myTracks].fDCAZ = dca2[1];
827 fTempStruct[myTracks].fDCA = dca3d;
828 fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
829 fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
830 fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
831
832
833 if(aodtrack->Charge() > 0) positiveTracks++;
834 else negativeTracks++;
835
836
837 myTracks++;
838 }
839 }else {// ESDs
840
841 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fESD->GetNumberOfTracks());
842 PrimaryVertexESD = fESD->GetPrimaryVertex();
843 if(!PrimaryVertexESD) return;
844
845 primaryVtx[0]=PrimaryVertexESD->GetX(); primaryVtx[1]=PrimaryVertexESD->GetY(); primaryVtx[2]=PrimaryVertexESD->GetZ();
846 ((TH3F*)fOutputList->FindObject("fVertexDist1"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
847
848 if(fMCcase){
849 /////////////////////////////////////////////////
850 // Lam mc input
851 /////////////////////////////////////////////////
852 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
853 TParticle *mcInputTrack = (TParticle*)mcstack->Particle(it);
854 if (!mcInputTrack) {
855 Error("UserExec", "Could not receive track %d", it);
856 continue;
857 }
858
859
860 // Xi
861 if(mcInputTrack->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
862 if(mcInputTrack->GetPdgCode() == -kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
863
864 // XiStar
865 if(mcInputTrack->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
866 if(mcInputTrack->GetPdgCode() == -kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
867
868 }
869
870
871 }
872
873
874
875 if(fabs(primaryVtx[2]) > 10.) return; // Z-Vertex Cut
876 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fESD->GetNumberOfTracks());
877
878
879 if(fESD->IsPileupFromSPD()) return; // Reject Pile-up events
880
881 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fESD->GetNumberOfTracks());
882 ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
883
884 if(PrimaryVertexESD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fESD->GetNumberOfTracks());
885
886 Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
887 // fNtracks Cut
888 if(fESD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
889 bField = fESD->GetMagneticField();
890
891
892 // Track loop
893 for (Int_t i = 0; i < fESD->GetNumberOfTracks(); i++) {
894 AliESDtrack* esdtrack = fESD->GetTrack(i);
895 if (!esdtrack) continue;
896
897 status=esdtrack->GetStatus();
898
899 if(!fTrackCut->AcceptTrack(esdtrack)) continue;
900
901 Bool_t goodMomentum = esdtrack->GetPxPyPz( fTempStruct[myTracks].fP);
902 if(!goodMomentum) continue;
903 esdtrack->GetXYZ( fTempStruct[myTracks].fX);
904
905
906 esdtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
907 //esdtrack->GetImpactParameters(dca2, cov);
908 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
909 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
910 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
911
912 ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fESD->GetNumberOfTracks(), dca3d);
913 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(esdtrack->Charge(), esdtrack->Phi(), esdtrack->Pt());
914 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(esdtrack->Charge(), esdtrack->Pt(), esdtrack->Eta());
915
916
917
918 fTempStruct[myTracks].fStatus = status;
919 fTempStruct[myTracks].fID = esdtrack->GetID();
920 fTempStruct[myTracks].fLabel = esdtrack->GetLabel();
921 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
922 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
923 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
924 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
925 fTempStruct[myTracks].fEta = esdtrack->Eta();
926 fTempStruct[myTracks].fCharge = esdtrack->Charge();
927 fTempStruct[myTracks].fDCAXY = dca2[0];
928 fTempStruct[myTracks].fDCAZ = dca2[1];
929 fTempStruct[myTracks].fDCA = dca3d;
930 fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kPion));
931 fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kKaon));
932 fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kProton));
fe45e21b 933 fTempStruct[myTracks].fNclusTPC = esdtrack->GetTPCNcls();
934
935
3311ad64 936 if(esdtrack->Charge() > 0) positiveTracks++;
937 else negativeTracks++;
938
939 myTracks++;
940 }
941
942 }// end of ESD case
943
944
945
946 if(myTracks >= 1) {
947 ((TH1F*)fOutputList->FindObject("fMultDist5"))->Fill(myTracks);
948 ((TH3F*)fOutputList->FindObject("fMultDist3d"))->Fill(positiveTracks+negativeTracks, positiveTracks, negativeTracks);
949 }
950
951
952 cout<<"There are "<<myTracks<<" myTracks"<<endl;
953
954 // set Z Vertex bin
955 for(Int_t i=0; i<fZvertexBins; i++){
956 if( (primaryVtx[2] > zStart+i*zStep) && (primaryVtx[2] < zStart+(i+1)*zStep) ){
957 zBin=i;
958 break;
959 }
960 }
961
962 // set Multiplicity bin
963 for(Int_t i=0; i<fMultBins; i++){
964 if( ( myTracks > fMultLimits[i]) && ( myTracks <= fMultLimits[i+1]) ) { mBin=i; break;}
965 }
966
967
968
969 ////////////////////////////////////
970 // Add event to buffer if > 0 tracks
971 if(myTracks > 0){
972 fEC[zBin][mBin]->FIFOShift();
973 (fEvt) = fEC[zBin][mBin]->fEvtStr;
974 (fEvt)->fNTracks = myTracks;
975 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
976 }
977
978
979
980
981 if(fMCcase && fAODcase){// get Input MC information for AOD case
982
983 /////////////////////////////////////////////////
984 // Xi mc input
985 /////////////////////////////////////////////////
986 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
987 AliAODMCParticle *mcInputTrackXi = (AliAODMCParticle*)mcArray->At(it);
988
989 if (!mcInputTrackXi) {
990 Error("UserExec", "Could not receive track %d", it);
991 continue;
992 }
993
994 if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
995 if(!mcInputTrackXi->IsPhysicalPrimary()) continue;
996
997
998 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
999 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1000
1001
1002
1003 AliAODMCParticle *mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(0)));
1004 AliAODMCParticle *mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(1)));
1005
1006 if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1007 if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1008
1009
1010 // At this point we have the right Xi decay channel
1011
1012 AliAODMCParticle *mcInputTrackLamD1;
1013 AliAODMCParticle *mcInputTrackLamD2;
1014
1015 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1016 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1017 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1018 }
1019 else{
1020 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1021 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1022 }
1023
1024
1025 if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1026 if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1027
1028 // At this point we have the right Lambda decay channel
1029
1030 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputXi"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1031 else ((TH3F*)fOutputList->FindObject("fMCinputXibar"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1032
1033 }
1034
1035
1036
1037 /////////////////////////////////////////////////
1038 // XiStar mc input
1039 /////////////////////////////////////////////////
1040 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
1041 AliAODMCParticle *mcInputTrackXiStar = (AliAODMCParticle*)mcArray->At(it);
1042 if (!mcInputTrackXiStar) {
1043 Error("UserExec", "Could not receive track %d", it);
1044 continue;
1045 }
1046
1047 if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1048 if(!mcInputTrackXiStar->IsPhysicalPrimary()) continue;
1049
1050 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1051 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1052
1053
1054 AliAODMCParticle *mcInputTrackXiStarD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(0)));
1055 AliAODMCParticle *mcInputTrackXiStarD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(1)));
1056
1057 if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kXiCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kXiCode) continue;
1058 if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kPionCode) continue;
1059
1060
1061 // At this point we have the right Xi(1530) decay channel
1062
1063 AliAODMCParticle *mcInputTrackXiD1;
1064 AliAODMCParticle *mcInputTrackXiD2;
1065 if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1066 mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(0)));
1067 mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(1)));
1068 }
1069 else{
1070 mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(0)));
1071 mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(1)));
1072 }
1073
1074 if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1075 if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1076
1077
1078 // At this point we have the right Xi decay channel
1079
1080 AliAODMCParticle *mcInputTrackLamD1;
1081 AliAODMCParticle *mcInputTrackLamD2;
1082
1083 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1084 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1085 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1086 }
1087 else{
1088 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1089 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1090 }
1091
1092 if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1093 if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1094
1095 // At this point we the right Lambda decay channel
1096
1097 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputXiStar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1098 else ((TH3F*)fOutputList->FindObject("fMCinputXiStarbar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1099
1100 if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1101 ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1102 ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1103 }else{
1104 ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1105 ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1106 }
1107 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode){
1108 ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1109 ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1110 }else{
1111 ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1112 ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1113 }
1114 if(abs(mcInputTrackLamD1->GetPdgCode())==kProtonCode){
1115 ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1116 ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1117 }else {
1118 ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1119 ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1120 }
1121
1122
1123 }
1124 }
1125
1126
1127
1128 if(fMCcase && !fAODcase){// get Input MC information for ESD case
1129
1130 /////////////////////////////////////////////////
1131 // Xi mc input
1132 /////////////////////////////////////////////////
1133 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1134 TParticle *mcInputTrackXi = (TParticle*)mcstack->Particle(it);
1135 if (!mcInputTrackXi) {
1136 Error("UserExec", "Could not receive track %d", it);
1137 continue;
1138 }
1139
1140 //if(!mcstack->IsPhysicalPrimary(it)) continue;
1141 if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
1142
1143
1144 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1145 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1146
1147 }
1148
1149
1150 /////////////////////////////////////////////////
1151 // XiStar mc input
1152 /////////////////////////////////////////////////
1153 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1154 TParticle *mcInputTrackXiStar = (TParticle*)mcstack->Particle(it);
1155 if (!mcInputTrackXiStar) {
1156 Error("UserExec", "Could not receive track %d", it);
1157 continue;
1158 }
1159
1160 //if(!mcstack->IsPhysicalPrimary(it)) continue;
1161 if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1162
1163 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1164 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1165
1166 }
1167 }
1168
1169
1170
1171 if(fAODcase) {cout<<"AOD XiVertexer not currently supported! Exiting event"<<endl; return;}
1172
1173 ////////////////////////////////////////////////
1174 // Reconstruction
fe45e21b 1175
3311ad64 1176 for(Int_t i=0; i<fESD->GetNumberOfCascades(); i++){
1177
1178 AliESDcascade *Xicandidate = fESD->GetCascade(i);
1179
1180 if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetNindex())) continue;
1181 if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1182 if(TMath::Abs( Xicandidate->GetNindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1183
1184 AliESDtrack *pTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetPindex()));
1185 AliESDtrack *nTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetNindex()));
1186 AliESDtrack *bTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetBindex()));
1187
1188 // Standard track QA cuts
1189 if(!fTrackCut->AcceptTrack(pTrackXi)) continue;
1190 if(!fTrackCut->AcceptTrack(nTrackXi)) continue;
1191 if(!fTrackCut->AcceptTrack(bTrackXi)) continue;
fe45e21b 1192
1193 //////////////////////
1194 // DecayParameters Key (number represents array index)
1195 // NclustersTPC: 0=proton, 1=pion first, 2=pion second, 3=pion third
1196 // DCAVtx: 4=proton, 5=pion first, 6=pion second, 7=lambda, 8=pion third
1197 // 9 = DCA proton-pion
1198 // 10 = DCA Lambda-pion
1199 // 11 = Rxy Lambda
1200 // 12 = Rxy Xi
1201 // 13 = Cos PA Lambda
1202 // 14 = Cos PA Xi
1203
1204
1205 fDecayParameters[2] = bTrackXi->GetTPCNcls();
3311ad64 1206
fe45e21b 1207
3311ad64 1208 if(Xicandidate->Charge() == -1){
fe45e21b 1209 fDecayParameters[0] = pTrackXi->GetTPCNcls();
1210 fDecayParameters[1] = nTrackXi->GetTPCNcls();
1211 fDecayParameters[4] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1212 fDecayParameters[5] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
3311ad64 1213 }else{
fe45e21b 1214 fDecayParameters[0] = nTrackXi->GetTPCNcls();
1215 fDecayParameters[1] = pTrackXi->GetTPCNcls();
1216 fDecayParameters[4] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1217 fDecayParameters[5] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
3311ad64 1218 }
1219
1220
fe45e21b 1221 fDecayParameters[6] = fabs(bTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion second
1222 fDecayParameters[7] = fabs(Xicandidate->GetD(primaryVtx[0],primaryVtx[1],primaryVtx[2]));// DCA Vtx Lambda
1223 fDecayParameters[9] = fabs(Xicandidate->GetDcaV0Daughters());// DCA proton-pion
1224 fDecayParameters[10] = fabs(Xicandidate->GetDcaXiDaughters());// DCA Lambda-pion
1225
1226
3311ad64 1227
1228 Double_t tempX[3]={0};
1229 Xicandidate->GetXYZ(tempX[0], tempX[1], tempX[2]);
fe45e21b 1230 fDecayParameters[11] = sqrt( pow(tempX[0],2) + pow(tempX[1],2));// Rxy Lambda
3311ad64 1231 if(sqrt( pow(tempX[0],2) + pow(tempX[1],2) ) > fMaxDecayLength) continue;
1232
1233
fe45e21b 1234 fDecayParameters[13] = Xicandidate->GetV0CosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Lambda
1235 fDecayParameters[14] = Xicandidate->GetCascadeCosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Xi
3311ad64 1236
3311ad64 1237 xiP[0] = Xicandidate->Px();
1238 xiP[1] = Xicandidate->Py();
1239 xiP[2] = Xicandidate->Pz();
1240 xiVtx[0] = Xicandidate->Xv();
1241 xiVtx[1] = Xicandidate->Yv();
1242 xiVtx[2] = Xicandidate->Zv();
1243 xiPt = Xicandidate->Pt();
1244 xiY = Xicandidate->RapXi();
1245 xiMass = Xicandidate->M();
1246 xiCharge = Xicandidate->Charge();
1247
1248 decayLengthXY = sqrt( pow(xiVtx[0]-primaryVtx[0],2) + pow(xiVtx[1]-primaryVtx[1],2) );
fe45e21b 1249 fDecayParameters[12] = decayLengthXY;// Rxy Xi
3311ad64 1250 if(decayLengthXY > fMaxDecayLength) continue;// 2d version
1251
fe45e21b 1252 Bool_t StandardXi=kTRUE;
1253 if(fDecayParameters[0] < fCutValues[0][0]) StandardXi=kFALSE;// Nclus proton
1254 if(fDecayParameters[1] < fCutValues[0][1]) StandardXi=kFALSE;// Nclus pion first
1255 if(fDecayParameters[2] < fCutValues[0][2]) StandardXi=kFALSE;// Nclus pion second
1256 //
1257 if(fDecayParameters[4] < fCutValues[0][4]) StandardXi=kFALSE;// DCAVtx proton
1258 if(fDecayParameters[5] < fCutValues[0][5]) StandardXi=kFALSE;// DCAVtx pion first
1259 if(fDecayParameters[6] < fCutValues[0][6]) StandardXi=kFALSE;// DCAVtx pion second
1260 if(fDecayParameters[7] < fCutValues[0][7]) StandardXi=kFALSE;// DCAVtx Lambda
1261 //
1262 if(fDecayParameters[9] > fCutValues[0][9]) StandardXi=kFALSE;// DCAV proton-pion
1263 if(fDecayParameters[10] > fCutValues[0][10]) StandardXi=kFALSE;// DCAV Lambda-pion
1264 //
1265 if(fDecayParameters[11] < fCutValues[0][11]) StandardXi=kFALSE;// Rxy Lambda
1266 if(fDecayParameters[12] < fCutValues[0][12]) StandardXi=kFALSE;// Rxy Xi
1267 //
1268 if(fDecayParameters[13] < fCutValues[0][13]) StandardXi=kFALSE;// Cos PA Lambda
1269 if(fDecayParameters[14] < fCutValues[0][14]) StandardXi=kFALSE;// Cos PA Xi
3311ad64 1270
fe45e21b 1271 if(StandardXi){
1272 if(xiCharge == -1) CutVar[0].fXi->Fill(xiPt, xiY, xiMass);
1273 else CutVar[0].fXibar->Fill(xiPt, xiY, xiMass);
1274 }
3311ad64 1275
1276 // MC associaton
1277 mcXiFilled = kFALSE;
1278 if(fMCcase && !fAODcase){
1279
1280 MCXiD2esd = (TParticle*)mcstack->Particle(abs(bTrackXi->GetLabel()));
1281
1282 if(abs(MCXiD2esd->GetPdgCode())==kPionCode){
1283
1284 MCLamD1esd = (TParticle*)mcstack->Particle(abs(pTrackXi->GetLabel()));
1285 MCLamD2esd = (TParticle*)mcstack->Particle(abs(nTrackXi->GetLabel()));
1286
1287 if(MCLamD1esd->GetMother(0) == MCLamD2esd->GetMother(0)){
1288 if(abs(MCLamD1esd->GetPdgCode())==kProtonCode || abs(MCLamD2esd->GetPdgCode())==kProtonCode) {
1289 if(abs(MCLamD1esd->GetPdgCode())==kPionCode || abs(MCLamD2esd->GetPdgCode())==kPionCode) {
1290
1291 MCLamesd = (TParticle*)mcstack->Particle(abs(MCLamD1esd->GetMother(0)));
1292 if(abs(MCLamesd->GetPdgCode())==kLambdaCode) {
1293
1294 if(MCLamesd->GetMother(0) == MCXiD2esd->GetMother(0)){
1295 MCXiesd = (TParticle*)mcstack->Particle(abs(MCLamesd->GetMother(0)));
1296 if(abs(MCXiesd->GetPdgCode())==kXiCode) {
3311ad64 1297 mcXiFilled = kTRUE;
fe45e21b 1298
1299 if(StandardXi){
1300 if(Xicandidate->Charge() == -1) {
1301 CutVar[0].fMCrecXi->Fill(xiPt, xiY, xiMass);
1302 }else {
1303 CutVar[0].fMCrecXibar->Fill(xiPt, xiY, xiMass);
1304 }
1305 }
1306
3311ad64 1307 }
1308 }
1309 }
1310 }
1311 }
1312 }
1313 }
1314 }// MC association
1315
1316
1317 if(fabs(xiMass-fTrueMassXi) > fMassWindow) continue;
1318
1319
1320
1321 fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
1322
1323
1324 //////////////////////////////////////////////////////////
1325 // Reconstruct Xi(1530)
1326 for(Int_t EN=0; EN<fEventsToMix+1; EN++){// Event buffer loop
1327
1328 for(Int_t l=0; l<(fEvt+EN)->fNTracks; l++){// Present(EN=0) and Past(EN from 1 to fEventsToMix) event track loop
1329
1330 if(EN==0) {
1331 if((fEvt+EN)->fTracks[l].fID == pTrackXi->GetID()) continue;
1332 if((fEvt+EN)->fTracks[l].fID == nTrackXi->GetID()) continue;
1333 if((fEvt+EN)->fTracks[l].fID == bTrackXi->GetID()) continue;
1334 }
1335
1336 fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
fe45e21b 1337
3311ad64 1338 if(!fESDTrack4) continue;
92650a3e 1339 fESDTrack4->Set((fEvt+EN)->fTracks[l].fX, (fEvt+EN)->fTracks[l].fP, (fEvt+EN)->fTracks[l].fCov, (fEvt+EN)->fTracks[l].fCharge);
3311ad64 1340 if(fAODcase){
1341 if((Bool_t)(((1<<5) & (fEvt+EN)->fTracks[l].fFilterMap) == 0)) continue;// AOD filterbit cut, "Standard cuts with tight dca"
1342 }else{
fe45e21b 1343 fDecayParameters[8] = (fEvt+EN)->fTracks[l].fDCAXY;// DCA Vtx pion third
3311ad64 1344 if((fEvt+EN)->fTracks[l].fDCAZ > 2) continue;
1345 if( (((fEvt+EN)->fTracks[l].fStatus)&AliESDtrack::kITSrefit)==0) continue;// Require itsrefit
1346 // no Chi^2 cut applied for ESDs. Info not available in my track structure.
1347 }
1348
3311ad64 1349 if(fabs((fEvt+EN)->fTracks[l].fEta) > 0.8) continue;
1350
fe45e21b 1351 fDecayParameters[3] = (fEvt+EN)->fTracks[l].fNclusTPC;
3311ad64 1352
1353 AliVertex *XiStarVtx = new AliVertex((fEvt+EN)->fTracks[l].fX,0,0);
1354 //fESDTrack4->PropagateToDCA(fXiTrack, bField);// Propagate tracks to dca, both tracks are budged
1355 if(!(fXiTrack->PropagateToDCA(XiStarVtx, bField, 3))) continue;// Propagate tracks to dca, version which assumes fESDTrack4 is already primary
1356 /////////////
1357 fXiTrack->GetPxPyPz(pDaughter1);
1358 fXiTrack->GetXYZ(xDaughter1);
1359 fESDTrack4->GetPxPyPz(pDaughter2);
1360 fESDTrack4->GetXYZ(xDaughter2);
1361 //////////////////////////
1362
1363
1364
fe45e21b 1365 //xiStarVtx[0] = (xDaughter1[0]+xDaughter2[0])/2.;
1366 //xiStarVtx[1] = (xDaughter1[1]+xDaughter2[1])/2.;
1367 //xiStarVtx[2] = (xDaughter1[2]+xDaughter2[2])/2.;
1368 //decayLength = sqrt(pow(xiStarVtx[0]-primaryVtx[0],2)+pow(xiStarVtx[1]-primaryVtx[1],2)+pow(xiStarVtx[2]-primaryVtx[2],2));
3311ad64 1369
1370 px1=pDaughter1[0];
1371 py1=pDaughter1[1];
1372 pz1=pDaughter1[2];
1373 px2=pDaughter2[0];
1374 py2=pDaughter2[1];
1375 pz2=pDaughter2[2];
1376
1377 p1sq=px1*px1+py1*py1+pz1*pz1;
1378 p2sq=px2*px2+py2*py2+pz2*pz2;
1379 if(p1sq <=0 || p2sq <=0) continue;
1380
1381 e1=sqrt(p1sq+fTrueMassXi*fTrueMassXi);
1382 e2=sqrt(p2sq+fTrueMassPi*fTrueMassPi);
1383 angle=px1*px2+py1*py2+pz1*pz2;
1384 xiStarMass=fTrueMassXi*fTrueMassXi+fTrueMassPi*fTrueMassPi+2.*e1*e2-2.*angle;
1385 if(xiStarMass<0.) xiStarMass=1.e-8;
1386 xiStarMass=sqrt(xiStarMass);
1387
1388
1389 xiStarP[0] = px1+px2;
1390 xiStarP[1] = py1+py2;
1391 xiStarP[2] = pz1+pz2;
1392 xiStarMom = sqrt(pow(xiStarP[0],2)+pow(xiStarP[1],2)+pow(xiStarP[2],2));
1393 if(xiStarMom==0) continue; // So one of the following lines doesnt break
1394 xiStarPt = sqrt(xiStarP[0]*xiStarP[0] + xiStarP[1]*xiStarP[1]);
1395 xiStarY = .5*log( ((e1+e2) + xiStarP[2])/((e1+e2) - xiStarP[2]));
fe45e21b 1396 //xiStarE = e1 + e2;
3311ad64 1397
1398
fe45e21b 1399 //if( (xiStarP[0]*(xiStarVtx[0]-primaryVtx[0]) + xiStarP[1]*(xiStarVtx[1]-primaryVtx[1]) + xiStarP[2]*(xiStarVtx[2]-primaryVtx[2]))/xiStarMom/decayLength < fXiStarCosTheta) continue;
1400
1401 for(int cv=0; cv<kNCutVariations; cv++){
1402 if(fDecayParameters[0] < fCutValues[cv][0]) continue;// Nclus proton
1403 if(fDecayParameters[1] < fCutValues[cv][1]) continue;// Nclus pion first
1404 if(fDecayParameters[2] < fCutValues[cv][2]) continue;// Nclus pion second
1405 if(fDecayParameters[3] < fCutValues[cv][3]) continue;// Nclus pion third
1406 //
1407 if(fDecayParameters[4] < fCutValues[cv][4]) continue;// DCAVtx proton
1408 if(fDecayParameters[5] < fCutValues[cv][5]) continue;// DCAVtx pion first
1409 if(fDecayParameters[6] < fCutValues[cv][6]) continue;// DCAVtx pion second
1410 if(fDecayParameters[7] < fCutValues[cv][7]) continue;// DCAVtx Lambda
1411 if(cv!=8) {if(fDecayParameters[8] > (0.0182 + 0.035/pow((fEvt+EN)->fTracks[l].fPt,1.01))) continue;}// DCAVtx pion third
1412 else {if(fDecayParameters[8] > fCutValues[cv][8]) continue;}// DCAVtx pion third
1413 //
1414 if(fDecayParameters[9] > fCutValues[cv][9]) continue;// DCAV proton-pion
1415 if(fDecayParameters[10] > fCutValues[cv][10]) continue;// DCAV Lambda-pion
1416 //
1417 if(fDecayParameters[11] < fCutValues[cv][11]) continue;// Rxy Lambda
1418 if(fDecayParameters[12] < fCutValues[cv][12]) continue;// Rxy Xi
1419 //
1420 if(fDecayParameters[13] < fCutValues[cv][13]) continue;// Cos PA Lambda
1421 if(fDecayParameters[14] < fCutValues[cv][14]) continue;// Cos PA Xi
1422
3311ad64 1423
fe45e21b 1424 if(EN==0){
1425 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1426 else if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1427 else if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1428 else CutVar[cv].fXiPlusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1429 }else {
1430 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1431 else if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1432 else if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1433 else CutVar[cv].fXiPlusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1434 }
1435
3311ad64 1436
1437 /*
fe45e21b 1438 // MC associaton AOD
3311ad64 1439 if(fMCcase && mcXiFilled && EN==0 && fAODcase){// AOD MC's
1440
1441 MCXiStarD2 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[l].fLabel));
1442
1443 if(abs(MCXiStarD2->GetPdgCode())==kPionCode){
1444 if(MCXi->GetMother() == MCXiStarD2->GetMother()){
1445 MCXiStar = (AliAODMCParticle*)mcArray->At(MCXi->GetMother());
1446 if(abs(MCXiStar->GetPdgCode())==kXiStarCode) {
1447
1448 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1449 if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1450
1451 }
1452 }
1453 }
1454 }
1455 */
1456
fe45e21b 1457 // MC associaton ESD
1458 if(fMCcase && mcXiFilled && EN==0 && !fAODcase){// ESD MC's
1459 MCXiStarD2esd = (TParticle*)mcstack->Particle(abs((fEvt)->fTracks[l].fLabel));
1460
1461 if(abs(MCXiStarD2esd->GetPdgCode())==kPionCode){
1462 if(MCXiesd->GetMother(0) == MCXiStarD2esd->GetMother(0)){
1463
1464 MCXiStaresd = (TParticle*)mcstack->Particle(abs(MCXiesd->GetMother(0)));
1465 if(abs(MCXiStaresd->GetPdgCode())==kXiStarCode) {
1466
1467 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1468 if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1469
1470 }
3311ad64 1471 }
1472 }
1473 }
3311ad64 1474
fe45e21b 1475 }// Cut Variation loop
3311ad64 1476 }// 3rd pion loop
1477 }// Event mixing loop
1478
1479 }// Xi loop
1480
1481
1482
1483 // Post output data.
1484 PostData(1, fOutputList);
1485
1486}
1487//________________________________________________________________________
1488void AliXiStar::Terminate(Option_t *)
1489{
1490 cout<<"Done"<<endl;
1491}
1492//________________________________________________________________________
1493Double_t AliXiStar::LinearPropagateToDCA(AliESDtrack *v, AliESDtrack *t, Double_t b) {// Adapted from AliCascadeVertexer.cxx
1494 //--------------------------------------------------------------------
1495 // This function returns the DCA between the V0 and the track
1496 //--------------------------------------------------------------------
1497
1498 Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
1499 Double_t r[3]; t->GetXYZ(r);
1500 Double_t x1=r[0], y1=r[1], z1=r[2];
1501 Double_t p[3]; t->GetPxPyPz(p);
1502 Double_t px1=p[0], py1=p[1], pz1=p[2];
1503
1504 Double_t x2[3]={0};
1505 Double_t p2[3]={0};
1506 Double_t vx2,vy2,vz2; // position and momentum of V0
1507 Double_t px2,py2,pz2;
1508
1509 v->GetXYZ(x2);
1510 v->GetPxPyPz(p2);
1511 vx2=x2[0], vy2=x2[1], vz2=x2[2];
1512 px2=p2[0], py2=p2[1], pz2=p2[2];
1513
1514// calculation dca
1515
1516 Double_t dd= Det(vx2-x1,vy2-y1,vz2-z1,px1,py1,pz1,px2,py2,pz2);
1517 Double_t ax= Det(py1,pz1,py2,pz2);
1518 Double_t ay=-Det(px1,pz1,px2,pz2);
1519 Double_t az= Det(px1,py1,px2,py2);
1520
1521 Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
1522
1523//points of the DCA
1524 Double_t t1 = Det(vx2-x1,vy2-y1,vz2-z1,px2,py2,pz2,ax,ay,az)/
1525 Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
1526
1527 x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1;
1528
1529
1530 //propagate track to the points of DCA
1531
1532 x1=x1*cs1 + y1*sn1;
1533 if (!t->PropagateTo(x1,b)) {
1534 Error("PropagateToDCA","Propagation failed !");
1535 return 1.e+33;
1536 }
1537
1538 return dca;
1539}
1540
1541
1542//________________________________________________________________________
1543Double_t AliXiStar::Det(Double_t a00, Double_t a01, Double_t a10, Double_t a11) const {// Taken from AliCascadeVertexer
1544 //--------------------------------------------------------------------
1545 // This function calculates locally a 2x2 determinant
1546 //--------------------------------------------------------------------
1547 return a00*a11 - a01*a10;
1548}
1549//________________________________________________________________________
1550Double_t AliXiStar::Det(Double_t a00,Double_t a01,Double_t a02,
1551 Double_t a10,Double_t a11,Double_t a12,
1552 Double_t a20,Double_t a21,Double_t a22) const {// Taken from AliCascadeVertexer
1553 //--------------------------------------------------------------------
1554 // This function calculates locally a 3x3 determinant
1555 //--------------------------------------------------------------------
1556 return a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);
1557}
1558
1559
1560