1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ////////////////////////////////////////////////////////////////////////////////
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.
22 // authors: Dhevan Gangadharan (dhevan.raja.gangadharan@cern.ch)
24 ////////////////////////////////////////////////////////////////////////////////
34 #include "TObjString.h"
44 #include "AliAnalysisTask.h"
45 #include "AliAnalysisManager.h"
48 #include "AliESDEvent.h"
49 #include "AliESDInputHandler.h"
50 #include "AliESDtrackCuts.h"
51 #include "AliMCEventHandler.h"
52 #include "AliMCEvent.h"
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"
63 #include "AliXiStar.h"
68 // Author: Dhevan Gangadharan
72 //________________________________________________________________________
73 AliXiStar::AliXiStar():
101 // Default Constructor
102 for (Int_t i=0; i<21; i++){
103 fCovMatrix[i]=-99999.;
104 if (i<12) fMultLimits[i] = 0;
106 for (Int_t i=0; i<kNCuts; i++){
107 fDecayParameters[i]=0;
108 for (Int_t j=0; j<kNCutVariations; j++){
113 for (Int_t cv=0; cv<kNCutVariations; cv++){
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;
121 CutVar[cv].fXiMinusPiPlusbkg=0x0;
122 CutVar[cv].fXiMinusPiMinusbkg=0x0;
123 CutVar[cv].fXiPlusPiPlusbkg=0x0;
124 CutVar[cv].fXiPlusPiMinusbkg=0x0;
126 CutVar[cv].fMCrecXi=0x0;
127 CutVar[cv].fMCrecXibar=0x0;
128 CutVar[cv].fMCrecXiMinusPiPlus=0x0;
129 CutVar[cv].fMCrecXiPlusPiMinus=0x0;
133 //________________________________________________________________________
134 AliXiStar::AliXiStar(const char *name, Bool_t AODdecision, Bool_t MCdecision, Int_t CutListOption)
135 : AliAnalysisTaskSE(name),
149 fAODcase(AODdecision),
160 fCutList(CutListOption)
164 for (Int_t i=0; i<21; i++){
165 fCovMatrix[i]=-99999.;
166 if (i<12) fMultLimits[i] = 0;
168 for (Int_t i=0; i<kNCuts; i++){
169 fDecayParameters[i]=0;
170 for (Int_t j=0; j<kNCutVariations; j++){
175 for (Int_t cv=0; cv<kNCutVariations; cv++){
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;
183 CutVar[cv].fXiMinusPiPlusbkg=0x0;
184 CutVar[cv].fXiMinusPiMinusbkg=0x0;
185 CutVar[cv].fXiPlusPiPlusbkg=0x0;
186 CutVar[cv].fXiPlusPiMinusbkg=0x0;
188 CutVar[cv].fMCrecXi=0x0;
189 CutVar[cv].fMCrecXibar=0x0;
190 CutVar[cv].fMCrecXiMinusPiPlus=0x0;
191 CutVar[cv].fMCrecXiPlusPiMinus=0x0;
195 // Define output slots here
197 DefineOutput(1, TList::Class());
200 //________________________________________________________________________
201 AliXiStar::AliXiStar(const AliXiStar &obj)
202 : AliAnalysisTaskSE(obj.fname),
206 fOutputList(obj.fOutputList),
207 fTrackCut(obj.fTrackCut),
208 fPIDResponse(obj.fPIDResponse),
211 fTempStruct(obj.fTempStruct),
212 fZvertexBins(obj.fZvertexBins),
213 fEventsToMix(obj.fEventsToMix),
214 fMultBins(obj.fMultBins),
216 fMCcase(obj.fMCcase),
217 fAODcase(obj.fAODcase),
218 fEventCounter(obj.fEventCounter),
219 fMaxDecayLength(obj.fMaxDecayLength),
220 fMassWindow(obj.fMassWindow),
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)
232 for (Int_t i=0; i<21; i++){
233 fCovMatrix[i]=obj.fCovMatrix[i];
234 if (i<12) fMultLimits[i]=obj.fMultLimits[i];
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];
244 //________________________________________________________________________
245 AliXiStar &AliXiStar::operator=(const AliXiStar &obj)
247 // Assignment operator
254 fOutputList = obj.fOutputList;
255 fTrackCut = obj.fTrackCut;
256 fPIDResponse = obj.fPIDResponse;
259 fTempStruct = obj.fTempStruct;
260 fZvertexBins = obj.fZvertexBins;
261 fEventsToMix = obj.fEventsToMix;
262 fMultBins = obj.fMultBins;
263 for (Int_t i=0; i<12; i++){
264 fMultLimits[i]=obj.fMultLimits[i];
266 fMCcase = obj.fMCcase;
267 fAODcase = obj.fAODcase;
268 fEventCounter = obj.fEventCounter;
269 fMaxDecayLength = obj.fMaxDecayLength;
270 fMassWindow = obj.fMassWindow;
271 for (Int_t i=0; i<21; i++){
272 fCovMatrix[i]=obj.fCovMatrix[i];
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;
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];
293 //________________________________________________________________________
294 AliXiStar::~AliXiStar()
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;
303 if(fEvt) delete fEvt;
304 if(fTempStruct) delete fTempStruct;
305 if(fESDTrack4) delete fESDTrack4;
306 if(fXiTrack) delete fXiTrack;
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;
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;
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;
328 //________________________________________________________________________
329 void AliXiStar::XiStarInit()
332 //Inits cuts and analysis settings
335 fEventCounter=0;// event counter initialization
336 cout<<"AliXiStar XiStarInit() call"<<endl;
339 ///////////////////////////////////////////////
340 // Track Cuts for ESD analysis
341 fTrackCut = new AliESDtrackCuts();
342 fTrackCut->SetPtRange(.15,1000);
343 fTrackCut->SetAcceptKinkDaughters(kFALSE);
344 //fTrackCut->SetMinNClustersTPC(70);
345 fTrackCut->SetRequireTPCRefit(kTRUE);
346 ////////////////////////////////////////////////
349 fMultBins = 11;// This must also be set in AliXiStar.h
350 if(fMCcase) fEventsToMix = 0;
351 else fEventsToMix = 40;
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;
358 fEC = new AliXiStarEventCollection **[fZvertexBins];
359 for(unsigned short i=0; i<fZvertexBins; i++){
361 fEC[i] = new AliXiStarEventCollection *[fMultBins];
363 for(unsigned short j=0; j<fMultBins; j++){
365 fEC[i][j] = new AliXiStarEventCollection(fEventsToMix+1);
370 fTempStruct = new AliXiStarTrackStruct[kNbinsM];
372 fESDTrack4 = new AliESDtrack();
373 fXiTrack = new AliESDtrack();
376 fMaxDecayLength = 100.;
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
388 // 13 = Cos PA Lambda
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];
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
424 fTrueMassPr=.93827, fTrueMassPi=.13957, fTrueMassK=.493677, fTrueMassLam=1.11568, fTrueMassXi=1.32171;
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;
432 //________________________________________________________________________
433 void AliXiStar::UserCreateOutputObjects()
436 XiStarInit();// Initialize settings
439 fOutputList = new TList();
440 fOutputList->SetOwner();
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);
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);
454 TH2F *fDCADist = new TH2F("fDCADist","DCA distribution",kNbinsM,-.5,kNbinsM-.5, 100,0,10);
455 fOutputList->Add(fDCADist);
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);
466 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
467 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
468 fOutputList->Add(fMultDist1);
470 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
471 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
472 fOutputList->Add(fMultDist2);
474 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
475 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
476 fOutputList->Add(fMultDist3);
478 TH1F *fMultDist4 = new TH1F("fMultDist4","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
479 fMultDist4->GetXaxis()->SetTitle("Multiplicity");
480 fOutputList->Add(fMultDist4);
482 TH1F *fMultDist5 = new TH1F("fMultDist5","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
483 fMultDist5->GetXaxis()->SetTitle("Multiplicity");
484 fOutputList->Add(fMultDist5);
487 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","PtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
488 fOutputList->Add(fPtEtaDist);
490 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","PhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
491 fOutputList->Add(fPhiPtDist);
494 for(Int_t cv=0; cv<kNCutVariations; cv++){
497 TString *nameXi=new TString("fXi_");
498 TString *nameXibar=new TString("fXibar_");
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);
506 TString *nameMCrecXi = new TString("fMCrecXi_");
507 TString *nameMCrecXibar = new TString("fMCrecXi_");
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);
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);
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);
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);
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);
581 //////////////////////
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);
585 fOutputList->Add(fMCinputXiStar);
586 fOutputList->Add(fMCinputXiStarbar);
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);
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);
597 fOutputList->Add(fMCinputTotalXiStar1);
598 fOutputList->Add(fMCinputTotalXiStarbar1);
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);
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);
609 fOutputList->Add(fMCinputTotalXiStar3);
610 fOutputList->Add(fMCinputTotalXiStarbar3);
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);
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);
633 ///////////////////////////////////
635 PostData(1, fOutputList);
638 //________________________________________________________________________
639 void AliXiStar::Exec(Option_t *)
642 // Called for each event
643 cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
646 if(fAODcase) {cout<<"AODs not fully supported! Exiting event."<<endl; return;}
648 if(fAODcase) fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
649 else fESD = dynamic_cast<AliESDEvent*> (InputEvent());
651 if(fAODcase) {if (!fAOD) {Printf("ERROR: fAOD not available"); return;}}
652 else {if (!fESD) {Printf("ERROR: fESD not available"); return;}}
657 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected())) {cout<<"Event Rejected"<<endl; return;}
660 ///////////////////////////////////////////////////////////
661 const AliAODVertex *PrimaryVertexAOD;
662 const AliESDVertex *PrimaryVertexESD;
664 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
665 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
666 fPIDResponse = inputHandler->GetPIDResponse();
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;
683 Double_t px1,py1,pz1, px2,py2,pz2;
684 Double_t p1sq,p2sq,e1,e2,angle;
687 Double_t xiVtx[3];//, xiStarVtx[3];
688 Double_t xiP[3], xiStarP[3];
690 Double_t xiMass, xiStarMass;
691 Double_t xiPt, xiStarPt;
692 Double_t xiY, xiStarY;
694 Double_t decayLengthXY;
695 Double_t pDaughter1[3];
696 Double_t pDaughter2[3];
697 Double_t xDaughter1[3];
698 Double_t xDaughter2[3];
702 Int_t positiveTracks=0, negativeTracks=0;
705 Double_t primaryVtx[3]={0};
708 Double_t zStep=2*10/Double_t(fZvertexBins), zStart=-10.;
710 Bool_t mcXiFilled=kFALSE;// So that mctracks are never used uninitialized
714 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
716 cout<<"No MC particle branch found"<<endl;
721 if (!mcEvent) {Printf("ERROR: Could not retrieve MC event"); return;}
723 mcstack = mcEvent->Stack();
724 if (!mcstack) {Printf("ERROR: Could not retrieve the stack"); return;}
729 /////////////////////////////////////////////////
733 if(fAODcase){// AOD case
734 ////////////////////////////////
736 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
737 PrimaryVertexAOD = fAOD->GetPrimaryVertex();
738 if(!PrimaryVertexAOD) return;
741 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
742 AliAODMCParticle *mcInputTrack = (AliAODMCParticle*)mcArray->At(it);
744 Error("UserExec", "Could not receive track %d", it);
748 if(!mcInputTrack->IsPhysicalPrimary()) continue;
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());
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());
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());
768 if(fAOD->IsPileupFromSPD()) return; // Reject Pile-up events
770 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fAOD->GetNumberOfTracks());
771 ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
773 if(PrimaryVertexAOD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fAOD->GetNumberOfTracks());
777 Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
779 if(fAOD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
780 bField = fAOD->GetMagneticField();
783 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
784 AliAODTrack* aodtrack = fAOD->GetTrack(i);
785 if (!aodtrack) continue;
787 status=aodtrack->GetStatus();
790 if( (status&AliESDtrack::kTPCrefit)==0) continue;// Require tpcrefit
791 if( aodtrack->GetTPCNcls() < 70) continue;// TPC Ncluster cut
792 if(aodtrack->Pt() < 0.15) continue;
793 AliAODVertex *VtxType=aodtrack->GetProdVertex();
794 if((VtxType->GetType())==AliAODVertex::kKink) continue;// Reject Kinks
796 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
797 if(!goodMomentum) continue;
798 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
801 aodtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
803 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
804 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
805 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
807 ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fAOD->GetNumberOfTracks(), dca3d);
808 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
809 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
814 fTempStruct[myTracks].fStatus = status;
815 fTempStruct[myTracks].fFilterMap = aodtrack->GetFilterMap();
816 fTempStruct[myTracks].fID = aodtrack->GetID();
817 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
818 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
819 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
820 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
821 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
822 fTempStruct[myTracks].fEta = aodtrack->Eta();
823 fTempStruct[myTracks].fCharge = aodtrack->Charge();
824 fTempStruct[myTracks].fDCAXY = dca2[0];
825 fTempStruct[myTracks].fDCAZ = dca2[1];
826 fTempStruct[myTracks].fDCA = dca3d;
827 fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
828 fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
829 fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
832 if(aodtrack->Charge() > 0) positiveTracks++;
833 else negativeTracks++;
840 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fESD->GetNumberOfTracks());
841 PrimaryVertexESD = fESD->GetPrimaryVertex();
842 if(!PrimaryVertexESD) return;
844 primaryVtx[0]=PrimaryVertexESD->GetX(); primaryVtx[1]=PrimaryVertexESD->GetY(); primaryVtx[2]=PrimaryVertexESD->GetZ();
845 ((TH3F*)fOutputList->FindObject("fVertexDist1"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
848 /////////////////////////////////////////////////
850 /////////////////////////////////////////////////
851 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
852 TParticle *mcInputTrack = (TParticle*)mcstack->Particle(it);
854 Error("UserExec", "Could not receive track %d", it);
860 if(mcInputTrack->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
861 if(mcInputTrack->GetPdgCode() == -kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
864 if(mcInputTrack->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
865 if(mcInputTrack->GetPdgCode() == -kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
874 if(fabs(primaryVtx[2]) > 10.) return; // Z-Vertex Cut
875 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fESD->GetNumberOfTracks());
878 if(fESD->IsPileupFromSPD()) return; // Reject Pile-up events
880 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fESD->GetNumberOfTracks());
881 ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
883 if(PrimaryVertexESD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fESD->GetNumberOfTracks());
885 Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
887 if(fESD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
888 bField = fESD->GetMagneticField();
892 for (Int_t i = 0; i < fESD->GetNumberOfTracks(); i++) {
893 AliESDtrack* esdtrack = fESD->GetTrack(i);
894 if (!esdtrack) continue;
896 status=esdtrack->GetStatus();
898 if(!fTrackCut->AcceptTrack(esdtrack)) continue;
900 Bool_t goodMomentum = esdtrack->GetPxPyPz( fTempStruct[myTracks].fP);
901 if(!goodMomentum) continue;
902 esdtrack->GetXYZ( fTempStruct[myTracks].fX);
905 esdtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
906 //esdtrack->GetImpactParameters(dca2, cov);
907 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
908 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
909 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
911 ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fESD->GetNumberOfTracks(), dca3d);
912 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(esdtrack->Charge(), esdtrack->Phi(), esdtrack->Pt());
913 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(esdtrack->Charge(), esdtrack->Pt(), esdtrack->Eta());
917 fTempStruct[myTracks].fStatus = status;
918 fTempStruct[myTracks].fID = esdtrack->GetID();
919 fTempStruct[myTracks].fLabel = esdtrack->GetLabel();
920 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
921 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
922 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
923 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
924 fTempStruct[myTracks].fEta = esdtrack->Eta();
925 fTempStruct[myTracks].fCharge = esdtrack->Charge();
926 fTempStruct[myTracks].fDCAXY = dca2[0];
927 fTempStruct[myTracks].fDCAZ = dca2[1];
928 fTempStruct[myTracks].fDCA = dca3d;
929 fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kPion));
930 fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kKaon));
931 fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kProton));
932 fTempStruct[myTracks].fNclusTPC = esdtrack->GetTPCNcls();
935 if(esdtrack->Charge() > 0) positiveTracks++;
936 else negativeTracks++;
946 ((TH1F*)fOutputList->FindObject("fMultDist5"))->Fill(myTracks);
947 ((TH3F*)fOutputList->FindObject("fMultDist3d"))->Fill(positiveTracks+negativeTracks, positiveTracks, negativeTracks);
951 cout<<"There are "<<myTracks<<" myTracks"<<endl;
954 for(Int_t i=0; i<fZvertexBins; i++){
955 if( (primaryVtx[2] > zStart+i*zStep) && (primaryVtx[2] < zStart+(i+1)*zStep) ){
961 // set Multiplicity bin
962 for(Int_t i=0; i<fMultBins; i++){
963 if( ( myTracks > fMultLimits[i]) && ( myTracks <= fMultLimits[i+1]) ) { mBin=i; break;}
968 ////////////////////////////////////
969 // Add event to buffer if > 0 tracks
971 fEC[zBin][mBin]->FIFOShift();
972 (fEvt) = fEC[zBin][mBin]->fEvtStr;
973 (fEvt)->fNTracks = myTracks;
974 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
980 if(fMCcase && fAODcase){// get Input MC information for AOD case
982 /////////////////////////////////////////////////
984 /////////////////////////////////////////////////
985 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
986 AliAODMCParticle *mcInputTrackXi = (AliAODMCParticle*)mcArray->At(it);
988 if (!mcInputTrackXi) {
989 Error("UserExec", "Could not receive track %d", it);
993 if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
994 if(!mcInputTrackXi->IsPhysicalPrimary()) continue;
997 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
998 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1002 AliAODMCParticle *mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(0)));
1003 AliAODMCParticle *mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(1)));
1005 if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1006 if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1009 // At this point we have the right Xi decay channel
1011 AliAODMCParticle *mcInputTrackLamD1;
1012 AliAODMCParticle *mcInputTrackLamD2;
1014 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1015 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1016 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1019 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1020 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1024 if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1025 if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1027 // At this point we have the right Lambda decay channel
1029 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputXi"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1030 else ((TH3F*)fOutputList->FindObject("fMCinputXibar"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1036 /////////////////////////////////////////////////
1038 /////////////////////////////////////////////////
1039 for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
1040 AliAODMCParticle *mcInputTrackXiStar = (AliAODMCParticle*)mcArray->At(it);
1041 if (!mcInputTrackXiStar) {
1042 Error("UserExec", "Could not receive track %d", it);
1046 if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1047 if(!mcInputTrackXiStar->IsPhysicalPrimary()) continue;
1049 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1050 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1053 AliAODMCParticle *mcInputTrackXiStarD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(0)));
1054 AliAODMCParticle *mcInputTrackXiStarD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(1)));
1056 if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kXiCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kXiCode) continue;
1057 if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kPionCode) continue;
1060 // At this point we have the right Xi(1530) decay channel
1062 AliAODMCParticle *mcInputTrackXiD1;
1063 AliAODMCParticle *mcInputTrackXiD2;
1064 if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1065 mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(0)));
1066 mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(1)));
1069 mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(0)));
1070 mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(1)));
1073 if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1074 if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1077 // At this point we have the right Xi decay channel
1079 AliAODMCParticle *mcInputTrackLamD1;
1080 AliAODMCParticle *mcInputTrackLamD2;
1082 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1083 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1084 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1087 mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1088 mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1091 if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1092 if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1094 // At this point we the right Lambda decay channel
1096 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputXiStar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1097 else ((TH3F*)fOutputList->FindObject("fMCinputXiStarbar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1099 if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1100 ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1101 ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1103 ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1104 ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1106 if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode){
1107 ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1108 ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1110 ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1111 ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1113 if(abs(mcInputTrackLamD1->GetPdgCode())==kProtonCode){
1114 ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1115 ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1117 ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1118 ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1127 if(fMCcase && !fAODcase){// get Input MC information for ESD case
1129 /////////////////////////////////////////////////
1131 /////////////////////////////////////////////////
1132 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1133 TParticle *mcInputTrackXi = (TParticle*)mcstack->Particle(it);
1134 if (!mcInputTrackXi) {
1135 Error("UserExec", "Could not receive track %d", it);
1139 //if(!mcstack->IsPhysicalPrimary(it)) continue;
1140 if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
1143 if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1144 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1149 /////////////////////////////////////////////////
1151 /////////////////////////////////////////////////
1152 for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1153 TParticle *mcInputTrackXiStar = (TParticle*)mcstack->Particle(it);
1154 if (!mcInputTrackXiStar) {
1155 Error("UserExec", "Could not receive track %d", it);
1159 //if(!mcstack->IsPhysicalPrimary(it)) continue;
1160 if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1162 if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1163 else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1170 if(fAODcase) {cout<<"AOD XiVertexer not currently supported! Exiting event"<<endl; return;}
1172 ////////////////////////////////////////////////
1175 for(Int_t i=0; i<fESD->GetNumberOfCascades(); i++){
1177 AliESDcascade *Xicandidate = fESD->GetCascade(i);
1179 if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetNindex())) continue;
1180 if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1181 if(TMath::Abs( Xicandidate->GetNindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1183 AliESDtrack *pTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetPindex()));
1184 AliESDtrack *nTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetNindex()));
1185 AliESDtrack *bTrackXi = fESD->GetTrack(TMath::Abs( Xicandidate->GetBindex()));
1187 // Standard track QA cuts
1188 if(!fTrackCut->AcceptTrack(pTrackXi)) continue;
1189 if(!fTrackCut->AcceptTrack(nTrackXi)) continue;
1190 if(!fTrackCut->AcceptTrack(bTrackXi)) continue;
1192 //////////////////////
1193 // DecayParameters Key (number represents array index)
1194 // NclustersTPC: 0=proton, 1=pion first, 2=pion second, 3=pion third
1195 // DCAVtx: 4=proton, 5=pion first, 6=pion second, 7=lambda, 8=pion third
1196 // 9 = DCA proton-pion
1197 // 10 = DCA Lambda-pion
1200 // 13 = Cos PA Lambda
1204 fDecayParameters[2] = bTrackXi->GetTPCNcls();
1207 if(Xicandidate->Charge() == -1){
1208 fDecayParameters[0] = pTrackXi->GetTPCNcls();
1209 fDecayParameters[1] = nTrackXi->GetTPCNcls();
1210 fDecayParameters[4] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1211 fDecayParameters[5] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
1213 fDecayParameters[0] = nTrackXi->GetTPCNcls();
1214 fDecayParameters[1] = pTrackXi->GetTPCNcls();
1215 fDecayParameters[4] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1216 fDecayParameters[5] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
1220 fDecayParameters[6] = fabs(bTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion second
1221 fDecayParameters[7] = fabs(Xicandidate->GetD(primaryVtx[0],primaryVtx[1],primaryVtx[2]));// DCA Vtx Lambda
1222 fDecayParameters[9] = fabs(Xicandidate->GetDcaV0Daughters());// DCA proton-pion
1223 fDecayParameters[10] = fabs(Xicandidate->GetDcaXiDaughters());// DCA Lambda-pion
1227 Double_t tempX[3]={0};
1228 Xicandidate->GetXYZ(tempX[0], tempX[1], tempX[2]);
1229 fDecayParameters[11] = sqrt( pow(tempX[0],2) + pow(tempX[1],2));// Rxy Lambda
1230 if(sqrt( pow(tempX[0],2) + pow(tempX[1],2) ) > fMaxDecayLength) continue;
1233 fDecayParameters[13] = Xicandidate->GetV0CosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Lambda
1234 fDecayParameters[14] = Xicandidate->GetCascadeCosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Xi
1236 xiP[0] = Xicandidate->Px();
1237 xiP[1] = Xicandidate->Py();
1238 xiP[2] = Xicandidate->Pz();
1239 xiVtx[0] = Xicandidate->Xv();
1240 xiVtx[1] = Xicandidate->Yv();
1241 xiVtx[2] = Xicandidate->Zv();
1242 xiPt = Xicandidate->Pt();
1243 xiY = Xicandidate->RapXi();
1244 xiMass = Xicandidate->M();
1245 xiCharge = Xicandidate->Charge();
1247 decayLengthXY = sqrt( pow(xiVtx[0]-primaryVtx[0],2) + pow(xiVtx[1]-primaryVtx[1],2) );
1248 fDecayParameters[12] = decayLengthXY;// Rxy Xi
1249 if(decayLengthXY > fMaxDecayLength) continue;// 2d version
1251 Bool_t StandardXi=kTRUE;
1252 if(fDecayParameters[0] < fCutValues[0][0]) StandardXi=kFALSE;// Nclus proton
1253 if(fDecayParameters[1] < fCutValues[0][1]) StandardXi=kFALSE;// Nclus pion first
1254 if(fDecayParameters[2] < fCutValues[0][2]) StandardXi=kFALSE;// Nclus pion second
1256 if(fDecayParameters[4] < fCutValues[0][4]) StandardXi=kFALSE;// DCAVtx proton
1257 if(fDecayParameters[5] < fCutValues[0][5]) StandardXi=kFALSE;// DCAVtx pion first
1258 if(fDecayParameters[6] < fCutValues[0][6]) StandardXi=kFALSE;// DCAVtx pion second
1259 if(fDecayParameters[7] < fCutValues[0][7]) StandardXi=kFALSE;// DCAVtx Lambda
1261 if(fDecayParameters[9] > fCutValues[0][9]) StandardXi=kFALSE;// DCAV proton-pion
1262 if(fDecayParameters[10] > fCutValues[0][10]) StandardXi=kFALSE;// DCAV Lambda-pion
1264 if(fDecayParameters[11] < fCutValues[0][11]) StandardXi=kFALSE;// Rxy Lambda
1265 if(fDecayParameters[12] < fCutValues[0][12]) StandardXi=kFALSE;// Rxy Xi
1267 if(fDecayParameters[13] < fCutValues[0][13]) StandardXi=kFALSE;// Cos PA Lambda
1268 if(fDecayParameters[14] < fCutValues[0][14]) StandardXi=kFALSE;// Cos PA Xi
1271 if(xiCharge == -1) CutVar[0].fXi->Fill(xiPt, xiY, xiMass);
1272 else CutVar[0].fXibar->Fill(xiPt, xiY, xiMass);
1276 mcXiFilled = kFALSE;
1277 if(fMCcase && !fAODcase){
1279 MCXiD2esd = (TParticle*)mcstack->Particle(abs(bTrackXi->GetLabel()));
1281 if(abs(MCXiD2esd->GetPdgCode())==kPionCode){
1283 MCLamD1esd = (TParticle*)mcstack->Particle(abs(pTrackXi->GetLabel()));
1284 MCLamD2esd = (TParticle*)mcstack->Particle(abs(nTrackXi->GetLabel()));
1286 if(MCLamD1esd->GetMother(0) == MCLamD2esd->GetMother(0)){
1287 if(abs(MCLamD1esd->GetPdgCode())==kProtonCode || abs(MCLamD2esd->GetPdgCode())==kProtonCode) {
1288 if(abs(MCLamD1esd->GetPdgCode())==kPionCode || abs(MCLamD2esd->GetPdgCode())==kPionCode) {
1290 MCLamesd = (TParticle*)mcstack->Particle(abs(MCLamD1esd->GetMother(0)));
1291 if(abs(MCLamesd->GetPdgCode())==kLambdaCode) {
1293 if(MCLamesd->GetMother(0) == MCXiD2esd->GetMother(0)){
1294 MCXiesd = (TParticle*)mcstack->Particle(abs(MCLamesd->GetMother(0)));
1295 if(abs(MCXiesd->GetPdgCode())==kXiCode) {
1299 if(Xicandidate->Charge() == -1) {
1300 CutVar[0].fMCrecXi->Fill(xiPt, xiY, xiMass);
1302 CutVar[0].fMCrecXibar->Fill(xiPt, xiY, xiMass);
1316 if(fabs(xiMass-fTrueMassXi) > fMassWindow) continue;
1320 fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
1323 //////////////////////////////////////////////////////////
1324 // Reconstruct Xi(1530)
1325 for(Int_t EN=0; EN<fEventsToMix+1; EN++){// Event buffer loop
1327 for(Int_t l=0; l<(fEvt+EN)->fNTracks; l++){// Present(EN=0) and Past(EN from 1 to fEventsToMix) event track loop
1330 if((fEvt+EN)->fTracks[l].fID == pTrackXi->GetID()) continue;
1331 if((fEvt+EN)->fTracks[l].fID == nTrackXi->GetID()) continue;
1332 if((fEvt+EN)->fTracks[l].fID == bTrackXi->GetID()) continue;
1335 fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
1337 if(!fESDTrack4) continue;
1338 fESDTrack4->Set((fEvt+EN)->fTracks[l].fX, (fEvt+EN)->fTracks[l].fP, (fEvt+EN)->fTracks[l].fCov, (fEvt+EN)->fTracks[l].fCharge);
1340 if((Bool_t)(((1<<5) & (fEvt+EN)->fTracks[l].fFilterMap) == 0)) continue;// AOD filterbit cut, "Standard cuts with tight dca"
1342 fDecayParameters[8] = (fEvt+EN)->fTracks[l].fDCAXY;// DCA Vtx pion third
1343 if((fEvt+EN)->fTracks[l].fDCAZ > 2) continue;
1344 if( (((fEvt+EN)->fTracks[l].fStatus)&AliESDtrack::kITSrefit)==0) continue;// Require itsrefit
1345 // no Chi^2 cut applied for ESDs. Info not available in my track structure.
1348 if(fabs((fEvt+EN)->fTracks[l].fEta) > 0.8) continue;
1350 fDecayParameters[3] = (fEvt+EN)->fTracks[l].fNclusTPC;
1352 AliVertex *XiStarVtx = new AliVertex((fEvt+EN)->fTracks[l].fX,0,0);
1353 //fESDTrack4->PropagateToDCA(fXiTrack, bField);// Propagate tracks to dca, both tracks are budged
1354 if(!(fXiTrack->PropagateToDCA(XiStarVtx, bField, 3))) continue;// Propagate tracks to dca, version which assumes fESDTrack4 is already primary
1356 fXiTrack->GetPxPyPz(pDaughter1);
1357 fXiTrack->GetXYZ(xDaughter1);
1358 fESDTrack4->GetPxPyPz(pDaughter2);
1359 fESDTrack4->GetXYZ(xDaughter2);
1360 //////////////////////////
1364 //xiStarVtx[0] = (xDaughter1[0]+xDaughter2[0])/2.;
1365 //xiStarVtx[1] = (xDaughter1[1]+xDaughter2[1])/2.;
1366 //xiStarVtx[2] = (xDaughter1[2]+xDaughter2[2])/2.;
1367 //decayLength = sqrt(pow(xiStarVtx[0]-primaryVtx[0],2)+pow(xiStarVtx[1]-primaryVtx[1],2)+pow(xiStarVtx[2]-primaryVtx[2],2));
1376 p1sq=px1*px1+py1*py1+pz1*pz1;
1377 p2sq=px2*px2+py2*py2+pz2*pz2;
1378 if(p1sq <=0 || p2sq <=0) continue;
1380 e1=sqrt(p1sq+fTrueMassXi*fTrueMassXi);
1381 e2=sqrt(p2sq+fTrueMassPi*fTrueMassPi);
1382 angle=px1*px2+py1*py2+pz1*pz2;
1383 xiStarMass=fTrueMassXi*fTrueMassXi+fTrueMassPi*fTrueMassPi+2.*e1*e2-2.*angle;
1384 if(xiStarMass<0.) xiStarMass=1.e-8;
1385 xiStarMass=sqrt(xiStarMass);
1388 xiStarP[0] = px1+px2;
1389 xiStarP[1] = py1+py2;
1390 xiStarP[2] = pz1+pz2;
1391 xiStarMom = sqrt(pow(xiStarP[0],2)+pow(xiStarP[1],2)+pow(xiStarP[2],2));
1392 if(xiStarMom==0) continue; // So one of the following lines doesnt break
1393 xiStarPt = sqrt(xiStarP[0]*xiStarP[0] + xiStarP[1]*xiStarP[1]);
1394 xiStarY = .5*log( ((e1+e2) + xiStarP[2])/((e1+e2) - xiStarP[2]));
1395 //xiStarE = e1 + e2;
1398 //if( (xiStarP[0]*(xiStarVtx[0]-primaryVtx[0]) + xiStarP[1]*(xiStarVtx[1]-primaryVtx[1]) + xiStarP[2]*(xiStarVtx[2]-primaryVtx[2]))/xiStarMom/decayLength < fXiStarCosTheta) continue;
1400 for(int cv=0; cv<kNCutVariations; cv++){
1401 if(fDecayParameters[0] < fCutValues[cv][0]) continue;// Nclus proton
1402 if(fDecayParameters[1] < fCutValues[cv][1]) continue;// Nclus pion first
1403 if(fDecayParameters[2] < fCutValues[cv][2]) continue;// Nclus pion second
1404 if(fDecayParameters[3] < fCutValues[cv][3]) continue;// Nclus pion third
1406 if(fDecayParameters[4] < fCutValues[cv][4]) continue;// DCAVtx proton
1407 if(fDecayParameters[5] < fCutValues[cv][5]) continue;// DCAVtx pion first
1408 if(fDecayParameters[6] < fCutValues[cv][6]) continue;// DCAVtx pion second
1409 if(fDecayParameters[7] < fCutValues[cv][7]) continue;// DCAVtx Lambda
1410 if(cv!=8) {if(fDecayParameters[8] > (0.0182 + 0.035/pow((fEvt+EN)->fTracks[l].fPt,1.01))) continue;}// DCAVtx pion third
1411 else {if(fDecayParameters[8] > fCutValues[cv][8]) continue;}// DCAVtx pion third
1413 if(fDecayParameters[9] > fCutValues[cv][9]) continue;// DCAV proton-pion
1414 if(fDecayParameters[10] > fCutValues[cv][10]) continue;// DCAV Lambda-pion
1416 if(fDecayParameters[11] < fCutValues[cv][11]) continue;// Rxy Lambda
1417 if(fDecayParameters[12] < fCutValues[cv][12]) continue;// Rxy Xi
1419 if(fDecayParameters[13] < fCutValues[cv][13]) continue;// Cos PA Lambda
1420 if(fDecayParameters[14] < fCutValues[cv][14]) continue;// Cos PA Xi
1424 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1425 else if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1426 else if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1427 else CutVar[cv].fXiPlusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1429 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1430 else if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1431 else if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1432 else CutVar[cv].fXiPlusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1437 // MC associaton AOD
1438 if(fMCcase && mcXiFilled && EN==0 && fAODcase){// AOD MC's
1440 MCXiStarD2 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[l].fLabel));
1442 if(abs(MCXiStarD2->GetPdgCode())==kPionCode){
1443 if(MCXi->GetMother() == MCXiStarD2->GetMother()){
1444 MCXiStar = (AliAODMCParticle*)mcArray->At(MCXi->GetMother());
1445 if(abs(MCXiStar->GetPdgCode())==kXiStarCode) {
1447 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1448 if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1456 // MC associaton ESD
1457 if(fMCcase && mcXiFilled && EN==0 && !fAODcase){// ESD MC's
1458 MCXiStarD2esd = (TParticle*)mcstack->Particle(abs((fEvt)->fTracks[l].fLabel));
1460 if(abs(MCXiStarD2esd->GetPdgCode())==kPionCode){
1461 if(MCXiesd->GetMother(0) == MCXiStarD2esd->GetMother(0)){
1463 MCXiStaresd = (TParticle*)mcstack->Particle(abs(MCXiesd->GetMother(0)));
1464 if(abs(MCXiStaresd->GetPdgCode())==kXiStarCode) {
1466 if(fXiTrack->Charge() == -1 && fESDTrack4->Charge() == +1) CutVar[cv].fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1467 if(fXiTrack->Charge() == +1 && fESDTrack4->Charge() == -1) CutVar[cv].fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1474 }// Cut Variation loop
1476 }// Event mixing loop
1482 // Post output data.
1483 PostData(1, fOutputList);
1486 //________________________________________________________________________
1487 void AliXiStar::Terminate(Option_t *)
1491 //________________________________________________________________________
1492 Double_t AliXiStar::LinearPropagateToDCA(AliESDtrack *v, AliESDtrack *t, Double_t b) {// Adapted from AliCascadeVertexer.cxx
1493 //--------------------------------------------------------------------
1494 // This function returns the DCA between the V0 and the track
1495 //--------------------------------------------------------------------
1497 Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
1498 Double_t r[3]; t->GetXYZ(r);
1499 Double_t x1=r[0], y1=r[1], z1=r[2];
1500 Double_t p[3]; t->GetPxPyPz(p);
1501 Double_t px1=p[0], py1=p[1], pz1=p[2];
1505 Double_t vx2,vy2,vz2; // position and momentum of V0
1506 Double_t px2,py2,pz2;
1510 vx2=x2[0], vy2=x2[1], vz2=x2[2];
1511 px2=p2[0], py2=p2[1], pz2=p2[2];
1515 Double_t dd= Det(vx2-x1,vy2-y1,vz2-z1,px1,py1,pz1,px2,py2,pz2);
1516 Double_t ax= Det(py1,pz1,py2,pz2);
1517 Double_t ay=-Det(px1,pz1,px2,pz2);
1518 Double_t az= Det(px1,py1,px2,py2);
1520 Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
1523 Double_t t1 = Det(vx2-x1,vy2-y1,vz2-z1,px2,py2,pz2,ax,ay,az)/
1524 Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
1526 x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1;
1529 //propagate track to the points of DCA
1532 if (!t->PropagateTo(x1,b)) {
1533 Error("PropagateToDCA","Propagation failed !");
1541 //________________________________________________________________________
1542 Double_t AliXiStar::Det(Double_t a00, Double_t a01, Double_t a10, Double_t a11) const {// Taken from AliCascadeVertexer
1543 //--------------------------------------------------------------------
1544 // This function calculates locally a 2x2 determinant
1545 //--------------------------------------------------------------------
1546 return a00*a11 - a01*a10;
1548 //________________________________________________________________________
1549 Double_t AliXiStar::Det(Double_t a00,Double_t a01,Double_t a02,
1550 Double_t a10,Double_t a11,Double_t a12,
1551 Double_t a20,Double_t a21,Double_t a22) const {// Taken from AliCascadeVertexer
1552 //--------------------------------------------------------------------
1553 // This function calculates locally a 3x3 determinant
1554 //--------------------------------------------------------------------
1555 return a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);